This is related to the issue https://github.com/JuliaLang/julia/issues/41959, but I didn’t want to derail the issue so thought I would continue the conversation on discourse.
Fundamentally, CartesianIndices represent an iterator product of unitranges or stepranges. The inconsistency arises as ranges are conventionally 1D, however CartesianIndices are also interpreted as a multidiemnsional range. This conflation of an iterator product with a range leads to several unexpected – and perhaps undesired – outcomes.
To start with:
julia> C = CartesianIndices((1:2, 1:3))
2×3 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
CartesianIndex(1, 1) CartesianIndex(1, 2) CartesianIndex(1, 3)
CartesianIndex(2, 1) CartesianIndex(2, 2) CartesianIndex(2, 3)
julia> step(C)
CartesianIndex(1, 1)
The step here represents an iterator product of the step of each range, and not the step of the CartesianIndices object. This is inconsistent with conventional julia usage, where one would use step.(C) to obtain the step along each range. The notion of step(C) doesn’t make sense, as an iterator product does not have a fixed step size. In particular this doesn’t correspond to the step along any axis, but along a diagonal. To see that this may not make much sense, we may see
julia> r = CartesianIndex(1,1):CartesianIndex(1,5);
julia> step(r)
CartesianIndex(1, 1)
julia> r[1] + step(r) == r[2]
false
where (1,1) is clearly not a step along the line that joins the points.
Secondly, this conflation of 1D and nD iterators leads to issues in indexing, in particular the result of a function becomes dependent on the IndexStyle of arguments.
julia> h(a, I = eachindex(a)) = a[first(I):2step(I):last(I)]
h (generic function with 2 methods)
julia> h(a[1:end, 1:end])
2-element Vector{Int64}:
1
2
julia> h(@view a[1:end, 1:end])
1×1 Matrix{Int64}:
1
This breaks the fundamental idea that slices and views are equivalent, and should be treated as a bug.
The question therefore is – what should a range of CartesianIndexes return? My view is that if a step size is specified, the result should unambiguously be a StepRangeLen. A CartesianIndex may be thought of as a point in space, so choosing a start, stop and a step clearly indicates a directed path through space. Choosing a start, step and length also indicates a path through space, as does choosing a stop, step and length.
What should be the result of a:b for CartesianIndex arguments? I’m not totally sure about this, as both arguments seem convincing. On one hand, a:b for numbers produces a unitrange, however there is no equivalent to 1 as a step size in Cartesian space. A vector needs a direction as well as a magnitude. Currently
julia> oneunit(CartesianIndex{2})
CartesianIndex(1, 1)
julia> one(CartesianIndex{2})
CartesianIndex(1, 1)
seems somewhat arbitrary, especially as multiplication is not defined for CartesianIndexes.
julia> b = CartesianIndex(2,2);
julia> b * one(b)
ERROR: MethodError: no method matching *(::CartesianIndex{2}, ::CartesianIndex{2})
This is certainly not a unique choice as a step either, as the direction is arbitrarily chosen along a diagonal for square matrices. On the other hand, if we accept that there’s no notion of one in a vector space, the method a:b seems ideally placed to return a CartesianIndices that spans the entire volume between the two points, so the present behavior actually might make sense here.
julia> CartesianIndex(1,1):CartesianIndex(2,5) # present behavior
2×5 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
CartesianIndex(1, 1) CartesianIndex(1, 2) CartesianIndex(1, 3) CartesianIndex(1, 4) CartesianIndex(1, 5)
CartesianIndex(2, 1) CartesianIndex(2, 2) CartesianIndex(2, 3) CartesianIndex(2, 4) CartesianIndex(2, 5)
This might mean that the following indexing behavior is consistent with now eachindex behaves when used in indexing:
julia> g(a, I = eachindex(a)) = a[first(I):last(I)]
g (generic function with 2 methods)
julia> g(a[1:end, 1:end])
4-element Vector{Int64}:
1
3
2
4
julia> g(@view a[1:end, 1:end]) # view
2×2 Matrix{Int64}:
1 2
3 4
julia> h(a, I = eachindex(a)) = a[I]
h (generic function with 2 methods)
julia> h(a[1:end, 1:end])
4-element Vector{Int64}:
1
3
2
4
julia> h(@view a[1:end, 1:end])
2×2 Matrix{Int64}:
1 2
3 4
While not ideal, the results using views and slices at least have the same elements. One may reshape the results to 1D to get around the changing number of dimensions. The situation with a step specified is, however, worse, as the elements themselves may be different. I’m not sure if there is a way to resolve this without actually requiring that the result be a 1D range of CartesianIndexes and not a CartesianIndices.