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

scipy-svn@scip... scipy-svn@scip...
Wed Oct 3 16:38:00 CDT 2007


Author: wnbell
Date: 2007-10-03 16:37:54 -0500 (Wed, 03 Oct 2007)
New Revision: 3395

Added:
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
Removed:
   trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/coarsen.py
   trunk/scipy/sandbox/multigrid/multilevel.py
Log:
separated SA and RS codes



Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-03 21:37:54 UTC (rev 3395)
@@ -4,36 +4,36 @@
 
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
-from coarsen import sa_constant_interpolation
+from sa import sa_constant_interpolation
 #from utils import infinity_norm
 from utils import approximate_spectral_radius
 
-def fit_candidate(I,x):
-    """
-    For each aggregate in I (i.e. each column of I) compute vector R and 
-    sparse matrix Q (having the sparsity of I) such that the following holds:
+##def fit_candidate(I,x):
+##    """
+##    For each aggregate in I (i.e. each column of I) compute vector R and 
+##    sparse matrix Q (having the sparsity of I) such that the following holds:
+##
+##        Q*R = x     and   Q^T*Q = I
+##
+##    In otherwords, find a prolongator Q with orthonormal columns so that
+##    x is represented exactly on the coarser level by R.
+##    """
+##    x = asarray(x)
+##    Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
+##    R = sqrt(ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0)))  #column 2-norms  
+##
+##    Q.data *= (1.0/R)[Q.indices]  #normalize columns of Q
+##   
+##    #print "norm(R)",scipy.linalg.norm(R)
+##    #print "min(R),max(R)",min(R),max(R)
+##    #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
+##    #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
+##    #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
+##    return Q,R
 
-        Q*R = x     and   Q^T*Q = I
 
