Inconsistency in interpreting a range of CartesianIndexes as a CartesianIndices

[jishnub]

I'm not sure if this is the intended behavior, but the a:b:c method involving CartesianIndexes that results in a CartesianIndices seems inconsistent with conventional usage of range. This is because CartesianIndices behaves as a product of ranges. To illustrate the behavior:

julia> a = CartesianIndex(1,1);

julia> C = range(a, step = a, stop = 2a)
2×2 CartesianIndices{2, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)

julia> first(C) + step(C) == C[2]
false

While a:b:c is a convenient constructor for CartesianIndices, it is inconsistent with the idea of ranges. Should this simply return a StepRangeLen instead? What's worse is that the display of such a StepRangeLen is corrupted by this method

julia> StepRangeLen(a, a, 2)
CartesianIndex(1, 1):CartesianIndex(1, 1):CartesianIndex(2, 2)

julia> StepRangeLen(a, a, 2) == CartesianIndex(1, 1):CartesianIndex(1, 1):CartesianIndex(2, 2)
false

[Seelengrab]

This is because CartesianIndices behaves as a product of ranges.

I think this is intended - a single step is taken as a step in any of the given cartesian directions, resulting in the product of the involved possibilities.

---

What version are you on? StepRangeLen(a, a, 2) fails for me:

julia> versioninfo()                                                                              
Julia Version 1.7.0-beta3                                                                         
Commit e76c9dad42 (2021-07-07 08:12 UTC)                                                          
Platform Info:                                                                                    
  OS: Linux (x86_64-linux-gnu)                                                                    
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz                                                   
  WORD_SIZE: 64                                                                                   
  LIBM: libopenlibm                                                                               
  LLVM: libLLVM-12.0.0 (ORCJIT, skylake)                                                          
Environment:                                                                                      
  JULIA_PKG_SERVER =                                                                              
  JULIA_NUM_THREADS = 4                                                                           
                                                                                                  
julia> StepRangeLen(a, a, 2)                                                                      
Error showing value of type StepRangeLen{CartesianIndex{2}, CartesianIndex{2}, CartesianIndex{2}}:
ERROR: MethodError: no method matching CartesianIndex{2}(::CartesianIndex{2})                     

[johnnychen94]

This inconsistency is due to the fact that you're trying to project the N-dimensional space into 1-dimensional space and thus lose the step information. If you stick to the N-dimensional range concept, which might be a little bit weird at first sight, it's consistent:

# N-dimensional range
C = range(CartesianIndex(1, 1), step=CartesianIndex(2, 2), stop=CartesianIndex(10, 10))
R = CartesianIndices(C)
all(R .- one(eltype(R))) do i
    C[first(R)] + CartesianIndex(i.I .* step(C).I) == C[first(R)+i]
end # true

# 1-dimensional range
C = range(1, step=2, stop=10)
R = LinearIndices(C)
all(R .- one(eltype(R))) do i
    C[first(R)] + i * step(C) == C[first(R)+i]
end # true

[jishnub]

This is because CartesianIndices behaves as a product of ranges.

I think this is intended - a single step is taken as a step in any of the given cartesian directions, resulting in the product of the involved possibilities.

What version are you on? StepRangeLen(a, a, 2) fails for me:

julia> versioninfo()                                                                              
Julia Version 1.7.0-beta3                                                                         
Commit e76c9dad42 (2021-07-07 08:12 UTC)                                                          
Platform Info:                                                                                    
  OS: Linux (x86_64-linux-gnu)                                                                    
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz                                                   
  WORD_SIZE: 64                                                                                   
  LIBM: libopenlibm                                                                               
  LLVM: libLLVM-12.0.0 (ORCJIT, skylake)                                                          
Environment:                                                                                      
  JULIA_PKG_SERVER =                                                                              
  JULIA_NUM_THREADS = 4                                                                           
                                                                                                  
julia> StepRangeLen(a, a, 2)                                                                      
Error showing value of type StepRangeLen{CartesianIndex{2}, CartesianIndex{2}, CartesianIndex{2}}:
ERROR: MethodError: no method matching CartesianIndex{2}(::CartesianIndex{2})                     

Ah I'm sorry, I seem to have added the missing CartesianIndex constructor here. StepRangeLen does indeed fail without that. IMO this is simply a missing method here:

Base.CartesianIndex{N}(c::CartesianIndex{N}) where {N} = c

With this added,

julia> StepRangeLen(a, a, 2)
CartesianIndex(1, 1):CartesianIndex(1, 1):CartesianIndex(2, 2)

[jishnub]

This inconsistency is due to the fact that you're trying to project the N-dimensional space into 1-dimensional space and thus lose the step information. If you stick to the N-dimensional range concept, which might be a little bit weird at first sight, it's consistent:

