Changeset 484

Timestamp:

Nov 4, 2004, 2:13:02 PM (20 years ago)

Author:

ole

Message:

First stab at using quad trees in least_squares.
Unit tests pass and least squares produce results
much faster now.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• cell.py (deleted)
• doit.py (modified) (1 diff)
• general_mesh.py (modified) (2 diffs)
• least_squares.py (modified) (7 diffs)
• mesh.py (modified) (2 diffs)
• quad.py (added)
• test_cell.py (deleted)
• test_least_squares.py (modified) (2 diffs)
• test_quad.py (added)

inundation/ga/storm_surge/pyvolution/doit.py

@@ -8 +8 @@
 
 
-os.system('python compile.py')
+#os.system('python compile.py')
 os.system('python test_all.py')
 os.system('python run_profile.py')

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -122 +122 @@
             #Edgelengths
             self.edgelengths[i, :] = [l0, l1, l2]
+
+        #Build vertex list
+        self.build_vertexlist()        
             
     def __len__(self):
@@ -197 +200 @@
             M = 3*m        #Total number of unique vertices
             return reshape(array(range(M)).astype(Int), (m,3))
+
+    
+    def build_vertexlist(self):
+        """Build vertexlist index by vertex ids and for each entry (point id)
+        build a list of (triangles, vertex_id) pairs that use the point
+        as vertex.
+
+        Preconditions:
+          self.coordinates and self.triangles are defined  
+        
+        Postcondition:
+          self.vertexlist is built
+        """
+
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(self.number_of_elements):
+
+            a = self.triangles[i, 0]
+            b = self.triangles[i, 1]
+            c = self.triangles[i, 2]            
+
+            #Register the vertices v as lists of
+            #(triangle_id, vertex_id) tuples associated with them
+            #This is used for smoothing
+            for vertex_id, v in enumerate([a,b,c]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )    
+
+        self.vertexlist = vertexlist
+    
+

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -19 +19 @@
 
 
-#FIXME: Current implementation uses full matrices and a general solver.
-#Later on we may consider using sparse techniques
-
 import exceptions
 class ShapeError(exceptions.Exception): pass
@@ -28 +25 @@
 from Numeric import zeros, array, Float, Int, dot, transpose
 from LinearAlgebra import solve_linear_equations
-#from scipy import sparse
 from sparse import Sparse
 from cg_solve import conjugate_gradient, VectorShapeError
@@ -200 +196 @@
                 self.B = self.AtA
 
