python类eigh()的实例源码

def _eig_decomposition(self, A, largest=True):
        n_features = A.shape[0]

        if (largest):
            eig_range = (n_features - self.n_components, n_features - 1)
        else:
            eig_range = (0, self.n_components - 1)

        eigvals, eigvecs = spla.eigh(A, eigvals=eig_range)

        if (largest):
            eigvals = eigvals[::-1]
            eigvecs = sp.fliplr(eigvecs)

        return eigvals, eigvecs

def calc_W(self, S_b, S_w):
        # W is the array of eigenvectors, where each column is a vector.
        eigenvalues, W = eigh(S_b, S_w)

        return W

def test_optimized_W(cls):
        tolerance = 1e-100

        vals, W = eigh(cls.model.S_b, cls.model.S_w)
        cls.assert_same(cls.model.W, W, tolerance=tolerance)

def test_calc_W(self):
        tolerance = 1e-100
        S_b = [[17.70840444, 17.96889098, 18.19513973],
               [17.96889098, 18.24564939, 18.46561872],
               [18.19513973, 18.46561872, 18.69940039]]
        S_w = [[0.94088804,  -0.05751511,  0.01467744],
               [-0.05751511,  1.01617648, -0.03831551],
               [0.01467744,  -0.03831551,  0.88440609]]

        W_model = self.model.calc_W(S_b, S_w)
        _, W_truth = eigh(S_b, S_w)

        self.assert_same(W_model, W_truth, tolerance=tolerance)

def _calculate_whitening_transformation_matrix(self, sample_from_prior):
        samples_vec = sp.asarray([self._dict_to_to_vect(x)
                                  for x in sample_from_prior])
        # samples_vec is an array of shape nr_samples x nr_features
        means = samples_vec.mean(axis=0)
        centered = samples_vec - means
        covariance = centered.T.dot(centered)
        w, v = la.eigh(covariance)
        self._whitening_transformation_matrix = (
            v.dot(sp.diag(1. / sp.sqrt(w))).dot(v.T))

def _logdet(A):
    """Compute the log det of a symmetric matrix."""
    vals = linalg.eigh(A)[0]
    # avoid negative (numerical errors) or zero (semi-definite matrix) values
    tol = vals.max() * vals.size * np.finfo(np.float64).eps
    vals = np.where(vals > tol, vals, tol)
    return np.sum(np.log(vals))

def set(self,gse):
        '''
        Set the Lambda matrix, Q matrix and QT matrix of the block.

        Parameters
        ----------
        gse : number
            The groundstate energy.
        '''
        sign=self.controllers['sign']
        if self.method=='S':
            lczs,Qs=self.controllers['lczs'],self.controllers['Qs']
            self.data['niters']=np.zeros(Qs.shape[0],dtype=np.int32)
            self.data['Lambdas']=np.zeros((Qs.shape[0],Qs.shape[2]),dtype=np.float64)
            self.data['Qs']=np.zeros(Qs.shape,dtype=Qs.dtype)
            self.data['QTs']=np.zeros((Qs.shape[0],Qs.shape[2]),dtype=Qs.dtype)
            for i,(lcz,Q) in enumerate(zip(lczs,Qs)):
                E,V=sl.eigh(lcz.T,eigvals_only=False)
                self.data['niters'][i]=lcz.niter
                self.data['Lambdas'][i,0:lcz.niter]=sign*(E-gse)
                self.data['Qs'][i,:,0:lcz.niter]=Q[:,0:lcz.niter].dot(V)
                self.data['QTs'][i,0:lcz.niter]=lcz.P[0,0]*V[0,:].conjugate()
        else:
            lanczos=self.controllers['lanczos']
            E,V=sl.eigh(lanczos.T,eigvals_only=False)
            self.data['Lambda']=sign*(E-gse)
            self.data['Q']=lanczos.P[:min(lanczos.nv0,lanczos.niter),:].T.conjugate().dot(V[:min(lanczos.nv0,lanczos.niter),:])
            self.data['QT']=HM.dagger(self.data['Q'])

