_polybase.py 32 KB

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