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

Last change on this file since 6702 was 6456, checked in by rwilson, 16 years ago

More testing for ensure_numeric()
File size: 20.6 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        # check default num.array type, which is supposed to be num.int32
100        A = num.array((1,2,3,4))
101        assert isinstance(A, num.ndarray)
102        msg = "Expected dtype.char='l', got '%s'" % A.dtype.char
103        assert A.dtype.char == 'l', msg
104
105        A = num.array([1,2,3,4])
106        B = ensure_numeric(A, num.float)
107        assert isinstance(B, num.ndarray)
108        assert A.dtype.char == 'l'       
109        assert B.dtype.char == 'd'       
110        assert num.alltrue(A == B)   
111        assert A is not B   # Not the same object
112
113        # Check scalars
114        A = 1
115        B = ensure_numeric(A, num.float)
116        assert num.alltrue(A == B)
117
118        B = ensure_numeric(A, num.int)       
119        assert num.alltrue(A == B)
120
121#        # try to simulate getting (x,0) shape
122#        data_points = [[ 413634. ],]
123#        array_data_points = ensure_numeric(data_points)
124#        if not (0,) == array_data_points.shape:
125#            assert len(array_data_points.shape) == 2
126#            assert array_data_points.shape[1] == 2
127
128        # strings input should raise exception
129        self.failUnlessRaises(Exception, ensure_numeric(['abc',]))
130        self.failUnlessRaises(Exception, ensure_numeric(('abc',)))
131        self.failUnlessRaises(Exception, ensure_numeric(num.array(('abc',))))
132
133    def NO_test_ensure_numeric_char(self):
134        '''numpy can't handle this'''
135
136        # Error situation
137        B = ensure_numeric('hello', num.int)               
138        assert num.allclose(B, [104, 101, 108, 108, 111])
139
140
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
152    def test_gradient_more(self):
153        x0 = 2.0/3; y0 = 2.0/3
154        x1=  8.0/3; y1 = 2.0/3
155        x2 = 2.0/3; y2 = 8.0/3
156
157        q0 = 2.0+2.0/3
158        q1 = 8.0+2.0/3
159        q2 = 2.0+8.0/3
160
161        #Gradient of fitted pwl surface
162        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
163
164        assert abs(a - 3.0) < epsilon
165        assert abs(b - 1.0) < epsilon
166
167
168    def test_gradient2(self):
169        """Test two-point gradient
170        """
171       
172        x0 = 5.0; y0 = 5.0; z0 = 10.0
173        x1 = 8.0; y1 = 2.0; z1 = 1.0
174        x2 = 8.0; y2 = 8.0; z2 = 10.0
175
176        #Reference
177        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
178        a, b = gradient2(x0, y0, x1, y1, z0, z1)
179
180        assert zx == a
181        assert zy == b
182
183        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
184        assert z2_computed == z2
185
186       
187    def test_gradient2_more(self):
188        """Test two-point gradient more
189        """
190        x0 = 2.0; y0 = 2.0
191        x1 = 8.0; y1 = 3.0
192        x2 = 1.0; y2 = 8.0
193
194        q0 = 2.0
195        q1 = 8.0
196        q2 = q0
197
198        #Gradient of fitted pwl surface
199        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
200        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
201
202        assert a == a_ref
203        assert b == b_ref
204
205
206    def test_machine_precision(self):
207        """test_machine_precision(self):
208        Test the function that calculates epsilon. As this varies on
209        different machines, this is only an indication.
210        """
211
212        eps = get_machine_precision()
213
214        assert eps < 1.0e-12, 'Machine precision should be better than 1.0e-12'
215        assert eps > 0.0
216        assert 1.0+eps/2 == 1.0
217       
218       
219    def test_histogram(self):
220        """Test histogram with different bin boundaries
221        """
222       
223        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
224
225        #There are four elements greater than or equal to 3
226        bins = [3]
227        assert num.allclose(histogram(a, bins), [4])
228
229        bins = [ min(a) ]
230        assert num.allclose(histogram(a, bins), [len(a)])
231
232        bins = [ max(a)+0.00001 ]
233        assert num.allclose(histogram(a, bins), [0])       
234       
235        bins = [1,2,3,4]
236        assert num.allclose(histogram(a, bins), [8,3,3,1])
237
238        bins = [1.1,2,3.1,4]
239        #print histogram(a, bins)
240        assert num.allclose(histogram(a, bins), [0,6,0,1])
241
242        bins = [0,1.5,2,3]
243        assert num.allclose(histogram(a, bins), [8,0,3,4])
244        assert num.allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
245
246        # Check situation with #bins >= #datapoints
247        a = [1.7]
248        bins = [0,1.5,2,3]
249        assert num.allclose(histogram(a, bins), [0,1,0,0])
250
251        a = [1.7]
252        bins = [0]
253        assert num.allclose(histogram(a, bins), [1])
254
255        a = [-1.7]
256        bins = [0]
257        assert num.allclose(histogram(a, bins), [0])
258
259        a = [-1.7]
260        bins = [-1.7]
261        assert num.allclose(histogram(a, bins), [1])       
262       
263
264    def test_that_C_extension_compiles(self):
265        FN = 'util_ext.c'
266        try:
267            import util_ext
268        except:
269            from compile import compile
270
271            try:
272                compile(FN)
273            except:
274                raise 'Could not compile %s' %FN
275            else:
276                import util_ext
277
278
279    def test_gradient_C_extension(self):
280        from util_ext import gradient as gradient_c
281
282        x0 = 2.0/3; y0 = 2.0/3
283        x1=  8.0/3; y1 = 2.0/3
284        x2 = 2.0/3; y2 = 8.0/3
285
286        q0 = 2.0+2.0/3
287        q1 = 8.0+2.0/3
288        q2 = 2.0+8.0/3
289
290        #Gradient of fitted pwl surface
291        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
292
293        assert abs(a - 3.0) < epsilon
294        assert abs(b - 1.0) < epsilon
295
296
297    def test_gradient_C_extension3(self):
298        from util_ext import gradient as gradient_c
299
300        from RandomArray import uniform, seed
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        # Make 'points' a numeric array of float coordinates.
614        # Should be changed by ensure_numeric(, num.float96).
615        #####
616        points = num.array([[1.,2.], [3.,4.], [5.,6.]], num.float)
617        points_id = id(points)
618
619        points_new = ensure_numeric(points, num.float96)
620        points_new_id = id(points_new)
621
622        msg = 'ensure_numeric() should return a copy of the input'
623        self.failUnless(points_new_id != points_id, msg)
624
625        # should never change it's input parameter
626        msg = "ensure_numeric() changed it's input parameter"
627        self.failUnless(points_id == id(points), msg)
628
629################################################################################
630
631if __name__ == "__main__":
632    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
633    runner = unittest.TextTestRunner()
634    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.