source: storm_surge/pyvolution/test_least_squares.py @ 1

Last change on this file since 1 was 1, checked in by duncan, 20 years ago

initial import
File size: 15.8 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt
5
6from least_squares import *
7from Numeric import allclose, array
8
9def distance(x, y):
10    return sqrt( sum( (array(x)-array(y))**2 ))         
11
12def linear_function(point):
13    point = array(point)
14    return point[:,0]+point[:,1]
15   
16       
17class TestCase(unittest.TestCase):
18
19    def setUp(self):
20        pass
21       
22    def tearDown(self):
23        pass
24
25
26    def test_datapoint_at_centroid(self):
27        a = [0.0, 0.0]
28        b = [0.0, 2.0]
29        c = [2.0,0.0]
30        points = [a, b, c]
31        vertices = [ [1,0,2] ]   #bac
32
33        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
34       
35        interp = Interpolation(points, vertices, data)
36        assert allclose(interp.A, [[1./3, 1./3, 1./3]]) 
37       
38
39
40    def test_datapoints_at_vertices(self):
41        """Test that data points coinciding with vertices yield a diagonal matrix
42        """
43       
44        a = [0.0, 0.0]
45        b = [0.0, 2.0]
46        c = [2.0,0.0]
47        points = [a, b, c]
48        vertices = [ [1,0,2] ]   #bac
49
50        data = points #Use data at vertices       
51       
52        interp = Interpolation(points, vertices, data)
53        assert allclose(interp.A, [[1., 0., 0.],
54                                   [0., 1., 0.],
55                                   [0., 0., 1.]]) 
56       
57
58
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       
64        a = [0.0, 0.0]
65        b = [0.0, 2.0]
66        c = [2.0,0.0]
67        points = [a, b, c]
68        vertices = [ [1,0,2] ]   #bac
69
70        data = [ [0., 1.], [1., 0.], [1., 1.] ]
71       
72        interp = Interpolation(points, vertices, data)
73       
74        assert allclose(interp.A, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
75                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
76                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
77
78
79    def test_datapoints_on_edges(self):
80        """Try datapoints on edges -
81        each point should affect two matrix entries in proportion
82        """
83       
84        a = [0.0, 0.0]
85        b = [0.0, 2.0]
86        c = [2.0,0.0]
87        points = [a, b, c]
88        vertices = [ [1,0,2] ]   #bac
89
90        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
91       
92        interp = Interpolation(points, vertices, data)
93       
94        assert allclose(interp.A, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
95                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
96                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
97
98    def test_arbitrary_datapoints(self):
99        """Try arbitrary datapoints
100        """
101
102        from Numeric import sum
103       
104        a = [0.0, 0.0]
105        b = [0.0, 2.0]
106        c = [2.0,0.0]
107        points = [a, b, c]
108        vertices = [ [1,0,2] ]   #bac
109
110        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
111       
112        interp = Interpolation(points, vertices, data)
113
114        assert allclose(sum(interp.A, axis=1), 1.0)
115       
116
117           
118    def test_more_triangles(self):
119        a = [-1.0, 0.0]
120        b = [3.0, 4.0]
121        c = [4.0,1.0]
122        d = [-3.0, 2.0] #3
123        e = [-1.0,-2.0]
124        f = [1.0, -2.0] #5
125
126        points = [a, b, c, d,e,f]
127        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
128
129        #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       
134        answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],   #Affects point d     
135                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
136                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
137                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
138                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
139                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f
140        assert allclose(interp.A, answer)
141
142
143
144
145    def test_smooth_attributes_to_mesh(self):
146        a = [0.0, 0.0]
147        b = [0.0, 5.0]
148        c = [5.0, 0.0]
149        points = [a, b, c]
150        triangles = [ [1,0,2] ]   #bac
151         
152        d1 = [1.0, 1.0]
153        d2 = [1.0, 3.0]
154        d3 = [3.0,1.0]
155        z1 = 2
156        z2 = 4
157        z3 = 4
158        data_coords = [d1, d2, d3] 
159       
160        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
161        z = [z1, z2, z3]
162        f = interp.fit(z)
163        answer = [0, 5., 5.]
164       
165        assert allclose(f, answer)
166       
167       
168    def test_smooth_att_to_meshII(self):
169        a = [0.0, 0.0]
170        b = [0.0, 5.0]
171        c = [5.0, 0.0]
172        points = [a, b, c]
173        triangles = [ [1,0,2] ]   #bac
174         
175        d1 = [1.0, 1.0]
176        d2 = [1.0, 2.0]
177        d3 = [3.0,1.0]
178        data_coords = [d1, d2, d3]       
179        z = linear_function(data_coords)
180        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
181        f = interp.fit(z)
182        answer = linear_function(points)
183       
184        assert allclose(f, answer)
185   
186    def test_smooth_attributes_to_meshIII(self):
187        a = [-1.0, 0.0]
188        b = [3.0, 4.0]
189        c = [4.0,1.0]
190        d = [-3.0, 2.0] #3
191        e = [-1.0,-2.0]
192        f = [1.0, -2.0] #5
193
194        vertices = [a, b, c, d,e,f]
195        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
196
197        point_coords = [[-2.0, 2.0],
198                        [-1.0, 1.0],
199                        [0.0,2.0],
200                        [1.0, 1.0],
201                        [2.0, 1.0],
202                        [0.0,0.0],
203                        [1.0, 0.0],
204                        [0.0, -1.0],
205                        [-0.2,-0.5],
206                        [-0.9, -1.5],
207                        [0.5, -1.9],
208                        [3.0,1.0]]
209       
210        z = linear_function(point_coords)
211        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
212        f = interp.fit(z)
213        answer = linear_function(vertices)
214        assert allclose(f, answer)
215       
216
217    def test_smooth_attributes_to_meshIV(self):
218        """ Testing 2 attributes smoothed to the mesh
219        """
220        a = [0.0, 0.0]
221        b = [0.0, 5.0]
222        c = [5.0, 0.0]
223        points = [a, b, c]
224        triangles = [ [1,0,2] ]   #bac
225         
226        d1 = [1.0, 1.0]
227        d2 = [1.0, 3.0]
228        d3 = [3.0,1.0]
229        z1 = [2,4]
230        z2 = [4,8]
231        z3 = [4,8]
232        data_coords = [ d1, d2, d3] 
233       
234        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
235        z = [z1, z2, z3]
236        f =  interp.fit(z)
237        answer = [[0,0], [5., 10.], [5., 10.]]
238        assert allclose(f, answer)
239       
240    def test_interpolate_attributes_to_points(self):
241        v0 = [0.0, 0.0]
242        v1 = [0.0, 5.0]
243        v2 = [5.0, 0.0]
244       
245        vertices = [v0, v1, v2]
246        triangles = [ [1,0,2] ]   #bac
247         
248        d0 = [1.0, 1.0]
249        d1 = [1.0, 2.0]
250        d2 = [3.0,1.0]
251        point_coords = [ d0, d1, d2]
252       
253        interp = Interpolation(vertices, triangles, point_coords)
254        f = linear_function(vertices)
255        z = interp.interpolate(f)
256        answer = linear_function(point_coords)
257       
258
259        assert allclose(z, answer)
260
261       
262    def test_interpolate_attributes_to_pointsII(self):
263        a = [-1.0, 0.0]
264        b = [3.0, 4.0]
265        c = [4.0,1.0]
266        d = [-3.0, 2.0] #3
267        e = [-1.0,-2.0]
268        f = [1.0, -2.0] #5
269
270        vertices = [a, b, c, d,e,f]
271        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
272
273         
274        point_coords = [[-2.0, 2.0],
275                        [-1.0, 1.0],
276                        [0.0,2.0],
277                        [1.0, 1.0],
278                        [2.0, 1.0],
279                        [0.0,0.0],
280                        [1.0, 0.0],
281                        [0.0, -1.0],
282                        [-0.2,-0.5],
283                        [-0.9, -1.5],
284                        [0.5, -1.9],
285                        [3.0,1.0]]
286       
287        interp = Interpolation(vertices, triangles, point_coords)
288        f = linear_function(vertices)
289        z = interp.interpolate(f)
290        answer = linear_function(point_coords)
291        #print "z",z
292        #print "answer",answer
293        assert allclose(z, answer)
294
295       
296       
297    def test_smooth_attributes_to_mesh_function(self):
298        """ Testing 2 attributes smoothed to the mesh
299        """
300        """Test multiple attributes
301        """
302       
303        a = [0.0, 0.0]
304        b = [0.0, 5.0]
305        c = [5.0, 0.0]
306        points = [a, b, c]
307        triangles = [ [1,0,2] ]   #bac
308         
309        d1 = [1.0, 1.0]
310        d2 = [1.0, 3.0]
311        d3 = [3.0,1.0]
312        z1 = [2,4]
313        z2 = [4,8]
314        z3 = [4,8]
315        data_coords = [ d1, d2, d3]
316        z = [z1, z2, z3]
317       
318        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
319        answer = [[0, 0], [5., 10.], [5., 10.]]
320       
321        assert allclose(f, answer)
322
323
324
325    #Tests of smoothing matrix   
326    def test_smoothing_matrix_one_triangle(self):
327        from Numeric import dot
328        a = [0.0, 0.0]
329        b = [0.0, 2.0]
330        c = [2.0,0.0]
331        points = [a, b, c]
332       
333        vertices = [ [1,0,2] ]   #bac
334
335        interp = Interpolation(points, vertices)
336
337        assert allclose(interp.D, [[1, -0.5, -0.5],
338                                   [-0.5, 0.5, 0],
339                                   [-0.5, 0, 0.5]])
340
341        #Define f(x,y) = x
342        f = array([0,0,2]) #Value at global vertex 2
343
344        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
345        #           int 1 dx dy = area = 2
346        assert dot(dot(f, interp.D), f) == 2
347       
348        #Define f(x,y) = y
349        f = array([0,2,0])  #Value at global vertex 1
350
351        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
352        #           int 1 dx dy = area = 2
353        assert dot(dot(f, interp.D), f) == 2
354
355        #Define f(x,y) = x+y
356        f = array([0,2,2])  #Values at global vertex 1 and 2
357
358        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
359        #           int 2 dx dy = 2*area = 4
360        assert dot(dot(f, interp.D), f) == 4       
361       
362
363
364    def test_smoothing_matrix_more_triangles(self):
365        from Numeric import dot
366       
367        a = [0.0, 0.0]
368        b = [0.0, 2.0]
369        c = [2.0,0.0]
370        d = [0.0, 4.0]
371        e = [2.0, 2.0]
372        f = [4.0,0.0]
373
374        points = [a, b, c, d, e, f]
375        #bac, bce, ecf, dbe, daf, dae
376        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
377
378        interp = Interpolation(points, vertices)
379
380
381        #assert allclose(interp.D, [[1, -0.5, -0.5],
382        #                           [-0.5, 0.5, 0],
383        #                           [-0.5, 0, 0.5]])
384
385        #Define f(x,y) = x
386        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
387
388        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
389        #           int 1 dx dy = total area = 8
390        assert dot(dot(f, interp.D), f) == 8
391       
392        #Define f(x,y) = y
393        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
394
395        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
396        #           int 1 dx dy = area = 8
397        assert dot(dot(f, interp.D), f) == 8
398
399        #Define f(x,y) = x+y
400        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
401
402        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
403        #           int 2 dx dy = 2*area = 16
404        assert dot(dot(f, interp.D), f) == 16       
405
406
407    def test_fit_and_interpolation(self):
408        from mesh import Mesh
409       
410        a = [0.0, 0.0]
411        b = [0.0, 2.0]
412        c = [2.0, 0.0]
413        d = [0.0, 4.0]
414        e = [2.0, 2.0]
415        f = [4.0, 0.0]
416
417        points = [a, b, c, d, e, f]
418        #bac, bce, ecf, dbe, daf, dae
419        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
420
421        #Get (enough) datapoints
422        data_points = [[ 0.66666667, 0.66666667],
423                       [ 1.33333333, 1.33333333],
424                       [ 2.66666667, 0.66666667],
425                       [ 0.66666667, 2.66666667],
426                       [ 0.0, 1.0],
427                       [ 0.0, 3.0],
428                       [ 1.0, 0.0],
429                       [ 1.0, 1.0],
430                       [ 1.0, 2.0],
431                       [ 1.0, 3.0],                                         
432                       [ 2.0, 1.0],
433                       [ 3.0, 0.0],
434                       [ 3.0, 1.0]]                                         
435       
436        interp = Interpolation(points, triangles, data_points, alpha=0.0)
437
438        z = linear_function(data_points)
439        answer = linear_function(points)
440
441        f = interp.fit(z)
442
443        assert allclose(f, answer)
444
445        #Map back
446        z1 = interp.interpolate(f)
447        assert allclose(z, z1)
448
449
450    def test_smoothing_and_interpolation(self):
451       
452        a = [0.0, 0.0]
453        b = [0.0, 2.0]
454        c = [2.0, 0.0]
455        d = [0.0, 4.0]
456        e = [2.0, 2.0]
457        f = [4.0, 0.0]
458
459        points = [a, b, c, d, e, f]
460        #bac, bce, ecf, dbe, daf, dae
461        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
462
463        #Get (too few!) datapoints
464        data_points = [[ 0.66666667, 0.66666667],
465                       [ 1.33333333, 1.33333333],
466                       [ 2.66666667, 0.66666667],
467                       [ 0.66666667, 2.66666667]]
468
469        z = linear_function(data_points)
470        answer = linear_function(points)       
471
472        #Make interpolator with too few data points and no smoothing
473        interp = Interpolation(points, triangles, data_points, alpha=0.0)
474        #Must raise an exception       
475        try:
476            f = interp.fit(z)
477        except:
478            pass
479
480        #Now try with smoothing parameter
481        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
482       
483        f = interp.fit(z)
484
485        #f will be different from answerr due to smoothing
486        #assert allclose(f, answer)
487
488        #Map back
489        z1 = interp.interpolate(f)
490        assert allclose(z, z1)
491
492    def test_fit_to_mesh_file(self):
493        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
494             export_trianglulation_file   
495        import tempfile
496        import os
497       
498        # create a .tsh file, no user outline
499        mesh_dic = {}
500        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
501                                          [0.0, 5.0],
502                                          [5.0, 0.0]]
503        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
504        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
505        mesh_dic['generatedtriangleattributelist'] = [['']]
506        mesh_dic['generatedpointattributelist'] = [[], [], []]
507        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
508        mesh_dic['generatedsegmentmarkerlist'] = ['external',
509                                                  'external',
510                                                  'external']       
511        mesh_file = tempfile.mktemp(".tsh")
512        export_trianglulation_file(mesh_file,mesh_dic)
513       
514        # create an .xya file
515        point_file = tempfile.mktemp(".xya")
516        fd = open(point_file,'w')
517        fd.write("# demo \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
518        fd.close()
519       
520        mesh_output_file = "new_trianlge.tsh"
521        fit_to_mesh_file(mesh_file,
522                         point_file,
523                         mesh_output_file,
524                         alpha = 0.0) 
525        # load in the .tsh file we jusdt wrote
526        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
527
528        assert allclose(mesh_dic['generatedpointattributelist'],
529                        [[0.0, 0.0], [5.0, 10.0], [5.0,10.0]])
530       
531        #clean up
532        os.remove(mesh_file)
533        os.remove(point_file)
534       
535#-------------------------------------------------------------
536if __name__ == "__main__":
537    suite = unittest.makeSuite(TestCase,'test')
538    runner = unittest.TextTestRunner()
539    runner.run(suite)
540
541   
542   
543
544
545
Note: See TracBrowser for help on using the repository browser.