Random resampling of the doublegauss distribution
357596
def resample(self,N,**kwargs):
        
        lovals = self.mu - np.absolute(rand.normal(size=N)*self.siglo)
        hivals = self.mu + np.absolute(rand.normal(size=N)*self.sighi)

        u = rand.random(size=N)
        hi = (u < float(self.sighi)/(self.sighi + self.siglo))
        lo = (u >= float(self.sighi)/(self.sighi + self.siglo))

        vals = np.zeros(N)
        vals[hi] = hivals[hi]
        vals[lo] = lovals[lo]
        return vals
1340858
def samples(self, gp, Y_metadata):
        
        N1, N2 = gp.shape
        Ysim = np.zeros((N1,N2))
        ind = Y_metadata['output_index'].flatten()
        for j in np.unique(ind):
            flt = ind==j
            gp_filtered = gp[flt,:]
            n1 = gp_filtered.shape[0]
            lik = self.likelihoods_list[j]
            _ysim = np.array([np.random.normal(lik.gp_link.transf(gpj), scale=np.sqrt(lik.variance), size=1) for gpj in gp_filtered.flatten()])
            Ysim[flt,:] = _ysim.reshape(n1,N2)
        return Ysim
1146495
def normalizedGraylevelVariance(img):
    '''''''
    mean, stdev = cv2.meanStdDev(img)
    s = stdev[0]**2 / mean[0]
    return s[0]
1244856
private void createSamples( double sampleRadius , Point2D_F64 workPoint , Point2D_F64 originalPoint ) {

		workPoint.x = originalPoint.x + sampleRadius;
		if( estimator.process(work,false) ) {
			samples.grow().set( estimator.getWorldToCamera() );
		}
		workPoint.x = originalPoint.x - sampleRadius;
		if( estimator.process(work,false) ) {
			samples.grow().set( estimator.getWorldToCamera() );
		}
		workPoint.x = originalPoint.x;

		workPoint.y = originalPoint.y + sampleRadius;
		if( estimator.process(work,false) ) {
			samples.grow().set( estimator.getWorldToCamera() );
		}
		workPoint.y = originalPoint.y - sampleRadius;
		if( estimator.process(work,false) ) {
			samples.grow().set( estimator.getWorldToCamera() );
		}
		workPoint.set(originalPoint);
	}
459381
def sample_image(d, data, u, v, w, i=-1, verbose=0, imager='xy', wres=100):
    

    if i == -1:
        i = len(data)/2

    if imager == 'xy':
        image = rtlib.imgonefullxy(n.outer(u, d['freq']/d['freq_orig'][0]), n.outer(v, d['freq']/d['freq_orig'][0]), data[i], d['npixx'], d['npixy'], d['uvres'], verbose=verbose)
    elif imager == 'w':
        npix = max(d['npixx'], d['npixy'])
        bls, uvkers = rtlib.genuvkernels(w, wres, npix, d['uvres'], ksize=21, oversample=1)
        image = rtlib.imgonefullw(n.outer(u, d['freq']/d['freq_orig'][0]), n.outer(v, d['freq']/d['freq_orig'][0]), data[i], npix, d['uvres'], bls, uvkers, verbose=verbose)

#        bls, lmkers = rtlib.genlmkernels(w, wres, npix, d['uvres'])
#        image = rtlib.imgonefullw(n.outer(u, d['freq']/d['freq_orig'][0]), n.outer(v, d['freq']/d['freq_orig'][0]), data[i], npix, d['uvres'], [bls[0]], [lmkers[0]], verbose=verbose)

    return image
899992
def samples(self, gp, Y_metadata=None, **kw):
        
        orig_shape = gp.shape
        gp = gp.flatten()
        N = Y_metadata['trials']
        Ysim = np.random.binomial(N, self.gp_link.transf(gp))
        return Ysim.reshape(orig_shape)
697893
def sample_W(self, nSamples, raw_samples=False):
        
        assert isinstance(self.kern, kern.Linear)
        from ..util.linalg import pdinv
        N, D = self.Y.shape
        Q = self.X.shape[1]
        noise_var = self.likelihood.variance.values
        
        # Draw samples for X
        Xs = np.random.randn(*((nSamples,)+self.X.shape))*np.sqrt(self.X.variance.values)+self.X.mean.values
        b = np.random.rand(*((nSamples,)+self.X.shape))
        Xs[b>self.X.gamma.values] = 0
        
        invcov = (Xs[:,:,:,None]*Xs[:,:,None,:]).sum(1)/noise_var+np.eye(Q)
        cov = np.array([pdinv(invcov[s_idx])[0] for s_idx in range(invcov.shape[0])])
        Ws = np.empty((nSamples, Q, D))
        tmp = (np.transpose(Xs, (0,2,1)).reshape(nSamples*Q,N).dot(self.Y)).reshape(nSamples,Q,D)
        mean = (cov[:,:,:,None]*tmp[:,None,:,:]).sum(2)/noise_var
        zeros = np.zeros((Q,))
        for s_idx in range(Xs.shape[0]):
            Ws[s_idx] = (np.random.multivariate_normal(mean=zeros,cov=cov[s_idx],size=(D,))).T+mean[s_idx]
        
        if raw_samples:
            return Ws
        else:
            return Ws.mean(0), Ws.std(0)
