Python scipy.sparse.linalg 模块，lsqr() 实例源码

def solveLaplacianMesh(mesh, anchors, anchorsIdx, cotangent=True):
    n = mesh.VPos.shape[0] # N x 3
    k = anchorsIdx.shape[0]

    operator = (getLaplacianMatrixUmbrella, getLaplacianMatrixCotangent)

    L = operator[1](mesh, anchorsIdx) if cotangent else operator[0](mesh, anchorsIdx)
    delta = np.array(L.dot(mesh.VPos))

    # augment delta solution matrix with weighted anchors
    for i in range(k):
        delta[n + i, :] = WEIGHT * anchors[i, :]

    # update mesh vertices with least-squares solution
    for i in range(3):
        mesh.VPos[:, i] = lsqr(L, delta[:, i])[0]

    return mesh
#Purpose: Given a few RGB colors on a mesh, smoothly interpolate those colors
#by using their values as anchors and 
#Inputs: mesh (polygon mesh object), anchors (a K x 3 numpy array of anchor
#coordinates), colorsIdx (a parallel array of the indices of the RGB anchor indices)

def smoothColors(mesh, anchors, colorsIdx):
    colorsIdx = np.array(colorsIdx)
    n = mesh.VPos.shape[0] 
    k = anchors.shape[0]
    colors = np.zeros((n, 3))
    delta = np.zeros((n + k, 3))

    L = getLaplacianMatrixUmbrella(mesh, colorsIdx);

     # augment delta solution matrix with weighted anchors
    for i in range(k):
        delta[n + i, :] = WEIGHT * anchors[i, :]

    # update RGB values with least-squares solution
    for i in range(3):
        colors[:, i] = lsqr(L, delta[:, i])[0]

    return colors

#Purpose: Given a mesh, to smooth it by subtracting off the delta coordinates
#from each vertex, normalized by the degree of that vertex
#Inputs: mesh (polygon mesh object)
#Returns: Nothing (should update mesh.VPos)

def get_scores_and_p_values(self, tdm, category):
        '''
        Parameters
        ----------
        tdm: TermDocMatrix
        category: str, category name

        Returns
        -------
        pd.DataFrame(['coef', 'p-val'])
        '''
        X = tdm._X
        y = self._make_response_variable_1_or_negative_1(category, tdm)
        pX = X / X.sum(axis=1)
        ansX = self._anscombe_transform(pX.copy())
        B, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var\
            = lsqr(A=ansX, b=y, calc_var=True)

def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
    n_samples, n_features = X.shape
    coefs = np.empty((y.shape[1], n_features))
    n_iter = np.empty(y.shape[1], dtype=np.int32)

    # According to the lsqr documentation, alpha = damp^2.
    sqrt_alpha = np.sqrt(alpha)

    for i in range(y.shape[1]):
        y_column = y[:, i]
        info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
                              atol=tol, btol=tol, iter_lim=max_iter)
        coefs[i] = info[0]
        n_iter[i] = info[2]

    return coefs, n_iter

def makeMinimalSurface(mesh, anchors, anchorsIdx):
    n = mesh.VPos.shape[0] # N x 3
    k = anchorsIdx.shape[0]
    I = []
    J = []
    V = []

    # Build sparse Laplacian Matrix coordinates and values
    for i in range(n):
        neighbors = mesh.vertices[i].getVertexNeighbors()
        indices = map(lambda x: x.ID, neighbors)

        if i in anchorsIdx:
            I = I + [i]
            J = J + [i]
            V = V + [WEIGHT]
        else:    
            z = len(indices)
            I = I + ([i] * (z + 1)) # repeated row
            J = J + indices + [i] # column indices and this row
            V = V + ([-1 / float(z)] * z) + [1] # negative weights divided by degree and row degree

    L = sparse.coo_matrix((V, (I, J)), shape=(n, n)).tocsr()
    delta = np.zeros((n, 3))
    delta[np.array(anchorsIdx), :] = WEIGHT * anchors

    # update mesh vertices with least-squares solution
    for i in range(3):
        mesh.VPos[:, i] = lsqr(L, delta[:, i])[0]

    return mesh


##############################################################
##        Spectral Representations / Heat Flow              ##
##############################################################

#Purpose: Given a mesh, to compute first K eigenvectors of its Laplacian
#and the corresponding eigenvalues
#Inputs: mesh (polygon mesh object), K (number of eigenvalues/eigenvectors)
#Returns: (eigvalues, eigvectors): a tuple of the eigenvalues and eigenvectors

def getTexCoords(mesh, quadIdx):
    n = mesh.VPos.shape[0] # N x 3
    k = np.array(quadIdx).shape[0]
    I = []
    J = []
    V = []
    anchors = np.array([[0, 0, 0],
                        [0, 1, 0],
                        [1, 1, 0],
                        [1, 0, 0]])
    U = np.zeros((n, 2))

    # Build sparse Laplacian Matrix coordinates and values
    for i in range(n):
        vertex = mesh.vertices[i]
        neighbors = mesh.vertices[i].getVertexNeighbors()
        indices = map(lambda x: x.ID, neighbors)

        if i in quadIdx:
            I = I + [i]
            J = J + [i]
            V = V + [WEIGHT]
        else:    
            z = len(indices)
            I = I + ([i] * (z + 1)) # repeated row
            J = J + indices + [i] # column indices and this row
            V = V + ([-1 / max(float(z), 1)] * z) + [1] # negative weights divided by degree and row degree

    L = sparse.coo_matrix((V, (I, J)), shape=(n, n)).tocsr()
    delta = np.zeros((n, 3))
    delta[np.array(quadIdx), :] = WEIGHT * anchors

    # update mesh vertices with least-squares solution
    for i in range(2): #(only X and Y cooridinates)
        U[:, i] = lsqr(L, delta[:, i])[0]

    return U

