Extract inputs for TP from the state of temporal memory
    three information are extracted
    1. correctly predicted cells
    2. all active cells
    3. bursting cells (unpredicted input)
320331
def extractInputForTP(self, tm):
    

    # bursting cells in layer 4
    burstingColumns = tm.activeState["t"].sum(axis=1)
    burstingColumns[ burstingColumns < tm.cellsPerColumn ] = 0
    burstingColumns[ burstingColumns == tm.cellsPerColumn ] = 1
    # print "Bursting column indices=",burstingColumns.nonzero()[0]

    # correctly predicted cells in layer 4
    correctlyPredictedCells = numpy.zeros(self._inputDimensions).astype(realDType)
    idx = (tm.predictedState["t-1"] + tm.activeState["t"]) == 2
    idx = idx.reshape(self._inputDimensions)
    correctlyPredictedCells[idx] = 1.0
    # print "Predicted->active cell indices=",correctlyPredictedCells.nonzero()[0]

    # all currently active cells in layer 4
    spInputVector = tm.learnState["t"].reshape(self._inputDimensions)
    # spInputVector = tm.activeState["t"].reshape(self._inputDimensions)

    return (correctlyPredictedCells, spInputVector, burstingColumns)
883485
def _learnTransition(self, prevActiveCells, deltaLocation, newLocation):
    

    prevLocationPotentialOverlaps = self.internalConnections.computeActivity(
      prevActiveCells)
    deltaPotentialOverlaps = self.deltaConnections.computeActivity(
      deltaLocation)

    matchingDeltaSegments = np.where(
      (prevLocationPotentialOverlaps >= self.learningThreshold) &
      (deltaPotentialOverlaps >= self.learningThreshold))[0]

    # Cells with a active segment pair: reinforce the segment
    cellsForActiveSegments = self.internalConnections.mapSegmentsToCells(
      self.activeDeltaSegments)
    learningActiveDeltaSegments = self.activeDeltaSegments[
      np.in1d(cellsForActiveSegments, newLocation)]
    remainingCells = np.setdiff1d(newLocation, cellsForActiveSegments)

    # Remaining cells with a matching segment pair: reinforce the best matching
    # segment pair.
    candidateSegments = self.internalConnections.filterSegmentsByCell(
      matchingDeltaSegments, remainingCells)
    cellsForCandidateSegments = self.internalConnections.mapSegmentsToCells(
      candidateSegments)
    candidateSegments = matchingDeltaSegments[
      np.in1d(cellsForCandidateSegments, remainingCells)]
    onePerCellFilter = np2.argmaxMulti(
      prevLocationPotentialOverlaps[candidateSegments] +
      deltaPotentialOverlaps[candidateSegments],
      cellsForCandidateSegments)
    learningMatchingDeltaSegments = candidateSegments[onePerCellFilter]

    newDeltaSegmentCells = np.setdiff1d(remainingCells, cellsForCandidateSegments)

    for learningSegments in (learningActiveDeltaSegments,
                             learningMatchingDeltaSegments):
      self._learn(self.internalConnections, self.rng, learningSegments,
                  prevActiveCells, prevActiveCells,
                  prevLocationPotentialOverlaps,
                  self.initialPermanence, self.sampleSize,
                  self.permanenceIncrement, self.permanenceDecrement,
                  self.maxSynapsesPerSegment)
      self._learn(self.deltaConnections, self.rng, learningSegments,
                  deltaLocation, deltaLocation, deltaPotentialOverlaps,
                  self.initialPermanence, self.sampleSize,
                  self.permanenceIncrement, self.permanenceDecrement,
                  self.maxSynapsesPerSegment)

    numNewLocationSynapses = len(prevActiveCells)
    numNewDeltaSynapses = len(deltaLocation)

    if self.sampleSize != -1:
      numNewLocationSynapses = min(numNewLocationSynapses, self.sampleSize)
      numNewDeltaSynapses = min(numNewDeltaSynapses, self.sampleSize)

    if self.maxSynapsesPerSegment != -1:
      numNewLocationSynapses = min(numNewLocationSynapses,
                                   self.maxSynapsesPerSegment)
      numNewDeltaSynapses = min(numNewLocationSynapses,
                                self.maxSynapsesPerSegment)

    newPrevLocationSegments = self.internalConnections.createSegments(
      newDeltaSegmentCells)
    newDeltaSegments = self.deltaConnections.createSegments(
      newDeltaSegmentCells)

    assert np.array_equal(newPrevLocationSegments, newDeltaSegments)

    self.internalConnections.growSynapsesToSample(
      newPrevLocationSegments, prevActiveCells, numNewLocationSynapses,
      self.initialPermanence, self.rng)
    self.deltaConnections.growSynapsesToSample(
      newDeltaSegments, deltaLocation, numNewDeltaSynapses,
      self.initialPermanence, self.rng)
