richardson lucy deconvolution of data with psf h
    using spatial convolutions (h should be small then)

    mode_conv = "fft" or "spatial"
760081
def deconv_rl(data, h, Niter=10, mode_conv="fft", log_iter=False):
    

    if isinstance(data, np.ndarray):
        mode_data = "np"
    elif isinstance(data, OCLArray):
        mode_data = "gpu"
    else:
        raise TypeError("array argument (1) has bad type: %s"%type(arr_obj))

    if mode_conv=="fft":
        res = _deconv_rl_np_fft(data, h, Niter) if mode_data=="np" else _deconv_rl_gpu_fft(data, h, Niter)
        return res
    elif mode_conv=="spatial":
        res = _deconv_rl_np_conv(data, h, Niter) if mode_data=="np" else gpu_conv(data, h, Niter)
        return res
    else:
        raise KeyError("mode_conv = %s not known, either 'fft' or 'spatial'"%mode_conv)
311443
def make_psf_kernel(psf, npix, cdelt, xpix, ypix, psf_scale_fn=None, normalize=False):
    

    egy = psf.energies
    x = make_pixel_distance(npix, xpix, ypix)
    x *= cdelt

    k = np.zeros((len(egy), npix, npix))
    for i in range(len(egy)):
        k[i] = psf.eval(i, x, scale_fn=psf_scale_fn)

    if normalize:
        k /= (np.sum(k, axis=0)[np.newaxis, ...] * np.radians(cdelt) ** 2)

    return k
1371562
function fft2d(re, im) {
    var tre = [],
        tim = [],
        i = 0;
    // x-axis
    for (var y = 0; y < _n; y++) {
        i = y * _n;
        for (var x1 = 0; x1 < _n; x1++) {
            tre[x1] = re[x1 + i];
            tim[x1] = im[x1 + i];
        }
        fft1d(tre, tim);
        for (var x2 = 0; x2 < _n; x2++) {
            re[x2 + i] = tre[x2];
            im[x2 + i] = tim[x2];
        }
    }
    // y-axis
    for (var x = 0; x < _n; x++) {
        for (var y1 = 0; y1 < _n; y1++) {
            i = x + y1 * _n;
            tre[y1] = re[i];
            tim[y1] = im[i];
        }
        fft1d(tre, tim);
        for (var y2 = 0; y2 < _n; y2++) {
            i = x + y2 * _n;
            re[i] = tre[y2];
            im[i] = tim[y2];
        }
    }
}
596719
def sheardec2D(X, shearletsystem):
    """"""
    coeffs = np.zeros(shearletsystem.shearlets.shape, dtype=complex)
    Xfreq = fftshift(fft2(ifftshift(X)))
    for i in range(shearletsystem.nShearlets):
        coeffs[:, :, i] = fftshift(ifft2(ifftshift(Xfreq * np.conj(
                                   shearletsystem.shearlets[:, :, i]))))
    return coeffs.real
640617
def corr2(X):
    ''' 
    '''
    X = np.atleast_3d(X)
    N = X.shape[0]
    D = X.shape[2]

    if D == 1:
        return np.zeros(N, dtype=np.float)

    trii = np.triu_indices(D, k=1)
    DD = len(trii[0])
    r = np.zeros((N, DD))
    for i in np.arange(N):
        rmat = np.corrcoef(X[i])  # get the ith window from each signal, result will be DxD
        r[i] = rmat[trii]
    return r
345730
def calc_signal_spatial(self):
        
        # Calculate the surface intensity
        self.surface_intensity_sparse = self.calc_surface_intensity()

        # Calculate the probability per object-by-object level
        self.surface_intensity_object = self.kernel.pdf(self.catalog.lon,
                                                        self.catalog.lat)
        
        # Spatial component of signal probability
        u_spatial = self.surface_intensity_object
        return u_spatial
