rm(list=ls(all=TRUE))
isBatch = FALSE
if(isBatch) { options(echo  =  FALSE); } else { options(echo  =  TRUE); }

library("openxlsx")
library("e1071")
library("effsize")
library("nnet")
library("neuralnet")
library("randomForest")

# ----------------------- constants and parameters ---------------------------
theSeed=9295

VERBOSE=FALSE  # set this variable to TRUE to get messages concerning elaboration progress
OUT_TO_FILE=TRUE  # set this variable to TRUE to save results to files

EXCLUDE_OUTLIER=FALSE

# for K-times K-fold analysis:
Ktimes=10
Kfold=10

MODEL_TYPE="NN"   #   "SVR"   #   "RF"   #  
if(MODEL_TYPE=="SVR"){
  # selection of SRV kernel type
  krnlType="linear"  #  "radial"   #   "polynomial"   #  
}
if(MODEL_TYPE=="NN"){
  NNsize=2
}

NUM_TUNING_ITERATIONS=10 # how many time the hyperparameter tuning function is applied

dependentVariableName="Time"
measuresToBeUsed=c("CoCo", "McCC") # independent variables (features) to be used for model building

# ------------------ functions ----------------------------------------

#  compose a formula, given the dependent variable name and the array of indep. var. names
makeFmla<-function(depVarName, listOfVars, intercept=TRUE){
  for(j in 1:length(listOfVars)){
    if(j==1){
      fmla=paste(depVarName, " ~ ", listOfVars[j], sep="")
    } else {
      fmla=paste(fmla, " + ", listOfVars[j], sep="")
    }
  }
  if(!intercept){
    fmla=paste(fmla, " + 0", sep="")
  }
  return(as.formula(fmla))
}

# perform K-times K-fold with any type of model
doNfoldXval<-function(mdlType, fmla, theData, dependentVariableName, verbose=FALSE){
  thePred=c()
  theActual=c()
  n=nrow(theData)
  for(jj in 1:Ktimes){
#    cat("Iteration" , jj, "\n")
    nn=floor(nrow(theData)/Kfold)
    permut=sample(1:n, n)
    for(ii in 1:Kfold){
      if(ii==Kfold){
        theRange=((ii-1)*nn+1):n
      } else {
        theRange=((ii-1)*nn+1):(ii*nn)
      }
      testRows=permut[theRange]
      testSet=theData[testRows,]
      trainingSet=theData[-testRows,]
      if(mdlType=="SVR"){
        tuneResult=doTheTuning("SVR", fmla, trainingSet, nn=NUM_TUNING_ITERATIONS, verbose=verbose)
        if(!is.null(tuneResult)){
          svmfit=tuneResult$best.model
          if(!is.null(svmfit)){
            pred=predict(svmfit, newdata=testSet)
            actuals=testSet[,dependentVariableName]
            thePred=c(thePred, pred)
            theActual=c(theActual, actuals)
          } else {
            pred=c()
          }
        }
      }
      if(mdlType=="RF"){
        tuneResult=doTheTuning("RF", fmla, trainingSet, nn=NUM_TUNING_ITERATIONS, verbose=verbose)
        if(!is.null(tuneResult)){
          rffit=tuneResult$best.model
          if(!is.null(rffit)){
            pred=predict(rffit, newdata=testSet)
            actuals=testSet[,dependentVariableName]
            thePred=c(thePred, pred)
            theActual=c(theActual, actuals)
          }
        }
      }
      if(mdlType=="NN"){
        tuneResult=doTheTuning("NN", fmla, trainingSet, nn=NUM_TUNING_ITERATIONS, verbose=verbose)
        if(!is.null(tuneResult)){
          bestMdl=nnet(fmla, trainingSet, size=tuneResult$best.parameters$size,
                     linout=TRUE, rang=0.1, maxit = 100, decay=0, trace=TRUE)
          if(!is.null(bestMdl)){
            pred=predict(bestMdl, newdata=testSet)
            actuals=conv_x(testSet[,dependentVariableName])
            pred=conv_x(pred)
            thePred=c(thePred, pred)
            theActual=c(theActual, actuals)
          }
        }
      }
      vMARtest<<-c(vMARtest, mean(abs(actuals-pred)))
    }
  }
  return(list(theActual, thePred))
}

