Allow multidimensional array indexing with any eltype

[mbauman]

This enables the use of the multidimensional abstract array indexing fallbacks for all element types. In particular, this fixes #12271: A[[]]. It also allows indexing with A[Any[1,2,3]]. I think this is a good direction in general, since we simply index through the array and pass along the elements to a scalar indexing method. If the AbstractArray subtype doesn't support indexing with those elements, then it'll throw an error there. But I think this starts allowing fancier capabilities in packages, like indexing with an array of Dual numbers.

[mbauman]

I suppose doing it this way makes for some less intelligible error messages. We could be a little less lenient and keep the error, and add to_index(i::AbstractArray) = i.

[mbauman]

Ref #10525 (comment) and cc @timholy, @tlycken

[tomasaschan]

Nice! I had some test breakage in Interpolations.jl only yesterday because I tried to do something like itp[1:.1:15], which didn't work out of the box - I suppose this might fix that for me. Will have to take it for a spin... :)

[tomasaschan]

Nope, same error still:

julia> Pkg.checkout("Interpolations") # #10525 isn't utilized in the latest tagged version
julia> using Interpolations
julia> A = sin(1:20)
julia> itp = interpolate(A, BSpline(Quadratic(Flat)), OnGrid)
julia> itp[3:.4:15]
ERROR: MethodError: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
  *{T<:Union{Float64,Complex{Float32},Float32,Complex{Float64}},S}(::Union{DenseArray{T<:Union{Float64,Complex{Float32},Float32,Complex{Float64}},2},SubArray{T<:Union{Float64,Complex{Float32},Float32,Complex{Float64}},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},DenseArray{S,1}})
  *{TA,TB}(::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{DenseArray{TB,2},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},DenseArray{TB,1}})
  ...
 in getindex at /Users/tlycken/.julia/v0.4/Interpolations/src/b-splines/indexing.jl:46

I think the problem is that since we allow indexing with Any, our indexing for ranges don't fall back to the AbstractArray defaults from #10525, and thus tries to treat the index range as a scalar. We don't seem to have done anything fishy, though:

julia> methods(getindex, Tuple{Interpolations.AbstractInterpolation})
3-element Array{Any,1}:
 getindex{T,N}(itp::Interpolations.BSplineInterpolation{T,N,...}, xs...) at /Users/tlycken/.julia/v0.4/Interpolations/src/b-splines/indexing.jl:27
 getindex{T,N,...}(exp::Interpolations.AbstractExtrapolation{T,N,...}, xs...) at /Users/tlycken/.julia/v0.4/Interpolations/src/extrapolation/indexing.jl:9
 getindex(A::AbstractArray{T,N}, I...) at abstractarray.jl:463

What would be the best way to fix this? Do we need to add a method in Interpolations.jl, or should something rather be changed in base?

[mbauman]

I think the problem is that since we allow indexing with Any, our indexing for ranges don't fall back to the AbstractArray defaults

Ah, shucks. Of course. This would require dispatch on complements of Union. I suppose Interpolations could opt in to the fallbacks for Union{Int,AbstractArray,Colon} by explicitly dispatching to the internal undocumented methods... But that's intentionally not a stable API.

The other thing I realized is that you may run into trouble with element types since interpolating doesn't always return the eltype, but the array indexing methods preallocate an output array using eltype.

[tomasaschan]

I suppose Interpolations could opt in to the fallbacks for Union{Int,AbstractArray,Colon} by explicitly dispatching to the internal undocumented methods... But that's intentionally not a stable API.

This doesn't seem like a viable long-term solution, so I'll just hold off on it for now. But it would be awesome if base could learn to handle these things :)

The other thing I realized is that you may run into trouble with element types since interpolating doesn't always return the eltype, but the array indexing methods preallocate an output array using eltype.

We've gone quite some length to make sure that the T in the AbstractArray{T} type each interpolation type inherits actually corresponds to the type returned by indexing into the interpolation object if you're indexing with "reasonable" number types. We've also provided methods for interpolate where you can provide the T yourself, if you know you're going to index with some esoteric number type and want to make sure things Just Work(TM), so if users run into problems here I think we've catered for their needs too, even if it will require a small effort on their part.

[timholy]

I think the problem is that since we allow indexing with Any, our indexing for ranges don't fall back to the AbstractArray defaults

In Interpolations, if we want to make use of the fallbacks in #10525 we should probably restrict the index type to Number.

I approve of this PR.

[mbauman]

I want to simplify things here a bit further before this gets merged.

[mbauman]

Ok, this is much better now:

julia> A[[]]
0-element Array{Int64,1}

julia> A[im]
ERROR: indexing Array{Int64,3} with types Tuple{Complex{Bool}} is not supported
 in error at ./error.jl:21
 in getindex at abstractarray.jl:463

julia> A[[im]]
ERROR: ArgumentError: unable to check bounds for indices of type Complex{Bool}
 in getindex at abstractarray.jl:463

julia> A[Vector[[1,2],[3,4]]]
ERROR: ArgumentError: invalid index: Array{T,1}[[1,2],[3,4]]
 in getindex at abstractarray.jl:463

We're still not quite at the point where we can support any element types for non-scalar indexing. The goals of scalar and non-scalar indexing fallbacks are slightly different.

Scalar indexing:

1. Convert given indices to Int with to_index
2. checkbounds
3. Use sub2ind or ind2sub and call getindex with all Ints.

Non-scalar indexing:

1. Call checkbounds first so ranges and logical vectors are optimized
2. Iteratively call scalar getindex with the elements of the index arrays

For non-scalar indexing, we don't really care what the element types of the index arrays are... except for that hoisted checkbounds call. I'm still not sure how to best deal with that, but for now custom arrays could potentially override checkbounds for checking custom indexing behaviors.

Anyhow, that's an aside and a project for another time... and really just notes for my future self.

[timholy]

Doing the check before calling to_indexes is better-optimized for logical indexing.

[mbauman]

Yup, but this is scalar indexing so I<:Tuple{Vararg{Real}}. It's something that I think we can get rid of once float indices are no longer supported, but that will take some thought about what to do with non-Int integers. We'll still want to convert small unsigned integers to Int before doing sub2ind arithmetic... and ideally also allow larger integers, too.

[timholy]

Good point---I was still remembering the original thinking behind this code and not noticing what it had become. Great work as always.

[JeffBezanson]

Ready to merge?

[mbauman]

Yup!