Cache sigma calculate for use in tandem with the FilterBank class
171199
def sigma_cached(self, psd):
    
    if not hasattr(self, '_sigmasq'):
        from pycbc.opt import LimitedSizeDict
        self._sigmasq = LimitedSizeDict(size_limit=2**5)

    key = id(psd)
    if not hasattr(psd, '_sigma_cached_key'):
        psd._sigma_cached_key = {}

    if key not in self._sigmasq or id(self) not in psd._sigma_cached_key:
        psd._sigma_cached_key[id(self)] = True
        # If possible, we precalculate the sigmasq vector for all possible waveforms
        if pycbc.waveform.waveform_norm_exists(self.approximant):
            if not hasattr(psd, 'sigmasq_vec'):
                psd.sigmasq_vec = {}

            if self.approximant not in psd.sigmasq_vec:
                psd.sigmasq_vec[self.approximant] = pycbc.waveform.get_waveform_filter_norm(
                     self.approximant, psd, len(psd), psd.delta_f, self.f_lower)

            if not hasattr(self, 'sigma_scale'):
                # Get an amplitude normalization (mass dependant constant norm)
                amp_norm = pycbc.waveform.get_template_amplitude_norm(
                                     self.params, approximant=self.approximant)
                amp_norm = 1 if amp_norm is None else amp_norm
                self.sigma_scale = (DYN_RANGE_FAC * amp_norm) ** 2.0


            self._sigmasq[key] = self.sigma_scale * \
                psd.sigmasq_vec[self.approximant][self.end_idx-1]

        else:
            if not hasattr(self, 'sigma_view'):
                from pycbc.filter.matchedfilter import get_cutoff_indices
                N = (len(self) -1) * 2
                kmin, kmax = get_cutoff_indices(
                        self.min_f_lower or self.f_lower, self.end_frequency,
                        self.delta_f, N)
                self.sslice = slice(kmin, kmax)
                self.sigma_view = self[self.sslice].squared_norm() * 4.0 * self.delta_f

            if not hasattr(psd, 'invsqrt'):
                psd.invsqrt = 1.0 / psd[self.sslice]

            self._sigmasq[key] = self.sigma_view.inner(psd.invsqrt)
    return self._sigmasq[key]
992565
private long[] calcBucketOffset(long timestampMs) {
        long bucketOffset = toTimeSeriesPoint(timestampMs);
        long delta = bucketOffset % (RESOLUTION_MS * BUCKET_SIZE);
        long bucketId = bucketOffset - delta;
        return new long[] { bucketId, bucketOffset };
    }
427806
def smear(self, sigma):
        
        diff = [self.x[i + 1] - self.x[i] for i in range(len(self.x) - 1)]
        avg_x_per_step = np.sum(diff) / len(diff)
        if len(self.ydim) == 1:
            self.y = gaussian_filter1d(self.y, sigma / avg_x_per_step)
        else:
            self.y = np.array([
                gaussian_filter1d(self.y[:, k], sigma / avg_x_per_step)
                for k in range(self.ydim[1])]).T
864005
def _sigmainf(N, h, m, dW, Km0, Pm0):
    """"""
    M = m*(m-1)//2
    Im = broadcast_to(np.eye(m), (N, m, m))
    IM = broadcast_to(np.eye(M), (N, M, M))
    Ims0 = np.eye(m**2)
    factor1 = broadcast_to((2.0/h)*np.dot(Km0, Ims0 - Pm0), (N, M, m**2))
    factor2 = _kp2(Im, _dot(dW, _t(dW)))
    factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
    return 2*IM + _dot(_dot(factor1, factor2), factor3)
