Gradient evaluation

[tomasaschan]

I got it working! =)

The interface is currently a little clunky, especially in 1D, but I couldn't think of a way to define an overload of gradient for vector interpolations, that wouldn't be recursive. This is how it works:

using Interpolations
xg = 1:15
xf = 1:.1:15
f1(x) = sin(2pi/15 * (x-1))
g1(x)= 2pi/15 * cos(2pi/15*(x-1))

f1g = f1(xg)
g1g = g1(xg)

itp1 = Interpolation(f1g, Quadratic(Free(),OnCell()), ExtrapNaN())
f1i = Float64[itp1[x] for x in xf]
g1i = Float64[gradient(itp1, x)[1] for x in xf]

using Gadfly
plot(
    layer(x=xg,y=fg,Geom.point),
    layer(x=xg,y=gg,Geom.point,Theme(default_color=color("red"))),
    layer(x=xf,y=f1i,Geom.path),
    layer(x=xf,y=g1i,Geom.path,Theme(default_color=color("red"))),
)

It works for linear and constant interpolations too, although the results are of course less useful:

The meta-programming mumbo-jumbo used to define the correct coefficients for evaluation became more complicated, and as a result (other than the code being a little more opaque), the loading time of the package is bumped up even further; tic(); using Interpolations; toc() now reports 31 seconds on my machine. For comparison, tic(); using Gadfly; toc() reports 25... It might be possible to reduce the load time slightly by organizing the code into smaller modules where things like symm, sym and symp could be defined as constants, but I guess much of this would be alleviated with staged functions, so I'm not very concerned about it at the moment.

[tomasaschan]

There, some tests and multidimensional working too =)

Since the interpolated gradient is pretty far from the exact one, I had a hard time writing good tests for more than one dimension. I've verified that it works visually, though, using f(x,y) = sin(2pi/nx * (x-1)) * cos(4pi/ny * (y-1)):

Analytical gradient in the x-direction:

Interpolated gradient in x-direction:

I'm pretty happy with that result =)

[tomasaschan]

@timholy, do you have a magic way to avoid having to do gradient(itp, x)[1] for 1-dimensional interpolations?

[timholy]

The only method I can think of is to write a special method for 1d. But would the fact that it doesn't always return an array be a hindrance to generic programming?

Also, you want to write the basic method as gradient!(g::Array{T,1}, ...) for a pre-allocated g, then add a little stub function

gradient(itp::Interpolation, ...) = gradient!(Array(gradient_eltype(itp), ndims(itp)), itp, ...)

That way users who need to evaluate many gradients can avoid allocation by calling gradient! directly.

[tomasaschan]

Good points - we probably do want gradient(itp{T,1}, x) to return an 1-element array.

I've now made the bulk method gradient!, which calculates into a pre-allocated array, and created a stub for gradient just as you suggested. It doesn't handle types very gracefully (it just takes T from the interpolation object...) but I'm thinking it's easier to fix that all at once when I take care of #17.