Vector-valued evaluation

[tomasaschan]

Looking at #5, this is one of very few things left to do before we might call ourselves ready for a version 0.1.

I still think this makes sense only in a multi-linear context; the following is probably a good method definition for that:

getindex{T,N}(itp::Interpolation{T,N}, NTuple{N,AbstractVector}...)

Of course, it's going to have to be more specific in order to avoid ambiguity errors. I'm thinking we might not have to specify the type of the vector elements at all - we'll get method errors in the next stage anyway, if the user supplies some nonsensical index vector.

[timholy]

Agreed.

This is really exciting!

[tomasaschan]

Is there any way to do this into a pre-allocated array? (Do we want to?) If so, what should the method signature look like?

[timholy]

It could look like this:

getindex(itp, I::AbstractVector...) = getindex!(Array(eltype(itp), map(length, I)), itp, I...)
function getindex!(dest, itp, I::AbstractVector...)
    # real function goes here
end

There are some models of this in base/multidimensional.jl.

[tomasaschan]

Nice! I ended up using a promotion rule like this to define the array:

eval(ngenerate(:N, :(Array{T, N}), :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Any}...)),
    N->quote
        Tout = promote_type(T, map(eltype, @ntuple $N x)...)
        dest = Array(T, map(length, @ntuple $N x)...)
        getindex!(dest, itp, x...)
    end
))

However, I couldn't use the same expression to define the return type (that's why the above says Array{T,N} rather than Array{<expr for TOut>,N}) because N doesn't seem to be defined there. Does it matter that the declared return type for ngenerate isn't exactly correct?

Also, when trying to resolve all the ambiguity errors, I want to ngenerate methods for N>=2 to disallow linear indexing into multi-dimensional interpolation objects, but I can't figure out a way to accomplish that. Basically, I'd like to do something like

@nsplat 2:Base.Cartesian.CARTESIAN_DIMS getindex{T,N,...}(itp::Interpolation{T,N,...}, x::Union(Real,AbstractVector)) = error("...")

but of course @nsplat requires that the last argument is a splat (and it doesn't support that range either; I have to say explicitly 2:4, but I can live with that for now...). Any suggestions?

[timholy]

You can cascade calls:

getindex(A, I...) = getindex_typed(A, compute_return_type(A, I), I...)
@ngenerate N Array{S,N} getindex_typed{T,S,N}(A::Array{T,N}, ::Type{S}, ...)

[tomasaschan]

Leaving that here (as a mental note to myself, mostly, so I can see if I can take a stab at this if I get some quality Julia time this weekend :smile: ):

With julialang/julia#10525 this is basically done - the only thing (except for tests!) missing is range indexing, which can probably be worked around by requiring that indices inherit Number (as you mentioned in https://github.com/JuliaLang/julia/pull/12273#issuecomment-124060763).

[deszoeke]

+1 for vector indexing (and stopgap warnings that it's not supported yet). In Julia 0.5.0 or 0.5.1 it seems to work for an interpolator from interpolate(), but not one from extrapolate(itp, missingvalue). Is this possible with some kind of broadcasting of xi within getindex?

Here's the test:

x=[0:10:100 200:10:300][:]
y=[2:2:22 2:2:22][:]
itp=interpolate((x,),y,Gridded(Linear()))
xi=-55:10:355
yi=itp[xi] # works
itpne=extrapolate(itp,-99) # non-extrapolating
yi=itpne[xi] # doesn't work: yi=getindex(itpne,xi)
[yi[i] = itpne[xi[i]] for i=1:length(xi)] # works