Menu
▾
▴
[Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
|
From: Chris M. <chr...@gm...> - 2009-08-28 17:43:31
|
Recently, I had a need for a monotonic piece-wise cubic Hermite
interpolator. Matlab provides the function "pchip" (Piecewise Cubic
Hermite Interpolator), but when I Googled I didn't find any Python
equivalent. I tried "interp1d()" from scipy.interpolate but this was
a standard cubic spline using all of the data - not a piece-wise cubic
spline.
I had access to Matlab documentation, so I spent a some time tracing
through the code to figure out how I might write a Python duplicate.
This was an massive exercise in frustration and a potent reminder on
why I love Python and use Matlab only under duress. I find typical
Matlab code is poorly documented (if at all) and that apparently
includes the code included in their official releases. I also find
Matlab syntax “dated” and the code very difficult to “read”.
Wikipedia to the rescue.
Not to be deterred, I found a couple of very well written Wikipedia
entries, which explained in simple language how to compute the
interpolant. Hats off to whoever wrote these entries – they are
excellent. The result was a surprising small amount of code
considering the Matlab code was approaching 10 pages of
incomprehensible code. Again - strong evidence that things are just
better in Python...
Offered for those who might have the same need – a Python pchip()
equivalent ==> pypchip(). Since I'm not sure how attachments work (or
if they work at all...), I copied the code I used below, followed by a
PNG showing "success":
#
# pychip.py
# Michalski
# 20090818
#
# Piecewise cubic Hermite interpolation (monotonic...) in Python
#
# References:
#
# Wikipedia: Monotone cubic interpolation
# Cubic Hermite spline
#
# A cubic Hermte spline is a third degree spline with each
polynomial of the spline
# in Hermite form. The Hermite form consists of two control points
and two control
# tangents for each polynomial. Each interpolation is performed on
one sub-interval
# at a time (piece-wise). A monotone cubic interpolation is a
variant of cubic
# interpolation that preserves monotonicity of the data to be
interpolated (in other
# words, it controls overshoot). Monotonicity is preserved by
linear interpolation
# but not by cubic interpolation.
#
# Use:
#
# There are two separate calls, the first call, pchip_init(),
computes the slopes that
# the interpolator needs. If there are a large number of points to
compute,
# it is more efficient to compute the slopes once, rather than for
every point
# being evaluated. The second call, pchip_eval(), takes the slopes
computed by
# pchip_init() along with X, Y, and a vector of desired "xnew"s and
computes a vector
# of "ynew"s. If only a handful of points is needed, pchip() is a
third function
# which combines a call to pchip_init() followed by pchip_eval().
#
import pylab as P
#=========================================================
def pchip(x, y, xnew):
# Compute the slopes used by the piecewise cubic Hermite
interpolator
m = pchip_init(x, y)
# Use these slopes (along with the Hermite basis function) to
interpolate
ynew = pchip_eval(x, y, xnew)
return ynew
#=========================================================
def x_is_okay(x,xvec):
# Make sure "x" and "xvec" satisfy the conditions for
# running the pchip interpolator
n = len(x)
m = len(xvec)
# Make sure "x" is in sorted order (brute force, but works...)
xx = x.copy()
xx.sort()
total_matches = (xx == x).sum()
if total_matches != n:
print "*" * 50
print "x_is_okay()"
print "x values weren't in sorted order --- aborting"
return False
# Make sure 'x' doesn't have any repeated values
delta = x[1:] - x[:-1]
if (delta == 0.0).any():
print "*" * 50
print "x_is_okay()"
print "x values weren't monotonic--- aborting"
return False
# Check for in-range xvec values (beyond upper edge)
check = xvec > x[-1]
if check.any():
print "*" * 50
print "x_is_okay()"
print "Certain 'xvec' values are beyond the upper end of 'x'"
print "x_max = ", x[-1]
indices = P.compress(check, range(m))
print "out-of-range xvec's = ", xvec[indices]
print "out-of-range xvec indices = ", indices
return False
# Second - check for in-range xvec values (beyond lower edge)
check = xvec< x[0]
if check.any():
print "*" * 50
print "x_is_okay()"
print "Certain 'xvec' values are beyond the lower end of 'x'"
print "x_min = ", x[0]
indices = P.compress(check, range(m))
print "out-of-range xvec's = ", xvec[indices]
print "out-of-range xvec indices = ", indices
return False
return True
#=========================================================
def pchip_eval(x, y, m, xvec):
# Evaluate the piecewise cubic Hermite interpolant with
monoticity preserved
#
# x = array containing the x-data
# y = array containing the y-data
# m = slopes at each (x,y) point [computed to preserve
monotonicity]
# xnew = new "x" value where the interpolation is desired
#
# x must be sorted low to high... (no repeats)
# y can have repeated values
#
# This works with either a scalar or vector of "xvec"
n = len(x)
mm = len(xvec)
############################
# Make sure there aren't problems with the input data
############################
if not x_is_okay(x, xvec):
print "pchip_eval2() - ill formed 'x' vector!!!!!!!!!!!!!"
# Cause a hard crash...
STOP_pchip_eval2
# Find the indices "k" such that x[k] < xvec < x[k+1]
# Create "copies" of "x" as rows in a mxn 2-dimensional vector
xx = P.resize(x,(mm,n)).transpose()
xxx = xx > xvec
# Compute column by column differences
z = xxx[:-1,:] - xxx[1:,:]
# Collapse over rows...
k = z.argmax(axis=0)
# Create the Hermite coefficients
h = x[k+1] - x[k]
t = (xvec - x[k]) / h[k]
# Hermite basis functions
h00 = (2 * t**3) - (3 * t**2) + 1
h10 = t**3 - (2 * t**2) + t
h01 = (-2* t**3) + (3 * t**2)
h11 = t**3 - t**2
# Compute the interpolated value of "y"
ynew = h00*y[k] + h10*h*m[k] + h01*y[k+1] + h11*h*m[k+1]
return ynew
#=========================================================
def pchip_init(x,y):
# Evaluate the piecewise cubic Hermite interpolant with
monoticity preserved
#
# x = array containing the x-data
# y = array containing the y-data
#
# x must be sorted low to high... (no repeats)
# y can have repeated values
#
# x input conditioning is assumed but not checked
n = len(x)
# Compute the slopes of the secant lines between successive points
delta = (y[1:] - y[:-1]) / (x[1:] - x[:-1])
# Initialize the tangents at every points as the average of the
secants
m = P.zeros(n, dtype='d')
# At the endpoints - use one-sided differences
m[0] = delta[0]
m[n-1] = delta[-1]
# In the middle - use the average of the secants
m[1:-1] = (delta[:-1] + delta[1:]) / 2.0
# Special case: intervals where y[k] == y[k+1]
# Setting these slopes to zero guarantees the spline connecting
# these points will be flat which preserves monotonicity
indices_to_fix = P.compress((delta == 0.0), range(n))
# print "zero slope indices to fix = ", indices_to_fix
for ii in indices_to_fix:
m[ii] = 0.0
m[ii+1] = 0.0
alpha = m[:-1]/delta
beta = m[1:]/delta
dist = alpha**2 + beta**2
tau = 3.0 / P.sqrt(dist)
# To prevent overshoot or undershoot, restrict the position vector
# (alpha, beta) to a circle of radius 3. If (alpha**2 +
beta**2)>9,
# then set m[k] = tau[k]alpha[k]delta[k] and m[k+1] =
tau[k]beta[b]delta[k]
# where tau = 3/sqrt(alpha**2 + beta**2).
# Find the indices that need adjustment
over = (dist > 9.0)
indices_to_fix = P.compress(over, range(n))
# print "overshoot indices to fix... = ", indices_to_fix
for ii in indices_to_fix:
m[ii] = tau[ii] * alpha[ii] * delta[ii]
m[ii+1] = tau[ii] * beta[ii] * delta[ii]
return m
#=
=======================================================================
def CubicHermiteSpline(x, y, x_new):
# Piecewise Cubic Hermite Interpolation using Catmull-Rom
# method for computing the slopes.
#
# Note - this only works if delta-x is uniform?
# Find the two points which "bracket" "x_new"
found_it = False
for ii in range(len(x)-1):
if (x[ii] <= x_new) and (x[ii+1] > x_new):
found_it = True
break
if not found_it:
print
print "requested x=<%f> outside X range[%f,%f]" % (x_new,
x[0], x[-1])
STOP_CubicHermiteSpline()
# Starting and ending data points
x0 = x[ii]
x1 = x[ii+1]
y0 = y[ii]
y1 = y[ii+1]
# Starting and ending tangents (using Catmull-Rom spline method)
# Handle special cases (hit one of the endpoints...)
if ii == 0:
# Hit lower endpoint
m0 = (y[1] - y[0])
m1 = (y[2] - y[0]) / 2.0
elif ii == (len(x) - 2):
# Hit upper endpoints
m0 = (y[ii+1] - y[ii-1]) / 2.0
m1 = (y[ii+1] - y[ii])
else:
# Inside the field...
m0 = (y[ii+1] - y[ii-1])/ 2.0
m1 = (y[ii+2] - y[ii]) / 2.0
# Normalize to x_new to [0,1] interval
h = (x1 - x0)
t = (x_new - x0) / h
# Compute the four Hermite basis functions
h00 = ( 2.0 * t**3) - (3.0 * t**2) + 1.0
h10 = ( 1.0 * t**3) - (2.0 * t**2) + t
h01 = (-2.0 * t**3) + (3.0 * t**2)
h11 = ( 1.0 * t**3) - (1.0 * t**2)
h = 1
y_new = (h00 * y0) + (h10 * h * m0) + (h01 * y1) + (h11 * h * m1)
return y_new
#==============================================================
def main():
############################################################
# Sine wave test
############################################################
# Create a example vector containing a sine wave.
x = P.arange(30.0)/10.
y = P.sin(x)
# Interpolate the data above to the grid defined by "xvec"
xvec = P.arange(250.)/100.
# Initialize the interpolator slopes
m = pchip_init(x,y)
# Call the monotonic piece-wise Hermite cubic interpolator
yvec2 = pchip_eval(x, y, m, xvec)
P.figure(1)
P.plot(x,y, 'ro')
P.title("pchip() Sin test code")
# Plot the interpolated points
P.plot(xvec, yvec2, 'b')
#########################################################################
# Step function test...
#########################################################################
P.figure(2)
P.title("pchip() step function test")
# Create a step function (will demonstrate monotonicity)
x = P.arange(7.0) - 3.0
y = P.array([-1.0, -1,-1,0,1,1,1])
# Interpolate using monotonic piecewise Hermite cubic spline
xvec = P.arange(599.)/100. - 3.0
# Create the pchip slopes slopes
m = pchip_init(x,y)
# Interpolate...
yvec = pchip_eval(x, y, m, xvec)
# Call the Scipy cubic spline interpolator
from scipy.interpolate import interpolate
function = interpolate.interp1d(x, y, kind='cubic')
yvec2 = function(xvec)
# Non-montonic cubic Hermite spline interpolator using
# Catmul-Rom method for computing slopes...
yvec3 = []
for xx in xvec:
yvec3.append(CubicHermiteSpline(x,y,xx))
yvec3 = P.array(yvec3)
# Plot the results
P.plot(x, y, 'ro')
P.plot(xvec, yvec, 'b')
P.plot(xvec, yvec2, 'k')
P.plot(xvec, yvec3, 'g')
P.xlabel("X")
P.ylabel("Y")
P.title("Comparing pypchip() vs. Scipy interp1d() vs. non-
monotonic CHS")
legends = ["Data", "pypchip()", "interp1d","CHS"]
P.legend(legends, loc="upper left")
P.show()
###################################################################
main()
|