Low level interface help/documentation

[iancze]

I'm extremely happy to have found Grid.jl, it seems like the perfect tool for doing the type of interpolation that I need. However, I am having a bit of difficulty getting started for my specific example. While combing the source of interp.jl has gotten me pretty far, I still can't quite get the behavior I want.

My main data unit is a 1d spectrum, with length of about 200,000 pixels. There are three parameters which describe a spectrum (call them x, y, and z), and I have about 500 spectra which are produced for regularly spaced combinations of these three parameters. My problem is analogous to the image/RGB example in the documentation, but in this case instead of the two spatial dimensions x and y and a third color dimension c (with size 3), I have three dimensions (x, y, z) and my color dimension is size 200,000.

Therefore, what I actually want to interpolate is a 1d vector with length 200,000 (i.e., a spectrum) which has been interpolated over the x, y, z axes but not the color axis. My attempt is the following, but I cannot figure out where to properly limit the dimensions to avoid interpolation over the color axis.

grid #my array of spectra on a grid, with size (11,8,6,216631), ie axes (x, y, z, c)
interp_invert!(grid, BCnan, InterpCubic, 1:3)   # solve for generalized interp. coefficients
ic = InterpGridCoefs(grid, InterpCubic) #gives error

#Then, suppose I want to interpolate a spectrum with values x=1.5, y=1.5, z=1.5
set_position(ic, BCnil, false, false, [1.5, 1.5, 1.5])  # set position to a point in the grid, computes the coefs

spec = Array(Float64, (Npix,))
for i=1:Npix
    spec[i] = interp(ic, grid[:,:,:,i])
end

However, different combinations of setting dimensions and limiting the interpolation have been unsuccessful for me. Although the README alludes to the ability to do interpolation like this with RGB channels, it would be helpful if the documentation included a more detailed example. If someone would gratefully help me solve my interpolation problem, I would be very happy to submit a pull request with a detailed example.

Thank you for your help!

[timholy]

Hmm, you're right that this isn't well documented. I suspect that the construct you need is probably something like

index = 1
for i = 1:Npix
    spec[i] = interp(ic, grid, index)
    index += strd
end

where strd is the stride along the dimension you want to interpolate, e.g., strd = stride(grid, 4) in your case.

Can you test and see if this produces correct answers? (Linear interpolation might be the easiest case to check.) I'd be grateful for a detailed example in the README!

[iancze]

Thank you for the quick response! My current stumbling block actually seems to be the choice of dimension indexes in the first few setup calls, and my general confusion of when to pass the full 4d array and when to pass only 3d slices.

interp_invert!(grid, BCnan, InterpLinear, 1:3)   # solve for generalized interp. coefficients
ic = InterpGridCoefs(grid, InterpLinear)    # prepare for interpolation on this grid

My script currently exits on the ic = call with ERROR: Dimensionality mismatch. I noticed in interp.jl that there are additional dims and strides arguments to the InterpGridCoefs function, and so I also tried

ic = InterpGridCoefs(grid, InterpLinear, [11, 8, 6], [1, 11, 88])

where [11, 8, 6] and [1, 11, 88] are the dims and strides in Vector form, respectively, but I get the error

ERROR: `InterpGridCoefs{T,IT<:InterpType}` has no method matching InterpGridCoefs{T,IT<:InterpType}(::Array{Float64,4}, ::Type{InterpLinear}, ::Array{Int64,1}, ::Array{Int64,1}) 

This seems like a dispatch error. I think I am calling the function correctly, but I might be making a silly error. Apologies if this is something trivial that I am missing. Thanks!

[timholy]

The first argument is a type, not an array. Try

ic = InterpGridCoefs(eltype(grid), InterpLinear, [11, 8, 6], [1, 11, 88])

[iancze]

That did the trick! Thanks for the help. I added some text to the README in #45

[iancze]

Now I have the basic functionality to accomplish the sort of task that I need to do, and I am very thankful for that! However, my curiosity has gotten the better of me... I started digging into Grid.jl a bit more, and now I have some more questions :)

