library("CMA")

# The CMA-function 'scdaCMA' performing Nearest Shrunken Centroids doesn't return a 'model-object', 
# so independent observations cannot be predicted.
# Because of this the function was modified to return only a function to
# predict new observations:

scdaCMApred <- function (X, y, f, learnind, delta = 0.5, models = FALSE, ...) 
{
    nrx <- nrow(X)
    ly <- length(y)
    if (nrx != length(y)) 
        stop("Number of rows of 'X' must agree with length of y \n")
    if (missing(learnind)) 
        learnind <- 1:nrx
    if (length(learnind) > nrx) 
        stop("length of 'learnind' must be smaller than the number of observations. \n")
    y <- as.factor(y)
    K <- nlevels(y)
    levels(y) <- 1:K
    if (K > 2) 
        mode <- "multiclass"
    else mode <- "binary"
    y <- as.numeric(y) - 1
    Ylearn <- y[learnind]
    Xlearn <- X[learnind, ]
    centroids <- matrix(nrow = ncol(X), ncol = K)
    variances <- matrix(nrow = ncol(X), ncol = K)
    mk <- as.numeric(sqrt(1/table(Ylearn) + 1/length(learnind)))
    for (k in 1:K) {
        indk <- (Ylearn == (k - 1))
        centroids[, k] <- muk <- colMeans(Xlearn[indk, , drop = FALSE])
        variances[, k] <- colSums((sweep(Xlearn[indk, , drop = FALSE], 
            2, muk, "-"))^2)
    }
    variances <- rowSums(variances)/(length(learnind) - K)
    overallmu <- colMeans(Xlearn)
    D <- sweep(sweep(sweep(centroids, 1, overallmu, "-"), 2, 
        mk, "/"), 1, variances^(-0.5), "*")
    shrunkD <- sign(D) * pmax(abs(D) - delta, 0)
    shrunkcentroids <- sweep(sweep(sweep(shrunkD, 2, mk, "*"), 
        1, variances^0.5, "*"), 1, overallmu, "+")
    priors <- as.numeric(-2 * log(table(Ylearn)/sum(table(Ylearn))))
    Xtest <- X[-learnind, , drop = FALSE]
    if (nrow(Xtest) == 0) {
        Xtest <- Xlearn
        y <- Ylearn
    }
    else y <- y[-learnind]
    dist <- t(apply(Xtest, 1, function(z) colSums(sweep(sweep(-shrunkcentroids, 
        1, z, "+")^2, 1, variances, "/")) + priors))
    dist <- sweep(dist, 1, rowMeans(dist), "-")
    prob <- safeexp(-0.5 * dist)
    prob <- sweep(prob, 1, rowSums(prob), "/")
    yhat <- apply(prob, 1, which.max) - 1
    if (models == TRUE) 
        modd <- list(NULL)
    if (models == FALSE) 
        modd <- list(NULL)
    new("cloutput", yhat = yhat, y = y, learnind = learnind, 
        prob = prob, method = "scDA", mode = mode, model = modd)
		
		
  # This function returned by 'scdaCMApred' can be used to classify new observations:
  function(xnew) {
    dist <- t(apply(xnew, 1, function(z)   colSums(sweep(sweep(-shrunkcentroids, 1, z, "+")^2, 1, variances, "/")) + priors))
    dist <- sweep(dist, 1, rowMeans(dist), "-")
    prob <- safeexp(-0.5 * dist)
    prob <- sweep(prob, 1, rowSums(prob), "/")
    yhat <- apply(prob, 1, which.max)
    yhat
  }		
		
}

# To make the functions available inside 'scdaCMApred', we have to use the same environment
# as 'scdaCMA', in which these functions are contained:
environment(scdaCMApred) <- environment(scdaCMA)





# This function downloads the raw datasets.
# It takes the whole number 'z', which corresponds to the z-th line of 'scenariogrid',
# which contains the accession number of the respective dataset and a seed.

DownloadData <- function(z) {

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

set.seed(seed)
datasetsfolder <- "./ICSVCodeData/Datasets/"

if(!file.exists(datasetsfolder))
  dir.create(datasetsfolder)

datasetfolder <- paste(datasetsfolder, accnumber, sep="")

if(!file.exists(datasetfolder))
  dir.create(datasetfolder)

library("ArrayExpress")
getAE(accnumber, path = datasetfolder, type = "raw")

}