# N-dimensional range
C = range(CartesianIndex(1, 1), step=CartesianIndex(2, 2), stop=CartesianIndex(10, 10))
R = CartesianIndices(C)
all(R .- one(eltype(R))) do i
    C[first(R)] + CartesianIndex(i.I .* step(C).I) == C[first(R)+i]
end # true

# 1-dimensional range
C = range(1, step=2, stop=10)
R = LinearIndices(C)
all(R .- one(eltype(R))) do i
    C[first(R)] + i * step(C) == C[first(R)+i]
end # true

Yes this does indeed make sense as an N-dimensional iterator, however perhaps there should be a different name given to this? For example, this results in

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])
2×2 Matrix{Int64}:
 1  2
 3  4

Ideally views should not lead to different behavior from slices.

[johnnychen94]

This is another issue, though:

julia> IndexStyle(a[1:end, 1:end])
IndexLinear()

julia> IndexStyle(@view a[1:end, 1:end])
IndexCartesian()

When clarification is needed, I would try to avoid using eachindex and instead use either CartesianIndices or LinearIndices here.

I agree that having a different symbol to : to construct CartesianIndices would make things clearer, but I don't know how generic this symbol can be used for other cases. CartesianIndex might just be one very special case that we ran into this confusion (all because of the fact that CartesianIndex has an order and associated distance.)

[jishnub]

Indeed the IndexStyle differs, but views and slices are only interchangeable if all functions may accept both equivalently. The returned value should not depend on the IndexStyle of the arguments, as views and slices would often have different IndexStyles.

Also, isn't eachindex recommended here, as it will return the optimal iterator? Choosing either LinearIndices or CartesianIndices appears to be sub-optimal in general.

[Seelengrab]

eachindex is correct here, but first(I):last(I) is not - there may not be a last. The returned iterator is (as far as I can tell) not required to be finite. A somewhat contrived example: consider an infinite array, where elements outside a certain index are not stored explicitly. All indices exist, but there is never a last. Cartesian/N-dimensional indexing still makes sense though and if you consider CartesianIndex(0,0) as the origin, even linear indexing may make sense.

If you really have to filter out certain indices based on the index alone, you have to choose whether to filter based on cartesian or linear indices. If you want to be on the save side and know you have an array-like, cartesianindices is a better choice, because it's cheaper to convert a cartesian index to a linear one than it is the other way around (which requires division).

---

The index style for views may be specialized based on what indices are passed in - the information is already part of the type of the view, so I think this could be a good addition. It currently falls back to a generic IndexStyle for all SubArrays. This is however a different issue to how : for cartesian indices works.

[jishnub]

I've created a discourse post to continue the discussion. My view now is that the current behavior of the range a:b is consistent, however that of a:b:c is not. Used in indexing, linear and cartesian indexing lead to different results altogether. Perhaps I should clarify that the issue now is only with a:b:c where the arguments are CartesianIndexes, as once we specify a start, stop and a step, we are definitely talking about a path in the vector space (ie a vector, as we have a magnitude and direction fully specified), and not an iterator product of ranges. In my view this should produce a StepRangeLen to avoid issues like

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

where slices and views behave very differently. This should really be treated as a bug.

Perhaps we may continue further discussion on discourse to avoid cluttering this issue?

[johnnychen94]

My whole point here is that if you choose to look from a 1D perspective, then you're losing step information and you get unwanted things. LinearIndices and CartesianIndices live in intrinsically different worlds, it happens to work well the same when you're doing getindex/setindex!, but for all the other cases, they are quite different.

When you really need to differentiate between an N-dimensional array and a linear array, use CartesianIndices explicitly. eachindex is used to choose the fastest loop index style, and to loop over the array, you don't care about whether it's an N-dimensional array or a linear array.

The g example you provide is an unfair example, even if first(I):2step(I):last(I) produces a linear 1D range as you proposed in #41960, it won't be consistent either; a linear range from CartesianIndex(1, 1):CartesianIndex(2,2):CartesianIndex(4, 4) produces 3 elements [CartesianIndex(1,1), CartesianIndex(2,2), CartesianIndex(4,4)], while that from 1:2:16 produces 8 elements.

I've created a discourse post to continue the discussion.

seems somewhat arbitrary, especially as multiplication is not defined for CartesianIndex es.

one(::CartesianIndex) is deprecated already since Julia 1.5 and to remove it we need to wait for Julia 2.0

julia> one(CartesianIndex{2})
┌ Warning: `one(I::Type{CartesianIndex{N}}) where N` is deprecated, use `oneunit(I)` instead.
│   caller = top-level scope at REPL[1]:1
└ @ Core REPL[1]:1
CartesianIndex(1, 1)

[jishnub]

You're right, the linear range of CartesianIndexes does not resolve it either, I was mistaken there. This would probably require some more thinking. Perhaps a 1D view of CartesianIndices is the right solution here, instead of a 1D range of CartesianIndexes?

I'll close this issue for the time being and reopen if I can think of a better approach.