Changeset 424
- Timestamp:
- Oct 19, 2004, 6:19:43 PM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/least_squares.py
r405 r424 24 24 from Numeric import zeros, array, Float, Int, dot, transpose 25 25 from LinearAlgebra import solve_linear_equations 26 from scipy import sparse 27 from cg_solve import conjugate_gradient 26 28 27 29 try: … … 176 178 """Build final coefficient matrix 177 179 """ 180 #self.alpha = 0 178 181 if self.alpha <> 0: 179 182 self.build_smoothing_matrix_D() … … 182 185 if point_coordinates: 183 186 self.build_interpolation_matrix_A(point_coordinates) 184 AtA = dot(self.At, self.A) 187 #self.At = self.At.tocsc() 188 #print "self.At ",self.At 189 #print "self.A",self.A 185 190 191 #AtA = self.At.todense() * self.A.todense() 192 #print "AtA dense",AtA 193 #AtA = sparse.dok_matrix(AtA) 194 AtA = self.At * self.A 195 #print "AtA",AtA 186 196 if self.alpha <> 0: 187 197 self.B = AtA + self.alpha*self.D 188 198 else: 189 199 self.B = AtA 190 200 #print "end of building b" 191 201 192 202 def build_interpolation_matrix_A(self, point_coordinates): … … 195 205 m is the number of basis functions phi_k (one per vertex) 196 206 """ 197 198 207 #Convert input to Numeric arrays 199 208 point_coordinates = array(point_coordinates).astype(Float) … … 203 212 n = point_coordinates.shape[0] #Nbr of data points 204 213 205 self.A = zeros((n,m), Float) 214 #self.A = zeros((n,m), Float) 215 self.A = sparse.dok_matrix() 206 216 207 217 #Compute matrix elements … … 237 247 238 248 #Assign values to matrix A 239 j = self.mesh.triangles[k,0] #Global vertex id 240 self.A[i, j] = sigma0 241 242 j = self.mesh.triangles[k,1] #Global vertex id 243 self.A[i, j] = sigma1 244 245 j = self.mesh.triangles[k,2] #Global vertex id 246 self.A[i, j] = sigma2 247 248 249 250 j0 = self.mesh.triangles[k,0] #Global vertex id 251 self.A[i, j0] = sigma0 252 253 j1 = self.mesh.triangles[k,1] #Global vertex id 254 self.A[i, j1] = sigma1 255 256 j2 = self.mesh.triangles[k,2] #Global vertex id 257 self.A[i, j2] = sigma2 258 259 260 #self.A = self.A.todense() 249 261 #Precompute 250 self.At = transpose(self.A) 251 252 262 self.At = self.A.transp() 263 264 def get_A(self): 265 return self.A.todense() 266 267 def get_B(self): 268 return self.B.todense() 269 270 def get_D(self): 271 return self.D.todense() 272 253 273 #FIXME: Remember to re-introduce the 1/n factor in the 254 274 #interpolation term … … 286 306 m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) 287 307 288 self.D = zeros((m,m), Float) 308 #self.D = zeros((m,m), Float) 309 self.D = sparse.dok_matrix() 289 310 290 311 #For each triangle compute contributions to D = D1+D2 … … 351 372 352 373 #Compute right hand side based on data 353 Atz = dot(self.At, z)374 Atz = self.At * z 354 375 355 376 #Check sanity … … 363 384 364 385 #Solve and return 365 return solve_linear_equations(self.B, Atz) 366 386 #return solve_linear_equations(self.get_B(), Atz) 387 388 # caused errors 389 #return sparse.solve(self.B,Atz) 390 391 return conjugate_gradient(self.B, Atz, Atz) 367 392 #FIXME: Should we store the result here for later use? (ON) 368 393 … … 383 408 """ 384 409 385 return dot(self.A, f)410 return self.A * f 386 411 387 412 -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r393 r424 34 34 35 35 interp = Interpolation(points, vertices, data) 36 assert allclose(interp. A, [[1./3, 1./3, 1./3]])36 assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 37 37 38 38 … … 41 41 """Test that data points coinciding with vertices yield a diagonal matrix 42 42 """ 43 43 print "aa" 44 44 a = [0.0, 0.0] 45 45 b = [0.0, 2.0] … … 51 51 52 52 interp = Interpolation(points, vertices, data) 53 assert allclose(interp. A, [[1., 0., 0.],53 assert allclose(interp.get_A(), [[1., 0., 0.], 54 54 [0., 1., 0.], 55 55 [0., 0., 1.]]) … … 61 61 each point should affect two matrix entries equally 62 62 """ 63 print "bb" 63 64 64 65 a = [0.0, 0.0] … … 72 73 interp = Interpolation(points, vertices, data) 73 74 74 assert allclose(interp. A, [[0.5, 0.5, 0.0], #Affects vertex 1 and 075 assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 75 76 [0.5, 0.0, 0.5], #Affects vertex 0 and 2 76 77 [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2 … … 81 82 each point should affect two matrix entries in proportion 82 83 """ 84 print "cc" 83 85 84 86 a = [0.0, 0.0] … … 92 94 interp = Interpolation(points, vertices, data) 93 95 94 assert allclose(interp. A, [[0.25, 0.75, 0.0], #Affects vertex 1 and 096 assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0 95 97 [0.25, 0.0, 0.75], #Affects vertex 0 and 2 96 98 [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2 … … 112 114 interp = Interpolation(points, vertices, data) 113 115 114 assert allclose(sum(interp. A, axis=1), 1.0)116 assert allclose(sum(interp.get_A(), axis=1), 1.0) 115 117 116 118 117 119 118 def test_more_triangles(self): 120 def teszt_more_triangles(self): 121 print "dd" 119 122 a = [-1.0, 0.0] 120 123 b = [3.0, 4.0] … … 128 131 129 132 #Data points 130 data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ] 131 132 interp = Interpolation(points, triangles, data) 133 data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ] 134 interp = Interpolation(points, triangles, data_points) 133 135 134 136 answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d … … 138 140 [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b 139 141 [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f 140 assert allclose(interp.A, answer) 142 143 assert allclose(interp.get_A(), answer) 141 144 142 145 … … 144 147 145 148 def test_smooth_attributes_to_mesh(self): 149 print "ee" 146 150 a = [0.0, 0.0] 147 151 b = [0.0, 5.0] … … 158 162 data_coords = [d1, d2, d3] 159 163 160 interp = Interpolation(points, triangles, data_coords, alpha=0.0 )164 interp = Interpolation(points, triangles, data_coords, alpha=0.0000005) 161 165 z = [z1, z2, z3] 162 166 f = interp.fit(z) … … 167 171 168 172 def test_smooth_att_to_meshII(self): 173 print "ff" 169 174 a = [0.0, 0.0] 170 175 b = [0.0, 5.0] … … 185 190 186 191 def test_smooth_attributes_to_meshIII(self): 192 print "gg" 187 193 a = [-1.0, 0.0] 188 194 b = [3.0, 4.0] … … 216 222 217 223 def test_smooth_attributes_to_meshIV(self): 224 print "hh" 218 225 """ Testing 2 attributes smoothed to the mesh 219 226 """ 227 print "test_smooth_attributes_to_meshIV" 220 228 a = [0.0, 0.0] 221 229 b = [0.0, 5.0] … … 398 406 interp = Interpolation(points, vertices) 399 407 400 assert allclose(interp. D, [[1, -0.5, -0.5],408 assert allclose(interp.get_D(), [[1, -0.5, -0.5], 401 409 [-0.5, 0.5, 0], 402 410 [-0.5, 0, 0.5]]) … … 407 415 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 408 416 # int 1 dx dy = area = 2 409 assert dot(dot(f, interp. D), f) == 2417 assert dot(dot(f, interp.get_D()), f) == 2 410 418 411 419 #Define f(x,y) = y … … 414 422 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 415 423 # int 1 dx dy = area = 2 416 assert dot(dot(f, interp. D), f) == 2424 assert dot(dot(f, interp.get_D()), f) == 2 417 425 418 426 #Define f(x,y) = x+y … … 421 429 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 422 430 # int 2 dx dy = 2*area = 4 423 assert dot(dot(f, interp. D), f) == 4431 assert dot(dot(f, interp.get_D()), f) == 4 424 432 425 433 … … 442 450 443 451 444 #assert allclose(interp. D, [[1, -0.5, -0.5],452 #assert allclose(interp.get_D(), [[1, -0.5, -0.5], 445 453 # [-0.5, 0.5, 0], 446 454 # [-0.5, 0, 0.5]]) … … 451 459 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 452 460 # int 1 dx dy = total area = 8 453 assert dot(dot(f, interp. D), f) == 8461 assert dot(dot(f, interp.get_D()), f) == 8 454 462 455 463 #Define f(x,y) = y … … 458 466 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 459 467 # int 1 dx dy = area = 8 460 assert dot(dot(f, interp. D), f) == 8468 assert dot(dot(f, interp.get_D()), f) == 8 461 469 462 470 #Define f(x,y) = x+y … … 465 473 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 466 474 # int 2 dx dy = 2*area = 16 467 assert dot(dot(f, interp. D), f) == 16475 assert dot(dot(f, interp.get_D()), f) == 16 468 476 469 477
Note: See TracChangeset
for help on using the changeset viewer.