Implements http://arxiv.org/abs/1401.0118, and uses the
    local reparameterization trick from http://arxiv.org/abs/1506.02557
324454
def black_box_variational_inference(logprob, D, num_samples):
    """"""

    def unpack_params(params):
        # Variational dist is a diagonal Gaussian.
        mean, log_std = params[:D], params[D:]
        return mean, log_std

    def gaussian_entropy(log_std):
        return 0.5 * D * (1.0 + np.log(2*np.pi)) + np.sum(log_std)

    rs = npr.RandomState(0)
    def variational_objective(params, t):
        """Provides a stochastic estimate of the variational lower bound."""
        mean, log_std = unpack_params(params)
        samples = rs.randn(num_samples, D) * np.exp(log_std) + mean
        lower_bound = gaussian_entropy(log_std) + np.mean(logprob(samples, t))
        return -lower_bound

    gradient = grad(variational_objective)

    return variational_objective, gradient, unpack_params
1244574
def locality_preserving_projections(self, coordinates, num_dims=None):
    ''''''
    X = np.atleast_2d(coordinates)  # n x d
    L = self.laplacian(normed=True)  # n x n
    u,s,_ = np.linalg.svd(X.T.dot(X))
    Fplus = np.linalg.pinv(u * np.sqrt(s))  # d x d
    n, d = X.shape
    if n >= d:  # optimized order: F(X'LX)F'
      T = Fplus.dot(X.T.dot(L.dot(X))).dot(Fplus.T)
    else:  # optimized order: (FX')L(XF')
      T = Fplus.dot(X.T).dot(L.dot(X.dot(Fplus.T)))
    L = 0.5*(T+T.T)
    return _null_space(L, num_vecs=num_dims, overwrite=True)
436715
def bosonic_constraints(a):
    
    substitutions = {}
    for i, ai in enumerate(a):
        substitutions[ai * Dagger(ai)] = 1.0 + Dagger(ai) * ai
        for aj in a[i+1:]:
            # substitutions[ai*Dagger(aj)] = -Dagger(ai)*aj
            substitutions[ai*Dagger(aj)] = Dagger(aj)*ai
            substitutions[Dagger(ai)*aj] = aj*Dagger(ai)
            substitutions[ai*aj] = aj*ai
            substitutions[Dagger(ai) * Dagger(aj)] = Dagger(aj) * Dagger(ai)

    return substitutions
913112
function MinMaxScalerTransforms(vector = [], nan_value = -1, return_nan=false) {
              const average = avg$2(vector);
              const standard_dev = sd(vector);
              const maximum = max$1(vector);
              const minimum = min$1(vector);
              const scale = (z) => {
                const scaledValue = (z - average) / (maximum - minimum);
                if (isNaN(scaledValue) && return_nan) return scaledValue;
                else if (isNaN(scaledValue) && return_nan === false) return (isNaN(standard_dev)) ? z : standard_dev;
                else return scaledValue;
              }; // equivalent to MinMaxScaler(z)
              const descale = (scaledZ) => {
                const descaledValue = (scaledZ * (maximum - minimum)) + average;
                if (isNaN(descaledValue) && return_nan) return descaledValue;
                else if (isNaN(descaledValue) && return_nan === false) return (isNaN(standard_dev)) ? scaledZ : standard_dev;
                else return descaledValue;
              };
              const values = vector.map(scale)
                .map(val => {
                  if (isNaN(val)) return nan_value;
                  else return val;
                });
              return {
                components: {
                  average,
                  standard_dev,
                  maximum,
                  minimum,
                },
                scale,
                descale,
                values,
              };
            }
1316092
def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
    '''
    '''
    if not numpy.all(numpy.array(deltas) > epsilon):
        raise ArgumentError('all deltas have to be greater than epsilon')

    bound = []
    for delta in deltas:
        # get boundary paths
        paths = pseudo.contour_paths(delta)

        # get vertices on boundary
        vertices = paths.vertices()

        # evaluate polynomial
        supremum = numpy.max(numpy.abs(p(vertices)))

        # compute bound
        bound.append(epsilon/(delta-epsilon)
                     * paths.length()/(2*numpy.pi*delta)
                     * supremum)
    return bound