50401
function parse_RgceArea(blob, length, opts) {
	var w = 2;
	if(opts) {
		if(opts.biff >= 2 && opts.biff <= 5) return parse_RgceArea_BIFF2(blob, length, opts);
		else if(opts.biff == 12) w = 4;
	}
	var r=blob.read_shift(w), R=blob.read_shift(w);
	var c=parse_ColRelU(blob, 2);
	var C=parse_ColRelU(blob, 2);
	return { s:{r:r, c:c[0], cRel:c[1], rRel:c[2]}, e:{r:R, c:C[0], cRel:C[1], rRel:C[2]} };
}
1273455
def _ncc_c_3dim(x, y):
    
    den = norm(x, axis=1)[:, None] * norm(y, axis=1)
    den[den == 0] = np.Inf
    x_len = x.shape[-1]
    fft_size = 1 << (2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size))[:, None])
    cc = np.concatenate((cc[:,:,-(x_len-1):], cc[:,:,:x_len]), axis=2)
    return np.real(cc) / den.T[:, :, None]
836242
def sepconv_relu_sepconv(inputs,
                         filter_size,
                         output_size,
                         first_kernel_size=(1, 1),
                         second_kernel_size=(1, 1),
                         padding="LEFT",
                         nonpadding_mask=None,
                         dropout=0.0,
                         name=None):
  """"""
  with tf.variable_scope(name, "sepconv_relu_sepconv", [inputs]):
    inputs = maybe_zero_out_padding(inputs, first_kernel_size, nonpadding_mask)
    if inputs.get_shape().ndims == 3:
      is_3d = True
      inputs = tf.expand_dims(inputs, 2)
    else:
      is_3d = False
    h = separable_conv(
        inputs,
        filter_size,
        first_kernel_size,
        activation=tf.nn.relu,
        padding=padding,
        name="conv1")
    if dropout != 0.0:
      h = tf.nn.dropout(h, 1.0 - dropout)
    h = maybe_zero_out_padding(h, second_kernel_size, nonpadding_mask)
    ret = separable_conv(
        h, output_size, second_kernel_size, padding=padding, name="conv2")
    if is_3d:
      ret = tf.squeeze(ret, 2)
    return ret
