git issue 335 test passing; concentrate all direct calls to ForwardDiff to one file

Author: ajwheeler

src/Interpolations.jl

@@ -447,9 +447,29 @@ maybe_clamp(itp, xs) = maybe_clamp(BoundsCheckStyle(itp), itp, xs)
 maybe_clamp(::NeedsCheck, itp, xs) = map(clamp, xs, lbounds(itp), ubounds(itp))
 maybe_clamp(::CheckWillPass, itp, xs) = xs
 
+using ForwardDiff
+# TODO
+function maybe_clamp(::NeedsCheck, itp, xs::Tuple{Vararg{ForwardDiff.Dual}})
+    xs_values = just_dual_value.(xs)
+    clamped_vals = maybe_clamp(NeedsCheck(), itp, xs_values)
+    apply_partials.(xs, clamped_vals)
+end
+
 Base.hash(x::AbstractInterpolation, h::UInt) = Base.hash_uint(3h - objectid(x))
 Base.hash(x::AbstractExtrapolation, h::UInt) = Base.hash_uint(3h - objectid(x))
 
+# TODO use this, and define a method in a ForwardDiff package extension
+# stip off arbitrary layers of ForwardDiff.Dual, returning the innermost value
+just_dual_value(x::Number) = x
+just_dual_value(x::ForwardDiff.Dual) = just_dual_value(ForwardDiff.value(x))
+
+# apply partials from arbitrarily nested ForwardDiff.Dual to a value
+function apply_partials(x_dual::D, val::Number) where D <: ForwardDiff.Dual
+    ∂s = ForwardDiff.partials(x_dual)
+    apply_partials(ForwardDiff.value(x_dual), D(val, ∂s))
+end
+apply_partials(x_dual::Number, val::Number) = val
+
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")

src/b-splines/constant.jl

@@ -67,35 +67,35 @@ to `A[ceil(Int,x)]` without scaling.
 Constant
 
 function positions(c::Constant{Previous}, ax, x)  # discontinuity occurs at integer locations
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     xm = floorbounds(x_value, ax)
     δx = x - xm
     fast_trunc(Int, xm), δx
 end
 function positions(c::Constant{Next}, ax, x)  # discontinuity occurs at integer locations
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     xm = ceilbounds(x_value, ax)
     δx = x - xm
     fast_trunc(Int, xm), δx
 end
 function positions(c::Constant{Nearest}, ax, x)  # discontinuity occurs at half-integer locations
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     xm = roundbounds(x_value, ax)
     δx = x - xm
     i = fast_trunc(Int, xm)
     i, δx
 end
 
 function positions(c::Constant{Previous,Periodic{OnCell}}, ax, x)
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     # We do not use floorbounds because we do not want to add a half at
     # the lowerbound to round up.
     xm = floor(x_value)
     δx = x - xm
     modrange(fast_trunc(Int, xm), ax), δx
 end
 function positions(c::Constant{Next,Periodic{OnCell}}, ax, x)  # discontinuity occurs at integer locations
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     xm = ceilbounds(x_value, ax)
     δx = x - xm
     modrange(fast_trunc(Int, xm), ax), δx

src/b-splines/indexing.jl

@@ -3,7 +3,7 @@
 itpinfo(itp) = (tcollect(itpflag, itp), axes(itp))
 
 @inline function (itp::BSplineInterpolation{T,N})(x::Vararg{Number,N}) where {T,N}
-    @boundscheck (checkbounds(Bool, itp, ForwardDiff.value.(ForwardDiff.value.(x))...) || Base.throw_boundserror(itp, x))
+    @boundscheck (checkbounds(Bool, itp, just_dual_value.(x)...) || Base.throw_boundserror(itp, x))
     wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
     InterpGetindex(itp)[wis...]
 end

src/b-splines/linear.jl

@@ -1,5 +1,3 @@
-using ForwardDiff # TODO
-
 struct Linear{BC<:Union{Throw{OnGrid},Periodic{OnCell}}} <: DegreeBC{1}
     bc::BC
     function Linear{BC}(bc::BC=BC()) where BC<:Union{Throw{OnGrid},Periodic{OnCell}}
@@ -43,7 +41,7 @@ a piecewise linear function connecting each pair of neighboring data points.
 Linear
 
 function positions(deg::Linear, ax::AbstractUnitRange{<:Integer}, x)
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     f = floor(x_value)
     # When x == last(ax) we want to use the x-1, x pair
     f = ifelse(x_value == last(ax), f - oneunit(f), f)

src/monotonic/monotonic.jl

@@ -207,7 +207,7 @@ function interpolate(
 end
 
 function (itp::MonotonicInterpolation)(x::Number)
-    x_value = ForwardDiff.value(ForwardDiff.value(x))
+    x_value = just_dual_value.(x)
     @boundscheck (checkbounds(Bool, itp, x_value) || Base.throw_boundserror(itp, (x_value,)))
     k = searchsortedfirst(itp.knots, x_value)
     if k > 1

src/scaling/scaling.jl

@@ -76,7 +76,8 @@ ubound(ax::AbstractRange, ::DegreeBC, ::OnGrid) = last(ax)
 
 # For (), we scale the evaluation point
 @propagate_inbounds function (sitp::ScaledInterpolation{T,N})(xs::Vararg{Number,N}) where {T,N}
-    @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs))
+    xs_values = just_dual_value.(xs)
+    @boundscheck (checkbounds(Bool, sitp, xs_values...) || Base.throw_boundserror(sitp, xs_values))
     xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs))
     @inbounds sitp.itp(xl...)
 end