Broadcasting pointwise indexing operator

[toivoh]

I would like to suggest to add a broadcasting pointwise indexing operator to Julia. This a very expressive form of indexing, and quite handy, as remarked in #2247.

Let us call the new operation pointref(A, inds...) = A.[inds...], where A and the elements of inds are arrays. Then we can compute

X = pointref(A, inds...) = A.[inds...]

as follows:

• Broadcast all elements of inds to the same shape: Extrude each singleton
  axis in each ind array to the common non-singleton size of that
  axis among them. It is an error if the indices cannot be broadcast together.

• The shape of X is the common broadcast shape of the indices.

• Each entry of X is calculated by indexing A with the elements of inds at the corresponding position. Suppose that inds have been broadcast to the same shape. Then
  for each element of the result X[ks...],

  X[ks...] = A[[ind[ks...] for ind in inds]...]
  

Example:

A = [1 3
     2 4]

P = [2 1 1]
Q = [1 2 1
     2 1 2]

X = pointref(A, P, Q) = A.[P, Q]

Since size(P) = (1,3) and size(Q) = (2,3), their broadcast shape is (2,3). The broadcast index arrays are

Pb = [2 1 1
      2 1 1]
Qb = [1 2 1
      2 1 2]

where P has been extruded along the first axis. Then

X = [A[2,1] A[1,2] A[1,1]
     A[2,2] A[1,1] A[1,2]]

  = [2 3 1
     4 1 3]

Relationship to current indexing

Current Julia indexing allows

A[inds...]

where the elements of inds are vectors and/or scalars. This axiswise indexing can be expressed using pointwise indexing as

A[inds...] = A.[[swapdims(ind,1,k) for (k,ind) in enumerate(inds)]...]

where swapdims(A,j,k) swaps axes j and k of an array. We see that pointwise and axiswise indexing coincide when there is only one index, since swapdims(A,1,1) is a noop.

Julia also supports linear indexing

A[ind] = A[:][ind]

Analogously, one could consider

A.[ind] = A[:].[ind]

i.e. A.[ind] looks does a linear-index lookup in A for each element in the index array ind.
This seems to be the same behavior as suggested in #2247 for A[ind].

Wrap-up

• I think that pointwise indexing is useful enough to deserve its own operator.
• It's also different enough from regular indexing that it seems impossible to overload the two in the same operator, except for the single-index case A[B], where the two can be said to coincide. (As described in RFC: Linear indexing with multidimensional index. #2247 (comment), numpy does this also in the multi-index case, which becomes messy.)
• It also happens to be a natural way to express the kind of indexing that a gpu kernel is capable of, into an array stored in gpu memory.

[JeffBezanson]

Very good write-up. Quite compelling.

[ViralBShah]

This is indeed quite compelling. Keeping syntax aside, I don't like the idea of having a whole new operator for this. Will sleep on this for a bit.

[andyferris]

Can't this already be performed with

X[f(inds...)]

for suitable f? I'm trying to figure out why this needs to be syntax with different meaning to getindex.(X, inds...) which currently gets the elements from the nested arrays inside X (e.g. X::Vector{Vector{T}})

[mbauman]

Yes, it's possible with X[CartesianIndex.(inds...)].

[andyferris]

Would a Generator be perfect to do the kinds of things you could achieve with pointwise indexing?

The shape of X would be given by the for dimensions. The generator should create a tuple for multidimensional indexing

diagonals = matrix[i,i for i = 1:min(size(matrix))]

or else do linear indexing

mycopy = matrix[sub2ind(size(matrix), i, j) for i = 1:end, j = 1:end]

(actually my usage of end here is a bit unusual, not sure how that should work...)

[GunnarFarneback]

This proposal is growing on me, thanks to the analogy with dotted function calls and the potential of dot fusion. If I understand it correctly, the original example can in Julia 0.6 be implemented as

getindex.((A,), P, Q)

and it would participate in dot fusion. It does break down a bit for scalar P and Q though.
The proposed syntax A.[P, Q] would read a lot better, especially as part of a bigger expression.

[mbauman]

Perhaps we should also deprecate the use of the getindex syntax within @.?

Personally, I still find the interplay between non-scalar non-dotted indexing and dot-indexin quite confusing, but perhaps it's something that we'll get used to.

[tkf]

Re #27563 (comment):

using a . to enable non-scalar fusion through to particular dimensions of a non-scalar indexing expression. This idea can actually stand completely on its own and isn't simply an extension of this particular & proposal.

You are thinking something like f.(A[g(idx[:, .])]) to do [f(@view A[:, g(@view idx[:, k])]) for k in 1:size(idx, 2)]? (Edit: oops, that was a gibberish) Interesting. I have to think about it.