Gaussianize X: [samples x phenotypes]
	- each phentoype is converted to ranks and transformed back to normal using the inverse CDF
505882
def rankStandardizeNormal(X):
	
	Is = X.argsort(axis=0)
	RV = SP.zeros_like(X)
	rank = SP.zeros_like(X)
	for i in range(X.shape[1]):
		x =  X[:,i]
		i_nan = SP.isnan(x)
		if 0:
			Is = x.argsort()
			rank = SP.zeros_like(x)
			rank[Is] = SP.arange(X.shape[0])
			#add one to ensure nothing = 0
			rank +=1
		else:
			rank = st.rankdata(x[~i_nan])
		#devide by (N+1) which yields uniform [0,1]
		rank /= ((~i_nan).sum()+1)
		#apply inverse gaussian cdf
		RV[~i_nan,i] = SP.sqrt(2) * special.erfinv(2*rank-1)
		RV[i_nan,i] = x[i_nan]
	return RV
404655
def uniform_binning_correction(x, n_bits=8):
  
  n_bins = 2**n_bits
  batch_size, height, width, n_channels = common_layers.shape_list(x)
  hwc = float(height * width * n_channels)

  x = x + tf.random_uniform(
      shape=(batch_size, height, width, n_channels),
      minval=0.0, maxval=1.0/n_bins)
  objective = -np.log(n_bins) * hwc * tf.ones(batch_size)
  return x, objective
314955
def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None, f_hat=None, W=None, Ki_W_i=None):
        
        Ki_f_init = np.zeros_like(Y)

        if K is None:
            K = kern.K(X)

        if f_hat is None:
            f_hat, _ = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)

        if W is None:
            W = -likelihood.d2logpdf_df2(f_hat, Y, Y_metadata=Y_metadata)

        if Ki_W_i is None:
            _, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave)

        logpdf_dfhat = likelihood.dlogpdf_df(f_hat, Y, Y_metadata=Y_metadata)

        if W.shape[1] == 1:
            W = np.diagflat(W)

        #Eq 14, and 16
        var_site = 1./np.diag(W)[:, None]
        mu_site = f_hat + var_site*logpdf_dfhat
        prec_site = 1./var_site
        #Eq 19
        marginal_cov = Ki_W_i
        marginal_mu = marginal_cov.dot(np.diagflat(prec_site)).dot(mu_site)
        marginal_var = np.diag(marginal_cov)[:, None]
        #Eq 30 with using site parameters instead of Gaussian site parameters
        #(var_site instead of sigma^{2} )
        posterior_cav_var = 1./(1./marginal_var - 1./var_site)
        posterior_cav_mean = posterior_cav_var*((1./marginal_var)*marginal_mu - (1./var_site)*Y)

        flat_y = Y.flatten()
        flat_mu = posterior_cav_mean.flatten()
        flat_var = posterior_cav_var.flatten()

        if Y_metadata is not None:
            #Need to zip individual elements of Y_metadata aswell
            Y_metadata_flat = {}
            if Y_metadata is not None:
                for key, val in Y_metadata.items():
                    Y_metadata_flat[key] = np.atleast_1d(val).reshape(-1, 1)

            zipped_values = []

            for i in range(Y.shape[0]):
                y_m = {}
                for key, val in Y_metadata_flat.items():
                    if np.isscalar(val) or val.shape[0] == 1:
                        y_m[key] = val
                    else:
                        #Won't broadcast yet
                        y_m[key] = val[i]
                zipped_values.append((flat_y[i], flat_mu[i], flat_var[i], y_m))
        else:
            #Otherwise just pass along None's
            zipped_values = zip(flat_y, flat_mu, flat_var, [None]*Y.shape[0])

        def integral_generator(yi, mi, vi, yi_m):
            def f(fi_star):
                #More stable in the log space
                p_fi = np.exp(likelihood.logpdf(fi_star, yi, yi_m)
                              - 0.5*np.log(2*np.pi*vi)
                              - 0.5*np.square(mi-fi_star)/vi)
                return p_fi
            return f

        #Eq 30
        p_ystar, _ = zip(*[quad(integral_generator(y, m, v, yi_m), -np.inf, np.inf)
                           for y, m, v, yi_m in zipped_values])
        p_ystar = np.array(p_ystar).reshape(-1, 1)
        return np.log(p_ystar)
