Description
Since ForwardDiff 1.0 introduced a breaking change (requiring that the partials be zero to return true when==ing a real and a dual), linear_interpolation no longer reliably works when the knots are dual numbers. I realize that this is a bit of a weird case that the maintainers of ForwardDiff may justifiably not choose to prioritize, but I want to document the situation for others because it wasn't straightforward to figure out.
Here is a minimal(ish) working example.
ys = zeros(length(xs))
itp = interpolate(ys, BSpline(Linear()))
# these are OK, obviously
itp(2)
itp(2.5)
# the problem is in the scaling
xs = LinRange(
ForwardDiff.Dual(1e-20, 13.7985425771311066e-8),
ForwardDiff.Dual(2e-20, 3.916023275392894e-8),
4)
sitp = scale(itp, xs)
# triggers error
# can also be indirectly called when showing or printing sitp
sitp(xs[end])When running with --check-bounds=yes, this throws an error:
stack trace
BoundsError: attempt to access 4-element Vector{Float64} at index [5]
Stacktrace:
[1] throw_boundserror(A::Vector{Float64}, I::Tuple{Int64})
@ Base ./essentials.jl:14
[2] getindex
@ ./essentials.jl:916 [inlined]
[3] interp_getindex
@ [~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl:308](http://localhost:8888/lab/tree/~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl#line=307) [inlined]
[4] interp_getindex
@ [~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl:302](http://localhost:8888/lab/tree/~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl#line=301) [inlined]
[5] getindex
@ [~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl:286](http://localhost:8888/lab/tree/~/.julia/packages/Interpolations/91PhN/src/Interpolations.jl#line=285) [inlined]
[6] BSplineInterpolation
@ [~/.julia/packages/Interpolations/91PhN/src/b-splines/indexing.jl:8](http://localhost:8888/lab/tree/~/.julia/packages/Interpolations/91PhN/src/b-splines/indexing.jl#line=7) [inlined]
[7] (::ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{LinRange{ForwardDiff.Dual{Nothing, Float64, 1}, Int64}}})(xs::ForwardDiff.Dual{Nothing, Float64, 1})
@ Interpolations [~/.julia/packages/Interpolations/91PhN/src/scaling/scaling.jl:74](http://localhost:8888/lab/tree/~/.julia/packages/Interpolations/91PhN/src/scaling/scaling.jl#line=73)
[8] top-level scope
@ In[118]:18
The offending code is marked with @inbounds.
I first saw this behavior in my test suite (Pkg.test() runs with bounds checking on), but couldn't replicate it an outside environment (where --check-bounds=auto is the default). Instead you get inconsistent behavior (as you would expect when reading outside array bounds).
After some investigation, I think the problem fundamentally boils down to how the positions treats edge cases with non-zero partials. Here's an example.
Interpolations.positions(Linear(), Base.OneTo(4), 4.0)
# returns (3, 1.0)
Interpolations.positions(Linear(), Base.OneTo(4), ForwardDiff.Dual(4.0, 0.0))
# returns (3, Dual{Nothing}(1.0,0.0)), as expected
Interpolations.positions(Linear(), Base.OneTo(4), ForwardDiff.Dual(4.0, 1.0))
# returns (4, Dual{Nothing}(0.0,1.0)), oh no!Because ForwardDiff is a dependency of Interpolations, I think fixing this is as simple as calling ForwardDiff.value in the right place. I'm happy to work up a PR.