_polybase.py 38 KB

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