test_laguerre.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535
  1. """Tests for laguerre module.
  2. """
  3. from functools import reduce
  4. import numpy as np
  5. import numpy.polynomial.laguerre as lag
  6. from numpy.polynomial.polynomial import polyval
  7. from numpy.testing import assert_, assert_almost_equal, assert_equal, assert_raises
  8. L0 = np.array([1]) / 1
  9. L1 = np.array([1, -1]) / 1
  10. L2 = np.array([2, -4, 1]) / 2
  11. L3 = np.array([6, -18, 9, -1]) / 6
  12. L4 = np.array([24, -96, 72, -16, 1]) / 24
  13. L5 = np.array([120, -600, 600, -200, 25, -1]) / 120
  14. L6 = np.array([720, -4320, 5400, -2400, 450, -36, 1]) / 720
  15. Llist = [L0, L1, L2, L3, L4, L5, L6]
  16. def trim(x):
  17. return lag.lagtrim(x, tol=1e-6)
  18. class TestConstants:
  19. def test_lagdomain(self):
  20. assert_equal(lag.lagdomain, [0, 1])
  21. def test_lagzero(self):
  22. assert_equal(lag.lagzero, [0])
  23. def test_lagone(self):
  24. assert_equal(lag.lagone, [1])
  25. def test_lagx(self):
  26. assert_equal(lag.lagx, [1, -1])
  27. class TestArithmetic:
  28. x = np.linspace(-3, 3, 100)
  29. def test_lagadd(self):
  30. for i in range(5):
  31. for j in range(5):
  32. msg = f"At i={i}, j={j}"
  33. tgt = np.zeros(max(i, j) + 1)
  34. tgt[i] += 1
  35. tgt[j] += 1
  36. res = lag.lagadd([0] * i + [1], [0] * j + [1])
  37. assert_equal(trim(res), trim(tgt), err_msg=msg)
  38. def test_lagsub(self):
  39. for i in range(5):
  40. for j in range(5):
  41. msg = f"At i={i}, j={j}"
  42. tgt = np.zeros(max(i, j) + 1)
  43. tgt[i] += 1
  44. tgt[j] -= 1
  45. res = lag.lagsub([0] * i + [1], [0] * j + [1])
  46. assert_equal(trim(res), trim(tgt), err_msg=msg)
  47. def test_lagmulx(self):
  48. assert_equal(lag.lagmulx([0]), [0])
  49. assert_equal(lag.lagmulx([1]), [1, -1])
  50. for i in range(1, 5):
  51. ser = [0] * i + [1]
  52. tgt = [0] * (i - 1) + [-i, 2 * i + 1, -(i + 1)]
  53. assert_almost_equal(lag.lagmulx(ser), tgt)
  54. def test_lagmul(self):
  55. # check values of result
  56. for i in range(5):
  57. pol1 = [0] * i + [1]
  58. val1 = lag.lagval(self.x, pol1)
  59. for j in range(5):
  60. msg = f"At i={i}, j={j}"
  61. pol2 = [0] * j + [1]
  62. val2 = lag.lagval(self.x, pol2)
  63. pol3 = lag.lagmul(pol1, pol2)
  64. val3 = lag.lagval(self.x, pol3)
  65. assert_(len(pol3) == i + j + 1, msg)
  66. assert_almost_equal(val3, val1 * val2, err_msg=msg)
  67. def test_lagdiv(self):
  68. for i in range(5):
  69. for j in range(5):
  70. msg = f"At i={i}, j={j}"
  71. ci = [0] * i + [1]
  72. cj = [0] * j + [1]
  73. tgt = lag.lagadd(ci, cj)
  74. quo, rem = lag.lagdiv(tgt, ci)
  75. res = lag.lagadd(lag.lagmul(quo, ci), rem)
  76. assert_almost_equal(trim(res), trim(tgt), err_msg=msg)
  77. def test_lagpow(self):
  78. for i in range(5):
  79. for j in range(5):
  80. msg = f"At i={i}, j={j}"
  81. c = np.arange(i + 1)
  82. tgt = reduce(lag.lagmul, [c] * j, np.array([1]))
  83. res = lag.lagpow(c, j)
  84. assert_equal(trim(res), trim(tgt), err_msg=msg)
  85. class TestEvaluation:
  86. # coefficients of 1 + 2*x + 3*x**2
  87. c1d = np.array([9., -14., 6.])
  88. c2d = np.einsum('i,j->ij', c1d, c1d)
  89. c3d = np.einsum('i,j,k->ijk', c1d, c1d, c1d)
  90. # some random values in [-1, 1)
  91. x = np.random.random((3, 5)) * 2 - 1
  92. y = polyval(x, [1., 2., 3.])
  93. def test_lagval(self):
  94. # check empty input
  95. assert_equal(lag.lagval([], [1]).size, 0)
  96. # check normal input)
  97. x = np.linspace(-1, 1)
  98. y = [polyval(x, c) for c in Llist]
  99. for i in range(7):
  100. msg = f"At i={i}"
  101. tgt = y[i]
  102. res = lag.lagval(x, [0] * i + [1])
  103. assert_almost_equal(res, tgt, err_msg=msg)
  104. # check that shape is preserved
  105. for i in range(3):
  106. dims = [2] * i
  107. x = np.zeros(dims)
  108. assert_equal(lag.lagval(x, [1]).shape, dims)
  109. assert_equal(lag.lagval(x, [1, 0]).shape, dims)
  110. assert_equal(lag.lagval(x, [1, 0, 0]).shape, dims)
  111. def test_lagval2d(self):
  112. x1, x2, x3 = self.x
  113. y1, y2, y3 = self.y
  114. # test exceptions
  115. assert_raises(ValueError, lag.lagval2d, x1, x2[:2], self.c2d)
  116. # test values
  117. tgt = y1 * y2
  118. res = lag.lagval2d(x1, x2, self.c2d)
  119. assert_almost_equal(res, tgt)
  120. # test shape
  121. z = np.ones((2, 3))
  122. res = lag.lagval2d(z, z, self.c2d)
  123. assert_(res.shape == (2, 3))
  124. def test_lagval3d(self):
  125. x1, x2, x3 = self.x
  126. y1, y2, y3 = self.y
  127. # test exceptions
  128. assert_raises(ValueError, lag.lagval3d, x1, x2, x3[:2], self.c3d)
  129. # test values
  130. tgt = y1 * y2 * y3
  131. res = lag.lagval3d(x1, x2, x3, self.c3d)
  132. assert_almost_equal(res, tgt)
  133. # test shape
  134. z = np.ones((2, 3))
  135. res = lag.lagval3d(z, z, z, self.c3d)
  136. assert_(res.shape == (2, 3))
  137. def test_laggrid2d(self):
  138. x1, x2, x3 = self.x
  139. y1, y2, y3 = self.y
  140. # test values
  141. tgt = np.einsum('i,j->ij', y1, y2)
  142. res = lag.laggrid2d(x1, x2, self.c2d)
  143. assert_almost_equal(res, tgt)
  144. # test shape
  145. z = np.ones((2, 3))
  146. res = lag.laggrid2d(z, z, self.c2d)
  147. assert_(res.shape == (2, 3) * 2)
  148. def test_laggrid3d(self):
  149. x1, x2, x3 = self.x
  150. y1, y2, y3 = self.y
  151. # test values
  152. tgt = np.einsum('i,j,k->ijk', y1, y2, y3)
  153. res = lag.laggrid3d(x1, x2, x3, self.c3d)
  154. assert_almost_equal(res, tgt)
  155. # test shape
  156. z = np.ones((2, 3))
  157. res = lag.laggrid3d(z, z, z, self.c3d)
  158. assert_(res.shape == (2, 3) * 3)
  159. class TestIntegral:
  160. def test_lagint(self):
  161. # check exceptions
  162. assert_raises(TypeError, lag.lagint, [0], .5)
  163. assert_raises(ValueError, lag.lagint, [0], -1)
  164. assert_raises(ValueError, lag.lagint, [0], 1, [0, 0])
  165. assert_raises(ValueError, lag.lagint, [0], lbnd=[0])
  166. assert_raises(ValueError, lag.lagint, [0], scl=[0])
  167. assert_raises(TypeError, lag.lagint, [0], axis=.5)
  168. # test integration of zero polynomial
  169. for i in range(2, 5):
  170. k = [0] * (i - 2) + [1]
  171. res = lag.lagint([0], m=i, k=k)
  172. assert_almost_equal(res, [1, -1])
  173. # check single integration with integration constant
  174. for i in range(5):
  175. scl = i + 1
  176. pol = [0] * i + [1]
  177. tgt = [i] + [0] * i + [1 / scl]
  178. lagpol = lag.poly2lag(pol)
  179. lagint = lag.lagint(lagpol, m=1, k=[i])
  180. res = lag.lag2poly(lagint)
  181. assert_almost_equal(trim(res), trim(tgt))
  182. # check single integration with integration constant and lbnd
  183. for i in range(5):
  184. scl = i + 1
  185. pol = [0] * i + [1]
  186. lagpol = lag.poly2lag(pol)
  187. lagint = lag.lagint(lagpol, m=1, k=[i], lbnd=-1)
  188. assert_almost_equal(lag.lagval(-1, lagint), i)
  189. # check single integration with integration constant and scaling
  190. for i in range(5):
  191. scl = i + 1
  192. pol = [0] * i + [1]
  193. tgt = [i] + [0] * i + [2 / scl]
  194. lagpol = lag.poly2lag(pol)
  195. lagint = lag.lagint(lagpol, m=1, k=[i], scl=2)
  196. res = lag.lag2poly(lagint)
  197. assert_almost_equal(trim(res), trim(tgt))
  198. # check multiple integrations with default k
  199. for i in range(5):
  200. for j in range(2, 5):
  201. pol = [0] * i + [1]
  202. tgt = pol[:]
  203. for k in range(j):
  204. tgt = lag.lagint(tgt, m=1)
  205. res = lag.lagint(pol, m=j)
  206. assert_almost_equal(trim(res), trim(tgt))
  207. # check multiple integrations with defined k
  208. for i in range(5):
  209. for j in range(2, 5):
  210. pol = [0] * i + [1]
  211. tgt = pol[:]
  212. for k in range(j):
  213. tgt = lag.lagint(tgt, m=1, k=[k])
  214. res = lag.lagint(pol, m=j, k=list(range(j)))
  215. assert_almost_equal(trim(res), trim(tgt))
  216. # check multiple integrations with lbnd
  217. for i in range(5):
  218. for j in range(2, 5):
  219. pol = [0] * i + [1]
  220. tgt = pol[:]
  221. for k in range(j):
  222. tgt = lag.lagint(tgt, m=1, k=[k], lbnd=-1)
  223. res = lag.lagint(pol, m=j, k=list(range(j)), lbnd=-1)
  224. assert_almost_equal(trim(res), trim(tgt))
  225. # check multiple integrations with scaling
  226. for i in range(5):
  227. for j in range(2, 5):
  228. pol = [0] * i + [1]
  229. tgt = pol[:]
  230. for k in range(j):
  231. tgt = lag.lagint(tgt, m=1, k=[k], scl=2)
  232. res = lag.lagint(pol, m=j, k=list(range(j)), scl=2)
  233. assert_almost_equal(trim(res), trim(tgt))
  234. def test_lagint_axis(self):
  235. # check that axis keyword works
  236. c2d = np.random.random((3, 4))
  237. tgt = np.vstack([lag.lagint(c) for c in c2d.T]).T
  238. res = lag.lagint(c2d, axis=0)
  239. assert_almost_equal(res, tgt)
  240. tgt = np.vstack([lag.lagint(c) for c in c2d])
  241. res = lag.lagint(c2d, axis=1)
  242. assert_almost_equal(res, tgt)
  243. tgt = np.vstack([lag.lagint(c, k=3) for c in c2d])
  244. res = lag.lagint(c2d, k=3, axis=1)
  245. assert_almost_equal(res, tgt)
  246. class TestDerivative:
  247. def test_lagder(self):
  248. # check exceptions
  249. assert_raises(TypeError, lag.lagder, [0], .5)
  250. assert_raises(ValueError, lag.lagder, [0], -1)
  251. # check that zeroth derivative does nothing
  252. for i in range(5):
  253. tgt = [0] * i + [1]
  254. res = lag.lagder(tgt, m=0)
  255. assert_equal(trim(res), trim(tgt))
  256. # check that derivation is the inverse of integration
  257. for i in range(5):
  258. for j in range(2, 5):
  259. tgt = [0] * i + [1]
  260. res = lag.lagder(lag.lagint(tgt, m=j), m=j)
  261. assert_almost_equal(trim(res), trim(tgt))
  262. # check derivation with scaling
  263. for i in range(5):
  264. for j in range(2, 5):
  265. tgt = [0] * i + [1]
  266. res = lag.lagder(lag.lagint(tgt, m=j, scl=2), m=j, scl=.5)
  267. assert_almost_equal(trim(res), trim(tgt))
  268. def test_lagder_axis(self):
  269. # check that axis keyword works
  270. c2d = np.random.random((3, 4))
  271. tgt = np.vstack([lag.lagder(c) for c in c2d.T]).T
  272. res = lag.lagder(c2d, axis=0)
  273. assert_almost_equal(res, tgt)
  274. tgt = np.vstack([lag.lagder(c) for c in c2d])
  275. res = lag.lagder(c2d, axis=1)
  276. assert_almost_equal(res, tgt)
  277. class TestVander:
  278. # some random values in [-1, 1)
  279. x = np.random.random((3, 5)) * 2 - 1
  280. def test_lagvander(self):
  281. # check for 1d x
  282. x = np.arange(3)
  283. v = lag.lagvander(x, 3)
  284. assert_(v.shape == (3, 4))
  285. for i in range(4):
  286. coef = [0] * i + [1]
  287. assert_almost_equal(v[..., i], lag.lagval(x, coef))
  288. # check for 2d x
  289. x = np.array([[1, 2], [3, 4], [5, 6]])
  290. v = lag.lagvander(x, 3)
  291. assert_(v.shape == (3, 2, 4))
  292. for i in range(4):
  293. coef = [0] * i + [1]
  294. assert_almost_equal(v[..., i], lag.lagval(x, coef))
  295. def test_lagvander2d(self):
  296. # also tests lagval2d for non-square coefficient array
  297. x1, x2, x3 = self.x
  298. c = np.random.random((2, 3))
  299. van = lag.lagvander2d(x1, x2, [1, 2])
  300. tgt = lag.lagval2d(x1, x2, c)
  301. res = np.dot(van, c.flat)
  302. assert_almost_equal(res, tgt)
  303. # check shape
  304. van = lag.lagvander2d([x1], [x2], [1, 2])
  305. assert_(van.shape == (1, 5, 6))
  306. def test_lagvander3d(self):
  307. # also tests lagval3d for non-square coefficient array
  308. x1, x2, x3 = self.x
  309. c = np.random.random((2, 3, 4))
  310. van = lag.lagvander3d(x1, x2, x3, [1, 2, 3])
  311. tgt = lag.lagval3d(x1, x2, x3, c)
  312. res = np.dot(van, c.flat)
  313. assert_almost_equal(res, tgt)
  314. # check shape
  315. van = lag.lagvander3d([x1], [x2], [x3], [1, 2, 3])
  316. assert_(van.shape == (1, 5, 24))
  317. class TestFitting:
  318. def test_lagfit(self):
  319. def f(x):
  320. return x * (x - 1) * (x - 2)
  321. # Test exceptions
  322. assert_raises(ValueError, lag.lagfit, [1], [1], -1)
  323. assert_raises(TypeError, lag.lagfit, [[1]], [1], 0)
  324. assert_raises(TypeError, lag.lagfit, [], [1], 0)
  325. assert_raises(TypeError, lag.lagfit, [1], [[[1]]], 0)
  326. assert_raises(TypeError, lag.lagfit, [1, 2], [1], 0)
  327. assert_raises(TypeError, lag.lagfit, [1], [1, 2], 0)
  328. assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[[1]])
  329. assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[1, 1])
  330. assert_raises(ValueError, lag.lagfit, [1], [1], [-1,])
  331. assert_raises(ValueError, lag.lagfit, [1], [1], [2, -1, 6])
  332. assert_raises(TypeError, lag.lagfit, [1], [1], [])
  333. # Test fit
  334. x = np.linspace(0, 2)
  335. y = f(x)
  336. #
  337. coef3 = lag.lagfit(x, y, 3)
  338. assert_equal(len(coef3), 4)
  339. assert_almost_equal(lag.lagval(x, coef3), y)
  340. coef3 = lag.lagfit(x, y, [0, 1, 2, 3])
  341. assert_equal(len(coef3), 4)
  342. assert_almost_equal(lag.lagval(x, coef3), y)
  343. #
  344. coef4 = lag.lagfit(x, y, 4)
  345. assert_equal(len(coef4), 5)
  346. assert_almost_equal(lag.lagval(x, coef4), y)
  347. coef4 = lag.lagfit(x, y, [0, 1, 2, 3, 4])
  348. assert_equal(len(coef4), 5)
  349. assert_almost_equal(lag.lagval(x, coef4), y)
  350. #
  351. coef2d = lag.lagfit(x, np.array([y, y]).T, 3)
  352. assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
  353. coef2d = lag.lagfit(x, np.array([y, y]).T, [0, 1, 2, 3])
  354. assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
  355. # test weighting
  356. w = np.zeros_like(x)
  357. yw = y.copy()
  358. w[1::2] = 1
  359. y[0::2] = 0
  360. wcoef3 = lag.lagfit(x, yw, 3, w=w)
  361. assert_almost_equal(wcoef3, coef3)
  362. wcoef3 = lag.lagfit(x, yw, [0, 1, 2, 3], w=w)
  363. assert_almost_equal(wcoef3, coef3)
  364. #
  365. wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, 3, w=w)
  366. assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
  367. wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
  368. assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
  369. # test scaling with complex values x points whose square
  370. # is zero when summed.
  371. x = [1, 1j, -1, -1j]
  372. assert_almost_equal(lag.lagfit(x, x, 1), [1, -1])
  373. assert_almost_equal(lag.lagfit(x, x, [0, 1]), [1, -1])
  374. class TestCompanion:
  375. def test_raises(self):
  376. assert_raises(ValueError, lag.lagcompanion, [])
  377. assert_raises(ValueError, lag.lagcompanion, [1])
  378. def test_dimensions(self):
  379. for i in range(1, 5):
  380. coef = [0] * i + [1]
  381. assert_(lag.lagcompanion(coef).shape == (i, i))
  382. def test_linear_root(self):
  383. assert_(lag.lagcompanion([1, 2])[0, 0] == 1.5)
  384. class TestGauss:
  385. def test_100(self):
  386. x, w = lag.laggauss(100)
  387. # test orthogonality. Note that the results need to be normalized,
  388. # otherwise the huge values that can arise from fast growing
  389. # functions like Laguerre can be very confusing.
  390. v = lag.lagvander(x, 99)
  391. vv = np.dot(v.T * w, v)
  392. vd = 1 / np.sqrt(vv.diagonal())
  393. vv = vd[:, None] * vv * vd
  394. assert_almost_equal(vv, np.eye(100))
  395. # check that the integral of 1 is correct
  396. tgt = 1.0
  397. assert_almost_equal(w.sum(), tgt)
  398. class TestMisc:
  399. def test_lagfromroots(self):
  400. res = lag.lagfromroots([])
  401. assert_almost_equal(trim(res), [1])
  402. for i in range(1, 5):
  403. roots = np.cos(np.linspace(-np.pi, 0, 2 * i + 1)[1::2])
  404. pol = lag.lagfromroots(roots)
  405. res = lag.lagval(roots, pol)
  406. tgt = 0
  407. assert_(len(pol) == i + 1)
  408. assert_almost_equal(lag.lag2poly(pol)[-1], 1)
  409. assert_almost_equal(res, tgt)
  410. def test_lagroots(self):
  411. assert_almost_equal(lag.lagroots([1]), [])
  412. assert_almost_equal(lag.lagroots([0, 1]), [1])
  413. for i in range(2, 5):
  414. tgt = np.linspace(0, 3, i)
  415. res = lag.lagroots(lag.lagfromroots(tgt))
  416. assert_almost_equal(trim(res), trim(tgt))
  417. def test_lagtrim(self):
  418. coef = [2, -1, 1, 0]
  419. # Test exceptions
  420. assert_raises(ValueError, lag.lagtrim, coef, -1)
  421. # Test results
  422. assert_equal(lag.lagtrim(coef), coef[:-1])
  423. assert_equal(lag.lagtrim(coef, 1), coef[:-3])
  424. assert_equal(lag.lagtrim(coef, 2), [0])
  425. def test_lagline(self):
  426. assert_equal(lag.lagline(3, 4), [7, -4])
  427. def test_lag2poly(self):
  428. for i in range(7):
  429. assert_almost_equal(lag.lag2poly([0] * i + [1]), Llist[i])
  430. def test_poly2lag(self):
  431. for i in range(7):
  432. assert_almost_equal(lag.poly2lag(Llist[i]), [0] * i + [1])
  433. def test_weight(self):
  434. x = np.linspace(0, 10, 11)
  435. tgt = np.exp(-x)
  436. res = lag.lagweight(x)
  437. assert_almost_equal(res, tgt)