##This script demonstrates the use of stacking based ensembles
##based on a gas sensor data set.
##Please note, that some packages are required to run this script.
##Install CEGO (and dependencies) from CRAN with:
#install.packages("CEGO")
##SPOT2 can be installed from the provided .tar.gz file,
##which should be part of the supplemental files.
##The simplest way to do so is by using RStudio
##Go to: Tools -> Install Packages
##then select "Install from package archive file"
##the select the provided .tar.gz file
##Note: this may require that Rtools are installed
##if you use windows. See the CRAN homepage for details.

## load packages into workspace
require(SPOT2)
require(CEGO)

##define an adequate MLE optimizer:
optimizer <- function(x,fun,lower,upper,control,...){
  CEGO::optimInterface(x,fun,lower,upper,control=list(
    method="NLOPT_GN_DIRECT_L",funEvals=400,reltol=1e-8),...)
}

## select which set of sensors to use
sensor = 1
#sensor = 2

if(sensor==1){
  train <- dataGasSensor[dataGasSensor$Sensor==1,1:8]
  test <- dataGasSensor[dataGasSensor$Sensor==3,1:8]
}else{
  train <- dataGasSensor[dataGasSensor$Sensor==2,1:8]
  test <- dataGasSensor[dataGasSensor$Sensor==4,1:8]
}

## main experiment loop over ten repeats:
error <- NULL
for(i in 1:10){
  ## set random number generator seed (reproducibility)
  set.seed(i)
  
  ## build and evaluate Kriging model
  fitK <- buildKriging(data.matrix(
                        train[,c("X1","X2","X3","X4","X5","X6","X7")]),
                      data.matrix(train$Y),
                      control=list(algTheta=optimizer,useLambda=TRUE,
                                   reinterpolate=FALSE))
  predK <- predict(fitK,
                   data.matrix(test[,c("X1","X2","X3","X4","X5","X6","X7")]))$y
  
  ## build and evaluate linear model
  fitLM <- buildLM(data.matrix(
                      train[,c("X1","X2","X3","X4","X5","X6","X7")]),
                  data.matrix(train$Y),
                  control=list(useStep=F,formula="y~X1+X2+X3+X4+X5+X6+X7"))
  predLM <- predict(fitLM,
                    data.matrix(test[,c("X1","X2","X3","X4","X5","X6","X7")]))$y
  
  ## build and evaluate random forest model
  fitRF <- buildRandomForest(
                  data.matrix(train[,c("X1","X2","X3","X4","X5","X6","X7")]),
                  data.matrix(train$Y))
  predRF <-predict(fitRF,
                   data.matrix(test[,c("X1","X2","X3","X4","X5","X6","X7")]))$y
  
  ## build and evaluate ensemble (stacking)
  fitE <- buildEnsembleStack(
                  data.matrix(train[,c("X1","X2","X3","X4","X5","X6","X7")]),
                  data.matrix(train$Y),
                  control=list(modelL1Control=list(formula="y~V1+V2+V3"),
                     modelL0=c(buildLM,buildRandomForest,buildKriging),
                     modelL0Control=list(list(useStep=F,
                                              formula="y~X1+X2+X3+X4+X5+X6+X7"),
                                         list(),
                                         list(algTheta=optimizer,
                                              useLambda=TRUE,
                                              reinterpolate=FALSE)))) 
  predE <- predict(fitE,
                   data.matrix(test[,c("X1","X2","X3","X4","X5","X6","X7")]))$y
  
  ## compute mean square error for each model
  error <- rbind(error, data.frame(
  Kriging=mean(abs(predK - test$Y)^2),
  LinearModel=mean(abs(predLM - test$Y)^2),
  RandomForest=mean(abs(predRF - test$Y)^2),
  Ensemble=mean(abs(predE - test$Y)^2)))
}
## visualize results
boxplot(error)

## compute the respective mean and standard deviations
## which are presented in the respective table in the
## article
s<-round(apply(error,2,sd),digits=3)
y<-round(apply(error,2,mean),digits=3)

## print
s
y

## Note, that results may vary (in comparison to depicted table)
## due to numeric issues (e.g., precision, random number generation)