source: inundation/fit_interpolate/test_interpolate.py @ 2263

Last change on this file since 2263 was 2201, checked in by duncan, 19 years ago

comments

File size: 16.6 KB
RevLine 
[2187]1#!/usr/bin/env python
2
3#TEST
4import sys
5import unittest
6from math import sqrt
7
8
9from interpolate 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_Interpolate(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 = Interpolate(points, vertices)
40        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
41                        [[1./3, 1./3, 1./3]])
42
43
44    def test_quad_tree(self):
45        p0 = [-10.0, -10.0]
46        p1 = [20.0, -10.0]
47        p2 = [-10.0, 20.0]
48        p3 = [10.0, 50.0]
49        p4 = [30.0, 30.0]
50        p5 = [50.0, 10.0]
51        p6 = [40.0, 60.0]
52        p7 = [60.0, 40.0]
53        p8 = [-66.0, 20.0]
54        p9 = [10.0, -66.0]
55
56        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
57        triangles = [ [0, 1, 2],
58                      [3, 2, 4],
59                      [4, 2, 1],
60                      [4, 1, 5],
61                      [3, 4, 6],
62                      [6, 4, 7],
63                      [7, 4, 5],
64                      [8, 0, 2],
65                      [0, 9, 1]]
66
67        data = [ [4,4] ]
68        interp = Interpolate(points, triangles,
[2201]69                               max_vertices_per_cell = 4)
[2187]70        #print "PDSG - interp.get_A()", interp.get_A()
71        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
72                      0., 0. , 0., 0., 0., 0.]]
73
74       
75        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
76                        answer)
77        #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
78        #print "PDSG - interp.get_A()", interp.get_A()
79        data = [[-30, -30]]
80        answer =  [ [ 0.0,  0.0,  0.0,  0.,
81                      0., 0. , 0., 0., 0., 0.]]
82       
83        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
84                        answer)
85
86
87        #point outside of quad tree root cell
88        #interp.set_point_coordinates([[-70, -70]])
89        #print "PDSG - interp.get_A()", interp.get_A()
90        data = [[-70, -70]]
91        answer =  [ [ 0.0,  0.0,  0.0,  0.,
92                      0., 0. , 0., 0., 0., 0.]]
93        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
94                        answer)
95
96    def test_datapoints_at_vertices(self):
97        """Test that data points coinciding with vertices yield a diagonal matrix
98        """
99
100        a = [0.0, 0.0]
101        b = [0.0, 2.0]
102        c = [2.0,0.0]
103        points = [a, b, c]
104        vertices = [ [1,0,2] ]   #bac
105
106        data = points #Use data at vertices
107
108        interp = Interpolate(points, vertices)
109        answer = [[1., 0., 0.],
110                   [0., 1., 0.],
111                   [0., 0., 1.]]
112        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
113                        answer)
114
115
116    def test_datapoints_on_edge_midpoints(self):
117        """Try datapoints midway on edges -
118        each point should affect two matrix entries equally
119        """
120
121        a = [0.0, 0.0]
122        b = [0.0, 2.0]
123        c = [2.0,0.0]
124        points = [a, b, c]
125        vertices = [ [1,0,2] ]   #bac
126
127        data = [ [0., 1.], [1., 0.], [1., 1.] ]
128        answer =  [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
129                    [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
130                    [0.0, 0.5, 0.5]]
131        interp = Interpolate(points, vertices, data)
132
133        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
134                        answer)
135
136    def test_datapoints_on_edges(self):
137        """Try datapoints on edges -
138        each point should affect two matrix entries in proportion
139        """
140
141        a = [0.0, 0.0]
142        b = [0.0, 2.0]
143        c = [2.0,0.0]
144        points = [a, b, c]
145        vertices = [ [1,0,2] ]   #bac
146
147        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
148        answer =  [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
149                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
150                   [0.0, 0.25, 0.75]]
151
152        interp = Interpolate(points, vertices, data)
153
154        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
155                        answer)
156
157
158    def test_arbitrary_datapoints(self):
159        """Try arbitrary datapoints
160        """
161
162        from Numeric import sum
163
164        a = [0.0, 0.0]
165        b = [0.0, 2.0]
166        c = [2.0,0.0]
167        points = [a, b, c]
168        vertices = [ [1,0,2] ]   #bac
169
170        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
171
172        interp = Interpolate(points, vertices, data)
173        #print "interp.get_A()", interp.get_A()
174        results = interp._build_interpolation_matrix_A(data).todense()
175        assert allclose(sum(results, axis=1), 1.0)
176
177        #FIXME - have to change this test to check default info
178    def NO_test_arbitrary_datapoints_some_outside(self):
179        """Try arbitrary datapoints one outside the triangle.
180        That one should be ignored
181        """
182
183        from Numeric import sum
184
185        a = [0.0, 0.0]
186        b = [0.0, 2.0]
187        c = [2.0,0.0]
188        points = [a, b, c]
189        vertices = [ [1,0,2] ]   #bac
190
191        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
192
193
194        interp = Interpolate(points, vertices, data, precrop = True)
195       
196        results = interp._build_interpolation_matrix_A(data).todense()
197        assert allclose(sum(results, axis=1), 1.0)
198
199        interp = Interpolate(points, vertices, data, precrop = False)
200        results = interp._build_interpolation_matrix_A(data).todense()
201        assert allclose(sum(results, axis=1), [1,1,1,0])
202
203
204
205    # this causes a memory error in scipy.sparse
206    def test_more_triangles(self):
207
208        a = [-1.0, 0.0]
209        b = [3.0, 4.0]
210        c = [4.0,1.0]
211        d = [-3.0, 2.0] #3
212        e = [-1.0,-2.0]
213        f = [1.0, -2.0] #5
214
215        points = [a, b, c, d,e,f]
216        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
217
218        #Data points
219        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
220        interp = Interpolate(points, triangles)
221
222        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
223                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
224                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
225                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
226                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
227                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
228
229
230        A = interp._build_interpolation_matrix_A(data).todense()
231        for i in range(A.shape[0]):
232            for j in range(A.shape[1]):
233                if not allclose(A[i,j], answer[i][j]):
234                    print i,j,':',A[i,j], answer[i][j]
235
236
[2190]237        #results = interp._build_interpolation_matrix_A(data).todense()
[2187]238
[2190]239        assert allclose(A, answer)
[2187]240
[2190]241
[2189]242    def test_interpolate_attributes_to_points(self):
243        v0 = [0.0, 0.0]
244        v1 = [0.0, 5.0]
245        v2 = [5.0, 0.0]
[2187]246
[2189]247        vertices = [v0, v1, v2]
248        triangles = [ [1,0,2] ]   #bac
[2187]249
[2189]250        d0 = [1.0, 1.0]
251        d1 = [1.0, 2.0]
252        d2 = [3.0, 1.0]
253        point_coords = [ d0, d1, d2]
254
255        interp = Interpolate(vertices, triangles, point_coords)
256        f = linear_function(vertices)
257        z = interp.interpolate(f, point_coords)
258        answer = linear_function(point_coords)
259
260
261        assert allclose(z, answer)
262
263
264
265    def test_interpolate_attributes_to_pointsII(self):
266        a = [-1.0, 0.0]
267        b = [3.0, 4.0]
268        c = [4.0, 1.0]
269        d = [-3.0, 2.0] #3
270        e = [-1.0, -2.0]
271        f = [1.0, -2.0] #5
272
273        vertices = [a, b, c, d,e,f]
274        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
275
276
277        point_coords = [[-2.0, 2.0],
278                        [-1.0, 1.0],
279                        [0.0, 2.0],
280                        [1.0, 1.0],
281                        [2.0, 1.0],
282                        [0.0, 0.0],
283                        [1.0, 0.0],
284                        [0.0, -1.0],
285                        [-0.2, -0.5],
286                        [-0.9, -1.5],
287                        [0.5, -1.9],
288                        [3.0, 1.0]]
289
290        interp = Interpolate(vertices, triangles)
291        f = linear_function(vertices)
292        z = interp.interpolate(f, point_coords)
293        answer = linear_function(point_coords)
294        #print "z",z
295        #print "answer",answer
296        assert allclose(z, answer)
297
298    def test_interpolate_attributes_to_pointsIII(self):
299        """Test linear interpolation of known values at vertices to
300        new points inside a triangle
301        """
302        a = [0.0, 0.0]
303        b = [0.0, 5.0]
304        c = [5.0, 0.0]
305        d = [5.0, 5.0]
306
307        vertices = [a, b, c, d]
308        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
309
310        #Points within triangle 1
311        d0 = [1.0, 1.0]
312        d1 = [1.0, 2.0]
313        d2 = [3.0, 1.0]
314
315        #Point within triangle 2
316        d3 = [4.0, 3.0]
317
318        #Points on common edge
319        d4 = [2.5, 2.5]
320        d5 = [4.0, 1.0]
321
322        #Point on common vertex
323        d6 = [0., 5.]
324       
325        point_coords = [d0, d1, d2, d3, d4, d5, d6]
326
327        interp = Interpolate(vertices, triangles)
328
329        #Known values at vertices
330        #Functions are x+y, x+2y, 2x+y, x-y-5
331        f = [ [0., 0., 0., -5.],        # (0,0)
332              [5., 10., 5., -10.],      # (0,5)
333              [5., 5., 10.0, 0.],       # (5,0)
334              [10., 15., 15., -5.]]     # (5,5)
335
336        z = interp.interpolate(f, point_coords)
337        answer = [ [2., 3., 3., -5.],   # (1,1)
338                   [3., 5., 4., -6.],   # (1,2)
339                   [4., 5., 7., -3.],   # (3,1)
340                   [7., 10., 11., -4.], # (4,3)
341                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
342                   [5., 6., 9., -2.],   # (4,1)
343                   [5., 10., 5., -10.]]  # (0,5)
344
345        #print "***********"
346        #print "z",z
347        #print "answer",answer
348        #print "***********"
349
350        #Should an error message be returned if points are outside
351        # of the mesh? Not currently. 
352
353        assert allclose(z, answer)
354
355
356    def test_interpolate_point_outside_of_mesh(self):
357        """Test linear interpolation of known values at vertices to
358        new points inside a triangle
359        """
360        a = [0.0, 0.0]
361        b = [0.0, 5.0]
362        c = [5.0, 0.0]
363        d = [5.0, 5.0]
364
365        vertices = [a, b, c, d]
366        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
367
368        #Far away point
369        d7 = [-1., -1.]
370       
371        point_coords = [ d7]
372        interp = Interpolate(vertices, triangles)
373
374        #Known values at vertices
375        #Functions are x+y, x+2y, 2x+y, x-y-5
376        f = [ [0., 0., 0., -5.],        # (0,0)
377              [5., 10., 5., -10.],      # (0,5)
378              [5., 5., 10.0, 0.],       # (5,0)
379              [10., 15., 15., -5.]]     # (5,5)
380
381        z = interp.interpolate(f, point_coords)
382        answer = [ [0., 0., 0., 0.]] # (-1,-1)
383
384        #print "***********"
385        #print "z",z
386        #print "answer",answer
387        #print "***********"
388
389        #Should an error message be returned if points are outside
390        # of the mesh? Not currently. 
391
392        assert allclose(z, answer)
393
394    def test_interpolate_attributes_to_pointsIV(self):
395        a = [-1.0, 0.0]
396        b = [3.0, 4.0]
397        c = [4.0, 1.0]
398        d = [-3.0, 2.0] #3
399        e = [-1.0, -2.0]
400        f = [1.0, -2.0] #5
401
402        vertices = [a, b, c, d,e,f]
403        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
404
405
406        point_coords = [[-2.0, 2.0],
407                        [-1.0, 1.0],
408                        [0.0, 2.0],
409                        [1.0, 1.0],
410                        [2.0, 1.0],
411                        [0.0, 0.0],
412                        [1.0, 0.0],
413                        [0.0, -1.0],
414                        [-0.2, -0.5],
415                        [-0.9, -1.5],
416                        [0.5, -1.9],
417                        [3.0, 1.0]]
418
419        interp = Interpolate(vertices, triangles)
420        f = array([linear_function(vertices),2*linear_function(vertices) ])
421        f = transpose(f)
422        #print "f",f
423        z = interp.interpolate(f, point_coords)
424        answer = [linear_function(point_coords),
425                  2*linear_function(point_coords) ]
426        answer = transpose(answer)
427        #print "z",z
428        #print "answer",answer
429        assert allclose(z, answer)
430
431
432    def test_interpolate_blocking(self):
433        a = [-1.0, 0.0]
434        b = [3.0, 4.0]
435        c = [4.0, 1.0]
436        d = [-3.0, 2.0] #3
437        e = [-1.0, -2.0]
438        f = [1.0, -2.0] #5
439
440        vertices = [a, b, c, d,e,f]
441        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
442
443
444        point_coords = [[-2.0, 2.0],
445                        [-1.0, 1.0],
446                        [0.0, 2.0],
447                        [1.0, 1.0],
448                        [2.0, 1.0],
449                        [0.0, 0.0],
450                        [1.0, 0.0],
451                        [0.0, -1.0],
452                        [-0.2, -0.5],
453                        [-0.9, -1.5],
454                        [0.5, -1.9],
455                        [3.0, 1.0]]
456
457        interp = Interpolate(vertices, triangles)
458        f = array([linear_function(vertices),2*linear_function(vertices) ])
459        f = transpose(f)
460        #print "f",f
461        for blocking_max in range(len(point_coords)+2):
462        #if True:
463         #   blocking_max = 5
464            z = interp.interpolate(f, point_coords,
465                                   start_blocking_len=blocking_max)
466            answer = [linear_function(point_coords),
467                      2*linear_function(point_coords) ]
468            answer = transpose(answer)
469            #print "z",z
470            #print "answer",answer
471            assert allclose(z, answer)
472
473    def test_interpolate_reuse(self):
474        a = [-1.0, 0.0]
475        b = [3.0, 4.0]
476        c = [4.0, 1.0]
477        d = [-3.0, 2.0] #3
478        e = [-1.0, -2.0]
479        f = [1.0, -2.0] #5
480
481        vertices = [a, b, c, d,e,f]
482        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
483
484
485        point_coords = [[-2.0, 2.0],
486                        [-1.0, 1.0],
487                        [0.0, 2.0],
488                        [1.0, 1.0],
489                        [2.0, 1.0],
490                        [0.0, 0.0],
491                        [1.0, 0.0],
492                        [0.0, -1.0],
493                        [-0.2, -0.5],
494                        [-0.9, -1.5],
495                        [0.5, -1.9],
496                        [3.0, 1.0]]
497
498        interp = Interpolate(vertices, triangles)
499        f = array([linear_function(vertices),2*linear_function(vertices) ])
500        f = transpose(f)
501        z = interp.interpolate(f, point_coords,
502                               start_blocking_len=20)
503        answer = [linear_function(point_coords),
504                  2*linear_function(point_coords) ]
505        answer = transpose(answer)
506        #print "z",z
507        #print "answer",answer
508        assert allclose(z, answer)
[2201]509        assert allclose(interp._A_can_be_reused, True)
[2189]510
511        z = interp.interpolate(f)
512        assert allclose(z, answer)
513       
514        # This causes blocking to occur.
515        z = interp.interpolate(f, start_blocking_len=10)
516        assert allclose(z, answer)
[2201]517        assert allclose(interp._A_can_be_reused, False)
[2189]518
519        #A is recalculated
520        z = interp.interpolate(f)
521        assert allclose(z, answer)
[2201]522        assert allclose(interp._A_can_be_reused, True)
523       
[2189]524        interp = Interpolate(vertices, triangles)
525        #Must raise an exception, no points specified
526        try:
527            z = interp.interpolate(f)
528        except:
529            pass
530       
531
[2187]532#-------------------------------------------------------------
533if __name__ == "__main__":
534    suite = unittest.makeSuite(Test_Interpolate,'test')
535    runner = unittest.TextTestRunner(verbosity=1)
536    runner.run(suite)
537
538
539
540
541
Note: See TracBrowser for help on using the repository browser.