def doHodge(R, W, Y, verbose = False):
    """
    Given 
    :param R: NEdges x 2 matrix specfiying comparisons that have been made
    :param W: A flat array of NEdges weights parallel to the rows of R
    :param Y: A flat array of NEdges specifying preferences
    :returns: (s, I, H): s is scalar function, I is local inconsistency vector,
        H is global inconsistency vector
    """
    #Step 1: Get s
    if verbose:
        print("Making Delta0...")
    tic = time.time()
    D0 = makeDelta0(R)
    toc = time.time()
    if verbose:
        print("Elapsed Time: %g seconds"%(toc-tic))
    wSqrt = np.sqrt(W).flatten()
    WSqrt = scipy.sparse.spdiags(wSqrt, 0, len(W), len(W))
    WSqrtRecip = scipy.sparse.spdiags(1/wSqrt, 0, len(W), len(W))
    A = WSqrt*D0
    b = WSqrt.dot(Y)
    s = lsqr(A, b)[0]

    #Step 2: Get local inconsistencies
    if verbose:
        print("Making Delta1...")
    tic = time.time()
    D1 = makeDelta1(R)
    toc = time.time()
    if verbose:
        print("Elapsed Time: %g seconds"%(toc-tic))
    B = WSqrtRecip*D1.T
    resid = Y - D0.dot(s)  #This has been verified to be orthogonal under <resid, D0*s>_W

    u = wSqrt*resid
    if verbose:
        print("Solving for Phi...")
    tic = time.time()
    Phi = lsqr(B, u)[0]
    toc = time.time()
    if verbose:
        print("Elapsed Time: %g seconds"%(toc-tic))
    I = WSqrtRecip.dot(B.dot(Phi)) #Delta1* dot Phi, since Delta1* = (1/W) Delta1^T

    #Step 3: Get harmonic cocycle
    H = resid - I
    return (s, I, H)

def doFlattening(mesh, quadIdx):
    n = mesh.VPos.shape[0] # N x 3
    k = np.array(quadIdx).shape[0]
    I = []
    J = []
    V = []
    anchors = np.array([[0, 0, 0],
                        [0, 1, 0],
                        [1, 1, 0],
                        [1, 0, 0]])
    # Build sparse Laplacian Matrix coordinates and values
    for i in range(n):
        vertex = mesh.vertices[i]
        neighbors = mesh.vertices[i].getVertexNeighbors()
        indices = map(lambda x: x.ID, neighbors)

        if i in quadIdx:
            I = I + [i]
            J = J + [i]
            V = V + [WEIGHT]
        else:    
            z = len(indices)
            I = I + ([i] * (z + 1)) # repeated row
            J = J + indices + [i] # column indices and this row
            V = V + ([-1 / float(z)] * z) + [1] # negative weights divided by degree and row degree

    L = sparse.coo_matrix((V, (I, J)), shape=(n, n)).tocsr()
    delta = np.zeros((n, 3))
    delta[np.array(quadIdx), :] = WEIGHT * anchors

    # update mesh vertices with least-squares solution
    for i in range(3):
        mesh.VPos[:, i] = lsqr(L, delta[:, i])[0]

    return mesh

#Purpose: Given 4 vertex indices on a quadrilateral, to anchor them to the 
#square and flatten the rest of the mesh inside of that square.  Then, to 
#return these to be used as texture coordinates
#Inputs: mesh (polygon mesh object), quadIdxs (a length 4 array of indices
#into the mesh of the four points that are to be anchored, in CCW order)
#Returns: U (an N x 2 matrix of texture coordinates)

def solve(A, b, method, tol=1e-3):
    """ General sparse solver interface.

    method can be one of
    - spsolve_umfpack_mmd_ata
    - spsolve_umfpack_colamd
    - spsolve_superlu_mmd_ata
    - spsolve_superlu_colamd
    - bicg
    - bicgstab
    - cg
    - cgs
    - gmres
    - lgmres
    - minres
    - qmr
    - lsqr
    - lsmr
    """

    if method == 'spsolve_umfpack_mmd_ata':
        return spla.spsolve(A,b,use_umfpack=True, permc_spec='MMD_ATA')
    elif method == 'spsolve_umfpack_colamd':
        return spla.spsolve(A,b,use_umfpack=True, permc_spec='COLAMD')
    elif method == 'spsolve_superlu_mmd_ata':
        return spla.spsolve(A,b,use_umfpack=False, permc_spec='MMD_ATA')
    elif method == 'spsolve_superlu_colamd':
        return spla.spsolve(A,b,use_umfpack=False, permc_spec='COLAMD')
    elif method == 'bicg':
        res = spla.bicg(A,b,tol=tol)
        return res[0]
    elif method == 'bicgstab':
        res = spla.bicgstab(A,b,tol=tol)
        return res[0]
    elif method == 'cg':
        res = spla.cg(A,b,tol=tol)
        return res[0]
    elif method == 'cgs':
        res = spla.cgs(A,b,tol=tol)
        return res[0]
    elif method == 'gmres':
        res = spla.gmres(A,b,tol=tol)
        return res[0]
    elif method == 'lgmres':
        res = spla.lgmres(A,b,tol=tol)
        return res[0]
    elif method == 'minres':
        res = spla.minres(A,b,tol=tol)
        return res[0]
    elif method == 'qmr':
        res = spla.qmr(A,b,tol=tol)
        return res[0]
    elif method == 'lsqr':
        res = spla.lsqr(A,b,atol=tol,btol=tol)
        return res[0]
    elif method == 'lsmr':
        res = spla.lsmr(A,b,atol=tol,btol=tol)
        return res[0]
    else:
        raise Exception('UnknownSolverType')