Merge pull request #99 from tlycken/tl/cubic-gradient

gradient for cubic b-splines, and hessian for all b-splines

Author: Tomas Lycken

src/Interpolations.jl

@@ -7,6 +7,10 @@ export
     scale,
 
     gradient!,
+    gradient1,
+    hessian!,
+    hessian,
+    hessian1,
 
     AbstractInterpolation,
     AbstractExtrapolation,
@@ -74,17 +78,11 @@ 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")
 include("extrapolation/extrapolation.jl")
 include("scaling/scaling.jl")
+include("utils.jl")
 
 end # module

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

@@ -81,6 +81,56 @@ function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
     end
 end
 
+"""
+In this function we assume that `fx_d = x-ix_d` and we produce `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 for `Cubic`.
+"""
+function gradient_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)
+    symfx_sqr = symbol("fx_sqr_", d)
+    sym_1m_fx_sqr = symbol("one_m_fx_sqr_", d)
+    quote
+        $symfx_sqr = sqr($symfx)
+        $sym_1m_fx_sqr = sqr(1 - $symfx)
+
+        $symm  = -SimpleRatio(1,2) * $sym_1m_fx_sqr
+        $sym   =  SimpleRatio(3,2) * $symfx_sqr     - 2 * $symfx
+        $symp  = -SimpleRatio(3,2) * $sym_1m_fx_sqr + 2 * (1 - $symfx)
+        $sympp =  SimpleRatio(1,2) * $symfx_sqr
+    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
@@ -155,7 +205,7 @@ condition gives:
 """
 function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                    ::Type{Cubic{Line}}, ::Type{OnCell})
-    dl,d,du = inner_system_diags(T,n,Cubic{Flat})
+    dl,d,du = inner_system_diags(T,n,Cubic{Line})
     d[1] = d[end] = 3
     du[1] = dl[end] = -7
 
@@ -181,18 +231,29 @@ condition gives:
 This is the same system as `Quadratic{Line}`, `OnGrid` so we reuse the
 implementation
 """
-function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
-                                                ::Type{Cubic{Line}}, ::Type{GT})
-    prefiltering_system(T, TC, n, Quadratic{Line}, OnGrid)
+function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
+                                   ::Type{Cubic{Line}}, ::Type{OnGrid})
+    dl,d,du = inner_system_diags(T,n,Cubic{Line})
+    d[1] = d[end] = 1
+    du[1] = dl[end] = -2
+
+    # now need Woodbury correction to set :
+    #    - [1, 3] and [n, n-2] ==> 1
+    specs = _build_woodbury_specs(T, n,
+                                  (1, 3, one(T)),
+                                  (n, n-2, one(T)),
+                                  )
+
+    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
 end
 
 function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                 ::Type{Cubic{Periodic}}, ::Type{GT})
     dl, d, du = inner_system_diags(T,n,Cubic{Periodic})
 
     specs = _build_woodbury_specs(T, n,
-                                  (1, n, SimpleRatio(1, 6)),
-                                  (n, 1, SimpleRatio(1, 6))
+                                  (1, n, du[1]),
+                                  (n, 1, dl[end])
                                   )
 
     Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
@@ -210,38 +271,12 @@ This is the same system as `Quadratic{Free}` so we reuse the implementation
 """
 function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                 ::Type{Cubic{Free}}, ::Type{GT})
-    prefiltering_system(T, TC, n, Quadratic{Free}, GT)
-end
-
-@inline cub(x) = x*x*x
-
-"""
-Build `rowspec`, `valspec`, `colspec` such that the product
-
-`out = rowspec * valspec * colspec` will be equivalent to:
-
-```julia
-out = zeros(n, n)
-
-for (i, j, v) in args
-    out[i, j] = v
-end
-```
-
-"""
-function _build_woodbury_specs{T}(::Type{T}, n::Int, args::Tuple{Int, Int, Any}...)
-    m = length(args)
-    rowspec = spzeros(T, n, m)
-    colspec = spzeros(T, m, n)
-    valspec = zeros(T, m, m)
+    dl, d, du = inner_system_diags(T,n,Cubic{Periodic})
 
-    ix = 1
-    for (i, (row, col, val)) in enumerate(args)
-        rowspec[row, ix] = 1
-        colspec[ix, col] = 1
-        valspec[ix, ix] = val
-        ix += 1
-    end
+    specs = _build_woodbury_specs(T, n,
+                                  (1, n, du[1]),
+                                  (n, 1, dl[end])
+                                  )
 
-    rowspec, valspec, colspec
+    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
 end

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...)