More boundary conditions for quadratic B-splines

[tomasaschan]

I made an attempt at building linear extrapolation and boundar conditions as well, but I can't seem to get it right.

The math behind is quite straightforward; in the middle of the (1D) domain, the relevant in the tridiagonal system matrix is [1/8 3/4 1/8]. In order to specify the f' at the edge, I introduce a ghost point y0 at x=0 and give it a value such that it is on the line between y1 and y2, yielding y0 = 2y1 - y2. Inserting in the equation y0/8 + 3y1/4 + y2/8 = A[1] yields y1 = 0. Thus, I modify the system (here) so that the first row is [1, 0, ...] and the last row is [ ... 0, 1 ].

However, this doesn't seem to give me enough to get where I want - I suspect I have to do something more (like what I did for flat BC here) but I can't figure out what. I just stole that piece from Grid, but there is no corresponding BC implemented in Grid to look at for this.

[timholy]

Unless I'm confused, I get that the "edge" equation is y1-y2/8 = A[1], not y1=0. So I'd set d[1] = 1 (as you have done) but du[1] = dl[1] = -1/8. Of course, you also have to fix up the tail.

[tomasaschan]

Hm, I meant to say y1=A[1], but I'm still missing the second term so that may very well be it. I'll try it out when I've recharged my laptop :p

[tomasaschan]

Nope, that didn't work either. I'm thinking maybe we're overlooking something here...

I wrote a small latex document outlining the math behind what I've tried - I'd be very grateful if you had time to look at it and see if I'm overlooking something fundamental.

[timholy]

I'll definitely try to take a look. Your formula is the right one; I had the sign of c0 wrong.

It occurs to me that the source of the problem might be that you want to constrain the slope of the interpolated function, and that c0 = c1 - (c2 - c1) is probably not quite right.

[tomasaschan]

To make this even more complicated, I just reminded myself about our discussion about on-grid/on-cell interpolations; the governing equations for these things will of course need to be different depending on which type we're talking about.

I've started patching up the work I've already committed to take this into account, using the following assumptions and consequences:

1. Regardless of OnGrid or OnCell, the interpolation is defined so that if itp = Interpolation(A, ...), then itp[i] == A[i] for all integers i between 1 and length(A), inclusive. (The corresponding statements hold true for all higher dimensions as well, of course.) In other words, regardless of if grid cell boundaries are defined on the data point or between them, only the grid cells are translated, but never the data.
2. For OnGrid, the domain is naturally limited by the end points of the data, such that itp[x] is interpolation (as opposed to extrapolation) if and only if 1 <= x <= lenght(A) (and corresponding for higher dimensions).
3. Here comes the tricky part: for OnCell, the outermost data points are specified as midpoints in the grid cells. Therefore, itp[x] is actually interpolation in the interval 0.5 <= x <= length(A) + 0.5 (etc), i.e. also outside the given data, since the first cell is necessarily defined to be from x=0.5 to x=1.5.
4. To make life easier for ourselves, we could also make the following rather bold assumption about our users' intentions when specifying interpolation degree and grid representation: if they care about discontinuities in some order of the derivative, they choose an interpolation degree that is continuous under so many derivations - grid representation is only chosen to specify if the domain edges are on or outside the data points.

Following from definition 3, we should specify the boundary conditions for OnCell at x = 0.5, not at x=1. For example, this yields for the flat BC that c0 == c1, not c0 == c2 as it is for OnGrid.

From 4, we can also conclude that we don't need to worry too much about where the implementation makes the switch between cells inside the domain, as long as we get the boundary conditions and extrapolation behavior correct under both grid representations, it's perfectly fine e.g. for Quadratic to actually switch between interpolating polynomials between data points even when OnGrid is specified.

So the code I've shown here is really only valid for OnGrid, and not for OnCell. However, when we get OnGrid working it's probably quite trivial to figure out what changes are needed for OnCell to work as well.

[timholy]

OK, I looked at this a bit. I think we have to introduce two ghost points, at 0 and -1, because the interpolating function is quadratic and so we have two equations to solve: one that the coefficient of the quadratic component is zero (gives linear behavior), and the other that the slope matches the slope we want.

If you have a PDF document I believe I sent quite some months ago, this refers to Eqs. 1 and 7-10 (if I didn't edit the document in the meantime 😄). The equation I get for the vanishing of the quadratic term (on the interval from -0.5 to 0.5) is c_{-1} - 2c_0 +c_1 = 0, and the slope is just (c_1-c_{-1})/2, which you'd match to y_2-y_1, I guess.

[tomasaschan]

There is another way too: just introduce one equation, that specifies that the quadratic component vanishes, and then use whatever slope that turns out to be for the linear extrapolation (i.e. the linear term at the edge determines the slope of the extrapolation). I'm not sure which I think is more correct (or useful).

[timholy]

I like the "free floating" slope version, because you can use a grid of the same size as the original array. If you use my version, it seems that technically c_0 is not a ghost point: you'd have to pad by 1 along each edge (similar to the BCfill case in Grid). Unless memory is tight, padding probably isn't a big deal performance-wise, but there is at least a small advantage in avoiding it when possible. (To take real advantage of this, one would have to be using the Interpolation! variant.)

[tomasaschan]

I think I like the free slope better too, given that I can figure out how to evaluate it ;)

