# This script (1) allows the update of the classification of trends based on the new thresholds
# This script (2) allows the creation of an alluvial plot based on switch in trend classification
# It is based on changes on 2 functions : classifyLinearTrends and plotLinearClassification
# It updates : 
# - the figures obtained showing all species classification trends (points with error bar or barplot)
# - the tables obtained especially the column "trend"
# NB : the script works if the working directory is the root of the indicator routine

rm(list = ls()) 

# Load libraries
library(tidyverse)
library(ggalluvial)

# Result path
result_path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/"

# FUNCTIONS ----
# Load all functions required from the routine
source("R/identifyConvIssue.R")
source("R/formatTrendEstimates.R")
source("R/plotLinearClassification.R")
source("R/fillAbsence.R")

# Updated version of the function 'classifyLinearTrends' 
classifyLinearTrends_bis = function(data, distribution, threshold_inf = 0.964, threshold_sup = 1.036){
  data$trend = NA
  
  # Erase NA values
  dataNA = data[is.na(data$estimate),]
  data = data[!is.na(data$estimate),]
  
  # Turn to exponential estimates and confidence intervals
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = exp(data[,c("estimate", "infIC", "supIC")])
    limit = 1
  }else{
    limit = 0
  }
  
  # Decrease
  data$trend[data$pval < 0.05 & data$estimate < limit] = "Déclin modéré"
  data$trend[data$pval < 0.05 & data$estimate < threshold_inf] = "Déclin modéré à fort"
  data$trend[data$pval < 0.05 & data$supIC < threshold_inf] = "Déclin fort"
  
  # Increase
  data$trend[data$pval < 0.05 & data$estimate > limit] = "Augmentation modérée"
  data$trend[data$pval < 0.05 & data$estimate > threshold_sup] = "Augmentation modérée à forte"
  data$trend[data$pval < 0.05 & data$infIC > threshold_sup] = "Augmentation forte"
  
  # Stability
  data$trend[data$pval > 0.05 & (data$infIC < threshold_inf | data$supIC > threshold_sup)] = "Incertain"
  data$trend[data$pval > 0.05 & data$infIC >= threshold_inf & data$supIC <= threshold_sup] = "Stable"
  
  # Turn back coefficients
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = log(data[,c("estimate", "infIC", "supIC")])
  }
  
  # Bind all dataframes
  data = rbind(data, dataNA)
  
  return(data)
}

