Convenience constructors for Interpolations.jl

[jlperla]

I find the construction methods for Interpolations.jl to be almost impenetrable. It would be really nice to have some convenience constructors to make things easier to use. Coupled with the () access, I think this makes Interpolations.jl feasible for teaching.

I can't figure out whether we would usually want https://github.com/JuliaMath/Interpolations.jl#gridded-interpolation or https://github.com/JuliaMath/Interpolations.jl#scaled-bsplines ?

The big scenarios I can think of to make easy in 1 dimension are:

• Linear interpolation on a regular grid
• Linear interpolation on an irregular grid
• Cubic spline interpolation on a regular grid
• Cubic spline interpolation on an irregular grid.
• Look into the same for 2 or more dimensions

I think it would make sense to figure out:

• Document the canonical way to do each of them in the current Interpolations.jl
  • Where possible, we can discuss reasonable default parameters
• Do we need to get Cubic spline interpolation working in the gridded setup in order to support cubic splines on an irregular grid?
• See what sort of convenience constructor would make sense so these are easier to use.

[jlperla]

A much saner interface for construction (if not calling) is https://github.com/kbarbary/Dierckx.jl

[jlperla]

Note that I found a PR in progress for #48 which we might be able to help with. Sounds like gridded cubic splines are not supported yet. But what about cubic splines on a uniform grid? This is where the scaling comes in that I just don't understand.

[jlperla]

Saw this: https://discourse.julialang.org/t/julia-version-of-matlab-function-interp1/8540 Not sure it is the best idea to do it that way, but gives you a sense of the usage patters.

[chiyahn]

I made the following tutorial file that lists all possible scenarios and their solutions: https://github.com/chiyahn/notes/blob/master/interpolations/interpolation-scenarios.md including some suggestions for convenient constructors as well!

For cubic splines (or other nonlinear splines) on irregular grids, I can think of the following two solutions as I explained in the markdown file:

1. Construct interpolations using a recursive algorithm: Slow, but trivial. I tried implementation a bit by myself but it seems like integration with the current code is not trivial and requires substantial changes.

   • ref: https://github.com/floswald/ApproXD.jl

2. Construct gridded interpolation algorithm dedicated for cubic splines: Fast, but requires substantial amount of verification.

   • ref: Help with gridded cubic splines #48 (PR you mentioned)

[jlperla]

Great, thanks! But I didn't see the convenience constructors?

Also, see the final comments on #26 to get it ready to merge.

[jlperla]

Sorry, now I see the constructors. I don't love the interp1 approach especially much, since I am not sure how well it supports uniform vs. irregular grids, andCubic(Line()) is still pretty ugly. Also, I am not 100% sure that it needs to be written to be written specific to one dimension, since we can use Julia dispatching.

One thing that could help is that we can dispatch based on the StepRangeLen type, i.e.

julia> typeof(1.0:0.1:2.0)
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

julia> typeof(1:10)
UnitRange{Int64}

Is that a way to have a single interface support both irregular and regular grids? Maybe the convenience constructor should have the name of the interpolation method rather than the dimension? That is,

f(x) = sin(x) / 2
xs = linspace(1.0, 2.0 ,10) #Note, this does not return a vector!
xs_irregular = collect(xs); #Now this is an an AbstractVector  and no longer a StepRangeLen

 #This knows it is a regular grid since xs is StepRangeLen or a StepRange
LinearInterp(xs, f.(xs))

 #This could use the fact that typeof(1:10) is a UnitRange to realize that rescaling is not needed?
LinearInterp(1:10, f.(1:10))

#Cubic Spline knowing that xs is a StepRangeLen
CubicSplineInterp(xs, f.(xs)) #This knows it is a regular grid since xs is StepRangeLen

#Now, this uses the irregular grids, since it dispatches on the AbstractVector rather than the UnitRange or StepRangeLen?
LinearInterp(xs_irregular , f.(xs_irregular ))

#This isn't functional yet..., but could be added later.
CubicSplineInterp(xs_irregular , f.(xs_irregular ))

