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

This reworks bounds-checking around BoundsCheckStyle trait.

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/gridded.jl

@@ -69,6 +69,8 @@ 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")

src/scaling/scaling.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

test/extrapolation/runtests.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

test/scaling/withextrap.jl

@@ -5,15 +5,15 @@ using Interpolations, Test
     xs = range(-5, stop=5, length=10)
     ys = map(sin, xs)
 
-    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp) where {T,N,IT}
+    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp, ::OnGrid) where {T,N,IT}
         for x in xs
             @test ≈(sut(x),sin(x),atol=sqrt(eps(sin(x))))
         end
         @test sut(-5) == sut(-5.1) == sut(-15.8) == sut(-Inf) == itp(1)
         @test sut(5) == sut(5.1) == sut(15.8) == sut(Inf) == itp[end]
     end
 
-    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp) where {T,N,IT}
+    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp, ::OnCell) where {T,N,IT}
         halfcell = (xs[2] - xs[1]) / 2
         itps, axs = Interpolations.itpinfo(itp)
         for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
@@ -22,6 +22,8 @@ using Interpolations, Test
         end
     end
 
+    run_tests(sut, itp) = run_tests(sut, itp, Interpolations.degree(Interpolations.itpflag(itp)).bc.gt)
+
     for GT in (OnGrid, OnCell)
         itp = interpolate(ys, BSpline(Quadratic(Flat(GT()))))