_polybase.py 32 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041
  1. """
  2. Abstract base class for the various polynomial Classes.
  3. The ABCPolyBase class provides the methods needed to implement the common API
  4. for the various polynomial classes. It operates as a mixin, but uses the
  5. abc module from the stdlib, hence it is only available for Python >= 2.6.
  6. """
  7. import abc
  8. import numbers
  9. import numpy as np
  10. from . import polyutils as pu
  11. __all__ = ['ABCPolyBase']
  12. class ABCPolyBase(abc.ABC):
  13. """An abstract base class for immutable series classes.
  14. ABCPolyBase provides the standard Python numerical methods
  15. '+', '-', '*', '//', '%', 'divmod', '**', and '()' along with the
  16. methods listed below.
  17. .. versionadded:: 1.9.0
  18. Parameters
  19. ----------
  20. coef : array_like
  21. Series coefficients in order of increasing degree, i.e.,
  22. ``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``, where
  23. ``P_i`` is the basis polynomials of degree ``i``.
  24. domain : (2,) array_like, optional
  25. Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
  26. to the interval ``[window[0], window[1]]`` by shifting and scaling.
  27. The default value is the derived class domain.
  28. window : (2,) array_like, optional
  29. Window, see domain for its use. The default value is the
  30. derived class window.
  31. Attributes
  32. ----------
  33. coef : (N,) ndarray
  34. Series coefficients in order of increasing degree.
  35. domain : (2,) ndarray
  36. Domain that is mapped to window.
  37. window : (2,) ndarray
  38. Window that domain is mapped to.
  39. Class Attributes
  40. ----------------
  41. maxpower : int
  42. Maximum power allowed, i.e., the largest number ``n`` such that
  43. ``p(x)**n`` is allowed. This is to limit runaway polynomial size.
  44. domain : (2,) ndarray
  45. Default domain of the class.
  46. window : (2,) ndarray
  47. Default window of the class.
  48. """
  49. # Not hashable
  50. __hash__ = None
  51. # Opt out of numpy ufuncs and Python ops with ndarray subclasses.
  52. __array_ufunc__ = None
  53. # Limit runaway size. T_n^m has degree n*m
  54. maxpower = 100
  55. @property
  56. @abc.abstractmethod
  57. def domain(self):
  58. pass
  59. @property
  60. @abc.abstractmethod
  61. def window(self):
  62. pass
  63. @property
  64. @abc.abstractmethod
  65. def nickname(self):
  66. pass
  67. @property
  68. @abc.abstractmethod
  69. def basis_name(self):
  70. pass
  71. @staticmethod
  72. @abc.abstractmethod
  73. def _add(c1, c2):
  74. pass
  75. @staticmethod
  76. @abc.abstractmethod
  77. def _sub(c1, c2):
  78. pass
  79. @staticmethod
  80. @abc.abstractmethod
  81. def _mul(c1, c2):
  82. pass
  83. @staticmethod
  84. @abc.abstractmethod
  85. def _div(c1, c2):
  86. pass
  87. @staticmethod
  88. @abc.abstractmethod
  89. def _pow(c, pow, maxpower=None):
  90. pass
  91. @staticmethod
  92. @abc.abstractmethod
  93. def _val(x, c):
  94. pass
  95. @staticmethod
  96. @abc.abstractmethod
  97. def _int(c, m, k, lbnd, scl):
  98. pass
  99. @staticmethod
  100. @abc.abstractmethod
  101. def _der(c, m, scl):
  102. pass
  103. @staticmethod
  104. @abc.abstractmethod
  105. def _fit(x, y, deg, rcond, full):
  106. pass
  107. @staticmethod
  108. @abc.abstractmethod
  109. def _line(off, scl):
  110. pass
  111. @staticmethod
  112. @abc.abstractmethod
  113. def _roots(c):
  114. pass
  115. @staticmethod
  116. @abc.abstractmethod
  117. def _fromroots(r):
  118. pass
  119. def has_samecoef(self, other):
  120. """Check if coefficients match.
  121. .. versionadded:: 1.6.0
  122. Parameters
  123. ----------
  124. other : class instance
  125. The other class must have the ``coef`` attribute.
  126. Returns
  127. -------
  128. bool : boolean
  129. True if the coefficients are the same, False otherwise.
  130. """
  131. if len(self.coef) != len(other.coef):
  132. return False
  133. elif not np.all(self.coef == other.coef):
  134. return False
  135. else:
  136. return True
  137. def has_samedomain(self, other):
  138. """Check if domains match.
  139. .. versionadded:: 1.6.0
  140. Parameters
  141. ----------
  142. other : class instance
  143. The other class must have the ``domain`` attribute.
  144. Returns
  145. -------
  146. bool : boolean
  147. True if the domains are the same, False otherwise.
  148. """
  149. return np.all(self.domain == other.domain)
  150. def has_samewindow(self, other):
  151. """Check if windows match.
  152. .. versionadded:: 1.6.0
  153. Parameters
  154. ----------
  155. other : class instance
  156. The other class must have the ``window`` attribute.
  157. Returns
  158. -------
  159. bool : boolean
  160. True if the windows are the same, False otherwise.
  161. """
  162. return np.all(self.window == other.window)
  163. def has_sametype(self, other):
  164. """Check if types match.
  165. .. versionadded:: 1.7.0
  166. Parameters
  167. ----------
  168. other : object
  169. Class instance.
  170. Returns
  171. -------
  172. bool : boolean
  173. True if other is same class as self
  174. """
  175. return isinstance(other, self.__class__)
  176. def _get_coefficients(self, other):
  177. """Interpret other as polynomial coefficients.
  178. The `other` argument is checked to see if it is of the same
  179. class as self with identical domain and window. If so,
  180. return its coefficients, otherwise return `other`.
  181. .. versionadded:: 1.9.0
  182. Parameters
  183. ----------
  184. other : anything
  185. Object to be checked.
  186. Returns
  187. -------
  188. coef
  189. The coefficients of`other` if it is a compatible instance,
  190. of ABCPolyBase, otherwise `other`.
  191. Raises
  192. ------
  193. TypeError
  194. When `other` is an incompatible instance of ABCPolyBase.
  195. """
  196. if isinstance(other, ABCPolyBase):
  197. if not isinstance(other, self.__class__):
  198. raise TypeError("Polynomial types differ")
  199. elif not np.all(self.domain == other.domain):
  200. raise TypeError("Domains differ")
  201. elif not np.all(self.window == other.window):
  202. raise TypeError("Windows differ")
  203. return other.coef
  204. return other
  205. def __init__(self, coef, domain=None, window=None):
  206. [coef] = pu.as_series([coef], trim=False)
  207. self.coef = coef
  208. if domain is not None:
  209. [domain] = pu.as_series([domain], trim=False)
  210. if len(domain) != 2:
  211. raise ValueError("Domain has wrong number of elements.")
  212. self.domain = domain
  213. if window is not None:
  214. [window] = pu.as_series([window], trim=False)
  215. if len(window) != 2:
  216. raise ValueError("Window has wrong number of elements.")
  217. self.window = window
  218. def __repr__(self):
  219. coef = repr(self.coef)[6:-1]
  220. domain = repr(self.domain)[6:-1]
  221. window = repr(self.window)[6:-1]
  222. name = self.__class__.__name__
  223. return f"{name}({coef}, domain={domain}, window={window})"
  224. def __str__(self):
  225. coef = str(self.coef)
  226. name = self.nickname
  227. return f"{name}({coef})"
  228. @classmethod
  229. def _repr_latex_term(cls, i, arg_str, needs_parens):
  230. if cls.basis_name is None:
  231. raise NotImplementedError(
  232. "Subclasses must define either a basis name, or override "
  233. "_repr_latex_term(i, arg_str, needs_parens)")
  234. # since we always add parens, we don't care if the expression needs them
  235. return f"{{{cls.basis_name}}}_{{{i}}}({arg_str})"
  236. @staticmethod
  237. def _repr_latex_scalar(x):
  238. # TODO: we're stuck with disabling math formatting until we handle
  239. # exponents in this function
  240. return r'\text{{{}}}'.format(x)
  241. def _repr_latex_(self):
  242. # get the scaled argument string to the basis functions
  243. off, scale = self.mapparms()
  244. if off == 0 and scale == 1:
  245. term = 'x'
  246. needs_parens = False
  247. elif scale == 1:
  248. term = f"{self._repr_latex_scalar(off)} + x"
  249. needs_parens = True
  250. elif off == 0:
  251. term = f"{self._repr_latex_scalar(scale)}x"
  252. needs_parens = True
  253. else:
  254. term = (
  255. f"{self._repr_latex_scalar(off)} + "
  256. f"{self._repr_latex_scalar(scale)}x"
  257. )
  258. needs_parens = True
  259. mute = r"\color{{LightGray}}{{{}}}".format
  260. parts = []
  261. for i, c in enumerate(self.coef):
  262. # prevent duplication of + and - signs
  263. if i == 0:
  264. coef_str = f"{self._repr_latex_scalar(c)}"
  265. elif not isinstance(c, numbers.Real):
  266. coef_str = f" + ({self._repr_latex_scalar(c)})"
  267. elif not np.signbit(c):
  268. coef_str = f" + {self._repr_latex_scalar(c)}"
  269. else:
  270. coef_str = f" - {self._repr_latex_scalar(-c)}"
  271. # produce the string for the term
  272. term_str = self._repr_latex_term(i, term, needs_parens)
  273. if term_str == '1':
  274. part = coef_str
  275. else:
  276. part = rf"{coef_str}\,{term_str}"
  277. if c == 0:
  278. part = mute(part)
  279. parts.append(part)
  280. if parts:
  281. body = ''.join(parts)
  282. else:
  283. # in case somehow there are no coefficients at all
  284. body = '0'
  285. return rf"$x \mapsto {body}$"
  286. # Pickle and copy
  287. def __getstate__(self):
  288. ret = self.__dict__.copy()
  289. ret['coef'] = self.coef.copy()
  290. ret['domain'] = self.domain.copy()
  291. ret['window'] = self.window.copy()
  292. return ret
  293. def __setstate__(self, dict):
  294. self.__dict__ = dict
  295. # Call
  296. def __call__(self, arg):
  297. off, scl = pu.mapparms(self.domain, self.window)
  298. arg = off + scl*arg
  299. return self._val(arg, self.coef)
  300. def __iter__(self):
  301. return iter(self.coef)
  302. def __len__(self):
  303. return len(self.coef)
  304. # Numeric properties.
  305. def __neg__(self):
  306. return self.__class__(-self.coef, self.domain, self.window)
  307. def __pos__(self):
  308. return self
  309. def __add__(self, other):
  310. othercoef = self._get_coefficients(other)
  311. try:
  312. coef = self._add(self.coef, othercoef)
  313. except Exception:
  314. return NotImplemented
  315. return self.__class__(coef, self.domain, self.window)
  316. def __sub__(self, other):
  317. othercoef = self._get_coefficients(other)
  318. try:
  319. coef = self._sub(self.coef, othercoef)
  320. except Exception:
  321. return NotImplemented
  322. return self.__class__(coef, self.domain, self.window)
  323. def __mul__(self, other):
  324. othercoef = self._get_coefficients(other)
  325. try:
  326. coef = self._mul(self.coef, othercoef)
  327. except Exception:
  328. return NotImplemented
  329. return self.__class__(coef, self.domain, self.window)
  330. def __truediv__(self, other):
  331. # there is no true divide if the rhs is not a Number, although it
  332. # could return the first n elements of an infinite series.
  333. # It is hard to see where n would come from, though.
  334. if not isinstance(other, numbers.Number) or isinstance(other, bool):
  335. raise TypeError(
  336. f"unsupported types for true division: "
  337. f"'{type(self)}', '{type(other)}'"
  338. )
  339. return self.__floordiv__(other)
  340. def __floordiv__(self, other):
  341. res = self.__divmod__(other)
  342. if res is NotImplemented:
  343. return res
  344. return res[0]
  345. def __mod__(self, other):
  346. res = self.__divmod__(other)
  347. if res is NotImplemented:
  348. return res
  349. return res[1]
  350. def __divmod__(self, other):
  351. othercoef = self._get_coefficients(other)
  352. try:
  353. quo, rem = self._div(self.coef, othercoef)
  354. except ZeroDivisionError as e:
  355. raise e
  356. except Exception:
  357. return NotImplemented
  358. quo = self.__class__(quo, self.domain, self.window)
  359. rem = self.__class__(rem, self.domain, self.window)
  360. return quo, rem
  361. def __pow__(self, other):
  362. coef = self._pow(self.coef, other, maxpower=self.maxpower)
  363. res = self.__class__(coef, self.domain, self.window)
  364. return res
  365. def __radd__(self, other):
  366. try:
  367. coef = self._add(other, self.coef)
  368. except Exception:
  369. return NotImplemented
  370. return self.__class__(coef, self.domain, self.window)
  371. def __rsub__(self, other):
  372. try:
  373. coef = self._sub(other, self.coef)
  374. except Exception:
  375. return NotImplemented
  376. return self.__class__(coef, self.domain, self.window)
  377. def __rmul__(self, other):
  378. try:
  379. coef = self._mul(other, self.coef)
  380. except Exception:
  381. return NotImplemented
  382. return self.__class__(coef, self.domain, self.window)
  383. def __rdiv__(self, other):
  384. # set to __floordiv__ /.
  385. return self.__rfloordiv__(other)
  386. def __rtruediv__(self, other):
  387. # An instance of ABCPolyBase is not considered a
  388. # Number.
  389. return NotImplemented
  390. def __rfloordiv__(self, other):
  391. res = self.__rdivmod__(other)
  392. if res is NotImplemented:
  393. return res
  394. return res[0]
  395. def __rmod__(self, other):
  396. res = self.__rdivmod__(other)
  397. if res is NotImplemented:
  398. return res
  399. return res[1]
  400. def __rdivmod__(self, other):
  401. try:
  402. quo, rem = self._div(other, self.coef)
  403. except ZeroDivisionError as e:
  404. raise e
  405. except Exception:
  406. return NotImplemented
  407. quo = self.__class__(quo, self.domain, self.window)
  408. rem = self.__class__(rem, self.domain, self.window)
  409. return quo, rem
  410. def __eq__(self, other):
  411. res = (isinstance(other, self.__class__) and
  412. np.all(self.domain == other.domain) and
  413. np.all(self.window == other.window) and
  414. (self.coef.shape == other.coef.shape) and
  415. np.all(self.coef == other.coef))
  416. return res
  417. def __ne__(self, other):
  418. return not self.__eq__(other)
  419. #
  420. # Extra methods.
  421. #
  422. def copy(self):
  423. """Return a copy.
  424. Returns
  425. -------
  426. new_series : series
  427. Copy of self.
  428. """
  429. return self.__class__(self.coef, self.domain, self.window)
  430. def degree(self):
  431. """The degree of the series.
  432. .. versionadded:: 1.5.0
  433. Returns
  434. -------
  435. degree : int
  436. Degree of the series, one less than the number of coefficients.
  437. """
  438. return len(self) - 1
  439. def cutdeg(self, deg):
  440. """Truncate series to the given degree.
  441. Reduce the degree of the series to `deg` by discarding the
  442. high order terms. If `deg` is greater than the current degree a
  443. copy of the current series is returned. This can be useful in least
  444. squares where the coefficients of the high degree terms may be very
  445. small.
  446. .. versionadded:: 1.5.0
  447. Parameters
  448. ----------
  449. deg : non-negative int
  450. The series is reduced to degree `deg` by discarding the high
  451. order terms. The value of `deg` must be a non-negative integer.
  452. Returns
  453. -------
  454. new_series : series
  455. New instance of series with reduced degree.
  456. """
  457. return self.truncate(deg + 1)
  458. def trim(self, tol=0):
  459. """Remove trailing coefficients
  460. Remove trailing coefficients until a coefficient is reached whose
  461. absolute value greater than `tol` or the beginning of the series is
  462. reached. If all the coefficients would be removed the series is set
  463. to ``[0]``. A new series instance is returned with the new
  464. coefficients. The current instance remains unchanged.
  465. Parameters
  466. ----------
  467. tol : non-negative number.
  468. All trailing coefficients less than `tol` will be removed.
  469. Returns
  470. -------
  471. new_series : series
  472. Contains the new set of coefficients.
  473. """
  474. coef = pu.trimcoef(self.coef, tol)
  475. return self.__class__(coef, self.domain, self.window)
  476. def truncate(self, size):
  477. """Truncate series to length `size`.
  478. Reduce the series to length `size` by discarding the high
  479. degree terms. The value of `size` must be a positive integer. This
  480. can be useful in least squares where the coefficients of the
  481. high degree terms may be very small.
  482. Parameters
  483. ----------
  484. size : positive int
  485. The series is reduced to length `size` by discarding the high
  486. degree terms. The value of `size` must be a positive integer.
  487. Returns
  488. -------
  489. new_series : series
  490. New instance of series with truncated coefficients.
  491. """
  492. isize = int(size)
  493. if isize != size or isize < 1:
  494. raise ValueError("size must be a positive integer")
  495. if isize >= len(self.coef):
  496. coef = self.coef
  497. else:
  498. coef = self.coef[:isize]
  499. return self.__class__(coef, self.domain, self.window)
  500. def convert(self, domain=None, kind=None, window=None):
  501. """Convert series to a different kind and/or domain and/or window.
  502. Parameters
  503. ----------
  504. domain : array_like, optional
  505. The domain of the converted series. If the value is None,
  506. the default domain of `kind` is used.
  507. kind : class, optional
  508. The polynomial series type class to which the current instance
  509. should be converted. If kind is None, then the class of the
  510. current instance is used.
  511. window : array_like, optional
  512. The window of the converted series. If the value is None,
  513. the default window of `kind` is used.
  514. Returns
  515. -------
  516. new_series : series
  517. The returned class can be of different type than the current
  518. instance and/or have a different domain and/or different
  519. window.
  520. Notes
  521. -----
  522. Conversion between domains and class types can result in
  523. numerically ill defined series.
  524. Examples
  525. --------
  526. """
  527. if kind is None:
  528. kind = self.__class__
  529. if domain is None:
  530. domain = kind.domain
  531. if window is None:
  532. window = kind.window
  533. return self(kind.identity(domain, window=window))
  534. def mapparms(self):
  535. """Return the mapping parameters.
  536. The returned values define a linear map ``off + scl*x`` that is
  537. applied to the input arguments before the series is evaluated. The
  538. map depends on the ``domain`` and ``window``; if the current
  539. ``domain`` is equal to the ``window`` the resulting map is the
  540. identity. If the coefficients of the series instance are to be
  541. used by themselves outside this class, then the linear function
  542. must be substituted for the ``x`` in the standard representation of
  543. the base polynomials.
  544. Returns
  545. -------
  546. off, scl : float or complex
  547. The mapping function is defined by ``off + scl*x``.
  548. Notes
  549. -----
  550. If the current domain is the interval ``[l1, r1]`` and the window
  551. is ``[l2, r2]``, then the linear mapping function ``L`` is
  552. defined by the equations::
  553. L(l1) = l2
  554. L(r1) = r2
  555. """
  556. return pu.mapparms(self.domain, self.window)
  557. def integ(self, m=1, k=[], lbnd=None):
  558. """Integrate.
  559. Return a series instance that is the definite integral of the
  560. current series.
  561. Parameters
  562. ----------
  563. m : non-negative int
  564. The number of integrations to perform.
  565. k : array_like
  566. Integration constants. The first constant is applied to the
  567. first integration, the second to the second, and so on. The
  568. list of values must less than or equal to `m` in length and any
  569. missing values are set to zero.
  570. lbnd : Scalar
  571. The lower bound of the definite integral.
  572. Returns
  573. -------
  574. new_series : series
  575. A new series representing the integral. The domain is the same
  576. as the domain of the integrated series.
  577. """
  578. off, scl = self.mapparms()
  579. if lbnd is None:
  580. lbnd = 0
  581. else:
  582. lbnd = off + scl*lbnd
  583. coef = self._int(self.coef, m, k, lbnd, 1./scl)
  584. return self.__class__(coef, self.domain, self.window)
  585. def deriv(self, m=1):
  586. """Differentiate.
  587. Return a series instance of that is the derivative of the current
  588. series.
  589. Parameters
  590. ----------
  591. m : non-negative int
  592. Find the derivative of order `m`.
  593. Returns
  594. -------
  595. new_series : series
  596. A new series representing the derivative. The domain is the same
  597. as the domain of the differentiated series.
  598. """
  599. off, scl = self.mapparms()
  600. coef = self._der(self.coef, m, scl)
  601. return self.__class__(coef, self.domain, self.window)
  602. def roots(self):
  603. """Return the roots of the series polynomial.
  604. Compute the roots for the series. Note that the accuracy of the
  605. roots decrease the further outside the domain they lie.
  606. Returns
  607. -------
  608. roots : ndarray
  609. Array containing the roots of the series.
  610. """
  611. roots = self._roots(self.coef)
  612. return pu.mapdomain(roots, self.window, self.domain)
  613. def linspace(self, n=100, domain=None):
  614. """Return x, y values at equally spaced points in domain.
  615. Returns the x, y values at `n` linearly spaced points across the
  616. domain. Here y is the value of the polynomial at the points x. By
  617. default the domain is the same as that of the series instance.
  618. This method is intended mostly as a plotting aid.
  619. .. versionadded:: 1.5.0
  620. Parameters
  621. ----------
  622. n : int, optional
  623. Number of point pairs to return. The default value is 100.
  624. domain : {None, array_like}, optional
  625. If not None, the specified domain is used instead of that of
  626. the calling instance. It should be of the form ``[beg,end]``.
  627. The default is None which case the class domain is used.
  628. Returns
  629. -------
  630. x, y : ndarray
  631. x is equal to linspace(self.domain[0], self.domain[1], n) and
  632. y is the series evaluated at element of x.
  633. """
  634. if domain is None:
  635. domain = self.domain
  636. x = np.linspace(domain[0], domain[1], n)
  637. y = self(x)
  638. return x, y
  639. @classmethod
  640. def fit(cls, x, y, deg, domain=None, rcond=None, full=False, w=None,
  641. window=None):
  642. """Least squares fit to data.
  643. Return a series instance that is the least squares fit to the data
  644. `y` sampled at `x`. The domain of the returned instance can be
  645. specified and this will often result in a superior fit with less
  646. chance of ill conditioning.
  647. Parameters
  648. ----------
  649. x : array_like, shape (M,)
  650. x-coordinates of the M sample points ``(x[i], y[i])``.
  651. y : array_like, shape (M,) or (M, K)
  652. y-coordinates of the sample points. Several data sets of sample
  653. points sharing the same x-coordinates can be fitted at once by
  654. passing in a 2D-array that contains one dataset per column.
  655. deg : int or 1-D array_like
  656. Degree(s) of the fitting polynomials. If `deg` is a single integer
  657. all terms up to and including the `deg`'th term are included in the
  658. fit. For NumPy versions >= 1.11.0 a list of integers specifying the
  659. degrees of the terms to include may be used instead.
  660. domain : {None, [beg, end], []}, optional
  661. Domain to use for the returned series. If ``None``,
  662. then a minimal domain that covers the points `x` is chosen. If
  663. ``[]`` the class domain is used. The default value was the
  664. class domain in NumPy 1.4 and ``None`` in later versions.
  665. The ``[]`` option was added in numpy 1.5.0.
  666. rcond : float, optional
  667. Relative condition number of the fit. Singular values smaller
  668. than this relative to the largest singular value will be
  669. ignored. The default value is len(x)*eps, where eps is the
  670. relative precision of the float type, about 2e-16 in most
  671. cases.
  672. full : bool, optional
  673. Switch determining nature of return value. When it is False
  674. (the default) just the coefficients are returned, when True
  675. diagnostic information from the singular value decomposition is
  676. also returned.
  677. w : array_like, shape (M,), optional
  678. Weights. If not None the contribution of each point
  679. ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
  680. weights are chosen so that the errors of the products
  681. ``w[i]*y[i]`` all have the same variance. The default value is
  682. None.
  683. .. versionadded:: 1.5.0
  684. window : {[beg, end]}, optional
  685. Window to use for the returned series. The default
  686. value is the default class domain
  687. .. versionadded:: 1.6.0
  688. Returns
  689. -------
  690. new_series : series
  691. A series that represents the least squares fit to the data and
  692. has the domain and window specified in the call. If the
  693. coefficients for the unscaled and unshifted basis polynomials are
  694. of interest, do ``new_series.convert().coef``.
  695. [resid, rank, sv, rcond] : list
  696. These values are only returned if `full` = True
  697. resid -- sum of squared residuals of the least squares fit
  698. rank -- the numerical rank of the scaled Vandermonde matrix
  699. sv -- singular values of the scaled Vandermonde matrix
  700. rcond -- value of `rcond`.
  701. For more details, see `linalg.lstsq`.
  702. """
  703. if domain is None:
  704. domain = pu.getdomain(x)
  705. elif type(domain) is list and len(domain) == 0:
  706. domain = cls.domain
  707. if window is None:
  708. window = cls.window
  709. xnew = pu.mapdomain(x, domain, window)
  710. res = cls._fit(xnew, y, deg, w=w, rcond=rcond, full=full)
  711. if full:
  712. [coef, status] = res
  713. return cls(coef, domain=domain, window=window), status
  714. else:
  715. coef = res
  716. return cls(coef, domain=domain, window=window)
  717. @classmethod
  718. def fromroots(cls, roots, domain=[], window=None):
  719. """Return series instance that has the specified roots.
  720. Returns a series representing the product
  721. ``(x - r[0])*(x - r[1])*...*(x - r[n-1])``, where ``r`` is a
  722. list of roots.
  723. Parameters
  724. ----------
  725. roots : array_like
  726. List of roots.
  727. domain : {[], None, array_like}, optional
  728. Domain for the resulting series. If None the domain is the
  729. interval from the smallest root to the largest. If [] the
  730. domain is the class domain. The default is [].
  731. window : {None, array_like}, optional
  732. Window for the returned series. If None the class window is
  733. used. The default is None.
  734. Returns
  735. -------
  736. new_series : series
  737. Series with the specified roots.
  738. """
  739. [roots] = pu.as_series([roots], trim=False)
  740. if domain is None:
  741. domain = pu.getdomain(roots)
  742. elif type(domain) is list and len(domain) == 0:
  743. domain = cls.domain
  744. if window is None:
  745. window = cls.window
  746. deg = len(roots)
  747. off, scl = pu.mapparms(domain, window)
  748. rnew = off + scl*roots
  749. coef = cls._fromroots(rnew) / scl**deg
  750. return cls(coef, domain=domain, window=window)
  751. @classmethod
  752. def identity(cls, domain=None, window=None):
  753. """Identity function.
  754. If ``p`` is the returned series, then ``p(x) == x`` for all
  755. values of x.
  756. Parameters
  757. ----------
  758. domain : {None, array_like}, optional
  759. If given, the array must be of the form ``[beg, end]``, where
  760. ``beg`` and ``end`` are the endpoints of the domain. If None is
  761. given then the class domain is used. The default is None.
  762. window : {None, array_like}, optional
  763. If given, the resulting array must be if the form
  764. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  765. the window. If None is given then the class window is used. The
  766. default is None.
  767. Returns
  768. -------
  769. new_series : series
  770. Series of representing the identity.
  771. """
  772. if domain is None:
  773. domain = cls.domain
  774. if window is None:
  775. window = cls.window
  776. off, scl = pu.mapparms(window, domain)
  777. coef = cls._line(off, scl)
  778. return cls(coef, domain, window)
  779. @classmethod
  780. def basis(cls, deg, domain=None, window=None):
  781. """Series basis polynomial of degree `deg`.
  782. Returns the series representing the basis polynomial of degree `deg`.
  783. .. versionadded:: 1.7.0
  784. Parameters
  785. ----------
  786. deg : int
  787. Degree of the basis polynomial for the series. Must be >= 0.
  788. domain : {None, array_like}, optional
  789. If given, the array must be of the form ``[beg, end]``, where
  790. ``beg`` and ``end`` are the endpoints of the domain. If None is
  791. given then the class domain is used. The default is None.
  792. window : {None, array_like}, optional
  793. If given, the resulting array must be if the form
  794. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  795. the window. If None is given then the class window is used. The
  796. default is None.
  797. Returns
  798. -------
  799. new_series : series
  800. A series with the coefficient of the `deg` term set to one and
  801. all others zero.
  802. """
  803. if domain is None:
  804. domain = cls.domain
  805. if window is None:
  806. window = cls.window
  807. ideg = int(deg)
  808. if ideg != deg or ideg < 0:
  809. raise ValueError("deg must be non-negative integer")
  810. return cls([0]*ideg + [1], domain, window)
  811. @classmethod
  812. def cast(cls, series, domain=None, window=None):
  813. """Convert series to series of this class.
  814. The `series` is expected to be an instance of some polynomial
  815. series of one of the types supported by by the numpy.polynomial
  816. module, but could be some other class that supports the convert
  817. method.
  818. .. versionadded:: 1.7.0
  819. Parameters
  820. ----------
  821. series : series
  822. The series instance to be converted.
  823. domain : {None, array_like}, optional
  824. If given, the array must be of the form ``[beg, end]``, where
  825. ``beg`` and ``end`` are the endpoints of the domain. If None is
  826. given then the class domain is used. The default is None.
  827. window : {None, array_like}, optional
  828. If given, the resulting array must be if the form
  829. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  830. the window. If None is given then the class window is used. The
  831. default is None.
  832. Returns
  833. -------
  834. new_series : series
  835. A series of the same kind as the calling class and equal to
  836. `series` when evaluated.
  837. See Also
  838. --------
  839. convert : similar instance method
  840. """
  841. if domain is None:
  842. domain = cls.domain
  843. if window is None:
  844. window = cls.window
  845. return series.convert(domain, cls, window)