Implement and test hessian computation

Author: timholy

src/Interpolations.jl

@@ -187,16 +187,16 @@ Base.to_index(::AbstractInterpolation, x::Number) = x
 #     itp(to_indices(itp, x))
 # end
 function gradient(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
-    gradient(itp, to_indices(itp, x))
+    gradient(itp, to_indices(itp, x)...)
 end
 function gradient!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
-    gradient!(dest, itp, to_indices(itp, x))
+    gradient!(dest, itp, to_indices(itp, x)...)
 end
 function hessian(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
-    hessian(itp, to_indices(itp, x))
+    hessian(itp, to_indices(itp, x)...)
 end
 function hessian!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
-    hessian!(dest, itp, to_indices(itp, x))
+    hessian!(dest, itp, to_indices(itp, x)...)
 end
 
 # @inline function (itp::AbstractInterpolation)(x::Vararg{ExpandedIndexTypes})

src/b-splines/indexing.jl

@@ -17,10 +17,15 @@ end
     @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
     expand_gradient!(dest, itp, x)
 end
+
 @inline function hessian(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
     @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
     expand_hessian(itp, x)
 end
+@inline function hessian!(dest, itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
+    @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
+    expand_hessian(itp, x)
+end
 
 checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{Number,N}) where N =
     checklubounds(lbounds(itp), ubounds(itp), x)
@@ -77,7 +82,7 @@ Calculate the interpolated hessian of `itp` at `x`.
 function expand_hessian(itp::AbstractInterpolation, x::Tuple)
     coefs = coefficients(itp)
     degree = interpdegree(itp)
-    ixs, rxs = expand_indices_resid(degree, axes(coefs), x)
+    ixs, rxs = expand_indices_resid(degree, bounds(itp), x)
     cxs = expand_weights(value_weights, degree, rxs)
     gxs = expand_weights(gradient_weights, degree, rxs)
     hxs = expand_weights(hessian_weights, degree, rxs)
@@ -194,14 +199,76 @@ function expand!(dest, coefs, (vweights, gweights)::Tuple{HasNoInterp{N},HasNoIn
     i = 0
     for d = 1:N
         w = substitute(vweights, d, gweights)
-        w isa Weights || continue   # must have a NoInterp in it
+        w isa Weights || continue   # if this isn't true, it must have a NoInterp in it
         dest[i+=1] = expand(coefs, w, ixs)
     end
     dest
 end
 
-function expand(coefs, (vweights, gweights, hweights)::NTuple{3,Weights{N}}, ixs::Indexes{N}) where N
-    error("not yet implemented")
+# Expansion of the hessian
+# To handle the immutability of SMatrix we build static methods that visit just the entries we need,
+# which due to symmetry is just the upper triangular part
+ntuple_sym(f, ::Val{0}) = ()
+ntuple_sym(f, ::Val{1}) = (f(1,1),)
+ntuple_sym(f, ::Val{2}) = (f(1,1), f(1,2), f(2,2))
+ntuple_sym(f, ::Val{3}) = (f(1,1), f(1,2), f(2,2), f(1,3), f(2,3), f(3,3))
+ntuple_sym(f, ::Val{4}) = (f(1,1), f(1,2), f(2,2), f(1,3), f(2,3), f(3,3), f(1,4), f(2,4), f(3,4), f(4,4))
+@inline function ntuple_sym(f, ::Val{N}) where N
+    (ntuple_sym(f, Val(N-1))..., ntuple(i->f(i,N), Val(N))...)
+end
+
+sym2dense(t::Tuple{})              = t
+sym2dense(t::NTuple{1,T}) where T  = t
+sym2dense(t::NTuple{3,T}) where T  = (t[1], t[2], t[2], t[3])
+sym2dense(t::NTuple{6,T}) where T  = (t[1], t[2], t[4], t[2], t[3], t[5], t[4], t[5], t[6])
+sym2dense(t::NTuple{10,T}) where T = (t[1], t[2], t[4], t[7], t[2], t[3], t[5], t[8], t[4], t[5], t[6], t[9], t[7], t[8], t[9], t[10])
+function sym2dense(t::NTuple{L,T}) where {L,T}
+    # Warning: non-inferrable unless we make this @generated.
+    # Above 4 dims one might anyway prefer an Array, and use hessian!
+    N = ceil(Int, sqrt(2*L))
+    @assert (N*(N+1))÷2 == L
+    a = Vector{T}(undef, N*N)
+    idx = 0
+    for j = 1:N, i=1:N
+        iu, ju = ifelse(i>=j, (j, i), (i, j))  # index in the upper triangular
+        k = (ju*(ju+1))÷2 + iu
+        a[idx+=1] = t[k]
+    end
+    tuple(a...)
+end
+
+squarematrix(t::NTuple{1,T}) where T  = SMatrix{1,1,T}(t)
+squarematrix(t::NTuple{4,T}) where T  = SMatrix{2,2,T}(t)
+squarematrix(t::NTuple{9,T}) where T  = SMatrix{3,3,T}(t)
+squarematrix(t::NTuple{16,T}) where T = SMatrix{4,4,T}(t)
+function squarematrix(t::NTuple{L,T}) where {L,T}
+    # Warning: non-inferrable unless we make this @generated.
+    # Above 4 dims one might anyway prefer an Array, and use hessian!
+    N = floor(Int, sqrt(L))
+    @assert N*N == L
+    SMatrix{N,N,T}(t)
+end
+
+cumweights(w) = _cumweights(0, w...)
+_cumweights(c, w1, w...) = (c+1, _cumweights(c+1, w...)...)
+_cumweights(c, ::NoInterp, w...) = (c, _cumweights(c, w...)...)
+_cumweights(c) = ()
+
+function expand(coefs, (vweights, gweights, hweights)::NTuple{3,HasNoInterp{N}}, ixs::Indexes{N}) where N
+    coefs = ntuple_sym((i,j)->expand(coefs, substitute(vweights, i, j, gweights, hweights), ixs), Val(N))
+    squarematrix(sym2dense(skip_nointerp(coefs...)))
+end
+
+function expand!(dest, coefs, (vweights, gweights, hweights)::NTuple{3,HasNoInterp{N}}, ixs::Indexes{N}) where N
+    # The Hessian is nominally N × N, but if there are K NoInterp dims then it's N-K × N-K
+    indlookup = cumweights(hweights)   # for d in 1:N, indlookup[d] returns the appropriate index in 1:N-K
+    for d2 = 1:N, d1 = 1:d2
+        w = substitute(vweights, d1, d2, gweights, hweights)
+        w isa Weights || continue   # if this isn't true, it must have a NoInterp in it
+        i, j = indlookup[d1], indlookup[d2]
+        dest[i, j] = dest[j, i] = expand(coefs, w, ixs)
+    end
+    dest
 end
 
 function expand_indices_resid(degree, bounds, x)
@@ -255,7 +322,7 @@ roundbounds(x, bounds::Tuple{Integer,Integer}) = round(x)
 roundbounds(x, (l, u)) = ifelse(x == l, ceil(l), ifelse(x == u, floor(u), round(x)))
 
 floorbounds(x, bounds::Tuple{Integer,Integer}) = floor(x)
-function floorbounds(x, (l, u))
+function floorbounds(x, (l, u)::Tuple{Real,Real})
     ceill = ceil(l)
     ifelse(l <= x <= ceill, ceill, floor(x))
 end

src/utils.jl

@@ -14,13 +14,19 @@ fast_trunc(::Type{Int}, x::Rational) = x.num ÷ x.den
 iextract(f::Flag, d) = f
 iextract(t::Tuple, d) = t[d]
 
+# Substitution for gradient components
 function substitute(default::NTuple{N,Any}, d::Integer, subst::NTuple{N,Any}) where N
     ntuple(i->ifelse(i==d, subst[i], default[i]), Val(N))
 end
 function substitute(default::NTuple{N,Any}, d::Integer, val) where N
     ntuple(i->ifelse(i==d, val, default[i]), Val(N))
 end
 
+# Substitution for hessian components
+function substitute(default::NTuple{N,Any}, d1::Integer, d2::Integer, subst1::NTuple{N,Any}, subst2::NTuple{N,Any}) where N
+    ntuple(i->ifelse(i==d1==d2, subst2[i], ifelse(i==d1, subst1[i], ifelse(i==d2, subst1[i], default[i]))), Val(N))
+end
+
 @inline skip_nointerp(x, rest...) = (x, skip_nointerp(rest...)...)
 @inline skip_nointerp(::NoInterp, rest...) = skip_nointerp(rest...)
 skip_nointerp() = ()

test/InterpolationTestUtils.jl

@@ -1,10 +1,11 @@
 module InterpolationTestUtils
 
-using Test, Interpolations
+using Test, Interpolations, ForwardDiff, StaticArrays
 using Interpolations: degree, itpflag, bounds, lbounds, ubounds
 using Interpolations: substitute
 
-export check_axes, check_inbounds_values, check_oob, can_eval_near_boundaries
+export check_axes, check_inbounds_values, check_oob, can_eval_near_boundaries,
+       check_gradient, check_hessian
 export MyPair
 
 const failstore = Ref{Any}(nothing)   # stash the inputs to failing tests here
@@ -105,7 +106,34 @@ function can_eval_near_boundaries(itp::AbstractInterpolation)
     end
 end
 
-# Used for multi-valued tests
+# Generate a grid of points [1.0, 1.3333, 1.6667, 2.0, 2.3333, ...] along each coordinate
+thirds(axs) = Iterators.product(_thirds(axs...)...)
+
+_thirds(a, axs...) =
+    (sort(Float64[a; (first(a):last(a)-1) .+ 1/3; (first(a)+1:last(a)) .- 1/3]), _thirds(axs...)...)
+_thirds() = ()
+
+function check_gradient(itp::AbstractInterpolation, gtmp)
+    val(x) = itp(Tuple(x)...)
+    g!(gstore, x) = ForwardDiff.gradient!(gstore, val, x)
+    gtmp2 = similar(gtmp)
+    for i in thirds(axes(itp))
+        @test Interpolations.gradient(itp, i...) ≈ g!(gtmp, SVector(i))
+        @test Interpolations.gradient!(gtmp2, itp, i...) ≈ gtmp
+    end
+end
+
+function check_hessian(itp::AbstractInterpolation, htmp)
+    val(x) = itp(Tuple(x)...)
+    h!(hstore, x) = ForwardDiff.hessian!(hstore, val, x)
+    htmp2 = similar(htmp)
+    for i in thirds(axes(itp))
+        @test Interpolations.hessian(itp, i...) ≈ h!(htmp, SVector(i))
+        @test Interpolations.hessian!(htmp2, itp, i...) ≈ htmp
+    end
+end
+
+## A type used for multi-valued tests
 import Base: +, -, *, /, ≈
 
 struct MyPair{T}

test/gradient.jl

@@ -3,18 +3,40 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
 @testset "Gradients" begin
     nx = 10
     f1(x) = sin((x-3)*2pi/(nx-1) - 1)
-    g1(x) = 2pi/(nx-1) * cos((x-3)*2pi/(nx-1) - 1)
+    g1gt(x) = 2pi/(nx-1) * cos((x-3)*2pi/(nx-1) - 1)
+    A1 = Float64[f1(x) for x in 1:nx]
+    g1 = Array{Float64}(undef, 1)
+    A2 = rand(Float64, nx, nx) * 100
+    g2 = Array{Float64}(undef, 2)
+
+    for (A, g) in ((A1, g1), (A2, g2))
+        # Gradient of Constant should always be 0
+        itp = interpolate(A, BSpline(Constant()), OnGrid())
+        for x in InterpolationTestUtils.thirds(axes(A))
+            @test all(iszero, Interpolations.gradient(itp, x...))
+            @test all(iszero, Interpolations.gradient!(g, itp, x...))
+        end
 
-    # Gradient of Constant should always be 0
-    itp1 = interpolate(Float64[f1(x) for x in 1:nx],
-                       BSpline(Constant()), OnGrid())
+        for GT in (OnGrid, OnCell)
+            itp = interpolate(A, BSpline(Linear()), GT())
+            check_gradient(itp, g)
+            I = first(eachindex(itp))
+            @test Interpolations.gradient(itp, I) == Interpolations.gradient(itp, Tuple(I)...)
+        end
 
-    g = Array{Float64}(undef, 1)
+        for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
+            itp = interpolate(A, BSpline(Quadratic(BC())), GT())
+            check_gradient(itp, g)
+            I = first(eachindex(itp))
+            @test Interpolations.gradient(itp, I) == Interpolations.gradient(itp, Tuple(I)...)
+        end
 
-    for x in 1:nx
-        @test Interpolations.gradient(itp1, x)[1] == 0
-        @test Interpolations.gradient!(g, itp1, x)[1] == 0
-        @test g[1] == 0
+        for BC in (Line, Flat, Free, Periodic), GT in (OnGrid, OnCell)
+            itp = interpolate(A, BSpline(Cubic(BC())), GT())
+            check_gradient(itp, g)
+            I = first(eachindex(itp))
+            @test Interpolations.gradient(itp, I) == Interpolations.gradient(itp, Tuple(I)...)
+        end
     end
 
     # Since Linear is OnGrid in the domain, check the gradients between grid points
@@ -24,9 +46,9 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
     #             Gridded(Linear()))
     for itp in (itp1, )#itp2)
         for x in 2.5:nx-1.5
-            @test ≈(g1(x),(Interpolations.gradient(itp,x))[1],atol=abs(0.1 * g1(x)))
-            @test ≈(g1(x),(Interpolations.gradient!(g,itp,x))[1],atol=abs(0.1 * g1(x)))
-            @test ≈(g1(x),g[1],atol=abs(0.1 * g1(x)))
+            @test ≈(g1gt(x),(Interpolations.gradient(itp,x))[1],atol=abs(0.1 * g1gt(x)))
+            @test ≈(g1gt(x),(Interpolations.gradient!(g1,itp,x))[1],atol=abs(0.1 * g1gt(x)))
+            @test ≈(g1gt(x),g1[1],atol=abs(0.1 * g1gt(x)))
         end
 
         for i = 1:10
@@ -52,9 +74,9 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
     itp1 = interpolate(Float64[f1(x) for x in 1:nx-1],
                        BSpline(Quadratic(Periodic())), OnCell())
     for x in 2:nx-1
-        @test ≈(g1(x),(Interpolations.gradient(itp1,x))[1],atol=abs(0.05 * g1(x)))
-        @test ≈(g1(x),(Interpolations.gradient!(g,itp1,x))[1],atol=abs(0.05 * g1(x)))
-        @test ≈(g1(x),g[1],atol=abs(0.1 * g1(x)))
+        @test ≈(g1gt(x),(Interpolations.gradient(itp1,x))[1],atol=abs(0.05 * g1gt(x)))
+        @test ≈(g1gt(x),(Interpolations.gradient!(g1,itp1,x))[1],atol=abs(0.05 * g1gt(x)))
+        @test ≈(g1gt(x),g1[1],atol=abs(0.1 * g1gt(x)))
     end
 
     for i = 1:10

test/hessian.jl

@@ -0,0 +1,44 @@
+using Test, Interpolations, LinearAlgebra
+
+@testset "Hessians" begin
+    nx = 5
+    k = 2pi/(nx-1)
+    f1(x) = sin(k*(x-3) - 1)
+    A1 = Float64[f1(x) for x in 1:nx]
+    h1 = Array{Float64}(undef, 1, 1)
+    A2 = rand(Float64, nx, nx) * 100
+    h2 = Array{Float64}(undef, 2, 2)
+
+    for (A, h) in ((A1, h1), (A2, h2))
+        for GT in (OnGrid, OnCell)
+            for itp in (interpolate(A, BSpline(Constant()), GT()),
+                        interpolate(A, BSpline(Linear()), GT()))
+                if ndims(A) == 1
+                    # Hessian of Constant and Linear should always be 0 in 1d
+                    for x in InterpolationTestUtils.thirds(axes(A))
+                        @test all(iszero, Interpolations.hessian(itp, x...))
+                        @test all(iszero, Interpolations.hessian!(h, itp, x...))
+                    end
+                else
+                    for x in InterpolationTestUtils.thirds(axes(A))
+                        check_hessian(itp, h)
+                    end
+                end
+            end
+        end
+
+        for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
+            itp = interpolate(A, BSpline(Quadratic(BC())), GT())
+            check_hessian(itp, h)
+            I = first(eachindex(itp))
+            @test Interpolations.hessian(itp, I) == Interpolations.hessian(itp, Tuple(I)...)
+        end
+
+        for BC in (Line, Flat, Free, Periodic), GT in (OnGrid, OnCell)
+            itp = interpolate(A, BSpline(Cubic(BC())), GT())
+            check_hessian(itp, h)
+        end
+    end
+
+    # TODO: mixed interpolation (see gradient.jl)
+end

test/runtests.jl

@@ -23,6 +23,8 @@ using Interpolations
 
     # # test gradient evaluation
     include("gradient.jl")
+    # # test hessian evaluation
+    include("hessian.jl")
 
     # # gridded interpolation tests
     # include("gridded/runtests.jl")