Not wrapping CartesianIndices incorrectly processed by to_indices

[sethaxen]

Indexing a matrix with Not(::CartesianIndices) results in Not only being applied to the row selector, not the column selector:

julia> x = reshape(1:12, (4, 3))
4×3 reshape(::UnitRange{Int64}, 4, 3) with eltype Int64:
 1  5   9
 2  6  10
 3  7  11
 4  8  12

julia> inds = CartesianIndices(x)[2:3, 1:2]
CartesianIndices((2:3, 1:2))

julia> x[inds]
2×2 Matrix{Int64}:
 2  6
 3  7

julia> x[Not(inds)]  # expected setdiff(x, x[inds])
2×2 Matrix{Int64}:
 1  5
 4  8

julia> x[Not(vec(inds))] == setdiff(x, x[inds])  # fine if we reshape to a vector
true

julia> x[Not(Matrix(inds))] == setdiff(x, x[inds])  # fine if we convert to a Matrix
true

julia> lin_inds = LinearIndices(x)[2:3, 1:2]
2×2 Matrix{Int64}:
 2  6
 3  7

julia> x[Not(lin_inds)] == setdiff(x, x[inds])  # fine if we convert to linear indices
true

This happens because to_indices only applies the Not to the first selector being returned:

julia> Base.to_indices(x, (inds,))
(2:3, 1:2)

julia> Base.to_indices(x, (Not(inds),))  # Not only applied to row selection
(Not(2:3), 1:2)

julia> Base.to_indices(x, (vec(inds),))  # calling vec on the indices makes to_indices a no-op
(CartesianIndex{2}[CartesianIndex(2, 1), CartesianIndex(3, 1), CartesianIndex(2, 2), CartesianIndex(3, 2)],)

julia> Base.to_indices(x, (Matrix(inds),))  # same for converting to matrix
(CartesianIndex{2}[CartesianIndex(2, 1) CartesianIndex(2, 2); CartesianIndex(3, 1) CartesianIndex(3, 2)],)

AFAICT this is because these lines assumesthat Not's skipped indices all are resolved by to_indices to the first entry of the resulting tuple, but for CartesIndices this is (at least somtimes?) false.

InvertedIndices.jl/src/InvertedIndices.jl

Lines 181 to 183 in 1253c44

	new_indices = to_indices(A, inds, (I[1].skip, tail(I)...))
	skips = uniquesort(new_indices[1])
	picks = spanned_indices(inds, skips)[1]

[mbauman]

There's a question here — what should Not(CartesianIndex()) and Not(CartesianIndices(()) do?

Of course as the only index into a zero-dimensional array, they're not selecting the only value. That I can wrap my head around. But what about (1:3)[Not(CartesianIndex()), 1]? The behavior that "falls out" of #47 happens to be:

julia> (1:3)[CartesianIndex(()), 1]
1

julia> (1:3)[Not(CartesianIndex((1,))), 1]
1

julia> (1:3)[Not(CartesianIndex(())), 1]
Int64[]

julia> (1:3)[Not(CartesianIndex((1,))), 1]
2-element Vector{Int64}:
 2
 3

julia> (1:3)[CartesianIndices(()), 1]
0-dimensional Array{Int64, 0}:
1

julia> (1:3)[CartesianIndices((1,)), 1]
1-element Vector{Int64}:
 1

julia> (1:3)[Not(CartesianIndices(())), 1]
Int64[]

julia> (1:3)[Not(CartesianIndices((1,))), 1]
2-element Vector{Int64}:
 2
 3

I guess you're Not selecting no index at all? It's weird, a bit like saying let's not take that tree in the forest I can't hear falling.

[sethaxen]

I didn't know CartesianIndex() was even possible, but it looks like x[CartesianIndex()] is the same as x[] so is allowed for any length one n-dimensional array. In that case I'd say x[Not(CartesianIndex())] should behave the same as x[Not(:)] for these arrays, which already returns an empty array. This seems to be what already happens.

But it looks like one can insert an arbitrary number of CartesianIndex() into an indexing expression, and these are effectively removed:

julia> x = reshape(1:6, 2, 3);

julia> x[2, 3]
6

julia> x[CartesianIndex(), 2, CartesianIndex(), 3, CartesianIndex()]
6

I can't think of a reasonable thing to do if any of those CartesianIndex() are wrapped in a Not, except for them to be a no-op or to raise an error. No-op seems reasonable, since the insertion of a CartesianIndex has no effect on the indexing, wrapping it in a Not shouldn't either.

Similar with CartesianIndices; if its inclusion doesn't change which values are returned, then I think it's reasonable for wrapping it in a Not to have no effect. But these are all edge cases I haven't worked with or thought through extensively.