Viewing file:  polynomial.py (19.33 KB) -rw-r--r--
Select action/file-type:
 (+) |  (+) |  (+) | Code (+) | Session (+) |  (+) | SDB (+) |  (+) |  (+) |  (+) |  (+) |  (+) |
"""Functions for dealing with polynomials.
This module provides a number of functions that are useful in dealing with
polynomials as well as a ``Polynomial`` class that encapsuletes the usual
arithmetic operations. All arrays of polynomial coefficients are assumed to
be ordered from low to high degree, thus `array([1,2,3])` will be treated
as the polynomial ``1 + 2*x + 3*x**2``
Constants
---------
- polydomain -- Polynomial default domain
- polyzero -- Polynomial that evaluates to 0.
- polyone -- Polynomial that evaluates to 1.
- polyx -- Polynomial of the identity map (x).
Arithmetic
----------
- polyadd -- add a polynomial to another.
- polysub -- subtract a polynomial from another.
- polymul -- multiply a polynomial by another
- polydiv -- divide one polynomial by another.
- polyval -- evaluate a polynomial at given points.
Calculus
--------
- polyder -- differentiate a polynomial.
- polyint -- integrate a polynomial.
Misc Functions
--------------
- polyfromroots -- create a polynomial with specified roots.
- polyroots -- find the roots of a polynomial.
- polyvander -- Vandermode like matrix for powers.
- polyfit -- least squares fit returning a polynomial.
- polytrim -- trim leading coefficients from a polynomial.
- polyline -- Polynomial of given straight line
Classes
-------
- Polynomial -- polynomial class.
"""
from __future__ import division
__all__ = ['polyzero', 'polyone', 'polyx', 'polydomain',
'polyline','polyadd', 'polysub', 'polymul', 'polydiv', 'polyval',
'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
'polytrim', 'polyroots', 'Polynomial']
import numpy as np
import numpy.linalg as la
import polyutils as pu
import warnings
from polytemplate import polytemplate
polytrim = pu.trimcoef
#
# These are constant arrays are of integer type so as to be compatible
# with the widest range of other types, such as Decimal.
#
# Polynomial default domain.
polydomain = np.array([-1,1])
# Polynomial coefficients representing zero.
polyzero = np.array([0])
# Polynomial coefficients representing one.
polyone = np.array([1])
# Polynomial coefficients representing the identity x.
polyx = np.array([0,1])
#
# Polynomial series functions
#
def polyline(off, scl) :
"""Polynomial whose graph is a straight line.
The line has the formula ``off + scl*x``
Parameters:
-----------
off, scl : scalars
The specified line is given by ``off + scl*x``.
Returns:
--------
series : 1d ndarray
The polynomial equal to ``off + scl*x``.
"""
if scl != 0 :
return np.array([off,scl])
else :
return np.array([off])
def polyfromroots(roots) :
"""Generate a polynomial with given roots.
Generate a polynomial whose roots are given by `roots`. The resulting
series is the produet `(x - roots[0])*(x - roots[1])*...`
Inputs
------
roots : array_like
1-d array containing the roots in sorted order.
Returns
-------
series : ndarray
1-d array containing the coefficients of the Chebeshev series
ordered from low to high.
See Also
--------
polyroots
"""
if len(roots) == 0 :
return np.ones(1)
else :
[roots] = pu.as_series([roots], trim=False)
prd = np.zeros(len(roots) + 1, dtype=roots.dtype)
prd[-1] = 1
for i in range(len(roots)) :
prd[-(i+2):-1] -= roots[i]*prd[-(i+1):]
return prd
def polyadd(c1, c2):
"""Add one polynomial to another.
Returns the sum of two polynomials `c1` + `c2`. The arguments are
sequences of coefficients ordered from low to high, i.e., [1,2,3] is
the polynomial ``1 + 2*x + 3*x**2"``.
Parameters
----------
c1, c2 : array_like
1d arrays of polynomial coefficients ordered from low to
high.
Returns
-------
out : ndarray
polynomial of the sum.
See Also
--------
polysub, polymul, polydiv, polypow
"""
# 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 polysub(c1, c2):
"""Subtract one polynomial from another.
Returns the difference of two polynomials `c1` - `c2`. The arguments
are sequences of coefficients ordered from low to high, i.e., [1,2,3]
is the polynomial ``1 + 2*x + 3*x**2``.
Parameters
----------
c1, c2 : array_like
1d arrays of polynomial coefficients ordered from low to
high.
Returns
-------
out : ndarray
polynomial of the difference.
See Also
--------
polyadd, polymul, polydiv, polypow
Examples
--------
"""
# 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 polymul(c1, c2):
"""Multiply one polynomial by another.
Returns the product of two polynomials `c1` * `c2`. The arguments
are sequences of coefficients ordered from low to high, i.e., [1,2,3]
is the polynomial ``1 + 2*x + 3*x**2.``
Parameters
----------
c1, c2 : array_like
1d arrays of polyyshev series coefficients ordered from low to
high.
Returns
-------
out : ndarray
polynomial of the product.
See Also
--------
polyadd, polysub, polydiv, polypow
"""
# c1, c2 are trimmed copies
[c1, c2] = pu.as_series([c1, c2])
ret = np.convolve(c1, c2)
return pu.trimseq(ret)
def polydiv(c1, c2):
"""Divide one polynomial by another.
Returns the quotient of two polynomials `c1` / `c2`. The arguments are
sequences of coefficients ordered from low to high, i.e., [1,2,3] is
the series ``1 + 2*x + 3*x**2.``
Parameters
----------
c1, c2 : array_like
1d arrays of chebyshev series coefficients ordered from low to
high.
Returns
-------
[quo, rem] : ndarray
polynomial of the quotient and remainder.
See Also
--------
polyadd, polysub, polymul, polypow
Examples
--------
"""
# c1, c2 are trimmed copies
[c1, c2] = pu.as_series([c1, c2])
if c2[-1] == 0 :
raise ZeroDivisionError()
len1 = len(c1)
len2 = len(c2)
if len2 == 1 :
return c1/c2[-1], c1[:1]*0
elif len1 < len2 :
return c1[:1]*0, c1
else :
dlen = len1 - len2
scl = c2[-1]
c2 = c2[:-1]/scl
i = dlen
j = len1 - 1
while i >= 0 :
c1[i:j] -= c2*c1[j]
i -= 1
j -= 1
return c1[j+1:]/scl, pu.trimseq(c1[:j+1])
def polypow(cs, pow, maxpower=None) :
"""Raise a polynomial to a power.
Returns the polynomial `cs` raised to the power `pow`. The argument
`cs` is a sequence of coefficients ordered from low to high. i.e.,
[1,2,3] is the series ``1 + 2*x + 3*x**2.``
Parameters
----------
cs : array_like
1d array of chebyshev 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
Chebyshev series of power.
See Also
--------
chebadd, chebsub, chebmul, chebdiv
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 = np.convolve(prd, cs)
return prd
def polyder(cs, m=1, scl=1) :
"""Differentiate a polynomial.
Returns the polynomial `cs` differentiated `m` times. The argument `cs`
is a sequence of coefficients ordered from low to high. i.e., [1,2,3]
is the series ``1 + 2*x + 3*x**2.``
Parameters
----------
cs: array_like
1d array of chebyshev series coefficients ordered from low to high.
m : int, optional
Order of differentiation, must be non-negative. (default: 1)
scl : scalar, optional
The result of each derivation 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
polynomial of the derivative.
See Also
--------
polyint
Examples
--------
"""
# cs is a trimmed copy
[cs] = pu.as_series([cs])
if m < 0 :
raise ValueError, "The order of derivation must be positive"
if not np.isscalar(scl) :
raise ValueError, "The scl parameter must be a scalar"
if m == 0 :
return cs
elif m >= len(cs) :
return cs[:1]*0
else :
n = len(cs)
d = np.arange(n)*scl
for i in range(m) :
cs[i:] *= d[:n-i]
return cs[i+1:].copy()
def polyint(cs, m=1, k=[], lbnd=0, scl=1) :
"""Integrate a polynomial.
Returns the polynomial `cs` integrated from `lbnd` to x `m` times. At
each iteration the resulting series is multiplied by `scl` and an
integration constant specified by `k` is added. The scaling factor is
for use in a linear change of variable. The argument `cs` is a sequence
of coefficients ordered from low to high. i.e., [1,2,3] is the
polynomial ``1 + 2*x + 3*x**2``.
Parameters
----------
cs : array_like
1d array of chebyshev series coefficients ordered from low to high.
m : int, optional
Order of integration, must be positeve. (default: 1)
k : {[], list, scalar}, optional
Integration constants. The value of the first integral at zero is
the first value in the list, the value of the second integral at
zero is the second value in the list, and so on. If ``[]``
(default), all constants are set 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
-------
der : ndarray
polynomial of the integral.
Raises
------
ValueError
See Also
--------
polyder
Examples
--------
"""
if np.isscalar(k) :
k = [k]
if m < 1 :
raise ValueError, "The order of integration must be positive"
if len(k) > m :
raise ValueError, "Too many integration constants"
if not np.isscalar(lbnd) :
raise ValueError, "The lbnd parameter must be a scalar"
if not np.isscalar(scl) :
raise ValueError, "The scl parameter must be a scalar"
# cs is a trimmed copy
[cs] = pu.as_series([cs])
k = list(k) + [0]*(m - len(k))
fac = np.arange(1, len(cs) + m)/scl
ret = np.zeros(len(cs) + m, dtype=cs.dtype)
ret[m:] = cs
for i in range(m) :
ret[m - i:] /= fac[:len(cs) + i]
ret[m - i - 1] += k[i] - polyval(lbnd, ret[m - i - 1:])
return ret
def polyval(x, cs):
"""Evaluate a polynomial.
If `cs` is of length `n`, this function returns :
``p(x) = cs[0] + cs[1]*x + ... + cs[n-1]*x**(n-1)``
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
If x is a list or tuple, it is converted to an ndarray. Otherwise
it must support addition and multiplication with itself and the
elements of `cs`.
cs : array_like
1-d array of Chebyshev coefficients ordered from low to high.
Returns
-------
values : ndarray
The return array has the same shape as `x`.
See Also
--------
polyfit
Examples
--------
Notes
-----
The evaluation uses Horner's method.
Examples
--------
"""
# cs is a trimmed copy
[cs] = pu.as_series([cs])
if isinstance(x, tuple) or isinstance(x, list) :
x = np.asarray(x)
c0 = cs[-1] + x*0
for i in range(2, len(cs) + 1) :
c0 = cs[-i] + c0*x
return c0
def polyvander(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. If ``V`` is the returned matrix and `x` is a 2d array, then
the elements of ``V`` are ``V[i,j,k] = x[i,j]**k``
Parameters
----------
x : array_like
Array of points. The values are converted to double or complex doubles.
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.
"""
x = np.asarray(x) + 0.0
order = int(deg) + 1
v = np.ones(x.shape + (order,), dtype=x.dtype)
if order > 1 :
v[...,1] = x
for i in range(2, order) :
v[...,i] = x*v[...,i-1]
return v
def polyfit(x, y, deg, rcond=None, full=False):
"""Least squares fit of polynomial to data.
Fit a polynomial ``p(x) = p[0] * T_{deq}(x) + ... + p[deg] *
T_{0}(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.
Returns
-------
coef : ndarray, shape (M,) or (M, K)
Polynomial 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
--------
polyval : Evaluates a polynomial.
polyvander : Vandermonde matrix for powers.
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 polynomial ``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 double precision and polynomials tend to fail at about
degree 20. Fits using Chebyshev series are generally better
conditioned, but much can still 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 x.shape[0] != y.shape[0] :
raise TypeError, "expected x and y to have same length"
# set rcond
if rcond is None :
rcond = len(x)*np.finfo(x.dtype).eps
# set up the design matrix and solve the least squares equation
A = polyvander(x, deg)
scl = np.sqrt((A*A).sum(0))
c, resids, rank, s = la.lstsq(A/scl, y, 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 polyroots(cs):
"""Roots of a polynomial.
Compute the roots of the Chebyshev series `cs`. The argument `cs` is a
sequence of coefficients ordered from low to high. i.e., [1,2,3] is the
polynomial ``1 + 2*x + 3*x**2``.
Parameters
----------
cs : array_like of shape(M,)
1D array of polynomial coefficients ordered from low to high.
Returns
-------
out : ndarray
An array containing the complex roots of the polynomial series.
Examples
--------
"""
# 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([-cs[0]/cs[1]])
n = len(cs) - 1
cmat = np.zeros((n,n), dtype=cs.dtype)
cmat.flat[n::n+1] = 1
cmat[:,-1] -= cs[:-1]/cs[-1]
roots = la.eigvals(cmat)
roots.sort()
return roots
#
# polynomial class
#
exec polytemplate.substitute(name='Polynomial', nick='poly', domain='[-1,1]')
|