Implement Hessian

Author: Tomas Lycken

src/Interpolations.jl

@@ -7,6 +7,10 @@ export
     scale,
 
     gradient!,
+    gradient1,
+    hessian!,
+    hessian,
+    hessian1,
 
     AbstractInterpolation,
     AbstractExtrapolation,
@@ -74,13 +78,6 @@ gridtype{ITP<:AbstractInterpolation}(::Type{ITP}) = gridtype(super(ITP))
 gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
 count_interp_dims{T<:AbstractInterpolation}(::Type{T}, N) = N
 
-@generated function gradient{T,N}(itp::AbstractInterpolation{T,N}, xs...)
-    n = count_interp_dims(itp, N)
-    Tg = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
-    xargs = [:(xs[$d]) for d in 1:length(xs)]
-    :(gradient!(Array($Tg,$n), itp, $(xargs...)))
-end
-
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")

src/b-splines/constant.jl

@@ -15,6 +15,11 @@ function gradient_coefficients(::Type{BSpline{Constant}}, d)
     :($sym = 0)
 end
 
+function hessian_coefficients(::Type{BSpline{Constant}}, d)
+    sym = symbol("c_",d)
+    :($sym = 0)
+end
+
 function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Constant}}, ::Type{IT}, N::Integer, offsets...)
     if (length(offsets) < N)
         d = length(offsets)+1

src/b-splines/cubic.jl

@@ -108,6 +108,29 @@ function gradient_coefficients{C<:Cubic}(::Type{BSpline{C}}, d)
     end
 end
 
+"""
+In `hessian_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
+
+    cm_d  = p''(fx_d)
+    c_d   = q''(fx_d)
+    cp_d  = q''(1-fx_d)
+    cpp_d = p''(1-fx_d)
+
+where `p` and `q` are defined in the docstring entry for `Cubic`, and
+`fx_d` in the docstring entry for `define_indices_d`.
+"""
+function hessian_coefficients{C<:Cubic}(::Type{BSpline{C}}, d)
+    symm, sym, symp, sympp = symbol("cm_",d), symbol("c_",d), symbol("cp_",d), symbol("cpp_",d)
+    symfx = symbol("fx_",d)
+    quote
+        $symm  = 1 - $symfx
+        $sym   = 3 * $symfx - 2
+        $symp  = 1 - 3 * $symfx
+        $sympp = $symfx
+    end
+end
+
 function index_gen{C<:Cubic,IT<:DimSpec{BSpline}}(::Type{BSpline{C}}, ::Type{IT}, N::Integer, offsets...)
     if length(offsets) < N
         d = length(offsets)+1

src/b-splines/indexing.jl

@@ -13,6 +13,16 @@ function gradient_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim)
                         coefficients(iextract(IT, d), N, d) for d = 1:N]
     Expr(:block, exs...)
 end
+function hessian_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim1, dim2)
+    exs = if dim1 == dim2
+        Expr[d==dim1==dim2 ? hessian_coefficients(iextract(IT, dim), d) :
+                             coefficients(iextract(IT, d), N, d) for d in 1:N]
+    else
+        Expr[d==dim1 || d==dim2 ? gradient_coefficients(iextract(IT, dim), d) :
+                                  coefficients(iextract(IT, d), N, d) for d in 1:N]
+    end
+    Expr(:block, exs...)
+end
 
 index_gen{IT}(::Type{IT}, N::Integer, offsets...) = index_gen(iextract(IT, min(length(offsets)+1, N)), IT, N, offsets...)
 
@@ -88,6 +98,63 @@ end
     :(gradient!(g, itp, $(args...)))
 end
 
