Changeset 6553 for branches/numpy/anuga/fit_interpolate
- Timestamp:
- Mar 19, 2009, 1:43:34 PM (16 years ago)
- Location:
- branches/numpy/anuga/fit_interpolate
- Files:
-
- 4 edited
-
fit.py (modified) (10 diffs)
-
interpolate.py (modified) (1 diff)
-
search_functions.py (modified) (7 diffs)
-
test_search_functions.py (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/fit_interpolate/fit.py
r6533 r6553 280 280 281 281 element_found, sigma0, sigma1, sigma2, k = \ 282 search_tree_of_vertices(self.root, self.mesh,x)282 search_tree_of_vertices(self.root, x) 283 283 284 284 if element_found is True: … … 301 301 AtA[j,k] += sigmas[j]*sigmas[k] 302 302 else: 303 # FIXME(Ole): This is the message referred to in ticket:314304 305 303 flag = is_inside_polygon(x, 306 304 self.mesh_boundary_polygon, … … 310 308 assert flag is True, msg 311 309 310 # FIXME(Ole): This is the message referred to in ticket:314 312 311 minx = min(self.mesh_boundary_polygon[:,0]) 313 312 maxx = max(self.mesh_boundary_polygon[:,0]) … … 317 316 msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\ 318 317 % (minx, maxx, miny, maxy) 319 320 321 raise Exception(msg) 318 raise RuntimeError, msg 319 322 320 self.AtA = AtA 323 321 … … 375 373 self.build_fit_subset(points, z, verbose=verbose) 376 374 375 # FIXME(Ole): I thought this test would make sense here 376 # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py 377 # Committed 11 March 2009 378 msg = 'Matrix AtA was not built' 379 assert self.AtA is not None, msg 380 381 #print 'Matrix was built OK' 382 377 383 378 384 point_coordinates = None … … 382 388 if point_coordinates is None: 383 389 if verbose: print 'Warning: no data points in fit' 384 #assert self.AtA != None, 'no interpolation matrix' 385 #assert self.Atz != None 386 assert not self.AtA is None, 'no interpolation matrix' 387 assert not self.Atz is None 390 msg = 'No interpolation matrix.' 391 assert self.AtA is not None, msg 392 assert self.Atz is not None 388 393 389 394 # FIXME (DSG) - do a message … … 471 476 alpha=DEFAULT_ALPHA, 472 477 verbose=False, 473 acceptable_overshoot=1.01,474 # FIXME: Move to config - this value is assumed in caching test475 # FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.476 478 mesh_origin=None, 477 479 data_origin=None, … … 490 492 'alpha': alpha, 491 493 'verbose': verbose, 492 'acceptable_overshoot': acceptable_overshoot,493 494 'mesh_origin': mesh_origin, 494 495 'data_origin': data_origin, … … 546 547 alpha=DEFAULT_ALPHA, 547 548 verbose=False, 548 acceptable_overshoot=1.01,549 549 mesh_origin=None, 550 550 data_origin=None, … … 572 572 alpha: Smoothing parameter. 573 573 574 acceptable overshoot: NOT IMPLEMENTED575 controls the allowed factor by which576 fitted values577 may exceed the value of input data. The lower limit is defined578 as min(z) - acceptable_overshoot*delta z and upper limit579 as max(z) + acceptable_overshoot*delta z580 581 582 574 mesh_origin: A geo_reference object or 3-tuples consisting of 583 575 UTM zone, easting and northing. 584 576 If specified vertex coordinates are assumed to be 585 577 relative to their respective origins. 586 587 578 588 579 point_attributes: Vector or array of data at the -
branches/numpy/anuga/fit_interpolate/interpolate.py
r6428 r6553 484 484 x = point_coordinates[i] 485 485 element_found, sigma0, sigma1, sigma2, k = \ 486 search_tree_of_vertices(self.root, self.mesh,x)487 488 486 search_tree_of_vertices(self.root, x) 487 488 # Update interpolation matrix A if necessary 489 489 if element_found is True: 490 490 # Assign values to matrix A -
branches/numpy/anuga/fit_interpolate/search_functions.py
r6304 r6553 8 8 import time 9 9 10 from anuga.utilities.polygon import is_inside_triangle 10 11 from anuga.utilities.numerical_tools import get_machine_precision 11 12 from anuga.config import max_float … … 18 19 search_more_cells_time = initial_search_value 19 20 20 #FIXME test what happens if a 21 LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]), 22 num.array([max_float,max_float]), 23 num.array([max_float,max_float])), 24 (num.array([1,1]), 25 num.array([0,0]), 26 num.array([-1.1,-1.1]))]]] 21 # FIXME(Ole): Could we come up with a less confusing structure? 22 LAST_TRIANGLE = [[-10, 23 (num.array([[max_float, max_float], 24 [max_float, max_float], 25 [max_float, max_float]]), 26 (num.array([1.,1.]), 27 num.array([0.,0.]), 28 num.array([-1.1,-1.1])))]] 27 29 28 def search_tree_of_vertices(root, mesh,x):30 def search_tree_of_vertices(root, x): 29 31 """ 30 32 Find the triangle (element) that the point x is in. … … 32 34 Inputs: 33 35 root: A quad tree of the vertices 34 mesh: The underlying mesh35 36 x: The point being placed 36 37 … … 47 48 global search_more_cells_time 48 49 49 #Find triangle containing x:50 element_found = False51 52 # This will be returned if element_found = False53 sigma2 = -10.054 sigma0 = -10.055 sigma1 = -10.056 k = -10.057 58 50 # Search the last triangle first 59 try: 60 element_found, sigma0, sigma1, sigma2, k = \ 61 _search_triangles_of_vertices(mesh, last_triangle, x) 62 except: 63 element_found = False 51 element_found, sigma0, sigma1, sigma2, k = \ 52 _search_triangles_of_vertices(last_triangle, x) 64 53 65 #print "last_triangle", last_triangle66 54 if element_found is True: 67 #print "last_triangle", last_triangle68 55 return element_found, sigma0, sigma1, sigma2, k 69 56 70 # This was only slightly faster than just checking the 71 # last triangle and it significantly slowed down 72 # non-gridded fitting 73 # If the last element was a dud, search its neighbours 74 #print "last_triangle[0][0]", last_triangle[0][0] 75 #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0]) 76 #print "neighbours", neighbours 77 #neighbours = [] 78 # for k in neighbours: 79 # if k >= 0: 80 # tri = mesh.get_vertex_coordinates(k, 81 # absolute=True) 82 # n0 = mesh.get_normal(k, 0) 83 # n1 = mesh.get_normal(k, 1) 84 # n2 = mesh.get_normal(k, 2) 85 # triangle =[[k,(tri, (n0, n1, n2))]] 86 # element_found, sigma0, sigma1, sigma2, k = \ 87 # _search_triangles_of_vertices(mesh, 88 # triangle, x) 89 # if element_found is True: 90 # return element_found, sigma0, sigma1, sigma2, k 91 92 #t0 = time.time() 57 93 58 # Get triangles in the cell that the point is in. 94 59 # Triangle is a list, first element triangle_id, 95 60 # second element the triangle 96 61 triangles = root.search(x[0], x[1]) 62 element_found, sigma0, sigma1, sigma2, k = \ 63 _search_triangles_of_vertices(triangles, x) 64 97 65 is_more_elements = True 98 66 99 element_found, sigma0, sigma1, sigma2, k = \100 _search_triangles_of_vertices(mesh,101 triangles, x)102 #search_one_cell_time += time.time()-t0103 #print "search_one_cell_time",search_one_cell_time104 #t0 = time.time()105 67 while not element_found and is_more_elements: 106 68 triangles, branch = root.expand_search() … … 109 71 # been searched. This is the last try 110 72 element_found, sigma0, sigma1, sigma2, k = \ 111 _search_triangles_of_vertices( mesh,triangles, x)73 _search_triangles_of_vertices(triangles, x) 112 74 is_more_elements = False 113 75 else: 114 76 element_found, sigma0, sigma1, sigma2, k = \ 115 _search_triangles_of_vertices(mesh, triangles, x) 116 #search_more_cells_time += time.time()-t0 117 #print "search_more_cells_time", search_more_cells_time 77 _search_triangles_of_vertices(triangles, x) 78 118 79 119 80 return element_found, sigma0, sigma1, sigma2, k 120 81 121 82 122 def _search_triangles_of_vertices( mesh,triangles, x):123 """Search for triangle containing x amongs candidate_vertices in mesh83 def _search_triangles_of_vertices(triangles, x): 84 """Search for triangle containing x amongs candidate_vertices in triangles 124 85 125 86 This is called by search_tree_of_vertices once the appropriate node … … 132 93 global last_triangle 133 94 134 # these statments are needed if triangles is empty 135 #Find triangle containing x: 136 element_found = False 137 138 # This will be returned if element_found = False 95 # These statments are needed if triangles is empty 139 96 sigma2 = -10.0 140 97 sigma0 = -10.0 141 98 sigma1 = -10.0 142 99 k = -10 143 144 #For all vertices in same cell as point x100 # For all vertices in same cell as point x 101 element_found = False 145 102 for k, tri_verts_norms in triangles: 146 103 tri = tri_verts_norms[0] 147 n0, n1, n2 = tri_verts_norms[1]148 104 # k is the triangle index 149 105 # tri is a list of verts (x, y), representing a tringle 150 106 # Find triangle that contains x (if any) and interpolate 151 element_found, sigma0, sigma1, sigma2 =\ 152 find_triangle_compute_interpolation(tri, n0, n1, n2, x) 153 if element_found is True: 107 108 # Input check disabled to speed things up. 109 if is_inside_triangle(x, tri, 110 closed=True, 111 check_inputs=False): 112 113 n0, n1, n2 = tri_verts_norms[1] 114 sigma0, sigma1, sigma2 =\ 115 compute_interpolation_values(tri, n0, n1, n2, x) 116 117 element_found = True 118 154 119 # Don't look for any other triangles in the triangle list 155 last_triangle = [[k,tri_verts_norms]] 156 break 120 last_triangle = [[k, tri_verts_norms]] 121 break 122 157 123 return element_found, sigma0, sigma1, sigma2, k 158 124 159 125 160 126 161 def find_triangle_compute_interpolation(triangle, n0, n1, n2, x): 162 """Compute linear interpolation of point x and triangle k in mesh. 163 It is assumed that x belongs to triangle k.max_float 127 def compute_interpolation_values(triangle, n0, n1, n2, x): 128 """Compute linear interpolation of point x and triangle. 129 130 n0, n1, n2 are normal to the tree edges. 164 131 """ 165 132 … … 167 134 xi0, xi1, xi2 = triangle 168 135 169 # this is where we can call some fast c code.170 171 # Integrity check - machine precision is too hard172 # so we use hardwired single precision173 epsilon = 1.0e-6174 175 if x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:176 # print "max(xi0[0], xi1[0], xi2[0])", max(xi0[0], xi1[0], xi2[0])177 return False,0,0,0178 if x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:179 return False,0,0,0180 if x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:181 return False,0,0,0182 if x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:183 return False,0,0,0184 185 # machine precision on some machines (e.g. nautilus)186 epsilon = get_machine_precision() * 2187 188 # Compute interpolation - return as soon as possible189 # print "(xi0-xi1)", (xi0-xi1)190 # print "n0", n0191 # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)192 193 136 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) 194 if sigma0 < -epsilon:195 return False,0,0,0196 137 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) 197 if sigma1 < -epsilon:198 return False,0,0,0199 138 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) 200 if sigma2 < -epsilon:201 return False,0,0,0202 203 # epsilon = 1.0e-6204 # we want to speed this up, so don't do assertions205 #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero206 #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\207 # %(delta, epsilon)208 #assert delta < epsilon, msg209 139 210 211 # Check that this triangle contains the data point 212 # Sigmas are allowed to get negative within 213 # machine precision on some machines (e.g. nautilus) 214 #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon: 215 # element_found = True 216 #else: 217 # element_found = False 218 return True, sigma0, sigma1, sigma2 140 return sigma0, sigma1, sigma2 219 141 220 142 def set_last_triangle(): -
branches/numpy/anuga/fit_interpolate/test_search_functions.py
r6360 r6553 5 5 from search_functions import search_tree_of_vertices, set_last_triangle 6 6 from search_functions import _search_triangles_of_vertices 7 from search_functions import compute_interpolation_values 7 8 8 9 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh … … 10 11 from anuga.utilities.polygon import is_inside_polygon 11 12 from anuga.utilities.quad import build_quadtree, Cell 12 13 from anuga.utilities.numerical_tools import ensure_numeric 14 from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle 15 16 import numpy as num 13 17 14 18 class Test_search_functions(unittest.TestCase): … … 49 53 50 54 x = [0.2, 0.7] 51 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,x)55 found, s0, s1, s2, k = search_tree_of_vertices(root, x) 52 56 assert k == 1 # Triangle one 53 57 assert found is True 54 58 55 59 def test_bigger(self): 56 """test_larger mesh 60 """test_bigger 61 62 test larger mesh 57 63 """ 58 64 … … 69 75 [0.1,0.9], [0.4,0.6], [0.9,0.1], 70 76 [10, 3]]: 71 72 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) 77 78 found, s0, s1, s2, k = search_tree_of_vertices(root, 79 ensure_numeric(x)) 73 80 74 81 if k >= 0: … … 102 109 [10, 3]]: 103 110 104 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,x)111 found, s0, s1, s2, k = search_tree_of_vertices(root, x) 105 112 106 113 if k >= 0: … … 111 118 assert found is False 112 119 113 120 121 if m == k == 0: return 114 122 115 123 def test_underlying_function(self): … … 124 132 125 133 # One point 126 x = [0.5, 0.5]134 x = ensure_numeric([0.5, 0.5]) 127 135 candidate_vertices = root.search(x[0], x[1]) 128 136 129 137 #print x, candidate_vertices 130 138 found, sigma0, sigma1, sigma2, k = \ 131 _search_triangles_of_vertices(mesh, 132 candidate_vertices, 139 _search_triangles_of_vertices(candidate_vertices, 133 140 x) 134 141 … … 151 158 #print x, candidate_vertices 152 159 found, sigma0, sigma1, sigma2, k = \ 153 _search_triangles_of_vertices(mesh, 154 candidate_vertices, 155 x) 160 _search_triangles_of_vertices(candidate_vertices, 161 ensure_numeric(x)) 156 162 if k >= 0: 157 163 V = mesh.get_vertex_coordinates(k) # nodes for triangle k … … 199 205 x = [2.5, 1.5] 200 206 element_found, sigma0, sigma1, sigma2, k = \ 201 search_tree_of_vertices(root, mesh,x)207 search_tree_of_vertices(root, x) 202 208 # One point 203 209 x = [3.00005, 2.999994] 204 210 element_found, sigma0, sigma1, sigma2, k = \ 205 search_tree_of_vertices(root, mesh,x)211 search_tree_of_vertices(root, x) 206 212 assert element_found is True 207 213 assert k == 1 208 214 209 ################################################################################ 210 215 216 def test_compute_interpolation_values(self): 217 """test_compute_interpolation_values 218 219 Test that interpolation values are correc 220 221 This test used to check element_found as output from 222 find_triangle_compute_interpolation(triangle, n0, n1, n2, x) 223 and that failed before 18th March 2009. 224 225 Now this function no longer returns this flag, so the test 226 is merely checknig the sigmas. 227 228 """ 229 230 triangle = num.array([[306951.77151059, 6194462.14596986], 231 [306952.58403545, 6194459.65001246], 232 [306953.55109034, 6194462.0041216]]) 233 234 235 n0 = [0.92499377, -0.37998227] 236 n1 = [0.07945684, 0.99683831] 237 n2 = [-0.95088404, -0.30954732] 238 239 x = [306953.344, 6194461.5] 240 241 # Test that point is indeed inside triangle 242 assert is_inside_polygon(x, triangle, 243 closed=True, verbose=False) 244 assert is_inside_triangle(x, triangle, 245 closed=True, verbose=False) 246 247 sigma0, sigma1, sigma2 = \ 248 compute_interpolation_values(triangle, n0, n1, n2, x) 249 250 msg = 'Point which is clearly inside triangle was not found' 251 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg 252 253 #------------------------------------------------------------- 211 254 if __name__ == "__main__": 212 suite = unittest.makeSuite(Test_search_functions, 'test')213 runner = unittest.TextTestRunner( ) #verbosity=1)255 suite = unittest.makeSuite(Test_search_functions, 'test') 256 runner = unittest.TextTestRunner(verbosity=1) 214 257 runner.run(suite) 215 258
Note: See TracChangeset
for help on using the changeset viewer.
