Extract stats from neurons
355198
def extract_stats(neurons, config):
    ''''''

    stats = defaultdict(dict)
    for ns, modes in config['neurite'].items():
        for n in config['neurite_type']:
            n = _NEURITE_MAP[n]
            for mode in modes:
                stat_name = _stat_name(ns, mode)
                stat = eval_stats(nm.get(ns, neurons, neurite_type=n), mode)

                if stat is None or not stat.shape:
                    stats[n.name][stat_name] = stat
                else:
                    assert stat.shape in ((3, ), ), \
                        'Statistic must create a 1x3 result'

                    for i, suffix in enumerate('XYZ'):
                        compound_stat_name = stat_name + '_' + suffix
                        stats[n.name][compound_stat_name] = stat[i]

    for ns, modes in config['neuron'].items():
        for mode in modes:
            stat_name = _stat_name(ns, mode)
            stats[stat_name] = eval_stats(nm.get(ns, neurons), mode)

    return stats
255662
def extract_images(filename):
    """"""
    with gzip.open(filename) as bytestream:
        magic = _read32(bytestream)
        if magic != 2051:
            raise ValueError(
                'Invalid magic number %d in MNIST image file: %s' %
                (magic, filename))
        num_images = _read32(bytestream)
        rows = _read32(bytestream)
        cols = _read32(bytestream)
        buf = bytestream.read(rows * cols * num_images)
        data = numpy.frombuffer(buf, dtype=numpy.uint8)
        data = data.reshape(num_images, rows, cols, 1)
        data = data.astype('float32') / 255.0
        return data
1217030
def mlp(feature, hparams, name="mlp"):
  """"""
  with tf.variable_scope(name, "mlp", values=[feature]):
    num_mlp_layers = hparams.num_mlp_layers
    mlp_dim = hparams.mlp_dim
    for _ in range(num_mlp_layers):
      feature = common_layers.dense(feature, mlp_dim, activation=tf.nn.relu)
      feature = tf.nn.dropout(feature, keep_prob=1.-hparams.dropout)
    return feature
425806
def calculate_reduce_output_shapes(operator):
    '''
    
    '''
    check_input_and_output_numbers(operator, input_count_range=1, output_count_range=1)
    check_input_and_output_types(operator, good_input_types=[FloatTensorType])

    output_shape = copy.deepcopy(operator.inputs[0].type.shape)
    params = operator.raw_operator.reduce

    from coremltools.proto.NeuralNetwork_pb2 import ReduceLayerParams as Params
    # Adjust C-axis
    if params.axis in [Params.CHW, Params.C]:
        output_shape[1] = 1
    # Adjust H-axis
    if params.axis in [Params.CHW, Params.HW, Params.H]:
        output_shape[2] = 1
    # Adjust W-axis
    if params.axis in [Params.CHW, Params.HW, Params.W]:
        output_shape[3] = 1

    operator.outputs[0].type.shape = output_shape
976868
def get_nonmonotonic_neurites(neuron, tol=1e-6):
    '''
    '''
    return [n for n in neuron.neurites if not is_monotonic(n, tol)]
315222
function() {
      var input_depth = this.data[0].w.length;
      var num_classes = this.unique_labels.length;

      // sample network topology and hyperparameters
      var layer_defs = [];
      layer_defs.push({type:'input', out_sx:1, out_sy:1, out_depth: input_depth});
      var nl = weightedSample([0,1,2,3], [0.2, 0.3, 0.3, 0.2]); // prefer nets with 1,2 hidden layers
      for(var q=0;q<nl;q++) {
        var ni = randi(this.neurons_min, this.neurons_max);
        var act = ['tanh','maxout','relu'][randi(0,3)];
        if(randf(0,1)<0.5) {
          var dp = Math.random();
          layer_defs.push({type:'fc', num_neurons: ni, activation: act, drop_prob: dp});
        } else {
          layer_defs.push({type:'fc', num_neurons: ni, activation: act});
        }
      }
      layer_defs.push({type:'softmax', num_classes: num_classes});
      var net = new Net();
      net.makeLayers(layer_defs);

      // sample training hyperparameters
      var bs = randi(this.batch_size_min, this.batch_size_max); // batch size
      var l2 = Math.pow(10, randf(this.l2_decay_min, this.l2_decay_max)); // l2 weight decay
      var lr = Math.pow(10, randf(this.learning_rate_min, this.learning_rate_max)); // learning rate
      var mom = randf(this.momentum_min, this.momentum_max); // momentum. Lets just use 0.9, works okay usually ;p
      var tp = randf(0,1); // trainer type
      var trainer_def;
      if(tp<0.33) {
        trainer_def = {method:'adadelta', batch_size:bs, l2_decay:l2};
      } else if(tp<0.66) {
        trainer_def = {method:'adagrad', learning_rate: lr, batch_size:bs, l2_decay:l2};
      } else {
        trainer_def = {method:'sgd', learning_rate: lr, momentum: mom, batch_size:bs, l2_decay:l2};
      }
      
      var trainer = new Trainer(net, trainer_def);

      var cand = {};
      cand.acc = [];
      cand.accv = 0; // this will maintained as sum(acc) for convenience
      cand.layer_defs = layer_defs;
      cand.trainer_def = trainer_def;
      cand.net = net;
      cand.trainer = trainer;
      return cand;
    }
