Hi, Suppose I have a multidimensional array: tmp <- array(1:8, c(2,2,2)) is there a function out there that, given a one-dimensional array index, will return the separate indices for each array dimension? for instance, tmp[8] is equivalent to tmp[2,2,2]. I'd like to derive the vector (2,2,2) from the index 8. thanks, Brad Buchsbaum
I believe that the function below satisfies the request:
rev.index <-
function (x, dims)
{
ld <- length(dims)
cumdim <- cumprod(dims)
ans <- array(NA, c(length(x), ld))
ans[, ld] <- (x-1) %/% cumdim[ld - 1] + 1
ans[, ld-1] <- (x-1) %% cumdim[ld - 1] + 1
if(ld > 2) {
for(i in (ld-1):2) {
y <- ans[, i]
ans[, i] <- (y-1) %/% cumdim[i - 1] + 1
ans[, i-1] <- (y-1) %% cumdim[i - 1] + 1
}
}
ans
}
Patrick Burns
Burns Statistics
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")
Buchsbaum, Brad (NIH/NIMH) wrote:
>Hi,
>
>Suppose I have a multidimensional array:
>
>tmp <- array(1:8, c(2,2,2))
>
>is there a function out there that, given a one-dimensional array index,
>will
>return the separate indices for each array dimension?
>
>for instance, tmp[8] is equivalent to tmp[2,2,2]. I'd like to derive the
>vector (2,2,2)
>from the index 8.
>
>thanks,
>
>Brad Buchsbaum
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
>
>
>
Try this... Perhaps this will help you do what you want. The key idea is to use which() with option arr.ind = TRUE. arr <- array(rnorm(27),c(3,3,3)) dimarr <- dim(arr) tmparr <- array(1:prod(dimarr),dimarr) sapply(c(3),function(x,tmparr) which(tmparr==x,T),tmparr=tmparr) sapply(c(3,17,13,5),function(x,tmparr) which(tmparr==x,T),tmparr=tmparr) Jerome On July 30, 2003 01:42 pm, Buchsbaum, Brad (NIH/NIMH) wrote:> Hi, > > Suppose I have a multidimensional array: > > tmp <- array(1:8, c(2,2,2)) > > is there a function out there that, given a one-dimensional array index, > will > return the separate indices for each array dimension? > > for instance, tmp[8] is equivalent to tmp[2,2,2]. I'd like to derive the > vector (2,2,2) > from the index 8. > > thanks, > > Brad Buchsbaum > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help
"Buchsbaum, Brad (NIH/NIMH)" <BuchsbaB at intra.nimh.nih.gov> writes:> Hi, > > Suppose I have a multidimensional array: > > tmp <- array(1:8, c(2,2,2)) > > is there a function out there that, given a one-dimensional array index, > will > return the separate indices for each array dimension? > > for instance, tmp[8] is equivalent to tmp[2,2,2]. I'd like to derive the > vector (2,2,2) > from the index 8.One not very efficient way is this:> foo <- seq(along=tmp) > dim(foo)<-dim(tmp) > which(foo==8,arr.ind=TRUE)dim1 dim2 dim3 [1,] 2 2 2 For more efficient solutions peek inside the code of which() and see how it gets from ordinary indices to array indices using %% and %/% operations. (Not very easy to grasp, I know). -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
Jerome Asselin <jerome at hivnet.ubc.ca> suggests this:
arr <- array(rnorm(27),c(3,3,3))
dimarr <- dim(arr)
tmparr <- array(1:prod(dimarr),dimarr)
sapply(c(3),function(x,tmparr) which(tmparr==x,T),tmparr=tmparr)
sapply(c(3,17,13,5),function(x,tmparr) which(tmparr==x,T),tmparr=tmparr)
Of course, in R we can simplify the last two lines to
sapply(<<argument goes here>>, function(x) which(tmparr==x,T))
However, wearing my "computer scientist" hat, I have to wonder about
costs.
This is basically the equivalent of the APL "decode" operator.
Let's define
index.decode <- function (index, array) {
dimarr <- dim(arr)
tmparr <- array(1:prod(dimarr), dimarr)
sapply(index, function(x) which(tmparr == x, T))
}
The result is a matrix with C=length(index) columns
and R=length(dim(array)) rows. This has to take time O(R*C), because
the result occupies O(R*C) space and all of it has to be defined.
Now it is possible to implement index.decode so that it does take O(R*C)
time. Here's an outline, which I shan't bother to finish. (I'd do
ndims==4 and the general case if I were going to finish it. I'd also
have a drop= argument to handle the case where length(index)==1.)
index.decode <- function (index, array) {
jndex <- index - 1
dimarr <- dim(arr)
ndims <- length(dimarr)
if (ndims == 1) {
rbind(index)
} else
if (ndims == 2) {
rbind(jndex %% dimarr[1] + 1, jndex %/% dimarr[1] + 1)
} else
if (ndims == 3) {
rbind(jndex %% dimarr[1] + 1,
(jndex %/% dimarr[1]) %% dimarr[2] + 1,
jndex %/% (dimarr[1]*dimarr[2]) + 1)
} else {
stop("length(dims(array)) > 3 not yet implemented")
}
}
This is clearly O(R*C). What about the
sapply(index, function(x) which(tmparr==x, T))
approach?
tmparr is of size prod(dimarr); call that P. The expression tmparr==x
has to examine each element of tmparr, so that's O(P). This is done
for each element of index (C times), so the total is O(P*C).
Consider
mega <- array(1:1000000, c(100,100,100))
inxs <- as.integer(runif(10000, min=1, max=1000000))
Here C = length(inxs) = 10000, R = length(dim(mega)) = 3,
P = prod(dim(mega)) = 1000000. O(R*C) is looking *really* good
compared with O(P*C).
> system.time(index.decode(inxs, mega))
[1] 0.03 0.00 0.03 0.00 0.00> system.time(slow.decode(inxs, mega))
[1] 3.51 0.79 4.33 0.00 0.00
Mind you, on a 500MHz UltraSPARC, I had to use big arrays to get any
measurable time at all...
Dear Philippe, Perhaps you could try a different graphics device (maybe postscript). On my machine, the time differences were all 1 second rather than the 3 you reported. If 300s is really too long for you, you could get a new computer or run your script on a faster one. Andrew C. Ward CAPE Centre Department of Chemical Engineering The University of Queensland Brisbane Qld 4072 Australia andreww at cheque.uq.edu.au Quoting Philippe Hup? <Philippe.Hupe at curie.fr>:> hello, > > I use R1.7.1 under winXP and I am running the following > script example : > > > for (i in 1:10) > { > x <- rnorm(100) > png( paste("D:/essai",i,".png",sep="")) > plot(x) > t1 <- Sys.time() > dev.off() > t2 <- Sys.time() > print(t2-t1) > > } > > at each step, it takes about 3 seconds to shut down the > graphic device. > I want to generate about one hundred of image and of > course it takes too > much time. Is there any trick ? > > Philippe > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help >
Andrew C. Ward wrote:> Dear Philippe, > > Perhaps you could try a different graphics device (maybe > postscript). On my machine, the time differences were all 1 > second rather than the 3 you reported. If 300s is really > too long for you, you could get a new computer or run your > script on a faster one.Let me point out that the "shut down" of a device means that the final calculations and writing of the file is done at that time. So we have to expect it consumes some time. Uwe Ligges> Andrew C. Ward > > CAPE Centre > Department of Chemical Engineering > The University of Queensland > Brisbane Qld 4072 Australia > andreww at cheque.uq.edu.au > > > Quoting Philippe Hup? <Philippe.Hupe at curie.fr>: > > >>hello, >> >>I use R1.7.1 under winXP and I am running the following >>script example : >> >> >>for (i in 1:10) >>{ >> x <- rnorm(100) >> png( paste("D:/essai",i,".png",sep="")) >> plot(x) >> t1 <- Sys.time() >> dev.off() >> t2 <- Sys.time() >> print(t2-t1) >> >>} >> >>at each step, it takes about 3 seconds to shut down the >>graphic device. >>I want to generate about one hundred of image and of >>course it takes too >>much time. Is there any trick ? >> >>Philippe >> >>______________________________________________ >>R-help at stat.math.ethz.ch mailing list >>https://www.stat.math.ethz.ch/mailman/listinfo/r-help >> > > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help