B1 ratio for noise identification
    
        ratio of Standard Variace to AVAR
189413
def b1_noise_id(x, af, rate):
    
    (taus,devs,errs,ns) = at.adev(x,taus=[af*rate],data_type="phase", rate=rate) 
    oadev_x = devs[0]
    y = np.diff(x)
    y_cut = np.array( y[:len(y)-(len(y)%af)] ) # cut to length
    assert len(y_cut)%af == 0
    y_shaped = y_cut.reshape( ( int(len(y_cut)/af), af) )
    y_averaged = np.average(y_shaped,axis=1) # average
        
    var = np.var(y_averaged, ddof=1)
    return var/pow(oadev_x,2.0)
1374902
def rn_theory(af, b):
    
    # From IEEE1139-2008
    #   alpha   beta    ADEV_mu MDEV_mu Rn_mu
    #   -2      -4       1       1       0      Random Walk FM
    #   -1      -3       0       0       0      Flicker FM
    #    0      -2      -1      -1       0      White FM
    #    1      -1      -2      -2       0      Flicker PM
    #    2      0       -2      -3      -1      White PM
    
    # (a=-3 flicker walk FM)
    # (a=-4 random run FM)
    if b==0:
        return pow(af,-1)
    elif b==-1:
        # f_h = 0.5/tau0  (assumed!)
        # af = tau/tau0
        # so f_h*tau = 0.5/tau0 * af*tau0 = 0.5*af
        avar = (1.038+3*np.log(2*np.pi*0.5*af)) / (4.0*pow(np.pi,2))
        mvar = 3*np.log(256.0/27.0)/(8.0*pow(np.pi,2))
        return mvar/avar
    else:
        return pow(af,0)
1174182
def Kdiag(self, X):
        """"""
        X_prod = self._comp_prod(X)
        return self.variance*four_over_tau*np.arcsin(X_prod/(X_prod+1.))
506863
def obfn_g0var(self):
        

        return self.var_y0() if self.opt['AuxVarObj'] else \
            self.cnst_A0(self.X) - self.cnst_c0()
1237628
def obfn_reg(self):
        
        """
        l1 = np.sum(mp_wl1*np.abs(self.obfn_gvar()))
        return (self.lmbda*l1, l1)
628295
def _get_a1_value(bbar, dbar, slip, beta, mmax):
        
        return ((dbar - bbar) / dbar) * (slip / beta) *\
            np.exp(-(dbar / 2.) * mmax)
25605
def eval_objfn(self):
        
        """

        gvr = self.obfn_gvar()
        dfd = np.sum(np.abs(self.Wdf * gvr))
        reg = 0.5*np.linalg.norm(
            sl.idctii(self.Alpha*sl.dctii(self.X, axes=self.axes),
                      axes=self.axes))**2
        obj = dfd + self.lmbda*reg
        return (obj, dfd, reg)
910968
function autocorrelation(signal) {
        
        var n=signal.length;
        
        var i;
        
        var n2=nextpow2(n);

	if (n===0){
	    n2=0;
	}
        
        var paddedSignal=[];
        
        if (n==n2){
            paddedSignal=signal;
        }
        else{
            for (i=0;i<n2;i++){
                if (i<n){
                    paddedSignal[i]=signal[i];
                }
                else{
                    paddedSignal[i]=0;
                }
            }
        }

        var acv=autocovariance(paddedSignal);
        
        var variance=acv[0];
        
        var acf=[];
        
        var i;
        
        if (variance>0){
            for (i=0; i<n; i++){
                acf[i]=acv[i]/variance;    
            }
        }
        else{
            for (i=0; i<n; i++){
                acf[i]=0;    
            }            
        }
        return acf;
    
        function autocovariance(signal) {
            
            var n=signal.length;
            
            if (n===0)
                return [];
            
            var mean=0;
            var i;
            for (i=0; i<n; i++){
                mean+=signal[i];
            }
            mean/=n;
                      
            var padded_signal=[];
            
            for (i=0; i<n; i++){
                padded_signal[i]=signal[i]-mean;
                padded_signal[n+i]=0;
            }
                        
            var ft_signal = fft(padded_signal);
                       
            var pseudo_powerSpectralDensity=[];
            
            for (i=0; i<ft_signal.length; i++){
                pseudo_powerSpectralDensity[i]=complexMultiply(ft_signal[i], complexConjugate(ft_signal[i]));
            }
           
            var pseudo_autocovariance = ifft(pseudo_powerSpectralDensity);
    
            var mask=[];
            
            for (i=0; i<n; i++){
                mask[i]=1;
                mask[i+n]=0;
            }
            
            var ft_mask = fft(mask);
    
            var mask_powerSpectralDensity=[];
            
            for (i=0; i<ft_mask.length; i++){
                mask_powerSpectralDensity[i]=complexMultiply(ft_mask[i], complexConjugate(ft_mask[i]));
            }
            
            var mask_correction_factors = ifft(mask_powerSpectralDensity);
    
            var acv=[];
            
            for (i=0; i<n; i++){
                acv[i]=complexDivideRealPart(pseudo_autocovariance[i],mask_correction_factors[i]);
            }    
            
            return acv;
            /*
             * ComplexMultiply
             * 
             * From complex.js in fft-js in https://github.com/vail-systems/node-fft
             *
             * (c) Vail Systems. Joshua Jung and Ben Bryan. 2015
             */
            function complexMultiply(a, b){
                return [(a[0] * b[0] - a[1] * b[1]), 
                        (a[0] * b[1] + a[1] * b[0])];
            }
            
            function complexConjugate(a){
                return [a[0], -a[1]];
            }
            /*function complexDivide(a,b){
                return [ (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]),
                         (a[1]*b[0]-a[0]*b[1])/(b[0]*b[0]+b[1]*b[1])];
            }*/
            //We only need the real part here
            function complexDivideRealPart(a,b){
                return (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]);
            }           
            
        }
    }