-    In otherwords, find a prolongator Q with orthonormal columns so that
-    x is represented exactly on the coarser level by R.
-    """
-    x = asarray(x)
-    Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
-    R = sqrt(ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0)))  #column 2-norms  
 
-    Q.data *= (1.0/R)[Q.indices]  #normalize columns of Q
-   
-    #print "norm(R)",scipy.linalg.norm(R)
-    #print "min(R),max(R)",min(R),max(R)
-    #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
-    #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
-    #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
-    return Q,R
 
-
-
-
 ##def orthonormalize_candidate(I,x,basis):
 ##    Px = csr_matrix((x,I.indices,I.indptr),dims=I.shape,check=False) 
 ##    Rs = []
@@ -62,7 +62,6 @@
     V = concatenate((A.data,B.data))
     return coo_matrix((V,(I,J)),dims=(A.shape[0],A.shape[1]+B.shape[1])).tocsr()
 
-
 def vstack_csr(A,B):
     #OPTIMIZE THIS
     assert(A.shape[1] == B.shape[1])
@@ -145,7 +144,6 @@
     Ps = []
 
     for W in Ws:
-        #P,x = fit_candidate(W,x)
         P,x = fit_candidates(W,x)
         I   = smoothed_prolongator(P,A)  
         A   = I.T.tocsr() * A * I

Modified: trunk/scipy/sandbox/multigrid/coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/coarsen.py	2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/coarsen.py	2007-10-03 21:37:54 UTC (rev 3395)
@@ -1,163 +1,163 @@
-import scipy
-import numpy
-
-from numpy import arange,ones,zeros,sqrt,isinf,asarray,empty
-from scipy.sparse import csr_matrix,isspmatrix_csr
-
-from utils import diag_sparse,approximate_spectral_radius
-import multigridtools
-
-
-def rs_strong_connections(A,theta):
-    """
-    Return a strength of connection matrix using the method of Ruge and Stuben
-
-        An off-diagonal entry A[i.j] is a strong connection iff
-
-                -A[i,j] >= theta * max( -A[i,k] )   where k != i
-    """
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
-    Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
-    return csr_matrix((Sx,Sj,Sp),A.shape)
-
-
-def rs_interpolation(A,theta=0.25):
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-    
-    S = rs_strong_connections(A,theta)
-
-    T = S.T.tocsr()  #transpose S for efficient column access
-
-    Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
-                                               A.indptr,A.indices,A.data,\
-                                               S.indptr,S.indices,S.data,\
-                                               T.indptr,T.indices,T.data)
-
-    return csr_matrix((Ix,Ij,Ip))
-
-
-def sa_strong_connections(A,epsilon):
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
-    Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
-
-    return csr_matrix((Sx,Sj,Sp),A.shape)
-
-def sa_constant_interpolation(A,epsilon,blocks=None):
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-    
-    #handle epsilon = 0 case without creating strength of connection matrix?
-    
-    if blocks is not None:
-        num_dofs   = A.shape[0]
-        num_blocks = blocks.max()
-        
-        if num_dofs != len(blocks):
-            raise ValueError,'improper block specification'
-        
-        # for non-scalar problems, use pre-defined blocks in aggregation
-        # the strength of connection matrix is based on the Frobenius norms of the blocks
-        
-        B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
-        Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B #Frobenius norms of blocks entries of A
-
-        S = sa_strong_connections(Block_Frob,epsilon)
-    
-        Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
-        Pj = Pj[blocks] #expand block aggregates into constituent dofs
-        Pp = B.indptr
-        Px = B.data
-    else:
-        S = sa_strong_connections(A,epsilon)
-        
-        Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
-        Pp = numpy.arange(len(Pj)+1)
-        Px = numpy.ones(len(Pj))
-    
-    return csr_matrix((Px,Pj,Pp))
-
-
-def fit_candidates(AggOp,candidates):
-    K = len(candidates)
-
-    N_fine,N_coarse = AggOp.shape
-
-    if K > 1 and len(candidates[0]) == K*N_fine:
-        #see if fine space has been expanded (all levels except for first)
-        AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
-        N_fine = K*N_fine
-
-    R = zeros((K*N_coarse,K))
-
-    candidate_matrices = []
-    for i,c in enumerate(candidates):
-        X = csr_matrix((c.copy(),AggOp.indices,AggOp.indptr),dims=AggOp.shape)
-       
-        #TODO optimize this  
-
-        #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
-            X.data -= D_AtX[X.indices] * A.data
-
-            #AtX = csr_matrix(A.T.tocsr() * X
-            #R[j::K,i] = AtX.data
-            #X = X - A * AtX 
-    
-        #normalize X
-        XtX = X.T.tocsr() * X
-        col_norms = sqrt(asarray(XtX.sum(axis=0)).flatten())
-        R[i::K,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)
-
-
-    Q_indptr  = K*AggOp.indptr
-    Q_indices = (K*AggOp.indices).repeat(K)
-    for i in range(K):
-        Q_indices[i::K] += i
-    Q_data = empty(N_fine * K)
-    for i,X in enumerate(candidate_matrices):
-        Q_data[i::K] = X.data
-    Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
-
-    coarse_candidates = [R[:,i] for i in range(K)]
-
-    return Q,coarse_candidates
-    
-##    S = sa_strong_connections(A,epsilon)
+##import scipy
+##import numpy
 ##
-##    #tentative (non-smooth) interpolation operator I
-##    Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
-##    Pp = numpy.arange(len(Pj)+1)
-##    Px = numpy.ones(len(Pj))
+##from numpy import arange,ones,zeros,sqrt,isinf,asarray,empty
+##from scipy.sparse import csr_matrix,isspmatrix_csr
+##
+##from utils import diag_sparse,approximate_spectral_radius
+##import multigridtools
+##
+##
+##def rs_strong_connections(A,theta):
+##    """
+##    Return a strength of connection matrix using the method of Ruge and Stuben
+##
+##        An off-diagonal entry A[i.j] is a strong connection iff
+##
+##                -A[i,j] >= theta * max( -A[i,k] )   where k != i
+##    """
+##    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+##    Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
+##    return csr_matrix((Sx,Sj,Sp),A.shape)
+##
+##
+##def rs_interpolation(A,theta=0.25):
+##    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
 ##    
-##    return scipy.sparse.csr_matrix((Px,Pj,Pp))
-
-##def sa_smoother(A,S,omega):
-##    Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+##    S = rs_strong_connections(A,theta)
 ##
