source: anuga_core/source/anuga/utilities/test_numerical_tools.py @ 4551

Last change on this file since 4551 was 4551, checked in by ole, 17 years ago

Introduced maximal_inundation_elevation (runup) for sww files

File size: 8.6 KB
RevLine 
[1910]1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
[3514]7from anuga.config import epsilon
[1910]8
9from numerical_tools import *
10
11def test_function(x, y):
12    return x+y
13
14class Test_Numerical_Tools(unittest.TestCase):
15    def setUp(self):
16        pass
17
18    def tearDown(self):
19        pass
20
21
[2704]22    def test_angle1(self):
23        """Test angles between one vector and the x-axis
24        """
25        assert allclose(angle([1.0, 0.0])/pi*180, 0.0)     
[2526]26        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
[2704]27        assert allclose(angle([0.0, 1.0])/pi*180, 90.0)         
28        assert allclose(angle([-1.0, 1.0])/pi*180, 135.0)               
29        assert allclose(angle([-1.0, 0.0])/pi*180, 180.0)
30        assert allclose(angle([-1.0, -1.0])/pi*180, 225.0)
31        assert allclose(angle([0.0, -1.0])/pi*180, 270.0)
32        assert allclose(angle([1.0, -1.0])/pi*180, 315.0)               
33               
34                                                         
35    def test_angle2(self):
36        """Test angles between two arbitrary vectors
37        """   
38       
39        assert allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
40        assert allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
41               
42        assert allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)   
43        assert allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0) 
44       
45        assert allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
46        assert allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)   
47       
48        assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
[2709]49        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
50
51        #From test_get_boundary_polygon_V
52        v_prev = [-0.5, -0.5]
53        vc = [ 0.0,  -0.5]
54        assert allclose(angle(vc, v_prev)/pi*180, 45.0)
55
56        vc = [ 0.5,  0.0]
57        assert allclose(angle(vc, v_prev)/pi*180, 135.0)
58
59        vc = [ -0.5,  0.5]
60        assert allclose(angle(vc, v_prev)/pi*180, 270.0)               
61
62
[2704]63               
64       
65                               
66               
67       
68               
[1910]69
[2526]70
71    def test_anglediff(self):
72        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
73
74
75
[1910]76       
77    def test_ensure_numeric(self):
78        from numerical_tools import ensure_numeric
[4551]79        from Numeric import ArrayType, Float, Int, array
[1910]80
81        A = [1,2,3,4]
82        B = ensure_numeric(A)
83        assert type(B) == ArrayType
84        assert B.typecode() == 'l'
85        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
86
87
88        A = [1,2,3.14,4]
89        B = ensure_numeric(A)
90        assert type(B) == ArrayType
91        assert B.typecode() == 'd'
92        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
93
94
95        A = [1,2,3,4]
96        B = ensure_numeric(A, Float)
97        assert type(B) == ArrayType
98        assert B.typecode() == 'd'
99        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
100
101
102        A = [1,2,3,4]
103        B = ensure_numeric(A, Float)
104        assert type(B) == ArrayType
105        assert B.typecode() == 'd'
106        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
107
108
109        A = array([1,2,3,4])
110        B = ensure_numeric(A)
111        assert type(B) == ArrayType
112        assert B.typecode() == 'l'       
113        assert A == B   
114        assert A is B   #Same object
115
116
117        A = array([1,2,3,4])
118        B = ensure_numeric(A, Float)
119        assert type(B) == ArrayType
120        assert B.typecode() == 'd'       
121        assert A == B   
[2509]122        assert A is not B   #Not the same object
[1910]123
[4551]124        # Check scalars
125        A = 1
126        B = ensure_numeric(A, Float)
127        #print A, B[0], len(B), type(B)
128        #print B.shape
129        assert A == B
[1910]130
[4551]131        B = ensure_numeric(A, Int)       
132        #print A, B
133        #print B.shape
134        assert A == B
135
136        # Error situation
137
138        B = ensure_numeric('hello', Int)               
139        assert allclose(B, [104, 101, 108, 108, 111])
140
[2509]141    def test_gradient(self):
142        x0 = 0.0; y0 = 0.0; z0 = 0.0
143        x1 = 1.0; y1 = 0.0; z1 = -1.0
144        x2 = 0.0; y2 = 1.0; z2 = 0.0
145
146        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
147
148        assert zx == -1.0
149        assert zy == 0.0
150
151    def test_gradient_more(self):
152        x0 = 2.0/3; y0 = 2.0/3
153        x1=  8.0/3; y1 = 2.0/3
154        x2 = 2.0/3; y2 = 8.0/3
155
156        q0 = 2.0+2.0/3
157        q1 = 8.0+2.0/3
158        q2 = 2.0+8.0/3
159
160        #Gradient of fitted pwl surface
161        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
162
163        assert abs(a - 3.0) < epsilon
164        assert abs(b - 1.0) < epsilon
165
166
167    def test_gradient2(self):
168        """Test two-point gradient
169        """
170       
171        x0 = 5.0; y0 = 5.0; z0 = 10.0
172        x1 = 8.0; y1 = 2.0; z1 = 1.0
173        x2 = 8.0; y2 = 8.0; z2 = 10.0
174
175        #Reference
176        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
177        a, b = gradient2(x0, y0, x1, y1, z0, z1)
178
179        assert zx == a
180        assert zy == b
181
182        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
183        assert z2_computed == z2
184       
185    def test_gradient2_more(self):
186        """Test two-point gradient more
187        """
188        x0 = 2.0; y0 = 2.0
189        x1 = 8.0; y1 = 3.0
190        x2 = 1.0; y2 = 8.0
191
192        q0 = 2.0
193        q1 = 8.0
194        q2 = q0
195
196        #Gradient of fitted pwl surface
197        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
198        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
199
200        assert a == a_ref
201        assert b == b_ref
202
203
[3849]204    def test_machine_precision(self):
205        """test_machine_precision(self):
206        Test the function that calculates epsilon. As this varies on
207        different machines, this is only an indication.
208        """
209
210        eps = get_machine_precision()
211
212        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
213        assert eps > 0.0
214        assert 1.0+eps/2 == 1.0
215       
216       
[2533]217    def test_histogram(self):
218        """Test histogram with different bin boundaries
219        """
220       
221        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
222
223
224        #There are four elements greater than or equal to 3
225        bins = [3]
226        assert allclose(histogram(a, bins), [4])
227
228
229        bins = [ min(a) ]
230        assert allclose(histogram(a, bins), [len(a)])
231
232
233        bins = [ max(a)+0.00001 ]
234        assert allclose(histogram(a, bins), [0])       
235
236       
237        bins = [1,2,3,4]
238        assert allclose(histogram(a, bins), [8,3,3,1])
239
240
[2534]241        bins = [1.1,2,3.1,4]
242        #print histogram(a, bins)
243        assert allclose(histogram(a, bins), [0,6,0,1])
244
245
246        bins = [0,1.5,2,3]
[2533]247        assert allclose(histogram(a, bins), [8,0,3,4])
[2971]248        assert allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
[2533]249
[2971]250        # Check situation with #bins >= #datapoints
251        a = [1.7]
252        bins = [0,1.5,2,3]
253        assert allclose(histogram(a, bins), [0,1,0,0])
[2533]254
[2971]255        a = [1.7]
256        bins = [0]
257        assert allclose(histogram(a, bins), [1])
[2534]258
[2971]259        a = [-1.7]
260        bins = [0]
261        assert allclose(histogram(a, bins), [0])
[2534]262
[2971]263        a = [-1.7]
264        bins = [-1.7]
265        assert allclose(histogram(a, bins), [1])
[2533]266       
[2971]267       
[2533]268
[2509]269    def test_that_C_extension_compiles(self):
270        FN = 'util_ext.c'
271        try:
272            import util_ext
273        except:
274            from compile import compile
275
276            try:
277                compile(FN)
278            except:
279                raise 'Could not compile %s' %FN
280            else:
281                import util_ext
282
283
284    def test_gradient_C_extension(self):
285        from util_ext import gradient as gradient_c
286
287        x0 = 2.0/3; y0 = 2.0/3
288        x1=  8.0/3; y1 = 2.0/3
289        x2 = 2.0/3; y2 = 8.0/3
290
291        q0 = 2.0+2.0/3
292        q1 = 8.0+2.0/3
293        q2 = 2.0+8.0/3
294
295        #Gradient of fitted pwl surface
296        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
297
298        assert abs(a - 3.0) < epsilon
299        assert abs(b - 1.0) < epsilon
300
301
302    def test_gradient_C_extension3(self):
303        from util_ext import gradient as gradient_c
304
305        from RandomArray import uniform, seed
306        seed(17, 53)
307
308        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
309
310        q0 = uniform(0.0, 10.0, 4)
311        q1 = uniform(1.0, 3.0, 4)
312        q2 = uniform(7.0, 20.0, 4)
313
314
315        for i in range(4):
316            #Gradient of fitted pwl surface
[2516]317            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
318                                           q0[i], q1[i], q2[i])
[2509]319
320            #print a_ref, b_ref
321            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
322                              q0[i], q1[i], q2[i])
323
324            #print a, a_ref, b, b_ref
325            assert abs(a - a_ref) < epsilon
326            assert abs(b - b_ref) < epsilon
327
328       
329
330
[1910]331#-------------------------------------------------------------
332if __name__ == "__main__":
333    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
334    runner = unittest.TextTestRunner()
335    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.