Changeset 5775

Timestamp:

Sep 20, 2008, 12:14:46 AM (15 years ago)

Author:

ole

Message:

Rationalised and optimised interpolation

Location:

anuga_core/source/anuga/fit_interpolate

Files:

• interpolate.py (modified) (8 diffs)
• search_functions.py (modified) (2 diffs)
• test_interpolate.py (modified) (13 diffs)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -46 +46 @@
 
 
+def interpolate(vertex_coordinates,
+                triangles,
+                vertex_values,
+                interpolation_points,
+                mesh_origin=None,
+                max_vertices_per_cell=None,
+                start_blocking_len=500000,
+                use_cache=False,             
+                verbose=False):                 
+    """Interpolate vertex_values to interpolation points.
+    
+    Inputs (mandatory):
+
+   
+    vertex_coordinates: List of coordinate pairs [xi, eta] of
+                        points constituting a mesh 
+                        (or an m x 2 Numeric array or
+                        a geospatial object)
+                        Points may appear multiple times
+                        (e.g. if vertices have discontinuities)
+
+    triangles: List of 3-tuples (or a Numeric array) of
+               integers representing indices of all vertices 
+               in the mesh.
+               
+    vertex_values: Vector or array of data at the mesh vertices.
+                   If array, interpolation will be done for each column as
+                   per underlying matrix-matrix multiplication
+              
+    interpolation_points: Interpolate mesh data to these positions.
+                          List of coordinate pairs [x, y] of
+                          data points or an nx2 Numeric array or a Geospatial_data object
+               
+    Inputs (optional)
+               
+    mesh_origin: A geo_reference object or 3-tuples consisting of
+                 UTM zone, easting and northing.
+                 If specified vertex coordinates are assumed to be
+                 relative to their respective origins.
+
+    max_vertices_per_cell: Number of vertices in a quad tree cell
+                           at which the cell is split into 4.
+
+                           Note: Don't supply a vertex coords as a geospatial object and
+                           a mesh origin, since geospatial has its own mesh origin.
+                           
+    start_blocking_len: If the # of points is more or greater than this,
+                        start blocking 
+                        
+    use_cache: True or False
+                            
+
+    Output:
+    
+    Interpolated values at specified point_coordinates                           
+    
+    
+    Note: This function is a simple shortcut for case where interpolation matrix is unnecessary
+    Note: This function does not take blocking into account, but allows caching.
+    
+    """
+    
+    from anuga.caching import cache
+
+    # Create interpolation object with matrix
+    args = (vertex_coordinates, triangles)
+    kwargs = {'mesh_origin': mesh_origin,
+              'max_vertices_per_cell': max_vertices_per_cell,
+              'verbose': verbose}
+              
+    if use_cache is True:
+        I = cache(Interpolate, args, kwargs,
+                  verbose=verbose)
+    else:
+        I = apply(Interpolate, args, kwargs)
+                  
+    
+    # Call interpolate method with interpolation points
+    result = I.interpolate_block(vertex_values, interpolation_points,
+                                 use_cache=use_cache,
+                                 verbose=verbose)
+                           
+    return result
+    
+
 
 class Interpolate (FitInterpolate):
@@ -95 +180 @@
                                 verbose=verbose,
                                 max_vertices_per_cell=max_vertices_per_cell)
