_polybase.py 36 KB

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