LAPV' algorithm (Pech2000)
326553
def varianceOfLaplacian(img):
    '''''''
    lap = cv2.Laplacian(img, ddepth=-1)#cv2.cv.CV_64F)
    stdev = cv2.meanStdDev(lap)[1]
    s = stdev[0]**2
    return s[0]
1320211
def _stieltjes_approx(dist, order, accuracy, normed, **kws):
    """"""
    kws["rule"] = kws.get("rule", "C")
    assert kws["rule"].upper() != "G"

    absisas, weights = chaospy.quad.generate_quadrature(
        accuracy, dist.range(), **kws)
    weights = weights*dist.pdf(absisas)

    poly = chaospy.poly.variable(len(dist))
    orth = [poly*0, poly**0]

    inner = numpy.sum(absisas*weights, -1)
    norms = [numpy.ones(len(dist)), numpy.ones(len(dist))]
    coeff1 = []
    coeff2 = []

    for _ in range(order):

        coeff1.append(inner/norms[-1])
        coeff2.append(norms[-1]/norms[-2])
        orth.append((poly-coeff1[-1])*orth[-1] - orth[-2]*coeff2[-1])

        raw_nodes = orth[-1](*absisas)**2*weights
        inner = numpy.sum(absisas*raw_nodes, -1)
        norms.append(numpy.sum(raw_nodes, -1))

        if normed:
            orth[-1] = orth[-1]/numpy.sqrt(norms[-1])

    coeff1.append(inner/norms[-1])
    coeff2.append(norms[-1]/norms[-2])
    coeff1 = numpy.transpose(coeff1)
    coeff2 = numpy.transpose(coeff2)
    norms = numpy.array(norms[1:]).T
    orth = orth[1:]
    return orth, norms, coeff1, coeff2
557750
def sph_harmonic_tp(nrows, ncols, n, m):
    
    nuvec = pysphi.ynunm(n, m, n + 1)
    num = nuvec[1:] * (-1) ** m
    num = num[::-1]
    yvec = np.concatenate([num, nuvec])
    
    theta = np.linspace(0.0, np.pi, nrows)
    phi = np.linspace(0.0, 2.0 * np.pi - 2.0 * np.pi / ncols)
    nu = np.array(range(-n, n + 1), dtype=np.complex128)
    
    out = np.zeros([nrows, ncols], dtype=np.complex128)

    for nn in xrange(0, nrows):
        exp_vec = np.exp(1j * nu * theta[nn])
        for mm in xrange(0, ncols):          
            vl = np.sum(yvec * exp_vec)
            out[nn, mm] = 1j ** m * np.exp(1j * m * phi[mm]) * vl

    return out
314745
def fit(self,Params0=None,grad_threshold=1e-2):
        

        # GPVD initialization
        lik = limix.CLikNormalNULL()
        # Initial Params
        if Params0==None:
            n_params = self.C1.getNumberParams()
            n_params+= self.C2.getNumberParams()
            Params0 = SP.rand(n_params)
        # MultiGP Framework
        covar = []
        gp    = []
        mean  = []
        for i in range(self.N):
            covar.append(limix.CLinCombCF())
            covar[i].addCovariance(self.C1)
            covar[i].addCovariance(self.C2)
            coeff = SP.array([self.eigen[i],1])
            covar[i].setCoeff(coeff)
            mean.append(limix.CLinearMean(self.Yt[:,i],SP.eye(self.P)))
        gpVD = limix.CGPvarDecomp(covar[0],lik,mean[0],SP.ones(self.N),self.P,self.Yt,Params0)
        for i in range(self.N):
            gp.append(limix.CGPbase(covar[i],lik,mean[i]))
            gp[i].setY(self.Yt[:,i])
            gpVD.addGP(gp[i])
                
        # Optimization
        gpVD.initGPs()
        gpopt = limix.CGPopt(gpVD)
        LML0=-1.0*gpVD.LML()
        start_time = time.time()
        conv = gpopt.opt()
        time_train = time.time() - start_time
        LML=-1.0*gpVD.LML()
        LMLgrad = SP.linalg.norm(gpVD.LMLgrad()['covar'])
        Params = gpVD.getParams()['covar']
    
        # Check whether limix::CVarianceDecomposition.train() has converged
        if conv!=True or LMLgrad>grad_threshold or Params.max()>10:
            print('limix::CVarianceDecomposition::train has not converged')
            res=None
        else:
            res = {
                'Params0':          Params0,
                'Params':           Params,
                'LML':              SP.array([LML]),
                'LML0':             SP.array([LML0]),
                'LMLgrad':          SP.array([LMLgrad]),
                'time_train':       SP.array([time_train]),
                }
        return res
        pass
