Python linalg.LinearOperator方法代碼示例

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def si_c2(self,ww):
    """
    This computes the correlation part of the screened interaction using LinearOpt and lgmres
    lgmres method is much slower than np.linalg.solve !!
    """
    import numpy as np
    from scipy.sparse.linalg import lgmres
    from scipy.sparse.linalg import LinearOperator
    rf0 = si0 = self.rf0(ww)    
    for iw,w in enumerate(ww):                                
      k_c = np.dot(self.kernel_sq, rf0[iw,:,:])                                         
      b = np.dot(k_c, self.kernel_sq)               
      self.comega_current = w
      k_c_opt = LinearOperator((self.nprod,self.nprod), matvec=self.gw_vext2veffmatvec, dtype=self.dtypeComplex)  
      for m in range(self.nprod): 
         si0[iw,m,:],exitCode = lgmres(k_c_opt, b[m,:], atol=self.gw_iter_tol, maxiter=self.maxiter)   
      if exitCode != 0: print("LGMRES has not achieved convergence: exitCode = {}".format(exitCode))
      #np.allclose(np.dot(k_c, si0), b, atol=1e-05) == True  #Test   
    return si0 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def gw_comp_veff(self, vext, comega=1j*0.0):
    """
    This computes an effective field (scalar potential) given the external
    scalar potential as follows:
        (1-v\chi_{0})V_{eff} = V_{ext} = X_{a}^{n}V_{\mu}^{ab}X_{b}^{m} * 
                                         v\chi_{0}v * X_{a}^{n}V_{nu}^{ab}X_{b}^{m}
    
    returns V_{eff} as list for all n states(self.nn[s]).
    """
    
    from scipy.sparse.linalg import LinearOperator
    self.comega_current = comega
    veff_op = LinearOperator((self.nprod,self.nprod),
                             matvec=self.gw_vext2veffmatvec,
                             dtype=self.dtypeComplex)

    from scipy.sparse.linalg import lgmres
    resgm, info = lgmres(veff_op,
                         np.require(vext, dtype=self.dtypeComplex, requirements='C'),
                         atol=self.gw_iter_tol, maxiter=self.maxiter)
    if info != 0:
      print("LGMRES has not achieved convergence: exitCode = {}".format(info))
    return resgm 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def build_full(self, shape, fill_val):
        m, n = shape
        if fill_val == 0:
            return shape
        else:
            def matvec(v):
                return v.sum() * fill_val * np.ones(m)

            def rmatvec(v):
                return v.sum() * fill_val * np.ones(n)

            def matmat(M):
                return M.sum(axis=0) * fill_val * np.ones((m, M.shape[1]))

            return sla.LinearOperator(shape=shape,
                                      matvec=matvec,
                                      rmatvec=rmatvec,
                                      matmat=matmat,
                                      dtype=self.dtype) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def _makeOperator(operatorInput, expectedShape):
    """Takes a dense numpy array or a sparse matrix or
    a function and makes an operator performing matrix * blockvector
    products.

    Examples
    --------
    >>> A = _makeOperator( arrayA, (n, n) )
    >>> vectorB = A( vectorX )

    """
    if operatorInput is None:
        def ident(x):
            return x
        operator = LinearOperator(expectedShape, ident, matmat=ident)
    else:
        operator = aslinearoperator(operatorInput)

    if operator.shape != expectedShape:
        raise ValueError('operator has invalid shape')

    return operator 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def regularized_lsq_operator(J, diag):
    """Return a matrix arising in regularized least squares as LinearOperator.
    
    The matrix is
        [ J ]
        [ D ]
    where D is diagonal matrix with elements from `diag`.
    """
    J = aslinearoperator(J)
    m, n = J.shape

    def matvec(x):
        return np.hstack((J.matvec(x), diag * x))

    def rmatvec(x):
        x1 = x[:m]
        x2 = x[m:]
        return J.rmatvec(x1) + diag * x2

    return LinearOperator((m + n, n), matvec=matvec, rmatvec=rmatvec) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def right_multiply(J, d, copy=True):
    """Compute J diag(d).
    
    If `copy` is False, `J` is modified in place (unless being LinearOperator).
    """
    if copy and not isinstance(J, LinearOperator):
        J = J.copy()

    if issparse(J):
        J.data *= d.take(J.indices, mode='clip')  # scikit-learn recipe.
    elif isinstance(J, LinearOperator):
        J = right_multiplied_operator(J, d)
    else:
        J *= d

    return J 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def left_multiply(J, d, copy=True):
    """Compute diag(d) J.
    
    If `copy` is False, `J` is modified in place (unless being LinearOperator).
    """
    if copy and not isinstance(J, LinearOperator):
        J = J.copy()

    if issparse(J):
        J.data *= np.repeat(d, np.diff(J.indptr))  # scikit-learn recipe.
    elif isinstance(J, LinearOperator):
        J = left_multiplied_operator(J, d)
    else:
        J *= d[:, np.newaxis]

    return J 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def lsmr_operator(Jop, d, active_set):
    """Compute LinearOperator to use in LSMR by dogbox algorithm.

    `active_set` mask is used to excluded active variables from computations
    of matrix-vector products.
    """
    m, n = Jop.shape

    def matvec(x):
        x_free = x.ravel().copy()
        x_free[active_set] = 0
        return Jop.matvec(x * d)

    def rmatvec(x):
        r = d * Jop.rmatvec(x)
        r[active_set] = 0
        return r

    return LinearOperator((m, n), matvec=matvec, rmatvec=rmatvec, dtype=float) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def todense(self):
        """Return a dense array representation of this operator.

        Returns
        -------
        arr : ndarray, shape=(n, n)
            An array with the same shape and containing
            the same data represented by this `LinearOperator`.

        """
        s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho
        I = np.eye(*self.shape, dtype=self.dtype)
        Hk = I

        for i in range(n_corrs):
            A1 = I - s[i][:, np.newaxis] * y[i][np.newaxis, :] * rho[i]
            A2 = I - y[i][:, np.newaxis] * s[i][np.newaxis, :] * rho[i]

            Hk = np.dot(A1, np.dot(Hk, A2)) + (rho[i] * s[i][:, np.newaxis] *
                                                        s[i][np.newaxis, :])
        return Hk 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def lagrangian_hessian(self, z, v):
        """Returns scaled Lagrangian Hessian"""
        # Compute Hessian in relation to x and s
        Hx = self.lagrangian_hessian_x(z, v)
        if self.n_ineq > 0:
            S_Hs_S = self.lagrangian_hessian_s(z, v)

        # The scaled Lagragian Hessian is:
        #     [ Hx    0    ]
        #     [ 0   S Hs S ]
        def matvec(vec):
            vec_x = self.get_variables(vec)
            vec_s = self.get_slack(vec)
            if self.n_ineq > 0:
                return np.hstack((Hx.dot(vec_x), S_Hs_S*vec_s))
            else:
                return Hx.dot(vec_x)
        return LinearOperator((self.n_vars+self.n_ineq,
                               self.n_vars+self.n_ineq),
                              matvec) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def _check_reentrancy(solver, is_reentrant):
    def matvec(x):
        A = np.array([[1.0, 0, 0], [0, 2.0, 0], [0, 0, 3.0]])
        y, info = solver(A, x)
        assert_equal(info, 0)
        return y
    b = np.array([1, 1./2, 1./3])
    op = LinearOperator((3, 3), matvec=matvec, rmatvec=matvec,
                        dtype=b.dtype)

    if not is_reentrant:
        assert_raises(RuntimeError, solver, op, b)
    else:
        y, info = solver(op, b)
        assert_equal(info, 0)
        assert_allclose(y, [1, 1, 1])


