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

Last change on this file since 999 was 997, checked in by ole, 20 years ago

Implemented precropping facility in least squares
Also optimised inside_polygon

File size: 36.0 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    def test_datapoint_at_centroid(self):
29        a = [0.0, 0.0]
30        b = [0.0, 2.0]
31        c = [2.0,0.0]
32        points = [a, b, c]
33        vertices = [ [1,0,2] ]   #bac
34
35        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
36       
37        interp = Interpolation(points, vertices, data)
38        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 
39       
40
41    def test_quad_tree(self):
42        p0 = [-10.0, -10.0]
43        p1 = [20.0, -10.0]
44        p2 = [-10.0, 20.0]
45        p3 = [10.0, 50.0]
46        p4 = [30.0, 30.0]
47        p5 = [50.0, 10.0]
48        p6 = [40.0, 60.0]
49        p7 = [60.0, 40.0]
50        p8 = [-66.0, 20.0]
51        p9 = [10.0, -66.0]
52       
53        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
54        triangles = [ [0, 1, 2],
55                      [3, 2, 4],
56                      [4, 2, 1],
57                      [4, 1, 5],
58                      [3, 4, 6],
59                      [6, 4, 7],
60                      [7, 4, 5],
61                      [8, 0, 2],
62                      [0, 9, 1]]   
63
64        data = [ [4,4] ] 
65        interp = Interpolation(points, triangles, data, alpha = 0.0,
66                               max_points_per_cell = 4)
67        #print "PDSG - interp.get_A()", interp.get_A()
68        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
69                      0., 0. , 0., 0., 0., 0.]]
70        assert allclose(interp.get_A(), answer)
71        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
72        #print "PDSG - interp.get_A()", interp.get_A()
73        answer =  [ [ 0.0,  0.0,  0.0,  0.,
74                      0., 0. , 0., 0., 0., 0.]]
75        assert allclose(interp.get_A(), answer)
76
77       
78        #point outside of quad tree root cell
79        interp.set_point_coordinates([[-70, -70]])
80        #print "PDSG - interp.get_A()", interp.get_A()
81        answer =  [ [ 0.0,  0.0,  0.0,  0.,
82                      0., 0. , 0., 0., 0., 0.]]
83        assert allclose(interp.get_A(), answer)
84     
85    def test_expand_search(self):
86        p0 = [-10.0, -10.0]
87        p1 = [20.0, -10.0]
88        p2 = [-10.0, 20.0]
89        p3 = [10.0, 50.0]
90        p4 = [30.0, 30.0]
91        p5 = [50.0, 10.0]
92        p6 = [40.0, 60.0]
93        p7 = [60.0, 40.0]
94        p8 = [-66.0, 20.0]
95        p9 = [10.0, -66.0]
96       
97        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
98        triangles = [ [0, 1, 2],
99                      [3, 2, 4],
100                      [4, 2, 1],
101                      [4, 1, 5],
102                      [3, 4, 6],
103                      [6, 4, 7],
104                      [7, 4, 5],
105                      [8, 0, 2],
106                      [0, 9, 1]]   
107
108        data = [ [4,4],
109                 [-30,10],
110                 [-20,0],
111                 [-20,10],
112                 [0,30],
113                 [10,-40],
114                 [10,-30],
115                 [10,-20],
116                 [10,10],
117                 [10,20],
118                 [10,30],
119                 [10,40],
120                 [20,10],
121                 [25,45],
122                 [30,0],
123                 [30,10],
124                 [30,30],
125                 [30,40],
126                 [30,50],
127                 [40,10],
128                 [40,30],
129                 [40,40],
130                 [40,50],
131                 [50,20],
132                 [50,30],
133                 [50,40],
134                 [50,50],
135                 [30,0],
136                 [-20,-20]] 
137        point_attributes = [ -400000,
138                     10,
139                     10,
140                     10,
141                    10,
142                     10,
143                     10,
144                     10,
145                     10,
146                     10,
147                     10,
148                     10,
149                     10,
150                     10,
151                    10,
152                     10,
153                     10,
154                     10,
155                     10,
156                     10,
157                     10,
158                     10,
159                     10,
160                     10,
161                     10,
162                     10,
163                     10,
164                   10,
165                   99] 
166       
167        interp = Interpolation(points, triangles, data,
168                               alpha=0.0, expand_search=False, #verbose = True, #False,
169                               max_points_per_cell = 4)
170        calc = interp.fit_points(point_attributes, )
171        #print "calc",calc
172       
173        # the point at 4,4 is ignored.  An expanded search has to be done
174        # to fine which triangel it's in.
175        # An expanded search isn't done to find that the last point
176        # isn't in the mesh.  But this isn't tested.
177        answer= [ 10,
178                     10,
179                     10,
180                     10,
181                    10,
182                     10,
183                     10,
184                     10,
185                     10,
186                     10]
187        assert allclose(calc, answer)
188
189    def test_quad_treeII(self):
190        p0 = [-66.0, 14.0]
191        p1 = [14.0, -66.0]
192        p2 = [14.0, 14.0]
193        p3 = [60.0, 20.0]
194        p4 = [10.0, 60.0]
195        p5 = [60.0, 60.0]
196       
197        points = [p0, p1, p2, p3, p4, p5]
198        triangles = [ [0, 1, 2],
199                      [3, 2, 1],
200                      [0, 2, 4],
201                      [0, 2, 4],
202                      [4, 2, 5],
203                      [5, 2, 3]]   
204
205        data = [ [-26.0,-26.0] ]
206        interp = Interpolation(points, triangles, data, alpha = 0.0,
207                 max_points_per_cell = 4)
208        #print "PDSG - interp.get_A()", interp.get_A()
209        answer =  [ [ 0.5,  0.5,  0.0,  0.,
210                      0., 0.]]
211        assert allclose(interp.get_A(), answer)
212        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
213        #print "PDSG -30,-30 - interp.get_A()", interp.get_A()
214        answer =  [ [ 0.0,  0.0,  0.0,  0.,
215                      0., 0.]]
216        assert allclose(interp.get_A(), answer)
217
218       
219        #point outside of quad tree root cell
220        interp.set_point_coordinates([[-70, -70]])
221        #print "PDSG -70,-70 interp.get_A()", interp.get_A()
222        answer =  [ [ 0.0,  0.0,  0.0,  0.,
223                      0., 0. ]]
224        assert allclose(interp.get_A(), answer)
225           
226
227    def test_datapoints_at_vertices(self):
228        """Test that data points coinciding with vertices yield a diagonal matrix
229        """
230       
231        a = [0.0, 0.0]
232        b = [0.0, 2.0]
233        c = [2.0,0.0]
234        points = [a, b, c]
235        vertices = [ [1,0,2] ]   #bac
236
237        data = points #Use data at vertices       
238       
239        interp = Interpolation(points, vertices, data)
240        assert allclose(interp.get_A(), [[1., 0., 0.],
241                                   [0., 1., 0.],
242                                   [0., 0., 1.]]) 
243       
244
245
246    def test_datapoints_on_edge_midpoints(self):
247        """Try datapoints midway on edges -
248        each point should affect two matrix entries equally
249        """
250       
251        a = [0.0, 0.0]
252        b = [0.0, 2.0]
253        c = [2.0,0.0]
254        points = [a, b, c]
255        vertices = [ [1,0,2] ]   #bac
256
257        data = [ [0., 1.], [1., 0.], [1., 1.] ]
258       
259        interp = Interpolation(points, vertices, data)
260       
261        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
262                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
263                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
264
265
266    def test_datapoints_on_edges(self):
267        """Try datapoints on edges -
268        each point should affect two matrix entries in proportion
269        """
270       
271        a = [0.0, 0.0]
272        b = [0.0, 2.0]
273        c = [2.0,0.0]
274        points = [a, b, c]
275        vertices = [ [1,0,2] ]   #bac
276
277        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
278       
279        interp = Interpolation(points, vertices, data)
280       
281        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
282                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
283                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
284
285    def test_arbitrary_datapoints(self):
286        """Try arbitrary datapoints
287        """
288
289        from Numeric import sum
290       
291        a = [0.0, 0.0]
292        b = [0.0, 2.0]
293        c = [2.0,0.0]
294        points = [a, b, c]
295        vertices = [ [1,0,2] ]   #bac
296
297        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
298       
299        interp = Interpolation(points, vertices, data)
300        #print "interp.get_A()", interp.get_A()
301        assert allclose(sum(interp.get_A(), axis=1), 1.0)
302       
303    def test_arbitrary_datapoints_some_outside(self):
304        """Try arbitrary datapoints one outside the triangle.
305        That one should be ignored
306        """
307
308        from Numeric import sum
309       
310        a = [0.0, 0.0]
311        b = [0.0, 2.0]
312        c = [2.0,0.0]
313        points = [a, b, c]
314        vertices = [ [1,0,2] ]   #bac
315
316        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
317
318
319        interp = Interpolation(points, vertices, data, precrop = True)
320        assert allclose(sum(interp.get_A(), axis=1), 1.0)
321
322        interp = Interpolation(points, vertices, data, precrop = False)
323        assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
324       
325       
326
327    # this causes a memory error in scipy.sparse       
328    def test_more_triangles(self):
329       
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        points = [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        #Data points
341        data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
342        interp = Interpolation(points, triangles, data_points) 
343
344        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d     
345                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
346                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
347                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
348                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
349                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
350
351
352        A = interp.get_A()
353        for i in range(A.shape[0]):
354            for j in range(A.shape[1]):
355                if not allclose(A[i,j], answer[i][j]):
356                    print i,j,':',A[i,j], answer[i][j] 
357
358
359        assert allclose(interp.get_A(), answer)
360
361
362
363
364    def test_smooth_attributes_to_mesh(self):
365        a = [0.0, 0.0]
366        b = [0.0, 5.0]
367        c = [5.0, 0.0]
368        points = [a, b, c]
369        triangles = [ [1,0,2] ]   #bac
370         
371        d1 = [1.0, 1.0]
372        d2 = [1.0, 3.0]
373        d3 = [3.0,1.0]
374        z1 = 2
375        z2 = 4
376        z3 = 4
377        data_coords = [d1, d2, d3] 
378       
379        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
380        z = [z1, z2, z3]
381        f = interp.fit(z)
382        answer = [0, 5., 5.]
383
384        #print "f\n",f
385        #print "answer\n",answer
386       
387        assert allclose(f, answer, atol=1e-7)
388       
389       
390    def test_smooth_att_to_meshII(self):
391       
392        a = [0.0, 0.0]
393        b = [0.0, 5.0]
394        c = [5.0, 0.0]
395        points = [a, b, c]
396        triangles = [ [1,0,2] ]   #bac
397         
398        d1 = [1.0, 1.0]
399        d2 = [1.0, 2.0]
400        d3 = [3.0,1.0]
401        data_coords = [d1, d2, d3]       
402        z = linear_function(data_coords)
403        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
404        f = interp.fit(z)
405        answer = linear_function(points)
406       
407        assert allclose(f, answer)
408   
409    def test_smooth_attributes_to_meshIII(self):
410       
411        a = [-1.0, 0.0]
412        b = [3.0, 4.0]
413        c = [4.0,1.0]
414        d = [-3.0, 2.0] #3
415        e = [-1.0,-2.0]
416        f = [1.0, -2.0] #5
417
418        vertices = [a, b, c, d,e,f]
419        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
420
421        point_coords = [[-2.0, 2.0],
422                        [-1.0, 1.0],
423                        [0.0,2.0],
424                        [1.0, 1.0],
425                        [2.0, 1.0],
426                        [0.0,0.0],
427                        [1.0, 0.0],
428                        [0.0, -1.0],
429                        [-0.2,-0.5],
430                        [-0.9, -1.5],
431                        [0.5, -1.9],
432                        [3.0,1.0]]
433
434        z = linear_function(point_coords)
435        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
436
437        #print 'z',z
438        f = interp.fit(z)
439        answer = linear_function(vertices)
440        #print "f\n",f
441        #print "answer\n",answer
442        assert allclose(f, answer)
443       
444
445    def test_smooth_attributes_to_meshIV(self):
446        """ Testing 2 attributes smoothed to the mesh
447        """
448       
449        a = [0.0, 0.0]
450        b = [0.0, 5.0]
451        c = [5.0, 0.0]
452        points = [a, b, c]
453        triangles = [ [1,0,2] ]   #bac
454         
455        d1 = [1.0, 1.0]
456        d2 = [1.0, 3.0]
457        d3 = [3.0, 1.0]
458        z1 = [2, 4]
459        z2 = [4, 8]
460        z3 = [4, 8]
461        data_coords = [d1, d2, d3] 
462       
463        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
464        z = [z1, z2, z3]
465        f =  interp.fit_points(z)
466        answer = [[0,0], [5., 10.], [5., 10.]]
467        assert allclose(f, answer)
468       
469    def test_interpolate_attributes_to_points(self):
470        v0 = [0.0, 0.0]
471        v1 = [0.0, 5.0]
472        v2 = [5.0, 0.0]
473       
474        vertices = [v0, v1, v2]
475        triangles = [ [1,0,2] ]   #bac
476         
477        d0 = [1.0, 1.0]
478        d1 = [1.0, 2.0]
479        d2 = [3.0, 1.0]
480        point_coords = [ d0, d1, d2]
481       
482        interp = Interpolation(vertices, triangles, point_coords)
483        f = linear_function(vertices)
484        z = interp.interpolate(f)
485        answer = linear_function(point_coords)
486       
487
488        assert allclose(z, answer)
489
490       
491    def test_interpolate_attributes_to_pointsII(self):
492        a = [-1.0, 0.0]
493        b = [3.0, 4.0]
494        c = [4.0, 1.0]
495        d = [-3.0, 2.0] #3
496        e = [-1.0, -2.0]
497        f = [1.0, -2.0] #5
498
499        vertices = [a, b, c, d,e,f]
500        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
501
502         
503        point_coords = [[-2.0, 2.0],
504                        [-1.0, 1.0],
505                        [0.0, 2.0],
506                        [1.0, 1.0],
507                        [2.0, 1.0],
508                        [0.0, 0.0],
509                        [1.0, 0.0],
510                        [0.0, -1.0],
511                        [-0.2, -0.5],
512                        [-0.9, -1.5],
513                        [0.5, -1.9],
514                        [3.0, 1.0]]
515       
516        interp = Interpolation(vertices, triangles, point_coords)
517        f = linear_function(vertices)
518        z = interp.interpolate(f)
519        answer = linear_function(point_coords)
520        #print "z",z
521        #print "answer",answer
522        assert allclose(z, answer)
523   
524    def test_interpolate_attributes_to_pointsIII(self):
525        """Test linear interpolation of known values at vertices to
526        new points inside a triangle
527        """
528        a = [0.0, 0.0]
529        b = [0.0, 5.0]
530        c = [5.0, 0.0]
531        d = [5.0, 5.0]
532       
533        vertices = [a, b, c, d]
534        triangles = [ [1,0,2], [2,3,0] ]   #bac, cdb
535
536        #Points within triangle 1
537        d0 = [1.0, 1.0]
538        d1 = [1.0, 2.0]
539        d2 = [3.0, 1.0]
540
541        #Point within triangle 2
542        d3 = [4.0, 3.0]
543
544        #Points on common edge
545        d4 = [2.5, 2.5]
546        d5 = [4.0, 1.0]
547         
548        #Point on common vertex
549        d6 = [0., 5.]
550
551       
552        point_coords = [d0, d1, d2, d3, d4, d5, d6]
553       
554        interp = Interpolation(vertices, triangles, point_coords)
555
556        #Known values at vertices
557        #Functions are x+y, x+2y, 2x+y, x-y-5       
558        f = [ [0., 0., 0., -5.],        # (0,0)   
559              [5., 10., 5., -10.],      # (0,5)
560              [5., 5., 10.0, 0.],       # (5,0)
561              [10., 15., 15., -5.]]     # (5,5)
562         
563        z = interp.interpolate(f)
564        answer = [ [2., 3., 3., -5.],   # (1,1)
565                   [3., 5., 4., -6.],   # (1,2)
566                   [4., 5., 7., -3.],   # (3,1)
567                   [7., 10., 11., -4.], # (4,3)
568                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
569                   [5., 6., 9., -2.],   # (4,1)
570                   [5., 10., 5., -10.]] # (0,5)
571
572        #print "***********"
573        #print "z",z
574        #print "answer",answer
575        #print "***********"
576
577        assert allclose(z, answer)
578   
579    def test_interpolate_attributes_to_pointsIV(self):
580        a = [-1.0, 0.0]
581        b = [3.0, 4.0]
582        c = [4.0, 1.0]
583        d = [-3.0, 2.0] #3
584        e = [-1.0, -2.0]
585        f = [1.0, -2.0] #5
586
587        vertices = [a, b, c, d,e,f]
588        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
589
590         
591        point_coords = [[-2.0, 2.0],
592                        [-1.0, 1.0],
593                        [0.0, 2.0],
594                        [1.0, 1.0],
595                        [2.0, 1.0],
596                        [0.0, 0.0],
597                        [1.0, 0.0],
598                        [0.0, -1.0],
599                        [-0.2, -0.5],
600                        [-0.9, -1.5],
601                        [0.5, -1.9],
602                        [3.0, 1.0]]
603       
604        interp = Interpolation(vertices, triangles, point_coords)
605        f = array([linear_function(vertices),2*linear_function(vertices) ])
606        f = transpose(f)
607        #print "f",f
608        z = interp.interpolate(f)
609        answer = [linear_function(point_coords),             
610                  2*linear_function(point_coords) ]
611        answer = transpose(answer)
612        #print "z",z
613        #print "answer",answer
614        assert allclose(z, answer)
615       
616    def test_smooth_attributes_to_mesh_function(self):
617        """ Testing 2 attributes smoothed to the mesh
618        """
619       
620        a = [0.0, 0.0]
621        b = [0.0, 5.0]
622        c = [5.0, 0.0]
623        points = [a, b, c]
624        triangles = [ [1,0,2] ]   #bac
625         
626        d1 = [1.0, 1.0]
627        d2 = [1.0, 3.0]
628        d3 = [3.0, 1.0]
629        z1 = [2, 4]
630        z2 = [4, 8]
631        z3 = [4, 8]
632        data_coords = [d1, d2, d3]
633        z = [z1, z2, z3]
634       
635        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
636        answer = [[0, 0], [5., 10.], [5., 10.]]
637       
638        assert allclose(f, answer)
639
640
641       
642    def test_pts2rectangular(self):
643
644        import time, os
645        FN = 'xyatest' + str(time.time()) + '.xya'
646        fid = open(FN, 'w')
647        fid.write(' %s \n' %('elevation'))       
648        fid.write('%f %f %f\n' %(1,1,2) )
649        fid.write('%f %f %f\n' %(1,3,4) )
650        fid.write('%f %f %f\n' %(3,1,4) )               
651        fid.close()
652       
653        points, triangles, boundary, attributes =\
654                pts2rectangular(FN, 4, 4, format = 'asc')
655       
656
657        data_coords = [ [1,1], [1,3], [3,1] ]
658        z = [2, 4, 4]
659
660        ref = fit_to_mesh(points, triangles, data_coords, z)           
661
662        #print attributes
663        #print ref
664        assert allclose(attributes, ref)
665
666        os.remove(FN)
667       
668
669    #Tests of smoothing matrix   
670    def test_smoothing_matrix_one_triangle(self):
671        from Numeric import dot
672        a = [0.0, 0.0]
673        b = [0.0, 2.0]
674        c = [2.0,0.0]
675        points = [a, b, c]
676       
677        vertices = [ [1,0,2] ]   #bac
678
679        interp = Interpolation(points, vertices)
680
681        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
682                                   [-0.5, 0.5, 0],
683                                   [-0.5, 0, 0.5]])
684
685        #Define f(x,y) = x
686        f = array([0,0,2]) #Value at global vertex 2
687
688        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
689        #           int 1 dx dy = area = 2
690        assert dot(dot(f, interp.get_D()), f) == 2
691       
692        #Define f(x,y) = y
693        f = array([0,2,0])  #Value at global vertex 1
694
695        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
696        #           int 1 dx dy = area = 2
697        assert dot(dot(f, interp.get_D()), f) == 2
698
699        #Define f(x,y) = x+y
700        f = array([0,2,2])  #Values at global vertex 1 and 2
701
702        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
703        #           int 2 dx dy = 2*area = 4
704        assert dot(dot(f, interp.get_D()), f) == 4       
705       
706
707
708    def test_smoothing_matrix_more_triangles(self):
709        from Numeric import dot
710       
711        a = [0.0, 0.0]
712        b = [0.0, 2.0]
713        c = [2.0,0.0]
714        d = [0.0, 4.0]
715        e = [2.0, 2.0]
716        f = [4.0,0.0]
717
718        points = [a, b, c, d, e, f]
719        #bac, bce, ecf, dbe, daf, dae
720        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
721
722        interp = Interpolation(points, vertices)
723
724
725        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
726        #                           [-0.5, 0.5, 0],
727        #                           [-0.5, 0, 0.5]])
728
729        #Define f(x,y) = x
730        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
731
732        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
733        #           int 1 dx dy = total area = 8
734        assert dot(dot(f, interp.get_D()), f) == 8
735       
736        #Define f(x,y) = y
737        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
738
739        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
740        #           int 1 dx dy = area = 8
741        assert dot(dot(f, interp.get_D()), f) == 8
742
743        #Define f(x,y) = x+y
744        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
745
746        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
747        #           int 2 dx dy = 2*area = 16
748        assert dot(dot(f, interp.get_D()), f) == 16       
749
750
751    def test_fit_and_interpolation(self):
752        from mesh import Mesh
753       
754        a = [0.0, 0.0]
755        b = [0.0, 2.0]
756        c = [2.0, 0.0]
757        d = [0.0, 4.0]
758        e = [2.0, 2.0]
759        f = [4.0, 0.0]
760
761        points = [a, b, c, d, e, f]
762        #bac, bce, ecf, dbe, daf, dae
763        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
764
765        #Get (enough) datapoints
766        data_points = [[ 0.66666667, 0.66666667],
767                       [ 1.33333333, 1.33333333],
768                       [ 2.66666667, 0.66666667],
769                       [ 0.66666667, 2.66666667],
770                       [ 0.0, 1.0],
771                       [ 0.0, 3.0],
772                       [ 1.0, 0.0],
773                       [ 1.0, 1.0],
774                       [ 1.0, 2.0],
775                       [ 1.0, 3.0],                                         
776                       [ 2.0, 1.0],
777                       [ 3.0, 0.0],
778                       [ 3.0, 1.0]]                                         
779       
780        interp = Interpolation(points, triangles, data_points, alpha=0.0)
781
782        z = linear_function(data_points)
783        answer = linear_function(points)
784
785        f = interp.fit(z)
786
787        #print "f",f
788        #print "answer",answer
789        assert allclose(f, answer)
790
791        #Map back
792        z1 = interp.interpolate(f)
793        #print "z1\n", z1
794        #print "z\n",z
795        assert allclose(z, z1)
796
797
798    def test_smoothing_and_interpolation(self):
799       
800        a = [0.0, 0.0]
801        b = [0.0, 2.0]
802        c = [2.0, 0.0]
803        d = [0.0, 4.0]
804        e = [2.0, 2.0]
805        f = [4.0, 0.0]
806
807        points = [a, b, c, d, e, f]
808        #bac, bce, ecf, dbe, daf, dae
809        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
810
811        #Get (too few!) datapoints
812        data_points = [[ 0.66666667, 0.66666667],
813                       [ 1.33333333, 1.33333333],
814                       [ 2.66666667, 0.66666667],
815                       [ 0.66666667, 2.66666667]]
816
817        z = linear_function(data_points)
818        answer = linear_function(points)       
819
820        #Make interpolator with too few data points and no smoothing
821        interp = Interpolation(points, triangles, data_points, alpha=0.0)
822        #Must raise an exception       
823        try:
824            f = interp.fit(z)
825        except:
826            pass
827
828        #Now try with smoothing parameter
829        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
830       
831        f = interp.fit(z)
832        #f will be different from answerr due to smoothing
833        assert allclose(f, answer,atol=5)
834
835        #Map back
836        z1 = interp.interpolate(f)
837        assert allclose(z, z1)
838
839
840
841    def test_fit_and_interpolation_with_new_points(self):
842        """Fit a surface to one set of points. Then interpolate that surface
843        using another set of points.
844        """
845        from mesh import Mesh
846
847
848        #Setup mesh used to represent fitted function
849        a = [0.0, 0.0]
850        b = [0.0, 2.0]
851        c = [2.0, 0.0]
852        d = [0.0, 4.0]
853        e = [2.0, 2.0]
854        f = [4.0, 0.0]
855
856        points = [a, b, c, d, e, f]
857        #bac, bce, ecf, dbe, daf, dae
858        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
859
860        #Datapoints to fit from
861        data_points1 = [[ 0.66666667, 0.66666667],
862                        [ 1.33333333, 1.33333333],
863                        [ 2.66666667, 0.66666667],
864                        [ 0.66666667, 2.66666667],
865                        [ 0.0, 1.0],
866                        [ 0.0, 3.0],
867                        [ 1.0, 0.0],
868                        [ 1.0, 1.0],
869                        [ 15, -17],   #Outside mesh
870                        [ 1.0, 2.0],
871                        [ 1.0, 3.0],                                         
872                        [ 2.0, 1.0],
873                        [ 3.0, 0.0],
874                        [ 3.0, 1.0]]
875
876        #Fit surface to mesh
877        interp = Interpolation(points, triangles, data_points1, alpha=0.0,
878                               precrop = True)
879        z = linear_function(data_points1) #Example z-values
880        f = interp.fit(z)                 #Fitted values at vertices
881
882
883
884        #New datapoints where interpolated values are sought
885        data_points2 = [[ 0.0, 0.0],
886                        [ 0.5, 0.5],
887                        [ 0.7, 0.7],
888                        [-13, 65],  #Outside
889                        [ 1.0, 0.5],
890                        [ 2.0, 0.4],
891                        [ 2.8, 1.2]]                                         
892       
893
894
895        #Build new A matrix based on new points (without precrop)
896        interp.build_interpolation_matrix_A(data_points2, precrop = False)
897
898        #Interpolate using fitted surface
899        z1 = interp.interpolate(f)
900
901        #import Numeric
902        #data_points2 = Numeric.take(data_points2, interp.point_indices)
903
904        #Desired result (OK for points inside)
905
906        answer = linear_function(data_points2)
907        import Numeric
908        z1 = Numeric.take(z1, [0,1,2,4,5,6])
909        answer = Numeric.take(answer, [0,1,2,4,5,6])               
910        assert allclose(z1, answer)
911
912        #Build new A matrix based on new points (with precrop)
913        interp.build_interpolation_matrix_A(data_points2, precrop = True)
914
915        #Interpolate using fitted surface
916        z1 = interp.interpolate(f)
917
918        import Numeric
919        data_points2 = Numeric.take(data_points2, interp.point_indices)
920
921        #Desired result
922        answer = linear_function(data_points2)
923        assert allclose(z1, answer)
924
925
926
927    def test_fit_and_interpolation_with_different_origins(self):
928        """Fit a surface to one set of points. Then interpolate that surface
929        using another set of points.
930        This test tests situtaion where points and mesh belong to a different
931        coordinate system as defined by origin.
932        """
933        from mesh import Mesh
934
935        #Setup mesh used to represent fitted function
936        a = [0.0, 0.0]
937        b = [0.0, 2.0]
938        c = [2.0, 0.0]
939        d = [0.0, 4.0]
940        e = [2.0, 2.0]
941        f = [4.0, 0.0]
942
943        points = [a, b, c, d, e, f]
944        #bac, bce, ecf, dbe, daf, dae
945        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
946
947        #Datapoints to fit from
948        data_points1 = [[ 0.66666667, 0.66666667],
949                        [ 1.33333333, 1.33333333],
950                        [ 2.66666667, 0.66666667],
951                        [ 0.66666667, 2.66666667],
952                        [ 0.0, 1.0],
953                        [ 0.0, 3.0],
954                        [ 1.0, 0.0],
955                        [ 1.0, 1.0],
956                        [ 1.0, 2.0],
957                        [ 1.0, 3.0],                                         
958                        [ 2.0, 1.0],
959                        [ 3.0, 0.0],
960                        [ 3.0, 1.0]]
961
962
963        #First check that things are OK when using same origin
964        mesh_origin = (56, 290000, 618000) #zone, easting, northing
965        data_origin = (56, 290000, 618000) #zone, easting, northing
966
967
968        #Fit surface to mesh
969        interp = Interpolation(points, triangles, data_points1,
970                               alpha=0.0,
971                               data_origin = data_origin,
972                               mesh_origin = mesh_origin)
973       
974        z = linear_function(data_points1) #Example z-values
975        f = interp.fit(z)                 #Fitted values at vertices
976
977
978        #New datapoints where interpolated values are sought
979        data_points2 = [[ 0.0, 0.0],
980                        [ 0.5, 0.5],
981                        [ 0.7, 0.7],
982                        [ 1.0, 0.5],
983                        [ 2.0, 0.4],
984                        [ 2.8, 1.2]]                                         
985       
986
987        #Build new A matrix based on new points
988        interp.build_interpolation_matrix_A(data_points2)
989
990        #Interpolate using fitted surface
991        z1 = interp.interpolate(f)
992
993        #Desired result
994        answer = linear_function(data_points2)
995        assert allclose(z1, answer)
996
997
998        ##############################################
999
1000        #Then check situtaion where points are relative to a different
1001        #origin (same zone, though, until we figure that out (FIXME))
1002       
1003        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1004        data_origin = (56, 10000, 10000) #zone, easting, northing
1005
1006        #Shift datapoints according to new origin
1007
1008        for k in range(len(data_points1)):
1009            data_points1[k][0] += mesh_origin[1] - data_origin[1]
1010            data_points1[k][1] += mesh_origin[2] - data_origin[2]
1011
1012        for k in range(len(data_points2)):
1013            data_points2[k][0] += mesh_origin[1] - data_origin[1]
1014            data_points2[k][1] += mesh_origin[2] - data_origin[2]           
1015
1016
1017       
1018        #Fit surface to mesh
1019        interp = Interpolation(points, triangles, data_points1,
1020                               alpha=0.0,
1021                               data_origin = data_origin,
1022                               mesh_origin = mesh_origin)
1023       
1024        f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
1025
1026        assert allclose(f,f1), 'Fit should have been unaltered'
1027
1028
1029        #Build new A matrix based on new points
1030        interp.build_interpolation_matrix_A(data_points2)
1031
1032        #Interpolate using fitted surface
1033        z1 = interp.interpolate(f)
1034        assert allclose(z1, answer)
1035
1036
1037        #########################################################
1038        #Finally try to relate data_points2 to new origin without
1039        #rebuilding matrix
1040
1041        data_origin = (56, 2000, 2000) #zone, easting, northing
1042        for k in range(len(data_points2)):
1043            data_points2[k][0] += 8000
1044            data_points2[k][1] += 8000
1045
1046        #Build new A matrix based on new points
1047        interp.build_interpolation_matrix_A(data_points2,
1048                                            data_origin = data_origin)
1049
1050        #Interpolate using fitted surface
1051        z1 = interp.interpolate(f)
1052        assert allclose(z1, answer)
1053
1054
1055
1056
1057    def test_fit_to_mesh_file(self):
1058        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
1059             export_triangulation_file   
1060        import tempfile
1061        import os
1062       
1063        # create a .tsh file, no user outline
1064        mesh_dic = {}
1065        mesh_dic['vertices'] = [[0.0, 0.0],
1066                                          [0.0, 5.0],
1067                                          [5.0, 0.0]]
1068        mesh_dic['triangles'] =  [[0, 2, 1]]
1069        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1070        mesh_dic['triangle_tags'] = [['']]
1071        mesh_dic['vertex_attributes'] = [[], [], []]
1072        mesh_dic['vertiex_attribute_titles'] = []
1073        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1074        mesh_dic['segment_tags'] = ['external',
1075                                                  'external',
1076                                                  'external']       
1077        mesh_file = tempfile.mktemp(".tsh")
1078        export_triangulation_file(mesh_file,mesh_dic)
1079       
1080        # create an .xya file
1081        point_file = tempfile.mktemp(".xya")
1082        fd = open(point_file,'w')
1083        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")
1084        fd.close()
1085       
1086        mesh_output_file = "new_trianlge.tsh"
1087        fit_to_mesh_file(mesh_file,
1088                         point_file,
1089                         mesh_output_file,
1090                         alpha = 0.0) 
1091        # load in the .tsh file we just wrote
1092        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
1093        #print "mesh_dic",mesh_dic
1094        ans =[[0.0, 0.0],
1095              [5.0, 10.0],
1096              [5.0,10.0]]
1097        assert allclose(mesh_dic['vertex_attributes'],ans)
1098       
1099        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1100                        ['elevation','stage'],
1101                        'test_fit_to_mesh_file failed')
1102       
1103        #clean up
1104        os.remove(mesh_file)
1105        os.remove(point_file)
1106        os.remove(mesh_output_file)       
1107       
1108    def test_fit_to_mesh_fileII(self):
1109        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
1110             export_triangulation_file   
1111        import tempfile
1112        import os
1113       
1114        # create a .tsh file, no user outline
1115        mesh_dic = {}
1116        mesh_dic['vertices'] = [[0.0, 0.0],
1117                                          [0.0, 5.0],
1118                                          [5.0, 0.0]]
1119        mesh_dic['triangles'] =  [[0, 2, 1]]
1120        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1121        mesh_dic['triangle_tags'] = [['']]
1122        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1123        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1124        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1125        mesh_dic['segment_tags'] = ['external',
1126                                                  'external',
1127                                                  'external']       
1128        mesh_file = tempfile.mktemp(".tsh")
1129        export_triangulation_file(mesh_file,mesh_dic)
1130       
1131        # create an .xya file
1132        point_file = tempfile.mktemp(".xya")
1133        fd = open(point_file,'w')
1134        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")
1135        fd.close()
1136       
1137        mesh_output_file = "new_triangle.tsh"
1138        fit_to_mesh_file(mesh_file,
1139                         point_file,
1140                         mesh_output_file,
1141                         alpha = 0.0) 
1142        # load in the .tsh file we just wrote
1143        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
1144
1145        assert allclose(mesh_dic['vertex_attributes'],
1146                        [[1.0, 2.0,0.0, 0.0],
1147                         [1.0, 2.0,5.0, 10.0],
1148                         [1.0, 2.0,5.0,10.0]])
1149       
1150        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1151                        ['density', 'temp','elevation','stage'],
1152                        'test_fit_to_mesh_file failed')
1153       
1154        #clean up
1155        os.remove(mesh_file)
1156        os.remove(mesh_output_file)       
1157        os.remove(point_file)
1158       
1159#-------------------------------------------------------------
1160if __name__ == "__main__":
1161    suite = unittest.makeSuite(TestCase,'test')
1162
1163    #suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
1164    runner = unittest.TextTestRunner(verbosity=1)
1165    runner.run(suite)
1166
1167   
1168   
1169
1170
1171
Note: See TracBrowser for help on using the repository browser.