+@generated function gradient{T,N}(itp::AbstractInterpolation{T,N}, xs...)
+    n = count_interp_dims(itp, N)
+    Tg = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
+    xargs = [:(xs[$d]) for d in 1:length(xs)]
+    :(gradient!(Array($Tg,$n), itp, $(xargs...)))
+end
+
+gradient1{T}(itp::AbstractInterpolation{T,1}, x) = gradient(itp, x)[1]
+
+function hessian_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
+    meta = Expr(:meta, :inline)
+    # For each component of the hessian, alternately calculate
+    # coefficients and set component
+    n = count_interp_dims(IT, N)
+    exs = Expr[]
+    cntr = 0
+    for d1 in 1:N, d2 in 1:N
+        if count_interp_dims(iextract(IT,d1), 1) > 0 && count_interp_dims(iextract(IT,d2),1) > 0
+            cntr += 1
+            push!(exs, hessian_coefficients(IT, N, d1, d2))
+            push!(exs, :(@inbounds H[$cntr] = $(index_gen(IT, N))))
+        end
+    end
+    hessian_exprs = Expr(:block, exs...)
+
+    quote
+        $meta
+        size(H) == ($n,$n) || throw(ArgumentError(string("The size of the provided Hessian matrix wasn't a square matrix of size ", size(H))))
+        @nexprs $N d->(x_d = xs[d])
+
+        $(define_indices(IT, N, Pad))
+
+        $hessian_exprs
+
+        H
+    end
+end
+
+@generated function hessian!{T,N}(H::AbstractMatrix, itp::BSplineInterpolation{T,N}, xs::Number...)
+    length(xs) == N || throw(ArgumentError("Can only be called with $N indexes"))
+    hessian_impl(itp)
+end
+
+@generated function hessian!{T,N}(H::AbstractMatrix, itp::BSplineInterpolation{T,N}, index::CartesianIndex{N})
+    args = [:(index[$d]) for d in 1:N]
+    :(hessian!(H, itp, $(args...)))
+end
+
+@generated function hessian{T,N}(itp::AbstractInterpolation{T,N}, xs...)
+    n = count_interp_dims(itp,N)
+    TH = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
+    xargs = [:(xs[$d]) for d in 1:length(xs)]
+    :(hessian!(Array($TH,$n,$n), itp, $(xargs...)))
+end
+
+hessian1{T}(itp::AbstractInterpolation{T,1}, x) = hessian(itp, x)[1,1]
+
 offsetsym(off, d) = off == -1 ? symbol("ixm_", d) :
                     off ==  0 ? symbol("ix_", d) :
                     off ==  1 ? symbol("ixp_", d) :

src/b-splines/linear.jl

@@ -25,6 +25,13 @@ function gradient_coefficients(::Type{BSpline{Linear}}, d)
     end
 end
 
+function hessian_coefficients(::Type{BSpline{Linear}}, d)
+    sym, symp = symbol("c_",d), symbol("cp_",d)
+    quote
+        $sym = $symp = 0
+    end
+end
+
 # This assumes fractional values 0 <= fx_d <= 1, integral values ix_d and ixp_d (typically ixp_d = ix_d+1,
 #except at boundaries), and an array itp.coefs
 function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Linear}}, ::Type{IT}, N::Integer, offsets...)

src/b-splines/quadratic.jl

@@ -59,6 +59,26 @@ function gradient_coefficients{Q<:Quadratic}(::Type{BSpline{Q}}, d)
     end
 end
 
+"""
+In `hessian_coefficients` for a quadratic b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {m, _, p}` such that
+
+    cm_d  = p''(fx_d)
+    c_d   = q''(fx_d)
+    cp_d  = p''(1-fx_d)
+
+where `p` and `q` are defined in the docstring entry for `Quadratic`, and
+`fx_d` in the docstring entry for `define_indices_d`.
+"""
+function hessian_coefficients{Q<:Quadratic}(::Type{BSpline{Q}}, d)
+    symm, sym, symp = symbol("cm_",d), symbol("c_",d), symbol("cp_",d)
+    quote
+        $symm = 1
+        $sym  = -2
+        $symp = 1
+    end
+end
+
 # This assumes integral values ixm_d, ix_d, and ixp_d,
 # coefficients cm_d, c_d, and cp_d, and an array itp.coefs
 function index_gen{Q<:Quadratic,IT<:DimSpec{BSpline}}(::Type{BSpline{Q}}, ::Type{IT}, N::Integer, offsets...)

src/cubic.jl

@@ -1,4 +1,4 @@
-# module CubicTests
+module CubicTests
 
 using Base.Test
 using Interpolations
@@ -50,4 +50,36 @@ for (constructor, copier) in ((interpolate, identity), (interpolate!, copy))
     end
 end
 
-# end
+end
+
+module CubicGradientTests
+
+using Interpolations, Base.Test
+
+ix = 1:15
+f(x) = cos((x-1)*2pi/(length(ix)-1))
+g(x) = -2pi/14 * sin((x-1)*2pi/(length(ix)-1))
+
+A = f(ix)
+
+for (constructor, copier) in ((interpolate, identity), (interpolate!, copy))
+
+    for BC in (Line, Flat, Free, Periodic), GT in (OnGrid,OnCell)
+
+        itp = constructor(copier(A), BSpline(Cubic(BC())), GT())
+        # test that inner region is close to data
+        for x in linspace(ix[5],ix[end-4],100)
+            @test_approx_eq_eps g(x) gradient(itp,x)[1] cbrt(cbrt(eps(g(x))))
+        end
+    end
+end
+itp_flat_g = interpolate(A, BSpline(Cubic(Flat())), OnGrid())
+@test_approx_eq_eps gradient(itp_flat_g, 1)[1] 0 eps()
+@test_approx_eq_eps gradient(itp_flat_g, ix[end])[1] 0 eps()
+
+itp_flat_c = interpolate(A, BSpline(Cubic(Flat())), OnCell())
+@test_approx_eq_eps gradient(itp_flat_c, .5)[1] 0 eps()
+@test_approx_eq_eps gradient(itp_flat_c, ix[end]+.5)[1] 0 eps()
+
+
+end