Changeset 5775


Ignore:
Timestamp:
Sep 20, 2008, 12:14:46 AM (14 years ago)
Author:
ole
Message:

Rationalised and optimised interpolation

Location:
anuga_core/source/anuga/fit_interpolate
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5723 r5775  
    4646
    4747
     48def interpolate(vertex_coordinates,
     49                triangles,
     50                vertex_values,
     51                interpolation_points,
     52                mesh_origin=None,
     53                max_vertices_per_cell=None,
     54                start_blocking_len=500000,
     55                use_cache=False,             
     56                verbose=False):                 
     57    """Interpolate vertex_values to interpolation points.
     58   
     59    Inputs (mandatory):
     60
     61   
     62    vertex_coordinates: List of coordinate pairs [xi, eta] of
     63                        points constituting a mesh
     64                        (or an m x 2 Numeric array or
     65                        a geospatial object)
     66                        Points may appear multiple times
     67                        (e.g. if vertices have discontinuities)
     68
     69    triangles: List of 3-tuples (or a Numeric array) of
     70               integers representing indices of all vertices
     71               in the mesh.
     72               
     73    vertex_values: Vector or array of data at the mesh vertices.
     74                   If array, interpolation will be done for each column as
     75                   per underlying matrix-matrix multiplication
     76             
     77    interpolation_points: Interpolate mesh data to these positions.
     78                          List of coordinate pairs [x, y] of
     79                          data points or an nx2 Numeric array or a Geospatial_data object
     80               
     81    Inputs (optional)
     82               
     83    mesh_origin: A geo_reference object or 3-tuples consisting of
     84                 UTM zone, easting and northing.
     85                 If specified vertex coordinates are assumed to be
     86                 relative to their respective origins.
     87
     88    max_vertices_per_cell: Number of vertices in a quad tree cell
     89                           at which the cell is split into 4.
     90
     91                           Note: Don't supply a vertex coords as a geospatial object and
     92                           a mesh origin, since geospatial has its own mesh origin.
     93                           
     94    start_blocking_len: If the # of points is more or greater than this,
     95                        start blocking
     96                       
     97    use_cache: True or False
     98                           
     99
     100    Output:
     101   
     102    Interpolated values at specified point_coordinates                           
     103   
     104   
     105    Note: This function is a simple shortcut for case where interpolation matrix is unnecessary
     106    Note: This function does not take blocking into account, but allows caching.
     107   
     108    """
     109   
     110    from anuga.caching import cache
     111
     112    # Create interpolation object with matrix
     113    args = (vertex_coordinates, triangles)
     114    kwargs = {'mesh_origin': mesh_origin,
     115              'max_vertices_per_cell': max_vertices_per_cell,
     116              'verbose': verbose}
     117             
     118    if use_cache is True:
     119        I = cache(Interpolate, args, kwargs,
     120                  verbose=verbose)
     121    else:
     122        I = apply(Interpolate, args, kwargs)
     123                 
     124   
     125    # Call interpolate method with interpolation points
     126    result = I.interpolate_block(vertex_values, interpolation_points,
     127                                 use_cache=use_cache,
     128                                 verbose=verbose)
     129                           
     130    return result
     131   
     132
    48133
    49134class Interpolate (FitInterpolate):
     
    95180                                verbose=verbose,
    96181                                max_vertices_per_cell=max_vertices_per_cell)
     182                               
     183                               
    97184
    98185    def interpolate_polyline(self,
     
    208295        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
    209296        # method is called even if interpolation points are unchanged.
     297        # This really should use some kind of caching in cases where
     298        # interpolation points are reused.
    210299
    211300        #print "point_coordinates interpolate.interpolate", point_coordinates
     
    270359   
    271360
    272     def interpolate_block(self, f, point_coordinates, verbose=False):
     361    def interpolate_block(self, f, point_coordinates,
     362                          use_cache=False, verbose=False):
    273363        """
    274364        Call this if you want to control the blocking or make sure blocking
     
    287377        f = ensure_numeric(f, Float)       
    288378           
    289         self._A = self._build_interpolation_matrix_A(point_coordinates,
    290                                                      verbose=verbose)
    291 
     379        if use_cache is True:
     380            X = cache(self._build_interpolation_matrix_A,
     381                      (point_coordinates),
     382                      {'verbose': verbose},                       
     383                      verbose=verbose)       
     384        else:
     385            X = self._build_interpolation_matrix_A(point_coordinates,
     386                                                   verbose=verbose)   
     387       
     388        # Unpack result                                       
     389        self._A, self.inside_poly_indices, self.outside_poly_indices = X
     390       
    292391
    293392        # Check that input dimensions are compatible
     
    353452       
    354453        if verbose: print 'Getting indices inside mesh boundary'
    355         self.inside_poly_indices, self.outside_poly_indices  = \
    356                      in_and_outside_polygon(point_coordinates,
    357                                             self.mesh.get_boundary_polygon(),
    358                                             closed = True, verbose = verbose)
    359        
    360         #Build n x m interpolation matrix
     454        inside_poly_indices, outside_poly_indices  =\
     455            in_and_outside_polygon(point_coordinates,
     456                                   self.mesh.get_boundary_polygon(),
     457                                   closed = True, verbose = verbose)
     458       
     459        # Build n x m interpolation matrix
    361460        if verbose and len(self.outside_poly_indices) > 0:
    362461            print '\n WARNING: Points outside mesh boundary. \n'
     462           
    363463        # Since you can block, throw a warning, not an error.
    364         if verbose and 0 == len(self.inside_poly_indices):
     464        if verbose and 0 == len(inside_poly_indices):
    365465            print '\n WARNING: No points within the mesh! \n'
    366466           
     
    373473        A = Sparse(n,m)
    374474
    375         n = len(self.inside_poly_indices)
    376         #Compute matrix elements for points inside the mesh
     475        n = len(inside_poly_indices)
     476       
     477        # Compute matrix elements for points inside the mesh
    377478        if verbose: print 'Building interpolation matrix from %d points' %n
    378         for d, i in enumerate(self.inside_poly_indices):
     479        for d, i in enumerate(inside_poly_indices):
    379480            # For each data_coordinate point
    380481            if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
     
    399500                msg = 'Could not find triangle for point', x
    400501                raise Exception(msg)
    401         return A
     502        return A, inside_poly_indices, outside_poly_indices
    402503
    403504def benchmark_interpolate(vertices,
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r4872 r5775  
    5555
    5656    # Search the last triangle first
    57     element_found, sigma0, sigma1, sigma2, k = \
    58                    _search_triangles_of_vertices(mesh,last_triangle, x)
     57    try:
     58        element_found, sigma0, sigma1, sigma2, k = \
     59            _search_triangles_of_vertices(mesh, last_triangle, x)
     60    except:
     61        element_found = False
     62                   
    5963    #print "last_triangle", last_triangle
    6064    if element_found is True:
     
    103107            # been searched.  This is the last try
    104108            element_found, sigma0, sigma1, sigma2, k = \
    105                            _search_triangles_of_vertices(mesh,triangles, x)
     109                           _search_triangles_of_vertices(mesh, triangles, x)
    106110            is_more_elements = False
    107111        else:
    108112            element_found, sigma0, sigma1, sigma2, k = \
    109                        _search_triangles_of_vertices(mesh,triangles, x)
     113                       _search_triangles_of_vertices(mesh, triangles, x)
    110114        #search_more_cells_time += time.time()-t0
    111115    #print "search_more_cells_time", search_more_cells_time
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r5442 r5775  
    103103
    104104        interp = Interpolate(points, vertices)
    105         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
     105        A, _, _ = interp._build_interpolation_matrix_A(data)
     106        assert allclose(A.todense(),
    106107                        [[1./3, 1./3, 1./3]])
    107108
     
    115116        from abstract_2d_finite_volumes.quantity import Quantity
    116117
    117         #Create basic mesh
     118        # Create basic mesh
    118119        points, vertices, boundary = rectangular(1, 3)
    119120
    120         #Create shallow water domain
     121        # Create shallow water domain
    121122        domain = Domain(points, vertices, boundary)
    122123
    123         #---------------
     124        #----------------
    124125        #Constant values
    125         #---------------       
     126        #----------------       
    126127        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
    127128                                    [4,4,4],[5,5,5]])
     
    142143
    143144
    144         #---------------
    145         #Variable values
    146         #---------------
     145        #----------------
     146        # Variable values
     147        #----------------
    147148        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
    148149                                    [1,4,-9],[2,5,0]])
     
    162163       
    163164
     165    def test_simple_interpolation_example_using_direct_interface(self):
     166       
     167        from mesh_factory import rectangular
     168        from shallow_water import Domain
     169        from Numeric import zeros, Float
     170        from abstract_2d_finite_volumes.quantity import Quantity
     171
     172        # Create basic mesh
     173        points, vertices, boundary = rectangular(1, 3)
     174
     175        # Create shallow water domain
     176        domain = Domain(points, vertices, boundary)
     177
     178        #----------------
     179        # Constant values
     180        #----------------       
     181        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
     182                                    [4,4,4],[5,5,5]])
     183
     184
     185        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
     186        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
     187        # FIXME: This concat should roll into get_vertex_values
     188
     189
     190        # Get interpolated values at centroids
     191        interpolation_points = domain.get_centroid_coordinates()
     192        answer = quantity.get_values(location='centroids')
     193
     194        result = interpolate(vertex_coordinates, triangles, vertex_values, interpolation_points)
     195        assert allclose(result, answer)
     196
     197
     198        #----------------
     199        # Variable values
     200        #----------------
     201        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
     202                                    [1,4,-9],[2,5,0]])
     203       
     204        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
     205        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
     206        # FIXME: This concat should roll into get_vertex_values
     207
     208
     209        # Get interpolated values at centroids
     210        interpolation_points = domain.get_centroid_coordinates()
     211        answer = quantity.get_values(location='centroids')
     212
     213        result = interpolate(vertex_coordinates, triangles,
     214                             vertex_values, interpolation_points)
     215        assert allclose(result, answer)       
     216       
     217       
     218    def test_simple_interpolation_example_using_direct_interface_and_caching(self):
     219       
     220        from mesh_factory import rectangular
     221        from shallow_water import Domain
     222        from Numeric import zeros, Float
     223        from abstract_2d_finite_volumes.quantity import Quantity
     224
     225        # Create basic mesh
     226        points, vertices, boundary = rectangular(1, 3)
     227
     228        # Create shallow water domain
     229        domain = Domain(points, vertices, boundary)
     230
     231        #----------------
     232        # First call
     233        #----------------
     234        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
     235                                    [1,4,-9],[2,5,0]])
     236       
     237        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
     238        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
     239        # FIXME: This concat should roll into get_vertex_values
     240
     241
     242        # Get interpolated values at centroids
     243        interpolation_points = domain.get_centroid_coordinates()
     244        answer = quantity.get_values(location='centroids')
     245
     246        result = interpolate(vertex_coordinates, triangles,
     247                             vertex_values, interpolation_points,
     248                             use_cache=True,
     249                             verbose=False)
     250        assert allclose(result, answer)               
     251       
     252        # Second call using the cache
     253        result = interpolate(vertex_coordinates, triangles,
     254                             vertex_values, interpolation_points,
     255                             use_cache=True,
     256                             verbose=False)
     257        assert allclose(result, answer)                       
     258       
     259       
    164260    def test_quad_tree(self):
    165261        p0 = [-10.0, -10.0]
     
    192288                      0., 0. , 0., 0., 0., 0.]]
    193289
    194        
    195         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    196                         answer)
     290        A,_,_ = interp._build_interpolation_matrix_A(data)
     291        assert allclose(A.todense(), answer)
     292       
    197293        #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
    198294        #print "PDSG - interp.get_A()", interp.get_A()
     
    201297                      0., 0. , 0., 0., 0., 0.]]
    202298       
    203         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    204                         answer)
     299        A,_,_ = interp._build_interpolation_matrix_A(data)       
     300        assert allclose(A.todense(), answer)
    205301
    206302
     
    211307        answer =  [ [ 0.0,  0.0,  0.0,  0.,
    212308                      0., 0. , 0., 0., 0., 0.]]
    213         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    214                         answer)
     309                     
     310        A,_,_ = interp._build_interpolation_matrix_A(data)       
     311        assert allclose(A.todense(), answer)
     312
    215313
    216314    def test_datapoints_at_vertices(self):
     
    230328                   [0., 1., 0.],
    231329                   [0., 0., 1.]]
    232         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    233                         answer)
     330                   
     331        A,_,_ = interp._build_interpolation_matrix_A(data)
     332        assert allclose(A.todense(), answer)
    234333
    235334
     
    251350        interp = Interpolate(points, vertices)
    252351
    253         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    254                         answer)
     352        A,_,_ = interp._build_interpolation_matrix_A(data)
     353        assert allclose(A.todense(), answer)
    255354
    256355    def test_datapoints_on_edges(self):
     
    272371        interp = Interpolate(points, vertices)
    273372
    274         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
    275                         answer)
     373        A,_,_ = interp._build_interpolation_matrix_A(data)
     374        assert allclose(A.todense(), answer)
    276375
    277376
     
    292391        interp = Interpolate(points, vertices)
    293392        #print "interp.get_A()", interp.get_A()
    294         results = interp._build_interpolation_matrix_A(data).todense()
     393       
     394        A,_,_ = interp._build_interpolation_matrix_A(data)
     395        results = A.todense()
    295396        assert allclose(sum(results, axis=1), 1.0)
    296397
     
    311412
    312413        interp = Interpolate(points, vertices)
    313         results = interp._build_interpolation_matrix_A(data).todense()
     414       
     415        A,_,_ = interp._build_interpolation_matrix_A(data)
     416        results = A.todense()
    314417        assert allclose(sum(results, axis=1), [1,1,1,0])
    315418
     
    341444
    342445
    343         A = interp._build_interpolation_matrix_A(data).todense()
     446        A,_,_ = interp._build_interpolation_matrix_A(data)
     447        A = A.todense()
    344448        for i in range(A.shape[0]):
    345449            for j in range(A.shape[1]):
Note: See TracChangeset for help on using the changeset viewer.