225466
def _parse_statistic(data, scale):
        """"""
        i = 0
        for byte in bytearray(data):
            i = (i << 8) + byte

        if i == 32768:
            return None

        if scale == 0:
            return i

        return float(i) / (scale * 10)
743429
def get(self, spike_ids, channels=None):
        """"""
        if isinstance(spike_ids, slice):
            spike_ids = _range_from_slice(spike_ids,
                                          start=0,
                                          stop=self.n_spikes,
                                          )
        if not hasattr(spike_ids, '__len__'):
            spike_ids = [spike_ids]
        if channels is None:
            channels = slice(None, None, None)
            nc = self.n_channels
        else:
            channels = np.asarray(channels, dtype=np.int32)
            assert np.all(channels < self.n_channels)
            nc = len(channels)

        # Ensure a list of time samples are being requested.
        spike_ids = _as_array(spike_ids)
        n_spikes = len(spike_ids)

        # Initialize the array.
        # NOTE: last dimension is time to simplify things.
        shape = (n_spikes, nc, self._n_samples_extract)
        waveforms = np.zeros(shape, dtype=np.float32)

        # No traces: return null arrays.
        if self.n_samples_trace == 0:
            return np.transpose(waveforms, (0, 2, 1))

        # Load all spikes.
        for i, spike_id in enumerate(spike_ids):
            assert 0 <= spike_id < self.n_spikes
            time = self._spike_samples[spike_id]

            # Extract the waveforms on the unmasked channels.
            try:
                w = self._load_at(time, channels)
            except ValueError as e:  # pragma: no cover
                logger.warn("Error while loading waveform: %s", str(e))
                continue

            assert w.shape == (self._n_samples_extract, nc)

            waveforms[i, :, :] = w.T

        # Filter the waveforms.
        waveforms_f = waveforms.reshape((-1, self._n_samples_extract))
        # Only filter the non-zero waveforms.
        unmasked = waveforms_f.max(axis=1) != 0
        waveforms_f[unmasked] = self._filter(waveforms_f[unmasked], axis=1)
        waveforms_f = waveforms_f.reshape((n_spikes, nc,
                                           self._n_samples_extract))

        # Remove the margin.
        margin_before, margin_after = self._filter_margin
        if margin_after > 0:
            assert margin_before >= 0
            waveforms_f = waveforms_f[:, :, margin_before:-margin_after]

        assert waveforms_f.shape == (n_spikes,
                                     nc,
                                     self.n_samples_waveforms,
                                     )

        # NOTE: we transpose before returning the array.
        return np.transpose(waveforms_f, (0, 2, 1))
566310
def initialize_weights_nn(data, means, lognorm=True):
    
    # TODO
    genes, cells = data.shape
    k = means.shape[1]
    if lognorm:
        data = log1p(cell_normalize(data))
    for i in range(cells):
        for j in range(k):
            pass
317876
def _extract_mnist_images(filename, num_images):
  
  with gzip.open(filename) as bytestream:
    bytestream.read(16)
    buf = bytestream.read(_MNIST_IMAGE_SIZE * _MNIST_IMAGE_SIZE * num_images)
    data = np.frombuffer(buf, dtype=np.uint8)
    data = data.reshape(num_images, _MNIST_IMAGE_SIZE, _MNIST_IMAGE_SIZE, 1)
  return data