Non-interpolating B-splines: I've been boning up on my interpolating functions (Numerical Recipes, A Practical Guide to Splines by de Boor, etc) and it seems like there is a bit of overloaded terminology that is confusing me. For example, I've seen reference to "cubic splines" (that can be interpolating, ie, go through grid points, such as scipy's version ) and "cubic B-splines" that apparently do not necessarily go through grid points. For example, check out my interpolation test here. This is the interpolation of the value at a single pixel of my spectrum, for many different values of x1 (keeping x2 and x3 fixed). The black dots are the grid values and the blue are the interpolated values. To me, it seems like Grid.jl does a pretty good job interpolating through grid points! What is the distinction between non-interpolating v. interpolating here, and might there ever be functionality for the user to choose the behaviour (such as with scipy's other spline function)?

The second part of my question is comes up when I try to interpolate near the edges of the grid, but still within the domain. Same idea as before but now I've moved closer to the grid edges here.

Boundary Conditions: Looking through the previous issues (#31, #29) it seems like the BC's for InterpCubic are something to be mindful of. As a rule, should I ignore any interpolation that is less than three grid points from the edge of the domain?

If there's a major rewrite occurring I can to try and contribute with what limited knowledge/ code foo I can. Thanks again.

[tomasaschan]

@iancze Regarding interpolating and non-interpolating splines, I found this paper (unfortunately not available for free) very enlightening (thanks Tim for showing it to me!). One of the main conclusions of the paper can be summarized as follows:

When considering basis functions for interpolation, some basis functions (like linear splines, i.e. linear interpolation) are "intrinsically" interpolating, i.e. they satisfy

where f_k are the datapoints. For f(x) to pass through all data points, we must require that the basis functions vanish for all integer values of x-k except when x-k=0, where the basis function must be one.

In cases of basis functions that don't fulfil this requirement, we can either a) accept that the interpolated function f(x) does not pass through the data points, or b) solve a linear system of equations for a new set of coefficients c_k, and instead use the interpolation expression

where c_k, as mentioned above, can be found by solving a linear system of equations only dependent on f_k (the given data points) and phi(k) (the value of the basis function for integer arguments).

Grid.jl takes the second approach, and solves the linear system, so there's currently no (easy) way to choose a "non-interpolating" version (but why would you want to? :wink:); all interpolations done by Grid.jl are interpolating.

Regarding boundary conditions, I haven't had time to debug the behavior of InterpCubic to see what BC:s are working correctly and what are not, but the behavior in your image could very well be intended - it depends on what boundary condition you chose (a quadratic termination could show that behavior, for example). In its current incarnation, boundary conditions are also lumped together with extrapolation behavior, which makes these things a little more difficult to reason about efficiently.

There is in fact a major rewrite under way (see #38), but it will probably end up in a different package (at which point we will, eventually, deprecate Grid.jl). Grid.jl will definitely be around for at least as long as Julia 0.3 though, so if you want to contribute (e.g. by trying to fix #29 and/or #31) it would be most welcome =) Just ask questions in those issues (or new ones) if/when you need to, and Tim and I will try to answer as good as we can.

[tomasaschan]

Hm... the images with equations worked in the preview, but apparently not when posting...

The first one was supposed to say f(x) = \sum_k f_k \varphi(x-k) and the second one f(x) = \sum_k c_k \varphi(x-k). I hope you can compile LaTeX in your head ;)

[iancze]

Thank you very much for the help and quick turn around! Your explanation and the paper you provided cleared up a lot of the confusion in my head. I agree with you, every time I reach for an interpolating function, generally my first requirement is that it actually reproduce my grid :)

Also, thank you for the advice on boundary conditions. The behavior I most desire would be that the function connects to the end grid point with minimal oscillatory behavior. If this requires a linear interpolation for that segment then that's fine with me. I think I am able to complete my current project with the current behavior and simply by staying far enough away from the grid edges, but I may find myself digging into #31 :) It seems as though understanding the guts of Grid.jl might be a prerequisite for working on Interpolations.jl?

[timholy]

Interpolations.jl will have a very different set of guts, so as long as you know the math, there's probably only minor advantage to knowing the internals of Grid.

[tomasaschan]

@iancze I actually started doing some work on Interpolations.jl today, so you can take a look at the repo and see for yourself what the code looks like =)

As the readme says, it is very much a work in progress, but I'm hoping it will have feature-parity with the interpolation capacities of Grid.jl and be ready for release in time for Julia 0.4. And of course, all help is welcome =)