[Scipy-svn] r3459 - in trunk/scipy/sandbox/multigrid: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 24 00:20:15 EDT 2007
Author: wnbell
Date: 2007-10-23 23:19:51 -0500 (Tue, 23 Oct 2007)
New Revision: 3459
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
aSA works reasonably well now
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-23 04:08:10 UTC (rev 3458)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-24 04:19:51 UTC (rev 3459)
@@ -1,6 +1,7 @@
import numpy,scipy,scipy.sparse
from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
- asarray, hstack
+ asarray, hstack, ascontiguousarray, isinf
+from numpy.random import randn
from scipy.sparse import csr_matrix,coo_matrix
from relaxation import gauss_seidel
@@ -10,11 +11,12 @@
-def augment_candidates(AggOp, old_Q, old_R, new_candidate, level):
+def augment_candidates(AggOp, old_Q, old_R, new_candidate):
K = old_R.shape[1]
#determine blocksizes
- if level == 0:
+ if new_candidate.shape[0] == old_Q.shape[0]:
+ #then this is the first prolongator
old_bs = (1,K)
new_bs = (1,K+1)
else:
@@ -22,42 +24,51 @@
new_bs = (K+1,K+1)
AggOp = expand_into_blocks(AggOp,new_bs[0],1).tocsr() #TODO switch to block matrix
-
+
+
# tentative prolongator
#TODO USE BSR
Q_indptr = (K+1)*AggOp.indptr
Q_indices = ((K+1)*AggOp.indices).repeat(K+1)
for i in range(K+1):
Q_indices[i::K+1] += i
- Q_data = zeros((AggOp.shape[0]/new_bs[0],) + new_bs)
- Q_data[:,:new_bs[0],:new_bs[1]] = old_Q.data.reshape((-1,) + old_bs) #TODO BSR change
+ Q_data = zeros((AggOp.indptr[-1]/new_bs[0],) + new_bs)
+ Q_data[:,:old_bs[0],:old_bs[1]] = old_Q.data.reshape((-1,) + old_bs) #TODO BSR change
# coarse candidates
- R = zeros(((K+1)*AggOp.shape[1],K+1))
- for i in range(K):
- R[i::K+1,:K] = old_R[i::K,:]
+ #R = zeros(((K+1)*AggOp.shape[1],K+1))
+ #for i in range(K):
+ # R[i::K+1,:K] = old_R[i::K,:]
+ R = zeros((AggOp.shape[1],K+1,K+1))
+ R[:,:K,:K] = old_R.reshape(-1,K,K)
c = new_candidate.reshape(-1)[diff(AggOp.indptr) == 1] #eliminate DOFs that aggregation misses
+ threshold = 1e-10 / abs(c).max() # cutoff for small basis functions
+
X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=AggOp.shape)
#orthogonalize X against previous
for i in range(K):
old_c = ascontiguousarray(Q_data[:,:,i].reshape(-1))
D_AtX = csr_matrix((old_c*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
- R[i::K+1,-1] = D_AtX
+ R[:,i,K] = D_AtX
X.data -= D_AtX[X.indices] * old_c
#normalize X
D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
col_norms = sqrt(D_XtX)
- R[K::K+1,-1] = col_norms
-
+ mask = col_norms < threshold # find small basis functions
+ col_norms[mask] = 0 # and set them to zero
+
+ R[:,K,K] = col_norms # store diagonal entry into R
+
col_norms = 1.0/col_norms
- col_norms[isinf(col_norms)] = 0
+ col_norms[mask] = 0
X.data *= col_norms[X.indices]
- last_column = Q_data[:,:,-1].reshape(-1)
- last_column = X.data.reshape(-1)
+ Q_data[:,:,-1] = X.data.reshape(-1,1)
+
Q_data = Q_data.reshape(-1) #TODO BSR change
+ R = R.reshape(-1,K+1)
Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(AggOp.shape[0],(K+1)*AggOp.shape[1]))
@@ -171,8 +182,8 @@
max_coarse = max_coarse, \
mu = mu, epsilon = epsilon)
- Ws = AggOps
-
+ #Ws = AggOps
+ x /= abs(x).max()
fine_candidates = x
#create SA using x here
@@ -186,6 +197,7 @@
#TODO which is faster?
x = self.__develop_candidate(As,Ps,Ts,AggOps,self.candidates,mu=mu)
+ x /= abs(x).max()
fine_candidates = hstack((fine_candidates,x))
As,Ps,Ts,self.candidates = sa_hierarchy(A,AggOps,fine_candidates)
@@ -204,7 +216,7 @@
#step 1
A_l = A
- x = scipy.rand(A_l.shape[0],1)
+ x = randn(A_l.shape[0],1)
skip_f_to_i = False
#step 2
@@ -251,7 +263,7 @@
def __develop_candidate(self,As,Ps,Ts,AggOps,candidates,mu):
A = As[0]
- x = A*scipy.rand(A.shape[0],1)
+ x = randn(A.shape[0],1)
b = zeros_like(x)
x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
@@ -270,22 +282,22 @@
for i in range(len(As) - 2):
#TODO test augment_candidates against fit candidates
- if i == 0:
- temp = hstack((candidates[0],x))
- else:
- K = candidates[i].shape[1]
- temp = zeros((x.shape[0]/(K+1),K+1,K + 1))
- temp[:,:-1,:-1] = candidates[i].reshape(-1,K,K)
- temp[:,:,-1] = x.reshape(-1,K+1,1)
- temp = temp.reshape(-1,K+1)
- T_,R_ = sa_fit_candidates(AggOps[i],temp)
- #print "T - T_",(T - T_).data.max()
- #assert((T - T_).data.max() < 1e-10)
- #assert((R - R_).data.max() < 1e-10)
- T,R = T_,R_
+ T,R = augment_candidates(AggOps[i], Ts[i], candidates[i+1], x)
+ #if i == 0:
+ # temp = hstack((candidates[0],x))
+ #else:
+ # K = candidates[i].shape[1]
+ # temp = zeros((x.shape[0]/(K+1),K+1,K + 1))
+ # temp[:,:-1,:-1] = candidates[i].reshape(-1,K,K)
+ # temp[:,:,-1] = x.reshape(-1,K+1,1)
+ # temp = temp.reshape(-1,K+1)
+ #T_,R_ = sa_fit_candidates(AggOps[i],temp)
+ #print "T - T_",abs(T - T_).sum()
+ #assert(abs(T - T_).sum() < 1e-10)
+ #assert((R - R_).sum() < 1e-10)
+ #T,R = T_,R_
#TODO end test
- #T,R = augment_candidates(AggOps[i], Ts[i], candidates[i+1], x, i)
P = smoothed_prolongator(T,A)
A = P.T.tocsr() * A * P
@@ -385,76 +397,83 @@
return new_As,new_Ps,new_Ts,new_Ws
+if __name__ == '__main__':
+ from scipy import *
+ from utils import diag_sparse
+ from multilevel import poisson_problem1D,poisson_problem2D
-from scipy import *
-from utils import diag_sparse
-from multilevel import poisson_problem1D,poisson_problem2D
+ blocks = None
+
+ A = poisson_problem2D(100)
+ #A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
+ #A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
+ #A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
+ #A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
+
+
+ #D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+ #A = D * A * D
+
+ A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+ blocks = arange(A.shape[0]/2).repeat(2)
+
+ asa = adaptive_sa_solver(A,max_candidates=3,mu=10)
+ #scipy.random.seed(0) #make tests repeatable
+ x = randn(A.shape[0])
+ b = A*randn(A.shape[0])
+ #b = zeros(A.shape[0])
+
-blocks = None
-
-A = poisson_problem2D(100)
-#A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
-#A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
-#A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
-#A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
-
-
-#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-#A = D * A * D
-
-A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
-blocks = arange(A.shape[0]/2).repeat(2)
-
-asa = adaptive_sa_solver(A,max_candidates=4,mu=5)
-scipy.random.seed(0) #make tests repeatable
-x = rand(A.shape[0])
-#b = A*rand(A.shape[0])
-b = zeros(A.shape[0])
-
-
-print "solving"
-if True:
- x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=25,tol=1e-7,return_residuals=True)
-else:
- residuals = []
- def add_resid(x):
- residuals.append(linalg.norm(b - A*x))
- A.psolve = asa.solver.psolve
- x_sol = linalg.cg(A,b,x0=x,maxiter=20,tol=1e-12,callback=add_resid)[0]
-
-residuals = array(residuals)/residuals[0]
-
-print "residuals ",residuals
-print "mean convergence factor",(residuals[-1]/residuals[0])**(1.0/len(residuals))
-print "last convergence factor",residuals[-1]/residuals[-2]
-
-print
-print asa.solver
-
-print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
-
-for c in asa.candidates[0].T:
- print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
-
-
-
-##W = asa.AggOps[0]*asa.AggOps[1]
-##pcolor((W * rand(W.shape[1])).reshape((200,200)))
-
-def plot2d_arrows(x):
- from pylab import quiver
- x = x.reshape(-1)
- N = (len(x)/2)**0.5
- assert(2 * N * N == len(x))
- X = linspace(-1,1,N).reshape(1,N).repeat(N,0).reshape(-1)
- Y = linspace(-1,1,N).reshape(1,N).repeat(N,0).T.reshape(-1)
-
- dX = x[0::2]
- dY = x[1::2]
-
- quiver(X,Y,dX,dY)
-
-def plot2d(x):
- from pylab import pcolor
- pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
+ print "solving"
+ if True:
+ x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+ else:
+ residuals = []
+ def add_resid(x):
+ residuals.append(linalg.norm(b - A*x))
+ A.psolve = asa.solver.psolve
+ x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
+ residuals = array(residuals)/residuals[0]
+
+ print "residuals ",residuals
+ print "mean convergence factor",(residuals[-1]/residuals[0])**(1.0/len(residuals))
+ print "last convergence factor",residuals[-1]/residuals[-2]
+
+ print
+ print asa.solver
+
+ print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
+
+
+ def plot2d_arrows(x):
+ from pylab import figure,quiver,show
+ x = x.reshape(-1)
+ N = (len(x)/2)**0.5
+ assert(2 * N * N == len(x))
+ X = linspace(-1,1,N).reshape(1,N).repeat(N,0).reshape(-1)
+ Y = linspace(-1,1,N).reshape(1,N).repeat(N,0).T.reshape(-1)
+
+ dX = x[0::2]
+ dY = x[1::2]
+
+ figure()
+ quiver(X,Y,dX,dY)
+ show()
+
+ def plot2d(x):
+ from pylab import pcolor,figure,show
+ figure()
+ pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
+ show()
+
+ for c in asa.candidates[0].T:
+ plot2d_arrows(c)
+ print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+
+
+
+ ##W = asa.AggOps[0]*asa.AggOps[1]
+ ##pcolor((W * rand(W.shape[1])).reshape((200,200)))
+
+
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-23 04:08:10 UTC (rev 3458)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-24 04:19:51 UTC (rev 3459)
@@ -4,7 +4,7 @@
import scipy
import numpy
-from numpy import ones,zeros,zeros_like,array
+from numpy import ones,zeros,zeros_like,array,asarray
from numpy.linalg import norm
from scipy.linsolve import spsolve
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2007-10-23 04:08:10 UTC (rev 3458)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-10-24 04:19:51 UTC (rev 3459)
@@ -104,7 +104,7 @@
def sa_fit_candidates(AggOp,candidates):
K = candidates.shape[1] # num candidates
-
+
N_fine,N_coarse = AggOp.shape
if K > 1 and candidates.shape[0] == K*N_fine:
@@ -112,27 +112,32 @@
AggOp = expand_into_blocks(AggOp,K,1).tocsr()
N_fine = K*N_fine
- R = zeros((K*N_coarse,K)) #storage for coarse candidates
+ R = zeros((N_coarse,K,K)) #storage for coarse candidates
candidate_matrices = []
for i in range(K):
c = candidates[:,i]
- c = c[diff(AggOp.indptr) == 1] #eliminate DOFs that aggregation misses
+ c = c[diff(AggOp.indptr) == 1] # eliminate DOFs that aggregation misses
+
+ threshold = 1e-10 / abs(c).max() # cutoff for small basis functions
+
X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=AggOp.shape)
#orthogonalize X against previous
for j,A in enumerate(candidate_matrices):
D_AtX = csr_matrix((A.data*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
- R[j::K,i] = D_AtX
+ R[:,j,i] = D_AtX
X.data -= D_AtX[X.indices] * A.data
#normalize X
D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
col_norms = sqrt(D_XtX)
- R[i::K,i] = col_norms
+ mask = col_norms < threshold # set small basis functions to 0
+ col_norms[mask] = 0
+ R[:,i,i] = col_norms
col_norms = 1.0/col_norms
- col_norms[isinf(col_norms)] = 0
+ col_norms[mask] = 0
X.data *= col_norms[X.indices]
candidate_matrices.append(X)
@@ -147,6 +152,8 @@
Q_data[i::K] = X.data
Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
+ R = R.reshape(-1,K)
+
return Q,R
def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2007-10-23 04:08:10 UTC (rev 3458)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2007-10-24 04:19:51 UTC (rev 3459)
@@ -1,17 +1,70 @@
from numpy.testing import *
from scipy.sparse import csr_matrix
-from scipy import arange,ones,zeros,array,eye
+from numpy import arange,ones,zeros,array,eye,vstack,diff
set_package_path()
-pass
+from scipy.sandbox.multigrid.sa import sa_fit_candidates
+from scipy.sandbox.multigrid.adaptive import augment_candidates
restore_path()
+#import pdb; pdb.set_trace()
class TestAdaptiveSA(NumpyTestCase):
def setUp(self):
pass
+
+class TestAugmentCandidates(NumpyTestCase):
+ def setUp(self):
+ self.cases = []
+
+ #two candidates
+
+ ##block candidates
+ ##self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)), vstack((array([1]*9 + [0]*9),arange(2*9))).T ))
+
+
+ def check_first_level(self):
+ cases = []
+
+ ## tests where AggOp includes all DOFs
+ cases.append((csr_matrix((ones(4),array([0,0,1,1]),arange(5)),dims=(4,2)), vstack((ones(4),arange(4))).T ))
+ cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)), vstack((ones(9),arange(9))).T ))
+ cases.append((csr_matrix((ones(9),array([0,0,1,1,2,2,3,3,3]),arange(10)),dims=(9,4)), vstack((ones(9),arange(9))).T ))
+
+ ## tests where AggOp excludes some DOFs
+ cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), vstack((ones(5),arange(5))).T ))
+
+ # overdetermined blocks
+ cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), vstack((ones(5),arange(5),arange(5)**2)).T ))
+ cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)), vstack((ones(9),arange(9),arange(9)**2)).T ))
+ cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)), vstack((ones(9),arange(9))).T ))
+
+ def mask_candidate(AggOp,candidates):
+ #mask out all DOFs that are not included in the aggregation
+ candidates[diff(AggOp.indptr) == 0,:] = 0
+
+ for AggOp,fine_candidates in cases:
+
+ mask_candidate(AggOp,fine_candidates)
+
+ for i in range(1,fine_candidates.shape[1]):
+ Q_expected,R_expected = sa_fit_candidates(AggOp,fine_candidates[:,:i+1])
+
+ old_Q, old_R = sa_fit_candidates(AggOp,fine_candidates[:,:i])
+
+ Q_result,R_result = augment_candidates(AggOp, old_Q, old_R, fine_candidates[:,[i]])
+
+ # compare against SA method (which is assumed to be correct)
+ assert_almost_equal(Q_expected.todense(),Q_result.todense())
+ assert_almost_equal(R_expected,R_result)
+
+ #each fine level candidate should be fit exactly
+ assert_almost_equal(fine_candidates[:,:i+1],Q_result*R_result)
+ assert_almost_equal(Q_result*(Q_result.T*fine_candidates[:,:i+1]),fine_candidates[:,:i+1])
+
+
if __name__ == '__main__':
NumpyTest().run()
Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-23 04:08:10 UTC (rev 3458)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-24 04:19:51 UTC (rev 3459)
@@ -143,6 +143,10 @@
## tests where AggOp excludes some DOFs
self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), ones((5,1)) ))
self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), vstack((ones(5),arange(5))).T ))
+
+ # overdetermined blocks
+ self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), vstack((ones(5),arange(5),arange(5)**2)).T ))
+ self.cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)), vstack((ones(9),arange(9),arange(9)**2)).T ))
self.cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)), vstack((ones(9),arange(9))).T ))
def check_all_cases(self):
@@ -153,7 +157,6 @@
candidates[diff(AggOp.indptr) == 0,:] = 0
for AggOp,fine_candidates in self.cases:
-
mask_candidate(AggOp,fine_candidates)
Q,coarse_candidates = sa_fit_candidates(AggOp,fine_candidates)
More information about the Scipy-svn
mailing list