# This functions performs the normalization on the datasets
# and stores the normalized datasets as well as the information
# necessary for addon normalization.
# It takes the whole number 'z', which corresponds to the z-th line 
# of 'scenariogrid', which contains the necessary information
# on the datasets to normalize.

NormalizeTraining <- function(z) {

# Information on the dataset to normalize:

accnumber <- scenariogrid$accnumber[z] 
seed <- scenariogrid$seed[z]
size <- scenariogrid$size[z]

set.seed(seed)
datasetsfolder <- "./ICSVCodeData/Datasets/"

datasetfolder <- paste(datasetsfolder, accnumber, sep="")


# Load data info:
  
load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber, "_datainfo.Rda", sep=""))

celfiles <- paste(datasetfolder, "/", datainfo$'Array Data File', sep="")



# Setting with large training/test dataset:

if(size=="large") {

# Folder in which the information necessary for addon normalization is stored:

RMAparamfolderLarge <- "./ICSVCodeData/Datasets/RMAparamfilesLarge"
if(!file.exists(RMAparamfolderLarge))
  dir.create(RMAparamfolderLarge) 
  
# Folder in which the normalized datasets are stored::
  
DatanormfolderLarge <- "./ICSVCodeData/Datasets/DatanormfilesLarge"
if(!file.exists(DatanormfolderLarge))
  dir.create(DatanormfolderLarge) 

  
# Subset if more than 70 observations:
if(length(celfiles) > 70) {
  largeind <- sort(sample(1:length(celfiles), size=70))
}
else {
  largeind <- 1:length(celfiles)
}

# Subset dataset:
celfileslarge <- celfiles[largeind]
datainfolarge <- datainfo[largeind,]

# Read CEL-files:
library("affy")
aboall <- ReadAffy(filenames=celfileslarge)

# Perform normalization:
library("bapred")
rmaparams <- rmatrain(aboall)
y <- factor(as.numeric(datainfolarge$sexvar=="female")+1, levels=c(1,2))

save(rmaparams, y, file=paste(RMAparamfolderLarge, "/", accnumber, "_rmaparam.Rda", sep=""))

X <- rmaparams$xnorm

save(X, y, file=paste(DatanormfolderLarge, "/", accnumber, "_datanorm.Rda", sep=""))
  
}



# Setting with small training/test dataset:

if(size=="small") {

RMAparamfolderSmall <- "./ICSVCodeData/Datasets/RMAparamfilesSmall"
if(!file.exists(RMAparamfolderSmall))
  dir.create(RMAparamfolderSmall) 

DatanormfolderSmall <- "./ICSVCodeData/Datasets/DatanormfilesSmall"
if(!file.exists(DatanormfolderSmall))
  dir.create(DatanormfolderSmall) 


# Draw 20 observations under the restriction that the smaller
# class must have at least five observations:

tabsexvar <- table(datainfo$sexvar)
smallclass <- names(tabsexvar)[which.min(tabsexvar)]

indsmall <- which(datainfo$sexvar==smallclass)
obligfive <- sample(indsmall, size=5)

remainind <- sample(setdiff(1:nrow(datainfo), obligfive), size=20-length(obligfive))

smallind <- sort(c(obligfive, remainind))

celfilessmall <- celfiles[smallind]
datainfosmall <- datainfo[smallind,]


library("affy")
aboall <- ReadAffy(filenames=celfilessmall)


library("bapred")
rmaparams <- rmatrain(aboall)
y <- factor(as.numeric(datainfosmall$sexvar=="female")+1, levels=c(1,2))

save(rmaparams, y, file=paste(RMAparamfolderSmall, "/", accnumber, "_rmaparam.Rda", sep=""))

X <- rmaparams$xnorm

save(X, y, file=paste(DatanormfolderSmall, "/", accnumber, "_datanorm.Rda", sep=""))

}



# Setting with very small test dataset:

