Pothole Detection With Machine Learning

rm(list=ls())

#set.seed(4)
set.seed(7)


pacman::p_load(tidyverse,janitor,DataExplorer,skimr,ggpubr,viridis,kableExtra,caret,lightgbm,recipes,rsample,yardstick,pROC,xgboost,doParallel,mlr,readxl,stringr,parallel,VIM,GGally,DALEX,MLmetrics)


#Load  variable importance plot
source("Varplot.R")

source("EvaluationMetrics.R")


theme_set(theme_pubclean())

potholedata<- rio::import("/POTHOLEDETECTION/pothole.xlsx",sheet="Sheet6")

potholedata<-potholedata%>%mutate_if(is.character,as_factor)%>%rename(Label_Combined=Feature,Label=Feature2)
  
  


Date=paste(substr(20180709, 1, 4),substr(20180709, 5, 6),substr(20180709, 7, 8),sep = "-")

#Equivalently
Date<-str_c(substr(20180709, 1, 4),"-",substr(20180709, 5, 6),"-",substr(20180709, 7, 8))

extracdate<-function(x){
d<-str_c(str_sub(x, 1, 4),str_sub(x, 5, 6),str_sub(x, 7, 8),sep = "-")  
 return(d) 
}



potholedata$Date<-as.Date(as_vector(purrr::map(potholedata$Date,extracdate)) )    


#sapply(potholedata$Date, extracdate)

 
potholedata%>%head()%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

plot_str(potholedata)

glimpse(potholedata)

table(potholedata$Label)

#check classes distribution
prop.table(table(potholedata$Label))*100

aggr(potholedata , col=c('navyblue','yellow'),
                    numbers=TRUE, sortVars=TRUE,
                    labels=names(potholedata), cex.axis=.7,
                    gap=3, ylab=c("Missing data","Pattern"))

plot_missing(potholedata)

plot_correlation(potholedata,type = "continuous",theme_config = list(legend.position = "bottom", axis.text.x =
  element_text(angle = 90)))

potholedata%>%select(Label_Combined ,Peak.Duration, Vehicle.Speed..mph., Peak.Z.Accel,Peak.Lat.Accel, Peak.Long.Accel,Front.Slip, Right.Slip, Left.Slip, Direction)%>%ggpairs(title = "")+
  theme(legend.position = "top")

pn<-potholedata%>%select(Label_Combined ,Peak.Duration, Vehicle.Speed..mph., Peak.Z.Accel,Peak.Lat.Accel, Peak.Long.Accel,Front.Slip, Right.Slip, Left.Slip, Direction)%>%ggpairs(title = "",mapping = aes(color = Label_Combined ))+
  theme(legend.position = "top")

pn

pn[1,1]

pn[1,2]

pn[1,3]

pn[1,4]

pn[1,5]

pn[1,6]

pn[1,7]

pn[1,8]

pn[1,9]

pn[1,10]

pn[3,2]

pn[6,5]

pn[5,3]

pn[6,4]

pn[6,5]

pn[6,3]

skimmed <-skim_to_wide(potholedata)
skimmed%>%
  kable() %>%
  kable_styling()

#skimr::skim(potholedata)

  mlr::summarizeColumns(potholedata)%>%  
  
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

library(rsample)

pothole_data<-potholedata%>%dplyr::select(-Date,-Label,-File,-Date,-MLT2.Feature)

data_split <- initial_split(pothole_data, strata = "Label_Combined", prop = 0.8)

pothole_train <- training(data_split)
pothole_test  <- testing(data_split)