+
         
     def build_interpolation_matrix_A(self, point_coordinates):
         """Build n x m interpolation matrix, where
         n is the number of data points and
-        m is the number of basis functions phi_k (one per vertex)"""
+        m is the number of basis functions phi_k (one per vertex)
+
+        This algorithm uses a quad tree data structure for fast binning of data points
+        """
+
+        from quad import build_quadtree
         
         #Convert input to Numeric arrays
@@ -213 +215 @@
         n = point_coordinates.shape[0]     #Nbr of data points
         
-        #self.A = zeros((n,m), Float)
         self.A = Sparse(n,m)
         self.AtA = Sparse(m,m)
+
+        #Build quad tree of vertices
+        #print 'Building quad tree'
+        root = build_quadtree(self.mesh)
+        #print 'Quad tree built'
+        #root.show()
 
         #Compute matrix elements
@@ -221 +228 @@
             #For each data_coordinate point
 
-            #print 'Doing %d/%d' %(i, n)
+            #print 'Doing %d of %d' %(i, n)
+
+            x = point_coordinates[i]
+
+            #Find vertices near x
+            candidate_vertices = root.search(x[0], x[1])
+
+            #Find triangle containing x
+            element_found = False            
+            for v in candidate_vertices:
+                for triangle_id, vertex_id in self.mesh.vertexlist[v]:
+
+                    k = triangle_id
+                    
+                    #Get the three vertex_points of candidate triangle
+                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
+                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
+                    xi2 = self.mesh.get_vertex_coordinate(k, 2)                 
+
+                    #Get the three normals
+                    n0 = self.mesh.get_normal(k, 0)
+                    n1 = self.mesh.get_normal(k, 1)
+                    n2 = self.mesh.get_normal(k, 2)                
+
+                    #Compute interpolation
+                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
+                    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
+                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
+
+                    #FIXME: Maybe move out to test or something
+                    epsilon = 1.0e-6
+                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
+
+                    #Check that this triangle contains the data point
+                    if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
+                        element_found = True
+                        break
+
+                if element_found is True:
+                    #Don't look for any other triangle
+                    break
+
+            if element_found is True:        
+                #Assign values to matrix A
+
+                j0 = self.mesh.triangles[k,0] #Global vertex id
+                #self.A[i, j0] = sigma0
+
+                j1 = self.mesh.triangles[k,1] #Global vertex id
+                #self.A[i, j1] = sigma1
+
+                j2 = self.mesh.triangles[k,2] #Global vertex id
+                #self.A[i, j2] = sigma2
+
+                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
+                js     = [j0,j1,j2]
+
+                for j in js:
+                    self.A[i,j] = sigmas[j]
+                    for k in js:
+                        self.AtA[j,k] += sigmas[j]*sigmas[k]
+            else:
+                pass
+                #Ok if there is no triangle for datapoint (as in brute force version)
+                #raise 'Could not find triangle for point', x
+
+
+        
+    def build_interpolation_matrix_A_brute(self, point_coordinates):
+        """Build n x m interpolation matrix, where
+        n is the number of data points and
+        m is the number of basis functions phi_k (one per vertex)
+
+
+        This is the brute force which is too slow for large problems, but good for testing
+        """
+
+
+        
+        #Convert input to Numeric arrays
+        point_coordinates = array(point_coordinates).astype(Float)
+        
+        #Build n x m interpolation matrix        
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        n = point_coordinates.shape[0]     #Nbr of data points
+        
+        self.A = Sparse(n,m)
+        self.AtA = Sparse(m,m)
+
+        #Compute matrix elements
+        for i in range(n):
+            #For each data_coordinate point
 
             x = point_coordinates[i]
@@ -254 +352 @@
                     #Assign values to matrix A
 
-
                     j0 = self.mesh.triangles[k,0] #Global vertex id
                     #self.A[i, j0] = sigma0
@@ -273 +370 @@
                 k = k+1
         
-##         self.A   = (self.A).tocsc()
-##         self.AtA = (self.AtA).tocsc()
-##         self.At  =  self.A.transp()
-
 
         

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -120 +120 @@
         #Build neighbour structure    
         self.build_neighbour_structure()
-
-        #Build vertex list
-        self.build_vertexlist()        
 
         #Build surrogate neighbour structure    
@@ -194 +191 @@
                 self.number_of_boundaries[i] -= 1              
 
-
-    
-    def build_vertexlist(self):
-        """Build vertexlist index by vertex ids and for each entry (point id)
-        build a list of (triangles, vertex_id) pairs that use the point
-        as vertex.
-
-        Preconditions:
-          self.coordinates and self.triangles are defined  
-        
-        Postcondition:
-          self.vertexlist is built
-        """
-
-        vertexlist = [None]*len(self.coordinates)
-        for i in range(self.number_of_elements):
-
-            a = self.triangles[i, 0]
-            b = self.triangles[i, 1]
-            c = self.triangles[i, 2]            
-
-            #Register the vertices v as lists of
-            #(triangle_id, vertex_id) tuples associated with them
-            #This is used for smoothing
-            for vertex_id, v in enumerate([a,b,c]):
-                if vertexlist[v] is None:
-                    vertexlist[v] = []
-
-                vertexlist[v].append( (i, vertex_id) )    
-
-        self.vertexlist = vertexlist
-    
 
     def build_surrogate_neighbour_structure(self):

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -647 +647 @@
         fd.close()
         
-        mesh_output_file = "new_trianlge.tsh"
+        mesh_output_file = "new_triangle.tsh"
         fit_to_mesh_file(mesh_file,
                          point_file,
@@ -666 +666 @@
         #clean up
         os.remove(mesh_file)
+        os.remove(mesh_output_file)        
         os.remove(point_file)