-##    return csr_matrix((Bx,Bj,Bp),dims=A.shape)
-    
-def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-   
-    AggOp  = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
-    T,coarse_candidates = fit_candidates(AggOp,candidates)
-
-    D_inv = diag_sparse(1.0/diag_sparse(A))       
-    
-    D_inv_A  = D_inv * A
-    D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
-
-    P = T - (D_inv_A*T)  #same as I=S*P, (faster?)
-           
-    return P,coarse_candidates
-
-
-
+##    T = S.T.tocsr()  #transpose S for efficient column access
+##
+##    Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
+##                                               A.indptr,A.indices,A.data,\
+##                                               S.indptr,S.indices,S.data,\
+##                                               T.indptr,T.indices,T.data)
+##
+##    return csr_matrix((Ix,Ij,Ip))
+##
+##
+##def sa_strong_connections(A,epsilon):
+##    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+##    Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+##
+##    return csr_matrix((Sx,Sj,Sp),A.shape)
+##
+##def sa_constant_interpolation(A,epsilon,blocks=None):
+##    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##    
+##    #handle epsilon = 0 case without creating strength of connection matrix?
+##    
+##    if blocks is not None:
+##        num_dofs   = A.shape[0]
+##        num_blocks = blocks.max()
+##        
+##        if num_dofs != len(blocks):
+##            raise ValueError,'improper block specification'
+##        
+##        # for non-scalar problems, use pre-defined blocks in aggregation
+##        # the strength of connection matrix is based on the Frobenius norms of the blocks
+##        
+##        B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
+##        Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B #Frobenius norms of blocks entries of A
+##
+##        S = sa_strong_connections(Block_Frob,epsilon)
+##    
+##        Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+##        Pj = Pj[blocks] #expand block aggregates into constituent dofs
+##        Pp = B.indptr
+##        Px = B.data
+##    else:
+##        S = sa_strong_connections(A,epsilon)
+##        
+##        Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+##        Pp = numpy.arange(len(Pj)+1)
+##        Px = numpy.ones(len(Pj))
+##    
+##    return csr_matrix((Px,Pj,Pp))
+##
+##
+##def fit_candidates(AggOp,candidates):
+##    K = len(candidates)
+##
+##    N_fine,N_coarse = AggOp.shape
+##
+##    if K > 1 and len(candidates[0]) == K*N_fine:
+##        #see if fine space has been expanded (all levels except for first)
+##        AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
+##        N_fine = K*N_fine
+##
+##    R = zeros((K*N_coarse,K))
+##
+##    candidate_matrices = []
+##    for i,c in enumerate(candidates):
+##        X = csr_matrix((c.copy(),AggOp.indices,AggOp.indptr),dims=AggOp.shape)
+##       
+##        #TODO optimize this  
+##
+##        #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
+##            X.data -= D_AtX[X.indices] * A.data
+##
+##            #AtX = csr_matrix(A.T.tocsr() * X
+##            #R[j::K,i] = AtX.data
+##            #X = X - A * AtX 
+##    
+##        #normalize X
+##        XtX = X.T.tocsr() * X
+##        col_norms = sqrt(asarray(XtX.sum(axis=0)).flatten())
+##        R[i::K,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)
+##
+##
+##    Q_indptr  = K*AggOp.indptr
+##    Q_indices = (K*AggOp.indices).repeat(K)
+##    for i in range(K):
+##        Q_indices[i::K] += i
+##    Q_data = empty(N_fine * K)
+##    for i,X in enumerate(candidate_matrices):
+##        Q_data[i::K] = X.data
+##    Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
+##
+##    coarse_candidates = [R[:,i] for i in range(K)]
+##
+##    return Q,coarse_candidates
+##    
+####    S = sa_strong_connections(A,epsilon)
+####
+####    #tentative (non-smooth) interpolation operator I
+####    Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+####    Pp = numpy.arange(len(Pj)+1)
+####    Px = numpy.ones(len(Pj))
+####    
+####    return scipy.sparse.csr_matrix((Px,Pj,Pp))
+##
+####def sa_smoother(A,S,omega):
+####    Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+####
+####    return csr_matrix((Bx,Bj,Bp),dims=A.shape)
+##    
+##def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
+##    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##   
+##    AggOp  = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
+##    T,coarse_candidates = fit_candidates(AggOp,candidates)
+##
+##    D_inv = diag_sparse(1.0/diag_sparse(A))       
+##    
+##    D_inv_A  = D_inv * A
+##    D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
+##
+##    P = T - (D_inv_A*T)  #same as I=S*P, (faster?)
+##           
+##    return P,coarse_candidates
+##
+##
+##

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-03 21:37:54 UTC (rev 3395)
@@ -2,12 +2,13 @@
            'ruge_stuben_solver','smoothed_aggregation_solver',
            'multilevel_solver']
 
-from numpy.linalg import norm
-from numpy import zeros,zeros_like,array
 import scipy
 import numpy
+from numpy import zeros,zeros_like,array
+from numpy.linalg import norm
 
-from coarsen import sa_interpolation,rs_interpolation
+from sa import sa_interpolation
+from rs import rs_interpolation
 from relaxation import gauss_seidel,jacobi,sor
 from utils import infinity_norm
 
@@ -188,11 +189,12 @@
     #A = io.mmread("9pt-100x100.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()
+    
     #A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
     #candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
     #candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
     
-    ml = smoothed_aggregation_solver(A,candidates,max_coarse=100,max_levels=3)
+    ml = smoothed_aggregation_solver(A,candidates,max_coarse=10,max_levels=10)
     #ml = ruge_stuben_solver(A)
 
     x = rand(A.shape[0])

