Inconsistent ranges of CartesianIndex

[Seelengrab]

I'd have expected these two to be equivalent, but the second one errors:

julia> const ci = CartesianIndex
CartesianIndex

julia> ci(1,1):ci(3,1) |> collect
3×1 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 1)
 CartesianIndex(3, 1)

julia> ci(1,1):ci(1,0):ci(3,1) |> collect
ERROR: ArgumentError: step cannot be zero
Stacktrace:
 [1] steprange_last(start::Int64, step::Int64, stop::Int64)
   @ Base ./range.jl:325
 [2] StepRange
   @ ./range.jl:312 [inlined]
 [3] StepRange
   @ ./range.jl:367 [inlined]
 [4] _colon
   @ ./range.jl:24 [inlined]
 [5] Colon
   @ ./range.jl:22 [inlined]
 [6] #7
   @ ./multidimensional.jl:318 [inlined]
 [7] map (repeats 2 times)
   @ ./tuple.jl:266 [inlined]
 [8] (::Colon)(I::CartesianIndex{2}, S::CartesianIndex{2}, J::CartesianIndex{2})
   @ Base.IteratorsMD ./multidimensional.jl:317
 [9] top-level scope
   @ REPL[3]:1

Unintuitively, the "correct" way to write this would be ci(1,1):ci(1,x):ci(3,1) for any value x other than 0:

julia> ci(1,1):ci(1,1):ci(3,1) |> collect
3×1 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 1)
 CartesianIndex(3, 1)

julia> ci(1,1):ci(1,10):ci(3,1) |> collect
3×1 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 1)
 CartesianIndex(3, 1)

I think this is a consequence of how the cartesian product is formed right now - maybe dimensions where the step is zero can be just skipped during product expansion, provided both start and stop of the to-be-constructed range have the same coordinate in that point? Since no dimensions are lost, it should still be type stable.

[stev47]

The default step-size for ranges is one, even if start and end agree. In your example the step-size in the second dimension should be one, not zero.
The error boils down to step-ranges not allowing a zero step-size (try 1:0:1), which is a useful convention (think e.g. of length calculation). See also this discussion for further details.

[johnnychen94]

See also the discussions in #41959 and #41960

There are two valid interpretations of CartesianIndex{N} here:

• a point (x¹, x², ..., xᴺ)` in N-dimensional space: pⁱ = p₀ⁱ + xⁱ * Δⁱ where pⁱ are real numbers.
• a point (x,) in 1-dimensional space: p = p₀ + x * Δ where x is real number.

ci(1,1):ci(3,1) does not construct a 1-dimensional range; it goes the first interpretation and thus the output is a 2-D matrix instead of a 1-D vector.

[Seelengrab]

@stev47 Can you please not just close an issue where I already mention that this error is due to the assumption that the cartesian product of each dimension "just works" for the expansion over a nonzero step? The discourse discussion you're linking is also not relevant here - CartesianIndex(1,0) is decidedly not zero(CartesianIndex{N}) for any N.

@johnnychen94 Yes, I've read & participated in both of those discussions. I don't think they apply here. My complaint is not with a:b, but with a:b:c behaving wrongly in the case that b has a step of 0 in one direction, but otherwise does create elements. Specifically I think it would be nice if CartesianIndex(1,1):CartesianIndex(1,0):CartesianIndex(10,1) == CartesianIndex(1,1):CartesianIndex(10,1) or CartesianIndex(1,1):CartesianIndex(0,1):CartesianIndex(1,10) == CartesianIndex(1,1):CartesianIndex(1,10) would both hold.

In short, I think this should only error if the step is 0 in all dimensions (i.e. the step is zero(CartesianIndex{N})). In other words, dimensions in the step that are 0 must be equal in both start & stop, otherwise build the cartesian product over the nonzero dimensions with extent in start & stop. This does not loose dimension information and should be consistent with the current implementation.

[Seelengrab]

Maybe let me expand on the motivation here - creating a cartesian product is not my main goal. I'm mainly interested in getting steps like CartesianIndex(1,0,0,0) or CartesianIndex(0,0,5,0) to work, since they're very useful when wanting to create a line along a given axis in multidimensional arrays. I'm fully aware that this will not and cannot create a line that's off-axis (e.g. with step like CartesianIndex(1,0,3,0)). That is not the goal here.

[johnnychen94]

@Seelengrab I see your point now. We should revive the discussion #15218 first; once 1:x:1 is supported then it's immediately followed that CartesianIndex(1, 1):CartesianIndex(1,0):CartesianIndex(4, 1) is also supported. For the moment let me reopen this issue.

[Seelengrab]

I'm not sure we have to solve the general case mentioned in that issue for this to be possible - simply ignoring the dimensions with 0 as their step during the product would already be consistent with the cartesian-ness of the CartesianIndex range, no? After all, the dimension isn't lost since all points are then required to have the same coordinate in that dimension. It just errors right now because the code assumes that this is the same as a cartesian product over the individual stepranges, which imo is a different case and would only matter when the step is zero(CartesianIndex{N}).

[johnnychen94]

simply ignoring the dimensions with 0 as their step during the product would already be consistent with the cartesian-ness of the CartesianIndex range, no?

Yes, that's a shortcut but I think it's essentially no different to the general 1-d range case. I understand your point that CartesianIndex is more useful a case that we may want to support; maybe you can put up a PR and get more serious feedback from others? It should be quite easy, just needs some effort on docstrings and tests.

# something like this :)
# https://github.com/JuliaLang/julia/blob/ca11c2fc24a64a083c0b06a4abcdfb41c0c4f2cb/base/multidimensional.jl#L317-L318
function Base.:(:)(I::CartesianIndex{N}, S::CartesianIndex{N}, J::CartesianIndex{N}) where N
    R = map(Tuple(I), Tuple(S), Tuple(J)) do i, s, j
        if iszero(s) && i == j
            i:one(i):i
        else
            i:s:j
        end
    end
    CartesianIndices(R)
end

[stev47]

@Seelengrab sorry for being hasty but no harm was done to you.
@johnnychen94 thank you for linking #15218, that discussion is interesting and I can see some argument for allowing zero-steps in general.
I still think that both issues are linked, since one would expect 1:0:1 and CartesianIndex(1):CartesianIndex(0):CartesianIndex(1) to behave similarly.

[Seelengrab]

sorry for being hasty but no harm was done to you.

No worries! Just felt very dismissive, as I've explicitly mentioned that I'm concerned with nonzero steps. Maybe I could have phrased that more clearly in the initial post.

I still think that both issues are linked, since one would expect 1:0:1 and CartesianIndex(1):CartesianIndex(0):CartesianIndex(1) to behave similarly.

That is orthogonal to this issue though, as CartesianIndex(0) truly is zero(CartesianIndex{1}). I'm explicitly talking about cases where at least one dimension isn't zero, making sure the step is never zero of the relevant type :)

[stev47]

I see, then it is indeed more viable.
It would however mean, that they no longer respect the tensor product structure, i.e. for tuples a, b and s:

CartesianIndex(a):CartesianIndex(s):CartesianIndex(b) ==
  CartesianIndex.(Iterators.product(map((a, s, b) -> a:s:b, a, s, b)...))

Not sure how valuable that consistency is.