LinearIndices for non-1-based indices

General Usage
,
1

I’m not sure if I understand LinearIndices and CartesianIndices correctly, but I think they provide a mapping between those index types, and are inverses of each other?

This seems to be the case for normal 1-based indices, but fails for non-1-based indices:

Julia-1.1.1> CI = CartesianIndices((0:5, 0:3))
6×4 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 CartesianIndex(0, 0)  CartesianIndex(0, 1)  CartesianIndex(0, 2)  CartesianIndex(0, 3)
 CartesianIndex(1, 0)  CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 0)  CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 0)  CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)
 CartesianIndex(4, 0)  CartesianIndex(4, 1)  CartesianIndex(4, 2)  CartesianIndex(4, 3)
 CartesianIndex(5, 0)  CartesianIndex(5, 1)  CartesianIndex(5, 2)  CartesianIndex(5, 3)

Julia-1.1.1> LI = LinearIndices(CI)
6×4 LinearIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 1   7  13  19
 2   8  14  20
 3   9  15  21
 4  10  16  22
 5  11  17  23
 6  12  18  24

Julia-1.1.1> CI[1]
CartesianIndex(0, 0)

Julia-1.1.1> LI[0, 0]
ERROR: BoundsError: attempt to access 6×4 LinearIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}} at index [0, 0]
Stacktrace:
 [1] throw_boundserror(::LinearIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}, ::Tuple{Int64,Int64}) at .\abstractarray.jl:484
 [2] checkbounds at .\abstractarray.jl:449 [inlined]
 [3] _getindex at .\abstractarray.jl:949 [inlined]
 [4] getindex(::LinearIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}, ::Int64, ::Int64) at .\abstractarray.jl:927
 [5] top-level scope at none:0

Is this the intended behavior?

In an effort to understand this better, I went down a rabbit hole into the labyrinth of array indexing and nearly didn’t come out.

However, specializing these functions for LinearIndices seems to give the behavior I expected:

function Base.checkbounds(::Type{Bool}, A::LinearIndices, I...)
    Base.@_inline_meta
    Base.checkbounds_indices(Bool, A.indices, I)
end

function Base._sub2ind(A::LinearIndices, I...)
    Base.@_inline_meta
    Base._sub2ind(A.indices, I...)
end

Julia-1.1.1> CI[1]
CartesianIndex(0, 0)

Julia-1.1.1> LI[0, 0]
1

Julia-1.1.1> CI[20]
CartesianIndex(1, 3)

Julia-1.1.1> LI[1, 3]
20

Can anyone comment on this?

1 Like
2

You’ve got almost all of it right, you’re just missing one tiny (but important) fact: as you’ve constructed it, while CI encodes indices that span 0:5 and 0:3, it is accessed by indices that range from 1:6 and 1:4. Whether you realized it or not, this is actually what you asked for. The range 0:5 has its own indices, which are 1:6:

julia> r = 0:5
0:5

julia> r[0]
ERROR: BoundsError: attempt to access 6-element UnitRange{Int64} at index [0]
Stacktrace:
 [1] throw_boundserror(::UnitRange{Int64}, ::Int64) at ./abstractarray.jl:538
 [2] getindex(::UnitRange{Int64}, ::Int64) at ./range.jl:617
 [3] top-level scope at REPL[19]:1

julia> r[1]
0

julia> r[6]
5

All it did was give you back an array using the indices of the component parts you supplied.

Use Base.IdentityUnitRange when you want it accessed by the same indices it encodes:

julia> CI = CartesianIndices(Base.IdentityUnitRange.((0:5, 0:3)))
6×4 CartesianIndices{2,Tuple{Base.IdentityUnitRange{UnitRange{Int64}},Base.IdentityUnitRange{UnitRange{Int64}}}} with indices 0:5×0:3:
 CartesianIndex(0, 0)  CartesianIndex(0, 1)  CartesianIndex(0, 2)  CartesianIndex(0, 3)
 CartesianIndex(1, 0)  CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 0)  CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 0)  CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)
 CartesianIndex(4, 0)  CartesianIndex(4, 1)  CartesianIndex(4, 2)  CartesianIndex(4, 3)
 CartesianIndex(5, 0)  CartesianIndex(5, 1)  CartesianIndex(5, 2)  CartesianIndex(5, 3)

Note in particular the difference in the summary line: mine adds with indices 0:5×0:3 to the end, whereas yours does not (which means it’s 1-indexed).

With this one change, it should work as you were expecting.

I recognize this may be initially confusing. However, in my opinion there’s a good reason for all this power (and consequent confusion). Our indexing is basically governed by a single key axiom, which can be summarized with the following example: if a is an AbstractVector, inds is a vector-of-indices, b = a[inds], and j is another index (integer or vector), then b[j] = a[inds[j]]. Written this way, one realizes that anything else would be nuts. But the implication is that the indices of inds matter—so learn to pay attention to the indices of your indices!

As a developer, to me it’s actually quite surprising how far this single axiom goes towards answering what might otherwise seem like a huge pile of ambiguity about how one should design the API when indexing with arbitrary indices. It became the bedrock that repeatedly saved me from more stupid mistakes than you’d care to know :smile:.

11 Likes
3

Wow, this was super confusing, but it makes sense to me now.
Thank you for taking the time to provide such a great explanation!

1 Like
4

It seems IdentityUnitRange was introduced in Julia v1.1, (PR #28941). Is this correct?

Is there an equivalent way to get this functionality in the LTS version (v1.0.4)?
Could this be backported?
Is there a simple subset of that PR that I could use to patch v1.0.4, or is it more entwined?

5

That PR describes its predecessor: Base.Slice. Both still exist but now Base.Slice should be used solely for internal operations in SubArrays.

The easiest approach that’s agnostic to Julia version is to use OffsetArrays.IdentityUnitRange everywhere, see https://github.com/JuliaArrays/OffsetArrays.jl/blob/950bb889753bfbdcd5313b6f1eb814b28f1a26c2/src/OffsetArrays.jl#L6-L10.

2 Likes
6

Thank you!
This is very helpful.