# Number of CPUs to use:

ncpus <- 150


# Make scenario grid:

scenariogrid <- expand.grid(iteration=1:50, nbf=1:3)[,2:1]

# Generate seeds:

set.seed(1234)
seeds <- sample(1000:10000000, size=nrow(scenariogrid))

scenariogrid$seed <- seeds



# This function performs the CV error estimation.
# It takes the covariate matrix 'X' and the target variable
# 'y' and performs LDA using PLS components
# as covariates tuning the number of components by 3-fold CV
# repeated two times.

myclassificator <- function(X, y) {

require("CMA")

K <- 5
niter <- 2

cvdatasets <- GenerateLearningsets(y = y, method = "CV", fold = K, niter = niter)

comp.grid <- 1:10

tunenestedcv <- tune(X, y, learningsets = cvdatasets, classifier = pls_ldaCMA, 
  grids = list(comp = comp.grid))

nestedcrossval <- classification(X, y, learningsets = cvdatasets, 
  classifier = pls_ldaCMA, tuneres = tunenestedcv)
misclassnestedcv <- evaluation(nestedcrossval, measure = "misclassification")

errortune <- matrix(misclassnestedcv@score, nrow=K, ncol=niter)

mean(errortune)

}



# This function simulates one dataset, performs SVA and FAbatch and
# estimates the CV-error using the results.

performAnalysis <- function(z) {

require("bapred")

seed <- scenariogrid$seed[z]
nbf <- scenariogrid$nbf[z]

set.seed(seed)

X <- matrix(nrow=40, ncol=1000, data=rnorm(40*1000))
batch <- factor(rep(1:2, times=c(20,20)))
y <- factor(sample(1:2, size=40, replace=TRUE))

svaparams <- ba(x=X, y=y, batch=batch, method = "sva", nbf=nbf)
fabatchparams <- ba(x=X, y=y, batch=batch, method = "fabatch", nbf=nbf)

misclasswithoutbr <- myclassificator(X, y)
misclassaftersva <- myclassificator(svaparams$xadj, y)
misclassafterfabatch <- myclassificator(fabatchparams$xadj, y)

return(list(withoutbr=misclasswithoutbr, sva=misclassaftersva, fabatch=misclassafterfabatch))

}



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


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

# Export all objects to the CPUs:
sfExportAll()

# Perform the simulation:
Results <- sfLapply(1:nrow(scenariogrid), function(z) try({performAnalysis(z)}))

# Stop the cluster:
sfStop()

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

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