# 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)
}