Python linalg.eigvals函数代码示例

def modularity_spectrum(G):
    """Return eigenvalues of the modularity matrix of G.

    Parameters
    ----------
    G : Graph
       A NetworkX Graph or DiGraph

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    See Also
    --------
    modularity_matrix

    References
    ----------
    .. [1] M. E. J. Newman, "Modularity and community structure in networks",
       Proc. Natl. Acad. Sci. USA, vol. 103, pp. 8577-8582, 2006.
    """
    from scipy.linalg import eigvals
    if G.is_directed():
        return eigvals(nx.directed_modularity_matrix(G))
    else:
        return eigvals(nx.modularity_matrix(G))

def get_graph_props(g):
    travgl = g.transitivity_avglocal_undirected()
    tru = g.transitivity_undirected()
    d = g.diameter()
    asd = g.assortativity_degree()
    apl = g.average_path_length()
    omega = g.omega()
    density = g.density()
    maxd = g.maxdegree()
    medd = np.median(g.degree())
    plaw = ig.power_law_fit(g.degree())
    spnorm = max(np.abs(spl.eigvals(g.get_adjacency().data)))
    leigs = spl.eigvals(g.laplacian())
    algc = abs(sorted([e for e in leigs if e > 1e-10])[1])
    
    return [travgl,     # avg. local transitivity
            tru,        # global transitivity
            d,          # diameter
            asd,        # degree assortativity
            apl,        # avg. path length
            omega,      # max. clique
            density,
            maxd,       # max. degree
            medd,       # median degree
            plaw.alpha, # power law exponent
            spnorm,     # largest eigenvalue of adj. matrix
            algc,       # 2nd smallest non-zero eigenvalue of laplacian
            ]

    def testMinrealBrute(self):
        for n, m, p in permutations(range(1,6), 3):
            s = matlab.rss(n, p, m)
            sr = s.minreal()
            if s.states > sr.states:
                self.nreductions += 1
            else:
                np.testing.assert_array_almost_equal(
                    np.sort(eigvals(s.A)), np.sort(eigvals(sr.A)))
                for i in range(m):
                    for j in range(p):
                        ht1 = matlab.tf(
                            matlab.ss(s.A, s.B[:,i], s.C[j,:], s.D[j,i]))
                        ht2 = matlab.tf(
                            matlab.ss(sr.A, sr.B[:,i], sr.C[j,:], sr.D[j,i]))
                        try:
                            self.assert_numden_almost_equal(
                                ht1.num[0][0], ht2.num[0][0],
                                ht1.den[0][0], ht2.den[0][0])
                        except Exception as e:
                            # for larger systems, the tf minreal's
                            # the original rss, but not the balanced one
                            if n < 6:
                                raise e

        self.assertEqual(self.nreductions, 2)

def format_eig_svd():
    def format_cplx(z):
        if z.imag < 1e-300:
            return '{0:.4f}'.format(z.real)
        return '{0:.4f}+{1:.4f}i'.format(z.real, z.imag)

    eig12 = sp.eigvals(generate_matrix(12))
    svd12 = sp.svdvals(generate_matrix(12))

    eig25 = sp.eigvals(generate_matrix(25))
    svd25 = sp.svdvals(generate_matrix(25))

    result12 = r'\begin{tabular}{cc}' + '\n'
    result12 += r'    Eigenvalues&Singular values\\' + '\n'
    result12 += '     \\hline\n'
    result25 = copy.copy(result12)
    for k in range(25):
        if k < 12:
            result12 += r'    ${0}$&${1:.4f}$\\'.format(format_cplx(eig12[k]), svd12[k]) + '\n'
        result25 += r'    ${0}$&${1:.4f}$\\'.format(format_cplx(eig25[k]), svd25[k]) + '\n'

    result12 += '\\end{tabular}\n'
    result25 += '\\end{tabular}\n'

    print(result12)

    print(result25)

    def test_dare(self):
        A = matrix([[-0.6, 0],[-0.1, -0.4]])
        Q = matrix([[2, 1],[1, 0]])
        B = matrix([[2, 1],[0, 1]])
        R = matrix([[1, 0],[0, 1]])

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T * X * A - X -
            A.T * X * B * solve(B.T * X * B + R, B.T * X * A) + Q, zeros((2,2)))
        assert_array_almost_equal(solve(B.T * X * B + R, B.T * X * A), G)
        # check for stable closed loop
        lam = eigvals(A - B * G)
        assert_array_less(abs(lam), 1.0)

        A = matrix([[1, 0],[-1, 1]])
        Q = matrix([[0, 1],[1, 1]])
        B = matrix([[1],[0]])
        R = 2

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T * X * A - X -
            A.T * X * B * solve(B.T *  X * B + R, B.T * X * A) + Q, zeros((2,2)))
        assert_array_almost_equal(B.T * X * A / (B.T * X * B + R), G)
        # check for stable closed loop
        lam = eigvals(A - B * G)
        assert_array_less(abs(lam), 1.0)

