Changeset 6152
- Timestamp:
- Jan 13, 2009, 2:20:31 PM (16 years ago)
- Location:
- anuga_core/source/anuga/fit_interpolate
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py
r4872 r6152 23 23 import profile , pstats 24 24 from math import sqrt 25 from Numeric import array26 25 27 26 from anuga.fit_interpolate.search_functions import search_times, \ … … 38 37 from anuga.fit_interpolate.general_fit_interpolate import \ 39 38 get_build_quadtree_time 39 40 40 41 41 """ -
anuga_core/source/anuga/fit_interpolate/fit.py
r5855 r6152 28 28 import types 29 29 30 from Numeric import zeros, Float, ArrayType,take, Int31 32 30 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 33 31 from anuga.caching import cache … … 46 44 class TooFewPointsError(exceptions.Exception): pass 47 45 class VertsWithNoTrianglesError(exceptions.Exception): pass 46 47 import Numeric as num 48 48 49 49 50 #DEFAULT_ALPHA = 0.001 … … 252 253 if len(z.shape) > 1: 253 254 att_num = z.shape[1] 254 self.Atz = zeros((m,att_num),Float)255 self.Atz = num.zeros((m,att_num), num.Float) 255 256 else: 256 257 att_num = 1 257 self.Atz = zeros((m,),Float)258 self.Atz = num.zeros((m,), num.Float) 258 259 assert z.shape[0] == point_coordinates.shape[0] 259 260 … … 437 438 # Convert input to Numeric arrays 438 439 if z is not None: 439 z = ensure_numeric(z, Float)440 z = ensure_numeric(z, num.Float) 440 441 else: 441 442 msg = 'z not specified' … … 443 444 z = point_coordinates.get_attributes(attribute_name) 444 445 445 point_coordinates = ensure_numeric(point_coordinates, Float)446 point_coordinates = ensure_numeric(point_coordinates, num.Float) 446 447 self._build_matrix_AtA_Atz(point_coordinates, z, verbose) 447 448 … … 584 585 585 586 #Convert input to Numeric arrays 586 triangles = ensure_numeric(triangles, Int)587 triangles = ensure_numeric(triangles, num.Int) 587 588 vertex_coordinates = ensure_absolute(vertex_coordinates, 588 589 geo_reference = mesh_origin) … … 651 652 vertex_coordinates = mesh_dict['vertices'] 652 653 triangles = mesh_dict['triangles'] 653 if type(mesh_dict['vertex_attributes']) == ArrayType:654 if type(mesh_dict['vertex_attributes']) == num.ArrayType: 654 655 old_point_attributes = mesh_dict['vertex_attributes'].tolist() 655 656 else: 656 657 old_point_attributes = mesh_dict['vertex_attributes'] 657 658 658 if type(mesh_dict['vertex_attribute_titles']) == ArrayType:659 if type(mesh_dict['vertex_attribute_titles']) == num.ArrayType: 659 660 old_title_list = mesh_dict['vertex_attribute_titles'].tolist() 660 661 else: -
anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
r5361 r6152 21 21 from warnings import warn 22 22 23 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \24 ArrayType, allclose, take, NewAxis, arange25 26 23 from anuga.caching.caching import cache 27 24 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh … … 37 34 ensure_absolute 38 35 from anuga.fit_interpolate.search_functions import set_last_triangle 36 37 import Numeric as num 38 39 39 40 40 # tests fail if 2 is used … … 97 97 98 98 #Convert input to Numeric arrays 99 triangles = ensure_numeric(triangles, Int)99 triangles = ensure_numeric(triangles, num.Int) 100 100 vertex_coordinates = ensure_absolute(vertex_coordinates, 101 101 geo_reference = mesh_origin) -
anuga_core/source/anuga/fit_interpolate/interpolate.py
r6086 r6152 24 24 from math import sqrt 25 25 from csv import writer, DictWriter 26 27 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \28 ArrayType, allclose, take, NewAxis, arange29 26 30 27 from anuga.caching.caching import cache … … 41 38 from anuga.abstract_2d_finite_volumes.util import file_function 42 39 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 40 41 import Numeric as num 42 43 43 44 44 # Interpolation specific exceptions … … 130 130 131 131 # Create interpolation object with matrix 132 args = (ensure_numeric(vertex_coordinates, Float),132 args = (ensure_numeric(vertex_coordinates, num.Float), 133 133 ensure_numeric(triangles)) 134 134 kwargs = {'mesh_origin': mesh_origin, … … 253 253 254 254 from utilities.polygon import point_on_line 255 from Numeric import ones 256 257 z = ones(len(point_coordinates), Float) 255 256 z = num.ones(len(point_coordinates), num.Float) 258 257 259 258 # input sanity check … … 379 378 # creating a dummy array to concatenate to. 380 379 381 f = ensure_numeric(f, Float)380 f = ensure_numeric(f, num.Float) 382 381 if len(f.shape) > 1: 383 z = zeros((0, f.shape[1]))382 z = num.zeros((0, f.shape[1])) 384 383 else: 385 z = zeros((0,))384 z = num.zeros((0,)) 386 385 387 386 for end in range(start_blocking_len, … … 390 389 t = self.interpolate_block(f, point_coordinates[start:end], 391 390 verbose=verbose) 392 z = concatenate((z, t))391 z = num.concatenate((z, t)) 393 392 start = end 394 393 … … 396 395 t = self.interpolate_block(f, point_coordinates[start:end], 397 396 verbose=verbose) 398 z = concatenate((z, t))397 z = num.concatenate((z, t)) 399 398 return z 400 399 … … 426 425 427 426 # Convert lists to Numeric arrays if necessary 428 point_coordinates = ensure_numeric(point_coordinates, Float)429 f = ensure_numeric(f, Float)427 point_coordinates = ensure_numeric(point_coordinates, num.Float) 428 f = ensure_numeric(f, num.Float) 430 429 431 430 from anuga.caching import myhash 432 from Numeric import alltrue433 431 import sys 434 432 … … 452 450 if self.interpolation_matrices.has_key(key): 453 451 X, stored_points = self.interpolation_matrices[key] 454 if alltrue(stored_points == point_coordinates):452 if num.alltrue(stored_points == point_coordinates): 455 453 reuse_A = True # Reuse interpolation matrix 456 454 … … 534 532 535 533 # Convert point_coordinates to Numeric arrays, in case it was a list. 536 point_coordinates = ensure_numeric(point_coordinates, Float)534 point_coordinates = ensure_numeric(point_coordinates, num.Float) 537 535 538 536 if verbose: print 'Getting indices inside mesh boundary' … … 838 836 """ 839 837 840 from Numeric import array, zeros, Float, alltrue, concatenate,\841 reshape, ArrayType842 843 838 from anuga.config import time_format 844 839 import types … … 847 842 time = ensure_numeric(time) 848 843 msg = 'Time must be a monotonuosly increasing sequence %s' % time 849 assert alltrue(time[1:] - time[:-1] >= 0), msg844 assert num.alltrue(time[1:] - time[:-1] >= 0), msg 850 845 851 846 # Check if quantities is a single array only … … 875 870 # Thin timesteps if needed 876 871 # Note array() is used to make the thinned arrays contiguous in memory 877 self.time = array(time[::time_thinning])872 self.time = num.array(time[::time_thinning]) 878 873 for name in quantity_names: 879 874 if len(quantities[name].shape) == 2: 880 quantities[name] = array(quantities[name][::time_thinning,:])875 quantities[name] = num.array(quantities[name][::time_thinning,:]) 881 876 882 877 # Save for use with statistics … … 941 936 # FIXME (Ole): Why only Windoze? 942 937 from anuga.utilities.polygon import plot_polygons 943 #out_interp_pts = take(interpolation_points,[indices])938 #out_interp_pts = num.take(interpolation_points,[indices]) 944 939 title = ('Interpolation points fall ' 945 940 'outside specified mesh') … … 986 981 987 982 for name in quantity_names: 988 self.precomputed_values[name] = zeros((p, m),Float)983 self.precomputed_values[name] = num.zeros((p, m), num.Float) 989 984 990 985 # Build interpolator … … 1086 1081 1087 1082 from math import pi, cos, sin, sqrt 1088 from Numeric import zeros, Float1089 1083 from anuga.abstract_2d_finite_volumes.util import mean 1090 1084 … … 1122 1116 1123 1117 # Compute interpolated values 1124 q = zeros(len(self.quantity_names),Float)1118 q = num.zeros(len(self.quantity_names), num.Float) 1125 1119 for i, name in enumerate(self.quantity_names): 1126 1120 Q = self.precomputed_values[name] … … 1166 1160 return q 1167 1161 else: 1168 from Numeric import ones, Float1169 1162 # x is a vector - Create one constant column for each value 1170 1163 N = len(x) … … 1172 1165 res = [] 1173 1166 for col in q: 1174 res.append(col* ones(N,Float))1167 res.append(col*num.ones(N, num.Float)) 1175 1168 1176 1169 return res … … 1257 1250 1258 1251 #Add the x and y together 1259 vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)1252 vertex_coordinates = num.concatenate((x[:,num.NewAxis], y[:,num.NewAxis]),axis=1) 1260 1253 1261 1254 #Will return the quantity values at the specified times and locations -
anuga_core/source/anuga/fit_interpolate/search_functions.py
r5775 r6152 6 6 7 7 """ 8 from Numeric import dot9 8 import time 10 11 from Numeric import array12 9 13 10 from anuga.utilities.numerical_tools import get_machine_precision 14 11 from anuga.config import max_float 12 13 import Numeric as num 14 15 15 16 16 initial_search_value = 'uncomment search_functions code first'#0 … … 19 19 20 20 #FIXME test what happens if a 21 LAST_TRIANGLE = [[-10,[( array([max_float,max_float]),22 array([max_float,max_float]),23 array([max_float,max_float])),24 ( array([1,1]),array([0,0]),array([-1.1,-1.1]))]]]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]),num.array([0,0]),num.array([-1.1,-1.1]))]]] 25 25 26 26 def search_tree_of_vertices(root, mesh, x): … … 189 189 # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0) 190 190 191 sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)191 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) 192 192 if sigma0 < -epsilon: 193 193 return False,0,0,0 194 sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)194 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) 195 195 if sigma1 < -epsilon: 196 196 return False,0,0,0 197 sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)197 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) 198 198 if sigma2 < -epsilon: 199 199 return False,0,0,0 -
anuga_core/source/anuga/fit_interpolate/test_fit.py
r5349 r6152 12 12 import tempfile 13 13 import os 14 from Numeric import zeros, take, compress, Float, Int, dot, concatenate, \15 ArrayType, allclose, array16 14 17 15 from fit import * … … 23 21 from anuga.shallow_water import Domain 24 22 23 import Numeric as num 24 25 25 26 def distance(x, y): 26 return sqrt( sum( ( array(x)-array(y))**2 ))27 return sqrt( sum( (num.array(x)-num.array(y))**2 )) 27 28 28 29 def linear_function(point): 29 point = array(point)30 point = num.array(point) 30 31 return point[:,0]+point[:,1] 31 32 … … 64 65 #print "z",z 65 66 66 assert allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7)67 assert num.allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7) 67 68 68 69 f = fit.fit() … … 73 74 #print "answer\n",answer 74 75 75 assert allclose(f, answer, atol=1e-7)76 assert num.allclose(f, answer, atol=1e-7) 76 77 77 78 def test_smooth_att_to_meshII(self): … … 96 97 #print "answer\n",answer 97 98 98 assert allclose(f, answer)99 assert num.allclose(f, answer) 99 100 100 101 def test_smooth_attributes_to_meshIII(self): … … 132 133 #print "f\n",f 133 134 #print "answer\n",answer 134 assert allclose(f, answer)135 assert num.allclose(f, answer) 135 136 136 137 … … 158 159 f = fit.fit(data_coords,z) 159 160 answer = [[0,0], [5., 10.], [5., 10.]] 160 assert allclose(f, answer)161 assert num.allclose(f, answer) 161 162 162 163 def test_smooth_attributes_to_mesh_build_fit_subset(self): … … 204 205 #print "f\n",f 205 206 #print "answer\n",answer 206 assert allclose(f, answer)207 assert num.allclose(f, answer) 207 208 208 209 # test fit 2 mesh as well. … … 244 245 #print "f\n",f 245 246 #print "answer\n",answer 246 assert allclose(f, answer)247 assert num.allclose(f, answer) 247 248 os.remove(fileName) 248 249 … … 280 281 #print "f",f 281 282 #print "answer",answer 282 assert allclose(f, answer)283 assert num.allclose(f, answer) 283 284 284 285 # Delete file! … … 323 324 #print "f\n",f 324 325 #print "answer\n",answer 325 assert allclose(f, answer)326 assert num.allclose(f, answer) 326 327 os.remove(fileName) 327 328 … … 362 363 #print "f\n",f 363 364 #print "answer\n",answer 364 assert allclose(f, answer)365 assert num.allclose(f, answer) 365 366 os.remove(fileName) 366 367 os.remove(fileName_pts) … … 402 403 #print "f\n",f 403 404 #print "answer\n",answer 404 assert allclose(f, answer)405 assert num.allclose(f, answer) 405 406 os.remove(fileName) 406 407 os.remove(fileName_pts) … … 442 443 #print "f\n",f 443 444 #print "answer\n",answer 444 assert allclose(f, answer)445 assert num.allclose(f, answer) 445 446 446 447 os.remove(fileName) … … 483 484 #print "f\n",f 484 485 #print "answer\n",answer 485 assert allclose(f, answer)486 assert num.allclose(f, answer) 486 487 os.remove(fileName) 487 488 … … 524 525 #print "f",f 525 526 #print "answer",answer 526 assert allclose(f, answer)527 assert num.allclose(f, answer) 527 528 528 529 … … 562 563 f = interp.fit(data_points,z) 563 564 #f will be different from answer due to smoothing 564 assert allclose(f, answer,atol=5)565 assert num.allclose(f, answer,atol=5) 565 566 566 567 567 568 #Tests of smoothing matrix 568 569 def test_smoothing_matrix_one_triangle(self): 569 from Numeric import dot570 570 a = [0.0, 0.0] 571 571 b = [0.0, 2.0] … … 577 577 interp = Fit(points, vertices) 578 578 579 assert allclose(interp.get_D(), [[1, -0.5, -0.5],580 [-0.5, 0.5, 0],581 [-0.5, 0, 0.5]])579 assert num.allclose(interp.get_D(), [[1, -0.5, -0.5], 580 [-0.5, 0.5, 0], 581 [-0.5, 0, 0.5]]) 582 582 583 583 #Define f(x,y) = x 584 f = array([0,0,2]) #Value at global vertex 2584 f = num.array([0,0,2]) #Value at global vertex 2 585 585 586 586 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 587 587 # int 1 dx dy = area = 2 588 assert dot(dot(f, interp.get_D()), f) == 2588 assert num.dot(num.dot(f, interp.get_D()), f) == 2 589 589 590 590 #Define f(x,y) = y 591 f = array([0,2,0]) #Value at global vertex 1591 f = num.array([0,2,0]) #Value at global vertex 1 592 592 593 593 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 594 594 # int 1 dx dy = area = 2 595 assert dot(dot(f, interp.get_D()), f) == 2595 assert num.dot(num.dot(f, interp.get_D()), f) == 2 596 596 597 597 #Define f(x,y) = x+y 598 f = array([0,2,2]) #Values at global vertex 1 and 2598 f = num.array([0,2,2]) #Values at global vertex 1 and 2 599 599 600 600 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 601 601 # int 2 dx dy = 2*area = 4 602 assert dot(dot(f, interp.get_D()), f) == 4602 assert num.dot(num.dot(f, interp.get_D()), f) == 4 603 603 604 604 605 605 def test_smoothing_matrix_more_triangles(self): 606 from Numeric import dot607 608 606 a = [0.0, 0.0] 609 607 b = [0.0, 2.0] … … 625 623 626 624 #Define f(x,y) = x 627 f = array([0,0,2,0,2,4]) #f evaluated at points a-f625 f = num.array([0,0,2,0,2,4]) #f evaluated at points a-f 628 626 629 627 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 630 628 # int 1 dx dy = total area = 8 631 assert dot(dot(f, interp.get_D()), f) == 8629 assert num.dot(num.dot(f, interp.get_D()), f) == 8 632 630 633 631 #Define f(x,y) = y 634 f = array([0,2,0,4,2,0]) #f evaluated at points a-f632 f = num.array([0,2,0,4,2,0]) #f evaluated at points a-f 635 633 636 634 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 637 635 # int 1 dx dy = area = 8 638 assert dot(dot(f, interp.get_D()), f) == 8636 assert num.dot(num.dot(f, interp.get_D()), f) == 8 639 637 640 638 #Define f(x,y) = x+y 641 f = array([0,2,2,4,4,4]) #f evaluated at points a-f639 f = num.array([0,2,2,4,4,4]) #f evaluated at points a-f 642 640 643 641 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 644 642 # int 2 dx dy = 2*area = 16 645 assert dot(dot(f, interp.get_D()), f) == 16643 assert num.dot(num.dot(f, interp.get_D()), f) == 16 646 644 647 645 … … 709 707 f1 = interp.fit(data_points1,z) 710 708 711 assert allclose(f,f1), 'Fit should have been unaltered'709 assert num.allclose(f,f1), 'Fit should have been unaltered' 712 710 713 711 … … 736 734 answer = [[0, 0], [5., 10.], [5., 10.]] 737 735 738 assert allclose(f, answer)736 assert num.allclose(f, answer) 739 737 740 738 def test_fit_to_mesh_w_georef(self): … … 773 771 mesh_origin = mesh_geo.get_origin(), 774 772 alpha = 0) 775 assert allclose( zz, [0,5,5] )773 assert num.allclose( zz, [0,5,5] ) 776 774 777 775 … … 819 817 [5.0, 10.0], 820 818 [5.0,10.0]] 821 assert allclose(mesh_dic['vertex_attributes'],ans)819 assert num.allclose(mesh_dic['vertex_attributes'],ans) 822 820 823 821 self.failUnless(mesh_dic['vertex_attribute_titles'] == … … 827 825 828 826 answer = [0., 5., 5.] 829 assert allclose(domain.quantities['elevation'].vertex_values,830 answer)827 assert num.allclose(domain.quantities['elevation'].vertex_values, 828 answer) 831 829 #clean up 832 830 os.remove(mesh_file) … … 878 876 [5.0, 10.0], 879 877 [5.0,10.0]] 880 assert allclose(mesh_dic['vertex_attributes'],ans)878 assert num.allclose(mesh_dic['vertex_attributes'],ans) 881 879 882 880 self.failUnless(mesh_dic['vertex_attribute_titles'] == … … 929 927 mesh_dic = import_mesh_file(mesh_output_file) 930 928 931 assert allclose(mesh_dic['vertex_attributes'],932 [[1.0, 2.0,0.0, 0.0],933 [1.0, 2.0,5.0, 10.0],934 [1.0, 2.0,5.0,10.0]])929 assert num.allclose(mesh_dic['vertex_attributes'], 930 [[1.0, 2.0,0.0, 0.0], 931 [1.0, 2.0,5.0, 10.0], 932 [1.0, 2.0,5.0,10.0]]) 935 933 936 934 self.failUnless(mesh_dic['vertex_attribute_titles'] == -
anuga_core/source/anuga/fit_interpolate/test_interpolate.py
r5775 r6152 14 14 15 15 from Scientific.IO.NetCDF import NetCDFFile 16 from Numeric import allclose, array, transpose, zeros, Float, sometrue, \ 17 alltrue, take, where 16 17 import Numeric as num 18 18 19 19 20 … … 28 29 29 30 def distance(x, y): 30 return sqrt( sum( (array(x)-array(y))**2 ))31 return sqrt( num.sum( (num.array(x)-num.array(y))**2 )) 31 32 32 33 def linear_function(point): 33 point = array(point)34 point = num.array(point) 34 35 return point[:,0]+point[:,1] 35 36 … … 66 67 67 68 bed = domain.quantities['elevation'].vertex_values 68 stage = zeros(bed.shape,Float)69 stage = num.zeros(bed.shape, num.Float) 69 70 70 71 h = 0.3 … … 104 105 interp = Interpolate(points, vertices) 105 106 A, _, _ = interp._build_interpolation_matrix_A(data) 106 assert allclose(A.todense(), 107 [[1./3, 1./3, 1./3]]) 107 assert num.allclose(A.todense(), [[1./3, 1./3, 1./3]]) 108 108 109 109 … … 113 113 from mesh_factory import rectangular 114 114 from shallow_water import Domain 115 from Numeric import zeros, Float116 115 from abstract_2d_finite_volumes.quantity import Quantity 117 116 … … 130 129 131 130 x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False) 132 vertex_coordinates = concatenate( (x[:, NewAxis], y[:,NewAxis]), axis=1 )131 vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 ) 133 132 # FIXME: This concat should roll into get_vertex_values 134 133 … … 140 139 I = Interpolate(vertex_coordinates, triangles) 141 140 result = I.interpolate(vertex_values, interpolation_points) 142 assert allclose(result, answer)141 assert num.allclose(result, answer) 143 142 144 143 … … 150 149 151 150 x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False) 152 vertex_coordinates = concatenate( (x[:, NewAxis], y[:,NewAxis]), axis=1 )151 vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 ) 153 152 # FIXME: This concat should roll into get_vertex_values 154 153 … … 160 159 I = Interpolate(vertex_coordinates, triangles) 161 160 result = I.interpolate(vertex_values, interpolation_points) 162 assert allclose(result, answer)161 assert num.allclose(result, answer) 163 162 164 163 … … 167 166 from mesh_factory import rectangular 168 167 from shallow_water import Domain 169 from Numeric import zeros, Float170 168 from abstract_2d_finite_volumes.quantity import Quantity 171 169 … … 184 182 185 183 x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False) 186 vertex_coordinates = concatenate( (x[:, NewAxis], y[:,NewAxis]), axis=1 )184 vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 ) 187 185 # FIXME: This concat should roll into get_vertex_values 188 186 … … 193 191 194 192 result = interpolate(vertex_coordinates, triangles, vertex_values, interpolation_points) 195 assert allclose(result, answer)193 assert num.allclose(result, answer) 196 194 197 195 … … 203 201 204 202 x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False) 205 vertex_coordinates = concatenate( (x[:, NewAxis], y[:,NewAxis]), axis=1 )203 vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 ) 206 204 # FIXME: This concat should roll into get_vertex_values 207 205 … … 213 211 result = interpolate(vertex_coordinates, triangles, 214 212 vertex_values, interpolation_points) 215 assert allclose(result, answer)213 assert num.allclose(result, answer) 216 214 217 215 … … 220 218 from mesh_factory import rectangular 221 219 from shallow_water import Domain 222 from Numeric import zeros, Float223 220 from abstract_2d_finite_volumes.quantity import Quantity 224 221 … … 236 233 237 234 x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False) 238 vertex_coordinates = concatenate( (x[:, NewAxis], y[:,NewAxis]), axis=1 )235 vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 ) 239 236 # FIXME: This concat should roll into get_vertex_values 240 237 … … 248 245 use_cache=True, 249 246 verbose=False) 250 assert allclose(result, answer)247 assert num.allclose(result, answer) 251 248 252 249 # Second call using the cache … … 255 252 use_cache=True, 256 253 verbose=False) 257 assert allclose(result, answer)254 assert num.allclose(result, answer) 258 255 259 256 … … 289 286 290 287 A,_,_ = interp._build_interpolation_matrix_A(data) 291 assert allclose(A.todense(), answer)288 assert num.allclose(A.todense(), answer) 292 289 293 290 #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh … … 298 295 299 296 A,_,_ = interp._build_interpolation_matrix_A(data) 300 assert allclose(A.todense(), answer)297 assert num.allclose(A.todense(), answer) 301 298 302 299 … … 309 306 310 307 A,_,_ = interp._build_interpolation_matrix_A(data) 311 assert allclose(A.todense(), answer)308 assert num.allclose(A.todense(), answer) 312 309 313 310 … … 330 327 331 328 A,_,_ = interp._build_interpolation_matrix_A(data) 332 assert allclose(A.todense(), answer)329 assert num.allclose(A.todense(), answer) 333 330 334 331 … … 351 348 352 349 A,_,_ = interp._build_interpolation_matrix_A(data) 353 assert allclose(A.todense(), answer)350 assert num.allclose(A.todense(), answer) 354 351 355 352 def test_datapoints_on_edges(self): … … 372 369 373 370 A,_,_ = interp._build_interpolation_matrix_A(data) 374 assert allclose(A.todense(), answer)371 assert num.allclose(A.todense(), answer) 375 372 376 373 … … 378 375 #Try arbitrary datapoints 379 376 380 381 from Numeric import sum382 377 383 378 a = [0.0, 0.0] … … 394 389 A,_,_ = interp._build_interpolation_matrix_A(data) 395 390 results = A.todense() 396 assert allclose(sum(results, axis=1), 1.0)391 assert num.allclose(num.sum(results, axis=1), 1.0) 397 392 398 393 def test_arbitrary_datapoints_some_outside(self): … … 400 395 #That one should be ignored 401 396 402 403 from Numeric import sum404 397 405 398 a = [0.0, 0.0] … … 415 408 A,_,_ = interp._build_interpolation_matrix_A(data) 416 409 results = A.todense() 417 assert allclose(sum(results, axis=1), [1,1,1,0])410 assert num.allclose(num.sum(results, axis=1), [1,1,1,0]) 418 411 419 412 … … 448 441 for i in range(A.shape[0]): 449 442 for j in range(A.shape[1]): 450 if not allclose(A[i,j], answer[i][j]):443 if not num.allclose(A[i,j], answer[i][j]): 451 444 print i,j,':',A[i,j], answer[i][j] 452 445 … … 454 447 #results = interp._build_interpolation_matrix_A(data).todense() 455 448 456 assert allclose(A, answer)449 assert num.allclose(A, answer) 457 450 458 451 def test_geo_ref(self): … … 481 474 #print "z",z 482 475 #print "answer",answer 483 assert allclose(z, answer)476 assert num.allclose(z, answer) 484 477 485 478 … … 489 482 #print "z",z 490 483 #print "answer",answer 491 assert allclose(z, answer)484 assert num.allclose(z, answer) 492 485 493 486 … … 516 509 #print "z",z 517 510 #print "answer",answer 518 assert allclose(z, answer)511 assert num.allclose(z, answer) 519 512 520 513 … … 524 517 #print "z",z 525 518 #print "answer",answer 526 assert allclose(z, answer)519 assert num.allclose(z, answer) 527 520 528 521 … … 552 545 #print "z",z 553 546 #print "answer",answer 554 assert allclose(z, answer)547 assert num.allclose(z, answer) 555 548 556 549 z = interp.interpolate(f, point_coords, start_blocking_len = 2) … … 559 552 #print "z",z 560 553 #print "answer",answer 561 assert allclose(z, answer)554 assert num.allclose(z, answer) 562 555 563 556 def test_interpolate_attributes_to_points(self): … … 581 574 #print "z",z 582 575 #print "answer",answer 583 assert allclose(z, answer)576 assert num.allclose(z, answer) 584 577 585 578 … … 589 582 #print "z",z 590 583 #print "answer",answer 591 assert allclose(z, answer)584 assert num.allclose(z, answer) 592 585 593 586 def test_interpolate_attributes_to_pointsII(self): … … 622 615 #print "z",z 623 616 #print "answer",answer 624 assert allclose(z, answer)617 assert num.allclose(z, answer) 625 618 626 619 z = interp.interpolate(f, point_coords, start_blocking_len = 2) … … 629 622 #print "z",z 630 623 #print "answer",answer 631 assert allclose(z, answer)624 assert num.allclose(z, answer) 632 625 633 626 def test_interpolate_attributes_to_pointsIII(self): … … 683 676 #print "***********" 684 677 685 assert allclose(z, answer)678 assert num.allclose(z, answer) 686 679 687 680 … … 690 683 #print "z",z 691 684 #print "answer",answer 692 assert allclose(z, answer)685 assert num.allclose(z, answer) 693 686 694 687 def test_interpolate_point_outside_of_mesh(self): … … 718 711 719 712 z = interp.interpolate(f, point_coords) #, verbose=True) 720 answer = array([ [NAN, NAN, NAN, NAN]]) # (-1,-1)713 answer = num.array([ [NAN, NAN, NAN, NAN]]) # (-1,-1) 721 714 722 715 #print "***********" … … 767 760 768 761 interp = Interpolate(vertices, triangles) 769 f = array([linear_function(vertices),2*linear_function(vertices) ])770 f = transpose(f)762 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 763 f = num.transpose(f) 771 764 #print "f",f 772 765 z = interp.interpolate(f, point_coords) 773 766 answer = [linear_function(point_coords), 774 767 2*linear_function(point_coords) ] 775 answer = transpose(answer)768 answer = num.transpose(answer) 776 769 #print "z",z 777 770 #print "answer",answer 778 assert allclose(z, answer)771 assert num.allclose(z, answer) 779 772 780 773 z = interp.interpolate(f, point_coords, start_blocking_len = 2) … … 782 775 #print "z",z 783 776 #print "answer",answer 784 assert allclose(z, answer)777 assert num.allclose(z, answer) 785 778 786 779 def test_interpolate_blocking(self): … … 810 803 811 804 interp = Interpolate(vertices, triangles) 812 f = array([linear_function(vertices),2*linear_function(vertices) ])813 f = transpose(f)805 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 806 f = num.transpose(f) 814 807 #print "f",f 815 808 for blocking_max in range(len(point_coords)+2): … … 820 813 answer = [linear_function(point_coords), 821 814 2*linear_function(point_coords) ] 822 answer = transpose(answer)815 answer = num.transpose(answer) 823 816 #print "z",z 824 817 #print "answer",answer 825 assert allclose(z, answer)818 assert num.allclose(z, answer) 826 819 827 f = array([linear_function(vertices),2*linear_function(vertices),828 2*linear_function(vertices) - 100 ])829 f = transpose(f)820 f = num.array([linear_function(vertices),2*linear_function(vertices), 821 2*linear_function(vertices) - 100 ]) 822 f = num.transpose(f) 830 823 #print "f",f 831 824 for blocking_max in range(len(point_coords)+2): … … 834 827 z = interp.interpolate(f, point_coords, 835 828 start_blocking_len=blocking_max) 836 answer = array([linear_function(point_coords),837 2*linear_function(point_coords) ,838 2*linear_function(point_coords)-100 ])839 z = transpose(z)829 answer = num.array([linear_function(point_coords), 830 2*linear_function(point_coords) , 831 2*linear_function(point_coords)-100 ]) 832 z = num.transpose(z) 840 833 #print "z",z 841 834 #print "answer",answer 842 assert allclose(z, answer)835 assert num.allclose(z, answer) 843 836 844 837 def test_interpolate_geo_spatial(self): … … 872 865 873 866 interp = Interpolate(vertices, triangles) 874 f = array([linear_function(vertices),2*linear_function(vertices) ])875 f = transpose(f)867 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 868 f = num.transpose(f) 876 869 #print "f",f 877 870 for blocking_max in range(14): … … 884 877 2*linear_function(point_coords.get_data_points( \ 885 878 absolute = True)) ] 886 answer = transpose(answer)879 answer = num.transpose(answer) 887 880 #print "z",z 888 881 #print "answer",answer 889 assert allclose(z, answer)882 assert num.allclose(z, answer) 890 883 891 f = array([linear_function(vertices),2*linear_function(vertices),892 2*linear_function(vertices) - 100 ])893 f = transpose(f)884 f = num.array([linear_function(vertices),2*linear_function(vertices), 885 2*linear_function(vertices) - 100 ]) 886 f = num.transpose(f) 894 887 #print "f",f 895 888 for blocking_max in range(14): … … 898 891 z = interp.interpolate(f, point_coords, 899 892 start_blocking_len=blocking_max) 900 answer = array([linear_function(point_coords.get_data_points( \901 absolute = True)),902 2*linear_function(point_coords.get_data_points( \903 absolute = True)) ,904 2*linear_function(point_coords.get_data_points( \905 absolute = True))-100 ])906 z = transpose(z)893 answer = num.array([linear_function(point_coords.get_data_points( \ 894 absolute = True)), 895 2*linear_function(point_coords.get_data_points( \ 896 absolute = True)) , 897 2*linear_function(point_coords.get_data_points( \ 898 absolute = True))-100 ]) 899 z = num.transpose(z) 907 900 #print "z",z 908 901 #print "answer",answer 909 assert allclose(z, answer)902 assert num.allclose(z, answer) 910 903 911 904 z = interp.interpolate(f, point_coords, start_blocking_len = 2) … … 913 906 #print "z",z 914 907 #print "answer",answer 915 assert allclose(z, answer)908 assert num.allclose(z, answer) 916 909 917 910 def test_interpolate_geo_spatial(self): … … 945 938 946 939 interp = Interpolate(vertices, triangles) 947 f = array([linear_function(vertices),2*linear_function(vertices) ])948 f = transpose(f)940 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 941 f = num.transpose(f) 949 942 #print "f",f 950 943 z = interp.interpolate_block(f, point_coords) … … 953 946 2*linear_function(point_coords.get_data_points( \ 954 947 absolute = True)) ] 955 answer = transpose(answer)948 answer = num.transpose(answer) 956 949 #print "z",z 957 950 #print "answer",answer 958 assert allclose(z, answer)951 assert num.allclose(z, answer) 959 952 960 953 z = interp.interpolate(f, point_coords, start_blocking_len = 2) … … 962 955 #print "z",z 963 956 #print "answer",answer 964 assert allclose(z, answer)957 assert num.allclose(z, answer) 965 958 966 959 … … 991 984 992 985 interp = Interpolate(vertices, triangles) 993 f = array([linear_function(vertices),2*linear_function(vertices) ])994 f = transpose(f)986 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 987 f = num.transpose(f) 995 988 z = interp.interpolate(f, point_coords, 996 989 start_blocking_len=20) 997 990 answer = [linear_function(point_coords), 998 991 2*linear_function(point_coords) ] 999 answer = transpose(answer)992 answer = num.transpose(answer) 1000 993 #print "z",z 1001 994 #print "answer",answer 1002 assert allclose(z, answer)1003 assert allclose(interp._A_can_be_reused, True)995 assert num.allclose(z, answer) 996 assert num.allclose(interp._A_can_be_reused, True) 1004 997 1005 998 z = interp.interpolate(f) 1006 assert allclose(z, answer)999 assert num.allclose(z, answer) 1007 1000 1008 1001 # This causes blocking to occur. 1009 1002 z = interp.interpolate(f, start_blocking_len=10) 1010 assert allclose(z, answer)1011 assert allclose(interp._A_can_be_reused, False)1003 assert num.allclose(z, answer) 1004 assert num.allclose(interp._A_can_be_reused, False) 1012 1005 1013 1006 #A is recalculated 1014 1007 z = interp.interpolate(f) 1015 assert allclose(z, answer)1016 assert allclose(interp._A_can_be_reused, True)1008 assert num.allclose(z, answer) 1009 assert num.allclose(interp._A_can_be_reused, True) 1017 1010 1018 1011 interp = Interpolate(vertices, triangles) … … 1054 1047 1055 1048 interp = Interpolate(vertices, triangles) 1056 f = array([linear_function(vertices),2*linear_function(vertices) ])1057 f = transpose(f)1049 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 1050 f = num.transpose(f) 1058 1051 z = interp.interpolate(f, point_coords) 1059 1052 answer = [linear_function(point_coords), 1060 1053 2*linear_function(point_coords) ] 1061 answer = transpose(answer)1062 1063 assert allclose(z, answer)1064 assert allclose(interp._A_can_be_reused, True)1054 answer = num.transpose(answer) 1055 1056 assert num.allclose(z, answer) 1057 assert num.allclose(interp._A_can_be_reused, True) 1065 1058 1066 1059 1067 1060 z = interp.interpolate(f) # None 1068 assert allclose(z, answer)1061 assert num.allclose(z, answer) 1069 1062 z = interp.interpolate(f, point_coords) # Repeated (not really a test) 1070 assert allclose(z, answer)1063 assert num.allclose(z, answer) 1071 1064 1072 1065 … … 1085 1078 1086 1079 #One quantity 1087 Q = zeros( (3,6),Float )1080 Q = num.zeros( (3,6), num.Float ) 1088 1081 1089 1082 #Linear in time and space … … 1111 1104 #Check temporal interpolation 1112 1105 for i in [0,1,2]: 1113 assert allclose(I(time[i]), mean(Q[i,:]))1106 assert num.allclose(I(time[i]), mean(Q[i,:])) 1114 1107 1115 1108 #Midway 1116 assert allclose(I( (time[0] + time[1])/2 ),1117 (I(time[0]) + I(time[1]))/2 )1118 1119 assert allclose(I( (time[1] + time[2])/2 ),1120 (I(time[1]) + I(time[2]))/2 )1121 1122 assert allclose(I( (time[0] + time[2])/2 ),1123 (I(time[0]) + I(time[2]))/2 )1109 assert num.allclose(I( (time[0] + time[1])/2 ), 1110 (I(time[0]) + I(time[1]))/2 ) 1111 1112 assert num.allclose(I( (time[1] + time[2])/2 ), 1113 (I(time[1]) + I(time[2]))/2 ) 1114 1115 assert num.allclose(I( (time[0] + time[2])/2 ), 1116 (I(time[0]) + I(time[2]))/2 ) 1124 1117 1125 1118 #1/3 1126 assert allclose(I( (time[0] + time[2])/3 ),1127 (I(time[0]) + I(time[2]))/3 )1119 assert num.allclose(I( (time[0] + time[2])/3 ), 1120 (I(time[0]) + I(time[2]))/3 ) 1128 1121 1129 1122 … … 1191 1184 for j in range(50): #t in [1, 6] 1192 1185 for id in range(len(interpolation_points)): 1193 assert allclose(I(t, id), answer[id])1186 assert num.allclose(I(t, id), answer[id]) 1194 1187 t += 0.1 1195 1188 … … 1232 1225 1233 1226 #One quantity 1234 Q = zeros( (3,6),Float )1227 Q = num.zeros( (3,6), num.Float ) 1235 1228 1236 1229 #Linear in time and space … … 1250 1243 for j in range(50): #t in [1, 6] 1251 1244 for id in range(len(interpolation_points)): 1252 assert allclose(I(t, id), t*answer[id])1245 assert num.allclose(I(t, id), t*answer[id]) 1253 1246 t += 0.1 1254 1247 … … 1292 1285 1293 1286 #One quantity 1294 Q = zeros( (8,6),Float )1287 Q = num.zeros( (8,6), num.Float ) 1295 1288 1296 1289 #Linear in time and space … … 1312 1305 for j in range(50): #t in [1, 6] 1313 1306 for id in range(len(interpolation_points)): 1314 assert allclose(I(t, id), t*answer[id])1307 assert num.allclose(I(t, id), t*answer[id]) 1315 1308 t += 0.1 1316 1309 … … 1326 1319 1327 1320 assert len(I.time) == 4 1328 assert( allclose(I.time, [1.0, 4.0, 7.0, 9.0] ))1321 assert( num.allclose(I.time, [1.0, 4.0, 7.0, 9.0] )) 1329 1322 1330 1323 answer = linear_function(interpolation_points) … … 1333 1326 for j in range(50): #t in [1, 6] 1334 1327 for id in range(len(interpolation_points)): 1335 assert allclose(I(t, id), t*answer[id])1328 assert num.allclose(I(t, id), t*answer[id]) 1336 1329 t += 0.1 1337 1330 … … 1361 1354 1362 1355 #One quantity 1363 Q = zeros( (2,6),Float )1356 Q = num.zeros( (2,6), num.Float ) 1364 1357 1365 1358 #Linear in time and space … … 1371 1364 1372 1365 interp = Interpolate(points, triangles) 1373 f = array([linear_function(points),2*linear_function(points) ])1374 f = transpose(f)1366 f = num.array([linear_function(points),2*linear_function(points) ]) 1367 f = num.transpose(f) 1375 1368 #print "f",f 1376 1369 z = interp.interpolate(f, interpolation_points) 1377 1370 answer = [linear_function(interpolation_points), 1378 1371 2*linear_function(interpolation_points) ] 1379 answer = transpose(answer)1372 answer = num.transpose(answer) 1380 1373 #print "z",z 1381 1374 #print "answer",answer 1382 assert allclose(z, answer)1375 assert num.allclose(z, answer) 1383 1376 1384 1377 … … 1393 1386 1394 1387 msg = 'Interpolation failed' 1395 assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg1388 assert num.allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg 1396 1389 #self.failUnless( I.precomputed_values['Attribute'][1] == 60.0, 1397 1390 # ' failed') … … 1427 1420 1428 1421 # One quantity 1429 Q = zeros( (3,6),Float )1422 Q = num.zeros( (3,6), num.Float ) 1430 1423 1431 1424 # Linear in time and space … … 1442 1435 1443 1436 1444 assert alltrue(I.precomputed_values['Attribute'][:,4] != NAN)1445 assert sometrue(I.precomputed_values['Attribute'][:,5] == NAN)1437 assert num.alltrue(I.precomputed_values['Attribute'][:,4] != NAN) 1438 assert num.sometrue(I.precomputed_values['Attribute'][:,5] == NAN) 1446 1439 1447 1440 #X = I.precomputed_values['Attribute'][1,:] … … 1455 1448 for j in range(50): #t in [1, 6] 1456 1449 for id in range(len(interpolation_points)-1): 1457 assert allclose(I(t, id), t*answer[id])1450 assert num.allclose(I(t, id), t*answer[id]) 1458 1451 t += 0.1 1459 1452 … … 1477 1470 1478 1471 1479 time = array(\1472 time = num.array(\ 1480 1473 [0.00000000e+00, 5.00000000e-02, 1.00000000e-01, 1.50000000e-01, 1481 1474 2.00000000e-01, 2.50000000e-01, 3.00000000e-01, 3.50000000e-01, … … 1616 1609 1617 1610 #One quantity 1618 Q = zeros( (len(time),6),Float )1611 Q = num.zeros( (len(time),6), num.Float ) 1619 1612 1620 1613 #Linear in time and space … … 1661 1654 1662 1655 interp = Interpolate(vertices, triangles) 1663 f = array([linear_function(vertices),2*linear_function(vertices) ])1664 f = transpose(f)1656 f = num.array([linear_function(vertices),2*linear_function(vertices) ]) 1657 f = num.transpose(f) 1665 1658 #print "f",f 1666 1659 z = interp.interpolate(f, geo_data) … … 1668 1661 answer = [linear_function(point_coords), 1669 1662 2*linear_function(point_coords) ] 1670 answer = transpose(answer)1663 answer = num.transpose(answer) 1671 1664 answer[2,:] = [NAN, NAN] 1672 1665 answer[3,:] = [NAN, NAN] … … 1674 1667 #print "z",z 1675 1668 #print "answer _ fixed",answer 1676 assert allclose(z[0:1], answer[0:1])1677 assert allclose(z[4:10], answer[4:10])1669 assert num.allclose(z[0:1], answer[0:1]) 1670 assert num.allclose(z[4:10], answer[4:10]) 1678 1671 for i in [2,3,11]: 1679 1672 self.failUnless( z[i,1] == answer[11,1], 'Fail!') … … 1765 1758 #print "velocity_answers_array", velocity_answers[i] 1766 1759 msg = 'Interpolation failed' 1767 assert allclose(float(depths[i]), depth_answers[i]), msg1768 assert allclose(float(velocitys[i]), velocity_answers[i]), msg1760 assert num.allclose(float(depths[i]), depth_answers[i]), msg 1761 assert num.allclose(float(velocitys[i]), velocity_answers[i]), msg 1769 1762 1770 1763 velocity_y_file_handle = file(velocity_y_file) … … 1780 1773 #print "velocity_answers_array", velocity_answers[i] 1781 1774 msg = 'Interpolation failed' 1782 assert allclose(float(depths[i]), depth_answers[i]), msg1783 assert allclose(float(velocitys[i]), velocity_answers[i]), msg1775 assert num.allclose(float(depths[i]), depth_answers[i]), msg 1776 assert num.allclose(float(velocitys[i]), velocity_answers[i]), msg 1784 1777 1785 1778 # clean up … … 1852 1845 #print "z",z 1853 1846 #print "answer",answer 1854 assert allclose(z, answer)1847 assert num.allclose(z, answer) 1855 1848 1856 1849 #------------------------------------------------------------- -
anuga_core/source/anuga/fit_interpolate/test_search_functions.py
r4859 r6152 5 5 from search_functions import search_tree_of_vertices, set_last_triangle 6 6 from search_functions import _search_triangles_of_vertices 7 8 9 from Numeric import zeros, array, allclose10 7 11 8 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
Note: See TracChangeset
for help on using the changeset viewer.