npplus.pwpoly module¶

Piecewise polynomial interpolation and curve fitting.

This module simplifies, replaces, or enhances many scipy.interpolate programming interfaces. The interface consists of two classes, PwPoly and its periodic subclass PerPwPoly, and four functions, pline, spline, plfit, splfit. The pline and plfit functions are special cases of the general spline and splfit functions returning piecewise linear functions; spline and splfit return piecewise cubics by default, but can be used to construct piecewise polynomials of any degree.

The spline and pline functions return the interpolating function passing exactly through a given set of points. The splfit and plfit functions return the least squares best fit to a cloud of given points with given knot points (the x values separating the polynomial pieces). All four functions work for curves in multidimensional space. A quick tour of PwPoly functionality:

f = pline(x, y)  # polyline through points (x, y)
g = spline(x, y)  # natural cubic spline through (x, y)
f = plfit(x, xdat, ydat)  # best fit polyline with knots at x
g = splfit(x, xdat, ydat)  # best fit cubic spline with knots at x
f(xp)  # if f from pline, same as interp(xp, x, y)
g(xp)  # evaluate cubic spline at points xp
h = 2*f + g  # operators +, -, * work for piecewise polynomials
xv = h.roots(value)  # all xv such that h(xv) = value
dydx = h.deriv()  # piecewise poly of one less degree than h
iydx = h.integ()  # piecewise poly of one greater degree than h

The only common use cases for the PwPoly constructor are for making histograms, that is, piecewise constant functions, which you do by supplying one more or less y value than x value, and for constructing piecewise cubic functions passing through (x,y) with given dy/dx:

histo = PwPoly(x, y)  # y.size = x.size+1 or x.size-1
cubic = PwPoly(x, y, dydx)  # x, y, dydx same size

The spline interpolation and fitting routines here use only standard numpy plus the scipy solve_banded function. The scipy.interpolate spline routines are wrappers around the old Fortran fitpack package. A key advantage of the pwpoly module is that you can mine the code here and adapt it to slightly different problems, whereas compiled fitpack code is not adaptable.

class PerPwPoly(xk=None, *args)[source]¶

Bases: PwPoly

Periodic piecewise polynomial function.

Like PwPoly, except first and last knot points xk are the same point, but advanced by one period. Input arguments to the function will be reduced to this period before evaluation, so there is no extrapolation.

Note that when you initialize a PerPwPoly, you must explicitly supply the duplicate final point at the end of the period.

class PwPoly(xk=None, *args)[source]¶

Bases: object

Piecewise polynomial function.

Typical usage:

pwp = PwPoly(xk, yk, dydxk)  # define a piecewise cubic function
y = pwp(x)                   # evaluate function at x
y, dydx, d2ydx2 = pwp(x, 2)  # evaluate function and 2 derivatives

The points xk are called “knot points” of the piecewise polynomial. Outside the endpoints of xk, the function will have the degree to which it is continuous at the other endpoints.

Parameters

xk (array_like) – List of abcissa values bounding the pieces, in increasing order. As a convenience, strictly decreasing xk are also accepted.
yk (array_like) –
Each is a list of function and derivative values yk, dydxk, d2ydx2k/2, etc. corresponding to the points xk. This produces an odd degree polynomial in each interval, continuous to the specified derivatives at the xk points.

Note that when you specify derivatives higher than first, you need to divide the n-th derivative by n!. That is, you are actually furnishing the Taylor series coefficients, not the derivatives.

To get even degree polynomials, specify the final derivative with one fewer point than xk – this produces a polynomial of degree P=2*N continuous at xk up to the N-1 derivatives you specified at all xk points, and with the P-th (not N-th) derivative/p! equal to the final values specified.

Finally, for complete generality, you can specify every argument with one more point than xk. In this case, the arguments simply become the coefficients of the polynomials in each interval with no continuity guaranteed at the xk. Note that the origin for coefficients in each interval is the first point of the interval (the smaller), except for the semi-infinite interval before the first knot point, which is relative to the first knot.

