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)