250904
def local_attention2d_spatial_decoder(x, kv_dim, heads_dim,
                                      feedforward_dim, hparams):
  """"""
  batch_dim, length_dim, model_dim = x.shape.dims
  blocks_h_dim = mtf.Dimension("blocksh", hparams.block_height)
  blocks_w_dim = mtf.Dimension("blocksw", hparams.block_width)
  num_h_blocks_dim = mtf.Dimension("num_h_blocks",
                                   hparams.img_len // hparams.block_height)
  num_w_blocks_dim = mtf.Dimension(
      "num_w_blocks",
      hparams.img_len * hparams.num_channels // hparams.block_width)
  x = mtf.transpose(
      mtf.reshape(
          x,
          mtf.Shape([
              batch_dim, num_h_blocks_dim, blocks_h_dim,
              num_w_blocks_dim, blocks_w_dim, model_dim
          ])),
      mtf.Shape([
          batch_dim, num_h_blocks_dim, num_w_blocks_dim,
          blocks_h_dim, blocks_w_dim, model_dim
      ]))
  # Image Transformer Decoder
  # [ self attention - ffn - residual + dropout] x n
  for layer in range(hparams.num_decoder_layers):
    layer_name = "decoder_layer_%d" % layer
    with tf.variable_scope(layer_name):
      # Self attention layer
      x += layer_prepostprocess_dropout(
          mtf.layers.local_2d_self_attention_spatial_blocks(
              mtf.layers.layer_norm(x, model_dim, name="layer_norm_att"),
              kv_dim,
              heads_dim,
              memory_h_dim=num_h_blocks_dim,
              memory_w_dim=num_w_blocks_dim,
              name="self_att"), hparams)
      # ffn layer
      x += layer_prepostprocess_dropout(
          mtf.layers.dense_relu_dense(
              mtf.layers.layer_norm(x, model_dim, name="layer_norm_ffn"),
              feedforward_dim,
              hparams.dropout,
              dropout_broadcast_dims=[length_dim]), hparams)

  output = mtf.layers.layer_norm(x, model_dim, name="final_layer_norm")
  return output
401210
def ks_synth(freq):
  
  ks_mem = (sum(lz.sinusoid(x * freq) for x in [1, 3, 9]) +
            lz.white_noise() + lz.Stream(-1, 1)) / 5
  return lz.karplus_strong(freq, memory=ks_mem)
1358570
def algorithm_u(ns, m):
    
    def visit(nv, av):
        ps = [[] for _ in range(m)]
        for jv in range(nv):
            ps[av[jv + 1]].append(ns[jv])
        return ps

    def f(mu, nu, sigma, n, a):
        if mu == 2:
            yield visit(n, a)
        else:
            for v in f(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
                yield v
        if nu == mu + 1:
            a[mu] = mu - 1
            yield visit(n, a)
            while a[nu] > 0:
                a[nu] = a[nu] - 1
                yield visit(n, a)
        elif nu > mu + 1:
            if (mu + sigma) % 2 == 1:
                a[nu - 1] = mu - 1
            else:
                a[mu] = mu - 1
            if (a[nu] + sigma) % 2 == 1:
                for v in b(mu, nu - 1, 0, n, a):
                    yield v
            else:
                for v in f(mu, nu - 1, 0, n, a):
                    yield v
            while a[nu] > 0:
                a[nu] = a[nu] - 1
                if (a[nu] + sigma) % 2 == 1:
                    for v in b(mu, nu - 1, 0, n, a):
                        yield v
                else:
                    for v in f(mu, nu - 1, 0, n, a):
                        yield v

    def b(mu, nu, sigma, n, a):
        if nu == mu + 1:
            while a[nu] < mu - 1:
                yield visit(n, a)
                a[nu] = a[nu] + 1
            yield visit(n, a)
            a[mu] = 0
        elif nu > mu + 1:
            if (a[nu] + sigma) % 2 == 1:
                for v in f(mu, nu - 1, 0, n, a):
                    yield v
            else:
                for v in b(mu, nu - 1, 0, n, a):
                    yield v
            while a[nu] < mu - 1:
                a[nu] = a[nu] + 1
                if (a[nu] + sigma) % 2 == 1:
                    for v in f(mu, nu - 1, 0, n, a):
                        yield v
                else:
                    for v in b(mu, nu - 1, 0, n, a):
                        yield v
            if (mu + sigma) % 2 == 1:
                a[nu - 1] = 0
            else:
                a[mu] = 0
        if mu == 2:
            yield visit(n, a)
        else:
            for v in b(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
                yield v

    ng = len(ns)
    ag = [0] * (ng + 1)
    for j in range(1, m + 1):
        ag[ng - m + j] = j - 1
    return f(m, ng, 0, ng, ag) if m > 1 else [[ns]]
1099446
def get_sigma(imt):
    
    if imt.period < 0.2:
        return np.log(10**0.23)
    elif imt.period > 1.0:
        return np.log(10**0.27)
    else:
        return np.log(10**(0.23 + (imt.period - 0.2)/0.8 * 0.04))
1153751
def errors(self):
        

        n = self.n
        dt = self.dt
        order = self._order
        sigma = self.sigma

        error = np.zeros(order + 1)
        std = np.zeros(order + 1)

        if n == 0:
            return (error, std)

        if order == 0:
            error[0] = sigma/sqrt(n)
            std[0] = sigma/sqrt(n)

        elif order == 1:
            if n > 1:
                error[0] = sigma * sqrt(2*(2*n-1) / (n*(n+1)))
                error[1] = sigma * sqrt(12. / (n*(n*n-1)*dt*dt))
            std[0] = sigma * sqrt((2*(2*n-1)) / (n*(n+1)))
            std[1] = (sigma/dt) * sqrt(12. / (n*(n*n-1)))

        elif order == 2:
            dt2 = dt * dt

            if n >= 3:
                error[0] = sigma * sqrt(3*(3*n*n-3*n+2) / (n*(n+1)*(n+2)))
                error[1] = sigma * sqrt(12*(16*n*n-30*n+11) /
                                        (n*(n*n-1)*(n*n-4)*dt2))
                error[2] = sigma * sqrt(720/(n*(n*n-1)*(n*n-4)*dt2*dt2))

            std[0] = sigma * sqrt((3*(3*n*n - 3*n + 2)) / (n*(n+1)*(n+2)))
            std[1] = (sigma/dt) * sqrt((12*(16*n*n - 30*n + 11)) /
                                       (n*(n*n - 1)*(n*n - 4)))
            std[2] = (sigma/dt2) * sqrt(720 / (n*(n*n-1)*(n*n-4)))

        return error, std
1238127
function convolute2DI(ftSignal, ftFilter, ftRows, ftCols) {
    var re, im;
    for (var iRow = 0; iRow < ftRows / 2; iRow++) {
        for (var iCol = 0; iCol < ftCols; iCol++) {
            //
            re = ftSignal[(iRow * 2) * ftCols + iCol]
                * ftFilter[(iRow * 2) * ftCols + iCol]
                - ftSignal[(iRow * 2 + 1) * ftCols + iCol]
                * ftFilter[(iRow * 2 + 1) * ftCols + iCol];
            im = ftSignal[(iRow * 2) * ftCols + iCol]
                * ftFilter[(iRow * 2 + 1) * ftCols + iCol]
                + ftSignal[(iRow * 2 + 1) * ftCols + iCol]
                * ftFilter[(iRow * 2) * ftCols + iCol];
            //
            ftSignal[(iRow * 2) * ftCols + iCol] = re;
            ftSignal[(iRow * 2 + 1) * ftCols + iCol] = im;
        }
    }
}
1186593
def filter_data(data, filter_name, filter_file, bits, filterbank_off=False,
                    swstat_channel_name=None):
    '''
    
    '''

    # if filterbank is off then return a time series of zeroes
    if filterbank_off:
        return numpy.zeros(len(data))

    # loop over the 10 filters in the filterbank
    for i in range(10):

        # read the filter
        filter = Filter(filter_file[filter_name][i])

        # if bit is on then filter the data
        bit = int(bits[-(i+1)])
        if bit:
            logging.info('filtering with filter module %d', i)

            # if there are second-order sections then filter with them
            if len(filter.sections):
                data = filter.apply(data)

            # else it is a filter with only gain so apply the gain
            else:
                coeffs = iir2z(filter_file[filter_name][i])
                if len(coeffs) > 1:
                    logging.info('Gain-only filter module return more than one number')
                    sys.exit()
                gain = coeffs[0]
                data = gain * data

    return  data
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")
480571
public static <T extends ImageBase<T>>
	boolean invokeNativeGaussian(T input, T output, double sigma , int radius, T storage) {
		boolean processed = false;
		if( BOverrideBlurImageOps.gaussian != null ) {
			try {
				BOverrideBlurImageOps.gaussian.processGaussian(input,output,sigma,radius,storage);
				processed = true;
			} catch( RuntimeException ignore ) {}
		}
		return processed;
	}
938054
void estimateTranslation( DMatrixRMaj R , List<AssociatedPair> points , Vector3D_F64 T )
	{
		final int N = points.size();
		W.reshape(N*2,3);
		y.reshape(N*2,1);
		Wty.reshape(3,1);

		DMatrix3x3 Rtmp = new DMatrix3x3();
		ConvertDMatrixStruct.convert(R,Rtmp);

		int indexY = 0,indexW = 0;
		for (int i = 0; i < N; i++) {
			AssociatedPair p = points.get(i);

			// rotate into camera frame
			double u1 = Rtmp.a11*p.p1.x + Rtmp.a12*p.p1.y;
			double u2 = Rtmp.a21*p.p1.x + Rtmp.a22*p.p1.y;
			double u3 = Rtmp.a31*p.p1.x + Rtmp.a32*p.p1.y;

			W.data[indexW++] = 1;
			W.data[indexW++] = 0;
			W.data[indexW++] = -p.p2.x;
			W.data[indexW++] = 0;
			W.data[indexW++] = 1;
			W.data[indexW++] = -p.p2.y;

			y.data[indexY++] = p.p2.x*u3 - u1;
			y.data[indexY++] = p.p2.y*u3 - u2;
		}

		//======= Compute Pseudo Inverse
		// WW = inv(W^T*W)
		CommonOps_DDRM.multTransA(W,W,WW);
		CommonOps_DDRM.invert(WW);

		// W^T*y
		CommonOps_DDRM.multTransA(W,y,Wty);

		// translation = inv(W^T*W)*W^T*y
		W.reshape(3,1);
		CommonOps_DDRM.mult(WW,Wty,W);

		T.x = W.data[0];
		T.y = W.data[1];
		T.z = W.data[2];
	}
267366
def precompute(self, cache_dir=None, swath_usage=0, **kwargs):
        """"""
        if kwargs.get('mask') is not None:
            LOG.warning("'mask' parameter has no affect during EWA "
                        "resampling")

        del kwargs
        source_geo_def = self.source_geo_def
        target_geo_def = self.target_geo_def

        if cache_dir:
            LOG.warning("'cache_dir' is not used by EWA resampling")

        # Satpy/PyResample don't support dynamic grids out of the box yet
        lons, lats = source_geo_def.get_lonlats()
        if isinstance(lons, xr.DataArray):
            # get dask arrays
            lons = lons.data
            lats = lats.data
        # we are remapping to a static unchanging grid/area with all of
        # its parameters specified
        chunks = (2,) + lons.chunks
        res = da.map_blocks(self._call_ll2cr, lons, lats,
                            target_geo_def, swath_usage,
                            dtype=lons.dtype, chunks=chunks, new_axis=[0])
        cols = res[0]
        rows = res[1]

        # save the dask arrays in the class instance cache
        # the on-disk cache will store the numpy arrays
        self.cache = {
            "rows": rows,
            "cols": cols,
        }

        return None
1091418
function decompressSamples (o, receivedDeltas) {
  // add the delta to the previous value
  for (let i = 1; i < 3; i++) {
    for (let j = 0; j < 4; j++) {
      o.decompressedSamples[i][j] = o.decompressedSamples[i - 1][j] - receivedDeltas[i - 1][j];
    }
  }
}
1205698
private double[] findAndEvaluateThreshold(DoubleDynamicHistogram histogram) {
    int n = histogram.getNumBins();
    double[] p1 = new double[n];
    double[] p2 = new double[n];
    double[] mu1 = new double[n];
    double[] mu2 = new double[n];
    double[] sigma1 = new double[n];
    double[] sigma2 = new double[n];
    double[] jt = new double[n];
    // Forward pass
    {
      MeanVariance mv = new MeanVariance();
      DoubleHistogram.Iter forward = histogram.iter();
      for(int i = 0; forward.valid(); i++, forward.advance()) {
        p1[i] = forward.getValue() + ((i > 0) ? p1[i - 1] : 0);
        mv.put(i, forward.getValue());
        mu1[i] = mv.getMean();
        sigma1[i] = mv.getNaiveStddev();
      }
    }
    // Backwards pass
    {
      MeanVariance mv = new MeanVariance();
      DoubleHistogram.Iter backwards = histogram.iter();
      backwards.seek(histogram.getNumBins() - 1); // Seek to last

      for(int j = n - 1; backwards.valid(); j--, backwards.retract()) {
        p2[j] = backwards.getValue() + ((j + 1 < n) ? p2[j + 1] : 0);
        mv.put(j, backwards.getValue());
        mu2[j] = mv.getMean();
        sigma2[j] = mv.getNaiveStddev();
      }
    }

    for(int i = 0; i < n; i++) {
      jt[i] = 1.0 + 2 * (p1[i] * (FastMath.log(sigma1[i]) - FastMath.log(p1[i])) + p2[i] * (FastMath.log(sigma2[i]) - FastMath.log(p2[i])));
    }

    int bestpos = -1;
    double bestgoodness = Double.NEGATIVE_INFINITY;

    double devPrev = jt[1] - jt[0];
    for(int i = 1; i < jt.length - 1; i++) {
      double devCur = jt[i + 1] - jt[i];
      // Local minimum found - calculate depth
      if(devCur >= 0 && devPrev <= 0) {
        double lowestMaxima = Double.POSITIVE_INFINITY;
        for(int j = i - 1; j > 0; j--) {
          if(jt[j - 1] < jt[j]) {
            lowestMaxima = Math.min(lowestMaxima, jt[j]);
            break;
          }
        }
        for(int j = i + 1; j < n - 2; j++) {
          if(jt[j + 1] < jt[j]) {
            lowestMaxima = Math.min(lowestMaxima, jt[j]);
            break;
          }
        }
        double localDepth = lowestMaxima - jt[i];

        final double mud = mu1[i] - mu2[i];
        double discriminability = mud * mud / (sigma1[i] * sigma1[i] + sigma2[i] * sigma2[i]);
        if(Double.isNaN(discriminability)) {
          discriminability = -1;
        }
        double goodness = localDepth * discriminability;
        if(goodness > bestgoodness) {
          bestgoodness = goodness;
          bestpos = i;
        }
      }
      devPrev = devCur;
    }
    Iter iter = histogram.iter();
    iter.seek(bestpos);
    return new double[] { iter.getRight(), bestgoodness };
  }