159233
public static void go_upstream_a( int[] p, RandomIter flowRandomIter, RandomIter tcaRandomIter, RandomIter lRandomIter,
            int[] param ) {
        double area = 0, lenght = 0;
        int[] point = new int[2];
        int kk = 0, count = 0;

        point[0] = p[0];
        point[1] = p[1];
        // check how many pixels are draining in the considered pixel and select
        // the pixel with maximun tca
        for( int k = 1; k <= 8; k++ ) {
            if (flowRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0) == dirIn[k][2]) {
                // counts how many pixels are draining in the considere
                count++;
                if (tcaRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0) >= area) {
                    // if two pixels has the same tca select the pixel with the
                    // maximum vale of hacklength
                    if (tcaRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0) == area) {
                        if (lRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0) > lenght) {
                            kk = k;
                            area = tcaRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0);
                            lenght = lRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0);
                            point[0] = p[0] + dirIn[k][1];
                            point[1] = p[1] + dirIn[k][0];
                        }
                    } else {
                        kk = k;
                        area = tcaRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0);
                        lenght = lRandomIter.getSampleDouble(p[0] + dirIn[k][1], p[1] + dirIn[k][0], 0);
                        point[0] = p[0] + dirIn[k][1];
                        point[1] = p[1] + dirIn[k][0];
                    }
                }
            }

        }
        p[0] = point[0];
        p[1] = point[1];
        param[0] = kk;
        param[1] = count;
    }
552428
def cublasCgemv(handle, trans, m, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasCgemv_v2(handle,
                                       _CUBLAS_OP[trans], m, n,
                                       ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
892631
def delta_e_cie2000(lab_color_vector, lab_color_matrix, Kl=1, Kc=1, Kh=1):
    
    L, a, b = lab_color_vector

    avg_Lp = (L + lab_color_matrix[:, 0]) / 2.0

    C1 = numpy.sqrt(numpy.sum(numpy.power(lab_color_vector[1:], 2)))
    C2 = numpy.sqrt(numpy.sum(numpy.power(lab_color_matrix[:, 1:], 2), axis=1))

    avg_C1_C2 = (C1 + C2) / 2.0

    G = 0.5 * (1 - numpy.sqrt(numpy.power(avg_C1_C2, 7.0) / (numpy.power(avg_C1_C2, 7.0) + numpy.power(25.0, 7.0))))

    a1p = (1.0 + G) * a
    a2p = (1.0 + G) * lab_color_matrix[:, 1]

    C1p = numpy.sqrt(numpy.power(a1p, 2) + numpy.power(b, 2))
    C2p = numpy.sqrt(numpy.power(a2p, 2) + numpy.power(lab_color_matrix[:, 2], 2))

    avg_C1p_C2p = (C1p + C2p) / 2.0

    h1p = numpy.degrees(numpy.arctan2(b, a1p))
    h1p += (h1p < 0) * 360

    h2p = numpy.degrees(numpy.arctan2(lab_color_matrix[:, 2], a2p))
    h2p += (h2p < 0) * 360

    avg_Hp = (((numpy.fabs(h1p - h2p) > 180) * 360) + h1p + h2p) / 2.0

    T = 1 - 0.17 * numpy.cos(numpy.radians(avg_Hp - 30)) + \
        0.24 * numpy.cos(numpy.radians(2 * avg_Hp)) + \
        0.32 * numpy.cos(numpy.radians(3 * avg_Hp + 6)) - \
        0.2 * numpy.cos(numpy.radians(4 * avg_Hp - 63))

    diff_h2p_h1p = h2p - h1p
    delta_hp = diff_h2p_h1p + (numpy.fabs(diff_h2p_h1p) > 180) * 360
    delta_hp -= (h2p > h1p) * 720

    delta_Lp = lab_color_matrix[:, 0] - L
    delta_Cp = C2p - C1p
    delta_Hp = 2 * numpy.sqrt(C2p * C1p) * numpy.sin(numpy.radians(delta_hp) / 2.0)

    S_L = 1 + ((0.015 * numpy.power(avg_Lp - 50, 2)) / numpy.sqrt(20 + numpy.power(avg_Lp - 50, 2.0)))
    S_C = 1 + 0.045 * avg_C1p_C2p
    S_H = 1 + 0.015 * avg_C1p_C2p * T

    delta_ro = 30 * numpy.exp(-(numpy.power(((avg_Hp - 275) / 25), 2.0)))
    R_C = numpy.sqrt((numpy.power(avg_C1p_C2p, 7.0)) / (numpy.power(avg_C1p_C2p, 7.0) + numpy.power(25.0, 7.0)))
    R_T = -2 * R_C * numpy.sin(2 * numpy.radians(delta_ro))

    return numpy.sqrt(
        numpy.power(delta_Lp / (S_L * Kl), 2) +
        numpy.power(delta_Cp / (S_C * Kc), 2) +
        numpy.power(delta_Hp / (S_H * Kh), 2) +
        R_T * (delta_Cp / (S_C * Kc)) * (delta_Hp / (S_H * Kh)))
1131832
def TP0(dv, u):
    ''''''
    return np.linalg.norm(np.array(dv)) + np.linalg.norm(np.array(u))
764672
def gpp(V,E):
    
    model = Model("gpp")

    x = {}
    y = {}
    for i in V:
        x[i] = model.addVar(vtype="B", name="x(%s)"%i)
    for (i,j) in E:
        y[i,j] = model.addVar(vtype="B", name="y(%s,%s)"%(i,j))

    model.addCons(quicksum(x[i] for i in V) == len(V)/2, "Partition")

    for (i,j) in E:
        model.addCons(x[i] - x[j] <= y[i,j], "Edge(%s,%s)"%(i,j))
        model.addCons(x[j] - x[i] <= y[i,j], "Edge(%s,%s)"%(j,i))

    model.setObjective(quicksum(y[i,j] for (i,j) in E), "minimize")

    model.data = x
    return model
483708
def run_riemannian_relaxation(laplacian, initial_guess,
                              intrinsic_dim, relaxation_kwds):
    """"""
    n, s = initial_guess.shape
    relaxation_kwds = initialize_kwds(relaxation_kwds, n, s, intrinsic_dim)
    if relaxation_kwds['save_init']:
        directory = relaxation_kwds['backup_dir']
        np.save(os.path.join(directory, 'Y0.npy'),initial_guess)
        sp.io.mmwrite(os.path.join(directory, 'L_used.mtx'),
                      sp.sparse.csc_matrix(laplacian))

    lossf = relaxation_kwds['lossf']
    return RiemannianRelaxation.init(lossf, laplacian, initial_guess,
                                     intrinsic_dim, relaxation_kwds)
988274
def cublasCher2k(handle, uplo, trans, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasCher2k_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], 
                                        _CUBLAS_OP[trans], 
                                        n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,                 
                                                                               alpha.imag)),
                                        int(A), lda, int(B), ldb, 
                                        ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                         beta.imag)),
                                        int(C), ldc)
    cublasCheckStatus(status)