if(size=="verysmall") {

RMAparamfolderVerysmall <- "./ICSVCodeData/Datasets/RMAparamfilesVerysmall"
if(!file.exists(RMAparamfolderVerysmall))
  dir.create(RMAparamfolderVerysmall) 

DatanormfolderVerysmall <- "./ICSVCodeData/Datasets/DatanormfilesVerysmall"
if(!file.exists(DatanormfolderVerysmall))
  dir.create(DatanormfolderVerysmall) 


# Draw 5 observations:

verysmallind <- sample(1:nrow(datainfo), size=5)

celfilesVerysmall <- celfiles[verysmallind]
datainfoVerysmall <- datainfo[verysmallind,]


library("affy")
aboall <- ReadAffy(filenames=celfilesVerysmall)

library("bapred")
rmaparams <- rmatrain(aboall)
y <- factor(as.numeric(datainfoVerysmall$sexvar=="female")+1, levels=c(1,2))

save(rmaparams, y, file=paste(RMAparamfolderVerysmall, "/", accnumber, "_rmaparam.Rda", sep=""))

X <- rmaparams$xnorm

save(X, y, file=paste(DatanormfolderVerysmall, "/", accnumber, "_datanorm.Rda", sep=""))

}


}






# This functions performs the addon normalization
# It takes the whole number 'z', which corresponds to the z-th line 
# of 'scenariogrid', which contains the necessary information
# on the datasets to addon normalize.

NormalizeAddonTest <- function(z) {

# Information on the dataset to normalize:

datasettrainid <- scenariogrid$accnumbertrain[z] 
datasettestid <- scenariogrid$accnumbertest[z] 
seed <- scenariogrid$seed[z]
size <- scenariogrid$size[z]

datasetsfolder <- "./ICSVCodeData/Datasets/"

datasettestfolder <- paste(datasetsfolder, datasettestid, sep="")



library("affy")


RMAparamfolderLarge <- "./ICSVCodeData/Datasets/RMAparamfilesLarge"
RMAparamfolderSmall <- "./ICSVCodeData/Datasets/RMAparamfilesSmall"



# Setting with large test dataset:

if(size=="large") {


# Large training dataset:

set.seed(seed)

# Folders in which the data after addon normalization is stored:
DatanormfoldertestaddonLargeTrainLargeTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonLargeTrainLargeTest"
if(!file.exists(DatanormfoldertestaddonLargeTrainLargeTest))
  dir.create(DatanormfoldertestaddonLargeTrainLargeTest) 
  

# Load information necessary for addon normalization:

load(paste(RMAparamfolderLarge, "/", datasettrainid, "_rmaparam.Rda", sep=""))


# Subset test dataset:

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))
celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")

if(length(celfiles) > 70) {
  largeind <- sort(sample(1:length(celfiles), size=70))
}
else {
  largeind <- 1:length(celfiles)
}

celfileslarge <- celfiles[largeind]
datainfolarge <- datainfo[largeind,]


# Read CEL-files of test data:

abotest <- ReadAffy(filenames=celfileslarge)


# Perform addon normalization:

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfolarge$sexvar=="female")+1, levels=c(1,2))

save(Xtest, ytest, file=paste(DatanormfoldertestaddonLargeTrainLargeTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))



# Small training dataset:

set.seed(seed)

DatanormfoldertestaddonSmallTrainLargeTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonSmallTrainLargeTest"
if(!file.exists(DatanormfoldertestaddonSmallTrainLargeTest))
  dir.create(DatanormfoldertestaddonSmallTrainLargeTest)  


load(paste(RMAparamfolderSmall, "/", datasettrainid, "_rmaparam.Rda", sep=""))

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))
celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")

if(length(celfiles) > 70) {
  largeind <- sort(sample(1:length(celfiles), size=70))
}
else {
  largeind <- 1:length(celfiles)
}

celfileslarge <- celfiles[largeind]
datainfolarge <- datainfo[largeind,]


abotest <- ReadAffy(filenames=celfileslarge)

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfolarge$sexvar=="female")+1, levels=c(1,2))

save(Xtest, ytest, file=paste(DatanormfoldertestaddonSmallTrainLargeTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))

}
  
  
if(size=="small") {

set.seed(seed)

DatanormfoldertestaddonLargeTrainSmallTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonLargeTrainSmallTest"
if(!file.exists(DatanormfoldertestaddonLargeTrainSmallTest))
  dir.create(DatanormfoldertestaddonLargeTrainSmallTest) 
 

load(paste(RMAparamfolderLarge, "/", datasettrainid, "_rmaparam.Rda", sep=""))

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))

celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")


smallind <- sort(sample(1:length(celfiles), size=20))
while(min(table(datainfo$sexvar[smallind])) < 2)
  smallind <- sort(sample(1:length(celfiles), size=20))

celfilessmall <- celfiles[smallind]
datainfosmall <- datainfo[smallind,]


abotest <- ReadAffy(filenames=celfilessmall)

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfosmall$sexvar=="female")+1, levels=c(1,2))


save(Xtest, ytest, file=paste(DatanormfoldertestaddonLargeTrainSmallTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))


set.seed(seed)

DatanormfoldertestaddonSmallTrainSmallTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonSmallTrainSmallTest"
if(!file.exists(DatanormfoldertestaddonSmallTrainSmallTest))
  dir.create(DatanormfoldertestaddonSmallTrainSmallTest)  


load(paste(RMAparamfolderSmall, "/", datasettrainid, "_rmaparam.Rda", sep=""))

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))
celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")

smallind <- sort(sample(1:length(celfiles), size=20))
while(min(table(datainfo$sexvar[smallind])) < 2)
  smallind <- sort(sample(1:length(celfiles), size=20))

celfilessmall <- celfiles[smallind]
datainfosmall <- datainfo[smallind,]


abotest <- ReadAffy(filenames=celfilessmall)

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfosmall$sexvar=="female")+1, levels=c(1,2))

save(Xtest, ytest, file=paste(DatanormfoldertestaddonSmallTrainSmallTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))
  
}
  
  

if(size=="verysmall") {

set.seed(seed)

DatanormfoldertestaddonLargeTrainVerysmallTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonLargeTrainVerysmallTest"
if(!file.exists(DatanormfoldertestaddonLargeTrainVerysmallTest))
  dir.create(DatanormfoldertestaddonLargeTrainVerysmallTest) 
 

load(paste(RMAparamfolderLarge, "/", datasettrainid, "_rmaparam.Rda", sep=""))

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))

celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")



verysmallind <- sample(1:nrow(datainfo), size=5)

celfilesVerysmall <- celfiles[verysmallind]
datainfoVerysmall <- datainfo[verysmallind,]


abotest <- ReadAffy(filenames=celfilesVerysmall)

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfoVerysmall$sexvar=="female")+1, levels=c(1,2))


save(Xtest, ytest, file=paste(DatanormfoldertestaddonLargeTrainVerysmallTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))



set.seed(seed)

DatanormfoldertestaddonSmallTrainVerysmallTest <- "./ICSVCodeData/Datasets/DatanormfilestestaddonSmallTrainVerysmallTest"
if(!file.exists(DatanormfoldertestaddonSmallTrainVerysmallTest))
  dir.create(DatanormfoldertestaddonSmallTrainVerysmallTest)    

  
load(paste(RMAparamfolderSmall, "/", datasettrainid, "_rmaparam.Rda", sep=""))

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", datasettestid, "_datainfo.Rda", sep=""))
celfiles <- paste(datasettestfolder, "/", datainfo$'Array Data File', sep="")

verysmallind <- sample(1:nrow(datainfo), size=5)

celfilesVerysmall <- celfiles[verysmallind]
datainfoVerysmall <- datainfo[verysmallind,]


abotest <- ReadAffy(filenames=celfilesVerysmall)

library("bapred")
Xtest <- rmaaddon(rmaparams, abotest)

ytest <- factor(as.numeric(datainfoVerysmall$sexvar=="female")+1, levels=c(1,2))

save(Xtest, ytest, file=paste(DatanormfoldertestaddonSmallTrainVerysmallTest, "/", datasettrainid, "train_", datasettestid, "test_normaddon.Rda", sep=""))
  
}
 
}
















# This function is called by "evaluatesetting()", see below. 
# It performs the classification using one of the considered classification
# methods and returns the cross-table between the predictions of the target variable and the true
# values of the target variable as well as the predictions of the target variable.

# Input parameters:

# Xtrain      - covariate matrix of training dataset
# ytrain      - values of the target variable in the training dataset
# Xtest       - covariate matrix of test dataset
# ytest       - values of the target variable in the test dataset
# classmethod - one of the following: "boosting", "boostingvarsel",
#               "plslda", "nsc", "plsldavarsel", "rf", "knnvarsel".
  
performValidation <- function(Xtrain, ytrain, Xtest, ytest, classmethod) {

  require("CMA")

  # Boosting:
  
  if(classmethod=="boosting") {
  
    # Grid used for tuning the number of boosting iterations: 
    mstop.grid <- c(50, 100, 200, 500, 1000)
    
	# Tuning by 3-fold CV is repeated 10 times to obtain more stable results:
    mstopbestvec <- 0
    for(i in 1:10) {
      tunecv <- tune(Xtrain, ytrain, classifier = compBoostCMA, grids = list(mstop = mstop.grid), trace=FALSE)
      mstopbestvec[i] <- mstop.grid[which.min(tunecv@tuneres[[1]])]
    }
    mstopbest <- median(mstopbestvec)
    mstopbest <- mstop.grid[max(which(mstop.grid <= mstopbest))]

    # Learn the prediction rule and evaluate it using the test dataset:
	boostingobject <- compBoostCMA(rbind(Xtrain, Xtest), c(ytrain, ytest), learnind=1:nrow(Xtrain),
                mstop = mstopbest, models=TRUE)
    
	ytesthat <- boostingobject@yhat + 1
	ytesthat <- factor(ytesthat, levels=1:2)

  }

  
  # Boosting with preliminary variable selection:  
  
  if(classmethod=="boostingvarsel") {
  
    # Pre-select the 2000 variables with smallest p-values out of two-sample t-tests:
    pvals <- apply(Xtrain, 2, function(x, y) t.test(x[y==1], x[y==2])$p.value, y = ytrain)
    selind <- order(pvals)[1:2000]

	# Subset variables in trainng and test data correspondingly:
    Xtrain <- Xtrain[,selind]
    Xtest <- Xtest[,selind]
	
	mstop.grid <- c(50, 100, 200, 500, 1000)

    mstopbestvec <- 0
    for(i in 1:10) {
      tunecv <- tune(Xtrain, ytrain, classifier = compBoostCMA, grids = list(mstop = mstop.grid), trace=FALSE)
      mstopbestvec[i] <- mstop.grid[which.min(tunecv@tuneres[[1]])]
    }
    mstopbest <- median(mstopbestvec)
    mstopbest <- mstop.grid[max(which(mstop.grid <= mstopbest))]

	boostingobject <- compBoostCMA(rbind(Xtrain, Xtest), c(ytrain, ytest), learnind=1:nrow(Xtrain),
                mstop = mstopbest, models=TRUE)
    
	ytesthat <- boostingobject@yhat + 1
	ytesthat <- factor(ytesthat, levels=1:2)
	
  }
  
  
  # PLS-LDA:
  
  if(classmethod=="plslda") {
  
    # Grid used for tuning the number of components: 
    comp.grid <- 1:10

    compbestvec <- 0
    for(i in 1:10) {
      tunecv <- tune(Xtrain, ytrain, classifier = pls_ldaCMA, grids = list(comp = comp.grid))
      compbestvec[i] <- comp.grid[which.min(tunecv@tuneres[[1]])]
    }
    compbest <- floor(median(compbestvec))

    plsldaobject <- pls_ldaCMA(rbind(Xtrain, Xtest), c(ytrain, ytest), 
      learnind=1:length(ytrain), comp = compbest)

    ytesthat <- plsldaobject@yhat + 1
    ytesthat <- factor(ytesthat, levels=c(1,2))
  
  }
  
  
  # NSC:
  
  if(classmethod=="nsc") {
    
    # Grid used for tuning the shrinkage parameter: 
    delta.grid <- c(0.05, 0.1, 0.25, 0.5, 1, 1.5)
    
    deltabestvec <- 0
    for(i in 1:10) {
      tunecv <- tune(Xtrain, ytrain, classifier = scdaCMA, grids = list(delta = delta.grid), trace=FALSE)
      deltabestvec[i] <- delta.grid[which.min(tunecv@tuneres[[1]])]
    }
    deltabest <- floor(median(deltabestvec))
    
    classifier <- scdaCMApred(Xtrain, ytrain, delta = deltabest, models=TRUE)
    
    ytesthat <- classifier(Xtest)
	names(ytesthat) <- NULL
    ytesthat <- factor(ytesthat, levels=c(1,2))
    
  }
  
  
  # PLS-LDA with preliminary variable selection:
  
  if(classmethod=="plsldavarsel") {
    
    pvals <- apply(Xtrain, 2, function(x, y) t.test(x[y==1], x[y==2])$p.value, y = ytrain)
    selind <- order(pvals)[1:2000]
    
    Xtrain <- Xtrain[,selind]
    Xtest <- Xtest[,selind]
    
    comp.grid <- 1:10
    
    compbestvec <- 0
    for(i in 1:10) {
      tunecv <- tune(Xtrain, ytrain, classifier = pls_ldaCMA, grids = list(comp = comp.grid))
      compbestvec[i] <- comp.grid[which.min(tunecv@tuneres[[1]])]
    }
    compbest <- floor(median(compbestvec))
    
    plsldaobject <- pls_ldaCMA(rbind(Xtrain, Xtest), c(ytrain, ytest), 
                               learnind=1:length(ytrain), comp = compbest)
    
    ytesthat <- plsldaobject@yhat + 1
    ytesthat <- factor(ytesthat, levels=c(1,2))
    
  }

  
  # Random Forest:
  
  if(classmethod=="rf") {
    
	# Grid used for tuning "mtry": 
    mtry.grid <- ceiling(c(0.25, 1, 2)*sqrt(ncol(Xtrain)))
    
	# For Random Forest we do not repeat the tuning 10 times, but perform it only once
	# in order to reduce computational costs:
    tunecv <- tune(Xtrain, ytrain, classifier = rfCMA, grids = list(mtry = mtry.grid))
    mtrybest <- mtry.grid[which.min(tunecv@tuneres[[1]])]
    
    classifiercma <- rfCMA(Xtrain, ytrain, mtry = mtrybest, ntree=500, models=TRUE)
    randfor <- classifiercma@model[[1]]
    
    ytesthat <- predict(object=randfor, newdata=Xtest)
    names(ytesthat) <- NULL
    
  }
  
  
  # kNN with preliminary variable selection:
  
  if(classmethod=="knnvarsel") {
   
	pvals <- apply(Xtrain, 2, function(x, y) t.test(x[y==1], x[y==2])$p.value, y = ytrain)
	selind <- order(pvals)[1:2000]

	Xtrain <- Xtrain[,selind]
	Xtest <- Xtest[,selind]

    # Grid used for tuning the number of nearest neighbours to be used: 
	k.grid <- 1:10

	kbestvec <- 0
	for(i in 1:10) {
 		 tunecv <- tune(Xtrain, ytrain, classifier = knnCMA, grids = list(k = k.grid), trace=FALSE)
  		kbestvec[i] <- k.grid[which.min(tunecv@tuneres[[1]])]
	}
	kbest <- median(kbestvec)
	kbest <- k.grid[max(which(k.grid <= kbest))]

	knnobject <- knnCMA(rbind(Xtrain, Xtest), c(ytrain, ytest), learnind=1:nrow(Xtrain),
                    k = kbest, models=TRUE)
  
  	ytesthat <- knnobject@yhat + 1
	ytesthat <- factor(ytesthat, levels=1:2)
  
  }
    
  # cross-table between predictions of the target variable
  # and its true values:
  confmatrix <- table(ytesthat, ytest)
  
  # Return cross-table from above, and the predictions of the 
  # values of the target variable:
  return(list(confmatrix=confmatrix, ytesthat=ytesthat))
  
}