#------------------------------------------------------------------------------ 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def get_subset_lin_op(lin_op, sub_idx):
    """ Subset a linear operator to the indices in `sub_idx`. Equivalent to A' = A[sub_idx, :]
    :param LinearOperator lin_op: input linear operator
    :param np.ndarray[int] sub_idx: subset index
    :return: the subset linear operator
    :rtype: LinearOperator
    """
    if lin_op is None:
        return None
    if type(lin_op) is IndexOperator:
        # subsetting IndexOperator yields a new IndexOperator
        return IndexOperator(lin_op.index_map[sub_idx], dim_x=lin_op.dim_x)
    elif isinstance(lin_op, MatrixLinearOperator):
        # subsetting a matrix multiplication operation yields a new matrix
        return MatrixLinearOperator(lin_op.A[sub_idx, :])
    # in the general case, append a sub-indexing operator
    return IndexOperator(sub_idx, dim_x=lin_op.shape[0]) * lin_op 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def rmatvec_nd(lin_op, x):
    """
    Project a 1D or 2D numpy or sparse array using rmatvec. This is different from rmatvec
    because it applies rmatvec to each row and column. If x is n x n and lin_op is n x k,
    the result will be k x k.

    :param LinearOperator lin_op: The linear operator to apply to x
    :param np.ndarray|sp.spmatrix x: array/matrix to be projected
    :return: the projected array
    :rtype: np.ndarray|sp.spmatrix
    """
    if x is None or lin_op is None:
        return x
    if isinstance(x, sp.spmatrix):
        y = x.toarray()
    elif np.isscalar(x):
        y = np.array(x, ndmin=1)
    else:
        y = np.copy(x)
    proj_func = lambda z: lin_op.rmatvec(z)
    for j in range(y.ndim):
        if y.shape[j] == lin_op.shape[0]:
            y = np.apply_along_axis(proj_func, j, y)
    return y 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def _check_reentrancy(solver, is_reentrant):
    def matvec(x):
        A = np.array([[1.0, 0, 0], [0, 2.0, 0], [0, 0, 3.0]])
        y, info = solver(A, x)
        assert_equal(info, 0)
        return y
    b = np.array([1, 1./2, 1./3])
    op = LinearOperator((3, 3), matvec=matvec, rmatvec=matvec,
                        dtype=b.dtype)

    if not is_reentrant:
        assert_raises(RuntimeError, solver, op, b)
    else:
        y, info = solver(op, b)
        assert_equal(info, 0)
        assert_allclose(y, [1, 1, 1]) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import LinearOperator [as 別名]
def test_eigs_for_k_greater():
    # Test eigs() for k beyond limits.
    A_sparse = diags([1, -2, 1], [-1, 0, 1], shape=(4, 4))  # sparse
    A = generate_matrix(4, sparse=False)
    M_dense = np.random.random((4, 4))
    M_sparse = generate_matrix(4, sparse=True)
    M_linop = aslinearoperator(M_dense)
    eig_tuple1 = eig(A, b=M_dense)
    eig_tuple2 = eig(A, b=M_sparse)

    with suppress_warnings() as sup:
        sup.filter(RuntimeWarning)

        assert_equal(eigs(A, M=M_dense, k=3), eig_tuple1)
        assert_equal(eigs(A, M=M_dense, k=4), eig_tuple1)
        assert_equal(eigs(A, M=M_dense, k=5), eig_tuple1)
        assert_equal(eigs(A, M=M_sparse, k=5), eig_tuple2)

        # M as LinearOperator
        assert_raises(TypeError, eigs, A, M=M_linop, k=3)

        # Test 'A' for different types
        assert_raises(TypeError, eigs, aslinearoperator(A), k=3)
        assert_raises(TypeError, eigs, A_sparse, k=3)