Implement hessian

[tomasaschan]

From Wikipedia:

In mathematics, the Hessian matrix or Hessian is a square matrix of second-order partial derivatives of a scalar-valued function, or scalar field. It describes the local curvature of a function of many variables.

It is basically a (symmetric) matrix of second-order partial derivatives. For scalar-valued function of a scalar, f(x), it is simply the second derivative; for a vector-valued function of a vector, it is

          |  f''_xx(v)     f''_xy(v) |
H(f)(v) = |                          |
          |  f''_yx(v)     f''_yy(v) |

or on component form, H_{i,j} = \frac{\partial^2 f}{\partial x_i \partial x_j} (yes, until #46 is in place you have to render LaTeX in your head. Or read the wiki entry...).

Implementations for lower degrees than cubic are trivial - they should just return 0 everywhere.

[timholy]

Actually, for Quadratic it should be piecewise constant. Similar to this:

julia> using Interpolations

julia> a = rand(5)
5-element Array{Float64,1}:
 0.139006
 0.984092
 0.827687
 0.822077
 0.179375

julia> itp = interpolate(a, BSpline(Linear()), OnGrid())
5-element Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0}:
 0.139006
 0.984092
 0.827687
 0.822077
 0.179375

julia> gradient(itp, 2.2)
1-element Array{Float64,1}:
 -0.156405

julia> gradient(itp, 2.8)
1-element Array{Float64,1}:
 -0.156405

julia> gradient(itp, 3.1)
1-element Array{Float64,1}:
 -0.0056092

[tomasaschan]

You're right, of course.

[tomasaschan]

Fixed in #99.