Changeset 441
- Timestamp:
- Oct 24, 2004, 10:08:40 AM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
-
least_squares.py (modified) (8 diffs)
-
test_least_squares.py (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/least_squares.py
r436 r441 20 20 #FIXME: Current implementation uses full matrices and a general solver. 21 21 #Later on we may consider using sparse techniques 22 23 import exceptions 24 class ShapeError(exceptions.Exception): pass 22 25 23 26 from general_mesh import General_Mesh … … 178 181 179 182 def build_coefficient_matrix_B(self, point_coordinates=None): 180 """Build final coefficient matrix 181 """182 #self.alpha = 0 183 """Build final coefficient matrix""" 184 185 183 186 if self.alpha <> 0: 184 187 self.build_smoothing_matrix_D() 185 188 186 187 189 if point_coordinates: 190 188 191 self.build_interpolation_matrix_A(point_coordinates) 189 #self.At = self.At.tocsc() 190 #print "self.At ",self.At 191 #print "self.A",self.A 192 193 #AtA = self.At.todense() * self.A.todense() 194 #print "AtA dense",AtA 195 #AtA = sparse.dok_matrix(AtA) 196 AtA = self.At * self.A 197 #print "AtA",AtA 192 198 193 if self.alpha <> 0: 199 self.B = AtA + self.alpha*self.D194 self.B = self.AtA + self.alpha*self.D 200 195 else: 201 self.B = AtA202 #print "end of building b" 196 self.B = self.AtA 197 203 198 204 199 def build_interpolation_matrix_A(self, point_coordinates): 205 200 """Build n x m interpolation matrix, where 206 201 n is the number of data points and 207 m is the number of basis functions phi_k (one per vertex) 208 """202 m is the number of basis functions phi_k (one per vertex)""" 203 209 204 #Convert input to Numeric arrays 210 205 point_coordinates = array(point_coordinates).astype(Float) … … 216 211 #self.A = zeros((n,m), Float) 217 212 self.A = sparse.dok_matrix() 213 self.AtA = sparse.dok_matrix() 218 214 219 215 #Compute matrix elements … … 249 245 250 246 #Assign values to matrix A 251 247 248 252 249 j0 = self.mesh.triangles[k,0] #Global vertex id 253 self.A[i, j0] = sigma0250 #self.A[i, j0] = sigma0 254 251 255 252 j1 = self.mesh.triangles[k,1] #Global vertex id 256 self.A[i, j1] = sigma1253 #self.A[i, j1] = sigma1 257 254 258 255 j2 = self.mesh.triangles[k,2] #Global vertex id 259 self.A[i, j2] = sigma2 260 261 262 #self.A = self.A.todense() 263 #Precompute 264 self.At = self.A.transp() 256 #self.A[i, j2] = sigma2 257 258 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} 259 js = [j0,j1,j2] 260 261 for j in js: 262 self.A[i,j] = sigmas[j] 263 for k in js: 264 self.AtA[j,k] += sigmas[j]*sigmas[k] 265 266 267 self.A = (self.A).tocsc() 268 self.AtA = (self.AtA).tocsc() 269 self.At = self.A.transp() 265 270 266 271 def get_A(self): … … 354 359 self.D[v2,v0] += e20 355 360 self.D[v0,v2] += e20 356 361 362 self.D = (self.D).tocsc() 357 363 358 364 def fit(self, z): 359 """Fit a smooth surface to given data points z.365 """Fit a smooth surface to given 1d array of data points z. 360 366 361 367 The smooth surface is computed at each vertex in the underlying … … 366 372 367 373 Inputs: 368 z: Vector or array of data at the point_coordinates. 369 If z is an array, smoothing will be done for each column 374 z: Single 1d vector or array of data at the point_coordinates. 370 375 """ 371 376 372 377 #Convert input to Numeric arrays 373 378 z = array(z).astype(Float) 379 380 381 if len(z.shape) > 1 : 382 raise VectorShapeError, 'Can only deal with 1d data vector' 374 383 375 384 #Compute right hand side based on data … … 392 401 393 402 return conjugate_gradient(self.B, Atz, Atz) 394 #FIXME: Should we store the result here for later use? (ON) 403 #FIXME: Should we store the result here for later use? (ON) 395 404 396 405 def fit_points(self, z): … … 432 441 If f is an array, interpolation will be done for each column 433 442 """ 434 435 return self.A * f 436 443 444 try: 445 return self.A * f 446 except ValueError: 447 # We are here, so it probalby means that f is 2 dimensional 448 # and so will do each column separately due to problems in 449 # sparse matrix, 2d array multiplication 450 451 (N , M) = self.A.shape 452 f = array(f).astype(Float) 453 (m , n) = f.shape 454 if m != M : 455 raise VectorShapeError, 'Mismatch between A and f dimensions' 456 457 y = zeros( (N,n), Float) 458 459 for i in range(n): 460 y[:,i] = self.A * f[:,i] 461 462 return y 437 463 438 464 #FIXME: We will need a method 'evaluate(self):' that will interpolate -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r436 r441 17 17 class TestCase(unittest.TestCase): 18 18 19 def setUp(self):20 pass21 22 def tearDown(self):23 pass19 ## def setUp(self): 20 ## pass 21 22 ## def tearDown(self): 23 ## pass 24 24 25 25 26 26 def test_datapoint_at_centroid(self): 27 27 28 a = [0.0, 0.0] 28 29 b = [0.0, 2.0] … … 39 40 40 41 def test_datapoints_at_vertices(self): 41 """Test that data points coinciding with vertices yield a diagonal matrix 42 """ 43 print "aa" 42 43 44 44 a = [0.0, 0.0] 45 45 b = [0.0, 2.0] … … 58 58 59 59 def test_datapoints_on_edge_midpoints(self): 60 """Try datapoints midway on edges - 61 each point should affect two matrix entries equally 62 """ 63 print "bb" 64 60 65 61 a = [0.0, 0.0] 66 62 b = [0.0, 2.0] … … 79 75 80 76 def test_datapoints_on_edges(self): 81 """Try datapoints on edges -82 each point should affect two matrix entries in proportion83 """84 print "cc"85 77 86 78 a = [0.0, 0.0] … … 99 91 100 92 def test_arbitrary_datapoints(self): 101 """Try arbitrary datapoints 102 """ 93 103 94 104 95 from Numeric import sum … … 117 108 118 109 119 # this causes a memory error in scipy.sparce 120 def bad_test_more_triangles(self): 121 print "dd" 110 def test_more_triangles(self): 111 122 112 a = [-1.0, 0.0] 123 113 b = [3.0, 4.0] … … 131 121 132 122 #Data points 133 data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] 123 data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3], [-1.0,-1.5 ], [1.0,-1.0]] 134 124 interp = Interpolation(points, triangles, data_points) 135 125 … … 147 137 148 138 def test_smooth_attributes_to_mesh(self): 149 print "ee"139 150 140 a = [0.0, 0.0] 151 141 b = [0.0, 5.0] … … 170 160 171 161 def test_smooth_att_to_meshII(self): 172 print "ff"162 173 163 a = [0.0, 0.0] 174 164 b = [0.0, 5.0] … … 189 179 190 180 def test_smooth_attributes_to_meshIII(self): 191 print "gg" 181 192 182 a = [-1.0, 0.0] 193 183 b = [3.0, 4.0] … … 221 211 222 212 def test_smooth_attributes_to_meshIV(self): 223 print "hh" 224 """ Testing 2 attributes smoothed to the mesh 225 """ 226 print "test_smooth_attributes_to_meshIV" 213 227 214 a = [0.0, 0.0] 228 215 b = [0.0, 5.0] … … 246 233 247 234 def test_interpolate_attributes_to_points(self): 235 248 236 v0 = [0.0, 0.0] 249 237 v1 = [0.0, 5.0] … … 268 256 269 257 def test_interpolate_attributes_to_pointsII(self): 258 270 259 a = [-1.0, 0.0] 271 260 b = [3.0, 4.0] … … 301 290 302 291 def test_interpolate_attributes_to_pointsIII(self): 292 303 293 v0 = [0.0, 0.0] 304 294 v1 = [0.0, 5.0] … … 328 318 329 319 def test_interpolate_attributes_to_pointsIV(self): 320 330 321 a = [-1.0, 0.0] 331 322 b = [3.0, 4.0] … … 366 357 367 358 def test_smooth_attributes_to_mesh_function(self): 368 """ Testing 2 attributes smoothed to the mesh369 """370 """Test multiple attributes371 """372 359 373 360 a = [0.0, 0.0] … … 433 420 434 421 def test_smoothing_matrix_more_triangles(self): 422 423 435 424 from Numeric import dot 436 425 … … 476 465 477 466 def test_fit_and_interpolation(self): 467 478 468 from mesh import Mesh 479 469 … … 663 653 if __name__ == "__main__": 664 654 suite = unittest.makeSuite(TestCase,'test') 665 runner = unittest.TextTestRunner( )655 runner = unittest.TextTestRunner(verbosity=2) 666 656 runner.run(suite) 667 657
Note: See TracChangeset
for help on using the changeset viewer.
