Changeset 1904

Timestamp:

Oct 13, 2005, 11:38:46 AM (18 years ago)

Author:

ole

Message:

Modified Interpolation_function to allow for time-independent quantities (such as elevation) and wrote test

Location:

inundation/pyvolution

Files:

• least_squares.py (modified) (5 diffs)
• test_least_squares.py (modified) (2 diffs)

inundation/pyvolution/least_squares.py

@@ -920 +920 @@
     Mandatory input
         time:               px1 array of monotonously increasing times (Float)
-        quantities:         Dictionary of pxm arrays or 1 pxm array (Float)
+        quantities:         Dictionary of arrays or 1 array (Float) 
+                            The arrays must either have dimensions pxm or mx1.
+                            The resulting function will be time dependent in
+                            the former case while it will be constan with
+                            respect to time in the latter case.
         
     Optional input:
@@ -926 +930 @@
         vertex_coordinates: mx2 array of coordinates (Float)
         triangles:          nx3 array of indices into vertex_coordinates (Int)
-        interpolation_points: array of coordinates to be interpolated to 
+        interpolation_points: Nx2 array of coordinates to be interpolated to 
         verbose:            Level of reporting
     
@@ -1002 +1006 @@
         self.vertex_coordinates = vertex_coordinates 
         self.interpolation_points = interpolation_points
-        self.T = time[:]  #Time assumed to be relative to starttime
-        self.index = 0    #Initial time index
+        self.T = time[:]  # Time assumed to be relative to starttime
+        self.index = 0    # Initial time index
         self.precomputed_values = {}
             
@@ -1013 +1017 @@
                 raise 'Triangles and vertex_coordinates must be specified'
             
-
             try:
-                self.interpolation_points =\
-                                          ensure_numeric(self.interpolation_points)
+                self.interpolation_points = ensure_numeric(interpolation_points)
             except:
                 msg = 'Interpolation points must be an N x 2 Numeric array '+\
                       'or a list of points\n'
-                msg += 'I got: %s.' %( str(self.interpolation_points)[:60] + '...')
+                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
+                                      '...')
                 raise msg
 
 
+            m = len(self.interpolation_points)
+            p = len(self.T)
+            
             for name in quantity_names:
-                self.precomputed_values[name] =\
-                                              zeros((len(self.T),
-                                                     len(self.interpolation_points)),
-                                                    Float)
+                self.precomputed_values[name] = zeros((p, m), Float)
 
             #Build interpolator
             interpol = Interpolation(vertex_coordinates,
                                      triangles,
-                                     point_coordinates = self.interpolation_points,
+                                     point_coordinates = \
+                                     self.interpolation_points,
                                      alpha = 0,
                                      precrop = False, 
@@ -1039 +1043 @@
 
             if verbose: print 'Interpolate'
-            n = len(self.T)
             for i, t in enumerate(self.T):
                 #Interpolate quantities at this timestep
-                if verbose and i%((n+10)/10)==0:
-                    print ' time step %d of %d' %(i, n)
+                if verbose and i%((p+10)/10)==0:
+                    print ' time step %d of %d' %(i, p)
                     
                 for name in quantity_names:
-                    self.precomputed_values[name][i, :] =\
-                    interpol.interpolate(quantities[name][i,:])
+                    if len(quantities[name].shape) == 2:
+                        result = interpol.interpolate(quantities[name][i,:])
+                    else:
+                       #Assume no time dependency 
+                       result = interpol.interpolate(quantities[name][:])
+                       
+                    self.precomputed_values[name][i, :] = result
+                   
+                        
 
             #Report

inundation/pyvolution/test_least_squares.py

@@ -1129 +1129 @@
         
 
+    def test_interpolation_function_spatial_only(self):
+        """Test spatio-temporal interpolation with constant time
+        """
+
+        #Three timesteps
+        time = [1.0, 5.0, 6.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 linear in space
+        Q = linear_function(points)
+
+
+        #Check interpolation of one quantity using interpolaton points
+        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), answer[id])
+
+            t += 0.1    
+
+
+        try:    
+            I(1)
+        except:
+            pass
+        else:
+            raise 'Should raise exception'
+            
+
+
     def test_interpolation_function(self):
         """Test spatio-temporal interpolation
@@ -1218 +1279 @@
 
 
+
     def test_fit_and_interpolation_with_different_origins(self):
         """Fit a surface to one set of points. Then interpolate that surface