def VCATEB(engine,app):
    '''
    This method calculates the topological Hamiltonian's spectrum.
    '''
    engine.rundependences(app.name)
    engine.gf(omega=app.mu)
    H=lambda kx,ky: -inv(engine.gf(k=[kx,ky]))
    result=np.zeros((app.path.rank('k'),engine.nopt+1))
    for i,paras in enumerate(app.path('+')):
        result[i,0]=i
        result[i,1:]=eigh(H(paras['k'][0],paras['k'][1]),eigvals_only=True)
    name='%s_%s'%(engine,app.name)
    if app.savedata: np.savetxt('%s/%s.dat'%(engine.dout,name),result)
    if app.plot: app.figure('L',result,'%s/%s'%(engine.dout,name))
    if app.returndata: return result

def berry_phase(H,path,ns):
    '''
    Calculate the Berry phase of some bands of a Hamiltonian along a certain path.

    Parameters
    ----------
    H : callable
        Input function which returns the Hamiltonian as a 2D array.
    path : iterable
        The path along which to calculate the Berry phase.
    ns : iterable of int
        The sequences of bands whose Berry phases are wanted.

    Returns
    -------
    1d ndarray
        The wanted Berry phase of the selected bands.
    '''
    ns=np.array(ns)
    for i,parameters in enumerate(path):
        new=eigh(H(**parameters))[1][:,ns]
        if i==0:
            result=np.ones(len(ns),new.dtype)
            evs=new
        else:
            for j in xrange(len(ns)):
                result[j]*=np.vdot(old[:,j],new[:,j])
        old=new
    else:
        for j in xrange(len(ns)):
            result[j]*=np.vdot(old[:,j],evs[:,j])
    return np.angle(result)/np.pi

def set(self,matrix):
        '''
        Set the basis.

        Parameters
        ----------
        matrix : callable
            The function to get the single particle matrix.
        '''
        Eup,Uup,Edw,Udw=[],[],[],[]
        for k in [()] if self.BZ is None else self.BZ.mesh('k'):
            m=matrix(k)
            es,us=sl.eigh(m[:m.shape[0]/2,:m.shape[0]/2])
            Edw.append(es)
            Udw.append(us)
            es,us=sl.eigh(m[m.shape[0]/2:,m.shape[0]/2:])
            Eup.append(es)
            Uup.append(us)
        Eup,Uup=np.asarray(Eup),np.asarray(Uup).transpose((1,0,2))
        Edw,Udw=np.asarray(Edw),np.asarray(Udw).transpose((1,0,2))
        if self.polarization=='up':
            self._E1_=Edw
            self._E2_=Eup
            self._U1_=Udw
            self._U2_=Uup
        else:
            self._E1_=Eup
            self._E2_=Edw
            self._U1_=Uup
            self._U2_=Udw

def __init__(self,ne=6,method='eigh',**karg):
        '''
        Constructor.

        Parameters
        ----------
        ne : integer, optional
            The number of energy spectrums.
        method : 'eigh','eigsh'
            The function used to calculate the spectrums. 'eigh' for `sl.eigh` and 'eigsh' for `HM.eigsh`.
        '''
        super(EB,self).__init__(**karg)
        self.ne=ne
        self.method=method

def eigs(self):
        '''
        The eigenvalues of the Lanczos matrix.

        Returns
        -------
        1d ndarray
            The eigenvalues of the Lanczos matrix.
        '''
        return sl.eigh(self.T,eigvals_only=True)

def test_lanczos():
    print 'test_lanczos'
    N,Nv,Ne,Niter=4000,15,1,750

    np.random.seed(1)
    m=np.random.random((N,N))+1j*np.random.random((N,N))
    m=m+m.T.conjugate()
    v0=np.random.random((Nv,N))
    exacteigs=sl.eigh(m,eigvals_only=True)[:Ne] if N<1000 else eigsh(m,k=Ne,return_eigenvectors=False,which='SA')[::-1]
    print 'Exact eigs:',exacteigs
    print

    a=Lanczos(m,deepcopy(v0),maxiter=Niter,keepstate=True)
    for i in xrange(Niter): a.iter()
    print 'diff from Hermitian: %s.'%sl.norm(a.T-a.T.T.conjugate())
    h=np.zeros((Niter,Niter),dtype=m.dtype)
    for i in xrange(Niter):
        hi=a.matrix.dot(a.vectors[i]).conjugate()
        for j in xrange(Niter):
            h[i,j]=hi.dot(a.vectors[j])
    print 'diff from h: %s.'%sl.norm(a.T-h)
    V=np.asarray(a.vectors)
    print 'vecotrs input diff: %s.'%sl.norm(V[0:Nv,:].T.dot(a.P)-v0.T)
    print 'vectors orthonormal diff: %s.'%sl.norm(np.identity(Niter)-V.dot(V.T.conjugate()))
    Leigs=a.eigs()[:Ne]
    print 'Lanczos eigs:',Leigs
    print 'eigenvalues diff: %s.'%sl.norm(exacteigs-Leigs)
    print

