# Packages needed for parallel computation
# library(doParallel)
# library(foreach) <-- should be installed by doParallel
# library(parallel) <-- should be installed by doParallel
# Useful commands
# detectCores() <-- number of available workers
# cl <- makeCluster(4) <-- make a cluster of 4 workers
# registerDoParallel(cl) <-- register the cluster
# ptime <- system.time({ CODE HERE })[3] <-- a wrapper for timing code
# getDoParWorkers() <-- verify number of parallel workers do Parallel will use
# stopCluster(cl) <-- obvious
############ Parallel Subset Selection with CRN ############
SubsetParallel <- function(k, alpha, n, seed){
# function to do subset selection with CRN in parallel
# k = number of systems
# n = number of replications (equal)
# 1-alpha = confidence level
# seed = random number seed
# simulate:
MySim <- MySim # doParallel needs function to be local
# parallel simulation
Yall <- foreach(x=1:k, .combine=cbind) %dopar% {MySim(x, n, seed)}
# subset selection
Ybar <- apply(Yall, 2, mean)
S2 <- cov(Yall)/n
tval <- qt(1-alpha/(k-1), df = n-1)
Subset <- 1:k
for (i in 1:k){
for (j in 1:k){
if (Ybar[i] < (Ybar[j]-tval*sqrt(S2[i,i] + S2[j,j] - 2*S2[i,j]))){
Subset[i] <- 0
break
}
}
}
list(Subset = Subset[Subset != 0], Ybar = Ybar, S2 = S2, corr=cor(Yall))
}
############ Parallel bi-PASS ############
bipassParallel <- function(k, c, n0, Nmax){
# synchronized biPASS with pooled variance
# k = number of systems
# c = constant needed to guarantee EFER
# n0 = first-stage sample size
# Nmax = maximum number of replications
MySim <- MySim # make the simulation local
g <- function(t, calpha=c){sqrt((calpha + log(t+1))*(t+1))}
II <- 1:k
Active <- rep(TRUE, k) # systems not eliminated
Elim <- rep(0, k) # rep on which elimination occurs
# get n0 reps from each system and compute pooled variance
Yn0 <- foreach(i=1:k, .combine=rbind) %dopar% {MySim(i, n0)}
S2 <- mean(apply(Yn0,1,var))
# start sequential elimination
Ysum <- apply(Yn0, 1, sum)
r <- n0
N <- n0*k
# main elimination loop
while(sum(Active) > 1 && N < Nmax){
r <- r + 1
N <- N + sum(Active)
Ynew <- foreach(i=II[Active], .combine=rbind) %dopar% {MySim(i)}
Ysum[Active] <- Ysum[Active] + Ynew
rmuhat <- sum(Ysum[Active])/sum(Active)
Survivors <- foreach(l=II[Active], .combine=rbind) %dopar% {Ysum[l] - rmuhat > -g(r/S2)*S2}
Elim[Active][Survivors == FALSE] <- r
Active[Active] <- Survivors
}
list(Best = II[Active], n = r, Elim=Elim, Means = Ysum[Active]/r)
}
############ Better Parallel bi-PASS ############
bipassFast <- function(k, c, n0, dn, Nmax){
# synchronized biPASS with pooled variance
# k = number of systems
# c = constant needed to guarantee EFER
# n0 = first-stage sample size
# dn = batch size per run
# Nmax = maximum number of replications
MySim <- MySim # make the simulation local
g <- function(t, calpha=c){sqrt((calpha + log(t+1))*(t+1))}
II <- 1:k
Active <- rep(TRUE, k) # systems not eliminated
Elim <- rep(0, k) # rep on which elimination occurs
# get n0 reps from each system and compute pooled variance
Yn0 <- foreach(i=1:k, .combine=rbind) %dopar% {MySim(i, n0)}
S2 <- mean(apply(Yn0,1,var))
# start sequential elimination
Ysum <- apply(Yn0, 1, sum)
r <- n0
N <- n0*k
# main elimination loop
while(sum(Active) > 1 && N < Nmax){
r <- r + dn
N <- N + dn*sum(Active)
Ynew <- foreach(i=II[Active], .combine=rbind) %dopar% {MySim(i, dn)}
Ysum[Active] <- Ysum[Active] + apply(Ynew, 1, sum)
rmuhat <- sum(Ysum[Active])/sum(Active)
for(l in II[Active]){
if(Ysum[l] - rmuhat <= -g(r/S2)*S2){
Active[l] <- FALSE
Elim[l] <- r
}
}
}
list(Best = II[Active], n = r, Elim=Elim, Means = Ysum[Active]/r)
}
############ Better Non-Parallel bi-PASS ############
bipassSlow <- function(k, c, n0, dn, Nmax){
# synchronized biPASS with pooled variance
# k = number of systems
# c = constant needed to guarantee EFER
# n0 = first-stage sample size
# dn = batch size per run
# Nmax = maximum number of replications
MySim <- MySim # make the simulation local
g <- function(t, calpha=c){sqrt((calpha + log(t+1))*(t+1))}
II <- 1:k
Active <- rep(TRUE, k) # systems not eliminated
Elim <- rep(0, k) # rep on which elimination occurs
# get n0 reps from each system and compute pooled variance
Yn0 <- NULL
for (i in 1:k){
Yn0 <- rbind(Yn0, MySim(i, n0))}
S2 <- mean(apply(Yn0,1,var))
# start sequential elimination
Ysum <- apply(Yn0, 1, sum)
r <- n0
N <- n0*k
# main elimination loop
while(sum(Active) > 1 && N < Nmax){
r <- r + dn
N <- N + dn*sum(Active)
Ynew <- NULL
for (i in II[Active]){
Ynew <- rbind(Ynew, MySim(i, dn))}
Ysum[Active] <- Ysum[Active] + apply(Ynew, 1, sum)
rmuhat <- sum(Ysum[Active])/sum(Active)
for(l in II[Active]){
if(Ysum[l] - rmuhat <= -g(r/S2)*S2){
Active[l] <- FALSE
Elim[l] <- r
}
}
}
list(Best = II[Active], n = r, Elim=Elim, Means = Ysum[Active]/r)
}