Interpolating over a vector of non-integer values raises "ERROR: InexactError()"

[btengels]

Is there a reason why I can vectorize interpolation when using a Gridded Interpolation object and not a plain vanilla B-spline?

using Interpolations

#----------------------------------------------------------------------- 
# B-spline (doesn't work)
#-----------------------------------------------------------------------
# make an array of floats and an interpolation object
J = broadcast(+, collect(linspace(2,5,10)), collect(linspace(1,5,10))' )
J_itp = interpolate(J, (BSpline(Linear()), NoInterp()), OnGrid())

# attempt to interpolate over a vector 
J_itp[linspace(1,2,10),1]
# J_itp[1.5,1]  # a single point works fine 

#-----------------------------------------------------------------------
# works on gridded interpolation just fine though...
#-----------------------------------------------------------------------
knots = (collect(linspace(2,5,10)), [1:10])
J_itp = interpolate(knots, J, Gridded(Linear()))
J_itp[linspace(2,3,10),1]

The interpolation call gives me the following error:

ERROR: InexactError()
in to_index at deprecated.jl:454
in _unsafe_getindex at multidimensional.jl:192
in getindex at abstractarray.jl:488

Great project btw. I hope to help out somehow though I'm not really a julia expert.

[sglyon]

@btengels great to see you here!

The only reason is that it hasn't been implemented.

For gridded interpolation there is a significant efficiency gain if we can cache the indices. The gain isn't as large with with cardinal ("vanilla") B-splines.

I have plans to implement vector evaluation for cardinal B-splines because you can actually get decent speedups if you make certain assumptions. One assumption is that you assume the points you want to evaluate at are sorted or otherwise close to one another. Then you can avoid a whole binary search across the vector by starting with a guess for the index at point i being near the index at point i-1.

[timholy]

With a Gridded object, there's a big performance savings to vectorization, because you have to do a grid-lookup for each evaluation and you can effectively cache the results. With a plain BSpline, there's no performance advantage in vectorization, and it seems that we were hoping to use julia's fallbacks: for most AbstractArray types, you only have to write routines for scalar indexes, and julia takes care of the rest. However, that seems to require that the indexes are integers, which obviously doesn't work with Interpolations.

So, someone presumably needs to implement the vectorization within Interpolations. It should be pretty simple, e.g., the one-liner

getindex(itp, index1, index2) = [itp[i1,i2] for i1 in index1, i2 in index2]

work for 2d. But obviously it has to be done for arbitrary dimensionality (using generated functions, presumably).

[tomasaschan]

If it can be done, I'd really like to "force" the vectorized fallbacks in Julia to become compatible with this* instead of implementing it here. For one thing, the code here is going to be very similar to the one in base Julia, and for another I think this package is a good use case for testing that the fallbacks become flexible enough to actually be useful.

There might be some perfomance optimizations for B-splines that are still possible after that, but I doubt that any of them are worth the added complexity, and I definitely think we should wait with implementing until we can measure that it's needed.

In the meantime, @btengels, I just do [J_itp[x,1] for x in linspace(1,2,10)] manually whenever I need something like that.

*) although it's a little over my head to implement it, so I'm just hoping @timholy and @ mbauman (who doesn't need to be pinged here...) will sort it out 😄

For reference, see also #24, #54, #55, #96 for other discussions related to vector-valued evaluation.

[timholy]

*) although it's a little over my head to implement it, so I'm just hoping @timholy and @ mbauman (who doesn't need to be pinged here...) will sort it out 

We will. But it will only arrive in julia 0.5.

[tomasaschan]

Not supporting itp[linspace(1,10,100), 3.5] until Julia 0.5 arrives is fine by me - after all, it's really only nicer syntax for functionality that already exists in more verbose forms (list comprehensions). If others disagree strongly there might be merit in doing something about it, but MHO is that the added value isn't worth the added complexity. Interpolations.jl's source is already pretty complex as it is.

[tomasaschan]

Closing this for housekeeping - feel free to re-open if there is more to discuss here.