#!/usr/bin/env python #TEST import unittest from math import sqrt from least_squares import * from Numeric import allclose, array, transpose def distance(x, y): return sqrt( sum( (array(x)-array(y))**2 )) def linear_function(point): point = array(point) return point[:,0]+point[:,1] class TestCase(unittest.TestCase): def setUp(self): pass def tearDown(self): pass def test_datapoint_at_centroid(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point interp = Interpolation(points, vertices, data) assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) def test_quad_tree(self): p0 = [-10.0, -10.0] p1 = [20.0, -10.0] p2 = [-10.0, 20.0] p3 = [10.0, 50.0] p4 = [30.0, 30.0] p5 = [50.0, 10.0] p6 = [40.0, 60.0] p7 = [60.0, 40.0] p8 = [-66.0, 20.0] p9 = [10.0, -66.0] points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9] triangles = [ [0, 1, 2], [3, 2, 4], [4, 2, 1], [4, 1, 5], [3, 4, 6], [6, 4, 7], [7, 4, 5], [8, 0, 2], [0, 9, 1]] data = [ [4,4] ] interp = Interpolation(points, triangles, data, alpha = 0.0) #print "PDSG - interp.get_A()", interp.get_A() answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0., 0., 0. , 0., 0., 0., 0.]] assert allclose(interp.get_A(), answer) interp.set_point_coordinates([[-30, -30]]) #point outside of mesh #print "PDSG - interp.get_A()", interp.get_A() answer = [ [ 0.0, 0.0, 0.0, 0., 0., 0. , 0., 0., 0., 0.]] assert allclose(interp.get_A(), answer) #point outside of quad tree root cell interp.set_point_coordinates([[-70, -70]]) #print "PDSG - interp.get_A()", interp.get_A() answer = [ [ 0.0, 0.0, 0.0, 0., 0., 0. , 0., 0., 0., 0.]] assert allclose(interp.get_A(), answer) def test_datapoints_at_vertices(self): """Test that data points coinciding with vertices yield a diagonal matrix """ a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac data = points #Use data at vertices interp = Interpolation(points, vertices, data) assert allclose(interp.get_A(), [[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) def test_datapoints_on_edge_midpoints(self): """Try datapoints midway on edges - each point should affect two matrix entries equally """ a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac data = [ [0., 1.], [1., 0.], [1., 1.] ] interp = Interpolation(points, vertices, data) assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 [0.5, 0.0, 0.5], #Affects vertex 0 and 2 [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2 def test_datapoints_on_edges(self): """Try datapoints on edges - each point should affect two matrix entries in proportion """ a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ] interp = Interpolation(points, vertices, data) assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0 [0.25, 0.0, 0.75], #Affects vertex 0 and 2 [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2 def test_arbitrary_datapoints(self): """Try arbitrary datapoints """ from Numeric import sum a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ] interp = Interpolation(points, vertices, data) assert allclose(sum(interp.get_A(), axis=1), 1.0) # this causes a memory error in scipy.sparce def test_more_triangles(self): a = [-1.0, 0.0] b = [3.0, 4.0] c = [4.0,1.0] d = [-3.0, 2.0] #3 e = [-1.0,-2.0] f = [1.0, -2.0] #5 points = [a, b, c, d,e,f] triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc #Data points data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ] interp = Interpolation(points, triangles, data_points) answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d [0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d [0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b [0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f #FIXME: This test failed on our Linux box, but not on Windows. A = interp.get_A() for i in range(A.shape[0]): for j in range(A.shape[1]): if not allclose(A[i,j], answer[i][j]): print i,j,':',A[i,j], answer[i][j] assert allclose(interp.get_A(), answer) def test_smooth_attributes_to_mesh(self): a = [0.0, 0.0] b = [0.0, 5.0] c = [5.0, 0.0] points = [a, b, c] triangles = [ [1,0,2] ] #bac d1 = [1.0, 1.0] d2 = [1.0, 3.0] d3 = [3.0,1.0] z1 = 2 z2 = 4 z3 = 4 data_coords = [d1, d2, d3] interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20) z = [z1, z2, z3] f = interp.fit(z) answer = [0, 5., 5.] #print "f\n",f #print "answer\n",answer assert allclose(f, answer, atol=1e-7) def test_smooth_att_to_meshII(self): a = [0.0, 0.0] b = [0.0, 5.0] c = [5.0, 0.0] points = [a, b, c] triangles = [ [1,0,2] ] #bac d1 = [1.0, 1.0] d2 = [1.0, 2.0] d3 = [3.0,1.0] data_coords = [d1, d2, d3] z = linear_function(data_coords) interp = Interpolation(points, triangles, data_coords, alpha=0.0) f = interp.fit(z) answer = linear_function(points) assert allclose(f, answer) def test_smooth_attributes_to_meshIII(self): a = [-1.0, 0.0] b = [3.0, 4.0] c = [4.0,1.0] d = [-3.0, 2.0] #3 e = [-1.0,-2.0] f = [1.0, -2.0] #5 vertices = [a, b, c, d,e,f] triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc point_coords = [[-2.0, 2.0], [-1.0, 1.0], [0.0,2.0], [1.0, 1.0], [2.0, 1.0], [0.0,0.0], [1.0, 0.0], [0.0, -1.0], [-0.2,-0.5], [-0.9, -1.5], [0.5, -1.9], [3.0,1.0]] z = linear_function(point_coords) interp = Interpolation(vertices, triangles, point_coords, alpha=0.0) #print 'z',z f = interp.fit(z) answer = linear_function(vertices) #print "f\n",f #print "answer\n",answer assert allclose(f, answer) def test_smooth_attributes_to_meshIV(self): """ Testing 2 attributes smoothed to the mesh """ a = [0.0, 0.0] b = [0.0, 5.0] c = [5.0, 0.0] points = [a, b, c] triangles = [ [1,0,2] ] #bac d1 = [1.0, 1.0] d2 = [1.0, 3.0] d3 = [3.0, 1.0] z1 = [2, 4] z2 = [4, 8] z3 = [4, 8] data_coords = [d1, d2, d3] interp = Interpolation(points, triangles, data_coords, alpha=0.0) z = [z1, z2, z3] f = interp.fit_points(z) answer = [[0,0], [5., 10.], [5., 10.]] assert allclose(f, answer) def test_interpolate_attributes_to_points(self): v0 = [0.0, 0.0] v1 = [0.0, 5.0] v2 = [5.0, 0.0] vertices = [v0, v1, v2] triangles = [ [1,0,2] ] #bac d0 = [1.0, 1.0] d1 = [1.0, 2.0] d2 = [3.0, 1.0] point_coords = [ d0, d1, d2] interp = Interpolation(vertices, triangles, point_coords) f = linear_function(vertices) z = interp.interpolate(f) answer = linear_function(point_coords) assert allclose(z, answer) def test_interpolate_attributes_to_pointsII(self): a = [-1.0, 0.0] b = [3.0, 4.0] c = [4.0, 1.0] d = [-3.0, 2.0] #3 e = [-1.0, -2.0] f = [1.0, -2.0] #5 vertices = [a, b, c, d,e,f] triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc point_coords = [[-2.0, 2.0], [-1.0, 1.0], [0.0, 2.0], [1.0, 1.0], [2.0, 1.0], [0.0, 0.0], [1.0, 0.0], [0.0, -1.0], [-0.2, -0.5], [-0.9, -1.5], [0.5, -1.9], [3.0, 1.0]] interp = Interpolation(vertices, triangles, point_coords) f = linear_function(vertices) z = interp.interpolate(f) answer = linear_function(point_coords) #print "z",z #print "answer",answer assert allclose(z, answer) def test_interpolate_attributes_to_pointsIII(self): v0 = [0.0, 0.0] v1 = [0.0, 5.0] v2 = [5.0, 0.0] vertices = [v0, v1, v2] triangles = [ [1,0,2] ] #bac d0 = [1.0, 1.0] d1 = [1.0, 2.0] d2 = [3.0,1.0] point_coords = [ d0, d1, d2] interp = Interpolation(vertices, triangles, point_coords) f = [ [0., 0.], [5., 10.], [5., 10.]] #linear_function(vertices) #print "f",f z = interp.interpolate(f) answer = [ [2., 4.], [3., 6.], [4., 8.]] #print "***********" #print "z",z #print "answer",answer #print "***********" assert allclose(z, answer) def test_interpolate_attributes_to_pointsIV(self): a = [-1.0, 0.0] b = [3.0, 4.0] c = [4.0, 1.0] d = [-3.0, 2.0] #3 e = [-1.0, -2.0] f = [1.0, -2.0] #5 vertices = [a, b, c, d,e,f] triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc point_coords = [[-2.0, 2.0], [-1.0, 1.0], [0.0, 2.0], [1.0, 1.0], [2.0, 1.0], [0.0, 0.0], [1.0, 0.0], [0.0, -1.0], [-0.2, -0.5], [-0.9, -1.5], [0.5, -1.9], [3.0, 1.0]] interp = Interpolation(vertices, triangles, point_coords) f = array([linear_function(vertices),2*linear_function(vertices) ]) f = transpose(f) #print "f",f z = interp.interpolate(f) answer = [linear_function(point_coords), 2*linear_function(point_coords) ] answer = transpose(answer) #print "z",z #print "answer",answer assert allclose(z, answer) def test_smooth_attributes_to_mesh_function(self): """ Testing 2 attributes smoothed to the mesh """ """Test multiple attributes """ a = [0.0, 0.0] b = [0.0, 5.0] c = [5.0, 0.0] points = [a, b, c] triangles = [ [1,0,2] ] #bac d1 = [1.0, 1.0] d2 = [1.0, 3.0] d3 = [3.0, 1.0] z1 = [2, 4] z2 = [4, 8] z3 = [4, 8] data_coords = [d1, d2, d3] z = [z1, z2, z3] f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0) answer = [[0, 0], [5., 10.], [5., 10.]] assert allclose(f, answer) #Tests of smoothing matrix def test_smoothing_matrix_one_triangle(self): from Numeric import dot a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] points = [a, b, c] vertices = [ [1,0,2] ] #bac interp = Interpolation(points, vertices) assert allclose(interp.get_D(), [[1, -0.5, -0.5], [-0.5, 0.5, 0], [-0.5, 0, 0.5]]) #Define f(x,y) = x f = array([0,0,2]) #Value at global vertex 2 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 1 dx dy = area = 2 assert dot(dot(f, interp.get_D()), f) == 2 #Define f(x,y) = y f = array([0,2,0]) #Value at global vertex 1 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 1 dx dy = area = 2 assert dot(dot(f, interp.get_D()), f) == 2 #Define f(x,y) = x+y f = array([0,2,2]) #Values at global vertex 1 and 2 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 2 dx dy = 2*area = 4 assert dot(dot(f, interp.get_D()), f) == 4 def test_smoothing_matrix_more_triangles(self): from Numeric import dot a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe, daf, dae vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] interp = Interpolation(points, vertices) #assert allclose(interp.get_D(), [[1, -0.5, -0.5], # [-0.5, 0.5, 0], # [-0.5, 0, 0.5]]) #Define f(x,y) = x f = array([0,0,2,0,2,4]) #f evaluated at points a-f #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 1 dx dy = total area = 8 assert dot(dot(f, interp.get_D()), f) == 8 #Define f(x,y) = y f = array([0,2,0,4,2,0]) #f evaluated at points a-f #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 1 dx dy = area = 8 assert dot(dot(f, interp.get_D()), f) == 8 #Define f(x,y) = x+y f = array([0,2,2,4,4,4]) #f evaluated at points a-f #Check that int (df/dx)**2 + (df/dy)**2 dx dy = # int 2 dx dy = 2*area = 16 assert dot(dot(f, interp.get_D()), f) == 16 def test_fit_and_interpolation(self): from mesh import Mesh a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0, 0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0, 0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe, daf, dae triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] #Get (enough) datapoints data_points = [[ 0.66666667, 0.66666667], [ 1.33333333, 1.33333333], [ 2.66666667, 0.66666667], [ 0.66666667, 2.66666667], [ 0.0, 1.0], [ 0.0, 3.0], [ 1.0, 0.0], [ 1.0, 1.0], [ 1.0, 2.0], [ 1.0, 3.0], [ 2.0, 1.0], [ 3.0, 0.0], [ 3.0, 1.0]] interp = Interpolation(points, triangles, data_points, alpha=0.0) z = linear_function(data_points) answer = linear_function(points) f = interp.fit(z) #print "f",f #print "answer",answer assert allclose(f, answer) #Map back z1 = interp.interpolate(f) #print "z1\n", z1 #print "z\n",z assert allclose(z, z1) def test_smoothing_and_interpolation(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0, 0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0, 0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe, daf, dae triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] #Get (too few!) datapoints data_points = [[ 0.66666667, 0.66666667], [ 1.33333333, 1.33333333], [ 2.66666667, 0.66666667], [ 0.66666667, 2.66666667]] z = linear_function(data_points) answer = linear_function(points) #Make interpolator with too few data points and no smoothing interp = Interpolation(points, triangles, data_points, alpha=0.0) #Must raise an exception try: f = interp.fit(z) except: pass #Now try with smoothing parameter interp = Interpolation(points, triangles, data_points, alpha=1.0e-13) f = interp.fit(z) #f will be different from answerr due to smoothing assert allclose(f, answer,atol=5) #Map back z1 = interp.interpolate(f) assert allclose(z, z1) def test_fit_to_mesh_file(self): from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \ export_trianglulation_file import tempfile import os # create a .tsh file, no user outline mesh_dic = {} mesh_dic['generatedpointlist'] = [[0.0, 0.0], [0.0, 5.0], [5.0, 0.0]] mesh_dic['generatedtrianglelist'] = [[0, 2, 1]] mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]] mesh_dic['generatedtriangleattributelist'] = [['']] mesh_dic['generatedpointattributelist'] = [[], [], []] mesh_dic['generatedpointattributetitlelist'] = [] mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]] mesh_dic['generatedsegmentmarkerlist'] = ['external', 'external', 'external'] mesh_file = tempfile.mktemp(".tsh") export_trianglulation_file(mesh_file,mesh_dic) # create an .xya file point_file = tempfile.mktemp(".xya") fd = open(point_file,'w') fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n") fd.close() mesh_output_file = "new_trianlge.tsh" fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha = 0.0) # load in the .tsh file we just wrote mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file) assert allclose(mesh_dic['generatedpointattributelist'], [[0.0, 0.0], [5.0, 10.0], [5.0,10.0]]) self.failUnless(mesh_dic['generatedpointattributetitlelist'] == ['elevation','stage'], 'test_fit_to_mesh_file failed') #clean up os.remove(mesh_file) os.remove(point_file) def test_fit_to_mesh_fileII(self): from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \ export_trianglulation_file import tempfile import os # create a .tsh file, no user outline mesh_dic = {} mesh_dic['generatedpointlist'] = [[0.0, 0.0], [0.0, 5.0], [5.0, 0.0]] mesh_dic['generatedtrianglelist'] = [[0, 2, 1]] mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]] mesh_dic['generatedtriangleattributelist'] = [['']] mesh_dic['generatedpointattributelist'] = [[1,2], [1,2], [1,2]] mesh_dic['generatedpointattributetitlelist'] = ['density', 'temp'] mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]] mesh_dic['generatedsegmentmarkerlist'] = ['external', 'external', 'external'] mesh_file = tempfile.mktemp(".tsh") export_trianglulation_file(mesh_file,mesh_dic) # create an .xya file point_file = tempfile.mktemp(".xya") fd = open(point_file,'w') fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n") fd.close() mesh_output_file = "new_triangle.tsh" fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha = 0.0) # load in the .tsh file we just wrote mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file) assert allclose(mesh_dic['generatedpointattributelist'], [[1.0, 2.0,0.0, 0.0], [1.0, 2.0,5.0, 10.0], [1.0, 2.0,5.0,10.0]]) self.failUnless(mesh_dic['generatedpointattributetitlelist'] == ['density', 'temp','elevation','stage'], 'test_fit_to_mesh_file failed') #clean up os.remove(mesh_file) os.remove(mesh_output_file) os.remove(point_file) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') runner = unittest.TextTestRunner(verbosity=1) runner.run(suite)