Changeset 330 for inundation/ga/storm_surge/pyvolution/least_squares.py

Timestamp:

Sep 21, 2004, 4:57:47 PM (20 years ago)

Author:

ole

Message:

Resolved conflicts

File:

• inundation/ga/storm_surge/pyvolution/least_squares.py (modified) (8 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -1 +1 @@
 """Least squares smooting and interpolation.
-   Very first cut - just to check a few things.
-   Architecture is not set in concrete at all.
-
-   O
+
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia, 2004.   
 """
 
-
-"""Classes implementing general 2D geometrical mesh of triangles.
-
-   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-   Geoscience Australia, 2004   
-"""
-
-EPSILON = 1.0e-13
+#FIXME: Current implementation uses full matrices and a general solver.
+#Later on we may consider using sparse techniques
 
 from general_mesh import General_Mesh
-from Numeric import zeros, array, Float, Int, dot,transpose
+from Numeric import zeros, array, Float, Int, dot, transpose
 from LinearAlgebra import solve_linear_equations
-
 
 def smooth_attributes_to_mesh(vertex_coordinates,
@@ -38 +30 @@
                  vertex_coordinates,
                  triangles,
-                 point_coordinates = None):
+                 point_coordinates = None,
+                 alpha = 0.001):
 
         
@@ -54 +47 @@
           point_coordinates: List of coordinate pairs [x, y] of data points
           (or an nx2 Numeric array)
