Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

More boundary conditions for quadratic B-splines #9

Merged
merged 24 commits into from
Dec 30, 2014
Merged

Conversation

tomasaschan
Copy link
Contributor

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 ].

linear-bc

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.

@tomasaschan tomasaschan mentioned this pull request Nov 24, 2014
17 tasks
@timholy
Copy link
Member

timholy commented Nov 24, 2014

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Nov 25, 2014

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 2, 2014

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 2, 2014

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 2, 2014

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 2, 2014

Thanks for the link.

@tomasaschan tomasaschan changed the title First attempt at linear bc/extrap for B2 More boundary conditions for quadratic B-splines Dec 26, 2014
Tomas Lycken added 7 commits December 26, 2014 06:12
The extra equations from the boundary conditions will not always allow
indexing as we did before, by e.g. doing cm_d = cp_d at the lower boundary
for some BC:s. In order to generalize better, we instead pad the equation
system with one line for each extra equation to allow more general BC specs.

This had the advantage of significantly simplifying some of the code (and,
even better, some of the concepts) required to implement new functionality,
so even though some of the underlying infrastructure is now more complicated
than before, new contributors will hopefully have a lower threshold to start
adding features.
@tomasaschan
Copy link
Contributor Author

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 ci-1 = ci+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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 27, 2014

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
Copy link
Member

timholy commented Dec 27, 2014

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 29, 2014

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Dec 30, 2014

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
Copy link
Contributor Author

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
Copy link
Contributor Author

@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.

tomasaschan pushed a commit that referenced this pull request Dec 30, 2014
More boundary conditions for quadratic B-splines
@tomasaschan tomasaschan merged commit 9820504 into master Dec 30, 2014
@timholy
Copy link
Member

timholy commented Dec 30, 2014

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
Copy link
Contributor Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants