R indexing arrays. How to index 3 dimensional array by using a matrix for the 3rd dimension

Question

I'm having a question about indexing 3 dim arrays.

Say I have a 3 dimensional array

x<- c(1:36)
dim(x) <- c(3,4,3) 

Now I want to extract values out of this array according to a matrix holding the 3rd dimension indices for all [i,j] positions.

y <- c(rep(1,4),rep(2,4),rep(3,4))
dim(y) <- c(3,4)

y
      [,1] [,2] [,3] [,4]
[1,]    1    1    2    3
[2,]    1    2    2    3
[3,]    1    2    3    3

So the result should be giving this:

     [,1] [,2] [,3] [,4]
[1,]    1    4   19   34
[2,]    2   17   20   35
[3,]    3   18   33   36

Is there some elegant way to do this? I know how to use two for loops to go over the array, but this is too slow for my data.

Answer 1

help("[") tells us this:

Matrices and arrays

[...]

A third form of indexing is via a numeric matrix with the one column for each dimension: each row of the index matrix then selects a single element of the array, and the result is a vector.

Thus, we transform your y matrix to a shape that conforms with this.

library(reshape2)
z <- x[as.matrix(melt(y))]
dim(z) <- dim(y)
#     [,1] [,2] [,3] [,4]
#[1,]    1    4   19   34
#[2,]    2   17   20   35
#[3,]    3   18   33   36

Answer 2

I'm looking at this as an opportunity for some code golf. It's definitely possible to do this as a one-liner:

> `dim<-`(x[cbind(c(row(y)), c(col(y)), c(y))], dim(y))
     [,1] [,2] [,3] [,4]
[1,]    1    4   19   34
[2,]    2   17   20   35
[3,]    3   18   33   36

As @Roland's answer shows, matrix/array indexing involves creating an n-column matrix and setting the columns equal to row, column, etc. position of each dimension of an n-dimensional array. We can use the row() and col() functions to extract the row and column positions of each element in y:

> row(y)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3
> col(y)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4

and y itself gives the third-dimension positions. wrapping each of those in c() turns them into a vector, so that they can be cbind-ed together to create an extraction matrix.

Then, there's just some fun use of dim<-() to fit it all on one line.