507755
def advance(self):
        """"""
        if self.iter >= self.maxiter:
            raise ArgumentError('Maximum number of iterations reached.')
        if self.invariant:
            raise ArgumentError('Krylov subspace was found to be invariant '
                                'in the previous iteration.')

        N = self.V.shape[0]
        k = self.iter

        # the matrix-vector multiplication
        Av = self.A * self.V[:, [k]]

        if self.ortho == 'house':
            # Householder
            for j in range(k+1):
                Av[j:] = self.houses[j].apply(Av[j:])
                Av[j] *= numpy.conj(self.houses[j].alpha)
            if k+1 < N:
                house = House(Av[k+1:])
                self.houses.append(house)
                Av[k+1:] = house.apply(Av[k+1:]) * numpy.conj(house.alpha)
                self.H[:k+2, [k]] = Av[:k+2]
            else:
                self.H[:k+1, [k]] = Av[:k+1]
            # next line is safe due to the multiplications with alpha
            self.H[k+1, k] = numpy.abs(self.H[k+1, k])
            if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                    <= 1e-14:
                self.invariant = True
            else:
                vnew = numpy.zeros((N, 1), dtype=self.dtype)
                vnew[k+1] = 1
                for j in range(k+1, -1, -1):
                    vnew[j:] = self.houses[j].apply(vnew[j:])
                self.V[:, [k+1]] = vnew * self.houses[-1].alpha
        else:
            # determine vectors for orthogonalization
            start = 0

            # Lanczos?
            if self.ortho == 'lanczos':
                start = k
                if k > 0:
                    self.H[k-1, k] = self.H[k, k-1]
                    if self.M is not None \
                            and not isinstance(self.M, IdentityLinearOperator):
                        Av -= self.H[k, k-1] * self.P[:, [k-1]]
                    else:
                        Av -= self.H[k, k-1] * self.V[:, [k-1]]

            # (double) modified Gram-Schmidt
            for reortho in range(self.reorthos+1):
                # orthogonalize
                for j in range(start, k+1):
                    alpha = inner(self.V[:, [j]], Av, ip_B=self.ip_B)[0, 0]
                    if self.ortho == 'lanczos':
                        # check if alpha is real
                        if abs(alpha.imag) > 1e-10:
                            warnings.warn(
                                'Iter {0}: abs(alpha.imag) = {1} > 1e-10. '
                                'Is your operator self-adjoint in the '
                                'provided inner product?'
                                .format(self.iter, abs(alpha.imag)))
                        alpha = alpha.real
                    self.H[j, k] += alpha
                    if self.M is not None:
                        Av -= alpha * self.P[:, [j]]
                    else:
                        Av -= alpha * self.V[:, [j]]
            if self.M is not None:
                MAv = self.M * Av
                self.H[k+1, k] = norm(Av, MAv, ip_B=self.ip_B)
            else:
                self.H[k+1, k] = norm(Av, ip_B=self.ip_B)
            if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                    <= 1e-14:
                self.invariant = True
            else:
                if self.M is not None:
                    self.P[:, [k+1]] = Av / self.H[k+1, k]
                    self.V[:, [k+1]] = MAv / self.H[k+1, k]
                else:
                    self.V[:, [k+1]] = Av / self.H[k+1, k]

        # increase iteration counter
        self.iter += 1
287749
def despeckle_local(M, stds=2, width=2):
    

    N = np.array(M, dtype=np.float64)
    n, m = M.shape
    for i, j in itertools.product(range(width, n - width),
                                  range(width, m - width)):
        square = M[i - width:i + width, j - width:j + width]
        avg = np.average(square)
        std = np.std(square)
        if M[i, j] >= avg + stds * std:
            N[i, j] = avg
    return (N + N.T) / 2
