Changeset 997

Timestamp:

Mar 3, 2005, 6:40:01 PM (20 years ago)

Author:

ole

Message:

Implemented precropping facility in least squares
Also optimised inside_polygon

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (13 diffs)
• mesh.py (modified) (5 diffs)
• test_least_squares.py (modified) (6 diffs)
• util.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -31 +31 @@
 
 from general_mesh import General_mesh
-from Numeric import zeros, array, Float, Int, dot, transpose, concatenate
+from mesh import Mesh
+
+from Numeric import zeros, take, array, Float, Int, dot, transpose, concatenate
 from sparse import Sparse, Sparse_CSR
 from cg_solve import conjugate_gradient, VectorShapeError
@@ -61 +63 @@
                      expand_search = False,
                      data_origin = None,
-                     mesh_origin = None):
+                     mesh_origin = None,
+                     precrop = False):
     """
     Given a mesh file (tsh) and a point attribute file (xya), fit
@@ -111 +114 @@
                     expand_search = expand_search,
                     data_origin = data_origin,
-                    mesh_origin = mesh_origin)
-    if verbose:print "finished fitting to mesh"
+                    mesh_origin = mesh_origin,
+                    precrop = precrop)
+    if verbose: print "finished fitting to mesh"
     
     # convert array to list of lists
@@ -141 +145 @@
                 expand_search = False,
                 data_origin = None,
-                mesh_origin = None):
+                mesh_origin = None,
+                precrop = False):
     """
     Fit a smooth surface to a triangulation,
@@ -175 +180 @@
                            expand_search = expand_search,
                            data_origin = data_origin,
-                           mesh_origin = mesh_origin)
-                           # expand_search = False)
+                           mesh_origin = mesh_origin,
+                           precrop = precrop)
     
     vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
@@ -248 +253 @@
                  max_points_per_cell = 30,
                  mesh_origin = None, 
-                 data_origin = None):
+                 data_origin = None,
+                 precrop = False):
 
         
@@ -282 +288 @@
         #Build underlying mesh
         if verbose: print 'Building mesh' 
-        self.mesh = General_mesh(vertex_coordinates, triangles,
-                                 origin = mesh_origin)
+        #self.mesh = General_mesh(vertex_coordinates, triangles,
+        #FIXME: Trying the normal mesh while testing precrop
+        self.mesh = Mesh(vertex_coordinates, triangles,
+                         origin = mesh_origin)
         self.data_origin = data_origin
+
+        self.point_indices = None
 
         #Smoothing parameter
@@ -295 +305 @@
                                         max_points_per_cell =\
                                         max_points_per_cell,
-                                        data_origin = data_origin)
+                                        data_origin = data_origin,
+                                        precrop = precrop)
         
 
@@ -308 +319 @@
                                    verbose = False, expand_search = True,
                                    max_points_per_cell=30,
-                                   data_origin = None):
+                                   data_origin = None,
+                                   precrop = False):
         """Build final coefficient matrix"""
         
@@ -324 +336 @@
                                               max_points_per_cell =\
                                               max_points_per_cell,
-                                              data_origin = data_origin)
+                                              data_origin = data_origin,
+                                              precrop = precrop)
 
             if self.alpha <> 0:
@@ -337 +350 @@
                                      verbose = False, expand_search = True,
                                      max_points_per_cell=30,
