图片解析应用
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

600 lines
20 KiB

  1. """Tests for polynomial module.
  2. """
  3. from functools import reduce
  4. import numpy as np
  5. import numpy.polynomial.polynomial as poly
  6. from numpy.testing import (
  7. assert_almost_equal, assert_raises, assert_equal, assert_,
  8. assert_warns, assert_array_equal, assert_raises_regex)
  9. def trim(x):
  10. return poly.polytrim(x, tol=1e-6)
  11. T0 = [1]
  12. T1 = [0, 1]
  13. T2 = [-1, 0, 2]
  14. T3 = [0, -3, 0, 4]
  15. T4 = [1, 0, -8, 0, 8]
  16. T5 = [0, 5, 0, -20, 0, 16]
  17. T6 = [-1, 0, 18, 0, -48, 0, 32]
  18. T7 = [0, -7, 0, 56, 0, -112, 0, 64]
  19. T8 = [1, 0, -32, 0, 160, 0, -256, 0, 128]
  20. T9 = [0, 9, 0, -120, 0, 432, 0, -576, 0, 256]
  21. Tlist = [T0, T1, T2, T3, T4, T5, T6, T7, T8, T9]
  22. class TestConstants:
  23. def test_polydomain(self):
  24. assert_equal(poly.polydomain, [-1, 1])
  25. def test_polyzero(self):
  26. assert_equal(poly.polyzero, [0])
  27. def test_polyone(self):
  28. assert_equal(poly.polyone, [1])
  29. def test_polyx(self):
  30. assert_equal(poly.polyx, [0, 1])
  31. class TestArithmetic:
  32. def test_polyadd(self):
  33. for i in range(5):
  34. for j in range(5):
  35. msg = f"At i={i}, j={j}"
  36. tgt = np.zeros(max(i, j) + 1)
  37. tgt[i] += 1
  38. tgt[j] += 1
  39. res = poly.polyadd([0]*i + [1], [0]*j + [1])
  40. assert_equal(trim(res), trim(tgt), err_msg=msg)
  41. def test_polysub(self):
  42. for i in range(5):
  43. for j in range(5):
  44. msg = f"At i={i}, j={j}"
  45. tgt = np.zeros(max(i, j) + 1)
  46. tgt[i] += 1
  47. tgt[j] -= 1
  48. res = poly.polysub([0]*i + [1], [0]*j + [1])
  49. assert_equal(trim(res), trim(tgt), err_msg=msg)
  50. def test_polymulx(self):
  51. assert_equal(poly.polymulx([0]), [0])
  52. assert_equal(poly.polymulx([1]), [0, 1])
  53. for i in range(1, 5):
  54. ser = [0]*i + [1]
  55. tgt = [0]*(i + 1) + [1]
  56. assert_equal(poly.polymulx(ser), tgt)
  57. def test_polymul(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(i + j + 1)
  62. tgt[i + j] += 1
  63. res = poly.polymul([0]*i + [1], [0]*j + [1])
  64. assert_equal(trim(res), trim(tgt), err_msg=msg)
  65. def test_polydiv(self):
  66. # check zero division
  67. assert_raises(ZeroDivisionError, poly.polydiv, [1], [0])
  68. # check scalar division
  69. quo, rem = poly.polydiv([2], [2])
  70. assert_equal((quo, rem), (1, 0))
  71. quo, rem = poly.polydiv([2, 2], [2])
  72. assert_equal((quo, rem), ((1, 1), 0))
  73. # check rest.
  74. for i in range(5):
  75. for j in range(5):
  76. msg = f"At i={i}, j={j}"
  77. ci = [0]*i + [1, 2]
  78. cj = [0]*j + [1, 2]
  79. tgt = poly.polyadd(ci, cj)
  80. quo, rem = poly.polydiv(tgt, ci)
  81. res = poly.polyadd(poly.polymul(quo, ci), rem)
  82. assert_equal(res, tgt, err_msg=msg)
  83. def test_polypow(self):
  84. for i in range(5):
  85. for j in range(5):
  86. msg = f"At i={i}, j={j}"
  87. c = np.arange(i + 1)
  88. tgt = reduce(poly.polymul, [c]*j, np.array([1]))
  89. res = poly.polypow(c, j)
  90. assert_equal(trim(res), trim(tgt), err_msg=msg)
  91. class TestEvaluation:
  92. # coefficients of 1 + 2*x + 3*x**2
  93. c1d = np.array([1., 2., 3.])
  94. c2d = np.einsum('i,j->ij', c1d, c1d)
  95. c3d = np.einsum('i,j,k->ijk', c1d, c1d, c1d)
  96. # some random values in [-1, 1)
  97. x = np.random.random((3, 5))*2 - 1
  98. y = poly.polyval(x, [1., 2., 3.])
  99. def test_polyval(self):
  100. #check empty input
  101. assert_equal(poly.polyval([], [1]).size, 0)
  102. #check normal input)
  103. x = np.linspace(-1, 1)
  104. y = [x**i for i in range(5)]
  105. for i in range(5):
  106. tgt = y[i]
  107. res = poly.polyval(x, [0]*i + [1])
  108. assert_almost_equal(res, tgt)
  109. tgt = x*(x**2 - 1)
  110. res = poly.polyval(x, [0, -1, 0, 1])
  111. assert_almost_equal(res, tgt)
  112. #check that shape is preserved
  113. for i in range(3):
  114. dims = [2]*i
  115. x = np.zeros(dims)
  116. assert_equal(poly.polyval(x, [1]).shape, dims)
  117. assert_equal(poly.polyval(x, [1, 0]).shape, dims)
  118. assert_equal(poly.polyval(x, [1, 0, 0]).shape, dims)
  119. #check masked arrays are processed correctly
  120. mask = [False, True, False]
  121. mx = np.ma.array([1, 2, 3], mask=mask)
  122. res = np.polyval([7, 5, 3], mx)
  123. assert_array_equal(res.mask, mask)
  124. #check subtypes of ndarray are preserved
  125. class C(np.ndarray):
  126. pass
  127. cx = np.array([1, 2, 3]).view(C)
  128. assert_equal(type(np.polyval([2, 3, 4], cx)), C)
  129. def test_polyvalfromroots(self):
  130. # check exception for broadcasting x values over root array with
  131. # too few dimensions
  132. assert_raises(ValueError, poly.polyvalfromroots,
  133. [1], [1], tensor=False)
  134. # check empty input
  135. assert_equal(poly.polyvalfromroots([], [1]).size, 0)
  136. assert_(poly.polyvalfromroots([], [1]).shape == (0,))
  137. # check empty input + multidimensional roots
  138. assert_equal(poly.polyvalfromroots([], [[1] * 5]).size, 0)
  139. assert_(poly.polyvalfromroots([], [[1] * 5]).shape == (5, 0))
  140. # check scalar input
  141. assert_equal(poly.polyvalfromroots(1, 1), 0)
  142. assert_(poly.polyvalfromroots(1, np.ones((3, 3))).shape == (3,))
  143. # check normal input)
  144. x = np.linspace(-1, 1)
  145. y = [x**i for i in range(5)]
  146. for i in range(1, 5):
  147. tgt = y[i]
  148. res = poly.polyvalfromroots(x, [0]*i)
  149. assert_almost_equal(res, tgt)
  150. tgt = x*(x - 1)*(x + 1)
  151. res = poly.polyvalfromroots(x, [-1, 0, 1])
  152. assert_almost_equal(res, tgt)
  153. # check that shape is preserved
  154. for i in range(3):
  155. dims = [2]*i
  156. x = np.zeros(dims)
  157. assert_equal(poly.polyvalfromroots(x, [1]).shape, dims)
  158. assert_equal(poly.polyvalfromroots(x, [1, 0]).shape, dims)
  159. assert_equal(poly.polyvalfromroots(x, [1, 0, 0]).shape, dims)
  160. # check compatibility with factorization
  161. ptest = [15, 2, -16, -2, 1]
  162. r = poly.polyroots(ptest)
  163. x = np.linspace(-1, 1)
  164. assert_almost_equal(poly.polyval(x, ptest),
  165. poly.polyvalfromroots(x, r))
  166. # check multidimensional arrays of roots and values
  167. # check tensor=False
  168. rshape = (3, 5)
  169. x = np.arange(-3, 2)
  170. r = np.random.randint(-5, 5, size=rshape)
  171. res = poly.polyvalfromroots(x, r, tensor=False)
  172. tgt = np.empty(r.shape[1:])
  173. for ii in range(tgt.size):
  174. tgt[ii] = poly.polyvalfromroots(x[ii], r[:, ii])
  175. assert_equal(res, tgt)
  176. # check tensor=True
  177. x = np.vstack([x, 2*x])
  178. res = poly.polyvalfromroots(x, r, tensor=True)
  179. tgt = np.empty(r.shape[1:] + x.shape)
  180. for ii in range(r.shape[1]):
  181. for jj in range(x.shape[0]):
  182. tgt[ii, jj, :] = poly.polyvalfromroots(x[jj], r[:, ii])
  183. assert_equal(res, tgt)
  184. def test_polyval2d(self):
  185. x1, x2, x3 = self.x
  186. y1, y2, y3 = self.y
  187. #test exceptions
  188. assert_raises_regex(ValueError, 'incompatible',
  189. poly.polyval2d, x1, x2[:2], self.c2d)
  190. #test values
  191. tgt = y1*y2
  192. res = poly.polyval2d(x1, x2, self.c2d)
  193. assert_almost_equal(res, tgt)
  194. #test shape
  195. z = np.ones((2, 3))
  196. res = poly.polyval2d(z, z, self.c2d)
  197. assert_(res.shape == (2, 3))
  198. def test_polyval3d(self):
  199. x1, x2, x3 = self.x
  200. y1, y2, y3 = self.y
  201. #test exceptions
  202. assert_raises_regex(ValueError, 'incompatible',
  203. poly.polyval3d, x1, x2, x3[:2], self.c3d)
  204. #test values
  205. tgt = y1*y2*y3
  206. res = poly.polyval3d(x1, x2, x3, self.c3d)
  207. assert_almost_equal(res, tgt)
  208. #test shape
  209. z = np.ones((2, 3))
  210. res = poly.polyval3d(z, z, z, self.c3d)
  211. assert_(res.shape == (2, 3))
  212. def test_polygrid2d(self):
  213. x1, x2, x3 = self.x
  214. y1, y2, y3 = self.y
  215. #test values
  216. tgt = np.einsum('i,j->ij', y1, y2)
  217. res = poly.polygrid2d(x1, x2, self.c2d)
  218. assert_almost_equal(res, tgt)
  219. #test shape
  220. z = np.ones((2, 3))
  221. res = poly.polygrid2d(z, z, self.c2d)
  222. assert_(res.shape == (2, 3)*2)
  223. def test_polygrid3d(self):
  224. x1, x2, x3 = self.x
  225. y1, y2, y3 = self.y
  226. #test values
  227. tgt = np.einsum('i,j,k->ijk', y1, y2, y3)
  228. res = poly.polygrid3d(x1, x2, x3, self.c3d)
  229. assert_almost_equal(res, tgt)
  230. #test shape
  231. z = np.ones((2, 3))
  232. res = poly.polygrid3d(z, z, z, self.c3d)
  233. assert_(res.shape == (2, 3)*3)
  234. class TestIntegral:
  235. def test_polyint(self):
  236. # check exceptions
  237. assert_raises(TypeError, poly.polyint, [0], .5)
  238. assert_raises(ValueError, poly.polyint, [0], -1)
  239. assert_raises(ValueError, poly.polyint, [0], 1, [0, 0])
  240. assert_raises(ValueError, poly.polyint, [0], lbnd=[0])
  241. assert_raises(ValueError, poly.polyint, [0], scl=[0])
  242. assert_raises(TypeError, poly.polyint, [0], axis=.5)
  243. with assert_warns(DeprecationWarning):
  244. poly.polyint([1, 1], 1.)
  245. # test integration of zero polynomial
  246. for i in range(2, 5):
  247. k = [0]*(i - 2) + [1]
  248. res = poly.polyint([0], m=i, k=k)
  249. assert_almost_equal(res, [0, 1])
  250. # check single integration with integration constant
  251. for i in range(5):
  252. scl = i + 1
  253. pol = [0]*i + [1]
  254. tgt = [i] + [0]*i + [1/scl]
  255. res = poly.polyint(pol, m=1, k=[i])
  256. assert_almost_equal(trim(res), trim(tgt))
  257. # check single integration with integration constant and lbnd
  258. for i in range(5):
  259. scl = i + 1
  260. pol = [0]*i + [1]
  261. res = poly.polyint(pol, m=1, k=[i], lbnd=-1)
  262. assert_almost_equal(poly.polyval(-1, res), i)
  263. # check single integration with integration constant and scaling
  264. for i in range(5):
  265. scl = i + 1
  266. pol = [0]*i + [1]
  267. tgt = [i] + [0]*i + [2/scl]
  268. res = poly.polyint(pol, m=1, k=[i], scl=2)
  269. assert_almost_equal(trim(res), trim(tgt))
  270. # check multiple integrations with default k
  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 = poly.polyint(tgt, m=1)
  277. res = poly.polyint(pol, m=j)
  278. assert_almost_equal(trim(res), trim(tgt))
  279. # check multiple integrations with defined k
  280. for i in range(5):
  281. for j in range(2, 5):
  282. pol = [0]*i + [1]
  283. tgt = pol[:]
  284. for k in range(j):
  285. tgt = poly.polyint(tgt, m=1, k=[k])
  286. res = poly.polyint(pol, m=j, k=list(range(j)))
  287. assert_almost_equal(trim(res), trim(tgt))
  288. # check multiple integrations with lbnd
  289. for i in range(5):
  290. for j in range(2, 5):
  291. pol = [0]*i + [1]
  292. tgt = pol[:]
  293. for k in range(j):
  294. tgt = poly.polyint(tgt, m=1, k=[k], lbnd=-1)
  295. res = poly.polyint(pol, m=j, k=list(range(j)), lbnd=-1)
  296. assert_almost_equal(trim(res), trim(tgt))
  297. # check multiple integrations with scaling
  298. for i in range(5):
  299. for j in range(2, 5):
  300. pol = [0]*i + [1]
  301. tgt = pol[:]
  302. for k in range(j):
  303. tgt = poly.polyint(tgt, m=1, k=[k], scl=2)
  304. res = poly.polyint(pol, m=j, k=list(range(j)), scl=2)
  305. assert_almost_equal(trim(res), trim(tgt))
  306. def test_polyint_axis(self):
  307. # check that axis keyword works
  308. c2d = np.random.random((3, 4))
  309. tgt = np.vstack([poly.polyint(c) for c in c2d.T]).T
  310. res = poly.polyint(c2d, axis=0)
  311. assert_almost_equal(res, tgt)
  312. tgt = np.vstack([poly.polyint(c) for c in c2d])
  313. res = poly.polyint(c2d, axis=1)
  314. assert_almost_equal(res, tgt)
  315. tgt = np.vstack([poly.polyint(c, k=3) for c in c2d])
  316. res = poly.polyint(c2d, k=3, axis=1)
  317. assert_almost_equal(res, tgt)
  318. class TestDerivative:
  319. def test_polyder(self):
  320. # check exceptions
  321. assert_raises(TypeError, poly.polyder, [0], .5)
  322. assert_raises(ValueError, poly.polyder, [0], -1)
  323. # check that zeroth derivative does nothing
  324. for i in range(5):
  325. tgt = [0]*i + [1]
  326. res = poly.polyder(tgt, m=0)
  327. assert_equal(trim(res), trim(tgt))
  328. # check that derivation is the inverse of integration
  329. for i in range(5):
  330. for j in range(2, 5):
  331. tgt = [0]*i + [1]
  332. res = poly.polyder(poly.polyint(tgt, m=j), m=j)
  333. assert_almost_equal(trim(res), trim(tgt))
  334. # check derivation with scaling
  335. for i in range(5):
  336. for j in range(2, 5):
  337. tgt = [0]*i + [1]
  338. res = poly.polyder(poly.polyint(tgt, m=j, scl=2), m=j, scl=.5)
  339. assert_almost_equal(trim(res), trim(tgt))
  340. def test_polyder_axis(self):
  341. # check that axis keyword works
  342. c2d = np.random.random((3, 4))
  343. tgt = np.vstack([poly.polyder(c) for c in c2d.T]).T
  344. res = poly.polyder(c2d, axis=0)
  345. assert_almost_equal(res, tgt)
  346. tgt = np.vstack([poly.polyder(c) for c in c2d])
  347. res = poly.polyder(c2d, axis=1)
  348. assert_almost_equal(res, tgt)
  349. class TestVander:
  350. # some random values in [-1, 1)
  351. x = np.random.random((3, 5))*2 - 1
  352. def test_polyvander(self):
  353. # check for 1d x
  354. x = np.arange(3)
  355. v = poly.polyvander(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], poly.polyval(x, coef))
  360. # check for 2d x
  361. x = np.array([[1, 2], [3, 4], [5, 6]])
  362. v = poly.polyvander(x, 3)
  363. assert_(v.shape == (3, 2, 4))
  364. for i in range(4):
  365. coef = [0]*i + [1]
  366. assert_almost_equal(v[..., i], poly.polyval(x, coef))
  367. def test_polyvander2d(self):
  368. # also tests polyval2d for non-square coefficient array
  369. x1, x2, x3 = self.x
  370. c = np.random.random((2, 3))
  371. van = poly.polyvander2d(x1, x2, [1, 2])
  372. tgt = poly.polyval2d(x1, x2, c)
  373. res = np.dot(van, c.flat)
  374. assert_almost_equal(res, tgt)
  375. # check shape
  376. van = poly.polyvander2d([x1], [x2], [1, 2])
  377. assert_(van.shape == (1, 5, 6))
  378. def test_polyvander3d(self):
  379. # also tests polyval3d for non-square coefficient array
  380. x1, x2, x3 = self.x
  381. c = np.random.random((2, 3, 4))
  382. van = poly.polyvander3d(x1, x2, x3, [1, 2, 3])
  383. tgt = poly.polyval3d(x1, x2, x3, c)
  384. res = np.dot(van, c.flat)
  385. assert_almost_equal(res, tgt)
  386. # check shape
  387. van = poly.polyvander3d([x1], [x2], [x3], [1, 2, 3])
  388. assert_(van.shape == (1, 5, 24))
  389. def test_polyvandernegdeg(self):
  390. x = np.arange(3)
  391. assert_raises(ValueError, poly.polyvander, x, -1)
  392. class TestCompanion:
  393. def test_raises(self):
  394. assert_raises(ValueError, poly.polycompanion, [])
  395. assert_raises(ValueError, poly.polycompanion, [1])
  396. def test_dimensions(self):
  397. for i in range(1, 5):
  398. coef = [0]*i + [1]
  399. assert_(poly.polycompanion(coef).shape == (i, i))
  400. def test_linear_root(self):
  401. assert_(poly.polycompanion([1, 2])[0, 0] == -.5)
  402. class TestMisc:
  403. def test_polyfromroots(self):
  404. res = poly.polyfromroots([])
  405. assert_almost_equal(trim(res), [1])
  406. for i in range(1, 5):
  407. roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
  408. tgt = Tlist[i]
  409. res = poly.polyfromroots(roots)*2**(i-1)
  410. assert_almost_equal(trim(res), trim(tgt))
  411. def test_polyroots(self):
  412. assert_almost_equal(poly.polyroots([1]), [])
  413. assert_almost_equal(poly.polyroots([1, 2]), [-.5])
  414. for i in range(2, 5):
  415. tgt = np.linspace(-1, 1, i)
  416. res = poly.polyroots(poly.polyfromroots(tgt))
  417. assert_almost_equal(trim(res), trim(tgt))
  418. def test_polyfit(self):
  419. def f(x):
  420. return x*(x - 1)*(x - 2)
  421. def f2(x):
  422. return x**4 + x**2 + 1
  423. # Test exceptions
  424. assert_raises(ValueError, poly.polyfit, [1], [1], -1)
  425. assert_raises(TypeError, poly.polyfit, [[1]], [1], 0)
  426. assert_raises(TypeError, poly.polyfit, [], [1], 0)
  427. assert_raises(TypeError, poly.polyfit, [1], [[[1]]], 0)
  428. assert_raises(TypeError, poly.polyfit, [1, 2], [1], 0)
  429. assert_raises(TypeError, poly.polyfit, [1], [1, 2], 0)
  430. assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[[1]])
  431. assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[1, 1])
  432. assert_raises(ValueError, poly.polyfit, [1], [1], [-1,])
  433. assert_raises(ValueError, poly.polyfit, [1], [1], [2, -1, 6])
  434. assert_raises(TypeError, poly.polyfit, [1], [1], [])
  435. # Test fit
  436. x = np.linspace(0, 2)
  437. y = f(x)
  438. #
  439. coef3 = poly.polyfit(x, y, 3)
  440. assert_equal(len(coef3), 4)
  441. assert_almost_equal(poly.polyval(x, coef3), y)
  442. coef3 = poly.polyfit(x, y, [0, 1, 2, 3])
  443. assert_equal(len(coef3), 4)
  444. assert_almost_equal(poly.polyval(x, coef3), y)
  445. #
  446. coef4 = poly.polyfit(x, y, 4)
  447. assert_equal(len(coef4), 5)
  448. assert_almost_equal(poly.polyval(x, coef4), y)
  449. coef4 = poly.polyfit(x, y, [0, 1, 2, 3, 4])
  450. assert_equal(len(coef4), 5)
  451. assert_almost_equal(poly.polyval(x, coef4), y)
  452. #
  453. coef2d = poly.polyfit(x, np.array([y, y]).T, 3)
  454. assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
  455. coef2d = poly.polyfit(x, np.array([y, y]).T, [0, 1, 2, 3])
  456. assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
  457. # test weighting
  458. w = np.zeros_like(x)
  459. yw = y.copy()
  460. w[1::2] = 1
  461. yw[0::2] = 0
  462. wcoef3 = poly.polyfit(x, yw, 3, w=w)
  463. assert_almost_equal(wcoef3, coef3)
  464. wcoef3 = poly.polyfit(x, yw, [0, 1, 2, 3], w=w)
  465. assert_almost_equal(wcoef3, coef3)
  466. #
  467. wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, 3, w=w)
  468. assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
  469. wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
  470. assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
  471. # test scaling with complex values x points whose square
  472. # is zero when summed.
  473. x = [1, 1j, -1, -1j]
  474. assert_almost_equal(poly.polyfit(x, x, 1), [0, 1])
  475. assert_almost_equal(poly.polyfit(x, x, [0, 1]), [0, 1])
  476. # test fitting only even Polyendre polynomials
  477. x = np.linspace(-1, 1)
  478. y = f2(x)
  479. coef1 = poly.polyfit(x, y, 4)
  480. assert_almost_equal(poly.polyval(x, coef1), y)
  481. coef2 = poly.polyfit(x, y, [0, 2, 4])
  482. assert_almost_equal(poly.polyval(x, coef2), y)
  483. assert_almost_equal(coef1, coef2)
  484. def test_polytrim(self):
  485. coef = [2, -1, 1, 0]
  486. # Test exceptions
  487. assert_raises(ValueError, poly.polytrim, coef, -1)
  488. # Test results
  489. assert_equal(poly.polytrim(coef), coef[:-1])
  490. assert_equal(poly.polytrim(coef, 1), coef[:-3])
  491. assert_equal(poly.polytrim(coef, 2), [0])
  492. def test_polyline(self):
  493. assert_equal(poly.polyline(3, 4), [3, 4])
  494. def test_polyline_zero(self):
  495. assert_equal(poly.polyline(3, 0), [3])