812508
def copy(self):
        
        np = Particle()
        memmove(byref(np), byref(self), sizeof(self))
        return np
1102406
def local_constraints(self):
        
        cs = []
        for feat, val in self._avm.items():
            try:
                if val.supertypes and not val._avm:
                    cs.append((feat, val))
                else:
                    for subfeat, subval in val.features():
                        cs.append(('{}.{}'.format(feat, subfeat), subval))
            except AttributeError:
                cs.append((feat, val))
        return cs
640977
def refine(pm, node, optimizations):
    
    # Sanitize input
    pm.apply(ExpandGlobals, node)
    pm.apply(ExpandImportAll, node)
    pm.apply(NormalizeTuples, node)
    pm.apply(ExpandBuiltins, node)
    pm.apply(ExpandImports, node)
    pm.apply(NormalizeMethodCalls, node)
    pm.apply(NormalizeIsNone, node)
    pm.apply(SplitStaticExpression, node)
    pm.apply(NormalizeStaticIf, node)
    pm.apply(NormalizeTuples, node)
    pm.apply(NormalizeException, node)
    pm.apply(NormalizeMethodCalls, node)

    # Some early optimizations
    pm.apply(ComprehensionPatterns, node)

    pm.apply(RemoveLambdas, node)
    pm.apply(RemoveNestedFunctions, node)
    pm.apply(NormalizeCompare, node)

    pm.gather(ExtendedSyntaxCheck, node)

    pm.apply(ListCompToGenexp, node)
    pm.apply(RemoveComprehension, node)
    pm.apply(RemoveNamedArguments, node)

    # sanitize input
    pm.apply(NormalizeReturn, node)
    pm.apply(UnshadowParameters, node)
    pm.apply(FalsePolymorphism, node)

    # some extra optimizations
    apply_optimisation = True
    while apply_optimisation:
        apply_optimisation = False
        for optimization in optimizations:
            apply_optimisation |= pm.apply(optimization, node)[0]
556852
def _opbend_transform_mean(rs, fn_low, deriv=0):
    
    v = 0.0
    d = np.zeros((4,3), float)
    dd = np.zeros((4,3,4,3), float)
    #loop over the 3 cyclic permutations
    for p in np.array([[0,1,2], [2,0,1], [1,2,0]]):
        opbend = _opbend_transform([rs[p[0]], rs[p[1]], rs[p[2]], rs[3]], fn_low, deriv)
        v += opbend[0]/3
        index0 = np.where(p==0)[0][0] #index0 is the index of the 0th atom (rs[0])
        index1 = np.where(p==1)[0][0]
        index2 = np.where(p==2)[0][0]
        index3 = 3
        if deriv>0:
            d[0] += opbend[1][index0]/3
            d[1] += opbend[1][index1]/3
            d[2] += opbend[1][index2]/3
            d[3] += opbend[1][index3]/3
        if deriv>1:
            dd[0, :, 0, :] += opbend[2][index0, :, index0, :]/3
            dd[0, :, 1, :] += opbend[2][index0, :, index1, :]/3
            dd[0, :, 2, :] += opbend[2][index0, :, index2, :]/3
            dd[0, :, 3, :] += opbend[2][index0, :, index3, :]/3

            dd[1, :, 0, :] += opbend[2][index1, :, index0, :]/3
            dd[1, :, 1, :] += opbend[2][index1, :, index1, :]/3
            dd[1, :, 2, :] += opbend[2][index1, :, index2, :]/3
            dd[1, :, 3, :] += opbend[2][index1, :, index3, :]/3

            dd[2, :, 0, :] += opbend[2][index2, :, index0, :]/3
            dd[2, :, 1, :] += opbend[2][index2, :, index1, :]/3
            dd[2, :, 2, :] += opbend[2][index2, :, index2, :]/3
            dd[2, :, 3, :] += opbend[2][index2, :, index3, :]/3

            dd[3, :, 0, :] += opbend[2][index3, :, index0, :]/3
            dd[3, :, 1, :] += opbend[2][index3, :, index1, :]/3
            dd[3, :, 2, :] += opbend[2][index3, :, index2, :]/3
            dd[3, :, 3, :] += opbend[2][index3, :, index3, :]/3
    if deriv==0:
        return v,
    elif deriv==1:
        return v, d
    elif deriv==2:
        return v, d, dd
    else:
        raise ValueError("deriv must be 0, 1 or 2.")
