Make CartesianRange an AbstractArray

[barche]

The use case for this change is easy conversion between linear and cartesian indices, which can be done by linearly indexing into the CartesianRange now.

[timholy]

An interesting and creative idea; I don't think this has ever occurred to me. With this, it seems like one might even be able to get rid ind2sub and sub2ind.

If I try to look into my crystal ball (something I'm not always great at), I think the biggest downside might be the risk of more ambiguities. That doesn't seem likely to be serious enough to block this. So I'm (ever so slightly nervously) 👍 on this.

[barche]

OK, I made the calculation slightly more elegant (imho) by moving everything into the CartesianIndex constructor.
I had forgotten about ind2sub and its cousin, but I have no ambition to replace those. For completeness, replacing ind2sub is easy:

CR = CartesianRange((1:3, 1:5))
CR[4] # CartesianIndex(1, 2)
ind2sub((3,5), 4) # (1,2)

The inverse is more cumbersome, however:

rs = reshape(linearindices(CR), size(CR))
rs[1,2] # 4
sub2ind((3,5), 1, 2) # 4

On the other hand, once you have CR and rs and keep them around, the indexing syntax is nice.

[nalimilan]

It would be nice to be able to get rid of ind2sub and sub2ind in 0.7.

[barche]

If there is consensus that this is a good idea, one approach would be as follows:

1. Replacing ind2sub just needs the current PR change and documentation changes
2. For sub2ind having to call reshape seems a bit cumbersome, so maybe a method taking a CartesianRange would be better, e.g:

linearmapping(iter::CartesianRange) = reshape(linearindices(iter), size(iter))

Then linearmapping (better names are possible, no doubt) would be exported and ind2sub and sub2ind no longer exported. Their implementations must remain however, since they are called internally anyway using the new method.

I'd be willing to extend this PR if this seems feasible for 0.7, and if someone can help with carefully checking that there are no performance regressions.

[mbauman]

This is something that I've considered myself, but I had been putting it in the same "performance trap" bucket as we have done for iterating over dimensions of individual CartesianIndexes. This is significantly different, though, since we have a way of stating the optimal way of indexing into arrays. Heck, if you need lots of ind2sub conversions, you can just reshape the array (or use non-scalar indexing) to take advantage of ReshapedArray's multiplicative inverses. That's cool.

Maybe we should change what indices(A) returns? Instead of a tuple of unit ranges, they could be wrapped in a CartesianRange. Then ind2sub(A, 3) is even more simply indices(A)[3].

The logical conclusion, then, I think would be to have linearindices return a reshaped UnitRange, such that sub2ind(A, (3,4,5)) becomes linearindices(A)[3,4,5].

Even better, this would mean that A == A[indices(A)] == A[linearindices(A)].

Just a brainstorm here. I like this, though.

[barche]

@mbauman I like the elegance of those ideas, but wouldn't this conflict with the custom indices specification, where different indexing types can be chosen along each dimension? Also, I tried changing indices as you suggested, but can't get sysimage to build or even get a successful evaluation using Revise.jl, so I guess this would be a very far-reaching change.

On the other hand, I noticed eachindex is currently undefined for CartesianRange, so maybe that could be overridden to return the remapped array?

[barche]

I have added the eachindex implementation and a constructor that takes an array, and also updated docs and NEWS.md.

I'll look into deprecating ind2sub and sub2ind, but that seems more involved.

[StefanKarpinski]

Triage finds this to be a lovely idea and would love to see it merged :)

[JeffBezanson]

I'm wondering if eachindex and linearindices should really be different here --- I think of eachindex as just being for performance, so we might not want to attach this subtle semantic difference to it.

[barche]

I was just working on the ind2sub and sub2ind deprecation, and I think instead of making a specialization on eachindex I should add a type CartesianToLinear that does the sub2ind conversion. The main advantage over just using reshaped linear indices is that this will also work with non-conventional indices, just like sub2ind does, so the following would be equivalent:

sub2ind((0:3,3:5), 0, 3)
CartesianToLinear(0:3,3:5)[0,3]

ind2sub((0:3,3:5), 1)
CartesianRange(0:3,3:5)[1]

This also has a nice symmetry to it.

[barche]

I've added the commit that deprecates sub2ind and ind2sub. Internally, I replaced the name by prefixing it with _. For now, I just replaced uses in Base with the renamed functions, but it may be more elegant to limit that to the file abstractarray.jl and use CartesianRange and CartesianToLinear elsewhere.

[nalimilan]

It feels a bit strange to me to have a type called CartesianToLinear: this sounds like the name of a conversion function rather than that of a type. I would find it more natural to do something like LinearIndex(dims::NTuple{Int}, x::CartesianIndex)::Int, i.e. similar to sub2ind but with a clearer name and accepting CartesianIndex.

[barche]

The reason it is a type is for symmetry with CartesianRange: this is an array that when indexed with a linear index yields the equivalent cartesian index. Similarly, an object of type CartesianToLinear yields a linear index when indexed with a cartesian index. Note that usually this will be called in a loop, so it's possible to use it like this:

linear = CartesianToLinear(dims)
for i in ...
  for j in ...
    k = linear[i,j]
    ...
  end
end

I agree it's possible to find a better name for the type, however :)

[timholy]

Sorry I lost track of this. I fixed some merge conflicts. This seems ready to merge when CI passes again.

The CI queue seems likely to be stalled out today, but that should not prevent this from being merged for 1.0.

[timholy]

Updated version in #25113