comparison DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/polynomial/tests/test_chebyshev.py @ 87:2a2c65a20a8b

Add Python libs and headers
author Chris Cannam
date Wed, 25 Feb 2015 14:05:22 +0000
parents
children
comparison
equal deleted inserted replaced
86:413a9d26189e 87:2a2c65a20a8b
1 """Tests for chebyshev module.
2
3 """
4 from __future__ import division, absolute_import, print_function
5
6 import numpy as np
7 import numpy.polynomial.chebyshev as cheb
8 from numpy.polynomial.polynomial import polyval
9 from numpy.testing import (
10 TestCase, assert_almost_equal, assert_raises,
11 assert_equal, assert_, run_module_suite)
12
13
14 def trim(x):
15 return cheb.chebtrim(x, tol=1e-6)
16
17 T0 = [1]
18 T1 = [0, 1]
19 T2 = [-1, 0, 2]
20 T3 = [0, -3, 0, 4]
21 T4 = [1, 0, -8, 0, 8]
22 T5 = [0, 5, 0, -20, 0, 16]
23 T6 = [-1, 0, 18, 0, -48, 0, 32]
24 T7 = [0, -7, 0, 56, 0, -112, 0, 64]
25 T8 = [1, 0, -32, 0, 160, 0, -256, 0, 128]
26 T9 = [0, 9, 0, -120, 0, 432, 0, -576, 0, 256]
27
28 Tlist = [T0, T1, T2, T3, T4, T5, T6, T7, T8, T9]
29
30
31 class TestPrivate(TestCase):
32
33 def test__cseries_to_zseries(self):
34 for i in range(5):
35 inp = np.array([2] + [1]*i, np.double)
36 tgt = np.array([.5]*i + [2] + [.5]*i, np.double)
37 res = cheb._cseries_to_zseries(inp)
38 assert_equal(res, tgt)
39
40 def test__zseries_to_cseries(self):
41 for i in range(5):
42 inp = np.array([.5]*i + [2] + [.5]*i, np.double)
43 tgt = np.array([2] + [1]*i, np.double)
44 res = cheb._zseries_to_cseries(inp)
45 assert_equal(res, tgt)
46
47
48 class TestConstants(TestCase):
49
50 def test_chebdomain(self):
51 assert_equal(cheb.chebdomain, [-1, 1])
52
53 def test_chebzero(self):
54 assert_equal(cheb.chebzero, [0])
55
56 def test_chebone(self):
57 assert_equal(cheb.chebone, [1])
58
59 def test_chebx(self):
60 assert_equal(cheb.chebx, [0, 1])
61
62
63 class TestArithmetic(TestCase):
64
65 def test_chebadd(self):
66 for i in range(5):
67 for j in range(5):
68 msg = "At i=%d, j=%d" % (i, j)
69 tgt = np.zeros(max(i, j) + 1)
70 tgt[i] += 1
71 tgt[j] += 1
72 res = cheb.chebadd([0]*i + [1], [0]*j + [1])
73 assert_equal(trim(res), trim(tgt), err_msg=msg)
74
75 def test_chebsub(self):
76 for i in range(5):
77 for j in range(5):
78 msg = "At i=%d, j=%d" % (i, j)
79 tgt = np.zeros(max(i, j) + 1)
80 tgt[i] += 1
81 tgt[j] -= 1
82 res = cheb.chebsub([0]*i + [1], [0]*j + [1])
83 assert_equal(trim(res), trim(tgt), err_msg=msg)
84
85 def test_chebmulx(self):
86 assert_equal(cheb.chebmulx([0]), [0])
87 assert_equal(cheb.chebmulx([1]), [0, 1])
88 for i in range(1, 5):
89 ser = [0]*i + [1]
90 tgt = [0]*(i - 1) + [.5, 0, .5]
91 assert_equal(cheb.chebmulx(ser), tgt)
92
93 def test_chebmul(self):
94 for i in range(5):
95 for j in range(5):
96 msg = "At i=%d, j=%d" % (i, j)
97 tgt = np.zeros(i + j + 1)
98 tgt[i + j] += .5
99 tgt[abs(i - j)] += .5
100 res = cheb.chebmul([0]*i + [1], [0]*j + [1])
101 assert_equal(trim(res), trim(tgt), err_msg=msg)
102
103 def test_chebdiv(self):
104 for i in range(5):
105 for j in range(5):
106 msg = "At i=%d, j=%d" % (i, j)
107 ci = [0]*i + [1]
108 cj = [0]*j + [1]
109 tgt = cheb.chebadd(ci, cj)
110 quo, rem = cheb.chebdiv(tgt, ci)
111 res = cheb.chebadd(cheb.chebmul(quo, ci), rem)
112 assert_equal(trim(res), trim(tgt), err_msg=msg)
113
114
115 class TestEvaluation(TestCase):
116 # coefficients of 1 + 2*x + 3*x**2
117 c1d = np.array([2.5, 2., 1.5])
118 c2d = np.einsum('i,j->ij', c1d, c1d)
119 c3d = np.einsum('i,j,k->ijk', c1d, c1d, c1d)
120
121 # some random values in [-1, 1)
122 x = np.random.random((3, 5))*2 - 1
123 y = polyval(x, [1., 2., 3.])
124
125 def test_chebval(self):
126 #check empty input
127 assert_equal(cheb.chebval([], [1]).size, 0)
128
129 #check normal input)
130 x = np.linspace(-1, 1)
131 y = [polyval(x, c) for c in Tlist]
132 for i in range(10):
133 msg = "At i=%d" % i
134 tgt = y[i]
135 res = cheb.chebval(x, [0]*i + [1])
136 assert_almost_equal(res, tgt, err_msg=msg)
137
138 #check that shape is preserved
139 for i in range(3):
140 dims = [2]*i
141 x = np.zeros(dims)
142 assert_equal(cheb.chebval(x, [1]).shape, dims)
143 assert_equal(cheb.chebval(x, [1, 0]).shape, dims)
144 assert_equal(cheb.chebval(x, [1, 0, 0]).shape, dims)
145
146 def test_chebval2d(self):
147 x1, x2, x3 = self.x
148 y1, y2, y3 = self.y
149
150 #test exceptions
151 assert_raises(ValueError, cheb.chebval2d, x1, x2[:2], self.c2d)
152
153 #test values
154 tgt = y1*y2
155 res = cheb.chebval2d(x1, x2, self.c2d)
156 assert_almost_equal(res, tgt)
157
158 #test shape
159 z = np.ones((2, 3))
160 res = cheb.chebval2d(z, z, self.c2d)
161 assert_(res.shape == (2, 3))
162
163 def test_chebval3d(self):
164 x1, x2, x3 = self.x
165 y1, y2, y3 = self.y
166
167 #test exceptions
168 assert_raises(ValueError, cheb.chebval3d, x1, x2, x3[:2], self.c3d)
169
170 #test values
171 tgt = y1*y2*y3
172 res = cheb.chebval3d(x1, x2, x3, self.c3d)
173 assert_almost_equal(res, tgt)
174
175 #test shape
176 z = np.ones((2, 3))
177 res = cheb.chebval3d(z, z, z, self.c3d)
178 assert_(res.shape == (2, 3))
179
180 def test_chebgrid2d(self):
181 x1, x2, x3 = self.x
182 y1, y2, y3 = self.y
183
184 #test values
185 tgt = np.einsum('i,j->ij', y1, y2)
186 res = cheb.chebgrid2d(x1, x2, self.c2d)
187 assert_almost_equal(res, tgt)
188
189 #test shape
190 z = np.ones((2, 3))
191 res = cheb.chebgrid2d(z, z, self.c2d)
192 assert_(res.shape == (2, 3)*2)
193
194 def test_chebgrid3d(self):
195 x1, x2, x3 = self.x
196 y1, y2, y3 = self.y
197
198 #test values
199 tgt = np.einsum('i,j,k->ijk', y1, y2, y3)
200 res = cheb.chebgrid3d(x1, x2, x3, self.c3d)
201 assert_almost_equal(res, tgt)
202
203 #test shape
204 z = np.ones((2, 3))
205 res = cheb.chebgrid3d(z, z, z, self.c3d)
206 assert_(res.shape == (2, 3)*3)
207
208
209 class TestIntegral(TestCase):
210
211 def test_chebint(self):
212 # check exceptions
213 assert_raises(ValueError, cheb.chebint, [0], .5)
214 assert_raises(ValueError, cheb.chebint, [0], -1)
215 assert_raises(ValueError, cheb.chebint, [0], 1, [0, 0])
216
217 # test integration of zero polynomial
218 for i in range(2, 5):
219 k = [0]*(i - 2) + [1]
220 res = cheb.chebint([0], m=i, k=k)
221 assert_almost_equal(res, [0, 1])
222
223 # check single integration with integration constant
224 for i in range(5):
225 scl = i + 1
226 pol = [0]*i + [1]
227 tgt = [i] + [0]*i + [1/scl]
228 chebpol = cheb.poly2cheb(pol)
229 chebint = cheb.chebint(chebpol, m=1, k=[i])
230 res = cheb.cheb2poly(chebint)
231 assert_almost_equal(trim(res), trim(tgt))
232
233 # check single integration with integration constant and lbnd
234 for i in range(5):
235 scl = i + 1
236 pol = [0]*i + [1]
237 chebpol = cheb.poly2cheb(pol)
238 chebint = cheb.chebint(chebpol, m=1, k=[i], lbnd=-1)
239 assert_almost_equal(cheb.chebval(-1, chebint), i)
240
241 # check single integration with integration constant and scaling
242 for i in range(5):
243 scl = i + 1
244 pol = [0]*i + [1]
245 tgt = [i] + [0]*i + [2/scl]
246 chebpol = cheb.poly2cheb(pol)
247 chebint = cheb.chebint(chebpol, m=1, k=[i], scl=2)
248 res = cheb.cheb2poly(chebint)
249 assert_almost_equal(trim(res), trim(tgt))
250
251 # check multiple integrations with default k
252 for i in range(5):
253 for j in range(2, 5):
254 pol = [0]*i + [1]
255 tgt = pol[:]
256 for k in range(j):
257 tgt = cheb.chebint(tgt, m=1)
258 res = cheb.chebint(pol, m=j)
259 assert_almost_equal(trim(res), trim(tgt))
260
261 # check multiple integrations with defined k
262 for i in range(5):
263 for j in range(2, 5):
264 pol = [0]*i + [1]
265 tgt = pol[:]
266 for k in range(j):
267 tgt = cheb.chebint(tgt, m=1, k=[k])
268 res = cheb.chebint(pol, m=j, k=list(range(j)))
269 assert_almost_equal(trim(res), trim(tgt))
270
271 # check multiple integrations with lbnd
272 for i in range(5):
273 for j in range(2, 5):
274 pol = [0]*i + [1]
275 tgt = pol[:]
276 for k in range(j):
277 tgt = cheb.chebint(tgt, m=1, k=[k], lbnd=-1)
278 res = cheb.chebint(pol, m=j, k=list(range(j)), lbnd=-1)
279 assert_almost_equal(trim(res), trim(tgt))
280
281 # check multiple integrations with scaling
282 for i in range(5):
283 for j in range(2, 5):
284 pol = [0]*i + [1]
285 tgt = pol[:]
286 for k in range(j):
287 tgt = cheb.chebint(tgt, m=1, k=[k], scl=2)
288 res = cheb.chebint(pol, m=j, k=list(range(j)), scl=2)
289 assert_almost_equal(trim(res), trim(tgt))
290
291 def test_chebint_axis(self):
292 # check that axis keyword works
293 c2d = np.random.random((3, 4))
294
295 tgt = np.vstack([cheb.chebint(c) for c in c2d.T]).T
296 res = cheb.chebint(c2d, axis=0)
297 assert_almost_equal(res, tgt)
298
299 tgt = np.vstack([cheb.chebint(c) for c in c2d])
300 res = cheb.chebint(c2d, axis=1)
301 assert_almost_equal(res, tgt)
302
303 tgt = np.vstack([cheb.chebint(c, k=3) for c in c2d])
304 res = cheb.chebint(c2d, k=3, axis=1)
305 assert_almost_equal(res, tgt)
306
307
308 class TestDerivative(TestCase):
309
310 def test_chebder(self):
311 # check exceptions
312 assert_raises(ValueError, cheb.chebder, [0], .5)
313 assert_raises(ValueError, cheb.chebder, [0], -1)
314
315 # check that zeroth deriviative does nothing
316 for i in range(5):
317 tgt = [0]*i + [1]
318 res = cheb.chebder(tgt, m=0)
319 assert_equal(trim(res), trim(tgt))
320
321 # check that derivation is the inverse of integration
322 for i in range(5):
323 for j in range(2, 5):
324 tgt = [0]*i + [1]
325 res = cheb.chebder(cheb.chebint(tgt, m=j), m=j)
326 assert_almost_equal(trim(res), trim(tgt))
327
328 # check derivation with scaling
329 for i in range(5):
330 for j in range(2, 5):
331 tgt = [0]*i + [1]
332 res = cheb.chebder(cheb.chebint(tgt, m=j, scl=2), m=j, scl=.5)
333 assert_almost_equal(trim(res), trim(tgt))
334
335 def test_chebder_axis(self):
336 # check that axis keyword works
337 c2d = np.random.random((3, 4))
338
339 tgt = np.vstack([cheb.chebder(c) for c in c2d.T]).T
340 res = cheb.chebder(c2d, axis=0)
341 assert_almost_equal(res, tgt)
342
343 tgt = np.vstack([cheb.chebder(c) for c in c2d])
344 res = cheb.chebder(c2d, axis=1)
345 assert_almost_equal(res, tgt)
346
347
348 class TestVander(TestCase):
349 # some random values in [-1, 1)
350 x = np.random.random((3, 5))*2 - 1
351
352 def test_chebvander(self):
353 # check for 1d x
354 x = np.arange(3)
355 v = cheb.chebvander(x, 3)
356 assert_(v.shape == (3, 4))
357 for i in range(4):
358 coef = [0]*i + [1]
359 assert_almost_equal(v[..., i], cheb.chebval(x, coef))
360
361 # check for 2d x
362 x = np.array([[1, 2], [3, 4], [5, 6]])
363 v = cheb.chebvander(x, 3)
364 assert_(v.shape == (3, 2, 4))
365 for i in range(4):
366 coef = [0]*i + [1]
367 assert_almost_equal(v[..., i], cheb.chebval(x, coef))
368
369 def test_chebvander2d(self):
370 # also tests chebval2d for non-square coefficient array
371 x1, x2, x3 = self.x
372 c = np.random.random((2, 3))
373 van = cheb.chebvander2d(x1, x2, [1, 2])
374 tgt = cheb.chebval2d(x1, x2, c)
375 res = np.dot(van, c.flat)
376 assert_almost_equal(res, tgt)
377
378 # check shape
379 van = cheb.chebvander2d([x1], [x2], [1, 2])
380 assert_(van.shape == (1, 5, 6))
381
382 def test_chebvander3d(self):
383 # also tests chebval3d for non-square coefficient array
384 x1, x2, x3 = self.x
385 c = np.random.random((2, 3, 4))
386 van = cheb.chebvander3d(x1, x2, x3, [1, 2, 3])
387 tgt = cheb.chebval3d(x1, x2, x3, c)
388 res = np.dot(van, c.flat)
389 assert_almost_equal(res, tgt)
390
391 # check shape
392 van = cheb.chebvander3d([x1], [x2], [x3], [1, 2, 3])
393 assert_(van.shape == (1, 5, 24))
394
395
396 class TestFitting(TestCase):
397
398 def test_chebfit(self):
399 def f(x):
400 return x*(x - 1)*(x - 2)
401
402 # Test exceptions
403 assert_raises(ValueError, cheb.chebfit, [1], [1], -1)
404 assert_raises(TypeError, cheb.chebfit, [[1]], [1], 0)
405 assert_raises(TypeError, cheb.chebfit, [], [1], 0)
406 assert_raises(TypeError, cheb.chebfit, [1], [[[1]]], 0)
407 assert_raises(TypeError, cheb.chebfit, [1, 2], [1], 0)
408 assert_raises(TypeError, cheb.chebfit, [1], [1, 2], 0)
409 assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[[1]])
410 assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[1, 1])
411
412 # Test fit
413 x = np.linspace(0, 2)
414 y = f(x)
415 #
416 coef3 = cheb.chebfit(x, y, 3)
417 assert_equal(len(coef3), 4)
418 assert_almost_equal(cheb.chebval(x, coef3), y)
419 #
420 coef4 = cheb.chebfit(x, y, 4)
421 assert_equal(len(coef4), 5)
422 assert_almost_equal(cheb.chebval(x, coef4), y)
423 #
424 coef2d = cheb.chebfit(x, np.array([y, y]).T, 3)
425 assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
426 # test weighting
427 w = np.zeros_like(x)
428 yw = y.copy()
429 w[1::2] = 1
430 y[0::2] = 0
431 wcoef3 = cheb.chebfit(x, yw, 3, w=w)
432 assert_almost_equal(wcoef3, coef3)
433 #
434 wcoef2d = cheb.chebfit(x, np.array([yw, yw]).T, 3, w=w)
435 assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
436 # test scaling with complex values x points whose square
437 # is zero when summed.
438 x = [1, 1j, -1, -1j]
439 assert_almost_equal(cheb.chebfit(x, x, 1), [0, 1])
440
441
442 class TestCompanion(TestCase):
443
444 def test_raises(self):
445 assert_raises(ValueError, cheb.chebcompanion, [])
446 assert_raises(ValueError, cheb.chebcompanion, [1])
447
448 def test_dimensions(self):
449 for i in range(1, 5):
450 coef = [0]*i + [1]
451 assert_(cheb.chebcompanion(coef).shape == (i, i))
452
453 def test_linear_root(self):
454 assert_(cheb.chebcompanion([1, 2])[0, 0] == -.5)
455
456
457 class TestGauss(TestCase):
458
459 def test_100(self):
460 x, w = cheb.chebgauss(100)
461
462 # test orthogonality. Note that the results need to be normalized,
463 # otherwise the huge values that can arise from fast growing
464 # functions like Laguerre can be very confusing.
465 v = cheb.chebvander(x, 99)
466 vv = np.dot(v.T * w, v)
467 vd = 1/np.sqrt(vv.diagonal())
468 vv = vd[:, None] * vv * vd
469 assert_almost_equal(vv, np.eye(100))
470
471 # check that the integral of 1 is correct
472 tgt = np.pi
473 assert_almost_equal(w.sum(), tgt)
474
475
476 class TestMisc(TestCase):
477
478 def test_chebfromroots(self):
479 res = cheb.chebfromroots([])
480 assert_almost_equal(trim(res), [1])
481 for i in range(1, 5):
482 roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
483 tgt = [0]*i + [1]
484 res = cheb.chebfromroots(roots)*2**(i-1)
485 assert_almost_equal(trim(res), trim(tgt))
486
487 def test_chebroots(self):
488 assert_almost_equal(cheb.chebroots([1]), [])
489 assert_almost_equal(cheb.chebroots([1, 2]), [-.5])
490 for i in range(2, 5):
491 tgt = np.linspace(-1, 1, i)
492 res = cheb.chebroots(cheb.chebfromroots(tgt))
493 assert_almost_equal(trim(res), trim(tgt))
494
495 def test_chebtrim(self):
496 coef = [2, -1, 1, 0]
497
498 # Test exceptions
499 assert_raises(ValueError, cheb.chebtrim, coef, -1)
500
501 # Test results
502 assert_equal(cheb.chebtrim(coef), coef[:-1])
503 assert_equal(cheb.chebtrim(coef, 1), coef[:-3])
504 assert_equal(cheb.chebtrim(coef, 2), [0])
505
506 def test_chebline(self):
507 assert_equal(cheb.chebline(3, 4), [3, 4])
508
509 def test_cheb2poly(self):
510 for i in range(10):
511 assert_almost_equal(cheb.cheb2poly([0]*i + [1]), Tlist[i])
512
513 def test_poly2cheb(self):
514 for i in range(10):
515 assert_almost_equal(cheb.poly2cheb(Tlist[i]), [0]*i + [1])
516
517 def test_weight(self):
518 x = np.linspace(-1, 1, 11)[1:-1]
519 tgt = 1./(np.sqrt(1 + x) * np.sqrt(1 - x))
520 res = cheb.chebweight(x)
521 assert_almost_equal(res, tgt)
522
523 def test_chebpts1(self):
524 #test exceptions
525 assert_raises(ValueError, cheb.chebpts1, 1.5)
526 assert_raises(ValueError, cheb.chebpts1, 0)
527
528 #test points
529 tgt = [0]
530 assert_almost_equal(cheb.chebpts1(1), tgt)
531 tgt = [-0.70710678118654746, 0.70710678118654746]
532 assert_almost_equal(cheb.chebpts1(2), tgt)
533 tgt = [-0.86602540378443871, 0, 0.86602540378443871]
534 assert_almost_equal(cheb.chebpts1(3), tgt)
535 tgt = [-0.9238795325, -0.3826834323, 0.3826834323, 0.9238795325]
536 assert_almost_equal(cheb.chebpts1(4), tgt)
537
538 def test_chebpts2(self):
539 #test exceptions
540 assert_raises(ValueError, cheb.chebpts2, 1.5)
541 assert_raises(ValueError, cheb.chebpts2, 1)
542
543 #test points
544 tgt = [-1, 1]
545 assert_almost_equal(cheb.chebpts2(2), tgt)
546 tgt = [-1, 0, 1]
547 assert_almost_equal(cheb.chebpts2(3), tgt)
548 tgt = [-1, -0.5, .5, 1]
549 assert_almost_equal(cheb.chebpts2(4), tgt)
550 tgt = [-1.0, -0.707106781187, 0, 0.707106781187, 1.0]
551 assert_almost_equal(cheb.chebpts2(5), tgt)
552
553 if __name__ == "__main__":
554 run_module_suite()