Given state dict and noise properties, estimate flux scale at the VLA 

    imstd and flagfrac are expected to be median (typical) values from sample in merged noise pkl.
353788
def calcfluxscale(d, imstd_med, flagfrac_med):
    

    # useful functions and VLA parameters
    sensitivity = lambda sefd, dt, bw, eta, nbl, npol: sefd/(eta*np.sqrt(nbl*2 * dt * bw * npol))
    nbl = lambda nant: nant*(nant-1)/2
    eta = {'L': 0.92, 'S': 0.92, 'C': 0.8, 'X': 0.8}   # correlator efficiency
    sefd = {'L': 420, 'S': 370, 'C': 310, 'X': 250}    # fixed to match exposure calculator int time to 100 microJy.

    bw = sum([d['spw_nchan_select'][i]*d['spw_chansize'][i] for i in range(len(d['spw_chansize']))])
    dt = d['inttime']
    npol = d['npol']
    nant = d['nants']
    freq = d['freq'][0]
    if (freq >= 1 and freq < 2):
        band = 'L'
    elif (freq >= 2 and freq < 4):
        band = 'S'
    elif (freq >= 4 and freq < 8):
        band = 'C'
    elif (freq >= 8 and freq < 12):
        band = 'X'
    else:
        logger.warn('first channel freq ({0}) not in bands L, S, C, or X. Assuming L band.'.format(freq))
        band = 'L'

    goodfrac = 1 - flagfrac_med # correct for flagged data
    slim_theory = sensitivity(sefd[band], dt, bw, eta[band], goodfrac*nbl(nant), npol)
    fluxscale = slim_theory/imstd_med

    return fluxscale
206592
def factorized_gaussian_noise(in_features, out_features, device):
    
    in_noise = scaled_noise(in_features, device=device)
    out_noise = scaled_noise(out_features, device=device)

    return out_noise.ger(in_noise), out_noise
676991
def estimate_lambda(pv):
    """"""
    LOD2 = sp.median(st.chi2.isf(pv,1))
    L = (LOD2/0.456)
    return (L)
156222
def bootstrap_noise(data, func, n=10000, std=1, symmetric=True):
    
    boot_dist = []
    arr = N.zeros(data.shape)
    for i in range(n):
        if symmetric:
            # Noise on all three axes
            arr = N.random.randn(*data.shape)*std
        else:
            # Only z-coordinate noise
            arr[:,-1] = N.random.randn(data.shape[0])*std
        boot_dist.append(func(data+arr))
    return N.array(boot_dist)
684323
def noise_uniform(self, lower_bound, upper_bound):
        """"""
        assert upper_bound > lower_bound
        nu = self.sym.sym('nu_{:d}'.format(len(self.scope['nu'])))
        self.scope['nu'].append(nu)
        return lower_bound + nu*(upper_bound - lower_bound)
1286802
def tf_smooth_eig_vec(self):
    """"""
    _, matrix_m = self.dual_object.get_full_psd_matrix()
    # Easier to think in terms of max so negating the matrix
    [eig_vals, eig_vectors] = tf.self_adjoint_eig(-matrix_m)
    exp_eig_vals = tf.exp(tf.divide(eig_vals, self.smooth_placeholder))
    scaling_factor = tf.reduce_sum(exp_eig_vals)
    # Multiplying each eig vector by exponential of corresponding eig value
    # Scaling factor normalizes the vector to be unit norm
    eig_vec_smooth = tf.divide(
        tf.matmul(eig_vectors, tf.diag(tf.sqrt(exp_eig_vals))),
        tf.sqrt(scaling_factor))
    return tf.reshape(
        tf.reduce_sum(eig_vec_smooth, axis=1),
        shape=[eig_vec_smooth.shape[0].value, 1])
1348892
def run_snr(self):
        

        if self.ecc:
            required_kwargs = {'dist_type': self.dist_type,
                               'initial_cond_type': self.initial_cond_type,
                               'ecc': True}
            input_args = [self.m1, self.m2, self.z_or_dist, self.initial_point,
                          self.eccentricity, self.observation_time]

        else:
            required_kwargs = {'dist_type': self.dist_type}
            input_args = [self.m1, self.m2, self.spin_1, self.spin_2,
                          self.z_or_dist, self.start_time, self.end_time]

        input_kwargs = {**required_kwargs,
                        **self.general,
                        **self.sensitivity_input,
                        **self.snr_input,
                        **self.parallel_input}

        self.final_dict = snr(*input_args, **input_kwargs)
        return
789568
def _cont_norm_running_quantile_mp(wl, fluxes, ivars, q, delta_lambda,
                                   n_proc=2, verbose=False):
    
    nStar = fluxes.shape[0]

    # start mp.Pool
    mp_results = []
    pool = mp.Pool(processes=n_proc)
    for i in xrange(nStar):
        mp_results.append(pool.apply_async(\
            _find_cont_running_quantile,
            (wl, fluxes[i, :].reshape((1, -1)), ivars[i, :].reshape((1, -1)),
             q, delta_lambda), {'verbose': False}))
        if verbose:
            print('@Bo Zhang: continuum normalizing star [%d/%d] ...'\
                  % (i + 1, nStar))
    # close mp.Pool
    pool.close()
    pool.join()

    # reshape results --> cont
    cont = np.zeros_like(fluxes)
    for i in xrange(nStar):
        cont[i, :] = mp_results[i].get() #.flatten()
    norm_fluxes = np.ones(fluxes.shape)
    norm_fluxes[cont!=0] = fluxes[cont!=0] / cont[cont!=0]
    norm_ivars = cont**2 * ivars

    print('@Bo Zhang: continuum normalization finished!')
    return norm_fluxes, norm_ivars
