RFC: Give AbstractArrays smart and performant indexing behaviors for free

[mbauman]

This is still a work in progress, but I'd like to get some feedback on the architecture and design here before applying this same sort of scheme to setindex!.

The basic idea is that the only getindex method defined in base for abstract arrays is getindex(::AbstractArray, I...). And the only methods that an AbstractArray subtype must define are size and just one getindex method:

getindex(::T, ::Int) # if linearindexing(T) == LinearFast()
getindex(::T, ::Int, ::Int, #=...ndims(A) indices...=#) if LinearSlow()

Unfortunately, it is currently impossible to express the latter method for arbitrary dimensionalities, but in practice it's not a big issue: most LinearSlow() arrays have a fixed dimension.

This is achieved through dispatch on an internal _getindex method, which recomputes the indices such that it can call the canonical getindex method that the user must define. If the user has not defined their canonical method, it will fall back to an error method in _getindex. I use similar scheme for unsafe_getindex, with the exception that we can fallback to the safe version if the subtype hasn't defined the canonical unsafe method. This enables fast vector indexing by checking bounds of the index vectors instead of on each element. And once @inbounds is extensible, AbstractArrays will be able to support it by default.

The difficulty with all this redirection is that an extra function call can wreck indexing performance, and it can be hard to avoid. I've had particular difficulty getting good performance with CartesianIndexes, and I still lag in performance there by 20x for big arrays. I think call site inline annotations would be a magic bullet, but there may be other tricks we can use, too. I've not looked into this very carefully yet, though. (Fixed with a more sensible inlining strategy)

TL/DR:

In my cursory performance tests hacked onto Tim's indexing perf suite from his reshape work (more tests are needed), I'm close to matching or outperforming master with Array with only these definitions:

julia> methods(getindex, (Array, Any...))
3-element Array{Any,1}:
 getindex(A::Array{T,N},i::Int) at array.jl:304
 getindex{T<:Real}(A::Array{T,N},I::Range{T<:Real}) at array.jl:347 # Needed for bootstrap
 getindex(A::AbstractArray{T,N},I...) at abstractarray.jl:492

(Of course, in places where we're not quite able to close the gap we can always reinstate the specialized methods. This is just a very useful stress-test of both functionality and performance.)

cc: @timholy

[mbauman]

(I was expecting the tests to fail, but I'm amazed at the variety of ways they did so. At a minimum, this needs #9607 (subsumed), #10337, and #10505).

[Jutho]

why is this size check necessary? Haven't you just created dest of that size?

[mbauman]

Not really - similar has. It could give us the wrong answer, and the check here should be very cheap compared to the allocations and many assignments (although I may need to hide it in a function; I haven't profiled yet). In general, my approach has been to not trust the output of similar. Master effectively does the same thing for arrays (although similar and the check are split across method boundaries): multidimensional.jl:219

[timholy]

I'm super-excited about this, and will try to review soon. I just have pretty limited time slots right now.

[johnmyleswhite]

This would be so awesome if it works out.

[timholy]

Do you need an @inline here to ensure good performance for higher dimensions?

[mbauman]

Yes, for N>5. See #10337 (comment).

[timholy]

This is awesome. In broad terms, I agree this is how indexing should work. You also have a lot of clever insights into how to implement this---it was fun to read your code.

You're basically building this on top of #8227? Perhaps we should just merge that and deal with the consequences later.

Unfortunately, it is currently impossible to express the latter method for arbitrary dimensionalities, but in practice it's not a big issue: most LinearSlow() arrays have a fixed dimension.

I'm not sure I agree with this. ReshapedArrays.jl, Interpolations.jl, and your own AxisArrays.jl (if you pass in a LinearSlow() parent) are good counter-examples. I fear we need to introduce getindexN.

Finally, scalar indexing with non-reals is a genuine possibility, see the remarks re DualNumbers.jl in #9586.

[mbauman]

Thanks!

You're basically building this on top of #8227?

Yes, in spirit. I just needed some way to express getindex without a bounds check; BitArrays and a few other data structures already use this idiom, too. I really want user-extensible @inbounds, too, but after working on this I'm not sure having two independent methods is the best solution. It adds a lot of indirection here, and it'd be nice if it could work for things like sub2ind.

I think your alternative proposal in #7799 (comment) might be more attractive, but that just punts the complexity up the chain to a system I don't know and only a few could implement. The method table would keep track of the methods compiled with and without bounds checks. When compiling without bounds checks, the compiler would simply elide everything within @boundscheck. I don't think @inbounds should propagate, so every performant getindex method would be written getindex(...) = (@boundscheck …; @inbounds …).

[mbauman]

The method table would keep track of the methods compiled with and without bounds checks.

Hunh, that sounds a lot like multiple dispatch. getindex(::Type{BoundsCheckOn}, x, I...)? That would actually solve the issues in #8227 with rewriting function names.

[simonster]

Using multiple dispatch has some elegance, but it has problems of its own.

The first problem is that it makes defining getindex uglier for cases where performance doesn't matter: if getindex(x, I...) is now getindex(::Type{BoundsCheckOn}, x, I...), then it seems that all existing definitions need to be changed. We could have a deprecation that gets thrown when an old getindex method is called, similar to @vtjnash's deprecation for convert with pointer types, but arguably users shouldn't need to think about bounds checks unless they need to turn them off for performance purposes.

The second problem is that, at least as far as I can tell, there's still a lowering issue. It's not possible to create an @inbounds macro that changes the first argument to getindex without modifying the frontend (or reproducing the lowering code in Julia), since ref is not lowered to getindex in the AST that the macro gets. Using multiple dispatch would solve @JeffBezanson's scoping concerns, but it seems we'd still need to add a new Expr or an argument to refto specify whether bounds checks are on or off.

[timholy]

Following up from #10507 (comment). @JeffBezanson, your input here would be helpful: how can we restrict dispatch to calling a particular method only for a specific number of varargs? I see a couple of options:

• Introduce getindexN as an indexing function that only gets called when length(indexes) == ndims(A)
• WIP: ReshapedArrays #10507 (comment)
• Introduce getindex{T,N}(A::AbstractArray{T,N}, I::NTuple{N}) (forces a tuple creation unless inlined)
• A new notation, getindex{T,N}(A::AbstractArray{T,N}, indexes...N).

To me the last seems most attractive, followed by the first or third.

[timholy]

One more option worth considering: rather than a tuple, use a CartesianIndex{N}. The main catch: in order to support things like InterpolationArray, the element type should be of any Number. This is a little awkward, given that the main mission for CartesianIndex is as a glorified counter.

[timholy]

Barring further suggestions, I'm going to (time permitting) start playing around with adding a new field to a jl_methlist_t: an index into the tvars field for the parameter specifying the length of the varargs list.
With some parser changes, that should allow us to support the syntax

getindex{T,N}(A::AbstractArray{T,N}, indexes...N)

[mbauman]

That would be awesome!