source: inundation/pyvolution/test_least_squares.py @ 2141

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

compiling triangle using setup.
File size: 54.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        #FIXME - make this test the georef in the points file as well.
1486        # create an .xya file
1487        point_file = tempfile.mktemp(".xya")
1488        fd = open(point_file,'w')
1489        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")
1490        fd.close()
1491
1492        mesh_output_file = tempfile.mktemp(".tsh")
1493        fit_to_mesh_file(mesh_file,
1494                         point_file,
1495                         mesh_output_file,
1496                         alpha = 0.0)
1497        # load in the .tsh file we just wrote
1498        mesh_dic = import_mesh_file(mesh_output_file)
1499        #print "mesh_dic",mesh_dic
1500        ans =[[0.0, 0.0],
1501              [5.0, 10.0],
1502              [5.0,10.0]]
1503        assert allclose(mesh_dic['vertex_attributes'],ans)
1504
1505        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1506                        ['elevation','stage'],
1507                        'test_fit_to_mesh_file failed')
1508
1509        #clean up
1510        os.remove(mesh_file)
1511        os.remove(point_file)
1512        os.remove(mesh_output_file)
1513
1514    def test_fit_to_mesh_fileII(self):
1515        from load_mesh.loadASCII import import_mesh_file, \
1516             export_mesh_file
1517        import tempfile
1518        import os
1519
1520        # create a .tsh file, no user outline
1521        mesh_dic = {}
1522        mesh_dic['vertices'] = [[0.0, 0.0],
1523                                [0.0, 5.0],
1524                                [5.0, 0.0]]
1525        mesh_dic['triangles'] =  [[0, 2, 1]]
1526        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1527        mesh_dic['triangle_tags'] = ['']
1528        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1529        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1530        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1531        mesh_dic['segment_tags'] = ['external',
1532                                                  'external',
1533                                                  'external']
1534        mesh_file = tempfile.mktemp(".tsh")
1535        export_mesh_file(mesh_file,mesh_dic)
1536
1537        # create an .xya file
1538        point_file = tempfile.mktemp(".xya")
1539        fd = open(point_file,'w')
1540        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")
1541        fd.close()
1542
1543        mesh_output_file = "new_triangle.tsh"
1544        fit_to_mesh_file(mesh_file,
1545                         point_file,
1546                         mesh_output_file,
1547                         alpha = 0.0)
1548        # load in the .tsh file we just wrote
1549        mesh_dic = import_mesh_file(mesh_output_file)
1550
1551        assert allclose(mesh_dic['vertex_attributes'],
1552                        [[1.0, 2.0,0.0, 0.0],
1553                         [1.0, 2.0,5.0, 10.0],
1554                         [1.0, 2.0,5.0,10.0]])
1555
1556        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1557                        ['density', 'temp','elevation','stage'],
1558                        'test_fit_to_mesh_file failed')
1559
1560        #clean up
1561        os.remove(mesh_file)
1562        os.remove(mesh_output_file)
1563        os.remove(point_file)
1564
1565    def test_fit_to_mesh_file_errors(self):
1566        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1567        import tempfile
1568        import os
1569
1570        # create a .tsh file, no user outline
1571        mesh_dic = {}
1572        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1573        mesh_dic['triangles'] =  [[0, 2, 1]]
1574        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1575        mesh_dic['triangle_tags'] = ['']
1576        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1577        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1578        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1579        mesh_dic['segment_tags'] = ['external', 'external','external']
1580        mesh_file = tempfile.mktemp(".tsh")
1581        export_mesh_file(mesh_file,mesh_dic)
1582
1583        # create an .xya file
1584        point_file = tempfile.mktemp(".xya")
1585        fd = open(point_file,'w')
1586        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")
1587        fd.close()
1588
1589        mesh_output_file = "new_triangle.tsh"
1590        try:
1591            fit_to_mesh_file(mesh_file, point_file,
1592                             mesh_output_file, display_errors = False)
1593        except IOError:
1594            pass
1595        else:
1596            #self.failUnless(0 ==1,  'Bad file did not raise error!')
1597            raise 'Bad file did not raise error!'
1598           
1599        #clean up
1600        os.remove(mesh_file)
1601        os.remove(point_file)
1602
1603    def test_fit_to_mesh_file_errorsII(self):
1604        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1605        import tempfile
1606        import os
1607
1608        # create a .tsh file, no user outline
1609        mesh_file = tempfile.mktemp(".tsh")
1610        fd = open(mesh_file,'w')
1611        fd.write("unit testing a bad .tsh file \n")
1612        fd.close()
1613
1614        # create an .xya file
1615        point_file = tempfile.mktemp(".xya")
1616        fd = open(point_file,'w')
1617        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")
1618        fd.close()
1619
1620        mesh_output_file = "new_triangle.tsh"
1621        try:
1622            fit_to_mesh_file(mesh_file, point_file,
1623                             mesh_output_file, display_errors = False)
1624        except IOError:
1625            pass
1626        else:
1627            raise 'Bad file did not raise error!'
1628           
1629        #clean up
1630        os.remove(mesh_file)
1631        os.remove(point_file)
1632
1633    def test_fit_to_mesh_file_errorsIII(self):
1634        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1635        import tempfile
1636        import os
1637
1638        # create a .tsh file, no user outline
1639        mesh_dic = {}
1640        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1641        mesh_dic['triangles'] =  [[0, 2, 1]]
1642        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1643        mesh_dic['triangle_tags'] = ['']
1644        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1645        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1646        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1647        mesh_dic['segment_tags'] = ['external', 'external','external']
1648        mesh_file = tempfile.mktemp(".tsh")
1649        export_mesh_file(mesh_file,mesh_dic)
1650
1651        # create an .xya file
1652        point_file = tempfile.mktemp(".xya")
1653        fd = open(point_file,'w')
1654        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")
1655        fd.close()
1656
1657        #This a deliberately illegal filename to invoke the error.
1658        mesh_output_file = ".../\z\z:ya.tsh"       
1659
1660        try:
1661            fit_to_mesh_file(mesh_file, point_file,
1662                             mesh_output_file, display_errors = False)
1663        except IOError:
1664            pass
1665        else:
1666            raise 'Bad file did not raise error!'
1667       
1668        #clean up
1669        os.remove(mesh_file)
1670        os.remove(point_file)
1671 
1672    def test_fit_to_mesh_file_errorsIV(self):
1673        import os
1674        command = '%s least_squares.py q q q e n 0.9 n'  %(sys.executable)
1675        status = os.system(command) 
1676        self.failUnless(status%255  == 1,
1677                        'command prompt least_squares.py failed.  Incorect exit status.')
1678       
1679    def test_fit_to_msh_netcdf_fileII(self):
1680        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1681        import tempfile
1682        import os
1683
1684        # create a .tsh file, no user outline
1685        mesh_dic = {}
1686        mesh_dic['vertices'] = [[0.0, 0.0],
1687                                [0.0, 5.0],
1688                                [5.0, 0.0]]
1689        mesh_dic['triangles'] =  [[0, 2, 1]]
1690        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1691        mesh_dic['triangle_tags'] = ['']
1692        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1693        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1694        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1695        mesh_dic['segment_tags'] = ['external',
1696                                                  'external',
1697                                                  'external']
1698        mesh_file = tempfile.mktemp(".msh")
1699        export_mesh_file(mesh_file,mesh_dic)
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.msh"
1708        fit_to_mesh_file(mesh_file,
1709                         point_file,
1710                         mesh_output_file,
1711                         alpha = 0.0)
1712        # load in the .tsh file we just wrote
1713        mesh_dic = import_mesh_file(mesh_output_file)
1714
1715        assert allclose(mesh_dic['vertex_attributes'],
1716                        [[1.0, 2.0,0.0, 0.0],
1717                         [1.0, 2.0,5.0, 10.0],
1718                         [1.0, 2.0,5.0,10.0]])
1719
1720        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1721                        ['density', 'temp','elevation','stage'],
1722                        'test_fit_to_mesh_file failed')
1723
1724        #clean up
1725        os.remove(mesh_file)
1726        os.remove(mesh_output_file)
1727        os.remove(point_file)
1728
1729#-------------------------------------------------------------
1730if __name__ == "__main__":
1731    suite = unittest.makeSuite(Test_Least_Squares,'test')
1732
1733    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
1734    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
1735    runner = unittest.TextTestRunner(verbosity=1)
1736    runner.run(suite)
1737
1738
1739
1740
1741
Note: See TracBrowser for help on using the repository browser.