Cubic interpolation

Interpolations.jl: The API may be clever, but it isn’t easy to use. How do I interpolate a set of points with a cubic curve?

I tried this:

 npts = size(xyz,1);
    s = fill(0.0, npts);
    for j in 2:npts
        s[j] = s[j-1]+norm(xyz[j,:]-xyz[j-1,:]);
    end 
    CubicSplineInterpolation((s,), xyz[:, 1])

so that the curve would be parameterized by the distance between the control points. No go:

julia>     CubicSplineInterpolation((s,), xyz[:, 1])                                                                                                                                                      
ERROR: MethodError: no method matching CubicSplineInterpolation(::Tuple{Array{Float64,1}}, ::Array{Float64,1})                                                                                            
Closest candidates are:                                                                                                                                                                                   
  CubicSplineInterpolation(::AbstractRange, ::AbstractArray{T,1} where T; bc, extrapolation_bc) at C:\Users\PKrysl\.julia\packages\Interpolations\9sq9m\src\convenience-constructors.jl:8                 
  CubicSplineInterpolation(::Tuple{Vararg{AbstractRange,N}}, ::AbstractArray{T,N}; bc, extrapolation_bc) where {N, T} at C:\Users\PKrysl\.julia\packages\Interpolations\9sq9m\src\convenience-constructors.jl:19                                                                                                                                                                                                    
Stacktrace:                                                                                                                                                                                               
 [1] top-level scope at REPL[25]:1  

If I replace cubic interp with LinearInterpolation, it works.

There is also GitHub - slabanja/SimplePCHIP: Simple monotone piecewise cubic interpolation.

I think Interpolations.jl only supports unequally spaced points for linear interpolation?

You might try https://github.com/kbarbary/Dierckx.jl

SimplePCHIP might be a good reference for a Julia Hermite interpolation package, but it is a bit out of date and the implementation looks rather Pythonic (not very type generic, and some potential performance problems).

Thanks for the hint, I tried to browse through the tests and it looks like you are right.
It doesn’t work for cubic curves.

Beautiful! Why can’t Interpolations.jl be so simple to use?

2 Likes

Thanks. This is good to know, but for my present purpose the PCHIP is too flat inbetween the control points. I really do prefer the regular spline in this case.

https://htmlpreview.github.io/?https://github.com/UMCTM/DataInterpolations.jl/blob/master/Example/DataInterpolations.html

3 Likes

Dierckx.jl? It is an easy to use package, but I don’t know about unequally spaced data.

I agree that Interpolations.jl should be reworked to be easier to use. It could really use a dedicated maintainer. It’s a nicely-isolated package (it doesn’t require you to maintain an entire ecosystem), mathematically interesting, and has huge implications for the community. Volunteers wanted!

9 Likes

I would like to contribute more to the Julia community but I’m not sure I would be qualified to be a maintainer. I wouldn’t mind focusing my extra time on contributing to the package rework though.

4 Likes

That would be awesome. I’ve been trying to review PRs to Interpolations when I can find the time; if you submit a few, not only will you be helping make it better but you’ll be in training for having more authority over the direction of the package.

2 Likes

Should I start by tackling issues or do you have a vision/roadmap for rework?

If there are some small issues easy to fix, that might be a great place to start.

My own vision looks something like this:

  • better documentation. Right now I think the documentation has been written by the same people who wrote the code, and that usually gives mixed results.
    The conundrum is that one could document the API and then it might change, so it’s a bit of a question of how much effort to throw in here.
    One fairly safe step might be to improve documenting the internals in a manner that teaches how the package works. (Obviously start by checking the current writeups and seeing whether you think they’re adequate or not; maybe it’s better than I think.) There are certain parts of the package that I think are solid, notably the WeightedIndex infrastructure and probably the prefiltering. The closer one gets to the user level, though, the more likely I think some of it will change in the long run.
  • check to see if the benchmarks work and whether we need more. See if work needs to be done to make running them easy, possibly as part of CI. This would allow us to detect whether future changes cause performance regressions.
  • start thinking about reworking the API. Can we take useful inspiration from other interpolation packages? Can we leverage inferrable keyword arguments (which Julia didn’t have when Interpolations was designed) to deliver a more convenient interface? The whole OnGrid/OnCell issue has added a lot of complexity, can that be simplified? https://github.com/JuliaMath/Interpolations.jl/issues/227 and https://github.com/JuliaMath/Interpolations.jl/issues/242 represent probably the biggest architectural changes but both serve as roadmaps for important changes.
4 Likes

It’s official: Interpolations.jl needs a new maintainer

2 Likes

[quote

I found this package is really neat. Is there any plan to add derivative and integration functions just like provided in Dierckx.jl? @ChrisRackauckas

1 Like
1 Like

@ChrisRackauckas. I have some questions. Liking DataInterpolations.jl a lot.

  1. Does it support extrapolation?
  2. Are there plans for multidimensional extrapolation/interpolation?
  3. Relative to other interpolation packages, why switch the positions of x,y (see below)?
using DataInterpolations
x = collect(1:0.1:5)
y = rand(length(x))
itp1 = LinearInterpolation(y, x)
@show itp1(4)
@show itp1(6)

I believe it does now.

Not really: we didn’t need it for our application.

It matches the differential equation solver, so the convention is “correct” :stuck_out_tongue:. Yeah… that’s how conventions go.

1 Like

Thanks @ChrisRackauckas. Great package.