continuous haar wavelet transform based on the paper
    "A practical guide to wavelet analysis" by Christopher Torrence and Gilbert P Compo
166895
def cwt_haar(x: np.ndarray, scale=10):
    
    next_power_two = 2 ** int(np.log2(len(x)))

    x = x[0:next_power_two]
    num_data = len(x)

    # get FFT of x (eq. (3) in paper)
    x_hat = np.fft.fft(x)

    # Get omega (eq. (5) in paper)
    f = (2.0 * np.pi / num_data)
    omega = f * np.concatenate((np.arange(0, num_data // 2), np.arange(num_data // 2, num_data) * -1))

    # get psi hat (eq. (6) in paper)
    psi_hat = np.sqrt(2.0 * np.pi * scale) * normalized_haar_wavelet(scale * omega, scale)

    # get W (eq. (4) in paper)
    W = np.fft.ifft(x_hat * psi_hat)

    return W[2 * scale:-2 * scale]
983745
function hslToRgb(h, s, l, a, out) {
		function hueToRgb(m1, m2, h) {
			h = h % 1;
			if (h < 0) {
				h += 1;
			}
			if (h < 1 / 6) {
				return m1 + (m2 - m1) * h * 6;
			}
			if (h < 1 / 2) {
				return m2;
			}
			if (h < 2 / 3) {
				return m1 + (m2 - m1) * (2/3 - h) * 6;
			}
			return m1;
		}

		var m1, m2;
		if (l < 0.5) {
			m2 = l * (s + 1);
		} else {
			m2 = l + s - l * s;
		}
		m1 = l * 2 - m2;

		if (!out) {
			out = [];
		}

		out[0] = hueToRgb(m1, m2, h + 1/3);
		out[1] = hueToRgb(m1, m2, h);
		out[2] = hueToRgb(m1, m2, h - 1/3);
		out[3] = a;

		return out;
	}
366004
def pywt_wavelet(wavelet):
    """"""
    if isinstance(wavelet, pywt.Wavelet):
        return wavelet
    else:
        return pywt.Wavelet(wavelet)
125572
function ifft2DArray(ft, ftRows, ftCols) {
    var tempTransform = new Array(ftRows * ftCols);
    var nRows = ftRows / 2;
    var nCols = (ftCols - 1) * 2;
    // reverse transform columns
    FFT.init(nRows);
    var tmpCols = {re: new Array(nRows), im: new Array(nRows)};
    var iRow, iCol;
    for (iCol = 0; iCol < ftCols; iCol++) {
        for (iRow = nRows - 1; iRow >= 0; iRow--) {
            tmpCols.re[iRow] = ft[(iRow * 2) * ftCols + iCol];
            tmpCols.im[iRow] = ft[(iRow * 2 + 1) * ftCols + iCol];
        }
        //Unnormalized inverse transform
        FFT.bt(tmpCols.re, tmpCols.im);
        for (iRow = nRows - 1; iRow >= 0; iRow--) {
            tempTransform[(iRow * 2) * ftCols + iCol] = tmpCols.re[iRow];
            tempTransform[(iRow * 2 + 1) * ftCols + iCol] = tmpCols.im[iRow];
        }
    }

    // reverse row transform
    var finalTransform = new Array(nRows * nCols);
    FFT.init(nCols);
    var tmpRows = {re: new Array(nCols), im: new Array(nCols)};
    var scale = nCols * nRows;
    for (iRow = 0; iRow < ftRows; iRow += 2) {
        tmpRows.re[0] = tempTransform[iRow * ftCols];
        tmpRows.im[0] = tempTransform[(iRow + 1) * ftCols];
        for (iCol = 1; iCol < ftCols; iCol++) {
            tmpRows.re[iCol] = tempTransform[iRow * ftCols + iCol];
            tmpRows.im[iCol] = tempTransform[(iRow + 1) * ftCols + iCol];
            tmpRows.re[nCols - iCol] = tempTransform[iRow * ftCols + iCol];
            tmpRows.im[nCols - iCol] = -tempTransform[(iRow + 1) * ftCols + iCol];
        }
        //Unnormalized inverse transform
        FFT.bt(tmpRows.re, tmpRows.im);

        var indexB = (iRow / 2) * nCols;
        for (iCol = nCols - 1; iCol >= 0; iCol--) {
            finalTransform[indexB + iCol] = tmpRows.re[iCol] / scale;
        }
    }
    return finalTransform;
}
1101605
def imagetransformer_base_10l_8h_big_uncond_dr03_dan_64_2d():
  """"""
  hparams = image_transformer2d_base()
  hparams.unconditional = True
  hparams.hidden_size = 512
  hparams.batch_size = 1
  hparams.img_len = 64
  hparams.num_heads = 8
  hparams.filter_size = 2048
  hparams.batch_size = 1
  hparams.max_length = 3075
  hparams.max_length = 14000
  hparams.layer_preprocess_sequence = "none"
  hparams.layer_postprocess_sequence = "dan"
  hparams.layer_prepostprocess_dropout = 0.1
  hparams.dec_attention_type = cia.AttentionType.LOCAL_2D
  hparams.query_shape = (16, 16)
  hparams.memory_flange = (8, 8)
  return hparams
684536
def Morlet (sampr, freq, width=7.0):
  
  dt = 1. / sampr
  sf = freq / width # sf in Hz
  st = 1. / (2. * np.pi * sf) # st in s
  t = np.arange(-3.5*st, 3.5*st, dt)
  A = 1. / (st * np.sqrt(2.*np.pi)) # A in 1 / s
  # units: 1/s * (exp (s**2 / s**2)) * exp( 1/ s * s)
  return A * np.exp(-t**2. / (2. * st**2.)) * np.exp(1.j * 2. * np.pi * freq * t)
388592
public static WaveletDescription<WlCoef_I32> biorthogonal_I32( int J ,
																   BorderType borderType ) {
		if( J != 5 ) {
			throw new IllegalArgumentException("Only 5 is currently supported");
		}

		WlCoef_I32 forward = new WlCoef_I32();

		forward.offsetScaling = -2;
		forward.offsetWavelet = 0;

		forward.scaling = new int[5];
		forward.wavelet = new int[3];

		forward.denominatorScaling = 8;
		forward.scaling[0] = -1;
		forward.scaling[1] = 2;
		forward.scaling[2] = 6;
		forward.scaling[3] = 2;
		forward.scaling[4] = -1;

		forward.denominatorWavelet = 2;
		forward.wavelet[0] = -1;
		forward.wavelet[1] = 2;
		forward.wavelet[2] = -1;

		BorderIndex1D border;
		WlBorderCoef<WlCoef_I32> inverse;

		if( borderType == BorderType.WRAP ) {
			WlCoef_I32 inner = computeInnerBiorthogonalInverse(forward);
			inverse = new WlBorderCoefStandard<>(inner);
			border = new BorderIndex1D_Wrap();
		} else if( borderType == BorderType.REFLECT ) {
			WlCoef_I32 inner = computeInnerBiorthogonalInverse(forward);
			inverse = convertToInt((WlBorderCoefFixed<WlCoef_F32>)biorthogonal_F32(J,borderType).getInverse(),inner);
			border = new BorderIndex1D_Reflect();
		} else {
			throw new IllegalArgumentException("Unsupported border type: "+borderType);
		}
		return new WaveletDescription<>(border, forward, inverse);

	}
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];
        }
    }
}
946167
def laplacian_of_gaussian(image, mask, size, sigma):
    '''
    '''
    half_size = size//2
    i,j = np.mgrid[-half_size:half_size+1,
                   -half_size:half_size+1].astype(float) / float(sigma)
    distance = (i**2 + j**2)/2
    gaussian = np.exp(-distance)
    #
    # Normalize the Gaussian
    #
    gaussian = gaussian / np.sum(gaussian)

    log = (distance - 1) * gaussian
    #
    # Normalize the kernel to have a sum of zero
    #
    log = log - np.mean(log)
    if mask is None:
        mask = np.ones(image.shape[:2], bool)
    masked_image = image.copy()
    masked_image[~mask] = 0
    output = convolve(masked_image, log, mode='constant', cval=0)
    #
    # Do the LoG of the inverse of the mask. This finds the magnitude of the
    # contribution of the masked pixels. We then fudge by multiplying by the
    # value at the pixel of interest - this effectively sets the value at a
    # masked pixel to that of the pixel of interest.
    #
    # It underestimates the LoG, that's not a terrible thing.
    #
    correction = convolve((~ mask).astype(float), log, mode='constant', cval = 1)
    output += correction * image
    output[~ mask] = image[~ mask]
    return output