1060347
def _extract_pynn_components_to_neuroml(nl_model, nml_doc=None):
    
    
    if nml_doc == None:
        from neuroml import NeuroMLDocument
        nml_doc = NeuroMLDocument(id="temp")
    
    for c in nl_model.cells:
        
        if c.pynn_cell:
            
            if nml_doc.get_by_id(c.id) == None:
                import pyNN.neuroml
                cell_params = c.parameters if c.parameters else {}
                #print('------- %s: %s' % (c, cell_params))
                for p in cell_params:
                    cell_params[p] = evaluate(cell_params[p], nl_model.parameters)
                #print('====== %s: %s' % (c, cell_params))
                for proj in nl_model.projections:

                    synapse = nl_model.get_child(proj.synapse, 'synapses')
                    post_pop = nl_model.get_child(proj.postsynaptic, 'populations')
                    if post_pop.component == c.id:
                        #print("--------- Cell %s in post pop %s of %s uses %s"%(c.id,post_pop.id, proj.id, synapse))

                        if synapse.pynn_receptor_type == 'excitatory':
                            post = '_E'
                        elif synapse.pynn_receptor_type == 'inhibitory':
                            post = '_I'
                        for p in synapse.parameters:
                            cell_params['%s%s' % (p, post)] = synapse.parameters[p]

                temp_cell = eval('pyNN.neuroml.%s(**cell_params)' % c.pynn_cell)
                if c.pynn_cell != 'SpikeSourcePoisson':
                    temp_cell.default_initial_values['v'] = temp_cell.parameter_space['v_rest'].base_value

                cell_id = temp_cell.add_to_nml_doc(nml_doc, None)
                cell = nml_doc.get_by_id(cell_id)
                cell.id = c.id
            
    for s in nl_model.synapses:
        if nml_doc.get_by_id(s.id) == None:
            
            if s.pynn_synapse_type and s.pynn_receptor_type:
                import neuroml
                
                if s.pynn_synapse_type == 'cond_exp':
                    syn = neuroml.ExpCondSynapse(id=s.id, tau_syn=s.parameters['tau_syn'], e_rev=s.parameters['e_rev'])
                    nml_doc.exp_cond_synapses.append(syn)
                elif s.pynn_synapse_type == 'cond_alpha':
                    syn = neuroml.AlphaCondSynapse(id=s.id, tau_syn=s.parameters['tau_syn'], e_rev=s.parameters['e_rev'])
                    nml_doc.alpha_cond_synapses.append(syn)
                elif s.pynn_synapse_type == 'curr_exp':
                    syn = neuroml.ExpCurrSynapse(id=s.id, tau_syn=s.parameters['tau_syn'])
                    nml_doc.exp_curr_synapses.append(syn)
                elif s.pynn_synapse_type == 'curr_alpha':
                    syn = neuroml.AlphaCurrSynapse(id=s.id, tau_syn=s.parameters['tau_syn'])
                    nml_doc.alpha_curr_synapses.append(syn)
    
    for i in nl_model.input_sources:
        
        #if nml_doc.get_by_id(i.id) == None:
      
        if i.pynn_input:
            import pyNN.neuroml
            input_params = i.parameters if i.parameters else {}
            exec('input__%s = pyNN.neuroml.%s(**input_params)' % (i.id, i.pynn_input))
            exec('temp_input = input__%s' % i.id)
            pg_id = temp_input.add_to_nml_doc(nml_doc, None)
            #for pp in nml_doc.pulse_generators:
            #    print('PG: %s: %s'%(pp,pp.id))
            pg = nml_doc.get_by_id(pg_id)
            pg.id = i.id
            
    return nml_doc
