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