# gives the tuned SVR model for the formula and data, performing nn iterations
doTheTuning<-function(mdlType, fmla, theData,  nn=8, verbose=TRUE){
  gammas = 2^(-1:3)
  costs = 2^(-5:5)
  epsilons = c(0.01, 0.005, 0.001)
  bestMAR=10000000
  bestIter=0
  theBest=NULL
  for(i in 1:nn){
    if(mdlType=="NN"){
      tuneResult <- tune.nnet(fmla, data = theData, size = NNsize, 
                linout=TRUE, rang=0.1, maxit = 100, 
                tunecontrol = tune.control(cross = 5, sampling = "cross"))
    } else {
      if(mdlType=="SVR"){
        if(krnlType=="linear"){
          tuneResult <- tune.svm(fmla, data = theData, kernel = krnlType,
                                 cost = costs, epsilon = epsilons,
                                 type = "eps-regression",
                                 tunecontrol = tune.control(sampling = "cross", cross = 5),
                                 scale=TRUE)
        } else {
          tuneResult <- tune.svm(fmla, data = theData, kernel = krnlType,
                                 cost = costs, epsilon = epsilons,
                                 gamma=gammas,
                                 type = "eps-regression",
                                 tunecontrol = tune.control(sampling = "cross", cross = 5),
                                 scale=TRUE)
        }
      } else if(mdlType=="RF"){
        tuneResult <- tune.randomForest(fmla, data = theData,
                                        tunecontrol = tune.control(cross = 5, sampling = "cross"))
      } else {
        stop("do the tuning: illegal argument")
      }
    }
    mdl <- tuneResult$best.model  
    if(verbose && mdlType=="SVR"){
      cat("Found mdl with cost=", mdl$cost, " epsilon=", mdl$epsilon, " gamma=", mdl$gamma, "\n", sep="")
    }
    pred=predict(mdl, newdata=theData)
    if(length(pred)>0){
      if(mdlType=="SVR" || mdlType=="RF"){
        abserr=abs(theData[,dependentVariableName]-pred)
      } else {
        if(mdlType=="NN"){
          abserr=abs((theData[,dependentVariableName]*maxminY+minY)-(pred*maxminY+minY))
        } else {
          stop("do the tuning: illegal argument")
        }
      }
      MAR=mean(abserr)
      if(MAR<bestMAR){
        bestMAR=MAR
        theBest=tuneResult
      }
    }
    if(verbose){
      cat("At end of iteration ", i, ", bestMAR=", bestMAR, "(previous was MAR=", MAR, ")\n")
    }
  }
  if(verbose){
    if(is.null(theBest)){
      cat(" DO THE TUNING FAILED !\n")
    } else {
      cat("End of tuning: best MAR=", bestMAR, ", obtained at iteration ", bestIter, "\n", sep="")
    }
  }
  vMARtraining<<-c(vMARtraining, bestMAR)
  return(theBest)
}

# normalize the dataset
min_max_norm <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}
# restore from normalization
conv_x <- function(x) {
  x * maxminY + minY
}

# save actuals and corresponding predictions to file
saveActPredToFile<-function(actuals, estimates, vMeasures, mtype){
  tdf=data.frame(actuals, estimates)
  theNames=""
  for(jj in 1:length(vMeasures)){
    theNames=paste(theNames, "_", vMeasures[jj], sep="")
  }
  fname=paste("./results_", mtype, "_", theNames, ".csv", sep="")
  write.csv(tdf, file=fname)
}

# ---------- end of functions ---------------------------------------------

# ----------- data preparation ----------------------------------------------

filename=paste("./measures.csv", sep="") # change if needed
theData=read.csv(filename)
if(EXCLUDE_OUTLIER){
  theData=subset(theData, theData$Time<=60)
}

par(mar=c(4.4, 4.0, 0.8, 0.4))  # bottom, left, top and right margins

# --------------- analysis starts here ---------------------
set.seed(theSeed)

vMARtraining=c()
vMARtest=c()

if(MODEL_TYPE=="NN"){
  minY=min(theData[,dependentVariableName])
  maxminY=max(theData[,dependentVariableName])-min(theData[,dependentVariableName])
  theData<- theData[, c(dependentVariableName, measuresToBeUsed)]
  theData <- as.data.frame(lapply(theData, min_max_norm))
}

fmla=makeFmla(dependentVariableName, measuresToBeUsed, intercept = FALSE)

ttt=doNfoldXval(MODEL_TYPE, fmla, theData, dependentVariableName, verbose=VERBOSE)
tfmla=as.character(fmla)
myTitle=paste("Results with ", tfmla[2], " ~ ", tfmla[3], sep="")
if(OUT_TO_FILE){
  saveActPredToFile(ttt[[1]], ttt[[2]], measuresToBeUsed, MODEL_TYPE)
}

