Rewrite with clear separation of concerns (closes #2)

[tomasaschan]

I just started working on some documentation for how I envision the pieces of this library to fit together. Most of this is obviously just sketches and fairy tales at the moment, but at least it might give us something to discuss and improve on. I'm hoping that discussions around these documents can help finalize some decisions that are probably easier to make before writing macro-heavy code than after =)

[tomasaschan]

CC @timholy I don't know if you get notifications on this without me mentioning you, but I'd love to hear your thougths on this =)

[timholy]

I do get notifications (you can see who it watching by clicking on the 5 up at the top of this page).

Don't have time to read it right now, but hope to soon.

[tomasaschan]

Check! No hurries =)

In the meantime, I'll see if I can come up with a skeleton implementation using this architecture, so there is some code to look at as well.

[tomasaschan]

OK, there's a very simple example of the architecture I have in mind over in the rewrite branch

[timholy]

It's looking more put-together, which is nice! But I'm still not sure I understand the mathematics you're proposing for extrapolation. If that's easier to write up as a LaTeX document, go for it.

Presumably you also saw timholy/Grid.jl#53, which could be turned into a well-defined proposal for how to extrapolate.

[tomasaschan]

There - there's now a working implementation of constant and linear interpolation, as well as quadratic interpolation in 1D. The results look OK to me: this is one period of sin(x) over 10 datapoints:

When looking at the prefiltering implementations in Grid.jl, I didn't really like all the custom slicing and iteration there, mainly because they make the code more opaque to new readers. What, exactly, should the prefiltering do for higher dimensions? Could we use some smart metaprogramming approach to get there?

Also, I don't know if my thoughts about extrapolation became clearer now that there is some substantial code to look at as well. If not, please let me know and I'll try to clarify further.

[tomasaschan]

Oh, there's also an IJulia notebook which shows off the way to reproduce that graph.

[timholy]

Yes, I agree that the custom iteration has to go. Now that Julia has Base.Cartesian (and in 0.4, easy-to-use iteration types), I suspect this should be quite feasible.

By prefiltering you mean the whole interp_invert! stuff?

[tomasaschan]

Yeah, exactly. I basically copied the relevant interp_invert! implementation to get the quadratic splines to work in 1D, but I didn't invest any time in ways to do it without the custom iteration (I'm still not that familiar with all the functionality in Base.Cartesian...).

interp_invert! in Grid.jl corresponds to prefilter/prefilter! in Interpolations.jl, and there is an example implementation in src/quadratic.jl.

[timholy]

I don't think you committed src/quadratic.jl.

[tomasaschan]

I didn't initially - sorry about that - but I fixed it later. So now it's here 😄

[timholy]

Ah, I was looking at the rewrite branch, not the doc branch.

While I admire the desire to make the inversions easier, I have to warn you it's nontrivial. Basically, the task in 2d is that you have to solve a tridiagonal system along each column, and then along each row. In 3d, think about a bundle of pencils all pointing in the same direction: each pencil represents a tridiagonal system of equations. Solve each pencil, then re-orient the entire bundle along a different axis and do it all over again. So if you have a KxMxN array, you solve M*N tridiagonal systems for the first axis, then K*N for the second axis, then K*M for the third.

The good news is that we have better tools now. Particularly, the combination of stagedfunctions and Base.Cartesian should make this more straightforward.

BUT...it's not just that supporting arbitrary dimensions is a bit tricky; if we want it to be really fast, we may be in for a certain amount of pain. The version in Grid basically "snips out" a pencil (the Slice1D) and works with that pencil in isolation---that allowed me to use the tridiagonal/Woodbury solvers as black boxes.

Unfortunately, in terms of the order in which it visits elements, that's analogous to the version of the algorithm described under "write cache-friendly codes" introduced as "It is natural to implement this as follows:". (You could have replaced the inner loop in Dahua's illustrative example with a black-box sum_over_row function, similar to how I'm using the tridiagonal/Woodbury solvers.) Instead, it would be desirable (from a performance perspective) to basically incorporate the tridiag/Woodbury solver into the loop itself, analogous to the algorithm that follows "...by changing the order of computation". (In that second example, there is no way one could replace a loop with a sum_over_row function, because it doesn't work row-by-row.) So rather than calling out to a tridiagonal/Woodbury solver, we might basically have to write one. Or probably better, write a stand-alone support package to abstract away the ugliness of this step.

I don't know how much of this is clear; I can try to explain any specific parts that are confusing.

[tomasaschan]

The pencils are a good analogy, so most of it is pretty clear. I too realize that this is very much nontrivial, so I guess I wasn't hoping for a quick "add a macro here" solution, but I think we could do this in several steps.