If things moved to multiple dimensional interpolation, then I think that this pattern for convenience constructors would work well, we would just need to dispatch on what xs looks like for a uniform grid.

Note that with this, we could write a 1D LinearInterp for the non-uniform case that completely moves out (and deprecates) https://github.com/QuantEcon/QuantEcon.jl/blob/master/src/interp.jl

[chiyahn]

I see your point! I was trying to imitate what MATLAB users are going to be more familiar with, but I agree with you that writing a constructor based on degrees (rather than dimension) is a more reasonable choice. The best way to handle it would be to create overloading constructors that are dependent on types of parameters so that (at least for linear case) it returns a scaled interpolation instance if the grids are regular and gridded one if they are not.

function LinearInterp(xs::AbstractArray,ys)
    itp = interpolate((xs,), ys, Gridded(Linear()))
    return(itp)
end

function LinearInterp(xs::Range,ys)
    itp = interpolate(ys, BSpline(Linear()), OnGrid())
    return(scale(itp, xs))
end

The following unit tests passed:

# test regular grid cases
xs = linspace(xmin,xmax,N)
ys = f(xs)
itp = LinearInterp(xs, ys)

# unit tests
@test itp[xs[1]] == f(xs[1])
@test itp[xs[2]] == f(xs[2])
@test itp[(xs[1] + xs[2]) / 2] == (f(xs[1]) + f(xs[2])) / 2

# test irregular grid cases
xs = [x^2 for x = 1:xmax]
ys = f(xs)
itp = LinearInterp(xs, ys)

# unit tests
@test itp[xs[1]] == f(xs[1])
@test itp[xs[2]] == f(xs[2])
@test itp[(xs[1] + xs[2]) / 2] == (f(xs[1]) + f(xs[2])) / 2

I will look into multi-dimensional cases tomorrow as soon as the survey conference (😭) is finally done.

[jlperla]

Great!

Though before multidimensional cases, I would try to finish JuliaMath/Interpolations.jl#206 (comment) so that we can tell people to merge your PR.

A few small things on this, which seems extremely close:

• The test failed for me for my f because of the exact floating point comparison, but I think you meant to use ≈ instead
• Probably worth adding LinearInterp(1:10, f.(1:10)) to the testsuite (i.e., a UnitRange) but it seems to work.
• This approach seems to work well with Cubic Splines as well. That is, I think it is as simple as

function CubicSplineInterp(xs::Range,ys)
    itp = interpolate(ys, BSpline(Cubic(Line())), OnGrid())
    return(scale(itp, xs))
end

• If an irregular grid version was eventually added, then CubicSplineInterp(xs::AbstractArray, ys) would be added.
• In the future: I think the more idiomatic way to return from small functions in Julia is to avoid the return, but I suspect it wouldn't change the compiling behavior in this case. That is,

CubicSplineInterp(xs::Range, ys) = scale(interpolate(ys, BSpline(Cubic(Line())), OnGrid()), xs)

• What are your thoughts about extrapolation in the convenience constructor?
  • Mine is that it should be avoided so there are no silent errors from misuse (and that people should use the full interface if they like). I have already made that mistake in testing
  • I think throwing is done with something like

CubicSplineInterp(xs::Range, ys) = extrapolate(scale(interpolate(ys, BSpline(Cubic(Line())), OnGrid()), xs),  Interpolations.Throw())

but we might want to profile it to see if this slows things down considerably.

[jlperla]

I did a basic profiling with using BenchmarkTools and @btime. For a simple lookup, it seems like itp4 takes 72ns with extrapolation and 92ns with no extrapolation. I think this is a reasonable price to pay for safety with a convenience constructor, but we can get other peoples feedback

[chiyahn]

• For multidimensional cases, it turns out I have to change the order of parameters for constructors from (ranges, values) to (values, ranges...) to allow arbitrary dimensions for ranges. Gridded interpolations (for irregular grids) work in the following version:

LinearInterp(vs, ranges::AbstractArray...) = interpolate(ranges, vs, Gridded(Linear()))

• But for scaled interpolations (for regular grids), there's a bit of problem; the following constructor seems to be a trivial choice:

