Merge pull request #245 from JuliaMath/teh/fix_#244

Avoid numerical-precision failures of bounds checking (fixes #244)

Author: timholy

src/Interpolations.jl

@@ -129,6 +129,18 @@ count_interp_dims(it::Type{IT}, n) where IT<:Tuple{Vararg{InterpolationType,N}}
     _count_interp_dims(c + count_interp_dims(IT1), args...)
 _count_interp_dims(c) = c
 
+"""
+    BoundsCheckStyle(itp)
+
+A trait to determine dispatch of bounds-checking for `itp`.
+Can return `NeedsCheck()`, in which case bounds-checking is performed, or `CheckWillPass()`
+in which case the check will return `true`.
+"""
+abstract type BoundsCheckStyle end
+struct NeedsCheck <: BoundsCheckStyle end
+struct CheckWillPass <: BoundsCheckStyle end
+
+BoundsCheckStyle(itp) = NeedsCheck()
 
 """
     wi = WeightedIndex(indexes, weights)
@@ -386,6 +398,23 @@ import Base: getindex
     itp(i)
 end
 
+@inline checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes,N}) where N =
+    _checkbounds(BoundsCheckStyle(itp), itp, x)
+
+_checkbounds(::CheckWillPass, itp, x) = true
+_checkbounds(::NeedsCheck, itp, x) = checklubounds(lbounds(itp), ubounds(itp), x)
+
+checklubounds(ls, us, xs) = _checklubounds(true, ls, us, xs)
+_checklubounds(tf::Bool, ls, us, xs::Tuple{Number, Vararg{Any}}) =
+    _checklubounds(tf & (ls[1] <= xs[1] <= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
+_checklubounds(tf::Bool, ls, us, xs::Tuple{AbstractVector, Vararg{Any}}) =
+    _checklubounds(tf & all(ls[1] .<= xs[1] .<= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
+_checklubounds(tf::Bool, ::Tuple{}, ::Tuple{}, ::Tuple{}) = tf
+
+maybe_clamp(itp, xs) = maybe_clamp(BoundsCheckStyle(itp), itp, xs)
+maybe_clamp(::NeedsCheck, itp, xs) = clamp.(xs, lbounds(itp), ubounds(itp))
+maybe_clamp(::CheckWillPass, itp, xs) = xs
+
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")

src/b-splines/indexing.jl

@@ -22,7 +22,7 @@ end
     reshape(ret, shape(wis...))
 end
 
-@inline function gradient(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
+@propagate_inbounds function gradient(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
     @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
     wis = weightedindexes((value_weights, gradient_weights), itpinfo(itp)..., x)
     SVector(map(inds->itp.coefs[inds...], wis))
@@ -31,7 +31,7 @@ end
     dest .= gradient(itp, x...)
 end
 
-@inline function hessian(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
+@propagate_inbounds function hessian(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
     @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
     wis = weightedindexes((value_weights, gradient_weights, hessian_weights), itpinfo(itp)..., x)
     symmatrix(map(inds->itp.coefs[inds...], wis))
@@ -40,16 +40,6 @@ end
     dest .= hessian(itp, x...)
 end
 
-checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes,N}) where N =
-    checklubounds(lbounds(itp), ubounds(itp), x)
-
-checklubounds(ls, us, xs) = _checklubounds(true, ls, us, xs)
-_checklubounds(tf::Bool, ls, us, xs::Tuple{Number, Vararg{Any}}) =
-    _checklubounds(tf & (ls[1] <= xs[1] <= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
-_checklubounds(tf::Bool, ls, us, xs::Tuple{AbstractVector, Vararg{Any}}) =
-    _checklubounds(tf & all(ls[1] .<= xs[1] .<= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
-_checklubounds(tf::Bool, ::Tuple{}, ::Tuple{}, ::Tuple{}) = tf
-
 # Leftovers from AbstractInterpolation
 @inline function (itp::BSplineInterpolation)(x::Vararg{UnexpandedIndexTypes})
     itp(to_indices(itp, x)...)

src/extrapolation/extrapolation.jl

@@ -15,6 +15,8 @@ const ExtrapDimSpec = Union{BoundaryCondition,Tuple{Vararg{Union{BoundaryConditi
 etptype(::Extrapolation{T,N,ITPT,IT,ET}) where {T,N,ITPT,IT,ET} = ET
 etpflag(etp::Extrapolation{T,N,ITPT,IT,ET}) where {T,N,ITPT,IT,ET} = etp.et
 
+BoundsCheckStyle(etp::AbstractExtrapolation) = CheckWillPass()
+
 """
 `extrapolate(itp, scheme)` adds extrapolation behavior to an interpolation object, according to the provided scheme.
 
@@ -35,11 +37,12 @@ extrapolate(itp::AbstractInterpolation{T,N,IT}, et::ET) where {T,N,IT,ET<:Extrap
 
 count_interp_dims(::Type{<:Extrapolation{T,N,ITPT}}, n) where {T,N,ITPT} = count_interp_dims(ITPT, n)
 
-@inline function (etp::Extrapolation{T,N})(x::Vararg{Number,N}) where {T,N}
+@propagate_inbounds function (etp::Extrapolation{T,N})(x::Vararg{Number,N}) where {T,N}
     itp = parent(etp)
     eflag = etpflag(etp)
     xs = inbounds_position(eflag, bounds(itp), x, etp, x)
-    extrapolate_value(eflag, skip_flagged_nointerp(itp, x), skip_flagged_nointerp(itp, xs), Tuple(gradient(itp, xs...)), itp(xs...))
+    g = @inbounds gradient(itp, xs...)
+    extrapolate_value(eflag, skip_flagged_nointerp(itp, x), skip_flagged_nointerp(itp, xs), Tuple(g), @inbounds(itp(xs...)))
 end
 @inline function (etp::Extrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
     itp = parent(etp)
@@ -58,7 +61,7 @@ end
     else
         eflag = tcollect(etpflag, etp)
         xs = inbounds_position(eflag, bounds(itp), x, etp, x)
-        g = gradient(itp, xs...)
+        g = @inbounds gradient(itp, xs...)
         skipni = t->skip_flagged_nointerp(itp, t)
         SVector(extrapolate_gradient.(skipni(eflag), skipni(x), skipni(xs), Tuple(g)))
     end
@@ -127,7 +130,7 @@ end
 periodic(y, l, u) = mod(y-l, u-l) + l
 
 
-function extrapolate_value(eflag, x, xs, g, val)
+@inline function extrapolate_value(eflag, x, xs, g, val)
     val = extrapolate_axis(getfirst(eflag), x[1], xs[1], g[1], val)
     extrapolate_value(getrest(eflag), Base.tail(x), Base.tail(xs), Base.tail(g), val)
 end

src/gridded/constant.jl

@@ -69,7 +69,7 @@ end
 
 lbounds(itp::GriddedInterpolation) = first.(itp.knots)
 ubounds(itp::GriddedInterpolation) = last.(itp.knots)
+lbound(ax::AbstractVector, gr::Gridded)   = lbound(ax, degree(gr))
+ubound(ax::AbstractVector, gr::Gridded)   = ubound(ax, degree(gr))
 
-include("constant.jl")
-include("linear.jl")
 include("indexing.jl")

src/gridded/gridded.jl

@@ -9,6 +9,7 @@ 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)
+BoundsCheckStyle(sitp::ScaledInterpolation) = BoundsCheckStyle(sitp.itp)
 
 """
 `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.
@@ -58,9 +59,10 @@ lbound(ax::AbstractRange, ::DegreeBC, ::OnGrid) = first(ax)
 ubound(ax::AbstractRange, ::DegreeBC, ::OnGrid) = last(ax)
 
 # For (), we scale the evaluation point
-function (sitp::ScaledInterpolation{T,N})(xs::Vararg{Number,N}) where {T,N}
-    xl = coordslookup(itpflag(sitp.itp), sitp.ranges, xs)
-    sitp.itp(xl...)
+@propagate_inbounds function (sitp::ScaledInterpolation{T,N})(xs::Vararg{Number,N}) where {T,N}
+    @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs))
+    xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs))
+    @inbounds sitp.itp(xl...)
 end
 @inline function (sitp::ScaledInterpolation)(x::Vararg{UnexpandedIndexTypes})
     xis = to_indices(sitp, x)
@@ -71,7 +73,6 @@ end
 (sitp::ScaledInterpolation{T,1}, x::Number, y::Int) where {T} = y == 1 ? sitp(x) : Base.throw_boundserror(sitp, (x, y))
 
 @inline function (itp::ScaledInterpolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
-    # @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
     [itp(i...) for i in Iterators.product(x...)]
 end
 
@@ -96,8 +97,9 @@ basetype(::Type{ScaledInterpolation{T,N,ITPT,IT,RT}}) where {T,N,ITPT,IT,RT} = I
 basetype(sitp::ScaledInterpolation) = basetype(typeof(sitp))
 
 
-function gradient(sitp::ScaledInterpolation{T,N}, xs::Vararg{Number,N}) where {T,N}
-    xl = coordslookup(itpflag(sitp.itp), sitp.ranges, xs)
+@propagate_inbounds function gradient(sitp::ScaledInterpolation{T,N}, xs::Vararg{Number,N}) where {T,N}
+    @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs))
+    xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs))
     g = gradient(sitp.itp, xl...)
     SVector(rescale_gradient_components(itpflag(sitp.itp), sitp.ranges, Tuple(g)))
 end

src/gridded/linear.jl

@@ -1,35 +1,35 @@
-module ExtrapNon1
-
 using Test, Interpolations, OffsetArrays
 
-f(x) = sin((x-3)*2pi/9 - 1)
-xinds = -3:6
-A = OffsetArray(Float64[f(x) for x in xinds], xinds)
-itpg = interpolate(A, BSpline(Linear()))
+@testset "ExtrapNon1" begin
 
-schemes = (
-    Flat,
-    Line,
-    Reflect,
-    Periodic
-)
+    f(x) = sin((x-3)*2pi/9 - 1)
+    xinds = -3:6
+    A = OffsetArray(Float64[f(x) for x in xinds], xinds)
+    itpg = interpolate(A, BSpline(Linear()))
 
-for etp in map(E -> extrapolate(itpg, E()), schemes), x in xinds
-    @test parent(etp) === itpg
-    @test @inferred(getindex(etp, x)) ≈ A[x]
-end
+    schemes = (
+        Flat,
+        Line,
+        Reflect,
+        Periodic
+    )
 
-g(y) = (y/100)^3
-yinds = 2:5
-A = OffsetArray(Float64[f(x)*g(y) for x in xinds, y in yinds], xinds, yinds)
-itp2 = interpolate(A, BSpline(Linear()))
+    for etp in map(E -> extrapolate(itpg, E()), schemes), x in xinds
+        @test parent(etp) === itpg
+        @test @inferred(getindex(etp, x)) ≈ A[x]
+    end
 
-for (etp2,E) in map(E -> (extrapolate(itp2, E()), E), schemes)
-    @test parent(etp2) === itp2
-    E == Periodic && continue  # g isn't periodic
-    for y in yinds, x in xinds
-        @test @inferred(getindex(etp2, x, y)) ≈ A[x, y]
+    g(y) = (y/100)^3
+    yinds = 2:5
+    A = OffsetArray(Float64[f(x)*g(y) for x in xinds, y in yinds], xinds, yinds)
+    itp2 = interpolate(A, BSpline(Linear()))
+
+    for (etp2,E) in map(E -> (extrapolate(itp2, E()), E), schemes)
+        @test parent(etp2) === itp2
+        E == Periodic && continue  # g isn't periodic
+        for y in yinds, x in xinds
+            @test @inferred(getindex(etp2, x, y)) ≈ A[x, y]
+        end
     end
-end
 
 end

src/scaling/scaling.jl

@@ -95,6 +95,12 @@ using Test
     targeterr = ArgumentError("cannot create a filled extrapolation with a type Line, use a value of this type instead (e.g., Line())")
     @test_throws targeterr extrapolate(itp, Line)
 
+    # Issue #244
+    xs = range(1e-2, stop = 8.3, length = 3)
+    ys = sort(rand(3))
+    itp = LinearInterpolation(xs, ys, extrapolation_bc = Flat())
+    @test itp(8.3) ≈ ys[end]
+
     include("type-stability.jl")
     include("non1.jl")
 end