Changeset 2659

Timestamp:

Apr 5, 2006, 11:22:10 AM (18 years ago)

Author:

duncan

Message:

Added an extra test, changed interpolate function to handle points outside the mesh

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (2 diffs)
• test_interpolate.py (modified) (7 diffs)

inundation/fit_interpolate/interpolate.py

@@ -265 +265 @@
         """
         z = self._A * f
-        for i in self.outside_poly_indices: # Taking into account points outside the
-                                            # mesh. 
-            z[i,:] = INF
+        # Taking into account points outside the mesh.
+        #print "self.outside_poly_indices", self.outside_poly_indices
+        #print "self.inside_poly_indices", self.inside_poly_indices
+        #print "z", z
+        #print "z.size", z.size
+        for i in self.outside_poly_indices: 
+            z[i] = INF
         return z
 
@@ -511 +515 @@
                     #Use precomputed point
                     Q0 = Q[self.index, point_id]
-                    if ratio > 0: Q1 = Q[self.index+1, point_id]
-
+                    if ratio > 0:
+                        Q1 = Q[self.index+1, point_id]
+            
             #Linear temporal interpolation    
             if ratio > 0:
-                q[i] = Q0 + ratio*(Q1 - Q0)
+                if Q0 == INF and Q1 == INF:
+                    q[i]  = Q0
+                else:
+                    q[i] = Q0 + ratio*(Q1 - Q0)
             else:
                 q[i] = Q0

inundation/fit_interpolate/test_interpolate.py

@@ -727 +727 @@
         # interpolations in both time and space
     
-
         #Three timesteps
-        time = [1.0, 5.0, 6.0]
-        
+        time = [1.0, 5.0, 6.0]    
 
         #Setup mesh used to represent fitted function
@@ -753 +751 @@
                                 [ 2.8, 1.2]]
 
-
         #One quantity
         Q = zeros( (3,6), Float )
@@ -760 +757 @@
         for i, t in enumerate(time):
             Q[i, :] = t*linear_function(points)
-
 
         #Check interpolation of one quantity using interpolaton points)
@@ -769 +765 @@
                                    verbose = False)
 
-
         answer = linear_function(interpolation_points)
         
@@ -776 +771 @@
             for id in range(len(interpolation_points)):
                 assert allclose(I(t, id), t*answer[id])
-
             t += 0.1    
 
@@ -787 +781 @@
 
 
-    def test_points_outside_the_polygon(self):
-        a = [-1.0, 0.0]
-        b = [3.0, 4.0]
-        c = [4.0, 1.0]
-        d = [-3.0, 2.0] #3
-        e = [-1.0, -2.0]
-        f = [1.0, -2.0] #5
-
-        vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
-
-        point_coords = [[-2.0, 2.0],
-                        [-1.0, 1.0],
-                        [0.0, 2.0],
-                        [1.0, 1.0],
-                        [2.0, 1.0],
-                        [0.0, 0.0],
-                        [1.0, 0.0],
-                        [0.0, -1.0],
-                        [-0.2, -0.5],
-                        [-0.9, -1.5],
-                        [0.5, -1.9],
-                        [999999, 9999999]] # point Outside poly
-        geo_data = Geospatial_data(data_points = point_coords)
-
-        interp = Interpolate(vertices, triangles)
-        f = array([linear_function(vertices),2*linear_function(vertices) ])
-        f = transpose(f)
-        #print "f",f
-        z = interp.interpolate(f, geo_data)
-        #z = interp.interpolate(f, point_coords)
-        answer = [linear_function(point_coords),
-                  2*linear_function(point_coords) ]
-        answer = transpose(answer)
-        answer[11,:] = [INF, INF]
-        #print "z",z
-        #print "answer _ fixed",answer
-        assert allclose(z[0:10], answer[0:10])
-        self.failUnless( z[11,1] == answer[11,1], 'Fail!')
-        self.failUnless( z[11,0] == answer[11,0], 'Fail!')
-
-        
+    def test_interpolation_function_outside_point(self):
+        # Test spatio-temporal interpolation
+        # Test that spatio temporal function performs the correct
+        # interpolations in both time and space
+    
+        #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],
+                                [ 545354534, 4354354353]] # outside the mesh
+
+        #One quantity
+        Q = zeros( (3,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)
+        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)-1):
+                assert allclose(I(t, id), t*answer[id])
+            t += 0.1
+            
+        # Now test the point outside the mesh
+        t = time[0]
+        for j in range(50): #t in [1, 6]
+            self.failUnless(I(t, 5) == INF, 'Fail!')
+            t += 0.1  
+            
+        try:    
+            I(1)
+        except:
+            pass
+        else:
+            raise 'Should raise exception'
+
     def test_points_outside_the_polygon(self):
         a = [-1.0, 0.0]
@@ -877 +893 @@
 if __name__ == "__main__":
 
-    #suite = unittest.makeSuite(Test_Interpolate,'test')
     suite = unittest.makeSuite(Test_Interpolate,'test')
     #suite = unittest.makeSuite(Test_Interpolate,'test_points_outside_the_polygon')