+
+          alpha: Smoothing parameter
           
         """
+
 
         #Convert input to Numeric arrays
@@ -64 +60 @@
         self.mesh = General_Mesh(vertex_coordinates, triangles)
 
-        self.A = zeros((1,1), Float) #Not initialised
-
+        #Smoothing parameter
+        self.alpha = alpha
+
+        #Build coefficient matrices
+        self.build_coefficient_matrix_B(point_coordinates)    
+
+
+        
+    def build_coefficient_matrix_B(self, point_coordinates=None):
+        """Build final coefficient matrix
+        """
+
+        self.build_smoothing_matrix_D()
+        
         if point_coordinates:
-            self.build_A(point_coordinates)
-        
-    def build_A(self, point_coordinates):
-        """
-        "A" is the matrix of hat functions.  It has 1 row per data point and 1 column per vertex.
-
-        Post Condition
-        "A" has been initialised
-        """
+            self.build_interpolation_matrix_A(point_coordinates)
+
+            self.B = self.AtA + self.alpha*self.D
+        else:
+            self.B = self.D
+
+        
+    def build_interpolation_matrix_A(self, point_coordinates):
+        """Build n x m interpolation matrix, where
+        n is the number of data points and
+        m is the number of basis functions phi_k (one per vertex)
+        """ 
 
         #Convert input to Numeric arrays
@@ -81 +92 @@
         
         #Build n x m interpolation matrix        
-        m = self.mesh.coordinates.shape[0] #Number of basis functions (1/vertex)
-        n = point_coordinates.shape[0]               #Number of data points         
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        n = point_coordinates.shape[0]     #Nbr of data points         
 
         self.A = zeros((n,m), Float)
@@ -111 +122 @@
 
                 #FIXME: Maybe move out to test or something
-                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < EPSILON
-
-                
+                epsilon = 1.0e-13
+                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
+
                 #Check that this triangle contains data point
                 if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
 
-                    #Assign values to A
+                    #Assign values to matrix A
                     j = self.mesh.triangles[k,0] #Global vertex id
                     self.A[i, j] = sigma0
@@ -126 +137 @@
                     j = self.mesh.triangles[k,2] #Global vertex id
                     self.A[i, j] = sigma2
+
     
-
-    def smooth_attributes_to_mesh(self, z):
-        """ Linearly interpolate many points with an attribute to the vertices.
+        #Precompute shorthands
+        self.At = transpose(self.A)
+        self.AtA = dot(self.At, self.A)
+
+        
+    def build_smoothing_matrix_D(self):
+        """Build m x m smoothing matrix, where
+        m is the number of basis functions phi_k (one per vertex)
+
+        The smoothing matrix is defined as
+
+        D = D1 + D2
+
+        where
+
+        [D1]_{k,l} = \int_\Omega
+           \frac{\partial \phi_k}{\partial x}
+           \frac{\partial \phi_l}{\partial x}\,
+           dx dy
+
+        [D2]_{k,l} = \int_\Omega
+           \frac{\partial \phi_k}{\partial y}
+           \frac{\partial \phi_l}{\partial y}\,
+           dx dy
+
+
+        The derivatives \frac{\partial \phi_k}{\partial x},
+        \frac{\partial \phi_k}{\partial x} for a particular triangle
+        are obtained by computing the gradient a_k, b_k for basis function k
+        """
+
+        #FIXME: algorithm might be optimised by computing local 9x9
+        #"element stiffness matrices:
+
+        from util import gradient
+        
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+
+        self.D = zeros((m,m), Float)
+
+        #For each triangle compute contributions to D = D1+D2        
+        for i in range(len(self.mesh)):
+
+            #Get area
+            area = self.mesh.areas[i]
+
+            #Get global vertex indices
+            v0 = self.mesh.triangles[i,0]
+            v1 = self.mesh.triangles[i,1]
+            v2 = self.mesh.triangles[i,2]
+            
+            #Get the three vertex_points
+            xi0 = self.mesh.get_vertex_coordinate(i, 0)
+            xi1 = self.mesh.get_vertex_coordinate(i, 1)
+            xi2 = self.mesh.get_vertex_coordinate(i, 2)                 
+
+            #Compute gradients for each vertex
+            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              1, 0, 0)
+
+            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              0, 1, 0)
+
+            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              0, 0, 1)            
+
+            #Compute diagonal contributions
+            self.D[v0,v0] += (a0*a0 + b0*b0)*area
+            self.D[v1,v1] += (a1*a1 + b1*b1)*area
+            self.D[v2,v2] += (a2*a2 + b2*b2)*area            
+
+            #Compute contributions for basis functions sharing edges
+            e01 = (a0*a1 + b0*b1)*area
+            self.D[v0,v1] += e01
+            self.D[v1,v0] += e01
+
+            e12 = (a1*a2 + b1*b2)*area
+            self.D[v1,v2] += e12
+            self.D[v2,v1] += e12
+
+            e20 = (a2*a0 + b2*b0)*area
+            self.D[v2,v0] += e20
+            self.D[v0,v2] += e20             
+            
+            
+    def fit(self, z):
+        """Fit a smooth surface to given data points z.
+
+        The smooth surface is computed at each vertex in the underlying
+        mesh using the formula given in the module doc string. 
+
         Pre Condition:
-          "A" has been initialised
+          self.A, self.At and self.B have been initialised
           
         Inputs:
-          z: Vector of data at the point_coordinates.  One value per point.
-          
-        """
+          z: Vector or array of data at the point_coordinates.
+          If z is an array, smoothing will be done for each column
+        """
+
         #Convert input to Numeric arrays
         z = array(z).astype(Float)
-        At = transpose(self.A)
-        #print "z", z
-        #print "At",At
-        self.AtA = dot(At, self.A) 
-        Atz = dot(At, z)
-        f = solve_linear_equations(self.AtA,Atz)
-        return f
-
-        
-    def interpolate_attributes_to_points(self, f):
-        """ Linearly interpolate vertices with attributes to the points.
+        
+        #Compute right hand side based on data
+        Atz = dot(self.At, z)
+
+        #Check sanity
+        n, m = self.A.shape
+        if n<m and self.alpha == 0.0:
+            msg = 'ERROR (least_squares): Too few data points\n'
+            msg += 'There only %d data points. Need at least %d\n' %(n,m)
+            msg += 'Alternatively, increase smoothing parameter alpha' 
+            raise msg
+
+
+        #Solve and return
+        return solve_linear_equations(self.B, Atz)
+
+        #FIXME: Should we store the result here for later use? (ON)
+    
+
+    def interpolate(self, f):
+        """Compute predicted values at data points implied in self.A. 
+
+        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:
-          A has been initialised
+          self.A has been initialised
         
         Inputs:
-          f: Matrix of data at the  vertices. 
-               One attribute per colum. 
-          
-        """
-        return dot(self.A,f)
+          f: Vector or array of data at the mesh vertices.
+          If f is an array, interpolation will be done for each column
+        """
+        
+        return dot(self.A, f)
+
+     
+    #FIXME: We may need a method 'evaluate(self):' that will interpolate
+    #a computed surface living on the mesh onto a collection of
+    #arbitrary data points
+    #
+    #Precondition: self.fit(z) has stored its result in self.f.
+    #
+    #Input: data_points
+    #Algorithm:
+    #  1 Build a new temporary A matrix based on mesh and new data points
+    #  2 Apply it to self.f (return A*self.f)
+    #
+    # ON
+
+
 
 #-------------------------------------------------------------
@@ -191 +323 @@
         print f
         
-        
+
+
+
+
+
+
+
+
+
+
+