915354
def noise(params, amplitude=1, offset=0):
    '''
    
    '''
    amplitude = create_buffer(params, amplitude)
    offset = create_buffer(params, offset)
    output = offset + amplitude * (np.random.random(params.length)*2 - 1)
    return output
859239
def vtas2mach(tas, h):
    
    a = vvsound(h)
    M = tas / a
    return M
531223
function estimateBandwidth(array) {
  var n = array.length,
      q = quartiles(array),
      h = (q[2] - q[0]) / 1.34;
  return 1.06 * Math.min(Math.sqrt(variance(array)), h) * Math.pow(n, -0.2);
}
213456
def generateNoise(self):
        
        # Fill wfb array with white noise based on given discrete variance
        wfb = np.zeros(self.nr*2)
        wfb[:self.nr] = np.random.normal(0, np.sqrt(self.qd), self.nr)
        # Generate the hfb coefficients based on the noise type
        mhb = -self.b/2.0
        hfb = np.zeros(self.nr*2)
        hfb = np.zeros(self.nr*2)
        hfb[0] = 1.0
        indices = np.arange(self.nr-1)
        hfb[1:self.nr] = (mhb+indices)/(indices+1.0)
        hfb[:self.nr] = np.multiply.accumulate(hfb[:self.nr])
        # Perform discrete Fourier transform of wfb and hfb time series
        wfb_fft = np.fft.rfft(wfb)
        hfb_fft = np.fft.rfft(hfb)
        # Perform inverse Fourier transform of the product of wfb and hfb FFTs
        time_series = np.fft.irfft(wfb_fft*hfb_fft)[:self.nr]
        self.time_series = time_series
