polynomial.py 48 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556
  1. """
  2. Objects for dealing with polynomials.
  3. This module provides a number of objects (mostly functions) useful for
  4. dealing with polynomials, including a `Polynomial` class that
  5. encapsulates the usual arithmetic operations. (General information
  6. on how this module represents and works with polynomial objects is in
  7. the docstring for its "parent" sub-package, `numpy.polynomial`).
  8. Constants
  9. ---------
  10. - `polydomain` -- Polynomial default domain, [-1,1].
  11. - `polyzero` -- (Coefficients of the) "zero polynomial."
  12. - `polyone` -- (Coefficients of the) constant polynomial 1.
  13. - `polyx` -- (Coefficients of the) identity map polynomial, ``f(x) = x``.
  14. Arithmetic
  15. ----------
  16. - `polyadd` -- add two polynomials.
  17. - `polysub` -- subtract one polynomial from another.
  18. - `polymul` -- multiply two polynomials.
  19. - `polydiv` -- divide one polynomial by another.
  20. - `polypow` -- raise a polynomial to an positive integer power
  21. - `polyval` -- evaluate a polynomial at given points.
  22. - `polyval2d` -- evaluate a 2D polynomial at given points.
  23. - `polyval3d` -- evaluate a 3D polynomial at given points.
  24. - `polygrid2d` -- evaluate a 2D polynomial on a Cartesian product.
  25. - `polygrid3d` -- evaluate a 3D polynomial on a Cartesian product.
  26. Calculus
  27. --------
  28. - `polyder` -- differentiate a polynomial.
  29. - `polyint` -- integrate a polynomial.
  30. Misc Functions
  31. --------------
  32. - `polyfromroots` -- create a polynomial with specified roots.
  33. - `polyroots` -- find the roots of a polynomial.
  34. - `polyvander` -- Vandermonde-like matrix for powers.
  35. - `polyvander2d` -- Vandermonde-like matrix for 2D power series.
  36. - `polyvander3d` -- Vandermonde-like matrix for 3D power series.
  37. - `polycompanion` -- companion matrix in power series form.
  38. - `polyfit` -- least-squares fit returning a polynomial.
  39. - `polytrim` -- trim leading coefficients from a polynomial.
  40. - `polyline` -- polynomial representing given straight line.
  41. Classes
  42. -------
  43. - `Polynomial` -- polynomial class.
  44. See Also
  45. --------
  46. `numpy.polynomial`
  47. """
  48. from __future__ import division, absolute_import, print_function
  49. __all__ = [
  50. 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd',
  51. 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval',
  52. 'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
  53. 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
  54. 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d']
  55. import warnings
  56. import numpy as np
  57. import numpy.linalg as la
  58. from . import polyutils as pu
  59. from ._polybase import ABCPolyBase
  60. polytrim = pu.trimcoef
  61. #
  62. # These are constant arrays are of integer type so as to be compatible
  63. # with the widest range of other types, such as Decimal.
  64. #
  65. # Polynomial default domain.
  66. polydomain = np.array([-1, 1])
  67. # Polynomial coefficients representing zero.
  68. polyzero = np.array([0])
  69. # Polynomial coefficients representing one.
  70. polyone = np.array([1])
  71. # Polynomial coefficients representing the identity x.
  72. polyx = np.array([0, 1])
  73. #
  74. # Polynomial series functions
  75. #
  76. def polyline(off, scl):
  77. """
  78. Returns an array representing a linear polynomial.
  79. Parameters
  80. ----------
  81. off, scl : scalars
  82. The "y-intercept" and "slope" of the line, respectively.
  83. Returns
  84. -------
  85. y : ndarray
  86. This module's representation of the linear polynomial ``off +
  87. scl*x``.
  88. See Also
  89. --------
  90. chebline
  91. Examples
  92. --------
  93. >>> from numpy.polynomial import polynomial as P
  94. >>> P.polyline(1,-1)
  95. array([ 1, -1])
  96. >>> P.polyval(1, P.polyline(1,-1)) # should be 0
  97. 0.0
  98. """
  99. if scl != 0:
  100. return np.array([off, scl])
  101. else:
  102. return np.array([off])
  103. def polyfromroots(roots):
  104. """
  105. Generate a monic polynomial with given roots.
  106. Return the coefficients of the polynomial
  107. .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
  108. where the `r_n` are the roots specified in `roots`. If a zero has
  109. multiplicity n, then it must appear in `roots` n times. For instance,
  110. if 2 is a root of multiplicity three and 3 is a root of multiplicity 2,
  111. then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear
  112. in any order.
  113. If the returned coefficients are `c`, then
  114. .. math:: p(x) = c_0 + c_1 * x + ... + x^n
  115. The coefficient of the last term is 1 for monic polynomials in this
  116. form.
  117. Parameters
  118. ----------
  119. roots : array_like
  120. Sequence containing the roots.
  121. Returns
  122. -------
  123. out : ndarray
  124. 1-D array of the polynomial's coefficients If all the roots are
  125. real, then `out` is also real, otherwise it is complex. (see
  126. Examples below).
  127. See Also
  128. --------
  129. chebfromroots, legfromroots, lagfromroots, hermfromroots
  130. hermefromroots
  131. Notes
  132. -----
  133. The coefficients are determined by multiplying together linear factors
  134. of the form `(x - r_i)`, i.e.
  135. .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n)
  136. where ``n == len(roots) - 1``; note that this implies that `1` is always
  137. returned for :math:`a_n`.
  138. Examples
  139. --------
  140. >>> from numpy.polynomial import polynomial as P
  141. >>> P.polyfromroots((-1,0,1)) # x(x - 1)(x + 1) = x^3 - x
  142. array([ 0., -1., 0., 1.])
  143. >>> j = complex(0,1)
  144. >>> P.polyfromroots((-j,j)) # complex returned, though values are real
  145. array([ 1.+0.j, 0.+0.j, 1.+0.j])
  146. """
  147. if len(roots) == 0:
  148. return np.ones(1)
  149. else:
  150. [roots] = pu.as_series([roots], trim=False)
  151. roots.sort()
  152. p = [polyline(-r, 1) for r in roots]
  153. n = len(p)
  154. while n > 1:
  155. m, r = divmod(n, 2)
  156. tmp = [polymul(p[i], p[i+m]) for i in range(m)]
  157. if r:
  158. tmp[0] = polymul(tmp[0], p[-1])
  159. p = tmp
  160. n = m
  161. return p[0]
  162. def polyadd(c1, c2):
  163. """
  164. Add one polynomial to another.
  165. Returns the sum of two polynomials `c1` + `c2`. The arguments are
  166. sequences of coefficients from lowest order term to highest, i.e.,
  167. [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
  168. Parameters
  169. ----------
  170. c1, c2 : array_like
  171. 1-D arrays of polynomial coefficients ordered from low to high.
  172. Returns
  173. -------
  174. out : ndarray
  175. The coefficient array representing their sum.
  176. See Also
  177. --------
  178. polysub, polymul, polydiv, polypow
  179. Examples
  180. --------
  181. >>> from numpy.polynomial import polynomial as P
  182. >>> c1 = (1,2,3)
  183. >>> c2 = (3,2,1)
  184. >>> sum = P.polyadd(c1,c2); sum
  185. array([ 4., 4., 4.])
  186. >>> P.polyval(2, sum) # 4 + 4(2) + 4(2**2)
  187. 28.0
  188. """
  189. # c1, c2 are trimmed copies
  190. [c1, c2] = pu.as_series([c1, c2])
  191. if len(c1) > len(c2):
  192. c1[:c2.size] += c2
  193. ret = c1
  194. else:
  195. c2[:c1.size] += c1
  196. ret = c2
  197. return pu.trimseq(ret)
  198. def polysub(c1, c2):
  199. """
  200. Subtract one polynomial from another.
  201. Returns the difference of two polynomials `c1` - `c2`. The arguments
  202. are sequences of coefficients from lowest order term to highest, i.e.,
  203. [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
  204. Parameters
  205. ----------
  206. c1, c2 : array_like
  207. 1-D arrays of polynomial coefficients ordered from low to
  208. high.
  209. Returns
  210. -------
  211. out : ndarray
  212. Of coefficients representing their difference.
  213. See Also
  214. --------
  215. polyadd, polymul, polydiv, polypow
  216. Examples
  217. --------
  218. >>> from numpy.polynomial import polynomial as P
  219. >>> c1 = (1,2,3)
  220. >>> c2 = (3,2,1)
  221. >>> P.polysub(c1,c2)
  222. array([-2., 0., 2.])
  223. >>> P.polysub(c2,c1) # -P.polysub(c1,c2)
  224. array([ 2., 0., -2.])
  225. """
  226. # c1, c2 are trimmed copies
  227. [c1, c2] = pu.as_series([c1, c2])
  228. if len(c1) > len(c2):
  229. c1[:c2.size] -= c2
  230. ret = c1
  231. else:
  232. c2 = -c2
  233. c2[:c1.size] += c1
  234. ret = c2
  235. return pu.trimseq(ret)
  236. def polymulx(c):
  237. """Multiply a polynomial by x.
  238. Multiply the polynomial `c` by x, where x is the independent
  239. variable.
  240. Parameters
  241. ----------
  242. c : array_like
  243. 1-D array of polynomial coefficients ordered from low to
  244. high.
  245. Returns
  246. -------
  247. out : ndarray
  248. Array representing the result of the multiplication.
  249. Notes
  250. -----
  251. .. versionadded:: 1.5.0
  252. """
  253. # c is a trimmed copy
  254. [c] = pu.as_series([c])
  255. # The zero series needs special treatment
  256. if len(c) == 1 and c[0] == 0:
  257. return c
  258. prd = np.empty(len(c) + 1, dtype=c.dtype)
  259. prd[0] = c[0]*0
  260. prd[1:] = c
  261. return prd
  262. def polymul(c1, c2):
  263. """
  264. Multiply one polynomial by another.
  265. Returns the product of two polynomials `c1` * `c2`. The arguments are
  266. sequences of coefficients, from lowest order term to highest, e.g.,
  267. [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.``
  268. Parameters
  269. ----------
  270. c1, c2 : array_like
  271. 1-D arrays of coefficients representing a polynomial, relative to the
  272. "standard" basis, and ordered from lowest order term to highest.
  273. Returns
  274. -------
  275. out : ndarray
  276. Of the coefficients of their product.
  277. See Also
  278. --------
  279. polyadd, polysub, polydiv, polypow
  280. Examples
  281. --------
  282. >>> from numpy.polynomial import polynomial as P
  283. >>> c1 = (1,2,3)
  284. >>> c2 = (3,2,1)
  285. >>> P.polymul(c1,c2)
  286. array([ 3., 8., 14., 8., 3.])
  287. """
  288. # c1, c2 are trimmed copies
  289. [c1, c2] = pu.as_series([c1, c2])
  290. ret = np.convolve(c1, c2)
  291. return pu.trimseq(ret)
  292. def polydiv(c1, c2):
  293. """
  294. Divide one polynomial by another.
  295. Returns the quotient-with-remainder of two polynomials `c1` / `c2`.
  296. The arguments are sequences of coefficients, from lowest order term
  297. to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``.
  298. Parameters
  299. ----------
  300. c1, c2 : array_like
  301. 1-D arrays of polynomial coefficients ordered from low to high.
  302. Returns
  303. -------
  304. [quo, rem] : ndarrays
  305. Of coefficient series representing the quotient and remainder.
  306. See Also
  307. --------
  308. polyadd, polysub, polymul, polypow
  309. Examples
  310. --------
  311. >>> from numpy.polynomial import polynomial as P
  312. >>> c1 = (1,2,3)
  313. >>> c2 = (3,2,1)
  314. >>> P.polydiv(c1,c2)
  315. (array([ 3.]), array([-8., -4.]))
  316. >>> P.polydiv(c2,c1)
  317. (array([ 0.33333333]), array([ 2.66666667, 1.33333333]))
  318. """
  319. # c1, c2 are trimmed copies
  320. [c1, c2] = pu.as_series([c1, c2])
  321. if c2[-1] == 0:
  322. raise ZeroDivisionError()
  323. len1 = len(c1)
  324. len2 = len(c2)
  325. if len2 == 1:
  326. return c1/c2[-1], c1[:1]*0
  327. elif len1 < len2:
  328. return c1[:1]*0, c1
  329. else:
  330. dlen = len1 - len2
  331. scl = c2[-1]
  332. c2 = c2[:-1]/scl
  333. i = dlen
  334. j = len1 - 1
  335. while i >= 0:
  336. c1[i:j] -= c2*c1[j]
  337. i -= 1
  338. j -= 1
  339. return c1[j+1:]/scl, pu.trimseq(c1[:j+1])
  340. def polypow(c, pow, maxpower=None):
  341. """Raise a polynomial to a power.
  342. Returns the polynomial `c` raised to the power `pow`. The argument
  343. `c` is a sequence of coefficients ordered from low to high. i.e.,
  344. [1,2,3] is the series ``1 + 2*x + 3*x**2.``
  345. Parameters
  346. ----------
  347. c : array_like
  348. 1-D array of array of series coefficients ordered from low to
  349. high degree.
  350. pow : integer
  351. Power to which the series will be raised
  352. maxpower : integer, optional
  353. Maximum power allowed. This is mainly to limit growth of the series
  354. to unmanageable size. Default is 16
  355. Returns
  356. -------
  357. coef : ndarray
  358. Power series of power.
  359. See Also
  360. --------
  361. polyadd, polysub, polymul, polydiv
  362. Examples
  363. --------
  364. """
  365. # c is a trimmed copy
  366. [c] = pu.as_series([c])
  367. power = int(pow)
  368. if power != pow or power < 0:
  369. raise ValueError("Power must be a non-negative integer.")
  370. elif maxpower is not None and power > maxpower:
  371. raise ValueError("Power is too large")
  372. elif power == 0:
  373. return np.array([1], dtype=c.dtype)
  374. elif power == 1:
  375. return c
  376. else:
  377. # This can be made more efficient by using powers of two
  378. # in the usual way.
  379. prd = c
  380. for i in range(2, power + 1):
  381. prd = np.convolve(prd, c)
  382. return prd
  383. def polyder(c, m=1, scl=1, axis=0):
  384. """
  385. Differentiate a polynomial.
  386. Returns the polynomial coefficients `c` differentiated `m` times along
  387. `axis`. At each iteration the result is multiplied by `scl` (the
  388. scaling factor is for use in a linear change of variable). The
  389. argument `c` is an array of coefficients from low to high degree along
  390. each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``
  391. while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is
  392. ``x`` and axis=1 is ``y``.
  393. Parameters
  394. ----------
  395. c : array_like
  396. Array of polynomial coefficients. If c is multidimensional the
  397. different axis correspond to different variables with the degree
  398. in each axis given by the corresponding index.
  399. m : int, optional
  400. Number of derivatives taken, must be non-negative. (Default: 1)
  401. scl : scalar, optional
  402. Each differentiation is multiplied by `scl`. The end result is
  403. multiplication by ``scl**m``. This is for use in a linear change
  404. of variable. (Default: 1)
  405. axis : int, optional
  406. Axis over which the derivative is taken. (Default: 0).
  407. .. versionadded:: 1.7.0
  408. Returns
  409. -------
  410. der : ndarray
  411. Polynomial coefficients of the derivative.
  412. See Also
  413. --------
  414. polyint
  415. Examples
  416. --------
  417. >>> from numpy.polynomial import polynomial as P
  418. >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3
  419. >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2
  420. array([ 2., 6., 12.])
  421. >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24
  422. array([ 24.])
  423. >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2
  424. array([ -2., -6., -12.])
  425. >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x
  426. array([ 6., 24.])
  427. """
  428. c = np.array(c, ndmin=1, copy=1)
  429. if c.dtype.char in '?bBhHiIlLqQpP':
  430. # astype fails with NA
  431. c = c + 0.0
  432. cdt = c.dtype
  433. cnt, iaxis = [int(t) for t in [m, axis]]
  434. if cnt != m:
  435. raise ValueError("The order of derivation must be integer")
  436. if cnt < 0:
  437. raise ValueError("The order of derivation must be non-negative")
  438. if iaxis != axis:
  439. raise ValueError("The axis must be integer")
  440. if not -c.ndim <= iaxis < c.ndim:
  441. raise ValueError("The axis is out of range")
  442. if iaxis < 0:
  443. iaxis += c.ndim
  444. if cnt == 0:
  445. return c
  446. c = np.rollaxis(c, iaxis)
  447. n = len(c)
  448. if cnt >= n:
  449. c = c[:1]*0
  450. else:
  451. for i in range(cnt):
  452. n = n - 1
  453. c *= scl
  454. der = np.empty((n,) + c.shape[1:], dtype=cdt)
  455. for j in range(n, 0, -1):
  456. der[j - 1] = j*c[j]
  457. c = der
  458. c = np.rollaxis(c, 0, iaxis + 1)
  459. return c
  460. def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
  461. """
  462. Integrate a polynomial.
  463. Returns the polynomial coefficients `c` integrated `m` times from
  464. `lbnd` along `axis`. At each iteration the resulting series is
  465. **multiplied** by `scl` and an integration constant, `k`, is added.
  466. The scaling factor is for use in a linear change of variable. ("Buyer
  467. beware": note that, depending on what one is doing, one may want `scl`
  468. to be the reciprocal of what one might expect; for more information,
  469. see the Notes section below.) The argument `c` is an array of
  470. coefficients, from low to high degree along each axis, e.g., [1,2,3]
  471. represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]]
  472. represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is
  473. ``y``.
  474. Parameters
  475. ----------
  476. c : array_like
  477. 1-D array of polynomial coefficients, ordered from low to high.
  478. m : int, optional
  479. Order of integration, must be positive. (Default: 1)
  480. k : {[], list, scalar}, optional
  481. Integration constant(s). The value of the first integral at zero
  482. is the first value in the list, the value of the second integral
  483. at zero is the second value, etc. If ``k == []`` (the default),
  484. all constants are set to zero. If ``m == 1``, a single scalar can
  485. be given instead of a list.
  486. lbnd : scalar, optional
  487. The lower bound of the integral. (Default: 0)
  488. scl : scalar, optional
  489. Following each integration the result is *multiplied* by `scl`
  490. before the integration constant is added. (Default: 1)
  491. axis : int, optional
  492. Axis over which the integral is taken. (Default: 0).
  493. .. versionadded:: 1.7.0
  494. Returns
  495. -------
  496. S : ndarray
  497. Coefficient array of the integral.
  498. Raises
  499. ------
  500. ValueError
  501. If ``m < 1``, ``len(k) > m``.
  502. See Also
  503. --------
  504. polyder
  505. Notes
  506. -----
  507. Note that the result of each integration is *multiplied* by `scl`. Why
  508. is this important to note? Say one is making a linear change of
  509. variable :math:`u = ax + b` in an integral relative to `x`. Then
  510. .. math::`dx = du/a`, so one will need to set `scl` equal to
  511. :math:`1/a` - perhaps not what one would have first thought.
  512. Examples
  513. --------
  514. >>> from numpy.polynomial import polynomial as P
  515. >>> c = (1,2,3)
  516. >>> P.polyint(c) # should return array([0, 1, 1, 1])
  517. array([ 0., 1., 1., 1.])
  518. >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
  519. array([ 0. , 0. , 0. , 0.16666667, 0.08333333,
  520. 0.05 ])
  521. >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1])
  522. array([ 3., 1., 1., 1.])
  523. >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1])
  524. array([ 6., 1., 1., 1.])
  525. >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2])
  526. array([ 0., -2., -2., -2.])
  527. """
  528. c = np.array(c, ndmin=1, copy=1)
  529. if c.dtype.char in '?bBhHiIlLqQpP':
  530. # astype doesn't preserve mask attribute.
  531. c = c + 0.0
  532. cdt = c.dtype
  533. if not np.iterable(k):
  534. k = [k]
  535. cnt, iaxis = [int(t) for t in [m, axis]]
  536. if cnt != m:
  537. raise ValueError("The order of integration must be integer")
  538. if cnt < 0:
  539. raise ValueError("The order of integration must be non-negative")
  540. if len(k) > cnt:
  541. raise ValueError("Too many integration constants")
  542. if iaxis != axis:
  543. raise ValueError("The axis must be integer")
  544. if not -c.ndim <= iaxis < c.ndim:
  545. raise ValueError("The axis is out of range")
  546. if iaxis < 0:
  547. iaxis += c.ndim
  548. if cnt == 0:
  549. return c
  550. k = list(k) + [0]*(cnt - len(k))
  551. c = np.rollaxis(c, iaxis)
  552. for i in range(cnt):
  553. n = len(c)
  554. c *= scl
  555. if n == 1 and np.all(c[0] == 0):
  556. c[0] += k[i]
  557. else:
  558. tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt)
  559. tmp[0] = c[0]*0
  560. tmp[1] = c[0]
  561. for j in range(1, n):
  562. tmp[j + 1] = c[j]/(j + 1)
  563. tmp[0] += k[i] - polyval(lbnd, tmp)
  564. c = tmp
  565. c = np.rollaxis(c, 0, iaxis + 1)
  566. return c
  567. def polyval(x, c, tensor=True):
  568. """
  569. Evaluate a polynomial at points x.
  570. If `c` is of length `n + 1`, this function returns the value
  571. .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n
  572. The parameter `x` is converted to an array only if it is a tuple or a
  573. list, otherwise it is treated as a scalar. In either case, either `x`
  574. or its elements must support multiplication and addition both with
  575. themselves and with the elements of `c`.
  576. If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
  577. `c` is multidimensional, then the shape of the result depends on the
  578. value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
  579. x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
  580. scalars have shape (,).
  581. Trailing zeros in the coefficients will be used in the evaluation, so
  582. they should be avoided if efficiency is a concern.
  583. Parameters
  584. ----------
  585. x : array_like, compatible object
  586. If `x` is a list or tuple, it is converted to an ndarray, otherwise
  587. it is left unchanged and treated as a scalar. In either case, `x`
  588. or its elements must support addition and multiplication with
  589. with themselves and with the elements of `c`.
  590. c : array_like
  591. Array of coefficients ordered so that the coefficients for terms of
  592. degree n are contained in c[n]. If `c` is multidimensional the
  593. remaining indices enumerate multiple polynomials. In the two
  594. dimensional case the coefficients may be thought of as stored in
  595. the columns of `c`.
  596. tensor : boolean, optional
  597. If True, the shape of the coefficient array is extended with ones
  598. on the right, one for each dimension of `x`. Scalars have dimension 0
  599. for this action. The result is that every column of coefficients in
  600. `c` is evaluated for every element of `x`. If False, `x` is broadcast
  601. over the columns of `c` for the evaluation. This keyword is useful
  602. when `c` is multidimensional. The default value is True.
  603. .. versionadded:: 1.7.0
  604. Returns
  605. -------
  606. values : ndarray, compatible object
  607. The shape of the returned array is described above.
  608. See Also
  609. --------
  610. polyval2d, polygrid2d, polyval3d, polygrid3d
  611. Notes
  612. -----
  613. The evaluation uses Horner's method.
  614. Examples
  615. --------
  616. >>> from numpy.polynomial.polynomial import polyval
  617. >>> polyval(1, [1,2,3])
  618. 6.0
  619. >>> a = np.arange(4).reshape(2,2)
  620. >>> a
  621. array([[0, 1],
  622. [2, 3]])
  623. >>> polyval(a, [1,2,3])
  624. array([[ 1., 6.],
  625. [ 17., 34.]])
  626. >>> coef = np.arange(4).reshape(2,2) # multidimensional coefficients
  627. >>> coef
  628. array([[0, 1],
  629. [2, 3]])
  630. >>> polyval([1,2], coef, tensor=True)
  631. array([[ 2., 4.],
  632. [ 4., 7.]])
  633. >>> polyval([1,2], coef, tensor=False)
  634. array([ 2., 7.])
  635. """
  636. c = np.array(c, ndmin=1, copy=0)
  637. if c.dtype.char in '?bBhHiIlLqQpP':
  638. # astype fails with NA
  639. c = c + 0.0
  640. if isinstance(x, (tuple, list)):
  641. x = np.asarray(x)
  642. if isinstance(x, np.ndarray) and tensor:
  643. c = c.reshape(c.shape + (1,)*x.ndim)
  644. c0 = c[-1] + x*0
  645. for i in range(2, len(c) + 1):
  646. c0 = c[-i] + c0*x
  647. return c0
  648. def polyval2d(x, y, c):
  649. """
  650. Evaluate a 2-D polynomial at points (x, y).
  651. This function returns the value
  652. .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j
  653. The parameters `x` and `y` are converted to arrays only if they are
  654. tuples or a lists, otherwise they are treated as a scalars and they
  655. must have the same shape after conversion. In either case, either `x`
  656. and `y` or their elements must support multiplication and addition both
  657. with themselves and with the elements of `c`.
  658. If `c` has fewer than two dimensions, ones are implicitly appended to
  659. its shape to make it 2-D. The shape of the result will be c.shape[2:] +
  660. x.shape.
  661. Parameters
  662. ----------
  663. x, y : array_like, compatible objects
  664. The two dimensional series is evaluated at the points `(x, y)`,
  665. where `x` and `y` must have the same shape. If `x` or `y` is a list
  666. or tuple, it is first converted to an ndarray, otherwise it is left
  667. unchanged and, if it isn't an ndarray, it is treated as a scalar.
  668. c : array_like
  669. Array of coefficients ordered so that the coefficient of the term
  670. of multi-degree i,j is contained in `c[i,j]`. If `c` has
  671. dimension greater than two the remaining indices enumerate multiple
  672. sets of coefficients.
  673. Returns
  674. -------
  675. values : ndarray, compatible object
  676. The values of the two dimensional polynomial at points formed with
  677. pairs of corresponding values from `x` and `y`.
  678. See Also
  679. --------
  680. polyval, polygrid2d, polyval3d, polygrid3d
  681. Notes
  682. -----
  683. .. versionadded:: 1.7.0
  684. """
  685. try:
  686. x, y = np.array((x, y), copy=0)
  687. except:
  688. raise ValueError('x, y are incompatible')
  689. c = polyval(x, c)
  690. c = polyval(y, c, tensor=False)
  691. return c
  692. def polygrid2d(x, y, c):
  693. """
  694. Evaluate a 2-D polynomial on the Cartesian product of x and y.
  695. This function returns the values:
  696. .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j
  697. where the points `(a, b)` consist of all pairs formed by taking
  698. `a` from `x` and `b` from `y`. The resulting points form a grid with
  699. `x` in the first dimension and `y` in the second.
  700. The parameters `x` and `y` are converted to arrays only if they are
  701. tuples or a lists, otherwise they are treated as a scalars. In either
  702. case, either `x` and `y` or their elements must support multiplication
  703. and addition both with themselves and with the elements of `c`.
  704. If `c` has fewer than two dimensions, ones are implicitly appended to
  705. its shape to make it 2-D. The shape of the result will be c.shape[2:] +
  706. x.shape + y.shape.
  707. Parameters
  708. ----------
  709. x, y : array_like, compatible objects
  710. The two dimensional series is evaluated at the points in the
  711. Cartesian product of `x` and `y`. If `x` or `y` is a list or
  712. tuple, it is first converted to an ndarray, otherwise it is left
  713. unchanged and, if it isn't an ndarray, it is treated as a scalar.
  714. c : array_like
  715. Array of coefficients ordered so that the coefficients for terms of
  716. degree i,j are contained in ``c[i,j]``. If `c` has dimension
  717. greater than two the remaining indices enumerate multiple sets of
  718. coefficients.
  719. Returns
  720. -------
  721. values : ndarray, compatible object
  722. The values of the two dimensional polynomial at points in the Cartesian
  723. product of `x` and `y`.
  724. See Also
  725. --------
  726. polyval, polyval2d, polyval3d, polygrid3d
  727. Notes
  728. -----
  729. .. versionadded:: 1.7.0
  730. """
  731. c = polyval(x, c)
  732. c = polyval(y, c)
  733. return c
  734. def polyval3d(x, y, z, c):
  735. """
  736. Evaluate a 3-D polynomial at points (x, y, z).
  737. This function returns the values:
  738. .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k
  739. The parameters `x`, `y`, and `z` are converted to arrays only if
  740. they are tuples or a lists, otherwise they are treated as a scalars and
  741. they must have the same shape after conversion. In either case, either
  742. `x`, `y`, and `z` or their elements must support multiplication and
  743. addition both with themselves and with the elements of `c`.
  744. If `c` has fewer than 3 dimensions, ones are implicitly appended to its
  745. shape to make it 3-D. The shape of the result will be c.shape[3:] +
  746. x.shape.
  747. Parameters
  748. ----------
  749. x, y, z : array_like, compatible object
  750. The three dimensional series is evaluated at the points
  751. `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
  752. any of `x`, `y`, or `z` is a list or tuple, it is first converted
  753. to an ndarray, otherwise it is left unchanged and if it isn't an
  754. ndarray it is treated as a scalar.
  755. c : array_like
  756. Array of coefficients ordered so that the coefficient of the term of
  757. multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
  758. greater than 3 the remaining indices enumerate multiple sets of
  759. coefficients.
  760. Returns
  761. -------
  762. values : ndarray, compatible object
  763. The values of the multidimensional polynomial on points formed with
  764. triples of corresponding values from `x`, `y`, and `z`.
  765. See Also
  766. --------
  767. polyval, polyval2d, polygrid2d, polygrid3d
  768. Notes
  769. -----
  770. .. versionadded:: 1.7.0
  771. """
  772. try:
  773. x, y, z = np.array((x, y, z), copy=0)
  774. except:
  775. raise ValueError('x, y, z are incompatible')
  776. c = polyval(x, c)
  777. c = polyval(y, c, tensor=False)
  778. c = polyval(z, c, tensor=False)
  779. return c
  780. def polygrid3d(x, y, z, c):
  781. """
  782. Evaluate a 3-D polynomial on the Cartesian product of x, y and z.
  783. This function returns the values:
  784. .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k
  785. where the points `(a, b, c)` consist of all triples formed by taking
  786. `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
  787. a grid with `x` in the first dimension, `y` in the second, and `z` in
  788. the third.
  789. The parameters `x`, `y`, and `z` are converted to arrays only if they
  790. are tuples or a lists, otherwise they are treated as a scalars. In
  791. either case, either `x`, `y`, and `z` or their elements must support
  792. multiplication and addition both with themselves and with the elements
  793. of `c`.
  794. If `c` has fewer than three dimensions, ones are implicitly appended to
  795. its shape to make it 3-D. The shape of the result will be c.shape[3:] +
  796. x.shape + y.shape + z.shape.
  797. Parameters
  798. ----------
  799. x, y, z : array_like, compatible objects
  800. The three dimensional series is evaluated at the points in the
  801. Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
  802. list or tuple, it is first converted to an ndarray, otherwise it is
  803. left unchanged and, if it isn't an ndarray, it is treated as a
  804. scalar.
  805. c : array_like
  806. Array of coefficients ordered so that the coefficients for terms of
  807. degree i,j are contained in ``c[i,j]``. If `c` has dimension
  808. greater than two the remaining indices enumerate multiple sets of
  809. coefficients.
  810. Returns
  811. -------
  812. values : ndarray, compatible object
  813. The values of the two dimensional polynomial at points in the Cartesian
  814. product of `x` and `y`.
  815. See Also
  816. --------
  817. polyval, polyval2d, polygrid2d, polyval3d
  818. Notes
  819. -----
  820. .. versionadded:: 1.7.0
  821. """
  822. c = polyval(x, c)
  823. c = polyval(y, c)
  824. c = polyval(z, c)
  825. return c
  826. def polyvander(x, deg):
  827. """Vandermonde matrix of given degree.
  828. Returns the Vandermonde matrix of degree `deg` and sample points
  829. `x`. The Vandermonde matrix is defined by
  830. .. math:: V[..., i] = x^i,
  831. where `0 <= i <= deg`. The leading indices of `V` index the elements of
  832. `x` and the last index is the power of `x`.
  833. If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
  834. matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and
  835. ``polyval(x, c)`` are the same up to roundoff. This equivalence is
  836. useful both for least squares fitting and for the evaluation of a large
  837. number of polynomials of the same degree and sample points.
  838. Parameters
  839. ----------
  840. x : array_like
  841. Array of points. The dtype is converted to float64 or complex128
  842. depending on whether any of the elements are complex. If `x` is
  843. scalar it is converted to a 1-D array.
  844. deg : int
  845. Degree of the resulting matrix.
  846. Returns
  847. -------
  848. vander : ndarray.
  849. The Vandermonde matrix. The shape of the returned matrix is
  850. ``x.shape + (deg + 1,)``, where the last index is the power of `x`.
  851. The dtype will be the same as the converted `x`.
  852. See Also
  853. --------
  854. polyvander2d, polyvander3d
  855. """
  856. ideg = int(deg)
  857. if ideg != deg:
  858. raise ValueError("deg must be integer")
  859. if ideg < 0:
  860. raise ValueError("deg must be non-negative")
  861. x = np.array(x, copy=0, ndmin=1) + 0.0
  862. dims = (ideg + 1,) + x.shape
  863. dtyp = x.dtype
  864. v = np.empty(dims, dtype=dtyp)
  865. v[0] = x*0 + 1
  866. if ideg > 0:
  867. v[1] = x
  868. for i in range(2, ideg + 1):
  869. v[i] = v[i-1]*x
  870. return np.rollaxis(v, 0, v.ndim)
  871. def polyvander2d(x, y, deg):
  872. """Pseudo-Vandermonde matrix of given degrees.
  873. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  874. points `(x, y)`. The pseudo-Vandermonde matrix is defined by
  875. .. math:: V[..., deg[1]*i + j] = x^i * y^j,
  876. where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
  877. `V` index the points `(x, y)` and the last index encodes the powers of
  878. `x` and `y`.
  879. If ``V = polyvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
  880. correspond to the elements of a 2-D coefficient array `c` of shape
  881. (xdeg + 1, ydeg + 1) in the order
  882. .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
  883. and ``np.dot(V, c.flat)`` and ``polyval2d(x, y, c)`` will be the same
  884. up to roundoff. This equivalence is useful both for least squares
  885. fitting and for the evaluation of a large number of 2-D polynomials
  886. of the same degrees and sample points.
  887. Parameters
  888. ----------
  889. x, y : array_like
  890. Arrays of point coordinates, all of the same shape. The dtypes
  891. will be converted to either float64 or complex128 depending on
  892. whether any of the elements are complex. Scalars are converted to
  893. 1-D arrays.
  894. deg : list of ints
  895. List of maximum degrees of the form [x_deg, y_deg].
  896. Returns
  897. -------
  898. vander2d : ndarray
  899. The shape of the returned matrix is ``x.shape + (order,)``, where
  900. :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same
  901. as the converted `x` and `y`.
  902. See Also
  903. --------
  904. polyvander, polyvander3d. polyval2d, polyval3d
  905. """
  906. ideg = [int(d) for d in deg]
  907. is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
  908. if is_valid != [1, 1]:
  909. raise ValueError("degrees must be non-negative integers")
  910. degx, degy = ideg
  911. x, y = np.array((x, y), copy=0) + 0.0
  912. vx = polyvander(x, degx)
  913. vy = polyvander(y, degy)
  914. v = vx[..., None]*vy[..., None,:]
  915. # einsum bug
  916. #v = np.einsum("...i,...j->...ij", vx, vy)
  917. return v.reshape(v.shape[:-2] + (-1,))
  918. def polyvander3d(x, y, z, deg):
  919. """Pseudo-Vandermonde matrix of given degrees.
  920. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  921. points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
  922. then The pseudo-Vandermonde matrix is defined by
  923. .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k,
  924. where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
  925. indices of `V` index the points `(x, y, z)` and the last index encodes
  926. the powers of `x`, `y`, and `z`.
  927. If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
  928. of `V` correspond to the elements of a 3-D coefficient array `c` of
  929. shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
  930. .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
  931. and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the
  932. same up to roundoff. This equivalence is useful both for least squares
  933. fitting and for the evaluation of a large number of 3-D polynomials
  934. of the same degrees and sample points.
  935. Parameters
  936. ----------
  937. x, y, z : array_like
  938. Arrays of point coordinates, all of the same shape. The dtypes will
  939. be converted to either float64 or complex128 depending on whether
  940. any of the elements are complex. Scalars are converted to 1-D
  941. arrays.
  942. deg : list of ints
  943. List of maximum degrees of the form [x_deg, y_deg, z_deg].
  944. Returns
  945. -------
  946. vander3d : ndarray
  947. The shape of the returned matrix is ``x.shape + (order,)``, where
  948. :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
  949. be the same as the converted `x`, `y`, and `z`.
  950. See Also
  951. --------
  952. polyvander, polyvander3d. polyval2d, polyval3d
  953. Notes
  954. -----
  955. .. versionadded:: 1.7.0
  956. """
  957. ideg = [int(d) for d in deg]
  958. is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
  959. if is_valid != [1, 1, 1]:
  960. raise ValueError("degrees must be non-negative integers")
  961. degx, degy, degz = ideg
  962. x, y, z = np.array((x, y, z), copy=0) + 0.0
  963. vx = polyvander(x, degx)
  964. vy = polyvander(y, degy)
  965. vz = polyvander(z, degz)
  966. v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
  967. # einsum bug
  968. #v = np.einsum("...i, ...j, ...k->...ijk", vx, vy, vz)
  969. return v.reshape(v.shape[:-3] + (-1,))
  970. def polyfit(x, y, deg, rcond=None, full=False, w=None):
  971. """
  972. Least-squares fit of a polynomial to data.
  973. Return the coefficients of a polynomial of degree `deg` that is the
  974. least squares fit to the data values `y` given at points `x`. If `y` is
  975. 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
  976. fits are done, one for each column of `y`, and the resulting
  977. coefficients are stored in the corresponding columns of a 2-D return.
  978. The fitted polynomial(s) are in the form
  979. .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n,
  980. where `n` is `deg`.
  981. Parameters
  982. ----------
  983. x : array_like, shape (`M`,)
  984. x-coordinates of the `M` sample (data) points ``(x[i], y[i])``.
  985. y : array_like, shape (`M`,) or (`M`, `K`)
  986. y-coordinates of the sample points. Several sets of sample points
  987. sharing the same x-coordinates can be (independently) fit with one
  988. call to `polyfit` by passing in for `y` a 2-D array that contains
  989. one data set per column.
  990. deg : int or 1-D array_like
  991. Degree(s) of the fitting polynomials. If `deg` is a single integer
  992. all terms up to and including the `deg`'th term are included in the
  993. fit. For Numpy versions >= 1.11 a list of integers specifying the
  994. degrees of the terms to include may be used instead.
  995. rcond : float, optional
  996. Relative condition number of the fit. Singular values smaller
  997. than `rcond`, relative to the largest singular value, will be
  998. ignored. The default value is ``len(x)*eps``, where `eps` is the
  999. relative precision of the platform's float type, about 2e-16 in
  1000. most cases.
  1001. full : bool, optional
  1002. Switch determining the nature of the return value. When ``False``
  1003. (the default) just the coefficients are returned; when ``True``,
  1004. diagnostic information from the singular value decomposition (used
  1005. to solve the fit's matrix equation) is also returned.
  1006. w : array_like, shape (`M`,), optional
  1007. Weights. If not None, the contribution of each point
  1008. ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
  1009. weights are chosen so that the errors of the products ``w[i]*y[i]``
  1010. all have the same variance. The default value is None.
  1011. .. versionadded:: 1.5.0
  1012. Returns
  1013. -------
  1014. coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`)
  1015. Polynomial coefficients ordered from low to high. If `y` was 2-D,
  1016. the coefficients in column `k` of `coef` represent the polynomial
  1017. fit to the data in `y`'s `k`-th column.
  1018. [residuals, rank, singular_values, rcond] : list
  1019. These values are only returned if `full` = True
  1020. resid -- sum of squared residuals of the least squares fit
  1021. rank -- the numerical rank of the scaled Vandermonde matrix
  1022. sv -- singular values of the scaled Vandermonde matrix
  1023. rcond -- value of `rcond`.
  1024. For more details, see `linalg.lstsq`.
  1025. Raises
  1026. ------
  1027. RankWarning
  1028. Raised if the matrix in the least-squares fit is rank deficient.
  1029. The warning is only raised if `full` == False. The warnings can
  1030. be turned off by:
  1031. >>> import warnings
  1032. >>> warnings.simplefilter('ignore', RankWarning)
  1033. See Also
  1034. --------
  1035. chebfit, legfit, lagfit, hermfit, hermefit
  1036. polyval : Evaluates a polynomial.
  1037. polyvander : Vandermonde matrix for powers.
  1038. linalg.lstsq : Computes a least-squares fit from the matrix.
  1039. scipy.interpolate.UnivariateSpline : Computes spline fits.
  1040. Notes
  1041. -----
  1042. The solution is the coefficients of the polynomial `p` that minimizes
  1043. the sum of the weighted squared errors
  1044. .. math :: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
  1045. where the :math:`w_j` are the weights. This problem is solved by
  1046. setting up the (typically) over-determined matrix equation:
  1047. .. math :: V(x) * c = w * y,
  1048. where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
  1049. coefficients to be solved for, `w` are the weights, and `y` are the
  1050. observed values. This equation is then solved using the singular value
  1051. decomposition of `V`.
  1052. If some of the singular values of `V` are so small that they are
  1053. neglected (and `full` == ``False``), a `RankWarning` will be raised.
  1054. This means that the coefficient values may be poorly determined.
  1055. Fitting to a lower order polynomial will usually get rid of the warning
  1056. (but may not be what you want, of course; if you have independent
  1057. reason(s) for choosing the degree which isn't working, you may have to:
  1058. a) reconsider those reasons, and/or b) reconsider the quality of your
  1059. data). The `rcond` parameter can also be set to a value smaller than
  1060. its default, but the resulting fit may be spurious and have large
  1061. contributions from roundoff error.
  1062. Polynomial fits using double precision tend to "fail" at about
  1063. (polynomial) degree 20. Fits using Chebyshev or Legendre series are
  1064. generally better conditioned, but much can still depend on the
  1065. distribution of the sample points and the smoothness of the data. If
  1066. the quality of the fit is inadequate, splines may be a good
  1067. alternative.
  1068. Examples
  1069. --------
  1070. >>> from numpy.polynomial import polynomial as P
  1071. >>> x = np.linspace(-1,1,51) # x "data": [-1, -0.96, ..., 0.96, 1]
  1072. >>> y = x**3 - x + np.random.randn(len(x)) # x^3 - x + N(0,1) "noise"
  1073. >>> c, stats = P.polyfit(x,y,3,full=True)
  1074. >>> c # c[0], c[2] should be approx. 0, c[1] approx. -1, c[3] approx. 1
  1075. array([ 0.01909725, -1.30598256, -0.00577963, 1.02644286])
  1076. >>> stats # note the large SSR, explaining the rather poor results
  1077. [array([ 38.06116253]), 4, array([ 1.38446749, 1.32119158, 0.50443316,
  1078. 0.28853036]), 1.1324274851176597e-014]
  1079. Same thing without the added noise
  1080. >>> y = x**3 - x
  1081. >>> c, stats = P.polyfit(x,y,3,full=True)
  1082. >>> c # c[0], c[2] should be "very close to 0", c[1] ~= -1, c[3] ~= 1
  1083. array([ -1.73362882e-17, -1.00000000e+00, -2.67471909e-16,
  1084. 1.00000000e+00])
  1085. >>> stats # note the minuscule SSR
  1086. [array([ 7.46346754e-31]), 4, array([ 1.38446749, 1.32119158,
  1087. 0.50443316, 0.28853036]), 1.1324274851176597e-014]
  1088. """
  1089. x = np.asarray(x) + 0.0
  1090. y = np.asarray(y) + 0.0
  1091. deg = np.asarray(deg)
  1092. # check arguments.
  1093. if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
  1094. raise TypeError("deg must be an int or non-empty 1-D array of int")
  1095. if deg.min() < 0:
  1096. raise ValueError("expected deg >= 0")
  1097. if x.ndim != 1:
  1098. raise TypeError("expected 1D vector for x")
  1099. if x.size == 0:
  1100. raise TypeError("expected non-empty vector for x")
  1101. if y.ndim < 1 or y.ndim > 2:
  1102. raise TypeError("expected 1D or 2D array for y")
  1103. if len(x) != len(y):
  1104. raise TypeError("expected x and y to have same length")
  1105. if deg.ndim == 0:
  1106. lmax = deg
  1107. order = lmax + 1
  1108. van = polyvander(x, lmax)
  1109. else:
  1110. deg = np.sort(deg)
  1111. lmax = deg[-1]
  1112. order = len(deg)
  1113. van = polyvander(x, lmax)[:, deg]
  1114. # set up the least squares matrices in transposed form
  1115. lhs = van.T
  1116. rhs = y.T
  1117. if w is not None:
  1118. w = np.asarray(w) + 0.0
  1119. if w.ndim != 1:
  1120. raise TypeError("expected 1D vector for w")
  1121. if len(x) != len(w):
  1122. raise TypeError("expected x and w to have same length")
  1123. # apply weights. Don't use inplace operations as they
  1124. # can cause problems with NA.
  1125. lhs = lhs * w
  1126. rhs = rhs * w
  1127. # set rcond
  1128. if rcond is None:
  1129. rcond = len(x)*np.finfo(x.dtype).eps
  1130. # Determine the norms of the design matrix columns.
  1131. if issubclass(lhs.dtype.type, np.complexfloating):
  1132. scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
  1133. else:
  1134. scl = np.sqrt(np.square(lhs).sum(1))
  1135. scl[scl == 0] = 1
  1136. # Solve the least squares problem.
  1137. c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
  1138. c = (c.T/scl).T
  1139. # Expand c to include non-fitted coefficients which are set to zero
  1140. if deg.ndim == 1:
  1141. if c.ndim == 2:
  1142. cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
  1143. else:
  1144. cc = np.zeros(lmax + 1, dtype=c.dtype)
  1145. cc[deg] = c
  1146. c = cc
  1147. # warn on rank reduction
  1148. if rank != order and not full:
  1149. msg = "The fit may be poorly conditioned"
  1150. warnings.warn(msg, pu.RankWarning)
  1151. if full:
  1152. return c, [resids, rank, s, rcond]
  1153. else:
  1154. return c
  1155. def polycompanion(c):
  1156. """
  1157. Return the companion matrix of c.
  1158. The companion matrix for power series cannot be made symmetric by
  1159. scaling the basis, so this function differs from those for the
  1160. orthogonal polynomials.
  1161. Parameters
  1162. ----------
  1163. c : array_like
  1164. 1-D array of polynomial coefficients ordered from low to high
  1165. degree.
  1166. Returns
  1167. -------
  1168. mat : ndarray
  1169. Companion matrix of dimensions (deg, deg).
  1170. Notes
  1171. -----
  1172. .. versionadded:: 1.7.0
  1173. """
  1174. # c is a trimmed copy
  1175. [c] = pu.as_series([c])
  1176. if len(c) < 2:
  1177. raise ValueError('Series must have maximum degree of at least 1.')
  1178. if len(c) == 2:
  1179. return np.array([[-c[0]/c[1]]])
  1180. n = len(c) - 1
  1181. mat = np.zeros((n, n), dtype=c.dtype)
  1182. bot = mat.reshape(-1)[n::n+1]
  1183. bot[...] = 1
  1184. mat[:, -1] -= c[:-1]/c[-1]
  1185. return mat
  1186. def polyroots(c):
  1187. """
  1188. Compute the roots of a polynomial.
  1189. Return the roots (a.k.a. "zeros") of the polynomial
  1190. .. math:: p(x) = \\sum_i c[i] * x^i.
  1191. Parameters
  1192. ----------
  1193. c : 1-D array_like
  1194. 1-D array of polynomial coefficients.
  1195. Returns
  1196. -------
  1197. out : ndarray
  1198. Array of the roots of the polynomial. If all the roots are real,
  1199. then `out` is also real, otherwise it is complex.
  1200. See Also
  1201. --------
  1202. chebroots
  1203. Notes
  1204. -----
  1205. The root estimates are obtained as the eigenvalues of the companion
  1206. matrix, Roots far from the origin of the complex plane may have large
  1207. errors due to the numerical instability of the power series for such
  1208. values. Roots with multiplicity greater than 1 will also show larger
  1209. errors as the value of the series near such points is relatively
  1210. insensitive to errors in the roots. Isolated roots near the origin can
  1211. be improved by a few iterations of Newton's method.
  1212. Examples
  1213. --------
  1214. >>> import numpy.polynomial.polynomial as poly
  1215. >>> poly.polyroots(poly.polyfromroots((-1,0,1)))
  1216. array([-1., 0., 1.])
  1217. >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype
  1218. dtype('float64')
  1219. >>> j = complex(0,1)
  1220. >>> poly.polyroots(poly.polyfromroots((-j,0,j)))
  1221. array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j])
  1222. """
  1223. # c is a trimmed copy
  1224. [c] = pu.as_series([c])
  1225. if len(c) < 2:
  1226. return np.array([], dtype=c.dtype)
  1227. if len(c) == 2:
  1228. return np.array([-c[0]/c[1]])
  1229. m = polycompanion(c)
  1230. r = la.eigvals(m)
  1231. r.sort()
  1232. return r
  1233. #
  1234. # polynomial class
  1235. #
  1236. class Polynomial(ABCPolyBase):
  1237. """A power series class.
  1238. The Polynomial class provides the standard Python numerical methods
  1239. '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
  1240. attributes and methods listed in the `ABCPolyBase` documentation.
  1241. Parameters
  1242. ----------
  1243. coef : array_like
  1244. Polynomial coefficients in order of increasing degree, i.e.,
  1245. ``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``.
  1246. domain : (2,) array_like, optional
  1247. Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
  1248. to the interval ``[window[0], window[1]]`` by shifting and scaling.
  1249. The default value is [-1, 1].
  1250. window : (2,) array_like, optional
  1251. Window, see `domain` for its use. The default value is [-1, 1].
  1252. .. versionadded:: 1.6.0
  1253. """
  1254. # Virtual Functions
  1255. _add = staticmethod(polyadd)
  1256. _sub = staticmethod(polysub)
  1257. _mul = staticmethod(polymul)
  1258. _div = staticmethod(polydiv)
  1259. _pow = staticmethod(polypow)
  1260. _val = staticmethod(polyval)
  1261. _int = staticmethod(polyint)
  1262. _der = staticmethod(polyder)
  1263. _fit = staticmethod(polyfit)
  1264. _line = staticmethod(polyline)
  1265. _roots = staticmethod(polyroots)
  1266. _fromroots = staticmethod(polyfromroots)
  1267. # Virtual properties
  1268. nickname = 'poly'
  1269. domain = np.array(polydomain)
  1270. window = np.array(polydomain)