source: inundation/pyvolution/test_least_squares.py @ 2260

Last change on this file since 2260 was 2260, checked in by duncan, 18 years ago

add test with different geo-refs in mesh and points files

File size: 56.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    def test_fit_to_mesh_file(self):
1411        from load_mesh.loadASCII import import_mesh_file, \
1412             export_mesh_file
1413        import tempfile
1414        import os
1415
1416        # create a .tsh file, no user outline
1417        mesh_dic = {}
1418        mesh_dic['vertices'] = [[0.0, 0.0],
1419                                          [0.0, 5.0],
1420                                          [5.0, 0.0]]
1421        mesh_dic['triangles'] =  [[0, 2, 1]]
1422        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1423        mesh_dic['triangle_tags'] = ['']
1424        mesh_dic['vertex_attributes'] = [[], [], []]
1425        mesh_dic['vertiex_attribute_titles'] = []
1426        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1427        mesh_dic['segment_tags'] = ['external',
1428                                                  'external',
1429                                                  'external']
1430        mesh_file = tempfile.mktemp(".tsh")
1431        export_mesh_file(mesh_file,mesh_dic)
1432
1433        # create an .xya file
1434        point_file = tempfile.mktemp(".xya")
1435        fd = open(point_file,'w')
1436        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")
1437        fd.close()
1438
1439        mesh_output_file = tempfile.mktemp(".tsh") 
1440        fit_to_mesh_file(mesh_file,
1441                         point_file,
1442                         mesh_output_file,
1443                         alpha = 0.0)
1444        # load in the .tsh file we just wrote
1445        mesh_dic = import_mesh_file(mesh_output_file)
1446        #print "mesh_dic",mesh_dic
1447        ans =[[0.0, 0.0],
1448              [5.0, 10.0],
1449              [5.0,10.0]]
1450        assert allclose(mesh_dic['vertex_attributes'],ans)
1451
1452        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1453                        ['elevation','stage'],
1454                        'test_fit_to_mesh_file failed')
1455
1456        #clean up
1457        os.remove(mesh_file)
1458        os.remove(point_file)
1459        os.remove(mesh_output_file)
1460
1461    def test_fit_to_mesh_file3(self):
1462        from load_mesh.loadASCII import import_mesh_file, \
1463             export_mesh_file
1464        import tempfile
1465        import os
1466
1467        # create a .tsh file, no user outline
1468        mesh_dic = {}
1469        mesh_dic['vertices'] = [[0.76, 0.76],
1470                                          [0.76, 5.76],
1471                                          [5.76, 0.76]]
1472        mesh_dic['triangles'] =  [[0, 2, 1]]
1473        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1474        mesh_dic['triangle_tags'] = ['']
1475        mesh_dic['vertex_attributes'] = [[], [], []]
1476        mesh_dic['vertiex_attribute_titles'] = []
1477        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1478        mesh_dic['segment_tags'] = ['external',
1479                                                  'external',
1480                                                  'external']
1481        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1482        mesh_file = tempfile.mktemp(".tsh")
1483        export_mesh_file(mesh_file,mesh_dic)
1484
1485        # create an .xya file
1486        point_file = tempfile.mktemp(".xya")
1487        fd = open(point_file,'w')
1488        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")
1489        fd.close()
1490
1491        mesh_output_file = tempfile.mktemp(".tsh")
1492        fit_to_mesh_file(mesh_file,
1493                         point_file,
1494                         mesh_output_file,
1495                         alpha = 0.0)
1496        # load in the .tsh file we just wrote
1497        mesh_dic = import_mesh_file(mesh_output_file)
1498        #print "mesh_dic",mesh_dic
1499        ans =[[0.0, 0.0],
1500              [5.0, 10.0],
1501              [5.0,10.0]]
1502        assert allclose(mesh_dic['vertex_attributes'],ans)
1503
1504        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1505                        ['elevation','stage'],
1506                        'test_fit_to_mesh_file failed')
1507
1508        #clean up
1509        os.remove(mesh_file)
1510        os.remove(point_file)
1511        os.remove(mesh_output_file)
1512
1513    def test_fit_to_mesh_file4(self):
1514        from load_mesh.loadASCII import import_mesh_file, \
1515             export_mesh_file
1516        import tempfile
1517        import os
1518
1519        # create a .tsh file, no user outline
1520        mesh_dic = {}
1521        mesh_dic['vertices'] = [[0.76, 0.76],
1522                                          [0.76, 5.76],
1523                                          [5.76, 0.76]]
1524        mesh_dic['triangles'] =  [[0, 2, 1]]
1525        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1526        mesh_dic['triangle_tags'] = ['']
1527        mesh_dic['vertex_attributes'] = [[], [], []]
1528        mesh_dic['vertiex_attribute_titles'] = []
1529        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1530        mesh_dic['segment_tags'] = ['external',
1531                                                  'external',
1532                                                  'external']
1533        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1534        mesh_file = tempfile.mktemp(".tsh")
1535        export_mesh_file(mesh_file,mesh_dic)
1536
1537        geo_ref = Geo_reference(56,-200,-400)
1538        # create an .xya file
1539        point_file = tempfile.mktemp(".xya")
1540        fd = open(point_file,'w')
1541        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")
1542        geo_ref.write_ASCII(fd)
1543        fd.close()
1544
1545        mesh_output_file = tempfile.mktemp(".tsh")
1546        fit_to_mesh_file(mesh_file,
1547                         point_file,
1548                         mesh_output_file,
1549                         alpha = 0.0)
1550        # load in the .tsh file we just wrote
1551        mesh_dic = import_mesh_file(mesh_output_file)
1552        #print "mesh_dic",mesh_dic
1553        ans =[[0.0, 0.0],
1554              [5.0, 10.0],
1555              [5.0,10.0]]
1556        assert allclose(mesh_dic['vertex_attributes'],ans)
1557
1558        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1559                        ['elevation','stage'],
1560                        'test_fit_to_mesh_file failed')
1561
1562        #clean up
1563        os.remove(mesh_file)
1564        os.remove(point_file)
1565        os.remove(mesh_output_file)
1566
1567    def test_fit_to_mesh_fileII(self):
1568        from load_mesh.loadASCII import import_mesh_file, \
1569             export_mesh_file
1570        import tempfile
1571        import os
1572
1573        # create a .tsh file, no user outline
1574        mesh_dic = {}
1575        mesh_dic['vertices'] = [[0.0, 0.0],
1576                                [0.0, 5.0],
1577                                [5.0, 0.0]]
1578        mesh_dic['triangles'] =  [[0, 2, 1]]
1579        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1580        mesh_dic['triangle_tags'] = ['']
1581        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1582        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1583        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1584        mesh_dic['segment_tags'] = ['external',
1585                                                  'external',
1586                                                  'external']
1587        mesh_file = tempfile.mktemp(".tsh")
1588        export_mesh_file(mesh_file,mesh_dic)
1589
1590        # create an .xya file
1591        point_file = tempfile.mktemp(".xya")
1592        fd = open(point_file,'w')
1593        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")
1594        fd.close()
1595
1596        mesh_output_file = "new_triangle.tsh"
1597        fit_to_mesh_file(mesh_file,
1598                         point_file,
1599                         mesh_output_file,
1600                         alpha = 0.0)
1601        # load in the .tsh file we just wrote
1602        mesh_dic = import_mesh_file(mesh_output_file)
1603
1604        assert allclose(mesh_dic['vertex_attributes'],
1605                        [[1.0, 2.0,0.0, 0.0],
1606                         [1.0, 2.0,5.0, 10.0],
1607                         [1.0, 2.0,5.0,10.0]])
1608
1609        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1610                        ['density', 'temp','elevation','stage'],
1611                        'test_fit_to_mesh_file failed')
1612
1613        #clean up
1614        os.remove(mesh_file)
1615        os.remove(mesh_output_file)
1616        os.remove(point_file)
1617
1618    def test_fit_to_mesh_file_errors(self):
1619        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1620        import tempfile
1621        import os
1622
1623        # create a .tsh file, no user outline
1624        mesh_dic = {}
1625        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1626        mesh_dic['triangles'] =  [[0, 2, 1]]
1627        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1628        mesh_dic['triangle_tags'] = ['']
1629        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1630        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1631        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1632        mesh_dic['segment_tags'] = ['external', 'external','external']
1633        mesh_file = tempfile.mktemp(".tsh")
1634        export_mesh_file(mesh_file,mesh_dic)
1635
1636        # create an .xya file
1637        point_file = tempfile.mktemp(".xya")
1638        fd = open(point_file,'w')
1639        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")
1640        fd.close()
1641
1642        mesh_output_file = "new_triangle.tsh"
1643        try:
1644            fit_to_mesh_file(mesh_file, point_file,
1645                             mesh_output_file, display_errors = False)
1646        except IOError:
1647            pass
1648        else:
1649            #self.failUnless(0 ==1,  'Bad file did not raise error!')
1650            raise 'Bad file did not raise error!'
1651           
1652        #clean up
1653        os.remove(mesh_file)
1654        os.remove(point_file)
1655
1656    def test_fit_to_mesh_file_errorsII(self):
1657        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1658        import tempfile
1659        import os
1660
1661        # create a .tsh file, no user outline
1662        mesh_file = tempfile.mktemp(".tsh")
1663        fd = open(mesh_file,'w')
1664        fd.write("unit testing a bad .tsh file \n")
1665        fd.close()
1666
1667        # create an .xya file
1668        point_file = tempfile.mktemp(".xya")
1669        fd = open(point_file,'w')
1670        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")
1671        fd.close()
1672
1673        mesh_output_file = "new_triangle.tsh"
1674        try:
1675            fit_to_mesh_file(mesh_file, point_file,
1676                             mesh_output_file, display_errors = False)
1677        except IOError:
1678            pass
1679        else:
1680            raise 'Bad file did not raise error!'
1681           
1682        #clean up
1683        os.remove(mesh_file)
1684        os.remove(point_file)
1685
1686    def test_fit_to_mesh_file_errorsIII(self):
1687        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1688        import tempfile
1689        import os
1690
1691        # create a .tsh file, no user outline
1692        mesh_dic = {}
1693        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1694        mesh_dic['triangles'] =  [[0, 2, 1]]
1695        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1696        mesh_dic['triangle_tags'] = ['']
1697        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1698        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1699        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1700        mesh_dic['segment_tags'] = ['external', 'external','external']
1701        mesh_file = tempfile.mktemp(".tsh")
1702        export_mesh_file(mesh_file,mesh_dic)
1703
1704        # create an .xya file
1705        point_file = tempfile.mktemp(".xya")
1706        fd = open(point_file,'w')
1707        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")
1708        fd.close()
1709
1710        #This a deliberately illegal filename to invoke the error.
1711        mesh_output_file = ".../\z\z:ya.tsh"       
1712
1713        try:
1714            fit_to_mesh_file(mesh_file, point_file,
1715                             mesh_output_file, display_errors = False)
1716        except IOError:
1717            pass
1718        else:
1719            raise 'Bad file did not raise error!'
1720       
1721        #clean up
1722        os.remove(mesh_file)
1723        os.remove(point_file)
1724 
1725    def test_fit_to_mesh_file_errorsIV(self):
1726        import os
1727        command = '%s least_squares.py q q q e n 0.9 n'  %(sys.executable)
1728        status = os.system(command) 
1729        self.failUnless(status%255  == 1,
1730                        'command prompt least_squares.py failed.  Incorect exit status.')
1731       
1732    def test_fit_to_msh_netcdf_fileII(self):
1733        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1734        import tempfile
1735        import os
1736
1737        # create a .tsh file, no user outline
1738        mesh_dic = {}
1739        mesh_dic['vertices'] = [[0.0, 0.0],
1740                                [0.0, 5.0],
1741                                [5.0, 0.0]]
1742        mesh_dic['triangles'] =  [[0, 2, 1]]
1743        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1744        mesh_dic['triangle_tags'] = ['']
1745        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1746        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1747        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1748        mesh_dic['segment_tags'] = ['external',
1749                                                  'external',
1750                                                  'external']
1751        mesh_file = tempfile.mktemp(".msh")
1752        export_mesh_file(mesh_file,mesh_dic)
1753
1754        # create an .xya file
1755        point_file = tempfile.mktemp(".xya")
1756        fd = open(point_file,'w')
1757        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")
1758        fd.close()
1759
1760        mesh_output_file = "new_triangle.msh"
1761        fit_to_mesh_file(mesh_file,
1762                         point_file,
1763                         mesh_output_file,
1764                         alpha = 0.0)
1765        # load in the .tsh file we just wrote
1766        mesh_dic = import_mesh_file(mesh_output_file)
1767
1768        assert allclose(mesh_dic['vertex_attributes'],
1769                        [[1.0, 2.0,0.0, 0.0],
1770                         [1.0, 2.0,5.0, 10.0],
1771                         [1.0, 2.0,5.0,10.0]])
1772
1773        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1774                        ['density', 'temp','elevation','stage'],
1775                        'test_fit_to_mesh_file failed')
1776
1777        #clean up
1778        os.remove(mesh_file)
1779        os.remove(mesh_output_file)
1780        os.remove(point_file)
1781
1782#-------------------------------------------------------------
1783if __name__ == "__main__":
1784    suite = unittest.makeSuite(Test_Least_Squares,'test')
1785
1786    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
1787    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
1788    runner = unittest.TextTestRunner(verbosity=1)
1789    runner.run(suite)
1790
1791
1792
1793
1794
Note: See TracBrowser for help on using the repository browser.