420853
def _laplace_fit(self,obj_type):
        

        # Get Mode and Inverse Hessian information
        y = self.fit(method='PML',printer=False)

        if y.ihessian is None:
            raise Exception("No Hessian information - Laplace approximation cannot be performed")
        else:

            self.latent_variables.estimation_method = 'Laplace'

            theta, Y, scores, states, states_var, X_names = self._categorize_model_output(self.latent_variables.get_z_values())

            # Change this in future
            try:
                latent_variables_store = self.latent_variables.copy()
            except:
                latent_variables_store = self.latent_variables

            return LaplaceResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method='Laplace',ihessian=y.ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
1100694
def cublasStpsv(handle, uplo, trans, diag, n, AP, x, incx):
    
    
    status = _libcublas.cublasStpsv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(AP), int(x), incx)
    cublasCheckStatus(status)
980787
def neighborhood_preserving_embedding(self, X, num_dims=None, reweight=True):
    ''''''
    if reweight:
      W = self.barycenter_edge_weights(X).matrix()
    else:
      W = self.matrix()
    # compute M = (I-W)'(I-W) as in LLE
    M = W.T.dot(W) - W.T - W
    if issparse(M):
      M = M.toarray()
    M.flat[::M.shape[0] + 1] += 1
    # solve generalized eig problem: X'MXa = \lambda X'Xa
    vals, vecs = eig(X.T.dot(M).dot(X), X.T.dot(X), overwrite_a=True,
                     overwrite_b=True)
    if num_dims is None:
      return vecs
    return vecs[:,:num_dims]
6969
def cublasZhpmv(handle, uplo, n, alpha, AP, x, incx, beta, y, incy):
    
    
    status = _libcublas.cublasZhpmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                            alpha.imag)),
                                       int(AP), int(x), incx,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real, 
                                                                         beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
120968
def _solve(self, x0, A, l, u, xmin, xmax):
        
        s = pips(self._costfcn, x0, A, l, u, xmin, xmax,
                 self._consfcn, self._hessfcn, self.opt)
        return s
401581
def compute_lunisolar
      dp, de = 0.0, 0.0

      begin
        l  = compute_l_iers2003(@jc)
        lp = compute_lp_mhb2000(@jc)
        f  = compute_f_iers2003(@jc)
        d  = compute_d_mhb2000(@jc)
        om = compute_om_iers2003(@jc)
        Const::NUT_LS.map do |x|
          x[0, 5] + x[5..-1].map { |y| y.to_s.sub(/\./, "").to_i }
        end.reverse.each do |x|
          arg = (x[0] * l + x[1] * lp + x[2] * f +
                 x[3] * d + x[4] * om) % Const::PI2
          sarg, carg = Math.sin(arg), Math.cos(arg)
          dp += (x[5] + x[6] * @jc) * sarg + x[ 7] * carg
          de += (x[8] + x[9] * @jc) * carg + x[10] * sarg
        end
        return [dp * Const::U2R, de * Const::U2R]
      rescue => e
        raise
      end
    end