949733
function resample_hermite(canvas, W, H, W2, H2){
	var time1 = Date.now();
	W2 = Math.round(W2);
	H2 = Math.round(H2);
	var img = canvas.getContext("2d").getImageData(0, 0, W, H);
	var img2 = canvas.getContext("2d").getImageData(0, 0, W2, H2);
	var data = img.data;
	var data2 = img2.data;
	var ratio_w = W / W2;
	var ratio_h = H / H2;
	var ratio_w_half = Math.ceil(ratio_w/2);
	var ratio_h_half = Math.ceil(ratio_h/2);
	
	for(var j = 0; j < H2; j++){
		for(var i = 0; i < W2; i++){
			var x2 = (i + j*W2) * 4;
			var weight = 0;
			var weights = 0;
			var weights_alpha = 0;
			var gx_r = gx_g = gx_b = gx_a = 0;
			var center_y = (j + 0.5) * ratio_h;
			for(var yy = Math.floor(j * ratio_h); yy < (j + 1) * ratio_h; yy++){
				var dy = Math.abs(center_y - (yy + 0.5)) / ratio_h_half;
				var center_x = (i + 0.5) * ratio_w;
				var w0 = dy*dy //pre-calc part of w
				for(var xx = Math.floor(i * ratio_w); xx < (i + 1) * ratio_w; xx++){
					var dx = Math.abs(center_x - (xx + 0.5)) / ratio_w_half;
					var w = Math.sqrt(w0 + dx*dx);
					if(w >= -1 && w <= 1){
						//hermite filter
						weight = 2 * w*w*w - 3*w*w + 1;
						if(weight > 0){
							dx = 4*(xx + yy*W);
							//alpha
							gx_a += weight * data[dx + 3];
							weights_alpha += weight;
							//colors
							if(data[dx + 3] < 255)
								weight = weight * data[dx + 3] / 250;
							gx_r += weight * data[dx];
							gx_g += weight * data[dx + 1];
							gx_b += weight * data[dx + 2];
							weights += weight;
							}
						}
					}		
				}
			data2[x2]     = gx_r / weights;
			data2[x2 + 1] = gx_g / weights;
			data2[x2 + 2] = gx_b / weights;
			data2[x2 + 3] = gx_a / weights_alpha;
			}
		}
	log( "pong-io","hermite = "+(Math.round(Date.now() - time1)/1000)+" s");
	canvas.getContext("2d").clearRect(0, 0, Math.max(W, W2), Math.max(H, H2));
	canvas.width = W2;
	canvas.height = H2;
	canvas.getContext("2d").putImageData(img2, 0, 0);
}
46975
def psd(ts, nperseg=1500, noverlap=1200, plot=True):
    """"""
    ts = ts.squeeze()
    if ts.ndim is 1:
        ts = ts.reshape((-1, 1))
    fs = (len(ts) - 1.0) / (ts.tspan[-1] - ts.tspan[0])
    window = signal.hamming(nperseg, sym=False)
    nfft = max(256, 2**np.int(np.log2(nperseg) + 1))
    freqs, pxx = signal.welch(ts, fs, window, nperseg, noverlap, nfft,
                              detrend='linear', axis=0)
    # Discard estimates for freq bins that are too low for the window size.
    # (require two full cycles to fit within the window)
    index = np.nonzero(freqs >= 2.0*fs/nperseg)[0][0]
    if index > 0:
        freqs = freqs[index:]
        pxx = pxx[index:]
    # Discard estimate for last freq bin as too high for Nyquist frequency:
    freqs = freqs[:-1]
    pxx = pxx[:-1]
    if plot is True:
        _plot_psd(ts, freqs, pxx)
    return freqs, pxx
956899
function(distribution, coverageAlpha) {
        if (distribution.numTrials < 1000) {
            // don't even bother trying to optimize for small-ish sample sizes
            return [0, distribution.numTrials];
        } else {
            return [
                Math.floor(distribution.inverseCdf(coverageAlpha / 2)),
                Math.ceil(distribution.inverseSurvival(coverageAlpha / 2))
            ];
        }
    }
