Re: [Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
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.compres
Re: [Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
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.bar...@noaa.gov -- 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
Re: [Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
Thanks for the inputs... perhaps it will provide the impetus for future postings as well... chris On Aug 29, 2009, at 11:49 AM, John Hunter wrote: > On Sat, Aug 29, 2009 at 1:12 PM, Eric Firing > 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. 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. > > JDH -- 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
Re: [Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
On Sat, Aug 29, 2009 at 1:12 PM, Eric Firing 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. 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 -- 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
Re: [Matplotlib-users] "Piecewise Cubic Hermite Interpolating Polynomial" in python
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 -- 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
