def test_diag(self):
     a = (100*get_mat(5)).astype('f')
     b = a.copy()
     for k in range(5):
         for l in range(max((k-1,0)),5):
             b[l,k] = 0
     assert_equal(triu(a,k=2),b)
     b = a.copy()
     for k in range(5):
         for l in range(k+3,5):
             b[l,k] = 0
     assert_equal(triu(a,k=-2),b)

def make_edges(n):
    A = la.triu(np.random.randint(1,50,(n,n))*(np.random.rand(n,n)>.5))
    S = []
    for index, x in np.ndenumerate(A):
        if x != 0:
            S.append((str(index[0]), str(index[1]), x))
    return S

 def test_basic(self):
     a = (100*get_mat(5)).astype('l')
     b = a.copy()
     for k in range(5):
         for l in range(k+1,5):
             b[l,k] = 0
     assert_equal(triu(a),b)

def Problem1(n):
    """Use linalg.toeplitz() and linalg.triu() to generate matrices of arbitrary size"""

    from scipy.linalg import triu


    ut = triu([[0]*i+[x for x in range(1,(n+1)-i)] for i in range(n)])
    toep = la.toeplitz([1.0/(i+1) for i in range(n)])

    return ut, toep

 def yty2():
     invT = S
     #log("old invT = " + str(invT))
     if j == p_eff - 1:
         invT[:p_eff, :p_eff] = triu(X2[:p_eff, :m].dot(np.conj(X2)[:p_eff, :m].T))
         #log("invT = " + str(invT))
         for jj in range(p_eff):
             invT[jj,jj] = (invT[jj,jj] - 1.)/2.
     #log("invT = {}".format(invT))
     return invT

def spd(A):
    """
    spd(A) -> True if A is a symmetric positive definite matrix

    constructs an upper triangular matrix from A and then checks that each
    diagonal element is greater than zero (this comes from the test that a
    matrix is SPD iff all principal leading minors are nonzero, and that these
    determinants will carry through to the diagonalization.

    alternative test is testing eigs(A) > 0 but this is much faster as it
    doesn't involve root-finding
    """
    # test symmetric first, as it's easier
    if not symmetric(A):
        return False
    

    t = linalg.triu(A)
    return (np.diag(t) > 0).all()

def kruskal(A):
    size=A.shape
    minSpanTree=sp.zeros(size)
    nodesTree=sp.arange(size[0])
    D=sp.ones(size)*sp.arange(size[0])
    C=la.triu(A)
    Q=sp.concatenate([[C.flatten()],[D.T.flatten()],[D.flatten()]])
    W=Q[:,sp.nonzero(C.flatten())[0]]
    edges=W[:,W[0,:].argsort()]
    i=0
    j=0
    #while np.sum(nodesTree!=nodesTree[0])>0:
    while (j<(size[0]-1)):
        now=edges[:,i]
        i=i+1
        if nodesTree[now[1]]!=nodesTree[now[2]]:
            minSpanTree[now[1],now[2]]=now[0]
            nodesTree[nodesTree==nodesTree[now[2]]]=nodesTree[now[1]]
            j=j+1
    return minSpanTree+minSpanTree.T

def GaussSeidel(A, b, tolerance = 1.e-10, MaxSteps = 100):
    """Solve the linear system A x = b using the Gauss-Seidel method, 
    starting from the trivial initial guess."""
    
    x = np.zeros_like(b)
    
    Anorm = A.copy()
    bnorm = b.copy()
    n = len(b)
    
    for i in range(n):
        bnorm[i] /= A[i, i]
        Anorm[i, :] /= A[i, i]
    
    # Compute the split
    D = np.eye(n)
    AL = la.tril(D - Anorm)
    AU = la.triu(D - Anorm)
    N = np.eye(n) - AL
    P = AU
    
    # Compute the convergence matrix and check its spectral radius
    M = np.dot(la.inv(N), P)
    eigenvalues, eigenvectors = la.eig(M)
    rho = np.amax(np.absolute(eigenvalues))
    if (rho > 1):
        print("Gauss-Seidel will not converge as the"\
              " largest eigenvalue of the convergence matrix is {}".format(rho))
    
    for j in range(MaxSteps):
        x_old = x.copy()
        for i in range(n):
            x[i] = bnorm[i] + np.dot(AL[i, :], x) + np.dot(AU[i, :], x_old)
        if (la.norm(x - x_old) < tolerance):
            print("Gauss-Seidel converged in {} iterations.".format(j))
            break
    
    return x

