# Inconsistency in interpreting and indexing with ranges of CartesianIndexes if a step is specified while constructing them

This is related to the issue Inconsistency in interpreting a range of CartesianIndexes as a CartesianIndices · Issue #41959 · JuliaLang/julia · GitHub, 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.

``````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 `CartesianIndex`es 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 `CartesianIndex`es.

``````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 `CartesianIndex`es and not a `CartesianIndices`.

One of the points raised in the original issue was that `first(I):last(I)` may not be the right way to index an array, as there may not be a `last`. I feel this is beside the point here. We may limit the discussion to `Array`s that do have a well-defined `lastindex`.