Python linalg.LinearOperator类代码示例

 def __init__(self, u, v):
     assert (len(u) == len(v))
     self.u, self.v = u, v
     if LINOP_SUBCLASSING:
         LinearOperator.__init__(self, dtype = float, shape = (len(u), len(v)))
     else:
         LinearOperator.__init__(self, dtype = float, shape = (len(u), len(v)), matvec = lambda x : type(self)._matvec(self, x))

def makeOperator( operatorInput, expectedShape ):
    """Internal. 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')

    if sys.version_info[0] >= 3:
        # special methods are looked up on the class -- so make a new one
        operator.__class__ = CallableLinearOperator
    else:
        operator.__call__ = operator.matmat

    return operator

def makeOperator( operatorInput, expectedShape ):
    """Internal. Takes a dense numpy array or a sparse matrix or
    a function and makes an operator performing matrix * blockvector
    products.

    Example
    -------

    >>> 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')

    operator.__call__ = operator.matmat

    return operator

 def __init__(self, F, Delta, u):
     n = F.shape[0]
     assert (F.shape[1] == n) and (u.shape == (n,))
     self.F, self.Delta, self.u = F, Delta, u
     self.u_prime = self.Delta - self.u
     self.ones = np.ones(n, dtype=float)
     LinearOperator.__init__(self, dtype=float, shape=self.F.shape)

    def __init__(self,
            A,
            coarsening=pyamgcl_ext.coarsening.smoothed_aggregation,
            relaxation=pyamgcl_ext.relaxation.spai0,
            prm={}
            ):
        """
        Class constructor.

        Creates algebraic multigrid hierarchy.

        Parameters
        ----------
        A : the system matrix in scipy.sparse format
        coarsening : {ruge_stuben, aggregation, *smoothed_aggregation*, smoothed_aggr_emin}
            The coarsening type to use for construction of the multigrid
            hierarchy.
        relaxation : {damped_jacobi, gauss_seidel, chebyshev, *spai0*, ilu0}
            The relaxation scheme to use for multigrid cycles.
        prm : dictionary with amgcl parameters
        """
        Acsr = A.tocsr()

        self.P = pyamgcl_ext.make_preconditioner(
                coarsening, relaxation, prm,
                Acsr.indptr.astype(numpy.int32),
                Acsr.indices.astype(numpy.int32),
                Acsr.data.astype(numpy.float64)
                )

        LinearOperator.__init__(self, A.shape, self.P)

    def __init__(self, shape, matvec):

        self._shape = shape
        self._action = matvec

        LinearOperator.__init__(self,
                                shape=self._shape,
                                matvec=self.parallelDot,
                                dtype=np.complex128)

 def __init__(self, K, quad, **kwargs):
     assert len(K.shape) == 1
     shape = (K.shape[0]-4, K.shape[0]-4)
     N = shape[0]
     LinearOperator.__init__(self, shape, None, **kwargs)
     ck = ones(K.shape)
     if quad == "GC": ck[N-1] = 2
     self.dd = -4*pi*(K[:N]+1)/(K[:N]+3)*(K[:N]+2)**2
     self.ud = [2*pi*(K[:N-2]+1)*(K[:N-2]+2)]        
     self.ld = [2*pi*(K[2:N]-1)*(K[2:N]+2)]

def Operator(f, size=None):
	"""Create a stacked version of the function f, casted as a LinearOperator.
	
	size is the dimension of the operator f."""
	
	if size == None: size = f.Size
	
	operator = LinearOperator((size,size), matvec = StackFunction(f), dtype=float64)
	operator.__call__ = lambda x: operator * x
	
	return operator

    def __init__(self, mat, r=10, leaf_side=16):
        LinearOperator.__init__(self, dtype=mat.dtype, shape=mat.shape,
                                matvec=self._matvec,
                                rmatvec=self._rmatvec)
        self.mat = BlackBox(mat)
        self.r = r
        self.leaf_side = leaf_side
        self.leaf_size = leaf_side**2

        N = int(np.sqrt(self.mat.shape[0]))

        self.pattern = distance_matrix(N, format='coo')
        perm = hilbert_traverse(N)
        conjugate_sparse(self.pattern, perm)
        self.pattern = self.pattern.tocsr()
        self.mat.permutate(perm)
        self.root = hmat_node(self, tuple(zip((0, 0), self.mat.shape)))
        return

def _diagonal_operator(diag):
    """Creates an operator representing a 
    multiplication with a diagonal matrix"""
    diag = diag.ravel()[:, np.newaxis]

    def diag_matvec(vec):
        if vec.ndim > 1:
            return diag * vec
        else:
            return diag.ravel() * vec

    linop = LinearOperator(shape=(len(diag), len(diag)),
                           matvec=diag_matvec,
                           rmatvec=diag_matvec,
                           dtype=np.float64)
    linop.matvec = diag_matvec
    linop.rmatvec = diag_matvec

    return linop