564466
def random_rescale_to_mahalanobis(self, x):
        """"""
        x -= self.mean
        if any(x):
            x *= sum(self.randn(len(x))**2)**0.5 / self.mahalanobis_norm(x)
        x += self.mean
        return x
223877
def _initGP(self):
        
        if self._inference=='GP2KronSum': 
            signalPos = sp.where(sp.arange(self.n_randEffs)!=self.noisPos)[0][0]
            gp = GP2KronSum(Y=self.Y, F=self.sample_designs, A=self.trait_designs,
                            Cg=self.trait_covars[signalPos], Cn=self.trait_covars[self.noisPos],
                            R=self.sample_covars[signalPos])
        else:
            mean = MeanKronSum(self.Y, self.sample_designs, self.trait_designs) 
            Iok = vec(~sp.isnan(mean.Y))[:,0] 
            if Iok.all():   Iok = None
            covar = SumCov(*[KronCov(self.trait_covars[i], self.sample_covars[i], Iok=Iok) for i in range(self.n_randEffs)])
            gp = GP(covar = covar, mean = mean) 
        self.gp = gp
1088051
def gen_toy(f, nsample, bound, accuracy=10000, quiet=True, **kwd):
    
    # based on inverting cdf this is fast but you will need to give it a reasonable range
    # unlike roofit which is based on accept reject

    vnames = describe(f)
    if not quiet:
        print(vnames)
    my_arg = [kwd[v] for v in vnames[1:]]
    # random number
    # if accuracy is None: accuracy=10*numtoys
    r = npr.random_sample(nsample)
    x = np.linspace(bound[0], bound[1], accuracy)
    pdf = _vector_apply(f, x, tuple(my_arg))
    cdf = compute_cdf(pdf, x)
    if cdf[-1] < 0.01:
        warn(SmallIntegralWarning('Integral for given funcition is'
                                  ' really low. Did you give it a reasonable range?'))
    cdfnorm = cdf[-1]
    cdf /= cdfnorm

    # now convert that to toy
    ret = invert_cdf(r, cdf, x)

    if not quiet:
        # move this to plotting
        from matplotlib import pyplot as plt
        plt.figure()
        plt.title('comparison')
        numbin = 100
        h, e = np.histogram(ret, bins=numbin)
        mp = (e[1:] + e[:-1]) / 2.
        err = np.sqrt(h)
        plt.errorbar(mp, h, err, fmt='.b')
        bw = e[1] - e[0]
        y = pdf * len(ret) / cdfnorm * bw
        ylow = y + np.sqrt(y)
        yhigh = y - np.sqrt(y)
        plt.plot(x, y, label='pdf', color='r')
        plt.fill_between(x, yhigh, ylow, color='g', alpha=0.2)
        plt.grid(True)
        plt.xlim(bound)
        plt.ylim(ymin=0)
    return ret
640550
def filter_genes_cv_deprecated(X, Ecutoff, cvFilter):
    
    if issparse(X):
        raise ValueError('Not defined for sparse input. See `filter_genes_dispersion`.')
    mean_filter = np.mean(X, axis=0) > Ecutoff
    var_filter = np.std(X, axis=0) / (np.mean(X, axis=0) + .0001) > cvFilter
    gene_subset = np.nonzero(np.all([mean_filter, var_filter], axis=0))[0]
    return gene_subset