-                                     data_origin = None):
+                                     data_origin = None,
+                                     precrop = False):
         """Build n x m interpolation matrix, where
         n is the number of data points and
@@ -396 +410 @@
 
         #Remove points falling outside mesh boundary
-        #from util import inside_polygon
-        #P = self.mesh.get_boundary_polygon()
-        #FIXME: TODO
+        #This reduced one example from 1356 seconds to 825 seconds
+        #And more could be had by writing util.inside_polygon in C
+        if precrop is True:
+            from Numeric import take
+            from util import inside_polygon
+
+            if verbose: print 'Getting boundary polygon'
+            P = self.mesh.get_boundary_polygon()
+
+            if verbose: print 'Getting indices inside mesh boundary'           
+            indices = inside_polygon(point_coordinates, P, verbose = verbose)
+
+            if verbose:
+                if len(indices) != point_coordinates.shape[0]:
+                    print '%d points outside mesh have been cropped.'\
+                          %(point_coordinates.shape[0] - len(indices))
+            point_coordinates = take(point_coordinates, indices)
+            self.point_indices = indices
 
 
@@ -717 +746 @@
         z = array(z).astype(Float)
 
-
         if len(z.shape) > 1 :
             raise VectorShapeError, 'Can only deal with 1d data vector'
+
+        if self.point_indices is not None:
+            #Remove values for any points that were outside mesh
+            z = take(z, self.point_indices) 
         
         #Compute right hand side based on data

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -58 +58 @@
     
     def __init__(self, coordinates, triangles, boundary = None,
-                 tagged_elements = None):
+                 tagged_elements = None, origin = None):
         """
         Build triangles from x,y coordinates (sequence of 2-tuples or
@@ -67 +67 @@
         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
 
-        General_mesh.__init__(self, coordinates, triangles)
+        General_mesh.__init__(self, coordinates, triangles, origin)
 
         N = self.number_of_elements
@@ -350 +350 @@
         """Return bounding polygon as a list of points
         """
-        from Numeric import allclose
+        from Numeric import allclose, sqrt, array, minimum, maximum
+    
         
 
@@ -356 +357 @@
         segments = {}
 
-        pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
+        #pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
+        #pmax = (max(self.coordinates[:,0]), max(self.coordinates[:,1]))
+
+        #FIXME:Can this be written more compactly, e.g.
+        #using minimum and maximium?
+        pmin = array( [min(self.coordinates[:,0]),
+                       min(self.coordinates[:,1]) ] )
+
+        pmax = array( [max(self.coordinates[:,0]),
+                       max(self.coordinates[:,1]) ] )        
+        
+        mindist = sqrt(sum( (pmax-pmin)**2 ))
         for i, edge_id in self.boundary.keys():
             #Find vertex ids for boundary segment
@@ -366 +378 @@
             B = tuple(self.get_vertex_coordinate(i, b))
 
+            #Take the point closest to pmin as starting point
+            #Note: Could be arbitrary, but nice to have
+            #a unique way of selecting
+            dist_A = sqrt(sum( (A-pmin)**2 ))
+            dist_B = sqrt(sum( (B-pmin)**2 ))                         
+
             #Find minimal point
-            if allclose(A, pmin): p0 = A
-            if allclose(B, pmin): p0 = B                
+            if dist_A < mindist:
+                mindist = dist_A
+                p0 = A
+            if dist_B < mindist:
+                mindist = dist_B
+                p0 = B
+
+
+            if p0 is None:
+                raise 'Weird'
+                p0 = A #We need a starting point (FIXME)
+
+                print 'A', A
+                print 'B', B
+                print 'pmin', pmin
+                print
 
             segments[A] = B

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -301 +301 @@
         assert allclose(sum(interp.get_A(), axis=1), 1.0)
         
-
-    # this causes a memory error in scipy.sparce       
+    def test_arbitrary_datapoints_some_outside(self):
+        """Try arbitrary datapoints one outside the triangle.
+        That one should be ignored
+        """
+
+        from Numeric import sum
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        points = [a, b, c]
+        vertices = [ [1,0,2] ]   #bac
+
+        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
+
+
+        interp = Interpolation(points, vertices, data, precrop = True)
+        assert allclose(sum(interp.get_A(), axis=1), 1.0)
+
+        interp = Interpolation(points, vertices, data, precrop = False)
+        assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
+        
+        
+
+    # this causes a memory error in scipy.sparse       
     def test_more_triangles(self):
         
@@ -844 +867 @@
                         [ 1.0, 0.0],
                         [ 1.0, 1.0],
+                        [ 15, -17],   #Outside mesh
                         [ 1.0, 2.0],
                         [ 1.0, 3.0],                                          
@@ -851 +875 @@
 
         #Fit surface to mesh
-        interp = Interpolation(points, triangles, data_points1, alpha=0.0)
+        interp = Interpolation(points, triangles, data_points1, alpha=0.0,
+                               precrop = True)
         z = linear_function(data_points1) #Example z-values
         f = interp.fit(z)                 #Fitted values at vertices
@@ -861 +886 @@
                         [ 0.5, 0.5],
                         [ 0.7, 0.7],
+                        [-13, 65],  #Outside
                         [ 1.0, 0.5],
                         [ 2.0, 0.4],
@@ -866 +892 @@
         
 
-        #Build new A matrix based on new points
-        interp.build_interpolation_matrix_A(data_points2)
+
+        #Build new A matrix based on new points (without precrop)
+        interp.build_interpolation_matrix_A(data_points2, precrop = False)
 
         #Interpolate using fitted surface
         z1 = interp.interpolate(f)
+
+        #import Numeric
+        #data_points2 = Numeric.take(data_points2, interp.point_indices)
+
+        #Desired result (OK for points inside)
+
+        answer = linear_function(data_points2)
+        import Numeric
+        z1 = Numeric.take(z1, [0,1,2,4,5,6])
+        answer = Numeric.take(answer, [0,1,2,4,5,6])                
+        assert allclose(z1, answer)
+
+        #Build new A matrix based on new points (with precrop)
+        interp.build_interpolation_matrix_A(data_points2, precrop = True)
+
+        #Interpolate using fitted surface
+        z1 = interp.interpolate(f)
+
+        import Numeric
+        data_points2 = Numeric.take(data_points2, interp.point_indices)
 
         #Desired result
@@ -1112 +1159 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(TestCase,'test')
-
     suite = unittest.makeSuite(TestCase,'test')
+
     #suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
     runner = unittest.TextTestRunner(verbosity=1)

inundation/ga/storm_surge/pyvolution/util.py

@@ -901 +901 @@
 #FIXME: ALl these should be put into new module polygon.py
 
-def inside_polygon(point, polygon, closed = True):
+def inside_polygon(point, polygon, closed = True, verbose = False):
     """Determine whether points are inside or outside a polygon
     
@@ -969 +969 @@
 
     N = polygon.shape[0] #Number of vertices in polygon 
-
+    M = points.shape[0]  #Number of points
+    
     px = polygon[:,0]
     py = polygon[:,1]   
 
-    #Begin algorithm
+    #Used for an optimisation when points are far away from polgon
+    minpx = min(px); maxpx = max(px)
+    minpy = min(py); maxpy = max(py)    
+
+
+    #Begin main loop (FIXME: It'd be crunchy to have this written in C:-)
     indices = []
-    for k in range(points.shape[0]):
+    for k in range(M):
+
+        if verbose:
+            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
+        
         #for each point   
         x = points[k, 0]
@@ -981 +991 @@
 
         inside = False
+
+        #Optimisation
+        if x > maxpx or x < minpx: continue
+        if y > maxpy or y < minpy: continue        
+
+        #Check polygon 
         for i in range(N):
             j = (i+1)%N