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