Chris@87: """ Chris@87: Objects for dealing with polynomials. Chris@87: Chris@87: This module provides a number of objects (mostly functions) useful for Chris@87: dealing with polynomials, including a `Polynomial` class that Chris@87: encapsulates the usual arithmetic operations. (General information Chris@87: on how this module represents and works with polynomial objects is in Chris@87: the docstring for its "parent" sub-package, `numpy.polynomial`). Chris@87: Chris@87: Constants Chris@87: --------- Chris@87: - `polydomain` -- Polynomial default domain, [-1,1]. Chris@87: - `polyzero` -- (Coefficients of the) "zero polynomial." Chris@87: - `polyone` -- (Coefficients of the) constant polynomial 1. Chris@87: - `polyx` -- (Coefficients of the) identity map polynomial, ``f(x) = x``. Chris@87: Chris@87: Arithmetic Chris@87: ---------- Chris@87: - `polyadd` -- add two polynomials. Chris@87: - `polysub` -- subtract one polynomial from another. Chris@87: - `polymul` -- multiply two polynomials. Chris@87: - `polydiv` -- divide one polynomial by another. Chris@87: - `polypow` -- raise a polynomial to an positive integer power Chris@87: - `polyval` -- evaluate a polynomial at given points. Chris@87: - `polyval2d` -- evaluate a 2D polynomial at given points. Chris@87: - `polyval3d` -- evaluate a 3D polynomial at given points. Chris@87: - `polygrid2d` -- evaluate a 2D polynomial on a Cartesian product. Chris@87: - `polygrid3d` -- evaluate a 3D polynomial on a Cartesian product. Chris@87: Chris@87: Calculus Chris@87: -------- Chris@87: - `polyder` -- differentiate a polynomial. Chris@87: - `polyint` -- integrate a polynomial. Chris@87: Chris@87: Misc Functions Chris@87: -------------- Chris@87: - `polyfromroots` -- create a polynomial with specified roots. Chris@87: - `polyroots` -- find the roots of a polynomial. Chris@87: - `polyvander` -- Vandermonde-like matrix for powers. Chris@87: - `polyvander2d` -- Vandermonde-like matrix for 2D power series. Chris@87: - `polyvander3d` -- Vandermonde-like matrix for 3D power series. Chris@87: - `polycompanion` -- companion matrix in power series form. Chris@87: - `polyfit` -- least-squares fit returning a polynomial. Chris@87: - `polytrim` -- trim leading coefficients from a polynomial. Chris@87: - `polyline` -- polynomial representing given straight line. Chris@87: Chris@87: Classes Chris@87: ------- Chris@87: - `Polynomial` -- polynomial class. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: `numpy.polynomial` Chris@87: Chris@87: """ Chris@87: from __future__ import division, absolute_import, print_function Chris@87: Chris@87: __all__ = [ Chris@87: 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd', Chris@87: 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval', Chris@87: 'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit', Chris@87: 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d', Chris@87: 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d'] Chris@87: Chris@87: import warnings Chris@87: import numpy as np Chris@87: import numpy.linalg as la Chris@87: Chris@87: from . import polyutils as pu Chris@87: from ._polybase import ABCPolyBase Chris@87: Chris@87: polytrim = pu.trimcoef Chris@87: Chris@87: # Chris@87: # These are constant arrays are of integer type so as to be compatible Chris@87: # with the widest range of other types, such as Decimal. Chris@87: # Chris@87: Chris@87: # Polynomial default domain. Chris@87: polydomain = np.array([-1, 1]) Chris@87: Chris@87: # Polynomial coefficients representing zero. Chris@87: polyzero = np.array([0]) Chris@87: Chris@87: # Polynomial coefficients representing one. Chris@87: polyone = np.array([1]) Chris@87: Chris@87: # Polynomial coefficients representing the identity x. Chris@87: polyx = np.array([0, 1]) Chris@87: Chris@87: # Chris@87: # Polynomial series functions Chris@87: # Chris@87: Chris@87: Chris@87: def polyline(off, scl): Chris@87: """ Chris@87: Returns an array representing a linear polynomial. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: off, scl : scalars Chris@87: The "y-intercept" and "slope" of the line, respectively. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: y : ndarray Chris@87: This module's representation of the linear polynomial ``off + Chris@87: scl*x``. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: chebline Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> P.polyline(1,-1) Chris@87: array([ 1, -1]) Chris@87: >>> P.polyval(1, P.polyline(1,-1)) # should be 0 Chris@87: 0.0 Chris@87: Chris@87: """ Chris@87: if scl != 0: Chris@87: return np.array([off, scl]) Chris@87: else: Chris@87: return np.array([off]) Chris@87: Chris@87: Chris@87: def polyfromroots(roots): Chris@87: """ Chris@87: Generate a monic polynomial with given roots. Chris@87: Chris@87: Return the coefficients of the polynomial Chris@87: Chris@87: .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), Chris@87: Chris@87: where the `r_n` are the roots specified in `roots`. If a zero has Chris@87: multiplicity n, then it must appear in `roots` n times. For instance, Chris@87: if 2 is a root of multiplicity three and 3 is a root of multiplicity 2, Chris@87: then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear Chris@87: in any order. Chris@87: Chris@87: If the returned coefficients are `c`, then Chris@87: Chris@87: .. math:: p(x) = c_0 + c_1 * x + ... + x^n Chris@87: Chris@87: The coefficient of the last term is 1 for monic polynomials in this Chris@87: form. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: roots : array_like Chris@87: Sequence containing the roots. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: 1-D array of the polynomial's coefficients If all the roots are Chris@87: real, then `out` is also real, otherwise it is complex. (see Chris@87: Examples below). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: chebfromroots, legfromroots, lagfromroots, hermfromroots Chris@87: hermefromroots Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The coefficients are determined by multiplying together linear factors Chris@87: of the form `(x - r_i)`, i.e. Chris@87: Chris@87: .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n) Chris@87: Chris@87: where ``n == len(roots) - 1``; note that this implies that `1` is always Chris@87: returned for :math:`a_n`. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> P.polyfromroots((-1,0,1)) # x(x - 1)(x + 1) = x^3 - x Chris@87: array([ 0., -1., 0., 1.]) Chris@87: >>> j = complex(0,1) Chris@87: >>> P.polyfromroots((-j,j)) # complex returned, though values are real Chris@87: array([ 1.+0.j, 0.+0.j, 1.+0.j]) Chris@87: Chris@87: """ Chris@87: if len(roots) == 0: Chris@87: return np.ones(1) Chris@87: else: Chris@87: [roots] = pu.as_series([roots], trim=False) Chris@87: roots.sort() Chris@87: p = [polyline(-r, 1) for r in roots] Chris@87: n = len(p) Chris@87: while n > 1: Chris@87: m, r = divmod(n, 2) Chris@87: tmp = [polymul(p[i], p[i+m]) for i in range(m)] Chris@87: if r: Chris@87: tmp[0] = polymul(tmp[0], p[-1]) Chris@87: p = tmp Chris@87: n = m Chris@87: return p[0] Chris@87: Chris@87: Chris@87: def polyadd(c1, c2): Chris@87: """ Chris@87: Add one polynomial to another. Chris@87: Chris@87: Returns the sum of two polynomials `c1` + `c2`. The arguments are Chris@87: sequences of coefficients from lowest order term to highest, i.e., Chris@87: [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c1, c2 : array_like Chris@87: 1-D arrays of polynomial coefficients ordered from low to high. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The coefficient array representing their sum. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polysub, polymul, polydiv, polypow Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c1 = (1,2,3) Chris@87: >>> c2 = (3,2,1) Chris@87: >>> sum = P.polyadd(c1,c2); sum Chris@87: array([ 4., 4., 4.]) Chris@87: >>> P.polyval(2, sum) # 4 + 4(2) + 4(2**2) Chris@87: 28.0 Chris@87: Chris@87: """ Chris@87: # c1, c2 are trimmed copies Chris@87: [c1, c2] = pu.as_series([c1, c2]) Chris@87: if len(c1) > len(c2): Chris@87: c1[:c2.size] += c2 Chris@87: ret = c1 Chris@87: else: Chris@87: c2[:c1.size] += c1 Chris@87: ret = c2 Chris@87: return pu.trimseq(ret) Chris@87: Chris@87: Chris@87: def polysub(c1, c2): Chris@87: """ Chris@87: Subtract one polynomial from another. Chris@87: Chris@87: Returns the difference of two polynomials `c1` - `c2`. The arguments Chris@87: are sequences of coefficients from lowest order term to highest, i.e., Chris@87: [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c1, c2 : array_like Chris@87: 1-D arrays of polynomial coefficients ordered from low to Chris@87: high. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Of coefficients representing their difference. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyadd, polymul, polydiv, polypow Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c1 = (1,2,3) Chris@87: >>> c2 = (3,2,1) Chris@87: >>> P.polysub(c1,c2) Chris@87: array([-2., 0., 2.]) Chris@87: >>> P.polysub(c2,c1) # -P.polysub(c1,c2) Chris@87: array([ 2., 0., -2.]) Chris@87: Chris@87: """ Chris@87: # c1, c2 are trimmed copies Chris@87: [c1, c2] = pu.as_series([c1, c2]) Chris@87: if len(c1) > len(c2): Chris@87: c1[:c2.size] -= c2 Chris@87: ret = c1 Chris@87: else: Chris@87: c2 = -c2 Chris@87: c2[:c1.size] += c1 Chris@87: ret = c2 Chris@87: return pu.trimseq(ret) Chris@87: Chris@87: Chris@87: def polymulx(c): Chris@87: """Multiply a polynomial by x. Chris@87: Chris@87: Multiply the polynomial `c` by x, where x is the independent Chris@87: variable. Chris@87: Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : array_like Chris@87: 1-D array of polynomial coefficients ordered from low to Chris@87: high. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Array representing the result of the multiplication. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.5.0 Chris@87: Chris@87: """ Chris@87: # c is a trimmed copy Chris@87: [c] = pu.as_series([c]) Chris@87: # The zero series needs special treatment Chris@87: if len(c) == 1 and c[0] == 0: Chris@87: return c Chris@87: Chris@87: prd = np.empty(len(c) + 1, dtype=c.dtype) Chris@87: prd[0] = c[0]*0 Chris@87: prd[1:] = c Chris@87: return prd Chris@87: Chris@87: Chris@87: def polymul(c1, c2): Chris@87: """ Chris@87: Multiply one polynomial by another. Chris@87: Chris@87: Returns the product of two polynomials `c1` * `c2`. The arguments are Chris@87: sequences of coefficients, from lowest order term to highest, e.g., Chris@87: [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.`` Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c1, c2 : array_like Chris@87: 1-D arrays of coefficients representing a polynomial, relative to the Chris@87: "standard" basis, and ordered from lowest order term to highest. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Of the coefficients of their product. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyadd, polysub, polydiv, polypow Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c1 = (1,2,3) Chris@87: >>> c2 = (3,2,1) Chris@87: >>> P.polymul(c1,c2) Chris@87: array([ 3., 8., 14., 8., 3.]) Chris@87: Chris@87: """ Chris@87: # c1, c2 are trimmed copies Chris@87: [c1, c2] = pu.as_series([c1, c2]) Chris@87: ret = np.convolve(c1, c2) Chris@87: return pu.trimseq(ret) Chris@87: Chris@87: Chris@87: def polydiv(c1, c2): Chris@87: """ Chris@87: Divide one polynomial by another. Chris@87: Chris@87: Returns the quotient-with-remainder of two polynomials `c1` / `c2`. Chris@87: The arguments are sequences of coefficients, from lowest order term Chris@87: to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c1, c2 : array_like Chris@87: 1-D arrays of polynomial coefficients ordered from low to high. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: [quo, rem] : ndarrays Chris@87: Of coefficient series representing the quotient and remainder. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyadd, polysub, polymul, polypow Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c1 = (1,2,3) Chris@87: >>> c2 = (3,2,1) Chris@87: >>> P.polydiv(c1,c2) Chris@87: (array([ 3.]), array([-8., -4.])) Chris@87: >>> P.polydiv(c2,c1) Chris@87: (array([ 0.33333333]), array([ 2.66666667, 1.33333333])) Chris@87: Chris@87: """ Chris@87: # c1, c2 are trimmed copies Chris@87: [c1, c2] = pu.as_series([c1, c2]) Chris@87: if c2[-1] == 0: Chris@87: raise ZeroDivisionError() Chris@87: Chris@87: len1 = len(c1) Chris@87: len2 = len(c2) Chris@87: if len2 == 1: Chris@87: return c1/c2[-1], c1[:1]*0 Chris@87: elif len1 < len2: Chris@87: return c1[:1]*0, c1 Chris@87: else: Chris@87: dlen = len1 - len2 Chris@87: scl = c2[-1] Chris@87: c2 = c2[:-1]/scl Chris@87: i = dlen Chris@87: j = len1 - 1 Chris@87: while i >= 0: Chris@87: c1[i:j] -= c2*c1[j] Chris@87: i -= 1 Chris@87: j -= 1 Chris@87: return c1[j+1:]/scl, pu.trimseq(c1[:j+1]) Chris@87: Chris@87: Chris@87: def polypow(c, pow, maxpower=None): Chris@87: """Raise a polynomial to a power. Chris@87: Chris@87: Returns the polynomial `c` raised to the power `pow`. The argument Chris@87: `c` is a sequence of coefficients ordered from low to high. i.e., Chris@87: [1,2,3] is the series ``1 + 2*x + 3*x**2.`` Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : array_like Chris@87: 1-D array of array of series coefficients ordered from low to Chris@87: high degree. Chris@87: pow : integer Chris@87: Power to which the series will be raised Chris@87: maxpower : integer, optional Chris@87: Maximum power allowed. This is mainly to limit growth of the series Chris@87: to unmanageable size. Default is 16 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: coef : ndarray Chris@87: Power series of power. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyadd, polysub, polymul, polydiv Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Chris@87: """ Chris@87: # c is a trimmed copy Chris@87: [c] = pu.as_series([c]) Chris@87: power = int(pow) Chris@87: if power != pow or power < 0: Chris@87: raise ValueError("Power must be a non-negative integer.") Chris@87: elif maxpower is not None and power > maxpower: Chris@87: raise ValueError("Power is too large") Chris@87: elif power == 0: Chris@87: return np.array([1], dtype=c.dtype) Chris@87: elif power == 1: Chris@87: return c Chris@87: else: Chris@87: # This can be made more efficient by using powers of two Chris@87: # in the usual way. Chris@87: prd = c Chris@87: for i in range(2, power + 1): Chris@87: prd = np.convolve(prd, c) Chris@87: return prd Chris@87: Chris@87: Chris@87: def polyder(c, m=1, scl=1, axis=0): Chris@87: """ Chris@87: Differentiate a polynomial. Chris@87: Chris@87: Returns the polynomial coefficients `c` differentiated `m` times along Chris@87: `axis`. At each iteration the result is multiplied by `scl` (the Chris@87: scaling factor is for use in a linear change of variable). The Chris@87: argument `c` is an array of coefficients from low to high degree along Chris@87: each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2`` Chris@87: while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is Chris@87: ``x`` and axis=1 is ``y``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : array_like Chris@87: Array of polynomial coefficients. If c is multidimensional the Chris@87: different axis correspond to different variables with the degree Chris@87: in each axis given by the corresponding index. Chris@87: m : int, optional Chris@87: Number of derivatives taken, must be non-negative. (Default: 1) Chris@87: scl : scalar, optional Chris@87: Each differentiation is multiplied by `scl`. The end result is Chris@87: multiplication by ``scl**m``. This is for use in a linear change Chris@87: of variable. (Default: 1) Chris@87: axis : int, optional Chris@87: Axis over which the derivative is taken. (Default: 0). Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: der : ndarray Chris@87: Polynomial coefficients of the derivative. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyint Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3 Chris@87: >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2 Chris@87: array([ 2., 6., 12.]) Chris@87: >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24 Chris@87: array([ 24.]) Chris@87: >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2 Chris@87: array([ -2., -6., -12.]) Chris@87: >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x Chris@87: array([ 6., 24.]) Chris@87: Chris@87: """ Chris@87: c = np.array(c, ndmin=1, copy=1) Chris@87: if c.dtype.char in '?bBhHiIlLqQpP': Chris@87: # astype fails with NA Chris@87: c = c + 0.0 Chris@87: cdt = c.dtype Chris@87: cnt, iaxis = [int(t) for t in [m, axis]] Chris@87: Chris@87: if cnt != m: Chris@87: raise ValueError("The order of derivation must be integer") Chris@87: if cnt < 0: Chris@87: raise ValueError("The order of derivation must be non-negative") Chris@87: if iaxis != axis: Chris@87: raise ValueError("The axis must be integer") Chris@87: if not -c.ndim <= iaxis < c.ndim: Chris@87: raise ValueError("The axis is out of range") Chris@87: if iaxis < 0: Chris@87: iaxis += c.ndim Chris@87: Chris@87: if cnt == 0: Chris@87: return c Chris@87: Chris@87: c = np.rollaxis(c, iaxis) Chris@87: n = len(c) Chris@87: if cnt >= n: Chris@87: c = c[:1]*0 Chris@87: else: Chris@87: for i in range(cnt): Chris@87: n = n - 1 Chris@87: c *= scl Chris@87: der = np.empty((n,) + c.shape[1:], dtype=cdt) Chris@87: for j in range(n, 0, -1): Chris@87: der[j - 1] = j*c[j] Chris@87: c = der Chris@87: c = np.rollaxis(c, 0, iaxis + 1) Chris@87: return c Chris@87: Chris@87: Chris@87: def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Chris@87: """ Chris@87: Integrate a polynomial. Chris@87: Chris@87: Returns the polynomial coefficients `c` integrated `m` times from Chris@87: `lbnd` along `axis`. At each iteration the resulting series is Chris@87: **multiplied** by `scl` and an integration constant, `k`, is added. Chris@87: The scaling factor is for use in a linear change of variable. ("Buyer Chris@87: beware": note that, depending on what one is doing, one may want `scl` Chris@87: to be the reciprocal of what one might expect; for more information, Chris@87: see the Notes section below.) The argument `c` is an array of Chris@87: coefficients, from low to high degree along each axis, e.g., [1,2,3] Chris@87: represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]] Chris@87: represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is Chris@87: ``y``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : array_like Chris@87: 1-D array of polynomial coefficients, ordered from low to high. Chris@87: m : int, optional Chris@87: Order of integration, must be positive. (Default: 1) Chris@87: k : {[], list, scalar}, optional Chris@87: Integration constant(s). The value of the first integral at zero Chris@87: is the first value in the list, the value of the second integral Chris@87: at zero is the second value, etc. If ``k == []`` (the default), Chris@87: all constants are set to zero. If ``m == 1``, a single scalar can Chris@87: be given instead of a list. Chris@87: lbnd : scalar, optional Chris@87: The lower bound of the integral. (Default: 0) Chris@87: scl : scalar, optional Chris@87: Following each integration the result is *multiplied* by `scl` Chris@87: before the integration constant is added. (Default: 1) Chris@87: axis : int, optional Chris@87: Axis over which the integral is taken. (Default: 0). Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: S : ndarray Chris@87: Coefficient array of the integral. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: ValueError Chris@87: If ``m < 1``, ``len(k) > m``. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyder Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Note that the result of each integration is *multiplied* by `scl`. Why Chris@87: is this important to note? Say one is making a linear change of Chris@87: variable :math:`u = ax + b` in an integral relative to `x`. Then Chris@87: .. math::`dx = du/a`, so one will need to set `scl` equal to Chris@87: :math:`1/a` - perhaps not what one would have first thought. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> c = (1,2,3) Chris@87: >>> P.polyint(c) # should return array([0, 1, 1, 1]) Chris@87: array([ 0., 1., 1., 1.]) Chris@87: >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20]) Chris@87: array([ 0. , 0. , 0. , 0.16666667, 0.08333333, Chris@87: 0.05 ]) Chris@87: >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1]) Chris@87: array([ 3., 1., 1., 1.]) Chris@87: >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1]) Chris@87: array([ 6., 1., 1., 1.]) Chris@87: >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2]) Chris@87: array([ 0., -2., -2., -2.]) Chris@87: Chris@87: """ Chris@87: c = np.array(c, ndmin=1, copy=1) Chris@87: if c.dtype.char in '?bBhHiIlLqQpP': Chris@87: # astype doesn't preserve mask attribute. Chris@87: c = c + 0.0 Chris@87: cdt = c.dtype Chris@87: if not np.iterable(k): Chris@87: k = [k] Chris@87: cnt, iaxis = [int(t) for t in [m, axis]] Chris@87: Chris@87: if cnt != m: Chris@87: raise ValueError("The order of integration must be integer") Chris@87: if cnt < 0: Chris@87: raise ValueError("The order of integration must be non-negative") Chris@87: if len(k) > cnt: Chris@87: raise ValueError("Too many integration constants") Chris@87: if iaxis != axis: Chris@87: raise ValueError("The axis must be integer") Chris@87: if not -c.ndim <= iaxis < c.ndim: Chris@87: raise ValueError("The axis is out of range") Chris@87: if iaxis < 0: Chris@87: iaxis += c.ndim Chris@87: Chris@87: if cnt == 0: Chris@87: return c Chris@87: Chris@87: k = list(k) + [0]*(cnt - len(k)) Chris@87: c = np.rollaxis(c, iaxis) Chris@87: for i in range(cnt): Chris@87: n = len(c) Chris@87: c *= scl Chris@87: if n == 1 and np.all(c[0] == 0): Chris@87: c[0] += k[i] Chris@87: else: Chris@87: tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt) Chris@87: tmp[0] = c[0]*0 Chris@87: tmp[1] = c[0] Chris@87: for j in range(1, n): Chris@87: tmp[j + 1] = c[j]/(j + 1) Chris@87: tmp[0] += k[i] - polyval(lbnd, tmp) Chris@87: c = tmp Chris@87: c = np.rollaxis(c, 0, iaxis + 1) Chris@87: return c Chris@87: Chris@87: Chris@87: def polyval(x, c, tensor=True): Chris@87: """ Chris@87: Evaluate a polynomial at points x. Chris@87: Chris@87: If `c` is of length `n + 1`, this function returns the value Chris@87: Chris@87: .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n Chris@87: Chris@87: The parameter `x` is converted to an array only if it is a tuple or a Chris@87: list, otherwise it is treated as a scalar. In either case, either `x` Chris@87: or its elements must support multiplication and addition both with Chris@87: themselves and with the elements of `c`. Chris@87: Chris@87: If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If Chris@87: `c` is multidimensional, then the shape of the result depends on the Chris@87: value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + Chris@87: x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that Chris@87: scalars have shape (,). Chris@87: Chris@87: Trailing zeros in the coefficients will be used in the evaluation, so Chris@87: they should be avoided if efficiency is a concern. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like, compatible object Chris@87: If `x` is a list or tuple, it is converted to an ndarray, otherwise Chris@87: it is left unchanged and treated as a scalar. In either case, `x` Chris@87: or its elements must support addition and multiplication with Chris@87: with themselves and with the elements of `c`. Chris@87: c : array_like Chris@87: Array of coefficients ordered so that the coefficients for terms of Chris@87: degree n are contained in c[n]. If `c` is multidimensional the Chris@87: remaining indices enumerate multiple polynomials. In the two Chris@87: dimensional case the coefficients may be thought of as stored in Chris@87: the columns of `c`. Chris@87: tensor : boolean, optional Chris@87: If True, the shape of the coefficient array is extended with ones Chris@87: on the right, one for each dimension of `x`. Scalars have dimension 0 Chris@87: for this action. The result is that every column of coefficients in Chris@87: `c` is evaluated for every element of `x`. If False, `x` is broadcast Chris@87: over the columns of `c` for the evaluation. This keyword is useful Chris@87: when `c` is multidimensional. The default value is True. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: values : ndarray, compatible object Chris@87: The shape of the returned array is described above. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyval2d, polygrid2d, polyval3d, polygrid3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The evaluation uses Horner's method. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial.polynomial import polyval Chris@87: >>> polyval(1, [1,2,3]) Chris@87: 6.0 Chris@87: >>> a = np.arange(4).reshape(2,2) Chris@87: >>> a Chris@87: array([[0, 1], Chris@87: [2, 3]]) Chris@87: >>> polyval(a, [1,2,3]) Chris@87: array([[ 1., 6.], Chris@87: [ 17., 34.]]) Chris@87: >>> coef = np.arange(4).reshape(2,2) # multidimensional coefficients Chris@87: >>> coef Chris@87: array([[0, 1], Chris@87: [2, 3]]) Chris@87: >>> polyval([1,2], coef, tensor=True) Chris@87: array([[ 2., 4.], Chris@87: [ 4., 7.]]) Chris@87: >>> polyval([1,2], coef, tensor=False) Chris@87: array([ 2., 7.]) Chris@87: Chris@87: """ Chris@87: c = np.array(c, ndmin=1, copy=0) Chris@87: if c.dtype.char in '?bBhHiIlLqQpP': Chris@87: # astype fails with NA Chris@87: c = c + 0.0 Chris@87: if isinstance(x, (tuple, list)): Chris@87: x = np.asarray(x) Chris@87: if isinstance(x, np.ndarray) and tensor: Chris@87: c = c.reshape(c.shape + (1,)*x.ndim) Chris@87: Chris@87: c0 = c[-1] + x*0 Chris@87: for i in range(2, len(c) + 1): Chris@87: c0 = c[-i] + c0*x Chris@87: return c0 Chris@87: Chris@87: Chris@87: def polyval2d(x, y, c): Chris@87: """ Chris@87: Evaluate a 2-D polynomial at points (x, y). Chris@87: Chris@87: This function returns the value Chris@87: Chris@87: .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j Chris@87: Chris@87: The parameters `x` and `y` are converted to arrays only if they are Chris@87: tuples or a lists, otherwise they are treated as a scalars and they Chris@87: must have the same shape after conversion. In either case, either `x` Chris@87: and `y` or their elements must support multiplication and addition both Chris@87: with themselves and with the elements of `c`. Chris@87: Chris@87: If `c` has fewer than two dimensions, ones are implicitly appended to Chris@87: its shape to make it 2-D. The shape of the result will be c.shape[2:] + Chris@87: x.shape. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y : array_like, compatible objects Chris@87: The two dimensional series is evaluated at the points `(x, y)`, Chris@87: where `x` and `y` must have the same shape. If `x` or `y` is a list Chris@87: or tuple, it is first converted to an ndarray, otherwise it is left Chris@87: unchanged and, if it isn't an ndarray, it is treated as a scalar. Chris@87: c : array_like Chris@87: Array of coefficients ordered so that the coefficient of the term Chris@87: of multi-degree i,j is contained in `c[i,j]`. If `c` has Chris@87: dimension greater than two the remaining indices enumerate multiple Chris@87: sets of coefficients. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: values : ndarray, compatible object Chris@87: The values of the two dimensional polynomial at points formed with Chris@87: pairs of corresponding values from `x` and `y`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyval, polygrid2d, polyval3d, polygrid3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: try: Chris@87: x, y = np.array((x, y), copy=0) Chris@87: except: Chris@87: raise ValueError('x, y are incompatible') Chris@87: Chris@87: c = polyval(x, c) Chris@87: c = polyval(y, c, tensor=False) Chris@87: return c Chris@87: Chris@87: Chris@87: def polygrid2d(x, y, c): Chris@87: """ Chris@87: Evaluate a 2-D polynomial on the Cartesian product of x and y. Chris@87: Chris@87: This function returns the values: Chris@87: Chris@87: .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j Chris@87: Chris@87: where the points `(a, b)` consist of all pairs formed by taking Chris@87: `a` from `x` and `b` from `y`. The resulting points form a grid with Chris@87: `x` in the first dimension and `y` in the second. Chris@87: Chris@87: The parameters `x` and `y` are converted to arrays only if they are Chris@87: tuples or a lists, otherwise they are treated as a scalars. In either Chris@87: case, either `x` and `y` or their elements must support multiplication Chris@87: and addition both with themselves and with the elements of `c`. Chris@87: Chris@87: If `c` has fewer than two dimensions, ones are implicitly appended to Chris@87: its shape to make it 2-D. The shape of the result will be c.shape[2:] + Chris@87: x.shape + y.shape. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y : array_like, compatible objects Chris@87: The two dimensional series is evaluated at the points in the Chris@87: Cartesian product of `x` and `y`. If `x` or `y` is a list or Chris@87: tuple, it is first converted to an ndarray, otherwise it is left Chris@87: unchanged and, if it isn't an ndarray, it is treated as a scalar. Chris@87: c : array_like Chris@87: Array of coefficients ordered so that the coefficients for terms of Chris@87: degree i,j are contained in ``c[i,j]``. If `c` has dimension Chris@87: greater than two the remaining indices enumerate multiple sets of Chris@87: coefficients. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: values : ndarray, compatible object Chris@87: The values of the two dimensional polynomial at points in the Cartesian Chris@87: product of `x` and `y`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyval, polyval2d, polyval3d, polygrid3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: c = polyval(x, c) Chris@87: c = polyval(y, c) Chris@87: return c Chris@87: Chris@87: Chris@87: def polyval3d(x, y, z, c): Chris@87: """ Chris@87: Evaluate a 3-D polynomial at points (x, y, z). Chris@87: Chris@87: This function returns the values: Chris@87: Chris@87: .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k Chris@87: Chris@87: The parameters `x`, `y`, and `z` are converted to arrays only if Chris@87: they are tuples or a lists, otherwise they are treated as a scalars and Chris@87: they must have the same shape after conversion. In either case, either Chris@87: `x`, `y`, and `z` or their elements must support multiplication and Chris@87: addition both with themselves and with the elements of `c`. Chris@87: Chris@87: If `c` has fewer than 3 dimensions, ones are implicitly appended to its Chris@87: shape to make it 3-D. The shape of the result will be c.shape[3:] + Chris@87: x.shape. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y, z : array_like, compatible object Chris@87: The three dimensional series is evaluated at the points Chris@87: `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If Chris@87: any of `x`, `y`, or `z` is a list or tuple, it is first converted Chris@87: to an ndarray, otherwise it is left unchanged and if it isn't an Chris@87: ndarray it is treated as a scalar. Chris@87: c : array_like Chris@87: Array of coefficients ordered so that the coefficient of the term of Chris@87: multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension Chris@87: greater than 3 the remaining indices enumerate multiple sets of Chris@87: coefficients. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: values : ndarray, compatible object Chris@87: The values of the multidimensional polynomial on points formed with Chris@87: triples of corresponding values from `x`, `y`, and `z`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyval, polyval2d, polygrid2d, polygrid3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: try: Chris@87: x, y, z = np.array((x, y, z), copy=0) Chris@87: except: Chris@87: raise ValueError('x, y, z are incompatible') Chris@87: Chris@87: c = polyval(x, c) Chris@87: c = polyval(y, c, tensor=False) Chris@87: c = polyval(z, c, tensor=False) Chris@87: return c Chris@87: Chris@87: Chris@87: def polygrid3d(x, y, z, c): Chris@87: """ Chris@87: Evaluate a 3-D polynomial on the Cartesian product of x, y and z. Chris@87: Chris@87: This function returns the values: Chris@87: Chris@87: .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k Chris@87: Chris@87: where the points `(a, b, c)` consist of all triples formed by taking Chris@87: `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form Chris@87: a grid with `x` in the first dimension, `y` in the second, and `z` in Chris@87: the third. Chris@87: Chris@87: The parameters `x`, `y`, and `z` are converted to arrays only if they Chris@87: are tuples or a lists, otherwise they are treated as a scalars. In Chris@87: either case, either `x`, `y`, and `z` or their elements must support Chris@87: multiplication and addition both with themselves and with the elements Chris@87: of `c`. Chris@87: Chris@87: If `c` has fewer than three dimensions, ones are implicitly appended to Chris@87: its shape to make it 3-D. The shape of the result will be c.shape[3:] + Chris@87: x.shape + y.shape + z.shape. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y, z : array_like, compatible objects Chris@87: The three dimensional series is evaluated at the points in the Chris@87: Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a Chris@87: list or tuple, it is first converted to an ndarray, otherwise it is Chris@87: left unchanged and, if it isn't an ndarray, it is treated as a Chris@87: scalar. Chris@87: c : array_like Chris@87: Array of coefficients ordered so that the coefficients for terms of Chris@87: degree i,j are contained in ``c[i,j]``. If `c` has dimension Chris@87: greater than two the remaining indices enumerate multiple sets of Chris@87: coefficients. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: values : ndarray, compatible object Chris@87: The values of the two dimensional polynomial at points in the Cartesian Chris@87: product of `x` and `y`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyval, polyval2d, polygrid2d, polyval3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: c = polyval(x, c) Chris@87: c = polyval(y, c) Chris@87: c = polyval(z, c) Chris@87: return c Chris@87: Chris@87: Chris@87: def polyvander(x, deg): Chris@87: """Vandermonde matrix of given degree. Chris@87: Chris@87: Returns the Vandermonde matrix of degree `deg` and sample points Chris@87: `x`. The Vandermonde matrix is defined by Chris@87: Chris@87: .. math:: V[..., i] = x^i, Chris@87: Chris@87: where `0 <= i <= deg`. The leading indices of `V` index the elements of Chris@87: `x` and the last index is the power of `x`. Chris@87: Chris@87: If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the Chris@87: matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and Chris@87: ``polyval(x, c)`` are the same up to roundoff. This equivalence is Chris@87: useful both for least squares fitting and for the evaluation of a large Chris@87: number of polynomials of the same degree and sample points. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like Chris@87: Array of points. The dtype is converted to float64 or complex128 Chris@87: depending on whether any of the elements are complex. If `x` is Chris@87: scalar it is converted to a 1-D array. Chris@87: deg : int Chris@87: Degree of the resulting matrix. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: vander : ndarray. Chris@87: The Vandermonde matrix. The shape of the returned matrix is Chris@87: ``x.shape + (deg + 1,)``, where the last index is the power of `x`. Chris@87: The dtype will be the same as the converted `x`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyvander2d, polyvander3d Chris@87: Chris@87: """ Chris@87: ideg = int(deg) Chris@87: if ideg != deg: Chris@87: raise ValueError("deg must be integer") Chris@87: if ideg < 0: Chris@87: raise ValueError("deg must be non-negative") Chris@87: Chris@87: x = np.array(x, copy=0, ndmin=1) + 0.0 Chris@87: dims = (ideg + 1,) + x.shape Chris@87: dtyp = x.dtype Chris@87: v = np.empty(dims, dtype=dtyp) Chris@87: v[0] = x*0 + 1 Chris@87: if ideg > 0: Chris@87: v[1] = x Chris@87: for i in range(2, ideg + 1): Chris@87: v[i] = v[i-1]*x Chris@87: return np.rollaxis(v, 0, v.ndim) Chris@87: Chris@87: Chris@87: def polyvander2d(x, y, deg): Chris@87: """Pseudo-Vandermonde matrix of given degrees. Chris@87: Chris@87: Returns the pseudo-Vandermonde matrix of degrees `deg` and sample Chris@87: points `(x, y)`. The pseudo-Vandermonde matrix is defined by Chris@87: Chris@87: .. math:: V[..., deg[1]*i + j] = x^i * y^j, Chris@87: Chris@87: where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of Chris@87: `V` index the points `(x, y)` and the last index encodes the powers of Chris@87: `x` and `y`. Chris@87: Chris@87: If ``V = polyvander2d(x, y, [xdeg, ydeg])``, then the columns of `V` Chris@87: correspond to the elements of a 2-D coefficient array `c` of shape Chris@87: (xdeg + 1, ydeg + 1) in the order Chris@87: Chris@87: .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ... Chris@87: Chris@87: and ``np.dot(V, c.flat)`` and ``polyval2d(x, y, c)`` will be the same Chris@87: up to roundoff. This equivalence is useful both for least squares Chris@87: fitting and for the evaluation of a large number of 2-D polynomials Chris@87: of the same degrees and sample points. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y : array_like Chris@87: Arrays of point coordinates, all of the same shape. The dtypes Chris@87: will be converted to either float64 or complex128 depending on Chris@87: whether any of the elements are complex. Scalars are converted to Chris@87: 1-D arrays. Chris@87: deg : list of ints Chris@87: List of maximum degrees of the form [x_deg, y_deg]. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: vander2d : ndarray Chris@87: The shape of the returned matrix is ``x.shape + (order,)``, where Chris@87: :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same Chris@87: as the converted `x` and `y`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyvander, polyvander3d. polyval2d, polyval3d Chris@87: Chris@87: """ Chris@87: ideg = [int(d) for d in deg] Chris@87: is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] Chris@87: if is_valid != [1, 1]: Chris@87: raise ValueError("degrees must be non-negative integers") Chris@87: degx, degy = ideg Chris@87: x, y = np.array((x, y), copy=0) + 0.0 Chris@87: Chris@87: vx = polyvander(x, degx) Chris@87: vy = polyvander(y, degy) Chris@87: v = vx[..., None]*vy[..., None,:] Chris@87: # einsum bug Chris@87: #v = np.einsum("...i,...j->...ij", vx, vy) Chris@87: return v.reshape(v.shape[:-2] + (-1,)) Chris@87: Chris@87: Chris@87: def polyvander3d(x, y, z, deg): Chris@87: """Pseudo-Vandermonde matrix of given degrees. Chris@87: Chris@87: Returns the pseudo-Vandermonde matrix of degrees `deg` and sample Chris@87: points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`, Chris@87: then The pseudo-Vandermonde matrix is defined by Chris@87: Chris@87: .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k, Chris@87: Chris@87: where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading Chris@87: indices of `V` index the points `(x, y, z)` and the last index encodes Chris@87: the powers of `x`, `y`, and `z`. Chris@87: Chris@87: If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns Chris@87: of `V` correspond to the elements of a 3-D coefficient array `c` of Chris@87: shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order Chris@87: Chris@87: .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},... Chris@87: Chris@87: and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the Chris@87: same up to roundoff. This equivalence is useful both for least squares Chris@87: fitting and for the evaluation of a large number of 3-D polynomials Chris@87: of the same degrees and sample points. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x, y, z : array_like Chris@87: Arrays of point coordinates, all of the same shape. The dtypes will Chris@87: be converted to either float64 or complex128 depending on whether Chris@87: any of the elements are complex. Scalars are converted to 1-D Chris@87: arrays. Chris@87: deg : list of ints Chris@87: List of maximum degrees of the form [x_deg, y_deg, z_deg]. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: vander3d : ndarray Chris@87: The shape of the returned matrix is ``x.shape + (order,)``, where Chris@87: :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will Chris@87: be the same as the converted `x`, `y`, and `z`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: polyvander, polyvander3d. polyval2d, polyval3d Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: ideg = [int(d) for d in deg] Chris@87: is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] Chris@87: if is_valid != [1, 1, 1]: Chris@87: raise ValueError("degrees must be non-negative integers") Chris@87: degx, degy, degz = ideg Chris@87: x, y, z = np.array((x, y, z), copy=0) + 0.0 Chris@87: Chris@87: vx = polyvander(x, degx) Chris@87: vy = polyvander(y, degy) Chris@87: vz = polyvander(z, degz) Chris@87: v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:] Chris@87: # einsum bug Chris@87: #v = np.einsum("...i, ...j, ...k->...ijk", vx, vy, vz) Chris@87: return v.reshape(v.shape[:-3] + (-1,)) Chris@87: Chris@87: Chris@87: def polyfit(x, y, deg, rcond=None, full=False, w=None): Chris@87: """ Chris@87: Least-squares fit of a polynomial to data. Chris@87: Chris@87: Return the coefficients of a polynomial of degree `deg` that is the Chris@87: least squares fit to the data values `y` given at points `x`. If `y` is Chris@87: 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple Chris@87: fits are done, one for each column of `y`, and the resulting Chris@87: coefficients are stored in the corresponding columns of a 2-D return. Chris@87: The fitted polynomial(s) are in the form Chris@87: Chris@87: .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n, Chris@87: Chris@87: where `n` is `deg`. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like, shape (`M`,) Chris@87: x-coordinates of the `M` sample (data) points ``(x[i], y[i])``. Chris@87: y : array_like, shape (`M`,) or (`M`, `K`) Chris@87: y-coordinates of the sample points. Several sets of sample points Chris@87: sharing the same x-coordinates can be (independently) fit with one Chris@87: call to `polyfit` by passing in for `y` a 2-D array that contains Chris@87: one data set per column. Chris@87: deg : int Chris@87: Degree of the polynomial(s) to be fit. Chris@87: rcond : float, optional Chris@87: Relative condition number of the fit. Singular values smaller Chris@87: than `rcond`, relative to the largest singular value, will be Chris@87: ignored. The default value is ``len(x)*eps``, where `eps` is the Chris@87: relative precision of the platform's float type, about 2e-16 in Chris@87: most cases. Chris@87: full : bool, optional Chris@87: Switch determining the nature of the return value. When ``False`` Chris@87: (the default) just the coefficients are returned; when ``True``, Chris@87: diagnostic information from the singular value decomposition (used Chris@87: to solve the fit's matrix equation) is also returned. Chris@87: w : array_like, shape (`M`,), optional Chris@87: Weights. If not None, the contribution of each point Chris@87: ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the Chris@87: weights are chosen so that the errors of the products ``w[i]*y[i]`` Chris@87: all have the same variance. The default value is None. Chris@87: Chris@87: .. versionadded:: 1.5.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`) Chris@87: Polynomial coefficients ordered from low to high. If `y` was 2-D, Chris@87: the coefficients in column `k` of `coef` represent the polynomial Chris@87: fit to the data in `y`'s `k`-th column. Chris@87: Chris@87: [residuals, rank, singular_values, rcond] : list Chris@87: These values are only returned if `full` = True Chris@87: Chris@87: resid -- sum of squared residuals of the least squares fit Chris@87: rank -- the numerical rank of the scaled Vandermonde matrix Chris@87: sv -- singular values of the scaled Vandermonde matrix Chris@87: rcond -- value of `rcond`. Chris@87: Chris@87: For more details, see `linalg.lstsq`. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: RankWarning Chris@87: Raised if the matrix in the least-squares fit is rank deficient. Chris@87: The warning is only raised if `full` == False. The warnings can Chris@87: be turned off by: Chris@87: Chris@87: >>> import warnings Chris@87: >>> warnings.simplefilter('ignore', RankWarning) Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: chebfit, legfit, lagfit, hermfit, hermefit Chris@87: polyval : Evaluates a polynomial. Chris@87: polyvander : Vandermonde matrix for powers. Chris@87: linalg.lstsq : Computes a least-squares fit from the matrix. Chris@87: scipy.interpolate.UnivariateSpline : Computes spline fits. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The solution is the coefficients of the polynomial `p` that minimizes Chris@87: the sum of the weighted squared errors Chris@87: Chris@87: .. math :: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2, Chris@87: Chris@87: where the :math:`w_j` are the weights. This problem is solved by Chris@87: setting up the (typically) over-determined matrix equation: Chris@87: Chris@87: .. math :: V(x) * c = w * y, Chris@87: Chris@87: where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the Chris@87: coefficients to be solved for, `w` are the weights, and `y` are the Chris@87: observed values. This equation is then solved using the singular value Chris@87: decomposition of `V`. Chris@87: Chris@87: If some of the singular values of `V` are so small that they are Chris@87: neglected (and `full` == ``False``), a `RankWarning` will be raised. Chris@87: This means that the coefficient values may be poorly determined. Chris@87: Fitting to a lower order polynomial will usually get rid of the warning Chris@87: (but may not be what you want, of course; if you have independent Chris@87: reason(s) for choosing the degree which isn't working, you may have to: Chris@87: a) reconsider those reasons, and/or b) reconsider the quality of your Chris@87: data). The `rcond` parameter can also be set to a value smaller than Chris@87: its default, but the resulting fit may be spurious and have large Chris@87: contributions from roundoff error. Chris@87: Chris@87: Polynomial fits using double precision tend to "fail" at about Chris@87: (polynomial) degree 20. Fits using Chebyshev or Legendre series are Chris@87: generally better conditioned, but much can still depend on the Chris@87: distribution of the sample points and the smoothness of the data. If Chris@87: the quality of the fit is inadequate, splines may be a good Chris@87: alternative. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> from numpy.polynomial import polynomial as P Chris@87: >>> x = np.linspace(-1,1,51) # x "data": [-1, -0.96, ..., 0.96, 1] Chris@87: >>> y = x**3 - x + np.random.randn(len(x)) # x^3 - x + N(0,1) "noise" Chris@87: >>> c, stats = P.polyfit(x,y,3,full=True) Chris@87: >>> c # c[0], c[2] should be approx. 0, c[1] approx. -1, c[3] approx. 1 Chris@87: array([ 0.01909725, -1.30598256, -0.00577963, 1.02644286]) Chris@87: >>> stats # note the large SSR, explaining the rather poor results Chris@87: [array([ 38.06116253]), 4, array([ 1.38446749, 1.32119158, 0.50443316, Chris@87: 0.28853036]), 1.1324274851176597e-014] Chris@87: Chris@87: Same thing without the added noise Chris@87: Chris@87: >>> y = x**3 - x Chris@87: >>> c, stats = P.polyfit(x,y,3,full=True) Chris@87: >>> c # c[0], c[2] should be "very close to 0", c[1] ~= -1, c[3] ~= 1 Chris@87: array([ -1.73362882e-17, -1.00000000e+00, -2.67471909e-16, Chris@87: 1.00000000e+00]) Chris@87: >>> stats # note the minuscule SSR Chris@87: [array([ 7.46346754e-31]), 4, array([ 1.38446749, 1.32119158, Chris@87: 0.50443316, 0.28853036]), 1.1324274851176597e-014] Chris@87: Chris@87: """ Chris@87: order = int(deg) + 1 Chris@87: x = np.asarray(x) + 0.0 Chris@87: y = np.asarray(y) + 0.0 Chris@87: Chris@87: # check arguments. Chris@87: if deg < 0: Chris@87: raise ValueError("expected deg >= 0") Chris@87: if x.ndim != 1: Chris@87: raise TypeError("expected 1D vector for x") Chris@87: if x.size == 0: Chris@87: raise TypeError("expected non-empty vector for x") Chris@87: if y.ndim < 1 or y.ndim > 2: Chris@87: raise TypeError("expected 1D or 2D array for y") Chris@87: if len(x) != len(y): Chris@87: raise TypeError("expected x and y to have same length") Chris@87: Chris@87: # set up the least squares matrices in transposed form Chris@87: lhs = polyvander(x, deg).T Chris@87: rhs = y.T Chris@87: if w is not None: Chris@87: w = np.asarray(w) + 0.0 Chris@87: if w.ndim != 1: Chris@87: raise TypeError("expected 1D vector for w") Chris@87: if len(x) != len(w): Chris@87: raise TypeError("expected x and w to have same length") Chris@87: # apply weights. Don't use inplace operations as they Chris@87: # can cause problems with NA. Chris@87: lhs = lhs * w Chris@87: rhs = rhs * w Chris@87: Chris@87: # set rcond Chris@87: if rcond is None: Chris@87: rcond = len(x)*np.finfo(x.dtype).eps Chris@87: Chris@87: # Determine the norms of the design matrix columns. Chris@87: if issubclass(lhs.dtype.type, np.complexfloating): Chris@87: scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1)) Chris@87: else: Chris@87: scl = np.sqrt(np.square(lhs).sum(1)) Chris@87: scl[scl == 0] = 1 Chris@87: Chris@87: # Solve the least squares problem. Chris@87: c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) Chris@87: c = (c.T/scl).T Chris@87: Chris@87: # warn on rank reduction Chris@87: if rank != order and not full: Chris@87: msg = "The fit may be poorly conditioned" Chris@87: warnings.warn(msg, pu.RankWarning) Chris@87: Chris@87: if full: Chris@87: return c, [resids, rank, s, rcond] Chris@87: else: Chris@87: return c Chris@87: Chris@87: Chris@87: def polycompanion(c): Chris@87: """ Chris@87: Return the companion matrix of c. Chris@87: Chris@87: The companion matrix for power series cannot be made symmetric by Chris@87: scaling the basis, so this function differs from those for the Chris@87: orthogonal polynomials. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : array_like Chris@87: 1-D array of polynomial coefficients ordered from low to high Chris@87: degree. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: mat : ndarray Chris@87: Companion matrix of dimensions (deg, deg). Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: """ Chris@87: # c is a trimmed copy Chris@87: [c] = pu.as_series([c]) Chris@87: if len(c) < 2: Chris@87: raise ValueError('Series must have maximum degree of at least 1.') Chris@87: if len(c) == 2: Chris@87: return np.array([[-c[0]/c[1]]]) Chris@87: Chris@87: n = len(c) - 1 Chris@87: mat = np.zeros((n, n), dtype=c.dtype) Chris@87: bot = mat.reshape(-1)[n::n+1] Chris@87: bot[...] = 1 Chris@87: mat[:, -1] -= c[:-1]/c[-1] Chris@87: return mat Chris@87: Chris@87: Chris@87: def polyroots(c): Chris@87: """ Chris@87: Compute the roots of a polynomial. Chris@87: Chris@87: Return the roots (a.k.a. "zeros") of the polynomial Chris@87: Chris@87: .. math:: p(x) = \\sum_i c[i] * x^i. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: c : 1-D array_like Chris@87: 1-D array of polynomial coefficients. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Array of the roots of the polynomial. If all the roots are real, Chris@87: then `out` is also real, otherwise it is complex. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: chebroots Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The root estimates are obtained as the eigenvalues of the companion Chris@87: matrix, Roots far from the origin of the complex plane may have large Chris@87: errors due to the numerical instability of the power series for such Chris@87: values. Roots with multiplicity greater than 1 will also show larger Chris@87: errors as the value of the series near such points is relatively Chris@87: insensitive to errors in the roots. Isolated roots near the origin can Chris@87: be improved by a few iterations of Newton's method. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> import numpy.polynomial.polynomial as poly Chris@87: >>> poly.polyroots(poly.polyfromroots((-1,0,1))) Chris@87: array([-1., 0., 1.]) Chris@87: >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype Chris@87: dtype('float64') Chris@87: >>> j = complex(0,1) Chris@87: >>> poly.polyroots(poly.polyfromroots((-j,0,j))) Chris@87: array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j]) Chris@87: Chris@87: """ Chris@87: # c is a trimmed copy Chris@87: [c] = pu.as_series([c]) Chris@87: if len(c) < 2: Chris@87: return np.array([], dtype=c.dtype) Chris@87: if len(c) == 2: Chris@87: return np.array([-c[0]/c[1]]) Chris@87: Chris@87: m = polycompanion(c) Chris@87: r = la.eigvals(m) Chris@87: r.sort() Chris@87: return r Chris@87: Chris@87: Chris@87: # Chris@87: # polynomial class Chris@87: # Chris@87: Chris@87: class Polynomial(ABCPolyBase): Chris@87: """A power series class. Chris@87: Chris@87: The Polynomial class provides the standard Python numerical methods Chris@87: '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the Chris@87: attributes and methods listed in the `ABCPolyBase` documentation. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: coef : array_like Chris@87: Polynomial coefficients in order of increasing degree, i.e., Chris@87: ``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``. Chris@87: domain : (2,) array_like, optional Chris@87: Domain to use. The interval ``[domain[0], domain[1]]`` is mapped Chris@87: to the interval ``[window[0], window[1]]`` by shifting and scaling. Chris@87: The default value is [-1, 1]. Chris@87: window : (2,) array_like, optional Chris@87: Window, see `domain` for its use. The default value is [-1, 1]. Chris@87: Chris@87: .. versionadded:: 1.6.0 Chris@87: Chris@87: """ Chris@87: # Virtual Functions Chris@87: _add = staticmethod(polyadd) Chris@87: _sub = staticmethod(polysub) Chris@87: _mul = staticmethod(polymul) Chris@87: _div = staticmethod(polydiv) Chris@87: _pow = staticmethod(polypow) Chris@87: _val = staticmethod(polyval) Chris@87: _int = staticmethod(polyint) Chris@87: _der = staticmethod(polyder) Chris@87: _fit = staticmethod(polyfit) Chris@87: _line = staticmethod(polyline) Chris@87: _roots = staticmethod(polyroots) Chris@87: _fromroots = staticmethod(polyfromroots) Chris@87: Chris@87: # Virtual properties Chris@87: nickname = 'poly' Chris@87: domain = np.array(polydomain) Chris@87: window = np.array(polydomain)