source: inundation/pyvolution/test_least_squares.py @ 2267

Last change on this file since 2267 was 2262, checked in by ole, 19 years ago

Fixed missing geo reference in set_quantity and added tests
File size: 57.2 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4import sys
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_points_interp_only(self):
493        v0 = [0.0, 0.0]
494        v1 = [0.0, 5.0]
495        v2 = [5.0, 0.0]
496
497        vertices = [v0, v1, v2]
498        triangles = [ [1,0,2] ]   #bac
499
500        d0 = [1.0, 1.0]
501        d1 = [1.0, 2.0]
502        d2 = [3.0, 1.0]
503        point_coords = [ d0, d1, d2]
504
505        interp = Interpolation(vertices, triangles, point_coords,
506                               interp_only = True)
507       
508        f = linear_function(vertices)
509        z = interp.interpolate(f)
510        answer = linear_function(point_coords)
511        #print "answer", answer
512        #print "z", z
513
514        assert allclose(z, answer)
515
516    def test_interpolate_attributes_to_pointsII(self):
517        a = [-1.0, 0.0]
518        b = [3.0, 4.0]
519        c = [4.0, 1.0]
520        d = [-3.0, 2.0] #3
521        e = [-1.0, -2.0]
522        f = [1.0, -2.0] #5
523
524        vertices = [a, b, c, d,e,f]
525        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
526
527
528        point_coords = [[-2.0, 2.0],
529                        [-1.0, 1.0],
530                        [0.0, 2.0],
531                        [1.0, 1.0],
532                        [2.0, 1.0],
533                        [0.0, 0.0],
534                        [1.0, 0.0],
535                        [0.0, -1.0],
536                        [-0.2, -0.5],
537                        [-0.9, -1.5],
538                        [0.5, -1.9],
539                        [3.0, 1.0]]
540
541        interp = Interpolation(vertices, triangles, point_coords)
542        f = linear_function(vertices)
543        z = interp.interpolate(f)
544        answer = linear_function(point_coords)
545        #print "z",z
546        #print "answer",answer
547        assert allclose(z, answer)
548
549    def test_interpolate_attributes_to_pointsIII(self):
550        """Test linear interpolation of known values at vertices to
551        new points inside a triangle
552        """
553        a = [0.0, 0.0]
554        b = [0.0, 5.0]
555        c = [5.0, 0.0]
556        d = [5.0, 5.0]
557
558        vertices = [a, b, c, d]
559        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
560
561        #Points within triangle 1
562        d0 = [1.0, 1.0]
563        d1 = [1.0, 2.0]
564        d2 = [3.0, 1.0]
565
566        #Point within triangle 2
567        d3 = [4.0, 3.0]
568
569        #Points on common edge
570        d4 = [2.5, 2.5]
571        d5 = [4.0, 1.0]
572
573        #Point on common vertex
574        d6 = [0., 5.]
575       
576        point_coords = [d0, d1, d2, d3, d4, d5, d6]
577
578        interp = Interpolation(vertices, triangles, point_coords)
579
580        #Known values at vertices
581        #Functions are x+y, x+2y, 2x+y, x-y-5
582        f = [ [0., 0., 0., -5.],        # (0,0)
583              [5., 10., 5., -10.],      # (0,5)
584              [5., 5., 10.0, 0.],       # (5,0)
585              [10., 15., 15., -5.]]     # (5,5)
586
587        z = interp.interpolate(f)
588        answer = [ [2., 3., 3., -5.],   # (1,1)
589                   [3., 5., 4., -6.],   # (1,2)
590                   [4., 5., 7., -3.],   # (3,1)
591                   [7., 10., 11., -4.], # (4,3)
592                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
593                   [5., 6., 9., -2.],   # (4,1)
594                   [5., 10., 5., -10.]]  # (0,5)
595
596        #print "***********"
597        #print "z",z
598        #print "answer",answer
599        #print "***********"
600
601        #Should an error message be returned if points are outside
602        # of the mesh? Not currently. 
603
604        assert allclose(z, answer)
605
606
607    def test_interpolate_point_outside_of_mesh(self):
608        """Test linear interpolation of known values at vertices to
609        new points inside a triangle
610        """
611        a = [0.0, 0.0]
612        b = [0.0, 5.0]
613        c = [5.0, 0.0]
614        d = [5.0, 5.0]
615
616        vertices = [a, b, c, d]
617        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
618
619        #Far away point
620        d7 = [-1., -1.]
621       
622        point_coords = [ d7]
623
624        interp = Interpolation(vertices, triangles, point_coords)
625
626        #Known values at vertices
627        #Functions are x+y, x+2y, 2x+y, x-y-5
628        f = [ [0., 0., 0., -5.],        # (0,0)
629              [5., 10., 5., -10.],      # (0,5)
630              [5., 5., 10.0, 0.],       # (5,0)
631              [10., 15., 15., -5.]]     # (5,5)
632
633        z = interp.interpolate(f)
634        answer = [ [0., 0., 0., 0.]] # (-1,-1)
635
636        #print "***********"
637        #print "z",z
638        #print "answer",answer
639        #print "***********"
640
641        #Should an error message be returned if points are outside
642        # of the mesh? Not currently. 
643
644        assert allclose(z, answer)
645
646    def test_interpolate_attributes_to_pointsIV(self):
647        a = [-1.0, 0.0]
648        b = [3.0, 4.0]
649        c = [4.0, 1.0]
650        d = [-3.0, 2.0] #3
651        e = [-1.0, -2.0]
652        f = [1.0, -2.0] #5
653
654        vertices = [a, b, c, d,e,f]
655        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
656
657
658        point_coords = [[-2.0, 2.0],
659                        [-1.0, 1.0],
660                        [0.0, 2.0],
661                        [1.0, 1.0],
662                        [2.0, 1.0],
663                        [0.0, 0.0],
664                        [1.0, 0.0],
665                        [0.0, -1.0],
666                        [-0.2, -0.5],
667                        [-0.9, -1.5],
668                        [0.5, -1.9],
669                        [3.0, 1.0]]
670
671        interp = Interpolation(vertices, triangles, point_coords)
672        f = array([linear_function(vertices),2*linear_function(vertices) ])
673        f = transpose(f)
674        #print "f",f
675        z = interp.interpolate(f)
676        answer = [linear_function(point_coords),
677                  2*linear_function(point_coords) ]
678        answer = transpose(answer)
679        #print "z",z
680        #print "answer",answer
681        assert allclose(z, answer)
682
683    def test_smooth_attributes_to_mesh_function(self):
684        """ Testing 2 attributes smoothed to the mesh
685        """
686
687        a = [0.0, 0.0]
688        b = [0.0, 5.0]
689        c = [5.0, 0.0]
690        points = [a, b, c]
691        triangles = [ [1,0,2] ]   #bac
692
693        d1 = [1.0, 1.0]
694        d2 = [1.0, 3.0]
695        d3 = [3.0, 1.0]
696        z1 = [2, 4]
697        z2 = [4, 8]
698        z3 = [4, 8]
699        data_coords = [d1, d2, d3]
700        z = [z1, z2, z3]
701
702        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
703        answer = [[0, 0], [5., 10.], [5., 10.]]
704
705        assert allclose(f, answer)
706
707
708
709    def test_pts2rectangular(self):
710
711        import time, os
712        FN = 'xyatest' + str(time.time()) + '.xya'
713        fid = open(FN, 'w')
714        fid.write(' %s \n' %('elevation'))
715        fid.write('%f %f %f\n' %(1,1,2) )
716        fid.write('%f %f %f\n' %(1,3,4) )
717        fid.write('%f %f %f\n' %(3,1,4) )
718        fid.close()
719
720        points, triangles, boundary, attributes =\
721                pts2rectangular(FN, 4, 4)
722
723
724        data_coords = [ [1,1], [1,3], [3,1] ]
725        z = [2, 4, 4]
726
727        ref = fit_to_mesh(points, triangles, data_coords, z)
728
729        #print attributes
730        #print ref
731        assert allclose(attributes, ref)
732
733        os.remove(FN)
734
735
736    #Tests of smoothing matrix
737    def test_smoothing_matrix_one_triangle(self):
738        from Numeric import dot
739        a = [0.0, 0.0]
740        b = [0.0, 2.0]
741        c = [2.0,0.0]
742        points = [a, b, c]
743
744        vertices = [ [1,0,2] ]   #bac
745
746        interp = Interpolation(points, vertices)
747
748        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
749                                   [-0.5, 0.5, 0],
750                                   [-0.5, 0, 0.5]])
751
752        #Define f(x,y) = x
753        f = array([0,0,2]) #Value at global vertex 2
754
755        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
756        #           int 1 dx dy = area = 2
757        assert dot(dot(f, interp.get_D()), f) == 2
758
759        #Define f(x,y) = y
760        f = array([0,2,0])  #Value at global vertex 1
761
762        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
763        #           int 1 dx dy = area = 2
764        assert dot(dot(f, interp.get_D()), f) == 2
765
766        #Define f(x,y) = x+y
767        f = array([0,2,2])  #Values at global vertex 1 and 2
768
769        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
770        #           int 2 dx dy = 2*area = 4
771        assert dot(dot(f, interp.get_D()), f) == 4
772
773
774
775    def test_smoothing_matrix_more_triangles(self):
776        from Numeric import dot
777
778        a = [0.0, 0.0]
779        b = [0.0, 2.0]
780        c = [2.0,0.0]
781        d = [0.0, 4.0]
782        e = [2.0, 2.0]
783        f = [4.0,0.0]
784
785        points = [a, b, c, d, e, f]
786        #bac, bce, ecf, dbe, daf, dae
787        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
788
789        interp = Interpolation(points, vertices)
790
791
792        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
793        #                           [-0.5, 0.5, 0],
794        #                           [-0.5, 0, 0.5]])
795
796        #Define f(x,y) = x
797        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
798
799        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
800        #           int 1 dx dy = total area = 8
801        assert dot(dot(f, interp.get_D()), f) == 8
802
803        #Define f(x,y) = y
804        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
805
806        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
807        #           int 1 dx dy = area = 8
808        assert dot(dot(f, interp.get_D()), f) == 8
809
810        #Define f(x,y) = x+y
811        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
812
813        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
814        #           int 2 dx dy = 2*area = 16
815        assert dot(dot(f, interp.get_D()), f) == 16
816
817
818    def test_fit_and_interpolation(self):
819        from mesh import Mesh
820
821        a = [0.0, 0.0]
822        b = [0.0, 2.0]
823        c = [2.0, 0.0]
824        d = [0.0, 4.0]
825        e = [2.0, 2.0]
826        f = [4.0, 0.0]
827
828        points = [a, b, c, d, e, f]
829        #bac, bce, ecf, dbe, daf, dae
830        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
831
832        #Get (enough) datapoints
833        data_points = [[ 0.66666667, 0.66666667],
834                       [ 1.33333333, 1.33333333],
835                       [ 2.66666667, 0.66666667],
836                       [ 0.66666667, 2.66666667],
837                       [ 0.0, 1.0],
838                       [ 0.0, 3.0],
839                       [ 1.0, 0.0],
840                       [ 1.0, 1.0],
841                       [ 1.0, 2.0],
842                       [ 1.0, 3.0],
843                       [ 2.0, 1.0],
844                       [ 3.0, 0.0],
845                       [ 3.0, 1.0]]
846
847        interp = Interpolation(points, triangles, data_points, alpha=0.0)
848
849        z = linear_function(data_points)
850        answer = linear_function(points)
851
852        f = interp.fit(z)
853
854        #print "f",f
855        #print "answer",answer
856        assert allclose(f, answer)
857
858        #Map back
859        z1 = interp.interpolate(f)
860        #print "z1\n", z1
861        #print "z\n",z
862        assert allclose(z, z1)
863
864
865    def test_smoothing_and_interpolation(self):
866
867        a = [0.0, 0.0]
868        b = [0.0, 2.0]
869        c = [2.0, 0.0]
870        d = [0.0, 4.0]
871        e = [2.0, 2.0]
872        f = [4.0, 0.0]
873
874        points = [a, b, c, d, e, f]
875        #bac, bce, ecf, dbe, daf, dae
876        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
877
878        #Get (too few!) datapoints
879        data_points = [[ 0.66666667, 0.66666667],
880                       [ 1.33333333, 1.33333333],
881                       [ 2.66666667, 0.66666667],
882                       [ 0.66666667, 2.66666667]]
883
884        z = linear_function(data_points)
885        answer = linear_function(points)
886
887        #Make interpolator with too few data points and no smoothing
888        interp = Interpolation(points, triangles, data_points, alpha=0.0)
889        #Must raise an exception
890        try:
891            f = interp.fit(z)
892        except:
893            pass
894
895        #Now try with smoothing parameter
896        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
897
898        f = interp.fit(z)
899        #f will be different from answer due to smoothing
900        assert allclose(f, answer,atol=5)
901
902        #Map back
903        z1 = interp.interpolate(f)
904        assert allclose(z, z1)
905
906
907
908    def test_fit_and_interpolation_with_new_points(self):
909        """Fit a surface to one set of points. Then interpolate that surface
910        using another set of points.
911        """
912        from mesh import Mesh
913
914
915        #Setup mesh used to represent fitted function
916        a = [0.0, 0.0]
917        b = [0.0, 2.0]
918        c = [2.0, 0.0]
919        d = [0.0, 4.0]
920        e = [2.0, 2.0]
921        f = [4.0, 0.0]
922
923        points = [a, b, c, d, e, f]
924        #bac, bce, ecf, dbe, daf, dae
925        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
926
927        #Datapoints to fit from
928        data_points1 = [[ 0.66666667, 0.66666667],
929                        [ 1.33333333, 1.33333333],
930                        [ 2.66666667, 0.66666667],
931                        [ 0.66666667, 2.66666667],
932                        [ 0.0, 1.0],
933                        [ 0.0, 3.0],
934                        [ 1.0, 0.0],
935                        [ 1.0, 1.0],
936                        [ 15, -17],   #Outside mesh
937                        [ 1.0, 2.0],
938                        [ 1.0, 3.0],
939                        [ 2.0, 1.0],
940                        [ 3.0, 0.0],
941                        [ 3.0, 1.0]]
942
943        #Fit surface to mesh
944        interp = Interpolation(points, triangles, data_points1, alpha=0.0,
945                               precrop = True)
946        z = linear_function(data_points1) #Example z-values
947        f = interp.fit(z)                 #Fitted values at vertices
948
949
950
951        #New datapoints where interpolated values are sought
952        data_points2 = [[ 0.0, 0.0],
953                        [ 0.5, 0.5],
954                        [ 0.7, 0.7],
955                        [-13, 65],  #Outside
956                        [ 1.0, 0.5],
957                        [ 2.0, 0.4],
958                        [ 2.8, 1.2]]
959
960
961
962        #Build new A matrix based on new points (without precrop)
963        interp.build_interpolation_matrix_A(data_points2, precrop = False)
964
965        #Interpolate using fitted surface
966        z1 = interp.interpolate(f)
967
968        #import Numeric
969        #data_points2 = Numeric.take(data_points2, interp.point_indices)
970
971        #Desired result (OK for points inside)
972
973        answer = linear_function(data_points2)
974        import Numeric
975        z1 = Numeric.take(z1, [0,1,2,4,5,6])
976        answer = Numeric.take(answer, [0,1,2,4,5,6])
977        assert allclose(z1, answer)
978
979        #Build new A matrix based on new points (with precrop)
980        interp.build_interpolation_matrix_A(data_points2, precrop = True)
981
982        #Interpolate using fitted surface
983        z1 = interp.interpolate(f)
984
985        import Numeric
986        data_points2 = Numeric.take(data_points2, interp.point_indices)
987
988        #Desired result
989        answer = linear_function(data_points2)
990        assert allclose(z1, answer)
991
992
993
994
995
996
997    def test_interpolation_from_discontinuous_vertex_values(self):
998        """test_interpolation_from_discontinuous_vertex_values.
999        This will test the format used internally in pyvolution and also
1000        interpolation from sww files
1001        """
1002       
1003        from mesh import Mesh
1004
1005
1006        #Setup mesh used to represent discontinuous function
1007        a = [0.0, 0.0]
1008        b = [0.0, 2.0]
1009        c = [2.0, 0.0]
1010        d = [0.0, 4.0]
1011        e = [2.0, 2.0]
1012        f = [4.0, 0.0]
1013
1014        points = [b, a, c,
1015                  b, c, e,
1016                  e, c, f,
1017                  d, b, e]
1018       
1019        #bac, bce, ecf, dbe
1020        triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
1021
1022       
1023        vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
1024                         
1025       
1026
1027        #New datapoints where interpolated values are sought
1028        data_points = [[0.0, 0.0],  #T0
1029                       [0.5, 0.5],  #T0
1030                       [1.5, 1.5],  #T1
1031                       [2.5, 0.5],  #T2
1032                       [0.0, 3.0],  #T3
1033                       [1.0, 2.0],  #In between T1 and T3 (T1 is used) FIXME?
1034                       [2.0, 1.0],  #In between T1 and T2 (T1 is used) FIXME?
1035                       [1.0, 1.0]]  #In between T1 and T0 (T0 is used) FIXME?
1036       
1037
1038
1039
1040        #Build interpolation matrix
1041        interp = Interpolation(points, triangles, data_points)
1042                               #, alpha=0.0, precrop = True)
1043
1044        #print interp.A.todense()
1045        #print vertex_values
1046
1047        #Interpolate using fitted surface
1048        z = interp.interpolate(vertex_values)
1049
1050        #print z
1051
1052        assert allclose(z, [0,0,1,2,5,1,1,0])
1053
1054
1055
1056
1057    def test_interpolation_function_time_only(self):
1058        """Test spatio-temporal interpolation
1059        Test that spatio temporal function performs the correct
1060        interpolations in both time and space
1061        """
1062
1063
1064        #Three timesteps
1065        time = [1.0, 5.0, 6.0]
1066       
1067
1068        #One quantity
1069        Q = zeros( (3,6), Float )
1070
1071        #Linear in time and space
1072        a = [0.0, 0.0]
1073        b = [0.0, 2.0]
1074        c = [2.0, 0.0]
1075        d = [0.0, 4.0]
1076        e = [2.0, 2.0]
1077        f = [4.0, 0.0]
1078
1079        points = [a, b, c, d, e, f]
1080       
1081        for i, t in enumerate(time):
1082            Q[i, :] = t*linear_function(points)
1083
1084           
1085        #Check basic interpolation of one quantity using averaging
1086        #(no interpolation points or spatial info)
1087        from util import mean       
1088        I = Interpolation_function(time, [mean(Q[0,:]),
1089                                          mean(Q[1,:]),
1090                                          mean(Q[2,:])])
1091
1092
1093
1094        #Check temporal interpolation
1095        for i in [0,1,2]:
1096            assert allclose(I(time[i]), mean(Q[i,:]))
1097
1098        #Midway   
1099        assert allclose(I( (time[0] + time[1])/2 ),
1100                        (I(time[0]) + I(time[1]))/2 )
1101
1102        assert allclose(I( (time[1] + time[2])/2 ),
1103                        (I(time[1]) + I(time[2]))/2 )
1104
1105        assert allclose(I( (time[0] + time[2])/2 ),
1106                        (I(time[0]) + I(time[2]))/2 )                 
1107
1108        #1/3
1109        assert allclose(I( (time[0] + time[2])/3 ),
1110                        (I(time[0]) + I(time[2]))/3 )                         
1111
1112
1113        #Out of bounds checks
1114        try:
1115            I(time[0]-1) 
1116        except:
1117            pass
1118        else:
1119            raise 'Should raise exception'
1120
1121        try:
1122            I(time[-1]+1) 
1123        except:
1124            pass
1125        else:
1126            raise 'Should raise exception'       
1127
1128
1129       
1130
1131    def test_interpolation_function_spatial_only(self):
1132        """Test spatio-temporal interpolation with constant time
1133        """
1134
1135        #Three timesteps
1136        time = [1.0, 5.0, 6.0]
1137       
1138       
1139        #Setup mesh used to represent fitted function
1140        a = [0.0, 0.0]
1141        b = [0.0, 2.0]
1142        c = [2.0, 0.0]
1143        d = [0.0, 4.0]
1144        e = [2.0, 2.0]
1145        f = [4.0, 0.0]
1146
1147        points = [a, b, c, d, e, f]
1148        #bac, bce, ecf, dbe
1149        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1150
1151
1152        #New datapoints where interpolated values are sought
1153        interpolation_points = [[ 0.0, 0.0],
1154                                [ 0.5, 0.5],
1155                                [ 0.7, 0.7],
1156                                [ 1.0, 0.5],
1157                                [ 2.0, 0.4],
1158                                [ 2.8, 1.2]]
1159
1160
1161        #One quantity linear in space
1162        Q = linear_function(points)
1163
1164
1165        #Check interpolation of one quantity using interpolaton points
1166        I = Interpolation_function(time, Q,
1167                                   vertex_coordinates = points,
1168                                   triangles = triangles, 
1169                                   interpolation_points = interpolation_points,
1170                                   verbose = False)
1171
1172
1173        answer = linear_function(interpolation_points)
1174
1175        t = time[0]
1176        for j in range(50): #t in [1, 6]
1177            for id in range(len(interpolation_points)):
1178                assert allclose(I(t, id), answer[id])
1179
1180            t += 0.1   
1181
1182
1183        try:   
1184            I(1)
1185        except:
1186            pass
1187        else:
1188            raise 'Should raise exception'
1189           
1190
1191
1192    def test_interpolation_function(self):
1193        """Test spatio-temporal interpolation
1194        Test that spatio temporal function performs the correct
1195        interpolations in both time and space
1196        """
1197
1198
1199        #Three timesteps
1200        time = [1.0, 5.0, 6.0]
1201       
1202
1203        #Setup mesh used to represent fitted function
1204        a = [0.0, 0.0]
1205        b = [0.0, 2.0]
1206        c = [2.0, 0.0]
1207        d = [0.0, 4.0]
1208        e = [2.0, 2.0]
1209        f = [4.0, 0.0]
1210
1211        points = [a, b, c, d, e, f]
1212        #bac, bce, ecf, dbe
1213        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1214
1215
1216        #New datapoints where interpolated values are sought
1217        interpolation_points = [[ 0.0, 0.0],
1218                                [ 0.5, 0.5],
1219                                [ 0.7, 0.7],
1220                                [ 1.0, 0.5],
1221                                [ 2.0, 0.4],
1222                                [ 2.8, 1.2]]
1223
1224
1225        #One quantity
1226        Q = zeros( (3,6), Float )
1227
1228        #Linear in time and space
1229        for i, t in enumerate(time):
1230            Q[i, :] = t*linear_function(points)
1231
1232
1233        #Check interpolation of one quantity using interpolaton points)
1234        I = Interpolation_function(time, Q,
1235                                   vertex_coordinates = points,
1236                                   triangles = triangles, 
1237                                   interpolation_points = interpolation_points,
1238                                   verbose = False)
1239
1240
1241        answer = linear_function(interpolation_points)
1242
1243        t = time[0]
1244        for j in range(50): #t in [1, 6]
1245            for id in range(len(interpolation_points)):
1246                assert allclose(I(t, id), t*answer[id])
1247
1248            t += 0.1   
1249
1250
1251        try:   
1252            I(1)
1253        except:
1254            pass
1255        else:
1256            raise 'Should raise exception'
1257           
1258        #
1259        #interpolation_points = [[ 0.0, 0.0],
1260        #                        [ 0.5, 0.5],
1261        #                        [ 0.7, 0.7],
1262        #                        [-13, 65],  #Outside
1263        #                        [ 1.0, 0.5],
1264        #                        [ 2.0, 0.4],
1265        #                        [ 2.8, 1.2]]
1266        #
1267        #try:
1268        #    I = Interpolation_function(time, Q,
1269        #                               vertex_coordinates = points,
1270        #                               triangles = triangles,
1271        #                               interpolation_points = interpolation_points,
1272        #                               verbose = False)
1273        #except:
1274        #    pass
1275        #else:
1276        #    raise 'Should raise exception'
1277
1278
1279
1280
1281
1282    def test_fit_and_interpolation_with_different_origins(self):
1283        """Fit a surface to one set of points. Then interpolate that surface
1284        using another set of points.
1285        This test tests situtaion where points and mesh belong to a different
1286        coordinate system as defined by origin.
1287        """
1288        from mesh import Mesh
1289
1290        #Setup mesh used to represent fitted function
1291        a = [0.0, 0.0]
1292        b = [0.0, 2.0]
1293        c = [2.0, 0.0]
1294        d = [0.0, 4.0]
1295        e = [2.0, 2.0]
1296        f = [4.0, 0.0]
1297
1298        points = [a, b, c, d, e, f]
1299        #bac, bce, ecf, dbe, daf, dae
1300        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1301
1302        #Datapoints to fit from
1303        data_points1 = [[ 0.66666667, 0.66666667],
1304                        [ 1.33333333, 1.33333333],
1305                        [ 2.66666667, 0.66666667],
1306                        [ 0.66666667, 2.66666667],
1307                        [ 0.0, 1.0],
1308                        [ 0.0, 3.0],
1309                        [ 1.0, 0.0],
1310                        [ 1.0, 1.0],
1311                        [ 1.0, 2.0],
1312                        [ 1.0, 3.0],
1313                        [ 2.0, 1.0],
1314                        [ 3.0, 0.0],
1315                        [ 3.0, 1.0]]
1316
1317
1318        #First check that things are OK when using same origin
1319        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1320        data_origin = (56, 290000, 618000) #zone, easting, northing
1321
1322
1323        #Fit surface to mesh
1324        interp = Interpolation(points, triangles, data_points1,
1325                               alpha=0.0,
1326                               data_origin = data_origin,
1327                               mesh_origin = mesh_origin)
1328
1329        z = linear_function(data_points1) #Example z-values
1330        f = interp.fit(z)                 #Fitted values at vertices
1331
1332
1333        #New datapoints where interpolated values are sought
1334        data_points2 = [[ 0.0, 0.0],
1335                        [ 0.5, 0.5],
1336                        [ 0.7, 0.7],
1337                        [ 1.0, 0.5],
1338                        [ 2.0, 0.4],
1339                        [ 2.8, 1.2]]
1340
1341
1342        #Build new A matrix based on new points
1343        interp.build_interpolation_matrix_A(data_points2)
1344
1345        #Interpolate using fitted surface
1346        z1 = interp.interpolate(f)
1347
1348        #Desired result
1349        answer = linear_function(data_points2)
1350        assert allclose(z1, answer)
1351
1352
1353        ##############################################
1354
1355        #Then check situation where points are relative to a different
1356        #origin (same zone, though, until we figure that out (FIXME))
1357
1358        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1359        data_origin = (56, 10000, 10000) #zone, easting, northing
1360
1361        #Shift datapoints according to new origin
1362
1363        for k in range(len(data_points1)):
1364            data_points1[k][0] += mesh_origin[1] - data_origin[1]
1365            data_points1[k][1] += mesh_origin[2] - data_origin[2]
1366
1367        for k in range(len(data_points2)):
1368            data_points2[k][0] += mesh_origin[1] - data_origin[1]
1369            data_points2[k][1] += mesh_origin[2] - data_origin[2]
1370
1371
1372
1373        #Fit surface to mesh
1374        interp = Interpolation(points, triangles, data_points1,
1375                               alpha=0.0,
1376                               data_origin = data_origin,
1377                               mesh_origin = mesh_origin)
1378
1379        f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
1380
1381        assert allclose(f,f1), 'Fit should have been unaltered'
1382
1383
1384        #Build new A matrix based on new points
1385        interp.build_interpolation_matrix_A(data_points2)
1386
1387        #Interpolate using fitted surface
1388        z1 = interp.interpolate(f)
1389        assert allclose(z1, answer)
1390
1391
1392        #########################################################
1393        #Finally try to relate data_points2 to new origin without
1394        #rebuilding matrix
1395
1396        data_origin = (56, 2000, 2000) #zone, easting, northing
1397        for k in range(len(data_points2)):
1398            data_points2[k][0] += 8000
1399            data_points2[k][1] += 8000
1400
1401        #Build new A matrix based on new points
1402        interp.build_interpolation_matrix_A(data_points2,
1403                                            data_origin = data_origin)
1404
1405        #Interpolate using fitted surface
1406        z1 = interp.interpolate(f)
1407        assert allclose(z1, answer)
1408
1409
1410
1411    def test_fit_to_mesh_w_georef(self):
1412        """Simple check that georef works at the fit_to_mesh level
1413        """
1414       
1415        from coordinate_transforms.geo_reference import Geo_reference
1416
1417        #Mesh
1418        vertex_coordinates = [[0.76, 0.76],
1419                              [0.76, 5.76],
1420                              [5.76, 0.76]]
1421        triangles = [[0,2,1]]
1422
1423        mesh_geo = Geo_reference(56,-0.76,-0.76)
1424
1425
1426        #Data                       
1427        data_points = [[ 201.0, 401.0],
1428                       [ 201.0, 403.0],
1429                       [ 203.0, 401.0]]
1430
1431        z = [2, 4, 4]
1432       
1433        data_geo = Geo_reference(56,-200,-400)
1434       
1435        #Fit
1436        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
1437                         data_origin = data_geo.get_origin(),
1438                         mesh_origin = mesh_geo.get_origin(),
1439                         alpha = 0)
1440
1441        assert allclose( zz, [0,5,5] )
1442
1443
1444    def test_fit_to_mesh_file(self):
1445        from load_mesh.loadASCII import import_mesh_file, \
1446             export_mesh_file
1447        import tempfile
1448        import os
1449
1450        # create a .tsh file, no user outline
1451        mesh_dic = {}
1452        mesh_dic['vertices'] = [[0.0, 0.0],
1453                                          [0.0, 5.0],
1454                                          [5.0, 0.0]]
1455        mesh_dic['triangles'] =  [[0, 2, 1]]
1456        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1457        mesh_dic['triangle_tags'] = ['']
1458        mesh_dic['vertex_attributes'] = [[], [], []]
1459        mesh_dic['vertiex_attribute_titles'] = []
1460        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1461        mesh_dic['segment_tags'] = ['external',
1462                                                  'external',
1463                                                  'external']
1464        mesh_file = tempfile.mktemp(".tsh")
1465        export_mesh_file(mesh_file,mesh_dic)
1466
1467        # create an .xya file
1468        point_file = tempfile.mktemp(".xya")
1469        fd = open(point_file,'w')
1470        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")
1471        fd.close()
1472
1473        mesh_output_file = tempfile.mktemp(".tsh") 
1474        fit_to_mesh_file(mesh_file,
1475                         point_file,
1476                         mesh_output_file,
1477                         alpha = 0.0)
1478        # load in the .tsh file we just wrote
1479        mesh_dic = import_mesh_file(mesh_output_file)
1480        #print "mesh_dic",mesh_dic
1481        ans =[[0.0, 0.0],
1482              [5.0, 10.0],
1483              [5.0,10.0]]
1484        assert allclose(mesh_dic['vertex_attributes'],ans)
1485
1486        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1487                        ['elevation','stage'],
1488                        'test_fit_to_mesh_file failed')
1489
1490        #clean up
1491        os.remove(mesh_file)
1492        os.remove(point_file)
1493        os.remove(mesh_output_file)
1494
1495    def test_fit_to_mesh_file3(self):
1496        from load_mesh.loadASCII import import_mesh_file, \
1497             export_mesh_file
1498        import tempfile
1499        import os
1500
1501        # create a .tsh file, no user outline
1502        mesh_dic = {}
1503        mesh_dic['vertices'] = [[0.76, 0.76],
1504                                          [0.76, 5.76],
1505                                          [5.76, 0.76]]
1506        mesh_dic['triangles'] =  [[0, 2, 1]]
1507        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1508        mesh_dic['triangle_tags'] = ['']
1509        mesh_dic['vertex_attributes'] = [[], [], []]
1510        mesh_dic['vertiex_attribute_titles'] = []
1511        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1512        mesh_dic['segment_tags'] = ['external',
1513                                                  'external',
1514                                                  'external']
1515        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1516        mesh_file = tempfile.mktemp(".tsh")
1517        export_mesh_file(mesh_file,mesh_dic)
1518
1519        # create an .xya file
1520        point_file = tempfile.mktemp(".xya")
1521        fd = open(point_file,'w')
1522        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")
1523        fd.close()
1524
1525        mesh_output_file = tempfile.mktemp(".tsh")
1526        fit_to_mesh_file(mesh_file,
1527                         point_file,
1528                         mesh_output_file,
1529                         alpha = 0.0)
1530        # load in the .tsh file we just wrote
1531        mesh_dic = import_mesh_file(mesh_output_file)
1532        #print "mesh_dic",mesh_dic
1533        ans =[[0.0, 0.0],
1534              [5.0, 10.0],
1535              [5.0,10.0]]
1536        assert allclose(mesh_dic['vertex_attributes'],ans)
1537
1538        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1539                        ['elevation','stage'],
1540                        'test_fit_to_mesh_file failed')
1541
1542        #clean up
1543        os.remove(mesh_file)
1544        os.remove(point_file)
1545        os.remove(mesh_output_file)
1546
1547    def test_fit_to_mesh_file4(self):
1548        from load_mesh.loadASCII import import_mesh_file, \
1549             export_mesh_file
1550        import tempfile
1551        import os
1552
1553        # create a .tsh file, no user outline
1554        mesh_dic = {}
1555        mesh_dic['vertices'] = [[0.76, 0.76],
1556                                [0.76, 5.76],
1557                                [5.76, 0.76]]
1558        mesh_dic['triangles'] =  [[0, 2, 1]]
1559        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1560        mesh_dic['triangle_tags'] = ['']
1561        mesh_dic['vertex_attributes'] = [[], [], []]
1562        mesh_dic['vertiex_attribute_titles'] = []
1563        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1564        mesh_dic['segment_tags'] = ['external',
1565                                    'external',
1566                                    'external']
1567        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1568        mesh_file = tempfile.mktemp(".tsh")
1569        export_mesh_file(mesh_file,mesh_dic)
1570
1571        geo_ref = Geo_reference(56,-200,-400)
1572        # create an .xya file
1573        point_file = tempfile.mktemp(".xya")
1574        fd = open(point_file,'w')
1575        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")
1576        geo_ref.write_ASCII(fd)
1577        fd.close()
1578
1579        mesh_output_file = tempfile.mktemp(".tsh")
1580        fit_to_mesh_file(mesh_file,
1581                         point_file,
1582                         mesh_output_file,
1583                         alpha = 0.0)
1584        # load in the .tsh file we just wrote
1585        mesh_dic = import_mesh_file(mesh_output_file)
1586        #print "mesh_dic",mesh_dic
1587        ans =[[0.0, 0.0],
1588              [5.0, 10.0],
1589              [5.0, 10.0]]
1590        assert allclose(mesh_dic['vertex_attributes'],ans)
1591
1592        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1593                        ['elevation','stage'],
1594                        'test_fit_to_mesh_file failed')
1595
1596        #clean up
1597        os.remove(mesh_file)
1598        os.remove(point_file)
1599        os.remove(mesh_output_file)
1600
1601    def test_fit_to_mesh_fileII(self):
1602        from load_mesh.loadASCII import import_mesh_file, \
1603             export_mesh_file
1604        import tempfile
1605        import os
1606
1607        # create a .tsh file, no user outline
1608        mesh_dic = {}
1609        mesh_dic['vertices'] = [[0.0, 0.0],
1610                                [0.0, 5.0],
1611                                [5.0, 0.0]]
1612        mesh_dic['triangles'] =  [[0, 2, 1]]
1613        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1614        mesh_dic['triangle_tags'] = ['']
1615        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1616        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1617        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1618        mesh_dic['segment_tags'] = ['external',
1619                                                  'external',
1620                                                  'external']
1621        mesh_file = tempfile.mktemp(".tsh")
1622        export_mesh_file(mesh_file,mesh_dic)
1623
1624        # create an .xya file
1625        point_file = tempfile.mktemp(".xya")
1626        fd = open(point_file,'w')
1627        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")
1628        fd.close()
1629
1630        mesh_output_file = "new_triangle.tsh"
1631        fit_to_mesh_file(mesh_file,
1632                         point_file,
1633                         mesh_output_file,
1634                         alpha = 0.0)
1635        # load in the .tsh file we just wrote
1636        mesh_dic = import_mesh_file(mesh_output_file)
1637
1638        assert allclose(mesh_dic['vertex_attributes'],
1639                        [[1.0, 2.0,0.0, 0.0],
1640                         [1.0, 2.0,5.0, 10.0],
1641                         [1.0, 2.0,5.0,10.0]])
1642
1643        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1644                        ['density', 'temp','elevation','stage'],
1645                        'test_fit_to_mesh_file failed')
1646
1647        #clean up
1648        os.remove(mesh_file)
1649        os.remove(mesh_output_file)
1650        os.remove(point_file)
1651
1652    def test_fit_to_mesh_file_errors(self):
1653        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1654        import tempfile
1655        import os
1656
1657        # create a .tsh file, no user outline
1658        mesh_dic = {}
1659        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1660        mesh_dic['triangles'] =  [[0, 2, 1]]
1661        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1662        mesh_dic['triangle_tags'] = ['']
1663        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1664        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1665        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1666        mesh_dic['segment_tags'] = ['external', 'external','external']
1667        mesh_file = tempfile.mktemp(".tsh")
1668        export_mesh_file(mesh_file,mesh_dic)
1669
1670        # create an .xya file
1671        point_file = tempfile.mktemp(".xya")
1672        fd = open(point_file,'w')
1673        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")
1674        fd.close()
1675
1676        mesh_output_file = "new_triangle.tsh"
1677        try:
1678            fit_to_mesh_file(mesh_file, point_file,
1679                             mesh_output_file, display_errors = False)
1680        except IOError:
1681            pass
1682        else:
1683            #self.failUnless(0 ==1,  'Bad file did not raise error!')
1684            raise 'Bad file did not raise error!'
1685           
1686        #clean up
1687        os.remove(mesh_file)
1688        os.remove(point_file)
1689
1690    def test_fit_to_mesh_file_errorsII(self):
1691        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1692        import tempfile
1693        import os
1694
1695        # create a .tsh file, no user outline
1696        mesh_file = tempfile.mktemp(".tsh")
1697        fd = open(mesh_file,'w')
1698        fd.write("unit testing a bad .tsh file \n")
1699        fd.close()
1700
1701        # create an .xya file
1702        point_file = tempfile.mktemp(".xya")
1703        fd = open(point_file,'w')
1704        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")
1705        fd.close()
1706
1707        mesh_output_file = "new_triangle.tsh"
1708        try:
1709            fit_to_mesh_file(mesh_file, point_file,
1710                             mesh_output_file, display_errors = False)
1711        except IOError:
1712            pass
1713        else:
1714            raise 'Bad file did not raise error!'
1715           
1716        #clean up
1717        os.remove(mesh_file)
1718        os.remove(point_file)
1719
1720    def test_fit_to_mesh_file_errorsIII(self):
1721        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1722        import tempfile
1723        import os
1724
1725        # create a .tsh file, no user outline
1726        mesh_dic = {}
1727        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1728        mesh_dic['triangles'] =  [[0, 2, 1]]
1729        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1730        mesh_dic['triangle_tags'] = ['']
1731        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1732        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1733        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1734        mesh_dic['segment_tags'] = ['external', 'external','external']
1735        mesh_file = tempfile.mktemp(".tsh")
1736        export_mesh_file(mesh_file,mesh_dic)
1737
1738        # create an .xya file
1739        point_file = tempfile.mktemp(".xya")
1740        fd = open(point_file,'w')
1741        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")
1742        fd.close()
1743
1744        #This a deliberately illegal filename to invoke the error.
1745        mesh_output_file = ".../\z\z:ya.tsh"       
1746
1747        try:
1748            fit_to_mesh_file(mesh_file, point_file,
1749                             mesh_output_file, display_errors = False)
1750        except IOError:
1751            pass
1752        else:
1753            raise 'Bad file did not raise error!'
1754       
1755        #clean up
1756        os.remove(mesh_file)
1757        os.remove(point_file)
1758 
1759    def test_fit_to_mesh_file_errorsIV(self):
1760        import os
1761        command = '%s least_squares.py q q q e n 0.9 n'  %(sys.executable)
1762        status = os.system(command) 
1763        self.failUnless(status%255  == 1,
1764                        'command prompt least_squares.py failed.  Incorect exit status.')
1765       
1766    def test_fit_to_msh_netcdf_fileII(self):
1767        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1768        import tempfile
1769        import os
1770
1771        # create a .tsh file, no user outline
1772        mesh_dic = {}
1773        mesh_dic['vertices'] = [[0.0, 0.0],
1774                                [0.0, 5.0],
1775                                [5.0, 0.0]]
1776        mesh_dic['triangles'] =  [[0, 2, 1]]
1777        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1778        mesh_dic['triangle_tags'] = ['']
1779        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1780        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1781        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1782        mesh_dic['segment_tags'] = ['external',
1783                                                  'external',
1784                                                  'external']
1785        mesh_file = tempfile.mktemp(".msh")
1786        export_mesh_file(mesh_file,mesh_dic)
1787
1788        # create an .xya file
1789        point_file = tempfile.mktemp(".xya")
1790        fd = open(point_file,'w')
1791        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")
1792        fd.close()
1793
1794        mesh_output_file = "new_triangle.msh"
1795        fit_to_mesh_file(mesh_file,
1796                         point_file,
1797                         mesh_output_file,
1798                         alpha = 0.0)
1799        # load in the .tsh file we just wrote
1800        mesh_dic = import_mesh_file(mesh_output_file)
1801
1802        assert allclose(mesh_dic['vertex_attributes'],
1803                        [[1.0, 2.0,0.0, 0.0],
1804                         [1.0, 2.0,5.0, 10.0],
1805                         [1.0, 2.0,5.0,10.0]])
1806
1807        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1808                        ['density', 'temp','elevation','stage'],
1809                        'test_fit_to_mesh_file failed')
1810
1811        #clean up
1812        os.remove(mesh_file)
1813        os.remove(mesh_output_file)
1814        os.remove(point_file)
1815
1816#-------------------------------------------------------------
1817if __name__ == "__main__":
1818    suite = unittest.makeSuite(Test_Least_Squares,'test')
1819
1820    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
1821    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
1822    runner = unittest.TextTestRunner(verbosity=1)
1823    runner.run(suite)
1824
1825
1826
1827
1828
Note: See TracBrowser for help on using the repository browser.