make vector indexing work with all extrapolations

[marius311]

Makes e.g. this work,

julia> itp = extrapolate(interpolate([1,2],BSpline(Linear()),OnGrid()),0)
2-element Interpolations.FilledExtrapolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Float64}:
 1.0
 2.0

julia> itp[[0,1,2,3]]
4-element Array{Float64,1}:
 0.0
 1.0
 2.0
 0.0

with all extrapolation types (previously only some worked).

I just let getindex fall through to the default implementation, and added a checkbounds which always succeeds since we are extrapolating to no-matter-where. Tests seem to be passing including the new ones I've written, but I'm not at all familiar with this code-base, so someone please check me!

I think this also fully closes #128

[sglyon]

@marius311 thanks for being willing to open a pull request here.

Can you provide an example of what didn't work before? Is the no-op implementation of checkbounds necessary for your example to work?

[marius311]

Sure, here's the above output before this fix (on 0.6),

julia> itp = extrapolate(interpolate([1,2],BSpline(Linear()),OnGrid()),0)
2-element Interpolations.FilledExtrapolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Float64}:
 1.0
 2.0

julia> itp[[0,1,2,3]]
ERROR: BoundsError: attempt to access 2-element Interpolations.FilledExtrapolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Float64} at index [[0, 1, 2, 3]]
Stacktrace:
 [1] throw_boundserror(::Interpolations.FilledExtrapolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Float64}, ::Tuple{Array{Int64,1}}) at ./abstractarray.jl:426
 [2] checkbounds at ./abstractarray.jl:355 [inlined]
 [3] macro expansion at ./multidimensional.jl:441 [inlined]
 [4] _getindex at ./multidimensional.jl:438 [inlined]
 [5] getindex(::Interpolations.FilledExtrapolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Float64}, ::Array{Int64,1}) at ./abstractarray.jl:875

And yea, that checkbounds overrides the default one which throws the error you see above.

Also, I think I need to fix that 0.4 unittest problem, do you happen to know what's wrong with the @inferred usage there?

[sglyon]

Ahh I see. Make sense.

Let's remove the ::Numbertype restriction on the xs...

Also I think that 0.4 might be having a hard time recognizing the [[3]] as a function call to getindex given everything else that is going on in that line. To test that maybe you could actually call getindex(stuff, [3]) and see if it helps?

[marius311]

Let's remove the ::Numbertype restriction on the xs...

That has to be there actually, so that the call to getindex first goes through the default stuff in Base which handles the logic of expanding out vector indexing.

I think that 0.4 might be having a hard time recognizing the [[3]]

Good call, that could well be it. I'll give it a go.

[marius311]

Ok I think that fixed 0.4 but the 0.5 tests fail now, but maybe just randomly stalled? Its fine locally for me.

[sglyon]

I've restarted the tests on 0.5.

I would like to have the opinion of at least one other person regarding the addition of the ::Number type restriction. I have only used this package to interpolate with numbers, but I've seen examples floating around for interpolating things. cc @timholy @ararslan @mbauman

[marius311]

Ahh I understand the concern now. FWIW, glancing through the source code most of the custom getindex methods defined in Interpolations.jl accept only Number, so this makes this consistent with that. In fact, this is basically how I arrived at this solution I give here, I was trying to understand why e.g. interpolation worked with vector indexing but not extrapolation, and it was basically because the getindex'es for interpolation methods were defined with ::Number, which led to the call falling back to the logic in Base first, so I made extrapolation do the same.

[marius311]

I do agree it'd be nice to allow any generic object, although that's probably beyond this PR/me.

[timholy]

Definitely the right idea. I'm not certain the Base mechanisms will work entirely, but this seems to be a step in the right direction.

One category where we'd be better off writing our own vector-indexing methods is for irregular grids; we can save ourselves some lookup time in that case. But I vaguely remember we might have already done that.

[timholy]

We should change this now to xs::Vararg{Number,N} (insist that N xs are provided).

[timholy]

Actually, I'd take that back. By putting on the Vararg constraint, you'd suddenly allow trailing 1s and indexing with fewer indices than dimensions (Base machinery will adjust the inputs to produce N outputs). So by not constraining to N, you ensure that this method "intercepts" the call.

[timholy]

Does it work with a less trivial test? E.g.,

etp = extrapolate(interpolate([0,0],BSpline(Linear()),OnGrid()),E)
@test @inferred(getindex(etp, [1.2, 1.8, 3.1])) === [0.0, 0.0, 0.0]

?

[marius311]

It does, I added this as a better test.

[mbauman]

With regards to allowing non-numbers, there's a challenge in structuring dispatch in such a way that you can differentiate between scalar and non-scalar indexing.

In 0.6, you can use Base's non-scalar indexing fallbacks with custom scalar indices with a setup like the following:

getindex(::ArrayType, ::ScalarIndex...) = # scalar indexing method
to_index(::ArrayType, i::ScalarIndex) = fill(i)

It's a little oblique, but the idea is that to_index will only get called by the base callbacks when your scalar specialization didn't take priority. That is — it'll happen for non-scalar indexing. And for non-scalar indexing, we can "protect" those scalar types that Base doesn't know about by wrapping them in a 0-dimensional array in order to ensure that they are given the correct scalar-like meaning in those fallbacks (it really should be a custom immutable 0-d array to avoid allocations). Now you might need to teach Base how to check bounds with your scalar types, and then you're good… just so long as the scalar indexing expression still returns elements of the array's eltype(). That last part is the kicker for interpolations — and I have an open PR at JuliaLang/julia#12567 that aims to solve it. The test in that PR should be particularly instructive.

This becomes significantly harder if you want to use ::Any... as your scalar indexing method because you'll now completely shadow the base fallbacks. And it's not really ::Any... that you want — it's any arguments so long as none are arrays. It's still possible, but you'll need to some work up-front to determine if you want to use scalar or non-scalar indexing methods, and then you'll have to directly call the undocumented internal interfaces to hit the fallbacks.

[sglyon]

Thanks for the comments @timholy and @mbauman .

From what I understood (I'm still working on wrapping my head around the new 0.6 indexing fallbacks) it seems like this PR is in mergeable condition, modulo addressing @timholy question about the tests. Is that right?

[timholy]

Personally I'd be fine with dropping support for julia-0.4. Might be the most expedient way to get this passing tests on all supported versions.

[sglyon]

Thanks @marius311 !

[marius311]

Thanks for the help/comments!