def set_aw(sys,poles):
    """Divide in controller in input and feedback part
       for anti-windup

    Usage
    =====
    [sys_in,sys_fbk]=set_aw(sys,poles)

    Inputs
    ------

    sys: controller
    poles : poles for the anti-windup filter

    Outputs
    -------
    sys_in, sys_fbk: controller in input and feedback part
    """
    sys=ss(sys);
    den_old=poly(eigvals(sys.A))
    den = poly(poles)
    tmp= tf(den_old,den,sys.Tsamp)
    tmpss=tf2ss(tmp)
    sys_in=ss(tmp*sys)
    sys_in.Tsamp=sys.Tsamp
    sys_fbk=ss(1-tmp)
    sys_fbk.Tsamp=sys.Tsamp
    return sys_in, sys_fbk

def bb_step(sys,X0=None,Tf=None,Ts=0.001):
    """Plot the step response of the continous system sys

    Call:
    y=bb_step(sys [,Tf=final time] [,Ts=time step])

    Parameters
    ----------
    sys : Continous System in State Space form
    X0: Initial state vector (not used yet)
    Ts  : sympling time
    Tf  : Final simulation time
 
    Returns
    -------
    Nothing

    """
    if Tf==None:
        vals = eigvals(sys.A)
        r = min(abs(real(vals)))
        if r < 1e-10:
            r = 0.1
        Tf = 7.0 / r
    sysd=c2d(sys,Ts)
    dstep(sysd,Tf=Tf)

    def testMinreal(self, verbose=False):
        """Test a minreal model reduction"""
        #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
        A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
        #B = [0.3, -1.3; 0.1, 0; 1, 0]
        B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
        #C = [0, 0.1, 0; -0.3, -0.2, 0]
        C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
        #D = [0 -0.8; -0.3 0]
        D = [[0., -0.8], [-0.3, 0.]]
        # sys = ss(A, B, C, D)

        sys = ss(A, B, C, D)
        sysr = minreal(sys)
        self.assertEqual(sysr.states, 2)
        self.assertEqual(sysr.inputs, sys.inputs)
        self.assertEqual(sysr.outputs, sys.outputs)
        np.testing.assert_array_almost_equal(
            eigvals(sysr.A), [-2.136154, -0.1638459])

        s = tf([1, 0], [1])
        h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
        hm = minreal(h)
        hr = (s+1)/(s**2+s+1)
        np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0])
        np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0])

