### bart16eOptimizationViaMultimodelSimulation
### Version 2.0 2017 Aug 26
### Thomas Bartz-Beielstein and Martin Zaefferer
### Note: The file "bart16eAllData.csv" has to be in the 
### current working directory of R.
### 

require(SPOT)
require(DEoptim)

#initializes RNG seed
set.seed(123)
#read data file
df <- read.csv("bart16eAllData.csv")
#select relevant columns from data
df <- df[, c("E", "he","be","Dt","ht","hz","Da","Du")]
#manually remove some outliers
df <- df[ df$E < 100,]
df <- df[ df$E > 80,]
#split into input and output data
x <- as.matrix(df[,c("he","be","Dt","ht","hz","Da","Du")])
y <- as.matrix(df[,"E",drop=F])
#because we will do minimization (and E is maximized):
y <- -y 
#build stacked ensemble
st <- buildEnsembleStack(x, y)  
#generate a target function from the fit
fun <- evaluateModel(st)
#wrapper for optimization algorithm (DE uses vector inputs)
f <- function(x){fun(matrix(x,1))}
#determine search space bounds
lo <- apply(x,2,min) #lower bounds
up <- apply(x,2,max) #upper bounds
#optimize the target function. Note, that this may require some time to compute.
set.seed(123) #initializes RNG seed
optimum <-  DEoptim(fn = f,lower=lo,upper=up,
                      control=list(trace=TRUE,iter=400,reltol=1e-4,steptol=20))
#print results: best solution
optimum$optim$bestmem
#best solution scaled
optimum$optim$bestmem / 160
#fitness of best solution
optimum$optim$bestval
#compare resulting optimum with bounds....
res <- rbind(lo,optimum$optim$bestmem,up)
rownames(res) <- c("lower","optimum","upper")
pdf(file="de1.pdf")
barplot(res
        ,ylim = c(0, 121)
        ,main=""
        ,xlab="parameter"
        ,col = c(gray(.1), gray(.5), gray(.9))
        ,legend = rownames(res)
        ,args.legend=list(x=10, y= 120)
        ,beside = TRUE
        ,xpd = FALSE) 
dev.off()