Find simplex structure in eigenvectors to begin PCCA+.


    Parameters
    ----------
    right_eigenvectors :  ndarray
        Right eigenvectors of transition matrix

    Returns
    -------
    index : ndarray
        Indices of simplex
495244
def index_search(right_eigenvectors):
    

    num_micro, num_eigen = right_eigenvectors.shape

    index = np.zeros(num_eigen, 'int')

    # first vertex: row with largest norm
    index[0] = np.argmax(
            [norm(right_eigenvectors[i]) for i in range(num_micro)])

    ortho_sys = right_eigenvectors - np.outer(np.ones(num_micro),
                                              right_eigenvectors[index[0]])

    for j in range(1, num_eigen):
        temp = ortho_sys[index[j - 1]].copy()
        for l in range(num_micro):
            ortho_sys[l] -= temp * dot(ortho_sys[l], temp)

        dist_list = np.array([norm(ortho_sys[l]) for l in range(num_micro)])

        index[j] = np.argmax(dist_list)

        ortho_sys /= dist_list.max()

    return index
910674
def cublasZgerc(handle, m, n, alpha, x, incx, y, incy, A, lda):
    

    status = _libcublas.cublasZgerc_v2(handle,
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                               alpha.imag)),
                                       int(x), incx, int(y), incy, int(A), lda)
    cublasCheckStatus(status)
347355
private static <T> IntPair findBridgeIndexes(PSArray<Point2D<T>> earlyHull, PSArray<Point2D<T>> laterHull, int earlyStart, int laterStart, MultipliableNumberSystem<T> ns) {
        int early = earlyStart;
        int later = laterStart;
        while (true) {
            int nextEarly = getPreIndex(early, earlyHull.size());
            int nextLater = getNextIndex(later, laterHull.size());
            if (LeftTurn.is(laterHull.get(later), earlyHull.get(early), earlyHull.get(nextEarly), ns))
                early = nextEarly;
            else if (RightTurn.is(earlyHull.get(early), laterHull.get(later), laterHull.get(nextLater), ns))
                later = nextLater;
            else
                break;
        }
        return new IntPair(early, later);
    }
827195
def laplacian_pca(self, coordinates, num_dims=None, beta=0.5):
    ''''''
    X = np.atleast_2d(coordinates)
    L = self.laplacian(normed=True)
    kernel = X.dot(X.T)
    kernel /= eigsh(kernel, k=1, which='LM', return_eigenvectors=False)
    L /= eigsh(L, k=1, which='LM', return_eigenvectors=False)
    W = (1-beta)*(np.identity(kernel.shape[0]) - kernel) + beta*L
    if num_dims is None:
      vals, vecs = np.linalg.eigh(W)
    else:
      vals, vecs = eigh(W, eigvals=(0, num_dims-1), overwrite_a=True)
    return X.T.dot(vecs).dot(vecs.T).T
947785
def cublasSsyr(handle, uplo, n, alpha, x, incx, A, lda): 
    
   
    status = _libcublas.cublasSsyr_v2(handle,
                                      _CUBLAS_FILL_MODE[uplo], n, 
                                      ctypes.byref(ctypes.c_float(alpha)),
                                      int(x), incx, int(A), lda)
    cublasCheckStatus(status)
648212
def _find_line_start_index(self, index):
        
        indexes = self._line_start_indexes

        pos = bisect.bisect_right(indexes, index) - 1
        return pos, indexes[pos]
736779
def get_nodal_planes(self):
        
        # Convert reference frame to NED
        self.tensor, self.tensor_sigma = self._to_ned()
        self.ref_frame = 'NED'
        # Eigenvalue decomposition
        # Tensor
        _, evect = utils.eigendecompose(self.tensor)
        # Rotation matrix
        _, rot_vec = utils.eigendecompose(np.matrix([[0., 0., -1],
                                                     [0., 0., 0.],
                                                     [-1., 0., 0.]]))
        rotation_matrix = (np.matrix(evect * rot_vec.T)).T
        if np.linalg.det(rotation_matrix) < 0.:
            rotation_matrix *= -1.
        flip_dc = np.matrix([[0., 0., -1.],
                             [0., -1., 0.],
                             [-1., 0., 0.]])
        rotation_matrices = sorted(
            [rotation_matrix, flip_dc * rotation_matrix],
            cmp=cmp_mat)
        nodal_planes = GCMTNodalPlanes()
        dip, strike, rake = [(180. / pi) * angle
                             for angle in utils.matrix_to_euler(rotation_matrices[0])]
        # 1st Nodal Plane
        nodal_planes.nodal_plane_1 = {'strike': strike % 360,
                                      'dip': dip,
                                      'rake': -rake}

        # 2nd Nodal Plane
        dip, strike, rake = [(180. / pi) * angle
                             for angle in utils.matrix_to_euler(rotation_matrices[1])]
        nodal_planes.nodal_plane_2 = {'strike': strike % 360.,
                                      'dip': dip,
                                      'rake': -rake}
        return nodal_planes
