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

Last change on this file since 5897 was 5897, checked in by ole, 15 years ago

Reverted numpy changes to the trunk that should have been made to the branch.
The command was svn merge -r 5895:5890 .
File size: 10.0 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from Numeric import ArrayType, Float, Int, array, alltrue
7
8from math import sqrt, pi
9from anuga.config import epsilon
10from numerical_tools import *
11
12def test_function(x, y):
13    return x+y
14
15class Test_Numerical_Tools(unittest.TestCase):
16    def setUp(self):
17        pass
18
19    def tearDown(self):
20        pass
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    def test_anglediff(self):
64        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
65
66       
67    def test_ensure_numeric(self):
68        A = [1,2,3,4]
69        B = ensure_numeric(A)
70        assert type(B) == ArrayType
71        assert B.typecode() == 'l'
72        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
73
74        A = [1,2,3.14,4]
75        B = ensure_numeric(A)
76        assert type(B) == ArrayType
77        assert B.typecode() == 'd'
78        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
79
80        A = [1,2,3,4]
81        B = ensure_numeric(A, Float)
82        assert type(B) == ArrayType
83        assert B.typecode() == 'd'
84        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
85
86        A = [1,2,3,4]
87        B = ensure_numeric(A, Float)
88        assert type(B) == ArrayType
89        assert B.typecode() == 'd'
90        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
91
92        A = array([1,2,3,4])
93        B = ensure_numeric(A)
94        assert type(B) == ArrayType
95        assert B.typecode() == 'l'       
96        assert alltrue(A == B)   
97        assert A is B   #Same object
98
99        A = array([1,2,3,4])
100        B = ensure_numeric(A, Float)
101        assert type(B) == ArrayType
102        assert B.typecode() == 'd'       
103        assert alltrue(A == B)   
104        assert A is not B   #Not the same object
105
106        # Check scalars
107        A = 1
108        B = ensure_numeric(A, Float)
109        #print A, B[0], len(B), type(B)
110        #print B.shape
111        assert alltrue(A == B)
112
113        B = ensure_numeric(A, Int)       
114        #print A, B
115        #print B.shape
116        assert alltrue(A == B)
117
118        # Error situation
119
120        B = ensure_numeric('hello', Int)               
121        assert allclose(B, [104, 101, 108, 108, 111])
122
123
124    def test_gradient(self):
125        x0 = 0.0; y0 = 0.0; z0 = 0.0
126        x1 = 1.0; y1 = 0.0; z1 = -1.0
127        x2 = 0.0; y2 = 1.0; z2 = 0.0
128
129        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
130
131        assert zx == -1.0
132        assert zy == 0.0
133
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       
170    def test_gradient2_more(self):
171        """Test two-point gradient more
172        """
173        x0 = 2.0; y0 = 2.0
174        x1 = 8.0; y1 = 3.0
175        x2 = 1.0; y2 = 8.0
176
177        q0 = 2.0
178        q1 = 8.0
179        q2 = q0
180
181        #Gradient of fitted pwl surface
182        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
183        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
184
185        assert a == a_ref
186        assert b == b_ref
187
188
189    def test_machine_precision(self):
190        """test_machine_precision(self):
191        Test the function that calculates epsilon. As this varies on
192        different machines, this is only an indication.
193        """
194
195        eps = get_machine_precision()
196
197        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
198        assert eps > 0.0
199        assert 1.0+eps/2 == 1.0
200       
201       
202    def test_histogram(self):
203        """Test histogram with different bin boundaries
204        """
205       
206        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
207
208        #There are four elements greater than or equal to 3
209        bins = [3]
210        assert allclose(histogram(a, bins), [4])
211
212        bins = [ min(a) ]
213        assert allclose(histogram(a, bins), [len(a)])
214
215        bins = [ max(a)+0.00001 ]
216        assert allclose(histogram(a, bins), [0])       
217       
218        bins = [1,2,3,4]
219        assert allclose(histogram(a, bins), [8,3,3,1])
220
221        bins = [1.1,2,3.1,4]
222        #print histogram(a, bins)
223        assert allclose(histogram(a, bins), [0,6,0,1])
224
225        bins = [0,1.5,2,3]
226        assert allclose(histogram(a, bins), [8,0,3,4])
227        assert allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
228
229        # Check situation with #bins >= #datapoints
230        a = [1.7]
231        bins = [0,1.5,2,3]
232        assert allclose(histogram(a, bins), [0,1,0,0])
233
234        a = [1.7]
235        bins = [0]
236        assert allclose(histogram(a, bins), [1])
237
238        a = [-1.7]
239        bins = [0]
240        assert allclose(histogram(a, bins), [0])
241
242        a = [-1.7]
243        bins = [-1.7]
244        assert allclose(histogram(a, bins), [1])       
245       
246
247    def test_that_C_extension_compiles(self):
248        FN = 'util_ext.c'
249        try:
250            import util_ext
251        except:
252            from compile import compile
253
254            try:
255                compile(FN)
256            except:
257                raise 'Could not compile %s' %FN
258            else:
259                import util_ext
260
261
262    def test_gradient_C_extension(self):
263        from util_ext import gradient as gradient_c
264
265        x0 = 2.0/3; y0 = 2.0/3
266        x1=  8.0/3; y1 = 2.0/3
267        x2 = 2.0/3; y2 = 8.0/3
268
269        q0 = 2.0+2.0/3
270        q1 = 8.0+2.0/3
271        q2 = 2.0+8.0/3
272
273        #Gradient of fitted pwl surface
274        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
275
276        assert abs(a - 3.0) < epsilon
277        assert abs(b - 1.0) < epsilon
278
279
280    def test_gradient_C_extension3(self):
281        from util_ext import gradient as gradient_c
282
283        from RandomArray import uniform, seed
284        seed(17, 53)
285
286        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
287
288        q0 = uniform(0.0, 10.0, 4)
289        q1 = uniform(1.0, 3.0, 4)
290        q2 = uniform(7.0, 20.0, 4)
291
292        for i in range(4):
293            #Gradient of fitted pwl surface
294            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
295                                           q0[i], q1[i], q2[i])
296
297            #print a_ref, b_ref
298            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
299                              q0[i], q1[i], q2[i])
300
301            #print a, a_ref, b, b_ref
302            assert abs(a - a_ref) < epsilon
303            assert abs(b - b_ref) < epsilon
304
305     
306    def test_err(self):
307        x = [2,5] # diff at first position = 4, 4^2 = 16
308        y = [6,7] # diff at secnd position = 2, 2^2 = 4
309        # 16 + 4 = 20
310       
311        # If there is x and y, n=2 and relative=False, this will calc;
312        # sqrt(sum_over_x&y((xi - yi)^2))
313        err__1 = err(x,y,2,False)
314        assert err__1 == sqrt(20)
315        #print "err_", err_
316        #rmsd_1 = err__1*sqrt(1./len(x))
317        #print "err__1*sqrt(1./len(x))", err__1*sqrt(1./len(x))
318        #print "sqrt(10)", sqrt(10)
319       
320        x = [2,7,100]
321        y = [5,10,103]
322        err__2 = err(x,y,2,False)
323        assert err__2 == sqrt(27)
324        #rmsd_2 = err__2*sqrt(1./len(x))
325        #print "err__2*sqrt(1./len(x))", err__2*sqrt(1./len(x))
326
327        x = [2,5,2,7,100]
328        y = [6,7,5,10,103]
329        err_3 = err(x,y,2,False)
330        assert err_3 == sqrt(47)
331       
332        #rmsd_3 = err_3*sqrt(1./len(x))
333        #print "err__3*sqrt(1./len(x))", err__3*sqrt(1./len(x))
334        #print "rmsd_3", rmsd_3
335        #print "sqrt(err_1*err__1+err__2*err__2)/sqrt(5)", \
336        # sqrt(err__1*err__1+err__2*err__2)/sqrt(5)
337        #print "(rmsd_1 + rmsd_2)/2.", (rmsd_1 + rmsd_2)/2.
338        #print "sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.", \
339        #sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.
340       
341    def test_norm(self):
342        x = norm(ensure_numeric([3,4]))
343        assert x == 5.
344
345                                   
346
347#-------------------------------------------------------------
348if __name__ == "__main__":
349    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
350    #suite = unittest.makeSuite(Test_Numerical_Tools,'test_err')
351    runner = unittest.TextTestRunner()
352    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.