def plot_ellipse(splot, mean, cov, color):
    v, w = linalg.eigh(cov)
    u = w[0] / linalg.norm(w[0])
    angle = np.arctan(u[1] / u[0])
    angle = 180 * angle / np.pi  # convert to degrees
    # filled Gaussian at 2 standard deviation
    ell = mpl.patches.Ellipse(mean, 2 * v[0] ** 0.5, 2 * v[1] ** 0.5,
                              180 + angle, color=color)
    ell.set_clip_box(splot.bbox)
    ell.set_alpha(0.5)
    splot.add_artist(ell)
    splot.set_xticks(())
    splot.set_yticks(())

def _pre_compute(self, X, y):
        # even if X is very sparse, K is usually very dense
        K = safe_sparse_dot(X, X.T, dense_output=True)
        v, Q = linalg.eigh(K)
        QT_y = np.dot(Q.T, y)
        return v, Q, QT_y

def PCA(data, dims_rescaled_data=2):
    """
    returns: data transformed in 2 dims/columns + regenerated original data
    pass in: data as 2D NumPy array
    """
    import numpy as NP
    from scipy import linalg as LA
    m, n = data.shape
    # mean center the data
    data -= data.mean(axis=0)
    # calculate the covariance matrix
    R = NP.cov(data, rowvar=False)
    # calculate eigenvectors & eigenvalues of the covariance matrix
    # use 'eigh' rather than 'eig' since R is symmetric, 
    # the performance gain is substantial
    evals, evecs = LA.eigh(R)
    # sort eigenvalue in decreasing order
    idx = NP.argsort(evals)[::-1]
    evecs = evecs[:,idx]
    # sort eigenvectors according to same index
    evals = evals[idx]
    # select the first n eigenvectors (n is desired dimension
    # of rescaled data array, or dims_rescaled_data)
    evecs = evecs[:, :dims_rescaled_data]
    # carry out the transformation on the data using eigenvectors
    # and return the re-scaled data, eigenvalues, and eigenvectors
    return NP.dot(evecs.T, data.T).T, evals, evecs

def solve_eigen_problem(self, A, B, solver):
        N = A.testfunction[0].N
        s = A.testfunction[0].slice()
        self.V = np.zeros((N, N))
        self.lmbda = np.ones(N)
        if solver == 0:
            self.lmbda[s], self.V[s, s] = scipy_la.eigh(A.diags().toarray(),
                                                        B.diags().toarray())

        elif solver == 1:
            #self.lmbda[s], self.V[s, s] = scipy_la.eigh(B.diags().toarray())
            a = np.zeros((3, N-2))
            a[0, :] = B[0]
            a[2, :-2] = B[2]
            self.lmbda[s], self.V[s, s] = scipy_la.eig_banded(a, lower=True)

def fit(self, X, y):
        """Train xdawn spatial filters.

        Parameters
        ----------
        X : ndarray, shape (n_trials, n_channels, n_samples)
            ndarray of trials.
        y : ndarray shape (n_trials, 1)
            labels corresponding to each trial.

        Returns
        -------
        self : Xdawn instance
            The Xdawn instance.
        """
        Nt, Ne, Ns = X.shape

        self.classes_ = (numpy.unique(y) if self.classes is None else
                         self.classes)

        # FIXME : too many reshape operation
        tmp = X.transpose((1, 2, 0))
        Cx = numpy.matrix(self.estimator(tmp.reshape(Ne, Ns * Nt)))

        self.evokeds_ = []
        self.filters_ = []
        self.patterns_ = []
        for c in self.classes_:
            # Prototyped responce for each class
            P = numpy.mean(X[y == c, :, :], axis=0)

            # Covariance matrix of the prototyper response & signal
            C = numpy.matrix(self.estimator(P))

            # Spatial filters
            evals, evecs = eigh(C, Cx)
            evecs = evecs[:, numpy.argsort(evals)[::-1]]  # sort eigenvectors
            evecs /= numpy.apply_along_axis(numpy.linalg.norm, 0, evecs)
            V = evecs
            A = numpy.linalg.pinv(V.T)
            # create the reduced prototyped response
            self.filters_.append(V[:, 0:self.nfilter].T)
            self.patterns_.append(A[:, 0:self.nfilter].T)
            self.evokeds_.append(numpy.dot(V[:, 0:self.nfilter].T, P))

        self.evokeds_ = numpy.concatenate(self.evokeds_, axis=0)
        self.filters_ = numpy.concatenate(self.filters_, axis=0)
        self.patterns_ = numpy.concatenate(self.patterns_, axis=0)
        return self