def Jacobi(A, b, tolerance = 1.e-10, MaxSteps = 100):
    """Solve the linear system A x = b using Jacobi's method, 
    starting from the trivial initial guess."""
    
    x = np.zeros_like(b)
    
    Anorm = A.copy()
    bnorm = b.copy()
    n = len(b)
    
    for i in range(n):
        bnorm[i] /= A[i, i]
        Anorm[i, :] /= A[i, i]
    
    # Compute the split
    N = np.eye(n)
    P = N - Anorm
    AL = la.tril(P)
    AU = la.triu(P)
    
    # Compute the convergence matrix and check its spectral radius
    M = np.dot(la.inv(N), P)
    eigenvalues, eigenvectors = la.eig(M)
    rho = np.amax(np.absolute(eigenvalues))
    if (rho > 1):
        print("Jacobi will not converge as the"\
            " largest eigenvalue of the convergence matrix is {}".format(rho))
    
    for j in range(MaxSteps):
        x_old = x.copy()
        x = bnorm + np.dot(AL + AU, x)
        if (la.norm(x - x_old) < tolerance):
            print "Jacobi converged in ", j, " iterations."
            break
    
    return x

def pseudoSpect(A, npts=200, s=2., gridPointSelect=100, verbose=True,
                lstSqSolve=True):
    """ 
    original code from http://www.cs.ox.ac.uk/projects/pseudospectra/psa.m
    % psa.m - Simple code for 2-norm pseudospectra of given matrix A.
    %         Typically about N/4 times faster than the obvious SVD method.
    %         Comes with no guarantees!   - L. N. Trefethen, March 1999.
    
    parameter: A: the matrix to analyze
               npts: number of points at the grid
               s: axis limits (-s ... +s)
               gridPointSelect: ???
               verbose: prints progress messages
               lstSqSolve: if true, use least squares in algorithm where
                  solve could be used (probably) instead. (replacement for
                  ldivide in MatLab)
    """
    
    from scipy.linalg import schur, triu
    from pylab import (meshgrid, norm, dot, zeros, eye, diag, find,  linspace,                       
                       arange, isreal, inf, ones, lstsq, solve, sqrt, randn,
                       eig, all)

    ldiv = lambda M1,M2 :lstsq(M1,M2)[0] if lstSqSolve else lambda M1,M2: solve(M1,M2)

    def planerot(x):
        '''
        return (G,y)
        with a matrix G such that y = G*x with y[1] = 0    
        '''
        G = zeros((2,2))
        xn = x / norm(x)
        G[0,0] = xn[0]
        G[1,0] = -xn[1]
        G[0,1] = xn[1]
        G[1,1] = xn[0]
        return G, dot(G,x)

    xmin = -s
    xmax = s
    ymin = -s
    ymax = s;  
    x = linspace(xmin,xmax,npts,endpoint=False)
    y = linspace(ymin,ymax,npts,endpoint=False)
    xx,yy = meshgrid(x,y)
    zz = xx + 1j*yy
     
    #% Compute Schur form and plot eigenvalues:
    T,Z = schur(A,output='complex');
        
    T = triu(T)
    eigA = diag(T)
    
    # Reorder Schur decomposition and compress to interesting subspace:
    select = find( eigA.real > -250)           # % <- ALTER SUBSPACE SELECTION
    n = len(select)
    for i in arange(n):
        for k in arange(select[i]-1,i,-1): #:-1:i
            G = planerot([T[k,k+1],T[k,k]-T[k+1,k+1]] )[0].T[::-1,::-1]
            J = slice(k,k+2)
            T[:,J] = dot(T[:,J],G)
            T[J,:] = dot(G.T,T[J,:])
          
    T = triu(T[:n,:n])
    I = eye(n);
    
    # Compute resolvent norms by inverse Lanczos iteration and plot contours:
    sigmin = inf*ones((len(y),len(x)));
    #A = eye(5)
    niter = 0
    for i in arange(len(y)): # 1:length(y)        
        if all(isreal(A)) and (ymax == -ymin) and (i > len(y)/2):
            sigmin[i,:] = sigmin[len(y) - i,:]
        else:
            for jj in arange(len(x)):
                z = zz[i,jj]
                T1 = z * I - T 
                T2 = T1.conj().T
                if z.real < gridPointSelect:    # <- ALTER GRID POINT SELECTION
                    sigold = 0
                    qold = zeros((n,1))
                    beta = 0
                    H = zeros((100,100))                
                    q = randn(n,1) + 1j*randn(n,1)                
                    while norm(q) < 1e-8:
                        q = randn(n,1) + 1j*randn(n,1)                
                    q = q/norm(q)
                    for k in arange(99):
                        v = ldiv(T1,(ldiv(T2,q))) - dot(beta,qold)
                        #stop
                        alpha = dot(q.conj().T, v).real
                        v = v - alpha*q
                        beta = norm(v)
                        qold = q
                        q = v/beta
                        H[k+1,k] = beta
                        H[k,k+1] = beta
                        H[k,k] = alpha
                        if (alpha > 1e100):
                            sig = alpha 
