Feature wishlist

[tomasaschan]

Now that basic functionality is in place, we need to extend what we have with more cases, to make the library feature complete. That's a good case for a wishlist =) This issue is intended to keep track of overall progress toward feature parity with (interpolation parts of) Grid.jl - actual implementation details can be discussed in separate issues.

Boundary condtions (assuming a grid on 1:n)

• Flat, i.e. f'(x=1) = f'(x=n) = 0 (WIP features (so far: Flat BC, Constant extrap) #8)
• Line (or Natural), i.e. f''(x=1) = f''(x=n) = 0 (More boundary conditions for quadratic B-splines #9)
• Tangent, i.e. with specified derivatives at the edges (Flat is essentially Tangent(0,0))
• Reflect, i.e. all continuous derivatives are 0 at the edge, to preserve continuity properties when extrapolating (Added in Reflecting BC for quadratic itp #20)
• Periodic, i.e. function value and all derivatives are equal at both edges; with OnGrid, the data should include both cell boundaries (and A[1]==A[end]), while with OnCell this assumption is not necessary (but rather A[-.5] == A[end+.5]) (More boundary conditions for quadratic B-splines #9)
• Free, i.e. with extra continuous derivatives close to the edges, to close the system. (For e.g. cubic splines, which have continuous second derivatives, the third derivatives are also continuous at the second-to-last boundaries). (More boundary conditions for quadratic B-splines #9)

Extrapolation behaviors

• Constant, i.e. return the value at the closest data point (WIP features (so far: Flat BC, Constant extrap) #8)
• Linear (Linear extrapolation #23)
• Quadratic etc, valid up to the same degree as the interpolation. When the same degree, just continue the outermost polynomial; when lower degree, construct an extrapolation polynomial with as many continuous derivatives as necessary for a smooth transition. For example, linear extrapolation of a quadratic interpolation object, would have continuous value and first derivative (which conveniently specifies the extrapolation uniquely...).
• Periodic, i.e. f(x + T) = f(x) for a function with period T (More boundary conditions for quadratic B-splines #9)
• Reflec, i.e. f(end - 1 + T/3) = f(end - T/3 + 1) (More boundary conditions for quadratic B-splines #9)

Interpolation degrees

• Cubic b-splines (WIP on cubic #90 - thanks @spencerlyon2!)

Other functionality

• Derivative evaluation, up to the highest continuous degree supported by the interpolation object. E.g. gradient(itp::Quadratic, x) or hessian(itp::Cubic, x). (Gradient fixed in Gradient evaluation #18, higher orders will be added when higher order interpolations are added...)
• Multi-valued interpolation (i.e. interpolation that yields more than one value per coordinate), More of a decision: do we do anything special to make it easier, or just force users to shape their data accordingly? Discussed separately under Multi-valued interpolation #3
• Range/vector-valued indexing (i.e. indexing with multiple values at the same time, e.g. itp[1:.1:10]) (Vector-valued evaluation #24, Restrict scalar indexing to Numbers #54, Improve testing #55)
• Non-uniform interpolation (Non-uniform interpolation #11, Add GriddedInterpolation for non-uniform grids #56)
• Uniform interpolation on non-unit scales (the same problem that CoordInterpGrid in Grid.jl solves) (Something equivalent to CoordInterpGrid #25, RFC: Scaling of interpolation objects (fixes #25) #47)

[timholy]

Looks good! One point I'm not sure about is why you separate boundary conditions from extrapolation behavior; it basically seems like one follows from the other.

[timholy]

One more topic of importance: efficient multi-valued interpolation (e.g., timholy/Grid.jl#44). Obviously this will just "drop out" for people who use Arrays of immutables, but for other cases (like the one in that issue) one might need another approach.

My current thinking is that this is a case where it's best to wait for Julia to support fixed-size arrays efficiently, and not try to hack around this in Interpolations.jl. In the meantime, people who need that can use Grid. As a close facsimile, we could conceivably support arrays of tuples.

[tomasaschan]

Multi-valued interpolation is important indeed, but I'm more and more leaning toward supporting this simply by duck typing, and forcing the user to use other tools to get their data in an appropriate shape. As we've said previously somewhere, we currently only demand that the data type supports multiplication with floats and addition (i.e. are basically vectors in some vector space...), and I would argue that reshaping data that is in some other form (i.e. an image in an Array{T,3}, the third dimension being color spaces) is out of scope for Interpolations.jl.

Re: bc/extrap separation, there are basically two reasons: one is that they are basically independent in the implementation, so it seems like a bad idea to couple them more than necessary in the API, and the other is that we have a couple of extrapolation behaviors that don't map nicely to boundary conditions (e.g. ExtrapError), and it's therefore better to let the user specify them independently. For many BCs, we could probably have a constructor that selects a corresponding extrapolation behavior, so the users don't have to say e.g. Periodic twice.

[timholy]

but I'm more and more leaning toward supporting this simply by duck typing

that's basically what I'm saying, too. (Your version is clearer, though.)

Re the boundary conditions: got it. Carry on 😄.

[albop]

The wishlist is promising !
About extrapolation, I also think that it is nice to separate boundary conditions and extrapolation. For some applications, using linear extrapolation, at any order, improves the chances that the extrapolated behavior is not crazy, not returning to zero or wildly diverging.

[TheBB]

Typical boundary conditions used for cubic splines are:

• Free: continuous third derivatives in the second and second-to-last points.
• Natural: zero second derivative at the boundaries.
• Tangent: specific derivatives at the end. This is a special case of Flat. It would be good to have an inhomogeneous Flat type for this, and let Flat = Tangent(0, 0) or something along those lines.

[tomasaschan]

@TheBB Thanks for the input! I think the natural boundary conditions will be identical to a linear BC, since the way to implement the latter is also to set the second derivative to zero at the boundary. I'll update the wishlist with the other two.

[kzapfe]

I have some of a theorical/implementation question:
Why packages such as Pyplot have such a large set of interpolation routines to plot functions over 2D domains ("images") and on the other hand turns out to be so cumbersome to implement the same interpolations "just to the data"? What I am missing here?
I think I have the same problem as other people here: I need continous, reliable Hessians for my work. I see what can I understand of your package and see if I can contribute. Keep the splendid work!

[tomasaschan]

@kzapfe PyCall/PyPlot are of a very specific type of Julia packages - they're both wrappers for code in other languages that have had much more time than Julia to evolve and mature. Thus, comparison to pure Julia packages like Interpolations.jl isn't quite "fair" to either type of package. Since PyPlot is built specifically for showing things on-screen, it's maybe not so surprising that it focuses its efforts on 2D. Did you have any other packages with the same limitations in mind?

It's also a much larger problem to write general code that works in any number of dimensions than to write for just one case, so if you only need 1D and 2D there's no reason to build something more general.

That said, you can probably find a lot of good interpolation code in various languages that already has been/would be easy enough to wrap in Julia and that does what you need.

My ambition with this library is to, eventually, eliminate the need for calling out to other languages, at least if you don't need something very specific - but part of the reason for building it is to have an excuse to practice Julia programming :)

[abahm]

This looks great - thanks for all your work! I'm looking at a task that would require barycentric Lagrange interpolation, which might fit here. Is anyone already working on this?

[tomasaschan]

@abahm Not that I know of. If you feel up to it, a PR is always welcome :) The API should probably be something like that for gridded interpolation, e.g. interpolate(knots, data, Lagrange) or even interpolate(knots, data, Lagrange(Barycentric)) to make room for regular Lagrange interpolation too.

[sglyon]

I'm going to try to get my hands dirty with cubic splines! I might reach out for help if I get stuck, hope you don't mind.

[sglyon]

Also, as I go would you be OK with me trying my best to add docstrings to currently undocumented functions? That would help me learn what everything does.

[tomasaschan]

@spencerlyon2 Great! I'm really excited about this project :)

Docstrings are great (we have already discussed moving most of the documentation into them, and away from the Readme, see #46) - but you might also find the documentation in doc/latex useful when figuring out how the current B-spline functionality fits together. I'd love if you could add a few sections for the cubic spline interpolation as well (you can probably copy most of the text, and just update the matrix equations).

My main concerns with this library is that the API stays consistent, and keeps leaving room for incorporating new interpolation schemes, and that we ensure that the composability that we have between the current implementations of b-splines, extrapolation and scaling, is kept for other schemes as well. Together with Tim's focus on squeezing out every single ounce of performance I really think we've found a sweet spot :) I have the ambition of helping you shape your contributions to fit well into this framework. I'm telling you this to make it clear from the start that I do value your contributions - and then I hope you'll have some patience if I'm nit-picky sometimes in my code review... ;)

Feel free to file issues with questions whenever you need it - we'll do our best to answer them promptly.

[timholy]

Yeah, this sounds fantastic. Do shout for help if you find any of the metaprogramming-heavy stuff confusing.

[greimel]

Is there a way of performing cubic spline interpolation on a non-uniform grid with this package?

If no, is there a specific reason? It seems to me that the generalization from uniform grid to a non-uniform should be independent of the uniform interpolation scheme.

My workaround so far looks like

function itp_nonuniform(val::Real, itp_uniform, x, y)
    i = 2
    while val > x[i]
        i += 1
    end
    λ = (val - x[i-1])/(x[i] - x[i-1])

    itp_uniform[i - 1 + λ]
end

itp_nonuniform(vals::AbstractVector, itp_uniform, x, y) =
    broadcast(val -> itp_nonuniform(val, itp_uniform, x, y), vals)

Some illustration

using Interpolations, Plots

N = 40
xmax = 2*round(π, 4)

x = xmax * [0; sort(rand(N-2)); 1]
y = sin.(x)

itp_uniform = interpolate(y, BSpline(Cubic(Natural())), OnGrid())

scatter(x, y, label="true")
vals = linspace(0, xmax, 5*N)
plot!(vals, itp_nonuniform(vals, itp_uniform, x, y), label="interp")

EDIT: I probably should've read #131 before my comment. My workaround is most likely doing something different than gridded cubic spline interpolation.

[dagap]

One thing that would be nice to have is also sinc interpolation.

[sglyon]

Thanks @dagap I've opened a separate issue #145

[babaq]

i wish this package include smooth splines too. recently i tried to convert some R code to Julia which uses R smooth.spline function. The Dierckx.jl can do smooth spline, but the only parameters is a smoothness. R version accept degree of freedoms which is very useful to control not only smoothness but also the complexity of the spline. right now i am using RCall, but hope to replace the hack some day.

another thought is why not use callable overload to sugar coat the interpolation. Dierckx.jl makes the splines callable, so instead of s[newx], i found it more natural to consider the splines as a function s(newx) since most of the time, we use it to fit a smooth functions. of course, this can be added easily, is there any confusions were both be added?

[yakir12]

I'll second @babaq's wish/need for a smoothing functionality.

I will admit though, that it is arguable if smoothing an interpolating object is within the scope of what an interpolating library does. The main objection I can think of is that once you smooth the interpolation it isn't interpolating between the data points any more. The fact remains though that people come looking for some smoothing functionality when interpolating. As the main interpolation package for Julia, I think there should be some information about smoothing the interpolation/data:

1. why Interpolations.jl doesn't do it
2. some suggestions on where and how to do it somewhere else
   Ideally there would be some companion package that does it...

[PallHaraldsson]

@gerlero Should monotonic, i.e. https://github.com/gerlero/PCHIPInterpolation.jl be on the list here?

I used PCHIP from SciPy at the time when I needed it (and extrapolation) and it wasn't then available in Julia (when porting from MATLAB). Maybe it makes sense to add that package as a dependency (I believe its BSD license may be compatible with MIT here), not to reimplement. If not add it as a forth alternative in the README?

[mkitti]

be on the list here

No. Please create a new issue.

[PallHaraldsson]

I rather just added a doc to README #493, at least for now. Did you for sure mean to close since 2 of the tasks not yet completed?

[mkitti]

Necroposting on an 8-year old issue where most of the participants have moved on is simply not helpful. For the current maintainers it is also very hard to find and track.