source: inundation/ga/storm_surge/pyvolution/test_least_squares.py @ 353

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

added more least squares tests, investigating sww files
File size: 17.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, transpose
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    def test_interpolate_attributes_to_pointsIII(self):
296        v0 = [0.0, 0.0]
297        v1 = [0.0, 5.0]
298        v2 = [5.0, 0.0]
299       
300        vertices = [v0, v1, v2]
301        triangles = [ [1,0,2] ]   #bac
302         
303        d0 = [1.0, 1.0]
304        d1 = [1.0, 2.0]
305        d2 = [3.0,1.0]
306        point_coords = [ d0, d1, d2]
307       
308        interp = Interpolation(vertices, triangles, point_coords)
309        f = [ [0., 0.],  [5., 10.],  [5., 10.]] #linear_function(vertices)
310        #print "f",f
311       
312        z = interp.interpolate(f)
313        answer = [ [2., 4.],  [3., 6.],  [4., 8.]]
314
315        #print "***********"
316        #print "z",z
317        #print "answer",answer
318        #print "***********"
319
320        assert allclose(z, answer)
321   
322    def test_interpolate_attributes_to_pointsIV(self):
323        a = [-1.0, 0.0]
324        b = [3.0, 4.0]
325        c = [4.0,1.0]
326        d = [-3.0, 2.0] #3
327        e = [-1.0,-2.0]
328        f = [1.0, -2.0] #5
329
330        vertices = [a, b, c, d,e,f]
331        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
332
333         
334        point_coords = [[-2.0, 2.0],
335                        [-1.0, 1.0],
336                        [0.0,2.0],
337                        [1.0, 1.0],
338                        [2.0, 1.0],
339                        [0.0,0.0],
340                        [1.0, 0.0],
341                        [0.0, -1.0],
342                        [-0.2,-0.5],
343                        [-0.9, -1.5],
344                        [0.5, -1.9],
345                        [3.0,1.0]]
346       
347        interp = Interpolation(vertices, triangles, point_coords)
348        f = array([linear_function(vertices),2*linear_function(vertices) ])
349        f = transpose(f)
350        #print "f",f
351        z = interp.interpolate(f)
352        answer = [linear_function(point_coords),
353             
354     2*linear_function(point_coords) ]
355        answer = transpose(answer)
356        #print "z",z
357        #print "answer",answer
358        assert allclose(z, answer)
359       
360    def test_smooth_attributes_to_mesh_function(self):
361        """ Testing 2 attributes smoothed to the mesh
362        """
363        """Test multiple attributes
364        """
365       
366        a = [0.0, 0.0]
367        b = [0.0, 5.0]
368        c = [5.0, 0.0]
369        points = [a, b, c]
370        triangles = [ [1,0,2] ]   #bac
371         
372        d1 = [1.0, 1.0]
373        d2 = [1.0, 3.0]
374        d3 = [3.0,1.0]
375        z1 = [2,4]
376        z2 = [4,8]
377        z3 = [4,8]
378        data_coords = [ d1, d2, d3]
379        z = [z1, z2, z3]
380       
381        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
382        answer = [[0, 0], [5., 10.], [5., 10.]]
383       
384        assert allclose(f, answer)
385
386
387
388    #Tests of smoothing matrix   
389    def test_smoothing_matrix_one_triangle(self):
390        from Numeric import dot
391        a = [0.0, 0.0]
392        b = [0.0, 2.0]
393        c = [2.0,0.0]
394        points = [a, b, c]
395       
396        vertices = [ [1,0,2] ]   #bac
397
398        interp = Interpolation(points, vertices)
399
400        assert allclose(interp.D, [[1, -0.5, -0.5],
401                                   [-0.5, 0.5, 0],
402                                   [-0.5, 0, 0.5]])
403
404        #Define f(x,y) = x
405        f = array([0,0,2]) #Value at global vertex 2
406
407        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
408        #           int 1 dx dy = area = 2
409        assert dot(dot(f, interp.D), f) == 2
410       
411        #Define f(x,y) = y
412        f = array([0,2,0])  #Value at global vertex 1
413
414        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
415        #           int 1 dx dy = area = 2
416        assert dot(dot(f, interp.D), f) == 2
417
418        #Define f(x,y) = x+y
419        f = array([0,2,2])  #Values at global vertex 1 and 2
420
421        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
422        #           int 2 dx dy = 2*area = 4
423        assert dot(dot(f, interp.D), f) == 4       
424       
425
426
427    def test_smoothing_matrix_more_triangles(self):
428        from Numeric import dot
429       
430        a = [0.0, 0.0]
431        b = [0.0, 2.0]
432        c = [2.0,0.0]
433        d = [0.0, 4.0]
434        e = [2.0, 2.0]
435        f = [4.0,0.0]
436
437        points = [a, b, c, d, e, f]
438        #bac, bce, ecf, dbe, daf, dae
439        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
440
441        interp = Interpolation(points, vertices)
442
443
444        #assert allclose(interp.D, [[1, -0.5, -0.5],
445        #                           [-0.5, 0.5, 0],
446        #                           [-0.5, 0, 0.5]])
447
448        #Define f(x,y) = x
449        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
450
451        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
452        #           int 1 dx dy = total area = 8
453        assert dot(dot(f, interp.D), f) == 8
454       
455        #Define f(x,y) = y
456        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
457
458        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
459        #           int 1 dx dy = area = 8
460        assert dot(dot(f, interp.D), f) == 8
461
462        #Define f(x,y) = x+y
463        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
464
465        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
466        #           int 2 dx dy = 2*area = 16
467        assert dot(dot(f, interp.D), f) == 16       
468
469
470    def test_fit_and_interpolation(self):
471        from mesh import Mesh
472       
473        a = [0.0, 0.0]
474        b = [0.0, 2.0]
475        c = [2.0, 0.0]
476        d = [0.0, 4.0]
477        e = [2.0, 2.0]
478        f = [4.0, 0.0]
479
480        points = [a, b, c, d, e, f]
481        #bac, bce, ecf, dbe, daf, dae
482        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
483
484        #Get (enough) datapoints
485        data_points = [[ 0.66666667, 0.66666667],
486                       [ 1.33333333, 1.33333333],
487                       [ 2.66666667, 0.66666667],
488                       [ 0.66666667, 2.66666667],
489                       [ 0.0, 1.0],
490                       [ 0.0, 3.0],
491                       [ 1.0, 0.0],
492                       [ 1.0, 1.0],
493                       [ 1.0, 2.0],
494                       [ 1.0, 3.0],                                         
495                       [ 2.0, 1.0],
496                       [ 3.0, 0.0],
497                       [ 3.0, 1.0]]                                         
498       
499        interp = Interpolation(points, triangles, data_points, alpha=0.0)
500
501        z = linear_function(data_points)
502        answer = linear_function(points)
503
504        f = interp.fit(z)
505
506        assert allclose(f, answer)
507
508        #Map back
509        z1 = interp.interpolate(f)
510        assert allclose(z, z1)
511
512
513    def test_smoothing_and_interpolation(self):
514       
515        a = [0.0, 0.0]
516        b = [0.0, 2.0]
517        c = [2.0, 0.0]
518        d = [0.0, 4.0]
519        e = [2.0, 2.0]
520        f = [4.0, 0.0]
521
522        points = [a, b, c, d, e, f]
523        #bac, bce, ecf, dbe, daf, dae
524        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
525
526        #Get (too few!) datapoints
527        data_points = [[ 0.66666667, 0.66666667],
528                       [ 1.33333333, 1.33333333],
529                       [ 2.66666667, 0.66666667],
530                       [ 0.66666667, 2.66666667]]
531
532        z = linear_function(data_points)
533        answer = linear_function(points)       
534
535        #Make interpolator with too few data points and no smoothing
536        interp = Interpolation(points, triangles, data_points, alpha=0.0)
537        #Must raise an exception       
538        try:
539            f = interp.fit(z)
540        except:
541            pass
542
543        #Now try with smoothing parameter
544        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
545       
546        f = interp.fit(z)
547
548        #f will be different from answerr due to smoothing
549        #assert allclose(f, answer)
550
551        #Map back
552        z1 = interp.interpolate(f)
553        assert allclose(z, z1)
554
555    def test_fit_to_mesh_file(self):
556        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
557             export_trianglulation_file   
558        import tempfile
559        import os
560       
561        # create a .tsh file, no user outline
562        mesh_dic = {}
563        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
564                                          [0.0, 5.0],
565                                          [5.0, 0.0]]
566        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
567        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
568        mesh_dic['generatedtriangleattributelist'] = [['']]
569        mesh_dic['generatedpointattributelist'] = [[], [], []]
570        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
571        mesh_dic['generatedsegmentmarkerlist'] = ['external',
572                                                  'external',
573                                                  'external']       
574        mesh_file = tempfile.mktemp(".tsh")
575        export_trianglulation_file(mesh_file,mesh_dic)
576       
577        # create an .xya file
578        point_file = tempfile.mktemp(".xya")
579        fd = open(point_file,'w')
580        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")
581        fd.close()
582       
583        mesh_output_file = "new_trianlge.tsh"
584        fit_to_mesh_file(mesh_file,
585                         point_file,
586                         mesh_output_file,
587                         alpha = 0.0) 
588        # load in the .tsh file we just wrote
589        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
590
591        assert allclose(mesh_dic['generatedpointattributelist'],
592                        [[0.0, 0.0], [5.0, 10.0], [5.0,10.0]])
593       
594        #clean up
595        os.remove(mesh_file)
596        os.remove(point_file)
597       
598#-------------------------------------------------------------
599if __name__ == "__main__":
600    suite = unittest.makeSuite(TestCase,'test')
601    runner = unittest.TextTestRunner()
602    runner.run(suite)
603
604   
605   
606
607
608
Note: See TracBrowser for help on using the repository browser.