# Updated version of the function 'plotLinearClassification' 
plotLinearClassification_bis <- function(data, distribution, type = "barplot", threshold_inf = 0.964, threshold_sup = 1.036,
                                         path = ""){
  
  # Erase NA values
  data = data[!is.na(data$trend),]
  
  # Transform coefficients according to distribution
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = exp(data[,c("estimate", "infIC", "supIC")])
    zeroVal = 1
  }else{
    zeroVal = 0
  }
  
  # Order and colors for trends
  ## Names
  refOrder = c("Strong decline", "Moderate to strong decline", "Moderate decline", 
               "Stable", "Moderate increase", "Moderate to strong increase", 
               "Strong increase", "Uncertain")
  
  ## Colors
  refCols = c(rev(c("#274034", "#4d8067", "#81b29a", "#F2CC8F","#f57e56","#d8400d", "#d61a0d")),  "gray")
  
  ## Point types
  refTypes = c(16, 1, 3, 3, 3, 1, 16, 4)
  
  # Filter for trends that occur
  refCols = refCols[which(!is.na(match(refOrder,unique(data$trend))))]
  refTypes = refTypes[which(!is.na(match(refOrder,unique(data$trend))))]
  refOrder = refOrder[which(!is.na(match(refOrder,unique(data$trend))))]
  
  # Reorder data trends with refOrder
  data$trend = factor(data$trend, levels = refOrder)
  
  # Reorder data species with decreasing growthrate
  data = dplyr::arrange(data, estimate)
  data$species = factor(data$species, data$species)
  
  # Make barplot
  if(type == "barplot"){
    plot = ggplot2::ggplot(data, ggplot2::aes(x = trend, fill = trend)) +
      
      # Make barplot
      ggplot2::geom_bar(stat = "count") +
      
      # Adapt colors
      ggplot2::stat_count(geom = "text", colour = "white", size = 6.5,
                          ggplot2::aes(label = ggplot2::after_stat(count)), 
                          position = ggplot2::position_stack(vjust=0.5)) +
      
      # Change colors
      ggplot2::scale_fill_manual("Classification", values = refCols) +
      
      # Change theme
      ggplot2::theme_bw() + 
      
      # Set a maximum y value
      ggplot2::ylim(c(0,max(table(data$trend)))) + 
      
      # Change axis names
      ggplot2::xlab("") + ggplot2::ylab("Number of species") +
      
      # Format theme
      ggplot2::theme(panel.grid.major.x = ggplot2::element_blank(),
                     panel.grid.minor.y = ggplot2::element_blank(),
                     axis.text.x = ggplot2::element_blank(), 
                     axis.ticks.x = ggplot2::element_blank(),
                     axis.text.y = ggplot2::element_text(size = 18),
                     axis.title.y = ggplot2::element_text(size = 20),
                     legend.text = ggplot2::element_text(size = 12),
                     legend.title = ggplot2::element_text(size = 14),
                     legend.position = "none")
    
  }else if(type == "errorbar"){
    plot = ggplot2::ggplot(data, ggplot2::aes(y = species, x = estimate, col = trend, shape = trend)) + 
      # Add estimates and errorbars
      ggplot2::geom_pointrange(ggplot2::aes(xmin = infIC, xmax = supIC), linewidth = 0.8, fatten = 3) +
      
      # Add vertical line for stable/uncertain effect
      ggplot2::geom_vline(xintercept = zeroVal, linetype = "dashed", col = "gray") +
      
      # Add vertical line for strong negative trend
      ggplot2::geom_vline(xintercept = threshold_inf, linetype = "dashed", col = "#d61a0d") +
      
      # Add vertical line for strong positive effect
      ggplot2::geom_vline(xintercept = threshold_sup, linetype = "dashed", col = "#274034") +
      
      # Change color
      ggplot2::scale_color_manual("Classification", values = refCols) +
      
      # Change point types
      ggplot2::scale_shape_manual("Classification", values = refTypes) +
      
      # Change names of axis
      ggplot2::ylab("Species") + ggplot2::xlab("Growth rates") + 
      
      # Apply classic theme
      ggplot2::theme_bw() +
      
      # Format theme
      ggplot2::theme(axis.ticks.y = ggplot2::element_blank(), 
                     axis.text.y = ggplot2::element_blank(),
                     axis.title = ggplot2::element_text(size = 20),
                     legend.title = ggplot2::element_text(size = 14),
                     legend.text = ggplot2::element_text(size = 12),
                     axis.text.x = ggplot2::element_text(size = 18),
                     legend.position = "none")
    
  }
  
  # Save result
  ggplot2::ggsave(plot = plot, filename = paste0(type, ".png"), device = "png", path = path, 
                  width = 20, height = 15, units = "cm")
  
  return(plot)
}

# 1. UPDATE CLASSIFICATION FIGURES ----

# Medium & short-term trend
dataLT = data.frame(data.table::fread(file = "data/STOC_2023_extra/countingData.csv", encoding = "UTF-8", drop = "V1"))
dataMT = data.frame(data.table::fread(file = "data/STOC_2023/countingData.csv", encoding = "UTF-8", drop = "V1"))
speciesList = unique(dataMT$species)
speciesList = speciesList[!(speciesList %in% c("AYTFUL", "HIMHIM", "RALAQU", "TETRAX", "TRITOT", "CLAGLA"))]
dataMT$ID = paste0(dataMT$year, dataMT$site)
dataLT$ID = paste0(dataLT$year, dataLT$site)
dataMT = fillAbsence(dataMT[,c("ID","species", "year", "site", "abondance")], interestVar = "abondance", speciesList, "once")
dataLT = fillAbsence(dataLT[,c("ID","species", "year", "site", "abondance")], interestVar = "abondance", speciesList, "once")


dataST = dataMT[dataMT$year >= 2013,]

dataLongTerm = formatTrendEstimates(dataLT, speciesList, "STOC_2023_extra", "longTermTrend")
dataMediumTerm = formatTrendEstimates(dataMT, speciesList, "STOC_2023", "longTermTrend")
dataShortTerm = formatTrendEstimates(dataST, speciesList, "STOC_2023", "shortTermTrend")



medClass = classifyLinearTrends_bis(data = dataMediumTerm, distribution = "nbinom2", threshold_inf = 0.961, threshold_sup = 1.030)
shortClass = classifyLinearTrends_bis(data = dataShortTerm, distribution = "nbinom2", threshold_inf = 0.961, threshold_sup = 1.030)
longClass = classifyLinearTrends_bis(data = dataLongTerm, distribution = "nbinom2", threshold_inf = 0.961, threshold_sup = 1.030)

