source: branches/numpy/anuga/utilities/test_numerical_tools.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 14 years ago

Initial commit of numpy changes. Still a long way to go.
File size: 14.7 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import numpy as num
6
7from math import sqrt, pi
8from anuga.config import epsilon
9from numerical_tools import *
10
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 num.allclose(angle([1.0, 0.0])/pi*180, 0.0)         
26        assert num.allclose(angle([1.0, 1.0])/pi*180, 45.0)
27        assert num.allclose(angle([0.0, 1.0])/pi*180, 90.0)             
28        assert num.allclose(angle([-1.0, 1.0])/pi*180, 135.0)           
29        assert num.allclose(angle([-1.0, 0.0])/pi*180, 180.0)
30        assert num.allclose(angle([-1.0, -1.0])/pi*180, 225.0)
31        assert num.allclose(angle([0.0, -1.0])/pi*180, 270.0)
32        assert num.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 num.allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
40        assert num.allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
41               
42        assert num.allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)       
43        assert num.allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0)     
44       
45        assert num.allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
46        assert num.allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)       
47       
48        assert num.allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)       
49        assert num.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 num.allclose(angle(vc, v_prev)/pi*180, 45.0)
55
56        vc = [ 0.5,  0.0]
57        assert num.allclose(angle(vc, v_prev)/pi*180, 135.0)
58
59        vc = [ -0.5,  0.5]
60        assert num.allclose(angle(vc, v_prev)/pi*180, 270.0)               
61
62
63    def test_anglediff(self):
64        assert num.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 isinstance(B, num.ndarray)
71        assert B.dtype.char == '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 isinstance(B, num.ndarray)
77        assert B.dtype.char == '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, num.float)
82        assert isinstance(B, num.ndarray)
83        assert B.dtype.char == '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, num.float)
88        assert isinstance(B, num.ndarray)
89        assert B.dtype.char == '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 = num.array([1,2,3,4])
93        B = ensure_numeric(A)
94        assert isinstance(B, num.ndarray)
95        assert B.dtype.char == 'l'       
96        assert num.alltrue(A == B)   
97        assert A is B   #Same object
98
99# THIS FAILS!  ASSUMED TYPE IS num.int!?
100#        # check default num.array type, which is supposed to be num.float
101#        A = num.array((1,2,3,4))
102#        assert isinstance(A, num.ndarray)
103#        assert A.dtype.char == 'd', \
104#                "Expected dtype='d', got '%s'" % A.dtype.char
105
106        A = num.array([1,2,3,4])
107        B = ensure_numeric(A, num.float)
108        assert isinstance(B, num.ndarray)
109        assert A.dtype.char == 'l'       
110        assert B.dtype.char == 'd'       
111        assert num.alltrue(A == B)   
112        assert A is not B   # Not the same object
113
114        # Check scalars
115        A = 1
116        B = ensure_numeric(A, num.float)
117        assert num.alltrue(A == B)
118#        print 'A=%s' % str(A)
119#        print 'B=%s, B.shape=%s' % (str(B), str(B.shape))
120
121        B = ensure_numeric(A, num.int)       
122        assert num.alltrue(A == B)
123
124        # try to simulate getting (x,0) shape
125        data_points = [[ 413634. ],]
126        array_data_points = ensure_numeric(data_points)
127        #if not (0,) == array_data_points.shape:
128        #    assert len(array_data_points.shape) == 2
129        #    assert array_data_points.shape[1] == 2
130
131
132    def NO_test_ensure_numeric_char(self):
133        '''numpy can't handle this'''
134
135        # Error situation
136        B = ensure_numeric('hello', num.int)               
137        assert num.allclose(B, [104, 101, 108, 108, 111])
138
139
140    def test_gradient(self):
141        x0 = 0.0; y0 = 0.0; z0 = 0.0
142        x1 = 1.0; y1 = 0.0; z1 = -1.0
143        x2 = 0.0; y2 = 1.0; z2 = 0.0
144
145        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
146
147        assert zx == -1.0
148        assert zy == 0.0
149
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       
186    def test_gradient2_more(self):
187        """Test two-point gradient more
188        """
189        x0 = 2.0; y0 = 2.0
190        x1 = 8.0; y1 = 3.0
191        x2 = 1.0; y2 = 8.0
192
193        q0 = 2.0
194        q1 = 8.0
195        q2 = q0
196
197        #Gradient of fitted pwl surface
198        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
199        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
200
201        assert a == a_ref
202        assert b == b_ref
203
204
205    def test_machine_precision(self):
206        """test_machine_precision(self):
207        Test the function that calculates epsilon. As this varies on
208        different machines, this is only an indication.
209        """
210
211        eps = get_machine_precision()
212
213        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
214        assert eps > 0.0
215        assert 1.0+eps/2 == 1.0
216       
217       
218    def test_histogram(self):
219        """Test histogram with different bin boundaries
220        """
221       
222        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
223
224        #There are four elements greater than or equal to 3
225        bins = [3]
226        assert num.allclose(histogram(a, bins), [4])
227
228        bins = [ min(a) ]
229        assert num.allclose(histogram(a, bins), [len(a)])
230
231        bins = [ max(a)+0.00001 ]
232        assert num.allclose(histogram(a, bins), [0])       
233       
234        bins = [1,2,3,4]
235        assert num.allclose(histogram(a, bins), [8,3,3,1])
236
237        bins = [1.1,2,3.1,4]
238        #print histogram(a, bins)
239        assert num.allclose(histogram(a, bins), [0,6,0,1])
240
241        bins = [0,1.5,2,3]
242        assert num.allclose(histogram(a, bins), [8,0,3,4])
243        assert num.allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
244
245        # Check situation with #bins >= #datapoints
246        a = [1.7]
247        bins = [0,1.5,2,3]
248        assert num.allclose(histogram(a, bins), [0,1,0,0])
249
250        a = [1.7]
251        bins = [0]
252        assert num.allclose(histogram(a, bins), [1])
253
254        a = [-1.7]
255        bins = [0]
256        assert num.allclose(histogram(a, bins), [0])
257
258        a = [-1.7]
259        bins = [-1.7]
260        assert num.allclose(histogram(a, bins), [1])       
261       
262
263    def test_that_C_extension_compiles(self):
264        FN = 'util_ext.c'
265        try:
266            import util_ext
267        except:
268            from compile import compile
269
270            try:
271                compile(FN)
272            except:
273                raise 'Could not compile %s' %FN
274            else:
275                import util_ext
276
277
278    def test_gradient_C_extension(self):
279        from util_ext import gradient as gradient_c
280
281        x0 = 2.0/3; y0 = 2.0/3
282        x1=  8.0/3; y1 = 2.0/3
283        x2 = 2.0/3; y2 = 8.0/3
284
285        q0 = 2.0+2.0/3
286        q1 = 8.0+2.0/3
287        q2 = 2.0+8.0/3
288
289        #Gradient of fitted pwl surface
290        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
291
292        assert abs(a - 3.0) < epsilon
293        assert abs(b - 1.0) < epsilon
294
295
296    def test_gradient_C_extension3(self):
297        from util_ext import gradient as gradient_c
298
299        from RandomArray import uniform, seed
300        seed(17, 53)
301
302        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
303
304        q0 = uniform(0.0, 10.0, 4)
305        q1 = uniform(1.0, 3.0, 4)
306        q2 = uniform(7.0, 20.0, 4)
307
308        for i in range(4):
309            #Gradient of fitted pwl surface
310            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
311                                           q0[i], q1[i], q2[i])
312
313            #print a_ref, b_ref
314            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
315                              q0[i], q1[i], q2[i])
316
317            #print a, a_ref, b, b_ref
318            assert abs(a - a_ref) < epsilon
319            assert abs(b - b_ref) < epsilon
320
321     
322    def test_err(self):
323        x = [2,5] # diff at first position = 4, 4^2 = 16
324        y = [6,7] # diff at secnd position = 2, 2^2 = 4
325        # 16 + 4 = 20
326       
327        # If there is x and y, n=2 and relative=False, this will calc;
328        # sqrt(sum_over_x&y((xi - yi)^2))
329        err__1 = err(x,y,2,False)
330        assert err__1 == sqrt(20)
331        #print "err_", err_
332        #rmsd_1 = err__1*sqrt(1./len(x))
333        #print "err__1*sqrt(1./len(x))", err__1*sqrt(1./len(x))
334        #print "sqrt(10)", sqrt(10)
335       
336        x = [2,7,100]
337        y = [5,10,103]
338        err__2 = err(x,y,2,False)
339        assert err__2 == sqrt(27)
340        #rmsd_2 = err__2*sqrt(1./len(x))
341        #print "err__2*sqrt(1./len(x))", err__2*sqrt(1./len(x))
342
343        x = [2,5,2,7,100]
344        y = [6,7,5,10,103]
345        err_3 = err(x,y,2,False)
346        assert err_3 == sqrt(47)
347       
348        #rmsd_3 = err_3*sqrt(1./len(x))
349        #print "err__3*sqrt(1./len(x))", err__3*sqrt(1./len(x))
350        #print "rmsd_3", rmsd_3
351        #print "sqrt(err_1*err__1+err__2*err__2)/sqrt(5)", \
352        # sqrt(err__1*err__1+err__2*err__2)/sqrt(5)
353        #print "(rmsd_1 + rmsd_2)/2.", (rmsd_1 + rmsd_2)/2.
354        #print "sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.", \
355        #sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.
356       
357    def test_norm(self):
358        x = norm(ensure_numeric([3,4]))
359        assert x == 5.
360
361
362################################################################################
363# Test the is_num_????() functions.
364################################################################################
365
366    def test_is_float(self):
367        def t(val, expected):
368            if expected == True:
369                msg = 'should be num.float?'
370            else:
371                msg = 'should not be num.float?'
372            msg = '%s (%s) %s' % (str(val), type(val), msg)
373            assert is_num_float(val) == expected, msg
374
375        t(1, False)
376        t(1.0, False)
377        t('abc', False)
378        t(None, False)
379        t(num.array(None), False)
380        # can't create array(None, num.int)
381#        t(num.array(None, num.int), False)
382        t(num.array(None, num.float), True)
383        t(num.array(()), True)
384        t(num.array((), num.int), False)
385        t(num.array((), num.float), True)
386        t(num.array((1), num.int), False)
387        t(num.array((1), num.float), True)
388
389        t(num.array((1,2)), False)
390        t(num.array((1,2), num.int), False)
391        t(num.array((1,2), num.float), True)
392        t(num.array([1,2]), False)
393        t(num.array([1,2], num.int), False)
394        t(num.array([1,2], num.float), True)
395
396        t(num.array((1.0,2.0)), True)
397        t(num.array((1.0,2.0), num.int), False)
398        t(num.array((1.0,2.0), num.float), True)
399        t(num.array([1.0,2.0]), True)
400        t(num.array([1.0,2.0], num.int), False)
401        t(num.array([1.0,2.0], num.float), True)
402
403        t(num.array(((1.0,2.0),(3.0,4.0))), True)
404        t(num.array(((1.0,2.0),(3.0,4.0)), num.int), False)
405        t(num.array(((1.0,2.0),(3.0,4.0)), num.float), True)
406        t(num.array([[1.0,2.0],[3.0,4.0]]), True)
407        t(num.array([1.0,2.0], num.int), False)
408        t(num.array([1.0,2.0], num.float), True)
409
410        t(num.array('abc'), False)
411        t(num.array('abc', num.character), False)
412        # can't create array as int from string
413#        t(num.array('abc', num.int), False)
414        # can't create array as float from string
415#        t(num.array('abc', num.float), True)
416
417    def test_is_int(self):
418        def t(val, expected):
419            if expected == True:
420                msg = 'should be num.int?'
421            else:
422                msg = 'should not be num.int?'
423            msg = '%s (%s) %s' % (str(val), type(val), msg)
424            assert is_num_int(val) == expected, msg
425
426        t(1, False)
427        t(1.0, False)
428        t('abc', False)
429        t(None, False)
430        t(num.array(None), False)
431        # can't create array(None, num.int)
432#        t(num.array(None, num.int), True)
433        t(num.array(None, num.float), False)
434        t(num.array((), num.int), True)
435        t(num.array(()), False)
436        t(num.array((), num.float), False)
437        t(num.array((1), num.int), True)
438        t(num.array((1), num.float), False)
439
440        t(num.array((1,2)), True)
441        t(num.array((1,2), num.int), True)
442        t(num.array((1,2), num.float), False)
443        t(num.array([1,2]), True)
444        t(num.array([1,2], num.int), True)
445        t(num.array([1,2], num.float), False)
446
447        t(num.array((1.0,2.0)), False)
448        t(num.array((1.0,2.0), num.int), True)
449        t(num.array((1.0,2.0), num.float), False)
450        t(num.array([1.0,2.0]), False)
451        t(num.array([1.0,2.0], num.int), True)
452        t(num.array([1.0,2.0], num.float), False)
453
454        t(num.array(((1.0,2.0),(3.0,4.0))), False)
455        t(num.array(((1.0,2.0),(3.0,4.0)), num.int), True)
456        t(num.array(((1.0,2.0),(3.0,4.0)), num.float), False)
457        t(num.array([[1.0,2.0],[3.0,4.0]]), False)
458        t(num.array([1.0,2.0], num.int), True)
459        t(num.array([1.0,2.0], num.float), False)
460
461        t(num.array('abc'), False)
462        t(num.array('abc', num.character), False)
463        # can't create array as int from string
464#        t(num.array('abc', num.int), True)
465        # can't create array as float from string
466#        t(num.array('abc', num.float), False)
467
468
469#-------------------------------------------------------------
470if __name__ == "__main__":
471    #suite = unittest.makeSuite(Test_Numerical_Tools,'test')
472    suite = unittest.makeSuite(Test_Numerical_Tools,'test_is_')
473    runner = unittest.TextTestRunner()
474    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.