Deleted: trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py	2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/tests/test_coarsen.py	2007-10-03 21:37:54 UTC (rev 3395)
@@ -1,161 +0,0 @@
-from numpy.testing import *
-
-from numpy import sqrt,empty,ones,arange,array_split
-from scipy import rand
-from scipy.sparse import spdiags,csr_matrix,lil_matrix
-import numpy
-
-set_package_path()
-import scipy.sandbox.multigrid
-from scipy.sandbox.multigrid.coarsen import sa_strong_connections,sa_constant_interpolation
-from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D
-restore_path()
-
-
-def reference_sa_strong_connections(A,epsilon):
-    A_coo = A.tocoo()
-    S = lil_matrix(A.shape)
-
-    for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
-        if i == j: continue #skip diagonal
-
-        if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
-            S[i,j] = v
-
-    return S.tocsr()
-
-
-# note that this method only tests the current implementation, not
-# all possible implementations
-def reference_sa_constant_interpolation(A,epsilon):
-    S = sa_strong_connections(A,epsilon)
-    S = array_split(S.indices,S.indptr[1:-1])
-
-    n = A.shape[0]
-
-    R = set(range(n))
-    j = 0
-
-    aggregates = empty(n,dtype=A.indices.dtype)
-    aggregates[:] = -1
-
-    # Pass #1
-    for i,row in enumerate(S):
-        Ni = set(row) | set([i])
-
-        if Ni.issubset(R):
-            R -= Ni
-            for x in Ni:
-                aggregates[x] = j
-            j += 1
-
-    # Pass #2
-    Old_R = R.copy()
-    for i,row in enumerate(S):
-        if i not in R: continue
-        
-        for x in row:
-            if x not in Old_R:
-                aggregates[i] = aggregates[x]
-                R.remove(i)
-                break
-
-
-    # Pass #3
-    for i,row in enumerate(S):
-        if i not in R: continue
-        Ni = set(row) | set([i])
-
-        for x in Ni:
-            if x in R:
-                aggregates[x] = j
-            j += 1
-
-    assert(len(R) == 0)
-
-    Pj = aggregates
-    Pp = arange(n+1)
-    Px = ones(n)
- 
-    return csr_matrix((Px,Pj,Pp))
-
-class TestSaStrongConnections(NumpyTestCase):
-    def check_simple(self):
-        N = 4
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
-        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
-        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
-        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
-       
-        N = 100
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
-        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
-        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
-        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
-
-    def check_random(self):
-        numpy.random.seed(0)
-
-        for N in [2,3,5,10]:
-            A = csr_matrix(rand(N,N))
-            for epsilon in [0.0,0.1,0.5,0.8,1.0,10.0]:
-                S_result = sa_strong_connections(A,epsilon)
-                S_expected = reference_sa_strong_connections(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-    def check_poisson1D(self):
-        for N in [2,3,5,7,10,11,19]:
-            A = poisson_problem1D(N)
-            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
-                S_result   = sa_strong_connections(A,epsilon)
-                S_expected = reference_sa_strong_connections(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-    def check_poisson2D(self):
-        for N in [2,3,5,7,10,11,19]:
-            A = poisson_problem2D(N)
-            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
-                S_result   = sa_strong_connections(A,epsilon)
-                S_expected = reference_sa_strong_connections(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-
-class TestSaConstantInterpolation(NumpyTestCase):
-    def check_random(self):
-        numpy.random.seed(0)
-        for N in [2,3,5,10]:
-            A = csr_matrix(rand(N,N))
-            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
-                S_result   = sa_constant_interpolation(A,epsilon)
-                S_expected = reference_sa_constant_interpolation(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-    def check_poisson1D(self):
-        for N in [2,3,5,7,10,11,20,21,29,30]:
-            A = poisson_problem1D(N)
-            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
-                S_result   = sa_constant_interpolation(A,epsilon)
-                S_expected = reference_sa_constant_interpolation(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-    def check_poisson2D(self):
-        for N in [2,3,5,7,10,11,20,21,29,30]:
-            A = poisson_problem2D(N)
-            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
-                S_result   = sa_constant_interpolation(A,epsilon)
-                S_expected = reference_sa_constant_interpolation(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
-
-##    def check_sample_data(self):
-##        from examples import all_examples,read_matrix
-##
-##        for filename in all_examples:
-##            A = read_matrix(filename)            
-##            for epsilon in [0.0,0.08,0.51,1.0]:
-##                S_result   = sa_constant_interpolation(A,epsilon)
-##                S_expected = reference_sa_constant_interpolation(A,epsilon)
-##                assert_array_equal((S_result - S_expected).nnz,0)
-
-if __name__ == '__main__':
-    NumpyTest().run()
-

Copied: trunk/scipy/sandbox/multigrid/tests/test_sa.py (from rev 3393, trunk/scipy/sandbox/multigrid/tests/test_coarsen.py)
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py	2007-10-03 06:19:08 UTC (rev 3393)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-03 21:37:54 UTC (rev 3395)
@@ -0,0 +1,160 @@
+from numpy.testing import *
+
+from numpy import sqrt,empty,ones,arange,array_split
+from scipy import rand
+from scipy.sparse import spdiags,csr_matrix,lil_matrix
+import numpy
+
+set_package_path()
+import scipy.sandbox.multigrid
+from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
+                                        sa_interpolation, sa_fit_candidates
+from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D
+restore_path()
+
+
+def reference_sa_strong_connections(A,epsilon):
+    A_coo = A.tocoo()
+    S = lil_matrix(A.shape)
+
+    for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
+        if i == j: continue #skip diagonal
+
+        if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+            S[i,j] = v
+
+    return S.tocsr()
+
+class TestSAStrongConnections(NumpyTestCase):
+    def check_simple(self):
+        N = 4
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
+        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+       
+        N = 100
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
+        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+    def check_random(self):
+        numpy.random.seed(0)
+
+        for N in [2,3,5]:
+            A = csr_matrix(rand(N,N))
+            for epsilon in [0.0,0.1,0.5,1.0,10.0]:
+                S_result = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson1D(self):
+        for N in [2,3,5,7,10,11,19]:
+            A = poisson_problem1D(N)
+            for epsilon in [0.0,0.1,0.5,1.0]:
+                S_result   = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson2D(self):
+        for N in [2,3,5,7,10,11]:
+            A = poisson_problem2D(N)
+            for epsilon in [0.0,0.1,0.5,1.0]:
+                S_result   = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+
+# note that this method only tests the current implementation, not
+# all possible implementations
+def reference_sa_constant_interpolation(A,epsilon):
+    S = sa_strong_connections(A,epsilon)
+    S = array_split(S.indices,S.indptr[1:-1])
+
+    n = A.shape[0]
+
+    R = set(range(n))
+    j = 0
+
+    aggregates = empty(n,dtype=A.indices.dtype)
+    aggregates[:] = -1
+
+    # Pass #1
+    for i,row in enumerate(S):
+        Ni = set(row) | set([i])
+
+        if Ni.issubset(R):
+            R -= Ni
+            for x in Ni:
+                aggregates[x] = j
+            j += 1
+
+    # Pass #2
+    Old_R = R.copy()
+    for i,row in enumerate(S):
+        if i not in R: continue
+        
+        for x in row:
+            if x not in Old_R:
+                aggregates[i] = aggregates[x]
+                R.remove(i)
+                break
+
+    # Pass #3
+    for i,row in enumerate(S):
+        if i not in R: continue
+        Ni = set(row) | set([i])
+
+        for x in Ni:
+            if x in R:
+                aggregates[x] = j
+            j += 1
+
+    assert(len(R) == 0)
+
+    Pj = aggregates
+    Pp = arange(n+1)
+    Px = ones(n)
+ 
+    return csr_matrix((Px,Pj,Pp))
+
+class TestSAConstantInterpolation(NumpyTestCase):
+    def check_random(self):
+        numpy.random.seed(0)
+        for N in [2,3,5,10]:
+            A = csr_matrix(rand(N,N))
+            for epsilon in [0.0,0.1,0.5,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson1D(self):
+        for N in [2,3,5,7,10,11,20,21,29,30]:
+            A = poisson_problem1D(N)
+            for epsilon in [0.0,0.1,0.5,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson2D(self):
+        for N in [2,3,5,7,10,11]:
+            A = poisson_problem2D(N)
+            for epsilon in [0.0,0.1,0.5,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+##    def check_sample_data(self):
+##        from examples import all_examples,read_matrix
+##
+##        for filename in all_examples:
+##            A = read_matrix(filename)            
+##            for epsilon in [0.0,0.08,0.51,1.0]:
+##                S_result   = sa_constant_interpolation(A,epsilon)
+##                S_expected = reference_sa_constant_interpolation(A,epsilon)
+##                assert_array_equal((S_result - S_expected).nnz,0)
+
+if __name__ == '__main__':
+    NumpyTest().run()
+



More information about the Scipy-svn mailing list