12313
def generate_gaussian_profile(seeing_fwhm):
    """"""
    FWHM_G = 2 * math.sqrt(2 * math.log(2))
    sigma = seeing_fwhm / FWHM_G
    amplitude = 1.0 / (2 * math.pi * sigma * sigma)
    seeing_model = Gaussian2D(amplitude=amplitude,
                              x_mean=0.0,
                              y_mean=0.0,
                              x_stddev=sigma,
                              y_stddev=sigma)
    return seeing_model
383895
def samples(self, gp, Y_metadata=None):
        
        orig_shape = gp.shape
        gp = gp.flatten()
        #FIXME: Very slow as we are computing a new random variable per input!
        #Can't get it to sample all at the same time
        #student_t_samples = np.array([stats.t.rvs(self.v, self.gp_link.transf(gpj),scale=np.sqrt(self.sigma2), size=1) for gpj in gp])
        dfs = np.ones_like(gp)*self.v
        scales = np.ones_like(gp)*np.sqrt(self.sigma2)
        student_t_samples = stats.t.rvs(dfs, loc=self.gp_link.transf(gp),
                                        scale=scales)
        return student_t_samples.reshape(orig_shape)
1219160
def _resample_obspy(samples, sr, newsr, window='hanning', lowpass=True):
    # type: (np.ndarray, int, int, str, bool) -> np.ndarray
    
    from scipy.signal import resample
    from math import ceil
    factor = sr/float(newsr)
    if newsr < sr and lowpass:
        # be sure filter still behaves good
        if factor > 16:
            logger.info("Automatic filter design is unstable for resampling "
                        "factors (current sampling rate/new sampling rate) " 
                        "above 16. Manual resampling is necessary.")
        freq = min(sr, newsr) * 0.5 / float(factor)
        logger.debug(f"resample_obspy: lowpass {freq}")
        samples = lowpass_cheby2(samples, freq=freq, sr=sr, maxorder=12)
    num = int(ceil(len(samples) / factor))

    return _applyMultichan(samples, 
                           lambda S: resample(S, num, window=window))
1013889
def interpoled_resampling(W, x):
    
    N = W.shape[0] 
    idx = np.argsort(x)
    xs = x[idx] 
    ws = W[idx]
    cs = np.cumsum(avg_n_nplusone(ws)) 
    u = random.rand(N)
    xrs = np.empty(N)
    where = np.searchsorted(cs, u)
    # costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
    for n in range(N): 
        m = where[n]
        if m==0:
            xrs[n] = xs[0]
        elif m==N:
            xrs[n] = xs[-1]
        else:
            xrs[n] = interpol(cs[m-1], cs[m], xs[m-1], xs[m], u[n])
    return xrs
390811
def sample_double_norm(mean, std_upper, std_lower, size):
    """"""
    from scipy.special import erfinv

    # There's probably a better way to do this. We first draw percentiles
    # uniformly between 0 and 1. We want the peak of the distribution to occur
    # at `mean`. However, if we assign 50% of the samples to the lower half
    # and 50% to the upper half, the side with the smaller variance will be
    # overrepresented because of the 1/sigma normalization of the Gaussian
    # PDF. Therefore we need to divide points between the two halves with a
    # fraction `cutoff` (defined below) going to the lower half. Having
    # partitioned them this way, we can then use the standard Gaussian
    # quantile function to go from percentiles to sample values -- except that
    # we must remap from [0, cutoff] to [0, 0.5] and from [cutoff, 1] to [0.5,
    # 1].

    samples = np.empty(size)
    percentiles = np.random.uniform(0., 1., size)
    cutoff = std_lower / (std_lower + std_upper)

    w = (percentiles < cutoff)
    percentiles[w] *= 0.5 / cutoff
    samples[w] = mean + np.sqrt(2) * std_lower * erfinv(2 * percentiles[w] - 1)

    w = ~w
    percentiles[w] = 1 - (1 - percentiles[w]) * 0.5 / (1 - cutoff)
    samples[w] = mean + np.sqrt(2) * std_upper * erfinv(2 * percentiles[w] - 1)

    return samples
