# CODE SVM, RVM AND YOUD EMPIRICAL EQUATION

# Load Packages
library('e1071')    # SUPPORT VECTOR MACHINE
library('boot')     # FUNCTION corr
library('qualV')    # FOR MEASURES OF ERROR RMSE & MAE


#####################################
# Functions #########################
#####################################
# Pearson's product-moment correlation coefficient
CoeffCorr <- function(O,P){
  ComB <- cbind(O, P)
  r <- corr(ComB, w=rep(1, nrow(ComB))/nrow(ComB))
  r <- round(r*100, 2)
  return(r)
}
# (Nash & Sutcliffe , 1970) Coefficient of Efficiency
CoeffEffic <- function(O,P){
  E <- ((1 - ((sum((O - P)^2))/(sum((O - mean(O))^2))))*100)
  E <- round(E, 2)
  return(E)
}

############################################################################
### READING THE WHOLE DATA AND IDENTIFYING THE FREE FACE AND SLOPE INDEX ###
############################################################################

dat <- read.table("researchLateralFreeFaceData.csv",sep=",",header=TRUE)# Change the file name here
N <- nrow(dat)
SlopeIndex <- which(dat$W==1)                    # Index of the Slope data 
FreeFaceIndex <-  setdiff(x = 1:N, SlopeIndex)   # Index of free face data
FreeFaceData <- matrix(nrow = length(FreeFaceIndex), ncol = ncol(dat), data = 0)
FreeFaceData <- dat[FreeFaceIndex,]              # Free face data
SlopeData <- matrix(nrow = length(SlopeIndex), ncol = ncol(dat), data = 0)  
SlopeData <- dat[SlopeIndex,]                    # Slope Data

   

############################################################################
### PREPARING THE FREE FACE DATA FOR N FOLD CROSS VALIDATION ###############
############################################################################

n = 5 # YOU CAN CHANGE THE NUMBER OF CROSS VALIDATION SETS

# FREE FACE ANALYSIS

N = nrow(FreeFaceData)
nrows <- ceiling(N/n)
set.seed(1)
rf <- sample(N, size = N)
FreeFaceNewData <- as.matrix(FreeFaceData[rf, ])
FreeFaceStackData <- matrix(nrow = N, ncol = 11, data = 0)
 
  FFDataInd <- 1:219
  FreeFaceNewData <- cbind(FreeFaceNewData, FFDataInd)
  SVdata <- list()

########################### SVM FREE FACE ##################################

## RUN THE MODEL
   
   ypredictResultSVM <- matrix(nrow = N, ncol = 2, data = 0)
   ypredictCrossvalResultSVM <- matrix(nrow = (n+1), ncol = 4, data = 0)
   end = 0
     for(i in 1:n){
     start <- end +1
     end <- (N - (nrows*(n-i)))
     TestData <- FreeFaceNewData[(start:end),]
     TrainDataIndex <-  setdiff(x = 1:N, (start:end)) 
     TrainData <- FreeFaceNewData[TrainDataIndex,]
     # Training Data
     xtrain <- TrainData[ , c(1,2,3,5,6,7)]
     ytrain <- log(as.matrix(TrainData[ , 8]))
     #Testing Data
     xtest <- TestData[ , c(1,2,3,5,6,7)]
     yactual <- log(as.matrix(TestData[ ,8]))
     # SVM Parameters
     g <- .136         # Gamma
     c <-  9           # Cost
     eps <- 0.162      # Epsilon

     # SVM Model
     model <- svm(xtrain, ytrain, scale = TRUE, type = "eps-regression",
             kernel = "radial", cross = 10, cost = c, epsilon = eps,
             gamma = g)
     # Model Predictions
     ypredictSVMtest <- predict(model,xtest)
     count=1
       for (j in start:end){
       ypredictResultSVM[j,1] <- exp(yactual[count])
       ypredictResultSVM[j,2] <- exp(ypredictSVMtest[count])
       FreeFaceStackData[j,1:6] <- xtest[count,]
       FreeFaceStackData[j,7] <- exp(ypredictSVMtest[count])
       FreeFaceStackData[j,10] <- exp(yactual[count])
       count = count+1
       }
      SVdata[[i]] <- list(TrainData, TestData, SVindex=model$index, Ypredict=exp(ypredictSVMtest)) # Data saved for SV analysis

    # Measure of performance for the particular cross validation set
    ypredictCrossvalResultSVM[i,1] <- round(RMSE(exp(yactual),exp(ypredictSVMtest)),2)
    ypredictCrossvalResultSVM[i,2] <- CoeffCorr(exp(yactual),exp(ypredictSVMtest))
    ypredictCrossvalResultSVM[i,3] <- CoeffEffic(exp(yactual),exp(ypredictSVMtest))
    ypredictCrossvalResultSVM[i,4] <- nrow(model$SV)
    }
  # Average Measure of performance for all Runs
  ypredictCrossvalResultSVM[(i+1),1] <- round(RMSE(ypredictResultSVM[,1]
                                             ,ypredictResultSVM[,2]),2)
  ypredictCrossvalResultSVM[(i+1),2] <- CoeffCorr(ypredictResultSVM[,1]
                                             ,ypredictResultSVM[,2])
  ypredictCrossvalResultSVM[(i+1),3] <- CoeffEffic(ypredictResultSVM[,1]
                                             ,ypredictResultSVM[,2])
  ypredictCrossvalResultSVM

####### End #################