_polybase.py 38 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184
  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. """
  568. return len(self) - 1
  569. def cutdeg(self, deg):
  570. """Truncate series to the given degree.
  571. Reduce the degree of the series to `deg` by discarding the
  572. high order terms. If `deg` is greater than the current degree a
  573. copy of the current series is returned. This can be useful in least
  574. squares where the coefficients of the high degree terms may be very
  575. small.
  576. .. versionadded:: 1.5.0
  577. Parameters
  578. ----------
  579. deg : non-negative int
  580. The series is reduced to degree `deg` by discarding the high
  581. order terms. The value of `deg` must be a non-negative integer.
  582. Returns
  583. -------
  584. new_series : series
  585. New instance of series with reduced degree.
  586. """
  587. return self.truncate(deg + 1)
  588. def trim(self, tol=0):
  589. """Remove trailing coefficients
  590. Remove trailing coefficients until a coefficient is reached whose
  591. absolute value greater than `tol` or the beginning of the series is
  592. reached. If all the coefficients would be removed the series is set
  593. to ``[0]``. A new series instance is returned with the new
  594. coefficients. The current instance remains unchanged.
  595. Parameters
  596. ----------
  597. tol : non-negative number.
  598. All trailing coefficients less than `tol` will be removed.
  599. Returns
  600. -------
  601. new_series : series
  602. New instance of series with trimmed coefficients.
  603. """
  604. coef = pu.trimcoef(self.coef, tol)
  605. return self.__class__(coef, self.domain, self.window, self.symbol)
  606. def truncate(self, size):
  607. """Truncate series to length `size`.
  608. Reduce the series to length `size` by discarding the high
  609. degree terms. The value of `size` must be a positive integer. This
  610. can be useful in least squares where the coefficients of the
  611. high degree terms may be very small.
  612. Parameters
  613. ----------
  614. size : positive int
  615. The series is reduced to length `size` by discarding the high
  616. degree terms. The value of `size` must be a positive integer.
  617. Returns
  618. -------
  619. new_series : series
  620. New instance of series with truncated coefficients.
  621. """
  622. isize = int(size)
  623. if isize != size or isize < 1:
  624. raise ValueError("size must be a positive integer")
  625. if isize >= len(self.coef):
  626. coef = self.coef
  627. else:
  628. coef = self.coef[:isize]
  629. return self.__class__(coef, self.domain, self.window, self.symbol)
  630. def convert(self, domain=None, kind=None, window=None):
  631. """Convert series to a different kind and/or domain and/or window.
  632. Parameters
  633. ----------
  634. domain : array_like, optional
  635. The domain of the converted series. If the value is None,
  636. the default domain of `kind` is used.
  637. kind : class, optional
  638. The polynomial series type class to which the current instance
  639. should be converted. If kind is None, then the class of the
  640. current instance is used.
  641. window : array_like, optional
  642. The window of the converted series. If the value is None,
  643. the default window of `kind` is used.
  644. Returns
  645. -------
  646. new_series : series
  647. The returned class can be of different type than the current
  648. instance and/or have a different domain and/or different
  649. window.
  650. Notes
  651. -----
  652. Conversion between domains and class types can result in
  653. numerically ill defined series.
  654. """
  655. if kind is None:
  656. kind = self.__class__
  657. if domain is None:
  658. domain = kind.domain
  659. if window is None:
  660. window = kind.window
  661. return self(kind.identity(domain, window=window, symbol=self.symbol))
  662. def mapparms(self):
  663. """Return the mapping parameters.
  664. The returned values define a linear map ``off + scl*x`` that is
  665. applied to the input arguments before the series is evaluated. The
  666. map depends on the ``domain`` and ``window``; if the current
  667. ``domain`` is equal to the ``window`` the resulting map is the
  668. identity. If the coefficients of the series instance are to be
  669. used by themselves outside this class, then the linear function
  670. must be substituted for the ``x`` in the standard representation of
  671. the base polynomials.
  672. Returns
  673. -------
  674. off, scl : float or complex
  675. The mapping function is defined by ``off + scl*x``.
  676. Notes
  677. -----
  678. If the current domain is the interval ``[l1, r1]`` and the window
  679. is ``[l2, r2]``, then the linear mapping function ``L`` is
  680. defined by the equations::
  681. L(l1) = l2
  682. L(r1) = r2
  683. """
  684. return pu.mapparms(self.domain, self.window)
  685. def integ(self, m=1, k=[], lbnd=None):
  686. """Integrate.
  687. Return a series instance that is the definite integral of the
  688. current series.
  689. Parameters
  690. ----------
  691. m : non-negative int
  692. The number of integrations to perform.
  693. k : array_like
  694. Integration constants. The first constant is applied to the
  695. first integration, the second to the second, and so on. The
  696. list of values must less than or equal to `m` in length and any
  697. missing values are set to zero.
  698. lbnd : Scalar
  699. The lower bound of the definite integral.
  700. Returns
  701. -------
  702. new_series : series
  703. A new series representing the integral. The domain is the same
  704. as the domain of the integrated series.
  705. """
  706. off, scl = self.mapparms()
  707. if lbnd is None:
  708. lbnd = 0
  709. else:
  710. lbnd = off + scl*lbnd
  711. coef = self._int(self.coef, m, k, lbnd, 1./scl)
  712. return self.__class__(coef, self.domain, self.window, self.symbol)
  713. def deriv(self, m=1):
  714. """Differentiate.
  715. Return a series instance of that is the derivative of the current
  716. series.
  717. Parameters
  718. ----------
  719. m : non-negative int
  720. Find the derivative of order `m`.
  721. Returns
  722. -------
  723. new_series : series
  724. A new series representing the derivative. The domain is the same
  725. as the domain of the differentiated series.
  726. """
  727. off, scl = self.mapparms()
  728. coef = self._der(self.coef, m, scl)
  729. return self.__class__(coef, self.domain, self.window, self.symbol)
  730. def roots(self):
  731. """Return the roots of the series polynomial.
  732. Compute the roots for the series. Note that the accuracy of the
  733. roots decrease the further outside the domain they lie.
  734. Returns
  735. -------
  736. roots : ndarray
  737. Array containing the roots of the series.
  738. """
  739. roots = self._roots(self.coef)
  740. return pu.mapdomain(roots, self.window, self.domain)
  741. def linspace(self, n=100, domain=None):
  742. """Return x, y values at equally spaced points in domain.
  743. Returns the x, y values at `n` linearly spaced points across the
  744. domain. Here y is the value of the polynomial at the points x. By
  745. default the domain is the same as that of the series instance.
  746. This method is intended mostly as a plotting aid.
  747. .. versionadded:: 1.5.0
  748. Parameters
  749. ----------
  750. n : int, optional
  751. Number of point pairs to return. The default value is 100.
  752. domain : {None, array_like}, optional
  753. If not None, the specified domain is used instead of that of
  754. the calling instance. It should be of the form ``[beg,end]``.
  755. The default is None which case the class domain is used.
  756. Returns
  757. -------
  758. x, y : ndarray
  759. x is equal to linspace(self.domain[0], self.domain[1], n) and
  760. y is the series evaluated at element of x.
  761. """
  762. if domain is None:
  763. domain = self.domain
  764. x = np.linspace(domain[0], domain[1], n)
  765. y = self(x)
  766. return x, y
  767. @classmethod
  768. def fit(cls, x, y, deg, domain=None, rcond=None, full=False, w=None,
  769. window=None, symbol='x'):
  770. """Least squares fit to data.
  771. Return a series instance that is the least squares fit to the data
  772. `y` sampled at `x`. The domain of the returned instance can be
  773. specified and this will often result in a superior fit with less
  774. chance of ill conditioning.
  775. Parameters
  776. ----------
  777. x : array_like, shape (M,)
  778. x-coordinates of the M sample points ``(x[i], y[i])``.
  779. y : array_like, shape (M,)
  780. y-coordinates of the M sample points ``(x[i], y[i])``.
  781. deg : int or 1-D array_like
  782. Degree(s) of the fitting polynomials. If `deg` is a single integer
  783. all terms up to and including the `deg`'th term are included in the
  784. fit. For NumPy versions >= 1.11.0 a list of integers specifying the
  785. degrees of the terms to include may be used instead.
  786. domain : {None, [beg, end], []}, optional
  787. Domain to use for the returned series. If ``None``,
  788. then a minimal domain that covers the points `x` is chosen. If
  789. ``[]`` the class domain is used. The default value was the
  790. class domain in NumPy 1.4 and ``None`` in later versions.
  791. The ``[]`` option was added in numpy 1.5.0.
  792. rcond : float, optional
  793. Relative condition number of the fit. Singular values smaller
  794. than this relative to the largest singular value will be
  795. ignored. The default value is len(x)*eps, where eps is the
  796. relative precision of the float type, about 2e-16 in most
  797. cases.
  798. full : bool, optional
  799. Switch determining nature of return value. When it is False
  800. (the default) just the coefficients are returned, when True
  801. diagnostic information from the singular value decomposition is
  802. also returned.
  803. w : array_like, shape (M,), optional
  804. Weights. If not None, the weight ``w[i]`` applies to the unsquared
  805. residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
  806. chosen so that the errors of the products ``w[i]*y[i]`` all have
  807. the same variance. When using inverse-variance weighting, use
  808. ``w[i] = 1/sigma(y[i])``. The default value is None.
  809. .. versionadded:: 1.5.0
  810. window : {[beg, end]}, optional
  811. Window to use for the returned series. The default
  812. value is the default class domain
  813. .. versionadded:: 1.6.0
  814. symbol : str, optional
  815. Symbol representing the independent variable. Default is 'x'.
  816. Returns
  817. -------
  818. new_series : series
  819. A series that represents the least squares fit to the data and
  820. has the domain and window specified in the call. If the
  821. coefficients for the unscaled and unshifted basis polynomials are
  822. of interest, do ``new_series.convert().coef``.
  823. [resid, rank, sv, rcond] : list
  824. These values are only returned if ``full == True``
  825. - resid -- sum of squared residuals of the least squares fit
  826. - rank -- the numerical rank of the scaled Vandermonde matrix
  827. - sv -- singular values of the scaled Vandermonde matrix
  828. - rcond -- value of `rcond`.
  829. For more details, see `linalg.lstsq`.
  830. """
  831. if domain is None:
  832. domain = pu.getdomain(x)
  833. elif type(domain) is list and len(domain) == 0:
  834. domain = cls.domain
  835. if window is None:
  836. window = cls.window
  837. xnew = pu.mapdomain(x, domain, window)
  838. res = cls._fit(xnew, y, deg, w=w, rcond=rcond, full=full)
  839. if full:
  840. [coef, status] = res
  841. return (
  842. cls(coef, domain=domain, window=window, symbol=symbol), status
  843. )
  844. else:
  845. coef = res
  846. return cls(coef, domain=domain, window=window, symbol=symbol)
  847. @classmethod
  848. def fromroots(cls, roots, domain=[], window=None, symbol='x'):
  849. """Return series instance that has the specified roots.
  850. Returns a series representing the product
  851. ``(x - r[0])*(x - r[1])*...*(x - r[n-1])``, where ``r`` is a
  852. list of roots.
  853. Parameters
  854. ----------
  855. roots : array_like
  856. List of roots.
  857. domain : {[], None, array_like}, optional
  858. Domain for the resulting series. If None the domain is the
  859. interval from the smallest root to the largest. If [] the
  860. domain is the class domain. The default is [].
  861. window : {None, array_like}, optional
  862. Window for the returned series. If None the class window is
  863. used. The default is None.
  864. symbol : str, optional
  865. Symbol representing the independent variable. Default is 'x'.
  866. Returns
  867. -------
  868. new_series : series
  869. Series with the specified roots.
  870. """
  871. [roots] = pu.as_series([roots], trim=False)
  872. if domain is None:
  873. domain = pu.getdomain(roots)
  874. elif type(domain) is list and len(domain) == 0:
  875. domain = cls.domain
  876. if window is None:
  877. window = cls.window
  878. deg = len(roots)
  879. off, scl = pu.mapparms(domain, window)
  880. rnew = off + scl*roots
  881. coef = cls._fromroots(rnew) / scl**deg
  882. return cls(coef, domain=domain, window=window, symbol=symbol)
  883. @classmethod
  884. def identity(cls, domain=None, window=None, symbol='x'):
  885. """Identity function.
  886. If ``p`` is the returned series, then ``p(x) == x`` for all
  887. values of x.
  888. Parameters
  889. ----------
  890. domain : {None, array_like}, optional
  891. If given, the array must be of the form ``[beg, end]``, where
  892. ``beg`` and ``end`` are the endpoints of the domain. If None is
  893. given then the class domain is used. The default is None.
  894. window : {None, array_like}, optional
  895. If given, the resulting array must be if the form
  896. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  897. the window. If None is given then the class window is used. The
  898. default is None.
  899. symbol : str, optional
  900. Symbol representing the independent variable. Default is 'x'.
  901. Returns
  902. -------
  903. new_series : series
  904. Series of representing the identity.
  905. """
  906. if domain is None:
  907. domain = cls.domain
  908. if window is None:
  909. window = cls.window
  910. off, scl = pu.mapparms(window, domain)
  911. coef = cls._line(off, scl)
  912. return cls(coef, domain, window, symbol)
  913. @classmethod
  914. def basis(cls, deg, domain=None, window=None, symbol='x'):
  915. """Series basis polynomial of degree `deg`.
  916. Returns the series representing the basis polynomial of degree `deg`.
  917. .. versionadded:: 1.7.0
  918. Parameters
  919. ----------
  920. deg : int
  921. Degree of the basis polynomial for the series. Must be >= 0.
  922. domain : {None, array_like}, optional
  923. If given, the array must be of the form ``[beg, end]``, where
  924. ``beg`` and ``end`` are the endpoints of the domain. If None is
  925. given then the class domain is used. The default is None.
  926. window : {None, array_like}, optional
  927. If given, the resulting array must be if the form
  928. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  929. the window. If None is given then the class window is used. The
  930. default is None.
  931. symbol : str, optional
  932. Symbol representing the independent variable. Default is 'x'.
  933. Returns
  934. -------
  935. new_series : series
  936. A series with the coefficient of the `deg` term set to one and
  937. all others zero.
  938. """
  939. if domain is None:
  940. domain = cls.domain
  941. if window is None:
  942. window = cls.window
  943. ideg = int(deg)
  944. if ideg != deg or ideg < 0:
  945. raise ValueError("deg must be non-negative integer")
  946. return cls([0]*ideg + [1], domain, window, symbol)
  947. @classmethod
  948. def cast(cls, series, domain=None, window=None):
  949. """Convert series to series of this class.
  950. The `series` is expected to be an instance of some polynomial
  951. series of one of the types supported by by the numpy.polynomial
  952. module, but could be some other class that supports the convert
  953. method.
  954. .. versionadded:: 1.7.0
  955. Parameters
  956. ----------
  957. series : series
  958. The series instance to be converted.
  959. domain : {None, array_like}, optional
  960. If given, the array must be of the form ``[beg, end]``, where
  961. ``beg`` and ``end`` are the endpoints of the domain. If None is
  962. given then the class domain is used. The default is None.
  963. window : {None, array_like}, optional
  964. If given, the resulting array must be if the form
  965. ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
  966. the window. If None is given then the class window is used. The
  967. default is None.
  968. Returns
  969. -------
  970. new_series : series
  971. A series of the same kind as the calling class and equal to
  972. `series` when evaluated.
  973. See Also
  974. --------
  975. convert : similar instance method
  976. """
  977. if domain is None:
  978. domain = cls.domain
  979. if window is None:
  980. window = cls.window
  981. return series.convert(domain, cls, window)