1296679
def reset(self):
    
    self.bumpPhases = np.empty((2,0), dtype="float")
    self.phaseDisplacement = np.empty((0,2), dtype="float")
    self.cellsForActivePhases = np.empty(0, dtype="int")
    self.activeCells = np.empty(0, dtype="int")
    self.learningCells = np.empty(0, dtype="int")
    self.sensoryAssociatedCells = np.empty(0, dtype="int")
546460
def _preSynapticTRNCells(self, i, j):
    
    xmin = max(i - 1, 0)
    xmax = min(i + 2, self.trnWidth)
    ymin = max(j - 1, 0)
    ymax = min(j + 2, self.trnHeight)
    trnCells = [
      (x, y) for x in range(xmin, xmax) for y in range(ymin, ymax)
    ]

    return trnCells
758808
def _initializeTRNToRelayCellConnections(self):
    
    for x in range(self.relayWidth):
      for y in range(self.relayHeight):

        # Create one dendrite for each trn cell that projects to this relay cell
        # This dendrite contains one synapse corresponding to this TRN->relay
        # connection.
        relayCellIndex = self.relayCellIndex((x,y))
        trnCells = self._preSynapticTRNCells(x, y)
        for trnCell in trnCells:
          newSegment = self.relayConnections.createSegments([relayCellIndex])
          self.relayConnections.growSynapses(newSegment,
                                             [self.trnCellIndex(trnCell)], 1.0)
1061900
def getActiveSegment(self, c, i, timeStep):
    

    # todo: put back preference for sequence segments.

    nSegments = len(self.cells[c][i])

    bestActivation = self.activationThreshold
    which = -1

    for j,s in enumerate(self.cells[c][i]):

      activity = self.getSegmentActivityLevel(s, self.activeState[timeStep], connectedSynapsesOnly = True)

      if activity >= bestActivation:
        bestActivation = activity
        which = j

    if which != -1:
      return self.cells[c][i][which]
    else:
      return None
133807
def compute(self, inputs, outputs):
    

    # If there's a reset, don't call compute. In some implementations, an empty
    # input might cause unwanted effects.
    if "resetIn" in inputs:
      assert len(inputs["resetIn"]) == 1
      if inputs["resetIn"][0] != 0:
        # send empty output
        self._tm.reset()
        outputs["activeCells"][:] = 0
        outputs["predictedActiveCells"][:] = 0
        outputs["winnerCells"][:] = 0
        return

    activeColumns = inputs["activeColumns"].nonzero()[0]

    if "basalInput" in inputs:
      basalInput = inputs["basalInput"].nonzero()[0]
    else:
      basalInput = np.empty(0, dtype="uint32")

    if "apicalInput" in inputs:
      apicalInput = inputs["apicalInput"].nonzero()[0]
    else:
      apicalInput = np.empty(0, dtype="uint32")

    if "basalGrowthCandidates" in inputs:
      basalGrowthCandidates = inputs["basalGrowthCandidates"].nonzero()[0]
    else:
      basalGrowthCandidates = basalInput

    if "apicalGrowthCandidates" in inputs:
      apicalGrowthCandidates = inputs["apicalGrowthCandidates"].nonzero()[0]
    else:
      apicalGrowthCandidates = apicalInput

    self._tm.compute(activeColumns, basalInput, apicalInput,
                     basalGrowthCandidates, apicalGrowthCandidates, self.learn)

    # Extract the active / predicted cells and put them into binary arrays.
    outputs["activeCells"][:] = 0
    outputs["activeCells"][self._tm.getActiveCells()] = 1
    outputs["predictedCells"][:] = 0
    outputs["predictedCells"][
      self._tm.getPredictedCells()] = 1
    outputs["predictedActiveCells"][:] = (outputs["activeCells"] *
                                          outputs["predictedCells"])
    outputs["winnerCells"][:] = 0
    outputs["winnerCells"][self._tm.getWinnerCells()] = 1
