Require that singleton dimensions be NoInterp. Fixes #289.

[timholy]

CC @zygmuntszpak. I wonder if this could be breaking, so I'll leave it open a while for comment.

[codecov-io]

Codecov Report

Merging #290 into master will increase coverage by 0.2%.
The diff coverage is 75%.

@@           Coverage Diff           @@
##           master   #290     +/-   ##
=======================================
+ Coverage    47.8%    48%   +0.2%     
=======================================
  Files          21     21             
  Lines        1071   1079      +8     
=======================================
+ Hits          512    518      +6     
- Misses        559    561      +2

Impacted Files	Coverage Δ
src/b-splines/b-splines.jl	46.83% <75%> (+3.17%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0b3fb7a...c78fc96. Read the comment docs.

[zygmuntszpak]

Thank you! This will break imresize which uses interpolations under the hood in the extreme case where someone is trying to resize an image consisting of a single row into multiple rows.

using Images, Interpolations
z = reshape(collect(0:1:255),1,256)
# This will now throw the new error message that you created
imresize(img, (z, 256))

I suspect this is unlikely to break peoples code in the wild though. I only discovered this because I was working on test cases for a new CLAHE implementation.

I can get things to work using your proposed fix. One corner case that I have is a situation where I would like to create:

        # A is a 1x1 matrix, e.g. [5]
        itp = interpolate(A, (NoInterp(), NoInterp()))
        sitp = scale(itp, Base.OneTo(1), Base.OneTo(1))
        etp =  extrapolate(sitp, Flat())

and I would like the NoInterp() operation to just return the value stored in A and etp with Flat() to also return the value stored in A.

If this is not desirable, not possible or too much work then no worries. I would like this property to allow me to handle a particular corner case, but I can always exclude the corner case from consideration.

[timholy]

It's worth discussing some of the subtleties here. The way this typically gets evaluated can be revealed with Rebugger (note I'm running Julia with --check-bounds=yes, and I'm using master, not this branch):

julia> a = [i for i = 1:3, j = 1:1]
3×1 Array{Int64,2}:
 1
 2
 3

julia> itp = interpolate(a, BSpline(Linear()))
3×1 interpolate(::Array{Float64,2}, BSpline(Linear())) with element type Float64:
 #undef
 #undef
 #undef

# Note: the next comes from stepping into `itp(1.1,1)`:
(itp::Interpolations.BSplineInterpolation{T,N,TCoefs,IT,Axs} where Axs<:Tuple{Vararg{AbstractUnitRange,N}} where IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, BSpline},N} where N}, BSpline} where TCoefs<:AbstractArray)(x::Vararg{Number,N}) where {T, N} in Interpolations at /home/tim/.julia/dev/Interpolations/src/b-splines/indexing.jl:6
  itp = [#undef; #undef; #undef]
  x = (1.1, 1)
  T = Float64
  N = 2
rebug> @eval Interpolations let (itp, x, T, N) = Main.Rebugger.getstored("0a91928a-06bf-11e9-2d75-97032c03971d")
       begin
           @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
           wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
           @show wis   # this line inserted during Rebugging
           itp.coefs[wis...]
       end
       end
wis = (Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8999999999999999, 0.10000000000000009)), Interpolations.WeightedAdjIndex{2,Int64}(0, (0, 1)))
ERROR: BoundsError: attempt to access 3×1 Array{Float64,2} at index [2, 0]
Stacktrace:
 [1] getindex at ./array.jl:732 [inlined]
 [2] interp_getindex at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:298 [inlined]
 [3] macro expansion at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:282 [inlined]
 [4] interp_getindex at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:273 [inlined]
 [5] macro expansion at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:282 [inlined]
 [6] interp_getindex at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:273 [inlined]
 [7] _getindex at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:238 [inlined]
 [8] getindex(::Array{Float64,2}, ::Interpolations.WeightedAdjIndex{2,Float64}, ::Interpolations.WeightedAdjIndex{2,Int64}) at ./abstractarray.jl:905
 [9] top-level scope at REBUG:6
 [10] eval(::Module, ::Any) at ./boot.jl:319
 [11] top-level scope at none:0

Note the index for dimension 2 is

Interpolations.WeightedAdjIndex{2,Int64}(0, (0, 1))

which means it will use 0*a[i,0] + 1*a[i,1]. (In case you haven't seen it, explanation is here.)
You get the NaN in #289 when whatever a[i,0] points to happens to be NaN, otherwise the value at a[i,0] is irrelevant since it gets multiplied by 0.

Options for fixing it:

1. check this case at runtime and return valid index objects of different type (WeightedAdjIndex always uses adjacent indices, and hence can't represent a single index for orders higher than Contant). Ideally for a singleton dimension one would return an Int but that would not be inferrable, and since it could happen in any dimension union-splitting would be unlikely to come to the rescue in multidimensional problems.
2. eliminate WeightedAdjIndex and always use WeightedArbIndex, in this case checking at runtime and returning WeightedArbIndex((1,1),(0,1)) for that second dimension and thus computing 0*a[i,1] + 1*a[i,1]. I imagine this would prevent SIMD vectorization, though I can't guarantee it's being used now everywhere it could/should be.
3. automatically filling in Contant() or NoInterp() for singleton dimensions even when the user specifies a single interpolation order for the whole array. The downside is that interpolate would no longer be inferrable.
4. throwing an error, like in this PR. This preserves inferrability at the cost of forcing users to jump through more hoops.

These are the only options I can think of. I suspect 2 and 4 are the most viable, but thoughts from others would be welcome.

[zygmuntszpak]

Thank you for the detailed clarification and explaining the trade-off. I think option 4 is reasonable.