47765
def shrink(self):
        
        x = np.nan_to_num(self.X.values)

        # de-mean returns
        t, n = np.shape(x)
        meanx = x.mean(axis=0)
        x = x - np.tile(meanx, (t, 1))

        # compute sample covariance matrix
        sample = (1.0 / t) * np.dot(x.T, x)

        # compute prior
        var = np.diag(sample).reshape(-1, 1)
        sqrtvar = np.sqrt(var)
        _var = np.tile(var, (n,))
        _sqrtvar = np.tile(sqrtvar, (n,))
        r_bar = (np.sum(sample / (_sqrtvar * _sqrtvar.T)) - n) / (n * (n - 1))
        prior = r_bar * (_sqrtvar * _sqrtvar.T)
        prior[np.eye(n) == 1] = var.reshape(-1)

        # compute shrinkage parameters and constant
        if self.delta is None:

            # what we call pi-hat
            y = x ** 2.0
            phi_mat = np.dot(y.T, y) / t - 2 * np.dot(x.T, x) * sample / t + sample ** 2
            phi = np.sum(phi_mat)

            # what we call rho-hat
            term1 = np.dot((x ** 3).T, x) / t
            help_ = np.dot(x.T, x) / t
            help_diag = np.diag(help_)
            term2 = np.tile(help_diag, (n, 1)).T * sample
            term3 = help_ * _var
            term4 = _var * sample
            theta_mat = term1 - term2 - term3 + term4
            theta_mat[np.eye(n) == 1] = np.zeros(n)
            rho = sum(np.diag(phi_mat)) + r_bar * np.sum(
                np.dot((1.0 / sqrtvar), sqrtvar.T) * theta_mat
            )

            # what we call gamma-hat
            gamma = np.linalg.norm(sample - prior, "fro") ** 2

            # compute shrinkage constant
            kappa = (phi - rho) / gamma
            shrinkage = max(0.0, min(1.0, kappa / t))
            self.delta = shrinkage
        else:
            # use specified constant
            shrinkage = self.delta

        # compute the estimator
        sigma = shrinkage * prior + (1 - shrinkage) * sample
        return self.format_and_annualise(sigma)
1316080
def gauss_box_model(x, amplitude=1.0, mean=0.0, stddev=1.0, hpix=0.5):
    """"""
    z = (x - mean) / stddev
    z2 = z + hpix / stddev
    z1 = z - hpix / stddev
    return amplitude * (norm.cdf(z2) - norm.cdf(z1))
131998
def get_pca_ks_stats(self, maxrange=5):
        
        pc_to_phenotype_pairs = {}
        num_components = self.principal_observations_df.shape[1]
        if num_components < maxrange:
            maxrange = num_components

        phenotypes = self.principal_observations_df.phenotype.unique().tolist()
        for i in range(0, maxrange):
            phenotype_pair_to_ks = {}
            for p1, p2 in combinations(phenotypes, 2):
                p1_pc = self.principal_observations_df[self.principal_observations_df.phenotype == p1].iloc[:,i].as_matrix()
                p2_pc = self.principal_observations_df[self.principal_observations_df.phenotype == p2].iloc[:,i].as_matrix()
                phenotype_pair_to_ks[(p1, p2)] = ks_2samp(p1_pc, p2_pc)
            pc_to_phenotype_pairs[i + 1] = phenotype_pair_to_ks

        return pc_to_phenotype_pairs
236809
def sample_norm(x):
    
    original = x
    # square
    x = x ** 2
    # feature-map-wise sum
    x = torch.sum(x, dim=1)
    # scale by number of feature maps
    x *= 1.0 / original.shape[1]
    x += 10e-8
    x = torch.sqrt(x)
    return original / x.view(-1, 1, x.shape[-1])