medClass$trend = factor(medClass$trend, levels = c("Augmentation forte","Augmentation modérée à forte", "Augmentation modérée", "Stable",
                                                   "Déclin modéré", "Déclin modéré à fort", "Déclin fort", "Incertain"),
                        labels = c("Strong increase", "Moderate to strong increase", "Moderate increase",
                                   "Stable", "Moderate decline", "Moderate to strong decline", "Strong decline", 
                                   "Uncertain"))
plotLinearClassification_bis(medClass, "nbinom2", type = "barplot", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/")
plotLinearClassification_bis(medClass, "nbinom2", type = "errorbar", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/")

shortClass$trend = factor(shortClass$trend, levels = c("Augmentation forte","Augmentation modérée à forte", "Augmentation modérée", "Stable",
                                                       "Déclin modéré", "Déclin modéré à fort", "Déclin fort", "Incertain"),
                        labels = c("Strong increase", "Moderate to strong increase", "Moderate increase",
                                   "Stable", "Moderate decline", "Moderate to strong decline", "Strong decline", 
                                   "Uncertain"))
plotLinearClassification_bis(shortClass, "nbinom2", type = "barplot", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/")
plotLinearClassification_bis(shortClass, "nbinom2", type = "errorbar", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/")

# 2. UPDATE CLASSIFICATION TABLES ----
# Important tables with old classification
dataTrend_MT = read.csv("outputs/STOC_2023/tables/trends/longTermTrends.csv", row.names = 1)
dataTrend_ST = read.csv("outputs/STOC_2023/tables/trends/shortTermTrends.csv", row.names = 1)

# Erase columns that need updating
dataTrend_MT = dataTrend_MT[,-c(16:18)] 
dataTrend_ST = dataTrend_ST[,-c(16:18)] 

# Extract only new classification
medClass = medClass[,c("species", "trend")]
shortClass = shortClass[,c("species", "trend")]

# Update the classification column and save the resulting data.Frame
dataTrend_MT = merge(dataTrend_MT, medClass, by = "species")
write.csv(dataTrend_MT, paste0(result_path, "longTermTrends_V2.csv"))

dataTrend_ST = merge(dataTrend_ST, shortClass, by = "species")
write.csv(dataTrend_ST, paste0(result_path,"shortTermTrends_V2.csv"))

# 3. MAKE ALLUVIAL PLOT ----

# Adapt colnames for the alluvial plot
colnames(medClass)[2] = "trend_MT"
colnames(shortClass)[2] = "trend_ST"

# Merge both classifications
dataClass = merge(medClass, shortClass, by = "species")

# Count each pair of class
dataAlluvial = group_by(dataClass, trend_MT, trend_ST) %>%
  summarise(count = n())

# Erase all pairs that correspond to no estimated trends
dataAlluvial = dataAlluvial[!is.na(dataAlluvial$trend_MT) & !is.na(dataAlluvial$trend_ST),]
# Erase all pairs that correspond to both uncertain trends
dataAlluvial = dataAlluvial[!(dataAlluvial$trend_MT == "Uncertain" & dataAlluvial$trend_ST == "Uncertain"),]

# Specify labels and colors
lab = c("Strong increase", "Moderate to strong increase", "Moderate increase", 
        "Stable", "Moderate decline", "Moderate to strong decline", "Strong decline", 
        "Uncertain") 

colTrends = c("#274034", "#4d8067", "#81b29a", "#F2CC8F","#f57e56","#d8400d", "#d61a0d",  "gray")

# Make the alluvial
ggplot(dataAlluvial,
       aes(y = count, axis1 = trend_MT, axis2 = trend_ST)) +
  geom_alluvium(data = dataAlluvial, aes(fill = trend_MT), alpha = 0.9) + 
  geom_stratum(fill = rev(rep(colTrends,2)),
               width = 1/6, col = "gray10") +
  scale_fill_manual("Classification", values = colTrends) +
  theme_bw() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)),
            reverse = TRUE, size = 2.5, col = "black") +
  ylab("Number of species") +
  scale_x_discrete(limits = c("2001 - 2023\nperiod", "2013 - 2023\nperiod"), expand = c(.05, .05)) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
        axis.title.y = element_text(size = 16),
        axis.text.y = element_text(size = 14))

# Save the plot
ggsave("alluvial_plot_zoom_STOC_V2.png", last_plot(), path = result_path, device = "png", 
       width = 22, height = 18, units = "cm")
