Closed timholy closed 10 years ago
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?
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:
body_gen
that takes care of evaluation outside the data domain.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.(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 :smile:, 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.
@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
.
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...)
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 =)
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.
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.
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 https://github.com/timholy/Images.jl/pull/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.
Looks like it got sorted based on date, so you'll have to scroll up to find it.
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.
Woot! Very exciting. Thanks for doing this!
Now that Interpolations.jl exists, should we close this PR and continue discussion there?
Closing in favor of Interpolations.jl.
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 sayingIn 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.