Changeset 330
- Timestamp:
- Sep 21, 2004, 4:57:47 PM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/least_squares.py
r328 r330 1 1 """Least squares smooting and interpolation. 2 Very first cut - just to check a few things. 3 Architecture is not set in concrete at all. 4 5 O 2 3 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou 4 Geoscience Australia, 2004. 6 5 """ 7 6 8 9 """Classes implementing general 2D geometrical mesh of triangles. 10 11 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou 12 Geoscience Australia, 2004 13 """ 14 15 EPSILON = 1.0e-13 7 #FIXME: Current implementation uses full matrices and a general solver. 8 #Later on we may consider using sparse techniques 16 9 17 10 from general_mesh import General_Mesh 18 from Numeric import zeros, array, Float, Int, dot, transpose11 from Numeric import zeros, array, Float, Int, dot, transpose 19 12 from LinearAlgebra import solve_linear_equations 20 21 13 22 14 def smooth_attributes_to_mesh(vertex_coordinates, … … 38 30 vertex_coordinates, 39 31 triangles, 40 point_coordinates = None): 32 point_coordinates = None, 33 alpha = 0.001): 41 34 42 35 … … 54 47 point_coordinates: List of coordinate pairs [x, y] of data points 55 48 (or an nx2 Numeric array) 49 50 alpha: Smoothing parameter 56 51 57 52 """ 53 58 54 59 55 #Convert input to Numeric arrays … … 64 60 self.mesh = General_Mesh(vertex_coordinates, triangles) 65 61 66 self.A = zeros((1,1), Float) #Not initialised 67 62 #Smoothing parameter 63 self.alpha = alpha 64 65 #Build coefficient matrices 66 self.build_coefficient_matrix_B(point_coordinates) 67 68 69 70 def build_coefficient_matrix_B(self, point_coordinates=None): 71 """Build final coefficient matrix 72 """ 73 74 self.build_smoothing_matrix_D() 75 68 76 if point_coordinates: 69 self.build_A(point_coordinates) 70 71 def build_A(self, point_coordinates): 72 """ 73 "A" is the matrix of hat functions. It has 1 row per data point and 1 column per vertex. 74 75 Post Condition 76 "A" has been initialised 77 """ 77 self.build_interpolation_matrix_A(point_coordinates) 78 79 self.B = self.AtA + self.alpha*self.D 80 else: 81 self.B = self.D 82 83 84 def build_interpolation_matrix_A(self, point_coordinates): 85 """Build n x m interpolation matrix, where 86 n is the number of data points and 87 m is the number of basis functions phi_k (one per vertex) 88 """ 78 89 79 90 #Convert input to Numeric arrays … … 81 92 82 93 #Build n x m interpolation matrix 83 m = self.mesh.coordinates.shape[0] #N umber of basis functions (1/vertex)84 n = point_coordinates.shape[0] #Number of data points94 m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) 95 n = point_coordinates.shape[0] #Nbr of data points 85 96 86 97 self.A = zeros((n,m), Float) … … 111 122 112 123 #FIXME: Maybe move out to test or something 113 assert abs(sigma0 + sigma1 + sigma2 - 1.0) < EPSILON114 115 124 epsilon = 1.0e-13 125 assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon 126 116 127 #Check that this triangle contains data point 117 128 if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0: 118 129 119 #Assign values to A130 #Assign values to matrix A 120 131 j = self.mesh.triangles[k,0] #Global vertex id 121 132 self.A[i, j] = sigma0 … … 126 137 j = self.mesh.triangles[k,2] #Global vertex id 127 138 self.A[i, j] = sigma2 139 128 140 129 130 def smooth_attributes_to_mesh(self, z): 131 """ Linearly interpolate many points with an attribute to the vertices. 141 #Precompute shorthands 142 self.At = transpose(self.A) 143 self.AtA = dot(self.At, self.A) 144 145 146 def build_smoothing_matrix_D(self): 147 """Build m x m smoothing matrix, where 148 m is the number of basis functions phi_k (one per vertex) 149 150 The smoothing matrix is defined as 151 152 D = D1 + D2 153 154 where 155 156 [D1]_{k,l} = \int_\Omega 157 \frac{\partial \phi_k}{\partial x} 158 \frac{\partial \phi_l}{\partial x}\, 159 dx dy 160 161 [D2]_{k,l} = \int_\Omega 162 \frac{\partial \phi_k}{\partial y} 163 \frac{\partial \phi_l}{\partial y}\, 164 dx dy 165 166 167 The derivatives \frac{\partial \phi_k}{\partial x}, 168 \frac{\partial \phi_k}{\partial x} for a particular triangle 169 are obtained by computing the gradient a_k, b_k for basis function k 170 """ 171 172 #FIXME: algorithm might be optimised by computing local 9x9 173 #"element stiffness matrices: 174 175 from util import gradient 176 177 m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) 178 179 self.D = zeros((m,m), Float) 180 181 #For each triangle compute contributions to D = D1+D2 182 for i in range(len(self.mesh)): 183 184 #Get area 185 area = self.mesh.areas[i] 186 187 #Get global vertex indices 188 v0 = self.mesh.triangles[i,0] 189 v1 = self.mesh.triangles[i,1] 190 v2 = self.mesh.triangles[i,2] 191 192 #Get the three vertex_points 193 xi0 = self.mesh.get_vertex_coordinate(i, 0) 194 xi1 = self.mesh.get_vertex_coordinate(i, 1) 195 xi2 = self.mesh.get_vertex_coordinate(i, 2) 196 197 #Compute gradients for each vertex 198 a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1], 199 1, 0, 0) 200 201 a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1], 202 0, 1, 0) 203 204 a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1], 205 0, 0, 1) 206 207 #Compute diagonal contributions 208 self.D[v0,v0] += (a0*a0 + b0*b0)*area 209 self.D[v1,v1] += (a1*a1 + b1*b1)*area 210 self.D[v2,v2] += (a2*a2 + b2*b2)*area 211 212 #Compute contributions for basis functions sharing edges 213 e01 = (a0*a1 + b0*b1)*area 214 self.D[v0,v1] += e01 215 self.D[v1,v0] += e01 216 217 e12 = (a1*a2 + b1*b2)*area 218 self.D[v1,v2] += e12 219 self.D[v2,v1] += e12 220 221 e20 = (a2*a0 + b2*b0)*area 222 self.D[v2,v0] += e20 223 self.D[v0,v2] += e20 224 225 226 def fit(self, z): 227 """Fit a smooth surface to given data points z. 228 229 The smooth surface is computed at each vertex in the underlying 230 mesh using the formula given in the module doc string. 231 132 232 Pre Condition: 133 "A" hasbeen initialised233 self.A, self.At and self.B have been initialised 134 234 135 235 Inputs: 136 z: Vector of data at the point_coordinates. One value per point. 137 138 """ 236 z: Vector or array of data at the point_coordinates. 237 If z is an array, smoothing will be done for each column 238 """ 239 139 240 #Convert input to Numeric arrays 140 241 z = array(z).astype(Float) 141 At = transpose(self.A) 142 #print "z", z 143 #print "At",At 144 self.AtA = dot(At, self.A) 145 Atz = dot(At, z) 146 f = solve_linear_equations(self.AtA,Atz) 147 return f 148 149 150 def interpolate_attributes_to_points(self, f): 151 """ Linearly interpolate vertices with attributes to the points. 242 243 #Compute right hand side based on data 244 Atz = dot(self.At, z) 245 246 #Check sanity 247 n, m = self.A.shape 248 if n<m and self.alpha == 0.0: 249 msg = 'ERROR (least_squares): Too few data points\n' 250 msg += 'There only %d data points. Need at least %d\n' %(n,m) 251 msg += 'Alternatively, increase smoothing parameter alpha' 252 raise msg 253 254 255 #Solve and return 256 return solve_linear_equations(self.B, Atz) 257 258 #FIXME: Should we store the result here for later use? (ON) 259 260 261 def interpolate(self, f): 262 """Compute predicted values at data points implied in self.A. 263 264 The mesh values representing a smooth surface are 265 assumed to be specified in f. This argument could, 266 for example have been obtained from the method self.fit() 267 152 268 Pre Condition: 153 A has been initialised269 self.A has been initialised 154 270 155 271 Inputs: 156 f: Matrix of data at the vertices. 157 One attribute per colum. 158 159 """ 160 return dot(self.A,f) 272 f: Vector or array of data at the mesh vertices. 273 If f is an array, interpolation will be done for each column 274 """ 275 276 return dot(self.A, f) 277 278 279 #FIXME: We may need a method 'evaluate(self):' that will interpolate 280 #a computed surface living on the mesh onto a collection of 281 #arbitrary data points 282 # 283 #Precondition: self.fit(z) has stored its result in self.f. 284 # 285 #Input: data_points 286 #Algorithm: 287 # 1 Build a new temporary A matrix based on mesh and new data points 288 # 2 Apply it to self.f (return A*self.f) 289 # 290 # ON 291 292 161 293 162 294 #------------------------------------------------------------- … … 191 323 print f 192 324 193 325 326 327 328 329 330 331 332 333 334 335
Note: See TracChangeset
for help on using the changeset viewer.