Python linalg.lsqr方法代码示例

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
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) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
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) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_Regression(par):
    """Dot-test, inversion and apply for Regression operator
    """
    np.random.seed(10)
    order = 4
    t = np.arange(par['ny'], dtype=np.float32)
    LRop = Regression(t, order=order, dtype=par['dtype'])
    assert dottest(LRop, par['ny'], order+1)

    x = np.array([1., 2., 0., 3., -1.], dtype=np.float32)
    xlsqr = lsqr(LRop, LRop*x, damp=1e-10, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)

    y = LRop * x
    y1 = LRop.apply(t, x)
    assert_array_almost_equal(y, y1, decimal=3) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_MatrixMult_repeated(par):
    """Dot-test and inversion for test_MatrixMult operator repeated
    along another dimension
    """
    np.random.seed(10)
    G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
        par['imag'] * np.random.normal(0, 10, (par['ny'],
                                               par['nx'])).astype('float32')
    Gop = MatrixMult(G, dims=5, dtype=par['dtype'])
    assert dottest(Gop, par['ny']*5, par['nx']*5,
                   complexflag=0 if par['imag'] == 1 else 3)

    x = (np.ones((par['nx'], 5)) +
         par['imag'] * np.ones((par['nx'], 5))).flatten()
    xlsqr = lsqr(Gop, Gop*x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=4) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_Kroneker(par):
    """Dot-test, inversion and comparison with np.kron for Kronecker operator
    """
    np.random.seed(10)
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    x = np.ones(par['nx']**2) + par['imag']*np.ones(par['nx']**2)

    Kop = Kronecker(MatrixMult(G1, dtype=par['dtype']),
                    MatrixMult(G2, dtype=par['dtype']),
                    dtype=par['dtype'])
    assert dottest(Kop, par['ny']**2, par['nx']**2,
                   complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(Kop, Kop * x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=2)

    # Comparison with numpy
    assert_array_almost_equal(np.kron(G1, G2), Kop * np.eye(par['nx']**2),
                              decimal=3) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_Fredholm1(par):
    """Dot-test and inversion for Fredholm1 operator
    """
    np.random.seed(10)

    _F = np.arange(par['nsl'] * par['nx'] * par['ny']).reshape(par['nsl'],
                                                               par['nx'],
                                                               par['ny'])
    F = _F - par['imag'] * _F

    x = np.ones((par['nsl'], par['ny'], par['nz'])) + \
        par['imag'] * np.ones((par['nsl'], par['ny'], par['nz']))

    Fop = Fredholm1(F, nz=par['nz'], saveGt=par['saveGt'],
                    usematmul=par['usematmul'], dtype=par['dtype'])
    assert dottest(Fop, par['nsl']*par['nx']*par['nz'],
                   par['nsl']*par['ny']*par['nz'],
                   complexflag=0 if par['imag'] == 0 else 3)
    xlsqr = lsqr(Fop, Fop * x.flatten(), damp=1e-20,
                 iter_lim=30, show=0)[0]
    xlsqr = xlsqr.reshape(par['nsl'], par['ny'], par['nz'])
    assert_array_almost_equal(x, xlsqr, decimal=3) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_Diagonal_2dsignal(par):
    """Dot-test and inversion for Diagonal operator for 2d signal
    """
    for idim, ddim in enumerate((par['nx'], par['nt'])):
        d = np.arange(ddim) + 1. +\
            par['imag'] * (np.arange(ddim) + 1.)

        Dop = Diagonal(d, dims=(par['nx'], par['nt']),
                       dir=idim, dtype=par['dtype'])
        assert dottest(Dop, par['nx']*par['nt'], par['nx']*par['nt'],
                       complexflag=0 if par['imag'] == 0 else 3)

        x = np.ones((par['nx'], par['nt'])) + \
            par['imag']*np.ones((par['nx'], par['nt']))
        xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]

        assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_Diagonal_3dsignal(par):
    """Dot-test and inversion for Diagonal operator for 3d signal
    """
    for idim, ddim in enumerate((par['ny'], par['nx'], par['nt'])):
        d = np.arange(ddim) + 1. +\
            par['imag'] * (np.arange(ddim) + 1.)

        Dop = Diagonal(d, dims=(par['ny'], par['nx'], par['nt']),
                       dir=idim, dtype=par['dtype'])
        assert dottest(Dop, par['ny']*par['nx']*par['nt'],
                       par['ny']*par['nx']*par['nt'],
                       complexflag=0 if par['imag'] == 0 else 3)

        x = np.ones((par['ny'], par['nx'], par['nt'])) + \
            par['imag']*np.ones((par['ny'], par['nx'], par['nt']))
        xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]

        assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_DWT_1dsignal(par):
    """Dot-test and inversion for DWT operator for 1d signal
    """
    DWTop = DWT(dims=[par['nt']], dir=0, wavelet='haar', level=3)
    x = np.random.normal(0., 1., par['nt']) + \
        par['imag'] * np.random.normal(0., 1., par['nt'])

    assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                   complexflag=0 if par['imag'] == 0 else 3)

    y = DWTop * x
    xadj = DWTop.H * y  # adjoint is same as inverse for dwt
    xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

    assert_array_almost_equal(x, xadj, decimal=8)
    assert_array_almost_equal(x, xinv, decimal=8) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_DWT_2dsignal(par):
    """Dot-test and inversion for DWT operator for 2d signal
    """
    for dir in [0, 1]:
        DWTop = DWT(dims=(par['nt'], par['nx']),
                    dir=dir, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_DWT_3dsignal(par):
    """Dot-test and inversion for DWT operator for 3d signal
    """
    for dir in [0, 1, 2]:
        DWTop = DWT(dims=(par['nt'], par['nx'], par['ny']),
                    dir=dir, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_DWT2D_2dsignal(par):
    """Dot-test and inversion for DWT2D operator for 2d signal
    """
    DWTop = DWT2D(dims=(par['nt'], par['nx']),
                dirs=(0, 1), wavelet='haar', level=3)
    x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
        par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))

    assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                   complexflag=0 if par['imag'] == 0 else 3)

    y = DWTop * x.ravel()
    xadj = DWTop.H * y  # adjoint is same as inverse for dwt
    xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

    assert_array_almost_equal(x.ravel(), xadj, decimal=8)
    assert_array_almost_equal(x.ravel(), xinv, decimal=8) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
def test_DWT2D_3dsignal(par):
    """Dot-test and inversion for DWT operator for 3d signal
    """
    for dirs in [(0, 1), (0, 2), (1, 2)]:
        DWTop = DWT2D(dims=(par['nt'], par['nx'], par['ny']),
                      dirs=dirs, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
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) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import lsqr [as 别名]
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), dtype=X.dtype)
    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