702845
public function getEigenvectors()
    {
        $vectors = $this->V;

        // Always return the eigenvectors of length 1.0
        $vectors = new Matrix($vectors);
        $vectors = array_map(function ($vect) {
            $sum = 0;
            for ($i = 0; $i < count($vect); ++$i) {
                $sum += $vect[$i] ** 2;
            }

            $sum = sqrt($sum);
            for ($i = 0; $i < count($vect); ++$i) {
                $vect[$i] /= $sum;
            }

            return $vect;
        }, $vectors->transpose()->toArray());

        return $vectors;
    }
1318192
def get_stationary_distribution(self):
        
        # The stationary distribution is proportional to the left-eigenvector
        # associated with the largest eigenvalue (i.e., 1) of the transition
        # matrix.
        check_is_fitted(self, "transmat_")
        eigvals, eigvecs = np.linalg.eig(self.transmat_.T)
        eigvec = np.real_if_close(eigvecs[:, np.argmax(eigvals)])
        return eigvec / eigvec.sum()
521515
def tocvx(B):
    
    Bcoo = B.tocoo()
    return spmatrix(Bcoo.data, Bcoo.row.tolist(), Bcoo.col.tolist())
727315
def _csr_matrix_indices(S):
    
    m, n = S.shape

    for i in range(m):
        for j in range(S.indptr[i], S.indptr[i+1]):
            row_index, col_index = i, S.indices[j]
            yield row_index, col_index
1143022
def _init_matrix(self, data, index, columns, dtype=None):
        
        data = prep_ndarray(data, copy=False)
        index, columns = self._prep_index(data, index, columns)
        data = {idx: data[:, i] for i, idx in enumerate(columns)}
        return self._init_dict(data, index, columns, dtype)
219890
public IVector getVectorFromColumn(int index) {
        IVector result = new IVector(rows);
        for (int i = 0; i < rows; i++) {
            result.realvector[i] = realmatrix[i][index];
            result.imagvector[i] = imagmatrix[i][index];
        }
        return result;
    }
1269821
def S(self):
        
        if self.c is None:
            raise RuntimeError('Cross-spectral density requires noise '
                               'covariance matrix c.')
        H = self.H()
        # TODO: can we do that more efficiently?
        S = np.empty(H.shape, dtype=H.dtype)
        for f in range(H.shape[2]):
            S[:, :, f] = H[:, :, f].dot(self.c).dot(H[:, :, f].conj().T)
        return S
399161
def get_vectors(self, indices=None):
        ''''''
        H_ = self._deflated_solver.H
        (n_, n) = H_.shape
        coeffs = self.coeffs if indices is None else self.coeffs[:, indices]
        return numpy.c_[self._deflated_solver.V[:, :n],
                        self._deflated_solver.projection.U].dot(coeffs)
183600
def preprocess_with_pca(adata, n_pcs=None, random_state=0):
    
    if n_pcs == 0:
        logg.info('    using data matrix X directly (no PCA)')
        return adata.X
    elif n_pcs is None and 'X_pca' in adata.obsm_keys():
        logg.info('    using \'X_pca\' with n_pcs = {}'
                  .format(adata.obsm['X_pca'].shape[1]))
        return adata.obsm['X_pca']
    elif ('X_pca' in adata.obsm_keys()
          and adata.obsm['X_pca'].shape[1] >= n_pcs):
        logg.info('    using \'X_pca\' with n_pcs = {}'
                  .format(n_pcs))
        return adata.obsm['X_pca'][:, :n_pcs]
    else:
        n_pcs = N_PCS if n_pcs is None else n_pcs
        if adata.X.shape[1] > n_pcs:
            logg.info('    computing \'X_pca\' with n_pcs = {}'.format(n_pcs))
            logg.hint('avoid this by setting n_pcs = 0')
            X = pca(adata.X, n_comps=n_pcs, random_state=random_state)
            adata.obsm['X_pca'] = X
            return X
        else:
            logg.info('    using data matrix X directly (no PCA)')
            return adata.X
129561
def eigen_table(self):
        """"""
        idx = ["Eigenvalue", "Variability (%)", "Cumulative (%)"]
        table = pd.DataFrame(
            np.array(
                [self.eigenvalues, self.inertia, self.cumulative_inertia]
            ),
            columns=["F%s" % i for i in range(1, self.keep + 1)],
            index=idx,
        )

        return table