951950
def rfft2d_freqs(h, w):
    """"""

    fy = np.fft.fftfreq(h)[:, None]
    # when we have an odd input dimension we need to keep one additional
    # frequency and later cut off 1 pixel
    if w % 2 == 1:
        fx = np.fft.fftfreq(w)[: w // 2 + 2]
    else:
        fx = np.fft.fftfreq(w)[: w // 2 + 1]
    return np.sqrt(fx * fx + fy * fy)
1323895
def _l2ycbcr(self, mode):
        
        self._check_modes(("L", "LA"))

        luma = self.channels[0]
        zeros = np.ma.zeros(luma.shape)
        zeros.mask = luma.mask

        self.channels = [luma, zeros, zeros] + self.channels[1:]

        if self.fill_value is not None:
            self.fill_value = [self.fill_value[0], 0, 0] + self.fill_value[1:]

        self.mode = mode
356263
def convolve_fft(array, kernel):
    

    array = np.asarray(array, dtype=np.complex)
    kernel = np.asarray(kernel, dtype=np.complex)

    if array.ndim != kernel.ndim:
        raise ValueError("Image and kernel must have same number of "
                         "dimensions")

    array_shape = array.shape
    kernel_shape = kernel.shape
    new_shape = np.array(array_shape) + np.array(kernel_shape)

    array_slices = []
    kernel_slices = []
    for (new_dimsize, array_dimsize, kernel_dimsize) in zip(
            new_shape, array_shape, kernel_shape):
        center = new_dimsize - (new_dimsize + 1) // 2
        array_slices += [slice(center - array_dimsize // 2,
                         center + (array_dimsize + 1) // 2)]
        kernel_slices += [slice(center - kernel_dimsize // 2,
                          center + (kernel_dimsize + 1) // 2)]

    array_slices = tuple(array_slices)
    kernel_slices = tuple(kernel_slices)

    if not np.all(new_shape == array_shape):
        big_array = np.zeros(new_shape, dtype=np.complex)
        big_array[array_slices] = array
    else:
        big_array = array

    if not np.all(new_shape == kernel_shape):
        big_kernel = np.zeros(new_shape, dtype=np.complex)
        big_kernel[kernel_slices] = kernel
    else:
        big_kernel = kernel

    array_fft = np.fft.fftn(big_array)
    kernel_fft = np.fft.fftn(np.fft.ifftshift(big_kernel))

    rifft = np.fft.ifftn(array_fft * kernel_fft)

    return rifft[array_slices].real
1034572
def lowpass(ts, cutoff_hz, order=3):
    """"""
    orig_ndim = ts.ndim
    if ts.ndim is 1:
        ts = ts[:, np.newaxis]
    channels = ts.shape[1]
    fs = (len(ts) - 1.0) / (ts.tspan[-1] - ts.tspan[0])
    nyq = 0.5 * fs
    cutoff = cutoff_hz/nyq
    b, a = signal.butter(order, cutoff, btype='low')
    if not np.all(np.abs(np.roots(a)) < 1.0):
        raise ValueError('Filter will not be stable with these values.')
    dtype = ts.dtype
    output = np.zeros((len(ts), channels), dtype)
    for i in range(channels):
        output[:, i] = signal.filtfilt(b, a, ts[:, i])
    if orig_ndim is 1:
        output = output[:, 0]
    return Timeseries(output, ts.tspan, labels=ts.labels)
460274
def _cutout_psf(self, image, subgrid_res):
        
        self._init_mask_psf()
        return image[self._x_min_psf*subgrid_res:(self._x_max_psf+1)*subgrid_res, self._y_min_psf*subgrid_res:(self._y_max_psf+1)*subgrid_res]
201636
def _fft_convolve_numpy(data, h, plan = None,
                        kernel_is_fft = False,
                        kernel_is_fftshifted = False):
    

    if data.shape != h.shape:
        raise ValueError("data and kernel must have same size! %s vs %s "%(str(data.shape),str(h.shape)))

    
    data_g = OCLArray.from_array(data.astype(np.complex64))

    if not kernel_is_fftshifted:
        h = np.fft.fftshift(h)

    
    h_g = OCLArray.from_array(h.astype(np.complex64))
    res_g = OCLArray.empty_like(data_g)
    
    _fft_convolve_gpu(data_g,h_g,res_g = res_g,
                      plan = plan,
                      kernel_is_fft = kernel_is_fft)

    res =  abs(res_g.get())

    del data_g
    del h_g
    del res_g
    
    return res
1225092
def convolve_sep2(data, hx, hy, res_g=None, sub_blocks=None):
    

    if isinstance(data, np.ndarray):
        data = np.ascontiguousarray(data)

        if sub_blocks == (1, 1) or sub_blocks is None:
            return _convolve_sep2_numpy(data, hx, hy)
        else:
            # cut the image into tile and operate on every of them
            N_sub = [int(np.ceil(1. * n / s)) for n, s in zip(data.shape, sub_blocks)]
            Npads = [int(len(_h) / 2) for _h in [hy, hx]]
            res = np.empty(data.shape, np.float32)
            for data_tile, data_s_src, data_s_dest \
                    in tile_iterator(data, blocksize=N_sub,
                                     padsize=Npads,
                                     mode="constant"):
                res_tile = _convolve_sep2_numpy(data_tile.copy(),
                                                hx, hy)
                res[data_s_src] = res_tile[data_s_dest]
            return res
    elif isinstance(data, OCLArray):
        if not sub_blocks is None:
            raise NotImplementedError()
        return _convolve_sep2_gpu(data, hx, hy, res_g=res_g)
    else:
        raise TypeError("array argument (1) has bad type: %s" % type(data))
1000190
def local_attention1d_spatial_decoder(x, kv_dim, heads_dim,
                                      feedforward_dim, hparams):
  """"""
  batch_dim, length_dim, model_dim = x.shape.dims
  blocks_w_dim = mtf.Dimension("blocksw", hparams.block_length)
  num_w_blocks_dim = mtf.Dimension("num_wblocks",
                                   length_dim.size // blocks_w_dim.size)
  x = mtf.reshape(
      x, mtf.Shape([batch_dim, num_w_blocks_dim, blocks_w_dim, model_dim]))
  # [ 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_self_attention_spatial_blocks(
              mtf.layers.layer_norm(x, model_dim, name="layer_norm_att"),
              kv_dim,
              heads_dim,
              memory_w_dim=blocks_w_dim,
              mask_right=True,
              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