From 8a5ed09610e56f9f60089ca017b7e4abd3ce9264 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Wed, 15 Sep 2010 09:24:31 -0600 Subject: Remove unused function legtimesx, it has been replaced by legmulx. --- numpy/polynomial/legendre.py | 46 -------------------------------------------- 1 file changed, 46 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py index 9aec256cd..3ea68a5e6 100644 --- a/numpy/polynomial/legendre.py +++ b/numpy/polynomial/legendre.py @@ -227,52 +227,6 @@ def legline(off, scl) : return np.array([off]) -def legtimesx(cs): - """Multiply a Legendre series by x. - - Multiply the Legendre series `cs` by x, where x is the independent - variable. - - - Parameters - ---------- - cs : array_like - 1-d array of Legendre series coefficients ordered from low to - high. - - Returns - ------- - out : ndarray - Array representing the result of the multiplication. - - Notes - ----- - The multiplication uses the recursion relationship for Legendre - polynomials in the form - - .. math:: - - xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) - - """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs - - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1] = cs[0] - for i in range(1, len(cs)): - j = i + 1 - k = i - 1 - s = i + j - prd[j] = (cs[i]*j)/s - prd[k] += (cs[i]*i)/s - return prd - - def legfromroots(roots) : """ Generate a Legendre series with the given roots. -- cgit v1.2.1 From 96c4eea33d8d69655be9a847d35e2ff96b9d5ffd Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Wed, 22 Sep 2010 19:33:06 -0600 Subject: ENH: First commit of hermite and laguerre polynomials. The documentation and tests still need fixes. --- numpy/polynomial/hermite.py | 1141 ++++++++++++++++++++++++++++++ numpy/polynomial/hermite_e.py | 1138 +++++++++++++++++++++++++++++ numpy/polynomial/laguerre.py | 1141 ++++++++++++++++++++++++++++++ numpy/polynomial/tests/test_hermite.py | 536 ++++++++++++++ numpy/polynomial/tests/test_hermite_e.py | 536 ++++++++++++++ numpy/polynomial/tests/test_laguerre.py | 530 ++++++++++++++ 6 files changed, 5022 insertions(+) create mode 100644 numpy/polynomial/hermite.py create mode 100644 numpy/polynomial/hermite_e.py create mode 100644 numpy/polynomial/laguerre.py create mode 100644 numpy/polynomial/tests/test_hermite.py create mode 100644 numpy/polynomial/tests/test_hermite_e.py create mode 100644 numpy/polynomial/tests/test_laguerre.py (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py new file mode 100644 index 000000000..550e316de --- /dev/null +++ b/numpy/polynomial/hermite.py @@ -0,0 +1,1141 @@ +""" +Objects for dealing with Hermite series. + +This module provides a number of objects (mostly functions) useful for +dealing with Hermite series, including a `Hermite` class that +encapsulates the usual arithmetic operations. (General information +on how this module represents and works with such polynomials is in the +docstring for its "parent" sub-package, `numpy.polynomial`). + +Constants +--------- +- `hermdomain` -- Hermite series default domain, [-1,1]. +- `hermzero` -- Hermite series that evaluates identically to 0. +- `hermone` -- Hermite series that evaluates identically to 1. +- `hermx` -- Hermite series for the identity map, ``f(x) = x``. + +Arithmetic +---------- +- `hermmulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``. +- `hermadd` -- add two Hermite series. +- `hermsub` -- subtract one Hermite series from another. +- `hermmul` -- multiply two Hermite series. +- `hermdiv` -- divide one Hermite series by another. +- `hermval` -- evaluate a Hermite series at given points. + +Calculus +-------- +- `hermder` -- differentiate a Hermite series. +- `hermint` -- integrate a Hermite series. + +Misc Functions +-------------- +- `hermfromroots` -- create a Hermite series with specified roots. +- `hermroots` -- find the roots of a Hermite series. +- `hermvander` -- Vandermonde-like matrix for Hermite polynomials. +- `hermfit` -- least-squares fit returning a Hermite series. +- `hermtrim` -- trim leading coefficients from a Hermite series. +- `hermline` -- Hermite series of given straight line. +- `herm2poly` -- convert a Hermite series to a polynomial. +- `poly2herm` -- convert a polynomial to a Hermite series. + +Classes +------- +- `Hermite` -- A Hermite series class. + +See also +-------- +`numpy.polynomial` + +""" +from __future__ import division + +__all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline', + 'hermadd', 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermval', + 'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots', + 'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite'] + +import numpy as np +import numpy.linalg as la +import polyutils as pu +import warnings +from polytemplate import polytemplate + +hermtrim = pu.trimcoef + +def poly2herm(pol) : + """ + poly2herm(pol) + + Convert a polynomial to a Hermite series. + + Convert an array representing the coefficients of a polynomial (relative + to the "standard" basis) ordered from lowest degree to highest, to an + array of the coefficients of the equivalent Hermite series, ordered + from lowest to highest degree. + + Parameters + ---------- + pol : array_like + 1-d array containing the polynomial coefficients + + Returns + ------- + cs : ndarray + 1-d array containing the coefficients of the equivalent Hermite + series. + + See Also + -------- + herm2poly + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> from numpy import polynomial as P + >>> p = P.Polynomial(np.arange(4)) + >>> p + Polynomial([ 0., 1., 2., 3.], [-1., 1.]) + >>> c = P.Hermite(P.poly2herm(p.coef)) + >>> c + Hermite([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + + """ + [pol] = pu.as_series([pol]) + deg = len(pol) - 1 + res = 0 + for i in range(deg, -1, -1) : + res = hermadd(hermmulx(res), pol[i]) + return res + + +def herm2poly(cs) : + """ + Convert a Hermite series to a polynomial. + + Convert an array representing the coefficients of a Hermite series, + ordered from lowest degree to highest, to an array of the coefficients + of the equivalent polynomial (relative to the "standard" basis) ordered + from lowest to highest degree. + + Parameters + ---------- + cs : array_like + 1-d array containing the Hermite series coefficients, ordered + from lowest order term to highest. + + Returns + ------- + pol : ndarray + 1-d array containing the coefficients of the equivalent polynomial + (relative to the "standard" basis) ordered from lowest order term + to highest. + + See Also + -------- + poly2herm + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> c = P.Hermite(range(4)) + >>> c + Hermite([ 0., 1., 2., 3.], [-1., 1.]) + >>> p = c.convert(kind=P.Polynomial) + >>> p + Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) + >>> P.herm2poly(range(4)) + array([-1. , -3.5, 3. , 7.5]) + + + """ + from polynomial import polyadd, polysub, polymulx + + [cs] = pu.as_series([cs]) + n = len(cs) + if n == 1: + return cs + if n == 2: + cs[1] *= 2 + return cs + else: + c0 = cs[-2] + c1 = cs[-1] + # i is the current degree of c1 + for i in range(n - 1, 1, -1) : + tmp = c0 + c0 = polysub(cs[i - 2], c1*(2*(i - 1))) + c1 = polyadd(tmp, polymulx(c1)*2) + return polyadd(c0, polymulx(c1)*2) + +# +# These are constant arrays are of integer type so as to be compatible +# with the widest range of other types, such as Decimal. +# + +# Hermite +hermdomain = np.array([-1,1]) + +# Hermite coefficients representing zero. +hermzero = np.array([0]) + +# Hermite coefficients representing one. +hermone = np.array([1]) + +# Hermite coefficients representing the identity x. +hermx = np.array([0, 1/2]) + + +def hermline(off, scl) : + """ + Hermite series whose graph is a straight line. + + + + Parameters + ---------- + off, scl : scalars + The specified line is given by ``off + scl*x``. + + Returns + ------- + y : ndarray + This module's representation of the Hermite series for + ``off + scl*x``. + + See Also + -------- + polyline, chebline + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.hermline(3,2) + array([3, 2]) + >>> L.hermval(-3, L.hermline(3,2)) # should be -3 + -3.0 + + """ + if scl != 0 : + return np.array([off,scl/2]) + else : + return np.array([off]) + + +def hermfromroots(roots) : + """ + Generate a Hermite series with the given roots. + + Return the array of coefficients for the P-series whose roots (a.k.a. + "zeros") are given by *roots*. The returned array of coefficients is + ordered from lowest order "term" to highest, and zeros of multiplicity + greater than one must be included in *roots* a number of times equal + to their multiplicity (e.g., if `2` is a root of multiplicity three, + then [2,2,2] must be in *roots*). + + Parameters + ---------- + roots : array_like + Sequence containing the roots. + + Returns + ------- + out : ndarray + 1-d array of the Hermite series coefficients, ordered from low to + high. If all roots are real, ``out.dtype`` is a float type; + otherwise, ``out.dtype`` is a complex type, even if all the + coefficients in the result are real (see Examples below). + + See Also + -------- + polyfromroots, chebfromroots + + Notes + ----- + What is returned are the :math:`c_i` such that: + + .. math:: + + \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) + + where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite + (basis) polynomial over the domain `[-1,1]`. Note that, unlike + `polyfromroots`, due to the nature of the Hermite basis set, the + above identity *does not* imply :math:`c_n = 1` identically (see + Examples). + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.hermfromroots((-1,0,1)) # x^3 - x relative to the standard basis + array([ 0. , -0.4, 0. , 0.4]) + >>> j = complex(0,1) + >>> L.hermfromroots((-j,j)) # x^2 + 1 relative to the standard basis + array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + + """ + if len(roots) == 0 : + return np.ones(1) + else : + [roots] = pu.as_series([roots], trim=False) + prd = np.array([1], dtype=roots.dtype) + for r in roots: + prd = hermsub(hermmulx(prd), r*prd) + return prd + + +def hermadd(c1, c2): + """ + Add one Hermite series to another. + + Returns the sum of two Hermite series `c1` + `c2`. The arguments + are sequences of coefficients ordered from lowest order term to + highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the Hermite series of their sum. + + See Also + -------- + hermsub, hermmul, hermdiv, hermpow + + Notes + ----- + Unlike multiplication, division, etc., the sum of two Hermite series + is a Hermite series (without having to "reproject" the result onto + the basis set) so addition, just like that of "standard" polynomials, + is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermadd(c1,c2) + array([ 4., 4., 4.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] += c2 + ret = c1 + else : + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def hermsub(c1, c2): + """ + Subtract one Hermite series from another. + + Returns the difference of two Hermite series `c1` - `c2`. The + sequences of coefficients are from lowest order term to highest, i.e., + [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Hermite series coefficients representing their difference. + + See Also + -------- + hermadd, hermmul, hermdiv, hermpow + + Notes + ----- + Unlike multiplication, division, etc., the difference of two Hermite + series is a Hermite series (without having to "reproject" the result + onto the basis set) so subtraction, just like that of "standard" + polynomials, is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermsub(c1,c2) + array([-2., 0., 2.]) + >>> L.hermsub(c2,c1) # -C.hermsub(c1,c2) + array([ 2., 0., -2.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] -= c2 + ret = c1 + else : + c2 = -c2 + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def hermmulx(cs): + """Multiply a Hermite series by x. + + Multiply the Hermite series `cs` by x, where x is the independent + variable. + + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the result of the multiplication. + + Notes + ----- + The multiplication uses the recursion relationship for Hermite + polynomials in the form + + .. math:: + + xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + # The zero series needs special treatment + if len(cs) == 1 and cs[0] == 0: + return cs + + prd = np.empty(len(cs) + 1, dtype=cs.dtype) + prd[0] = cs[0]*0 + prd[1] = cs[0]/2 + for i in range(1, len(cs)): + prd[i + 1] = cs[i]/2 + prd[i - 1] += cs[i]*i + return prd + + +def hermmul(c1, c2): + """ + Multiply one Hermite series by another. + + Returns the product of two Hermite series `c1` * `c2`. The arguments + are sequences of coefficients, from lowest order "term" to highest, + e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Hermite series coefficients representing their product. + + See Also + -------- + hermadd, hermsub, hermdiv, hermpow + + Notes + ----- + In general, the (polynomial) product of two C-series results in terms + that are not in the Hermite polynomial basis set. Thus, to express + the product as a Hermite series, it is necessary to "re-project" the + product onto said basis set, which may produce "un-intuitive" (but + correct) results; see Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2) + >>> P.hermmul(c1,c2) # multiplication requires "reprojection" + array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + + """ + # s1, s2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + + if len(c1) > len(c2): + cs = c2 + xs = c1 + else: + cs = c1 + xs = c2 + + if len(cs) == 1: + c0 = cs[0]*xs + c1 = 0 + elif len(cs) == 2: + c0 = cs[0]*xs + c1 = cs[1]*xs + else : + nd = len(cs) + c0 = cs[-2]*xs + c1 = cs[-1]*xs + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = hermsub(cs[-i]*xs, c1*(2*(nd - 1))) + c1 = hermadd(tmp, hermmulx(c1)*2) + return hermadd(c0, hermmulx(c1)*2) + + +def hermdiv(c1, c2): + """ + Divide one Hermite series by another. + + Returns the quotient-with-remainder of two Hermite series + `c1` / `c2`. The arguments are sequences of coefficients from lowest + order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + [quo, rem] : ndarrays + Of Hermite series coefficients representing the quotient and + remainder. + + See Also + -------- + hermadd, hermsub, hermmul, hermpow + + Notes + ----- + In general, the (polynomial) division of one Hermite series by another + results in quotient and remainder terms that are not in the Hermite + polynomial basis set. Thus, to express these results as a Hermite + series, it is necessary to "re-project" the results onto the Hermite + basis set, which may produce "un-intuitive" (but correct) results; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermdiv(c1,c2) # quotient "intuitive," remainder not + (array([ 3.]), array([-8., -4.])) + >>> c2 = (0,1,2,3) + >>> L.hermdiv(c2,c1) # neither "intuitive" + (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if c2[-1] == 0 : + raise ZeroDivisionError() + + lc1 = len(c1) + lc2 = len(c2) + if lc1 < lc2 : + return c1[:1]*0, c1 + elif lc2 == 1 : + return c1/c2[-1], c1[:1]*0 + else : + quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) + rem = c1 + for i in range(lc1 - lc2, - 1, -1): + p = hermmul([0]*i + [1], c2) + q = rem[-1]/p[-1] + rem = rem[:-1] - q*p[:-1] + quo[i] = q + return quo, pu.trimseq(rem) + + +def hermpow(cs, pow, maxpower=16) : + """Raise a Hermite series to a power. + + Returns the Hermite series `cs` raised to the power `pow`. The + arguement `cs` is a sequence of coefficients ordered from low to high. + i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` + + Parameters + ---------- + cs : array_like + 1d array of Hermite series coefficients ordered from low to + high. + pow : integer + Power to which the series will be raised + maxpower : integer, optional + Maximum power allowed. This is mainly to limit growth of the series + to umanageable size. Default is 16 + + Returns + ------- + coef : ndarray + Hermite series of power. + + See Also + -------- + hermadd, hermsub, hermmul, hermdiv + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + power = int(pow) + if power != pow or power < 0 : + raise ValueError("Power must be a non-negative integer.") + elif maxpower is not None and power > maxpower : + raise ValueError("Power is too large") + elif power == 0 : + return np.array([1], dtype=cs.dtype) + elif power == 1 : + return cs + else : + # This can be made more efficient by using powers of two + # in the usual way. + prd = cs + for i in range(2, power + 1) : + prd = hermmul(prd, cs) + return prd + + +def hermder(cs, m=1, scl=1) : + """ + Differentiate a Hermite series. + + Returns the series `cs` differentiated `m` times. At each iteration the + result is multiplied by `scl` (the scaling factor is for use in a linear + change of variable). The argument `cs` is the sequence of coefficients + from lowest order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + cs: array_like + 1-d array of Hermite series coefficients ordered from low to high. + m : int, optional + Number of derivatives taken, must be non-negative. (Default: 1) + scl : scalar, optional + Each differentiation is multiplied by `scl`. The end result is + multiplication by ``scl**m``. This is for use in a linear change of + variable. (Default: 1) + + Returns + ------- + der : ndarray + Hermite series of the derivative. + + See Also + -------- + hermint + + Notes + ----- + In general, the result of differentiating a Hermite series does not + resemble the same operation on a power series. Thus the result of this + function may be "un-intuitive," albeit correct; see Examples section + below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3,4) + >>> L.hermder(cs) + array([ 6., 9., 20.]) + >>> L.hermder(cs,3) + array([ 60.]) + >>> L.hermder(cs,scl=-1) + array([ -6., -9., -20.]) + >>> L.hermder(cs,2,-1) + array([ 9., 60.]) + + """ + cnt = int(m) + + if cnt != m: + raise ValueError, "The order of derivation must be integer" + if cnt < 0 : + raise ValueError, "The order of derivation must be non-negative" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + elif cnt >= len(cs): + return cs[:1]*0 + else : + for i in range(cnt): + n = len(cs) - 1 + cs *= scl + der = np.empty(n, dtype=cs.dtype) + for j in range(n, 0, -1): + der[j - 1] = (2*j)*cs[j] + cs = der + return cs + + +def hermint(cs, m=1, k=[], lbnd=0, scl=1): + """ + Integrate a Hermite series. + + Returns a Hermite series that is the Hermite series `cs`, integrated + `m` times from `lbnd` to `x`. At each iteration the resulting series + is **multiplied** by `scl` and an integration constant, `k`, is added. + The scaling factor is for use in a linear change of variable. ("Buyer + beware": note that, depending on what one is doing, one may want `scl` + to be the reciprocal of what one might expect; for more information, + see the Notes section below.) The argument `cs` is a sequence of + coefficients, from lowest order Hermite series "term" to highest, + e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`. + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients, ordered from low to high. + m : int, optional + Order of integration, must be positive. (Default: 1) + k : {[], list, scalar}, optional + Integration constant(s). The value of the first integral at + ``lbnd`` is the first value in the list, the value of the second + integral at ``lbnd`` is the second value, etc. If ``k == []`` (the + default), all constants are set to zero. If ``m == 1``, a single + scalar can be given instead of a list. + lbnd : scalar, optional + The lower bound of the integral. (Default: 0) + scl : scalar, optional + Following each integration the result is *multiplied* by `scl` + before the integration constant is added. (Default: 1) + + Returns + ------- + S : ndarray + Hermite series coefficients of the integral. + + Raises + ------ + ValueError + If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or + ``np.isscalar(scl) == False``. + + See Also + -------- + hermder + + Notes + ----- + Note that the result of each integration is *multiplied* by `scl`. + Why is this important to note? Say one is making a linear change of + variable :math:`u = ax + b` in an integral relative to `x`. Then + :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` + - perhaps not what one would have first thought. + + Also note that, in general, the result of integrating a C-series needs + to be "re-projected" onto the C-series basis set. Thus, typically, + the result of this function is "un-intuitive," albeit correct; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3) + >>> L.hermint(cs) + array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermint(cs,3) + array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, + -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) + >>> L.hermint(cs, k=3) + array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermint(cs, lbnd=-2) + array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermint(cs, scl=2) + array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + + """ + cnt = int(m) + if np.isscalar(k) : + k = [k] + + if cnt != m: + raise ValueError, "The order of integration must be integer" + if cnt < 0 : + raise ValueError, "The order of integration must be non-negative" + if len(k) > cnt : + raise ValueError, "Too many integration constants" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + + k = list(k) + [0]*(cnt - len(k)) + for i in range(cnt) : + n = len(cs) + cs *= scl + if n == 1 and cs[0] == 0: + cs[0] += k[i] + else: + tmp = np.empty(n + 1, dtype=cs.dtype) + tmp[0] = cs[0]*0 + tmp[1] = cs[0]/2 + for j in range(1, n): + tmp[j + 1] = cs[j]/(2*(j + 1)) + tmp[0] += k[i] - hermval(lbnd, tmp) + cs = tmp + return cs + + +def hermval(x, cs): + """Evaluate a Hermite series. + + If `cs` is of length `n`, this function returns : + + ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)`` + + If x is a sequence or array then p(x) will have the same shape as x. + If r is a ring_like object that supports multiplication and addition + by the values in `cs`, then an object of the same type is returned. + + Parameters + ---------- + x : array_like, ring_like + Array of numbers or objects that support multiplication and + addition with themselves and with the elements of `cs`. + cs : array_like + 1-d array of Hermite coefficients ordered from low to high. + + Returns + ------- + values : ndarray, ring_like + If the return is an ndarray then it has the same shape as `x`. + + See Also + -------- + hermfit + + Examples + -------- + + Notes + ----- + The evaluation uses Clenshaw recursion, aka synthetic division. + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if isinstance(x, tuple) or isinstance(x, list) : + x = np.asarray(x) + + x2 = x*2 + if len(cs) == 1 : + c0 = cs[0] + c1 = 0 + elif len(cs) == 2 : + c0 = cs[0] + c1 = cs[1] + else : + nd = len(cs) + c0 = cs[-2] + c1 = cs[-1] + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = cs[-i] - c1*(2*(nd - 1)) + c1 = tmp + c1*x2 + return c0 + c1*x2 + + +def hermvander(x, deg) : + """Vandermonde matrix of given degree. + + Returns the Vandermonde matrix of degree `deg` and sample points `x`. + This isn't a true Vandermonde matrix because `x` can be an arbitrary + ndarray and the Hermite polynomials aren't powers. If ``V`` is the + returned matrix and `x` is a 2d array, then the elements of ``V`` are + ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Hermite polynomial + of degree ``k``. + + Parameters + ---------- + x : array_like + Array of points. The values are converted to double or complex + doubles. If x is scalar it is converted to a 1D array. + deg : integer + Degree of the resulting matrix. + + Returns + ------- + vander : Vandermonde matrix. + The shape of the returned matrix is ``x.shape + (deg+1,)``. The last + index is the degree. + + """ + ideg = int(deg) + if ideg != deg: + raise ValueError("deg must be integer") + if ideg < 0: + raise ValueError("deg must be non-negative") + + x = np.array(x, copy=0, ndmin=1) + 0.0 + v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype) + v[0] = x*0 + 1 + if ideg > 0 : + x2 = x*2 + v[1] = x2 + for i in range(2, ideg + 1) : + v[i] = (v[i-1]*x2 - v[i-2]*(2*(i - 1))) + return np.rollaxis(v, 0, v.ndim) + + +def hermfit(x, y, deg, rcond=None, full=False, w=None): + """ + Least squares fit of Hermite series to data. + + Fit a Hermite series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] * + P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of + coefficients `p` that minimises the squared error. + + Parameters + ---------- + x : array_like, shape (M,) + x-coordinates of the M sample points ``(x[i], y[i])``. + y : array_like, shape (M,) or (M, K) + y-coordinates of the sample points. Several data sets of sample + points sharing the same x-coordinates can be fitted at once by + passing in a 2D-array that contains one dataset per column. + deg : int + Degree of the fitting polynomial + rcond : float, optional + Relative condition number of the fit. Singular values smaller than + this relative to the largest singular value will be ignored. The + default value is len(x)*eps, where eps is the relative precision of + the float type, about 2e-16 in most cases. + full : bool, optional + Switch determining nature of return value. When it is False (the + default) just the coefficients are returned, when True diagnostic + information from the singular value decomposition is also returned. + w : array_like, shape (`M`,), optional + Weights. If not None, the contribution of each point + ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the + weights are chosen so that the errors of the products ``w[i]*y[i]`` + all have the same variance. The default value is None. + + Returns + ------- + coef : ndarray, shape (M,) or (M, K) + Hermite coefficients ordered from low to high. If `y` was 2-D, + the coefficients for the data in column k of `y` are in column + `k`. + + [residuals, rank, singular_values, rcond] : present when `full` = True + Residuals of the least-squares fit, the effective rank of the + scaled Vandermonde matrix and its singular values, and the + specified value of `rcond`. For more details, see `linalg.lstsq`. + + Warns + ----- + RankWarning + The rank of the coefficient matrix in the least-squares fit is + deficient. The warning is only raised if `full` = False. The + warnings can be turned off by + + >>> import warnings + >>> warnings.simplefilter('ignore', RankWarning) + + See Also + -------- + hermval : Evaluates a Hermite series. + hermvander : Vandermonde matrix of Hermite series. + polyfit : least squares fit using polynomials. + chebfit : least squares fit using Chebyshev series. + linalg.lstsq : Computes a least-squares fit from the matrix. + scipy.interpolate.UnivariateSpline : Computes spline fits. + + Notes + ----- + The solution are the coefficients ``c[i]`` of the Hermite series + ``P(x)`` that minimizes the squared error + + ``E = \\sum_j |y_j - P(x_j)|^2``. + + This problem is solved by setting up as the overdetermined matrix + equation + + ``V(x)*c = y``, + + where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are + the coefficients to be solved for, and the elements of `y` are the + observed values. This equation is then solved using the singular value + decomposition of ``V``. + + If some of the singular values of ``V`` are so small that they are + neglected, then a `RankWarning` will be issued. This means that the + coeficient values may be poorly determined. Using a lower order fit + will usually get rid of the warning. The `rcond` parameter can also be + set to a value smaller than its default, but the resulting fit may be + spurious and have large contributions from roundoff error. + + Fits using Hermite series are usually better conditioned than fits + using power series, but much can depend on the distribution of the + sample points and the smoothness of the data. If the quality of the fit + is inadequate splines may be a good alternative. + + References + ---------- + .. [1] Wikipedia, "Curve fitting", + http://en.wikipedia.org/wiki/Curve_fitting + + Examples + -------- + + """ + order = int(deg) + 1 + x = np.asarray(x) + 0.0 + y = np.asarray(y) + 0.0 + + # check arguments. + if deg < 0 : + raise ValueError, "expected deg >= 0" + if x.ndim != 1: + raise TypeError, "expected 1D vector for x" + if x.size == 0: + raise TypeError, "expected non-empty vector for x" + if y.ndim < 1 or y.ndim > 2 : + raise TypeError, "expected 1D or 2D array for y" + if len(x) != len(y): + raise TypeError, "expected x and y to have same length" + + # set up the least squares matrices + lhs = hermvander(x, deg) + rhs = y + if w is not None: + w = np.asarray(w) + 0.0 + if w.ndim != 1: + raise TypeError, "expected 1D vector for w" + if len(x) != len(w): + raise TypeError, "expected x and w to have same length" + # apply weights + if rhs.ndim == 2: + lhs *= w[:, np.newaxis] + rhs *= w[:, np.newaxis] + else: + lhs *= w[:, np.newaxis] + rhs *= w + + # set rcond + if rcond is None : + rcond = len(x)*np.finfo(x.dtype).eps + + # scale the design matrix and solve the least squares equation + scl = np.sqrt((lhs*lhs).sum(0)) + c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond) + c = (c.T/scl).T + + # warn on rank reduction + if rank != order and not full: + msg = "The fit may be poorly conditioned" + warnings.warn(msg, pu.RankWarning) + + if full : + return c, [resids, rank, s, rcond] + else : + return c + + +def hermroots(cs): + """ + Compute the roots of a Hermite series. + + Return the roots (a.k.a "zeros") of the Hermite series represented by + `cs`, which is the sequence of coefficients from lowest order "term" + to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients ordered from low to high. + + Returns + ------- + out : ndarray + Array of the roots. If all the roots are real, then so is the + dtype of ``out``; otherwise, ``out``'s dtype is complex. + + See Also + -------- + polyroots + chebroots + + Notes + ----- + Algorithm(s) used: + + Remember: because the Hermite series basis set is different from the + "standard" basis set, the results of this function *may* not be what + one is expecting. + + Examples + -------- + >>> import numpy.polynomial as P + >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots + array([-0.60582959+0.j , -0.07208521-0.63832674j, + -0.07208521+0.63832674j]) + >>> P.hermroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots + array([-0.85099543, -0.11407192, 0.51506735]) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if len(cs) <= 1 : + return np.array([], dtype=cs.dtype) + if len(cs) == 2 : + return np.array([-.5*cs[0]/cs[1]]) + + n = len(cs) - 1 + cs /= cs[-1] + cmat = np.zeros((n,n), dtype=cs.dtype) + cmat[1, 0] = .5 + for i in range(1, n): + cmat[i - 1, i] = i + if i != n - 1: + cmat[i + 1, i] = .5 + else: + cmat[:, i] -= cs[:-1]*.5 + roots = la.eigvals(cmat) + roots.sort() + return roots + + +# +# Hermite series class +# + +exec polytemplate.substitute(name='Hermite', nick='herm', domain='[-1,1]') diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py new file mode 100644 index 000000000..36e452074 --- /dev/null +++ b/numpy/polynomial/hermite_e.py @@ -0,0 +1,1138 @@ +""" +Objects for dealing with Hermite series. + +This module provides a number of objects (mostly functions) useful for +dealing with Hermite series, including a `Hermite` class that +encapsulates the usual arithmetic operations. (General information +on how this module represents and works with such polynomials is in the +docstring for its "parent" sub-package, `numpy.polynomial`). + +Constants +--------- +- `hermedomain` -- Hermite series default domain, [-1,1]. +- `hermezero` -- Hermite series that evaluates identically to 0. +- `hermeone` -- Hermite series that evaluates identically to 1. +- `hermex` -- Hermite series for the identity map, ``f(x) = x``. + +Arithmetic +---------- +- `hermemulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``. +- `hermeadd` -- add two Hermite series. +- `hermesub` -- subtract one Hermite series from another. +- `hermemul` -- multiply two Hermite series. +- `hermediv` -- divide one Hermite series by another. +- `hermeval` -- evaluate a Hermite series at given points. + +Calculus +-------- +- `hermeder` -- differentiate a Hermite series. +- `hermeint` -- integrate a Hermite series. + +Misc Functions +-------------- +- `hermefromroots` -- create a Hermite series with specified roots. +- `hermeroots` -- find the roots of a Hermite series. +- `hermevander` -- Vandermonde-like matrix for Hermite polynomials. +- `hermefit` -- least-squares fit returning a Hermite series. +- `hermetrim` -- trim leading coefficients from a Hermite series. +- `hermeline` -- Hermite series of given straight line. +- `herme2poly` -- convert a Hermite series to a polynomial. +- `poly2herme` -- convert a polynomial to a Hermite series. + +Classes +------- +- `Hermite` -- A Hermite series class. + +See also +-------- +`numpy.polynomial` + +""" +from __future__ import division + +__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', + 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval', + 'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots', + 'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'Hermite_e'] + +import numpy as np +import numpy.linalg as la +import polyutils as pu +import warnings +from polytemplate import polytemplate + +hermetrim = pu.trimcoef + +def poly2herme(pol) : + """ + poly2herme(pol) + + Convert a polynomial to a Hermite series. + + Convert an array representing the coefficients of a polynomial (relative + to the "standard" basis) ordered from lowest degree to highest, to an + array of the coefficients of the equivalent Hermite series, ordered + from lowest to highest degree. + + Parameters + ---------- + pol : array_like + 1-d array containing the polynomial coefficients + + Returns + ------- + cs : ndarray + 1-d array containing the coefficients of the equivalent Hermite + series. + + See Also + -------- + herme2poly + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> from numpy import polynomial as P + >>> p = P.Polynomial(np.arange(4)) + >>> p + Polynomial([ 0., 1., 2., 3.], [-1., 1.]) + >>> c = P.Hermite(P.poly2herme(p.coef)) + >>> c + Hermite([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + + """ + [pol] = pu.as_series([pol]) + deg = len(pol) - 1 + res = 0 + for i in range(deg, -1, -1) : + res = hermeadd(hermemulx(res), pol[i]) + return res + + +def herme2poly(cs) : + """ + Convert a Hermite series to a polynomial. + + Convert an array representing the coefficients of a Hermite series, + ordered from lowest degree to highest, to an array of the coefficients + of the equivalent polynomial (relative to the "standard" basis) ordered + from lowest to highest degree. + + Parameters + ---------- + cs : array_like + 1-d array containing the Hermite series coefficients, ordered + from lowest order term to highest. + + Returns + ------- + pol : ndarray + 1-d array containing the coefficients of the equivalent polynomial + (relative to the "standard" basis) ordered from lowest order term + to highest. + + See Also + -------- + poly2herme + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> c = P.Hermite(range(4)) + >>> c + Hermite([ 0., 1., 2., 3.], [-1., 1.]) + >>> p = c.convert(kind=P.Polynomial) + >>> p + Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) + >>> P.herme2poly(range(4)) + array([-1. , -3.5, 3. , 7.5]) + + + """ + from polynomial import polyadd, polysub, polymulx + + [cs] = pu.as_series([cs]) + n = len(cs) + if n == 1: + return cs + if n == 2: + return cs + else: + c0 = cs[-2] + c1 = cs[-1] + # i is the current degree of c1 + for i in range(n - 1, 1, -1) : + tmp = c0 + c0 = polysub(cs[i - 2], c1*(i - 1)) + c1 = polyadd(tmp, polymulx(c1)) + return polyadd(c0, polymulx(c1)) + +# +# These are constant arrays are of integer type so as to be compatible +# with the widest range of other types, such as Decimal. +# + +# Hermite +hermedomain = np.array([-1,1]) + +# Hermite coefficients representing zero. +hermezero = np.array([0]) + +# Hermite coefficients representing one. +hermeone = np.array([1]) + +# Hermite coefficients representing the identity x. +hermex = np.array([0, 1]) + + +def hermeline(off, scl) : + """ + Hermite series whose graph is a straight line. + + + + Parameters + ---------- + off, scl : scalars + The specified line is given by ``off + scl*x``. + + Returns + ------- + y : ndarray + This module's representation of the Hermite series for + ``off + scl*x``. + + See Also + -------- + polyline, chebline + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.hermeline(3,2) + array([3, 2]) + >>> L.hermeval(-3, L.hermeline(3,2)) # should be -3 + -3.0 + + """ + if scl != 0 : + return np.array([off,scl]) + else : + return np.array([off]) + + +def hermefromroots(roots) : + """ + Generate a Hermite series with the given roots. + + Return the array of coefficients for the P-series whose roots (a.k.a. + "zeros") are given by *roots*. The returned array of coefficients is + ordered from lowest order "term" to highest, and zeros of multiplicity + greater than one must be included in *roots* a number of times equal + to their multiplicity (e.g., if `2` is a root of multiplicity three, + then [2,2,2] must be in *roots*). + + Parameters + ---------- + roots : array_like + Sequence containing the roots. + + Returns + ------- + out : ndarray + 1-d array of the Hermite series coefficients, ordered from low to + high. If all roots are real, ``out.dtype`` is a float type; + otherwise, ``out.dtype`` is a complex type, even if all the + coefficients in the result are real (see Examples below). + + See Also + -------- + polyfromroots, chebfromroots + + Notes + ----- + What is returned are the :math:`c_i` such that: + + .. math:: + + \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) + + where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite + (basis) polynomial over the domain `[-1,1]`. Note that, unlike + `polyfromroots`, due to the nature of the Hermite basis set, the + above identity *does not* imply :math:`c_n = 1` identically (see + Examples). + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.hermefromroots((-1,0,1)) # x^3 - x relative to the standard basis + array([ 0. , -0.4, 0. , 0.4]) + >>> j = complex(0,1) + >>> L.hermefromroots((-j,j)) # x^2 + 1 relative to the standard basis + array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + + """ + if len(roots) == 0 : + return np.ones(1) + else : + [roots] = pu.as_series([roots], trim=False) + prd = np.array([1], dtype=roots.dtype) + for r in roots: + prd = hermesub(hermemulx(prd), r*prd) + return prd + + +def hermeadd(c1, c2): + """ + Add one Hermite series to another. + + Returns the sum of two Hermite series `c1` + `c2`. The arguments + are sequences of coefficients ordered from lowest order term to + highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the Hermite series of their sum. + + See Also + -------- + hermesub, hermemul, hermediv, hermepow + + Notes + ----- + Unlike multiplication, division, etc., the sum of two Hermite series + is a Hermite series (without having to "reproject" the result onto + the basis set) so addition, just like that of "standard" polynomials, + is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermeadd(c1,c2) + array([ 4., 4., 4.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] += c2 + ret = c1 + else : + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def hermesub(c1, c2): + """ + Subtract one Hermite series from another. + + Returns the difference of two Hermite series `c1` - `c2`. The + sequences of coefficients are from lowest order term to highest, i.e., + [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Hermite series coefficients representing their difference. + + See Also + -------- + hermeadd, hermemul, hermediv, hermepow + + Notes + ----- + Unlike multiplication, division, etc., the difference of two Hermite + series is a Hermite series (without having to "reproject" the result + onto the basis set) so subtraction, just like that of "standard" + polynomials, is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermesub(c1,c2) + array([-2., 0., 2.]) + >>> L.hermesub(c2,c1) # -C.hermesub(c1,c2) + array([ 2., 0., -2.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] -= c2 + ret = c1 + else : + c2 = -c2 + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def hermemulx(cs): + """Multiply a Hermite series by x. + + Multiply the Hermite series `cs` by x, where x is the independent + variable. + + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the result of the multiplication. + + Notes + ----- + The multiplication uses the recursion relationship for Hermite + polynomials in the form + + .. math:: + + xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + # The zero series needs special treatment + if len(cs) == 1 and cs[0] == 0: + return cs + + prd = np.empty(len(cs) + 1, dtype=cs.dtype) + prd[0] = cs[0]*0 + prd[1] = cs[0]/2 + for i in range(1, len(cs)): + prd[i + 1] = cs[i] + prd[i - 1] += cs[i]*i + return prd + + +def hermemul(c1, c2): + """ + Multiply one Hermite series by another. + + Returns the product of two Hermite series `c1` * `c2`. The arguments + are sequences of coefficients, from lowest order "term" to highest, + e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Hermite series coefficients representing their product. + + See Also + -------- + hermeadd, hermesub, hermediv, hermepow + + Notes + ----- + In general, the (polynomial) product of two C-series results in terms + that are not in the Hermite polynomial basis set. Thus, to express + the product as a Hermite series, it is necessary to "re-project" the + product onto said basis set, which may produce "un-intuitive" (but + correct) results; see Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2) + >>> P.hermemul(c1,c2) # multiplication requires "reprojection" + array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + + """ + # s1, s2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + + if len(c1) > len(c2): + cs = c2 + xs = c1 + else: + cs = c1 + xs = c2 + + if len(cs) == 1: + c0 = cs[0]*xs + c1 = 0 + elif len(cs) == 2: + c0 = cs[0]*xs + c1 = cs[1]*xs + else : + nd = len(cs) + c0 = cs[-2]*xs + c1 = cs[-1]*xs + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = hermesub(cs[-i]*xs, c1*(nd - 1)) + c1 = hermeadd(tmp, hermemulx(c1)) + return hermeadd(c0, hermemulx(c1)) + + +def hermediv(c1, c2): + """ + Divide one Hermite series by another. + + Returns the quotient-with-remainder of two Hermite series + `c1` / `c2`. The arguments are sequences of coefficients from lowest + order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Hermite series coefficients ordered from low to + high. + + Returns + ------- + [quo, rem] : ndarrays + Of Hermite series coefficients representing the quotient and + remainder. + + See Also + -------- + hermeadd, hermesub, hermemul, hermepow + + Notes + ----- + In general, the (polynomial) division of one Hermite series by another + results in quotient and remainder terms that are not in the Hermite + polynomial basis set. Thus, to express these results as a Hermite + series, it is necessary to "re-project" the results onto the Hermite + basis set, which may produce "un-intuitive" (but correct) results; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.hermediv(c1,c2) # quotient "intuitive," remainder not + (array([ 3.]), array([-8., -4.])) + >>> c2 = (0,1,2,3) + >>> L.hermediv(c2,c1) # neither "intuitive" + (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if c2[-1] == 0 : + raise ZeroDivisionError() + + lc1 = len(c1) + lc2 = len(c2) + if lc1 < lc2 : + return c1[:1]*0, c1 + elif lc2 == 1 : + return c1/c2[-1], c1[:1]*0 + else : + quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) + rem = c1 + for i in range(lc1 - lc2, - 1, -1): + p = hermemul([0]*i + [1], c2) + q = rem[-1]/p[-1] + rem = rem[:-1] - q*p[:-1] + quo[i] = q + return quo, pu.trimseq(rem) + + +def hermepow(cs, pow, maxpower=16) : + """Raise a Hermite series to a power. + + Returns the Hermite series `cs` raised to the power `pow`. The + arguement `cs` is a sequence of coefficients ordered from low to high. + i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` + + Parameters + ---------- + cs : array_like + 1d array of Hermite series coefficients ordered from low to + high. + pow : integer + Power to which the series will be raised + maxpower : integer, optional + Maximum power allowed. This is mainly to limit growth of the series + to umanageable size. Default is 16 + + Returns + ------- + coef : ndarray + Hermite series of power. + + See Also + -------- + hermeadd, hermesub, hermemul, hermediv + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + power = int(pow) + if power != pow or power < 0 : + raise ValueError("Power must be a non-negative integer.") + elif maxpower is not None and power > maxpower : + raise ValueError("Power is too large") + elif power == 0 : + return np.array([1], dtype=cs.dtype) + elif power == 1 : + return cs + else : + # This can be made more efficient by using powers of two + # in the usual way. + prd = cs + for i in range(2, power + 1) : + prd = hermemul(prd, cs) + return prd + + +def hermeder(cs, m=1, scl=1) : + """ + Differentiate a Hermite series. + + Returns the series `cs` differentiated `m` times. At each iteration the + result is multiplied by `scl` (the scaling factor is for use in a linear + change of variable). The argument `cs` is the sequence of coefficients + from lowest order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + cs: array_like + 1-d array of Hermite series coefficients ordered from low to high. + m : int, optional + Number of derivatives taken, must be non-negative. (Default: 1) + scl : scalar, optional + Each differentiation is multiplied by `scl`. The end result is + multiplication by ``scl**m``. This is for use in a linear change of + variable. (Default: 1) + + Returns + ------- + der : ndarray + Hermite series of the derivative. + + See Also + -------- + hermeint + + Notes + ----- + In general, the result of differentiating a Hermite series does not + resemble the same operation on a power series. Thus the result of this + function may be "un-intuitive," albeit correct; see Examples section + below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3,4) + >>> L.hermeder(cs) + array([ 6., 9., 20.]) + >>> L.hermeder(cs,3) + array([ 60.]) + >>> L.hermeder(cs,scl=-1) + array([ -6., -9., -20.]) + >>> L.hermeder(cs,2,-1) + array([ 9., 60.]) + + """ + cnt = int(m) + + if cnt != m: + raise ValueError, "The order of derivation must be integer" + if cnt < 0 : + raise ValueError, "The order of derivation must be non-negative" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + elif cnt >= len(cs): + return cs[:1]*0 + else : + for i in range(cnt): + n = len(cs) - 1 + cs *= scl + der = np.empty(n, dtype=cs.dtype) + for j in range(n, 0, -1): + der[j - 1] = j*cs[j] + cs = der + return cs + + +def hermeint(cs, m=1, k=[], lbnd=0, scl=1): + """ + Integrate a Hermite series. + + Returns a Hermite series that is the Hermite series `cs`, integrated + `m` times from `lbnd` to `x`. At each iteration the resulting series + is **multiplied** by `scl` and an integration constant, `k`, is added. + The scaling factor is for use in a linear change of variable. ("Buyer + beware": note that, depending on what one is doing, one may want `scl` + to be the reciprocal of what one might expect; for more information, + see the Notes section below.) The argument `cs` is a sequence of + coefficients, from lowest order Hermite series "term" to highest, + e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`. + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients, ordered from low to high. + m : int, optional + Order of integration, must be positive. (Default: 1) + k : {[], list, scalar}, optional + Integration constant(s). The value of the first integral at + ``lbnd`` is the first value in the list, the value of the second + integral at ``lbnd`` is the second value, etc. If ``k == []`` (the + default), all constants are set to zero. If ``m == 1``, a single + scalar can be given instead of a list. + lbnd : scalar, optional + The lower bound of the integral. (Default: 0) + scl : scalar, optional + Following each integration the result is *multiplied* by `scl` + before the integration constant is added. (Default: 1) + + Returns + ------- + S : ndarray + Hermite series coefficients of the integral. + + Raises + ------ + ValueError + If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or + ``np.isscalar(scl) == False``. + + See Also + -------- + hermeder + + Notes + ----- + Note that the result of each integration is *multiplied* by `scl`. + Why is this important to note? Say one is making a linear change of + variable :math:`u = ax + b` in an integral relative to `x`. Then + :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` + - perhaps not what one would have first thought. + + Also note that, in general, the result of integrating a C-series needs + to be "re-projected" onto the C-series basis set. Thus, typically, + the result of this function is "un-intuitive," albeit correct; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3) + >>> L.hermeint(cs) + array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermeint(cs,3) + array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, + -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) + >>> L.hermeint(cs, k=3) + array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermeint(cs, lbnd=-2) + array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.hermeint(cs, scl=2) + array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + + """ + cnt = int(m) + if np.isscalar(k) : + k = [k] + + if cnt != m: + raise ValueError, "The order of integration must be integer" + if cnt < 0 : + raise ValueError, "The order of integration must be non-negative" + if len(k) > cnt : + raise ValueError, "Too many integration constants" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + + k = list(k) + [0]*(cnt - len(k)) + for i in range(cnt) : + n = len(cs) + cs *= scl + if n == 1 and cs[0] == 0: + cs[0] += k[i] + else: + tmp = np.empty(n + 1, dtype=cs.dtype) + tmp[0] = cs[0]*0 + tmp[1] = cs[0] + for j in range(1, n): + tmp[j + 1] = cs[j]/(j + 1) + tmp[0] += k[i] - hermeval(lbnd, tmp) + cs = tmp + return cs + + +def hermeval(x, cs): + """Evaluate a Hermite series. + + If `cs` is of length `n`, this function returns : + + ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)`` + + If x is a sequence or array then p(x) will have the same shape as x. + If r is a ring_like object that supports multiplication and addition + by the values in `cs`, then an object of the same type is returned. + + Parameters + ---------- + x : array_like, ring_like + Array of numbers or objects that support multiplication and + addition with themselves and with the elements of `cs`. + cs : array_like + 1-d array of Hermite coefficients ordered from low to high. + + Returns + ------- + values : ndarray, ring_like + If the return is an ndarray then it has the same shape as `x`. + + See Also + -------- + hermefit + + Examples + -------- + + Notes + ----- + The evaluation uses Clenshaw recursion, aka synthetic division. + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if isinstance(x, tuple) or isinstance(x, list) : + x = np.asarray(x) + + if len(cs) == 1 : + c0 = cs[0] + c1 = 0 + elif len(cs) == 2 : + c0 = cs[0] + c1 = cs[1] + else : + nd = len(cs) + c0 = cs[-2] + c1 = cs[-1] + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = cs[-i] - c1*(nd - 1) + c1 = tmp + c1*x + return c0 + c1*x + + +def hermevander(x, deg) : + """Vandermonde matrix of given degree. + + Returns the Vandermonde matrix of degree `deg` and sample points `x`. + This isn't a true Vandermonde matrix because `x` can be an arbitrary + ndarray and the Hermite polynomials aren't powers. If ``V`` is the + returned matrix and `x` is a 2d array, then the elements of ``V`` are + ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Hermite polynomial + of degree ``k``. + + Parameters + ---------- + x : array_like + Array of points. The values are converted to double or complex + doubles. If x is scalar it is converted to a 1D array. + deg : integer + Degree of the resulting matrix. + + Returns + ------- + vander : Vandermonde matrix. + The shape of the returned matrix is ``x.shape + (deg+1,)``. The last + index is the degree. + + """ + ideg = int(deg) + if ideg != deg: + raise ValueError("deg must be integer") + if ideg < 0: + raise ValueError("deg must be non-negative") + + x = np.array(x, copy=0, ndmin=1) + 0.0 + v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype) + v[0] = x*0 + 1 + if ideg > 0 : + v[1] = x + for i in range(2, ideg + 1) : + v[i] = (v[i-1]*x - v[i-2]*(i - 1)) + return np.rollaxis(v, 0, v.ndim) + + +def hermefit(x, y, deg, rcond=None, full=False, w=None): + """ + Least squares fit of Hermite series to data. + + Fit a Hermite series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] * + P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of + coefficients `p` that minimises the squared error. + + Parameters + ---------- + x : array_like, shape (M,) + x-coordinates of the M sample points ``(x[i], y[i])``. + y : array_like, shape (M,) or (M, K) + y-coordinates of the sample points. Several data sets of sample + points sharing the same x-coordinates can be fitted at once by + passing in a 2D-array that contains one dataset per column. + deg : int + Degree of the fitting polynomial + rcond : float, optional + Relative condition number of the fit. Singular values smaller than + this relative to the largest singular value will be ignored. The + default value is len(x)*eps, where eps is the relative precision of + the float type, about 2e-16 in most cases. + full : bool, optional + Switch determining nature of return value. When it is False (the + default) just the coefficients are returned, when True diagnostic + information from the singular value decomposition is also returned. + w : array_like, shape (`M`,), optional + Weights. If not None, the contribution of each point + ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the + weights are chosen so that the errors of the products ``w[i]*y[i]`` + all have the same variance. The default value is None. + + Returns + ------- + coef : ndarray, shape (M,) or (M, K) + Hermite coefficients ordered from low to high. If `y` was 2-D, + the coefficients for the data in column k of `y` are in column + `k`. + + [residuals, rank, singular_values, rcond] : present when `full` = True + Residuals of the least-squares fit, the effective rank of the + scaled Vandermonde matrix and its singular values, and the + specified value of `rcond`. For more details, see `linalg.lstsq`. + + Warns + ----- + RankWarning + The rank of the coefficient matrix in the least-squares fit is + deficient. The warning is only raised if `full` = False. The + warnings can be turned off by + + >>> import warnings + >>> warnings.simplefilter('ignore', RankWarning) + + See Also + -------- + hermeval : Evaluates a Hermite series. + hermevander : Vandermonde matrix of Hermite series. + polyfit : least squares fit using polynomials. + chebfit : least squares fit using Chebyshev series. + linalg.lstsq : Computes a least-squares fit from the matrix. + scipy.interpolate.UnivariateSpline : Computes spline fits. + + Notes + ----- + The solution are the coefficients ``c[i]`` of the Hermite series + ``P(x)`` that minimizes the squared error + + ``E = \\sum_j |y_j - P(x_j)|^2``. + + This problem is solved by setting up as the overdetermined matrix + equation + + ``V(x)*c = y``, + + where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are + the coefficients to be solved for, and the elements of `y` are the + observed values. This equation is then solved using the singular value + decomposition of ``V``. + + If some of the singular values of ``V`` are so small that they are + neglected, then a `RankWarning` will be issued. This means that the + coeficient values may be poorly determined. Using a lower order fit + will usually get rid of the warning. The `rcond` parameter can also be + set to a value smaller than its default, but the resulting fit may be + spurious and have large contributions from roundoff error. + + Fits using Hermite series are usually better conditioned than fits + using power series, but much can depend on the distribution of the + sample points and the smoothness of the data. If the quality of the fit + is inadequate splines may be a good alternative. + + References + ---------- + .. [1] Wikipedia, "Curve fitting", + http://en.wikipedia.org/wiki/Curve_fitting + + Examples + -------- + + """ + order = int(deg) + 1 + x = np.asarray(x) + 0.0 + y = np.asarray(y) + 0.0 + + # check arguments. + if deg < 0 : + raise ValueError, "expected deg >= 0" + if x.ndim != 1: + raise TypeError, "expected 1D vector for x" + if x.size == 0: + raise TypeError, "expected non-empty vector for x" + if y.ndim < 1 or y.ndim > 2 : + raise TypeError, "expected 1D or 2D array for y" + if len(x) != len(y): + raise TypeError, "expected x and y to have same length" + + # set up the least squares matrices + lhs = hermevander(x, deg) + rhs = y + if w is not None: + w = np.asarray(w) + 0.0 + if w.ndim != 1: + raise TypeError, "expected 1D vector for w" + if len(x) != len(w): + raise TypeError, "expected x and w to have same length" + # apply weights + if rhs.ndim == 2: + lhs *= w[:, np.newaxis] + rhs *= w[:, np.newaxis] + else: + lhs *= w[:, np.newaxis] + rhs *= w + + # set rcond + if rcond is None : + rcond = len(x)*np.finfo(x.dtype).eps + + # scale the design matrix and solve the least squares equation + scl = np.sqrt((lhs*lhs).sum(0)) + c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond) + c = (c.T/scl).T + + # warn on rank reduction + if rank != order and not full: + msg = "The fit may be poorly conditioned" + warnings.warn(msg, pu.RankWarning) + + if full : + return c, [resids, rank, s, rcond] + else : + return c + + +def hermeroots(cs): + """ + Compute the roots of a Hermite series. + + Return the roots (a.k.a "zeros") of the Hermite series represented by + `cs`, which is the sequence of coefficients from lowest order "term" + to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + + Parameters + ---------- + cs : array_like + 1-d array of Hermite series coefficients ordered from low to high. + + Returns + ------- + out : ndarray + Array of the roots. If all the roots are real, then so is the + dtype of ``out``; otherwise, ``out``'s dtype is complex. + + See Also + -------- + polyroots + chebroots + + Notes + ----- + Algorithm(s) used: + + Remember: because the Hermite series basis set is different from the + "standard" basis set, the results of this function *may* not be what + one is expecting. + + Examples + -------- + >>> import numpy.polynomial as P + >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots + array([-0.60582959+0.j , -0.07208521-0.63832674j, + -0.07208521+0.63832674j]) + >>> P.hermeroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots + array([-0.85099543, -0.11407192, 0.51506735]) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if len(cs) <= 1 : + return np.array([], dtype=cs.dtype) + if len(cs) == 2 : + return np.array([-.5*cs[0]/cs[1]]) + + n = len(cs) - 1 + cs /= cs[-1] + cmat = np.zeros((n,n), dtype=cs.dtype) + cmat[1, 0] = 1 + for i in range(1, n): + cmat[i - 1, i] = i + if i != n - 1: + cmat[i + 1, i] = 1 + else: + cmat[:, i] -= cs[:-1] + roots = la.eigvals(cmat) + roots.sort() + return roots + + +# +# Hermite_e series class +# + +exec polytemplate.substitute(name='Hermite_e', nick='herme', domain='[-1,1]') diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py new file mode 100644 index 000000000..94c495deb --- /dev/null +++ b/numpy/polynomial/laguerre.py @@ -0,0 +1,1141 @@ +""" +Objects for dealing with Laguerre series. + +This module provides a number of objects (mostly functions) useful for +dealing with Laguerre series, including a `Laguerre` class that +encapsulates the usual arithmetic operations. (General information +on how this module represents and works with such polynomials is in the +docstring for its "parent" sub-package, `numpy.polynomial`). + +Constants +--------- +- `lagdomain` -- Laguerre series default domain, [-1,1]. +- `lagzero` -- Laguerre series that evaluates identically to 0. +- `lagone` -- Laguerre series that evaluates identically to 1. +- `lagx` -- Laguerre series for the identity map, ``f(x) = x``. + +Arithmetic +---------- +- `lagmulx` -- multiply a Laguerre series in ``P_i(x)`` by ``x``. +- `lagadd` -- add two Laguerre series. +- `lagsub` -- subtract one Laguerre series from another. +- `lagmul` -- multiply two Laguerre series. +- `lagdiv` -- divide one Laguerre series by another. +- `lagval` -- evaluate a Laguerre series at given points. + +Calculus +-------- +- `lagder` -- differentiate a Laguerre series. +- `lagint` -- integrate a Laguerre series. + +Misc Functions +-------------- +- `lagfromroots` -- create a Laguerre series with specified roots. +- `lagroots` -- find the roots of a Laguerre series. +- `lagvander` -- Vandermonde-like matrix for Laguerre polynomials. +- `lagfit` -- least-squares fit returning a Laguerre series. +- `lagtrim` -- trim leading coefficients from a Laguerre series. +- `lagline` -- Laguerre series of given straight line. +- `lag2poly` -- convert a Laguerre series to a polynomial. +- `poly2lag` -- convert a polynomial to a Laguerre series. + +Classes +------- +- `Laguerre` -- A Laguerre series class. + +See also +-------- +`numpy.polynomial` + +""" +from __future__ import division + +__all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline', + 'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagval', + 'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots', + 'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre'] + +import numpy as np +import numpy.linalg as la +import polyutils as pu +import warnings +from polytemplate import polytemplate + +lagtrim = pu.trimcoef + +def poly2lag(pol) : + """ + poly2lag(pol) + + Convert a polynomial to a Laguerre series. + + Convert an array representing the coefficients of a polynomial (relative + to the "standard" basis) ordered from lowest degree to highest, to an + array of the coefficients of the equivalent Laguerre series, ordered + from lowest to highest degree. + + Parameters + ---------- + pol : array_like + 1-d array containing the polynomial coefficients + + Returns + ------- + cs : ndarray + 1-d array containing the coefficients of the equivalent Laguerre + series. + + See Also + -------- + lag2poly + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> from numpy import polynomial as P + >>> p = P.Polynomial(np.arange(4)) + >>> p + Polynomial([ 0., 1., 2., 3.], [-1., 1.]) + >>> c = P.Laguerre(P.poly2lag(p.coef)) + >>> c + Laguerre([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + + """ + [pol] = pu.as_series([pol]) + deg = len(pol) - 1 + res = 0 + for i in range(deg, -1, -1) : + res = lagadd(lagmulx(res), pol[i]) + return res + + +def lag2poly(cs) : + """ + Convert a Laguerre series to a polynomial. + + Convert an array representing the coefficients of a Laguerre series, + ordered from lowest degree to highest, to an array of the coefficients + of the equivalent polynomial (relative to the "standard" basis) ordered + from lowest to highest degree. + + Parameters + ---------- + cs : array_like + 1-d array containing the Laguerre series coefficients, ordered + from lowest order term to highest. + + Returns + ------- + pol : ndarray + 1-d array containing the coefficients of the equivalent polynomial + (relative to the "standard" basis) ordered from lowest order term + to highest. + + See Also + -------- + poly2lag + + Notes + ----- + The easy way to do conversions between polynomial basis sets + is to use the convert method of a class instance. + + Examples + -------- + >>> c = P.Laguerre(range(4)) + >>> c + Laguerre([ 0., 1., 2., 3.], [-1., 1.]) + >>> p = c.convert(kind=P.Polynomial) + >>> p + Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) + >>> P.lag2poly(range(4)) + array([-1. , -3.5, 3. , 7.5]) + + + """ + from polynomial import polyadd, polysub, polymulx + + [cs] = pu.as_series([cs]) + n = len(cs) + if n == 1: + return cs + else: + c0 = cs[-2] + c1 = cs[-1] + # i is the current degree of c1 + for i in range(n - 1, 1, -1): + tmp = c0 + c0 = polysub(cs[i - 2], (c1*(i - 1))/i) + c1 = polyadd(tmp, polysub((2*i - 1)*c1, polymulx(c1))/i) + return polyadd(c0, polysub(c1, polymulx(c1))) + +# +# These are constant arrays are of integer type so as to be compatible +# with the widest range of other types, such as Decimal. +# + +# Laguerre +lagdomain = np.array([0,1]) + +# Laguerre coefficients representing zero. +lagzero = np.array([0]) + +# Laguerre coefficients representing one. +lagone = np.array([1]) + +# Laguerre coefficients representing the identity x. +lagx = np.array([1, -1]) + + +def lagline(off, scl) : + """ + Laguerre series whose graph is a straight line. + + + + Parameters + ---------- + off, scl : scalars + The specified line is given by ``off + scl*x``. + + Returns + ------- + y : ndarray + This module's representation of the Laguerre series for + ``off + scl*x``. + + See Also + -------- + polyline, chebline + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.lagline(3,2) + array([3, 2]) + >>> L.lagval(-3, L.lagline(3,2)) # should be -3 + -3.0 + + """ + if scl != 0 : + return np.array([off + scl, -scl]) + else : + return np.array([off]) + + +def lagfromroots(roots) : + """ + Generate a Laguerre series with the given roots. + + Return the array of coefficients for the P-series whose roots (a.k.a. + "zeros") are given by *roots*. The returned array of coefficients is + ordered from lowest order "term" to highest, and zeros of multiplicity + greater than one must be included in *roots* a number of times equal + to their multiplicity (e.g., if `2` is a root of multiplicity three, + then [2,2,2] must be in *roots*). + + Parameters + ---------- + roots : array_like + Sequence containing the roots. + + Returns + ------- + out : ndarray + 1-d array of the Laguerre series coefficients, ordered from low to + high. If all roots are real, ``out.dtype`` is a float type; + otherwise, ``out.dtype`` is a complex type, even if all the + coefficients in the result are real (see Examples below). + + See Also + -------- + polyfromroots, chebfromroots + + Notes + ----- + What is returned are the :math:`c_i` such that: + + .. math:: + + \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) + + where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Laguerre + (basis) polynomial over the domain `[-1,1]`. Note that, unlike + `polyfromroots`, due to the nature of the Laguerre basis set, the + above identity *does not* imply :math:`c_n = 1` identically (see + Examples). + + Examples + -------- + >>> import numpy.polynomial.legendre as L + >>> L.lagfromroots((-1,0,1)) # x^3 - x relative to the standard basis + array([ 0. , -0.4, 0. , 0.4]) + >>> j = complex(0,1) + >>> L.lagfromroots((-j,j)) # x^2 + 1 relative to the standard basis + array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + + """ + if len(roots) == 0 : + return np.ones(1) + else : + [roots] = pu.as_series([roots], trim=False) + prd = np.array([1], dtype=roots.dtype) + for r in roots: + prd = lagsub(lagmulx(prd), r*prd) + return prd + + +def lagadd(c1, c2): + """ + Add one Laguerre series to another. + + Returns the sum of two Laguerre series `c1` + `c2`. The arguments + are sequences of coefficients ordered from lowest order term to + highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Laguerre series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the Laguerre series of their sum. + + See Also + -------- + lagsub, lagmul, lagdiv, lagpow + + Notes + ----- + Unlike multiplication, division, etc., the sum of two Laguerre series + is a Laguerre series (without having to "reproject" the result onto + the basis set) so addition, just like that of "standard" polynomials, + is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.lagadd(c1,c2) + array([ 4., 4., 4.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] += c2 + ret = c1 + else : + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def lagsub(c1, c2): + """ + Subtract one Laguerre series from another. + + Returns the difference of two Laguerre series `c1` - `c2`. The + sequences of coefficients are from lowest order term to highest, i.e., + [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Laguerre series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Laguerre series coefficients representing their difference. + + See Also + -------- + lagadd, lagmul, lagdiv, lagpow + + Notes + ----- + Unlike multiplication, division, etc., the difference of two Laguerre + series is a Laguerre series (without having to "reproject" the result + onto the basis set) so subtraction, just like that of "standard" + polynomials, is simply "component-wise." + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.lagsub(c1,c2) + array([-2., 0., 2.]) + >>> L.lagsub(c2,c1) # -C.lagsub(c1,c2) + array([ 2., 0., -2.]) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if len(c1) > len(c2) : + c1[:c2.size] -= c2 + ret = c1 + else : + c2 = -c2 + c2[:c1.size] += c1 + ret = c2 + return pu.trimseq(ret) + + +def lagmulx(cs): + """Multiply a Laguerre series by x. + + Multiply the Laguerre series `cs` by x, where x is the independent + variable. + + + Parameters + ---------- + cs : array_like + 1-d array of Laguerre series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Array representing the result of the multiplication. + + Notes + ----- + The multiplication uses the recursion relationship for Laguerre + polynomials in the form + + .. math:: + + xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + # The zero series needs special treatment + if len(cs) == 1 and cs[0] == 0: + return cs + + prd = np.empty(len(cs) + 1, dtype=cs.dtype) + prd[0] = cs[0] + prd[1] = -cs[0] + for i in range(1, len(cs)): + prd[i + 1] = -cs[i]*(i + 1) + prd[i] += cs[i]*(2*i + 1) + prd[i - 1] -= cs[i]*i + return prd + + +def lagmul(c1, c2): + """ + Multiply one Laguerre series by another. + + Returns the product of two Laguerre series `c1` * `c2`. The arguments + are sequences of coefficients, from lowest order "term" to highest, + e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Laguerre series coefficients ordered from low to + high. + + Returns + ------- + out : ndarray + Of Laguerre series coefficients representing their product. + + See Also + -------- + lagadd, lagsub, lagdiv, lagpow + + Notes + ----- + In general, the (polynomial) product of two C-series results in terms + that are not in the Laguerre polynomial basis set. Thus, to express + the product as a Laguerre series, it is necessary to "re-project" the + product onto said basis set, which may produce "un-intuitive" (but + correct) results; see Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2) + >>> P.lagmul(c1,c2) # multiplication requires "reprojection" + array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + + """ + # s1, s2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + + if len(c1) > len(c2): + cs = c2 + xs = c1 + else: + cs = c1 + xs = c2 + + if len(cs) == 1: + c0 = cs[0]*xs + c1 = 0 + elif len(cs) == 2: + c0 = cs[0]*xs + c1 = cs[1]*xs + else : + nd = len(cs) + c0 = cs[-2]*xs + c1 = cs[-1]*xs + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = lagsub(cs[-i]*xs, (c1*(nd - 1))/nd) + c1 = lagadd(tmp, lagsub((2*nd - 1)*c1, lagmulx(c1))/nd) + return lagadd(c0, lagsub(c1, lagmulx(c1))) + + +def lagdiv(c1, c2): + """ + Divide one Laguerre series by another. + + Returns the quotient-with-remainder of two Laguerre series + `c1` / `c2`. The arguments are sequences of coefficients from lowest + order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + c1, c2 : array_like + 1-d arrays of Laguerre series coefficients ordered from low to + high. + + Returns + ------- + [quo, rem] : ndarrays + Of Laguerre series coefficients representing the quotient and + remainder. + + See Also + -------- + lagadd, lagsub, lagmul, lagpow + + Notes + ----- + In general, the (polynomial) division of one Laguerre series by another + results in quotient and remainder terms that are not in the Laguerre + polynomial basis set. Thus, to express these results as a Laguerre + series, it is necessary to "re-project" the results onto the Laguerre + basis set, which may produce "un-intuitive" (but correct) results; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> c1 = (1,2,3) + >>> c2 = (3,2,1) + >>> L.lagdiv(c1,c2) # quotient "intuitive," remainder not + (array([ 3.]), array([-8., -4.])) + >>> c2 = (0,1,2,3) + >>> L.lagdiv(c2,c1) # neither "intuitive" + (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + + """ + # c1, c2 are trimmed copies + [c1, c2] = pu.as_series([c1, c2]) + if c2[-1] == 0 : + raise ZeroDivisionError() + + lc1 = len(c1) + lc2 = len(c2) + if lc1 < lc2 : + return c1[:1]*0, c1 + elif lc2 == 1 : + return c1/c2[-1], c1[:1]*0 + else : + quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) + rem = c1 + for i in range(lc1 - lc2, - 1, -1): + p = lagmul([0]*i + [1], c2) + q = rem[-1]/p[-1] + rem = rem[:-1] - q*p[:-1] + quo[i] = q + return quo, pu.trimseq(rem) + + +def lagpow(cs, pow, maxpower=16) : + """Raise a Laguerre series to a power. + + Returns the Laguerre series `cs` raised to the power `pow`. The + arguement `cs` is a sequence of coefficients ordered from low to high. + i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` + + Parameters + ---------- + cs : array_like + 1d array of Laguerre series coefficients ordered from low to + high. + pow : integer + Power to which the series will be raised + maxpower : integer, optional + Maximum power allowed. This is mainly to limit growth of the series + to umanageable size. Default is 16 + + Returns + ------- + coef : ndarray + Laguerre series of power. + + See Also + -------- + lagadd, lagsub, lagmul, lagdiv + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + power = int(pow) + if power != pow or power < 0 : + raise ValueError("Power must be a non-negative integer.") + elif maxpower is not None and power > maxpower : + raise ValueError("Power is too large") + elif power == 0 : + return np.array([1], dtype=cs.dtype) + elif power == 1 : + return cs + else : + # This can be made more efficient by using powers of two + # in the usual way. + prd = cs + for i in range(2, power + 1) : + prd = lagmul(prd, cs) + return prd + + +def lagder(cs, m=1, scl=1) : + """ + Differentiate a Laguerre series. + + Returns the series `cs` differentiated `m` times. At each iteration the + result is multiplied by `scl` (the scaling factor is for use in a linear + change of variable). The argument `cs` is the sequence of coefficients + from lowest order "term" to highest, e.g., [1,2,3] represents the series + ``P_0 + 2*P_1 + 3*P_2``. + + Parameters + ---------- + cs: array_like + 1-d array of Laguerre series coefficients ordered from low to high. + m : int, optional + Number of derivatives taken, must be non-negative. (Default: 1) + scl : scalar, optional + Each differentiation is multiplied by `scl`. The end result is + multiplication by ``scl**m``. This is for use in a linear change of + variable. (Default: 1) + + Returns + ------- + der : ndarray + Laguerre series of the derivative. + + See Also + -------- + lagint + + Notes + ----- + In general, the result of differentiating a Laguerre series does not + resemble the same operation on a power series. Thus the result of this + function may be "un-intuitive," albeit correct; see Examples section + below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3,4) + >>> L.lagder(cs) + array([ 6., 9., 20.]) + >>> L.lagder(cs,3) + array([ 60.]) + >>> L.lagder(cs,scl=-1) + array([ -6., -9., -20.]) + >>> L.lagder(cs,2,-1) + array([ 9., 60.]) + + """ + cnt = int(m) + + if cnt != m: + raise ValueError, "The order of derivation must be integer" + if cnt < 0 : + raise ValueError, "The order of derivation must be non-negative" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + elif cnt >= len(cs): + return cs[:1]*0 + else : + for i in range(cnt): + n = len(cs) - 1 + cs *= scl + der = np.empty(n, dtype=cs.dtype) + for j in range(n, 0, -1): + der[j - 1] = -cs[j] + cs[j - 1] += cs[j] + cs = der + return cs + + +def lagint(cs, m=1, k=[], lbnd=0, scl=1): + """ + Integrate a Laguerre series. + + Returns a Laguerre series that is the Laguerre series `cs`, integrated + `m` times from `lbnd` to `x`. At each iteration the resulting series + is **multiplied** by `scl` and an integration constant, `k`, is added. + The scaling factor is for use in a linear change of variable. ("Buyer + beware": note that, depending on what one is doing, one may want `scl` + to be the reciprocal of what one might expect; for more information, + see the Notes section below.) The argument `cs` is a sequence of + coefficients, from lowest order Laguerre series "term" to highest, + e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`. + + Parameters + ---------- + cs : array_like + 1-d array of Laguerre series coefficients, ordered from low to high. + m : int, optional + Order of integration, must be positive. (Default: 1) + k : {[], list, scalar}, optional + Integration constant(s). The value of the first integral at + ``lbnd`` is the first value in the list, the value of the second + integral at ``lbnd`` is the second value, etc. If ``k == []`` (the + default), all constants are set to zero. If ``m == 1``, a single + scalar can be given instead of a list. + lbnd : scalar, optional + The lower bound of the integral. (Default: 0) + scl : scalar, optional + Following each integration the result is *multiplied* by `scl` + before the integration constant is added. (Default: 1) + + Returns + ------- + S : ndarray + Laguerre series coefficients of the integral. + + Raises + ------ + ValueError + If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or + ``np.isscalar(scl) == False``. + + See Also + -------- + lagder + + Notes + ----- + Note that the result of each integration is *multiplied* by `scl`. + Why is this important to note? Say one is making a linear change of + variable :math:`u = ax + b` in an integral relative to `x`. Then + :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` + - perhaps not what one would have first thought. + + Also note that, in general, the result of integrating a C-series needs + to be "re-projected" onto the C-series basis set. Thus, typically, + the result of this function is "un-intuitive," albeit correct; see + Examples section below. + + Examples + -------- + >>> from numpy.polynomial import legendre as L + >>> cs = (1,2,3) + >>> L.lagint(cs) + array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.lagint(cs,3) + array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, + -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) + >>> L.lagint(cs, k=3) + array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.lagint(cs, lbnd=-2) + array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) + >>> L.lagint(cs, scl=2) + array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + + """ + cnt = int(m) + if np.isscalar(k) : + k = [k] + + if cnt != m: + raise ValueError, "The order of integration must be integer" + if cnt < 0 : + raise ValueError, "The order of integration must be non-negative" + if len(k) > cnt : + raise ValueError, "Too many integration constants" + + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if cnt == 0: + return cs + + k = list(k) + [0]*(cnt - len(k)) + for i in range(cnt) : + n = len(cs) + cs *= scl + if n == 1 and cs[0] == 0: + cs[0] += k[i] + else: + tmp = np.empty(n + 1, dtype=cs.dtype) + tmp[0] = cs[0] + tmp[1] = -cs[0] + for j in range(1, n): + tmp[j] += cs[j] + tmp[j + 1] = -cs[j] + tmp[0] += k[i] - lagval(lbnd, tmp) + cs = tmp + return cs + + +def lagval(x, cs): + """Evaluate a Laguerre series. + + If `cs` is of length `n`, this function returns : + + ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)`` + + If x is a sequence or array then p(x) will have the same shape as x. + If r is a ring_like object that supports multiplication and addition + by the values in `cs`, then an object of the same type is returned. + + Parameters + ---------- + x : array_like, ring_like + Array of numbers or objects that support multiplication and + addition with themselves and with the elements of `cs`. + cs : array_like + 1-d array of Laguerre coefficients ordered from low to high. + + Returns + ------- + values : ndarray, ring_like + If the return is an ndarray then it has the same shape as `x`. + + See Also + -------- + lagfit + + Examples + -------- + + Notes + ----- + The evaluation uses Clenshaw recursion, aka synthetic division. + + Examples + -------- + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if isinstance(x, tuple) or isinstance(x, list) : + x = np.asarray(x) + + if len(cs) == 1 : + c0 = cs[0] + c1 = 0 + elif len(cs) == 2 : + c0 = cs[0] + c1 = cs[1] + else : + nd = len(cs) + c0 = cs[-2] + c1 = cs[-1] + for i in range(3, len(cs) + 1) : + tmp = c0 + nd = nd - 1 + c0 = cs[-i] - (c1*(nd - 1))/nd + c1 = tmp + (c1*((2*nd - 1) - x))/nd + return c0 + c1*(1 - x) + + +def lagvander(x, deg) : + """Vandermonde matrix of given degree. + + Returns the Vandermonde matrix of degree `deg` and sample points `x`. + This isn't a true Vandermonde matrix because `x` can be an arbitrary + ndarray and the Laguerre polynomials aren't powers. If ``V`` is the + returned matrix and `x` is a 2d array, then the elements of ``V`` are + ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Laguerre polynomial + of degree ``k``. + + Parameters + ---------- + x : array_like + Array of points. The values are converted to double or complex + doubles. If x is scalar it is converted to a 1D array. + deg : integer + Degree of the resulting matrix. + + Returns + ------- + vander : Vandermonde matrix. + The shape of the returned matrix is ``x.shape + (deg+1,)``. The last + index is the degree. + + """ + ideg = int(deg) + if ideg != deg: + raise ValueError("deg must be integer") + if ideg < 0: + raise ValueError("deg must be non-negative") + + x = np.array(x, copy=0, ndmin=1) + 0.0 + v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype) + v[0] = x*0 + 1 + if ideg > 0 : + v[1] = 1 - x + for i in range(2, ideg + 1) : + v[i] = (v[i-1]*(2*i - 1 - x) - v[i-2]*(i - 1))/i + return np.rollaxis(v, 0, v.ndim) + + +def lagfit(x, y, deg, rcond=None, full=False, w=None): + """ + Least squares fit of Laguerre series to data. + + Fit a Laguerre series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] * + P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of + coefficients `p` that minimises the squared error. + + Parameters + ---------- + x : array_like, shape (M,) + x-coordinates of the M sample points ``(x[i], y[i])``. + y : array_like, shape (M,) or (M, K) + y-coordinates of the sample points. Several data sets of sample + points sharing the same x-coordinates can be fitted at once by + passing in a 2D-array that contains one dataset per column. + deg : int + Degree of the fitting polynomial + rcond : float, optional + Relative condition number of the fit. Singular values smaller than + this relative to the largest singular value will be ignored. The + default value is len(x)*eps, where eps is the relative precision of + the float type, about 2e-16 in most cases. + full : bool, optional + Switch determining nature of return value. When it is False (the + default) just the coefficients are returned, when True diagnostic + information from the singular value decomposition is also returned. + w : array_like, shape (`M`,), optional + Weights. If not None, the contribution of each point + ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the + weights are chosen so that the errors of the products ``w[i]*y[i]`` + all have the same variance. The default value is None. + + Returns + ------- + coef : ndarray, shape (M,) or (M, K) + Laguerre coefficients ordered from low to high. If `y` was 2-D, + the coefficients for the data in column k of `y` are in column + `k`. + + [residuals, rank, singular_values, rcond] : present when `full` = True + Residuals of the least-squares fit, the effective rank of the + scaled Vandermonde matrix and its singular values, and the + specified value of `rcond`. For more details, see `linalg.lstsq`. + + Warns + ----- + RankWarning + The rank of the coefficient matrix in the least-squares fit is + deficient. The warning is only raised if `full` = False. The + warnings can be turned off by + + >>> import warnings + >>> warnings.simplefilter('ignore', RankWarning) + + See Also + -------- + lagval : Evaluates a Laguerre series. + lagvander : Vandermonde matrix of Laguerre series. + polyfit : least squares fit using polynomials. + chebfit : least squares fit using Chebyshev series. + linalg.lstsq : Computes a least-squares fit from the matrix. + scipy.interpolate.UnivariateSpline : Computes spline fits. + + Notes + ----- + The solution are the coefficients ``c[i]`` of the Laguerre series + ``P(x)`` that minimizes the squared error + + ``E = \\sum_j |y_j - P(x_j)|^2``. + + This problem is solved by setting up as the overdetermined matrix + equation + + ``V(x)*c = y``, + + where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are + the coefficients to be solved for, and the elements of `y` are the + observed values. This equation is then solved using the singular value + decomposition of ``V``. + + If some of the singular values of ``V`` are so small that they are + neglected, then a `RankWarning` will be issued. This means that the + coeficient values may be poorly determined. Using a lower order fit + will usually get rid of the warning. The `rcond` parameter can also be + set to a value smaller than its default, but the resulting fit may be + spurious and have large contributions from roundoff error. + + Fits using Laguerre series are usually better conditioned than fits + using power series, but much can depend on the distribution of the + sample points and the smoothness of the data. If the quality of the fit + is inadequate splines may be a good alternative. + + References + ---------- + .. [1] Wikipedia, "Curve fitting", + http://en.wikipedia.org/wiki/Curve_fitting + + Examples + -------- + + """ + order = int(deg) + 1 + x = np.asarray(x) + 0.0 + y = np.asarray(y) + 0.0 + + # check arguments. + if deg < 0 : + raise ValueError, "expected deg >= 0" + if x.ndim != 1: + raise TypeError, "expected 1D vector for x" + if x.size == 0: + raise TypeError, "expected non-empty vector for x" + if y.ndim < 1 or y.ndim > 2 : + raise TypeError, "expected 1D or 2D array for y" + if len(x) != len(y): + raise TypeError, "expected x and y to have same length" + + # set up the least squares matrices + lhs = lagvander(x, deg) + rhs = y + if w is not None: + w = np.asarray(w) + 0.0 + if w.ndim != 1: + raise TypeError, "expected 1D vector for w" + if len(x) != len(w): + raise TypeError, "expected x and w to have same length" + # apply weights + if rhs.ndim == 2: + lhs *= w[:, np.newaxis] + rhs *= w[:, np.newaxis] + else: + lhs *= w[:, np.newaxis] + rhs *= w + + # set rcond + if rcond is None : + rcond = len(x)*np.finfo(x.dtype).eps + + # scale the design matrix and solve the least squares equation + scl = np.sqrt((lhs*lhs).sum(0)) + c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond) + c = (c.T/scl).T + + # warn on rank reduction + if rank != order and not full: + msg = "The fit may be poorly conditioned" + warnings.warn(msg, pu.RankWarning) + + if full : + return c, [resids, rank, s, rcond] + else : + return c + + +def lagroots(cs): + """ + Compute the roots of a Laguerre series. + + Return the roots (a.k.a "zeros") of the Laguerre series represented by + `cs`, which is the sequence of coefficients from lowest order "term" + to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + + Parameters + ---------- + cs : array_like + 1-d array of Laguerre series coefficients ordered from low to high. + + Returns + ------- + out : ndarray + Array of the roots. If all the roots are real, then so is the + dtype of ``out``; otherwise, ``out``'s dtype is complex. + + See Also + -------- + polyroots + chebroots + + Notes + ----- + Algorithm(s) used: + + Remember: because the Laguerre series basis set is different from the + "standard" basis set, the results of this function *may* not be what + one is expecting. + + Examples + -------- + >>> import numpy.polynomial as P + >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots + array([-0.60582959+0.j , -0.07208521-0.63832674j, + -0.07208521+0.63832674j]) + >>> P.lagroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots + array([-0.85099543, -0.11407192, 0.51506735]) + + """ + # cs is a trimmed copy + [cs] = pu.as_series([cs]) + if len(cs) <= 1 : + return np.array([], dtype=cs.dtype) + if len(cs) == 2 : + return np.array([1 + cs[0]/cs[1]]) + + n = len(cs) - 1 + cs /= cs[-1] + cmat = np.zeros((n,n), dtype=cs.dtype) + cmat[0, 0] = 1 + cmat[1, 0] = -1 + for i in range(1, n): + cmat[i - 1, i] = -i + cmat[i, i] = 2*i + 1 + if i != n - 1: + cmat[i + 1, i] = -(i + 1) + else: + cmat[:, i] += cs[:-1]*(i + 1) + roots = la.eigvals(cmat) + roots.sort() + return roots + + +# +# Laguerre series class +# + +exec polytemplate.substitute(name='Laguerre', nick='lag', domain='[-1,1]') diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py new file mode 100644 index 000000000..7e361e804 --- /dev/null +++ b/numpy/polynomial/tests/test_hermite.py @@ -0,0 +1,536 @@ +"""Tests for hermendre module. + +""" +from __future__ import division + +import numpy as np +import numpy.polynomial.hermite as herm +import numpy.polynomial.polynomial as poly +from numpy.testing import * + +H0 = np.array([ 1]) +H1 = np.array([0, 2]) +H2 = np.array([ -2, 0, 4]) +H3 = np.array([0, -12, 0, 8]) +H4 = np.array([ 12, 0, -48, 0, 16]) +H5 = np.array([0, 120, 0, -160, 0, 32]) +H6 = np.array([-120, 0, 720, 0, -480, 0, 64]) +H7 = np.array([0, -1680, 0, 3360, 0, -1344, 0, 128]) +H8 = np.array([1680, 0, -13440, 0, 13440, 0, -3584, 0, 256]) +H9 = np.array([0, 30240, 0, -80640, 0, 48384, 0, -9216, 0, 512]) + +Hlist = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9] + +def trim(x) : + return herm.hermtrim(x, tol=1e-6) + + +class TestConstants(TestCase) : + + def test_hermdomain(self) : + assert_equal(herm.hermdomain, [-1, 1]) + + def test_hermzero(self) : + assert_equal(herm.hermzero, [0]) + + def test_hermone(self) : + assert_equal(herm.hermone, [1]) + + def test_hermx(self) : + assert_equal(herm.hermx, [0, .5]) + + +class TestArithmetic(TestCase) : + x = np.linspace(-3, 3, 100) + y0 = poly.polyval(x, H0) + y1 = poly.polyval(x, H1) + y2 = poly.polyval(x, H2) + y3 = poly.polyval(x, H3) + y4 = poly.polyval(x, H4) + y5 = poly.polyval(x, H5) + y6 = poly.polyval(x, H6) + y7 = poly.polyval(x, H7) + y8 = poly.polyval(x, H8) + y9 = poly.polyval(x, H9) + y = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9] + + def test_hermval(self) : + def f(x) : + return x*(x**2 - 1) + + #check empty input + assert_equal(herm.hermval([], [1]).size, 0) + + #check normal input) + for i in range(10) : + msg = "At i=%d" % i + ser = np.zeros + tgt = self.y[i] + res = herm.hermval(self.x, [0]*i + [1]) + assert_almost_equal(res, tgt, err_msg=msg) + + #check that shape is preserved + for i in range(3) : + dims = [2]*i + x = np.zeros(dims) + assert_equal(herm.hermval(x, [1]).shape, dims) + assert_equal(herm.hermval(x, [1,0]).shape, dims) + assert_equal(herm.hermval(x, [1,0,0]).shape, dims) + + def test_hermadd(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] += 1 + res = herm.hermadd([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_hermsub(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] -= 1 + res = herm.hermsub([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_hermmulx(self): + assert_equal(herm.hermmulx([0]), [0]) + assert_equal(herm.hermmulx([1]), [0,.5]) + for i in range(1, 5): + ser = [0]*i + [1] + tgt = [0]*(i - 1) + [i, 0, .5] + assert_equal(herm.hermmulx(ser), tgt) + + def test_hermmul(self) : + # check values of result + for i in range(5) : + pol1 = [0]*i + [1] + val1 = herm.hermval(self.x, pol1) + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + pol2 = [0]*j + [1] + val2 = herm.hermval(self.x, pol2) + pol3 = herm.hermmul(pol1, pol2) + val3 = herm.hermval(self.x, pol3) + assert_(len(pol3) == i + j + 1, msg) + assert_almost_equal(val3, val1*val2, err_msg=msg) + + def test_hermdiv(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + ci = [0]*i + [1] + cj = [0]*j + [1] + tgt = herm.hermadd(ci, cj) + quo, rem = herm.hermdiv(tgt, ci) + res = herm.hermadd(herm.hermmul(quo, ci), rem) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + +class TestCalculus(TestCase) : + + def test_hermint(self) : + # check exceptions + assert_raises(ValueError, herm.hermint, [0], .5) + assert_raises(ValueError, herm.hermint, [0], -1) + assert_raises(ValueError, herm.hermint, [0], 1, [0,0]) + + # test integration of zero polynomial + for i in range(2, 5): + k = [0]*(i - 2) + [1] + res = herm.hermint([0], m=i, k=k) + assert_almost_equal(res, [0, .5]) + + # check single integration with integration constant + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [1/scl] + hermpol = herm.poly2herm(pol) + hermint = herm.hermint(hermpol, m=1, k=[i]) + res = herm.herm2poly(hermint) + assert_almost_equal(trim(res), trim(tgt)) + + # check single integration with integration constant and lbnd + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + hermpol = herm.poly2herm(pol) + hermint = herm.hermint(hermpol, m=1, k=[i], lbnd=-1) + assert_almost_equal(herm.hermval(-1, hermint), i) + + # check single integration with integration constant and scaling + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [2/scl] + hermpol = herm.poly2herm(pol) + hermint = herm.hermint(hermpol, m=1, k=[i], scl=2) + res = herm.herm2poly(hermint) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with default k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herm.hermint(tgt, m=1) + res = herm.hermint(pol, m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with defined k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herm.hermint(tgt, m=1, k=[k]) + res = herm.hermint(pol, m=j, k=range(j)) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with lbnd + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herm.hermint(tgt, m=1, k=[k], lbnd=-1) + res = herm.hermint(pol, m=j, k=range(j), lbnd=-1) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with scaling + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herm.hermint(tgt, m=1, k=[k], scl=2) + res = herm.hermint(pol, m=j, k=range(j), scl=2) + assert_almost_equal(trim(res), trim(tgt)) + + def test_hermder(self) : + # check exceptions + assert_raises(ValueError, herm.hermder, [0], .5) + assert_raises(ValueError, herm.hermder, [0], -1) + + # check that zeroth deriviative does nothing + for i in range(5) : + tgt = [1] + [0]*i + res = herm.hermder(tgt, m=0) + assert_equal(trim(res), trim(tgt)) + + # check that derivation is the inverse of integration + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = herm.hermder(herm.hermint(tgt, m=j), m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check derivation with scaling + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = herm.hermder(herm.hermint(tgt, m=j, scl=2), m=j, scl=.5) + assert_almost_equal(trim(res), trim(tgt)) + + +class TestMisc(TestCase) : + + def test_hermfromroots(self) : + res = herm.hermfromroots([]) + assert_almost_equal(trim(res), [1]) + for i in range(1,5) : + roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2]) + pol = herm.hermfromroots(roots) + res = herm.hermval(roots, pol) + tgt = 0 + assert_(len(pol) == i + 1) + assert_almost_equal(herm.herm2poly(pol)[-1], 1) + assert_almost_equal(res, tgt) + + def test_hermroots(self) : + assert_almost_equal(herm.hermroots([1]), []) + assert_almost_equal(herm.hermroots([1, 1]), [-.5]) + for i in range(2,5) : + tgt = np.linspace(-1, 1, i) + res = herm.hermroots(herm.hermfromroots(tgt)) + assert_almost_equal(trim(res), trim(tgt)) + + def test_hermvander(self) : + # check for 1d x + x = np.arange(3) + v = herm.hermvander(x, 3) + assert_(v.shape == (3,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], herm.hermval(x, coef)) + + # check for 2d x + x = np.array([[1,2],[3,4],[5,6]]) + v = herm.hermvander(x, 3) + assert_(v.shape == (3,2,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], herm.hermval(x, coef)) + + def test_hermfit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + + # Test exceptions + assert_raises(ValueError, herm.hermfit, [1], [1], -1) + assert_raises(TypeError, herm.hermfit, [[1]], [1], 0) + assert_raises(TypeError, herm.hermfit, [], [1], 0) + assert_raises(TypeError, herm.hermfit, [1], [[[1]]], 0) + assert_raises(TypeError, herm.hermfit, [1, 2], [1], 0) + assert_raises(TypeError, herm.hermfit, [1], [1, 2], 0) + assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[[1]]) + assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[1,1]) + + # Test fit + x = np.linspace(0,2) + y = f(x) + # + coef3 = herm.hermfit(x, y, 3) + assert_equal(len(coef3), 4) + assert_almost_equal(herm.hermval(x, coef3), y) + # + coef4 = herm.hermfit(x, y, 4) + assert_equal(len(coef4), 5) + assert_almost_equal(herm.hermval(x, coef4), y) + # + coef2d = herm.hermfit(x, np.array([y,y]).T, 3) + assert_almost_equal(coef2d, np.array([coef3,coef3]).T) + # test weighting + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + y[0::2] = 0 + wcoef3 = herm.hermfit(x, yw, 3, w=w) + assert_almost_equal(wcoef3, coef3) + # + wcoef2d = herm.hermfit(x, np.array([yw,yw]).T, 3, w=w) + assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T) + + def test_hermtrim(self) : + coef = [2, -1, 1, 0] + + # Test exceptions + assert_raises(ValueError, herm.hermtrim, coef, -1) + + # Test results + assert_equal(herm.hermtrim(coef), coef[:-1]) + assert_equal(herm.hermtrim(coef, 1), coef[:-3]) + assert_equal(herm.hermtrim(coef, 2), [0]) + + def test_hermline(self) : + assert_equal(herm.hermline(3,4), [3, 2]) + + def test_herm2poly(self) : + for i in range(10) : + assert_almost_equal(herm.herm2poly([0]*i + [1]), Hlist[i]) + + def test_poly2herm(self) : + for i in range(10) : + assert_almost_equal(herm.poly2herm(Hlist[i]), [0]*i + [1]) + + +def assert_poly_almost_equal(p1, p2): + assert_almost_equal(p1.coef, p2.coef) + assert_equal(p1.domain, p2.domain) + + +class TestHermiteClass(TestCase) : + + p1 = herm.Hermite([1,2,3]) + p2 = herm.Hermite([1,2,3], [0,1]) + p3 = herm.Hermite([1,2]) + p4 = herm.Hermite([2,2,3]) + p5 = herm.Hermite([3,2,3]) + + def test_equal(self) : + assert_(self.p1 == self.p1) + assert_(self.p2 == self.p2) + assert_(not self.p1 == self.p2) + assert_(not self.p1 == self.p3) + assert_(not self.p1 == [1,2,3]) + + def test_not_equal(self) : + assert_(not self.p1 != self.p1) + assert_(not self.p2 != self.p2) + assert_(self.p1 != self.p2) + assert_(self.p1 != self.p3) + assert_(self.p1 != [1,2,3]) + + def test_add(self) : + tgt = herm.Hermite([2,4,6]) + assert_(self.p1 + self.p1 == tgt) + assert_(self.p1 + [1,2,3] == tgt) + assert_([1,2,3] + self.p1 == tgt) + + def test_sub(self) : + tgt = herm.Hermite([1]) + assert_(self.p4 - self.p1 == tgt) + assert_(self.p4 - [1,2,3] == tgt) + assert_([2,2,3] - self.p1 == tgt) + + def test_mul(self) : + tgt = herm.Hermite([ 81., 52., 82., 12., 9.]) + assert_poly_almost_equal(self.p1 * self.p1, tgt) + assert_poly_almost_equal(self.p1 * [1,2,3], tgt) + assert_poly_almost_equal([1,2,3] * self.p1, tgt) + + def test_floordiv(self) : + tgt = herm.Hermite([1]) + assert_(self.p4 // self.p1 == tgt) + assert_(self.p4 // [1,2,3] == tgt) + assert_([2,2,3] // self.p1 == tgt) + + def test_mod(self) : + tgt = herm.Hermite([1]) + assert_((self.p4 % self.p1) == tgt) + assert_((self.p4 % [1,2,3]) == tgt) + assert_(([2,2,3] % self.p1) == tgt) + + def test_divmod(self) : + tquo = herm.Hermite([1]) + trem = herm.Hermite([2]) + quo, rem = divmod(self.p5, self.p1) + assert_(quo == tquo and rem == trem) + quo, rem = divmod(self.p5, [1,2,3]) + assert_(quo == tquo and rem == trem) + quo, rem = divmod([3,2,3], self.p1) + assert_(quo == tquo and rem == trem) + + def test_pow(self) : + tgt = herm.Hermite([1]) + for i in range(5) : + res = self.p1**i + assert_(res == tgt) + tgt = tgt*self.p1 + + def test_call(self) : + # domain = [-1, 1] + x = np.linspace(-1, 1) + tgt = 3*(4*x**2 - 2) + 2*(2*x) + 1 + assert_almost_equal(self.p1(x), tgt) + + # domain = [0, 1] + x = np.linspace(0, 1) + xx = 2*x - 1 + assert_almost_equal(self.p2(x), self.p1(xx)) + + def test_degree(self) : + assert_equal(self.p1.degree(), 2) + + def test_trimdeg(self) : + assert_raises(ValueError, self.p1.cutdeg, .5) + assert_raises(ValueError, self.p1.cutdeg, -1) + assert_equal(len(self.p1.cutdeg(3)), 3) + assert_equal(len(self.p1.cutdeg(2)), 3) + assert_equal(len(self.p1.cutdeg(1)), 2) + assert_equal(len(self.p1.cutdeg(0)), 1) + + def test_convert(self) : + x = np.linspace(-1,1) + p = self.p1.convert(domain=[0,1]) + assert_almost_equal(p(x), self.p1(x)) + + def test_mapparms(self) : + parms = self.p2.mapparms() + assert_almost_equal(parms, [-1, 2]) + + def test_trim(self) : + coef = [1, 1e-6, 1e-12, 0] + p = herm.Hermite(coef) + assert_equal(p.trim().coef, coef[:3]) + assert_equal(p.trim(1e-10).coef, coef[:2]) + assert_equal(p.trim(1e-5).coef, coef[:1]) + + def test_truncate(self) : + assert_raises(ValueError, self.p1.truncate, .5) + assert_raises(ValueError, self.p1.truncate, 0) + assert_equal(len(self.p1.truncate(4)), 3) + assert_equal(len(self.p1.truncate(3)), 3) + assert_equal(len(self.p1.truncate(2)), 2) + assert_equal(len(self.p1.truncate(1)), 1) + + def test_copy(self) : + p = self.p1.copy() + assert_(self.p1 == p) + + def test_integ(self) : + p = self.p2.integ() + assert_almost_equal(p.coef, herm.hermint([1,2,3], 1, 0, scl=.5)) + p = self.p2.integ(lbnd=0) + assert_almost_equal(p(0), 0) + p = self.p2.integ(1, 1) + assert_almost_equal(p.coef, herm.hermint([1,2,3], 1, 1, scl=.5)) + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.coef, herm.hermint([1,2,3], 2, [1,2], scl=.5)) + + def test_deriv(self) : + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef) + assert_almost_equal(p.deriv(2).coef, self.p2.coef) + + def test_roots(self) : + p = herm.Hermite(herm.poly2herm([0, -1, 0, 1]), [0, 1]) + res = p.roots() + tgt = [0, .5, 1] + assert_almost_equal(res, tgt) + + def test_linspace(self): + xdes = np.linspace(0, 1, 20) + ydes = self.p2(xdes) + xres, yres = self.p2.linspace(20) + assert_almost_equal(xres, xdes) + assert_almost_equal(yres, ydes) + + def test_fromroots(self) : + roots = [0, .5, 1] + p = herm.Hermite.fromroots(roots, domain=[0, 1]) + res = p.coef + tgt = herm.poly2herm([0, -1, 0, 1]) + assert_almost_equal(res, tgt) + + def test_fit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + x = np.linspace(0,3) + y = f(x) + + # test default value of domain + p = herm.Hermite.fit(x, y, 3) + assert_almost_equal(p.domain, [0,3]) + + # test that fit works in given domains + p = herm.Hermite.fit(x, y, 3, None) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [0,3]) + p = herm.Hermite.fit(x, y, 3, []) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [-1, 1]) + # test that fit accepts weights. + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + yw[0::2] = 0 + p = herm.Hermite.fit(x, yw, 3, w=w) + assert_almost_equal(p(x), y) + + def test_identity(self) : + x = np.linspace(0,3) + p = herm.Hermite.identity() + assert_almost_equal(p(x), x) + p = herm.Hermite.identity([1,3]) + assert_almost_equal(p(x), x) +# + +if __name__ == "__main__": + run_module_suite() diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py new file mode 100644 index 000000000..aa01baf8e --- /dev/null +++ b/numpy/polynomial/tests/test_hermite_e.py @@ -0,0 +1,536 @@ +"""Tests for hermeendre module. + +""" +from __future__ import division + +import numpy as np +import numpy.polynomial.hermite_e as herme +import numpy.polynomial.polynomial as poly +from numpy.testing import * + +He0 = np.array([ 1 ]) +He1 = np.array([ 0 , 1 ]) +He2 = np.array([ -1 ,0 , 1 ]) +He3 = np.array([ 0 , -3 ,0 , 1 ]) +He4 = np.array([ 3 ,0 , -6 ,0 , 1 ]) +He5 = np.array([ 0 , 15 ,0 , -10 ,0 , 1 ]) +He6 = np.array([ -15 ,0 , 45 ,0 , -15 ,0 , 1 ]) +He7 = np.array([ 0 , -105 ,0 , 105 ,0 , -21 ,0 , 1 ]) +He8 = np.array([ 105 ,0 , -420 ,0 , 210 ,0 , -28 ,0 , 1 ]) +He9 = np.array([ 0 , 945 ,0 , -1260 ,0 , 378 ,0 , -36 ,0 , 1 ]) + +Helist = [He0, He1, He2, He3, He4, He5, He6, He7, He8, He9] + +def trim(x) : + return herme.hermetrim(x, tol=1e-6) + + +class TestConstants(TestCase) : + + def test_hermedomain(self) : + assert_equal(herme.hermedomain, [-1, 1]) + + def test_hermezero(self) : + assert_equal(herme.hermezero, [0]) + + def test_hermeone(self) : + assert_equal(herme.hermeone, [1]) + + def test_hermex(self) : + assert_equal(herme.hermex, [0, 1]) + + +class TestArithmetic(TestCase) : + x = np.linspace(-3, 3, 100) + y0 = poly.polyval(x, He0) + y1 = poly.polyval(x, He1) + y2 = poly.polyval(x, He2) + y3 = poly.polyval(x, He3) + y4 = poly.polyval(x, He4) + y5 = poly.polyval(x, He5) + y6 = poly.polyval(x, He6) + y7 = poly.polyval(x, He7) + y8 = poly.polyval(x, He8) + y9 = poly.polyval(x, He9) + y = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9] + + def test_hermeval(self) : + def f(x) : + return x*(x**2 - 1) + + #check empty input + assert_equal(herme.hermeval([], [1]).size, 0) + + #check normal input) + for i in range(10) : + msg = "At i=%d" % i + ser = np.zeros + tgt = self.y[i] + res = herme.hermeval(self.x, [0]*i + [1]) + assert_almost_equal(res, tgt, err_msg=msg) + + #check that shape is preserved + for i in range(3) : + dims = [2]*i + x = np.zeros(dims) + assert_equal(herme.hermeval(x, [1]).shape, dims) + assert_equal(herme.hermeval(x, [1,0]).shape, dims) + assert_equal(herme.hermeval(x, [1,0,0]).shape, dims) + + def test_hermeadd(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] += 1 + res = herme.hermeadd([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_hermesub(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] -= 1 + res = herme.hermesub([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_hermemulx(self): + assert_equal(herme.hermemulx([0]), [0]) + assert_equal(herme.hermemulx([1]), [0,.5]) + for i in range(1, 5): + ser = [0]*i + [1] + tgt = [0]*(i - 1) + [i, 0, .5] + assert_equal(herme.hermemulx(ser), tgt) + + def test_hermemul(self) : + # check values of result + for i in range(5) : + pol1 = [0]*i + [1] + val1 = herme.hermeval(self.x, pol1) + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + pol2 = [0]*j + [1] + val2 = herme.hermeval(self.x, pol2) + pol3 = herme.hermemul(pol1, pol2) + val3 = herme.hermeval(self.x, pol3) + assert_(len(pol3) == i + j + 1, msg) + assert_almost_equal(val3, val1*val2, err_msg=msg) + + def test_hermediv(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + ci = [0]*i + [1] + cj = [0]*j + [1] + tgt = herme.hermeadd(ci, cj) + quo, rem = herme.hermediv(tgt, ci) + res = herme.hermeadd(herme.hermemul(quo, ci), rem) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + +class TestCalculus(TestCase) : + + def test_hermeint(self) : + # check exceptions + assert_raises(ValueError, herme.hermeint, [0], .5) + assert_raises(ValueError, herme.hermeint, [0], -1) + assert_raises(ValueError, herme.hermeint, [0], 1, [0,0]) + + # test integration of zero polynomial + for i in range(2, 5): + k = [0]*(i - 2) + [1] + res = herme.hermeint([0], m=i, k=k) + assert_almost_equal(res, [0, .5]) + + # check single integration with integration constant + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [1/scl] + hermepol = herme.poly2herme(pol) + hermeint = herme.hermeint(hermepol, m=1, k=[i]) + res = herme.herme2poly(hermeint) + assert_almost_equal(trim(res), trim(tgt)) + + # check single integration with integration constant and lbnd + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + hermepol = herme.poly2herme(pol) + hermeint = herme.hermeint(hermepol, m=1, k=[i], lbnd=-1) + assert_almost_equal(herme.hermeval(-1, hermeint), i) + + # check single integration with integration constant and scaling + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [2/scl] + hermepol = herme.poly2herme(pol) + hermeint = herme.hermeint(hermepol, m=1, k=[i], scl=2) + res = herme.herme2poly(hermeint) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with default k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herme.hermeint(tgt, m=1) + res = herme.hermeint(pol, m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with defined k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herme.hermeint(tgt, m=1, k=[k]) + res = herme.hermeint(pol, m=j, k=range(j)) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with lbnd + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herme.hermeint(tgt, m=1, k=[k], lbnd=-1) + res = herme.hermeint(pol, m=j, k=range(j), lbnd=-1) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with scaling + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = herme.hermeint(tgt, m=1, k=[k], scl=2) + res = herme.hermeint(pol, m=j, k=range(j), scl=2) + assert_almost_equal(trim(res), trim(tgt)) + + def test_hermeder(self) : + # check exceptions + assert_raises(ValueError, herme.hermeder, [0], .5) + assert_raises(ValueError, herme.hermeder, [0], -1) + + # check that zeroth deriviative does nothing + for i in range(5) : + tgt = [1] + [0]*i + res = herme.hermeder(tgt, m=0) + assert_equal(trim(res), trim(tgt)) + + # check that derivation is the inverse of integration + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = herme.hermeder(herme.hermeint(tgt, m=j), m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check derivation with scaling + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = herme.hermeder(herme.hermeint(tgt, m=j, scl=2), m=j, scl=.5) + assert_almost_equal(trim(res), trim(tgt)) + + +class TestMisc(TestCase) : + + def test_hermefromroots(self) : + res = herme.hermefromroots([]) + assert_almost_equal(trim(res), [1]) + for i in range(1,5) : + roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2]) + pol = herme.hermefromroots(roots) + res = herme.hermeval(roots, pol) + tgt = 0 + assert_(len(pol) == i + 1) + assert_almost_equal(herme.herme2poly(pol)[-1], 1) + assert_almost_equal(res, tgt) + + def test_hermeroots(self) : + assert_almost_equal(herme.hermeroots([1]), []) + assert_almost_equal(herme.hermeroots([1, 1]), [-.5]) + for i in range(2,5) : + tgt = np.linspace(-1, 1, i) + res = herme.hermeroots(herme.hermefromroots(tgt)) + assert_almost_equal(trim(res), trim(tgt)) + + def test_hermevander(self) : + # check for 1d x + x = np.arange(3) + v = herme.hermevander(x, 3) + assert_(v.shape == (3,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], herme.hermeval(x, coef)) + + # check for 2d x + x = np.array([[1,2],[3,4],[5,6]]) + v = herme.hermevander(x, 3) + assert_(v.shape == (3,2,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], herme.hermeval(x, coef)) + + def test_hermefit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + + # Test exceptions + assert_raises(ValueError, herme.hermefit, [1], [1], -1) + assert_raises(TypeError, herme.hermefit, [[1]], [1], 0) + assert_raises(TypeError, herme.hermefit, [], [1], 0) + assert_raises(TypeError, herme.hermefit, [1], [[[1]]], 0) + assert_raises(TypeError, herme.hermefit, [1, 2], [1], 0) + assert_raises(TypeError, herme.hermefit, [1], [1, 2], 0) + assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[[1]]) + assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[1,1]) + + # Test fit + x = np.linspace(0,2) + y = f(x) + # + coef3 = herme.hermefit(x, y, 3) + assert_equal(len(coef3), 4) + assert_almost_equal(herme.hermeval(x, coef3), y) + # + coef4 = herme.hermefit(x, y, 4) + assert_equal(len(coef4), 5) + assert_almost_equal(herme.hermeval(x, coef4), y) + # + coef2d = herme.hermefit(x, np.array([y,y]).T, 3) + assert_almost_equal(coef2d, np.array([coef3,coef3]).T) + # test weighting + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + y[0::2] = 0 + wcoef3 = herme.hermefit(x, yw, 3, w=w) + assert_almost_equal(wcoef3, coef3) + # + wcoef2d = herme.hermefit(x, np.array([yw,yw]).T, 3, w=w) + assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T) + + def test_hermetrim(self) : + coef = [2, -1, 1, 0] + + # Test exceptions + assert_raises(ValueError, herme.hermetrim, coef, -1) + + # Test results + assert_equal(herme.hermetrim(coef), coef[:-1]) + assert_equal(herme.hermetrim(coef, 1), coef[:-3]) + assert_equal(herme.hermetrim(coef, 2), [0]) + + def test_hermeline(self) : + assert_equal(herme.hermeline(3,4), [3, 4]) + + def test_herme2poly(self) : + for i in range(10) : + assert_almost_equal(herme.herme2poly([0]*i + [1]), Hlist[i]) + + def test_poly2herme(self) : + for i in range(10) : + assert_almost_equal(herme.poly2herme(Hlist[i]), [0]*i + [1]) + + +def assert_poly_almost_equal(p1, p2): + assert_almost_equal(p1.coef, p2.coef) + assert_equal(p1.domain, p2.domain) + + +class TestHermite_eClass(TestCase) : + + p1 = herme.Hermite_e([1,2,3]) + p2 = herme.Hermite_e([1,2,3], [0,1]) + p3 = herme.Hermite_e([1,2]) + p4 = herme.Hermite_e([2,2,3]) + p5 = herme.Hermite_e([3,2,3]) + + def test_equal(self) : + assert_(self.p1 == self.p1) + assert_(self.p2 == self.p2) + assert_(not self.p1 == self.p2) + assert_(not self.p1 == self.p3) + assert_(not self.p1 == [1,2,3]) + + def test_not_equal(self) : + assert_(not self.p1 != self.p1) + assert_(not self.p2 != self.p2) + assert_(self.p1 != self.p2) + assert_(self.p1 != self.p3) + assert_(self.p1 != [1,2,3]) + + def test_add(self) : + tgt = herme.Hermite_e([2,4,6]) + assert_(self.p1 + self.p1 == tgt) + assert_(self.p1 + [1,2,3] == tgt) + assert_([1,2,3] + self.p1 == tgt) + + def test_sub(self) : + tgt = herme.Hermite_e([1]) + assert_(self.p4 - self.p1 == tgt) + assert_(self.p4 - [1,2,3] == tgt) + assert_([2,2,3] - self.p1 == tgt) + + def test_mul(self) : + tgt = herme.Hermite_e([ 81., 52., 82., 12., 9.]) + assert_poly_almost_equal(self.p1 * self.p1, tgt) + assert_poly_almost_equal(self.p1 * [1,2,3], tgt) + assert_poly_almost_equal([1,2,3] * self.p1, tgt) + + def test_floordiv(self) : + tgt = herme.Hermite_e([1]) + assert_(self.p4 // self.p1 == tgt) + assert_(self.p4 // [1,2,3] == tgt) + assert_([2,2,3] // self.p1 == tgt) + + def test_mod(self) : + tgt = herme.Hermite_e([1]) + assert_((self.p4 % self.p1) == tgt) + assert_((self.p4 % [1,2,3]) == tgt) + assert_(([2,2,3] % self.p1) == tgt) + + def test_divmod(self) : + tquo = herme.Hermite_e([1]) + trem = herme.Hermite_e([2]) + quo, rem = divmod(self.p5, self.p1) + assert_(quo == tquo and rem == trem) + quo, rem = divmod(self.p5, [1,2,3]) + assert_(quo == tquo and rem == trem) + quo, rem = divmod([3,2,3], self.p1) + assert_(quo == tquo and rem == trem) + + def test_pow(self) : + tgt = herme.Hermite_e([1]) + for i in range(5) : + res = self.p1**i + assert_(res == tgt) + tgt = tgt*self.p1 + + def test_call(self) : + # domain = [-1, 1] + x = np.linspace(-1, 1) + tgt = 3*(4*x**2 - 2) + 2*(2*x) + 1 + assert_almost_equal(self.p1(x), tgt) + + # domain = [0, 1] + x = np.linspace(0, 1) + xx = 2*x - 1 + assert_almost_equal(self.p2(x), self.p1(xx)) + + def test_degree(self) : + assert_equal(self.p1.degree(), 2) + + def test_trimdeg(self) : + assert_raises(ValueError, self.p1.cutdeg, .5) + assert_raises(ValueError, self.p1.cutdeg, -1) + assert_equal(len(self.p1.cutdeg(3)), 3) + assert_equal(len(self.p1.cutdeg(2)), 3) + assert_equal(len(self.p1.cutdeg(1)), 2) + assert_equal(len(self.p1.cutdeg(0)), 1) + + def test_convert(self) : + x = np.linspace(-1,1) + p = self.p1.convert(domain=[0,1]) + assert_almost_equal(p(x), self.p1(x)) + + def test_mapparms(self) : + parms = self.p2.mapparms() + assert_almost_equal(parms, [-1, 2]) + + def test_trim(self) : + coef = [1, 1e-6, 1e-12, 0] + p = herme.Hermite_e(coef) + assert_equal(p.trim().coef, coef[:3]) + assert_equal(p.trim(1e-10).coef, coef[:2]) + assert_equal(p.trim(1e-5).coef, coef[:1]) + + def test_truncate(self) : + assert_raises(ValueError, self.p1.truncate, .5) + assert_raises(ValueError, self.p1.truncate, 0) + assert_equal(len(self.p1.truncate(4)), 3) + assert_equal(len(self.p1.truncate(3)), 3) + assert_equal(len(self.p1.truncate(2)), 2) + assert_equal(len(self.p1.truncate(1)), 1) + + def test_copy(self) : + p = self.p1.copy() + assert_(self.p1 == p) + + def test_integ(self) : + p = self.p2.integ() + assert_almost_equal(p.coef, herme.hermeint([1,2,3], 1, 0, scl=.5)) + p = self.p2.integ(lbnd=0) + assert_almost_equal(p(0), 0) + p = self.p2.integ(1, 1) + assert_almost_equal(p.coef, herme.hermeint([1,2,3], 1, 1, scl=.5)) + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.coef, herme.hermeint([1,2,3], 2, [1,2], scl=.5)) + + def test_deriv(self) : + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef) + assert_almost_equal(p.deriv(2).coef, self.p2.coef) + + def test_roots(self) : + p = herme.Hermite_e(herme.poly2herme([0, -1, 0, 1]), [0, 1]) + res = p.roots() + tgt = [0, .5, 1] + assert_almost_equal(res, tgt) + + def test_linspace(self): + xdes = np.linspace(0, 1, 20) + ydes = self.p2(xdes) + xres, yres = self.p2.linspace(20) + assert_almost_equal(xres, xdes) + assert_almost_equal(yres, ydes) + + def test_fromroots(self) : + roots = [0, .5, 1] + p = herme.Hermite_e.fromroots(roots, domain=[0, 1]) + res = p.coef + tgt = herme.poly2herme([0, -1, 0, 1]) + assert_almost_equal(res, tgt) + + def test_fit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + x = np.linspace(0,3) + y = f(x) + + # test default value of domain + p = herme.Hermite_e.fit(x, y, 3) + assert_almost_equal(p.domain, [0,3]) + + # test that fit works in given domains + p = herme.Hermite_e.fit(x, y, 3, None) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [0,3]) + p = herme.Hermite_e.fit(x, y, 3, []) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [-1, 1]) + # test that fit accepts weights. + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + yw[0::2] = 0 + p = herme.Hermite_e.fit(x, yw, 3, w=w) + assert_almost_equal(p(x), y) + + def test_identity(self) : + x = np.linspace(0,3) + p = herme.Hermite_e.identity() + assert_almost_equal(p(x), x) + p = herme.Hermite_e.identity([1,3]) + assert_almost_equal(p(x), x) + + +if __name__ == "__main__": + run_module_suite() diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py new file mode 100644 index 000000000..4d1b80e4b --- /dev/null +++ b/numpy/polynomial/tests/test_laguerre.py @@ -0,0 +1,530 @@ +"""Tests for hermendre module. + +""" +from __future__ import division + +import numpy as np +import numpy.polynomial.laguerre as lag +import numpy.polynomial.polynomial as poly +from numpy.testing import * + +L0 = np.array([1 ])/1 +L1 = np.array([1 , -1 ])/1 +L2 = np.array([2 , -4 , 1 ])/2 +L3 = np.array([6 , -18 , 9 , -1 ])/6 +L4 = np.array([24 , -96 , 72 , -16 , 1 ])/24 +L5 = np.array([120 , -600 , 600 , -200 , 25 , -1 ])/120 +L6 = np.array([720 , -4320 , 5400 , -2400 , 450 , -36 , 1 ])/720 + +Llist = [L0, L1, L2, L3, L4, L5, L6] + +def trim(x) : + return lag.lagtrim(x, tol=1e-6) + + +class TestConstants(TestCase) : + + def test_lagdomain(self) : + assert_equal(lag.lagdomain, [0, 1]) + + def test_lagzero(self) : + assert_equal(lag.lagzero, [0]) + + def test_lagone(self) : + assert_equal(lag.lagone, [1]) + + def test_lagx(self) : + assert_equal(lag.lagx, [1, -1]) + + +class TestArithmetic(TestCase) : + x = np.linspace(-3, 3, 100) + y0 = poly.polyval(x, L0) + y1 = poly.polyval(x, L1) + y2 = poly.polyval(x, L2) + y3 = poly.polyval(x, L3) + y4 = poly.polyval(x, L4) + y5 = poly.polyval(x, L5) + y6 = poly.polyval(x, L6) + y = [y0, y1, y2, y3, y4, y5, y6] + + def test_lagval(self) : + def f(x) : + return x*(x**2 - 1) + + #check empty input + assert_equal(lag.lagval([], [1]).size, 0) + + #check normal input) + for i in range(7) : + msg = "At i=%d" % i + ser = np.zeros + tgt = self.y[i] + res = lag.lagval(self.x, [0]*i + [1]) + assert_almost_equal(res, tgt, err_msg=msg) + + #check that shape is preserved + for i in range(3) : + dims = [2]*i + x = np.zeros(dims) + assert_equal(lag.lagval(x, [1]).shape, dims) + assert_equal(lag.lagval(x, [1,0]).shape, dims) + assert_equal(lag.lagval(x, [1,0,0]).shape, dims) + + def test_lagadd(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] += 1 + res = lag.lagadd([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_lagsub(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + tgt = np.zeros(max(i,j) + 1) + tgt[i] += 1 + tgt[j] -= 1 + res = lag.lagsub([0]*i + [1], [0]*j + [1]) + assert_equal(trim(res), trim(tgt), err_msg=msg) + + def test_lagmulx(self): + assert_equal(lag.lagmulx([0]), [0]) + assert_equal(lag.lagmulx([1]), [1,-1]) + for i in range(1, 5): + ser = [0]*i + [1] + tgt = [0]*(i - 1) + [-i, 2*i + 1, -(i + 1)] + assert_almost_equal(lag.lagmulx(ser), tgt) + + def test_lagmul(self) : + # check values of result + for i in range(5) : + pol1 = [0]*i + [1] + val1 = lag.lagval(self.x, pol1) + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + pol2 = [0]*j + [1] + val2 = lag.lagval(self.x, pol2) + pol3 = lag.lagmul(pol1, pol2) + val3 = lag.lagval(self.x, pol3) + assert_(len(pol3) == i + j + 1, msg) + assert_almost_equal(val3, val1*val2, err_msg=msg) + + def test_lagdiv(self) : + for i in range(5) : + for j in range(5) : + msg = "At i=%d, j=%d" % (i,j) + ci = [0]*i + [1] + cj = [0]*j + [1] + tgt = lag.lagadd(ci, cj) + quo, rem = lag.lagdiv(tgt, ci) + res = lag.lagadd(lag.lagmul(quo, ci), rem) + assert_almost_equal(trim(res), trim(tgt), err_msg=msg) + + +class TestCalculus(TestCase) : + + def test_lagint(self) : + # check exceptions + assert_raises(ValueError, lag.lagint, [0], .5) + assert_raises(ValueError, lag.lagint, [0], -1) + assert_raises(ValueError, lag.lagint, [0], 1, [0,0]) + + # test integration of zero polynomial + for i in range(2, 5): + k = [0]*(i - 2) + [1] + res = lag.lagint([0], m=i, k=k) + assert_almost_equal(res, [1, -1]) + + # check single integration with integration constant + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [1/scl] + lagpol = lag.poly2lag(pol) + lagint = lag.lagint(lagpol, m=1, k=[i]) + res = lag.lag2poly(lagint) + assert_almost_equal(trim(res), trim(tgt)) + + # check single integration with integration constant and lbnd + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + lagpol = lag.poly2lag(pol) + lagint = lag.lagint(lagpol, m=1, k=[i], lbnd=-1) + assert_almost_equal(lag.lagval(-1, lagint), i) + + # check single integration with integration constant and scaling + for i in range(5) : + scl = i + 1 + pol = [0]*i + [1] + tgt = [i] + [0]*i + [2/scl] + lagpol = lag.poly2lag(pol) + lagint = lag.lagint(lagpol, m=1, k=[i], scl=2) + res = lag.lag2poly(lagint) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with default k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = lag.lagint(tgt, m=1) + res = lag.lagint(pol, m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with defined k + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = lag.lagint(tgt, m=1, k=[k]) + res = lag.lagint(pol, m=j, k=range(j)) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with lbnd + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = lag.lagint(tgt, m=1, k=[k], lbnd=-1) + res = lag.lagint(pol, m=j, k=range(j), lbnd=-1) + assert_almost_equal(trim(res), trim(tgt)) + + # check multiple integrations with scaling + for i in range(5) : + for j in range(2,5) : + pol = [0]*i + [1] + tgt = pol[:] + for k in range(j) : + tgt = lag.lagint(tgt, m=1, k=[k], scl=2) + res = lag.lagint(pol, m=j, k=range(j), scl=2) + assert_almost_equal(trim(res), trim(tgt)) + + def test_lagder(self) : + # check exceptions + assert_raises(ValueError, lag.lagder, [0], .5) + assert_raises(ValueError, lag.lagder, [0], -1) + + # check that zeroth deriviative does nothing + for i in range(5) : + tgt = [1] + [0]*i + res = lag.lagder(tgt, m=0) + assert_equal(trim(res), trim(tgt)) + + # check that derivation is the inverse of integration + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = lag.lagder(lag.lagint(tgt, m=j), m=j) + assert_almost_equal(trim(res), trim(tgt)) + + # check derivation with scaling + for i in range(5) : + for j in range(2,5) : + tgt = [1] + [0]*i + res = lag.lagder(lag.lagint(tgt, m=j, scl=2), m=j, scl=.5) + assert_almost_equal(trim(res), trim(tgt)) + + +class TestMisc(TestCase) : + + def test_lagfromroots(self) : + res = lag.lagfromroots([]) + assert_almost_equal(trim(res), [1]) + for i in range(1,5) : + roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2]) + pol = lag.lagfromroots(roots) + res = lag.lagval(roots, pol) + tgt = 0 + assert_(len(pol) == i + 1) + assert_almost_equal(lag.lag2poly(pol)[-1], 1) + assert_almost_equal(res, tgt) + + def test_lagroots(self) : + assert_almost_equal(lag.lagroots([1]), []) + assert_almost_equal(lag.lagroots([0, 1]), [1]) + for i in range(2,5) : + tgt = np.linspace(0, 3, i) + res = lag.lagroots(lag.lagfromroots(tgt)) + assert_almost_equal(trim(res), trim(tgt)) + + def test_lagvander(self) : + # check for 1d x + x = np.arange(3) + v = lag.lagvander(x, 3) + assert_(v.shape == (3,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], lag.lagval(x, coef)) + + # check for 2d x + x = np.array([[1,2],[3,4],[5,6]]) + v = lag.lagvander(x, 3) + assert_(v.shape == (3,2,4)) + for i in range(4) : + coef = [0]*i + [1] + assert_almost_equal(v[...,i], lag.lagval(x, coef)) + + def test_lagfit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + + # Test exceptions + assert_raises(ValueError, lag.lagfit, [1], [1], -1) + assert_raises(TypeError, lag.lagfit, [[1]], [1], 0) + assert_raises(TypeError, lag.lagfit, [], [1], 0) + assert_raises(TypeError, lag.lagfit, [1], [[[1]]], 0) + assert_raises(TypeError, lag.lagfit, [1, 2], [1], 0) + assert_raises(TypeError, lag.lagfit, [1], [1, 2], 0) + assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[[1]]) + assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[1,1]) + + # Test fit + x = np.linspace(0,2) + y = f(x) + # + coef3 = lag.lagfit(x, y, 3) + assert_equal(len(coef3), 4) + assert_almost_equal(lag.lagval(x, coef3), y) + # + coef4 = lag.lagfit(x, y, 4) + assert_equal(len(coef4), 5) + assert_almost_equal(lag.lagval(x, coef4), y) + # + coef2d = lag.lagfit(x, np.array([y,y]).T, 3) + assert_almost_equal(coef2d, np.array([coef3,coef3]).T) + # test weighting + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + y[0::2] = 0 + wcoef3 = lag.lagfit(x, yw, 3, w=w) + assert_almost_equal(wcoef3, coef3) + # + wcoef2d = lag.lagfit(x, np.array([yw,yw]).T, 3, w=w) + assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T) + + def test_lagtrim(self) : + coef = [2, -1, 1, 0] + + # Test exceptions + assert_raises(ValueError, lag.lagtrim, coef, -1) + + # Test results + assert_equal(lag.lagtrim(coef), coef[:-1]) + assert_equal(lag.lagtrim(coef, 1), coef[:-3]) + assert_equal(lag.lagtrim(coef, 2), [0]) + + def test_lagline(self) : + assert_equal(lag.lagline(3,4), [7, -4]) + + def test_lag2poly(self) : + for i in range(7) : + assert_almost_equal(lag.lag2poly([0]*i + [1]), Llist[i]) + + def test_poly2lag(self) : + for i in range(7) : + assert_almost_equal(lag.poly2lag(Llist[i]), [0]*i + [1]) + + +def assert_poly_almost_equal(p1, p2): + assert_almost_equal(p1.coef, p2.coef) + assert_equal(p1.domain, p2.domain) + + +class TestLaguerreClass(TestCase) : + + p1 = lag.Laguerre([1,2,3]) + p2 = lag.Laguerre([1,2,3], [0,1]) + p3 = lag.Laguerre([1,2]) + p4 = lag.Laguerre([2,2,3]) + p5 = lag.Laguerre([3,2,3]) + + def test_equal(self) : + assert_(self.p1 == self.p1) + assert_(self.p2 == self.p2) + assert_(not self.p1 == self.p2) + assert_(not self.p1 == self.p3) + assert_(not self.p1 == [1,2,3]) + + def test_not_equal(self) : + assert_(not self.p1 != self.p1) + assert_(not self.p2 != self.p2) + assert_(self.p1 != self.p2) + assert_(self.p1 != self.p3) + assert_(self.p1 != [1,2,3]) + + def test_add(self) : + tgt = lag.Laguerre([2,4,6]) + assert_(self.p1 + self.p1 == tgt) + assert_(self.p1 + [1,2,3] == tgt) + assert_([1,2,3] + self.p1 == tgt) + + def test_sub(self) : + tgt = lag.Laguerre([1]) + assert_(self.p4 - self.p1 == tgt) + assert_(self.p4 - [1,2,3] == tgt) + assert_([2,2,3] - self.p1 == tgt) + + def test_mul(self) : + tgt = lag.Laguerre([ 14., -16., 56., -72., 54.]) + assert_poly_almost_equal(self.p1 * self.p1, tgt) + assert_poly_almost_equal(self.p1 * [1,2,3], tgt) + assert_poly_almost_equal([1,2,3] * self.p1, tgt) + + def test_floordiv(self) : + tgt = lag.Laguerre([1]) + assert_(self.p4 // self.p1 == tgt) + assert_(self.p4 // [1,2,3] == tgt) + assert_([2,2,3] // self.p1 == tgt) + + def test_mod(self) : + tgt = lag.Laguerre([1]) + assert_((self.p4 % self.p1) == tgt) + assert_((self.p4 % [1,2,3]) == tgt) + assert_(([2,2,3] % self.p1) == tgt) + + def test_divmod(self) : + tquo = lag.Laguerre([1]) + trem = lag.Laguerre([2]) + quo, rem = divmod(self.p5, self.p1) + assert_(quo == tquo and rem == trem) + quo, rem = divmod(self.p5, [1,2,3]) + assert_(quo == tquo and rem == trem) + quo, rem = divmod([3,2,3], self.p1) + assert_(quo == tquo and rem == trem) + + def test_pow(self) : + tgt = lag.Laguerre([1]) + for i in range(5) : + res = self.p1**i + assert_(res == tgt) + tgt = tgt*self.p1 + + def test_call(self) : + # domain = [0, 1] + x = np.linspace(0, 1) + tgt = 3*(.5*x**2 - 2*x + 1) + 2*(-x + 1) + 1 + assert_almost_equal(self.p1(x), tgt) + + # domain = [0, 1] + x = np.linspace(.5, 1) + xx = 2*x - 1 + assert_almost_equal(self.p2(x), self.p1(xx)) + + def test_degree(self) : + assert_equal(self.p1.degree(), 2) + + def test_trimdeg(self) : + assert_raises(ValueError, self.p1.cutdeg, .5) + assert_raises(ValueError, self.p1.cutdeg, -1) + assert_equal(len(self.p1.cutdeg(3)), 3) + assert_equal(len(self.p1.cutdeg(2)), 3) + assert_equal(len(self.p1.cutdeg(1)), 2) + assert_equal(len(self.p1.cutdeg(0)), 1) + + def test_convert(self) : + x = np.linspace(-1,1) + p = self.p1.convert(domain=[0,1]) + assert_almost_equal(p(x), self.p1(x)) + + def test_mapparms(self) : + parms = self.p2.mapparms() + assert_almost_equal(parms, [-1, 2]) + + def test_trim(self) : + coef = [1, 1e-6, 1e-12, 0] + p = lag.Laguerre(coef) + assert_equal(p.trim().coef, coef[:3]) + assert_equal(p.trim(1e-10).coef, coef[:2]) + assert_equal(p.trim(1e-5).coef, coef[:1]) + + def test_truncate(self) : + assert_raises(ValueError, self.p1.truncate, .5) + assert_raises(ValueError, self.p1.truncate, 0) + assert_equal(len(self.p1.truncate(4)), 3) + assert_equal(len(self.p1.truncate(3)), 3) + assert_equal(len(self.p1.truncate(2)), 2) + assert_equal(len(self.p1.truncate(1)), 1) + + def test_copy(self) : + p = self.p1.copy() + assert_(self.p1 == p) + + def test_integ(self) : + p = self.p2.integ() + assert_almost_equal(p.coef, lag.lagint([1,2,3], 1, 0, scl=.5)) + p = self.p2.integ(lbnd=0) + assert_almost_equal(p(0), 0) + p = self.p2.integ(1, 1) + assert_almost_equal(p.coef, lag.lagint([1,2,3], 1, 1, scl=.5)) + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.coef, lag.lagint([1,2,3], 2, [1,2], scl=.5)) + + def test_deriv(self) : + p = self.p2.integ(2, [1, 2]) + assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef) + assert_almost_equal(p.deriv(2).coef, self.p2.coef) + + def test_roots(self) : + p = lag.Laguerre(lag.poly2lag([0, -1, 0, 1]), [0, 1]) + res = p.roots() + tgt = [0, .5, 1] + assert_almost_equal(res, tgt) + + def test_linspace(self): + xdes = np.linspace(0, 1, 20) + ydes = self.p2(xdes) + xres, yres = self.p2.linspace(20) + assert_almost_equal(xres, xdes) + assert_almost_equal(yres, ydes) + + def test_fromroots(self) : + roots = [0, .5, 1] + p = lag.Laguerre.fromroots(roots, domain=[0, 1]) + res = p.coef + tgt = lag.poly2lag([0, -1, 0, 1]) + assert_almost_equal(res, tgt) + + def test_fit(self) : + def f(x) : + return x*(x - 1)*(x - 2) + x = np.linspace(0,3) + y = f(x) + + # test default value of domain + p = lag.Laguerre.fit(x, y, 3) + assert_almost_equal(p.domain, [0,3]) + + # test that fit works in given domains + p = lag.Laguerre.fit(x, y, 3, None) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [0,3]) + p = lag.Laguerre.fit(x, y, 3, []) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, [-1, 1]) + # test that fit accepts weights. + w = np.zeros_like(x) + yw = y.copy() + w[1::2] = 1 + yw[0::2] = 0 + p = lag.Laguerre.fit(x, yw, 3, w=w) + assert_almost_equal(p(x), y) + + def test_identity(self) : + x = np.linspace(0,3) + p = lag.Laguerre.identity() + assert_almost_equal(p(x), x) + p = lag.Laguerre.identity([1,3]) + assert_almost_equal(p(x), x) +# + +if __name__ == "__main__": + run_module_suite() -- cgit v1.2.1 From 16041b5c012d60fe80da173059532ffb2d79c530 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Thu, 3 Mar 2011 21:13:13 -0700 Subject: ENH: A window attribute in polytemplate. This is helpful in defining the mappings for the Hermite and Laguerre polynomials where the domains have infinite bounds. The window allows one to specify the interval that the domain maps to instead of using the default domain as was done before. --- numpy/polynomial/polytemplate.py | 379 ++++++++++++++++++++++++++------------- 1 file changed, 250 insertions(+), 129 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/polytemplate.py b/numpy/polynomial/polytemplate.py index 2106ad84e..657b48508 100644 --- a/numpy/polynomial/polytemplate.py +++ b/numpy/polynomial/polytemplate.py @@ -31,16 +31,22 @@ class $name(pu.PolyBase) : $name coefficients, in increasing order. For example, ``(1, 2, 3)`` implies ``P_0 + 2P_1 + 3P_2`` where the ``P_i`` are a graded polynomial basis. - domain : (2,) array_like + domain : (2,) array_like, optional Domain to use. The interval ``[domain[0], domain[1]]`` is mapped to - the interval ``$domain`` by shifting and scaling. + the interval ``[window[0], window[1]]`` by shifting and scaling. + The default value is $domain. + window : (2,) array_like, optional + Window, see ``domain`` for its use. The default value is $domain. + .. versionadded:: 1.6.0 Attributes ---------- coef : (N,) array $name coefficients, from low to high. - domain : (2,) array_like - Domain that is mapped to ``$domain``. + domain : (2,) array + Domain that is mapped to ``window``. + window : (2,) array + Window that ``domain`` is mapped to. Class Attributes ---------------- @@ -49,6 +55,8 @@ class $name(pu.PolyBase) : ``p(x)**n`` is allowed. This is to limit runaway polynomial size. domain : (2,) ndarray Default domain of the class. + window : (2,) ndarray + Default window of the class. Notes ----- @@ -65,25 +73,115 @@ class $name(pu.PolyBase) : maxpower = 16 # Default domain domain = np.array($domain) + # Default window + window = np.array($domain) # Don't let participate in array operations. Value doesn't matter. __array_priority__ = 0 - def __init__(self, coef, domain=$domain) : - [coef, domain] = pu.as_series([coef, domain], trim=False) - if len(domain) != 2 : + def has_samecoef(self, other): + """Check if coefficients match. + + Parameters + ---------- + other : class instance + The other class must have the ``coef`` attribute. + + Returns + ------- + bool : boolean + True if the coefficients are the same, False otherwise. + + Notes + ----- + .. versionadded:: 1.6.0 + + """ + if len(self.coef) != len(other.coef): + return False + elif not np.all(self.coef == other.coef): + return False + else: + return True + + def has_samedomain(self, other): + """Check if domains match. + + Parameters + ---------- + other : class instance + The other class must have the ``domain`` attribute. + + Returns + ------- + bool : boolean + True if the domains are the same, False otherwise. + + Notes + ----- + .. versionadded:: 1.6.0 + + """ + return np.all(self.domain == other.domain) + + def has_samewindow(self, other): + """Check if windows match. + + Parameters + ---------- + other : class instance + The other class must have the ``window`` attribute. + + Returns + ------- + bool : boolean + True if the windows are the same, False otherwise. + + Notes + ----- + .. versionadded:: 1.6.0 + + """ + return np.all(self.window == other.window) + + def has_samewindow(self, other): + """Check if windows match. + + Parameters + ---------- + other : class instance + The other class must have the ``window`` attribute. + + Returns + ------- + bool : boolean + True if the windows are the same, False otherwise. + + """ + return np.all(self.window == other.window) + + def __init__(self, coef, domain=$domain, window=$domain) : + [coef, dom, win] = pu.as_series([coef, domain, window], trim=False) + if len(dom) != 2 : raise ValueError("Domain has wrong number of elements.") + if len(win) != 2 : + raise ValueError("Window has wrong number of elements.") self.coef = coef - self.domain = domain + self.domain = dom + self.window = win def __repr__(self): - format = "%s(%s, %s)" + format = "%s(%s, %s, %s)" coef = repr(self.coef)[6:-1] domain = repr(self.domain)[6:-1] - return format % ('$name', coef, domain) + window = repr(self.domain)[6:-1] + return format % ('$name', coef, domain, window) def __str__(self) : - format = "%s(%s, %s)" - return format % ('$nick', str(self.coef), str(self.domain)) + format = "%s(%s, %s, %s)" + coef = str(self.coef)[6:-1] + domain = str(self.domain)[6:-1] + window = str(self.domain)[6:-1] + return format % ('$nick', coef, domain, window) # Pickle and copy @@ -91,6 +189,7 @@ class $name(pu.PolyBase) : ret = self.__dict__.copy() ret['coef'] = self.coef.copy() ret['domain'] = self.domain.copy() + ret['window'] = self.window.copy() return ret def __setstate__(self, dict) : @@ -99,11 +198,10 @@ class $name(pu.PolyBase) : # Call def __call__(self, arg) : - off, scl = pu.mapparms(self.domain, $domain) + off, scl = pu.mapparms(self.domain, self.window) arg = off + scl*arg return ${nick}val(arg, self.coef) - def __iter__(self) : return iter(self.coef) @@ -112,9 +210,8 @@ class $name(pu.PolyBase) : # Numeric properties. - def __neg__(self) : - return self.__class__(-self.coef, self.domain) + return self.__class__(-self.coef, self.domain, self.window) def __pos__(self) : return self @@ -122,7 +219,7 @@ class $name(pu.PolyBase) : def __add__(self, other) : """Returns sum""" if isinstance(other, self.__class__) : - if np.all(self.domain == other.domain) : + if self.has_samedomain(other) and self.has_samewindow(other): coef = ${nick}add(self.coef, other.coef) else : raise PolyDomainError() @@ -131,12 +228,12 @@ class $name(pu.PolyBase) : coef = ${nick}add(self.coef, other) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __sub__(self, other) : """Returns difference""" if isinstance(other, self.__class__) : - if np.all(self.domain == other.domain) : + if self.has_samedomain(other) and self.has_samewindow(other): coef = ${nick}sub(self.coef, other.coef) else : raise PolyDomainError() @@ -145,12 +242,12 @@ class $name(pu.PolyBase) : coef = ${nick}sub(self.coef, other) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __mul__(self, other) : """Returns product""" if isinstance(other, self.__class__) : - if np.all(self.domain == other.domain) : + if self.has_samedomain(other) and self.has_samewindow(other): coef = ${nick}mul(self.coef, other.coef) else : raise PolyDomainError() @@ -159,7 +256,7 @@ class $name(pu.PolyBase) : coef = ${nick}mul(self.coef, other) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __div__(self, other): # set to __floordiv__ /. @@ -175,10 +272,11 @@ class $name(pu.PolyBase) : else : return NotImplemented elif np.isscalar(other) : + # this might be overly restrictive coef = self.coef/other else : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __floordiv__(self, other) : """Returns the quotient.""" @@ -192,12 +290,12 @@ class $name(pu.PolyBase) : quo, rem = ${nick}div(self.coef, other) except : return NotImplemented - return self.__class__(quo, self.domain) + return self.__class__(quo, self.domain, self.window) def __mod__(self, other) : """Returns the remainder.""" if isinstance(other, self.__class__) : - if np.all(self.domain == other.domain) : + if self.has_samedomain(other) and self.has_samewindow(other): quo, rem = ${nick}div(self.coef, other.coef) else : raise PolyDomainError() @@ -206,12 +304,12 @@ class $name(pu.PolyBase) : quo, rem = ${nick}div(self.coef, other) except : return NotImplemented - return self.__class__(rem, self.domain) + return self.__class__(rem, self.domain, self.window) def __divmod__(self, other) : """Returns quo, remainder""" if isinstance(other, self.__class__) : - if np.all(self.domain == other.domain) : + if self.has_samedomain(other) and self.has_samewindow(other): quo, rem = ${nick}div(self.coef, other.coef) else : raise PolyDomainError() @@ -220,35 +318,37 @@ class $name(pu.PolyBase) : quo, rem = ${nick}div(self.coef, other) except : return NotImplemented - return self.__class__(quo, self.domain), self.__class__(rem, self.domain) + quo = self.__class__(quo, self.domain, self.window) + rem = self.__class__(rem, self.domain, self.window) + return quo, rem def __pow__(self, other) : try : coef = ${nick}pow(self.coef, other, maxpower = self.maxpower) except : raise - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __radd__(self, other) : try : coef = ${nick}add(other, self.coef) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __rsub__(self, other): try : coef = ${nick}sub(other, self.coef) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __rmul__(self, other) : try : coef = ${nick}mul(other, self.coef) except : return NotImplemented - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def __rdiv__(self, other): # set to __floordiv__ /. @@ -263,46 +363,62 @@ class $name(pu.PolyBase) : quo, rem = ${nick}div(other, self.coef[0]) except : return NotImplemented - return self.__class__(quo, self.domain) + return self.__class__(quo, self.domain, self.window) def __rfloordiv__(self, other) : try : quo, rem = ${nick}div(other, self.coef) except : return NotImplemented - return self.__class__(quo, self.domain) + return self.__class__(quo, self.domain, self.window) def __rmod__(self, other) : try : quo, rem = ${nick}div(other, self.coef) except : return NotImplemented - return self.__class__(rem, self.domain) + return self.__class__(rem, self.domain, self.window) def __rdivmod__(self, other) : try : quo, rem = ${nick}div(other, self.coef) except : return NotImplemented - return self.__class__(quo, self.domain), self.__class__(rem, self.domain) + quo = self.__class__(quo, self.domain, self.window) + rem = self.__class__(rem, self.domain, self.window) + return quo, rem # Enhance me # some augmented arithmetic operations could be added here def __eq__(self, other) : res = isinstance(other, self.__class__) \ - and len(self.coef) == len(other.coef) \ - and np.all(self.domain == other.domain) \ - and np.all(self.coef == other.coef) + and self.has_samecoef(other) \ + and self.has_samedomain(other) \ + and self.has_samewindow(other) return res def __ne__(self, other) : return not self.__eq__(other) # - # Extra numeric functions. + # Extra methods. # + def copy(self) : + """Return a copy. + + A new instance of $name is returned that has the same + coefficients and domain as the current instance. + + Returns + ------- + new_instance : $name + New instance of $name with the same coefficients and domain. + + """ + return self.__class__(self.coef, self.domain, self.window) + def degree(self) : """The degree of the series. @@ -340,68 +456,6 @@ class $name(pu.PolyBase) : """ return self.truncate(deg + 1) - def convert(self, domain=None, kind=None) : - """Convert to different class and/or domain. - - Parameters - ---------- - domain : array_like, optional - The domain of the new series type instance. If the value is None, - then the default domain of `kind` is used. - kind : class, optional - The polynomial series type class to which the current instance - should be converted. If kind is None, then the class of the - current instance is used. - - Returns - ------- - new_series_instance : `kind` - The returned class can be of different type than the current - instance and/or have a different domain. - - Notes - ----- - Conversion between domains and class types can result in - numerically ill defined series. - - Examples - -------- - - """ - if kind is None : - kind = $name - if domain is None : - domain = kind.domain - return self(kind.identity(domain)) - - def mapparms(self) : - """Return the mapping parameters. - - The returned values define a linear map ``off + scl*x`` that is - applied to the input arguments before the series is evaluated. The - of the map depend on the domain; if the current domain is equal to - the default domain ``$domain`` the resulting map is the identity. - If the coeffients of the ``$name`` instance are to be used - separately, then the linear function must be substituted for the - ``x`` in the standard representation of the base polynomials. - - Returns - ------- - off, scl : floats or complex - The mapping function is defined by ``off + scl*x``. - - Notes - ----- - If the current domain is the interval ``[l_1, r_1]`` and the default - interval is ``[l_2, r_2]``, then the linear mapping function ``L`` is - defined by the equations:: - - L(l_1) = l_2 - L(r_1) = r_2 - - """ - return pu.mapparms(self.domain, $domain) - def trim(self, tol=0) : """Remove small leading coefficients @@ -422,7 +476,8 @@ class $name(pu.PolyBase) : Contains the new set of coefficients. """ - return self.__class__(pu.trimcoef(self.coef, tol), self.domain) + coef = pu.trimcoef(self.coef, tol) + return self.__class__(coef, self.domain, self.window) def truncate(self, size) : """Truncate series to length `size`. @@ -448,23 +503,75 @@ class $name(pu.PolyBase) : if isize != size or isize < 1 : raise ValueError("size must be a positive integer") if isize >= len(self.coef) : - return self.__class__(self.coef, self.domain) + coef = self.coef else : - return self.__class__(self.coef[:isize], self.domain) + coef = self.coef[:isize] + return self.__class__(coef, self.domain, self.window) - def copy(self) : - """Return a copy. + def convert(self, domain=None, kind=None, window=None) : + """Convert to different class and/or domain. - A new instance of $name is returned that has the same - coefficients and domain as the current instance. + Parameters + ---------- + domain : array_like, optional + The domain of the new series type instance. If the value is None, + then the default domain of `kind` is used. + kind : class, optional + The polynomial series type class to which the current instance + should be converted. If kind is None, then the class of the + current instance is used. Returns ------- - new_instance : $name - New instance of $name with the same coefficients and domain. + new_series_instance : `kind` + The returned class can be of different type than the current + instance and/or have a different domain. + + Notes + ----- + Conversion between domains and class types can result in + numerically ill defined series. + + Examples + -------- + + """ + if kind is None: + kind = $name + if domain is None: + domain = kind.domain + if window is None: + window = kind.window + return self(kind.identity(domain, window=window)) + + def mapparms(self) : + """Return the mapping parameters. + + The returned values define a linear map ``off + scl*x`` that is + applied to the input arguments before the series is evaluated. The + map depends on the ``domain`` and ``window``; if the current + ``domain`` is equal to the ``window`` the resulting map is the + identity. If the coeffients of the ``$name`` instance are to be + used by themselves outside this class, then the linear function + must be substituted for the ``x`` in the standard representation of + the base polynomials. + + Returns + ------- + off, scl : floats or complex + The mapping function is defined by ``off + scl*x``. + + Notes + ----- + If the current domain is the interval ``[l_1, r_1]`` and the window + is ``[l_2, r_2]``, then the linear mapping function ``L`` is + defined by the equations:: + + L(l_1) = l_2 + L(r_1) = r_2 """ - return self.__class__(self.coef, self.domain) + return pu.mapparms(self.domain, self.window) def integ(self, m=1, k=[], lbnd=None) : """Integrate. @@ -501,7 +608,7 @@ class $name(pu.PolyBase) : else : lbnd = off + scl*lbnd coef = ${nick}int(self.coef, m, k, lbnd, 1./scl) - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def deriv(self, m=1): """Differentiate. @@ -527,7 +634,7 @@ class $name(pu.PolyBase) : """ off, scl = self.mapparms() coef = ${nick}der(self.coef, m, scl) - return self.__class__(coef, self.domain) + return self.__class__(coef, self.domain, self.window) def roots(self) : """Return list of roots. @@ -543,9 +650,9 @@ class $name(pu.PolyBase) : """ roots = ${nick}roots(self.coef) - return pu.mapdomain(roots, $domain, self.domain) + return pu.mapdomain(roots, self.window, self.domain) - def linspace(self, n=100): + def linspace(self, n=100, domain=None): """Return x,y values at equally spaced points in domain. Returns x, y values at `n` equally spaced points across domain. @@ -566,14 +673,17 @@ class $name(pu.PolyBase) : .. versionadded:: 1.5.0 """ - x = np.linspace(self.domain[0], self.domain[1], n) + if domain is None: + domain = self.domain + x = np.linspace(domain[0], domain[1], n) y = self(x) return x, y @staticmethod - def fit(x, y, deg, domain=None, rcond=None, full=False, w=None) : + def fit(x, y, deg, domain=None, rcond=None, full=False, w=None, + window=$domain): """Least squares fit to data. Return a `$name` instance that is the least squares fit to the data @@ -616,6 +726,10 @@ class $name(pu.PolyBase) : ``w[i]*y[i]`` all have the same variance. The default value is None. .. versionadded:: 1.5.0 + window : {[beg, end]}, optional + Window to use for the returned $name instance. The default + value is ``$domain`` + .. versionadded:: 1.6.0 Returns ------- @@ -634,21 +748,25 @@ class $name(pu.PolyBase) : ${nick}fit : similar function """ - if domain is None : + if domain is None: domain = pu.getdomain(x) - elif domain == [] : + elif domain == []: domain = $domain - xnew = pu.mapdomain(x, domain, $domain) + + if window == []: + window = $domain + + xnew = pu.mapdomain(x, domain, window) res = ${nick}fit(xnew, y, deg, w=w, rcond=rcond, full=full) if full : [coef, status] = res - return $name(coef, domain=domain), status + return $name(coef, domain=domain, window=window), status else : coef = res - return $name(coef, domain=domain) + return $name(coef, domain=domain, window=window) @staticmethod - def fromroots(roots, domain=$domain) : + def fromroots(roots, domain=$domain, window=$domain) : """Return $name object with specified roots. See ${nick}fromroots for full documentation. @@ -660,12 +778,12 @@ class $name(pu.PolyBase) : """ if domain is None : domain = pu.getdomain(roots) - rnew = pu.mapdomain(roots, domain, $domain) + rnew = pu.mapdomain(roots, domain, window) coef = ${nick}fromroots(rnew) - return $name(coef, domain=domain) + return $name(coef, domain=domain, window=window) @staticmethod - def identity(domain=$domain) : + def identity(domain=$domain, window=$domain) : """Identity function. If ``p`` is the returned $name object, then ``p(x) == x`` for all @@ -676,13 +794,16 @@ class $name(pu.PolyBase) : domain : array_like The resulting array must be if the form ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of the domain. + window : array_like + The resulting array must be if the form ``[beg, end]``, where + ``beg`` and ``end`` are the endpoints of the window. Returns ------- identity : $name object """ - off, scl = pu.mapparms($domain, domain) + off, scl = pu.mapparms(window, domain) coef = ${nick}line(off, scl) - return $name(coef, domain) + return $name(coef, domain, window) '''.replace('REL_IMPORT', rel_import)) -- cgit v1.2.1 From bb70738f14caaa1bfd2478af283dd101d6931f17 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 19:38:54 -0600 Subject: ENH: Rename test_trimdeg to test_cutdeg to match method and add ability to run as script. --- numpy/polynomial/tests/test_chebyshev.py | 6 +++++- numpy/polynomial/tests/test_legendre.py | 6 +++++- numpy/polynomial/tests/test_polynomial.py | 6 +++++- 3 files changed, 15 insertions(+), 3 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py index 65bb877f4..21e4728bf 100644 --- a/numpy/polynomial/tests/test_chebyshev.py +++ b/numpy/polynomial/tests/test_chebyshev.py @@ -462,7 +462,7 @@ class TestChebyshevClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) @@ -564,3 +564,7 @@ class TestChebyshevClass(TestCase) : assert_almost_equal(p(x), x) p = ch.Chebyshev.identity([1,3]) assert_almost_equal(p(x), x) +# + +if __name__ == "__main__": + run_module_suite() diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py index f963bd2df..a35e6605e 100644 --- a/numpy/polynomial/tests/test_legendre.py +++ b/numpy/polynomial/tests/test_legendre.py @@ -429,7 +429,7 @@ class TestLegendreClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) @@ -531,3 +531,7 @@ class TestLegendreClass(TestCase) : assert_almost_equal(p(x), x) p = leg.Legendre.identity([1,3]) assert_almost_equal(p(x), x) +# + +if __name__ == "__main__": + run_module_suite() diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py index 5890ac13f..4b93ea118 100644 --- a/numpy/polynomial/tests/test_polynomial.py +++ b/numpy/polynomial/tests/test_polynomial.py @@ -400,7 +400,7 @@ class TestPolynomialClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) @@ -502,3 +502,7 @@ class TestPolynomialClass(TestCase) : assert_almost_equal(p(x), x) p = poly.Polynomial.identity([1,3]) assert_almost_equal(p(x), x) +# + +if __name__ == "__main__": + run_module_suite() -- cgit v1.2.1 From 8ddfe35f050a0b3a5daf0b320bec8f7f373a54a5 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 19:41:10 -0600 Subject: BUG: Fix hermemulx, rename class to HermiteE, and move __all__ after imports. --- numpy/polynomial/hermite_e.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index 36e452074..df0ff4e5e 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -50,17 +50,17 @@ See also """ from __future__ import division -__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', - 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval', - 'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots', - 'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'Hermite_e'] - import numpy as np import numpy.linalg as la import polyutils as pu import warnings from polytemplate import polytemplate +__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', + 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval', + 'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots', + 'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE'] + hermetrim = pu.trimcoef def poly2herme(pol) : @@ -430,7 +430,7 @@ def hermemulx(cs): prd = np.empty(len(cs) + 1, dtype=cs.dtype) prd[0] = cs[0]*0 - prd[1] = cs[0]/2 + prd[1] = cs[0] for i in range(1, len(cs)): prd[i + 1] = cs[i] prd[i - 1] += cs[i]*i @@ -1132,7 +1132,7 @@ def hermeroots(cs): # -# Hermite_e series class +# HermiteE series class # -exec polytemplate.substitute(name='Hermite_e', nick='herme', domain='[-1,1]') +exec polytemplate.substitute(name='HermiteE', nick='herme', domain='[-1,1]') -- cgit v1.2.1 From 6799e995afbcd790e6d746c16715384560cbd30d Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 19:43:46 -0600 Subject: BUG: Fix wrong target values. Change Hermite_e to HermiteE, rename test_trimdeg to test_cutdeg to match method name. --- numpy/polynomial/tests/test_hermite_e.py | 60 ++++++++++++++++---------------- 1 file changed, 30 insertions(+), 30 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py index aa01baf8e..6026e0e64 100644 --- a/numpy/polynomial/tests/test_hermite_e.py +++ b/numpy/polynomial/tests/test_hermite_e.py @@ -99,10 +99,10 @@ class TestArithmetic(TestCase) : def test_hermemulx(self): assert_equal(herme.hermemulx([0]), [0]) - assert_equal(herme.hermemulx([1]), [0,.5]) + assert_equal(herme.hermemulx([1]), [0,1]) for i in range(1, 5): ser = [0]*i + [1] - tgt = [0]*(i - 1) + [i, 0, .5] + tgt = [0]*(i - 1) + [i, 0, 1] assert_equal(herme.hermemulx(ser), tgt) def test_hermemul(self) : @@ -143,7 +143,7 @@ class TestCalculus(TestCase) : for i in range(2, 5): k = [0]*(i - 2) + [1] res = herme.hermeint([0], m=i, k=k) - assert_almost_equal(res, [0, .5]) + assert_almost_equal(res, [0, 1]) # check single integration with integration constant for i in range(5) : @@ -333,11 +333,11 @@ class TestMisc(TestCase) : def test_herme2poly(self) : for i in range(10) : - assert_almost_equal(herme.herme2poly([0]*i + [1]), Hlist[i]) + assert_almost_equal(herme.herme2poly([0]*i + [1]), Helist[i]) def test_poly2herme(self) : for i in range(10) : - assert_almost_equal(herme.poly2herme(Hlist[i]), [0]*i + [1]) + assert_almost_equal(herme.poly2herme(Helist[i]), [0]*i + [1]) def assert_poly_almost_equal(p1, p2): @@ -345,13 +345,13 @@ def assert_poly_almost_equal(p1, p2): assert_equal(p1.domain, p2.domain) -class TestHermite_eClass(TestCase) : +class TestHermiteEClass(TestCase) : - p1 = herme.Hermite_e([1,2,3]) - p2 = herme.Hermite_e([1,2,3], [0,1]) - p3 = herme.Hermite_e([1,2]) - p4 = herme.Hermite_e([2,2,3]) - p5 = herme.Hermite_e([3,2,3]) + p1 = herme.HermiteE([1,2,3]) + p2 = herme.HermiteE([1,2,3], [0,1]) + p3 = herme.HermiteE([1,2]) + p4 = herme.HermiteE([2,2,3]) + p5 = herme.HermiteE([3,2,3]) def test_equal(self) : assert_(self.p1 == self.p1) @@ -368,38 +368,38 @@ class TestHermite_eClass(TestCase) : assert_(self.p1 != [1,2,3]) def test_add(self) : - tgt = herme.Hermite_e([2,4,6]) + tgt = herme.HermiteE([2,4,6]) assert_(self.p1 + self.p1 == tgt) assert_(self.p1 + [1,2,3] == tgt) assert_([1,2,3] + self.p1 == tgt) def test_sub(self) : - tgt = herme.Hermite_e([1]) + tgt = herme.HermiteE([1]) assert_(self.p4 - self.p1 == tgt) assert_(self.p4 - [1,2,3] == tgt) assert_([2,2,3] - self.p1 == tgt) def test_mul(self) : - tgt = herme.Hermite_e([ 81., 52., 82., 12., 9.]) + tgt = herme.HermiteE([ 23., 28., 46., 12., 9.]) assert_poly_almost_equal(self.p1 * self.p1, tgt) assert_poly_almost_equal(self.p1 * [1,2,3], tgt) assert_poly_almost_equal([1,2,3] * self.p1, tgt) def test_floordiv(self) : - tgt = herme.Hermite_e([1]) + tgt = herme.HermiteE([1]) assert_(self.p4 // self.p1 == tgt) assert_(self.p4 // [1,2,3] == tgt) assert_([2,2,3] // self.p1 == tgt) def test_mod(self) : - tgt = herme.Hermite_e([1]) + tgt = herme.HermiteE([1]) assert_((self.p4 % self.p1) == tgt) assert_((self.p4 % [1,2,3]) == tgt) assert_(([2,2,3] % self.p1) == tgt) def test_divmod(self) : - tquo = herme.Hermite_e([1]) - trem = herme.Hermite_e([2]) + tquo = herme.HermiteE([1]) + trem = herme.HermiteE([2]) quo, rem = divmod(self.p5, self.p1) assert_(quo == tquo and rem == trem) quo, rem = divmod(self.p5, [1,2,3]) @@ -408,7 +408,7 @@ class TestHermite_eClass(TestCase) : assert_(quo == tquo and rem == trem) def test_pow(self) : - tgt = herme.Hermite_e([1]) + tgt = herme.HermiteE([1]) for i in range(5) : res = self.p1**i assert_(res == tgt) @@ -417,7 +417,7 @@ class TestHermite_eClass(TestCase) : def test_call(self) : # domain = [-1, 1] x = np.linspace(-1, 1) - tgt = 3*(4*x**2 - 2) + 2*(2*x) + 1 + tgt = 3*(x**2 - 1) + 2*(x) + 1 assert_almost_equal(self.p1(x), tgt) # domain = [0, 1] @@ -428,7 +428,7 @@ class TestHermite_eClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) @@ -447,7 +447,7 @@ class TestHermite_eClass(TestCase) : def test_trim(self) : coef = [1, 1e-6, 1e-12, 0] - p = herme.Hermite_e(coef) + p = herme.HermiteE(coef) assert_equal(p.trim().coef, coef[:3]) assert_equal(p.trim(1e-10).coef, coef[:2]) assert_equal(p.trim(1e-5).coef, coef[:1]) @@ -480,7 +480,7 @@ class TestHermite_eClass(TestCase) : assert_almost_equal(p.deriv(2).coef, self.p2.coef) def test_roots(self) : - p = herme.Hermite_e(herme.poly2herme([0, -1, 0, 1]), [0, 1]) + p = herme.HermiteE(herme.poly2herme([0, -1, 0, 1]), [0, 1]) res = p.roots() tgt = [0, .5, 1] assert_almost_equal(res, tgt) @@ -494,7 +494,7 @@ class TestHermite_eClass(TestCase) : def test_fromroots(self) : roots = [0, .5, 1] - p = herme.Hermite_e.fromroots(roots, domain=[0, 1]) + p = herme.HermiteE.fromroots(roots, domain=[0, 1]) res = p.coef tgt = herme.poly2herme([0, -1, 0, 1]) assert_almost_equal(res, tgt) @@ -506,14 +506,14 @@ class TestHermite_eClass(TestCase) : y = f(x) # test default value of domain - p = herme.Hermite_e.fit(x, y, 3) + p = herme.HermiteE.fit(x, y, 3) assert_almost_equal(p.domain, [0,3]) # test that fit works in given domains - p = herme.Hermite_e.fit(x, y, 3, None) + p = herme.HermiteE.fit(x, y, 3, None) assert_almost_equal(p(x), y) assert_almost_equal(p.domain, [0,3]) - p = herme.Hermite_e.fit(x, y, 3, []) + p = herme.HermiteE.fit(x, y, 3, []) assert_almost_equal(p(x), y) assert_almost_equal(p.domain, [-1, 1]) # test that fit accepts weights. @@ -521,14 +521,14 @@ class TestHermite_eClass(TestCase) : yw = y.copy() w[1::2] = 1 yw[0::2] = 0 - p = herme.Hermite_e.fit(x, yw, 3, w=w) + p = herme.HermiteE.fit(x, yw, 3, w=w) assert_almost_equal(p(x), y) def test_identity(self) : x = np.linspace(0,3) - p = herme.Hermite_e.identity() + p = herme.HermiteE.identity() assert_almost_equal(p(x), x) - p = herme.Hermite_e.identity([1,3]) + p = herme.HermiteE.identity([1,3]) assert_almost_equal(p(x), x) -- cgit v1.2.1 From 0a7870dd85f4dbbf8e6a3544856d8b752e02ac6d Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 19:46:02 -0600 Subject: ENH: Change test_trimdeg to test_cutdeg to match method name. --- numpy/polynomial/tests/test_hermite.py | 3 ++- numpy/polynomial/tests/test_laguerre.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py index 7e361e804..dea32f24a 100644 --- a/numpy/polynomial/tests/test_hermite.py +++ b/numpy/polynomial/tests/test_hermite.py @@ -21,6 +21,7 @@ H9 = np.array([0, 30240, 0, -80640, 0, 48384, 0, -9216, 0, 512]) Hlist = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9] + def trim(x) : return herm.hermtrim(x, tol=1e-6) @@ -428,7 +429,7 @@ class TestHermiteClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py index 4d1b80e4b..f3a4a930b 100644 --- a/numpy/polynomial/tests/test_laguerre.py +++ b/numpy/polynomial/tests/test_laguerre.py @@ -422,7 +422,7 @@ class TestLaguerreClass(TestCase) : def test_degree(self) : assert_equal(self.p1.degree(), 2) - def test_trimdeg(self) : + def test_cutdeg(self) : assert_raises(ValueError, self.p1.cutdeg, .5) assert_raises(ValueError, self.p1.cutdeg, -1) assert_equal(len(self.p1.cutdeg(3)), 3) -- cgit v1.2.1 From becc1257b87b4d4138e33f67a2bd3db872c11848 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 19:56:10 -0600 Subject: ENH: Import Hermite, HermiteE, and Laguerre into package namespace. --- numpy/polynomial/__init__.py | 3 +++ 1 file changed, 3 insertions(+) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/__init__.py b/numpy/polynomial/__init__.py index 6b004acc4..32b39aa98 100644 --- a/numpy/polynomial/__init__.py +++ b/numpy/polynomial/__init__.py @@ -17,6 +17,9 @@ from polynomial import * from chebyshev import * from legendre import * from polyutils import * +from hermite import Hermite +from hermite_e import HermiteE +from laguerre import Laguerre from numpy.testing import Tester test = Tester().test -- cgit v1.2.1 From a9f4f3ce9c4860b380294097b29af8efeeba6c97 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 13 Mar 2011 23:11:56 -0600 Subject: DOC: Add examples for hermite, hermite_e, and laguerre polynomials. --- numpy/polynomial/hermite.py | 170 +++++++++++++++++++++-------------------- numpy/polynomial/hermite_e.py | 171 ++++++++++++++++++++++-------------------- numpy/polynomial/laguerre.py | 169 +++++++++++++++++++++-------------------- 3 files changed, 263 insertions(+), 247 deletions(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 550e316de..d266a6453 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -96,13 +96,9 @@ def poly2herm(pol) : Examples -------- - >>> from numpy import polynomial as P - >>> p = P.Polynomial(np.arange(4)) - >>> p - Polynomial([ 0., 1., 2., 3.], [-1., 1.]) - >>> c = P.Hermite(P.poly2herm(p.coef)) - >>> c - Hermite([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + >>> from numpy.polynomial.hermite_e import poly2herme + >>> poly2herm(np.arange(4)) + array([ 1. , 2.75 , 0.5 , 0.375]) """ [pol] = pu.as_series([pol]) @@ -146,15 +142,9 @@ def herm2poly(cs) : Examples -------- - >>> c = P.Hermite(range(4)) - >>> c - Hermite([ 0., 1., 2., 3.], [-1., 1.]) - >>> p = c.convert(kind=P.Polynomial) - >>> p - Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) - >>> P.herm2poly(range(4)) - array([-1. , -3.5, 3. , 7.5]) - + >>> from numpy.polynomial.hermite import herm2poly + >>> herm2poly([ 1. , 2.75 , 0.5 , 0.375]) + array([ 0., 1., 2., 3.]) """ from polynomial import polyadd, polysub, polymulx @@ -217,11 +207,11 @@ def hermline(off, scl) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.hermline(3,2) - array([3, 2]) - >>> L.hermval(-3, L.hermline(3,2)) # should be -3 - -3.0 + >>> from numpy.polynomial.hermite import hermline, hermval + >>> hermval(0,hermline(3, 2)) + 3.0 + >>> hermval(1,hermline(3, 2)) + 5.0 """ if scl != 0 : @@ -274,12 +264,13 @@ def hermfromroots(roots) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.hermfromroots((-1,0,1)) # x^3 - x relative to the standard basis - array([ 0. , -0.4, 0. , 0.4]) - >>> j = complex(0,1) - >>> L.hermfromroots((-j,j)) # x^2 + 1 relative to the standard basis - array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + >>> from numpy.polynomial.hermite import hermfromroots, hermval + >>> coef = hermfromroots((-1, 0, 1)) + >>> hermval((-1, 0, 1), coef) + array([ 0., 0., 0.]) + >>> coef = hermfromroots((-1j, 1j)) + >>> hermval((-1j, 1j), coef) + array([ 0.+0.j, 0.+0.j]) """ if len(roots) == 0 : @@ -324,11 +315,9 @@ def hermadd(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermadd(c1,c2) - array([ 4., 4., 4.]) + >>> from numpy.polynomial.hermite import hermadd + >>> hermadd([1, 2, 3], [1, 2, 3, 4]) + array([ 2., 4., 6., 4.]) """ # c1, c2 are trimmed copies @@ -374,13 +363,9 @@ def hermsub(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermsub(c1,c2) - array([-2., 0., 2.]) - >>> L.hermsub(c2,c1) # -C.hermsub(c1,c2) - array([ 2., 0., -2.]) + >>> from numpy.polynomial.hermite import hermsub + >>> hermsub([1, 2, 3, 4], [1, 2, 3]) + array([ 0., 0., 0., 4.]) """ # c1, c2 are trimmed copies @@ -420,7 +405,13 @@ def hermmulx(cs): .. math:: - xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + xP_i(x) = (P_{i + 1}(x)/2 + i*P_{i - 1}(x)) + + Examples + -------- + >>> from numpy.polynomial.hermite import hermmulx + >>> hermmulx([1, 2, 3]) + array([ 2. , 6.5, 1. , 1.5]) """ # cs is a trimmed copy @@ -471,11 +462,9 @@ def hermmul(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2) - >>> P.hermmul(c1,c2) # multiplication requires "reprojection" - array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + >>> from numpy.polynomial.hermite import hermmul + >>> hermmul([1, 2, 3], [0, 1, 2]) + array([ 52., 29., 52., 7., 6.]) """ # s1, s2 are trimmed copies @@ -542,14 +531,13 @@ def hermdiv(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermdiv(c1,c2) # quotient "intuitive," remainder not - (array([ 3.]), array([-8., -4.])) - >>> c2 = (0,1,2,3) - >>> L.hermdiv(c2,c1) # neither "intuitive" - (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + >>> from numpy.polynomial.hermite import hermdiv + >>> hermdiv([ 52., 29., 52., 7., 6.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 0.])) + >>> hermdiv([ 54., 31., 52., 7., 6.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 2., 2.])) + >>> hermdiv([ 53., 30., 52., 7., 6.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 1., 1.])) """ # c1, c2 are trimmed copies @@ -603,6 +591,9 @@ def hermpow(cs, pow, maxpower=16) : Examples -------- + >>> from numpy.polynomial.hermite import hermpow + >>> hermpow([1, 2, 3], 2) + array([ 81., 52., 82., 12., 9.]) """ # cs is a trimmed copy @@ -664,16 +655,11 @@ def hermder(cs, m=1, scl=1) : Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3,4) - >>> L.hermder(cs) - array([ 6., 9., 20.]) - >>> L.hermder(cs,3) - array([ 60.]) - >>> L.hermder(cs,scl=-1) - array([ -6., -9., -20.]) - >>> L.hermder(cs,2,-1) - array([ 9., 60.]) + >>> from numpy.polynomial.hermite import hermder + >>> hermder([ 1. , 0.5, 0.5, 0.5]) + array([ 1., 2., 3.]) + >>> hermder([-0.5, 1./2., 1./8., 1./12., 1./16.], m=2) + array([ 1., 2., 3.]) """ cnt = int(m) @@ -762,19 +748,17 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3) - >>> L.hermint(cs) - array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermint(cs,3) - array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, - -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) - >>> L.hermint(cs, k=3) - array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermint(cs, lbnd=-2) - array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermint(cs, scl=2) - array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + >>> from numpy.polynomial.hermite import hermint + >>> hermint([1,2,3]) # integrate once, value 0 at 0. + array([ 1. , 0.5, 0.5, 0.5]) + >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0 + array([-0.5 , 0.5 , 0.125 , 0.08333333, 0.0625 ]) + >>> hermint([1,2,3], k=1) # integrate once, value 1 at 0. + array([ 2. , 0.5, 0.5, 0.5]) + >>> hermint([1,2,3], lbnd=-1) # integrate once, value 0 at -1 + array([-2. , 0.5, 0.5, 0.5]) + >>> hermint([1,2,3], m=2, k=[1,2], lbnd=-1) + array([ 1.66666667, -0.5 , 0.125 , 0.08333333, 0.0625 ]) """ cnt = int(m) @@ -847,6 +831,13 @@ def hermval(x, cs): Examples -------- + >>> from numpy.polynomial.hermite import hermval + >>> coef = [1,2,3] + >>> hermval(1, coef) + 11.0 + >>> hermval([[1,2],[3,4]], coef) + array([[ 11., 51.], + [ 115., 203.]]) """ # cs is a trimmed copy @@ -897,6 +888,15 @@ def hermvander(x, deg) : The shape of the returned matrix is ``x.shape + (deg+1,)``. The last index is the degree. + Examples + -------- + >>> from numpy.polynomial.hermite import hermvander + >>> x = np.array([-1, 0, 1]) + >>> hermvander(x, 3) + array([[ 1., -2., 2., 4.], + [ 1., 0., -2., -0.], + [ 1., 2., 2., -4.]]) + """ ideg = int(deg) if ideg != deg: @@ -1015,6 +1015,12 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): Examples -------- + >>> from numpy.polynomial.hermite import hermfit, hermval + >>> x = np.linspace(-10, 10) + >>> err = np.random.randn(len(x))/10 + >>> y = hermval(x, [1, 2, 3]) + err + >>> hermfit(x, y, 2) + array([ 0.97902637, 1.99849131, 3.00006 ]) """ order = int(deg) + 1 @@ -1104,12 +1110,12 @@ def hermroots(cs): Examples -------- - >>> import numpy.polynomial as P - >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots - array([-0.60582959+0.j , -0.07208521-0.63832674j, - -0.07208521+0.63832674j]) - >>> P.hermroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots - array([-0.85099543, -0.11407192, 0.51506735]) + >>> from numpy.polynomial.hermite import hermroots, hermfromroots + >>> coef = hermfromroots([-1, 0, 1]) + >>> coef + array([ 0. , 0.25 , 0. , 0.125]) + >>> hermroots(coef) + array([ -1.00000000e+00, -1.38777878e-17, 1.00000000e+00]) """ # cs is a trimmed copy diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index df0ff4e5e..ecb4bc3c3 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -96,13 +96,9 @@ def poly2herme(pol) : Examples -------- - >>> from numpy import polynomial as P - >>> p = P.Polynomial(np.arange(4)) - >>> p - Polynomial([ 0., 1., 2., 3.], [-1., 1.]) - >>> c = P.Hermite(P.poly2herme(p.coef)) - >>> c - Hermite([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + >>> from numpy.polynomial.hermite_e import poly2herme + >>> poly2herme(np.arange(4)) + array([ 2., 10., 2., 3.]) """ [pol] = pu.as_series([pol]) @@ -146,15 +142,9 @@ def herme2poly(cs) : Examples -------- - >>> c = P.Hermite(range(4)) - >>> c - Hermite([ 0., 1., 2., 3.], [-1., 1.]) - >>> p = c.convert(kind=P.Polynomial) - >>> p - Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) - >>> P.herme2poly(range(4)) - array([-1. , -3.5, 3. , 7.5]) - + >>> from numpy.polynomial.hermite_e import herme2poly + >>> herme2poly([ 2., 10., 2., 3.]) + array([ 0., 1., 2., 3.]) """ from polynomial import polyadd, polysub, polymulx @@ -216,11 +206,12 @@ def hermeline(off, scl) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.hermeline(3,2) - array([3, 2]) - >>> L.hermeval(-3, L.hermeline(3,2)) # should be -3 - -3.0 + >>> from numpy.polynomial.hermite_e import hermeline + >>> from numpy.polynomial.hermite_e import hermeline, hermeval + >>> hermeval(0,hermeline(3, 2)) + 3.0 + >>> hermeval(1,hermeline(3, 2)) + 5.0 """ if scl != 0 : @@ -273,12 +264,13 @@ def hermefromroots(roots) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.hermefromroots((-1,0,1)) # x^3 - x relative to the standard basis - array([ 0. , -0.4, 0. , 0.4]) - >>> j = complex(0,1) - >>> L.hermefromroots((-j,j)) # x^2 + 1 relative to the standard basis - array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + >>> from numpy.polynomial.hermite_e import hermefromroots, hermeval + >>> coef = hermefromroots((-1, 0, 1)) + >>> hermeval((-1, 0, 1), coef) + array([ 0., 0., 0.]) + >>> coef = hermefromroots((-1j, 1j)) + >>> hermeval((-1j, 1j), coef) + array([ 0.+0.j, 0.+0.j]) """ if len(roots) == 0 : @@ -323,11 +315,9 @@ def hermeadd(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermeadd(c1,c2) - array([ 4., 4., 4.]) + >>> from numpy.polynomial.hermite_e import hermeadd + >>> hermeadd([1, 2, 3], [1, 2, 3, 4]) + array([ 2., 4., 6., 4.]) """ # c1, c2 are trimmed copies @@ -373,13 +363,9 @@ def hermesub(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermesub(c1,c2) - array([-2., 0., 2.]) - >>> L.hermesub(c2,c1) # -C.hermesub(c1,c2) - array([ 2., 0., -2.]) + >>> from numpy.polynomial.hermite_e import hermesub + >>> hermesub([1, 2, 3, 4], [1, 2, 3]) + array([ 0., 0., 0., 4.]) """ # c1, c2 are trimmed copies @@ -419,7 +405,13 @@ def hermemulx(cs): .. math:: - xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + xP_i(x) = (P_{i + 1}(x) + iP_{i - 1}(x))) + + Examples + -------- + >>> from numpy.polynomial.hermite_e import hermemulx + >>> hermemulx([1, 2, 3]) + array([ 2., 7., 2., 3.]) """ # cs is a trimmed copy @@ -470,11 +462,9 @@ def hermemul(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2) - >>> P.hermemul(c1,c2) # multiplication requires "reprojection" - array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + >>> from numpy.polynomial.hermite_e import hermemul + >>> hermemul([1, 2, 3], [0, 1, 2]) + array([ 14., 15., 28., 7., 6.]) """ # s1, s2 are trimmed copies @@ -541,14 +531,11 @@ def hermediv(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.hermediv(c1,c2) # quotient "intuitive," remainder not - (array([ 3.]), array([-8., -4.])) - >>> c2 = (0,1,2,3) - >>> L.hermediv(c2,c1) # neither "intuitive" - (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + >>> from numpy.polynomial.hermite_e import hermediv + >>> hermediv([ 14., 15., 28., 7., 6.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 0.])) + >>> hermediv([ 15., 17., 28., 7., 6.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 1., 2.])) """ # c1, c2 are trimmed copies @@ -602,6 +589,9 @@ def hermepow(cs, pow, maxpower=16) : Examples -------- + >>> from numpy.polynomial.hermite_e import hermepow + >>> hermepow([1, 2, 3], 2) + array([ 23., 28., 46., 12., 9.]) """ # cs is a trimmed copy @@ -663,16 +653,11 @@ def hermeder(cs, m=1, scl=1) : Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3,4) - >>> L.hermeder(cs) - array([ 6., 9., 20.]) - >>> L.hermeder(cs,3) - array([ 60.]) - >>> L.hermeder(cs,scl=-1) - array([ -6., -9., -20.]) - >>> L.hermeder(cs,2,-1) - array([ 9., 60.]) + >>> from numpy.polynomial.hermite_e import hermeder + >>> hermeder([ 1., 1., 1., 1.]) + array([ 1., 2., 3.]) + >>> hermeder([-0.25, 1., 1./2., 1./3., 1./4 ], m=2) + array([ 1., 2., 3.]) """ cnt = int(m) @@ -761,19 +746,17 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3) - >>> L.hermeint(cs) - array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermeint(cs,3) - array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, - -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) - >>> L.hermeint(cs, k=3) - array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermeint(cs, lbnd=-2) - array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.hermeint(cs, scl=2) - array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + >>> from numpy.polynomial.hermite_e import hermeint + >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0. + array([ 1., 1., 1., 1.]) + >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0 + array([-0.25 , 1. , 0.5 , 0.33333333, 0.25 ]) + >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0. + array([ 2., 1., 1., 1.]) + >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1 + array([-1., 1., 1., 1.]) + >>> hermeint([1, 2, 3], m=2, k=[1,2], lbnd=-1) + array([ 1.83333333, 0. , 0.5 , 0.33333333, 0.25 ]) """ cnt = int(m) @@ -846,6 +829,13 @@ def hermeval(x, cs): Examples -------- + >>> from numpy.polynomial.hermite_e import hermeval + >>> coef = [1,2,3] + >>> hermeval(1, coef) + 3.0 + >>> hermeval([[1,2],[3,4]], coef) + array([[ 3., 14.], + [ 31., 54.]]) """ # cs is a trimmed copy @@ -895,6 +885,15 @@ def hermevander(x, deg) : The shape of the returned matrix is ``x.shape + (deg+1,)``. The last index is the degree. + Examples + -------- + >>> from numpy.polynomial.hermite_e import hermevander + >>> x = np.array([-1, 0, 1]) + >>> hermevander(x, 3) + array([[ 1., -1., 0., 2.], + [ 1., 0., -1., -0.], + [ 1., 1., 0., -2.]]) + """ ideg = int(deg) if ideg != deg: @@ -1012,6 +1011,12 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): Examples -------- + >>> from numpy.polynomial.hermite_e import hermefit, hermeval + >>> x = np.linspace(-10, 10) + >>> err = np.random.randn(len(x))/10 + >>> y = hermeval(x, [1, 2, 3]) + err + >>> hermefit(x, y, 2) + array([ 1.01690445, 1.99951418, 2.99948696]) """ order = int(deg) + 1 @@ -1020,7 +1025,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): # check arguments. if deg < 0 : - raise ValueError, "expected deg >= 0" + raise valueerror, "expected deg >= 0" if x.ndim != 1: raise TypeError, "expected 1D vector for x" if x.size == 0: @@ -1101,12 +1106,12 @@ def hermeroots(cs): Examples -------- - >>> import numpy.polynomial as P - >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots - array([-0.60582959+0.j , -0.07208521-0.63832674j, - -0.07208521+0.63832674j]) - >>> P.hermeroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots - array([-0.85099543, -0.11407192, 0.51506735]) + >>> from numpy.polynomial.hermite_e import hermeroots, hermefromroots + >>> coef = hermefromroots([-1, 0, 1]) + >>> coef + array([ 0., 2., 0., 1.]) + >>> hermeroots(coef) + array([-1., 0., 1.]) """ # cs is a trimmed copy diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py index 94c495deb..b6389bf63 100644 --- a/numpy/polynomial/laguerre.py +++ b/numpy/polynomial/laguerre.py @@ -96,13 +96,9 @@ def poly2lag(pol) : Examples -------- - >>> from numpy import polynomial as P - >>> p = P.Polynomial(np.arange(4)) - >>> p - Polynomial([ 0., 1., 2., 3.], [-1., 1.]) - >>> c = P.Laguerre(P.poly2lag(p.coef)) - >>> c - Laguerre([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + >>> from numpy.polynomial.laguerre import poly2lag + >>> poly2lag(np.arange(4)) + array([ 23., -63., 58., -18.]) """ [pol] = pu.as_series([pol]) @@ -146,15 +142,9 @@ def lag2poly(cs) : Examples -------- - >>> c = P.Laguerre(range(4)) - >>> c - Laguerre([ 0., 1., 2., 3.], [-1., 1.]) - >>> p = c.convert(kind=P.Polynomial) - >>> p - Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) - >>> P.lag2poly(range(4)) - array([-1. , -3.5, 3. , 7.5]) - + >>> from numpy.polynomial.laguerre import lag2poly + >>> lag2poly([ 23., -63., 58., -18.]) + array([ 0., 1., 2., 3.]) """ from polynomial import polyadd, polysub, polymulx @@ -214,11 +204,11 @@ def lagline(off, scl) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.lagline(3,2) - array([3, 2]) - >>> L.lagval(-3, L.lagline(3,2)) # should be -3 - -3.0 + >>> from numpy.polynomial.laguerre import lagline, lagval + >>> lagval(0,lagline(3, 2)) + 3.0 + >>> lagval(1,lagline(3, 2)) + 5.0 """ if scl != 0 : @@ -271,12 +261,13 @@ def lagfromroots(roots) : Examples -------- - >>> import numpy.polynomial.legendre as L - >>> L.lagfromroots((-1,0,1)) # x^3 - x relative to the standard basis - array([ 0. , -0.4, 0. , 0.4]) - >>> j = complex(0,1) - >>> L.lagfromroots((-j,j)) # x^2 + 1 relative to the standard basis - array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) + >>> from numpy.polynomial.laguerre import lagfromroots, lagval + >>> coef = lagfromroots((-1, 0, 1)) + >>> lagval((-1, 0, 1), coef) + array([ 0., 0., 0.]) + >>> coef = lagfromroots((-1j, 1j)) + >>> lagval((-1j, 1j), coef) + array([ 0.+0.j, 0.+0.j]) """ if len(roots) == 0 : @@ -321,11 +312,10 @@ def lagadd(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.lagadd(c1,c2) - array([ 4., 4., 4.]) + >>> from numpy.polynomial.laguerre import lagadd + >>> lagadd([1, 2, 3], [1, 2, 3, 4]) + array([ 2., 4., 6., 4.]) + """ # c1, c2 are trimmed copies @@ -371,13 +361,9 @@ def lagsub(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.lagsub(c1,c2) - array([-2., 0., 2.]) - >>> L.lagsub(c2,c1) # -C.lagsub(c1,c2) - array([ 2., 0., -2.]) + >>> from numpy.polynomial.laguerre import lagsub + >>> lagsub([1, 2, 3, 4], [1, 2, 3]) + array([ 0., 0., 0., 4.]) """ # c1, c2 are trimmed copies @@ -417,7 +403,13 @@ def lagmulx(cs): .. math:: - xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) + xP_i(x) = (-(i + 1)*P_{i + 1}(x) + (2i + 1)P_{i}(x) - iP_{i - 1}(x)) + + Examples + -------- + >>> from numpy.polynomial.laguerre import lagmulx + >>> lagmulx([1, 2, 3]) + array([ -1., -1., 11., -9.]) """ # cs is a trimmed copy @@ -469,11 +461,9 @@ def lagmul(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2) - >>> P.lagmul(c1,c2) # multiplication requires "reprojection" - array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) + >>> from numpy.polynomial.laguerre import lagmul + >>> lagmul([1, 2, 3], [0, 1, 2]) + array([ 8., -13., 38., -51., 36.]) """ # s1, s2 are trimmed copies @@ -540,14 +530,11 @@ def lagdiv(c1, c2): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> c1 = (1,2,3) - >>> c2 = (3,2,1) - >>> L.lagdiv(c1,c2) # quotient "intuitive," remainder not - (array([ 3.]), array([-8., -4.])) - >>> c2 = (0,1,2,3) - >>> L.lagdiv(c2,c1) # neither "intuitive" - (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) + >>> from numpy.polynomial.laguerre import lagdiv + >>> lagdiv([ 8., -13., 38., -51., 36.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 0.])) + >>> lagdiv([ 9., -12., 38., -51., 36.], [0, 1, 2]) + (array([ 1., 2., 3.]), array([ 1., 1.])) """ # c1, c2 are trimmed copies @@ -601,6 +588,9 @@ def lagpow(cs, pow, maxpower=16) : Examples -------- + >>> from numpy.polynomial.laguerre import lagpow + >>> lagpow([1, 2, 3], 2) + array([ 14., -16., 56., -72., 54.]) """ # cs is a trimmed copy @@ -662,16 +652,11 @@ def lagder(cs, m=1, scl=1) : Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3,4) - >>> L.lagder(cs) - array([ 6., 9., 20.]) - >>> L.lagder(cs,3) - array([ 60.]) - >>> L.lagder(cs,scl=-1) - array([ -6., -9., -20.]) - >>> L.lagder(cs,2,-1) - array([ 9., 60.]) + >>> from numpy.polynomial.laguerre import lagder + >>> lagder([ 1., 1., 1., -3.]) + array([ 1., 2., 3.]) + >>> lagder([ 1., 0., 0., -4., 3.], m=2) + array([ 1., 2., 3.]) """ cnt = int(m) @@ -761,19 +746,17 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1): Examples -------- - >>> from numpy.polynomial import legendre as L - >>> cs = (1,2,3) - >>> L.lagint(cs) - array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.lagint(cs,3) - array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, - -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) - >>> L.lagint(cs, k=3) - array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.lagint(cs, lbnd=-2) - array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) - >>> L.lagint(cs, scl=2) - array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) + >>> from numpy.polynomial.laguerre import lagint + >>> lagint([1,2,3]) + array([ 1., 1., 1., -3.]) + >>> lagint([1,2,3], m=2) + array([ 1., 0., 0., -4., 3.]) + >>> lagint([1,2,3], k=1) + array([ 2., 1., 1., -3.]) + >>> lagint([1,2,3], lbnd=-1) + array([ 11.5, 1. , 1. , -3. ]) + >>> lagint([1,2], m=2, k=[1,2], lbnd=-1) + array([ 11.16666667, -5. , -3. , 2. ]) """ cnt = int(m) @@ -847,6 +830,13 @@ def lagval(x, cs): Examples -------- + >>> from numpy.polynomial.laguerre import lagval + >>> coef = [1,2,3] + >>> lagval(1, coef) + -0.5 + >>> lagval([[1,2],[3,4]], coef) + array([[-0.5, -4. ], + [-4.5, -2. ]]) """ # cs is a trimmed copy @@ -896,6 +886,15 @@ def lagvander(x, deg) : The shape of the returned matrix is ``x.shape + (deg+1,)``. The last index is the degree. + Examples + -------- + >>> from numpy.polynomial.laguerre import lagvander + >>> x = np.array([0, 1, 2]) + >>> lagvander(x, 3) + array([[ 1. , 1. , 1. , 1. ], + [ 1. , 0. , -0.5 , -0.66666667], + [ 1. , -1. , -1. , -0.33333333]]) + """ ideg = int(deg) if ideg != deg: @@ -1013,6 +1012,12 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): Examples -------- + >>> from numpy.polynomial.laguerre import lagfit, lagval + >>> x = np.linspace(0, 10) + >>> err = np.random.randn(len(x))/10 + >>> y = lagval(x, [1, 2, 3]) + err + >>> lagfit(x, y, 2) + array([ 0.96971004, 2.00193749, 3.00288744]) """ order = int(deg) + 1 @@ -1102,12 +1107,12 @@ def lagroots(cs): Examples -------- - >>> import numpy.polynomial as P - >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots - array([-0.60582959+0.j , -0.07208521-0.63832674j, - -0.07208521+0.63832674j]) - >>> P.lagroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots - array([-0.85099543, -0.11407192, 0.51506735]) + >>> from numpy.polynomial.laguerre import lagroots, lagfromroots + >>> coef = lagfromroots([0, 1, 2]) + >>> coef + array([ 2., -8., 12., -6.]) + >>> lagroots(coef) + array([ -4.44089210e-16, 1.00000000e+00, 2.00000000e+00]) """ # cs is a trimmed copy -- cgit v1.2.1 From c2e4c9c034a3446de643fb4f8e96d21249abb7b9 Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Mon, 14 Mar 2011 10:18:36 -0600 Subject: BUG: Fix valueerror typo. --- numpy/polynomial/hermite_e.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'numpy/polynomial') diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index ecb4bc3c3..e644e345a 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -1025,7 +1025,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): # check arguments. if deg < 0 : - raise valueerror, "expected deg >= 0" + raise ValueError, "expected deg >= 0" if x.ndim != 1: raise TypeError, "expected 1D vector for x" if x.size == 0: -- cgit v1.2.1