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

Last change on this file since 4150 was 3849, checked in by ole, 18 years ago

Added function and test for calculating machine precision for Python Floats

File size: 8.2 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7from anuga.config import epsilon
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
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)     
26        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
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)   
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
63               
64       
65                               
66               
67       
68               
69
70
71    def test_anglediff(self):
72        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
73
74
75
76       
77    def test_ensure_numeric(self):
78        from numerical_tools import ensure_numeric
79        from Numeric import ArrayType, Float, array
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   
122        assert A is not B   #Not the same object
123
124
125    def test_gradient(self):
126        x0 = 0.0; y0 = 0.0; z0 = 0.0
127        x1 = 1.0; y1 = 0.0; z1 = -1.0
128        x2 = 0.0; y2 = 1.0; z2 = 0.0
129
130        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
131
132        assert zx == -1.0
133        assert zy == 0.0
134
135    def test_gradient_more(self):
136        x0 = 2.0/3; y0 = 2.0/3
137        x1=  8.0/3; y1 = 2.0/3
138        x2 = 2.0/3; y2 = 8.0/3
139
140        q0 = 2.0+2.0/3
141        q1 = 8.0+2.0/3
142        q2 = 2.0+8.0/3
143
144        #Gradient of fitted pwl surface
145        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
146
147        assert abs(a - 3.0) < epsilon
148        assert abs(b - 1.0) < epsilon
149
150
151    def test_gradient2(self):
152        """Test two-point gradient
153        """
154       
155        x0 = 5.0; y0 = 5.0; z0 = 10.0
156        x1 = 8.0; y1 = 2.0; z1 = 1.0
157        x2 = 8.0; y2 = 8.0; z2 = 10.0
158
159        #Reference
160        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
161        a, b = gradient2(x0, y0, x1, y1, z0, z1)
162
163        assert zx == a
164        assert zy == b
165
166        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
167        assert z2_computed == z2
168       
169    def test_gradient2_more(self):
170        """Test two-point gradient more
171        """
172        x0 = 2.0; y0 = 2.0
173        x1 = 8.0; y1 = 3.0
174        x2 = 1.0; y2 = 8.0
175
176        q0 = 2.0
177        q1 = 8.0
178        q2 = q0
179
180        #Gradient of fitted pwl surface
181        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
182        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
183
184        assert a == a_ref
185        assert b == b_ref
186
187
188    def test_machine_precision(self):
189        """test_machine_precision(self):
190        Test the function that calculates epsilon. As this varies on
191        different machines, this is only an indication.
192        """
193
194        eps = get_machine_precision()
195
196        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
197        assert eps > 0.0
198        assert 1.0+eps/2 == 1.0
199       
200       
201    def test_histogram(self):
202        """Test histogram with different bin boundaries
203        """
204       
205        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
206
207
208        #There are four elements greater than or equal to 3
209        bins = [3]
210        assert allclose(histogram(a, bins), [4])
211
212
213        bins = [ min(a) ]
214        assert allclose(histogram(a, bins), [len(a)])
215
216
217        bins = [ max(a)+0.00001 ]
218        assert allclose(histogram(a, bins), [0])       
219
220       
221        bins = [1,2,3,4]
222        assert allclose(histogram(a, bins), [8,3,3,1])
223
224
225        bins = [1.1,2,3.1,4]
226        #print histogram(a, bins)
227        assert allclose(histogram(a, bins), [0,6,0,1])
228
229
230        bins = [0,1.5,2,3]
231        assert allclose(histogram(a, bins), [8,0,3,4])
232        assert allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
233
234        # Check situation with #bins >= #datapoints
235        a = [1.7]
236        bins = [0,1.5,2,3]
237        assert allclose(histogram(a, bins), [0,1,0,0])
238
239        a = [1.7]
240        bins = [0]
241        assert allclose(histogram(a, bins), [1])
242
243        a = [-1.7]
244        bins = [0]
245        assert allclose(histogram(a, bins), [0])
246
247        a = [-1.7]
248        bins = [-1.7]
249        assert allclose(histogram(a, bins), [1])
250       
251       
252
253    def test_that_C_extension_compiles(self):
254        FN = 'util_ext.c'
255        try:
256            import util_ext
257        except:
258            from compile import compile
259
260            try:
261                compile(FN)
262            except:
263                raise 'Could not compile %s' %FN
264            else:
265                import util_ext
266
267
268    def test_gradient_C_extension(self):
269        from util_ext import gradient as gradient_c
270
271        x0 = 2.0/3; y0 = 2.0/3
272        x1=  8.0/3; y1 = 2.0/3
273        x2 = 2.0/3; y2 = 8.0/3
274
275        q0 = 2.0+2.0/3
276        q1 = 8.0+2.0/3
277        q2 = 2.0+8.0/3
278
279        #Gradient of fitted pwl surface
280        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
281
282        assert abs(a - 3.0) < epsilon
283        assert abs(b - 1.0) < epsilon
284
285
286    def test_gradient_C_extension3(self):
287        from util_ext import gradient as gradient_c
288
289        from RandomArray import uniform, seed
290        seed(17, 53)
291
292        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
293
294        q0 = uniform(0.0, 10.0, 4)
295        q1 = uniform(1.0, 3.0, 4)
296        q2 = uniform(7.0, 20.0, 4)
297
298
299        for i in range(4):
300            #Gradient of fitted pwl surface
301            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
302                                           q0[i], q1[i], q2[i])
303
304            #print a_ref, b_ref
305            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
306                              q0[i], q1[i], q2[i])
307
308            #print a, a_ref, b, b_ref
309            assert abs(a - a_ref) < epsilon
310            assert abs(b - b_ref) < epsilon
311
312       
313
314
315#-------------------------------------------------------------
316if __name__ == "__main__":
317    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
318    runner = unittest.TextTestRunner()
319    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.