def get_grad_linop(X, Y, invcovB, invcovN, alpha):
    """
    Linear operator implementing the gradient of the functional
    \frac{1}{2} \|Y - XB\|^2_{\Sigma_n} + \frac{1}{2} \|B\|^2_{\Sigma_s}

    which reads
    grad_B = X^T(XB - Y)\Sigma_n^{-1} + \lambda B\Sigma_s^{-1}
    """

    N, P = X.shape
    T = invcovB.shape[0]

    if P <= N:
        XTX = aslinearoperator(X.T.dot(X))
        XTYinvcovN = invcovN.rmatvec(Y.T.dot(X)).T

        def matvec(vecB):
            XTXB = XTX.matvec(vecB.reshape(T, P).T)
            XTXB_invcovN = invcovN.rmatvec(XTXB.T).T
            B_incovB = invcovB.rmatvec(vecB.reshape(T, P)).T
            result = XTXB_invcovN - XTYinvcovN + alpha * B_incovB
            return result.T.ravel()
    else:
        # raise(Exception)
        def matvec(vecB):
            XB_minus_Y_invcovN = invcovN.rmatvec(
                (X.dot(vecB.reshape(T, P).T) - Y).T).T
            XT_XB_minus_Y_invcovN = X.T.dot(XB_minus_Y_invcovN)
            B_incovB = invcovB.rmatvec(vecB.reshape(T, P)).T
            result = XT_XB_minus_Y_invcovN + alpha * B_incovB
            return result.T.ravel()

    linop = LinearOperator(shape=tuple([X.shape[1] * Y.shape[1]] * 2),
                           matvec=matvec,
                           rmatvec=matvec,
                           dtype=np.dtype('float64'))
    linop.matvec = matvec
    linop.rmatvec = matvec
    linop.dtype = np.dtype('float64')

    return linop

def _woodbury_inverse(Ainv, Cinv, U, V):
    """Uses Woodbury Matrix Identity to invert the Matrix
    (A + UCV) ^ (-1)
    See http://en.wikipedia.org/wiki/Woodbury_matrix_identity"""

    def matvec(x):
        # this is probably wildly suboptimal, but it works
        Ainv_x = Ainv.matvec(x)
        Cinv_mat = Cinv.matvec(np.eye(Cinv.shape[0]))
        VAinvU = V.dot(Ainv.matvec(U))
        inv_Cinv_plus_VAinvU = np.linalg.inv(Cinv_mat + VAinvU)
        VAinv_x = V.dot(Ainv_x)
        inv_blabla_VAinv_x = inv_Cinv_plus_VAinvU.dot(VAinv_x)
        whole_big_block = Ainv.matvec(
            U.dot(inv_blabla_VAinv_x))
        return Ainv_x - whole_big_block

    shape = Ainv.shape
    linop = LinearOperator(shape=shape, matvec=matvec)
    linop.matvec = matvec
    linop.rmatvec = matvec
    return linop

    def __init__(self, A, pmask, prm={}):
        """
        Class constructor.

        Parameters
        ----------
        A : the system matrix in scipy.sparse format
        prm : dictionary with amgcl parameters
        """
        Acsr = A.tocsr()

        self.P = pyamgcl_ext.make_simple(
                prm,
                Acsr.indptr.astype(numpy.int32),
                Acsr.indices.astype(numpy.int32),
                Acsr.data.astype(numpy.float64),
                pmask.astype(numpy.int32)
                )

        if [int(v) for v in scipy.__version__.split('.')] < [0, 16, 0]:
            LinearOperator.__init__(self, A.shape, self.P)
        else:
            LinearOperator.__init__(self, dtype=numpy.float64, shape=A.shape)

