# [Scipy-svn] r3459 - in trunk/scipy/sandbox/multigrid: . tests

scipy-svn@scip... scipy-svn@scip...
Tue Oct 23 23:20:15 CDT 2007

Author: wnbell
Date: 2007-10-23 23:19:51 -0500 (Tue, 23 Oct 2007)
New Revision: 3459

Modified:
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
aSA works reasonably well now

===================================================================
--- 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
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)
+
+
+    #D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+    #A = D * A * D
+
+    blocks = arange(A.shape[0]/2).repeat(2)
+
+    #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)
-
-
-#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-#A = D * A * D
-
-blocks = arange(A.shape[0]/2).repeat(2)
-
-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 = []
-        residuals.append(linalg.norm(b - A*x))
-    A.psolve = asa.solver.psolve
-
-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 = []
+            residuals.append(linalg.norm(b - A*x))
+        A.psolve = asa.solver.psolve

+    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
+        R[:,i,i] = col_norms
col_norms = 1.0/col_norms
-        col_norms[isinf(col_norms)] = 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):

===================================================================
--- 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
restore_path()

+#import pdb; pdb.set_trace()

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 ))
+
+            #mask out all DOFs that are not included in the aggregation
+            candidates[diff(AggOp.indptr) == 0,:] = 0
+
+        for AggOp,fine_candidates in cases:
+
+
+            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:
-