source: inundation/utilities/test_numerical_tools.py @ 2709

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

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

File size: 7.4 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7from pyvolution.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_histogram(self):
189        """Test histogram with different bin boundaries
190        """
191       
192        a = [1,1,1,1,1,2,1,3,2,3,1,2,3,4,1]
193
194
195        #There are four elements greater than or equal to 3
196        bins = [3]
197        assert allclose(histogram(a, bins), [4])
198
199
200        bins = [ min(a) ]
201        assert allclose(histogram(a, bins), [len(a)])
202
203
204        bins = [ max(a)+0.00001 ]
205        assert allclose(histogram(a, bins), [0])       
206
207       
208        bins = [1,2,3,4]
209        assert allclose(histogram(a, bins), [8,3,3,1])
210
211
212        bins = [1.1,2,3.1,4]
213        #print histogram(a, bins)
214        assert allclose(histogram(a, bins), [0,6,0,1])
215
216
217        bins = [0,1.5,2,3]
218        assert allclose(histogram(a, bins), [8,0,3,4])
219
220
221        assert allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
222
223
224       
225
226    def test_that_C_extension_compiles(self):
227        FN = 'util_ext.c'
228        try:
229            import util_ext
230        except:
231            from compile import compile
232
233            try:
234                compile(FN)
235            except:
236                raise 'Could not compile %s' %FN
237            else:
238                import util_ext
239
240
241    def test_gradient_C_extension(self):
242        from util_ext import gradient as gradient_c
243
244        x0 = 2.0/3; y0 = 2.0/3
245        x1=  8.0/3; y1 = 2.0/3
246        x2 = 2.0/3; y2 = 8.0/3
247
248        q0 = 2.0+2.0/3
249        q1 = 8.0+2.0/3
250        q2 = 2.0+8.0/3
251
252        #Gradient of fitted pwl surface
253        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
254
255        assert abs(a - 3.0) < epsilon
256        assert abs(b - 1.0) < epsilon
257
258
259    def test_gradient_C_extension3(self):
260        from util_ext import gradient as gradient_c
261
262        from RandomArray import uniform, seed
263        seed(17, 53)
264
265        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
266
267        q0 = uniform(0.0, 10.0, 4)
268        q1 = uniform(1.0, 3.0, 4)
269        q2 = uniform(7.0, 20.0, 4)
270
271
272        for i in range(4):
273            #Gradient of fitted pwl surface
274            a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2,
275                                           q0[i], q1[i], q2[i])
276
277            #print a_ref, b_ref
278            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
279                              q0[i], q1[i], q2[i])
280
281            #print a, a_ref, b, b_ref
282            assert abs(a - a_ref) < epsilon
283            assert abs(b - b_ref) < epsilon
284
285       
286
287
288#-------------------------------------------------------------
289if __name__ == "__main__":
290    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
291    runner = unittest.TextTestRunner()
292    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.