#.........这里部分代码省略.........

def get_stat_from_dynamics(singlecdt, tmin=None, tmax=None):
    """Computes stationary autocorrelation vector from autocorr matrix.

    This function uses the autocorrelation matrix, that stores two-point
    autocorrelation functions between the various acquisition time-points,
    to compute the autocorrelation function in case of stationary hypothesis.
    Computation is fast but result is mostly unreliable (depends very much
    on the accuracy of the dynamical autocorrelation estimates, which is
    usually quite low). Use set_stationary_autocorrelation instead.

    Parameters
    ----------
    singlecdt : UnivariateConditioned instance
    tmin : float (default None)
    tmax : float (default None)

    Returns
    -------
    dts, cts, res
    dts : array of floats
        time intervals
    cts : array of ints
        sample counts
    res : array of floats
        autocorrelation values

    Note
    ----
    The estimate of the autocorrelation function using this procedure gives
    very poor accuracy estimates, and should be used only for quick inspection
    when a Univariate has been created and computed.
    For a better autocorrelation function estimate, it is necessary to parse
    samples another time, using only the sample average estimated in
    Univariate conditioned instances.
    """
    times = singlecdt.time
    autocorr = singlecdt.autocorr
    cts = singlecdt.count_two
    # Resize matrices depending on time limits
    indexlow, indexup = 0, None
    if tmin is not None:
        while indexlow < len(times) and times[indexlow] < tmin:
            indexlow += 1
    if tmax is not None:
        indexup = indexlow
        while indexup < len(times) and times[indexup] < tmax:
            indexup += 1
    sl = slice(indexlow, indexup)
    times = times[sl]
    autocorr = autocorr[sl, sl]
    cts = cts[sl, sl]
    # how many time-points
    nframes = len(times)
    all_counts = np.zeros(nframes, dtype=np.int)
    res = np.zeros(nframes, dtype=np.float)
    dts = np.zeros(nframes, dtype=np.float)
    col = np.zeros(nframes, dtype=np.int16)
    col[-1] = 1  # initialisation
    for k in range(nframes):
        col[k] = 1
        col[k-1] -= 1
        forward = triu(toeplitz(col))
        all_counts[k] = np.sum(forward * cts)
        res[k] = np.sum(forward * cts * autocorr)/all_counts[k]
        dts[k] = times[k] - times[0]
    return dts, all_counts, res

def toepOne(n):
    return la.triu(la.toeplitz(sp.arange(1,n+1,1)))