library(caret)
library(tidyverse)
library(DMwR)
library(randomForest)
library(FSelector)
# library(randomForestExplainer)
# may need to: options(expressions = 5e5) to avoid stackoverflow for installing package

set.seed(42)

ngramfilegold<-"goldstandard01_anno_ml_synfeat_nstop1"
ngramfiletest<-"wikilist_cleanup_syn"
setwd(dirname(rstudioapi::getSourceEditorContext()$path))
stopwords <- readLines(con = "../data/stopwords.txt",encoding="UTF-8")
oringramme <- read.csv(paste("../data/",ngramfilegold,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
oringrammetest <- read.csv(paste("../data/",ngramfiletest,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
syfeaturenames <- read.csv("../data/syfeatures.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
# syfeaturenames$navalue<-sapply(syfeaturenames$navalue,as.numeric)

deleteStopwords =  function(wl, stopwords = NULL) {
  wl[!(wl %in% stopwords)]
}

# featuresets

covars <- c("CO_LL", "CO_Z", "CO_G", "CO_T", "CO_LOGDICE", "CO_PMI", "CO_MI3", "CO_DEREKO", "CO_SGT", "CO_WIN5_VEC","CO_WIN5_VEC_AUTOSEM")
syvars <- c(syfeaturenames$synames,"NSTOPW")
vars <- c(covars,syvars)

oringramme <- oringramme %>%
  mutate(CO_IDIOM = ifelse(is.na(CO_IDIOM),0,CO_IDIOM)) %>%
  filter(CO_IDIOM < 2)   # just two classes: 0 no idiom, 1 idiom

# Reduce number of classes, treat null values, add NSTOPW, change names for SY features

ngramme <- oringramme %>%
  add_column(NSTOPW = sapply(oringramme$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
  add_column(nstokens = sapply(oringramme$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
  # select(-matches("CO_TOKEN.*"), -tokens) %>%
  select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
  mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
  mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
  rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
  mutate(across(everything(), ~ replace_na(.x, 0))) %>%
  mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM == 1, "idiom", "no_idiom"))) %>% # just two classes: 0 no idiom, 1 idiom
  add_column(GOLD=1) %>%
  select(c("CO_IDIOM","tokens","nstokens","GOLD",all_of(syvars)))

# discard all ngrams with less than 2 non stopwords (no syntagmatic features possible)

ngramme <- ngramme %>%
  filter(NSTOPW > 1)

# remove duplicates after stopword exclusion

ngramme <- ngramme %>% distinct(nstokens,.keep_all=T)

# wiki ngrams

wikingramme <- oringrammetest %>%
  add_column(NSTOPW = sapply(oringrammetest$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
  add_column(nstokens = sapply(oringrammetest$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
  # select(-matches("CO_TOKEN.*"), -tokens) %>%
  select(-c(IDIOM,KERN)) %>% # keep tokens for interpretability
  mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
  mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
  rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
  mutate(across(everything(), ~ replace_na(.x, 0))) %>%
  add_column(CO_IDIOM = as.factor("idiom")) %>%
  add_column(GOLD = 0) %>%
  select(c("CO_IDIOM","tokens","nstokens","GOLD",all_of(syvars)))

# remove all with NSTOPW > 1 (no syntagmatic context available)

wikingramme <- wikingramme %>%
  filter(NSTOPW > 1)

# remove duplicates after stopword exclusion

wikingramme <- wikingramme %>% distinct(nstokens,.keep_all=T)

# find duplicates by lower cased tokens without stopwords

bothngramme <- merge(ngramme[,c("tokens","nstokens","CO_IDIOM")],wikingramme[,c("tokens","nstokens","CO_IDIOM")],by="nstokens")

# 100% agreement ;)

# combine

allngramme <- rbind(ngramme,wikingramme)


# and again remove duplicates

allngramme <- allngramme %>% distinct(nstokens,.keep_all=T)

ngramme <- allngramme %>% filter(GOLD==1)
wikingramme <- allngramme %>% filter(GOLD==0)

# formulae for training and testing rf

fmla <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
fmlaco <- as.formula(paste("CO_IDIOM ~ ", paste(covars, collapse= "+")))
fmlasy <- as.formula(paste("CO_IDIOM ~ ", paste(syvars, collapse= "+")))

# Train/Test split

# Training: 80% Gold Standard, no_idiom + 100% Gold Standard, idiom
# Test: 20% Gold Standard, no_idiom + 100 % Wiki, idiom

noidiomsgold <- ngramme %>% filter(CO_IDIOM=="no_idiom")
idiomsgold <- ngramme %>% filter(CO_IDIOM=="idiom")

trainRows <- sample(nrow(noidiomsgold), nrow(noidiomsgold)*0.8, replace = FALSE)
train <- rbind(noidiomsgold[trainRows,],idiomsgold)
test <- rbind(noidiomsgold[setdiff(1:nrow(noidiomsgold),trainRows),],wikingramme)

cat("Random Forest\n")

rf_classifier = randomForest(fmlasy, train, importance=TRUE)

prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM))

res <- confusionMatrix(prediction_for_table, test$CO_IDIOM, positive= "idiom")
print(res)
collected_results <- bind_cols("rf" = res$byClass)

# Sensitivity is recall of class 1
# Pos Pred Value is precision
varImpPlot(rf_classifier)

cat("Random Forest with cutoff\n")
prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.2, 0.8))
res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom")
collected_results <- bind_cols(collected_results, "rf with cutoff" = res$byClass)
print(res)

cat("With SMOTE resampled training data\n")
smoted.data <- SMOTE(fmlasy, subset(train, select = c("CO_IDIOM", syvars)), perc.over = 1200, perc.under = 100)
rf_classifier = randomForest(fmlasy, smoted.data, importance=TRUE)
prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM))
res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom")
collected_results <- bind_cols(collected_results, "rf with SMOTE" = res$byClass)
print(res)

cat("With SMOTE and cutoff\n")
prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.2, 0.8))
res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom")
collected_results <- bind_cols(collected_results, "rf with SMOTE and cutoff" = res$byClass)
print(res)

collected_results <- collected_results %>%
  round(3) %>%
  add_column(measure = names(res$byClass)) %>%
  column_to_rownames("measure")

View(collected_results)

# Analysing tradeoff between Fscore, Recall, Precision for various cutoffs
# full range from precision almost 100% to recall almost 100%
rf_classifier = randomForest(fmlasy, train, importance=TRUE)
cvalues<-tibble()
for (c in c(seq(from=0.4,to=0.99,by=0.025),0.999)) {
  prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM), cutoff = c(1-c, c))
  conf<-confusionMatrix(prediction_for_table, test$CO_IDIOM, positive = "idiom")
  cvalues <-bind_rows(cvalues, c(cutoff=c, conf$byClass))
}
cvalues %>%
  select(c("cutoff", "Recall", "Precision", "F1", "Specificity", "Balanced Accuracy")) %>%
  pivot_longer(!cutoff, names_to=c("measure")) %>%
  ggplot(aes(cutoff, value, colour=measure)) + geom_line() +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10))


