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

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

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