First, obviously, we should to build something that works at all. I guess a naïve, cache-unfriendly version that might be slow but gets the job done is definitely good enough. A lot of work from Grid can probably be re-used here, but using Base.Cartesian and other tools rather than custom-built slicing and iteration.

When we have that in place, we can start thinking about ways to make it faster. A package with general-dimensional, cache-friendly versions of Woodbury and friends would be super cool, but it's definitely outside of my comfort zone (not very far outside, but still outside - if I do it, it's going to take a while...), so maybe it's something we should post for help for on the Julia users list or something.

But since all of this matrix system solving is done only once for each interpolation object, I think it's more important to get indexing (and, eventually, also evaluation of derivatives) fast, and to implement as many useful interpolation types, boundary conditions etc as possible.

[tomasaschan]

After giving this another round of pondering on the subway ride to work, I think my main concern here is that I want to separate implementation details from the math as much as possible also in this part of the library, in the same way as we do with generation of the getindex methods and the various expressions for indexing, boundary conditions etc. Basically, we could separate the prefiltering steps into three distinct parts:

• Setting up the un-altered tridiagonal system matrix, that depends only on interpolation degree and grid representation
• Altering the tridiagonal system matrix to take boundary conditions into account (what Grid does using Woodbury matrices)
• Actually solving the system (in all the relevant directions)

I think the basic approach should be that when implementing a new interpolation type, the coder just provides an expression for the un-altered system (or re-uses one, if there is one already), and another expression for the alterations, and then there is some macro magic that takes those and, like for getindex, generates the relevant prefilter methods which do the actual solving. That way, for example implementing new boundary conditions for an existing interpolation degree entails as little copy-pasting of existing code as possible. That way, we can also improve the performance of the code across the board by giving the macro magic more unicorns and rainbows.

[tomasaschan]

There, now I've generalized this to use macro-ized prefiltering.

However, there is one problem: because of how I index into the coefficient matrix, the prefiltering isn't applied. If I understand how this works correctly, vec(A[1:stride(A,dim):size(A,dim)*stride(A,dim)]) from here makes a copy, so A_ldiv_B! doesn't actually change anything in A.

Is this behavior going to change in 0.4 (I'm testing on 0.3 atm...)? If not, how should I do this indexing to actually alter A in-place?

Also, I realize that this is still not cache-friendly. However, I don't think there is a way to do this cache-friendlier, at least not without copying a lot of data between runs, since we want to do the matrix solving in all dimensions (thus being cache-un-friendly in some dimension regardless of which ones we start in...). If the matrix solving algorithm has a lot of cache misses, it might be benefitial to copy out the vector that we want to solve with, solve, and then copy back the result, but that will impose a copying cost instead, so I'm not sure what will be fastest. If this doesn't cut it, we'll have to profile and decide then.

[tomasaschan]

...and of course I had simplified the problem too much. That code, when altered to have the correct copy behavior, only works in 1D - anything higher, and there is no iteration over each row and column anymore, so most of the array is never filtered. I'm looking into it...

[timholy]

On Julia 0.4, with the new SubArray work I just merged, I think it will become possible to rewrite julia's mapslices function in a way that doesn't have horrible performance. (It will still not be ideal, but much better than the disaster it is now.) This might provide a simple way forward, although likely quite imperfect from the standpoint of performance (likely worse than what we have now).

[tomasaschan]

Voila! As far as I can tell, this works now for at least 1D and 2D interpolation (and none of the implementation is dimension specific, so I would assume that it carries to larger dimensions). I also re-worked the notebook into a sort-of tutorial that shows off and explains the currently implemented functionality.

I haven't yet figured out a nice way to test all this stuff, so Travis will still complain. I'm thinking about stealing relevant parts of Grid.jl's test suite and adapting the interpolation constructors.

@timholy, how do you feel about merging this and start building in smaller chunks from here?

[timholy]

Why make a copy when, for nearest and linear, it's not necessary?

[tomasaschan]

For concistency.

Copying is necessary for higher orders (if we don't want to change the input array already on creation of the interpolation object...), so if we don't copy also for nearest and linear, interpolation objects of different degrees won't behave the same:

A = rand(10)
itp1 = Interpolation(A, Linear(OnGrid()), ExtrapError())
itp2 = Interpolation(A, Quadratic(ExtendInner(), OnCell()), ExtrapError())

A[3] = 5.0

Now, the coefficients of itp1 have changed, but not the coefficients of itp2...

[timholy]

I see your point. But for efficiency we might want to offer Interpolation! that avoids making a copy (and in the quadratic case destroys the input).

[timholy]

Overall, this is looking really promising!