Changeset 2189

Timestamp:

Jan 9, 2006, 4:09:06 PM (19 years ago)

Author:

duncan

Message:

refactoring, has old code commented out

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (8 diffs)
• test_interpolate.py (modified) (1 diff)

inundation/fit_interpolate/interpolate.py

@@ -33 +33 @@
                  vertex_coordinates,
                  triangles,
-                 point_coordinates = None,
+                 #point_coordinates = None,
                  mesh_origin = None,
                  verbose=False,
@@ -68 +68 @@
         from pyvolution.util import ensure_numeric
 
+        # Initialise variabels
+        self.A_can_be_reused = False
+        self.point_coordinates = None
+        
         #Convert input to Numeric arrays
         triangles = ensure_numeric(triangles, Int)
@@ -143 +147 @@
 
             if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
+
             x = point_coordinates[i]
-
-            #Find vertices near x
-            candidate_vertices = self.root.search(x[0], x[1])
-            is_more_elements = True
-
+            
             element_found, sigma0, sigma1, sigma2, k = \
-                self.search_triangles_of_vertices(candidate_vertices, x)
-            first_expansion = True
-            while not element_found and is_more_elements:
-                #if verbose: print 'Expanding search'
-                if first_expansion == True:
-                    #self.expanded_quad_searches.append(1)
-                    first_expansion = False
-                #else:
-                    #end = len(self.expanded_quad_searches) - 1
-                    #assert end >= 0
-                    #self.expanded_quad_searches[end] += 1
-                candidate_vertices, branch = self.root.expand_search()
-                if branch == []:
-                    # Searching all the verts from the root cell that haven't
-                    # been searched.  This is the last try
-                    element_found, sigma0, sigma1, sigma2, k = \
-                      self.search_triangles_of_vertices(candidate_vertices, x)
-                    is_more_elements = False
-                else:
-                    element_found, sigma0, sigma1, sigma2, k = \
-                      self.search_triangles_of_vertices(candidate_vertices, x)
-
-                
+                           self._search_tree_of_vertices(x)    
             #Update interpolation matrix A if necessary
             if element_found is True:
@@ -186 +165 @@
                 for j in js:
                     A[i,j] = sigmas[j]
-                    
             else:
-                print 'Could not find triangle for point', x
-
-
-
+                print 'Could not find triangle for point', x 
+                
         return A
 
