Python scipy.mat函数代码示例

    def __init__(self,m,n,l1_weight,eval_func, eval_func_for_test_set = None, output_func = None, tol = None):
        self.M = m
        self.N = n
        
        d = np.ones((n))
       
        self.diagH = sparse.diags(d,0)
        # to reference the element, use diag.data[0,0]

        self.iter = 0
        
        self.l1_weight = l1_weight
        self.eval_func = eval_func
        
        self.eval_func_for_test_set  = eval_func_for_test_set
        self.output_func = output_func

        self.x = scipy.mat([0]*n).T
        self.grad = scipy.mat([0]*n).T
        self.dir = scipy.mat([0]*n).T
        self.loss = 0
        
        if tol == None:
            self.tol = 1e-4
        else:
            self.tol = tol

def calcInvFisher(sigma, invSigma=None, factorSigma=None):
    """ Efficiently compute the exact inverse of the FIM of a Gaussian.
    Returns a list of the diagonal blocks. """
    if invSigma == None:
        invSigma = inv(sigma)
    if factorSigma == None:
        factorSigma = cholesky(sigma)
    dim = sigma.shape[0]

    invF = [mat(1 / (invSigma[-1, -1] + factorSigma[-1, -1] ** -2))]
    invD = 1 / invSigma[-1, -1]
    for k in reversed(list(range(dim - 1))):
        v = invSigma[k + 1:, k]
        w = invSigma[k, k]
        wr = w + factorSigma[k, k] ** -2
        u = dot(invD, v)
        s = dot(v, u)
        q = 1 / (w - s)
        qr = 1 / (wr - s)
        t = -(1 + q * s) / w
        tr = -(1 + qr * s) / wr
        invF.append(blockCombine([[qr, tr * u], [mat(tr * u).T, invD + qr * outer(u, u)]]))
        invD = blockCombine([[q , t * u], [mat(t * u).T, invD + q * outer(u, u)]])

    invF.append(sigma)
    invF.reverse()
    return invF

def steepest_descent(A, b, x0, tol=1e-8):
    """
    Uses the steepest descent method to find the x that satisfies Ax = b.

    Inputs:
        A: An m x n NumPy array
        b: An m x 1 NumPy array
        x0: An n x 1 NumPy array that represents the initial guess at a
            solution.
        tol (optional): The tolerance level for convergence. This is compared
                        against the norm(x_n+1 - x_n) each iteration.

    Outputs:
        x: The x that satisfies the equation.
    """
    A = sp.mat(A)
    b = sp.reshape(sp.mat(b),(b.size,1))


    def grad(A, b, x):
        """
        Find the gradient of ||Ax - b||
        Inputs:
            A: An m x n NumPy matrix.
            b: An m x 1 NumPy matrix.
            x: An n x a NumPy matrix.

        Outputs:
            grad: A NumPy matrix representing the gradient of ||Ax - b||
        """
        return np.mat(2  * A.T*(A*x - b))

    def solve_alpha_k(A, b, x):
        """
        Solves for alpha in the steepest descent algorithm
        x_n+1 = x_n - alpha * grad(x_n)

        Inputs:
            A: An m x n NumPy array
            b: An m x 1 NumPy array
            x: The x value where you want alpha to be defined for.

        Outputs:
            alpha: The alpha satisfying the algorithm above.
        """

        gradient = grad(A, b, x)
        return np.array(
            (gradient.T * gradient)/(2 * gradient.T * A.T * A * gradient))[0]



    xold = sp.reshape(sp.mat(x0),(x0.size,1))
    xnew = xold - grad(A, b, xold) * solve_alpha_k(A,b,xold)

    while la.norm(xold - xnew) > tol:
        xold = xnew
        xnew = xold - grad(A, b, xold) * solve_alpha_k(A,b,xold)

    return xnew