pothole_recipe2 <- recipe(Label_Combined ~ ., data = pothole_train ) %>%
  #Transform numeric skewed predictors
  step_YeoJohnson(all_numeric()) %>%
  # standardize the data 
  step_center(all_numeric(), -all_outcomes()) %>%
  #scale the data
  step_scale(all_numeric(), -all_outcomes()) %>%
  #step_kpca a specification of a recipe step that will convert numeric data into one or more principal components using a kernel basis expansion.
  #step_kpca(all_numeric(), num=6)%>%
  #step_log(Label, base = 10)
   # Lump factor levels that occur in <= 10% of data as "other"
  step_other(Direction, MLT1.Feature, threshold = 0.1) %>%
  # Create dummy variables for all nominal predictor factor variables except the response
  step_dummy(all_nominal(), -all_outcomes())%>%
  prep(data = pothole_train,retain = TRUE )

# split data into training and test set
test_tbl2  <- bake(pothole_recipe2, newdata = pothole_test)
train_tbl2 <- bake(pothole_recipe2, newdata = pothole_train) 
glimpse(train_tbl2)

library(ggrepel)
# Fit the multinomial logistic model
model_multlog  <- nnet::multinom(Label_Combined ~., data = train_tbl2)

#summary(model_multlog )

tidy(model_multlog)%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

z<-summary(model_multlog)$coefficients/summary(model_multlog)$standard.errors
# 2-tailed z test
 ((1 - pnorm(abs(z), 0, 1)) * 2) %>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

glance(model_multlog) %>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

# Make predictions
pred3 = predict(model_multlog , newdata = test_tbl2, type="class")

data1=data.frame(test_tbl2["Label_Combined"],predicted=pred3)

#check accuracy
yardstick::metrics(data1,truth = Label_Combined, estimate = predicted) 

# Make predictions
predicted.classes <- model_multlog  %>% predict(test_tbl2)
head(predicted.classes)

#Equivalently Model accuracy

mean(predicted.classes == test_tbl2$Label_Combined)

#plot variable importance
feature=factor(rownames(varImp(model_multlog ,scale=T)))
Importance=varImp(model_multlog ,scale=T)$Overall





Varplot(feature,Importance)

ggsave("/POTHOLEDETECTION/variable.png") 

#confusion matrix
table(predicted.classes,test_tbl2$Label_Combined)

model_multlog_pred<-predict(model_multlog , newdata = test_tbl2, type="probs")

#pred <- predict(model_multlog, newdata = test_tbl2, probability = TRUE)

#MLmetrics::MultiLogLoss(y_true = test_tbl2$Label_Combined, y_pred = attr(model_multlog_pred, "probabilities"))

l=caret::confusionMatrix(predicted.classes,test_tbl2$Label_Combined)

overall_classification_accuracy<-sum(diag(l$table))/sum(l$table)

l$overall%>%data.frame()%>%
  rename(`Macro Measure`=".")%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

l$byClass%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

macro_meaure=l$byClass%>%apply(.,2,mean)%>%
  data.frame()%>%
  rename(`Macro Measure`=".")

macro_meaure%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

l$table%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

pred_probs1<-predict(model_multlog , newdata = test_tbl2, type="probs")

#truth: integer vector with truth labels, values range from 0 to n - 1 classes
 # prob_matrix: predicted probs: column 1 => label 0, column 2 => label 1 and so on
  Multiclasslogloss = function(truth, pred_prob_matrix, eps = 1e-15){

     if(is.double(truth)=="FALSE"){
      
     truth=as.numeric(truth)-1 
      
    }
    
    
    if(max(truth) >= ncol(pred_prob_matrix) || min(truth) < 0){
      stop(cat('True labels should range from 0 to', ncol(pred_prob_matrix) - 1, '\n'))
    }
    
   
    
    
    

   pred_prob_matrix[pred_prob_matrix > 1 - eps] = 1 - eps
   pred_prob_matrix[pred_prob_matrix< eps] = eps
   pred_prob_matrix = t(apply(pred_prob_matrix, 1, function(r)r/sum(r)))
   truth_matrix = matrix(0, nrow = nrow(pred_prob_matrix), ncol = ncol(pred_prob_matrix))
   truth_matrix[matrix(c(1:nrow(pred_prob_matrix), truth + 1), ncol = 2)] = 1
   -sum(truth_matrix * log(pred_prob_matrix))/nrow(pred_prob_matrix)
  }
  

  
