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

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

fixed bug in least_squares algorithm.
File size: 25.8 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4
5import unittest
6from math import sqrt
7
8
9from least_squares import *
10from Numeric import allclose, array, transpose
11
12def distance(x, y):
13    return sqrt( sum( (array(x)-array(y))**2 ))         
14
15def linear_function(point):
16    point = array(point)
17    return point[:,0]+point[:,1]
18   
19       
20class TestCase(unittest.TestCase):
21
22    def setUp(self):
23        pass
24       
25    def tearDown(self):
26        pass
27
28
29    def test_datapoint_at_centroid(self):
30        a = [0.0, 0.0]
31        b = [0.0, 2.0]
32        c = [2.0,0.0]
33        points = [a, b, c]
34        vertices = [ [1,0,2] ]   #bac
35
36        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
37       
38        interp = Interpolation(points, vertices, data)
39        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 
40       
41
42    def test_quad_tree(self):
43        p0 = [-10.0, -10.0]
44        p1 = [20.0, -10.0]
45        p2 = [-10.0, 20.0]
46        p3 = [10.0, 50.0]
47        p4 = [30.0, 30.0]
48        p5 = [50.0, 10.0]
49        p6 = [40.0, 60.0]
50        p7 = [60.0, 40.0]
51        p8 = [-66.0, 20.0]
52        p9 = [10.0, -66.0]
53       
54        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
55        triangles = [ [0, 1, 2],
56                      [3, 2, 4],
57                      [4, 2, 1],
58                      [4, 1, 5],
59                      [3, 4, 6],
60                      [6, 4, 7],
61                      [7, 4, 5],
62                      [8, 0, 2],
63                      [0, 9, 1]]   
64
65        data = [ [4,4] ] 
66       
67        interp = Interpolation(points, triangles, data, alpha = 0.0)
68        #print "PDSG - interp.get_A()", interp.get_A()
69        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
70                      0., 0. , 0., 0., 0., 0.]]
71        assert allclose(interp.get_A(), answer)
72        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
73        #print "PDSG - interp.get_A()", interp.get_A()
74        answer =  [ [ 0.0,  0.0,  0.0,  0.,
75                      0., 0. , 0., 0., 0., 0.]]
76        assert allclose(interp.get_A(), answer)
77
78       
79        #point outside of quad tree root cell
80        interp.set_point_coordinates([[-70, -70]])
81        #print "PDSG - interp.get_A()", interp.get_A()
82        answer =  [ [ 0.0,  0.0,  0.0,  0.,
83                      0., 0. , 0., 0., 0., 0.]]
84        assert allclose(interp.get_A(), answer)
85     
86    def test_quad_treeII(self):
87        p0 = [-66.0, 14.0]
88        p1 = [14.0, -66.0]
89        p2 = [14.0, 14.0]
90        p3 = [60.0, 20.0]
91        p4 = [10.0, 60.0]
92        p5 = [60.0, 60.0]
93       
94        points = [p0, p1, p2, p3, p4, p5]
95        triangles = [ [0, 1, 2],
96                      [3, 2, 1],
97                      [0, 2, 4],
98                      [0, 2, 4],
99                      [4, 2, 5],
100                      [5, 2, 3]]   
101
102        data = [ [-26.0,-26.0] ] 
103       
104        interp = Interpolation(points, triangles, data, alpha = 0.0)
105        #print "PDSG - interp.get_A()", interp.get_A()
106        answer =  [ [ 0.5,  0.5,  0.0,  0.,
107                      0., 0.]]
108        assert allclose(interp.get_A(), answer)
109        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
110        #print "PDSG -30,-30 - interp.get_A()", interp.get_A()
111        answer =  [ [ 0.0,  0.0,  0.0,  0.,
112                      0., 0.]]
113        assert allclose(interp.get_A(), answer)
114
115       
116        #point outside of quad tree root cell
117        interp.set_point_coordinates([[-70, -70]])
118        #print "PDSG -70,-70 interp.get_A()", interp.get_A()
119        answer =  [ [ 0.0,  0.0,  0.0,  0.,
120                      0., 0. ]]
121        assert allclose(interp.get_A(), answer)
122           
123
124    def test_datapoints_at_vertices(self):
125        """Test that data points coinciding with vertices yield a diagonal matrix
126        """
127       
128        a = [0.0, 0.0]
129        b = [0.0, 2.0]
130        c = [2.0,0.0]
131        points = [a, b, c]
132        vertices = [ [1,0,2] ]   #bac
133
134        data = points #Use data at vertices       
135       
136        interp = Interpolation(points, vertices, data)
137        assert allclose(interp.get_A(), [[1., 0., 0.],
138                                   [0., 1., 0.],
139                                   [0., 0., 1.]]) 
140       
141
142
143    def test_datapoints_on_edge_midpoints(self):
144        """Try datapoints midway on edges -
145        each point should affect two matrix entries equally
146        """
147       
148        a = [0.0, 0.0]
149        b = [0.0, 2.0]
150        c = [2.0,0.0]
151        points = [a, b, c]
152        vertices = [ [1,0,2] ]   #bac
153
154        data = [ [0., 1.], [1., 0.], [1., 1.] ]
155       
156        interp = Interpolation(points, vertices, data)
157       
158        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
159                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
160                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
161
162
163    def test_datapoints_on_edges(self):
164        """Try datapoints on edges -
165        each point should affect two matrix entries in proportion
166        """
167       
168        a = [0.0, 0.0]
169        b = [0.0, 2.0]
170        c = [2.0,0.0]
171        points = [a, b, c]
172        vertices = [ [1,0,2] ]   #bac
173
174        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
175       
176        interp = Interpolation(points, vertices, data)
177       
178        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
179                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
180                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
181
182    def test_arbitrary_datapoints(self):
183        """Try arbitrary datapoints
184        """
185
186        from Numeric import sum
187       
188        a = [0.0, 0.0]
189        b = [0.0, 2.0]
190        c = [2.0,0.0]
191        points = [a, b, c]
192        vertices = [ [1,0,2] ]   #bac
193
194        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
195       
196        interp = Interpolation(points, vertices, data)
197        #print "interp.get_A()", interp.get_A()
198        assert allclose(sum(interp.get_A(), axis=1), 1.0)
199       
200
201    # this causes a memory error in scipy.sparce       
202    def test_more_triangles(self):
203       
204        a = [-1.0, 0.0]
205        b = [3.0, 4.0]
206        c = [4.0,1.0]
207        d = [-3.0, 2.0] #3
208        e = [-1.0,-2.0]
209        f = [1.0, -2.0] #5
210
211        points = [a, b, c, d,e,f]
212        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
213
214        #Data points
215        data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
216        interp = Interpolation(points, triangles, data_points) 
217
218        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d     
219                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
220                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
221                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
222                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
223                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
224
225
226        A = interp.get_A()
227        for i in range(A.shape[0]):
228            for j in range(A.shape[1]):
229                if not allclose(A[i,j], answer[i][j]):
230                    print i,j,':',A[i,j], answer[i][j] 
231
232
233        assert allclose(interp.get_A(), answer)
234
235
236
237
238    def test_smooth_attributes_to_mesh(self):
239        a = [0.0, 0.0]
240        b = [0.0, 5.0]
241        c = [5.0, 0.0]
242        points = [a, b, c]
243        triangles = [ [1,0,2] ]   #bac
244         
245        d1 = [1.0, 1.0]
246        d2 = [1.0, 3.0]
247        d3 = [3.0,1.0]
248        z1 = 2
249        z2 = 4
250        z3 = 4
251        data_coords = [d1, d2, d3] 
252       
253        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
254        z = [z1, z2, z3]
255        f = interp.fit(z)
256        answer = [0, 5., 5.]
257
258        #print "f\n",f
259        #print "answer\n",answer
260       
261        assert allclose(f, answer, atol=1e-7)
262       
263       
264    def test_smooth_att_to_meshII(self):
265       
266        a = [0.0, 0.0]
267        b = [0.0, 5.0]
268        c = [5.0, 0.0]
269        points = [a, b, c]
270        triangles = [ [1,0,2] ]   #bac
271         
272        d1 = [1.0, 1.0]
273        d2 = [1.0, 2.0]
274        d3 = [3.0,1.0]
275        data_coords = [d1, d2, d3]       
276        z = linear_function(data_coords)
277        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
278        f = interp.fit(z)
279        answer = linear_function(points)
280       
281        assert allclose(f, answer)
282   
283    def test_smooth_attributes_to_meshIII(self):
284       
285        a = [-1.0, 0.0]
286        b = [3.0, 4.0]
287        c = [4.0,1.0]
288        d = [-3.0, 2.0] #3
289        e = [-1.0,-2.0]
290        f = [1.0, -2.0] #5
291
292        vertices = [a, b, c, d,e,f]
293        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
294
295        point_coords = [[-2.0, 2.0],
296                        [-1.0, 1.0],
297                        [0.0,2.0],
298                        [1.0, 1.0],
299                        [2.0, 1.0],
300                        [0.0,0.0],
301                        [1.0, 0.0],
302                        [0.0, -1.0],
303                        [-0.2,-0.5],
304                        [-0.9, -1.5],
305                        [0.5, -1.9],
306                        [3.0,1.0]]
307
308        z = linear_function(point_coords)
309        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
310
311        #print 'z',z
312        f = interp.fit(z)
313        answer = linear_function(vertices)
314        #print "f\n",f
315        #print "answer\n",answer
316        assert allclose(f, answer)
317       
318
319    def test_smooth_attributes_to_meshIV(self):
320        """ Testing 2 attributes smoothed to the mesh
321        """
322       
323        a = [0.0, 0.0]
324        b = [0.0, 5.0]
325        c = [5.0, 0.0]
326        points = [a, b, c]
327        triangles = [ [1,0,2] ]   #bac
328         
329        d1 = [1.0, 1.0]
330        d2 = [1.0, 3.0]
331        d3 = [3.0, 1.0]
332        z1 = [2, 4]
333        z2 = [4, 8]
334        z3 = [4, 8]
335        data_coords = [d1, d2, d3] 
336       
337        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
338        z = [z1, z2, z3]
339        f =  interp.fit_points(z)
340        answer = [[0,0], [5., 10.], [5., 10.]]
341        assert allclose(f, answer)
342       
343    def test_interpolate_attributes_to_points(self):
344        v0 = [0.0, 0.0]
345        v1 = [0.0, 5.0]
346        v2 = [5.0, 0.0]
347       
348        vertices = [v0, v1, v2]
349        triangles = [ [1,0,2] ]   #bac
350         
351        d0 = [1.0, 1.0]
352        d1 = [1.0, 2.0]
353        d2 = [3.0, 1.0]
354        point_coords = [ d0, d1, d2]
355       
356        interp = Interpolation(vertices, triangles, point_coords)
357        f = linear_function(vertices)
358        z = interp.interpolate(f)
359        answer = linear_function(point_coords)
360       
361
362        assert allclose(z, answer)
363
364       
365    def test_interpolate_attributes_to_pointsII(self):
366        a = [-1.0, 0.0]
367        b = [3.0, 4.0]
368        c = [4.0, 1.0]
369        d = [-3.0, 2.0] #3
370        e = [-1.0, -2.0]
371        f = [1.0, -2.0] #5
372
373        vertices = [a, b, c, d,e,f]
374        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
375
376         
377        point_coords = [[-2.0, 2.0],
378                        [-1.0, 1.0],
379                        [0.0, 2.0],
380                        [1.0, 1.0],
381                        [2.0, 1.0],
382                        [0.0, 0.0],
383                        [1.0, 0.0],
384                        [0.0, -1.0],
385                        [-0.2, -0.5],
386                        [-0.9, -1.5],
387                        [0.5, -1.9],
388                        [3.0, 1.0]]
389       
390        interp = Interpolation(vertices, triangles, point_coords)
391        f = linear_function(vertices)
392        z = interp.interpolate(f)
393        answer = linear_function(point_coords)
394        #print "z",z
395        #print "answer",answer
396        assert allclose(z, answer)
397   
398    def test_interpolate_attributes_to_pointsIII(self):
399        v0 = [0.0, 0.0]
400        v1 = [0.0, 5.0]
401        v2 = [5.0, 0.0]
402       
403        vertices = [v0, v1, v2]
404        triangles = [ [1,0,2] ]   #bac
405         
406        d0 = [1.0, 1.0]
407        d1 = [1.0, 2.0]
408        d2 = [3.0,1.0]
409        point_coords = [ d0, d1, d2]
410       
411        interp = Interpolation(vertices, triangles, point_coords)
412        f = [ [0., 0.],  [5., 10.],  [5., 10.]] #linear_function(vertices)
413        #print "f",f
414       
415        z = interp.interpolate(f)
416        answer = [ [2., 4.],  [3., 6.],  [4., 8.]]
417
418        #print "***********"
419        #print "z",z
420        #print "answer",answer
421        #print "***********"
422
423        assert allclose(z, answer)
424   
425    def test_interpolate_attributes_to_pointsIV(self):
426        a = [-1.0, 0.0]
427        b = [3.0, 4.0]
428        c = [4.0, 1.0]
429        d = [-3.0, 2.0] #3
430        e = [-1.0, -2.0]
431        f = [1.0, -2.0] #5
432
433        vertices = [a, b, c, d,e,f]
434        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
435
436         
437        point_coords = [[-2.0, 2.0],
438                        [-1.0, 1.0],
439                        [0.0, 2.0],
440                        [1.0, 1.0],
441                        [2.0, 1.0],
442                        [0.0, 0.0],
443                        [1.0, 0.0],
444                        [0.0, -1.0],
445                        [-0.2, -0.5],
446                        [-0.9, -1.5],
447                        [0.5, -1.9],
448                        [3.0, 1.0]]
449       
450        interp = Interpolation(vertices, triangles, point_coords)
451        f = array([linear_function(vertices),2*linear_function(vertices) ])
452        f = transpose(f)
453        #print "f",f
454        z = interp.interpolate(f)
455        answer = [linear_function(point_coords),
456             
457     2*linear_function(point_coords) ]
458        answer = transpose(answer)
459        #print "z",z
460        #print "answer",answer
461        assert allclose(z, answer)
462       
463    def test_smooth_attributes_to_mesh_function(self):
464        """ Testing 2 attributes smoothed to the mesh
465        """
466        """Test multiple attributes
467        """
468       
469        a = [0.0, 0.0]
470        b = [0.0, 5.0]
471        c = [5.0, 0.0]
472        points = [a, b, c]
473        triangles = [ [1,0,2] ]   #bac
474         
475        d1 = [1.0, 1.0]
476        d2 = [1.0, 3.0]
477        d3 = [3.0, 1.0]
478        z1 = [2, 4]
479        z2 = [4, 8]
480        z3 = [4, 8]
481        data_coords = [d1, d2, d3]
482        z = [z1, z2, z3]
483       
484        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
485        answer = [[0, 0], [5., 10.], [5., 10.]]
486       
487        assert allclose(f, answer)
488
489
490
491    #Tests of smoothing matrix   
492    def test_smoothing_matrix_one_triangle(self):
493        from Numeric import dot
494        a = [0.0, 0.0]
495        b = [0.0, 2.0]
496        c = [2.0,0.0]
497        points = [a, b, c]
498       
499        vertices = [ [1,0,2] ]   #bac
500
501        interp = Interpolation(points, vertices)
502
503        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
504                                   [-0.5, 0.5, 0],
505                                   [-0.5, 0, 0.5]])
506
507        #Define f(x,y) = x
508        f = array([0,0,2]) #Value at global vertex 2
509
510        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
511        #           int 1 dx dy = area = 2
512        assert dot(dot(f, interp.get_D()), f) == 2
513       
514        #Define f(x,y) = y
515        f = array([0,2,0])  #Value at global vertex 1
516
517        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
518        #           int 1 dx dy = area = 2
519        assert dot(dot(f, interp.get_D()), f) == 2
520
521        #Define f(x,y) = x+y
522        f = array([0,2,2])  #Values at global vertex 1 and 2
523
524        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
525        #           int 2 dx dy = 2*area = 4
526        assert dot(dot(f, interp.get_D()), f) == 4       
527       
528
529
530    def test_smoothing_matrix_more_triangles(self):
531        from Numeric import dot
532       
533        a = [0.0, 0.0]
534        b = [0.0, 2.0]
535        c = [2.0,0.0]
536        d = [0.0, 4.0]
537        e = [2.0, 2.0]
538        f = [4.0,0.0]
539
540        points = [a, b, c, d, e, f]
541        #bac, bce, ecf, dbe, daf, dae
542        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
543
544        interp = Interpolation(points, vertices)
545
546
547        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
548        #                           [-0.5, 0.5, 0],
549        #                           [-0.5, 0, 0.5]])
550
551        #Define f(x,y) = x
552        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
553
554        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
555        #           int 1 dx dy = total area = 8
556        assert dot(dot(f, interp.get_D()), f) == 8
557       
558        #Define f(x,y) = y
559        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
560
561        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
562        #           int 1 dx dy = area = 8
563        assert dot(dot(f, interp.get_D()), f) == 8
564
565        #Define f(x,y) = x+y
566        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
567
568        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
569        #           int 2 dx dy = 2*area = 16
570        assert dot(dot(f, interp.get_D()), f) == 16       
571
572
573    def test_fit_and_interpolation(self):
574        from mesh import Mesh
575       
576        a = [0.0, 0.0]
577        b = [0.0, 2.0]
578        c = [2.0, 0.0]
579        d = [0.0, 4.0]
580        e = [2.0, 2.0]
581        f = [4.0, 0.0]
582
583        points = [a, b, c, d, e, f]
584        #bac, bce, ecf, dbe, daf, dae
585        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
586
587        #Get (enough) datapoints
588        data_points = [[ 0.66666667, 0.66666667],
589                       [ 1.33333333, 1.33333333],
590                       [ 2.66666667, 0.66666667],
591                       [ 0.66666667, 2.66666667],
592                       [ 0.0, 1.0],
593                       [ 0.0, 3.0],
594                       [ 1.0, 0.0],
595                       [ 1.0, 1.0],
596                       [ 1.0, 2.0],
597                       [ 1.0, 3.0],                                         
598                       [ 2.0, 1.0],
599                       [ 3.0, 0.0],
600                       [ 3.0, 1.0]]                                         
601       
602        interp = Interpolation(points, triangles, data_points, alpha=0.0)
603
604        z = linear_function(data_points)
605        answer = linear_function(points)
606
607        f = interp.fit(z)
608
609        #print "f",f
610        #print "answer",answer
611        assert allclose(f, answer)
612
613        #Map back
614        z1 = interp.interpolate(f)
615        #print "z1\n", z1
616        #print "z\n",z
617        assert allclose(z, z1)
618
619
620    def test_smoothing_and_interpolation(self):
621       
622        a = [0.0, 0.0]
623        b = [0.0, 2.0]
624        c = [2.0, 0.0]
625        d = [0.0, 4.0]
626        e = [2.0, 2.0]
627        f = [4.0, 0.0]
628
629        points = [a, b, c, d, e, f]
630        #bac, bce, ecf, dbe, daf, dae
631        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
632
633        #Get (too few!) datapoints
634        data_points = [[ 0.66666667, 0.66666667],
635                       [ 1.33333333, 1.33333333],
636                       [ 2.66666667, 0.66666667],
637                       [ 0.66666667, 2.66666667]]
638
639        z = linear_function(data_points)
640        answer = linear_function(points)       
641
642        #Make interpolator with too few data points and no smoothing
643        interp = Interpolation(points, triangles, data_points, alpha=0.0)
644        #Must raise an exception       
645        try:
646            f = interp.fit(z)
647        except:
648            pass
649
650        #Now try with smoothing parameter
651        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
652       
653        f = interp.fit(z)
654        #f will be different from answerr due to smoothing
655        assert allclose(f, answer,atol=5)
656
657        #Map back
658        z1 = interp.interpolate(f)
659        assert allclose(z, z1)
660
661
662
663    def test_fit_and_interpolation_with_new_points(self):
664        """Fit a surface to one set of points. Then interpolate that surface
665        using another set of points.
666        """
667        from mesh import Mesh
668
669
670        #Setup mesh used to represent fitted function
671        a = [0.0, 0.0]
672        b = [0.0, 2.0]
673        c = [2.0, 0.0]
674        d = [0.0, 4.0]
675        e = [2.0, 2.0]
676        f = [4.0, 0.0]
677
678        points = [a, b, c, d, e, f]
679        #bac, bce, ecf, dbe, daf, dae
680        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
681
682        #Datapoints to fit from
683        data_points1 = [[ 0.66666667, 0.66666667],
684                        [ 1.33333333, 1.33333333],
685                        [ 2.66666667, 0.66666667],
686                        [ 0.66666667, 2.66666667],
687                        [ 0.0, 1.0],
688                        [ 0.0, 3.0],
689                        [ 1.0, 0.0],
690                        [ 1.0, 1.0],
691                        [ 1.0, 2.0],
692                        [ 1.0, 3.0],                                         
693                        [ 2.0, 1.0],
694                        [ 3.0, 0.0],
695                        [ 3.0, 1.0]]
696
697        #Fit surface to mesh
698        interp = Interpolation(points, triangles, data_points1, alpha=0.0)
699        z = linear_function(data_points1) #Example z-values
700        f = interp.fit(z)                 #Fitted values at vertices
701
702
703
704        #New datapoints where interpolated values are sought
705        data_points2 = [[ 0.0, 0.0],
706                        [ 0.5, 0.5],
707                        [ 0.7, 0.7],
708                        [ 1.0, 0.5],
709                        [ 2.0, 0.4],
710                        [ 2.8, 1.2]]                                         
711       
712
713        #Build new A matrix based on new points
714        interp.build_interpolation_matrix_A(data_points2)
715
716        #Interpolate using fitted surface
717        z1 = interp.interpolate(f)
718
719        #Desired result
720        answer = linear_function(data_points2)
721        assert allclose(z1, answer)
722
723
724
725
726    def test_fit_to_mesh_file(self):
727        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
728             export_trianglulation_file   
729        import tempfile
730        import os
731       
732        # create a .tsh file, no user outline
733        mesh_dic = {}
734        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
735                                          [0.0, 5.0],
736                                          [5.0, 0.0]]
737        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
738        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
739        mesh_dic['generatedtriangleattributelist'] = [['']]
740        mesh_dic['generatedpointattributelist'] = [[], [], []]
741        mesh_dic['generatedpointattributetitlelist'] = []
742        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
743        mesh_dic['generatedsegmentmarkerlist'] = ['external',
744                                                  'external',
745                                                  'external']       
746        mesh_file = tempfile.mktemp(".tsh")
747        export_trianglulation_file(mesh_file,mesh_dic)
748       
749        # create an .xya file
750        point_file = tempfile.mktemp(".xya")
751        fd = open(point_file,'w')
752        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")
753        fd.close()
754       
755        mesh_output_file = "new_trianlge.tsh"
756        fit_to_mesh_file(mesh_file,
757                         point_file,
758                         mesh_output_file,
759                         alpha = 0.0) 
760        # load in the .tsh file we just wrote
761        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
762
763        assert allclose(mesh_dic['generatedpointattributelist'],
764                        [[0.0, 0.0],
765                         [5.0, 10.0],
766                         [5.0,10.0]])
767       
768        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
769                        ['elevation','stage'],
770                        'test_fit_to_mesh_file failed')
771       
772        #clean up
773        os.remove(mesh_file)
774        os.remove(point_file)
775       
776    def test_fit_to_mesh_fileII(self):
777        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
778             export_trianglulation_file   
779        import tempfile
780        import os
781       
782        # create a .tsh file, no user outline
783        mesh_dic = {}
784        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
785                                          [0.0, 5.0],
786                                          [5.0, 0.0]]
787        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
788        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
789        mesh_dic['generatedtriangleattributelist'] = [['']]
790        mesh_dic['generatedpointattributelist'] = [[1,2], [1,2], [1,2]]
791        mesh_dic['generatedpointattributetitlelist'] = ['density', 'temp']
792        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
793        mesh_dic['generatedsegmentmarkerlist'] = ['external',
794                                                  'external',
795                                                  'external']       
796        mesh_file = tempfile.mktemp(".tsh")
797        export_trianglulation_file(mesh_file,mesh_dic)
798       
799        # create an .xya file
800        point_file = tempfile.mktemp(".xya")
801        fd = open(point_file,'w')
802        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")
803        fd.close()
804       
805        mesh_output_file = "new_triangle.tsh"
806        fit_to_mesh_file(mesh_file,
807                         point_file,
808                         mesh_output_file,
809                         alpha = 0.0) 
810        # load in the .tsh file we just wrote
811        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
812
813        assert allclose(mesh_dic['generatedpointattributelist'],
814                        [[1.0, 2.0,0.0, 0.0],
815                         [1.0, 2.0,5.0, 10.0],
816                         [1.0, 2.0,5.0,10.0]])
817       
818        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
819                        ['density', 'temp','elevation','stage'],
820                        'test_fit_to_mesh_file failed')
821       
822        #clean up
823        os.remove(mesh_file)
824        os.remove(mesh_output_file)       
825        os.remove(point_file)
826       
827#-------------------------------------------------------------
828if __name__ == "__main__":
829    #suite = unittest.makeSuite(TestCase,'test')
830
831    suite = unittest.makeSuite(TestCase,'test')
832    #suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
833    runner = unittest.TextTestRunner(verbosity=1)
834    runner.run(suite)
835
836   
837   
838
839
840
Note: See TracBrowser for help on using the repository browser.