LinearInterp(vs, ranges::Range...) = scale(interpolate(vs, BSpline(Linear()), OnGrid()), ranges)

but running the following lines returns an error:

xs = 1:.1:10
ys = 1:.5:10
f(x, y) = log(x+y)
A = [f(x,y) for x in xs, y in ys]
LinearInterp(A, xs, ys)

The reason is that when xs and ys are passed as parameters for ranges, they are going to be provided as a single Tuple instance while scale does not have any overloading constructor with tuples. In fact, using implementation of scale, the following version works:

LinearInterp(vs, ranges::Range...) = ScaledInterpolation(interpolate(vs, BSpline(Linear()), OnGrid()), ranges)

• On the other hand, notice that this bypasses check functions (for dimension compatibility, etc.) in scale; I am not sure why maintainers decided to make these two different versions (for safety?), but it seems like the only difference between ScaledInterpolation and scale is coming from these check functions. So we can proceed with one of the followings:

  • Add these check functions in constructors manually; this is the simplest and safest way to go, but makes the constructor look complicated and there's a possibility of being inconsistent with check functions in scale when there are changes in future
  • Stick with the version above; this makes the code neat and consistent with the other scripts in the package, but might be unsafe as users might not be aware of what the errors are about.
    I personally prefer the first second option, but I want to hear your opinions too. 😃

• I agree that it should be avoided. Calling extrapolation is relatively easy too.

• I agree that it is a reasonable price to pay. My understanding is that the value of having a safe & convenient constructor outweighs its computational cost as long as it is reasonable (92ns is) since it is relatively too small compared to the computation time needed for interpolation evaluations anyway.

[jlperla]

Interesting. So that vs in LinearInterp(vs, ranges::Range...) would be on the grid, with the convention that the ranges.... represent the cartesian product. Are you sure it isn't possible to do LinearInterp with tuples and splat? That is something like (the untested)

LinearInterp(ranges, vs) = ScaledInterpolation(interpolate(vs, BSpline(Linear()), OnGrid()), ranges...)
xs = 1:.1:10
ys = 1:.5:10
f(x, y) = log(x+y)
A = [f(x,y) for x in xs, y in ys]
LinearInterp((xs, ys), A)

I think it is reasonable to bypass the checks for simplicity. I like clear code, and the whole thing is kind of unsafe.

[jlperla]

Another, fancier, way to do it is to look into ways to represent cartesian products with operator overloading. We can see if http://juliaapproximation.github.io/ApproxFun.jl/latest/usage/domains.html does anything interesting... For example, instead of a tuple it might help to be explicit.

f(x, y) = log(x+y)
A = [f(x,y) for x in xs, y in ys]
LinearInterp(ProductDomain(xs, ys), A)
#(or with overloading)
LinearInterp(xs ⊗, ys, A)

[jlperla]

Check out QuantEcon/Expectations.jl#9 for some discussion of the generic, abstract types for ranges.... this will matter when it comes to the interface to call with uniform grid methods.

[chiyahn]

Oh, I see your point; in fact the code you provided works by changing ScaledInterpolation to scale, which takes splat on the second parameter:

LinearInterp(ranges, vs) = scale(interpolate(vs, BSpline(Linear()), OnGrid()), ranges...)

but the challenge here is that Tuple is not a UnionAll type so that when we want to specify the type Range the following line is not going to work with 2-dimensional tuples:

LinearInterp(ranges::Tuple{T}, vs) where T <: Range = print(typeof(ranges))

in fact, one has to define it alternatively as LinearInterp(ranges::Tuple{T,T}, vs) where T <: Range = print(typeof(ranges)), specifying the dimension as well, which is not desirable.

We can, for instance, alternatively take Array instances that are UnionAll as an alternative choice for grids, but it is not trivial on how we should handle the cases when the length for grids in the first and second dimension differ. Looking into other ways to represent cartesian products will be a good idea as well, but after all I think having the parameter for grids located after the value matrices as I suggested still makes sense, except 1D cases (as it is trivial to have it in the order of (x, y)).

[jlperla]

On 0.6,

