###################################################################
# Code is used to run machine learning regression and prediction for snow depth/Thaw depth modeling and prediction

# SVM: Support Vector Machine regression
# RF: Random Forest regression
# ANN: Artificial Neural Network regression
# kNN: k Nearest Neighbor regression
# LM: linear regression
# Coded by: Caiyun Zhang, czhang3@fau.edu
# Florida Atlantic University
# Coded date: 4/5/2025
# updated: 6/12/2026 to training inputs of WV-2, Lidar CHM, and DEM using field measurements

# load in required libraries
library(randomForest)
library(caret)
library(e1071)
library(nnet)

# set up working directory where the input data and output will be saved
setwd("I:/Research/Proposals_Project/USACE_2018-2019/Phase_2/FAU_Invoice_Report/Reports/TutorialDataset/Input")
# Import the training samples saved in csv file. 
# The first column is y (dependent variable), and other columns are x independent variables (predictors)
###################################################################

inTrainingFile<-"Training_WV.csv" # input training data with first column as the target
inForPredictionFile<-"InputForPrediction.csv" # input file for final map prediction

k<-5 # set up of k-fold

outMetrix<-"Out/StatisticsMetrix_.csv" # output of model performance
outSample<-"Out/samplePrediction.csv" # output of reference samples from models of k-fold
outPrediction<-"Out/FinalPrediction.csv" # final prediction for mapping to be imported into ArcGIS

set.seed(10)
samples <- read.csv(inTrainingFile)
Data<-read.csv(inForPredictionFile) 
y<-samples[,1] # first column is y
x<-samples[,-1] # other column is x, drop first column

Actual<-y

trainctrl <- trainControl(method = "cv", number = k, allowParallel = TRUE, verboseIter = TRUE, savePredictions = "final")

SVM_model <- train(x,y,method = "svmRadial", preProcess = c("center", "scale"), trControl = trainctrl,trace=F,tuneLength = 10)
print("SVM model is toned")
 
RF_model<- train(x,y,method = "rf",preProcess = c("center", "scale"),ntree = 500,trControl = trainctrl,trace=F,tuneLength = 10)
print("RF model is tuned")
 
# Neural network needs to be specifically tuned
nnGrid <-  expand.grid(size = seq(1, 5, 1), decay = seq(0.1,0.5,0.05))
  # linout: TRUE for regression, FALSE for classification
ANN_model<- train(x, y, method = "nnet", preProcess = c("center", "scale"),trControl = trainctrl,tuneGrid = nnGrid,tuneLength = 5,trace=F,linout=TRUE) 
print("ANN model is tuned")
  
  # get the prediction in each k-fold modeling and order them based on rowIndex
p_SVM<-SVM_model$pred[order(SVM_model$pred$rowIndex),]
p_RF<-RF_model$pred[order(RF_model$pred$rowIndex),]
p_ANN<-ANN_model$pred[order(ANN_model$pred$rowIndex),]
 
SVM<-p_SVM$pred
RF<-p_RF$pred
ANN<-p_ANN$pred
EA<-(SVM+RF+ANN)/3.0  # simply get the mean of three ML algorithms outputs for ensemble prediction

# get the statistics metrics, RMSE, R2, and MAE of the model
res_SVM <- postResample(y, SVM) 
res_RF <- postResample(y, RF)
res_ANN <- postResample(y, ANN)
res_EA <- postResample(y, EA)
  
pred_Sample<-cbind(Actual,RF,SVM,ANN,EA)
res<-cbind(res_RF,res_SVM,res_ANN,res_EA)
  
# write out to files
write.csv(round(res,digits = 2), outMetrix)
write.csv(round(pred_Sample,digits=2),outSample)

# apply tuned model for prediction of xx inputs
# need to read in xx data set
pred_ANN<-predict(ANN_model, Data)
pred_SVM<-predict(SVM_model, Data)
pred_RF<-predict(RF_model, Data)
pANNOut<-round(pred_ANN,digits=2)
pSVMOut<-round(pred_SVM,digits=2)
pRFOut<-round(pred_RF,digits=2)

pred_EA<-(pred_ANN+pred_SVM+pred_RF)/3
nRowTotal<-nrow(Data) # get the total number of objects/rows of the input data set
Id<-0:(nRowTotal-1)
EA_Final<-as.numeric(pred_EA)
final_pred <- cbind(Id,pred_RF,pred_SVM,pred_ANN,EA_Final)

write.csv(final_pred,outPrediction , row.names = FALSE)

print ("Game over")
###################################################################