# This functions loads the previously stored datasets after (addon)
# normalization, performs the (addon) batch effect adjustment,
# applies the two classification methods and returns the information
# necessary for calculating the MCC-values.
# It takes the whole number 'z', which corresponds to the z-th line 
# of 'scenariogrid', which contains the necessary information
# for evaluating a specific setting for a particular training and test set
# pair.

evaluatesetting <- function(z) {

# Load "bapred"-package:

require("bapred")


# Extract information for the z-th iteration:

setting <- scenariogrid$setting[z]
accnumbertrain <- scenariogrid$accnumbertrain[z]
accnumbertest <- scenariogrid$accnumbertest[z]
seed <- scenariogrid$seed[z]
batchremmethod <- scenariogrid$batchremmethod[z]
traintestsetsize <- scenariogrid$traintestsetsize[z]


# Load training and test sets - those corresponding to 
# the z-th iteration:

datasetsfolder <- "./ICSVCodeData/Datasets"

set.seed(seed)

if(setting=="sepnorm") {

if(traintestsetsize=="trainsmalltestverysmall") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesVerysmall/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}  

if(traintestsetsize=="trainsmalltestsmall") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}  

if(traintestsetsize=="trainsmalltestlarge") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}

if(traintestsetsize=="trainlargetestverysmall") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesVerysmall/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}  