# Using estimates by random forest on entire dataset

library(randomForest)
rf_classifier_full = randomForest(fmlasy, data=allngramme, importance=TRUE)
rf_classifier_full
# class.error is 1 - recall
varImpPlot(rf_classifier_full)

# Feature ranking

# rf features as table

# correlated features seem to split their rankings

rfranks<-importance(rf_classifier_full)[,3:4]

# ttest

idioms<-ngramme %>% filter(CO_IDIOM == "idiom")
nonidioms<-ngramme %>% filter(CO_IDIOM != "idiom")

ttestPvalues<-sapply(syvars,
                     function(sel) t.test(idioms[sel],nonidioms[sel])$p.value)

# information gain
# multiply by 1000 to avoid undersized bins
# features are ranked individually no matter their correlation
igain<-information.gain(fmlasy, data=ngramme%>%mutate_at(syvars, ~ . * 1000),unit="log2")

featureRanks<-cbind(rfranks,igain,ttestPvalues)

#randomForestExplainer::explain_forest(rf_classifier )

# averate estimates and feature ranks over 10 runs

errrate<-0
conf<-matrix(0,2,3)
featureRanks<-matrix(0,4,length(vars))
for (i in 1:10) {
  rfc =randomForest(fmlasy, data=ngramme, importance=TRUE)
  #rfc =randomForest(fmlasy, data=ngramme, importance=TRUE, cutoff=c(0.2, 0.8))
  errrate<-errrate+rfc$err.rate[100,1]
  conf<-conf+rfc$confusion
  featureRanks<-featureRanks+
    cbind(importance(rfc)[,3:4],
          information.gain(fmlasy, data=ngramme%>%mutate_at(syvars, ~ . * 1000),unit="log2"),
          sapply(syvars,
                 function(sel) t.test(idioms[sel],nonidioms[sel])$p.value))
  print(errrate/i)
  conf1<-round(
    rbind(
      cbind(conf[,1:2]/i,(1-conf[,3]/i)*100),
      c(100*diag(conf[,1:2])/colSums(conf[,1:2]),NA),
      c(rowSums(conf[,1:2]/i),NA)),digits=2)
  colnames(conf1)<-c("1","0","rec")
  rownames(conf1)<-c("1","0","prec","sum")
  print(conf1)
}
featureRanks<-featureRanks/10
colnames(featureRanks)<-c("MeanDecreaseAccuracy","MeanDecreaseGini","InformationGain","Ttest")