658908
def _build_index_param(self, n_l, n_V, n_smooth):
        
        idx_param_sing = {'Cholesky': np.arange(n_l), 'a1': n_l}
        # for simplified fitting
        idx_param_fitU = {'Cholesky': np.arange(n_l),
                          'a1': np.arange(n_l, n_l + n_V)}
        # for the likelihood function when we fit U (the shared covariance).
        idx_param_fitV = {'log_SNR2': np.arange(n_V - 1),
                          'c_space': n_V - 1, 'c_inten': n_V,
                          'c_both': np.arange(n_V - 1, n_V - 1 + n_smooth)}
        # for the likelihood function when we fit V (reflected by SNR of
        # each voxel)
        return idx_param_sing, idx_param_fitU, idx_param_fitV
970897
def estimate_band_connection(prev_eigvecs, eigvecs, prev_band_order):
    
    metric = np.abs(np.dot(prev_eigvecs.conjugate().T, eigvecs))
    connection_order = []
    for overlaps in metric:
        maxval = 0
        for i in reversed(range(len(metric))):
            val = overlaps[i]
            if i in connection_order:
                continue
            if val > maxval:
                maxval = val
                maxindex = i
        connection_order.append(maxindex)

    band_order = [connection_order[x] for x in prev_band_order]

    return band_order
644642
def compute_vis_noise(msname, sefd, spw_id=0):
    """"""
    from pyrap.tables import table

    tab = table(msname)
    spwtab = table(msname + "/SPECTRAL_WINDOW")

    freq0 = spwtab.getcol("CHAN_FREQ")[spw_id, 0]
    wavelength = 300e+6/freq0
    bw = spwtab.getcol("CHAN_WIDTH")[spw_id, 0]
    dt = tab.getcol("EXPOSURE", 0, 1)[0]
    dtf = (tab.getcol("TIME", tab.nrows()-1, 1)-tab.getcol("TIME", 0, 1))[0]

    # close tables properly, else the calls below will hang waiting for a lock...
    tab.close()
    spwtab.close()

    print(">>> %s freq %.2f MHz (lambda=%.2fm), bandwidth %.2g kHz, %.2fs integrations, %.2fh synthesis"%(msname, freq0*1e-6, wavelength, bw*1e-3, dt, dtf/3600))
    noise = sefd/math.sqrt(abs(2*bw*dt))
    print(">>> SEFD of %.2f Jy gives per-visibility noise of %.2f mJy"%(sefd, noise*1000))

    return noise
379698
def from_dict(d):
        
        return NoiseModel(
            gates=[KrausModel.from_dict(t) for t in d["gates"]],
            assignment_probs={int(qid): np.array(a) for qid, a in d["assignment_probs"].items()},
        )
12058
def corruptSparseVector(sdr, noiseLevel):
  

  numNoiseBits = int(noiseLevel * np.sum(sdr))
  if numNoiseBits <= 0:
    return sdr
  activeBits = np.where(sdr > 0)[0]
  inActiveBits = np.where(sdr == 0)[0]

  turnOffBits = np.random.permutation(activeBits)
  turnOnBits = np.random.permutation(inActiveBits)
  turnOffBits = turnOffBits[:numNoiseBits]
  turnOnBits = turnOnBits[:numNoiseBits]

  sdr[turnOffBits] = 0
  sdr[turnOnBits] = 1
1046152
def _gaussian_membership_sigma(smoothness, eps = 0.0005): # 275us @ smothness=10
    
    """
    error = 0
    deltas = [0.1, 0.01, 0.001, 0.0001]
    sigma = smoothness * 0.3
    point = -1. * (smoothness + 0.5)
    for delta in deltas:
        while error < eps:
            sigma += delta
            error = scipy.stats.norm.cdf(0.5, point, sigma) - scipy.stats.norm.cdf(-0.5, point, sigma) # x, mu, sigma
        sigma -= delta
    return sigma
1141558
def time_col_turbulent(EnergyDis, ConcAl, ConcClay, coag, material,
                       DiamTarget, DIM_FRACTAL):
    
    return((1/6) * (6/np.pi)**(1/9) * EnergyDis**(-1/3) * DiamTarget**(2/3)
           * frac_vol_floc_initial(ConcAl, ConcClay, coag, material)**(-8/9)
           * (DiamTarget / material.Diameter)**((8*(DIM_FRACTAL-3)) / 9)
           )
517870
def median_fltr(dem, fsize=7, origmask=False):
    
    print("Applying median filter with size %s" % fsize)
    from scipy.ndimage.filters import median_filter
    dem_filt_med = median_filter(dem.filled(np.nan), fsize)
    #Now mask all nans
    out = np.ma.fix_invalid(dem_filt_med, copy=False, fill_value=dem.fill_value)
    if origmask:
        out = np.ma.array(out, mask=dem.mask, fill_value=dem.fill_value)
    out.set_fill_value(dem.fill_value)
    return out
805563
def video_pixel_noise_bottom(x, model_hparams, vocab_size):
  """"""
  input_noise = getattr(model_hparams, "video_modality_input_noise", 0.25)
  inputs = x
  if model_hparams.mode == tf.estimator.ModeKeys.TRAIN:
    background = tfp.stats.percentile(inputs, 50., axis=[0, 1, 2, 3])
    input_shape = common_layers.shape_list(inputs)
    input_size = tf.reduce_prod(input_shape[:-1])
    input_mask = tf.multinomial(
        tf.log([[input_noise, 1.-input_noise]]), input_size)
    input_mask = tf.reshape(tf.cast(input_mask, tf.int32),
                            input_shape[:-1]+[1])
    inputs = inputs * input_mask + background * (1 - input_mask)
  return video_bottom(inputs, model_hparams, vocab_size)