Confusion On Cartesian And Linear Index

I know cartesian index is to be of this nature: x[i, j] and linear index is of x[i].

With that said, let me explain my confusion. I have a 3D array:

julia> A = [1 2 3; 4 5 6;;; 7 8 9; 10 12 13]
2×3×2 Array{Int64, 3}:
[:, :, 1] =
 1  2  3
 4  5  6

[:, :, 2] =
  7   8   9
 10  12  13

What I’m confused about is that: I know I can retrieve the element of row 2, column 2 of block 1 in cartesian or linear index using this:

julia> A[2, 2, 1]
5

julia> A[4]
5

But then I see everywhere (in most source code and blog posts) people using:

julia> A[CartesianIndex((2, 2, 1))]
5

And they go on talking about the usefulness of it. I’ve searched through the docs and every blog posts I’ve read on this, but I can’t still see the benefits of the latter over the former. And then there is also LinearIndices and CartesiaIndices . I’ve assuringly thought they’re just constructors for the syntax, but I feel I’m missing something.

So my questions arising from the confusions are:

  1. When should I just use the syntax and not worry over the constructors?
  2. What are the benefits of using the constructors over the syntax? If so, then what are the use cases?

I’m asking those questions not from developing a package or creating a custom array, but just using pure julia. Thank you.

2 Likes

Please read Tim Holy’s blog post on Multidimensional algorithms and iteration, where we can learn that a key use for CartesianIndices and Julia iterators is in writing generic multidimensional algorithms.

A simple example is provided below, but more complex ones can be found in other Discourse posts.

Suppose we wanted to write dimension-agnostic code that sets the values of an arbitrary multidimensional array to: -1 if all the indices are odd, +1 if they are all even, and 0 otherwise.

We could achieve this with CartesianIndices as follows:

A = Array{Int8}(undef, 3, 3, 2)     # can set arbitrarily the number of dimensions

for ci in CartesianIndices(A)
    A[ci] = all(isodd, Tuple(ci)) ? -1 : all(iseven, Tuple(ci)) ? 1 : 0
end
A
3 Likes

Thank you for this. I’m beginning to understand A[i, j] should be used in normal casual cases, but then, when one wants to do advance concepts with multi-dimensional arrays, certainly CartesianIndex or CartesianIndices are needed.

For example, after I ask the questions as about 2 days ago and got no reply, I decided to think for use cases where I would actually need CartesianIndex. So I gave up with a problem to create a function that returns the left-to-right and right-to-left diagonal of a 2D array. It was then it done on me, that in most cases I would actually need this:

"""
    ldiag(M; step::Integer=0)

Returns the left-to-right diagonal of a matrix as a vector.

A `step` value can be specified, and the diagonal
below (if `step < 0`) or above (if `step > 0`) the
principal diagonal is returned.
"""
function ldiag(M; step::Integer=0)
    row, col = size(M)
    r, c = 1, 1
    if step < 0
        r = abs(step - 1)  # start from this row
    elseif step > 0
        c = abs(-step - 1) # start from this column
    end
    M[CartesianIndex.(zip(r:row, c:col))]
end


"""
    rdiag(M; step::Integer=0)

Returns the right-to-left diagonal of a matrix as a vector.

A `step` value can be specified, and the diagonal below
(if `step < 0`) or above (if `step > 0`) the principal
diagonal is returned.
"""
function rdiag(M; step::Integer=0)
    row, col = size(M)
    r, c = 1, 1
    if step < 0
        r = abs(step - 1)  # start from this row
    elseif step > 0
        col = col - step    # start from this column
    end
    M[CartesianIndex.(zip(r:row, col:-1:c))]
end

So yes CartesianIndex and its counterparts are certainly needed for most advanced use cases.

FWIW, such a CartesianIndex broadcast could be written alternatively as:

[M[t...] for t in zip(r:row, col:-1:c)]
1 Like
[M[t...] for t in zip(r:row, col:-1:c)]

That part never occurred to me. As the original one I wrote with:

M[CartesianIndex.(zip(r:row, col:-1:c))]

used more allocations. Indeed for loops are faster in Julia. When possible, operations should be devectorized. Thank you for this.