453948
def getDefaultTMParams(self, inputSize, numInputBits):
    
    sampleSize = int(1.5 * numInputBits)

    if numInputBits == 20:
      activationThreshold = 18
      minThreshold = 18
    elif numInputBits == 10:
      activationThreshold = 8
      minThreshold = 8
    else:
      activationThreshold = int(numInputBits * .6)
      minThreshold = activationThreshold

    return {
      "columnCount": inputSize,
      "cellsPerColumn": 16,
      "learn": True,
      "learnOnOneCell": False,
      "initialPermanence": 0.41,
      "connectedPermanence": 0.6,
      "permanenceIncrement": 0.1,
      "permanenceDecrement": 0.03,
      "minThreshold": minThreshold,
      "basalPredictedSegmentDecrement": 0.003,
      "apicalPredictedSegmentDecrement": 0.0,
      "reducedBasalThreshold": int(activationThreshold*0.6),
      "activationThreshold": activationThreshold,
      "sampleSize": sampleSize,
      "implementation": "ApicalTiebreak",
      "seed": self.seed
    }
16343
def predictor(enc_flat,
              action,
              lstm_states,
              pred_depth,
              reuse=False,
              scope_prefix='',
              hparams=None):
  """"""
  with tf.variable_scope(scope_prefix + 'predict', reuse=reuse):

    enc_final_size = enc_flat.get_shape().as_list()[1]
    action_size = action.get_shape().as_list()[1]
    initial_size = (enc_final_size + action_size)

    batch_size = tf.shape(enc_flat)[0]

    init_stddev = 1e-2

    pre_pred = tf.concat([enc_flat, action], 1)
    pre_pred = tf.layers.dense(
        pre_pred,
        initial_size,
        kernel_initializer=tf.truncated_normal_initializer(stddev=init_stddev))

    # This is only needed or the GAN version.
    if hparams.pred_noise_std > 0:
      # Add the noise like this so a pretrained model can be used.
      pred_noise = tf.random_normal(
          shape=[batch_size, 100], stddev=hparams.pred_noise_std)
      pre_pred += tf.layers.dense(
          pred_noise,
          initial_size,
          kernel_initializer=tf.truncated_normal_initializer(
              stddev=init_stddev),
          name='noise_dense')

    pre_pred = tf.nn.relu(pre_pred)

    if lstm_states[pred_depth - 2] is None:
      back_connect = tf.tile(
          tf.get_variable(
              'back_connect_init',
              shape=[1, initial_size * 2],
              initializer=tf.truncated_normal_initializer(stddev=init_stddev))
          , (batch_size, 1))
    else:
      back_connect = lstm_states[pred_depth - 2]

    lstm_init_stddev = 1e-4

    part_pred, lstm_states[0] = common_video.lstm_cell(
        tf.concat([pre_pred, back_connect], 1),
        lstm_states[0],
        initial_size,
        use_peepholes=True,
        initializer=tf.truncated_normal_initializer(stddev=lstm_init_stddev),
        num_proj=initial_size)
    part_pred = tf.contrib.layers.layer_norm(part_pred)
    pred = part_pred

    for pred_layer_num in range(1, pred_depth, 2):
      part_pred, lstm_states[pred_layer_num] = common_video.lstm_cell(
          pred,
          lstm_states[pred_layer_num],
          initial_size,
          use_peepholes=True,
          initializer=tf.truncated_normal_initializer(stddev=lstm_init_stddev),
          num_proj=initial_size)
      pred += part_pred

      part_pred, lstm_states[pred_layer_num + 1] = common_video.lstm_cell(
          tf.concat([pred, pre_pred], 1),
          lstm_states[pred_layer_num + 1],
          initial_size,
          use_peepholes=True,
          initializer=tf.truncated_normal_initializer(stddev=lstm_init_stddev),
          num_proj=initial_size)
      part_pred = tf.contrib.layers.layer_norm(part_pred)
      pred += part_pred

    pred = tf.layers.dense(
        pred,
        enc_final_size,
        kernel_initializer=tf.truncated_normal_initializer(stddev=init_stddev))

    if hparams.enc_pred_use_l2norm:
      pred = tf.nn.l2_normalize(pred, 1)

    return pred