357684
def get_stats(self):
        """"""

        try:
            ret = fcntl.ioctl(self.ins, BIOCGSTATS, struct.pack("2I", 0, 0))
            return struct.unpack("2I", ret)
        except IOError:
            warning("Unable to get stats from BPF !")
            return (None, None)
1035153
public static void step(long seed, Neurons[] neurons, DeepLearningModel.DeepLearningModelInfo minfo, boolean training, double[] responses) {
    try {
      for (int i=1; i<neurons.length-1; ++i) {
        neurons[i].fprop(seed, training);
      }
      if (minfo.get_params().autoencoder) {
        neurons[neurons.length - 1].fprop(seed, training);
        if (training) {
          for (int i=neurons.length-1; i>0; --i) {
            neurons[i].bprop();
          }
        }
      } else {
        if (minfo.get_params().classification) {
          ((Neurons.Softmax) neurons[neurons.length - 1]).fprop();
          if (training) {
            for (int i = 1; i < neurons.length - 1; i++)
              Arrays.fill(neurons[i]._e.raw(), 0);
            int target_label;
            if (Double.isNaN(responses[0])) { //missing response
              target_label = Neurons.missing_int_value;
            } else {
              assert ((double) (int) responses[0] == responses[0]); //classification -> integer labels expected
              target_label = (int) responses[0];
            }
            ((Neurons.Softmax) neurons[neurons.length - 1]).bprop(target_label);
          }
        } else {
          ((Neurons.Linear) neurons[neurons.length - 1]).fprop();
          if (training) {
            for (int i = 1; i < neurons.length - 1; i++)
              Arrays.fill(neurons[i]._e.raw(), 0);
            float target_value;
            if (Double.isNaN(responses[0])) { //missing response
              target_value = Neurons.missing_real_value;
            } else {
              target_value = (float) responses[0];
            }
            ((Neurons.Linear) neurons[neurons.length - 1]).bprop(target_value);
          }
        }
        if (training) {
          for (int i=neurons.length-2; i>0; --i)
            neurons[i].bprop();
        }
      }
    }
    catch(RuntimeException ex) {
      Log.warn(ex.getMessage());
      minfo.set_unstable();
      throw new Job.JobCancelledException("Canceling job due to numerical instability.");
    }
  }
1005747
def get_stats(self):
    """"""
    ostr = ''
    errtotal = self.deletions['total']+self.insertions['total']+self.mismatches
    ostr += "ALIGNMENT_COUNT\t"+str(self.alignment_count)+"\n"
    ostr += "ALIGNMENT_BASES\t"+str(self.alignment_length)+"\n"
    ostr += "ANY_ERROR\t"+str(errtotal)+"\n"
    ostr += "MISMATCHES\t"+str(self.mismatches)+"\n"
    ostr += "ANY_DELETION\t"+str(self.deletions['total'])+"\n"
    ostr += "COMPLETE_DELETION\t"+str(self.deletions['specific'])+"\n"
    ostr += "HOMOPOLYMER_DELETION\t"+str(self.deletions['homopolymer'])+"\n"
    ostr += "ANY_INSERTION\t"+str(self.insertions['total'])+"\n"
    ostr += "COMPLETE_INSERTION\t"+str(self.insertions['specific'])+"\n"
    ostr += "HOMOPOLYMER_INSERTION\t"+str(self.insertions['homopolymer'])+"\n"
    return ostr
