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

[tomasaschan]

This is a WIP implementation for #94.

I still have to add tests, and the gradients aren't well-behaved for all boundary conditions. Not sure if the problem is with the gradients or (more likely) with the implementations of the boundary conditions in the cubic prefiltering functions.

using Interpolations, Gadfly
f0(x) = (x.^3 + 3*x.^2 - 15*x - 8) ./ 4
g0(x) = (3*x.^2 + 6*x - 15) ./ 4
f(x) = f0(x-9)
g(x) = g0(x-9)

ix, xs = 1:15, linspace(1,15,100)
itp = interpolate(f(ix), BSpline(Cubic(Line())), OnGrid())

plot(
    layer(x=ix, y=f(ix), Geom.point),
    layer(x=ix, y=g(ix), Geom.point),
    layer(x=xs, y=[itpc[x] for x in xs], Geom.line),
    layer(x=xs, y=[gradient(itpc, x)[1] for x in xs], Geom.line),
)

itpc = interpolate(f(ix), BSpline(Cubic(Flat())), OnGrid());
plot(
    layer(x=ix, y=f(ix), Geom.point),
    layer(x=ix, y=g(ix), Geom.point),
    layer(x=xs, y=[itpc[x] for x in xs], Geom.line),
    layer(x=xs, y=[gradient(itpc, x)[1] for x in xs], Geom.line),
)

[tomasaschan]

cc @spencerlyon2

[sglyon]

Hmm something is definitely off here. I'll take another look at my derivations and make sure the implementation matches .

I liked my Mathematica notebook as I didn't actually do any derivations by hand, I just let Mathematica take derivatives and collect coefficients for me. I think the most likely place for an error in the math is at the very beginning when I derived the polynomial form of the cubic spline

[tomasaschan]

Yeah - that is definitely a great selling point for doing things like this in Mathematica. (An embarassingly large portion of the time it took me to implement this was spent taking derivatives of polynomials...)

I would definitely be interested if there was some place online to share (and possibly actually execute) Mathematica notebooks, perhaps similar to JuliaBox or even nbviewer.

[tomasaschan]

Re: where the problem lies; it seems in the plots like it's mainly the boundary conditions that are off. The interior domain seems fine.

[sglyon]

Hmm that's true.

I'm surprised by this because once I have the main polynomial I don't do any calculations by hand.

Anyway, I'll double check and report back

[tomasaschan]

👍

[sglyon]

Hmm, I'm not sure what the problem is here...

Let's consider Line, OnGrid. Here we should set y''(0) = 0

The polynomial for the interior is: y(x) = { lower order stuff} + 1/2 (cp - 2c + cm)x^2 + 1/6(cpp - 3cp + 3p - cm)x^3.

This means that

y''(x) = cp - 2c + cm+ (cpp - 3cp + 3p - cm)x ==>
y'(0) = cp - 2c + cm

Which is exactly the condition implemented here (actual implementation is found in quadratic.jl)

The one thing I can see that is wrong is that we shouldn't have a type parameter over GT in cubic.jl there as we only borrow the quadratic implementation for OnGrid, but that shouldn't be causing the problem we are seeing above...

Any thoughts?

[sglyon]

This has been merged into Woodbury is now (unexported) WoodburyMatrices. sparse_factors in the latest tagged release (0.1.4)

[tomasaschan]

See #104

[tomasaschan]

Ah, that explains it. The lower boundary for OnGrid is at x=1, not x=0. (For OnCell it's at x=0.5 as you'd expect.)

[sglyon]

Oh man... my mistake. Sorry about that. how do we know that it is at 1?

[tomasaschan]

Well, the underlying data is usually an array, with indices 1:N. Starting the interpolation at 1 is necessary to make itp[i] = A[i] for integers.

I don't know if it's explicitly documented anywhere, but it should be implied in the latex docs (and by lbound(itp) == 1 :)).

[tomasaschan]

Ha! I found it!

It's not the boundary conditions after all - although you did the derivations of boundary conditions with the wrong boundary in mind, the results are actually equivalent; you considered y_0(x) at x=0, while it would have been more correct to consider y_1(x) at x=1. However, any choice of i would give the same results when considering y_i(x) at x=i, so the resulting matrix systems are correct anyway (although we might want to update the docstrings).

However, you copied one thing too many from the implementation of quadratic b-splines :)

A quadratic b-spline is supported by three coefficients, and centering over them the boundaries end up at half-intervals. In other words, y_i(x) is valid for x ∈ [i-.5, i+.5]. Cubic b-splines are supported by four coefficients, and so centering over them the boundaries end up at the data points; y_i(x) is valid for x ∈ [i, i+1]. These two calls to round should have been floor. (What tipped me off was the discontinuities in the function value, which were located at half-points rather than integers...)

With round:

With floor:

I'll patch master momentarily.

[sglyon]

Good catch. So with #102 is this fixed?

