"Piecewise Cubic Hermite Interpolating Polynomial" in python

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_example.png

···

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 * t3) - (3 * t2) + 1

h10 = t3 - (2 * t2) + t

h01 = (-2* t3) + (3 * t2)

h11 = t3 - t2

# Compute the interpolated value of "y"

ynew = h00y[k] + h10hm[k] + h01y[k+1] + h11hm[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 = alpha2 + beta2

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 * t3) - (3.0 * t2) + 1.0

h10 = ( 1.0 * t3) - (2.0 * t2) + t

h01 = (-2.0 * t3) + (3.0 * t2)

h11 = ( 1.0 * t3) - (1.0 * t2)

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()

Chris Michalski wrote:

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":

Chris,

This looks interesting. I successfully ran your program by using copy and paste to get it into a file, but for the future I certainly recommend that you attach such a file directly--file attachments generally work very well these days, but bad things can happen to code included as inline text.

It would also be helpful if you include the actual URLs of the Wikipedia articles that you used. (Or did I miss them?)

Thanks for providing this code, example, and references.

Eric

The other thing I recommend is do not use the pylab namespace for any
of the numerics. pylab is getting all the numerical functions from
numpy, so if you

  import numpy as np

and then refer to any numerical functions you need as np.somefunc.

Finally, for the functions to be suitable for inclusion in a
production package like numpy or matplotlib.mlab, you should not use
any print statements in the function, but rather a combination of
warnings.warn or exceptions or if it for matplotlib, use the
verbose.report infrastructure. That way users can configure how much
verbosity they want, where the output should be directed, etc.

After a cleanup, you may want to check with numpy or scipy to see of
it could find a home there. There was a discussion at scipy on the
need to improve scipy.interpolate and this seems to go part of the way
toward that objective. So I would start there.

JDH

···

On Sat, Aug 29, 2009 at 1:12 PM, Eric Firing<efiring@...202...> wrote:

This looks interesting. I successfully ran your program by using copy
and paste to get it into a file, but for the future I certainly
recommend that you attach such a file directly--file attachments
generally work very well these days, but bad things can happen to code
included as inline text.

Thanks for the inputs... perhaps it will provide the impetus for future postings as well...

chris

This looks interesting. I successfully ran your program by using copy
and paste to get it into a file, but for the future I certainly
recommend that you attach such a file directly--file attachments
generally work very well these days, but bad things can happen to code
included as inline text.

I haven't contributed to matplotlib or numpy even though I've used them for some years now, so I wasn't sure about the "etiquette" of file attachments.

The other thing I recommend is do not use the pylab namespace for any

of the numerics. pylab is getting all the numerical functions from
numpy, so if you

import numpy as np

and then refer to any numerical functions you need as np.somefunc.

Point well taken. Since pylab exposes most of the numpy calls I use, I typically include pylab instead for nump.

Finally, for the functions to be suitable for inclusion in a
production package like numpy or matplotlib.mlab, you should not use
any print statements in the function, but rather a combination of
warnings.warn or exceptions or if it for matplotlib, use the
verbose.report infrastructure. That way users can configure how much
verbosity they want, where the output should be directed, etc.

Point also well taken. I figured out when there were problems, but even after 7 years of writing large Python package, I haven't found the best way to handle exceptions. Usually I purposely cause a "crash" so I don't miss the fact that the code had ill formed data.

After a cleanup, you may want to check with numpy or scipy to see of
it could find a home there. There was a discussion at scipy on the
need to improve scipy.interpolate and this seems to go part of the way
toward that objective. So I would start there.

I'll send it along to the scipy people. I figured since I figured out a relatively simple solution to a problem that is often encountered, it might find use even in its primitive form. I'll add the URLs to the WIkipedia references as well.

···

On Aug 29, 2009, at 11:49 AM, John Hunter wrote:

On Sat, Aug 29, 2009 at 1:12 PM, Eric Firing<efiring@...202...> > wrote:

JDH

Chris Michalski wrote:

Thanks for the inputs... perhaps it will provide the impetus for future postings as well...

I think this would be a great addition to scipy.interpolate. I encourage you to massage your code to fit the API and scipy standards and contribute it.

-Chris

···

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...259...

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/"Piecewise-Cubic-Hermite-Interpolating-Polynomial"-in-python-tp25204843p25530075.html
Sent from the matplotlib - users mailing list archive at Nabble.com.