All the yk, dydxk, etc. may have a set of leading axes to define a multidimensional piecewise polynomial function of xk. That is, the interpolation direction is the last axis of xk, yk, dydxk, etc.
dydxk (array_like) –
Each is a list of function and derivative values yk, dydxk, d2ydx2k/2, etc. corresponding to the points xk. This produces an odd degree polynomial in each interval, continuous to the specified derivatives at the xk points.

Note that when you specify derivatives higher than first, you need to divide the n-th derivative by n!. That is, you are actually furnishing the Taylor series coefficients, not the derivatives.

To get even degree polynomials, specify the final derivative with one fewer point than xk – this produces a polynomial of degree P=2*N continuous at xk up to the N-1 derivatives you specified at all xk points, and with the P-th (not N-th) derivative/p! equal to the final values specified.

Finally, for complete generality, you can specify every argument with one more point than xk. In this case, the arguments simply become the coefficients of the polynomials in each interval with no continuity guaranteed at the xk. Note that the origin for coefficients in each interval is the first point of the interval (the smaller), except for the semi-infinite interval before the first knot point, which is relative to the first knot.

All the yk, dydxk, etc. may have a set of leading axes to define a multidimensional piecewise polynomial function of xk. That is, the interpolation direction is the last axis of xk, yk, dydxk, etc.
d2ydx2k/2 (array_like) –
Each is a list of function and derivative values yk, dydxk, d2ydx2k/2, etc. corresponding to the points xk. This produces an odd degree polynomial in each interval, continuous to the specified derivatives at the xk points.

Note that when you specify derivatives higher than first, you need to divide the n-th derivative by n!. That is, you are actually furnishing the Taylor series coefficients, not the derivatives.

To get even degree polynomials, specify the final derivative with one fewer point than xk – this produces a polynomial of degree P=2*N continuous at xk up to the N-1 derivatives you specified at all xk points, and with the P-th (not N-th) derivative/p! equal to the final values specified.

Finally, for complete generality, you can specify every argument with one more point than xk. In this case, the arguments simply become the coefficients of the polynomials in each interval with no continuity guaranteed at the xk. Note that the origin for coefficients in each interval is the first point of the interval (the smaller), except for the semi-infinite interval before the first knot point, which is relative to the first knot.

All the yk, dydxk, etc. may have a set of leading axes to define a multidimensional piecewise polynomial function of xk. That is, the interpolation direction is the last axis of xk, yk, dydxk, etc.
... (array_like) –
Each is a list of function and derivative values yk, dydxk, d2ydx2k/2, etc. corresponding to the points xk. This produces an odd degree polynomial in each interval, continuous to the specified derivatives at the xk points.

Note that when you specify derivatives higher than first, you need to divide the n-th derivative by n!. That is, you are actually furnishing the Taylor series coefficients, not the derivatives.

To get even degree polynomials, specify the final derivative with one fewer point than xk – this produces a polynomial of degree P=2*N continuous at xk up to the N-1 derivatives you specified at all xk points, and with the P-th (not N-th) derivative/p! equal to the final values specified.

Finally, for complete generality, you can specify every argument with one more point than xk. In this case, the arguments simply become the coefficients of the polynomials in each interval with no continuity guaranteed at the xk. Note that the origin for coefficients in each interval is the first point of the interval (the smaller), except for the semi-infinite interval before the first knot point, which is relative to the first knot.

All the yk, dydxk, etc. may have a set of leading axes to define a multidimensional piecewise polynomial function of xk. That is, the interpolation direction is the last axis of xk, yk, dydxk, etc.

xk¶

1D array of knot points in strictly increasing order.

Type: ndarray

c¶