julia> LinearInterp(ranges::NTuple{N,T}, vs) where {N,T <: Range} = print(typeof(ranges))
LinearInterp (generic function with 1 method)

julia> LinearInterp((1:10,2:5), 1)
Tuple{UnitRange{Int64},UnitRange{Int64}}

julia> LinearInterp((1:0.1:10,2:0.1:5), 1)
Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}

Still not 100% sure of the right approach. Regardless of which argument comes first, my feeling is that the ranges really represents some coherent "thing" and that variable argument parameters to the function do not capture that. Having it as a single object allows us to change behavior later to be more in line with the spirit of http://juliaapproximation.github.io/ApproxFun.jl/latest/usage/domains.html

If what I suggested works, maybe we can get a tuple based approach and then ask for people's feedback.

[chiyahn]

Whoa, didn't think about NTuple. I agree that having grids as a single object is a better approach, but I thought it was not trivial. In fact, the following implementation works:

LinearInterp(ranges::NTuple{N,T}, vs) where {N,T <: Range} = scale(interpolate(vs, BSpline(Linear()), OnGrid()), ranges...)
LinearInterp(ranges::NTuple{N,T}, vs) where {N,T <: AbstractArray} = interpolate(ranges, vs, Gridded(Linear()))

I tried some unit tests and they have all passed: https://github.com/chiyahn/notes/blob/master/interpolations/multivariate-constructor.md; we can use them for unit tests to be included in PR. Do you want me to make a PR for linear cases first?

[jlperla]

Maybe add in the throwing extrapolation as well? Lets stick with the Ntuple for now,

Why don't we get the full version up and running in a PR first (maybe without docs) to get feedback. I think the key is

• Linear interpolation with univariate uniform grids (i.e. abstractranges)
• Linear interpolation with univariate non-uniform grids (i.e. abstractranges)
• Linear interpolation with multivariate uniform grids (i.e. abstractranges)
• Linear interpolation with multivariate non-uniform grids (i.e. abstractarrays)
• Cubic spline interpolation with uniform grids

And that is pretty much it? After we get this done, we are 90% of the way there to making the library usable. Then we can think about whether it is the right library to stick with in the middle-term.

The only other thing I would think about is whether we want to change the names? I asked Spencer and he doesn't see a reason to have this stuff using QuantEcon naming. My suggestion is LinearInterpolation and CubicSplineInterpolation? I tend to hate shortened words where it is tough to remember how it was condensed (e.g. remembering whether it was LinearInterp vs. LinearInter)

[chiyahn]

No problem, I'll add in throwing extrapolations as well and work on the version for cubic spline interpolation! Meanwhile, note that the current version supports both univariate and multivariate cases; unit tests I created support both cases.

And I agree -- I prefer longer names too.

[jlperla]

I think this is the right strategy, and lets not overthink the interface. These were intended to patch up Interpolations.jl to a point where it is usable for our needs, but ultimately it might make sense to build on new interpolations libraries. Tough to know when the timeline for that new library is feasible, though.

[chiyahn]

Okay, here I have made a pull request: JuliaMath/Interpolations.jl#214. Even though implementation itself didn't take that much time it took me a while to remind myself that I had to add new functions to export 😃..

[chiyahn]

FYI, here are unit tests I have run (included in test/convenience-constructors.jl in the PR as well): https://github.com/chiyahn/notes/blob/master/interpolations/multivariate-constructor.ipynb

[jlperla]

(edit: hit enter too early!)

Great!

The last thing I see is that I think we need a 1D version of the interfaces... the whole tuple of one element thing is ugly.

Maybe add in something like

LinearInterpolation(range::T, vs) where {T <: Range} = LinearInterpolation((range,), vs)

That way you can just go LinearInterpolation(xs, A)

I also agree with the commenter on the PR that it woudl be nice to be able to swap on/off extrapolation easily.  But can't think of any other features we need to implement.  

[chiyahn]

I think this can be closed now -- the PR (JuliaMath/Interpolations.jl#214) has been merged! 🎉 There are a few changes in README.md for quickstart guide that have to be done (JuliaMath/Interpolations.jl#216), but implementation is now basically complete.

[jlperla]

Great!