if(traintestsetsize=="trainlargetestsmall") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}  

if(traintestsetsize=="trainlargetestlarge") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertest, "_datanorm.Rda", sep=""))
  Xtest <- X; ytest <- y
  rm(X, y); gc()
}
  
}

if(setting=="addonnorm") {

if(traintestsetsize=="trainsmalltestverysmall") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonLargeTrainVerysmallTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
} 

if(traintestsetsize=="trainsmalltestsmall") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonSmallTrainSmallTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
}  

if(traintestsetsize=="trainsmalltestlarge") {
  load(paste(datasetsfolder, "/DatanormfilesSmall/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonSmallTrainLargeTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
}  

if(traintestsetsize=="trainlargetestverysmall") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonLargeTrainVerysmallTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
} 

if(traintestsetsize=="trainlargetestsmall") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonLargeTrainSmallTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
}  

if(traintestsetsize=="trainlargetestlarge") {
  load(paste(datasetsfolder, "/DatanormfilesLarge/", accnumbertrain, "_datanorm.Rda", sep=""))
  Xtrain <- X; ytrain <- y
  rm(X, y); gc()
  load(paste(datasetsfolder, "/DatanormfilestestaddonLargeTrainLargeTest/", accnumbertrain, "train_", accnumbertest, "test_normaddon.Rda", sep=""))
  ytest <- factor(ytest, levels=c(1,2))
}    

}


# Perform batch effect adjustment on the training
# data and addon batch effect adjustment on the test data:

if(!(batchremmethod %in% c("fsvafast", "fsvaexact"))) {
params <- ba(x=Xtrain, y=ytrain, batch=factor(rep(1, nrow(Xtrain))), method = batchremmethod)
}
else {
if(batchremmethod=="fsvafast")
  params <- ba(x=Xtrain, y=ytrain, batch=factor(rep(1, nrow(Xtrain))), method = "sva", algorithm = "fast")
if(batchremmethod=="fsvaexact")
  params <- ba(x=Xtrain, y=ytrain, batch=factor(rep(1, nrow(Xtrain))), method = "sva", algorithm = "exact")
}

Xtrain <- params$xadj
Xtest <- baaddon(params, x=Xtest, batch=factor(rep(1, nrow(Xtest))))



# Perform validation with Boosting:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "boosting"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1boosting <- confmatrix[1,1]
true1pred2boosting <- confmatrix[2,1]
true2pred1boosting <- confmatrix[1,2]
true2pred2boosting <- confmatrix[2,2]

# Test set predictions:
ytesthatboosting <- traintestval$ytesthat




# Perform validation with Boosting with preliminary variable selection:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "boostingvarsel"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1boostingvarsel <- confmatrix[1,1]
true1pred2boostingvarsel <- confmatrix[2,1]
true2pred1boostingvarsel <- confmatrix[1,2]
true2pred2boostingvarsel <- confmatrix[2,2]

# Test set predictions:
ytesthatboostingvarsel <- traintestval$ytesthat



# Perform validation with PLS-LDA:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "plslda"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1plslda <- confmatrix[1,1]
true1pred2plslda <- confmatrix[2,1]
true2pred1plslda <- confmatrix[1,2]
true2pred2plslda <- confmatrix[2,2]

