make hessians of scaled interpolants work (#628)

* make hessians of scaled interpolants work, and add a test that would have caught the bug

* make the test deterministic

* Bump version and news

---------

Co-authored-by: Mark Kittisopikul <markkitt@gmail.com>

Author: ajwheeler

NEWS.md

@@ -1,5 +1,10 @@
 For a comprehensive list of changes, see [Releases](https://github.com/JuliaMath/Interpolations.jl/releases).
 
+# v0.16.1
+
+Fixes:
+- Fix `hessian` with ForwardDiff.jl for scaled interpolants ([#628](https://github.com/JuliaMath/Interpolations.jl/pull/628))
+
 # v0.16
 
 Breaking changes:

Project.toml

@@ -1,6 +1,6 @@
 name = "Interpolations"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.16.0"
+version = "0.16.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

ext/InterpolationsForwardDiffExt.jl

@@ -14,10 +14,11 @@ end
 
 # apply partials from arbitrarily nested ForwardDiff.Dual to a value
 # used in maybe_clamp, above
-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))
+function apply_partials(dual_number::D, x::Number) where D <: ForwardDiff.Dual
+    ∂s = ForwardDiff.partials(dual_number)
+    val = ForwardDiff.value(dual_number)
+    D(apply_partials(val, x) , ∂s)
 end
-apply_partials(_::Number, val::Number) = val
+apply_partials(_::Number, x::Number) = x
 
 end

test/hessian.jl

@@ -72,4 +72,21 @@ using Test, Interpolations, LinearAlgebra, ForwardDiff
             end
         end
     end
+
+    @testset "scaled interpolators / convienience constructors (issue 627)" begin
+        N = 11
+        x = range(1.0, 3.0, N)
+        y = range(-3.0, 3.0, N)
+        u = reshape(sin.(1 : (N^2)), N, N)
+
+        for uitp in [
+            Interpolations.linear_interpolation((x, y), u),
+            Interpolations.cubic_spline_interpolation((x, y), u),
+        ]
+            f = x -> uitp(x[1], x[2])
+            x = [π/3, -π/2]
+
+            @test ForwardDiff.hessian(f, x) ≈ Interpolations.hessian(uitp, x[1], x[2])
+        end
+    end
 end