WIP: metaprogramming-based interpolation

[timholy]

OK, here's a prototype of the kind of thing I had in mind with #36. You can test by checking out the cartesian branch and then, from the test/ directory saying

julia> include("newinterp.jl")
BCnil
elapsed time: 0.832196552 seconds (16 bytes allocated)
elapsed time: 0.212113778 seconds (16 bytes allocated)
elapsed time: 1.359122945 seconds (16 bytes allocated)
elapsed time: 0.422922048 seconds (16 bytes allocated)
BCnan
elapsed time: 0.902566273 seconds (16 bytes allocated)
elapsed time: 0.248787426 seconds (16 bytes allocated)
elapsed time: 1.421108665 seconds (16 bytes allocated)
elapsed time: 0.461047829 seconds (16 bytes allocated)
BCna
elapsed time: 1.070271096 seconds (16 bytes allocated)
elapsed time: 0.197870166 seconds (16 bytes allocated)
elapsed time: 1.362531069 seconds (16 bytes allocated)
elapsed time: 0.334651534 seconds (16 bytes allocated)
BCperiodic
elapsed time: 1.345538871 seconds (16 bytes allocated)
elapsed time: 0.435039454 seconds (16 bytes allocated)
elapsed time: 2.308907927 seconds (16 bytes allocated)
elapsed time: 0.728555211 seconds (16 bytes allocated)
BCnearest
elapsed time: 0.812522817 seconds (16 bytes allocated)
elapsed time: 0.187141129 seconds (16 bytes allocated)
elapsed time: 1.222788721 seconds (16 bytes allocated)
elapsed time: 0.359194482 seconds (16 bytes allocated)

In each case, old is on top and new is on the bottom. The first pair is for 1d, the second for 2d.

The improvement is not as large as I expected, but it's still substantial. Bounds-checking takes more time than the actual computation, at least in 1 & 2 dimensions. (You can profile to see that.)

I didn't implement all the boundary conditions, but only because I'm out of time for working on this further. I don't foresee any obstacles, though.

[coveralls]

Coverage remained the same when pulling cf01938 on cartesian into a099b87 on master.

[timholy]

There's one important thing to discuss: interpolation of multi-valued functions. Let's say I had an RGB image but for some reason didn't want to use the RGB immutable. Then each spatial location is 3-valued. The original implementation somewhat awkwardly allowed you to set the position once and then calculate the interpolated value re-using the same set of coefficients. Since this version doesn't define any temporary storage, that becomes more awkward.

We could reintroduce the storage, but I also wonder if there is a nicer API for this. For example, G = InterpGrid(A, BC, IT; value_dims = 1) might specify that G[x_1, x_2] would effectively calculate A[:, x_1, x_2]. We could define a val! method that stores the results in a properly-sized array output, so you aren't forced to allocate memory to return multi-valued outputs.

Thoughts?

[coveralls]

Coverage remained the same when pulling fc8a380 on cartesian into a099b87 on master.

[tomasaschan]

Huh. Cool. And somewhat magical...

I must admit that this code is even more opaque to me at first sight than the previous version - especially since so much happens in contexts where it's not always straightforward to see what various variables represent. (Stepping through code in my head is difficult enough - doing it through code that I first had to generate in my head makes it more difficult. Not very surprising... :P) Documentation and some examples in comments, once we're done, will probably help this quite a lot.

However, I do find it very elegant how it seems like every combination of interpolation order and boundary condition can be expressed with a body_gen method only. As long as we make sure to document thoroughly what body_gen is supposed to do, this probably makes it quite easy to extend the library with more cases (e.g. higher interpolation orders or new boundary conditions).

A couple of comments:

• I still want to separate the concepts of boundary conditions and extrapolation behavior (see Decision: Separate the concepts of boundary condition and extrapolation behavior #37). I think a good way to do this might be to have an abstraction in body_gen that takes care of evaluation outside the data domain.
• Regarding multi-valued interpolation, I'm wondering if it's not possible to do via a type parameter instead; today T<:FloatingPoint - couldn't we allow also e.g. T<:Vector{S}; S<:FloatingPoint? We could probably duck-type it based on what operations are needed to evaluate the interpolation. My main concern is that a keyword argument will a) violate type stability (maybe avoidable, if we're careful), and/or b) seem like even more magic than what's already happening in the macros. Given the speed increase, it also seems like evaluating the interpolation three times will still be no significant performance loss to what we did before, so maybe that's also an option. I haven't really worked with images, or any other context where this type of interpolation is used, so my intuition for it is quite limited.

