source: branches/numpy/obsolete_code/test_least_squares.py @ 8353

Last change on this file since 8353 was 3564, checked in by ole, 18 years ago

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