Changeset 2781

Timestamp:

Apr 28, 2006, 4:54:26 PM (19 years ago)

Author:

duncan

Message:

investigating ticket #8

Location:

inundation/fit_interpolate

Files:

• spike_least_squares.py (modified) (11 diffs)
• spike_test_least_squares.py (modified) (2 diffs)

inundation/fit_interpolate/spike_least_squares.py

@@ -384 +384 @@
                  vertex_coordinates,
                  triangles,
+                 z,
                  point_coordinates = None,
                  alpha = None,
@@ -419 +420 @@
           point coordinates and vertex coordinates are assumed to be
           relative to their respective origins.
+          
+          z: Single 1d vector or array of data at the point_coordinates.
 
         """
@@ -461 +464 @@
             if verbose: print 'Building interpolation matrix'
             self.build_interpolation_matrix_A(point_coordinates,
+                                              z,
                                               verbose = verbose,
                                               expand_search = expand_search,
@@ -489 +493 @@
             if verbose: print 'Building interpolation matrix'
             self.build_interpolation_matrix_A(point_coordinates,
+                                              z,
                                               verbose = verbose,
                                               data_origin = data_origin,
@@ -516 +521 @@
 
     def build_interpolation_matrix_A(self, point_coordinates,
+                                     z,
                                      verbose = False, expand_search = True,
                                      max_points_per_cell=30,
@@ -537 +543 @@
         from utilities.polygon import inside_polygon
         
+        #Convert input to Numeric arrays
+        z = ensure_numeric(z, Float)
+            
         #Convert input to Numeric arrays just in case.
         point_coordinates = ensure_numeric(point_coordinates, Float)
@@ -543 +552 @@
         m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
         n = point_coordinates.shape[0]     #Nbr of data points
-
+        if len(z.shape) > 1:
+            att_num = z.shape[1]
+        else:
+            att_num = 0
+            
+        assert z.shape[0] == point_coordinates.shape[0] 
+        if verbose: print 'Number of attributes: %d' %att_num
         if verbose: print 'Number of datapoints: %d' %n
         if verbose: print 'Number of basis functions: %d' %m
@@ -549 +564 @@
         self.A = Sparse(n,m)
         self.AtA = Sparse(m,m)
+        self.Atz = zeros((m,att_num), Float)
 
         #Build quad tree of vertices
@@ -603 +619 @@
                 for j in js:
                     self.A[i,j] = sigmas[j]
+                    #self.Atz[n-i-1] +=  sigmas[j]*z[n-j-1]
+                    #print "i",i
+                    #print "j",j
+                    #print "self.Atz",self.Atz
+                    #print "z",z 
+                    #a=self.Atz[j]
+                    #b=z[i]
+                    self.Atz[j] +=  sigmas[j]*z[i]
+                    #print "m-i-1",m-i-1
+                    #print "n-j-1",n-j-1
+                    #print "self.A[i,j]",self.A[i,j]
+                    #print "z[n-j-1]",z[n-j-1]
+                    #print "z[i]",z[i]
+                    #print "", 
+                    #print "self.Atz",self.Atz                    
                     for k in js:
                         if interp_only == False:
@@ -688 +719 @@
 
         This is the brute force which is too slow for large problems,
-        but could be used for testing
+but could be used for testing
         """
 
@@ -873 +904 @@
         #FIXME (DSG-DsG): could Sparse_CSR be used here?  Use this format
         # after a matrix is built, before calcs.
+        Atz = self.Atz
+        #print "self.get_A()", self.get_A() 
+        #print "self.Atz",self.Atz 
+        #print 'z', z
         Atz = self.A.trans_mult(z)
+        #print "self.A.trans_mult(z)",self.A.trans_mult(z) 
 
 

inundation/fit_interpolate/spike_test_least_squares.py

@@ -7 +7 @@
 
 
-from least_squares import *
+from spike_least_squares import *
 from Numeric import allclose, array, transpose
+from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
+from utilities.sparse import Sparse, Sparse_CSR
 
 from coordinate_transforms.geo_reference import Geo_reference
@@ -43 +45 @@
         data_coords = [d1, d2, d3]
 
-        interp = Interpolation(points, triangles, data_coords, alpha=0)
         z = [z1, z2, z3]
-
-        A = interp.get_A()
-        print "A",A
+        interp = Interpolation(points, triangles, z,data_coords, alpha=0)
+        #print "interp.get_A()", interp.get_A()
+        A = interp.A
+        #print "A",A
+        #print "z",z 
         Atz = A.trans_mult(z)
-        print "Atz",Atz
+        #print "Atz",Atz
         
         f = interp.fit(z)
         answer = [0, 5., 5.]
 
