Changeset 3900

Timestamp:

Nov 1, 2006, 10:58:50 AM (17 years ago)

Author:

ole

Message:

Added time_thinning to Interpolation_function and improved diagnostics.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/util.py (modified) (9 diffs)
• fit_interpolate/interpolate.py (modified) (5 diffs)
• fit_interpolate/test_interpolate.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -13 +13 @@
 
 def file_function(filename,
-                  domain = None,
-                  quantities = None,
-                  interpolation_points = None,
-                  verbose = False,
-                  use_cache = False):
+                  domain=None,
+                  quantities=None,
+                  interpolation_points=None,
+                  time_thinning=1,
+                  verbose=False,
+                  use_cache=False):
     """Read time history of spatial data from NetCDF file and return
     a callable object.
@@ -58 +59 @@
 
 
+    # Build arguments and keyword arguments for use with caching or apply.
+    args = (filename,)
+    
+    kwargs = {'domain': domain,
+              'quantities': quantities,
+              'interpolation_points': interpolation_points,
+              'time_thinning': time_thinning,                   
+              'verbose': verbose}
+
+
+    # Call underlying engine with or without caching
     if use_cache is True:
         try:
@@ -66 +78 @@
             raise msg
 
-
         f = cache(_file_function,
-                  filename,
-                  {'domain': domain,
-                   'quantities': quantities,
-                   'interpolation_points': interpolation_points,
-                   'verbose': verbose},
-                  dependencies = [filename],
-                  compression = False,
-                  verbose = verbose)
-        #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
-        #structure
+                  args, kwargs,
+                  dependencies=[filename],
+                  compression=False,                  
+                  verbose=verbose)
+
+    else:
+        f = apply(_file_function,
+                  args, kwargs)
+
+
+    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
+    #structure
         
-    else:
-        f = _file_function(filename,
-                           domain,
-                           quantities,
-                           interpolation_points,
-                           verbose)            
 
     return f
@@ -91 +98 @@
 
 def _file_function(filename,
-                   domain = None,
-                   quantities = None,
-                   interpolation_points = None,
-                   verbose = False):
+                   domain=None,
+                   quantities=None,
+                   interpolation_points=None,
+                   time_thinning=1,                                                
+                   verbose=False):
     """Internal function
     
@@ -133 +141 @@
         return get_netcdf_file_function(filename, domain, quantities,
                                         interpolation_points,
+                                        time_thinning=time_thinning,
                                         verbose = verbose)
     else:
@@ -143 +152 @@
                              quantity_names=None,
                              interpolation_points=None,
+                             time_thinning=1,                             
                              verbose = False):
     """Read time history of spatial data from NetCDF sww file and
@@ -238 +248 @@
     time = fid.variables['time'][:]    
 
+    # Get time independent stuff
     if spatial:
         #Get origin
@@ -296 +307 @@
         
     
-    #Produce values for desired data points at
-    #each timestep
-
+    # Produce values for desired data points at
+    # each timestep for each quantity
     quantities = {}
     for i, name in enumerate(quantity_names):
@@ -318 +328 @@
                                   triangles,
                                   interpolation_points,
+                                  time_thinning=time_thinning,
                                   verbose=verbose)
 

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -405 +405 @@
     FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
     geospatial data objects
+
+    Time assumed to be relative to starttime    
     """
   
@@ -415 +417 @@
                  triangles=None,
                  interpolation_points=None,
+                 time_thinning=1,
                  verbose=False):
         """Initialise object and build spatial interpolation if required
+
+        Time_thinning_number controls how many timesteps to use. Only timesteps with
+        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
+        will cause the algorithm to use every third time step.
         """
 
@@ -458 +465 @@
             self.spatial = True            
 
+        # Thin timesteps if needed
+        # Note array() is used to make the thinned arrays contiguous in memory
+        self.time = array(time[::time_thinning])          
+        for name in quantity_names:
+            if len(quantities[name].shape) == 2:
+                quantities[name] = array(quantities[name][::time_thinning,:])
              
         # Save for use with statistics
@@ -468 +481 @@
         self.vertex_coordinates = vertex_coordinates 
         self.interpolation_points = interpolation_points
-        self.time = time[:]  # Time assumed to be relative to starttime
+        
+
         self.index = 0    # Initial time index
         self.precomputed_values = {}
@@ -495 +509 @@
 
             # Build interpolator
-            if verbose: print 'Building interpolation matrix'            
+            if verbose:
+                msg = 'Building interpolation matrix from source mesh '
+                msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
+                                                       triangles.shape[0])
+                print msg
+
+                
             interpol = Interpolate(vertex_coordinates,
                                    triangles,
-                                   #point_coordinates = \
-                                   #self.interpolation_points,
-                                   #alpha = 0,
                                    verbose=verbose)
 
-            if verbose: print 'Interpolate'
+            if verbose:
+                print 'Interpolating (%d interpolation points, %d timesteps).'\
+                      %(self.interpolation_points.shape[0], self.time.shape[0]),
+            
+                if time_thinning > 1:
+                    print 'Timesteps were thinned by a factor of %d' %time_thinning
+                else:
+                    print
+
             for i, t in enumerate(self.time):
                 # Interpolate quantities at this timestep

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -1124 +1124 @@
         else:
             raise 'Should raise exception'
+
+
+
+    def test_interpolation_interface_with_time_thinning(self):
+        # Test spatio-temporal interpolation
+        # Test that spatio temporal function performs the correct
+        # interpolations in both time and space
+    
+        #Three timesteps
+        time = [1.0, 2.0, 4.0, 5.0, 7.0, 8.0, 9.0, 10.0]    
+
+        #Setup mesh used to represent fitted function
+        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
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        #New datapoints where interpolated values are sought
+        interpolation_points = [[ 0.0, 0.0],
+                                [ 0.5, 0.5],
+                                [ 0.7, 0.7],
+                                [ 1.0, 0.5],
+                                [ 2.0, 0.4],
+                                [ 2.8, 1.2]]
+
+        #One quantity
+        Q = zeros( (8,6), Float )
+
+        #Linear in time and space
+        for i, t in enumerate(time):
+            Q[i, :] = t*linear_function(points)
+
+        # Check interpolation of one quantity using interpolaton points) using default
+        # time_thinning of 1
+        I = Interpolation_function(time, Q,
+                                   vertex_coordinates=points,
+                                   triangles=triangles, 
+                                   interpolation_points=interpolation_points,
+                                   verbose=False)
+
+        answer = linear_function(interpolation_points)
+
+        
+        t = time[0]
+        for j in range(50): #t in [1, 6]
+            for id in range(len(interpolation_points)):
+                assert allclose(I(t, id), t*answer[id])
+            t += 0.1    
+
+
+        # Now check time_thinning
+        I = Interpolation_function(time, Q,
+                                   vertex_coordinates=points,
+                                   triangles=triangles, 
+                                   interpolation_points=interpolation_points,
+                                   time_thinning=2,
+                                   verbose=False)
+
+
+        assert len(I.time) == 4
+        assert( allclose(I.time, [1.0, 4.0, 7.0, 9.0] ))    
+
+        answer = linear_function(interpolation_points)
+
+        t = time[0]
+        for j in range(50): #t in [1, 6]
+            for id in range(len(interpolation_points)):
+                assert allclose(I(t, id), t*answer[id])
+            t += 0.1    
+
+