source: trunk/anuga_core/source/anuga/utilities/test_numerical_tools.py @ 7967

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

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

File size: 20.0 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import numpy as num
6from numpy.random import uniform, seed
7
8from math import sqrt, pi
9from anuga.config import epsilon
10from numerical_tools import *
11
12
13def test_function(x, y):
14    return x+y
15
16class Test_Numerical_Tools(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23    def test_angle1(self):
24        """Test angles between one vector and the x-axis
25        """
26        assert num.allclose(angle([1.0, 0.0])/pi*180, 0.0)         
27        assert num.allclose(angle([1.0, 1.0])/pi*180, 45.0)
28        assert num.allclose(angle([0.0, 1.0])/pi*180, 90.0)             
29        assert num.allclose(angle([-1.0, 1.0])/pi*180, 135.0)           
30        assert num.allclose(angle([-1.0, 0.0])/pi*180, 180.0)
31        assert num.allclose(angle([-1.0, -1.0])/pi*180, 225.0)
32        assert num.allclose(angle([0.0, -1.0])/pi*180, 270.0)
33        assert num.allclose(angle([1.0, -1.0])/pi*180, 315.0)
34               
35                                                         
36    def test_angle2(self):
37        """Test angles between two arbitrary vectors
38        """   
39       
40        assert num.allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
41        assert num.allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
42               
43        assert num.allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)       
44        assert num.allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0)     
45       
46        assert num.allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
47        assert num.allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)       
48       
49        assert num.allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)       
50        assert num.allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
51
52        #From test_get_boundary_polygon_V
53        v_prev = [-0.5, -0.5]
54        vc = [ 0.0,  -0.5]
55        assert num.allclose(angle(vc, v_prev)/pi*180, 45.0)
56
57        vc = [ 0.5,  0.0]
58        assert num.allclose(angle(vc, v_prev)/pi*180, 135.0)
59
60        vc = [ -0.5,  0.5]
61        assert num.allclose(angle(vc, v_prev)/pi*180, 270.0)               
62
63
64    def test_anglediff(self):
65        assert num.allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
66
67       
68    def test_ensure_numeric(self):
69        A = [1,2,3,4]
70        B = ensure_numeric(A)
71        assert isinstance(B, num.ndarray)
72        assert B.dtype.char == 'l'
73        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
74
75        A = [1,2,3.14,4]
76        B = ensure_numeric(A)
77        assert isinstance(B, num.ndarray)
78        assert B.dtype.char == 'd'
79        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
80
81        A = [1,2,3,4]
82        B = ensure_numeric(A, num.float)
83        assert isinstance(B, num.ndarray)
84        assert B.dtype.char == 'd'
85        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
86
87        A = [1,2,3,4]
88        B = ensure_numeric(A, num.float)
89        assert isinstance(B, num.ndarray)
90        assert B.dtype.char == 'd'
91        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
92
93        A = num.array([1,2,3,4])
94        B = ensure_numeric(A)
95        assert isinstance(B, num.ndarray)
96        assert B.dtype.char == 'l'       
97        assert num.alltrue(A == B)   
98        assert A is B   #Same object
99
100        # check default num.array type, which is supposed to be num.int32
101        A = num.array((1,2,3,4))
102        assert isinstance(A, num.ndarray)
103        msg = "Expected dtype.char='l', got '%s'" % A.dtype.char
104        assert A.dtype.char == 'l', msg
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
119        B = ensure_numeric(A, num.int)       
120        assert num.alltrue(A == B)
121
122#        # try to simulate getting (x,0) shape
123#        data_points = [[ 413634. ],]
124#        array_data_points = ensure_numeric(data_points)
125#        if not (0,) == array_data_points.shape:
126#            assert len(array_data_points.shape) == 2
127#            assert array_data_points.shape[1] == 2
128
129        # strings input should raise exception
130        self.failUnlessRaises(Exception, ensure_numeric(['abc',]))
131        self.failUnlessRaises(Exception, ensure_numeric(('abc',)))
132        self.failUnlessRaises(Exception, ensure_numeric(num.array(('abc',))))
133
134    def NO_test_ensure_numeric_char(self):
135        '''numpy can't handle this'''
136
137        # Error situation
138        B = ensure_numeric('hello', num.int)               
139        assert num.allclose(B, [104, 101, 108, 108, 111])
140
141
142    def test_gradient(self):
143        x0 = 0.0; y0 = 0.0; z0 = 0.0
144        x1 = 1.0; y1 = 0.0; z1 = -1.0
145        x2 = 0.0; y2 = 1.0; z2 = 0.0
146
147        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
148
149        assert zx == -1.0
150        assert zy == 0.0
151
152
153    def test_gradient_more(self):
154        x0 = 2.0/3; y0 = 2.0/3
155        x1=  8.0/3; y1 = 2.0/3
156        x2 = 2.0/3; y2 = 8.0/3
157
158        q0 = 2.0+2.0/3
159        q1 = 8.0+2.0/3
160        q2 = 2.0+8.0/3
161
162        #Gradient of fitted pwl surface
163        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
164
165        assert abs(a - 3.0) < epsilon
166        assert abs(b - 1.0) < epsilon
167
168
169    def test_gradient2(self):
170        """Test two-point gradient
171        """
172       
173        x0 = 5.0; y0 = 5.0; z0 = 10.0
174        x1 = 8.0; y1 = 2.0; z1 = 1.0
175        x2 = 8.0; y2 = 8.0; z2 = 10.0
176
177        #Reference
178        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
179        a, b = gradient2(x0, y0, x1, y1, z0, z1)
180
181        assert zx == a
182        assert zy == b
183
184        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
185        assert z2_computed == z2
186
187       
188    def test_gradient2_more(self):
189        """Test two-point gradient more
190        """
191        x0 = 2.0; y0 = 2.0
192        x1 = 8.0; y1 = 3.0
193        x2 = 1.0; y2 = 8.0
194
195        q0 = 2.0
196        q1 = 8.0
197        q2 = q0
198
199        #Gradient of fitted pwl surface
200        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
201        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
202
203        assert a == a_ref
204        assert b == b_ref
205
206
207    def test_machine_precision(self):
208        """test_machine_precision(self):
209        Test the function that calculates epsilon. As this varies on
210        different machines, this is only an indication.
211        """
212
213        eps = get_machine_precision()
214
215        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
216        assert eps > 0.0
217        assert 1.0+eps/2 == 1.0
218       
219       
220    def test_histogram(self):
221        """Test histogram with different bin boundaries
222        """
223       
224        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
225
226        #There are four elements greater than or equal to 3
227        bins = [3]
228        assert num.allclose(histogram(a, bins), [4])
229
230        bins = [ min(a) ]
231        assert num.allclose(histogram(a, bins), [len(a)])
232
233        bins = [ max(a)+0.00001 ]
234        assert num.allclose(histogram(a, bins), [0])       
235       
236        bins = [1,2,3,4]
237        assert num.allclose(histogram(a, bins), [8,3,3,1])
238
239        bins = [1.1,2,3.1,4]
240        #print histogram(a, bins)
241        assert num.allclose(histogram(a, bins), [0,6,0,1])
242
243        bins = [0,1.5,2,3]
244        assert num.allclose(histogram(a, bins), [8,0,3,4])
245        assert num.allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
246
247        # Check situation with #bins >= #datapoints
248        a = [1.7]
249        bins = [0,1.5,2,3]
250        assert num.allclose(histogram(a, bins), [0,1,0,0])
251
252        a = [1.7]
253        bins = [0]
254        assert num.allclose(histogram(a, bins), [1])
255
256        a = [-1.7]
257        bins = [0]
258        assert num.allclose(histogram(a, bins), [0])
259
260        a = [-1.7]
261        bins = [-1.7]
262        assert num.allclose(histogram(a, bins), [1])       
263       
264
265    def test_that_C_extension_compiles(self):
266        FN = 'util_ext.c'
267        try:
268            import util_ext
269        except:
270            from compile import compile
271
272            try:
273                compile(FN)
274            except:
275                raise 'Could not compile %s' %FN
276            else:
277                import util_ext
278
279
280    def test_gradient_C_extension(self):
281        from util_ext import gradient as gradient_c
282
283        x0 = 2.0/3; y0 = 2.0/3
284        x1=  8.0/3; y1 = 2.0/3
285        x2 = 2.0/3; y2 = 8.0/3
286
287        q0 = 2.0+2.0/3
288        q1 = 8.0+2.0/3
289        q2 = 2.0+8.0/3
290
291        #Gradient of fitted pwl surface
292        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
293
294        assert abs(a - 3.0) < epsilon
295        assert abs(b - 1.0) < epsilon
296
297
298    def test_gradient_C_extension3(self):
299        from util_ext import gradient as gradient_c
300
301        seed((17, 53))
302
303        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
304
305        q0 = uniform(0.0, 10.0, 4)
306        q1 = uniform(1.0, 3.0, 4)
307        q2 = uniform(7.0, 20.0, 4)
308
309        for i in range(4):
310            #Gradient of fitted pwl surface
311            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
312                                           q0[i], q1[i], q2[i])
313
314            #print a_ref, b_ref
315            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
316                              q0[i], q1[i], q2[i])
317
318            #print a, a_ref, b, b_ref
319            assert abs(a - a_ref) < epsilon
320            assert abs(b - b_ref) < epsilon
321
322     
323    def test_err(self):
324        x = [2,5] # diff at first position = 4, 4^2 = 16
325        y = [6,7] # diff at secnd position = 2, 2^2 = 4
326        # 16 + 4 = 20
327       
328        # If there is x and y, n=2 and relative=False, this will calc;
329        # sqrt(sum_over_x&y((xi - yi)^2))
330        err__1 = err(x,y,2,False)
331        assert err__1 == sqrt(20)
332        #print "err_", err_
333        #rmsd_1 = err__1*sqrt(1./len(x))
334        #print "err__1*sqrt(1./len(x))", err__1*sqrt(1./len(x))
335        #print "sqrt(10)", sqrt(10)
336       
337        x = [2,7,100]
338        y = [5,10,103]
339        err__2 = err(x,y,2,False)
340        assert err__2 == sqrt(27)
341        #rmsd_2 = err__2*sqrt(1./len(x))
342        #print "err__2*sqrt(1./len(x))", err__2*sqrt(1./len(x))
343
344        x = [2,5,2,7,100]
345        y = [6,7,5,10,103]
346        err_3 = err(x,y,2,False)
347        assert err_3 == sqrt(47)
348       
349        #rmsd_3 = err_3*sqrt(1./len(x))
350        #print "err__3*sqrt(1./len(x))", err__3*sqrt(1./len(x))
351        #print "rmsd_3", rmsd_3
352        #print "sqrt(err_1*err__1+err__2*err__2)/sqrt(5)", \
353        # sqrt(err__1*err__1+err__2*err__2)/sqrt(5)
354        #print "(rmsd_1 + rmsd_2)/2.", (rmsd_1 + rmsd_2)/2.
355        #print "sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.", \
356        #sqrt((rmsd_1*rmsd_1 + rmsd_2*rmsd_2))/2.
357       
358    def test_norm(self):
359        x = norm(ensure_numeric([3,4]))
360        assert x == 5.
361
362
363################################################################################
364# Test the is_num_????() functions.
365################################################################################
366
367    def test_is_float(self):
368        def t(val, expected):
369            if expected == True:
370                msg = 'should be num.float?'
371            else:
372                msg = 'should not be num.float?'
373            msg = '%s (%s) %s' % (str(val), type(val), msg)
374            assert is_num_float(val) == expected, msg
375
376        t(1, False)
377        t(1.0, False)
378        t('abc', False)
379        t(None, False)
380        t(num.array(None), False)
381        # can't create array(None, num.int)
382#        t(num.array(None, num.int), False)
383        t(num.array(None, num.float), True)
384        t(num.array(()), True)
385        t(num.array((), num.int), False)
386        t(num.array((), num.float), True)
387        t(num.array((1), num.int), False)
388        t(num.array((1), num.float), True)
389
390        t(num.array((1,2)), False)
391        t(num.array((1,2), num.int), False)
392        t(num.array((1,2), num.float), True)
393        t(num.array([1,2]), False)
394        t(num.array([1,2], num.int), False)
395        t(num.array([1,2], num.float), True)
396
397        t(num.array((1.0,2.0)), True)
398        t(num.array((1.0,2.0), num.int), False)
399        t(num.array((1.0,2.0), num.float), True)
400        t(num.array([1.0,2.0]), True)
401        t(num.array([1.0,2.0], num.int), False)
402        t(num.array([1.0,2.0], num.float), True)
403
404        t(num.array(((1.0,2.0),(3.0,4.0))), True)
405        t(num.array(((1.0,2.0),(3.0,4.0)), num.int), False)
406        t(num.array(((1.0,2.0),(3.0,4.0)), num.float), True)
407        t(num.array([[1.0,2.0],[3.0,4.0]]), True)
408        t(num.array([1.0,2.0], num.int), False)
409        t(num.array([1.0,2.0], num.float), True)
410
411        t(num.array('abc'), False)
412        t(num.array('abc', num.character), False)
413        # can't create array as int from string
414#        t(num.array('abc', num.int), False)
415        # can't create array as float from string
416#        t(num.array('abc', num.float), True)
417
418    def test_is_int(self):
419        def t(val, expected):
420            if expected == True:
421                msg = 'should be num.int?'
422            else:
423                msg = 'should not be num.int?'
424            msg = '%s (%s) %s' % (str(val), type(val), msg)
425            assert is_num_int(val) == expected, msg
426
427        t(1, False)
428        t(1.0, False)
429        t('abc', False)
430        t(None, False)
431        t(num.array(None), False)
432        # can't create array(None, num.int)
433#        t(num.array(None, num.int), True)
434        t(num.array(None, num.float), False)
435        t(num.array((), num.int), True)
436        t(num.array(()), False)
437        t(num.array((), num.float), False)
438        t(num.array((1), num.int), True)
439        t(num.array((1), num.float), False)
440
441        t(num.array((1,2)), True)
442        t(num.array((1,2), num.int), True)
443        t(num.array((1,2), num.float), False)
444        t(num.array([1,2]), True)
445        t(num.array([1,2], num.int), True)
446        t(num.array([1,2], num.float), False)
447
448        t(num.array((1.0,2.0)), False)
449        t(num.array((1.0,2.0), num.int), True)
450        t(num.array((1.0,2.0), num.float), False)
451        t(num.array([1.0,2.0]), False)
452        t(num.array([1.0,2.0], num.int), True)
453        t(num.array([1.0,2.0], num.float), False)
454
455        t(num.array(((1.0,2.0),(3.0,4.0))), False)
456        t(num.array(((1.0,2.0),(3.0,4.0)), num.int), True)
457        t(num.array(((1.0,2.0),(3.0,4.0)), num.float), False)
458        t(num.array([[1.0,2.0],[3.0,4.0]]), False)
459        t(num.array([1.0,2.0], num.int), True)
460        t(num.array([1.0,2.0], num.float), False)
461
462        t(num.array('abc'), False)
463        t(num.array('abc', num.character), False)
464        # can't create array as int from string
465#        t(num.array('abc', num.int), True)
466        # can't create array as float from string
467#        t(num.array('abc', num.float), False)
468
469    ##
470    # @brief Test to see if ensure_numeric() behaves as we expect.
471    # @note Under Numeric ensure_numeric() *always* returned a copy (bug).
472    #       Under numpy it copies only when it has to.
473    def test_ensure_numeric_copy(self):
474        #####
475        # Make 'points' a _list_ of coordinates.
476        # Should be changed by ensure_numeric().
477        #####
478        points = [[1.,2.], [3.,4.], [5.,6.]]
479        points_id = id(points)
480
481        points_new = ensure_numeric(points, num.float)
482        points_new_id = id(points_new)
483
484        msg = 'ensure_numeric() should return a copy of a list'
485        self.failUnless(points_new_id != points_id, msg)
486
487        # should never change it's input parameter
488        msg = "ensure_numeric() changed it's input parameter"
489        self.failUnless(points_id == id(points), msg)
490
491        #####
492        # Make 'points' a _tuple_ of coordinates.
493        # Should be changed by ensure_numeric().
494        #####
495        points = ((1.,2.), (3.,4.), (5.,6.))
496        points_id = id(points)
497
498        points_new = ensure_numeric(points, num.int)
499        points_new_id = id(points_new)
500
501        msg = 'ensure_numeric() should return a copy of a list'
502        self.failUnless(points_new_id != points_id, msg)
503
504        # should never change it's input parameter
505        msg = "ensure_numeric() changed it's input parameter"
506        self.failUnless(points_id == id(points), msg)
507
508        #####
509        # Make 'points' a numeric array of float coordinates.
510        # Should NOT be changed by ensure_numeric().
511        #####
512        points = num.array([[1.,2.], [3.,4.], [5.,6.]], num.float)
513        points_id = id(points)
514
515        points_new = ensure_numeric(points, num.float)
516        points_new_id = id(points_new)
517
518        msg = 'ensure_numeric() should return the original input'
519        self.failUnless(points_new_id == points_id, msg)
520
521        # should never change it's input parameter
522        msg = "ensure_numeric() changed it's input parameter"
523        self.failUnless(points_id == id(points), msg)
524
525        #####
526        # Make 'points' a numeric array of int coordinates.
527        # Should be changed by ensure_numeric(, num.float).
528        #####
529        points = num.array([[1,2], [3,4], [5,6]], num.int)
530        points_id = id(points)
531
532        points_new = ensure_numeric(points, num.float)
533        points_new_id = id(points_new)
534
535        msg = 'ensure_numeric() should return a copy of the input'
536        self.failUnless(points_new_id != points_id, msg)
537
538        # should never change it's input parameter
539        msg = "ensure_numeric() changed it's input parameter"
540        self.failUnless(points_id == id(points), msg)
541
542        #####
543        # Make 'points' a numeric array of int coordinates.
544        # Should NOT be changed by ensure_numeric(, num.int).
545        #####
546        points = num.array([[1,2], [3,4], [5,6]], num.int)
547        points_id = id(points)
548
549        points_new = ensure_numeric(points, num.int)
550        points_new_id = id(points_new)
551
552        msg = 'ensure_numeric() should return the original input'
553        self.failUnless(points_new_id == points_id, msg)
554
555        # should never change it's input parameter
556        msg = "ensure_numeric() changed it's input parameter"
557        self.failUnless(points_id == id(points), msg)
558
559        #####
560        # Make 'points' a numeric array of float32 coordinates.
561        # Should NOT be changed by ensure_numeric(, num.float32).
562        #####
563        points = num.array([[1.,2.], [3.,4.], [5.,6.]], num.float32)
564        points_id = id(points)
565
566        points_new = ensure_numeric(points, num.float32)
567        points_new_id = id(points_new)
568
569        msg = 'ensure_numeric() should return the original input'
570        self.failUnless(points_new_id == points_id, msg)
571
572        # should never change it's input parameter
573        msg = "ensure_numeric() changed it's input parameter"
574        self.failUnless(points_id == id(points), msg)
575
576        #####
577        # Make 'points' a numeric array of float32 coordinates.
578        # Should be changed by ensure_numeric(, num.float64).
579        #####
580        points = num.array([[1.,2.], [3.,4.], [5.,6.]], num.float32)
581        points_id = id(points)
582
583        points_new = ensure_numeric(points, num.float64)
584        points_new_id = id(points_new)
585
586        msg = 'ensure_numeric() should return a copy of the input'
587        self.failUnless(points_new_id != points_id, msg)
588
589        # should never change it's input parameter
590        msg = "ensure_numeric() changed it's input parameter"
591        self.failUnless(points_id == id(points), msg)
592
593        #####
594        # Make 'points' a numeric array of float coordinates.
595        # Should NOT be changed by ensure_numeric(, num.float64).
596        #####
597        points = num.array([[1.,2.], [3.,4.], [5.,6.]], num.float)
598        points_id = id(points)
599
600        points_new = ensure_numeric(points, num.float64)
601        points_new_id = id(points_new)
602
603        msg = 'ensure_numeric() should return the original input'
604        self.failUnless(points_new_id == points_id, msg)
605        #msg = 'ensure_numeric() should return a copy of the input'
606        #self.failUnless(points_new_id != points_id, msg)
607
608        # should never change it's input parameter
609        msg = "ensure_numeric() changed it's input parameter"
610        self.failUnless(points_id == id(points), msg)
611
612################################################################################
613
614if __name__ == "__main__":
615    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
616    runner = unittest.TextTestRunner()
617    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.