Decision: Separate the concepts of boundary condition and extrapolation behavior

[tomasaschan]

After thinking a little about #32 as well as poking around in the code to see what I can do about the problems with InterpCubic (#29, #31, #33) I'm getting more and more frustrated with the lack of separation between the two concepts "boundary conditions" and "extrapolation behavior". As I see it, those are two separate concepts:

• Boundary conditions define how the interpolation behaves close to, but inside the edges of the domain. For example, the current behavior of quadratic interpolations is (IIUC) to use the coefficients for x in [2,3] also for x in [1,2]. An alternative way of terminating the interpolation would be to use a linear interpolation in the edge interval, and still require smoothness at the edge.
• Extrapolation behavior defines what happens if an x is provided which is outside the domain (e.g. any x < 1). Currently, there are a couple of behaviors: BCnil and BCnan both do something to indicate an undefined value, while BCfill defines a constant extrapolation independent of the function value at the edge, and BCnearest is in effect a simple constant extrapolation. But this is (could be) independent of what happens inside the domain - it is only concerned with behavior outside.

Some of today's boundary conditions capture both these concepts - BCperiodic and BCreflect intrinsically defines both the boundary condition and the extrapolation behavior - but that's more a special circumstance than anything else.

Separating these two concepts would probably be a quite substantial rewrite of parts of interp.jl, and breaking API changes are likely (if they are separate concepts, both need to be specified by the user...). This is unfortunate, since Grid is used by a lot of people, and in a lot of packages.

Before I start working on anything, I'd love to hear some input from others on whether this is a good idea or not, and (if it is) if there are any special considerations that I might not have thought of.

[timholy]

I think that's a good idea. The main concern will be the combinatorial increase of cases to handle.

#36, if that happens, might make some of this easier. (Depending on how we do it; I should say, it could conceivably also make this harder!) The basic idea I have in mind for #36 is that rather than writing a generic interp and getting all specific behavior through wrap functions, just implement each combination of interpolation order/boundary condition as a separate function. That makes it easier to omit irrelevant calculations, and hence will improve performance. But for this issue, it would also make it easier to introduce special treatment of edges/boundaries, if we find we need it.

It's a little hard to know how this would play out until we try. I wonder if I need to at least post a couple of examples of what #36 might look like before you make big decisions about this issue?

[tomasaschan]

Yeah, if we're going to make a rewrite anyway for performance reasons, this can very well wait to be included there.

I think we can still do quite a lot of stuff with generic functions and get this working well, though - the boundary conditions basically affect the equation system to solve (i.e. the various incarnations of interp_invert_matrix in the current implementation) while the extrapolation behavior defines what to do if x < 1 || x > len, which is something we don't really do today (for the few cases where it should be done, we cheat and set x=1...).

[tomasaschan]

Since this is handled in Interpolations.jl, I don't think Grid.jl needs to be concerned about this.