376700
function analyzeNumericArray( v, include_costly ) {
	var min = v[0], max = min, is_int = false, is_float = false, is_same = true,
		dmin=0, dmax=0, is_dfloat = false,
		mean=0, n_type=0, d_mode=0, f_type=[0, NUMTYPE.UNKNOWN], 
		c_same=0, same=0,
		s_min = [min,min,min,min,min], s_min_i=0, s_max = [min,min,min,min,min], s_max_i=0, samples,
		a, b, d_type, cd, cv, lv = v[0];

	// Anlyze array items
	for (var i=0, l=v.length; i<l; ++i) {
		cv = v[i];

		// Exit on non-numeric cases
		if (typeof cv !== 'number')
			return null;

		// Update mean
		mean += cv;

		// Include costly calculations if enabled
		if (include_costly) {

			// Update delta
			if (i !== 0) {
				cd = lv - cv; if (cd < 0) cd = -cd;
				if (i === 1) {
					dmin = cd;
					dmax = cd;
				} else {
					if (cd < dmin) dmin = cd;
					if (cd > dmax) dmax = cd;
				}

				// Check if delta is float
				if ((cd !== 0) && (cd % 1 !== 0))
					is_dfloat = true;
			}

			// Update bounds & Keep samples
			if (cv < min) {
				min = cv;
				s_min[s_min_i] = cv;
				if (++s_min_i>5) s_min_i=0;
			}
			if (cv > max) {
				max = cv;
				s_max[s_max_i] = cv;
				if (++s_max_i>5) s_max_i=0;
			}

		} else {

			// Update bounds
			if (cv < min) min = cv;
			if (cv > max) max = cv;

		}

		// Check for same values
		if (cv === lv) {
			c_same++;
		} else {
			if (c_same > same)
				same = c_same;
			is_same = false;
			c_same = 0;
		}

		// Keep last value
		lv = cv;

		// Skip zeros from further analysis
		if (cv === 0) continue;

		// Update integer/float
		if (cv % 1 !== 0) {
			is_float = true;
		} else {
			is_int = true;
		}

	}

	// Finalize counters
	if (c_same > same) same = c_same;
	mean /= v.length;

	// Guess numerical type
	n_type = getNumType( min, max, is_float );

	// Calculate delta-encoding details
	d_type = NUMTYPE.UNKNOWN;
	if (include_costly) {
		if (!is_float && is_int) {

			// INTEGERS : Use Delta-Encoding (d_mode=1)

			if ((dmin > INT8_MIN) && (dmax < INT8_MAX)) {
				d_type = NUMTYPE.INT8;
				d_mode = 1;
			} else if ((dmin > INT16_MIN) && (dmax < INT16_MAX)) {
				d_type = NUMTYPE.INT16;
				d_mode = 1;
			} else if ((dmin > INT32_MIN) && (dmax < INT32_MAX)) {
				d_type = NUMTYPE.INT32;
				d_mode = 1;
			}

		} else if (is_float) {

			// FLOATS : Use Rebase Encoding (d_mode=2)

			// Get a couple more samples
			samples = [].concat(
				s_min, s_max,
				[
					v[Math.floor(Math.random()*v.length)],
					v[Math.floor(Math.random()*v.length)],
					v[Math.floor(Math.random()*v.length)],
					v[Math.floor(Math.random()*v.length)],
					v[Math.floor(Math.random()*v.length)],
				]
			);

			// Calculate float scale
			f_type = getFloatScale( v, min, max, 0.01 );
			if (f_type[1] != NUMTYPE.UNKNOWN) {
				d_type = f_type[1];
				d_mode = 2;
			}

		}
	}

	// Based on information detected so far, populate
	// the analysis results
	return {

		// Get numeric type
		'type'		: n_type,
		'delta_type': d_type,

		// Percentage of same items
		'psame'		: same / v.length,

		// Log numeric bounds
		'min'		: min,
		'max'		: max,
		'mean'		: mean,

		// Log delta bounds
		'dmin'		: dmin,
		'dmax' 		: dmax,

		// Delta mode
		'dmode'		: d_mode,
		'fscale' 	: f_type[0],

		// Expose information details
		'integer' 	: is_int && !is_float,
		'float' 	: is_float && !is_int,
		'mixed' 	: is_float && is_int,
		'same'		: is_same && (v.length > 1),

	};

}
1095810
def _margtimedist_loglr(self, mf_snr, opt_snr):
        
        logl = special.logsumexp(mf_snr, b=self._deltat)
        logl_marg = logl/self._dist_array
        opt_snr_marg = opt_snr/self._dist_array**2
        return special.logsumexp(logl_marg - 0.5*opt_snr_marg,
                                 b=self._deltad*self.dist_prior)
579107
def obfn_g0var(self):
        

        # Use of self.block_sep0(self.AXnr) instead of self.cnst_A0(self.X)
        # reduces number of calls to self.cnst_A0
        return self.var_y0() if self.opt['gEvalY'] else \
            self.block_sep0(self.AXnr)
508718
def _batch_norm(self, name, x):
        """"""
        with tf.variable_scope(name):
            params_shape = [x.get_shape()[-1]]

            beta = tf.get_variable(
                "beta",
                params_shape,
                tf.float32,
                initializer=tf.constant_initializer(0.0, tf.float32))
            gamma = tf.get_variable(
                "gamma",
                params_shape,
                tf.float32,
                initializer=tf.constant_initializer(1.0, tf.float32))

            if self.mode == "train":
                mean, variance = tf.nn.moments(x, [0, 1, 2], name="moments")

                moving_mean = tf.get_variable(
                    "moving_mean",
                    params_shape,
                    tf.float32,
                    initializer=tf.constant_initializer(0.0, tf.float32),
                    trainable=False)
                moving_variance = tf.get_variable(
                    "moving_variance",
                    params_shape,
                    tf.float32,
                    initializer=tf.constant_initializer(1.0, tf.float32),
                    trainable=False)

                self._extra_train_ops.append(
                    moving_averages.assign_moving_average(
                        moving_mean, mean, 0.9))
                self._extra_train_ops.append(
                    moving_averages.assign_moving_average(
                        moving_variance, variance, 0.9))
            else:
                mean = tf.get_variable(
                    "moving_mean",
                    params_shape,
                    tf.float32,
                    initializer=tf.constant_initializer(0.0, tf.float32),
                    trainable=False)
                variance = tf.get_variable(
                    "moving_variance",
                    params_shape,
                    tf.float32,
                    initializer=tf.constant_initializer(1.0, tf.float32),
                    trainable=False)
                tf.summary.histogram(mean.op.name, mean)
                tf.summary.histogram(variance.op.name, variance)
            # elipson used to be 1e-5. Maybe 0.001 solves NaN problem in deeper
            # net.
            y = tf.nn.batch_normalization(x, mean, variance, beta, gamma,
                                          0.001)
            y.set_shape(x.get_shape())
            return y
741590
def Kdiag(self, X):
        """"""
        K_ff = np.ones(X.shape[0])
        for i,x in enumerate(X):
            for il,l in enumerate(self.lengthscale):
                idx = il*2
                K_ff[i] *= self.k_ff(x[idx],x[idx],l)
        return K_ff * self.variances[0]