Changeset 452

Timestamp:

Oct 26, 2004, 11:55:00 PM (20 years ago)

Author:

steve

Message:

Made some quick and dirty additions to sparse.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (7 diffs)
• sparse.py (modified) (8 diffs)
• test_least_squares.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -27 +27 @@
 from Numeric import zeros, array, Float, Int, dot, transpose
 from LinearAlgebra import solve_linear_equations
-from scipy import sparse
+#from scipy import sparse
+from sparse import Sparse
 from cg_solve import conjugate_gradient, VectorShapeError
 
@@ -210 +211 @@
         
         #self.A = zeros((n,m), Float)
-        self.A = sparse.dok_matrix()
-        self.AtA = sparse.dok_matrix()
+        self.A = Sparse(n,m)
+
+        #print 'n by m ',self.A
+        self.AtA = Sparse(m,m)
 
         #Compute matrix elements
@@ -265 +268 @@
 
         
-        self.A   = (self.A).tocsc()
-        self.AtA = (self.AtA).tocsc()
-        self.At  =  self.A.transp()
+##         self.A   = (self.A).tocsc()
+##         self.AtA = (self.AtA).tocsc()
+##         self.At  =  self.A.transp()
 
     def get_A(self):
@@ -314 +317 @@
 
         #self.D = zeros((m,m), Float)
-        self.D = sparse.dok_matrix()
+        self.D = Sparse(m,m)
 
         #For each triangle compute contributions to D = D1+D2        
@@ -360 +363 @@
             self.D[v0,v2] += e20             
 
-        self.D = (self.D).tocsc()
+        #self.D = (self.D).tocsc()
             
     def fit(self, z):
@@ -383 +386 @@
         
         #Compute right hand side based on data
-        Atz = self.At * z
-
+        Atz = self.A.trans_mult(z)
+
+
+        #print 'fit: Atz',Atz
+        
         #Check sanity
         n, m = self.A.shape
@@ -485 +491 @@
     """
     import os, sys
-    usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh alpha" %         os.path.basename(sys.argv[0])
+    usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh alpha" % os.path.basename(sys.argv[0])
 
     if len(sys.argv) < 4:

inundation/ga/storm_surge/pyvolution/sparse.py

@@ -1 +1 @@
 """Proof of concept sparse matrix code
 """
-
-
+from scipy_base import *
+from cg_solve import conjugate_gradient, VectorShapeError
 
 class Sparse:
@@ -12 +12 @@
         self.M = M
         self.N = N
+        self.shape = (M,N)
         self.A = {}
 
@@ -40 +41 @@
             return 0.0
 
+    def copy(self):
+
+        new = Sparse(self.M,self.N)
+
+        for key in self.A.keys():
+            i, j = key
+
+            new[i,j] = self.A[i,j]
+
+        return new
+
 
     def todense(self):
@@ -49 +61 @@
             for j in range(self.N):
                 if self.A.has_key( (i,j) ):                
-                    D[i, j] = A[ (i,j) ]
+                    D[i, j] = self.A[ (i,j) ]
         return D
 
@@ -67 +79 @@
 
         #Assume numeric types from now on
-        R = zeros(B.shape, Float) #Result
+        R = zeros((self.M,), Float) #Result
         
         if len(B.shape) == 1:
             #Vector
 
+##             print 'B.shape      ',B.shape
+##             print 'self.shape ',self.shape
+            
             assert B.shape[0] == self.N, 'Mismatching dimensions'
 
@@ -78 +93 @@
                 i, j = key
 
-                R[i] += A[key]*B[j]
-
-        else:
-            raise 'Numeric matrix not yet implemented'
+                R[i] += self.A[key]*B[j]
+
+        else:
+            raise ValueError, 'Numeric matrix not yet implemented'
 
         return R
-        
-
+
+    def __add__(self, other):
+        """Add this matrix onto 'other' 
+        """
+
+        from Numeric import array, zeros, Float
+        
+        new = other.copy()
+
+ #       print 'self.shape',self.shape
+ #       print 'other.shape',other.shape
+        
+        for key in self.A.keys():
+            i, j = key
+
+            new[i,j] = new[i,j] + self.A[key]
+
+        return new
+
+
+    def __rmul__(self, other):
+        """Right multiply this matrix with scalar
+        """
+
+        from Numeric import array, zeros, Float
+        
+        if isscalar(other):
+            new = self.copy()
+            #Multiply nonzero elements
+            for key in new.A.keys():
+                i, j = key
+
+                new.A[key] = other*new.A[key]
+        else:
+            raise 'only right multiple with scalar implemented'
+
+#        print 'new.shape',new.shape
+
+        return new
+
+
+    def trans_mult(self, other):
+        """Multiply the transpose of  matrix with 'other' which can be
+        a Numeric vector.
+        """
+
+        from Numeric import array, zeros, Float
+        
+        try:
+            B = array(other)
+        except:
+            print 'FIXME: Only Numeric types implemented so far'
+
+
+        #Assume numeric types from now on
+
+        
+        if len(B.shape) == 1:
+            #Vector
+
+            assert B.shape[0] == self.M, 'Mismatching dimensions'
+
+            R = zeros((self.N,), Float) #Result
+
+            #Multiply nonzero elements
+            for key in self.A.keys():
+                i, j = key
+
+                R[j] += self.A[key]*B[i]
+
+        else:
+            raise 'Can only multiply with 1d array'
+
+        return R
+
+
+    
+    
 if __name__ == '__main__':
 
@@ -122 +213 @@
     assert allclose(u, [6,14,4])
 
+    u = A.trans_mult(v)
+    print u
+    assert allclose(u, [6,6,10])
+
     #Right hand side column
     v = array([[2,4],[3,4],[4,4]])
@@ -130 +225 @@
     #u = A*v[:,1]
     #print u
+    print A.shape
+
+    B = 3*A
+    print B.todense()
+
+    B[1,0] = 2
+
+    C = A+B
+
+    print C.todense()

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -3 +3 @@
 import unittest
 from math import sqrt
+from scipy import mat
 
 from least_squares import *
@@ -121 +122 @@
 
         #Data points
-        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3], [-1.0,-1.5 ], [1.0,-1.0]]
+        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3],
+                        [-1.0/3,-4./3], [-1.0,-1.5 ], [1.0,-1.0]]
         interp = Interpolation(points, triangles, data_points) 
-        
+
+        #FIXME Which matrix does this refer to? AtA?
         answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],   #Affects point d     
                   [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
@@ -131 +134 @@
                   [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f
 
-        assert allclose(interp.get_A(), answer)
+
+        A = mat(interp.get_A())
+        At = transpose(A)
+        AtA = At * A
+        AtA = AtA.asarray()
+        #print AtA
+
+        #FIXME These two matrices are not correct, recalculate
+        assert allclose(AtA, answer)
 
 
@@ -202 +213 @@
                         [0.5, -1.9],
                         [3.0,1.0]]
-        
+
         z = linear_function(point_coords)
         interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
+
+        #print 'z',z
         f = interp.fit(z)
         answer = linear_function(vertices)