Coefficient array. Axes are (degree, dimensionality, knots+1) The polynomial coefficients are for the polynomial in (x-xk[i-1]), where xk[i-1] <= x < xk[i], except for c[...,0], which like c[...,1] is for the polynomial in (x-xk[0]). Coefficients are in order of increasing powers of x.

Type: ndarray

See also

spline: PwPoly interpolation by cubic and other degree splines
splfit: PwPoly fitting by cubic and other degree splines
pline: polyline interpolation
plfit: polyline fitting

Notes

Let p and q be PwPoly instances. The following operators work:

p(x): Evaluate p at x, returning array of same shape as x. If p is multidimensional, its additional dimensions are leading dimensions of the result.
p+q, p-q, -p, p*q: Return a new PwPoly instance. If p and q have different xk, result the union of the two xk. Either p or q may be scalars (or arrays same shape as p(0)). For array_like q, division p/q also works.
p[i]: Components of p, only for multidimensional p.
len(p): First dimension length of p, only for multidimensional p.

The PwPoly constructor produces a smooth fit that is local, in the sense that the function in the interval between consecutive knot points depends only on the given function and derivative values at the interval endpoints.

Use the spline function to construct a PwPoly which is smoother (for a given degree) by using the function values you provide at all the knot points to determine the function within each interval. The splfit function constructs a PwPoly that is the statistical best fit to a cloud of data points you provide.

The pline and plfit functions are piecewise linear variants of spline and splfit.

PwPoly is intended for low degree polynomials, usually degree 1 and 3. It should work tolerably well at or below degree 7, but roundoff errors will become significant at higher degree. To get a better fit with PwPoly, use more knots, not higher degree.

PwPoly(x, y) where y has one more or one less element than x produces a histogram, that is, a piecewise polynomial of degree zero. This is probably the only useful application of the PwPoly constructor to produce a PwPoly of even degree. For other even degree piecewise polynomials, you probably want to use spline or splfit.

addknots(xk, _nocheck=False)[source]¶

Return a new PwPoly with knot points xk that matches this one.

The xk must be a superset of the existing knot points for this to work as intended.

Parameters

xk (array_like) – 1D array of new knot points.
n (int, optional) – Degree for new piecewise polynomial, default same as this one.

Returns

Like this PwPoly but with knots xk.

Return type

PwPoly

allknots(*args, **kwargs)[source]¶

Return union of knot points among piecewise polynomials.

Parameters

pwp1 (PwPoly) – The input PwPoly instances.
pwp2 (PwPoly) – The input PwPoly instances.
... (PwPoly) – The input PwPoly instances.
tol (float, optional keyword) – Relative tolerance compared to interval between knots, default 1e-6.

Returns

xk – 1D union of knot points in increasing order.

Return type

ndarray

property deg¶: Highest power of the piecewise polynomial as a property.

degree()[source]¶: Highest power of the piecewise polynomial as a method.

deriv(m=1)[source]¶

Derivative as a new piecewise polynomial.

Parameters: m (int, optional) – The number of derivatives to perform, default 1.
Returns: The derivative of the PwPoly.
Return type: PwPoly

integ(m=1, k=(), lbnd=0)[source]¶

Integral as a new piecewise polynomial.

Parameters

m (int, optional) – The number of integrations to perform, default 1.
k (array_like, optional) – Integration constants. The first constant is applied to the first integration, the second to the second, and so on. The list of values must less than or equal to m in length and any missing values are set to zero.
lbnd (float, optional) – The lower bound of the definite integral, default 0. The integration constant k is the value of the integral at x = lbnd, where x is the independent variable of the piecewise polynomial.

Returns

The integral of the PwPoly.

Return type

PwPoly

jumps(n=None)[source]¶

Return jumps in polynomial coefficients at the knot points.