# Test set predictions:
ytesthatplslda <- traintestval$ytesthat



# Perform validation with NSC:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "nsc"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1nsc <- confmatrix[1,1]
true1pred2nsc <- confmatrix[2,1]
true2pred1nsc <- confmatrix[1,2]
true2pred2nsc <- confmatrix[2,2]

# Test set predictions:
ytesthatnsc <- traintestval$ytesthat



# Perform validation with PLS-LDA with preliminary variable selection:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "plsldavarsel"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1plsldavarsel <- confmatrix[1,1]
true1pred2plsldavarsel <- confmatrix[2,1]
true2pred1plsldavarsel <- confmatrix[1,2]
true2pred2plsldavarsel <- confmatrix[2,2]

# Test set predictions:
ytesthatplsldavarsel <- traintestval$ytesthat




# Perform validation with Random Forest:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "rf"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1rf <- confmatrix[1,1]
true1pred2rf <- confmatrix[2,1]
true2pred1rf <- confmatrix[1,2]
true2pred2rf <- confmatrix[2,2]

# Test set predictions:
ytesthatrf <- traintestval$ytesthat




# Perform validation with kNN with preliminary variable selection:

traintestval <- try(0 + "0")
count <- 0

while((class(traintestval)=="try-error") & (count < 5)) {
  traintestval <- try(performValidation(Xtrain, ytrain, Xtest, ytest, classmethod = "knnvarsel"))
  count <- count + 1
}

confmatrix <- traintestval$confmatrix

# Extract the informations necessary for calculating the MCC-values:
true1pred1knnvarsel <- confmatrix[1,1]
true1pred2knnvarsel <- confmatrix[2,1]
true2pred1knnvarsel <- confmatrix[1,2]
true2pred2knnvarsel <- confmatrix[2,2]

# Test set predictions:
ytesthatknnvarsel <- traintestval$ytesthat



# Return all informations for identifying the specific iteration together with all information 
# necessary for calculating the MCC-values:

res <- list(setting=setting, acctrain=accnumbertrain, acctest=accnumbertest,
  batchremmethod=batchremmethod, traintestsetsize=traintestsetsize, ntrain=nrow(Xtrain), propclass2=mean(ytrain==2), ytest=ytest, 
  true1pred1boosting=true1pred1boosting, true1pred2boosting=true1pred2boosting, true2pred1boosting=true2pred1boosting, 
  true2pred2boosting=true2pred2boosting, ytesthatboosting=ytesthatboosting, true1pred1boostingvarsel=
    true1pred1boostingvarsel, true1pred2boostingvarsel=true1pred2boostingvarsel, true2pred1boostingvarsel=
    true2pred1boostingvarsel, true2pred2boostingvarsel=true2pred2boostingvarsel, 
  ytesthatboostingvarsel=ytesthatboostingvarsel, true1pred1plslda=true1pred1plslda, 
  true1pred2plslda=true1pred2plslda, true2pred1plslda=true2pred1plslda, 
  true2pred2plslda=true2pred2plslda, ytesthatplslda=ytesthatplslda, true1pred1nsc=true1pred1nsc, 
  true1pred2nsc=true1pred2nsc, true2pred1nsc=true2pred1nsc, true2pred2nsc=true2pred2nsc, 
  ytesthatnsc=ytesthatnsc, true1pred1plsldavarsel=true1pred1plsldavarsel, 
  true1pred2plsldavarsel=true1pred2plsldavarsel, true2pred1plsldavarsel=true2pred1plsldavarsel, 
  true2pred2plsldavarsel=true2pred2plsldavarsel, ytesthatplsldavarsel=ytesthatplsldavarsel, 
  true1pred1rf=true1pred1rf, true1pred2rf=true1pred2rf, true2pred1rf=true2pred1rf, 
  true2pred2rf=true2pred2rf, ytesthatrf=ytesthatrf, 
  true1pred1knnvarsel=true1pred1knnvarsel, true1pred2knnvarsel=true1pred2knnvarsel, true2pred1knnvarsel=true2pred1knnvarsel, 
  true2pred2knnvarsel=true2pred2knnvarsel, ytesthatknnvarsel=ytesthatknnvarsel)
  
return(res)

}
 