source: inundation/pyvolution/test_least_squares.py @ 2891

Last change on this file since 2891 was 2778, checked in by steve, 19 years ago

On our 64 bit machine, ran into problem in pmesh/mesh.py where seg[0], seg[1], triangle,
and neighbor where not seen as integers. Had to wrap them in int(..)
File size: 58.9 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4import sys
5import unittest
6from math import sqrt
7
8
9from pyvolution.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       
40        interp = Interpolation(points, vertices, data)
41        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
42
43
44    def test_quad_tree(self):
45        p0 = [-10.0, -10.0]
46        p1 = [20.0, -10.0]
47        p2 = [-10.0, 20.0]
48        p3 = [10.0, 50.0]
49        p4 = [30.0, 30.0]
50        p5 = [50.0, 10.0]
51        p6 = [40.0, 60.0]
52        p7 = [60.0, 40.0]
53        p8 = [-66.0, 20.0]
54        p9 = [10.0, -66.0]
55
56        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
57        triangles = [ [0, 1, 2],
58                      [3, 2, 4],
59                      [4, 2, 1],
60                      [4, 1, 5],
61                      [3, 4, 6],
62                      [6, 4, 7],
63                      [7, 4, 5],
64                      [8, 0, 2],
65                      [0, 9, 1]]
66
67        data = [ [4,4] ]
68        interp = Interpolation(points, triangles, data, alpha = 0.0,
69                               max_points_per_cell = 4)
70        #print "PDSG - interp.get_A()", interp.get_A()
71        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
72                      0., 0. , 0., 0., 0., 0.]]
73        assert allclose(interp.get_A(), answer)
74        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
75        #print "PDSG - interp.get_A()", interp.get_A()
76        answer =  [ [ 0.0,  0.0,  0.0,  0.,
77                      0., 0. , 0., 0., 0., 0.]]
78        assert allclose(interp.get_A(), answer)
79
80
81        #point outside of quad tree root cell
82        interp.set_point_coordinates([[-70, -70]])
83        #print "PDSG - interp.get_A()", interp.get_A()
84        answer =  [ [ 0.0,  0.0,  0.0,  0.,
85                      0., 0. , 0., 0., 0., 0.]]
86        assert allclose(interp.get_A(), answer)
87
88    def test_expand_search(self):
89        p0 = [-10.0, -10.0]
90        p1 = [20.0, -10.0]
91        p2 = [-10.0, 20.0]
92        p3 = [10.0, 50.0]
93        p4 = [30.0, 30.0]
94        p5 = [50.0, 10.0]
95        p6 = [40.0, 60.0]
96        p7 = [60.0, 40.0]
97        p8 = [-66.0, 20.0]
98        p9 = [10.0, -66.0]
99
100        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
101        triangles = [ [0, 1, 2],
102                      [3, 2, 4],
103                      [4, 2, 1],
104                      [4, 1, 5],
105                      [3, 4, 6],
106                      [6, 4, 7],
107                      [7, 4, 5],
108                      [8, 0, 2],
109                      [0, 9, 1]]
110
111        data = [ [4,4],
112                 [-30,10],
113                 [-20,0],
114                 [-20,10],
115                 [0,30],
116                 [10,-40],
117                 [10,-30],
118                 [10,-20],
119                 [10,10],
120                 [10,20],
121                 [10,30],
122                 [10,40],
123                 [20,10],
124                 [25,45],
125                 [30,0],
126                 [30,10],
127                 [30,30],
128                 [30,40],
129                 [30,50],
130                 [40,10],
131                 [40,30],
132                 [40,40],
133                 [40,50],
134                 [50,20],
135                 [50,30],
136                 [50,40],
137                 [50,50],
138                 [30,0],
139                 [-20,-20]]
140        point_attributes = [ -400000,
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                   10,
168                   99]
169
170        interp = Interpolation(points, triangles, data,
171                               alpha=0.0, expand_search=False, #verbose = True, #False,
172                               max_points_per_cell = 4)
173        calc = interp.fit_points(point_attributes, )
174        #print "calc",calc
175
176        # the point at 4,4 is ignored.  An expanded search has to be done
177        # to fine which triangel it's in.
178        # An expanded search isn't done to find that the last point
179        # isn't in the mesh.  But this isn't tested.
180        answer= [ 10,
181                     10,
182                     10,
183                     10,
184                    10,
185                     10,
186                     10,
187                     10,
188                     10,
189                     10]
190        assert allclose(calc, answer)
191
192    def test_quad_treeII(self):
193        p0 = [-66.0, 14.0]
194        p1 = [14.0, -66.0]
195        p2 = [14.0, 14.0]
196        p3 = [60.0, 20.0]
197        p4 = [10.0, 60.0]
198        p5 = [60.0, 60.0]
199
200        points = [p0, p1, p2, p3, p4, p5]
201        triangles = [ [0, 1, 2],
202                      [3, 2, 1],
203                      [0, 2, 4],
204                      [4, 2, 5],
205                      [5, 2, 3]]
206
207        data = [ [-26.0,-26.0] ]
208        interp = Interpolation(points, triangles, data, alpha = 0.0,
209                 max_points_per_cell = 4)
210        #print "PDSG - interp.get_A()", interp.get_A()
211        answer =  [ [ 0.5,  0.5,  0.0,  0.,
212                      0., 0.]]
213        assert allclose(interp.get_A(), answer)
214        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
215        #print "PDSG -30,-30 - interp.get_A()", interp.get_A()
216        answer =  [ [ 0.0,  0.0,  0.0,  0.,
217                      0., 0.]]
218        assert allclose(interp.get_A(), answer)
219
220
221        #point outside of quad tree root cell
222        interp.set_point_coordinates([[-70, -70]])
223        #print "PDSG -70,-70 interp.get_A()", interp.get_A()
224        answer =  [ [ 0.0,  0.0,  0.0,  0.,
225                      0., 0. ]]
226        assert allclose(interp.get_A(), answer)
227
228
229    def test_datapoints_at_vertices(self):
230        """Test that data points coinciding with vertices yield a diagonal matrix
231        """
232
233        a = [0.0, 0.0]
234        b = [0.0, 2.0]
235        c = [2.0,0.0]
236        points = [a, b, c]
237        vertices = [ [1,0,2] ]   #bac
238
239        data = points #Use data at vertices
240
241        interp = Interpolation(points, vertices, data)
242        assert allclose(interp.get_A(), [[1., 0., 0.],
243                                   [0., 1., 0.],
244                                   [0., 0., 1.]])
245
246
247
248    def test_datapoints_on_edge_midpoints(self):
249        """Try datapoints midway on edges -
250        each point should affect two matrix entries equally
251        """
252
253        a = [0.0, 0.0]
254        b = [0.0, 2.0]
255        c = [2.0,0.0]
256        points = [a, b, c]
257        vertices = [ [1,0,2] ]   #bac
258
259        data = [ [0., 1.], [1., 0.], [1., 1.] ]
260
261        interp = Interpolation(points, vertices, data)
262
263        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
264                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
265                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
266
267
268    def test_datapoints_on_edges(self):
269        """Try datapoints on edges -
270        each point should affect two matrix entries in proportion
271        """
272
273        a = [0.0, 0.0]
274        b = [0.0, 2.0]
275        c = [2.0,0.0]
276        points = [a, b, c]
277        vertices = [ [1,0,2] ]   #bac
278
279        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
280
281        interp = Interpolation(points, vertices, data)
282
283        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
284                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
285                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
286
287    def test_arbitrary_datapoints(self):
288        """Try arbitrary datapoints
289        """
290
291        from Numeric import sum
292
293        a = [0.0, 0.0]
294        b = [0.0, 2.0]
295        c = [2.0,0.0]
296        points = [a, b, c]
297        vertices = [ [1,0,2] ]   #bac
298
299        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
300
301        interp = Interpolation(points, vertices, data)
302        #print "interp.get_A()", interp.get_A()
303        assert allclose(sum(interp.get_A(), axis=1), 1.0)
304
305    def test_arbitrary_datapoints_some_outside(self):
306        """Try arbitrary datapoints one outside the triangle.
307        That one should be ignored
308        """
309
310        from Numeric import sum
311
312        a = [0.0, 0.0]
313        b = [0.0, 2.0]
314        c = [2.0,0.0]
315        points = [a, b, c]
316        vertices = [ [1,0,2] ]   #bac
317
318        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
319
320
321        interp = Interpolation(points, vertices, data, precrop = True)
322        assert allclose(sum(interp.get_A(), axis=1), 1.0)
323
324        interp = Interpolation(points, vertices, data, precrop = False)
325        assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
326
327
328
329    # this causes a memory error in scipy.sparse
330    def test_more_triangles(self):
331
332        a = [-1.0, 0.0]
333        b = [3.0, 4.0]
334        c = [4.0,1.0]
335        d = [-3.0, 2.0] #3
336        e = [-1.0,-2.0]
337        f = [1.0, -2.0] #5
338
339        points = [a, b, c, d,e,f]
340        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
341
342        #Data points
343        data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
344        interp = Interpolation(points, triangles, data_points)
345
346        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
347                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
348                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
349                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
350                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
351                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
352
353
354        A = interp.get_A()
355        for i in range(A.shape[0]):
356            for j in range(A.shape[1]):
357                if not allclose(A[i,j], answer[i][j]):
358                    print i,j,':',A[i,j], answer[i][j]
359
360
361        assert allclose(interp.get_A(), answer)
362
363
364
365
366    def test_smooth_attributes_to_mesh(self):
367        a = [0.0, 0.0]
368        b = [0.0, 5.0]
369        c = [5.0, 0.0]
370        points = [a, b, c]
371        triangles = [ [1,0,2] ]   #bac
372
373        d1 = [1.0, 1.0]
374        d2 = [1.0, 3.0]
375        d3 = [3.0,1.0]
376        z1 = 2
377        z2 = 4
378        z3 = 4
379        data_coords = [d1, d2, d3]
380
381        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
382        z = [z1, z2, z3]
383        f = interp.fit(z)
384        answer = [0, 5., 5.]
385
386        #print "f\n",f
387        #print "answer\n",answer
388
389        assert allclose(f, answer, atol=1e-7)
390
391
392    def test_smooth_att_to_meshII(self):
393
394        a = [0.0, 0.0]
395        b = [0.0, 5.0]
396        c = [5.0, 0.0]
397        points = [a, b, c]
398        triangles = [ [1,0,2] ]   #bac
399
400        d1 = [1.0, 1.0]
401        d2 = [1.0, 2.0]
402        d3 = [3.0,1.0]
403        data_coords = [d1, d2, d3]
404        z = linear_function(data_coords)
405        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
406        f = interp.fit(z)
407        answer = linear_function(points)
408
409        assert allclose(f, answer)
410
411    def test_smooth_attributes_to_meshIII(self):
412
413        a = [-1.0, 0.0]
414        b = [3.0, 4.0]
415        c = [4.0,1.0]
416        d = [-3.0, 2.0] #3
417        e = [-1.0,-2.0]
418        f = [1.0, -2.0] #5
419
420        vertices = [a, b, c, d,e,f]
421        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
422
423        point_coords = [[-2.0, 2.0],
424                        [-1.0, 1.0],
425                        [0.0,2.0],
426                        [1.0, 1.0],
427                        [2.0, 1.0],
428                        [0.0,0.0],
429                        [1.0, 0.0],
430                        [0.0, -1.0],
431                        [-0.2,-0.5],
432                        [-0.9, -1.5],
433                        [0.5, -1.9],
434                        [3.0,1.0]]
435
436        z = linear_function(point_coords)
437        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
438
439        #print 'z',z
440        f = interp.fit(z)
441        answer = linear_function(vertices)
442        #print "f\n",f
443        #print "answer\n",answer
444        assert allclose(f, answer)
445
446
447    def test_smooth_attributes_to_meshIV(self):
448        """ Testing 2 attributes smoothed to the mesh
449        """
450
451        a = [0.0, 0.0]
452        b = [0.0, 5.0]
453        c = [5.0, 0.0]
454        points = [a, b, c]
455        triangles = [ [1,0,2] ]   #bac
456
457        d1 = [1.0, 1.0]
458        d2 = [1.0, 3.0]
459        d3 = [3.0, 1.0]
460        z1 = [2, 4]
461        z2 = [4, 8]
462        z3 = [4, 8]
463        data_coords = [d1, d2, d3]
464
465        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
466        z = [z1, z2, z3]
467        f =  interp.fit_points(z)
468        answer = [[0,0], [5., 10.], [5., 10.]]
469        assert allclose(f, answer)
470
471    def test_interpolate_attributes_to_points(self):
472        v0 = [0.0, 0.0]
473        v1 = [0.0, 5.0]
474        v2 = [5.0, 0.0]
475
476        vertices = [v0, v1, v2]
477        triangles = [ [1,0,2] ]   #bac
478
479        d0 = [1.0, 1.0]
480        d1 = [1.0, 2.0]
481        d2 = [3.0, 1.0]
482        point_coords = [ d0, d1, d2]
483
484        interp = Interpolation(vertices, triangles, point_coords)
485        f = linear_function(vertices)
486        #z = interp.interpolate(f)
487        answer = linear_function(point_coords)
488
489
490        #assert allclose(z, answer)
491
492
493    def test_interpolate_attributes_to_points_interp_only(self):
494        v0 = [0.0, 0.0]
495        v1 = [0.0, 5.0]
496        v2 = [5.0, 0.0]
497
498        vertices = [v0, v1, v2]
499        triangles = [ [1,0,2] ]   #bac
500
501        d0 = [1.0, 1.0]
502        d1 = [1.0, 2.0]
503        d2 = [3.0, 1.0]
504        point_coords = [ d0, d1, d2]
505
506        interp = Interpolation(vertices, triangles, point_coords,
507                               interp_only = True)
508       
509        f = linear_function(vertices)
510        #z = interp.interpolate(f)
511        answer = linear_function(point_coords)
512        #print "answer", answer
513        #print "z", z
514
515        #assert allclose(z, answer)
516
517    def test_interpolate_attributes_to_pointsII(self):
518        a = [-1.0, 0.0]
519        b = [3.0, 4.0]
520        c = [4.0, 1.0]
521        d = [-3.0, 2.0] #3
522        e = [-1.0, -2.0]
523        f = [1.0, -2.0] #5
524
525        vertices = [a, b, c, d,e,f]
526        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
527
528
529        point_coords = [[-2.0, 2.0],
530                        [-1.0, 1.0],
531                        [0.0, 2.0],
532                        [1.0, 1.0],
533                        [2.0, 1.0],
534                        [0.0, 0.0],
535                        [1.0, 0.0],
536                        [0.0, -1.0],
537                        [-0.2, -0.5],
538                        [-0.9, -1.5],
539                        [0.5, -1.9],
540                        [3.0, 1.0]]
541
542        interp = Interpolation(vertices, triangles, point_coords)
543        f = linear_function(vertices)
544        #z = interp.interpolate(f)
545        answer = linear_function(point_coords)
546        #print "z",z
547        #print "answer",answer
548        #assert allclose(z, answer)
549
550    def test_interpolate_attributes_to_pointsIII(self):
551        """Test linear interpolation of known values at vertices to
552        new points inside a triangle
553        """
554        a = [0.0, 0.0]
555        b = [0.0, 5.0]
556        c = [5.0, 0.0]
557        d = [5.0, 5.0]
558
559        vertices = [a, b, c, d]
560        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
561
562        #Points within triangle 1
563        d0 = [1.0, 1.0]
564        d1 = [1.0, 2.0]
565        d2 = [3.0, 1.0]
566
567        #Point within triangle 2
568        d3 = [4.0, 3.0]
569
570        #Points on common edge
571        d4 = [2.5, 2.5]
572        d5 = [4.0, 1.0]
573
574        #Point on common vertex
575        d6 = [0., 5.]
576       
577        point_coords = [d0, d1, d2, d3, d4, d5, d6]
578
579        interp = Interpolation(vertices, triangles, point_coords)
580
581        #Known values at vertices
582        #Functions are x+y, x+2y, 2x+y, x-y-5
583        f = [ [0., 0., 0., -5.],        # (0,0)
584              [5., 10., 5., -10.],      # (0,5)
585              [5., 5., 10.0, 0.],       # (5,0)
586              [10., 15., 15., -5.]]     # (5,5)
587
588        #z = interp.interpolate(f)
589        answer = [ [2., 3., 3., -5.],   # (1,1)
590                   [3., 5., 4., -6.],   # (1,2)
591                   [4., 5., 7., -3.],   # (3,1)
592                   [7., 10., 11., -4.], # (4,3)
593                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
594                   [5., 6., 9., -2.],   # (4,1)
595                   [5., 10., 5., -10.]]  # (0,5)
596
597        #print "***********"
598        #print "z",z
599        #print "answer",answer
600        #print "***********"
601
602        #Should an error message be returned if points are outside
603        # of the mesh? Not currently. 
604
605        #assert allclose(z, answer)
606
607
608    def test_interpolate_point_outside_of_mesh(self):
609        """Test linear interpolation of known values at vertices to
610        new points inside a triangle
611        """
612        a = [0.0, 0.0]
613        b = [0.0, 5.0]
614        c = [5.0, 0.0]
615        d = [5.0, 5.0]
616
617        vertices = [a, b, c, d]
618        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
619
620        #Far away point
621        d7 = [-1., -1.]
622       
623        point_coords = [ d7]
624
625        interp = Interpolation(vertices, triangles, point_coords)
626
627        #Known values at vertices
628        #Functions are x+y, x+2y, 2x+y, x-y-5
629        f = [ [0., 0., 0., -5.],        # (0,0)
630              [5., 10., 5., -10.],      # (0,5)
631              [5., 5., 10.0, 0.],       # (5,0)
632              [10., 15., 15., -5.]]     # (5,5)
633
634
635        #z = interp.interpolate(f)
636        #answer = [ [0., 0., 0., 0.]] # (-1,-1)
637
638        #print "***********"
639        #print "z",z
640        #print "answer",answer
641        #print "***********"
642
643        #Should an error message be returned if points are outside
644        # of the mesh? Not currently. 
645
646        #assert allclose(z, answer)
647
648    def test_interpolate_attributes_to_pointsIV(self):
649        a = [-1.0, 0.0]
650        b = [3.0, 4.0]
651        c = [4.0, 1.0]
652        d = [-3.0, 2.0] #3
653        e = [-1.0, -2.0]
654        f = [1.0, -2.0] #5
655
656        vertices = [a, b, c, d,e,f]
657        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
658
659
660        point_coords = [[-2.0, 2.0],
661                        [-1.0, 1.0],
662                        [0.0, 2.0],
663                        [1.0, 1.0],
664                        [2.0, 1.0],
665                        [0.0, 0.0],
666                        [1.0, 0.0],
667                        [0.0, -1.0],
668                        [-0.2, -0.5],
669                        [-0.9, -1.5],
670                        [0.5, -1.9],
671                        [3.0, 1.0]]
672
673        interp = Interpolation(vertices, triangles, point_coords)
674        f = array([linear_function(vertices),2*linear_function(vertices) ])
675        f = transpose(f)
676        #print "f",f
677        #z = interp.interpolate(f)
678        answer = [linear_function(point_coords),
679                  2*linear_function(point_coords) ]
680        answer = transpose(answer)
681        #print "z",z
682        #print "answer",answer
683        #assert allclose(z, answer)
684
685    def test_smooth_attributes_to_mesh_function(self):
686        """ Testing 2 attributes smoothed to the mesh
687        """
688
689        a = [0.0, 0.0]
690        b = [0.0, 5.0]
691        c = [5.0, 0.0]
692        points = [a, b, c]
693        triangles = [ [1,0,2] ]   #bac
694
695        d1 = [1.0, 1.0]
696        d2 = [1.0, 3.0]
697        d3 = [3.0, 1.0]
698        z1 = [2, 4]
699        z2 = [4, 8]
700        z3 = [4, 8]
701        data_coords = [d1, d2, d3]
702        z = [z1, z2, z3]
703
704        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
705        answer = [[0, 0], [5., 10.], [5., 10.]]
706
707        assert allclose(f, answer)
708
709
710
711    def test_pts2rectangular(self):
712
713        import time, os
714        FN = 'xyatest' + str(time.time()) + '.xya'
715        fid = open(FN, 'w')
716        fid.write(' %s \n' %('elevation'))
717        fid.write('%f %f %f\n' %(1,1,2) )
718        fid.write('%f %f %f\n' %(1,3,4) )
719        fid.write('%f %f %f\n' %(3,1,4) )
720        fid.close()
721
722        points, triangles, boundary, attributes =\
723                pts2rectangular(FN, 4, 4)
724
725
726        data_coords = [ [1,1], [1,3], [3,1] ]
727        z = [2, 4, 4]
728
729        ref = fit_to_mesh(points, triangles, data_coords, z, verbose=False)
730
731        #print attributes
732        #print ref
733        assert allclose(attributes, ref)
734
735        os.remove(FN)
736
737
738    #Tests of smoothing matrix
739    def test_smoothing_matrix_one_triangle(self):
740        from Numeric import dot
741        a = [0.0, 0.0]
742        b = [0.0, 2.0]
743        c = [2.0,0.0]
744        points = [a, b, c]
745
746        vertices = [ [1,0,2] ]   #bac
747
748        interp = Interpolation(points, vertices)
749
750        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
751                                   [-0.5, 0.5, 0],
752                                   [-0.5, 0, 0.5]])
753
754        #Define f(x,y) = x
755        f = array([0,0,2]) #Value at global vertex 2
756
757        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
758        #           int 1 dx dy = area = 2
759        assert dot(dot(f, interp.get_D()), f) == 2
760
761        #Define f(x,y) = y
762        f = array([0,2,0])  #Value at global vertex 1
763
764        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
765        #           int 1 dx dy = area = 2
766        assert dot(dot(f, interp.get_D()), f) == 2
767
768        #Define f(x,y) = x+y
769        f = array([0,2,2])  #Values at global vertex 1 and 2
770
771        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
772        #           int 2 dx dy = 2*area = 4
773        assert dot(dot(f, interp.get_D()), f) == 4
774
775
776
777    def test_smoothing_matrix_more_triangles(self):
778        from Numeric import dot
779
780        a = [0.0, 0.0]
781        b = [0.0, 2.0]
782        c = [2.0,0.0]
783        d = [0.0, 4.0]
784        e = [2.0, 2.0]
785        f = [4.0,0.0]
786
787        points = [a, b, c, d, e, f]
788        #bac, bce, ecf, dbe, daf, dae
789        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
790
791        interp = Interpolation(points, vertices)
792
793
794        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
795        #                           [-0.5, 0.5, 0],
796        #                           [-0.5, 0, 0.5]])
797
798        #Define f(x,y) = x
799        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
800
801        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
802        #           int 1 dx dy = total area = 8
803        assert dot(dot(f, interp.get_D()), f) == 8
804
805        #Define f(x,y) = y
806        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
807
808        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
809        #           int 1 dx dy = area = 8
810        assert dot(dot(f, interp.get_D()), f) == 8
811
812        #Define f(x,y) = x+y
813        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
814
815        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
816        #           int 2 dx dy = 2*area = 16
817        assert dot(dot(f, interp.get_D()), f) == 16
818
819
820    def test_fit_and_interpolation(self):
821        from mesh import Mesh
822
823        a = [0.0, 0.0]
824        b = [0.0, 2.0]
825        c = [2.0, 0.0]
826        d = [0.0, 4.0]
827        e = [2.0, 2.0]
828        f = [4.0, 0.0]
829
830        points = [a, b, c, d, e, f]
831        #bac, bce, ecf, dbe, daf, dae
832        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
833
834        #Get (enough) datapoints
835        data_points = [[ 0.66666667, 0.66666667],
836                       [ 1.33333333, 1.33333333],
837                       [ 2.66666667, 0.66666667],
838                       [ 0.66666667, 2.66666667],
839                       [ 0.0, 1.0],
840                       [ 0.0, 3.0],
841                       [ 1.0, 0.0],
842                       [ 1.0, 1.0],
843                       [ 1.0, 2.0],
844                       [ 1.0, 3.0],
845                       [ 2.0, 1.0],
846                       [ 3.0, 0.0],
847                       [ 3.0, 1.0]]
848
849        interp = Interpolation(points, triangles, data_points, alpha=0.0)
850
851        z = linear_function(data_points)
852        answer = linear_function(points)
853
854        f = interp.fit(z)
855
856        #print "f",f
857        #print "answer",answer
858        assert allclose(f, answer)
859
860        #Map back
861        #z1 = interp.interpolate(f)
862        #print "z1\n", z1
863        #print "z\n",z
864        #assert allclose(z, z1)
865
866
867    def test_smoothing_and_interpolation(self):
868
869        a = [0.0, 0.0]
870        b = [0.0, 2.0]
871        c = [2.0, 0.0]
872        d = [0.0, 4.0]
873        e = [2.0, 2.0]
874        f = [4.0, 0.0]
875
876        points = [a, b, c, d, e, f]
877        #bac, bce, ecf, dbe, daf, dae
878        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
879
880        #Get (too few!) datapoints
881        data_points = [[ 0.66666667, 0.66666667],
882                       [ 1.33333333, 1.33333333],
883                       [ 2.66666667, 0.66666667],
884                       [ 0.66666667, 2.66666667]]
885
886        z = linear_function(data_points)
887        answer = linear_function(points)
888
889        #Make interpolator with too few data points and no smoothing
890        interp = Interpolation(points, triangles, data_points, alpha=0.0)
891        #Must raise an exception
892        try:
893            f = interp.fit(z)
894        except:
895            pass
896
897        #Now try with smoothing parameter
898        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
899
900        f = interp.fit(z)
901        #f will be different from answer due to smoothing
902        assert allclose(f, answer,atol=5)
903
904        #Map back
905        #z1 = interp.interpolate(f)
906        #assert allclose(z, z1)
907
908
909
910    def test_fit_and_interpolation_with_new_points(self):
911        """Fit a surface to one set of points. Then interpolate that surface
912        using another set of points.
913        """
914        from mesh import Mesh
915
916
917        #Setup mesh used to represent fitted function
918        a = [0.0, 0.0]
919        b = [0.0, 2.0]
920        c = [2.0, 0.0]
921        d = [0.0, 4.0]
922        e = [2.0, 2.0]
923        f = [4.0, 0.0]
924
925        points = [a, b, c, d, e, f]
926        #bac, bce, ecf, dbe, daf, dae
927        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
928
929        #Datapoints to fit from
930        data_points1 = [[ 0.66666667, 0.66666667],
931                        [ 1.33333333, 1.33333333],
932                        [ 2.66666667, 0.66666667],
933                        [ 0.66666667, 2.66666667],
934                        [ 0.0, 1.0],
935                        [ 0.0, 3.0],
936                        [ 1.0, 0.0],
937                        [ 1.0, 1.0],
938                        [ 15, -17],   #Outside mesh
939                        [ 1.0, 2.0],
940                        [ 1.0, 3.0],
941                        [ 2.0, 1.0],
942                        [ 3.0, 0.0],
943                        [ 3.0, 1.0]]
944
945        #Fit surface to mesh
946        interp = Interpolation(points, triangles, data_points1, alpha=0.0,
947                               precrop = True, verbose=False)
948        z = linear_function(data_points1) #Example z-values
949        f = interp.fit(z)                 #Fitted values at vertices
950
951
952
953        #New datapoints where interpolated values are sought
954        data_points2 = [[ 0.0, 0.0],
955                        [ 0.5, 0.5],
956                        [ 0.7, 0.7],
957                        [-13, 65],  #Outside
958                        [ 1.0, 0.5],
959                        [ 2.0, 0.4],
960                        [ 2.8, 1.2]]
961
962
963
964        #Build new A matrix based on new points (without precrop)
965        interp.build_interpolation_matrix_A(data_points2, precrop = False)
966
967        #Interpolate using fitted surface
968        #z1 = interp.interpolate(f)
969
970        #import Numeric
971        #data_points2 = Numeric.take(data_points2, interp.point_indices)
972
973        #Desired result (OK for points inside)
974
975        answer = linear_function(data_points2)
976        import Numeric
977        #z1 = Numeric.take(z1, [0,1,2,4,5,6])
978        answer = Numeric.take(answer, [0,1,2,4,5,6])
979        #assert allclose(z1, answer)
980
981        #Build new A matrix based on new points (with precrop)
982        interp.build_interpolation_matrix_A(data_points2, precrop = True)
983
984        #Interpolate using fitted surface
985        #z1 = interp.interpolate(f)
986
987        import Numeric
988        data_points2 = Numeric.take(data_points2, interp.point_indices)
989
990        #Desired result
991        answer = linear_function(data_points2)
992        #assert allclose(z1, answer)
993
994
995
996
997
998
999    def test_interpolation_from_discontinuous_vertex_values(self):
1000        """test_interpolation_from_discontinuous_vertex_values.
1001        This will test the format used internally in pyvolution and also
1002        interpolation from sww files
1003        """
1004       
1005        from mesh import Mesh
1006
1007
1008        #Setup mesh used to represent discontinuous function
1009        a = [0.0, 0.0]
1010        b = [0.0, 2.0]
1011        c = [2.0, 0.0]
1012        d = [0.0, 4.0]
1013        e = [2.0, 2.0]
1014        f = [4.0, 0.0]
1015
1016        points = [b, a, c,
1017                  b, c, e,
1018                  e, c, f,
1019                  d, b, e]
1020       
1021        #bac, bce, ecf, dbe
1022        triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
1023
1024       
1025        vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
1026                         
1027       
1028
1029        #New datapoints where interpolated values are sought
1030        data_points = [[0.0, 0.0],  #T0
1031                       [0.5, 0.5],  #T0
1032                       [1.5, 1.5],  #T1
1033                       [2.5, 0.5],  #T2
1034                       [0.0, 3.0],  #T3
1035                       [1.0, 2.0],  #In between T1 and T3 (T1 is used) FIXME?
1036                       [2.0, 1.0],  #In between T1 and T2 (T1 is used) FIXME?
1037                       [1.0, 1.0]]  #In between T1 and T0 (T0 is used) FIXME?
1038       
1039
1040
1041
1042        #Build interpolation matrix
1043        interp = Interpolation(points, triangles, data_points)
1044                               #, alpha=0.0, precrop = True)
1045
1046        #print interp.A.todense()
1047        #print vertex_values
1048
1049        #Interpolate using fitted surface
1050        #z = interp.interpolate(vertex_values)
1051
1052        #print z
1053
1054        #assert allclose(z, [0,0,1,2,5,1,1,0])
1055
1056
1057
1058
1059    def test_interpolation_function_time_only(self):
1060        """Test spatio-temporal interpolation
1061        Test that spatio temporal function performs the correct
1062        interpolations in both time and space
1063        """
1064
1065
1066        #Three timesteps
1067        time = [1.0, 5.0, 6.0]
1068       
1069
1070        #One quantity
1071        Q = zeros( (3,6), Float )
1072
1073        #Linear in time and space
1074        a = [0.0, 0.0]
1075        b = [0.0, 2.0]
1076        c = [2.0, 0.0]
1077        d = [0.0, 4.0]
1078        e = [2.0, 2.0]
1079        f = [4.0, 0.0]
1080
1081        points = [a, b, c, d, e, f]
1082       
1083        for i, t in enumerate(time):
1084            Q[i, :] = t*linear_function(points)
1085
1086           
1087        #Check basic interpolation of one quantity using averaging
1088        #(no interpolation points or spatial info)
1089        from utilities.numerical_tools import mean       
1090        I = Interpolation_function(time, [mean(Q[0,:]),
1091                                          mean(Q[1,:]),
1092                                          mean(Q[2,:])])
1093
1094
1095
1096        #Check temporal interpolation
1097        for i in [0,1,2]:
1098            assert allclose(I(time[i]), mean(Q[i,:]))
1099
1100        #Midway   
1101        assert allclose(I( (time[0] + time[1])/2 ),
1102                        (I(time[0]) + I(time[1]))/2 )
1103
1104        assert allclose(I( (time[1] + time[2])/2 ),
1105                        (I(time[1]) + I(time[2]))/2 )
1106
1107        assert allclose(I( (time[0] + time[2])/2 ),
1108                        (I(time[0]) + I(time[2]))/2 )                 
1109
1110        #1/3
1111        assert allclose(I( (time[0] + time[2])/3 ),
1112                        (I(time[0]) + I(time[2]))/3 )                         
1113
1114
1115        #Out of bounds checks
1116        try:
1117            I(time[0]-1) 
1118        except:
1119            pass
1120        else:
1121            raise 'Should raise exception'
1122
1123        try:
1124            I(time[-1]+1) 
1125        except:
1126            pass
1127        else:
1128            raise 'Should raise exception'       
1129
1130
1131       
1132
1133    def interpolation_test_interpolation_function_spatial_only(self):
1134        """Test spatio-temporal interpolation with constant time
1135        """
1136
1137        #Three timesteps
1138        time = [1.0, 5.0, 6.0]
1139       
1140       
1141        #Setup mesh used to represent fitted function
1142        a = [0.0, 0.0]
1143        b = [0.0, 2.0]
1144        c = [2.0, 0.0]
1145        d = [0.0, 4.0]
1146        e = [2.0, 2.0]
1147        f = [4.0, 0.0]
1148
1149        points = [a, b, c, d, e, f]
1150        #bac, bce, ecf, dbe
1151        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1152
1153
1154        #New datapoints where interpolated values are sought
1155        interpolation_points = [[ 0.0, 0.0],
1156                                [ 0.5, 0.5],
1157                                [ 0.7, 0.7],
1158                                [ 1.0, 0.5],
1159                                [ 2.0, 0.4],
1160                                [ 2.8, 1.2]]
1161
1162
1163        #One quantity linear in space
1164        Q = linear_function(points)
1165
1166
1167        #Check interpolation of one quantity using interpolaton points
1168        I = Interpolation_function(time, Q,
1169                                   vertex_coordinates = points,
1170                                   triangles = triangles, 
1171                                   interpolation_points = interpolation_points,
1172                                   verbose = False)
1173
1174
1175        answer = linear_function(interpolation_points)
1176
1177        t = time[0]
1178        for j in range(50): #t in [1, 6]
1179            for id in range(len(interpolation_points)):
1180                assert allclose(I(t, id), answer[id])
1181
1182            t += 0.1   
1183
1184
1185        try:   
1186            I(1)
1187        except:
1188            pass
1189        else:
1190            raise 'Should raise exception'
1191           
1192
1193
1194    def interpolation_test_interpolation_function(self):
1195        """Test spatio-temporal interpolation
1196        Test that spatio temporal function performs the correct
1197        interpolations in both time and space
1198        """
1199
1200
1201        #Three timesteps
1202        time = [1.0, 5.0, 6.0]
1203       
1204
1205        #Setup mesh used to represent fitted function
1206        a = [0.0, 0.0]
1207        b = [0.0, 2.0]
1208        c = [2.0, 0.0]
1209        d = [0.0, 4.0]
1210        e = [2.0, 2.0]
1211        f = [4.0, 0.0]
1212
1213        points = [a, b, c, d, e, f]
1214        #bac, bce, ecf, dbe
1215        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1216
1217
1218        #New datapoints where interpolated values are sought
1219        interpolation_points = [[ 0.0, 0.0],
1220                                [ 0.5, 0.5],
1221                                [ 0.7, 0.7],
1222                                [ 1.0, 0.5],
1223                                [ 2.0, 0.4],
1224                                [ 2.8, 1.2]]
1225
1226
1227        #One quantity
1228        Q = zeros( (3,6), Float )
1229
1230        #Linear in time and space
1231        for i, t in enumerate(time):
1232            Q[i, :] = t*linear_function(points)
1233
1234
1235        #Check interpolation of one quantity using interpolaton points)
1236        I = Interpolation_function(time, Q,
1237                                   vertex_coordinates = points,
1238                                   triangles = triangles, 
1239                                   interpolation_points = interpolation_points,
1240                                   verbose = False)
1241
1242
1243        answer = linear_function(interpolation_points)
1244
1245        t = time[0]
1246        for j in range(50): #t in [1, 6]
1247            for id in range(len(interpolation_points)):
1248                assert allclose(I(t, id), t*answer[id])
1249
1250            t += 0.1   
1251
1252
1253        try:   
1254            I(1)
1255        except:
1256            pass
1257        else:
1258            raise 'Should raise exception'
1259           
1260        #
1261        #interpolation_points = [[ 0.0, 0.0],
1262        #                        [ 0.5, 0.5],
1263        #                        [ 0.7, 0.7],
1264        #                        [-13, 65],  #Outside
1265        #                        [ 1.0, 0.5],
1266        #                        [ 2.0, 0.4],
1267        #                        [ 2.8, 1.2]]
1268        #
1269        #try:
1270        #    I = Interpolation_function(time, Q,
1271        #                               vertex_coordinates = points,
1272        #                               triangles = triangles,
1273        #                               interpolation_points = interpolation_points,
1274        #                               verbose = False)
1275        #except:
1276        #    pass
1277        #else:
1278        #    raise 'Should raise exception'
1279
1280
1281
1282
1283
1284    def test_fit_and_interpolation_with_different_origins(self):
1285        """Fit a surface to one set of points. Then interpolate that surface
1286        using another set of points.
1287        This test tests situtaion where points and mesh belong to a different
1288        coordinate system as defined by origin.
1289        """
1290        from mesh import Mesh
1291
1292        #Setup mesh used to represent fitted function
1293        a = [0.0, 0.0]
1294        b = [0.0, 2.0]
1295        c = [2.0, 0.0]
1296        d = [0.0, 4.0]
1297        e = [2.0, 2.0]
1298        f = [4.0, 0.0]
1299
1300        points = [a, b, c, d, e, f]
1301        #bac, bce, ecf, dbe, daf, dae
1302        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1303
1304        #Datapoints to fit from
1305        data_points1 = [[ 0.66666667, 0.66666667],
1306                        [ 1.33333333, 1.33333333],
1307                        [ 2.66666667, 0.66666667],
1308                        [ 0.66666667, 2.66666667],
1309                        [ 0.0, 1.0],
1310                        [ 0.0, 3.0],
1311                        [ 1.0, 0.0],
1312                        [ 1.0, 1.0],
1313                        [ 1.0, 2.0],
1314                        [ 1.0, 3.0],
1315                        [ 2.0, 1.0],
1316                        [ 3.0, 0.0],
1317                        [ 3.0, 1.0]]
1318
1319
1320        #First check that things are OK when using same origin
1321        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1322        data_origin = (56, 290000, 618000) #zone, easting, northing
1323
1324
1325        #Fit surface to mesh
1326        interp = Interpolation(points, triangles, data_points1,
1327                               alpha=0.0,
1328                               data_origin = data_origin,
1329                               mesh_origin = mesh_origin)
1330
1331        z = linear_function(data_points1) #Example z-values
1332        f = interp.fit(z)                 #Fitted values at vertices
1333
1334
1335        #New datapoints where interpolated values are sought
1336        data_points2 = [[ 0.0, 0.0],
1337                        [ 0.5, 0.5],
1338                        [ 0.7, 0.7],
1339                        [ 1.0, 0.5],
1340                        [ 2.0, 0.4],
1341                        [ 2.8, 1.2]]
1342
1343
1344        #Build new A matrix based on new points
1345        interp.build_interpolation_matrix_A(data_points2)
1346
1347        #Interpolate using fitted surface
1348        #z1 = interp.interpolate(f)
1349
1350        #Desired result
1351        #answer = linear_function(data_points2)
1352        #assert allclose(z1, answer)
1353
1354
1355        ##############################################
1356
1357        #Then check situation where points are relative to a different
1358        #origin (same zone, though, until we figure that out (FIXME))
1359
1360        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1361        data_origin = (56, 10000, 10000) #zone, easting, northing
1362
1363        #Shift datapoints according to new origin
1364
1365        for k in range(len(data_points1)):
1366            data_points1[k][0] += mesh_origin[1] - data_origin[1]
1367            data_points1[k][1] += mesh_origin[2] - data_origin[2]
1368
1369        for k in range(len(data_points2)):
1370            data_points2[k][0] += mesh_origin[1] - data_origin[1]
1371            data_points2[k][1] += mesh_origin[2] - data_origin[2]
1372
1373
1374
1375        #Fit surface to mesh
1376        interp = Interpolation(points, triangles, data_points1,
1377                               alpha=0.0,
1378                               data_origin = data_origin,
1379                               mesh_origin = mesh_origin)
1380
1381        f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
1382
1383        assert allclose(f,f1), 'Fit should have been unaltered'
1384
1385
1386        #Build new A matrix based on new points
1387        interp.build_interpolation_matrix_A(data_points2)
1388
1389        #Interpolate using fitted surface
1390        #z1 = interp.interpolate(f)
1391        #assert allclose(z1, answer)
1392
1393
1394        #########################################################
1395        #Finally try to relate data_points2 to new origin without
1396        #rebuilding matrix
1397
1398        data_origin = (56, 2000, 2000) #zone, easting, northing
1399        for k in range(len(data_points2)):
1400            data_points2[k][0] += 8000
1401            data_points2[k][1] += 8000
1402
1403        #Build new A matrix based on new points
1404        interp.build_interpolation_matrix_A(data_points2,
1405                                            data_origin = data_origin)
1406
1407        #Interpolate using fitted surface
1408        #z1 = interp.interpolate(f)
1409        #assert allclose(z1, answer)
1410
1411
1412
1413    def test_fit_to_mesh_w_georef(self):
1414        """Simple check that georef works at the fit_to_mesh level
1415        """
1416       
1417        from coordinate_transforms.geo_reference import Geo_reference
1418
1419        #Mesh
1420        vertex_coordinates = [[0.76, 0.76],
1421                              [0.76, 5.76],
1422                              [5.76, 0.76]]
1423        triangles = [[0,2,1]]
1424
1425        mesh_geo = Geo_reference(56,-0.76,-0.76)
1426
1427
1428        #Data                       
1429        data_points = [[ 201.0, 401.0],
1430                       [ 201.0, 403.0],
1431                       [ 203.0, 401.0]]
1432
1433        z = [2, 4, 4]
1434       
1435        data_geo = Geo_reference(56,-200,-400)
1436       
1437        #Fit
1438        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
1439                         data_origin = data_geo.get_origin(),
1440                         mesh_origin = mesh_geo.get_origin(),
1441                         alpha = 0)
1442        assert allclose( zz, [0,5,5] )
1443
1444
1445    def test_fit_to_mesh_file(self):
1446        from load_mesh.loadASCII import import_mesh_file, \
1447             export_mesh_file
1448        import tempfile
1449        import os
1450
1451        # create a .tsh file, no user outline
1452        mesh_dic = {}
1453        mesh_dic['vertices'] = [[0.0, 0.0],
1454                                          [0.0, 5.0],
1455                                          [5.0, 0.0]]
1456        mesh_dic['triangles'] =  [[0, 2, 1]]
1457        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1458        mesh_dic['triangle_tags'] = ['']
1459        mesh_dic['vertex_attributes'] = [[], [], []]
1460        mesh_dic['vertiex_attribute_titles'] = []
1461        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1462        mesh_dic['segment_tags'] = ['external',
1463                                                  'external',
1464                                                  'external']
1465        mesh_file = tempfile.mktemp(".tsh")
1466        export_mesh_file(mesh_file,mesh_dic)
1467
1468        # create an .xya file
1469        point_file = tempfile.mktemp(".xya")
1470        fd = open(point_file,'w')
1471        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")
1472        fd.close()
1473
1474        mesh_output_file = tempfile.mktemp(".tsh") 
1475        fit_to_mesh_file(mesh_file,
1476                         point_file,
1477                         mesh_output_file,
1478                         alpha = 0.0)
1479        # load in the .tsh file we just wrote
1480        mesh_dic = import_mesh_file(mesh_output_file)
1481        #print "mesh_dic",mesh_dic
1482        ans =[[0.0, 0.0],
1483              [5.0, 10.0],
1484              [5.0,10.0]]
1485        assert allclose(mesh_dic['vertex_attributes'],ans)
1486
1487        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1488                        ['elevation','stage'],
1489                        'test_fit_to_mesh_file failed')
1490
1491        #clean up
1492        os.remove(mesh_file)
1493        os.remove(point_file)
1494        os.remove(mesh_output_file)
1495
1496    def test_fit_to_mesh_file3(self):
1497        from load_mesh.loadASCII import import_mesh_file, \
1498             export_mesh_file
1499        import tempfile
1500        import os
1501
1502        # create a .tsh file, no user outline
1503        mesh_dic = {}
1504        mesh_dic['vertices'] = [[0.76, 0.76],
1505                                          [0.76, 5.76],
1506                                          [5.76, 0.76]]
1507        mesh_dic['triangles'] =  [[0, 2, 1]]
1508        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1509        mesh_dic['triangle_tags'] = ['']
1510        mesh_dic['vertex_attributes'] = [[], [], []]
1511        mesh_dic['vertiex_attribute_titles'] = []
1512        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1513        mesh_dic['segment_tags'] = ['external',
1514                                                  'external',
1515                                                  'external']
1516        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1517        mesh_file = tempfile.mktemp(".tsh")
1518        export_mesh_file(mesh_file,mesh_dic)
1519
1520        # create an .xya file
1521        point_file = tempfile.mktemp(".xya")
1522        fd = open(point_file,'w')
1523        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")
1524        fd.close()
1525
1526        mesh_output_file = tempfile.mktemp(".tsh")
1527        fit_to_mesh_file(mesh_file,
1528                         point_file,
1529                         mesh_output_file,
1530                         alpha = 0.0)
1531        # load in the .tsh file we just wrote
1532        mesh_dic = import_mesh_file(mesh_output_file)
1533        #print "mesh_dic",mesh_dic
1534        ans =[[0.0, 0.0],
1535              [5.0, 10.0],
1536              [5.0,10.0]]
1537        assert allclose(mesh_dic['vertex_attributes'],ans)
1538
1539        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1540                        ['elevation','stage'],
1541                        'test_fit_to_mesh_file failed')
1542
1543        #clean up
1544        os.remove(mesh_file)
1545        os.remove(point_file)
1546        os.remove(mesh_output_file)
1547
1548    def test_fit_to_mesh_file4(self):
1549        from load_mesh.loadASCII import import_mesh_file, \
1550             export_mesh_file
1551        import tempfile
1552        import os
1553
1554        # create a .tsh file, no user outline
1555        mesh_dic = {}
1556        mesh_dic['vertices'] = [[0.76, 0.76],
1557                                [0.76, 5.76],
1558                                [5.76, 0.76]]
1559        mesh_dic['triangles'] =  [[0, 2, 1]]
1560        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1561        mesh_dic['triangle_tags'] = ['']
1562        mesh_dic['vertex_attributes'] = [[], [], []]
1563        mesh_dic['vertiex_attribute_titles'] = []
1564        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1565        mesh_dic['segment_tags'] = ['external',
1566                                    'external',
1567                                    'external']
1568        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1569        mesh_file = tempfile.mktemp(".tsh")
1570        export_mesh_file(mesh_file,mesh_dic)
1571
1572        geo_ref = Geo_reference(56,-200,-400)
1573        # create an .xya file
1574        point_file = tempfile.mktemp(".xya")
1575        fd = open(point_file,'w')
1576        fd.write("elevation, stage \n 201.0, 401.0,2.,4 \n 201.0, 403.0,4,8 \n 203.0, 401.0,4.,8 \n")
1577        geo_ref.write_ASCII(fd)
1578        fd.close()
1579
1580        mesh_output_file = tempfile.mktemp(".tsh")
1581        fit_to_mesh_file(mesh_file,
1582                         point_file,
1583                         mesh_output_file,
1584                         alpha = 0.0)
1585        # load in the .tsh file we just wrote
1586        mesh_dic = import_mesh_file(mesh_output_file)
1587        #print "mesh_dic",mesh_dic
1588        ans =[[0.0, 0.0],
1589              [5.0, 10.0],
1590              [5.0, 10.0]]
1591        assert allclose(mesh_dic['vertex_attributes'],ans)
1592
1593        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1594                        ['elevation','stage'],
1595                        'test_fit_to_mesh_file failed')
1596
1597        #clean up
1598        os.remove(mesh_file)
1599        os.remove(point_file)
1600        os.remove(mesh_output_file)
1601
1602    def test_fit_to_mesh_fileII(self):
1603        from load_mesh.loadASCII import import_mesh_file, \
1604             export_mesh_file
1605        import tempfile
1606        import os
1607
1608        # create a .tsh file, no user outline
1609        mesh_dic = {}
1610        mesh_dic['vertices'] = [[0.0, 0.0],
1611                                [0.0, 5.0],
1612                                [5.0, 0.0]]
1613        mesh_dic['triangles'] =  [[0, 2, 1]]
1614        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1615        mesh_dic['triangle_tags'] = ['']
1616        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1617        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1618        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1619        mesh_dic['segment_tags'] = ['external',
1620                                                  'external',
1621                                                  'external']
1622        mesh_file = tempfile.mktemp(".tsh")
1623        export_mesh_file(mesh_file,mesh_dic)
1624
1625        # create an .xya file
1626        point_file = tempfile.mktemp(".xya")
1627        fd = open(point_file,'w')
1628        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")
1629        fd.close()
1630
1631        mesh_output_file = "new_triangle.tsh"
1632        fit_to_mesh_file(mesh_file,
1633                         point_file,
1634                         mesh_output_file,
1635                         alpha = 0.0)
1636        # load in the .tsh file we just wrote
1637        mesh_dic = import_mesh_file(mesh_output_file)
1638
1639        assert allclose(mesh_dic['vertex_attributes'],
1640                        [[1.0, 2.0,0.0, 0.0],
1641                         [1.0, 2.0,5.0, 10.0],
1642                         [1.0, 2.0,5.0,10.0]])
1643
1644        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1645                        ['density', 'temp','elevation','stage'],
1646                        'test_fit_to_mesh_file failed')
1647
1648        #clean up
1649        os.remove(mesh_file)
1650        os.remove(mesh_output_file)
1651        os.remove(point_file)
1652
1653    def test_fit_to_mesh_file_errors(self):
1654        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1655        import tempfile
1656        import os
1657
1658        # create a .tsh file, no user outline
1659        mesh_dic = {}
1660        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1661        mesh_dic['triangles'] =  [[0, 2, 1]]
1662        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1663        mesh_dic['triangle_tags'] = ['']
1664        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1665        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1666        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1667        mesh_dic['segment_tags'] = ['external', 'external','external']
1668        mesh_file = tempfile.mktemp(".tsh")
1669        export_mesh_file(mesh_file,mesh_dic)
1670
1671        # create an .xya file
1672        point_file = tempfile.mktemp(".xya")
1673        fd = open(point_file,'w')
1674        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")
1675        fd.close()
1676
1677        mesh_output_file = "new_triangle.tsh"
1678        try:
1679            fit_to_mesh_file(mesh_file, point_file,
1680                             mesh_output_file, display_errors = False)
1681        except IOError:
1682            pass
1683        else:
1684            #self.failUnless(0 ==1,  'Bad file did not raise error!')
1685            raise 'Bad file did not raise error!'
1686           
1687        #clean up
1688        os.remove(mesh_file)
1689        os.remove(point_file)
1690
1691    def test_fit_to_mesh_file_errorsII(self):
1692        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1693        import tempfile
1694        import os
1695
1696        # create a .tsh file, no user outline
1697        mesh_file = tempfile.mktemp(".tsh")
1698        fd = open(mesh_file,'w')
1699        fd.write("unit testing a bad .tsh file \n")
1700        fd.close()
1701
1702        # create an .xya file
1703        point_file = tempfile.mktemp(".xya")
1704        fd = open(point_file,'w')
1705        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")
1706        fd.close()
1707
1708        mesh_output_file = "new_triangle.tsh"
1709        try:
1710            fit_to_mesh_file(mesh_file, point_file,
1711                             mesh_output_file, display_errors = False)
1712        except IOError:
1713            pass
1714        else:
1715            raise 'Bad file did not raise error!'
1716           
1717        #clean up
1718        os.remove(mesh_file)
1719        os.remove(point_file)
1720
1721    def test_fit_to_mesh_file_errorsIII(self):
1722        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1723        import tempfile
1724        import os
1725
1726        # create a .tsh file, no user outline
1727        mesh_dic = {}
1728        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1729        mesh_dic['triangles'] =  [[0, 2, 1]]
1730        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1731        mesh_dic['triangle_tags'] = ['']
1732        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1733        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1734        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1735        mesh_dic['segment_tags'] = ['external', 'external','external']
1736        mesh_file = tempfile.mktemp(".tsh")
1737        export_mesh_file(mesh_file,mesh_dic)
1738
1739        # create an .xya file
1740        point_file = tempfile.mktemp(".xya")
1741        fd = open(point_file,'w')
1742        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")
1743        fd.close()
1744
1745        #This a deliberately illegal filename to invoke the error.
1746        mesh_output_file = ".../\z\z:ya.tsh"       
1747
1748        try:
1749            fit_to_mesh_file(mesh_file, point_file,
1750                             mesh_output_file, display_errors = False)
1751        except IOError:
1752            pass
1753        else:
1754            raise 'Bad file did not raise error!'
1755       
1756        #clean up
1757        os.remove(mesh_file)
1758        os.remove(point_file)
1759 
1760## FIXME?  Running from the Comand line isn't in vogue these days
1761#  The test was breaking when test_all at the inundation level was running
1762# was running it.issue - not running the test in this directory
1763    def Bad_test_fit_to_mesh_file_errorsIV(self):
1764        import os
1765        command = '%s least_squares.py q q q e n 0.9 n'  %(sys.executable)
1766        status = os.system(command) 
1767        self.failUnless(status%255  == 1,
1768                        'command prompt least_squares.py failed.  Incorect exit status.')
1769       
1770    def test_fit_to_msh_netcdf_fileII(self):
1771        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1772        import tempfile
1773        import os
1774
1775        # create a .tsh file, no user outline
1776        mesh_dic = {}
1777        mesh_dic['vertices'] = [[0.0, 0.0],
1778                                [0.0, 5.0],
1779                                [5.0, 0.0]]
1780        mesh_dic['triangles'] =  [[0, 2, 1]]
1781        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1782        mesh_dic['triangle_tags'] = ['']
1783        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1784        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1785        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1786        mesh_dic['segment_tags'] = ['external',
1787                                                  'external',
1788                                                  'external']
1789        mesh_file = tempfile.mktemp(".msh")
1790        export_mesh_file(mesh_file,mesh_dic)
1791
1792        # create an .xya file
1793        point_file = tempfile.mktemp(".xya")
1794        fd = open(point_file,'w')
1795        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")
1796        fd.close()
1797
1798        mesh_output_file = "new_triangle.msh"
1799        fit_to_mesh_file(mesh_file,
1800                         point_file,
1801                         mesh_output_file,
1802                         alpha = 0.0)
1803        # load in the .tsh file we just wrote
1804        mesh_dic = import_mesh_file(mesh_output_file)
1805
1806        assert allclose(mesh_dic['vertex_attributes'],
1807                        [[1.0, 2.0,0.0, 0.0],
1808                         [1.0, 2.0,5.0, 10.0],
1809                         [1.0, 2.0,5.0,10.0]])
1810
1811        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1812                        ['density', 'temp','elevation','stage'],
1813                        'test_fit_to_mesh_file failed')
1814
1815        #clean up
1816        os.remove(mesh_file)
1817        os.remove(mesh_output_file)
1818        os.remove(point_file)
1819
1820       
1821       
1822    def test_fit_using_fit_to_mesh(self):
1823        """Fit a surface to one set of points. Then interpolate that surface
1824        using another set of points.
1825        """
1826        from mesh import Mesh
1827
1828
1829        #Setup mesh used to represent fitted function
1830        a = [0.0, 0.0]
1831        b = [0.0, 2.0]
1832        c = [2.0, 0.0]
1833        d = [0.0, 4.0]
1834        e = [2.0, 2.0]
1835        f = [4.0, 0.0]
1836
1837        points = [a, b, c, d, e, f]
1838        #bac, bce, ecf, dbe, daf, dae
1839        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1840
1841        #Datapoints to fit from
1842        data_points1 = [[ 0.66666667, 0.66666667],
1843                        [ 1.33333333, 1.33333333],
1844                        [ 2.66666667, 0.66666667],
1845                        [ 0.66666667, 2.66666667],
1846                        [ 0.0, 1.0],
1847                        [ 0.0, 3.0],
1848                        [ 1.0, 0.0],
1849                        [ 1.0, 1.0],
1850                        [ 15, -17],   #Outside mesh
1851                        [ 1.0, 2.0],
1852                        [ 1.0, 3.0],
1853                        [ 2.0, 1.0],
1854                        [ 3.0, 0.0],
1855                        [ 3.0, 1.0]]
1856
1857        #Fit surface to mesh
1858        z = linear_function(data_points1) #Example z-values
1859        v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
1860                        precrop=True, verbose=False)
1861
1862        assert allclose(linear_function(points), v)
1863
1864
1865       
1866    def test_acceptable_overshoot(self):
1867        """Fit a surface to one set of points. Then interpolate that surface
1868        using another set of points.
1869        Check that exceedance in fitted values are caught.
1870        """
1871        from mesh import Mesh
1872
1873
1874        #Setup mesh used to represent fitted function
1875        a = [0.0, 0.0]
1876        b = [0.0, 2.0]
1877        c = [2.0, 0.0]
1878        d = [0.0, 4.0]
1879        e = [2.0, 2.0]
1880        f = [4.0, 0.0]
1881
1882        points = [a, b, c, d, e, f]
1883        #bac, bce, ecf, dbe, daf, dae
1884        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1885
1886        #Datapoints to fit from
1887        data_points1 = [[ 0.66666667, 0.66666667],
1888                        [ 1.33333333, 1.33333333],
1889                        [ 2.66666667, 0.66666667],
1890                        [ 0.66666667, 2.66666667],
1891                        [ 0.0, 1.0],
1892                        [ 0.0, 3.0],
1893                        [ 1.0, 0.0],
1894                        [ 1.0, 1.0],
1895                        [ 15, -17],   #Outside mesh
1896                        [ 1.0, 2.0],
1897                        [ 1.0, 3.0],
1898                        [ 2.0, 1.0],
1899                        [ 3.0, 0.0],
1900                        [ 3.0, 1.0]]
1901
1902        #Fit surface to mesh
1903        z = linear_function(data_points1) #Example z-values
1904
1905        try:
1906            v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
1907                            acceptable_overshoot = 0.2,
1908                            precrop=True, verbose=False)
1909        except FittingError, e:
1910            pass
1911        else:
1912            raise 'Should have raised exception'
1913           
1914
1915        #assert allclose(linear_function(points), v)
1916       
1917
1918       
1919#-------------------------------------------------------------
1920if __name__ == "__main__":
1921    #suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_function')
1922    #suite = unittest.makeSuite(Test_Least_Squares,'test_datapoint_at_centroid')
1923    suite = unittest.makeSuite(Test_Least_Squares,'test')
1924
1925    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
1926    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
1927    runner = unittest.TextTestRunner(verbosity=1)
1928    runner.run(suite)
1929
1930
1931
1932
1933
Note: See TracBrowser for help on using the repository browser.