98485
function dither(cb) {
  const rgb565Matrix = [1, 9, 3, 11, 13, 5, 15, 7, 4, 12, 2, 10, 16, 8, 14, 6];
  this.scanQuiet(0, 0, this.bitmap.width, this.bitmap.height, function(
    x,
    y,
    idx
  ) {
    const thresholdId = ((y & 3) << 2) + (x % 4);
    const dither = rgb565Matrix[thresholdId];
    this.bitmap.data[idx] = Math.min(this.bitmap.data[idx] + dither, 0xff);
    this.bitmap.data[idx + 1] = Math.min(
      this.bitmap.data[idx + 1] + dither,
      0xff
    );
    this.bitmap.data[idx + 2] = Math.min(
      this.bitmap.data[idx + 2] + dither,
      0xff
    );
  });

  if (isNodePattern(cb)) {
    cb.call(this, null, this);
  }

  return this;
}
338540
function DFT(bufferSize, sampleRate) {
  FourierTransform.call(this, bufferSize, sampleRate);

  var N = bufferSize/2 * bufferSize;
  var TWO_PI = 2 * Math.PI;

  this.sinTable = new Float32Array(N);
  this.cosTable = new Float32Array(N);

  for (var i = 0; i < N; i++) {
    this.sinTable[i] = Math.sin(i * TWO_PI / bufferSize);
    this.cosTable[i] = Math.cos(i * TWO_PI / bufferSize);
  }
}
15609
public static void horizontalInverse(BorderIndex1D border , WlBorderCoef<WlCoef_I32> inverseCoef , GrayI input , GrayI output ) {

		UtilWavelet.checkShape(output,input);

		int []trends = new int[ output.width ];
		int []details = new int[ output.width ];

		boolean isLarger = input.width >= output.width;
		int paddedWidth = output.width + output.width%2;

		final int lowerBorder = inverseCoef.getLowerLength()*2;
		final int upperBorder = output.width - inverseCoef.getUpperLength()*2;

		border.setLength(output.width+output.width%2);

		WlCoef_I32 coefficients = inverseCoef.getInnerCoefficients();
		final int e = coefficients.denominatorScaling*2;
		final int f = coefficients.denominatorWavelet*2;
		final int ef = e*f;
		final int ef2 = ef/2;

		for( int y = 0; y < output.height; y++ ) {

			for( int i = 0; i < details.length; i++ ) {
				details[i] = 0;
				trends[i] = 0;
			}

			for( int x = 0; x < output.width; x += 2 ) {
				int a = input.get(x/2,y);
				int d = input.get(input.width/2+x/2,y);

				if( x < lowerBorder ) {
					coefficients = inverseCoef.getBorderCoefficients(x);
				} else if( x >= upperBorder ) {
					coefficients = inverseCoef.getBorderCoefficients(x-paddedWidth);
				} else {
					coefficients = inverseCoef.getInnerCoefficients();
				}

				final int offsetA = coefficients.offsetScaling;
				final int offsetB = coefficients.offsetWavelet;
				final int[] alpha = coefficients.scaling;
				final int[] beta = coefficients.wavelet;

				// add the trend
				for( int i = 0; i < alpha.length; i++ ) {
					// if an odd image don't update the outer edge
					int xx = border.getIndex(x+offsetA+i);
					if( isLarger && xx >= output.width )
						continue;
					trends[xx] += a*alpha[i];
				}

				// add the detail signal
				for( int i = 0; i < beta.length; i++ ) {
					int xx = border.getIndex(x+offsetB+i);
					if( isLarger && xx >= output.width )
						continue;
					details[xx] += d*beta[i];
				}
			}

			for( int x = 0; x < output.width; x++ ) {
				output.set(x,y, UtilWavelet.round(trends[x]*f + details[x]*e ,ef2,ef));
			}
		}
	}