def calculate_delta(gamma,     #3-dim array (Nz,Nx,Ny) of complex shear
                    P_kd,      #line-of-sight kappa-delta transform
                    P_gk,      #gamma-kappa transform
                    S_dd,      #expected signal covariance of delta
                               # (can be arbitrary form)
                    N,         #shear noise
                    alpha,     #weiner filter strength
                    M_factor = 1 #optional scaling of M before cg-method
                    ):    
    """
    Implement equation A3 from Simon & Taylor 2009
    Note: their notation Q -> our notation P_kd

    Also, here we factor N_d^-1 out of M: M should be Hermitian for
    the cg method. The Simon 09 expression is Hermitian only if N_d is
    proportional to the identity, which will not be the case for
    deweighted border pixels.
    """
    P_kd = as_Lens3D_matrix(P_kd)
    P_gk = as_Lens3D_matrix(P_gk)
    S_dd = as_Lens3D_matrix(S_dd)
    N = as_Lens3D_matrix(N)
    
    P_gk_cross = P_gk.conj_transpose()
    P_kd_T = P_kd.transpose()

    print "calculating delta:"

    print "  shape of P_gk:",P_gk.shape
    print "  shape of P_kd:",P_kd.shape

    print "constructing linear operator M"
    #define an operator which performs matrix-vector
    # multiplication representing M
    def matvec(v):
        v0 = P_gk_cross.view_as_Lens3D_vec(v)
        
        v1 = N.matvec(v0)
        v1 *= alpha

        v2 = P_gk_cross.matvec( v0 )
        v2 = P_kd_T.matvec( v2 )
        v2 = S_dd.matvec( v2 )
        v2 = P_kd.matvec( v2 )
        v2 = P_gk.matvec( v2 )

        ret = numpy.zeros(v.shape,dtype=complex)
        ret += v1.vec
        ret += v2.vec

        return P_gk_cross.view_as_same_type( ret * M_factor , v )

    M = LinearOperator(P_gk.shape,
                       matvec=matvec,
                       dtype=complex)

    v = numpy.random.random(M.shape[1])
    t0 = time()
    v2 = M.matvec(v)
    t = time()-t0
    print "  M multiplication: %.3g sec" % t

    #print M.matvec(numpy.ones(M.shape[1]))[:10]
    #exit()

    print "constructing preconditioner for M"

    #define an operator which can quickly approximate the inverse of
    # M using fourier-space inversions.  This inverse will be exact for
    # a noiseless reconstruction on an infinite field
    P_gk_I = P_gk.inverse(False)
    #P_kd_I = P_kd.inverse()
    S_dd_I = S_dd.inverse(False)
    #P_kd_I_T = P_kd_I.transpose()
    P_gk_I_cross = P_gk_I.conj_transpose()

    def matvec_pc(v):
        v0 = P_gk_I.view_as_Lens3D_vec(v)
        v0 = P_gk_I.matvec(v0)
        #v0 = P_kd_I.matvec(v0)
        v0 = S_dd_I.matvec(v0)
        #v0 = P_kd_I_T.matvec(v0)
        v0 = P_gk_I_cross.matvec(v0)
        return P_gk_I.view_as_same_type(v0,v)
    
    M_pc = LinearOperator( (M.shape[1],M.shape[0]), 
                           matvec = matvec_pc,
                           dtype = M.dtype )
    
    v = numpy.random.random(M_pc.shape[1])
    t0 = time()
    v3 = M_pc.matvec(v)
    t_pc = time()-t0
    print "  preconditioner multiplication: %.3g sec" % t_pc

    step1_vec = gamma.vec

    use_cg = True

    #---define callback function---
#.........这里部分代码省略.........

def estimate_condition_number(P_kd,
                              P_gk,
                              S_dd,
                              N,
                              alpha,
                              compute_exact = False):
    P_kd = as_Lens3D_matrix(P_kd)
    P_gk = as_Lens3D_matrix(P_gk)
    S_dd = as_Lens3D_matrix(S_dd)
    N = as_Lens3D_matrix(N)
    
    P_gk_cross = P_gk.conj_transpose()
    P_kd_T = P_kd.transpose()
    
    def matvec(v):
        v0 = P_gk_cross.view_as_Lens3D_vec(v)

        v2 = P_gk_cross.matvec( v0 )
        v2 = P_kd_T.matvec( v2 )
        v2 = S_dd.matvec( v2 )
        v2 = P_kd.matvec( v2 )
        v2 = P_gk.matvec( v2 )
        
        v1 = N.matvec(v0)
        v1 *= alpha

        ret = numpy.zeros(v.shape,dtype=complex)
        ret += v1.vec
        ret += v2.vec

        return P_gk_cross.view_as_same_type( ret , v )

    M = LinearOperator(P_gk.shape,
                       matvec=matvec,
                       dtype=complex)

    #compute the exact condition number
    if compute_exact:
        v = numpy.random.random(M.shape[1])
        t0 = time()
        v2 = M.matvec(v)
        t = time()-t0
        print " - constructing matrix representation (est. %s)" \
            % printtime(t*M.shape[0])
        t0 = time()
        M_rep = get_mat_rep(M)
        print "    time to get mat rep: %.2g sec" % (time()-t0)
        print " - computing SVD"
        t0 = time()
        sig = numpy.linalg.svd(M_rep,compute_uv=False)
        print "    time for SVD: %.2g sec" % (time()-t0)
        print 'true condition number:      %.2e / %.2e = %.2e' \
            % (sig[0],sig[-1],
               sig[0]/sig[-1])

    #estimate condition number, assuming the noiseless matrix
    # is rank-deficient.  This will be true if there are more
    # source lens-planes than mass lens-planes
    eval_max,evec_max = arpack.eigen(M,1)
    print 'estimated condition number: %.2e / %.2e = %.2e' \
        % ( abs(eval_max[0]) , numpy.min(N.data),
            abs(eval_max[0]) / numpy.min(N.data) )