# regresses response variable Y column against 
# all possible subsets of the Xi predictor variables, 
# with subset size up through k; returns the
# adjusted R-squared value for each subset

# this version computes X'X and X'Y first, and stores it at the workers

# scheduling methods:
#
#   static (clusterApply())
#   dynamic (clusterApplyLB())
#   reverse the order of the tasks
#   varying chunk size (in dynamic case)

# arguments:
#    cls:  cluster
#    x:  matrix of predictors, one per column
#    y:  vector of the response variable
#    k:  max size of predictor set
#    reverse:  True means reverse the order of the iterations
#    dyn:  True means dynamic scheduling
#    chunk:  chunk size 
# return value:
#    R matrix, showing adjusted R-squared values, 
#    indexed by predictor set

snowapr1 <- function(cls,x,y,k,reverse=F,dyn=F,chunk=1) {
   # add 1s column
   x <- cbind(1,x)
   xpx <- crossprod(x,x)
   xpy <- crossprod(x,y)
   p <- ncol(x) - 1
   # generate matrix of predictor subsets
   allcombs <- genallcombs(p,k)
   ncombs <- length(allcombs)
   clusterExport(cls,"do1pset1")
   clusterExport(cls,"linregadjr2")
   # set up task indices 
   tasks <- if (!reverse) seq(1,ncombs,chunk) else 
      seq(ncombs,1,-chunk) 
   if (!dyn) {
      out <- clusterApply(cls,tasks,dochunk1,
         x,y,xpx,xpy,allcombs,chunk)
   } else {
      out <- clusterApplyLB(cls,tasks,dochunk1,
         x,y,xpx,xpy,allcombs,chunk)
   }
   Reduce(rbind,out)
}

# generate all nonempty subsets of 1..p of size <= k; 
# returns a list, one element per predictor set
genallcombs <- function(p,k) {
   allcombs <- list()
   for (i in 1:k) {
      tmp <- combn(1:p,i)
      allcombs <- c(allcombs,matrixtolist(tmp,rc=2))
   }
   allcombs
}

# extracts rows (rc=1) or columns (rc=2) of a matrix, producing a list
matrixtolist <- function(rc,m) {
   if (rc == 1) {
      Map(function(rownum) m[rownum,],1:nrow(m))
   } else Map(function(colnum) m[,colnum],1:ncol(m))
}

# process all the predictor sets in the chunk 
# whose first index is psetstart
dochunk1 <- function(psetstart,x,y,xpx,xpy,allcombs,chunk) {
   ncombs <- length(allcombs)
   lasttask <- min(psetstart+chunk-1,ncombs)
   t(sapply(allcombs[psetstart:lasttask],do1pset1,x,y,xpx,xpy))
}

# find the adjusted R-squared values for the given 
# predictor set index
do1pset1 <- function(onepset,x,y,xpx,xpy) {
   ps <- c(1,onepset+1)  # account for constant term
   x1 <- x[,ps]
   xpx1 <- xpx[ps,ps]
   xpy1 <- xpy[ps]
   ar2 <- linregadjr2(x1,y,xpx1,xpy1)
   n0s <- ncol(x) - length(ps)
   # form the report for this predictor set; need trailngs 0s so as to
   # form matrices of uniform numbers of rows, to use rbind() in
   # snowapr()
   c(ar2,onepset,rep(0,n0s))
}

linregadjr2 <- function(x,y,xpx,xpy) {
   bhat <- solve(xpx,xpy)
   resids <- y - x %*% bhat
   r2 <- 1 - sum(resids^2)/sum((y-mean(y))^2)
   n <-nrow(x); p <- ncol(x) - 1
   1 - (1-r2) * (n-1) / (n-p-1)
}

# predictor set seems best
test <- function(cls,n,p,k,chunk=1,dyn=F,rvrs=F) {
   gendata(n,p)
   snowapr1(cls,x,y,k,rvrs,dyn,chunk)
}

gendata <- function(n,p) {
   x <<- matrix(rnorm(n*p),ncol=p)
   y <<- x%*%c(rep(0.5,p)) + rnorm(n)
}