634109
def _fused_batch_norm_op(self, input_batch, mean, variance, use_batch_stats):
    """"""
    # Store the original shape of the mean and variance.
    mean_shape = mean.get_shape()
    variance_shape = variance.get_shape()
    # The fused batch norm expects the mean, variance, gamma and beta
    # tensors to have dimension 1, so we flatten them to remove the
    # extra dimensions. In addition, it expects the input_batch to have
    # dimension 4, so we reshape it accordingly.
    gamma_flatten = tf.reshape(self._gamma, shape=(self._num_channels,))
    beta_flatten = tf.reshape(self._beta, shape=(self._num_channels,))
    flatten_mean = tf.reshape(mean, shape=(self._num_channels,))
    flatten_variance = tf.reshape(variance, shape=(self._num_channels,))
    use_batch_stats = tf.convert_to_tensor(use_batch_stats)

    input_shape = input_batch.get_shape()
    output_shape = [-1] + input_shape.as_list()[1:]

    flat_image_size = np.prod(self._image_shape, dtype=np.int32)

    if len(self._data_format) == 4:
      fusable_data_format = self._data_format
      fusable_batch = input_batch
    elif self._channel_index == 1 and self._image_shape:
      fusable_data_format = "NCHW"
      fusable_batch = tf.reshape(
          input_batch,
          shape=(-1, self._num_channels, 1, flat_image_size))
    else:
      # The CPU implementation of FusedBatchNorm only supports NHWC tensor
      # format for now.
      fusable_data_format = "NHWC"
      fusable_batch = tf.reshape(
          input_batch,
          shape=(-1, 1, flat_image_size, self._num_channels))

    common_args = {
        "scale": gamma_flatten,
        "offset": beta_flatten,
        "epsilon": self._eps,
        "data_format": fusable_data_format,
        "name": "batch_norm"
    }

    def use_batch_stats_fused_batch_norm():
      return tf.nn.fused_batch_norm(
          fusable_batch,
          mean=None,
          variance=None,
          is_training=True,
          **common_args)

    def moving_average_fused_batch_norm():
      return tf.nn.fused_batch_norm(
          fusable_batch,
          mean=flatten_mean,
          variance=flatten_variance,
          is_training=False,
          **common_args)

    batch_norm_op, mean, variance = utils.smart_cond(
        use_batch_stats, use_batch_stats_fused_batch_norm,
        moving_average_fused_batch_norm)

    if len(self._data_format) != 4:
      batch_norm_op = tf.reshape(batch_norm_op, output_shape)
    mean = tf.reshape(mean, mean_shape)
    variance = tf.reshape(variance, variance_shape)
    return batch_norm_op, mean, variance
948945
def inverse_categorical(x, mu=255):
    
    flat = x.reshape((-1, x.shape[-1]))
    indices = np.argmax(flat, axis=1).astype(np.float32)
    indices = (indices / mu) - 0.5
    inverted = inverse_mu_law(indices, mu=mu).reshape(x.shape[:-1])
    return ArrayWithUnits(inverted, x.dimensions[:2])
524300
def GroupNorm(x, group=32, gamma_initializer=tf.constant_initializer(1.)):
    
    shape = x.get_shape().as_list()
    ndims = len(shape)
    assert ndims == 4, shape
    chan = shape[1]
    assert chan % group == 0, chan
    group_size = chan // group

    orig_shape = tf.shape(x)
    h, w = orig_shape[2], orig_shape[3]

    x = tf.reshape(x, tf.stack([-1, group, group_size, h, w]))

    mean, var = tf.nn.moments(x, [2, 3, 4], keep_dims=True)

    new_shape = [1, group, group_size, 1, 1]

    beta = tf.get_variable('beta', [chan], initializer=tf.constant_initializer())
    beta = tf.reshape(beta, new_shape)

    gamma = tf.get_variable('gamma', [chan], initializer=gamma_initializer)
    gamma = tf.reshape(gamma, new_shape)

    out = tf.nn.batch_normalization(x, mean, var, beta, gamma, 1e-5, name='output')
    return tf.reshape(out, orig_shape, name='output')