1008775
def funcSmthSpt(aryFuncChnk, varSdSmthSpt):
    
    varNdim = aryFuncChnk.ndim

    # Number of time points in this chunk:
    varNumVol = aryFuncChnk.shape[-1]

    # Loop through volumes:
    if varNdim == 4:
        for idxVol in range(0, varNumVol):

            aryFuncChnk[:, :, :, idxVol] = gaussian_filter(
                aryFuncChnk[:, :, :, idxVol],
                varSdSmthSpt,
                order=0,
                mode='nearest',
                truncate=4.0)
    elif varNdim == 5:
        varNumMtnDrctns = aryFuncChnk.shape[3]
        for idxVol in range(0, varNumVol):
            for idxMtn in range(0, varNumMtnDrctns):
                aryFuncChnk[:, :, :, idxMtn, idxVol] = gaussian_filter(
                    aryFuncChnk[:, :, :, idxMtn, idxVol],
                    varSdSmthSpt,
                    order=0,
                    mode='nearest',
                    truncate=4.0)

    # Output list:
    return aryFuncChnk
230637
def GaussianBlur(X, ksize_width, ksize_height, sigma_x, sigma_y):
    
    return image_transform(
        X,
        cv2.GaussianBlur,
        ksize=(ksize_width, ksize_height),
        sigmaX=sigma_x,
        sigmaY=sigma_y
    )
