_polybase.py 38 KB

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