Multiclasslogloss(test_tbl2$Label_Combined, pred_probs1)

 n = sum(l$table) # number of instances
 nc = nrow(l$table) # number of classes
 diag = diag(l$table) # number of correctly classified instances per class 
 rowsums = apply(l$table, 1, sum) # number of instances per class
 colsums = apply(l$table, 2, sum) # number of predictions per class
 p = rowsums / n # distribution of instances over the actual classes
 q = colsums / n # distribution of instances over the predicted classes
 
precision = diag / colsums 
recall = diag / rowsums 
f1 = 2 * precision * recall / (precision + recall) 
data.frame(precision, recall, f1)

pacman::p_load(furrr)



### faster version using furrr package
oneVsAll =furrr::future_map(1 : nc,
                      function(i){
                        m = c(l$table[i,i],
                              rowsums[i] - l$table[i,i],
                              colsums[i] - l$table[i,i],
                              n-rowsums[i] - colsums[i] + l$table[i,i]);
                        return(matrix(m, nrow = 2, byrow = T))})







#sum matrices element wise
#s=rowSums(c, dims = 2)



#equivalenlty
#apply(c, c(1,2), sum)




## A general-purpose adder:
#this function adds the elements of the matrix
#in the list elementwise
add <- function(x) Reduce("+", x)
sum_all<-add(oneVsAll)

Average_Accuracy <- sum(diag(sum_all)) / sum(sum_all)

Average_Accuracy

Micro_Accuracy<-(diag(sum_all) / apply(sum_all,1, sum))[1];
Micro_Accuracy

fitControl <- trainControl(method="repeatedcv",
 #For repeated k-fold cross-validation only: the number of complete sets of folds to compute                          
                     repeats = 10,
                     savePredictions = "final",
                     classProbs=TRUE,
          #Either the number of folds or number of resampling iterations
                     number = 10,
                     summaryFunction=multiClassSummary,
                     allowParallel=TRUE) 

knn_grid <- expand.grid(
  kmax = 1:9,
  distance = 1:3,
  kernel = c("rectangular", "triangular", "gaussian")
  )



knnFit <- caret::train(Label_Combined ~ ., 
                       data = train_tbl2, 
                       method = "kknn",
                       trControl = fitControl,
                       tuneGrid = knn_grid,
                       summaryFunction=multiClassSummary,
                        ## Specify which metric to optimize
                        metric = "Accuracy")

#Output of kNN fit
knnFit

getTrainPerf(knnFit)

ggplot(knnFit)

saveRDS(knnFit , "/POTHOLEDETECTION/knnFit.rds")
knnFit <- readRDS("/POTHOLEDETECTION/knnFit.rds")

kknnpred1 = predict(knnFit, newdata = test_tbl2, type="prob")
kknnpred = predict(knnFit, newdata = test_tbl2, type="raw")


library(yardstick)

# Compute the summary statistics
d2=data.frame(test_tbl2["Label_Combined"],predicted=kknnpred)
yardstick::metrics(d2,truth = Label_Combined, estimate = predicted) 

yardstick::conf_mat(d2,truth = Label_Combined, estimate = predicted)

l=caret::confusionMatrix(kknnpred,test_tbl2$Label_Combined)



l$overall%>%data.frame()%>%
  rename(`Macro Measure`=".")%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

l$byClass%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

macro_meaure=l$byClass%>%apply(.,2,mean)%>%
  data.frame()%>%
  rename(`Macro Measure`=".")

macro_meaure%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

 Multiclasslogloss(test_tbl2$Label_Combined, kknnpred1)

library(h2o)
h2o.init()

pf=as.h2o(train_tbl2)

pftest=as.h2o(test_tbl2)

# Run GLM 
myX = setdiff(colnames(train_tbl2), "Label_Combined")

