Chris,
regarding the print statements in the function you may want to use something
like:
total_matches = (xx == x).sum()
if total_matches != n:
raise ValueError(x values weren't in sorted order)
I also have one question:
in function pchip_eval instead of t = (xvec - x[k]) / h[k] should not
be t = (xvec - x[k]) / h ?
Best regards,
Ricardo Bessa
Chris Michalski-3 wrote:
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 xnews and
computes a vector
#of ynews. 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 -