509936
def create_model(samples_x, samples_y_aggregation,
                 n_restarts_optimizer=250, is_white_kernel=False):
    '''
    
    '''
    kernel = gp.kernels.ConstantKernel(constant_value=1,
                                       constant_value_bounds=(1e-12, 1e12)) * \
                                                gp.kernels.Matern(nu=1.5)
    if is_white_kernel is True:
        kernel += gp.kernels.WhiteKernel(noise_level=1, noise_level_bounds=(1e-12, 1e12))
    regressor = gp.GaussianProcessRegressor(kernel=kernel,
                                            n_restarts_optimizer=n_restarts_optimizer,
                                            normalize_y=True,
                                            alpha=1e-10)
    regressor.fit(numpy.array(samples_x), numpy.array(samples_y_aggregation))

    model = {}
    model['model'] = regressor
    model['kernel_prior'] = str(kernel)
    model['kernel_posterior'] = str(regressor.kernel_)
    model['model_loglikelihood'] = regressor.log_marginal_likelihood(regressor.kernel_.theta)

    return model
17299
def spm(x, model, y=None, n_samples=None, dx_min=-0.1,
        dx_max=0.1, n_dxs=5, dy_min=-0.1, dy_max=0.1, n_dys=5,
        angle_min=-30, angle_max=30, n_angles=31, black_border_size=0):
  
  if y is None:
    preds = model.get_probs(x)
    # Using model predictions as ground truth to avoid label leaking
    preds_max = reduce_max(preds, 1, keepdims=True)
    y = tf.to_float(tf.equal(preds, preds_max))
    y = tf.stop_gradient(y)
    del preds
  y = y / reduce_sum(y, 1, keepdims=True)

  # Define the range of transformations
  dxs = np.linspace(dx_min, dx_max, n_dxs)
  dys = np.linspace(dy_min, dy_max, n_dys)
  angles = np.linspace(angle_min, angle_max, n_angles)

  if n_samples is None:
    import itertools
    transforms = list(itertools.product(*[dxs, dys, angles]))
  else:
    sampled_dxs = np.random.choice(dxs, n_samples)
    sampled_dys = np.random.choice(dys, n_samples)
    sampled_angles = np.random.choice(angles, n_samples)
    transforms = zip(sampled_dxs, sampled_dys, sampled_angles)
  transformed_ims = parallel_apply_transformations(
      x, transforms, black_border_size)

  def _compute_xent(x):
    preds = model.get_logits(x)
    return tf.nn.softmax_cross_entropy_with_logits_v2(
        labels=y, logits=preds)

  all_xents = tf.map_fn(
      _compute_xent,
      transformed_ims,
      parallel_iterations=1)  # Must be 1 to avoid keras race conditions

  # Return the adv_x with worst accuracy

  # all_xents is n_total_samples x batch_size (SB)
  all_xents = tf.stack(all_xents)  # SB

  # We want the worst case sample, with the largest xent_loss
  worst_sample_idx = tf.argmax(all_xents, axis=0)  # B

  batch_size = tf.shape(x)[0]
  keys = tf.stack([
      tf.range(batch_size, dtype=tf.int32),
      tf.cast(worst_sample_idx, tf.int32)
  ], axis=1)
  transformed_ims_bshwc = tf.einsum('sbhwc->bshwc', transformed_ims)
  after_lookup = tf.gather_nd(transformed_ims_bshwc, keys)  # BHWC
  return after_lookup
546764
def cumulative_distribution(self, X):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        if (V == 0).all() or (U == 0).all():
            return np.zeros(V.shape[0])

        else:
            cdfs = [
                np.power(
                    np.power(U[i], -self.theta) + np.power(V[i], -self.theta) - 1,
                    -1.0 / self.theta
                )
                if (U[i] > 0 and V[i] > 0) else 0
                for i in range(len(U))
            ]

            return np.array([max(x, 0) for x in cdfs])