Count gradient dimensions more accurately (fixes #156) (#189)

Author: timholy

src/extrapolation/extrapolation.jl

@@ -30,6 +30,8 @@ Finally, you can specify different extrapolation behavior in different direction
 extrapolate(itp::AbstractInterpolation{T,N,IT,GT}, ::ET) where {T,N,IT,GT,ET<:ExtrapDimSpec} =
     Extrapolation{T,N,typeof(itp),IT,GT,ET}(itp)
 
+count_interp_dims(::Type{<:Extrapolation{T,N,ITPT}}, n) where {T,N,ITPT} = count_interp_dims(ITPT, n)
+
 include("throw.jl")
 include("flat.jl")
 include("linear.jl")

src/scaling/scaling.jl

@@ -15,6 +15,7 @@ end
 end
 
 Base.parent(A::ScaledInterpolation) = A.itp
+count_interp_dims(::Type{<:ScaledInterpolation{T,N,ITPT}}, n) where {T,N,ITPT} = count_interp_dims(ITPT, n)
 
 """
 `scale(itp, xs, ys, ...)` scales an existing interpolation object to allow for indexing using other coordinate axes than unit ranges, by wrapping the interpolation object and transforming the indices from the provided axes onto unit ranges upon indexing.

test/extrapolation/runtests.jl

@@ -91,6 +91,14 @@ for E in [0,Flat(),Linear(),Periodic(),Reflect()]
     @test (@inferred(getindex(extrapolate(interpolate([0,0],BSpline(Linear()),OnGrid()),E),[1.2, 1.8, 3.1]))) == [0,0,0]
 end
 
+# Issue #156
+F     = *(collect(1.:10.), collect(1:4)')
+itp   = interpolate(F, (BSpline(Linear()), NoInterp()), OnGrid());
+itps   = scale(itp, 1:10, 1:4)
+itpe   = extrapolate(itps, (Linear(), Interpolations.Throw()))
+@test itpe[10.1, 1] ≈ 10.1
+@test_throws BoundsError itpe[9.9, 0]
+
 end
 
 include("type-stability.jl")