But setting the quadratic term to 0 is a good first step; for OnGrid interpolations, we want the quadratic term to vanish at x=1, yielding c_0 - 2c_1 + c_2 = 0, which is quite straightforward. I guess I can figure out how to get the linear term from here.
For OnCell, on the other hand, the BC should be applied at x = 0.5, yielding either c_0 - 2 c_{0.5} + c_1 or c_{-0.5} - 2 c_{0.5} + c_{1.5}, both of which are problematic because they require (at least) one data point which is not on the unit grid. Am I missing something obvious here?

[timholy]

I haven't read up on your math about OnCell, can you remind me where your notebook is stashed? I checked old PRs but the links I found were broken.

Bottom line, I suspect that you should insist that all "control points" are on the unit grid. Define the extrapolation function y(x) over each interval in terms of the coefficients c, just like we have y(x) = (x+0.5)^2/2*c_{-1} + (0.75-x^2)*c_0 + (x-0.5)^2/2*c_1 for the OnGrid case. Then derive the equation that gives y the desired properties (don't assume it looks like a discrete 2nd deriviative on c).

[tomasaschan]

Yeah, I thought some more about it and realized that the discrete 2nd derivative on c was just a coincidence - what I really want to do is to specify sharply that y'' = 0 at the boundary - wherever that is - and that will give the right expressions. Unfortunately, I took an invalid shortcut... I'll do the math more thoroughly and get back on this - thanks for kicking me in the right direction.

And for future reference, I put a link to the most relevant (although already slightly outdated) notebook in the readme =)

[timholy]

Thanks for the link.

[tomasaschan]

I feel this is merge-ready now, although there are still a couple of issues that I thought this would resolve when I started working on this PR. Most notably, the linear extrapolation doesn't work at all, because it gets its slope from some mathematically insane assumption. However, since I want to implement gradient and Hessian evaluation anyway, I figured it could wait until I can just do gradient(itp, 1) to get the slope at the lower boundary...

Also, there are a couple of other BC:s that aren't implemented yet, but the infrastructure changes here do make that stuff somewhat easier to deal with, so I feel it's better to merge this and start working on those separately.

There is quite possible some performance loss here, because the most significant change I've done is to introduce padding for all BC:s for quadratic interpolations. The rationale behind that is that I wanted to get away from the voodoo surrounding interpolation evaluation close to the boundary that we had before with the indices function; for some BC:s, for example, it pretended that c_i-1 = c_i+1 just to be able to evaluate the spline polynomial properly, but that didn't generalize to all boundary conditions we want (Tangent, for example, will never be able to work that way). Now, if the input matrix A has size(A,dim) = n, we solve a system of n+2 equations, where the right-hand side is [0, the_slice_of(A), 0] for most BC:s. I've written the code to try to minimize the number of allocate and/or copy operations, but I'm not sure I've done the best possible thing. The relevant code is, mainly, here and here.

[tomasaschan]

I got the periodic BC/extrapolation working too =)

Also, I have had a notebook with some Gadfly plots of the various configuraitons as a sanity check I figured I could just as well share it: here. Once it's merged, I'll link to it from the readme as well. (And all the code it runs is under test/visual.jl, so you can test it locally too!)

As you can see in the notebook, everything seems to work well except linear extrapolation, which - as noted above - I'll fix when gradient has been implemented.

[timholy]

