# Number of iterations and CPUs to use:

niter <- 500; ncpus <- 250

# Make scenario grid for simulated data:

simulgrid <- expand.grid(iteration=1:niter, indscen=c("b0scen", "bjscen", "btildescen", 
  "CommonCor", "BatchspecificCor", "BatchClassspecificCor"), stringsAsFactors = FALSE)
simulgrid <- simulgrid[,ncol(simulgrid):1]


# Different seed for each simulated dataset:

set.seed(1234)
seeds <- sample(1000:10000000, size=length(unique(simulgrid$indscen))*niter)


# Make scenario grid for simulated datasets:
  
sceniterind <- as.numeric(factor(paste(simulgrid$iteration, simulgrid$indscen, sep="_")))
simulgrid$seed <- sapply(sceniterind, function(x) seeds[unique(sceniterind)==x])
simulgrid$datasetind <- 1:nrow(simulgrid)
nsimul <- nrow(simulgrid)

save(simulgrid, file="./FAbatchPaper/InterimResults/simulgrid.Rda")


# Make scenario grid for batch effect removal:

batchremovalgrid <- expand.grid(batchremmethod=c("none", "fabatch", "combat", "sva", "meancenter", 
  "standardize", "ratioa", "ratiog"), datasetind=simulgrid$datasetind, stringsAsFactors = FALSE)
batchremovalgrid <- batchremovalgrid[,ncol(batchremovalgrid):1]
batchremovalgrid$datasetbrind <- 1:nrow(batchremovalgrid)
batchremovalgrid$iteration <- sapply(batchremovalgrid$datasetind, function(x) simulgrid$iteration[simulgrid$datasetind==x])
batchremovalgrid$indscen <- sapply(batchremovalgrid$datasetind, function(x) simulgrid$indscen[simulgrid$datasetind==x])
nbatchremoval <- nrow(batchremovalgrid)

save(batchremovalgrid, file="./FAbatchPaper/InterimResults/batchremovalgrid.Rda")


# Make scenario grid for applying different metrics:

metricgrid <- expand.grid(metric=c("sep", "avedist", "kldist", "skew", "pvca", 
  "diffexpr", "cor"), datasetbrind=batchremovalgrid$datasetbrind, stringsAsFactors = FALSE)
metricgrid <- metricgrid[,ncol(metricgrid):1]
nmetric <- nrow(metricgrid)

save(metricgrid, file="./FAbatchPaper/InterimResults/metricgrid.Rda")

rm(simulgrid, batchremovalgrid, metricgrid)
gc()


# Source functions for performing the analysis on the cluster:

source("./FAbatchPaper/Functions/SimulationSnowfallFunctions.R")


# Load 'snowfall'-Package for parallel computation:
require("snowfall")


# Initiate cluster:
sfInit(parallel=TRUE, cpus=ncpus, type = "MPI")

# Export all objects to the CPUs:
sfExportAll()

# Simulate datasets:
SimulationDatasets <- sfLapply(1:nsimul, function(z) simulateData(z))

# Perform batch effect removal:
SimulationBatcheffremoval <- sfLapply(1:nbatchremoval, function(z) performBatchremoval(z))

# Apply metrics:
SimulationResults <- sfLapply(1:nmetric, function(z) evaluateMetric(z))

# Stop the cluster:
sfStop()

# Save results:
save(SimulationResults, file="./FAbatchPaper/Results/SimulationResults.Rda")

# Clear the workspace:
rm(list=ls()); gc()
