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 "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()
>
>
>
>
>
>
>
>
>
>
> ------------------------------------------------------------------------------
> Let Crystal Reports handle the reporting - Free Crystal Reports 2008
> 30-Day
> trial. Simplify your report design, integration and deployment - and focus
> on
> what you do best, core application coding. Discover what's new with
> Crystal Reports now. http://p.sf.net/sfu/bobj-july
> _______________________________________________
> Matplotlib-users mailing list
> Matplotlib-users@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>
>
--
View this message in context:
http://www.nabble.com/%22Piecewise-Cubic-Hermite-Interpolating-Polynomial%22-in-python-tp25204843p25530075.html
Sent from the matplotlib - users mailing list archive at Nabble.com.
------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