-        print "f\n",f
-        print "answer\n",answer
+        #print "f\n",f
+        #print "answer\n",answer
 
         assert allclose(f, answer, atol=1e-7)
 
-
+       
+    def test_smooth_attributes_to_mesh_one_point(self):
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0,1.0]
+        z1 = 2
+        z2 = 4
+        z3 = 4
+        data_coords = [d1]
+
+        z = [z1]
+        interp = Interpolation(points, triangles, z,data_coords, alpha=0)
+        #print "interp.get_A()", interp.get_A()
+        A = interp.A
+        #print "A",A
+        #print "z",z 
+        Atz_actual = A.trans_mult(z)
+        #Atz = interp.Atz.todense()
+
+        
+        #print "Atz",Atz
+        #print "Atz_actual",Atz_actual 
+        
+
+        #assert allclose(Atz_actual, Atz, atol=1e-7)
+
+    def ytest_chewin_the_fat(self):
+        
+        A = Sparse(2,2)
+        A[0,0] = 1
+        A[1,0] = 1
+        A[0,1] = 0
+        A[1,1] = 1
+        z = [1,1]
+        m = 2
+        n = 2
+        Atz = zeros((n), Float)
+        r = A.trans_mult(z)
+        for i in range(m):
+            for j in range(n):
+                Atz[i] += A[m-i-1,n-j-1]*z[i]
+        print "A.trans_mult(z)",A.trans_mult(z)
+        print "Atz",Atz 
+        print "A*z", A*z
+        
+
+        #print "answer\n",answer
+
+        assert allclose(f, answer, atol=1e-7)
+
+
+    def test_smooth_att_to_meshII(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 2.0]
+        d3 = [3.0,1.0]
+        data_coords = [d1, d2, d3]
+        z = linear_function(data_coords)
+        #print "z",z 
+
+        interp = Interpolation(points, triangles, z,
+                               data_coords, alpha=0.0)
+        f = interp.fit(z)
+        answer = linear_function(points)
+        #print "f\n",f
+        #print "answer\n",answer
+
+        assert allclose(f, answer)
+
+    def test_smooth_attributes_to_meshIII(self):
+
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0,2.0],
+                        [1.0, 1.0],
+                        [2.0, 1.0],
+                        [0.0,0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2,-0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0,1.0]]
+
+        z = linear_function(point_coords)
+        interp = Interpolation(vertices, triangles, z,
+                               point_coords, alpha=0.0)
+
+        #print 'z',z
+        f = interp.fit(z)
+        answer = linear_function(vertices)
+        #print "f\n",f
+        #print "answer\n",answer
+        assert allclose(f, answer)
+
+    def test_smooth_attributes_to_meshIV(self):
+        """ Testing 2 attributes smoothed to the mesh
+        """
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0, 1.0]
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [d1, d2, d3]
+
+        z = [z1, z2, z3]
+        interp = Interpolation(points, triangles, z,
+                               data_coords, alpha=0.0)
+        f =  interp.fit_points(z)
+        answer = [[0,0], [5., 10.], [5., 10.]]
+        assert allclose(f, answer)
+
+
+    def test_fit_and_interpolation(self):
+        from mesh import Mesh
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Get (enough) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667],
+                       [ 0.0, 1.0],
+                       [ 0.0, 3.0],
+                       [ 1.0, 0.0],
+                       [ 1.0, 1.0],
+                       [ 1.0, 2.0],
+                       [ 1.0, 3.0],
+                       [ 2.0, 1.0],
+                       [ 3.0, 0.0],
+                       [ 3.0, 1.0]]
+
+        z = linear_function(data_points)
+        interp = Interpolation(points, triangles, z,
+                               data_points, alpha=0.0)
+
+        answer = linear_function(points)
+
+        f = interp.fit(z)
+
+        #print "f",f
+        #print "answer",answer
+        assert allclose(f, answer)
+
+        
+    def test_smoothing_and_interpolation(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Get (too few!) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667]]
+
+        z = linear_function(data_points)
+        answer = linear_function(points)
+
+        #Make interpolator with too few data points and no smoothing
+        interp = Interpolation(points, triangles, z,
+                               data_points, alpha=0.0)
+        #Must raise an exception
+        try:
+            f = interp.fit(z)
+        except:
+            pass
+
+        #Now try with smoothing parameter
+        interp = Interpolation(points, triangles, z,
+                               data_points, alpha=1.0e-13)
+
+        f = interp.fit(z)
+        #f will be different from answer due to smoothing
+        assert allclose(f, answer,atol=5)
+
+        #Map back
+        #z1 = interp.interpolate(f)
+        #assert allclose(z, z1)
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Least_Squares,'test')
+    #suite = unittest.makeSuite(Test_Least_Squares,'test_smoothing_and_interpolation')
+    #suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_one_point')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)