+                                
+                                
 
     def interpolate_polyline(self,
@@ -208 +295 @@
         # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
         # method is called even if interpolation points are unchanged.
+        # This really should use some kind of caching in cases where
+        # interpolation points are reused.
 
         #print "point_coordinates interpolate.interpolate", point_coordinates
@@ -270 +359 @@
     
 
-    def interpolate_block(self, f, point_coordinates, verbose=False):
+    def interpolate_block(self, f, point_coordinates, 
+                          use_cache=False, verbose=False):
         """
         Call this if you want to control the blocking or make sure blocking
@@ -287 +377 @@
         f = ensure_numeric(f, Float)        
             
-        self._A = self._build_interpolation_matrix_A(point_coordinates,
-                                                     verbose=verbose)
-
+        if use_cache is True:
+            X = cache(self._build_interpolation_matrix_A,
+                      (point_coordinates),
+                      {'verbose': verbose},                        
+                      verbose=verbose)        
+        else:
+            X = self._build_interpolation_matrix_A(point_coordinates,
+                                                   verbose=verbose)    
+        
+        # Unpack result                                       
+        self._A, self.inside_poly_indices, self.outside_poly_indices = X
+        
 
         # Check that input dimensions are compatible
@@ -353 +452 @@
         
         if verbose: print 'Getting indices inside mesh boundary'
-        self.inside_poly_indices, self.outside_poly_indices  = \
-                     in_and_outside_polygon(point_coordinates,
-                                            self.mesh.get_boundary_polygon(),
-                                            closed = True, verbose = verbose)
-        
-        #Build n x m interpolation matrix
+        inside_poly_indices, outside_poly_indices  =\
+            in_and_outside_polygon(point_coordinates,
+                                   self.mesh.get_boundary_polygon(),
+                                   closed = True, verbose = verbose)
+        
+        # Build n x m interpolation matrix
         if verbose and len(self.outside_poly_indices) > 0:
             print '\n WARNING: Points outside mesh boundary. \n'
+            
         # Since you can block, throw a warning, not an error.
-        if verbose and 0 == len(self.inside_poly_indices):
+        if verbose and 0 == len(inside_poly_indices):
             print '\n WARNING: No points within the mesh! \n'
             
@@ -373 +473 @@
         A = Sparse(n,m)
 
-        n = len(self.inside_poly_indices)
-        #Compute matrix elements for points inside the mesh
+        n = len(inside_poly_indices)
+        
+        # Compute matrix elements for points inside the mesh
         if verbose: print 'Building interpolation matrix from %d points' %n
-        for d, i in enumerate(self.inside_poly_indices):
+        for d, i in enumerate(inside_poly_indices):
             # For each data_coordinate point
             if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
@@ -399 +500 @@
                 msg = 'Could not find triangle for point', x 
                 raise Exception(msg)
-        return A
+        return A, inside_poly_indices, outside_poly_indices
 
 def benchmark_interpolate(vertices,

anuga_core/source/anuga/fit_interpolate/search_functions.py

@@ -55 +55 @@
 
     # Search the last triangle first
-    element_found, sigma0, sigma1, sigma2, k = \
-                   _search_triangles_of_vertices(mesh,last_triangle, x)
+    try:
+        element_found, sigma0, sigma1, sigma2, k = \
+            _search_triangles_of_vertices(mesh, last_triangle, x)
+    except:
+        element_found = False
+                   
     #print "last_triangle", last_triangle
     if element_found is True:
@@ -103 +107 @@
             # been searched.  This is the last try
             element_found, sigma0, sigma1, sigma2, k = \
-                           _search_triangles_of_vertices(mesh,triangles, x)
+                           _search_triangles_of_vertices(mesh, triangles, x)
             is_more_elements = False
         else:
             element_found, sigma0, sigma1, sigma2, k = \
-                       _search_triangles_of_vertices(mesh,triangles, x)
+                       _search_triangles_of_vertices(mesh, triangles, x)
         #search_more_cells_time += time.time()-t0
     #print "search_more_cells_time", search_more_cells_time

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -103 +103 @@
 
         interp = Interpolate(points, vertices)
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
+        A, _, _ = interp._build_interpolation_matrix_A(data)
+        assert allclose(A.todense(),
                         [[1./3, 1./3, 1./3]])
 
@@ -115 +116 @@
         from abstract_2d_finite_volumes.quantity import Quantity
 
-        #Create basic mesh
+        # Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
 
-        #Create shallow water domain
+        # Create shallow water domain
         domain = Domain(points, vertices, boundary)
 
-        #---------------
+        #----------------
         #Constant values
-        #---------------        
+        #----------------        
         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5]])
@@ -142 +143 @@
 
 
-        #---------------
-        #Variable values
-        #---------------
+        #----------------
+        # Variable values
+        #----------------
         quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
                                     [1,4,-9],[2,5,0]])
@@ -162 +163 @@
         
 
+    def test_simple_interpolation_example_using_direct_interface(self):
+        
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        from abstract_2d_finite_volumes.quantity import Quantity
+
+        # Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+
+        #----------------
+        # Constant values
+        #----------------        
+        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
+                                    [4,4,4],[5,5,5]])
+
+
+        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+        # FIXME: This concat should roll into get_vertex_values
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        result = interpolate(vertex_coordinates, triangles, vertex_values, interpolation_points)
+        assert allclose(result, answer)
+
+
+        #----------------
+        # Variable values
+        #----------------
+        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
+                                    [1,4,-9],[2,5,0]])
+        
+        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+        # FIXME: This concat should roll into get_vertex_values
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        result = interpolate(vertex_coordinates, triangles,
+                             vertex_values, interpolation_points)
+        assert allclose(result, answer)        
+        
+        
+    def test_simple_interpolation_example_using_direct_interface_and_caching(self):
+        
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        from abstract_2d_finite_volumes.quantity import Quantity
+
+        # Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+
+        #----------------
+        # First call
+        #----------------
+        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
+                                    [1,4,-9],[2,5,0]])
+        
+        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+        # FIXME: This concat should roll into get_vertex_values
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        result = interpolate(vertex_coordinates, triangles,
+                             vertex_values, interpolation_points,
+                             use_cache=True,
+                             verbose=False)
+        assert allclose(result, answer)                
+        
+        # Second call using the cache
+        result = interpolate(vertex_coordinates, triangles,
+                             vertex_values, interpolation_points,
+                             use_cache=True,
+                             verbose=False)
+        assert allclose(result, answer)                        
+        
+        
     def test_quad_tree(self):
         p0 = [-10.0, -10.0]
@@ -192 +288 @@
                       0., 0. , 0., 0., 0., 0.]]
 
-        
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        assert allclose(A.todense(), answer)
+        
         #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
         #print "PDSG - interp.get_A()", interp.get_A()
@@ -201 +297 @@
                       0., 0. , 0., 0., 0., 0.]]
         
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+        A,_,_ = interp._build_interpolation_matrix_A(data)        
+        assert allclose(A.todense(), answer)
 
 
@@ -211 +307 @@
         answer =  [ [ 0.0,  0.0,  0.0,  0.,
                       0., 0. , 0., 0., 0., 0.]]
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+                      
+        A,_,_ = interp._build_interpolation_matrix_A(data)        
+        assert allclose(A.todense(), answer)
+
 
     def test_datapoints_at_vertices(self):
@@ -230 +328 @@
                    [0., 1., 0.],
                    [0., 0., 1.]]
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+                   
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        assert allclose(A.todense(), answer)
 
 
@@ -251 +350 @@
         interp = Interpolate(points, vertices)
 
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        assert allclose(A.todense(), answer)
 
     def test_datapoints_on_edges(self):
@@ -272 +371 @@
         interp = Interpolate(points, vertices)
 
-        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
-                        answer)
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        assert allclose(A.todense(), answer)
 
 
@@ -292 +391 @@
         interp = Interpolate(points, vertices)
         #print "interp.get_A()", interp.get_A()
-        results = interp._build_interpolation_matrix_A(data).todense()
+        
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        results = A.todense()
         assert allclose(sum(results, axis=1), 1.0)
 
@@ -311 +412 @@
 
         interp = Interpolate(points, vertices)
-        results = interp._build_interpolation_matrix_A(data).todense()
+        
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        results = A.todense()
         assert allclose(sum(results, axis=1), [1,1,1,0])
 
@@ -341 +444 @@
 
 
-        A = interp._build_interpolation_matrix_A(data).todense()
+        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A = A.todense()
         for i in range(A.shape[0]):
             for j in range(A.shape[1]):