Python LinearOperator.matvec方法代码示例

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

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

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

# 需要导入模块: from scipy.sparse.linalg import LinearOperator [as 别名]
# 或者: from scipy.sparse.linalg.LinearOperator import matvec [as 别名]
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---
#.........这里部分代码省略.........

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

# 需要导入模块: from scipy.sparse.linalg import LinearOperator [as 别名]
# 或者: from scipy.sparse.linalg.LinearOperator import matvec [as 别名]
        pylab.colorbar()

    pylab.show()
    exit()

#test Sigma->gamma
if False:
    print "testing Sigma->gamma"
    def matvec(v):
        v1 = P_kd.matvec(v)
        return P_gk.matvec(v1)

    M = LinearOperator(P_kd.shape,matvec=matvec,dtype=complex)
    
    kappa_test = P_kd.matvec(Sigma)
    gamma_test = M.matvec(Sigma.vec)
    gamma_test = P_kd.view_as_Lens3D_vec(gamma_test)
    
    for i in range(gamma.Nz):
        pylab.figure( figsize=(6,8) )
        pylab.subplot(211)
        kappa.imshow_lens_plane(i)
        gamma.fieldplot_lens_plane(i)
        pylab.title(r"$\kappa,\gamma\ \rm{true}\ (z=%.2f)$" % z_gamma[i])
        pylab.subplot(212)
        kappa_test.imshow_lens_plane(i)
        gamma_test.fieldplot_lens_plane(i)
        pylab.title(r"$\kappa,\gamma\ \rm{test}\ (z=%.2f)$" % z_gamma[i])
    pylab.show()
    exit()

# 需要导入模块: from scipy.sparse.linalg import LinearOperator [as 别名]
# 或者: from scipy.sparse.linalg.LinearOperator import matvec [as 别名]
	def make_kle(self, KLE_sigma, KLE_L,
				KLE_num_eig,
				nrow, ncol, nlay, Lx, Ly):

		print "making KLE, this will take a long time for large dimensions"

		print KLE_sigma, KLE_L, KLE_num_eig, nrow, ncol, nlay, Lx, Ly


		# initialise swig wrapped cpp class

		C_matrix = correlation.C_matrix(KLE_sigma, KLE_L)
		C_matrix.set_dims(nrow, ncol, Lx, Ly)

		out_vec = np.zeros((nrow*ncol), 'd')
		def KLE_Av(v):
			C_matrix.av_no_C(nrow*ncol, v, out_vec)
			return out_vec

		KLE_A = LinearOperator((nrow*ncol,nrow*ncol), matvec=KLE_Av, dtype='d')

		t1 = time.time()
		eig_vals, eig_vecs = eigsh( KLE_A, k=KLE_num_eig)
		t2 = time.time()


		# sometimes one gets -v rather than v?
		for i in range(KLE_num_eig):
			print "NORM", np.linalg.norm(eig_vecs[:,i])
			assert np.allclose(KLE_A.matvec(eig_vecs[:,i]), eig_vals[i]*eig_vecs[:,i])

		print "=================================="
		print "SVD took ", t2-t1, "seconds and "
		print "they seem to indeed be eigen vectors"
		print "=================================="


# plot eigenvectors
		from mpl_toolkits.mplot3d import axes3d
		import matplotlib.pyplot as plt
		fig = plt.figure(figsize=plt.figaspect(0.2))
		for i in range(eig_vecs.shape[1]):
			ax = fig.add_subplot(1, eig_vecs.shape[1], i, projection='3d')
			x = np.arange(0, nrow*Lx, Lx)
			X = np.empty((nrow, ncol))
			for col in range(ncol):
				X[:,col] = x
			y = np.arange(0, ncol*Ly, Ly)
			Y = np.empty((nrow, ncol))
			for row in range(nrow):
				Y[row,:] = y
			Z = eig_vecs[:,i].reshape((nrow, ncol))
			ax.plot_wireframe(X, Y, Z)
		plt.show()

# plot eigenvals

		# import matplotlib.pyplot as plt
		plt.plot(eig_vals, 'o-')
		plt.show()

		print "eig_vals", eig_vals

		np.savez(working_directory+'kle.npz', eig_vals=eig_vals, eig_vecs=eig_vecs, KLE_L=KLE_L, KLE_sigma=KLE_sigma, KLE_num_eig=KLE_num_eig, nrow=nrow, ncol=ncol, Lx=Lx, Ly=Ly)

		return eig_vals, eig_vecs

# 需要导入模块: from scipy.sparse.linalg import LinearOperator [as 别名]
# 或者: from scipy.sparse.linalg.LinearOperator import matvec [as 别名]
# how you'd use it 
##########################################

# mu = np.zeros(nrow*ncol) # could be the mean of any data you have
# scale = 1.
# # evaluate KLE with given modes
# # other parameters (other than modes) are mu, scale, eig_vecs, eig_vals (or really sigma and L on which they depend)
# def KLE(modes):
#     coefs = np.sqrt(eig_vals) * modes # elementwise
#     truncated_M = mu + scale * np.dot( eig_vecs, coefs)
#     unflattened = truncated_M.reshape(nrow,ncol)
#     return unflattened

# modes = np.ones(num_eig)
# print KLE(modes)

############################################


print "=================================="
print "done in ", t1-t0, "seconds"
print "=================================="


for i in range(num_eig):
	assert np.allclose(A.matvec(eig_vecs[:,i]), eig_vals[i]*eig_vecs[:,i]) 

print "=================================="
print "they are indeed eigenvectors"
print "=================================="