h2o_mult<-h2o.glm(y = "Label_Combined", x = myX, training_frame = pf, family = "multinomial",
nfolds = 0, alpha = 0, lambda_search = FALSE)

h2o.performance(h2o_mult)

pftest=as.h2o(test_tbl2)

predh20mult=h2o.predict(h2o_mult,pftest)

predict_model.H2OMultinomialModel <- function(x, newdata, type, ...) {
    
  
# Function performs prediction and returns dataframe with Response
  
# x is h2o model
# newdata is data frame
# type is only setup for data frame

    pred <- h2o.predict(x, as.h2o(newdata))

    # return classification probabilities only
    return(as.data.frame(pred[,-1]))

}

h2ologisticpred=predict_model.H2OMultinomialModel(h2o_mult,pftest)

h2o_multdata<-data.frame(h2o.varimp(h2o_mult))



Varplot(feature=h2o_multdata$names,Importance=h2o_multdata$coefficients)

ggsave("/POTHOLEDETECTION/h2o_multcombined.png") 

# predicted probabilities for each class
head(h2ologisticpred)

#select the indices of the column on rows  of the class with the highest predicted probability as the class label prediction


x=apply(h2ologisticpred, 1, which.max)

prob=apply(h2ologisticpred,1,max)


## Equivalent  way to convert numeric values to factor variables with dplyr
# recode works better, it preserves levels in the factors
y=recode_factor(x, `1` = "B1", `2` = "B2", `3` = "B3",`4` = "B4", `5` = "B5")



probdataframe<-data.frame(x=x,prob=prob,predclass=y)

#accuracy

mean(test_tbl2$Label_Combined==probdataframe$predclass)

l=caret::confusionMatrix(probdataframe$predclass,test_tbl2$Label_Combined)



l$overall%>%data.frame()%>%
  rename(`Macro Measure`=".")%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

l$byClass%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

macro_meaure=l$byClass%>%apply(.,2,mean)%>%
  data.frame()%>%
  rename(`Macro Measure`=".")

macro_meaure%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

Multiclasslogloss(test_tbl2$Label_Combined, as.matrix(predh20mult[,-1]))

xgbmGrid <-  expand.grid(max_depth = c(1, 5, 8,10),
                         nrounds=300,
                        colsample_bytree = c(0.3,0.4,0.5,0.6,0.8,1), 
                        eta = seq.int(from=0.01, to=0.3, length.out=5),
                        min_child_weight = 1:3,
                        subsample=1,
                        
                        #rate_drop=0.1, 
                        #skip_drop=0.1
                        gamma=0.001
                        )

xgbm <-caret::train(Label_Combined ~ ., data = train_tbl2, 
                 method = "xgbTree", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 ## Now specify the exact models 
                 ## to evaluate:
                  #summaryFunction=multiClassSummary,
                 tuneGrid = xgbmGrid)



saveRDS(xgbm , "/POTHOLEDETECTION/xgbm.rds")
xgbm <- readRDS("/POTHOLEDETECTION/xgbm.rds")



getTrainPerf(xgbm)

ggplot(xgbm)

xgbm$bestTune

xgbm_pred = predict(xgbm, newdata = test_tbl2, type="raw")
xgbm_pred1 = predict(xgbm, newdata = test_tbl2, type="prob")


xgbmd=data.frame(test_tbl2["Label_Combined"],predicted=xgbm_pred)
yardstick::metrics(xgbmd,truth = Label_Combined, estimate = predicted) 

l=caret::confusionMatrix(xgbm_pred,test_tbl2$Label_Combined)



l$overall%>%data.frame()%>%
  rename(`Macro Measure`=".")%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

l$byClass%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

macro_meaure=l$byClass%>%apply(.,2,mean)%>%
  data.frame()%>%
  rename(`Macro Measure`=".")

macro_meaure%>%
  kable(escape = F, align = "c") %>%
  kable_styling(c("striped", "condensed"), full_width = F)

 Multiclasslogloss(test_tbl2$Label_Combined, xgbm_pred1)