210220
def getUnionTemporalPoolerInput(self):
    
    activeCells = numpy.zeros(self.tm.numberOfCells()).astype(realDType)
    activeCells[list(self.tm.activeCellsIndices())] = 1

    predictedActiveCells = numpy.zeros(self.tm.numberOfCells()).astype(
      realDType)
    predictedActiveCells[list(self.tm.predictedActiveCellsIndices())] = 1

    burstingColumns = numpy.zeros(self.tm.numberOfColumns()).astype(realDType)
    burstingColumns[list(self.tm.unpredictedActiveColumns)] = 1

    return activeCells, predictedActiveCells, burstingColumns
96893
def filtered(w=250):
  
  # prepare filter bank kernels
  kernels = []
  for theta in range(4):
    theta = theta / 4. * np.pi
    for sigma in (1, 3):
      for frequency in (0.05, 0.25):
        kernel = np.real(gabor_kernel(frequency, theta=theta,
                                      sigma_x=sigma, sigma_y=sigma))
        kernels.append(kernel)

  print("Initializing thalamus")

  t = Thalamus(
    trnCellShape=(w, w),
    relayCellShape=(w, w),
    inputShape=(w, w),
    l6CellCount=128*128,
    trnThreshold=15,
  )


  ff = loadImage(t)

  for i,k in enumerate(kernels):
    plotActivity(k, "kernel"+str(i)+".jpg", "Filter kernel", vmax=k.max(),
                 vmin=k.min())
    filtered0 = power(ff, k)
    ft = np.zeros((w, w))
    ft[filtered0 > filtered0.mean() + filtered0.std()] = 1.0
    plotActivity(ft, "filtered"+str(i)+".jpg", "Filtered image", vmax=1.0)

  encoder = createLocationEncoder(t, w=17)
  trainThalamusLocations(t, encoder)

  filtered0 = power(ff, kernels[3])
  ft = np.zeros((w, w))
  ft[filtered0 > filtered0.mean() + filtered0.std()] = 1.0

  # Get a salt and pepper burst ready image
  print("Getting unions")
  l6Code = list(getUnionLocations(encoder, 125, 125, 150, step=10))
  print("Num active cells in L6 union:", len(l6Code),"out of", t.l6CellCount)
  ffOutput = inferThalamus(t, l6Code, ft)
  plotActivity(t.burstReadyCells, "relay_burstReady_filtered.jpg",
               title="Burst-ready cells",
               )
  plotActivity(ffOutput, "cajal_relay_output_filtered.jpg",
               title="Filtered activity",
               cmap="Greys")

  # Get a more detailed filtered image
  print("Getting unions")
  l6Code = list(getUnionLocations(encoder, 125, 125, 150, step=3))
  print("Num active cells in L6 union:", len(l6Code),"out of", t.l6CellCount)
  ffOutput_all = inferThalamus(t, l6Code, ff)
  ffOutput_filtered = inferThalamus(t, l6Code, ft)
  ffOutput3 = ffOutput_all*0.4 + ffOutput_filtered
  plotActivity(t.burstReadyCells, "relay_burstReady_all.jpg",
               title="Burst-ready cells",
               )
  plotActivity(ffOutput3, "cajal_relay_output_filtered2.jpg",
               title="Filtered activity",
               cmap="Greys")
