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.