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

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

reducing the # of errors when testing.
Adding fit_points function that doesn't rely on broadcasting, since cg_solve doesn't broadcast.

File size: 20.5 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.get_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        print "aa"
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.get_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        print "bb"
64       
65        a = [0.0, 0.0]
66        b = [0.0, 2.0]
67        c = [2.0,0.0]
68        points = [a, b, c]
69        vertices = [ [1,0,2] ]   #bac
70
71        data = [ [0., 1.], [1., 0.], [1., 1.] ]
72       
73        interp = Interpolation(points, vertices, data)
74       
75        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
76                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
77                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
78
79
80    def test_datapoints_on_edges(self):
81        """Try datapoints on edges -
82        each point should affect two matrix entries in proportion
83        """
84        print "cc"
85       
86        a = [0.0, 0.0]
87        b = [0.0, 2.0]
88        c = [2.0,0.0]
89        points = [a, b, c]
90        vertices = [ [1,0,2] ]   #bac
91
92        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
93       
94        interp = Interpolation(points, vertices, data)
95       
96        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
97                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
98                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
99
100    def test_arbitrary_datapoints(self):
101        """Try arbitrary datapoints
102        """
103
104        from Numeric import sum
105       
106        a = [0.0, 0.0]
107        b = [0.0, 2.0]
108        c = [2.0,0.0]
109        points = [a, b, c]
110        vertices = [ [1,0,2] ]   #bac
111
112        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
113       
114        interp = Interpolation(points, vertices, data)
115
116        assert allclose(sum(interp.get_A(), axis=1), 1.0)
117       
118
119    # this causes a memory error in scipy.sparce       
120    def bad_test_more_triangles(self):
121        print "dd"
122        a = [-1.0, 0.0]
123        b = [3.0, 4.0]
124        c = [4.0,1.0]
125        d = [-3.0, 2.0] #3
126        e = [-1.0,-2.0]
127        f = [1.0, -2.0] #5
128
129        points = [a, b, c, d,e,f]
130        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
131
132        #Data points
133        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
134        interp = Interpolation(points, triangles, data_points) 
135       
136        answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],   #Affects point d     
137                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
138                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
139                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
140                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
141                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f
142
143        assert allclose(interp.get_A(), answer)
144
145
146
147
148    def test_smooth_attributes_to_mesh(self):
149        print "ee"
150        a = [0.0, 0.0]
151        b = [0.0, 5.0]
152        c = [5.0, 0.0]
153        points = [a, b, c]
154        triangles = [ [1,0,2] ]   #bac
155         
156        d1 = [1.0, 1.0]
157        d2 = [1.0, 3.0]
158        d3 = [3.0,1.0]
159        z1 = 2
160        z2 = 4
161        z3 = 4
162        data_coords = [d1, d2, d3] 
163       
164        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
165        z = [z1, z2, z3]
166        f = interp.fit(z)
167        answer = [0, 5., 5.]
168        assert allclose(f, answer)
169       
170       
171    def test_smooth_att_to_meshII(self):
172        print "ff"
173        a = [0.0, 0.0]
174        b = [0.0, 5.0]
175        c = [5.0, 0.0]
176        points = [a, b, c]
177        triangles = [ [1,0,2] ]   #bac
178         
179        d1 = [1.0, 1.0]
180        d2 = [1.0, 2.0]
181        d3 = [3.0,1.0]
182        data_coords = [d1, d2, d3]       
183        z = linear_function(data_coords)
184        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
185        f = interp.fit(z)
186        answer = linear_function(points)
187       
188        assert allclose(f, answer)
189   
190    def test_smooth_attributes_to_meshIII(self):
191        print "gg"
192        a = [-1.0, 0.0]
193        b = [3.0, 4.0]
194        c = [4.0,1.0]
195        d = [-3.0, 2.0] #3
196        e = [-1.0,-2.0]
197        f = [1.0, -2.0] #5
198
199        vertices = [a, b, c, d,e,f]
200        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
201
202        point_coords = [[-2.0, 2.0],
203                        [-1.0, 1.0],
204                        [0.0,2.0],
205                        [1.0, 1.0],
206                        [2.0, 1.0],
207                        [0.0,0.0],
208                        [1.0, 0.0],
209                        [0.0, -1.0],
210                        [-0.2,-0.5],
211                        [-0.9, -1.5],
212                        [0.5, -1.9],
213                        [3.0,1.0]]
214       
215        z = linear_function(point_coords)
216        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
217        f = interp.fit(z)
218        answer = linear_function(vertices)
219        assert allclose(f, answer)
220       
221
222    def test_smooth_attributes_to_meshIV(self):
223        print "hh"
224        """ Testing 2 attributes smoothed to the mesh
225        """
226        print "test_smooth_attributes_to_meshIV"
227        a = [0.0, 0.0]
228        b = [0.0, 5.0]
229        c = [5.0, 0.0]
230        points = [a, b, c]
231        triangles = [ [1,0,2] ]   #bac
232         
233        d1 = [1.0, 1.0]
234        d2 = [1.0, 3.0]
235        d3 = [3.0,1.0]
236        z1 = [2,4]
237        z2 = [4,8]
238        z3 = [4,8]
239        data_coords = [ d1, d2, d3] 
240       
241        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
242        z = [z1, z2, z3]
243        f =  interp.fit_points(z)
244        answer = [[0,0], [5., 10.], [5., 10.]]
245        assert allclose(f, answer)
246       
247    def test_interpolate_attributes_to_points(self):
248        v0 = [0.0, 0.0]
249        v1 = [0.0, 5.0]
250        v2 = [5.0, 0.0]
251       
252        vertices = [v0, v1, v2]
253        triangles = [ [1,0,2] ]   #bac
254         
255        d0 = [1.0, 1.0]
256        d1 = [1.0, 2.0]
257        d2 = [3.0,1.0]
258        point_coords = [ d0, d1, d2]
259       
260        interp = Interpolation(vertices, triangles, point_coords)
261        f = linear_function(vertices)
262        z = interp.interpolate(f)
263        answer = linear_function(point_coords)
264       
265
266        assert allclose(z, answer)
267
268       
269    def test_interpolate_attributes_to_pointsII(self):
270        a = [-1.0, 0.0]
271        b = [3.0, 4.0]
272        c = [4.0,1.0]
273        d = [-3.0, 2.0] #3
274        e = [-1.0,-2.0]
275        f = [1.0, -2.0] #5
276
277        vertices = [a, b, c, d,e,f]
278        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
279
280         
281        point_coords = [[-2.0, 2.0],
282                        [-1.0, 1.0],
283                        [0.0,2.0],
284                        [1.0, 1.0],
285                        [2.0, 1.0],
286                        [0.0,0.0],
287                        [1.0, 0.0],
288                        [0.0, -1.0],
289                        [-0.2,-0.5],
290                        [-0.9, -1.5],
291                        [0.5, -1.9],
292                        [3.0,1.0]]
293       
294        interp = Interpolation(vertices, triangles, point_coords)
295        f = linear_function(vertices)
296        z = interp.interpolate(f)
297        answer = linear_function(point_coords)
298        #print "z",z
299        #print "answer",answer
300        assert allclose(z, answer)
301   
302    def test_interpolate_attributes_to_pointsIII(self):
303        v0 = [0.0, 0.0]
304        v1 = [0.0, 5.0]
305        v2 = [5.0, 0.0]
306       
307        vertices = [v0, v1, v2]
308        triangles = [ [1,0,2] ]   #bac
309         
310        d0 = [1.0, 1.0]
311        d1 = [1.0, 2.0]
312        d2 = [3.0,1.0]
313        point_coords = [ d0, d1, d2]
314       
315        interp = Interpolation(vertices, triangles, point_coords)
316        f = [ [0., 0.],  [5., 10.],  [5., 10.]] #linear_function(vertices)
317        #print "f",f
318       
319        z = interp.interpolate(f)
320        answer = [ [2., 4.],  [3., 6.],  [4., 8.]]
321
322        #print "***********"
323        #print "z",z
324        #print "answer",answer
325        #print "***********"
326
327        assert allclose(z, answer)
328   
329    def test_interpolate_attributes_to_pointsIV(self):
330        a = [-1.0, 0.0]
331        b = [3.0, 4.0]
332        c = [4.0,1.0]
333        d = [-3.0, 2.0] #3
334        e = [-1.0,-2.0]
335        f = [1.0, -2.0] #5
336
337        vertices = [a, b, c, d,e,f]
338        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
339
340         
341        point_coords = [[-2.0, 2.0],
342                        [-1.0, 1.0],
343                        [0.0,2.0],
344                        [1.0, 1.0],
345                        [2.0, 1.0],
346                        [0.0,0.0],
347                        [1.0, 0.0],
348                        [0.0, -1.0],
349                        [-0.2,-0.5],
350                        [-0.9, -1.5],
351                        [0.5, -1.9],
352                        [3.0,1.0]]
353       
354        interp = Interpolation(vertices, triangles, point_coords)
355        f = array([linear_function(vertices),2*linear_function(vertices) ])
356        f = transpose(f)
357        #print "f",f
358        z = interp.interpolate(f)
359        answer = [linear_function(point_coords),
360             
361     2*linear_function(point_coords) ]
362        answer = transpose(answer)
363        #print "z",z
364        #print "answer",answer
365        assert allclose(z, answer)
366       
367    def test_smooth_attributes_to_mesh_function(self):
368        """ Testing 2 attributes smoothed to the mesh
369        """
370        """Test multiple attributes
371        """
372       
373        a = [0.0, 0.0]
374        b = [0.0, 5.0]
375        c = [5.0, 0.0]
376        points = [a, b, c]
377        triangles = [ [1,0,2] ]   #bac
378         
379        d1 = [1.0, 1.0]
380        d2 = [1.0, 3.0]
381        d3 = [3.0,1.0]
382        z1 = [2,4]
383        z2 = [4,8]
384        z3 = [4,8]
385        data_coords = [ d1, d2, d3]
386        z = [z1, z2, z3]
387       
388        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
389        answer = [[0, 0], [5., 10.], [5., 10.]]
390       
391        assert allclose(f, answer)
392
393
394
395    #Tests of smoothing matrix   
396    def test_smoothing_matrix_one_triangle(self):
397        from Numeric import dot
398        a = [0.0, 0.0]
399        b = [0.0, 2.0]
400        c = [2.0,0.0]
401        points = [a, b, c]
402       
403        vertices = [ [1,0,2] ]   #bac
404
405        interp = Interpolation(points, vertices)
406
407        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
408                                   [-0.5, 0.5, 0],
409                                   [-0.5, 0, 0.5]])
410
411        #Define f(x,y) = x
412        f = array([0,0,2]) #Value at global vertex 2
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.get_D()), f) == 2
417       
418        #Define f(x,y) = y
419        f = array([0,2,0])  #Value at global vertex 1
420
421        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
422        #           int 1 dx dy = area = 2
423        assert dot(dot(f, interp.get_D()), f) == 2
424
425        #Define f(x,y) = x+y
426        f = array([0,2,2])  #Values at global vertex 1 and 2
427
428        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
429        #           int 2 dx dy = 2*area = 4
430        assert dot(dot(f, interp.get_D()), f) == 4       
431       
432
433
434    def test_smoothing_matrix_more_triangles(self):
435        from Numeric import dot
436       
437        a = [0.0, 0.0]
438        b = [0.0, 2.0]
439        c = [2.0,0.0]
440        d = [0.0, 4.0]
441        e = [2.0, 2.0]
442        f = [4.0,0.0]
443
444        points = [a, b, c, d, e, f]
445        #bac, bce, ecf, dbe, daf, dae
446        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
447
448        interp = Interpolation(points, vertices)
449
450
451        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
452        #                           [-0.5, 0.5, 0],
453        #                           [-0.5, 0, 0.5]])
454
455        #Define f(x,y) = x
456        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
457
458        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
459        #           int 1 dx dy = total area = 8
460        assert dot(dot(f, interp.get_D()), f) == 8
461       
462        #Define f(x,y) = y
463        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
464
465        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
466        #           int 1 dx dy = area = 8
467        assert dot(dot(f, interp.get_D()), f) == 8
468
469        #Define f(x,y) = x+y
470        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
471
472        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
473        #           int 2 dx dy = 2*area = 16
474        assert dot(dot(f, interp.get_D()), f) == 16       
475
476
477    def test_fit_and_interpolation(self):
478        from mesh import Mesh
479       
480        a = [0.0, 0.0]
481        b = [0.0, 2.0]
482        c = [2.0, 0.0]
483        d = [0.0, 4.0]
484        e = [2.0, 2.0]
485        f = [4.0, 0.0]
486
487        points = [a, b, c, d, e, f]
488        #bac, bce, ecf, dbe, daf, dae
489        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
490
491        #Get (enough) datapoints
492        data_points = [[ 0.66666667, 0.66666667],
493                       [ 1.33333333, 1.33333333],
494                       [ 2.66666667, 0.66666667],
495                       [ 0.66666667, 2.66666667],
496                       [ 0.0, 1.0],
497                       [ 0.0, 3.0],
498                       [ 1.0, 0.0],
499                       [ 1.0, 1.0],
500                       [ 1.0, 2.0],
501                       [ 1.0, 3.0],                                         
502                       [ 2.0, 1.0],
503                       [ 3.0, 0.0],
504                       [ 3.0, 1.0]]                                         
505       
506        interp = Interpolation(points, triangles, data_points, alpha=0.0)
507
508        z = linear_function(data_points)
509        answer = linear_function(points)
510
511        f = interp.fit(z)
512
513        assert allclose(f, answer)
514
515        #Map back
516        z1 = interp.interpolate(f)
517        assert allclose(z, z1)
518
519
520    def test_smoothing_and_interpolation(self):
521       
522        a = [0.0, 0.0]
523        b = [0.0, 2.0]
524        c = [2.0, 0.0]
525        d = [0.0, 4.0]
526        e = [2.0, 2.0]
527        f = [4.0, 0.0]
528
529        points = [a, b, c, d, e, f]
530        #bac, bce, ecf, dbe, daf, dae
531        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
532
533        #Get (too few!) datapoints
534        data_points = [[ 0.66666667, 0.66666667],
535                       [ 1.33333333, 1.33333333],
536                       [ 2.66666667, 0.66666667],
537                       [ 0.66666667, 2.66666667]]
538
539        z = linear_function(data_points)
540        answer = linear_function(points)       
541
542        #Make interpolator with too few data points and no smoothing
543        interp = Interpolation(points, triangles, data_points, alpha=0.0)
544        #Must raise an exception       
545        try:
546            f = interp.fit(z)
547        except:
548            pass
549
550        #Now try with smoothing parameter
551        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
552       
553        f = interp.fit(z)
554
555        #f will be different from answerr due to smoothing
556        #assert allclose(f, answer)
557
558        #Map back
559        z1 = interp.interpolate(f)
560        assert allclose(z, z1)
561
562    def test_fit_to_mesh_file(self):
563        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
564             export_trianglulation_file   
565        import tempfile
566        import os
567       
568        # create a .tsh file, no user outline
569        mesh_dic = {}
570        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
571                                          [0.0, 5.0],
572                                          [5.0, 0.0]]
573        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
574        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
575        mesh_dic['generatedtriangleattributelist'] = [['']]
576        mesh_dic['generatedpointattributelist'] = [[], [], []]
577        mesh_dic['generatedpointattributetitlelist'] = []
578        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
579        mesh_dic['generatedsegmentmarkerlist'] = ['external',
580                                                  'external',
581                                                  'external']       
582        mesh_file = tempfile.mktemp(".tsh")
583        export_trianglulation_file(mesh_file,mesh_dic)
584       
585        # create an .xya file
586        point_file = tempfile.mktemp(".xya")
587        fd = open(point_file,'w')
588        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")
589        fd.close()
590       
591        mesh_output_file = "new_trianlge.tsh"
592        fit_to_mesh_file(mesh_file,
593                         point_file,
594                         mesh_output_file,
595                         alpha = 0.0) 
596        # load in the .tsh file we just wrote
597        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
598
599        assert allclose(mesh_dic['generatedpointattributelist'],
600                        [[0.0, 0.0],
601                         [5.0, 10.0],
602                         [5.0,10.0]])
603       
604        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
605                        ['elevation','stage'],
606                        'test_fit_to_mesh_file failed')
607       
608        #clean up
609        os.remove(mesh_file)
610        os.remove(point_file)
611       
612    def test_fit_to_mesh_fileII(self):
613        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
614             export_trianglulation_file   
615        import tempfile
616        import os
617       
618        # create a .tsh file, no user outline
619        mesh_dic = {}
620        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
621                                          [0.0, 5.0],
622                                          [5.0, 0.0]]
623        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
624        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
625        mesh_dic['generatedtriangleattributelist'] = [['']]
626        mesh_dic['generatedpointattributelist'] = [[1,2], [1,2], [1,2]]
627        mesh_dic['generatedpointattributetitlelist'] = ['density', 'temp']
628        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
629        mesh_dic['generatedsegmentmarkerlist'] = ['external',
630                                                  'external',
631                                                  'external']       
632        mesh_file = tempfile.mktemp(".tsh")
633        export_trianglulation_file(mesh_file,mesh_dic)
634       
635        # create an .xya file
636        point_file = tempfile.mktemp(".xya")
637        fd = open(point_file,'w')
638        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")
639        fd.close()
640       
641        mesh_output_file = "new_trianlge.tsh"
642        fit_to_mesh_file(mesh_file,
643                         point_file,
644                         mesh_output_file,
645                         alpha = 0.0) 
646        # load in the .tsh file we just wrote
647        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
648
649        assert allclose(mesh_dic['generatedpointattributelist'],
650                        [[1.0, 2.0,0.0, 0.0],
651                         [1.0, 2.0,5.0, 10.0],
652                         [1.0, 2.0,5.0,10.0]])
653       
654        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
655                        ['density', 'temp','elevation','stage'],
656                        'test_fit_to_mesh_file failed')
657       
658        #clean up
659        os.remove(mesh_file)
660        os.remove(point_file)
661       
662#-------------------------------------------------------------
663if __name__ == "__main__":
664    suite = unittest.makeSuite(TestCase,'test')
665    runner = unittest.TextTestRunner()
666    runner.run(suite)
667
668   
669   
670
671
672
Note: See TracBrowser for help on using the repository browser.