Parameters: n (int, optional) – Maximum degree of returned coefficients, all by default.
Returns: dc – Shape is (n, ..., K) where K is number of knot points. Degree is first axis, dc[0] is jump in function value, dc[1] is jump in derivative, dc[2] is half the jump in second derivative, and so on. Knot is final axis. Any additional axes are the dimensionality of the PwPoly.
Return type: ndarray

property ndim¶: Number of dimensions of result when evaluated at a scalar x.

classmethod new(xk, c)[source]¶

Make a new PwPoly with given knots and coefficients.

Parameters

xk (ndarray, shape (K) or (K+1)) –
c (ndarray, shape (N, ..., K+1)) –

Returns

p – The new PwPoly.

Return type

PwPoly or a subclass

Notes

If xk.size is K+1, xk will not be copied. If xk.size is K, xk[0] will be duplicated. The coefficient array c is always used uncopied. No error checking is done; this is a low level method for creating PwPoly results, not a user interface.

reknot(xk, n=None)[source]¶

Return a new PwPoly with knot points xk that approximates this one.

If n is the polynomial degree, the resulting piecewise polynomial will match (n-1)/2 derivatives of the old one at the points xk. If n is even, the n-th derivative at the midpoint of each xk interval will also match the original.

Parameters

xk (array_like) – 1D array of new knot points.
n (int, optional) – Degree for new piecewise polynomial, default same as this one.

Returns

Like this PwPoly but with knots xk.

Return type

PwPoly

roots(value=0.0, tol=1e-09)[source]¶

Return the values of x for which the piecewise polynomial is zero.

Only works for scalar PwPoly.

Parameters: value (float, optional) – Return values of x for which piecewise polynomial equals value. Defaults to 0.0.
Returns: 1D array of x values where this PwPoly evaluates to value.
Return type: ndarray

property shape¶: Shape of result when evaluated at a scalar x.

plfit(xk, x, y, sigy=None, lo=(), hi=(), per=False, extrap=None, solver=<function _plfitter>, **kwargs)[source]¶

Return best fit linear PwPoly with knots xk to given points (x,y).

Gives the same fit as splfit(xk, x, y, n=1, nc=0), that is, the best fit piecewise polynomial through (x,y) with knots xk. See splfit documentation for detailed description of the common parameters. If you need the cost of continuity constraints, use splfit instead.

See also

splfit: general PwPoly spline fitting function

pline(x, y, extrap=None, per=False)[source]¶

Return a piecewise linear function through given points (x,y).

Convenience shorthand for spline(x, y, n=1), intended as a more general version of the interp function with a slightly different interface:

pline(x, y)(xp)  <--->  interp(xp, x, y)

But pline permits y to have additional dimensions and x can be a decreasing series.

Note that pline(x, y) by default is constant before the first x and after the final x. Use pline(x, y, 1) to extrapolate using the linear function in the first and last intervals.

See also

spline: general PwPoly spline interpolation function

splfit(xk, x, y, sigy=None, n=3, nc=None, lo=(), hi=(), per=False, extrap=None, cost=None)[source]¶

Return the best fit PwPoly with knots xk to given points (x,y).

Parameters

xk (1D array_like) – Knot points of the PwPoly, sparse enough that there are some x values “near” each interval of xk. The xk must be strictly monotonic.
x (array_like) –
Cloud of points to fit. If f is the returned PwPoly,

chi2 = sum(((f(x) - y)/sigy)**2)

will be minimum for the given PwPoly knot points xk, degree n, and boundary constraints. The y array must have the same or conformable trailing dimensions as x. However, y may have additional leading dimensions to produce a multi-dimensional PwPoly curve. In this case, the chi2 sum extends over all the leading dimensions for each point. The order of the (x,y) points is irrelevant; indeed x need not be a 1D array.
y (array_like) –
Cloud of points to fit. If f is the returned PwPoly,

chi2 = sum(((f(x) - y)/sigy)**2)