def test_hamiltonian(self):
        # printing an example from a H2 file
        hfile = self._get_resource_path("H2Equilibrium.txt")
        hamiltonian = make_Hamiltonian(Hamiltonian_from_file(hfile))
        self.log.info(hamiltonian)
        # [[-0.24522469381221926 0 0 0.18093133934472627 ]
        # [0 -1.0636560168497590 0.18093133934472627 0]
        # [0 0.18093133934472627 -1.0636560168497592 0]
        # [0.18093133934472627 0 0 -1.8369675149908681]]
        self.assertSequenceEqual([str(i) for i in hamiltonian[0]],
                                 ['(-0.245224693812+0j)', '0j', '0j', '(0.180931339345+0j)'])
        self.assertSequenceEqual([str(i) for i in hamiltonian[1]],
                                 ['0j', '(-1.06365601685+0j)', '(0.180931339345+0j)', '0j'])
        self.assertSequenceEqual([str(i) for i in hamiltonian[2]],
                                 ['0j', '(0.180931339345+0j)', '(-1.06365601685+0j)', '0j'])
        self.assertSequenceEqual([str(i) for i in hamiltonian[3]],
                                 ['(0.180931339345+0j)', '0j', '0j', '(-1.83696751499+0j)'])

        # printing an example from a graph input
        n = 3
        v0 = np.zeros(n)
        v0[2] = 1
        v1 = np.zeros(n)
        v1[0] = 1
        v1[1] = 1
        v2 = np.zeros(n)
        v2[0] = 1
        v2[2] = 1
        v3 = np.zeros(n)
        v3[1] = 1
        v3[2] = 1

        pauli_list = [(1, Pauli(v0, np.zeros(n))), (1, Pauli(v1, np.zeros(n))),
                      (1, Pauli(v2, np.zeros(n))), (1, Pauli(v3, np.zeros(n)))]
        a = make_Hamiltonian(pauli_list)
        self.log.info(a)

        w, v = la.eigh(a, eigvals=(0, 0))
        self.log.info(w)
        self.log.info(v)

        data = {'000': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'001': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'010': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'011': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'100': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'101': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'110': 10}
        self.log.info(Energy_Estimate(data, pauli_list))
        data = {'111': 10}
        self.log.info(Energy_Estimate(data, pauli_list))

def nfiedler(A=None, tol=1e-12):
    """
    :param A: motif adjacency matrix A
    :param tol: no uesd in this program
    :return: the fiedler vector of the normalized laplacian of A
             (fiedler vector: eigenvector corresponding to the second smallest eigenvalue)

    """
    if A is None:
        print('Error! matrix A is None..')
        return None

    L = nlaplacian(A)
    # print('L\n', L, '\n')
    n = A.shape[0]
    # print(n)
    eigvalue, eigvector = SLA.eigh(L + np.eye(n))

    # print('eigvalue\n', eigvalue)
    # print('eigvector\n', eigvector[:, 0], '\n\n', eigvector[:, 1])
    # print(L + np.eye(n))
    # print('eigvalues\n', eigvalue)
    # print('eigvalue\n', eigvalue[:2])
    # print('eigvector\n', eigvector[0:1])
    # print(eigvalue)

    x = eigvector[:, 1]

    # print('\n\nx\n', x, '\n')
    # print(x)

    x = x / (np.sqrt(np.sum(A, 1)).T)

    # x = x.T
    # print('\nx\n', x.T)

    eigvalue = eigvalue[1] - 1

    # print('\neigvalue\n', eigvalue)
    # print('\nnp.argsort(e+igvector[1])\n', np.argsort(eigvector[1]), '\n')
    # print('eigvector[1]\n', eigvector[1])
    # print('\nx\n', x)
    # print('\neigvalue\n', eigvalue)

    return x, eigvalue