""" Extended functions to operate on polynomials """ # ------------------------------------------------------------------------- # Name: polynomial # Purpose: Functions to operate on polynomials. # # Author: pab # polyXXX functions are based on functions found in the matlab toolbox polyutil # written by # Author: Peter J. Acklam # E-mail: pjacklam@online.no # WWW URL: http://home.online.no/~pjacklam # # Created: 30.12.2008 # Copyright: (c) pab 2008 # Licence: LGPL # ------------------------------------------------------------------------- # !/usr/bin/env python import warnings # @UnusedImport from numpy.polynomial import polyutils as pu from plotbackend import plotbackend as plt import numpy as np from numpy import (zeros, asarray, newaxis, arange, logical_or, any, pi, cos, round, diff, all, exp, where, extract, linalg, sign, concatenate, floor, linspace, sum, meshgrid) from scipy.fftpack import dct, idct as _idct from numpy.lib.polynomial import * # @UnusedWildImport from scipy.misc.common import pade # @UnresolvedImport __all__ = np.lib.polynomial.__all__ __all__ = __all__ + ['pade', 'padefit', 'polyreloc', 'polyrescl', 'polytrim', 'poly2hstr', 'poly2str', 'polyshift', 'polyishift', 'map_from_intervall', 'map_to_intervall', 'cheb2poly', 'chebextr', 'chebroot', 'chebpoly', 'chebfit', 'chebval', 'chebder', 'chebint', 'Cheb1d', 'dct', 'idct'] def polyint(p, m=1, k=None): """ Return an antiderivative (indefinite integral) of a polynomial. The returned order `m` antiderivative `P` of polynomial `p` satisfies :math:`\\frac{d^m}{dx^m}P(x) = p(x)` and is defined up to `m - 1` integration constants `k`. The constants determine the low-order polynomial part .. math:: \\frac{k_{m-1}}{0!} x^0 + \\ldots + \\frac{k_0}{(m-1)!}x^{m-1} of `P` so that :math:`P^{(j)}(0) = k_{m-j-1}`. Parameters ---------- p : {array_like, poly1d} Polynomial to differentiate. A sequence is interpreted as polynomial coefficients, see `poly1d`. m : int, optional Order of the antiderivative. (Default: 1) k : {None, list of `m` scalars, scalar}, optional Integration constants. They are given in the order of integration: those corresponding to highest-order terms come first. If ``None`` (default), all constants are assumed to be zero. If `m = 1`, a single scalar can be given instead of a list. See Also -------- polyder : derivative of a polynomial poly1d.integ : equivalent method Examples -------- The defining property of the antiderivative: >>> p = np.poly1d([1,1,1]) >>> P = np.polyint(p) >>> P poly1d([ 0.33333333, 0.5 , 1. , 0. ]) >>> np.polyder(P) == p True The integration constants default to zero, but can be specified: >>> P = np.polyint(p, 3) >>> P(0) 0.0 >>> np.polyder(P)(0) 0.0 >>> np.polyder(P, 2)(0) 0.0 >>> P = np.polyint(p, 3, k=[6, 5, 3]) >>> P.coefficients.tolist() [0.016666666666666666, 0.041666666666666664, 0.16666666666666666, 3.0, 5.0, 3.0] Note that 3 = 6 / 2!, and that the constants are given in the order of integrations. Constant of the highest-order polynomial term comes first: >>> np.polyder(P, 2)(0) 6.0 >>> np.polyder(P, 1)(0) 5.0 >>> P(0) 3.0 """ m = int(m) if m < 0: raise ValueError("Order of integral must be positive (see polyder)") if k is None: k = zeros(m, float) k = np.atleast_1d(k) if len(k) == 1 and m > 1: k = k[0] * np.ones(m, float) if len(k) < m: raise ValueError( "k must be a scalar or a rank-1 array of length 1 or >m.") truepoly = isinstance(p, poly1d) p = asarray(p) if m == 0: if truepoly: return poly1d(p) return p else: ix = arange(len(p), 0, -1) if p.ndim > 1: ix = ix[..., newaxis] pieces = p.shape[-1] k0 = k[0] * np.ones((1, pieces), dtype=int) else: k0 = [k[0]] y = np.concatenate((p.__truediv__(ix), k0), axis=0) val = polyint(y, m - 1, k=k[1:]) if truepoly: return poly1d(val) return val def polyder(p, m=1): """ Return the derivative of the specified order of a polynomial. Parameters ---------- p : poly1d or sequence Polynomial to differentiate. A sequence is interpreted as polynomial coefficients, see `poly1d`. m : int, optional Order of differentiation (default: 1) Returns ------- der : poly1d A new polynomial representing the derivative. See Also -------- polyint : Anti-derivative of a polynomial. poly1d : Class for one-dimensional polynomials. Examples -------- The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is: >>> p = np.poly1d([1,1,1,1]) >>> p2 = np.polyder(p) >>> p2 poly1d([3, 2, 1]) which evaluates to: >>> p2(2.) 17.0 We can verify this, approximating the derivative with ``(f(x + h) - f(x))/h``: >>> (p(2. + 0.001) - p(2.)) / 0.001 17.007000999997857 The fourth-order derivative of a 3rd-order polynomial is zero: >>> np.polyder(p, 2) poly1d([6, 2]) >>> np.polyder(p, 3) poly1d([6]) >>> np.polyder(p, 4) poly1d([ 0.]) """ m = int(m) if m < 0: raise ValueError("Order of derivative must be positive (see polyint)") truepoly = isinstance(p, poly1d) p = asarray(p) if m == 0: if truepoly: return poly1d(p) return p else: n = len(p) - 1 ix = arange(n, 0, -1) if p.ndim > 1: ix = ix[..., newaxis] y = ix * p[:-1] val = polyder(y, m - 1) if truepoly: return poly1d(val) return val def polydeg(x, y): ''' Return optimal degree for polynomial fitting N = POLYDEG(X,Y) finds the optimal degree for polynomial fitting, according to the Akaike's information criterion. Assuming that you want to find the degree N of a polynomial that fits the data Y(X) best in a least-squares sense, the Akaike's information criterion is defined by: 2*(N + 1) + n * (log(2 * pi * RSS / n) + 1) where n is the number of points and RSS is the residual sum of squares. The optimal degree N is defined here as that which minimizes AIC: http://en.wikipedia.org/wiki/Akaike_Information_Criterion Notes: ----- If the number of data is small, POLYDEG may tend to return: N = (number of points)-1. ORTHOFIT is more appropriate than POLYFIT for polynomial fitting with relatively high degrees. Example: ------- >>> x = np.linspace(0,10,300) >>> y = np.sin(x ** 3 / 100) ** 2 + 0.05 * np.random.randn(x.size) >>> n = polydeg(x,y) >>> n 21 ys = orthofit(x,y,n); plt.plot(x, y, '.', x, ys, 'k') See also -------- polyfit, orthofit ''' x, y = np.atleast_1d(x, y) x = x.ravel() y = y.ravel() N = len(x) # Search the optimal degree minimizing the Akaike's information criterion # y(x) are fitted in a least-squares sense using a polynomial of degree n # developed in a series of orthogonal polynomials. ys = np.ones((N,)) * y.mean() # correction for small sample sizes logsum2 = (np.log(2 * pi * ((ys - y) ** 2).sum() / N) + 1) AIC = 2 + N * logsum2 + 4 / (N - 2) n = 1 nit = 0 # While-loop is stopped when a minimum is detected. 3 more steps are # required to take AIC noise into account and to ensure that this minimum # is a (likely) global minimum. while nit < 3: p = orthofit(x, y, n) ys = orthoval(p, x) # -- Akaike's Information Criterion aic = (2 * (n + 1) * (1 + (n + 2) / (N - n - 2)) + N * (np.log(2 * pi * sum((ys - y) ** 2) / N) + 1)) if aic >= AIC: nit += 1 else: nit = 0 AIC = aic n = n + 1 if n >= N: break n = n - nit - 1 return n def orthoval(p, x): ''' Evaluation of orthogonal polynomial Parameters ---------- p : array_like 2D array of polynomial coefficients (including coefficients equal to zero) from highest degree to the constant term. x : array_like A number or a 1D array of numbers "at" which to evaluate `p`. Returns ------- values : ndarray See Also -------- orthofit ''' p = np.atleast_2d(p) n = p.shape[1] - 1 xi = np.atleast_1d(x) shape0 = xi.shape if n == 0: return np.ones(shape0) * p[0] xi = xi.ravel() xn = np.ones((n + 1, len(xi))) xn[1] = xi - p[1, 1] for i in range(2, n + 1): xn[i, :] = (xi - p[1, i]) * xn[i - 1, :] - p[2, i] * xn[i - 2, :] ys = np.dot(p[0], xn) return ys.reshape(shape0) def ortho2poly(p): """ Converts orthogonal polynomial to ordinary polynomial coefficients Parameters ---------- p : array-like orthogonal polynomial coefficients Returns ------- p : ndarray ordinary polynomial coefficients It is not advised to do this for p.shape[1]>10 due to numerical cancellations. See also -------- orthoval orthofit Examples -------- >>> import numpy as np >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0]) >>> p = orthofit(x, y, 3) >>> p array([[ 0. , -0.30285714, -0.16071429, 0.08703704], [ 0. , 2.5 , 2.5 , 2.5 ], [ 0. , 0. , 2.91666667, 2.13333333]]) >>> ortho2poly(p) array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) >>> np.polyfit(x, y, 3) array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) References ---------- """ p = np.atleast_2d(p) n = p.shape[1] - 1 if n == 0: return p[0] x = [1, ] * (n + 1) x[1] = np.array([1, - p[1, 1]]) for i in range(2, n + 1): x[i] = polyadd(polymul([1, - p[1, i]], x[i - 1]), - p[2, i] * x[i - 2]) for i in range(n + 1): x[i] *= p[0, i] return reduce(polyadd, x) def orthofit(x, y, n): ''' Fit orthogonal polynomial to data. Parameters --------- x, y : arrays data Y(X) to fit to a polynomial. n : integer degree of fitted polynomial. Returns ------- p : array orthogonal polynomial Notes: ----- Orthofit smooths/fits data using a polynomial of degree N developed in a sequence of orthogonal polynomials. ORTHOFIT is more appropriate than polyfit for polynomial fitting and smoothing since this method does not involve any matrix linear system but a simple recursive procedure. Degrees much higher than 30 could be used with orthogonal polynomials, whereas badly conditioned matrices may appear with a classical polynomial fitting of degree typically higher than 10. To avoid using unnecessarily high degrees, you may let the function POLYDEG choose it for you. POLYDEG finds an optimal polynomial degree according to the Akaike's information criterion. Example: ------- >>> x = np.linspace(0,10,300); >>> y = np.sin(x**3/100)**2 + 0.05*np.random.randn(x.size) >>> p = orthofit(x, y, 25) >>> ys = orthoval(p, x) plot(x, y,'.',x, ys, 'k') See also -------- polydeg, polyfit, polyval Reference: --------- Methodes de calcul numerique 2. JP Nougier. Hermes Science Publications, 2001. Section 4.7 pp 116-121 ''' x, y = np.atleast_1d(x, y) x = x.ravel() y = y.ravel() # Particular case: n=0 if n == 0: return y.mean() # p = Coefficients of the orthogonal polynomials p = np.zeros((3, n + 1)) p[1, 1] = x.mean() N = len(x) PL = np.ones((n + 1, N)) PL[1] = x - p[1, 1] for i in range(2, n + 1): p[1, i] = np.dot(x, PL[i - 1] ** 2) / sum(PL[i - 1] ** 2) p[2, i] = np.dot(x, PL[i - 2] * PL[i - 1]) / sum(PL[i - 2] ** 2) PL[i] = (x - p[1, i]) * PL[i - 1] - p[2, i] * PL[i - 2] p[0, :] = np.dot(PL, y) / sum(PL ** 2, axis=1) return p # ys = np.dot(p[0, :], PL) # smoothed y def polyreloc(p, x, y=0.0): """ Relocate polynomial The polynomial `p` is relocated by "moving" it `x` units along the x-axis and `y` units along the y-axis. So the polynomial `r` is relative to the point (x,y) as the polynomial `p` is relative to the point (0,0). Parameters ---------- p : array-like, poly1d vector or matrix of column vectors of polynomial coefficients to relocate. (Polynomial coefficients are in decreasing order.) x : scalar distance to relocate P along x-axis y : scalar distance to relocate P along y-axis (default 0) Returns ------- r : ndarray, poly1d vector/matrix/poly1d of relocated polynomial coefficients. See also -------- polyrescl Example ------- >>> import numpy as np >>> p = np.arange(6); p.shape = (2,-1) >>> np.polyval(p,0) array([3, 4, 5]) >>> np.polyval(p,1) array([3, 5, 7]) >>> r = polyreloc(p,-1) # move to the left along x-axis >>> np.polyval(r,-1) # = polyval(p,0) array([3, 4, 5]) >>> np.polyval(r,0) # = polyval(p,1) array([3, 5, 7]) """ truepoly = isinstance(p, poly1d) r = np.atleast_1d(p).copy() n = r.shape[0] # Relocate polynomial using Horner's algorithm for ii in range(n, 1, -1): for i in range(1, ii): r[i] = r[i] - x * r[i - 1] r[-1] = r[-1] + y if r.ndim > 1 and r.shape[-1] == 1: r.shape = (r.size,) if truepoly: r = poly1d(r) return r def polyrescl(p, x, y=1.0): """ Rescale polynomial. Parameters ---------- p : array-like, poly1d vector or matrix of column vectors of polynomial coefficients to rescale. (Polynomial coefficients are in decreasing order.) x,y : scalars defining the factors to rescale the polynomial `p` in x-direction and y-direction, respectively. Returns ------- r : ndarray, poly1d vector/matrix/poly1d of rescaled polynomial coefficients. See also -------- polyreloc Example ------- >>> import numpy as np >>> p = np.arange(6); p.shape = (2,-1) >>> np.polyval(p,0) array([3, 4, 5]) >>> np.polyval(p,1) array([3, 5, 7]) >>> r = polyrescl(p,2) # scale by 2 along x-axis >>> np.polyval(r,0) # = polyval(p,0) array([ 3., 4., 5.]) >>> np.polyval(r,2) # = polyval(p,1) array([ 3., 5., 7.]) """ truepoly = isinstance(p, poly1d) r = np.atleast_1d(p) n = r.shape[0] xscale = (float(x) ** arange(1 - n, 1)) if r.ndim == 1: q = y * r * xscale else: q = y * r * xscale[:, newaxis] if truepoly: q = poly1d(q) return q def polytrim(p): """ Trim polynomial by stripping off leading zeros. Parameters ---------- p : array-like, poly1d vector or matrix of column vectors of polynomial coefficients in decreasing order. Returns ------- r : ndarray, poly1d vector/matrix/poly1d of trimmed polynomial coefficients. Example ------- >>> p = [0,1,2] >>> polytrim(p) array([1, 2]) >>> p1 = [[0,0],[1,2],[3,4]] >>> polytrim(p1) array([[1, 2], [3, 4]]) """ truepoly = isinstance(p, poly1d) if truepoly: return p else: r = np.atleast_1d(p).copy() # Remove leading zeros is_not_lead_zeros = logical_or.accumulate(r != 0, axis=0) if r.ndim == 1: r = r[is_not_lead_zeros] else: is_not_lead_zeros = any(is_not_lead_zeros, axis=1) r = r[is_not_lead_zeros, :] return r def poly2hstr(p, variable='x'): """ Return polynomial as a Horner represented string. Parameters ---------- p : array-like poly1d vector of polynomial coefficients in decreasing order. variable : string display character for variable Returns ------- p_str : string consisting of the polynomial coefficients in the vector P multiplied by powers of the given `variable`. Examples -------- >>> poly2hstr([1, 1, 2], 's' ) '(s + 1)*s + 2' See also -------- poly2str """ var = variable coefs = polytrim(np.atleast_1d(p)) order = len(coefs) - 1 # Order of polynomial. s = '' # Initialize output string. ix = 1 for expon in range(order, -1, -1): coef = coefs[order - expon] # There is no point in adding a zero term (except if it's the only # term, but we'll take care of that later). if coef == 0: ix += 1 else: # Append exponent if necessary. if ix > 1: exponstr = '%.0f' % ix s = '%s**%s' % (s, exponstr) ix = 1 # Is it the first term? isfirst = s == '' # We need the coefficient only if it is different from 1 or -1 or # when it is the constant term. needcoef = ( (abs(coef) != 1) | ( expon == 0) & isfirst) | 1 - isfirst # We need the variable except in the constant term. needvar = (expon != 0) # Add sign, but we don't need a leading plus-sign. if isfirst: if coef < 0: s = '-' # % Unary minus. else: if coef < 0: s = '%s - ' % s # % Binary minus (subtraction). else: s = '%s + ' % s # % Binary plus (addition). # Append the coefficient if it is different from one or when it is # the constant term. if needcoef: coefstr = '%.20g' % abs(coef) s = '%s%s' % (s, coefstr) # Append variable if necessary. if needvar: # Append a multiplication sign if necessary. if needcoef: if 1 - isfirst: s = '(%s)' % s s = '%s*' % s s = '%s%s' % (s, var) # Now treat the special case where the polynomial is zero. if s == '': s = '0' return s def poly2str(p, variable='x'): """ Return polynomial as a string. Parameters ---------- p : array-like poly1d vector of polynomial coefficients in decreasing order. variable : string display character for variable Returns ------- p_str : string consisting of the polynomial coefficients in the vector P multiplied by powers of the given `variable`. See also -------- poly2hstr Examples -------- >>> poly2str([1, 1, 2], 's' ) 's**2 + s + 2' """ thestr = "0" var = variable # Remove leading zeros coeffs = polytrim(np.atleast_1d(p)) N = len(coeffs) - 1 for k in range(len(coeffs)): coefstr = '%.4g' % abs(coeffs[k]) if coefstr[-4:] == '0000': coefstr = coefstr[:-5] power = (N - k) if power == 0: if coefstr != '0': newstr = '%s' % (coefstr,) else: if k == 0: newstr = '0' else: newstr = '' elif power == 1: if coefstr == '0': newstr = '' elif coefstr == 'b' or coefstr == '1': newstr = var else: newstr = '%s*%s' % (coefstr, var) else: if coefstr == '0': newstr = '' elif coefstr == 'b' or coefstr == '1': newstr = '%s**%d' % (var, power,) else: newstr = '%s*%s**%d' % (coefstr, var, power) if k > 0: if newstr != '': if coeffs[k] < 0: thestr = "%s - %s" % (thestr, newstr) else: thestr = "%s + %s" % (thestr, newstr) elif (k == 0) and (newstr != '') and (coeffs[k] < 0): thestr = "-%s" % (newstr,) else: thestr = newstr return thestr def polyshift(py, a=-1, b=1): """ Polynomial coefficient shift Polyshift shift the polynomial coefficients by a variable shift: Y = 2*(X-.5*(b+a))/(b-a) i.e., the interval -1 <= Y <= 1 is mapped to the interval a <= X <= b Parameters ---------- py : array-like polynomial coefficients for the variable y. a,b : scalars lower and upper limit. Returns ------- px : ndarray polynomial coefficients for the variable x. See also -------- polyishift Example ------- >>> py = [1, 0] >>> px = polyshift(py,0,5) >>> polyval(px,[0, 2.5, 5]) #% This is the same as the line below array([-1., 0., 1.]) >>> polyval(py,[-1, 0, 1 ]) array([-1, 0, 1]) """ if (a == -1) & (b == 1): return py L = b - a return polyishift(py, -(2. + b + a) / L, (2. - b - a) / L) def polyishift(px, a=-1, b=1): """ Inverse polynomial coefficient shift Polyishift does the inverse of Polyshift, shift the polynomial coefficients by a variable shift: Y = 2*(X-.5*(b+a)/(b-a) i.e., the interval a <= X <= b is mapped to the interval -1 <= Y <= 1 Parameters ---------- px : array-like polynomial coefficients for the variable x. a,b : scalars lower and upper limit. Returns ------- py : ndarray polynomial coefficients for the variable y. See also -------- polyishift Example ------- >>> px = [1, 0] >>> py = polyishift(px,0,5); >>> polyval(px,[0, 2.5, 5]) #% This is the same as the line below array([ 0. , 2.5, 5. ]) >>> polyval(py,[-1, 0, 1]) array([ 0. , 2.5, 5. ]) """ if (a == -1) & (b == 1): return px L = b - a xscale = 2. / L xloc = -float(a + b) / L return polyreloc(polyrescl(px, xscale), xloc) def map_from_interval(x, a, b): """F(x), where F: [a,b] -> [-1,1].""" return (x - (b + a) / 2.0) * (2.0 / (b - a)) def map_to_interval(x, a, b): """F(x), where F: [-1,1] -> [a,b].""" return (x * (b - a) + (b + a)) / 2.0 def poly2cheb(p, a=-1, b=1): """ Convert polynomial coefficients into Chebyshev coefficients Parameters ---------- p : array-like polynomial coefficients a,b : real scalars lower and upper limits (Default -1,1) Returns ------- ck : ndarray Chebychef coefficients POLY2CHEB do the inverse of CHEB2POLY: given a vector of polynomial coefficients AK, returns an equivalent vector of Chebyshev coefficients CK. This is useful for economization of power series. The steps for doing so: 1. Convert polynomial coefficients to Chebychev coefficients, CK. 2. Truncate the CK series to a smaller number of terms, using the coefficient of the first neglected Chebychev polynomial as an error estimate. 3 Convert back to a polynomial by CHEB2POLY See also -------- cheb2poly chebval chebfit Examples -------- >>> import numpy as np >>> p = np.arange(5) >>> ck = poly2cheb(p) >>> cheb2poly(ck) array([ 1., 2., 3., 4.]) Reference --------- William H. Press, Saul Teukolsky, William T. Wetterling and Brian P. Flannery (1997) "Numerical recipes in Fortran 77", Vol. 1, pp 184-194 """ f = poly1d(p) n = len(f.coeffs) return chebfit(f, n, a, b) def cheb2poly(ck, a=-1, b=1): """ Converts Chebyshev coefficients to polynomial coefficients Parameters ---------- ck : array-like Chebychef coefficients a,b : real, scalars lower and upper limits (Default -1,1) Returns ------- p : ndarray polynomial coefficients It is not advised to do this for len(ck)>10 due to numerical cancellations. See also -------- chebval chebfit Examples -------- >>> import numpy as np >>> p = np.arange(5) >>> ck = poly2cheb(p) >>> cheb2poly(ck) array([ 1., 2., 3., 4.]) References ---------- http://en.wikipedia.org/wiki/Chebyshev_polynomials http://en.wikipedia.org/wiki/Chebyshev_form http://en.wikipedia.org/wiki/Clenshaw_algorithm """ n = len(ck) b_Nmi = zeros(1) b_Nmip1 = zeros(1) y = np.r_[2 / (b - a), -(a + b) / (b - a)] y2 = 2. * y # Clenshaw recurence for ix in xrange(n - 1): tmp = b_Nmi b_Nmi = polymul(y2, b_Nmi) # polynomial multiplication nb = len(b_Nmip1) b_Nmip1[-1] = b_Nmip1[-1] - ck[ix] b_Nmi[-nb::] = b_Nmi[-nb::] - b_Nmip1 b_Nmip1 = tmp p = polymul(y, b_Nmi) # polynomial multiplication nb = len(b_Nmip1) b_Nmip1[-1] = b_Nmip1[-1] - ck[n - 1] p[-nb::] = p[-nb::] - b_Nmip1 return polytrim(p) def chebextr(n): """ Return roots of derivative of Chebychev polynomial of the first kind. All local extreme values of the polynomial are either -1 or 1. So, CHEBPOLY( N, CHEBEXTR(N) ) ) return the same as (-1).^(N:-1:0) except for the numerical noise in the former. Because the extreme values of Chebychev polynomials of the first kind are either -1 or 1, their roots are often used as starting values for the nodes in minimax approximations. Parameters ---------- n : scalar, integer degree of Chebychev polynomial. Examples -------- >>> x = chebextr(4) >>> chebpoly(4,x) array([ 1., -1., 1., -1., 1.]) Reference --------- http://en.wikipedia.org/wiki/Chebyshev_nodes http://en.wikipedia.org/wiki/Chebyshev_polynomials """ return - cos((pi * arange(n + 1)) / n) def chebroot(n, kind=1): """ Return roots of Chebychev polynomial of the first or second kind. The roots of the Chebychev polynomial of the first kind form a particularly good set of nodes for polynomial interpolation because the resulting interpolation polynomial minimizes the problem of Runge's phenomenon. Parameters ---------- n : scalar, integer degree of Chebychev polynomial. kind: 1 or 2, optional kind of Chebychev polynomial (default 1) Examples -------- >>> import numpy as np >>> x = chebroot(3) >>> np.abs(chebpoly(3,x))<1e-15 array([ True, True, True], dtype=bool) >>> chebpoly(3) array([ 4., 0., -3., 0.]) >>> x2 = chebroot(4,kind=2) >>> np.abs(chebpoly(4,x2,kind=2))<1e-15 array([ True, True, True, True], dtype=bool) >>> chebpoly(4,kind=2) array([ 16., 0., -12., 0., 1.]) Reference --------- http://en.wikipedia.org/wiki/Chebyshev_nodes http://en.wikipedia.org/wiki/Chebyshev_polynomials """ if kind not in (1, 2): raise ValueError('kind must be 1 or 2') return - cos(pi * (arange(n) + 0.5 * kind) / (n + kind - 1)) def chebpoly(n, x=None, kind=1): """ Return Chebyshev polynomial of the first or second kind. These polynomials are orthogonal on the interval [-1,1], with respect to the weight function w(x) = (1-x**2)**(-1/2+kind-1). chebpoly(n) returns coefficients of the Chebychev polynomial of degree N. chebpoly(n,x) returns the Chebychev polynomial of degree N evaluated at X. Parameters ---------- n : integer, scalar degree of Chebychev polynomial. x : array-like, optional evaluation points kind: 1 or 2, optional kind of Chebychev polynomial (default 1) Returns ------- p : ndarray polynomial coefficients if x is None. Chebyshev polynomial evaluated at x otherwise Examples -------- >>> import numpy as np >>> x = chebroot(3) >>> np.abs(chebpoly(3,x))<1e-15 array([ True, True, True], dtype=bool) >>> chebpoly(3) array([ 4., 0., -3., 0.]) >>> x2 = chebroot(4,kind=2) >>> np.abs(chebpoly(4,x2,kind=2))<1e-15 array([ True, True, True, True], dtype=bool) >>> chebpoly(4,kind=2) array([ 16., 0., -12., 0., 1.]) Reference --------- http://en.wikipedia.org/wiki/Chebyshev_polynomials """ if x is None: # Calculate coefficients. if n == 0: p = np.ones(1) else: p = round(pow(2, n - 2 + kind) * poly(chebroot(n, kind=kind))) p[1::2] = 0 return p else: # Evaluate polynomial in chebychev form ck = zeros(n + 1) ck[0] = 1. return _chebval(np.atleast_1d(x), ck, kind=kind) def chebfit(fun, n=10, a=-1, b=1, trace=False): """ Computes the Chebyshevs coefficients so that f(x) can be approximated by: n-1 f(x) = sum ck*Tk(x) k=0 where Tk is the k'th Chebyshev polynomial of the first kind. Parameters ---------- fun : callable function to approximate n : integer, scalar, optional number of base points (abscissas). Default n=10 (maximum 50) a,b : real, scalars, optional integration limits Returns ------- ck : ndarray polynomial coefficients in Chebychev form. Examples -------- Fit exp(x) >>> import matplotlib.pyplot as plt >>> a = 0; b = 2 >>> ck = chebfit(np.exp,7,a,b); >>> x = np.linspace(0,4); >>> h=plt.plot(x, np.exp(x), 'r', x, chebval(x,ck,a,b), 'g.') >>> x1 = chebroot(9)*(b-a)/2+(b+a)/2 >>> ck1 = chebfit(np.exp(x1)) >>> h = plt.plot(x,np.exp(x), 'r', x, chebval(x,ck1,a,b),'g.') >>> plt.close() See also -------- chebval Reference --------- http://en.wikipedia.org/wiki/Chebyshev_nodes http://mathworld.wolfram.com/ChebyshevApproximationFormula.html W. Fraser (1965) "A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable" Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314 """ if (n > 50): warnings.warn('CHEBFIT should only be used for n<50') if hasattr(fun, '__call__'): x = map_to_interval(chebroot(n), a, b) f = fun(x) if trace: plt.plot(x, f, '+') else: f = fun n = len(f) # N-1 # c[k] = (2/N) sum w[n] f[n]*cos(pi*k*(2n+1)/(2N)), 0 <= k < N. # n=0 # # w[0] = 0.5, w[n]=1 for n>0 ck = dct(f[::-1]) / n ck[0] = ck[0] / 2. return ck[::-1] def chebfit_dct(f, n=(10, ), domain=None): """ Fit Chebyshev series to N-dimensional function so that f(x1, x2,..., xn) can be approximated by: .. math:: f(x_1, x_2,...,x_n) = \\sum_{i,j,...k} c_i T_i(x_1)*...*c_k T_k(x_n) , where Tk is the k'th Chebyshev polynomial of the first kind. Parameters ---------- f : callable function to approximate n : list of integers, optional number of base points (abscissas) used for each dimension. Default n=10 (maximum 50) domain : vector [a1,b1,a2,b2 ,..., an, bn], optional defining the rectangle [a1, b1] x [a2, b2] x ...x [an, bn]. (default domain = (-1,1) * len(n)) Returns ------- ck : ndarray polynomial coefficients in Chebychev form. Examples -------- Fit exponential function >>> import matplotlib.pyplot as plt >>> domain = (0, 2) >>> ck = chebfit_dct(np.exp, 7, domain) >>> np.allclose(ck, [3.44152387e+00, 3.07252345e+00, 7.38000848e-01, ... 1.20520053e-01, 1.48805268e-02, 1.47579673e-03, ... 1.21719524e-04]) True >>> x1 = map_to_interval(chebroot(9), *domain) >>> ck1 = chebfit(np.exp(x1)) >>> np.allclose(ck1, [5.40019009e-07, 8.69418381e-06, 1.22261037e-04, ... 1.47582673e-03, 1.48805283e-02, 1.20520053e-01, ... 7.38000848e-01, 3.07252345e+00, 3.44152387e+00]) True x = np.linspace(0,4) h = plt.plot(x, np.exp(x), 'r', x, chebvalnd(ck, x,ck,a,b), 'g.') h = plt.plot(x, np.exp(x), 'r', x, chebvalnd(ck1, x,ck1,a,b),'b.') plt.close() See also -------- chebval, chebvalnd Reference --------- http://en.wikipedia.org/wiki/Chebyshev_nodes http://mathworld.wolfram.com/ChebyshevApproximationFormula.html W. Fraser (1965) "A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable" Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314 """ n = np.atleast_1d(n) if np.any(n > 50): warnings.warn('CHEBFIT should only be used for n<50') if hasattr(f, '__call__'): if domain is None: domain = (-1, 1) * len(n) domain = np.atleast_2d(domain).reshape((-1, 2)) xi = [map_to_interval(chebroot(ni), d[0], d[1]) for ni, d in zip(n, domain)] Xi = np.meshgrid(*xi) ck = f(*Xi) else: ck = f n = ck.shape ndim = len(n) for i in range(ndim): ck = dct(ck[..., ::-1]) ck[..., 0] = ck[..., 0] / 2. if i < ndim-1: ck = np.rollaxis(ck, axis=-1) return ck / np.product(n) def idct(x, n=None): """ Inverse Discrete Cosine Transform N-1 x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. n=0 w(0) = 1/2 w(n) = 1 for n>0 Examples -------- >>> import numpy as np >>> x = np.arange(5) >>> np.abs(x-idct(dct(x)))<1e-14 array([ True, True, True, True, True], dtype=bool) >>> np.abs(x-dct(idct(x)))<1e-14 array([ True, True, True, True, True], dtype=bool) Reference --------- http://en.wikipedia.org/wiki/Discrete_cosine_transform http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ """ return _idct(x, n=n, norm=None)*0.5/len(x) def _chebval(x, ck, kind=1): """ Evaluate polynomial in Chebyshev form. A polynomial of degree N in Chebyshev form is a polynomial p(x): N p(x) = sum ck*Tk(x) k=0 or N p(x) = sum ck*Uk(x) k=0 where Tk and Uk are the k'th Chebyshev polynomial of the first and second kind, respectively. References ---------- http://en.wikipedia.org/wiki/Clenshaw_algorithm http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html """ n = len(ck) b_Nmi = zeros(x.shape) # b_(N-i) b_Nmip1 = b_Nmi.copy() # b_(N-i+1) x2 = 2 * x # Clenshaw reccurence for ix in xrange(n - 1): tmp = b_Nmi b_Nmi = x2 * b_Nmi - b_Nmip1 + ck[ix] b_Nmip1 = tmp return kind * x * b_Nmi - b_Nmip1 + ck[n - 1] def chebval(x, ck, a=-1, b=1, kind=1, fill=None): """ Evaluate polynomial in Chebyshev form at X A polynomial of degree N in Chebyshev form is a polynomial p(x) of the form N p(x) = sum ck*Tk(x) k=0 where Tk is the k'th Chebyshev polynomial of the first or second kind. Paramaters ---------- x : array-like points to evaluate ck : array-like polynomial coefficients in Chebyshev form ordered from highest degree to zero a,b : real, scalars, optional limits for polynomial (Default -1,1) kind: 1 or 2, optional kind of Chebychev polynomial (default 1) fill : scalar, optional If provided, define value to return for `x < a` or `b < x`. Examples -------- Plot Chebychev polynomial of the first kind and order 4: >>> import matplotlib.pyplot as plt >>> x = np.linspace(-1,1) >>> ck = np.zeros(5); ck[-1]=1 >>> h = plt.plot(x,chebval(x,ck),x,chebpoly(4,x),'.') >>> plt.close() Fit exponential function: >>> import matplotlib.pyplot as plt >>> ck = chebfit(np.exp,7,0,2) >>> x = np.linspace(0,4); >>> h=plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x)) >>> plt.close() See also -------- chebfit References ---------- http://en.wikipedia.org/wiki/Clenshaw_algorithm http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html """ y = map_from_interval(atleast_1d(x), a, b) if fill is None: f = _chebval(y, ck, kind=kind) else: cond = (abs(y) <= 1) f = where(cond, 0, fill) if any(cond): yk = extract(cond, y) f[cond] = _chebval(yk, ck, kind=kind) return f def chebder(ck, a=-1, b=1): """ Differentiate Chebyshev polynomial Parameters ---------- ck : array-like polynomial coefficients in Chebyshev form of function to differentiate a,b : real, scalars limits for polynomial(Default -1,1) Return ------ cder : ndarray polynomial coefficients in Chebyshev form of the derivative Examples -------- Fit exponential function: >>> import matplotlib.pyplot as plt >>> ck = chebfit(np.exp,7,0,2) >>> x = np.linspace(0,4) >>> ck2 = chebder(ck,0,2); >>> h = plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x),'r') >>> plt.close() See also -------- chebint chebfit Reference --------- http://en.wikipedia.org/wiki/Chebyshev_polynomials W. Fraser (1965) "A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable" Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314 """ n = len(ck) - 1 cder = zeros(n, dtype=asarray(ck).dtype) cder[0] = 2 * n * ck[0] cder[1] = 2 * (n - 1) * ck[1] for j in xrange(2, n): cder[j] = cder[j - 2] + 2 * (n - j) * ck[j] return cder * 2. / (b - a) # Normalize to the interval b-a. def chebint(ck, a=-1, b=1): """ Integrate Chebyshev polynomial Parameters ---------- ck : array-like polynomial coefficients in Chebyshev form of function to integrate. a,b : real, scalars limits for polynomial(Default -1,1) Return ------ cint : ndarray polynomial coefficients in Chebyshev form of the integrated function Examples -------- Fit exponential function: >>> import matplotlib.pyplot as plt >>> ck = chebfit(np.exp,7,0,2) >>> x = np.linspace(0,4) >>> ck2 = chebint(ck,0,2); >>> h=plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x),'r.') >>> plt.close() See also -------- chebder chebfit Reference --------- http://en.wikipedia.org/wiki/Chebyshev_polynomials W. Fraser (1965) "A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable" Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314 """ # int T0(x) = T1(x)+1 # int T1(x) = 0.5*(T2(x)/2-T0/2) # int Tn(x) dx = 0.5*{Tn+1(x)/(n+1) - Tn-1(x)/(n-1)} # N # p(x) = sum cn*Tn(x) # n=0 # int p(x) dx = sum cn * int(Tn(x)dx) = # 0.5*sum cn *{Tn+1(x)/(n+1) - Tn-1(x)/(n-1)} = 0.5 sum (cn-1-cn+1)*Tn/n n>0 n = len(ck) cint = zeros(n) con = 0.25 * (b - a) dif1 = diff(ck[-1::-2]) ix1 = np.r_[1:n - 1:2] cint[ix1] = -(con * dif1) / ix1 if n > 3: dif2 = diff(ck[-2::-2]) ix2 = np.r_[2:n - 1:2] cint[ix2] = -(con * dif2) / ix2 cint = cint[::-1] # cint(n) is a special case cint[-1] = (con * ck[n - 2]) / (n - 1) # Set integration constant cint[0] = 2 * np.sum((-1) ** np.r_[0:n - 1] * cint[-2::-1]) return cint class Cheb1d(object): coeffs = None order = None a = None b = None kind = None def __init__(self, ck, a=-1, b=1, kind=1): if isinstance(ck, Cheb1d): for key in ck.__dict__.keys(): self.__dict__[key] = ck.__dict__[key] return cki = trim_zeros(atleast_1d(ck), 'b') if len(cki.shape) > 1: raise ValueError("Polynomial must be 1d only.") self.__dict__['coeffs'] = cki self.__dict__['order'] = len(cki) - 1 self.__dict__['a'] = a self.__dict__['b'] = b self.__dict__['kind'] = kind def __call__(self, x): return chebval(x, self.coeffs, self.a, self.b, self.kind) def __array__(self, t=None): if t: return asarray(self.coeffs, t) else: return asarray(self.coeffs) def __repr__(self): vals = repr(self.coeffs) vals = vals[6:-1] return "Cheb1d(%s)" % vals def __len__(self): return self.order def __str__(self): pass def __neg__(self): new = Cheb1d(self) new.coeffs = -self.coeffs return new def __pos__(self): return self def __add__(self, other): other = Cheb1d(other) new = Cheb1d(self) new.coeffs = polyadd(self.coeffs, other.coeffs) return new def __radd__(self, other): return self.__add__(other) def __sub__(self, other): other = Cheb1d(other) new = Cheb1d(self) new.coeffs = polysub(self.coeffs, other.coeffs) return new def __rsub__(self, other): other = Cheb1d(other) new = Cheb1d(self) new.coeffs = polysub(other.coeffs, new.coeffs) return new def __eq__(self, other): other = Cheb1d(other) return (all(self.coeffs == other.coeffs) and (self.a == other.a) and (self.b == other.b) and (self.kind == other.kind)) def __ne__(self, other): return any(self.coeffs != other.coeffs) or (self.a != other.a) or ( self.b != other.b) or (self.kind != other.kind) def __setattr__(self, key, val): raise ValueError("Attributes cannot be changed this way.") def __getattr__(self, key): if key in ['c', 'coef', 'coefficients']: return self.coeffs elif key in ['o']: return self.order elif key in ['a']: return self.a elif key in ['b']: return self.b elif key in ['k']: return self.kind else: try: return self.__dict__[key] except KeyError: raise AttributeError( "'%s' has no attribute '%s'" % (self.__class__, key)) def __getitem__(self, val): if val > self.order: return 0 if val < 0: return 0 return self.coeffs[val] def __setitem__(self, key, val): # ind = self.order - key if key < 0: raise ValueError("Does not support negative powers.") if key > self.order: zr = zeros(key - self.order, self.coeffs.dtype) self.__dict__['coeffs'] = concatenate((self.coeffs, zr)) self.__dict__['order'] = key self.__dict__['coeffs'][key] = val return def __iter__(self): return iter(self.coeffs) def integ(self, m=1): """ Return an antiderivative (indefinite integral) of this polynomial. Refer to `chebint` for full documentation. See Also -------- chebint : equivalent function """ integ = Cheb1d(self) integ.coeffs = chebint(self.coeffs, self.a, self.b) return integ def deriv(self, m=1): """ Return a derivative of this polynomial. Refer to `chebder` for full documentation. See Also -------- chebder : equivalent function """ der = Cheb1d(self) der.coeffs = chebder(self.coeffs, self.a, self.b) return der def padefit(c, m=None): """ Rational polynomial fitting from polynomial coefficients Parameters ---------- c : array-like coefficients of power series expansion from highest degree to zero. m : scalar integer order of denominator polynomial. (Default floor((len(c)-1)/2)) Returns ------- num, den : poly1d numerator and denominator polynomials for the pade approximation If the function is well approximated by M+N+1 f(x) = sum c(2*n+2-k)*x^k k=0 then the pade approximation is given by M sum c1(n-k+1)*x^k k=0 f(x) = ------------------------ N sum c2(n-k+1)*x^k k=0 Note: c must be ordered for direct use with polyval Example ------- Pade approximation to exp(x) >>> import scipy.special as sp >>> import matplotlib.pyplot as plt >>> c = poly1d(1./sp.gamma(np.r_[6+1:0:-1])) >>> [p, q] = padefit(c) >>> p; q poly1d([ 0.00277778, 0.03333333, 0.2 , 0.66666667, 1. ]) poly1d([ 0.03333333, -0.33333333, 1. ]) >>> x = np.linspace(0,4); >>> h = plt.plot(x,c(x),x,p(x)/q(x),'g-', x,np.exp(x),'r.') >>> plt.close() See also -------- scipy.misc.pade """ if not m: m = int(floor((len(c) - 1) * 0.5)) c = asarray(c) return pade(c[::-1], m) def test_pade(): cof = array(([1.0, 1.0, 1.0 / 2, 1. / 6, 1. / 24])) p, q = pade(cof, 2) t = arange(0, 2, 0.1) assert(all(abs(p(t) / q(t) - exp(t)) < 0.3)) def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True): """ Rational polynomial fitting. A minimax solution by least squares. Parameters ---------- fun : callable or or a two column matrix f=[x,f(x)] where length(x)>(m+k+1)*8. m, k : integer number of coefficients of the numerator and denominater, respectively. a, b : real scalars evaluation limits, (default a=-1,b=1) Returns ------- num, den : poly1d numerator and denominator polynomials for the pade approximation dev : ndarray maximum absolute deviation of the approximation The pade approximation is given by m sum c1[m-i]*x**i i=0 f(x) = ------------------------ k sum c2[k-i]*x**i i=0 If F is a two column matrix, [x f(x)], a good choice for x is: x = cos(pi/(N-1)*(N-1:-1:0))*(b-a)/2+ (a+b)/2, where N = (m+k+1)*8; Note: c1 and c2 are ordered for direct use with polyval Example ------- Pade approximation to exp(x) between 0 and 2 >>> import matplotlib.pyplot as plt >>> [c1, c2] = padefitlsq(np.exp,3,3,0,2) >>> c1; c2 poly1d([ 0.01443847, 0.128842 , 0.55284547, 0.99999962]) poly1d([-0.0049658 , 0.07610473, -0.44716929, 1. ]) >>> x = np.linspace(0,4) >>> h = plt.plot(x, polyval(c1,x)/polyval(c2,x),'g') >>> h = plt.plot(x, np.exp(x), 'r') See also -------- padefit Reference --------- William H. Press, Saul Teukolsky, William T. Wetterling and Brian P. Flannery (1997) "Numerical recipes in Fortran 77", Vol. 1, pp 197-20 """ NFAC = 8 BIG = 1e30 MAXIT = 5 smallest_devmax = BIG ncof = m + k + 1 # Number of points where function is evaluated, i.e. fineness of mesh npt = NFAC * ncof if x is None: if end_points: # Use the location of the local extreme values of # the Chebychev polynomial of the first kind of degree NPT-1. x = map_to_interval(chebextr(npt - 1), a, b) else: # Use the roots of the Chebychev polynomial of the first kind of # degree NPT. Note this is useful if there are singularities close # to the endpoints. x = map_to_interval(chebroot(npt, kind=1), a, b) if hasattr(fun, '__call__'): fs = fun(x) else: fs = fun n = len(fs) if n < npt: warnings.warn( 'Check the result! ' + 'Number of function values should be at least: %d' % npt) if trace: plt.plot(x, fs, '+') wt = np.ones((npt)) ee = np.ones((npt)) mad = 0 u = zeros((npt, ncof)) for ix in xrange(MAXIT): # Set up design matrix for least squares fit. pow1 = wt bb = pow1 * (fs + abs(mad) * sign(ee)) for jx in xrange(m + 1): u[:, jx] = pow1 pow1 = pow1 * x pow1 = -bb for jx in xrange(m + 1, ncof): pow1 = pow1 * x u[:, jx] = pow1 [u1, w, v] = linalg.svd(u, full_matrices=False) cof = where(w == 0, 0.0, np.dot(bb, u1) / w) cof = np.dot(cof, v) # Tabulate the deviations and revise the weights ee = polyval(cof[m::-1], x) / \ polyval(cof[ncof:m:-1].tolist() + [1, ], x) - fs wt = np.abs(ee) devmax = max(wt) mad = wt.mean() # mean absolute deviation # Save only the best coefficients found if (devmax <= smallest_devmax): smallest_devmax = devmax c1 = cof[m::-1] c2 = cof[ncof:m:-1].tolist() + [1, ] if trace: print('Iteration=%d, max error=%g' % (ix, devmax)) plt.plot(x, fs, x, ee + fs) return poly1d(c1), poly1d(c2) def main(): [c1, c2] = padefitlsq(exp, 3, 3, 0, 2) x = linspace(0, 4) plt.plot(x, polyval(c1, x) / polyval(c2, x), 'g') plt.plot(x, exp(x), 'r') import scipy.special as sp p = [[1, 1, 1], [2, 2, 2]] pi = polyint(p, 1) _pr = polyreloc(p, 2) _pd = polyder(p) _st = poly2str(p) c = poly1d( 1. / sp.gamma( np.r_[ 6 + 1:0:- 1])) # polynomial coeff exponential function [p, q] = padefit(c) x = linspace(0, 4) plt.plot(x, c(x), x, p(x) / q(x), 'g-', x, exp(x), 'r.') plt.close() x = arange(4) dx = dct(x) _idx = idct(dx) a = 0 b = 2 ck = chebfit(exp, 6, a, b) _t = chebval(0, ck, a, b) x = linspace(0, 2, 6) plt.plot(x, exp(x), 'r', x, chebval(x, ck, a, b), 'g.') # x1 = chebroot(9).'*(b-a)/2+(b+a)/2 ; # ck1 =chebfit([x1 exp(x1)],9,a,b); # plot(x,exp(x),'r'), hold on # plot(x,chebval(x,ck1,a,b),'g'), hold off _t = poly2hstr([1, 1, 2]) py = [1, 0] px = polyshift(py, 0, 5) _t1 = polyval(px, [0, 2.5, 5]) # % This is the same as the line below _t2 = polyval(py, [-1, 0, 1]) px = [1, 0] py = polyishift(px, 0, 5) t1 = polyval(px, [0, 2.5, 5]) # % This is the same as the line below t2 = polyval(py, [-1, 0, 1]) print(t1, t2) def test_polydeg(): x = np.linspace(0, 10, 300) y = np.sin(x ** 3 / 100) ** 2 + 0.05 * np.random.randn(x.size) n = polydeg(x, y) # n = 2 p = orthofit(x, y, n) xi = linspace(x.min(), x.max()) ys0 = orthoval(p, x) ys = orthoval(p, xi) ys2 = orthoval(p, xi) plt.plot(x, y, '.', x, ys0, 'k', xi, ys, 'r', xi, ys2, 'r.') p0 = ortho2poly(p) p1 = polyfit(x, ys0, n) plt.plot(xi, polyval(p0, xi), 'g-.', xi, polyval(p1, xi), 'go') plt.show('hold') def test_docstrings(): import doctest print('Testing docstrings in %s' % __file__) doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE) def chebvandernd(deg, *xi): """Pseudo-Vandermonde matrix of given degrees. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample points `(x1, x2, ..., xn)`. If `l, m, n` are the given degrees in `x1, x2, x3`, then The pseudo-Vandermonde matrix is defined by .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = T_i(x1)*T_j(x2)*T_k(x3), where `0 <= i <= l`, `0 <= j <= m`, and `0 <= k <= n`. The leading indices of `V` index the points `(x, y, z)` and the last index encodes the degrees of the Chebyshev polynomials. If ``V = chebvandernd([xdeg, ydeg, zdeg], x, y, z)``, then the columns of `V` correspond to the elements of a 3-D coefficient array `c` of shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},... and ``np.dot(V, c.flat)`` and ``chebvalnd(c, x, y, z)`` will be the same up to roundoff. This equivalence is useful both for least squares fitting and for the evaluation of a large number of N-D Chebyshev series of the same degrees and sample points. Parameters ---------- deg : list of ints List of maximum degrees of the form [x1_deg, x2_deg, ...,xn_deg]. x1, x2, ..., xn : array_like Arrays of point coordinates, all of the same shape. The dtypes will be converted to either float64 or complex128 depending on whether any of the elements are complex. Scalars are converted to 1-D arrays. Returns ------- vander : ndarray The shape of the returned matrix is ``x1.shape + (order,)``, where :math:`order = (deg[0]+1)*(deg([1]+1)*...*(deg[n-1]+1)`. The dtype will be the same as the converted `x1`, `x2`, ... `xn`. See Also -------- chebvander, chebvalnd, chebfitnd """ ideg = [int(d) for d in deg] is_valid = np.array([di == d and di >= 0 for di, d in zip(ideg, deg)]) if np.any(is_valid != 1): raise ValueError("degrees must be non-negative integers") ndim = len(xi) if len(ideg) != ndim: msg = 'length of deg must be the same as number of dimensions' raise ValueError(msg) xi = np.array(xi, copy=0) + 0.0 chebvander = np.polynomial.chebyshev.chebvander shape0 = xi[0].shape s0 = (1,) * ndim vxi = [chebvander(x, d).reshape(shape0 + s0[:i] + (-1,) + s0[i + 1::]) for i, (d, x) in enumerate(zip(ideg, xi))] v = reduce(np.multiply, vxi) return v.reshape(v.shape[:-ndim] + (-1,)) def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): """ Least squares fit of Chebyshev series to N-dimensional data. Return the coefficients of a Chebyshev series of degree `deg` that is the least squares fit to the data values `f` given at points `x1`, `x2`,..., `xn` The fitted polynomial(s) are in the form .. math:: p(x,y) = c_00 + c_11 * T_1(x)*T_1(y) + ..c_ij * T_i(x)*T_j(y). + c_nm * T_n(x)*T_m(y), where `n`, `m` is `deg`. Parameters ---------- xi: tuple x1-, x2-,....xn-coordinates of the sample points. f : array_like function values at the sample points ``(x1[i], x2[i], ..., xn[i])``. deg : list Degrees of the fitting series in the x1, x2, ..., xn directions, respectively. 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 size(x1)*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, optional Weights. If not None, the contribution of each point ``(x1[i], x2[i], ..., xn[i])`` to the fit is weighted by `w[i]`. Ideally the weights are chosen so that the errors of the products ``w[i]*f[i]`` all have the same variance. The default value is None. Returns ------- coef : ndarray, shape (M1, M2,..., Mn) Chebyshev coefficients ordered from low to high. [residuals, rank, singular_values, rcond] : list These values are only returned if `full` = True resid -- sum of squared residuals of the least squares fit rank -- the numerical rank of the scaled Vandermonde matrix sv -- singular values of the scaled Vandermonde matrix rcond -- 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 -------- chebvalnd, chebgridnd Notes ----- The solution is the coefficients of the Chebyshev series `p` that minimizes the sum of the weighted squared errors .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2, where :math:`w_j` are the weights. This problem is solved by setting up as the (typically) overdetermined matrix equation .. math:: V(x, y) * c = w * y, where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the coefficients to be solved for, `w` are the weights, and `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 coefficient 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 Chebyshev 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 -------- """ # xi = np.array(xi, copy=0) + 0.0 z = np.array(f) degrees = np.asarray(deg, dtype=int) orders = degrees + 1 order = np.product(orders) ndims = np.array([x.ndim for x in xi]) ndim = len(ndims) sizes = np.array([x.size for x in xi]) if np.any(ndims != ndim) or z.ndim != ndim: raise TypeError("expected %dD array for x1, x2,...,xn and f" % ndim) if np.any(sizes == 0): raise TypeError("expected non-empty vector for xi") lhs = chebvandernd(degrees, *xi).reshape((-1, order)) rhs = z.ravel() if w is not None: w = np.asarray(w).ravel() + 0.0 if len(lhs) != len(w): raise TypeError("expected x and w to have same length") lhs = lhs * w rhs = rhs * w if rcond is None: rcond = xi[0].size * np.finfo(x.dtype).eps if issubclass(lhs.dtype.type, np.complexfloating): scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(axis=0)) else: scl = np.sqrt(np.square(lhs).sum(axis=0)) scl[scl == 0] = 1 # Solve the least squares problem. c, resids, rank, s = np.linalg.lstsq(lhs/scl, rhs, rcond) c = (c/scl).reshape(orders) 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 chebvalnd(c, *xi): """ Evaluate a N-D Chebyshev series at points (x1, x2, ..., xn). This function returns the values: .. math:: p(x1,x2,...,xn) = \\sum_{i,j,...,k} c_{i,j,...,k} * T_i(x1) * T_j(x2)*...* T_k(xn) The parameters `x1`, `x2`, ...., `xn` are converted to arrays only if they are tuples or a lists, otherwise they are treated as a scalars and they must have the same shape after conversion. In either case, either `x1`, `x2`, ..., `xn` or their elements must support multiplication and addition both with themselves and with the elements of `c`. If `c` has fewer than N dimensions, ones are implicitly appended to its shape to make it N-D. The shape of the result will be c.shape[3:] + x1.shape. Parameters ---------- c : array_like Array of coefficients ordered so that the coefficient of the term of multi-degree i,j,...,k is contained in ``c[i,j,...,k]``. If `c` has dimension greater than N the remaining indices enumerate multiple sets of coefficients. x1, x2,..., xn : array_like, compatible object The N dimensional series is evaluated at the points `(x1, x2,...,xn)`, where `x1`, `x2`,..., `xn` must have the same shape. If any of `x1`, `x2`, ..., `xn` is a list or tuple, it is first converted to an ndarray, otherwise it is left unchanged and if it isn't an ndarray it is treated as a scalar. Returns ------- values : ndarray, compatible object The values of the multidimensional polynomial on points formed with triples of corresponding values from `x`, `y`, and `z`. See Also -------- chebval, chebgridnd, chebfitnd """ try: xi = np.array(xi, copy=0) except: raise ValueError('x, y, z are incompatible') chebval = np.polynomial.chebyshev.chebval c = chebval(xi[0], c) for x in xi[1:]: c = chebval(x, c, tensor=False) return c def chebgridnd(c, *xi): """ Evaluate a N-D Chebyshev series on the Cartesian product of x1, x2,..., xn. This function returns the values: .. math:: p(a,b,...) = \\sum_{i,j,...} c_{i,j,...} * T_i(a) * T_j(b) * ... where the points `(a, b, ...)` consist of all points formed by taking `a` from `x1`, `b` from `x2`, and so on. The resulting points form a grid with `x1` in the first dimension, `x2` in the second, and so on. The parameters `x1`, `x2`, ... and `xn` are converted to arrays only if they are tuples or a lists, otherwise they are treated as a scalars. In either case, either `x1`, `x2`,... and `xn` or their elements must support multiplication and addition both with themselves and with the elements of `c`. If `c` has fewer than N dimensions, ones are implicitly appended to its shape to make it N-D. The shape of the result will be c.shape[3:] + x1.shape + x2.shape + ... + xn.shape Parameters ---------- c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension greater than two the remaining indices enumerate multiple sets of coefficients. x1, x2,..., xn : ndarray, compatible object 1-D arrays representing the coordinates of a grid. The N dimensional series is evaluated at the points in the Cartesian product of `x1`, `x2`, ... and `xn`. If `xi`, is a list or tuple, it is first converted to an ndarray, otherwise it is left unchanged and, if it isn't an ndarray, it is treated as a scalar. Returns ------- values : ndarray, compatible object The values of the N dimensional polynomial at points in the Cartesian product of `x1`, `x2`, ... and `xn`. See Also -------- chebval, chebvalnd, chebfitnd """ chebval = np.polynomial.chebyshev.chebval for x in xi: c = chebval(x, c) return c def test_chebfit1d(): def f(x): return np.exp(-x**2) # x = chebroot(n=64, kind=1) # z = f(x) c = chebfit(f, n=64)[::-1] xi = np.linspace(-1, 1, 151) zi = np.polynomial.chebyshev.chebval(xi, c) # plt.plot(xi, zi,'.', xi, f(xi)) plt.semilogy(xi, np.abs(zi-f(xi))) plt.show('hold') def test_chebfit2d(): n = 3 xorder, yorder = n-1, n-1 x = chebroot(n=n, kind=1) xgrid, ygrid = meshgrid(x, x) def f(x, y): return np.exp(-x**2-6*y**2) zgrid = f(xgrid, ygrid) # v2d = np.polynomial.chebyshev.chebvander2d(xgrid, ygrid, # [xorder,yorder]).reshape((-1, (xorder+1)*(yorder+1))) # coeff, residuals, rank, s = np.linalg.lstsq(v2d, zgrid.ravel()) # doeff = coeff.reshape(xorder+1,yorder+1) _dcoeff2 = chebfitnd((xgrid, ygrid), zgrid, [xorder, yorder]) dcoeff = chebfit_dct(f, n=(xorder+1, yorder+1)) xi = np.linspace(-1, 1, 151) Xi, Yi = np.meshgrid(xi, xi) Zi = f(Xi, Yi) zzi = chebvalnd(dcoeff, Xi, Yi) _devi = Zi - zzi # plot residuals # zz = np.polynomial.chebyshev.chebval2d(xgrid, ygrid, dcoeff) zz = chebvalnd(dcoeff, xgrid, ygrid) dev = zgrid - zz # plt.spy(np.abs(dcoeff)>1e-13) plt.contourf(xgrid, ygrid, np.abs(dev)) # plt.contourf(Xi, Yi, np.abs(devi)) plt.colorbar() # plt.semilogy(np.abs(devi.ravel())) plt.show('hold') if __name__ == '__main__': if False: # True: # main() else: test_chebfit2d() # test_docstrings() # test_polydeg()