723015
def convert_instancenorm(node, **kwargs):
    
    name, input_nodes, attrs = get_inputs(node, kwargs)

    eps = float(attrs.get("eps", 0.001))

    node = onnx.helper.make_node(
        'InstanceNormalization',
        inputs=input_nodes,
        outputs=[name],
        name=name,
        epsilon=eps)

    return [node]
669327
def draw_variable_local(self, size): 
        
        return ss.norm.rvs(loc=self.mu0, scale=self.sigma0, size=size)
1222852
def ucamel_method(func):
    
    frame_locals = inspect.currentframe().f_back.f_locals
    frame_locals[snake2ucamel(func.__name__)] = func
    return func
1218739
def update_variables(X, Z, U, prox_f, step_f, prox_g, step_g, L):
    
    if not hasattr(prox_g, '__iter__'):
        if prox_g is not None:
            dX = step_f/step_g * L.T.dot(L.dot(X) - Z + U)
            X[:] = prox_f(X - dX, step_f)
            LX, R, S = do_the_mm(X, step_f, Z, U, prox_g, step_g, L)
        else:
            # fall back to simple fixed-point method for f
            # see do_the_mm for normal definitions of LX,Z,R,S
            S = -X.copy()
            X[:] = prox_f(X, step_f)
            LX = X
            Z[:] = X[:]
            R = np.zeros(X.shape, dtype=X.dtype)
            S += X

    else:
        M = len(prox_g)
        dX = np.sum([step_f/step_g[i] * L[i].T.dot(L[i].dot(X) - Z[i] + U[i]) for i in range(M)], axis=0)
        X[:] = prox_f(X - dX, step_f)
        LX = [None] * M
        R = [None] * M
        S = [None] * M
        for i in range(M):
            LX[i], R[i], S[i] = do_the_mm(X, step_f, Z[i], U[i], prox_g[i], step_g[i], L[i])
    return LX, R, S
777191
def folding_model_energy(rvec, rcut):
    """
    r = np.linalg.norm(rvec) - rcut
    rr = r ** 2
    if r < 0.0:
        return -2.5 * rr
    return 0.5 * (r - 2.0) * rr
975434
def _update_parameters(self):
        
        self.repr_points, self.repr_points_log = self.sampler.get_samples(self.num_repr_points, self.proposal_function, self.burn_in_steps)

        if np.any(np.isnan(self.repr_points_log)) or np.any(np.isposinf(self.repr_points_log)):
            raise RuntimeError("Sampler generated representer points with invalid log values: {}".format(self.repr_points_log))

        # Removing representer points that have 0 probability of being the minimum (corresponding to log probability being minus infinity)
        idx_to_remove = np.where(np.isneginf(self.repr_points_log))[0]
        if len(idx_to_remove) > 0:
            idx = list(set(range(self.num_repr_points)) - set(idx_to_remove))
            self.repr_points = self.repr_points[idx, :]
            self.repr_points_log = self.repr_points_log[idx]

        # We predict with the noise as we need to make sure that var is indeed positive definite.
        mu, _ = self.model.predict(self.repr_points)
        # we need a vector
        mu = np.ndarray.flatten(mu)
        var = self.model.predict_covariance(self.repr_points)
        
        self.logP, self.dlogPdMu, self.dlogPdSigma, self.dlogPdMudMu = epmgp.joint_min(mu, var, with_derivatives=True)
        # add a second dimension to the array
        self.logP = np.reshape(self.logP, (self.logP.shape[0], 1))