codekingpro/portable-devtools
114k
1"""Tests for legendre module.
2
3"""
4from functools import reduce
5
6import numpy as np
7import numpy.polynomial.legendre as leg
8from numpy.polynomial.polynomial import polyval
9from numpy.testing import assert_, assert_almost_equal, assert_equal, assert_raises
10
11L0 = np.array([1])
12L1 = np.array([0, 1])
13L2 = np.array([-1, 0, 3]) / 2
14L3 = np.array([0, -3, 0, 5]) / 2
15L4 = np.array([3, 0, -30, 0, 35]) / 8
16L5 = np.array([0, 15, 0, -70, 0, 63]) / 8
17L6 = np.array([-5, 0, 105, 0, -315, 0, 231]) / 16
18L7 = np.array([0, -35, 0, 315, 0, -693, 0, 429]) / 16
19L8 = np.array([35, 0, -1260, 0, 6930, 0, -12012, 0, 6435]) / 128
20L9 = np.array([0, 315, 0, -4620, 0, 18018, 0, -25740, 0, 12155]) / 128
21
22Llist = [L0, L1, L2, L3, L4, L5, L6, L7, L8, L9]
23
24
25def trim(x):
26 return leg.legtrim(x, tol=1e-6)
27
28
29class TestConstants:
30
31 def test_legdomain(self):
32 assert_equal(leg.legdomain, [-1, 1])
33
34 def test_legzero(self):
35 assert_equal(leg.legzero, [0])
36
37 def test_legone(self):
38 assert_equal(leg.legone, [1])
39
40 def test_legx(self):
41 assert_equal(leg.legx, [0, 1])
42
43
44class TestArithmetic:
45 x = np.linspace(-1, 1, 100)
46
47 def test_legadd(self):
48 for i in range(5):
49 for j in range(5):
50 msg = f"At i={i}, j={j}"
51 tgt = np.zeros(max(i, j) + 1)
52 tgt[i] += 1
53 tgt[j] += 1
54 res = leg.legadd([0] * i + [1], [0] * j + [1])
55 assert_equal(trim(res), trim(tgt), err_msg=msg)
56
57 def test_legsub(self):
58 for i in range(5):
59 for j in range(5):
60 msg = f"At i={i}, j={j}"
61 tgt = np.zeros(max(i, j) + 1)
62 tgt[i] += 1
63 tgt[j] -= 1
64 res = leg.legsub([0] * i + [1], [0] * j + [1])
65 assert_equal(trim(res), trim(tgt), err_msg=msg)
66
67 def test_legmulx(self):
68 assert_equal(leg.legmulx([0]), [0])
69 assert_equal(leg.legmulx([1]), [0, 1])
70 for i in range(1, 5):
71 tmp = 2 * i + 1
72 ser = [0] * i + [1]
73 tgt = [0] * (i - 1) + [i / tmp, 0, (i + 1) / tmp]
74 assert_equal(leg.legmulx(ser), tgt)
75
76 def test_legmul(self):
77 # check values of result
78 for i in range(5):
79 pol1 = [0] * i + [1]
80 val1 = leg.legval(self.x, pol1)
81 for j in range(5):
82 msg = f"At i={i}, j={j}"
83 pol2 = [0] * j + [1]
84 val2 = leg.legval(self.x, pol2)
85 pol3 = leg.legmul(pol1, pol2)
86 val3 = leg.legval(self.x, pol3)
87 assert_(len(pol3) == i + j + 1, msg)
88 assert_almost_equal(val3, val1 * val2, err_msg=msg)
89
90 def test_legdiv(self):
91 for i in range(5):
92 for j in range(5):
93 msg = f"At i={i}, j={j}"
94 ci = [0] * i + [1]
95 cj = [0] * j + [1]
96 tgt = leg.legadd(ci, cj)
97 quo, rem = leg.legdiv(tgt, ci)
98 res = leg.legadd(leg.legmul(quo, ci), rem)
99 assert_equal(trim(res), trim(tgt), err_msg=msg)
100
101 def test_legpow(self):
102 for i in range(5):
103 for j in range(5):
104 msg = f"At i={i}, j={j}"
105 c = np.arange(i + 1)
106 tgt = reduce(leg.legmul, [c] * j, np.array([1]))
107 res = leg.legpow(c, j)
108 assert_equal(trim(res), trim(tgt), err_msg=msg)
109
110
111class TestEvaluation:
112 # coefficients of 1 + 2*x + 3*x**2
113 c1d = np.array([2., 2., 2.])
114 c2d = np.einsum('i,j->ij', c1d, c1d)
115 c3d = np.einsum('i,j,k->ijk', c1d, c1d, c1d)
116
117 # some random values in [-1, 1)
118 x = np.random.random((3, 5)) * 2 - 1
119 y = polyval(x, [1., 2., 3.])
120
121 def test_legval(self):
122 # check empty input
123 assert_equal(leg.legval([], [1]).size, 0)
124
125 # check normal input)
126 x = np.linspace(-1, 1)
127 y = [polyval(x, c) for c in Llist]
128 for i in range(10):
129 msg = f"At i={i}"
130 tgt = y[i]
131 res = leg.legval(x, [0] * i + [1])
132 assert_almost_equal(res, tgt, err_msg=msg)
133
134 # check that shape is preserved
135 for i in range(3):
136 dims = [2] * i
137 x = np.zeros(dims)
138 assert_equal(leg.legval(x, [1]).shape, dims)
139 assert_equal(leg.legval(x, [1, 0]).shape, dims)
140 assert_equal(leg.legval(x, [1, 0, 0]).shape, dims)
141
142 def test_legval2d(self):
143 x1, x2, x3 = self.x
144 y1, y2, y3 = self.y
145
146 # test exceptions
147 assert_raises(ValueError, leg.legval2d, x1, x2[:2], self.c2d)
148
149 # test values
150 tgt = y1 * y2
151 res = leg.legval2d(x1, x2, self.c2d)
152 assert_almost_equal(res, tgt)
153
154 # test shape
155 z = np.ones((2, 3))
156 res = leg.legval2d(z, z, self.c2d)
157 assert_(res.shape == (2, 3))
158
159 def test_legval3d(self):
160 x1, x2, x3 = self.x
161 y1, y2, y3 = self.y
162
163 # test exceptions
164 assert_raises(ValueError, leg.legval3d, x1, x2, x3[:2], self.c3d)
165
166 # test values
167 tgt = y1 * y2 * y3
168 res = leg.legval3d(x1, x2, x3, self.c3d)
169 assert_almost_equal(res, tgt)
170
171 # test shape
172 z = np.ones((2, 3))
173 res = leg.legval3d(z, z, z, self.c3d)
174 assert_(res.shape == (2, 3))
175
176 def test_leggrid2d(self):
177 x1, x2, x3 = self.x
178 y1, y2, y3 = self.y
179
180 # test values
181 tgt = np.einsum('i,j->ij', y1, y2)
182 res = leg.leggrid2d(x1, x2, self.c2d)
183 assert_almost_equal(res, tgt)
184
185 # test shape
186 z = np.ones((2, 3))
187 res = leg.leggrid2d(z, z, self.c2d)
188 assert_(res.shape == (2, 3) * 2)
189
190 def test_leggrid3d(self):
191 x1, x2, x3 = self.x
192 y1, y2, y3 = self.y
193
194 # test values
195 tgt = np.einsum('i,j,k->ijk', y1, y2, y3)
196 res = leg.leggrid3d(x1, x2, x3, self.c3d)
197 assert_almost_equal(res, tgt)
198
199 # test shape
200 z = np.ones((2, 3))
201 res = leg.leggrid3d(z, z, z, self.c3d)
202 assert_(res.shape == (2, 3) * 3)
203
204
205class TestIntegral:
206
207 def test_legint(self):
208 # check exceptions
209 assert_raises(TypeError, leg.legint, [0], .5)
210 assert_raises(ValueError, leg.legint, [0], -1)
211 assert_raises(ValueError, leg.legint, [0], 1, [0, 0])
212 assert_raises(ValueError, leg.legint, [0], lbnd=[0])
213 assert_raises(ValueError, leg.legint, [0], scl=[0])
214 assert_raises(TypeError, leg.legint, [0], axis=.5)
215
216 # test integration of zero polynomial
217 for i in range(2, 5):
218 k = [0] * (i - 2) + [1]
219 res = leg.legint([0], m=i, k=k)
220 assert_almost_equal(res, [0, 1])
221
222 # check single integration with integration constant
223 for i in range(5):
224 scl = i + 1
225 pol = [0] * i + [1]
226 tgt = [i] + [0] * i + [1 / scl]
227 legpol = leg.poly2leg(pol)
228 legint = leg.legint(legpol, m=1, k=[i])
229 res = leg.leg2poly(legint)
230 assert_almost_equal(trim(res), trim(tgt))
231
232 # check single integration with integration constant and lbnd
233 for i in range(5):
234 scl = i + 1
235 pol = [0] * i + [1]
236 legpol = leg.poly2leg(pol)
237 legint = leg.legint(legpol, m=1, k=[i], lbnd=-1)
238 assert_almost_equal(leg.legval(-1, legint), i)
239
240 # check single integration with integration constant and scaling
241 for i in range(5):
242 scl = i + 1
243 pol = [0] * i + [1]
244 tgt = [i] + [0] * i + [2 / scl]
245 legpol = leg.poly2leg(pol)
246 legint = leg.legint(legpol, m=1, k=[i], scl=2)
247 res = leg.leg2poly(legint)
248 assert_almost_equal(trim(res), trim(tgt))
249
250 # check multiple integrations with default k
251 for i in range(5):
252 for j in range(2, 5):
253 pol = [0] * i + [1]
254 tgt = pol[:]
255 for k in range(j):
256 tgt = leg.legint(tgt, m=1)
257 res = leg.legint(pol, m=j)
258 assert_almost_equal(trim(res), trim(tgt))
259
260 # check multiple integrations with defined k
261 for i in range(5):
262 for j in range(2, 5):
263 pol = [0] * i + [1]
264 tgt = pol[:]
265 for k in range(j):
266 tgt = leg.legint(tgt, m=1, k=[k])
267 res = leg.legint(pol, m=j, k=list(range(j)))
268 assert_almost_equal(trim(res), trim(tgt))
269
270 # check multiple integrations with lbnd
271 for i in range(5):
272 for j in range(2, 5):
273 pol = [0] * i + [1]
274 tgt = pol[:]
275 for k in range(j):
276 tgt = leg.legint(tgt, m=1, k=[k], lbnd=-1)
277 res = leg.legint(pol, m=j, k=list(range(j)), lbnd=-1)
278 assert_almost_equal(trim(res), trim(tgt))
279
280 # check multiple integrations with scaling
281 for i in range(5):
282 for j in range(2, 5):
283 pol = [0] * i + [1]
284 tgt = pol[:]
285 for k in range(j):
286 tgt = leg.legint(tgt, m=1, k=[k], scl=2)
287 res = leg.legint(pol, m=j, k=list(range(j)), scl=2)
288 assert_almost_equal(trim(res), trim(tgt))
289
290 def test_legint_axis(self):
291 # check that axis keyword works
292 c2d = np.random.random((3, 4))
293
294 tgt = np.vstack([leg.legint(c) for c in c2d.T]).T
295 res = leg.legint(c2d, axis=0)
296 assert_almost_equal(res, tgt)
297
298 tgt = np.vstack([leg.legint(c) for c in c2d])
299 res = leg.legint(c2d, axis=1)
300 assert_almost_equal(res, tgt)
301
302 tgt = np.vstack([leg.legint(c, k=3) for c in c2d])
303 res = leg.legint(c2d, k=3, axis=1)
304 assert_almost_equal(res, tgt)
305
306 def test_legint_zerointord(self):
307 assert_equal(leg.legint((1, 2, 3), 0), (1, 2, 3))
308
309
310class TestDerivative:
311
312 def test_legder(self):
313 # check exceptions
314 assert_raises(TypeError, leg.legder, [0], .5)
315 assert_raises(ValueError, leg.legder, [0], -1)
316
317 # check that zeroth derivative does nothing
318 for i in range(5):
319 tgt = [0] * i + [1]
320 res = leg.legder(tgt, m=0)
321 assert_equal(trim(res), trim(tgt))
322
323 # check that derivation is the inverse of integration
324 for i in range(5):
325 for j in range(2, 5):
326 tgt = [0] * i + [1]
327 res = leg.legder(leg.legint(tgt, m=j), m=j)
328 assert_almost_equal(trim(res), trim(tgt))
329
330 # check derivation with scaling
331 for i in range(5):
332 for j in range(2, 5):
333 tgt = [0] * i + [1]
334 res = leg.legder(leg.legint(tgt, m=j, scl=2), m=j, scl=.5)
335 assert_almost_equal(trim(res), trim(tgt))
336
337 def test_legder_axis(self):
338 # check that axis keyword works
339 c2d = np.random.random((3, 4))
340
341 tgt = np.vstack([leg.legder(c) for c in c2d.T]).T
342 res = leg.legder(c2d, axis=0)
343 assert_almost_equal(res, tgt)
344
345 tgt = np.vstack([leg.legder(c) for c in c2d])
346 res = leg.legder(c2d, axis=1)
347 assert_almost_equal(res, tgt)
348
349 def test_legder_orderhigherthancoeff(self):
350 c = (1, 2, 3, 4)
351 assert_equal(leg.legder(c, 4), [0])
352
353class TestVander:
354 # some random values in [-1, 1)
355 x = np.random.random((3, 5)) * 2 - 1
356
357 def test_legvander(self):
358 # check for 1d x
359 x = np.arange(3)
360 v = leg.legvander(x, 3)
361 assert_(v.shape == (3, 4))
362 for i in range(4):
363 coef = [0] * i + [1]
364 assert_almost_equal(v[..., i], leg.legval(x, coef))
365
366 # check for 2d x
367 x = np.array([[1, 2], [3, 4], [5, 6]])
368 v = leg.legvander(x, 3)
369 assert_(v.shape == (3, 2, 4))
370 for i in range(4):
371 coef = [0] * i + [1]
372 assert_almost_equal(v[..., i], leg.legval(x, coef))
373
374 def test_legvander2d(self):
375 # also tests polyval2d for non-square coefficient array
376 x1, x2, x3 = self.x
377 c = np.random.random((2, 3))
378 van = leg.legvander2d(x1, x2, [1, 2])
379 tgt = leg.legval2d(x1, x2, c)
380 res = np.dot(van, c.flat)
381 assert_almost_equal(res, tgt)
382
383 # check shape
384 van = leg.legvander2d([x1], [x2], [1, 2])
385 assert_(van.shape == (1, 5, 6))
386
387 def test_legvander3d(self):
388 # also tests polyval3d for non-square coefficient array
389 x1, x2, x3 = self.x
390 c = np.random.random((2, 3, 4))
391 van = leg.legvander3d(x1, x2, x3, [1, 2, 3])
392 tgt = leg.legval3d(x1, x2, x3, c)
393 res = np.dot(van, c.flat)
394 assert_almost_equal(res, tgt)
395
396 # check shape
397 van = leg.legvander3d([x1], [x2], [x3], [1, 2, 3])
398 assert_(van.shape == (1, 5, 24))
399
400 def test_legvander_negdeg(self):
401 assert_raises(ValueError, leg.legvander, (1, 2, 3), -1)
402
403
404class TestFitting:
405
406 def test_legfit(self):
407 def f(x):
408 return x * (x - 1) * (x - 2)
409
410 def f2(x):
411 return x**4 + x**2 + 1
412
413 # Test exceptions
414 assert_raises(ValueError, leg.legfit, [1], [1], -1)
415 assert_raises(TypeError, leg.legfit, [[1]], [1], 0)
416 assert_raises(TypeError, leg.legfit, [], [1], 0)
417 assert_raises(TypeError, leg.legfit, [1], [[[1]]], 0)
418 assert_raises(TypeError, leg.legfit, [1, 2], [1], 0)
419 assert_raises(TypeError, leg.legfit, [1], [1, 2], 0)
420 assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[[1]])
421 assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[1, 1])
422 assert_raises(ValueError, leg.legfit, [1], [1], [-1,])
423 assert_raises(ValueError, leg.legfit, [1], [1], [2, -1, 6])
424 assert_raises(TypeError, leg.legfit, [1], [1], [])
425
426 # Test fit
427 x = np.linspace(0, 2)
428 y = f(x)
429 #
430 coef3 = leg.legfit(x, y, 3)
431 assert_equal(len(coef3), 4)
432 assert_almost_equal(leg.legval(x, coef3), y)
433 coef3 = leg.legfit(x, y, [0, 1, 2, 3])
434 assert_equal(len(coef3), 4)
435 assert_almost_equal(leg.legval(x, coef3), y)
436 #
437 coef4 = leg.legfit(x, y, 4)
438 assert_equal(len(coef4), 5)
439 assert_almost_equal(leg.legval(x, coef4), y)
440 coef4 = leg.legfit(x, y, [0, 1, 2, 3, 4])
441 assert_equal(len(coef4), 5)
442 assert_almost_equal(leg.legval(x, coef4), y)
443 # check things still work if deg is not in strict increasing
444 coef4 = leg.legfit(x, y, [2, 3, 4, 1, 0])
445 assert_equal(len(coef4), 5)
446 assert_almost_equal(leg.legval(x, coef4), y)
447 #
448 coef2d = leg.legfit(x, np.array([y, y]).T, 3)
449 assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
450 coef2d = leg.legfit(x, np.array([y, y]).T, [0, 1, 2, 3])
451 assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
452 # test weighting
453 w = np.zeros_like(x)
454 yw = y.copy()
455 w[1::2] = 1
456 y[0::2] = 0
457 wcoef3 = leg.legfit(x, yw, 3, w=w)
458 assert_almost_equal(wcoef3, coef3)
459 wcoef3 = leg.legfit(x, yw, [0, 1, 2, 3], w=w)
460 assert_almost_equal(wcoef3, coef3)
461 #
462 wcoef2d = leg.legfit(x, np.array([yw, yw]).T, 3, w=w)
463 assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
464 wcoef2d = leg.legfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
465 assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
466 # test scaling with complex values x points whose square
467 # is zero when summed.
468 x = [1, 1j, -1, -1j]
469 assert_almost_equal(leg.legfit(x, x, 1), [0, 1])
470 assert_almost_equal(leg.legfit(x, x, [0, 1]), [0, 1])
471 # test fitting only even Legendre polynomials
472 x = np.linspace(-1, 1)
473 y = f2(x)
474 coef1 = leg.legfit(x, y, 4)
475 assert_almost_equal(leg.legval(x, coef1), y)
476 coef2 = leg.legfit(x, y, [0, 2, 4])
477 assert_almost_equal(leg.legval(x, coef2), y)
478 assert_almost_equal(coef1, coef2)
479
480
481class TestCompanion:
482
483 def test_raises(self):
484 assert_raises(ValueError, leg.legcompanion, [])
485 assert_raises(ValueError, leg.legcompanion, [1])
486
487 def test_dimensions(self):
488 for i in range(1, 5):
489 coef = [0] * i + [1]
490 assert_(leg.legcompanion(coef).shape == (i, i))
491
492 def test_linear_root(self):
493 assert_(leg.legcompanion([1, 2])[0, 0] == -.5)
494
495
496class TestGauss:
497
498 def test_100(self):
499 x, w = leg.leggauss(100)
500
501 # test orthogonality. Note that the results need to be normalized,
502 # otherwise the huge values that can arise from fast growing
503 # functions like Laguerre can be very confusing.
504 v = leg.legvander(x, 99)
505 vv = np.dot(v.T * w, v)
506 vd = 1 / np.sqrt(vv.diagonal())
507 vv = vd[:, None] * vv * vd
508 assert_almost_equal(vv, np.eye(100))
509
510 # check that the integral of 1 is correct
511 tgt = 2.0
512 assert_almost_equal(w.sum(), tgt)
513
514
515class TestMisc:
516
517 def test_legfromroots(self):
518 res = leg.legfromroots([])
519 assert_almost_equal(trim(res), [1])
520 for i in range(1, 5):
521 roots = np.cos(np.linspace(-np.pi, 0, 2 * i + 1)[1::2])
522 pol = leg.legfromroots(roots)
523 res = leg.legval(roots, pol)
524 tgt = 0
525 assert_(len(pol) == i + 1)
526 assert_almost_equal(leg.leg2poly(pol)[-1], 1)
527 assert_almost_equal(res, tgt)
528
529 def test_legroots(self):
530 assert_almost_equal(leg.legroots([1]), [])
531 assert_almost_equal(leg.legroots([1, 2]), [-.5])
532 for i in range(2, 5):
533 tgt = np.linspace(-1, 1, i)
534 res = leg.legroots(leg.legfromroots(tgt))
535 assert_almost_equal(trim(res), trim(tgt))
536
537 def test_legtrim(self):
538 coef = [2, -1, 1, 0]
539
540 # Test exceptions
541 assert_raises(ValueError, leg.legtrim, coef, -1)
542
543 # Test results
544 assert_equal(leg.legtrim(coef), coef[:-1])
545 assert_equal(leg.legtrim(coef, 1), coef[:-3])
546 assert_equal(leg.legtrim(coef, 2), [0])
547
548 def test_legline(self):
549 assert_equal(leg.legline(3, 4), [3, 4])
550
551 def test_legline_zeroscl(self):
552 assert_equal(leg.legline(3, 0), [3])
553
554 def test_leg2poly(self):
555 for i in range(10):
556 assert_almost_equal(leg.leg2poly([0] * i + [1]), Llist[i])
557
558 def test_poly2leg(self):
559 for i in range(10):
560 assert_almost_equal(leg.poly2leg(Llist[i]), [0] * i + [1])
561
562 def test_weight(self):
563 x = np.linspace(-1, 1, 11)
564 tgt = 1.
565 res = leg.legweight(x)
566 assert_almost_equal(res, tgt)
567