611428
def getL4PredictedActiveCells(self):
    
    predictedActive = []
    for i in xrange(self.numColumns):
      region = self.network.regions["L4Column_" + str(i)]
      predictedActive.append(
        region.getOutputData("predictedActiveCells").nonzero()[0])
    return predictedActive
722172
def getMeanInpCurrents(params, numunits=100,
                       filepattern=os.path.join('simulation_output_default',
                                                'population_input_spikes*')):
    '''
    
    '''
    #convolution kernels
    x = np.arange(100) * params.dt
    kernel = np.exp(-x / params.model_params['tau_syn_ex'])


    #number of external inputs:
    K_bg = np.array(sum(params.K_bg, []))

    #compensate for DC CC connections if we're using that
    iDC =  K_bg * params.dc_amplitude * 1E-3 # unit ????
    
    data = {}
    
    #loop over network-populations
    for i, Y in enumerate(params.Y):
        if i % SIZE == RANK:
            #file to open
            fname = glob.glob(filepattern+'*' + Y + '*')[0]
            print fname
            #read in read data and assess units, up to numunits
            rawdata = np.array(helpers.read_gdf(fname))
            units = np.unique(rawdata[:, 0])
            if numunits > units.size:
                numcells = units.size
            else:
                numcells = numunits
            units = units[:numcells]
            
            #churn through data and extract the input currents per cell
            for j, unit in enumerate(units):
                slc = rawdata[:, 0] == unit
                #just the spikes:
                if j == 0:
                    dataslc = rawdata[slc, 2:]
                else:
                    dataslc = np.r_['0,3', dataslc, rawdata[slc, 2:]]
            
            #fix the datatype, it may be object
            dataslc = dataslc.astype(float)
                   
            #fill in data-structure
            data.update({
                Y : {
                    'E' : np.convolve(dataslc[:, :, 0].mean(axis=0),
                                      kernel, 'same')*1E-3 + float(iDC[i]),
                    'I' : np.convolve(dataslc[:, :, 1].mean(axis=0),
                                      kernel, 'same')*1E-3,
                    'tvec' : rawdata[slc, 1],
                    'numunits' : numunits,
                }
            })
    
    
    data = COMM.allgather(data)
    
    return {k: v for d in data for k, v in d.items()}
976868
def get_nonmonotonic_neurites(neuron, tol=1e-6):
    '''
    '''
    return [n for n in neuron.neurites if not is_monotonic(n, tol)]
1045169
def extract_data(data_path, feature):
    '''
    '''
    population = nm.load_neurons(data_path)

    feature_data = [nm.get(feature, n) for n in population]
    feature_data = list(chain(*feature_data))

    return stats.optimal_distribution(feature_data)
992782
def conv_lstm_2d(inputs, state, output_channels,
                 kernel_size=5, name=None, spatial_dims=None):
  """"""
  input_shape = common_layers.shape_list(inputs)
  batch_size, input_channels = input_shape[0], input_shape[-1]
  if spatial_dims is None:
    input_shape = input_shape[1:]
  else:
    input_shape = spatial_dims + [input_channels]

  cell = tf.contrib.rnn.ConvLSTMCell(
      2, input_shape, output_channels,
      [kernel_size, kernel_size], name=name)
  if state is None:
    state = cell.zero_state(batch_size, tf.float32)
  outputs, new_state = cell(inputs, state)
  return outputs, new_state