def dsimul(sys,u):
    """Simulate the discrete system sys
    Only for discrete systems!!!

    Call:
    y=dsimul(sys,u)

    Parameters
    ----------
    sys : Discrete System in State Space form
    u   : input vector
    Returns
    -------
    y: ndarray
    Simulation results

    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ns=shape(u)[1]
    xk=zeros((nx,1))
    for i in arange(0,ns):
        uk=u[:,i]
        xk_1=a*xk+b*uk
        yk=c*xk+d*uk
        xk=xk_1
        if i==0:
            y=yk
        else:
            y=hstack((y,yk))
    y=array(y).T
    return y

def qhull(V, qstring):
    """
    Use qhull to determine convex hull / volume / normals.
     V - [matrix] vertices
     qstring - [string] arguments to pass to qhull
    """
    try:
        qhullp = subprocess.Popen(["qhull", qstring],
                              stdin=subprocess.PIPE, stdout=subprocess.PIPE)
        Vc = qhullp.communicate(qhullstr(V))[0] #qhull output to Vc
        
        if qstring == "FS": #calc area and volume
            ks = Vc.split('\n')[-2]
            Vol = float(ks.split(' ')[-2]) #get volume of D-hull
            return Vol
        elif qstring == "Ft": #calc vertices and facets
            ks = Vc.split('\n')
            fms = int(ks[1].split(' ')[1]) #get size of facet matrix
            fmat = ks[-fms-1:-1]
            fmat = mat(';'.join(fmat)) #generate matrix
            fmatv = fmat[:, 1:] #vertices on facets
            return array(fmatv)
        elif qstring == "n": #calc convex hull and get normals
            ks = ';'.join(Vc.split('\n')[2:]) #remove leading dimension output
            k = mat(ks[:-1]) #convert to martrix with vertices
            return array(k)
        else:
            exit(1)
    except:
        raise NameError('QhullError')

def con2vert(A, b):
    """
    Convert sets of constraints to a list of vertices (of the feasible region).
    If the shape is open, con2vert returns False for the closed property.
    """
    # Python implementation of con2vert.m by Michael Kleder (July 2005),
    #  available: http://www.mathworks.com/matlabcentral/fileexchange/7894
    #  -con2vert-constraints-to-vertices
    # Author: Michael Kelder (Original)
    #         Andre Campher (Python implementation)
    c = linalg.lstsq(mat(A), mat(b))[0]
    btmp = mat(b)-mat(A)*c
    D = mat(A)/matlib.repmat(btmp, 1, A.shape[1])

    fmatv = qhull(D, "Ft") #vertices on facets

    G  = zeros((fmatv.shape[0], D.shape[1]))
    for ix in range(0, fmatv.shape[0]):
        F = D[fmatv[ix, :], :].squeeze()
        G[ix, :] = linalg.lstsq(F, ones((F.shape[0], 1)))[0].transpose()

    V = G + matlib.repmat(c.transpose(), G.shape[0], 1)
    ux = uniqm(V)

    eps = 1e-13
    Av = dot(A, ux.T)
    bv = tile(b, (1, ux.shape[0]))
    closed = sciall(Av - bv <= eps)

    return ux, closed

def qnwcheb1(n, a, b):
    """ Univariate Gauss-Chebyshev quadrature nodes and weights

    Parameters
    -----------
    n : int
        number of nodes
    a : float
        left endpoint
    b : float
        right endpoint

    Returns
    ---------
    x : array, shape (n,)
        nodes
    x : array, shape (n,)
        weights

    Notes
    ---------
    
    Port of the qnwcheb1 function in the compecon matlab toolbox.
    """
    x = ((b + a) / 2 - (b - a) / 2
         * sp.cos(sp.pi / n * sp.arange(0.5, n + 0.5, 1)))
    w2 =  sp.r_[1, -2. / (sp.r_[1:(n - 1):2] * sp.r_[3:(n + 1):2])]
    w1 = (sp.cos(sp.pi / n * sp.mat((sp.r_[0:n] + 0.5)).T *
                 sp.mat((sp.r_[0:n:2]))).A)
    w0 = (b - a) / n
    w = w0 * sp.dot(w1, w2)
    return x, w

def get_projection_matrix(u,z,v,k=TOP_NUM_SINGULAR_VALUES):
    """
        generate the projection matrix which contains a 
        score vector for the patterns for each word pair 
    """ 
    #determine the best patterns for comparison   
    column_indexes=get_top_k_column_indexes(z)
    
    #using the column indexes of the best patterns recreate the 
    #u & z matrices containing only those correspoding columns
    
    #creating the uk matrix
    uk=[]
    for r in range(len(u)):
        uk.append([])
        for index in column_indexes:
            uk[r].append(u[r][index])
        r+=1

    #creating the zk matrix
    zk=[]
    for index in column_indexes:
        zk.append([])
        for col in range(len(v)):
            if (col==index):
                zk[len(zk)-1].append(z[index])
            else: 
                zk[len(zk)-1].append(0)
    
    #calcualte the projecttion matrix by u.z
    return mat(uk)*mat(zk)

def arnoldi(A, v0, k):
    """
    Arnoldi algorithm (Krylov approximation of a matrix)
        input: 
            A: matrix to approximate
            v0: initial vector (should be in matrix form) 
            k: number of Krylov steps 
        output: 
            V: matrix (large, N*k) containing the orthogonal vectors
            H: matrix (small, k*k) containing the Krylov approximation of A

    Author: Vasile Gradinaru, 14.12.2007 (Rennes)
    """
    #print 'ARNOLDI METHOD'
    inputtype = A.dtype.type
    V = mat( v0.copy() / norm(v0), dtype=inputtype)
    H = mat( zeros((k+1,k), dtype=inputtype) )
    for m in xrange(k):
        vt = A*V[ :, m]
        for j in xrange( m+1):
            H[ j, m] = (V[ :, j].H * vt )[0,0]
            vt -= H[ j, m] * V[:, j]
        H[ m+1, m] = norm(vt);
        if m is not k-1:
            V =  hstack( (V, vt.copy() / H[ m+1, m] ) ) 
    return V,  H

 def genInitSigmaFactor(self):
     """ depending on the algorithm settings, we start out with in identity matrix, or perturb it """
     if self.perturbedInitSigma:
         res = mat(eye(self.xdim)*self.initSigmaCoeff+randn(self.xdim, self.xdim)*self.initSigmaRandCoeff)            
     else:
         res = mat(eye(self.xdim)*self.initSigmaCoeff)
     return res   

def minreal(sys):
    """Minimal representation for state space systems

    Usage
    =====
    [sysmin]=minreal[sys]

    Inputs
    ------

    sys: system in ss or tf form

    Outputs
    -------
    sysfin: system in state space form
    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ni=shape(b)[1]
    no=shape(c)[0]

    out=tb03ad(nx,no,ni,a,b,c,d,'R')

    nr=out[3]
    A=out[0][:nr,:nr]
    B=out[1][:nr,:ni]
    C=out[2][:no,:nr]
    sysf=ss(A,B,C,sys.D,sys.Tsamp)
    return sysf

