2

Let's say I have the array

TestArray=array(1:(3*3*4),c(3,3,4))

In the following I will refer to TestArray[i,,], TestArray[,j,] and TestArray[,,k] as the x=i, y=j and z=k subsets, respectively. In this specific example, the indices i and j can go from 1 to 3 and k from 1 to 4.

Now, I want to subset this 3 dimensional array so that I get the x=y subset. The output should be

do.call("cbind",
        list(TestArray[1,1,,drop=FALSE],
             TestArray[2,2,,drop=FALSE],
             TestArray[3,3,,drop=FALSE]
            )
        )

I have (naively) thought that such an operation should be possible by

library(Matrix)
TestArray[as.array(Diagonal(3,TRUE)),]

This works in 2 dimensions

matrix(1:9,3,3)[as.matrix(Diagonal(3,TRUE))]

However, in 3 dimensions it gives an error.

I know that I could produce an index array

IndexArray=outer(diag(1,3,3),c(1,1,1,1),"*")
mode(IndexArray)="logical"

and access the elements by

matrix(TestArray[IndexArray],nrow=4,ncol=3,byrow=TRUE)

But the first method would be much nicer and would need less memory as well. Do you know how I could fix TestArray[as.array(Diagonal(3,TRUE)),] so that it works as desired? Maybe I am just missing some syntactic sugar...

cryo111
  • 4,444
  • 1
  • 15
  • 37
  • would `?slice.index` be helpful? – Ben Bolker Jul 17 '13 at 18:00
  • based on http://stackoverflow.com/questions/16607232/r-array-slice-of-programmatically-determined-dimension/16609409#16609409 , maybe `abind::asub` is what you want instead. – Ben Bolker Jul 17 '13 at 18:02
  • Thank you for your suggestions regarding `slice.index` and `asub`. I did not know these functions yet. However, they do not seem to accomplish the task I have outlined above. – cryo111 Jul 17 '13 at 20:10
  • BTW: I think "irregular subsetting" is more appropriate than just "subsetting" in the title. I have changed that accordingly. – cryo111 Jul 17 '13 at 20:11

2 Answers2

1

I don't know if abind::asub will do what I (you) want. This uses a more efficient form of matrix indexing than you have above, but I still have to coerce the results into the right shape ...

indmat <- cbind(1:3,as.matrix(expand.grid(1:3,1:4)))
matrix(TestArray[indmat],nrow=4,ncol=3,byrow=TRUE)

Slightly more generally:

d <- dim(TestArray)[1]
d2 <- dim(TestArray)[3]
indmat <- cbind(1:d,as.matrix(expand.grid(1:d,1:d2))
matrix(TestArray[indmat],nrow=d2,ncol=d,byrow=TRUE)
cryo111
  • 4,444
  • 1
  • 15
  • 37
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Indeed, your indexing is more efficient. Thank you for your suggestion. As it seems, there is no way to do it as described in my original post. – cryo111 Jul 17 '13 at 20:13
1

In addition to Ben's answer, here a surprisingly simple modification to my original line of code that does the job.

matrix(TestArray[as.matrix(Diagonal(3,TRUE))],ncol=3,nrow=4,byrow=TRUE)

This works because as.matrix(Diagonal(3,TRUE)) gets recycled.

cryo111
  • 4,444
  • 1
  • 15
  • 37