will be minimum for the given PwPoly knot points xk, degree n, and boundary constraints. The y array must have the same or conformable trailing dimensions as x. However, y may have additional leading dimensions to produce a multi-dimensional PwPoly curve. In this case, the chi2 sum extends over all the leading dimensions for each point. The order of the (x,y) points is irrelevant; indeed x need not be a 1D array.
sigy (array_like, optional) – Standard deviations of y, must be conformable to y, default 1.
n (int, optional) – The degree of the PwPoly, default 3 (piecewise cubic).
nc (int, optional) – The degree to which the result PwPoly is continuous, default n-1 (that is, the PwPoly is a spline of degree n). Note that a smaller value of nc requires more (x,y) points per xk interval.
lo (tuple of array_like, optional) – Boundary conditions at the endpoints xk[0] and xk[-1]. Each tuple represents the values of (y, dydx, d2ydx2/2, d3ydx3/6, ...) at xk[0] for lo or xk[-1] for hi. Use None for a derivative to be not specified. See the docstring for spline for details; however, note that the lo and hi tuples for splfit begin with the function value (zero-th derivative), while the lo and hi tuples for spline begin with the first derivative.
hi (tuple of array_like, optional) – Boundary conditions at the endpoints xk[0] and xk[-1]. Each tuple represents the values of (y, dydx, d2ydx2/2, d3ydx3/6, ...) at xk[0] for lo or xk[-1] for hi. Use None for a derivative to be not specified. See the docstring for spline for details; however, note that the lo and hi tuples for splfit begin with the function value (zero-th derivative), while the lo and hi tuples for spline begin with the first derivative.
per (bool) – True to make all derivatives match at first and last points of xk. This produces a PerPwPoly function, which maps any inputs outside the first and last points of xk into that interval.
extrap (int or (int,int)) – Degree to maintain continuity beyond endpoints, can be set separately for before first and after last interval.
cost (bool, optional) – Set True to return cost array continuity of various degrees.

Returns

fit (PwPoly) – The best fit to the given (x,y) with the given knots xk.
cost (tuple of ndarray) – This return is only present with cost=1 parameter. The Lagrange multipliers used to enforce the chi2 minimization constraints associated with the continuous function values and nc derivative values, plus any additional boundary condition constraints you specified with the lo or hi keywords. The interior knot point constraints are cost[0] for function continuity, cost[1] for first derivative continuity, through cost[nc] for nc-derivative continuity. Then cost[nc+1:nc+1+nlo] are the costs of the lo constraints and cost[nc+1+nlo:] are the costs of the hi constraints, where nlo=len(lo).

The shape of each of the cost[:nc+1] interior constraints is y[...,1:-1].shape, that is, any leading dimensions of y followed by the number of interior knot points where the constraints apply. The shape of each of the cost[nc+1:] boundary constraints is y[...,0].shape, that is, any leading dimensions of y.

See also

plfit: special case of splfit for n=1, nc=0
spline: general PwPoly spline interpolation function

Notes

The cost values require a bit more explanation. Each interior constraint is schematically c[this] - l2r.dot(c[prev]) = 0 where l2r is the linear transformation mapping the polynomial coefficients c[prev] in the interval to the left to their values at xk[this] knot. That is, each interior constraint is that the jump in the value of the polynomial coefficient of some degree be zero. The cost is defined to be the partial derivative of chi2 with respect to that jump. In other words, if we relaxed just that one constraint and allowed a jump of dc, chi2 would change by cost*dc. We keep this sign convention for the lo and hi constraints as well – a positive cost means that if the coefficient changes by a positive amount for increasing x across the knot, then chi2 will increase.

spline(x, y, n=3, lo=(), hi=(), per=False, extrap=None)[source]¶

Construct a PwPoly as a spline through given points (x,y).

A spline is a piecewise polynomial of degree n (default 3) passing through all of the given points with continuous derivatives up to degree n-1 at every point. Before x[0] and after x[-1], the result will be of degree n-1, maintaining continuity of n-1 derivatives at every knot point, including the first and last.