def accuracy(p, colname='ClassLabel'):
    # find out how many classes are in the experiment
    numSamples = {}      # this is a list containing list of labels
    labels = p.column(col=2)[1]   # get the column of the actual labels
    for l in labels:
        if not(l in numSamples.keys()):
            numSamples.update({l:0})
    numLabels = len(numSamples.keys()) 
    # count number of samples per class
    labels = p.column(col=2)[1]
    for l in labels:
        numSamples[l] = numSamples[l] + 1  
    numLabels = len(numSamples.keys())
    confusionMatrix = scipy.mat( numpy.zeros( (numLabels,numLabels) ) ) 
    hdr_actual = '%s-actual' % colname
    hdr_predic = '%s-prediction' % colname
    mapLabels = {}
    cnt = 0
    for l in numSamples.keys():
       mapLabels.update({l:cnt})
       cnt = cnt + 1
    confusionMatrix = scipy.mat( numpy.zeros( (numLabels,numLabels) ) )  
    hdr_actual = '%s-actual' % colname
    hdr_predic = '%s-prediction' % colname
    actualLabels = p.column(hdr_actual)
    predictLabels = p.column(hdr_predic)
    for cnt in range(0,len(actualLabels[0])):
        pLabel = mapLabels[predictLabels[1][cnt]]     # predicted label
        aLabel = mapLabels[actualLabels[1][cnt]]      # actual Label
        if ( not(  predictLabels[0][cnt] ==   actualLabels[0][cnt]   ) ):     # it is just a sanity check, this event should happen ever. Just to make sure that it is comparing the same subject
              assert False, "This event should NOT happend ever !!! Are you sure you are using the correct Pyxel version??? I am comparing labels of two different subjects !! " 
        confusionMatrix[aLabel,pLabel] = confusionMatrix[aLabel,pLabel] + 1.0
    return confusionMatrix

def dcgain(sys):
    """Return the steady state value of the step response os sys

    Usage
    =====
    dcgain=dcgain(sys)

    Inputs
    ------

    sys: system

    Outputs
    -------
    dcgain : steady state value
    """

    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    if sys.Tsamp!=0.0:
        a=a-eye(nx,nx)
    r=rank(a)
    if r<nx:
        gm=[]
    else:
        gm=-c*inv(a)*b+d
    return array(gm)

def full_obs(sys,poles):
    """Full order observer of the system sys

    Call:
    obs=full_obs(sys,poles)

    Parameters
    ----------
    sys : System in State Space form
    poles: desired observer poles

    Returns
    -------
    obs: ss
    Observer

    """
    if isinstance(sys, TransferFunction):
        "System must be in state space form"
        return
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    L=place(a.T,c.T,poles)
    L=mat(L).T
    Ao=a-L*c
    Bo=hstack((b-L*d,L))
    n=shape(Ao)
    m=shape(Bo)
    Co=eye(n[0],n[1])
    Do=zeros((n[0],m[1]))
    obs=StateSpace(Ao,Bo,Co,Do,sys.dt)
    return obs

def ned2ecef(lat, lon, alt, n, e, d):
    X0, Y0, Z0 = coord.geodetic2ecef(lat, lon, alt)
    lat, lon = radians(lat), radians(lon)
    
    pitch = math.pi/2 + lat
    yaw = -lon 
    
    my = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(pitch), 0, -sin(pitch),
         0,1,0,
         sin(pitch), 0, cos(pitch)))
    
    mz = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(yaw), sin(yaw),0,
         -sin(yaw),cos(yaw),0,
         0,0,1))
    
    mr = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (-cos(lon)*sin(lat), -sin(lon), -cos(lat) * cos(lon), 
         -sin(lat)*sin(lon), cos(lon), -sin(lon)*cos(lat),
         cos(lat), 0, -sin(lat)))
    
    geo = mat('[%f; %f; %f]' % (X0, Y0, Z0))
    ned = mat('[%f; %f; %f]' % (n, e, d))
    res = mr*ned + geo
    return res[0], res[1], res[2]