496432
def parse_arguments(argv):
  """"""
  parser = DatalabParser(
      epilog=('Note that if using a DNN model, --hidden-layer-size1=NUM, '
              '--hidden-layer-size2=NUM, ..., is also required. '),
      datalab_epilog=("""
  Note that if using a DNN model,
  hidden-layer-size1: NUM
  hidden-layer-size2: NUM
  ...
  is also required. """))

  # HP parameters
  parser.add_argument(
      '--epsilon', type=float, default=0.0005, metavar='R',
      help='tf.train.AdamOptimizer epsilon. Only used in dnn models.')
  parser.add_argument(
      '--l1-regularization', type=float, default=0.0, metavar='R',
      help='L1 term for linear models.')
  parser.add_argument(
      '--l2-regularization', type=float, default=0.0, metavar='R',
      help='L2 term for linear models.')

  # Model parameters
  parser.add_argument(
    '--model', required=True,
    choices=['linear_classification', 'linear_regression', 'dnn_classification', 'dnn_regression'])
  parser.add_argument(
      '--top-n', type=int, default=0, metavar='N',
      help=('For classification problems, the output graph will contain the '
            'labels and scores for the top n classes, and results will be in the form of '
            '"predicted, predicted_2, ..., probability, probability_2, ...". '
            'If --top-n=0, then all labels and scores are returned in the form of '
            '"predicted, class_name1, class_name2,...".'))

  # HP parameters
  parser.add_argument(
      '--learning-rate', type=float, default=0.01, metavar='R',
      help='optimizer learning rate.')

  # Training input parameters
  parser.add_argument(
      '--max-steps', type=int, metavar='N',
      help='Maximum number of training steps to perform. If unspecified, will '
           'honor "max-epochs".')
  parser.add_argument(
      '--max-epochs', type=int, default=1000, metavar='N',
      help='Maximum number of training data epochs on which to train. If '
           'both "max-steps" and "max-epochs" are specified, the training '
           'job will run for "max-steps" or "num-epochs", whichever occurs '
           'first. If early stopping is enabled, training may also stop '
           'earlier.')
  parser.add_argument(
      '--train-batch-size', type=int, default=64, metavar='N',
      help='How many training examples are used per step. If num-epochs is '
           'used, the last batch may not be full.')
  parser.add_argument(
      '--eval-batch-size', type=int, default=64, metavar='N',
      help='Batch size during evaluation. Larger values increase performance '
           'but also increase peak memory usgae on the master node. One pass '
           'over the full eval set is performed per evaluation run.')
  parser.add_argument(
      '--min-eval-frequency', type=int, default=1000, metavar='N',
      help='Minimum number of training steps between evaluations. Evaluation '
           'does not occur if no new checkpoint is available, hence, this is '
           'the minimum. If 0, the evaluation will only happen after training. ')
  parser.add_argument(
      '--early-stopping-num_evals', type=int, default=3,
      help='Automatic training stop after results of specified number of evals '
           'in a row show the model performance does not improve. Set to 0 to '
           'disable early stopping.')
  parser.add_argument(
      '--logging-level', choices=['error', 'warning', 'info'],
      help='The TF logging level. If absent, use info for cloud training '
           'and warning for local training.')

  args, remaining_args = parser.parse_known_args(args=argv[1:])

  # All HP parambeters must be unique, so we need to support an unknown number
  # of --hidden-layer-size1=10 --lhidden-layer-size2=10 ...
  # Look at remaining_args for hidden-layer-size\d+ to get the layer info.

  # Get number of layers
  pattern = re.compile('hidden-layer-size(\d+)')
  num_layers = 0
  for other_arg in remaining_args:
    match = re.search(pattern, other_arg)
    if match:
      if int(match.group(1)) <= 0:
        raise ValueError('layer size must be a positive integer. Was given %s' % other_arg)
      num_layers = max(num_layers, int(match.group(1)))

  # Build a new parser so we catch unknown args and missing layer_sizes.
  parser = argparse.ArgumentParser()
  for i in range(num_layers):
    parser.add_argument('--hidden-layer-size%s' % str(i + 1), type=int, required=True)

  layer_args = vars(parser.parse_args(args=remaining_args))
  hidden_layer_sizes = []
  for i in range(num_layers):
    key = 'hidden_layer_size%s' % str(i + 1)
    hidden_layer_sizes.append(layer_args[key])

  assert len(hidden_layer_sizes) == num_layers
  args.hidden_layer_sizes = hidden_layer_sizes

  return args
246803
def inside_brain(stat_dset,atlas=None,p=0.001):
    ''''''
    atlas = find_atlas(atlas)
    if atlas==None:
        return None
    mask_dset = nl.suffix(stat_dset,'_atlasfrac')
    nl.run(['3dfractionize','-template',nl.strip_subbrick(stat_dset),'-input',nl.calc([atlas],'1+step(a-100)',datum='short'),'-preserve','-clip','0.2','-prefix',mask_dset],products=mask_dset,quiet=True,stderr=None)
    s = nl.roi_stats(mask_dset,nl.thresh(stat_dset,p))
    return 100.0 * s[2]['nzvoxels'] / (s[1]['nzvoxels'] + s[2]['nzvoxels'])
1118650
def n_neurites(neurites, neurite_type=NeuriteType.all):
    ''''''
    return sum(1 for _ in iter_neurites(neurites, filt=is_type(neurite_type)))
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")