421361
def imagetransformer_sep_channels_8l():
  """"""
  hparams = imagetransformer_base()
  hparams.num_heads = 4
  hparams.attention_key_channels = hparams.attention_value_channels = 0
  hparams.hidden_size = 256
  hparams.filter_size = 256
  hparams.num_hidden_layers = 8
  hparams.sampling_method = "random"
  return hparams
1346867
public static void horizontalInverse(BorderIndex1D border , WlBorderCoef<WlCoef_F32> inverseCoef , GrayF32 input , GrayF32 output ) {

		UtilWavelet.checkShape(output,input);

		float []trends = new float[ output.width ];
		float []details = new float[ output.width ];

		boolean isLarger = input.width >= output.width;
		int paddedWidth = output.width + output.width%2;

		final int lowerBorder = inverseCoef.getLowerLength()*2;
		final int upperBorder = output.width - inverseCoef.getUpperLength()*2;

		border.setLength(output.width+output.width%2);

		WlCoef_F32 coefficients;

		for( int y = 0; y < output.height; y++ ) {

			for( int i = 0; i < details.length; i++ ) {
				details[i] = 0;
				trends[i] = 0;
			}

			for( int x = 0; x < output.width; x += 2 ) {
				float a = input.get(x/2,y);
				float d = input.get(input.width/2+x/2,y);

				if( x < lowerBorder ) {
					coefficients = inverseCoef.getBorderCoefficients(x);
				} else if( x >= upperBorder ) {
					coefficients = inverseCoef.getBorderCoefficients(x-paddedWidth);
				} else {
					coefficients = inverseCoef.getInnerCoefficients();
				}

				final int offsetA = coefficients.offsetScaling;
				final int offsetB = coefficients.offsetWavelet;
				final float[] alpha = coefficients.scaling;
				final float[] beta = coefficients.wavelet;

				// add the trend
				for( int i = 0; i < alpha.length; i++ ) {
					// if an odd image don't update the outer edge
					int xx = border.getIndex(x+offsetA+i);
					if( isLarger && xx >= output.width )
						continue;
					trends[xx] += a*alpha[i];
				}

				// add the detail signal
				for( int i = 0; i < beta.length; i++ ) {
					int xx = border.getIndex(x+offsetB+i);
					if( isLarger && xx >= output.width )
						continue;
					details[xx] += d*beta[i];
				}
			}

			for( int x = 0; x < output.width; x++ ) {
				output.set(x,y, trends[x] + details[x]);
			}
		}
	}
1233241
def _hotellings_fourier(self):
        
        sigma_2 = self._catalog_object.sigma**2
        # compute residual
        df = self._df #degrees of freedom
        R = self._A - np.dot(self._X, self._Bhat)
        R = np.matrix(R)
        # residual covariance matrix
        Sigma_R = (R.H*R) * (1./df)
        # noise covariance matrix
        #raise TypeError("fourier domain hypothesis tests not implemented yet")
        Zft = np.matrix(self._basis_object._Z)
        Sigma_S = Zft.H*np.matrix(np.diag(sigma_2))*Zft
        # add covariances (sum of multivariate normals is normal)
        Sigma_Z = np.array(Sigma_R + Sigma_S)

        T_2_list = []
        p_value_list = []
        for i in np.arange(0,self._Bhat.shape[0]):

            Bstar       = self._Bhat[i,np.arange(0,np.shape(self._A)[1])]
            lstar       = float(np.shape(Bstar)[0])
            cx          = self._Cx[i,i]
            Einv        = np.linalg.pinv(Sigma_Z)
            Zs          = Bstar/np.sqrt(cx)
            last_part   = np.array(np.matrix(Zs)*np.matrix(Einv)*np.matrix(Zs).H).real
            T_2         = ((df - lstar + 1.)/(df*lstar))*last_part
            p_value     = 1. - stats.f.cdf(T_2, 2.*lstar, 2.*(df - lstar + 1.))
            p_value_list.append(p_value)
            T_2_list.append(T_2)
        # save pvalue results in a table
        results = [['Comparison','Hotellings T^2', "p-value"]]
        for i in np.arange(0,len(self._col_names)):
            results.append([self._col_names[i], T_2_list[i], p_value_list[i]])
        self._results = results