Hessian fails for function of two variables

[simonp0420]

In converting over some code that currently uses Grid.jl, I need to compute the curvature of a surface. The following REPL transcript illustrates the problem I'm having:

julia> A = rand(5,5)
5×5 Array{Float64,2}:
 0.264768   0.405079  0.253461  0.759207  0.673105
 0.0884208  0.528674  0.244601  0.545683  0.751446
 0.403174   0.68685   0.722238  0.746143  0.42529
 0.396175   0.651159  0.909525  0.402493  0.251807
 0.130116   0.986705  0.259471  0.51158   0.850608

julia> itp = interpolate(A, BSpline(Cubic(Natural())), OnGrid())
5×5 Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}:
 0.264768   0.405079  0.253461  0.759207  0.673105
 0.0884208  0.528674  0.244601  0.545683  0.751446
 0.403174   0.68685   0.722238  0.746143  0.42529
 0.396175   0.651159  0.909525  0.402493  0.251807
 0.130116   0.986705  0.259471  0.51158   0.850608

julia> h = hessian(itp,2.0,2.0)
ERROR: UndefVarError: dim not defined
 in hessian_coefficients(::Type{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}}}, ::Int64, ::Int64, ::Int64) at c:\home\simonp\.julia\v0.5\Interpolations\src\b-splines\indexing.jl:18
 in hessian_impl(::Type{Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}}) at c:\home\simonp\.julia\v0.5\Interpolations\src\b-splines\indexing.jl:123
 in hessian!(...) at c:\home\simonp\.julia\v0.5\Interpolations\src\b-splines\indexing.jl:144
 in hessian(::Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}, ::Float64, ::Float64) at c:\home\simonp\.julia\v0.5\Interpolations\src\b-splines\indexing.jl:153

This was run using Julia 0.5.0 and Interpolations 0.3.6. Note: I get the same error using OnCell() and when checking out the master version of Interpolations. Is hessian supposed to work for surfaces?

Thanks,
--Peter

[sglyon]

Error is in these lines

https://github.com/tlycken/Interpolations.jl/blob/ac5475dc6fe76cc28c0b88680460d63f640919ca/src/b-splines/indexing.jl#L18-L22

[aaronc8]

Recently using this package, I ran into the same error as well. Do we have to manually edit in the fix for these lines, and what is it? Or, would it be fixed in an update to the package?

[keatonmill]

I run into a similar, related problem:

julia> A = rand(5,5)
5×5 Array{Float64,2}:
 0.368021   0.440408  0.659274    0.629533  0.215752
 0.0122785  0.74314   0.192229    0.77579   0.937559
 0.0333598  0.757276  0.00512945  0.923899  0.312257
 0.479558   0.59737   0.238214    0.812611  0.882245
 0.25881    0.341119  0.384092    0.239488  0.797613

julia> itp = interpolate(A, BSpline(Cubic(Natural())), OnGrid())
5×5 interpolate(::Array{Float64,2}, BSpline(Cubic(Line())), OnGrid()) with element type Float64:
 0.368021   0.440408  0.659274    0.629533  0.215752
 0.0122785  0.74314   0.192229    0.77579   0.937559
 0.0333598  0.757276  0.00512945  0.923899  0.312257
 0.479558   0.59737   0.238214    0.812611  0.882245
 0.25881    0.341119  0.384092    0.239488  0.797613

julia> itp[2.5, 2.5]
0.38246048600552157

julia> gradient(itp, 2.5, 2.5)
2-element Array{Float64,1}:
 -0.158142
 -1.03777 

julia> hessian(itp, 2.5, 2.5)
ERROR: UndefVarError: TH not defined
 in hessian(::Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}, ::Float64, ::Float64) at /home/keaton/.julia/v0.5/Interpolations/src/b-splines/indexing.jl:160

[timholy]

The TH fix is easy, just change it to $TH. But the lines pointed out by @sglyon are much more problematic, those lines look quite borked to me. Someone will probably have to work through the math to figure out what really needs to be done.

[timholy]

Sorry. Use the $ inside the quoted expression, :(hessian(Array{$TH, 2}....).

[keatonmill]

Aha! Clearly I need to learn something about metaprogramming. Still get an error though. Maybe this is what you were referring to with @sglyon 's comment?

hessian(itp, 2.5 2.5)

UndefVarError: dim not defined
 in hessian_coefficients(::Type{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}}}, ::Int64, ::Int64, ::Int64) at indexing.jl:19
 in hessian_impl(::Type{Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}}) at indexing.jl:129
 in hessian!(...) at indexing.jl:151
 in hessian(::Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1}, ::Float64, ::Float64) at indexing.jl:160
 in include_string(::String, ::String) at loading.jl:441
 in eval(::Module, ::Any) at boot.jl:234
 in (::Atom.##65#68)() at eval.jl:102
 in withpath(::Atom.##65#68, ::Void) at utils.jl:30
 in withpath(::Function, ::Void) at eval.jl:38
 in macro expansion at eval.jl:101 [inlined]
 in (::Atom.##64#67{Dict{String,Any}})() at task.jl:60

I tried it on a one-dimensional interpolation and got this:

hessian(sitp, 3.1415)

MethodError: no method matching hessian!(::Array{Float64,2}, ::Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,Tuple{FloatRange{Float64}}}, ::Float64)
Closest candidates are:
  hessian!{T,N}(::AbstractArray{T,2}, !Matched::Interpolations.BSplineInterpolation{T,N,TCoefs<:AbstractArray,IT<:Union{Interpolations.BSpline,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.BSpline,Interpolations.NoInterp},N}}},GT<:Union{Interpolations.GridType,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.GridType,Interpolations.NoInterp},N}}},pad}, ::Number...) at /home/keaton/.julia/v0.5/Interpolations/src/b-splines/indexing.jl:150
  hessian!{T,N}(::AbstractArray{T,2}, !Matched::Interpolations.BSplineInterpolation{T,N,TCoefs<:AbstractArray,IT<:Union{Interpolations.BSpline,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.BSpline,Interpolations.NoInterp},N}}},GT<:Union{Interpolations.GridType,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.GridType,Interpolations.NoInterp},N}}},pad}, !Matched::CartesianIndex{N}) at /home/keaton/.julia/v0.5/Interpolations/src/b-splines/indexing.jl:155
 in hessian(::Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,1},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.OnGrid,Tuple{FloatRange{Float64}}}, ::Float64) at indexing.jl:160
 in include_string(::String, ::String) at loading.jl:441
 in eval(::Module, ::Any) at boot.jl:234
 in (::Atom.##65#68)() at eval.jl:102
 in withpath(::Atom.##65#68, ::Void) at utils.jl:30
 in withpath(::Function, ::Void) at eval.jl:38
 in macro expansion at eval.jl:101 [inlined]
 in (::Atom.##64#67{Dict{String,Any}})() at task.jl:60

[timholy]

Yep, that's what I was referring to.

Looks like the hessian stuff isn't as well tested as it needs to be. I've not used hessians myself because the optimization routines I use for my problems only exploit the gradient.