[tomasaschan]

So with #102 is this fixed?

Yep, should be.

I'm going to spend some time on this later tonight (in UTC+1...), hopefully outputting some tests for this PR and other things, as well as tagging and releasing a new version.

[sglyon]

Sounds great, thanks. Let me know if you need any help!

[tomasaschan]

Sigh. After spending way more time than I should on troubleshooting something that in the end turned out to be equivalent to this:

for GT in (OnGrid, OnCell)
    itp = interpolate(A, BSpline(Cubic(Line())), OnCell()) # note: OnCell(), not GT()
    #plot stuff
end

I now have a working implementation of both gradient and hessian for all degrees up to cubic, as well as a good start on migrating the documentation to docstrings for at least quadratic. It'll take a little time to clean up and put it all in neat commits, but I hope to get it in tomorrow.

[timholy]

Wow, great!

Sympathies for the misguided bug hunt. I had a remarkably similar thing happen to me last week, and it cost me a lot of time, too. Why is it always the bugs in the test suite that are the worst to track down? (I guess there's a good answer to that, but still...)

[tomasaschan]

I haven't written a test suite for these implementations yet, but given the plots below I feel fairly confident that they are correct (there are some tests for the cubic gradient). Unfortunately, I'm going to have less time to spend on Julia over the coming weeks/months, so I want to make sure the work I've done so far this time isn't left hanging.

My proposed plan of action is this:

1. Merge this, to make the functionality available
2. Merge a PR with documentation that I'll file shortly after finishing this comment

(and, later, or by someone else than me)

3. Re-vamp the `InPlace`/`InPlaceQ` stuff into an `Interior` boundary condition, which simply uses the outermost data points as boundary conditions, and is interpolating *exactly* on `2:N-1` for N data points (in 1D).
4. Re-vamp the test suite to test the interior domains using `Quadratic{Interior}` and `Cubic{Interior}`, and verifying the boundary conditions by their definitions, using `gradient` and `hessian`. For example, with an underlying cubic polynomial, we can test both that the interpolant is cubic, that the gradient is the corresponding quadratic and that the hessian is the corresponding linear function. This, together with some sanity tests for generalizations into multiple dimensions, should cover our bases pretty well.

Before I hit merge on this, I would like at least one other person to say that this is a sane plan, given that merging this in its current state would mean adding quite a lot of code that has basically no automated tests at all.

Plots for sanity tests (generating script at the bottom):

using Interpolations, Gadfly

ix = 1:15
k = 2pi/(length(ix)-1)
f(x) = sin((x-1) * k) 
g(x) = k * cos((x-1) * k)
h(x) = -k^2 * sin((x-1) * k)

b(::OnGrid) = 0
b(::OnCell) = .5
xs = linspace(.5,15.5,100)

fth = Theme(default_color=colorant"blue")
gth = Theme(default_color=colorant"green")
hth = Theme(default_color=colorant"orange")

for BC in (
        Constant(),
        Linear(),
        Quadratic(Line()),
        Cubic(Line()),
    ),
    GT in (
        OnGrid(),
        OnCell()
    )

    itp = interpolate(f(ix), BSpline(BC), GT)

    config = replace((@sprintf "%s, %s" typeof(BC) typeof(GT)), "Interpolations.", "")

    p = plot(
    layer(x=xs,y=[itp[x] for x in xs], Geom.line, fth),
    layer(x=ix,y=f(ix),Geom.point,fth),
    layer(x=xs,y=[gradient1(itp,x)[1] for x in xs], Geom.line, gth),
    layer(x=ix,y=g(ix),Geom.point,gth),
    layer(x=xs,y=[hessian1(itp,x) for x in xs], Geom.line, hth),
    layer(x=ix,y=h(ix),Geom.point,hth),
    Guide.title(config)
    )
    display(p)
    draw(PNG("$config.png", 6inch, 4inch), p)
end

[tomasaschan]

For reference: when this merges, it should fix #94 and #95.

[sglyon]

Hey @tlycken sorry for the delay.

Can you remind me which parts of this PR will introduce untested code?

Otherwise I think your proposal sounds great

[tomasaschan]

Can you remind me which parts of this PR will introduce untested code?

I've only written a few sanity tests for gradient on cubic b-splines - I'm not sure I covered all possible BC's, for example. And there are no tests at all for hessian.

I've also exported gradient1 and hessian1 as convenient methods for 1D interpolations; they just forward to gradient and hessian and then index into the result to return a scalar. These are also completely untested. (They also leave some room for improvement when it comes to performance - with specialized implementations it should be possible to avoid a whole bunch of array allocations and save a lot of memory...)

[tomasaschan]

I'm taking @timholy's silence as agreement with the plan, and hitting the button on this :)

[timholy]

I'm super-excited about this, but just haven't had two halves-of-an-hour to rub together. Good that you hit merge. I will check this out soon, I promise.