Parameters

x (array_like) – 1D array of knot points, strictly monotonic.
y (array_like) – The given (x,y) data points to be fit. The y data may have additional leading dimensions, in order to fit a curve in a multi-dimensional space.
n (int, optional) – The degree of the PwPoly. Defaults to 3, a piecewise cubic spline.
lo (tuple of array_like, optional) –
Boundary conditions at the endpoints x[0] and x[-1]. Each tuple represents the values of (dydx, d2ydx2/2, d3ydx3/6, ...) at x[0] for lo or x[-1] for hi. Use None for a derivative to be not specified. The derivatives are broadcast to any leading dimensions of y if necessary.

The highest derivative you can specify is n-1, that is d2ydx2/2 for a cubic spline, dydx for a quadratic spline, and so on. For the special case that you only wish to specify dydx at an endpoint, lo or hi need not be a tuple (dydx,); spline accepts simply dydx to mean (dydx, None, None, ...). The maximum number of derivatives you can specify in both lo and hi is n-1, that is, 2 for a cubic spline, 4 for a pentic, and so on. If you specify fewer than n-1 derivatives, spline will set the highest unspecified derivatives to zero, begining with lo and alternating lo, hi, lo, etc., until n-1 derivative are specified. For example, if you specify neither lo nor hi for a cubic spline, then spline sets d2ydx2=0 at both endpoints, and the natural spline with zero curvature at the endpoints is the result. This is equivalent to lo=(None,0.), hi=(None,0.).
hi (tuple of array_like, optional) –
Boundary conditions at the endpoints x[0] and x[-1]. Each tuple represents the values of (dydx, d2ydx2/2, d3ydx3/6, ...) at x[0] for lo or x[-1] for hi. Use None for a derivative to be not specified. The derivatives are broadcast to any leading dimensions of y if necessary.

The highest derivative you can specify is n-1, that is d2ydx2/2 for a cubic spline, dydx for a quadratic spline, and so on. For the special case that you only wish to specify dydx at an endpoint, lo or hi need not be a tuple (dydx,); spline accepts simply dydx to mean (dydx, None, None, ...). The maximum number of derivatives you can specify in both lo and hi is n-1, that is, 2 for a cubic spline, 4 for a pentic, and so on. If you specify fewer than n-1 derivatives, spline will set the highest unspecified derivatives to zero, begining with lo and alternating lo, hi, lo, etc., until n-1 derivative are specified. For example, if you specify neither lo nor hi for a cubic spline, then spline sets d2ydx2=0 at both endpoints, and the natural spline with zero curvature at the endpoints is the result. This is equivalent to lo=(None,0.), hi=(None,0.).
per (bool) – True to make all derivatives match at first and last points of x. This produces a PerPwPoly function, which maps any inputs outside the first and last points of x into that interval.
extrap (int or (int, int)) – Degree to maintain continuity beyond endpoints, can be set separately for before first and after last interval.

See also

pline: polyline constructor, use instead of spline for linear case
splfit: construct a pline by fitting points rather than interpolating

Notes

The K knots define K-1 intervals, requiring a total of (n+1)*(K-1)=n*K-n+K-1 polynomial coefficients. There are n*(K-2) continuity constraints at the interior knots and K y values, for a total of n*K-2*n+K constraint equations. We require n-1 additional constraints to determine a unique spline fit. When n is odd (linear, cubic, pentic, etc.), these additional constraints may be divided equally between the first and last knots points. When n is even, there is no way to treat the first and last knots symmetrically.

By default, spline forces the (n-1)//2 highest derivatives to zero at each endpoint, with the first endpoint getting an extra zero derivative in the case of even n. The lo and hi keywords allow you to split up the n-1 constraints among the possible derivatives and endpoints however you wish.

npplus

Navigation

Related Topics

npplus.pwpoly module¶