[timholy]

(Stepping through code in my head is difficult enough - doing it through code that I first had to generate in my head makes it more difficult. Not very surprising... :P) Documentation and some examples in comments, once we're done, will probably help this quite a lot.

Yeah, this is the big negative of this approach. We don't have to go this way; it's more important that this works for those nice people who help make Grid better 😄, even if we have to give up some on performance. I'm willing to accept your guidance on where the right balance is.

One tip, though: at least at first, instead of generating the code in your head, try calling (for example) body_gen(BCnan, InterpLinear, 2). It will generate the code for you, and print it out in the REPL. macroexpand is also really helpful for understanding what the Cartesian macros do.

I still want to separate the concepts of boundary conditions and extrapolation behavior

Absolutely. This was merely to give you a sense for what I was thinking about, upon which you might consider basing your own extensions.

My main concern is that a keyword argument will a) violate type stability (maybe avoidable, if we're careful)

Hm, that is interesting. I hadn't thought about whether one wants to encode the "spatial" dimensionality in the type. I'll have to think this through more carefully.

[coveralls]

Coverage remained the same when pulling 3844de1 on cartesian into a099b87 on master.

[timholy]

@tlycken, are we good to go with this approach? I'm thinking of punting on the multi-valued issue for now.

There's another approach we could take: since the low-level API will be quite different, and that low-level API is considered "public," we could just start a new Interpolate.jl package and slowly deprecate Grid. Or have Grid just be for restrict/prolong.

[tomasaschan]

Yeah, sure - I'm willing to learn what I need to understand the meta-programming approach, and it does look like it could provide a better framework.

However, I like the suggestion to deprecate Grid even better. We've already talked about before that restrict and prolong might be better suited for an image processing package (but Images.jl already has at least restrict, right? I know you mentioned somewhere that the algorithms were different and both had their merits, but I can't find where now...) and the first time I wanted interpolation stuff I had to ask on the mailing list because Grid wasn't an obvious name for it. Interpolate (or Interpolations) is a much better name.

If we don't introduce any big changes to Grid.jl now, we should be able to keep it stable and "good enough" to still be the go-to package (at least for now) when Julia 0.3 is released, and then have Interpolations be the go-to package for 0.4 and up. (It's not a good idea to deprecate Grid until a feature-worthy replacement is available...)

[tomasaschan]

Also, OT but maybe relevant for the implementation of this: I'm going to be afk with no internet connection for most of August, so don't hold your breath for replies from me if you don't have very big lungs =)

[timholy]

Sounds like a fun time!

OK, I'll start a new package. I'm immersed in other things at this moment, so this will unfold over time. Thanks for getting back to me.

[tomasaschan]

I'm back from the bush, and have started taking a look at this a little more closely. I think I'm getting the hang of how linear interpolation works in 1D, but I can't seem to wrap my head around the general principle for how multidimensional interpolation maps to recursion. (I've started reading this paper, but it's a long read...) For the linear case it's quite straightforward, but for e.g. quadratic or cubic B-splines it's not at all clear to me how (or even that) it works.

Anyway, it seems like a good approach, given that it can be generalized to higher-order interpolation, so I'm all for continuing with it. I think maybe the structure could be somewhat different to make new additions easier - and maybe that will help finding some more descriptive function names as well... - I'll see if I can come up with an example for (at least) linear interpolation.

[timholy]

Welcome back! Hope you had a blast.

To be honest, I haven't thought about higher-order interpolation at all. I've just assumed it will work similarly.

I'm currently utterly swamped by JuliaImages/Images.jl#135, so it will still be a little while before I can get to this. I did remember that I have a commit I never pushed (one that fleshes out the Linear case more completely). I'll push now.

[timholy]

Looks like it got sorted based on date, so you'll have to scroll up to find it.

[coveralls]

Coverage remained the same when pulling bfa96fb on cartesian into a099b87 on master.

[tomasaschan]

Nice!

I took the liberty of stealing some of your code, rewriting some of it quite a lot, and re-packaging it in a new (yet unregistered) package Interpolations.jl. It does almost exactly the same thing as yours, but it has the additional benefit of separation of boundary conditions and extrapolation behavior.

[timholy]

Woot! Very exciting. Thanks for doing this!

[tomasaschan]

Now that Interpolations.jl exists, should we close this PR and continue discussion there?

[timholy]

Closing in favor of Interpolations.jl.