2nd order derivative

[tomasaschan]

2nd order derivatives implemented through the Hessian matrix

Addition to high-level interface:

# example from test code
julia> y = rand(7,8);
julia> yi = InterpGrid(y, BCnil, InterpQuadratic);
julia> x = [2.2, 3.1];
julia> v,g,h = valgradhess(yi, x);  # new!
julia> h
  2x2 Array{Float64,2}:
   -0.0846558   0.148568 
    0.148568   -0.0781553

Addition to low-level interface:

julia> yc = copy(y)
julia> interp_invert!(yc, BCnan, InterpQuadratic)
julia> ic = InterpGridCoefs(yc, InterpQuadratic)
julia> set_hessian_coordinate(ic, 1, 2) # new! cross-derivative, d2/dxdy
julia> set_hessian_coordinate(ic, 2, 2) # new! second-order derivative d2/dy2

Note:

In order to get valgradhess working correctly, I had to introduce another parameter to set_position, as part of the low-level interface. Currently, the change I made is breaking (one needs to provide a second boolean argument, indicating wether to calculate hessian coefficients) but one might want to make that change non-breaking somehow. I don't know if that function is used a lot outside of the library (it is exported), but I couldn't make up a non-breaking way to change that without keeping its api consistent with itself.

A better (more future-proof), but still breaking, change, could be to make calc_grad and calc_hess keyword arguments, allowing to add further ones if desired.

Todo:

• [x] Documentation.

[tomasaschan]

The travis failure is for Julia 0.2 - it passes on 0.3 - and for something that I believe is unrelated. Should maybe REQUIRE.jl be updated to require julia 0.3?

[timholy]

The odd part is that the REQUIRE is updated for 0.3, and I've even made a release. Sorry to ping you yet again, @StefanKarpinski, but I'm confused about why this happened. Are we asking Travis to test this particular commit, no matter what?

[timholy]

Looks really good. The one change I'd suggest is that I don't think we actually need to make this breaking; we can introduce a version of set_position that has the original signature and just calls the new set_position with calc_hess=false. I agree it might be more flexible to make them keyword arguments, but I do worry about the performance implications. Grid's performance is something I should look into anyway---I suspect there are a number of places it can become better.

[tomasaschan]

That seems like a reasonable approach. I'll add a few lines.

[timholy]

Wanted to advertise that I just pushed another change (508d2d9774d3e7255490837b26867bd6bdf49006), which you may want if you're implementing CubicSpline. I noticed an opportunity for improvement when reviewing your changes, and it's had a major beneficial impact on performance. Sorry if this causes you any merging difficulties.

[tomasaschan]

No merging problems at all :)

I do have another question though; it doesn't seem to cause a problem - all interpolation tests pass - but there's something going on here that I can't seem to wrap my head around fully.

For the gradient, we visit each coordinate direction and save the 1D interpolation pattern there. For the hessian, it seems we only visit the diagonal elements. Should ic.hcoef1d really have been a Matrix{Vector{T}}, with hessian 1D coefs along the diagonal and gradient coefs off the diagonal? Could that maybe have been utilized so that some of the branching here could have been avoided?

[timholy]

Yes, that's right. There may not be substantial performance implications for doing it the way you're doing it currently, unless you're thinking about interpolation of multi-valued functions (an RGB image, for example). In such cases, you probably would be better off with a Matrix{Vector{T}}.