def adjacency_spectrum(G, weight='weight'):
    """Return eigenvalues of the adjacency matrix of G.

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    adjacency_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.adjacency_matrix(G,weight=weight).todense())

def ar_model_check_stable(A):
    """check if this AR model is stable

    :Parameters:
        A : ndarray
            The coefficient matrix of the model
    """

    # inits and checks
    m, p = A.shape
    p /= m
    if p != round(p):
        raise ValueError('bad inputs!')

    # check for stable model
    A1 = N.concatenate((
        A,
        N.concatenate((
            N.eye((p - 1) * m),
            N.zeros(((p - 1) * m, m))
        ), axis=1)
    ))
    lambdas = NL.eigvals(A1)
    rval = True
    if (N.absolute(lambdas) > 1).any():
        rval = False
    del A1, lambdas
    return rval

def compare_eigen_methods():
    """ timing of different eigensolver methods"""
    import scipy.linalg as linalg 
    print '\n *** diagonalize a real symmetric band matrix by different methods ***\n'
    mt=mytimer.create(10)    
    try:
        N=float(sys.argv[1])
        S=float(sys.argv[2])
    except:
        sys.exit('supply N S (command line arguments): matrix dimension N and and number of super-diagonals S ')

    np.random.seed(1)
    a=np.random.random((N,N))

    ab=GMatrix(np.random.random((S+1,N)),storage='upper_banded')
    mt[0].start('lapack symmetric upper banded')
    print ab.store+'\n',ab.eigvals()[:5]
    mt[0].stop()

    a=ab.restore('full')
    mt[1].start('lapack symmetric full')
    print a.store+'\n',a.eigvals()[:5]
    mt[1].stop()

    print 'lingalg'
    mt[2].start('linalg general')
    print np.sort(linalg.eigvals(a.m).real)[:5]
    mt[2].stop()

    mytimer.table()

    def test_timescales(self):
        P_dense=self.bdc.transition_matrix()
        P=self.bdc.transition_matrix_sparse()
        ev=eigvals(P_dense)
        """Sort with decreasing magnitude"""
        ev=ev[np.argsort(np.abs(ev))[::-1]]
        ts=-1.0/np.log(np.abs(ev))

        """k=None"""
        with self.assertRaises(ValueError):
            tsn=timescales(P)

        """k is not None"""
        tsn=timescales(P, k=self.k)
        self.assertTrue(np.allclose(ts[1:self.k], tsn[1:]))

        """k is not None, ncv is not None"""
        tsn=timescales(P, k=self.k, ncv=self.ncv)
        self.assertTrue(np.allclose(ts[1:self.k], tsn[1:]))
        

        """tau=7"""      

        """k is not None"""
        tsn=timescales(P, k=self.k, tau=7)
        self.assertTrue(np.allclose(7*ts[1:self.k], tsn[1:]))

def projection_shifts_init(A, E, B, shift_options):
    """Find starting shift parameters for low-rank ADI iteration using
    Galerkin projection on spaces spanned by LR-ADI iterates.

    See [PK16]_, pp. 92-95.

    Parameters
    ----------
    A
        The |Operator| A from the corresponding Lyapunov equation.
    E
        The |Operator| E from the corresponding Lyapunov equation.
    B
        The |VectorArray| B from the corresponding Lyapunov equation.
    shift_options
        The shift options to use (see :func:`lyap_solver_options`).

    Returns
    -------
    shifts
        A |NumPy array| containing a set of stable shift parameters.
    """
    for i in range(shift_options['init_maxiter']):
        Q = gram_schmidt(B, atol=0, rtol=0)
        shifts = spla.eigvals(A.apply2(Q, Q), E.apply2(Q, Q))
        shifts = shifts[np.real(shifts) < 0]
        if shifts.size == 0:
            # use random subspace instead of span{B} (with same dimensions)
            if shift_options['init_seed'] is not None:
                np.random.seed(shift_options['init_seed'])
                np.random.seed(np.random.random() + i)
            B = B.space.make_array(np.random.rand(len(B), B.space.dim))
        else:
            return shifts
    raise RuntimeError('Could not generate initial shifts for low-rank ADI iteration.')

	def __init__(self,para):	
		self.parms = Parameters()
		self.parms.N = para.N
		self.parms.Nin = para.Nin
		self.parms.Nout= para.Nout
		self.parms.inp_sc= para.inp_sc
		self.parms.spectr_rad= para.spectr_rad
		self.parms.leak_rate= para.leak_rate
		self.parms.inpBias = para.inpBias

		self.parms.Pscaling = para.Pscaling
		self.parms.OutScaling = para.OutScaling

		self.parms.alpha = para.alpha

		self.parms.t1 = para.t1
		self.parms.dt = para.dt

		self.parms.y_d = para.y_d

		self.P = np.eye(self.parms.N)*self.parms.Pscaling#100
		self.inpBias = self.parms.inpBias * np.random.randn(self.parms.N,1)

		wtemp = np.random.randn(self.parms.N,self.parms.N)
		self.W = self.parms.spectr_rad*wtemp/max(abs(linalg.eigvals(wtemp)))
		self.V = np.random.randn(self.parms.N,self.parms.Nin)*np.dot(np.ones((self.parms.N,1)),self.parms.inp_sc) 

		self.Woutp = np.random.randn(self.parms.Nout,self.parms.N)

		self.state_A = np.random.rand(self.parms.N,1)*2.-1.
		self.state_B = self.state_A
		self.orgstate = self.state_A

def _run_hamiltonian(verbose=True):
    c = classicalHamiltonian()
    if verbose:
        print(c.potential(array([-0.5, 0.5])))
        print(c.potential(array([-0.5, 0.0])))
        print(c.potential(array([0.0, 0.0])))

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)

    hessian = nd.Hessian(c.potential)

    H = hessian(xopt)
    true_H = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    error_estimate = np.NAN
    if verbose:
        print(xopt)
        print('H', H)
        print('H-true_H', np.abs(H - true_H))
        # print('error_estimate', info.error_estimate)

        eigenvalues = linalg.eigvals(H)
        normal_modes = c.normal_modes(eigenvalues)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)
    return H, error_estimate, true_H