Dear list members,

I would appreciate your advice on performance improvement of the following.

System used: intel core i7, 2.66GHz, 2 cores, 8GB RAM, MacOS 10.6.8, R64:
2.14.2

Starting with a binary array of n= 9e6 elements in rows versus 100 subjects
in cols, I would like to do resampling statistics using R=2000 reshuffles
on each row.
As result I want to store (a) the original statistic for each element
(row), (b) the R reshuffled statistics per element, for later calculation
of false discovery rate, and (c) the p-value of the original statistic
based on the null distribution observations.
If I did my homework right, for best performance I want to (1) use apply
instead of for loops, (2) parallelize the calculations, in (3) quite large
chunks of data.
The question I have mainly concerns (3) for big.matrix or ff objects.
Please find an example below, in which I used 9e2 instead of 9e6 elements
for convenience. The current code takes 21.9 s for n=9e2, 2638 s for n=9e4.
A rough estimation would be a ~70 hr run for 9e6 elements... Any ideas of
performance improvement are warmly appreciated.

I tried several options:
(1) using ff instead of big.matrix yields similar performance
(2) since foreach(i=1:n) seems a bottleneck, I tried the iblkcol example
from the ' Using the foreach package' documentation, intending to iterate
selected columns of AB, with simiar performance (see 'attempt 2' below).
(3) another possible bottleneck may be the  %dopar% { get.o(AB[,i]) }, in
which AB is distributed over the slaves at each iterations, but I can't
find a way to select sub.matrices for big.matrix or ff more efficiently.

Thanks for you time and efforts, Philip


### CODE  START ###

library(bigmemory)
library(foreach)
library(doMC)
registerDoMC(4)

# array AB containing 9e6 elements for 100 subjects into a char (0/1) array
# now using 9e2 in this code for convenience, some dummy 0/1's used here
n <- 9e2
AB <- big.matrix(nrow=100, ncol=n, type='char', init=0)
foreach(i=1:ncol(AB)) %dopar% { AB[,i] <- round(runif(nrow(AB))) })

# subject 1:50 belong to 'group 1' and subject 51:100 belong to 'group 2'
grp.org <- c(rep(1,50),rep(2,50))

# shuffle labels '1' and '2' R, say 2000, times
# result is grp.rnd holding R rows each containing reshuffle of 100 cols of
group labels 1/2
R <- 2000
grp.rnd <- array(dim=c(R, length(grp.org)))
for (i in 1:R) { grp.rnd[i,] <- sample(grp.org, length(grp.org))}

# some functions
get.o <- function(v,g=grp.org) {
  # function calculates [o]bservation, here sum of group 1 minus sum of
group 2 for a row of AB as vector
  # input: [v]ector, [g]roup labels, output: single value
  return(sum(v[which(g==1)]) - sum(v[which(g==2)]))
}
get.d <- function(v,s=grp.rnd) {
  # function calculates observations for reshuffled observations
  # input: [v]ector, re[s]huffled array Rx100, output: [d]istribution as
vector of R reshuffled observations
  return(apply(s, MARGIN=1, function(a) {get.o(v, a)}))
}
get.p <- function(o,d) {
  # function calculates onse-sided p-value
  # input: [o]riginal observation and reshuffled null [d]istribution,
output: single [p] value
  r <- sum(d <= o) / length(d)
  p <- ifelse(r < 0.5, r, 1-r)
  return(p)
}

# test on single cols
test.o <- get.o(AB[,1]); test.o
test.d <- get.d(AB[,1]); test.d
test.p <- get.p(test.o,test.d); test.p

# init big matrices for results
oBM <- big.matrix(nrow=n, ncol=1, type='integer', init=NULL)
dBM <- big.matrix(nrow=R, ncol=n, type='integer', init=NULL)
pBM <- big.matrix(nrow=n, ncol=1, type='double', init=NULL)

gc()
system.time({
  oBM <- foreach(i=1:n, .combine=c    ) %dopar% { get.o(AB[,i]) }
  dBM <- foreach(i=1:n, .combine=cbind) %dopar% { get.d(AB[,i]) }
  pBM <- foreach(i=1:n, .combine=c    ) %dopar% { get.p(oBM[i],dBM[,i]) }
})

### CODE END ###


### ATTEMPT 2: ###

iblkcol <- function(a, chunks) {
n <- ncol(a)
i <- 1
nextElem <- function() {
if (chunks <= 0 || n <= 0) stop('StopIteration')
m <- ceiling(n / chunks)
r <- seq(i, length=m)
i <<- i + m
n <<- n - m
chunks <<- chunks - 1
a[,r, drop=FALSE]
}
structure(list(nextElem=nextElem), class=c('iblkcol', 'iter'))
}
nextElem.iblkcol <- function(obj) obj$nextElem()

system.time({
  oBM <- foreach(these=iblkcol(AB,chunks=4), .combine=c) %dopar% {
    foreach(i=1:ncol(these), .combine=c) %do% { get.o(these[,i]) }}
  dBM <- foreach(these=iblkcol(AB,chunks=4), .combine=cbind) %dopar% {
   foreach(i=1:ncol(these), .combine=cbind) %do% { get.d(these[,i]) }}
  pBM <- foreach(i=1:n, .combine=c    ) %dopar% { get.p(oBM[i],dBM[,i]) }
})

	[[alternative HTML version deleted]]

