install.packages('Rmpi')
require(Rmpi)

#######################################################################
hello <- function(index){
  if(is.element(mpi.comm.rank(), index))
    print(paste("Hello ya'll from processsor ", mpi.comm.rank(),
                "!", sep = ""))
  else
    print(paste("Processor", mpi.comm.rank(), "has nothing to say!"))
}
#######################################################################
myRegTask <- function(mySubSets){
  n <- length(y)
  aicVec <- apply(mySubSets, MARGIN = 1, function(a){
        ind <- as.vector(a, mode = "logical")
        if(sum(a) == 0)
          extractAIC(lm(y ~ 1), k = log(n))
        else
          extractAIC(lm(y ~ X[, ind]), k = log(n))})
  location <- which.min(aicVec[2, ])
  list(mySubSets[location, ], aicVec[, location])
}

#create the data
set.seed(123)
n <- 1000
p <- 10
nDelete <- sample(1:p, size = 1)
indicesToDelete <- sample(1:p, nDelete)
beta <- rnorm(p, sd = 2)
beta[indicesToDelete] <- 0
beta
X <- matrix(rnorm(p*n), n, p)
y <- X%*%beta + rnorm(n)

#set up the subsets
nSubSets <- 2^{p - 2}
subSets <- expand.grid(rep(list(0:1), p))
indexList <- lapply(1:4, FUN = function(i)
subSets[((i - 1)*nSubSets + 1):(i *nSubSets), ])
for(i in 1:3) mpi.send.Robj(indexList[[i + 1]], i, 0)
mpi.bcast.cmd(mySubSets <- mpi.recv.Robj(mpi.any.source(),
mpi.any.tag()))
mpi.bcast.Robj2slave(y)
mpi.bcast.Robj2slave(X)

#send the function to the processors
mpi.bcast.Robj2slave(myRegTask)

#now everybody analyzes their data
mpi.remote.exec(myRegTask(mySubSets))
myRegTask(indexList[[1]])
#######################################################################
library(multicore)
timeVec <- vector(length = 2, mode = "numeric")
for(i in 1:100){
  timeVec[1] <- timeVec[1] +
        system.time(mclapply(indexList, myRegTask))[[3]]
  timeVec[2] <- timeVec[2] +
        system.time(myRegTask(subSets))[[3]]
}
#######################################################################
library(rstream)
library(Rmpi)
mpi.spawn.Rslaves(nslaves = 3)
#create a list of rstream.mrg32k3a objects
rngList <- vector(length = 4, mode = "list")
rngList <- c(new("rstream.mrg32k3a", seed = c(1, 2, 3, 4, 5, 6)),
replicate(3, new("rstream.mrg32k3a")))
#pack the objects that will be sent to the processors
for(i in 2:4) rstream.packed(rngList[[i]]) <- TRUE
#load the package and send the packed rstream objects
mpi.bcast.cmd(library(rstream))
for(i in 1:3) mpi.send.Robj(rngList[[i + 1]], i, 0)
mpi.bcast.cmd(myRng
<- mpi.recv.Robj(mpi.any.source(), mpi.any.tag()))
#unpack the objects on the processors
mpi.bcast.cmd(rstream.packed(myRng) <- FALSE)
#now generate the random numbers
mpi.remote.exec(rstream.sample(myRng, 1))
rstream.sample(rngList[[1]], 1)
#######################################################################
library(rstream)
library(multicore)
rngList <- vector(length = 4, mode = "list")
rngList <- c(new("rstream.mrg32k3a", seed = c(1, 2, 3, 4, 5, 6)),
replicate(3, new("rstream.mrg32k3a")))
mclapply(rngList, function(a) rstream.sample(a, 1))
#######################################################################