123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413 |
- """
- Utility classes and functions for the polynomial modules.
- This module provides: error and warning objects; a polynomial base class;
- and some routines used in both the `polynomial` and `chebyshev` modules.
- Error objects
- -------------
- .. autosummary::
- :toctree: generated/
- PolyError base class for this sub-package's errors.
- PolyDomainError raised when domains are mismatched.
- Warning objects
- ---------------
- .. autosummary::
- :toctree: generated/
- RankWarning raised in least-squares fit for rank-deficient matrix.
- Base class
- ----------
- .. autosummary::
- :toctree: generated/
- PolyBase Obsolete base class for the polynomial classes. Do not use.
- Functions
- ---------
- .. autosummary::
- :toctree: generated/
- as_series convert list of array_likes into 1-D arrays of common type.
- trimseq remove trailing zeros.
- trimcoef remove small trailing coefficients.
- getdomain return the domain appropriate for a given set of abscissae.
- mapdomain maps points between domains.
- mapparms parameters of the linear map between domains.
- """
- from __future__ import division, absolute_import, print_function
- import numpy as np
- __all__ = [
- 'RankWarning', 'PolyError', 'PolyDomainError', 'as_series', 'trimseq',
- 'trimcoef', 'getdomain', 'mapdomain', 'mapparms', 'PolyBase']
- #
- # Warnings and Exceptions
- #
- class RankWarning(UserWarning):
- """Issued by chebfit when the design matrix is rank deficient."""
- pass
- class PolyError(Exception):
- """Base class for errors in this module."""
- pass
- class PolyDomainError(PolyError):
- """Issued by the generic Poly class when two domains don't match.
- This is raised when an binary operation is passed Poly objects with
- different domains.
- """
- pass
- #
- # Base class for all polynomial types
- #
- class PolyBase(object):
- """
- Base class for all polynomial types.
- Deprecated in numpy 1.9.0, use the abstract
- ABCPolyBase class instead. Note that the latter
- requires a number of virtual functions to be
- implemented.
- """
- pass
- #
- # Helper functions to convert inputs to 1-D arrays
- #
- def trimseq(seq):
- """Remove small Poly series coefficients.
- Parameters
- ----------
- seq : sequence
- Sequence of Poly series coefficients. This routine fails for
- empty sequences.
- Returns
- -------
- series : sequence
- Subsequence with trailing zeros removed. If the resulting sequence
- would be empty, return the first element. The returned sequence may
- or may not be a view.
- Notes
- -----
- Do not lose the type info if the sequence contains unknown objects.
- """
- if len(seq) == 0:
- return seq
- else:
- for i in range(len(seq) - 1, -1, -1):
- if seq[i] != 0:
- break
- return seq[:i+1]
- def as_series(alist, trim=True):
- """
- Return argument as a list of 1-d arrays.
- The returned list contains array(s) of dtype double, complex double, or
- object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of
- size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays
- of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array
- raises a Value Error if it is not first reshaped into either a 1-d or 2-d
- array.
- Parameters
- ----------
- alist : array_like
- A 1- or 2-d array_like
- trim : boolean, optional
- When True, trailing zeros are removed from the inputs.
- When False, the inputs are passed through intact.
- Returns
- -------
- [a1, a2,...] : list of 1-D arrays
- A copy of the input data as a list of 1-d arrays.
- Raises
- ------
- ValueError
- Raised when `as_series` cannot convert its input to 1-d arrays, or at
- least one of the resulting arrays is empty.
- Examples
- --------
- >>> from numpy.polynomial import polyutils as pu
- >>> a = np.arange(4)
- >>> pu.as_series(a)
- [array([ 0.]), array([ 1.]), array([ 2.]), array([ 3.])]
- >>> b = np.arange(6).reshape((2,3))
- >>> pu.as_series(b)
- [array([ 0., 1., 2.]), array([ 3., 4., 5.])]
- >>> pu.as_series((1, np.arange(3), np.arange(2, dtype=np.float16)))
- [array([ 1.]), array([ 0., 1., 2.]), array([ 0., 1.])]
- >>> pu.as_series([2, [1.1, 0.]])
- [array([ 2.]), array([ 1.1])]
- >>> pu.as_series([2, [1.1, 0.]], trim=False)
- [array([ 2.]), array([ 1.1, 0. ])]
- """
- arrays = [np.array(a, ndmin=1, copy=0) for a in alist]
- if min([a.size for a in arrays]) == 0:
- raise ValueError("Coefficient array is empty")
- if any([a.ndim != 1 for a in arrays]):
- raise ValueError("Coefficient array is not 1-d")
- if trim:
- arrays = [trimseq(a) for a in arrays]
- if any([a.dtype == np.dtype(object) for a in arrays]):
- ret = []
- for a in arrays:
- if a.dtype != np.dtype(object):
- tmp = np.empty(len(a), dtype=np.dtype(object))
- tmp[:] = a[:]
- ret.append(tmp)
- else:
- ret.append(a.copy())
- else:
- try:
- dtype = np.common_type(*arrays)
- except Exception:
- raise ValueError("Coefficient arrays have no common type")
- ret = [np.array(a, copy=1, dtype=dtype) for a in arrays]
- return ret
- def trimcoef(c, tol=0):
- """
- Remove "small" "trailing" coefficients from a polynomial.
- "Small" means "small in absolute value" and is controlled by the
- parameter `tol`; "trailing" means highest order coefficient(s), e.g., in
- ``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``)
- both the 3-rd and 4-th order coefficients would be "trimmed."
- Parameters
- ----------
- c : array_like
- 1-d array of coefficients, ordered from lowest order to highest.
- tol : number, optional
- Trailing (i.e., highest order) elements with absolute value less
- than or equal to `tol` (default value is zero) are removed.
- Returns
- -------
- trimmed : ndarray
- 1-d array with trailing zeros removed. If the resulting series
- would be empty, a series containing a single zero is returned.
- Raises
- ------
- ValueError
- If `tol` < 0
- See Also
- --------
- trimseq
- Examples
- --------
- >>> from numpy.polynomial import polyutils as pu
- >>> pu.trimcoef((0,0,3,0,5,0,0))
- array([ 0., 0., 3., 0., 5.])
- >>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
- array([ 0.])
- >>> i = complex(0,1) # works for complex
- >>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
- array([ 0.0003+0.j , 0.0010-0.001j])
- """
- if tol < 0:
- raise ValueError("tol must be non-negative")
- [c] = as_series([c])
- [ind] = np.nonzero(np.abs(c) > tol)
- if len(ind) == 0:
- return c[:1]*0
- else:
- return c[:ind[-1] + 1].copy()
- def getdomain(x):
- """
- Return a domain suitable for given abscissae.
- Find a domain suitable for a polynomial or Chebyshev series
- defined at the values supplied.
- Parameters
- ----------
- x : array_like
- 1-d array of abscissae whose domain will be determined.
- Returns
- -------
- domain : ndarray
- 1-d array containing two values. If the inputs are complex, then
- the two returned points are the lower left and upper right corners
- of the smallest rectangle (aligned with the axes) in the complex
- plane containing the points `x`. If the inputs are real, then the
- two points are the ends of the smallest interval containing the
- points `x`.
- See Also
- --------
- mapparms, mapdomain
- Examples
- --------
- >>> from numpy.polynomial import polyutils as pu
- >>> points = np.arange(4)**2 - 5; points
- array([-5, -4, -1, 4])
- >>> pu.getdomain(points)
- array([-5., 4.])
- >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
- >>> pu.getdomain(c)
- array([-1.-1.j, 1.+1.j])
- """
- [x] = as_series([x], trim=False)
- if x.dtype.char in np.typecodes['Complex']:
- rmin, rmax = x.real.min(), x.real.max()
- imin, imax = x.imag.min(), x.imag.max()
- return np.array((complex(rmin, imin), complex(rmax, imax)))
- else:
- return np.array((x.min(), x.max()))
- def mapparms(old, new):
- """
- Linear map parameters between domains.
- Return the parameters of the linear map ``offset + scale*x`` that maps
- `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
- Parameters
- ----------
- old, new : array_like
- Domains. Each domain must (successfully) convert to a 1-d array
- containing precisely two values.
- Returns
- -------
- offset, scale : scalars
- The map ``L(x) = offset + scale*x`` maps the first domain to the
- second.
- See Also
- --------
- getdomain, mapdomain
- Notes
- -----
- Also works for complex numbers, and thus can be used to calculate the
- parameters required to map any line in the complex plane to any other
- line therein.
- Examples
- --------
- >>> from numpy.polynomial import polyutils as pu
- >>> pu.mapparms((-1,1),(-1,1))
- (0.0, 1.0)
- >>> pu.mapparms((1,-1),(-1,1))
- (0.0, -1.0)
- >>> i = complex(0,1)
- >>> pu.mapparms((-i,-1),(1,i))
- ((1+1j), (1+0j))
- """
- oldlen = old[1] - old[0]
- newlen = new[1] - new[0]
- off = (old[1]*new[0] - old[0]*new[1])/oldlen
- scl = newlen/oldlen
- return off, scl
- def mapdomain(x, old, new):
- """
- Apply linear map to input points.
- The linear map ``offset + scale*x`` that maps the domain `old` to
- the domain `new` is applied to the points `x`.
- Parameters
- ----------
- x : array_like
- Points to be mapped. If `x` is a subtype of ndarray the subtype
- will be preserved.
- old, new : array_like
- The two domains that determine the map. Each must (successfully)
- convert to 1-d arrays containing precisely two values.
- Returns
- -------
- x_out : ndarray
- Array of points of the same shape as `x`, after application of the
- linear map between the two domains.
- See Also
- --------
- getdomain, mapparms
- Notes
- -----
- Effectively, this implements:
- .. math ::
- x\\_out = new[0] + m(x - old[0])
- where
- .. math ::
- m = \\frac{new[1]-new[0]}{old[1]-old[0]}
- Examples
- --------
- >>> from numpy.polynomial import polyutils as pu
- >>> old_domain = (-1,1)
- >>> new_domain = (0,2*np.pi)
- >>> x = np.linspace(-1,1,6); x
- array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
- >>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out
- array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825,
- 6.28318531])
- >>> x - pu.mapdomain(x_out, new_domain, old_domain)
- array([ 0., 0., 0., 0., 0., 0.])
- Also works for complex numbers (and thus can be used to map any line in
- the complex plane to any other line therein).
- >>> i = complex(0,1)
- >>> old = (-1 - i, 1 + i)
- >>> new = (-1 + i, 1 - i)
- >>> z = np.linspace(old[0], old[1], 6); z
- array([-1.0-1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1.0+1.j ])
- >>> new_z = P.mapdomain(z, old, new); new_z
- array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j, 0.2-0.2j, 0.6-0.6j, 1.0-1.j ])
- """
- x = np.asanyarray(x)
- off, scl = mapparms(old, new)
- return off + scl*x
|