349006
function GraphicalEq(sampleRate) {
  this.FS = sampleRate;
  this.minFreq = 40.0;
  this.maxFreq = 16000.0;

  this.bandsPerOctave = 1.0;

  this.filters = [];
  this.freqzs = [];

  this.calculateFreqzs = true;

  this.recalculateFilters = function() {
    var bandCount = Math.round(Math.log(this.maxFreq/this.minFreq) * this.bandsPerOctave/ Math.LN2);

    this.filters = [];
    for (var i=0; i<bandCount; i++) {
      var freq = this.minFreq*(Math.pow(2, i/this.bandsPerOctave));
      var newFilter = new Biquad(DSP.PEAKING_EQ, this.FS);
      newFilter.setDbGain(0);
      newFilter.setBW(1/this.bandsPerOctave);
      newFilter.setF0(freq);
      this.filters[i] = newFilter;
      this.recalculateFreqz(i);
    }
  };

  this.setMinimumFrequency = function(freq) {
    this.minFreq = freq;
    this.recalculateFilters();
  };

  this.setMaximumFrequency = function(freq) {
    this.maxFreq = freq;
    this.recalculateFilters();
  };

  this.setBandsPerOctave = function(bands) {
    this.bandsPerOctave = bands;
    this.recalculateFilters();
  };

  this.setBandGain = function(bandIndex, gain) {
    if (bandIndex < 0 || bandIndex > (this.filters.length-1)) {
      throw "The band index of the graphical equalizer is out of bounds.";
    }

    if (!gain) {
      throw "A gain must be passed.";
    }
   
    this.filters[bandIndex].setDbGain(gain);
    this.recalculateFreqz(bandIndex);
  };
 
  this.recalculateFreqz = function(bandIndex) {
    if (!this.calculateFreqzs) {
      return;
    }

    if (bandIndex < 0 || bandIndex > (this.filters.length-1)) {
      throw "The band index of the graphical equalizer is out of bounds. " + bandIndex + " is out of [" + 0 + ", " + this.filters.length-1 + "]";
    }
       
    if (!this.w) {
      this.w = Float32Array(400);
      for (var i=0; i<this.w.length; i++) {
         this.w[i] = Math.PI/this.w.length * i;
      }
    }
   
    var b = [this.filters[bandIndex].b0, this.filters[bandIndex].b1, this.filters[bandIndex].b2];
    var a = [this.filters[bandIndex].a0, this.filters[bandIndex].a1, this.filters[bandIndex].a2];

    this.freqzs[bandIndex] = DSP.mag2db(DSP.freqz(b, a, this.w));
  };

  this.process = function(buffer) {
    var output = buffer;

    for (var i = 0; i < this.filters.length; i++) {
      output = this.filters[i].process(output);
    }

    return output;
  };

  this.processStereo = function(buffer) {
    var output = buffer;

    for (var i = 0; i < this.filters.length; i++) {
      output = this.filters[i].processStereo(output);
    }

    return output;
  };
}
1046122
def noise2(self, x, y):
		
		# Skew input space to determine which simplex (triangle) we are in
		s = (x + y) * _F2
		i = floor(x + s)
		j = floor(y + s)
		t = (i + j) * _G2
		x0 = x - (i - t) # "Unskewed" distances from cell origin
		y0 = y - (j - t)

		if x0 > y0:
			i1 = 1; j1 = 0 # Lower triangle, XY order: (0,0)->(1,0)->(1,1)
		else:
			i1 = 0; j1 = 1 # Upper triangle, YX order: (0,0)->(0,1)->(1,1)
		
		x1 = x0 - i1 + _G2 # Offsets for middle corner in (x,y) unskewed coords
		y1 = y0 - j1 + _G2
		x2 = x0 + _G2 * 2.0 - 1.0 # Offsets for last corner in (x,y) unskewed coords
		y2 = y0 + _G2 * 2.0 - 1.0

		# Determine hashed gradient indices of the three simplex corners
		perm = self.permutation
		ii = int(i) % self.period
		jj = int(j) % self.period
		gi0 = perm[ii + perm[jj]] % 12
		gi1 = perm[ii + i1 + perm[jj + j1]] % 12
		gi2 = perm[ii + 1 + perm[jj + 1]] % 12

		# Calculate the contribution from the three corners
		tt = 0.5 - x0**2 - y0**2
		if tt > 0:
			g = _GRAD3[gi0]
			noise = tt**4 * (g[0] * x0 + g[1] * y0)
		else:
			noise = 0.0
		
		tt = 0.5 - x1**2 - y1**2
		if tt > 0:
			g = _GRAD3[gi1]
			noise += tt**4 * (g[0] * x1 + g[1] * y1)
		
		tt = 0.5 - x2**2 - y2**2
		if tt > 0:
			g = _GRAD3[gi2]
			noise += tt**4 * (g[0] * x2 + g[1] * y2)

		return noise * 70.0