-
-    def search_triangles_of_vertices(self,
+    def _search_tree_of_vertices(self, x):
+        """
+        Find the triangle (element) that the point x is in.
+
+        Return the associated sigma and k values
+           (and if the element was found) .
+        """
+        #Find triangle containing x:
+        element_found = False
+
+        # This will be returned if element_found = False
+        sigma2 = -10.0
+        sigma0 = -10.0
+        sigma1 = -10.0
+        k = -10.0
+            
+        #Find vertices near x
+        candidate_vertices = self.root.search(x[0], x[1])
+        is_more_elements = True
+
+        element_found, sigma0, sigma1, sigma2, k = \
+                       self._search_triangles_of_vertices(candidate_vertices, x)
+        #FIXME Do we need this var?
+        # Exclude points outside the mesh before removing this. 
+        #first_expansion = True
+        while not element_found and is_more_elements:
+            #if verbose: print 'Expanding search'
+            #if first_expansion == True:
+            #self.expanded_quad_searches.append(1)
+            #first_expansion = False
+            #else:
+            #end = len(self.expanded_quad_searches) - 1
+            #assert end >= 0
+            #self.expanded_quad_searches[end] += 1
+            candidate_vertices, branch = self.root.expand_search()
+            if branch == []:
+                # Searching all the verts from the root cell that haven't
+                # been searched.  This is the last try
+                element_found, sigma0, sigma1, sigma2, k = \
+                      self._search_triangles_of_vertices(candidate_vertices, x)
+                is_more_elements = False
+            else:
+                element_found, sigma0, sigma1, sigma2, k = \
+                      self._search_triangles_of_vertices(candidate_vertices, x)
+
+        return element_found, sigma0, sigma1, sigma2, k
+
+    def _search_triangles_of_vertices(self,
                                  candidate_vertices, x):
             #Find triangle containing x:
@@ -215 +237 @@
                 #for each triangle id (k) which has v as a vertex
                 for k, _ in self.mesh.vertexlist[v]:
-
                     #Get the three vertex_points of candidate triangle
                     xi0 = self.mesh.get_vertex_coordinate(k, 0)
@@ -221 +242 @@
                     xi2 = self.mesh.get_vertex_coordinate(k, 2)
 
-                    #print "PDSG - k", k
-                    #print "PDSG - xi0", xi0
-                    #print "PDSG - xi1", xi1
-                    #print "PDSG - xi2", xi2
-                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
-                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
-
                     #Get the three normals
                     n0 = self.mesh.get_normal(k, 0)
@@ -233 +247 @@
                     n2 = self.mesh.get_normal(k, 2)
 
-
                     #Compute interpolation
                     sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
                     sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
                     sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
-
-                    #print "PDSG - sigma0", sigma0
-                    #print "PDSG - sigma1", sigma1
-                    #print "PDSG - sigma2", sigma2
 
                     #FIXME: Maybe move out to test or something
@@ -264 +273 @@
 
 
-    #def get_A(self):
-     #       return self.A.todense()
-    
-    def interpolate(self, f):
-        """Interpolate mesh data f to points.
+     # FIXME: What is a good start_blocking_count value?
+    def interpolate(self, f, point_coordinates = None,
+                    start_blocking_len = 500000, verbose=False):
+        """Interpolate mesh data f to determine values, z, at points.
 
         f is the data on the mesh vertices.
 
-
         The mesh values representing a smooth surface are
-        assumed to be specified in f. This argument could,
-        for example have been obtained from the method self.fit()
-
-        Pre Condition:
-          self.A has been initialised
+        assumed to be specified in f.
 
         Inputs:
           f: Vector or array of data at the mesh vertices.
-          If f is an array, interpolation will be done for each column as
-          per underlying matrix-matrix multiplication
+              If f is an array, interpolation will be done for each column as
+              per underlying matrix-matrix multiplication
+          point_coordinates: Interpolate mesh data to these positions.
+          start_blocking_len: If the # of points is more or greater than this,
+              start blocking 
 
         Output:
-          Interpolated values at data points implied in self.A
-
-        """
-
+          Interpolated values at inputted points (z).
+        """
+        
+        # Can I interpolate, based on previous point_coordinates? 
+        if point_coordinates is None:
+            if self.A_can_be_reused is True and \
+            len(self.point_coordinates) < start_blocking_len: 
+                z = self._get_point_data_z(f,
+                                           verbose=verbose)
+            elif self.point_coordinates is not None:
+                #     if verbose, give warning
+                if verbose:
+                    print 'WARNING: Recalculating A matrix, due to blocking.'
+                point_coordinates = self.point_coordinates
+            else:
+                #There are no good point_coordinates. import sys; sys.exit()
+                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
+                raise msg
+            
+        if point_coordinates is not None:
+            self.point_coordinates = point_coordinates
+            if len(point_coordinates) < start_blocking_len or \
+                   start_blocking_len == 0:
+                self.A_can_be_reused = True
+                z = self.interpolate_block(f, point_coordinates,
+                                           verbose=verbose)
+            else:
+                #Handle blocking
+                self.A_can_be_reused = False
+                start=0
+                z = self.interpolate_block(f, point_coordinates[0:0])
+                for end in range(start_blocking_len
+                                 ,len(point_coordinates)
+                                 ,start_blocking_len):
+                    t = self.interpolate_block(f, point_coordinates[start:end],
+                                               verbose=verbose)
+                    z = concatenate((z,t))
+                    start = end
+                end = len(point_coordinates)
+                t = self.interpolate_block(f, point_coordinates[start:end],
+                                           verbose=verbose)
+                z = concatenate((z,t))
+        return z
+
+    def interpolate_block(self, f, point_coordinates = None, verbose=False):
+        """
+        Call this if you want to control the blocking or make sure blocking
+        doesn't occur.
+
+        See interpolate for doc info.
+        """
+        if point_coordinates is not None:
+            self.A =self._build_interpolation_matrix_A(point_coordinates,
+                                                       verbose=verbose)
+        return self._get_point_data_z(f)
+
+    def _get_point_data_z(self, f, verbose=False):
         return self.A * f
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/fit_interpolate/test_interpolate.py

@@ -239 +239 @@
 
 
-
+    def test_interpolate_attributes_to_points(self):
+        v0 = [0.0, 0.0]
+        v1 = [0.0, 5.0]
+        v2 = [5.0, 0.0]
+
+        vertices = [v0, v1, v2]
+        triangles = [ [1,0,2] ]   #bac
+
+        d0 = [1.0, 1.0]
+        d1 = [1.0, 2.0]
+        d2 = [3.0, 1.0]
+        point_coords = [ d0, d1, d2]
+
+        interp = Interpolate(vertices, triangles, point_coords)
+        f = linear_function(vertices)
+        z = interp.interpolate(f, point_coords)
+        answer = linear_function(point_coords)
+
+
+        assert allclose(z, answer)
+
+
+
+    def test_interpolate_attributes_to_pointsII(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],
+                        [3.0, 1.0]]
+
+        interp = Interpolate(vertices, triangles)
+        f = linear_function(vertices)
+        z = interp.interpolate(f, point_coords)
+        answer = linear_function(point_coords)
+        #print "z",z
+        #print "answer",answer
+        assert allclose(z, answer)
+
+    def test_interpolate_attributes_to_pointsIII(self):
+        """Test linear interpolation of known values at vertices to
+        new points inside a triangle
+        """
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        d = [5.0, 5.0]
+
+        vertices = [a, b, c, d]
+        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
+
+        #Points within triangle 1
+        d0 = [1.0, 1.0]
+        d1 = [1.0, 2.0]
+        d2 = [3.0, 1.0]
+
+        #Point within triangle 2
+        d3 = [4.0, 3.0]
+
+        #Points on common edge
+        d4 = [2.5, 2.5]
+        d5 = [4.0, 1.0]
+
+        #Point on common vertex
+        d6 = [0., 5.]
+        
+        point_coords = [d0, d1, d2, d3, d4, d5, d6]
+
+        interp = Interpolate(vertices, triangles)
+
+        #Known values at vertices
+        #Functions are x+y, x+2y, 2x+y, x-y-5
+        f = [ [0., 0., 0., -5.],        # (0,0)
+              [5., 10., 5., -10.],      # (0,5)
+              [5., 5., 10.0, 0.],       # (5,0)
+              [10., 15., 15., -5.]]     # (5,5)
+
+        z = interp.interpolate(f, point_coords)
+        answer = [ [2., 3., 3., -5.],   # (1,1)
+                   [3., 5., 4., -6.],   # (1,2)
+                   [4., 5., 7., -3.],   # (3,1)
+                   [7., 10., 11., -4.], # (4,3)
+                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
+                   [5., 6., 9., -2.],   # (4,1)
+                   [5., 10., 5., -10.]]  # (0,5)
+
+        #print "***********"
+        #print "z",z
+        #print "answer",answer
+        #print "***********"
+
+        #Should an error message be returned if points are outside
+        # of the mesh? Not currently.  
+
+        assert allclose(z, answer)
+
+
+    def test_interpolate_point_outside_of_mesh(self):
+        """Test linear interpolation of known values at vertices to
+        new points inside a triangle
+        """
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        d = [5.0, 5.0]
+
+        vertices = [a, b, c, d]
+        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
+
+        #Far away point
+        d7 = [-1., -1.]
+        
+        point_coords = [ d7]
+        interp = Interpolate(vertices, triangles)
+
+        #Known values at vertices
+        #Functions are x+y, x+2y, 2x+y, x-y-5
+        f = [ [0., 0., 0., -5.],        # (0,0)
+              [5., 10., 5., -10.],      # (0,5)
+              [5., 5., 10.0, 0.],       # (5,0)
+              [10., 15., 15., -5.]]     # (5,5)
+
+        z = interp.interpolate(f, point_coords)
+        answer = [ [0., 0., 0., 0.]] # (-1,-1)
+
+        #print "***********"
+        #print "z",z
+        #print "answer",answer
+        #print "***********"
+
+        #Should an error message be returned if points are outside
+        # of the mesh? Not currently.  
+
+        assert allclose(z, answer)
+
+    def test_interpolate_attributes_to_pointsIV(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],
+                        [3.0, 1.0]]
+
+        interp = Interpolate(vertices, triangles)
+        f = array([linear_function(vertices),2*linear_function(vertices) ])
+        f = transpose(f)
+        #print "f",f
+        z = interp.interpolate(f, point_coords)
+        answer = [linear_function(point_coords),
+                  2*linear_function(point_coords) ]
+        answer = transpose(answer)
+        #print "z",z
+        #print "answer",answer
+        assert allclose(z, answer)
+
+
+    def test_interpolate_blocking(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],
+                        [3.0, 1.0]]
+
+        interp = Interpolate(vertices, triangles)
+        f = array([linear_function(vertices),2*linear_function(vertices) ])
+        f = transpose(f)
+        #print "f",f
+        for blocking_max in range(len(point_coords)+2):
+        #if True:
+         #   blocking_max = 5
+            z = interp.interpolate(f, point_coords,
+                                   start_blocking_len=blocking_max)
+            answer = [linear_function(point_coords),
+                      2*linear_function(point_coords) ]
+            answer = transpose(answer)
+            #print "z",z
+            #print "answer",answer
+            assert allclose(z, answer)
+
+    def test_interpolate_reuse(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],
+                        [3.0, 1.0]]
+
+        interp = Interpolate(vertices, triangles)
+        f = array([linear_function(vertices),2*linear_function(vertices) ])
+        f = transpose(f)
+        z = interp.interpolate(f, point_coords,
+                               start_blocking_len=20)
+        answer = [linear_function(point_coords),
+                  2*linear_function(point_coords) ]
+        answer = transpose(answer)
+        #print "z",z
+        #print "answer",answer
+        assert allclose(z, answer)
+        assert allclose(interp.A_can_be_reused, True)
+
+        z = interp.interpolate(f)
+        assert allclose(z, answer)
+        
+        # This causes blocking to occur. 
+        z = interp.interpolate(f, start_blocking_len=10)
+        assert allclose(z, answer)
+        assert allclose(interp.A_can_be_reused, False)
+
+        #A is recalculated
+        z = interp.interpolate(f)
+        assert allclose(z, answer)
+        assert allclose(interp.A_can_be_reused, True)
+
+        interp = Interpolate(vertices, triangles)
+        #Must raise an exception, no points specified
+        try:
+            z = interp.interpolate(f)
+        except:
+            pass
+        
 
 #-------------------------------------------------------------