I am not worried about the copy/padding operation; compared to solving the system of equations along slices, it will be trivial (for anything higher than 1d). However, right now the implementation of the filtering is sure to be very slow. I'd say, get started with a slow implementation and I'll try to write a standalone package to handle this step (I haven't had time to start it, however).

[timholy]

One really good way to debug 2d problems is to make one of the dimensions constant, i.e., just duplicate one column (or one row) many times. The advantage of this trick is then you can use the 1d functions for comparison, since the 1d functions applied to that same column (or row) should give the same answer.

[tomasaschan]

I've figured out what the problem is, and it turns out it's really simple: because I've padded the matrix that stores the coefficient, I have to compensate for that when indexing into it. I'll see if I can get a fix in later tonight (and if I can, I think I can pass all Grid's tests except the ones that make use of getindex(itp::Interpolation, x::AbstractVector) and friends).

[tomasaschan]

Hm. I got the indexing right now, but I still don't get correct results - it seems I still have some of the math related to the prefiltering step to grok. I'm all good with how it works in 1D, both mathematically and in code.

What I'm having problems with is how it works in the multi-dimensional case. Take 2D as an example: the input data to the interpolation is a matrix. When calculating the coefficients, we do it first column-wise, using each column of the matrix in the right-hand side of the same equation system as we did for 1D. Then, we do the same thing row-wise. (This is your "pencils aligned along all directions" analogy, if you remember...).

However, it seems that Grid.jl actually overwrites the column-wise solutions with the row-wise solutions. Have I understood this correctly? Why does this work? Is there something in particular that one has to do to make it work?

I ask mainly because now that I got the indexing working, I was fairly sure that I did almost exactly the same things as Grid.jl, except that I handle boundary conditions a little differently (but the result when evaluating should be the same). Also, my N-dimensional code works in 1D. However, I get completely nonsensical results in 2D, and since I'm not sure I understand what I'm aming for, it's difficult to know what might be wrong.

[timholy]

It's more accurate to say that Grid double-filters: it filters the column-wise solution over rows. I do everything in-place for performance reasons, but there's nothing magical about that---if you're taking the solution from column-wise filtering and feeding it into the row-wise filtering, it should work.

I haven't looked in detail at your filtering code, but I would recommend basing it on mapslices, unless you know your code is much faster. The only problem with mapslices is that it's currently slow, but now that we have fast subarrays and fast cartesian iteration I think that can change. (It will never be cache-efficient, however, so the best implementation will be hand-tuned.)

[tomasaschan]

OK - that at least answers the question on what I'm aming for. If you have any resources (now or at a more convenient, less vacation-y, time - just because I like to spend my holidays coding 24/7 doesn't mean you have to 😉) on why that works, I'd be happy to indulge in some reading!

I've started looking at mapslices, but I can't figure out how to write stuff back to the correct place in coefs without doing all the work I'm doing now with range indexing anyway; it seems to me that something like mapslices(A, dim) do v #= cool stuff =# end only gives me "read-only" access to v; for example, this leaves A unchanged:

for dim in 1:ndims(A)
    mapslices(A,dim) do v
        v[:] = 2v
    end
end
A

Quick searches through the docs, as well as the first page of Google results, don't give me any hints. Any ideas?

[tomasaschan]

Huzza! This PR now (usually, see below) passes all tests, including the ones I ported from Grid.jl that don't index with ranges =)

Unfortunately, it only does so most of the time - sometimes I get the weird NaN errors seen in the current Travis build. I have no idea what that comes from (except maybe some odd conditioning problem with the rand initialization of data arrays, possibly similar to timholy/Grid.jl#20 - but this fails much more often than that did...), since one or a few simple re-runs of Pkg.test("Interpolations") usually passes.

Any suggestions on how we can make the test suite less brittle?

[tomasaschan]

There - I've added a bunch of tests that do the same thing as the ones in the old Grid.jl test suite, but without relying on random number generation. However, the results are still not deterministic, as repeated runs of Pkg.test("Interpolations") might show.

A faster way of reproducing this is to open a REPL and do include(joinpath(Pkg.dir("Interpolations"), "test", "on-grid.jl"))) over and over again; then the Interpolations module won't have to be recompiled each time, so running the tests is almost instantaneous. On my machine it fails most of the time, but in different places, and some of the time (perhaps 30%) it doesn't fail at all.

The error message is always something like

ERROR: mismatch of non-finite elements: 
  itp[i,j,k] = NaN
  fg[i,j,k] = 0.0
 in test_approx_eq at test.jl:101
 in ThreeD at /home/tlycken/.julia/v0.3/Interpolations/test/on-grid.jl:69
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
while loading /home/tlycken/.julia/v0.3/Interpolations/test/on-grid.jl, in expression starting on line 85

In every single failure, fg[i...] == 0.0 (or -0.0), while itp[i...] == NaN, which makes me suspect that there is a division by zero somewhere (either in our code - unlikely, since we rarely, if ever, do division except in literals - or in the matrix inversion - no idea how likely...).

@timholy, do you have any idea what might cause this, and/or any suggestion on how to trouble-shoot it since it doesn't seem deterministic?

[timholy]

I haven't had time to look into this in detail. But it occurs to me that if you're getting random failures from what looks like deterministic code:

• You could have a memory error. If you sometimes use @inbounds in your code, try starting julia with julia --check-bounds=yes and run your tests. (I confess I actually did test this, and got a failure right away. So this is definitely a problem.) IIUC this will become the default for Travis tests with the new format of travis scripts.
• It's possible you might be using uninitialized memory, e.g., from Array(Float64, m, n). Occasionally that block of memory will have a NaN in it.

[tomasaschan]

Yep, julia --check-bounds=yes test/on-grid.jl throws a BoundsError first thing that happens. I'll look into it. (And thanks a lot for the tip - I would never have guessed that that could be the problem...)

[tomasaschan]

@timholy, you deserve a medal! The tip on @inbounds was gold - turns out I was indexing out of bounds when evaluating at the upper boundary in the linear extrapolation. I added some edge-case tests for boundary evaluation, to cover both OnCell and OnGrid, so hopefully we're decently covered now.

[timholy]

This is really exciting. Seems like this has gone from "promising package" to "something that's already useful." I'm going to have to start kicking the tires directly, and of course work on the efficiency of the filtering.

[tomasaschan]

Nice! Tire-kicking (and performance enhancements) are always welcome! In the meantime, I've started working on gradient evaluation; I'll file a separate PR for that soon.