# This script allows the creation of trend figures for the 4 birds species retained in the paper
# These figures are slightly modified compared to the ones obtained in the routine
# So that they are more readable (e.g, increased writing size)
# NB : the script work if the working directory is the root of the indicator routine

rm(list = ls()) 

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

# Load libraries
library(tidyverse)
library(ggalluvial)

# FUNCTIONS ----

# Load all functions required from the routine
source("R/formatTrendEstimates.R")
source("R/formatYearlyVariations.R")
source("R/identifyConvIssue.R")
source("R/makeSlopesFromVar.R")
source("R/supportFunctions.R")
source("R/fillAbsence.R")

# Updated version of the function 'percEvolutionToText' 
percEvolutionToText_bis <- function(data, distribution, language){
  
  # Extract percentage of change (and confidence interval)
  dataChange = apply(data[,c("estimate", "infIC", "supIC")], 2, function(x){
    
    # Calculate value for first year
    start = x * 0
    
    # Calculate Value for relative last year
    stop = x * (data$maxYear - data$minYear)
    
    # Turn to exponential if counting or occurrence data
    if(distribution != "gaussian"){
      start = exp(start)
      stop = exp(stop)
    }
    
    # Measure the percentage of difference
    perc = round(100 * (stop- start), 1)
    
    # Return the percentage of difference
    return(perc)
  })
  
  # Turn to text
  changeText = paste0(ifelse(dataChange["estimate"] > 0, "+", ""),
                      dataChange["estimate"], "% (",
                      
                      ifelse(dataChange["infIC"] > 0, "+", ""),
                      dataChange["infIC"], "% ; ",
                      
                      ifelse(dataChange["supIC"] > 0, "+", ""),
                      dataChange["supIC"], "%)")
  
  # Add a 
  if(language == "french"){ 
    preText = "Tendance linéaire "}
  else{preText = "Linear trend over "}
  
  changeText = paste0(preText, data$min, " - ", data$max, " : ", changeText)
  
  return(changeText)
}

# Updated version of the function 'plotLinearTrends' 
plotLinearTrends_bis <- function(speciesList, dataLongTerm, dataYearlyVariations, dataShortTerm = NULL,
                                 language = "french", dataName = NULL, distribution, plotGamm = TRUE, plotST = TRUE,
                                 uncertainty = TRUE, N = 100, weight = TRUE, save = TRUE, path = ""){
  
  plots <- sapply(speciesList, function(sp){ 
    # Filter tables for considered species
    dataLT_sp = dataLongTerm[dataLongTerm$species == sp,]
    dataVar_sp = dataYearlyVariations[dataYearlyVariations$species == sp,]
    
    if(nrow(dataLT_sp)==0|nrow(dataLT_sp)==0) {
      message('No plot for species ', sp, ': no model was found for either ',
              'long-term trend or yearly variations\n')
      plot = NULL
      
    }else if(any(is.na(c(dataLT_sp$estimate, dataVar_sp$estimate)))){
      message('No plot for species ', sp, ' due to no-convergence of either',
              " long-term trend or yearly variations\n")
      plot = NULL
      
    }else{
      # Extract values of the temporal series
      yearValues = dataLT_sp$minYear:dataLT_sp$maxYear
      
      ###################
      # GAMM VARIATIONS #
      ###################
      pathToGamm = here::here("outputs", repo, "models", "gammVariations", paste0(sp, ".rdata"))
      
      if(plotGamm & file.exists(pathToGamm)){
        # Load the corresponding file
        load(here::here("outputs", repo, "models", "gammVariations", paste0(sp, ".rdata")))
        
        if(is.null(gammVariations$warnings) & is.null(gammVariations$error)){
          # Extract values of gam
          valuesToPlotGAM = plot(gammVariations$value$gam, se = TRUE, n = 50 * length(yearValues))
          
          # Make data.frame
          dataGammVariations_sp = data.frame(year = valuesToPlotGAM[[1]]$x, 
                                             estimate = valuesToPlotGAM[[1]]$fit, 
                                             infIC = valuesToPlotGAM[[1]]$fit - 1.96 * valuesToPlotGAM[[1]]$se,
                                             supIC = valuesToPlotGAM[[1]]$fit + 1.96 * valuesToPlotGAM[[1]]$se,
                                             group = "GAM",
                                             rep = 1)
          
          # Centrer sur 0
          dataGammVariations_sp[,c("estimate", "infIC", "supIC")] = dataGammVariations_sp[,c("estimate", "infIC", "supIC")] - mean(dataGammVariations_sp$estimate)
          
          # Passer à l'exponentielle si données d'occurrence ou de comptage
          if(distribution != "gaussien"){
            dataGammVariations_sp[,c("estimate", "infIC", "supIC")] = exp(dataGammVariations_sp[,c("estimate", "infIC", "supIC")])
          }
        }else{
          plotGamm = FALSE
        }
        
      }else{
        plotGamm = FALSE
      }
      
      ####################
      # SHORT-TERM TREND #
      ####################    
      if(plotST & !is.null(dataShortTerm)){
        dataST_sp = dataShortTerm[dataShortTerm$species == sp,]
        if(!is.na(dataST_sp$estimate)){
          plotST = TRUE
          
          # Extract values range of year
          yearValuesST = dataST_sp$minYear : dataST_sp$maxYear
          
          if(uncertainty){
            # Make slopes from yearly variations
            STSlopes = makeSlopesFromVar(dataVar_sp, N = N, limits = yearValuesST, weight = weight, group =  "ST")
            
            # Aesthetics
            lwdST = .8
            alphaST = .1
          }else{
            # Values of abundance for the whole years range
            estST = dataST_sp$estimate * (yearValuesST - min(yearValuesST))
            
            # Center on the middle value of the yearly variations
            midVal = mean(dataVar_sp$estimate[dataVar_sp$year >= min(yearValuesST)])
            estST = estST - mean(estST) + midVal
            
            # Turn to exponential if required
            estST = sapply(estST, function(x) ifelse(distribution == "gaussian", x, exp(x)))
            # Make dataframe
            STSlopes = data.frame(year = yearValuesST,
                                  estimate = estST,
                                  infIC = NA,
                                  supIC = NA,
                                  group = "ST",
                                  rep = 1)
            
            # Aesthetics
            lwdST = 1.1
            alphaST = 1
          }
          
        }else{
          plotST = FALSE
        }
      }else{
        plotST = FALSE
      }
      
      
      ###################
      # LONG-TERM TREND #
      ###################
      if(uncertainty){
        LTSlopes = makeSlopesFromVar(dataVar_sp, N = N, limits = yearValues, weight = weight, group =  "LT")
        # Aesthetics
        lwdLT = .8
        alphaLT = .1
        
      }else{
        # Values of abundance for the whole years range
        estLT = dataLT_sp$estimate * (yearValues - min(yearValues))
        
        # Center on middle value
        estLT = estLT - mean(estLT)
        
        # Turn to exponential if required
        estLT = sapply(estLT, function(x) ifelse(distribution == "gaussian", x, exp(x)))
        
        # Make the data.frame
        LTSlopes = data.frame(year = yearValues,
                              estimate = estLT,
                              infIC = NA,
                              supIC = NA,
                              group = "LT",
                              rep = 1)
        
        # Aesthetics
        lwdLT = 1.1
        alphaLT = 1
        
      }
      
      #####################
      # YEARLY VARIATIONS #
      #####################
      # Centrer sur le milieu de la courbe
      dataVar_sp[,c("estimate", "infIC", "supIC")] = dataVar_sp[,c("estimate", "infIC", "supIC")] - mean(dataVar_sp$estimate)
      
      # Passer à l'exponentiel, si données d'occurrence ou de comptage
      if(distribution != "gaussian"){
        dataVar_sp[,c("estimate", "infIC", "supIC")] = exp(dataVar_sp[,c("estimate", "infIC", "supIC")])
      }
      
      # Garder les colonnes d'intérêt
      dataVar_sp = dataVar_sp[,c("year","estimate", "infIC", "supIC")]
      
      # Ajouter les
      dataVar_sp$group = "VAR"
      dataVar_sp$rep = 1
      
      # FINAL DATAFRAME #
      dataPlot = rbind(LTSlopes, dataVar_sp)
      
      if(plotST){
        dataPlot = rbind(dataPlot, STSlopes) 
      }
      if(plotGamm){
        dataPlot = rbind(dataPlot, dataGammVariations_sp)
      }
      
      ###################################
      #   TITLE, SUBTITLES & AXIS NAME  #
      ###################################
      
      if(is.null(dataName)){
        title <- sp
      }else{
        # Filter french names for considered species
        dataName_sp = dataName[dataName$species == sp,]
        sc_name = paste0("(",dataName_sp$scientific_name, ")")
        
        if(nrow(dataName_sp) >0){
          fr_name = dataName_sp$french_name
          title <- bquote(.(fr_name)~italic(.(sc_name)))
          if(language == "english"){
            en_name = dataName_sp$english_name
            title <- bquote(.(en_name)~italic(.(sc_name)))
          }
        }else{
          title <- sp
        }
        
      }
      
      # Subtitle
      legendLT = percEvolutionToText_bis(data = dataLT_sp, distribution = distribution, language = language)
      
      if(plotST){
        legendST = percEvolutionToText_bis(data = dataST_sp, distribution = distribution, language = language)
      }
      
      legendGAM = ifelse(language == "french", "Tendance non linéaire", "Non-linear trend")
      legendVar = ifelse(language == "french", "Estimations annuelles", "Annual estimates")
      
      # Name of y axis
      yName = ifelse(language == "french", "Abondance relative", "Relative abundance")
      
      ######################
      #   REFERENCE LINE   # 
      ######################
      
      # Identify reference intercept
      yRef = ifelse(distribution == "gaussian", 0, 1)
      
      ####################
      #   SIGNIFICANCE   #
      ####################
      
      # Turn sinificance of the trend to text ----
      signLT <- significanceToText(dataLT_sp$pval)
      
      # x-position of the significance
      xSignLT <- dataLT_sp$maxYear + 1 ; 
      
      # y-position of the significance
      ySignLT <- mean(LTSlopes$estimate[LTSlopes$year == dataLT_sp$maxYear]) 
      
      # If ST trend exists, do the same ----
      if(plotST){
        xSignST <- xSignLT
        signST <- significanceToText(dataST_sp$pval)
        ySignST <- mean(STSlopes$estimate[STSlopes$year == dataST_sp$maxYear ]) 
      }
      
      ###################
      #   GROWTH RATE   #
      ###################
      
      # Extract growth rate
      grLT = ifelse(distribution != "gaussian", exp(dataLT_sp$estimate), dataLT_sp$estimate)
      
      # Round it
      grLT = round(grLT, 3)
      if(grLT == 1){grLT = paste0(1, ".0")}
      
      # Extract x & y coordinate
      xGrLT = dataLT_sp$minYear
      yGrLT = mean(LTSlopes$estimate[LTSlopes$year == min(LTSlopes$year)]) + 0.15
      
      if(plotST){
        # Extract growth rate
        grST = ifelse(distribution != "gaussian", exp(dataST_sp$estimate), dataST_sp$estimate)
        # Round it
        grST = round(grST, 3)
        # Extract x & y coordinate
        xGrST = dataST_sp$minYear - 1.25
        yGrST = mean(STSlopes$estimate[STSlopes$year == min(STSlopes$year)])
        
        if(grST == 1){
          grST = paste0(grST, ".0")
          yGrST = yGrST - 0.15
          xGrST = xGrST + 1.25
          
        }
        
      }
      
      ##############
      #   COLORS   #
      ##############
      if (plotGamm & plotST){
        # Color values
        colValues = c( "#F8C457","#429EBD","#FF4359", "#F8C457")
        colValues_guide = c("#429EBD","#FF4359", "#F8C457", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendST, legendGAM, legendVar)
        
        # Other characteristics
        linetype = c(1, 1, 1, 0) ; shape = c(NA, NA, NA, 18)
        size = c(1.1, 1.1, 1.1, 3) ; alpha = c(1, 1, 1, 1)
        
      }else if(plotGamm & !plotST){
        # Color values
        colValues = c( "#F8C457","#429EBD", "#F8C457") 
        colValues_guide = c("#429EBD", "#F8C457", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendGAM, legendVar)
        
        # Other characteristics
        linetype = c(1, 1, 0) ; shape = c(NA, NA, 18)
        size = c(1.1, 1.1, 3) ; alpha = c(1, 1, 1)
        
      }else if(!plotGamm & plotST){
        # Color values
        colValues = colValues_guide = c("#429EBD","#FF4359", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendST, legendVar) 
        
        # Other characteristics
        linetype = c(1, 1, 0) ; shape = c(NA, NA, 18)
        size = c(1.1, 1.1, 3) ; alpha = c(1, 1, 1)
        
      }else{
        # Color values
        colValues = colValues_guide = c("#429EBD", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendVar) 
        
        # Other characteristics
        linetype = c(1, 0) ; shape = c(NA, 18)
        size = c(1.1, 3) ; alpha = c(1, 1)
      }
      
      ############
      #   PLOT   #
      ############
      if(all(dataPlot$estimate < 10^10)){
        
        # Initialiser le plot 
        plot <- ggplot2::ggplot()
        
        # ADD GAMM INFORMATION
        if(plotGamm){
          plot <- plot +  
            # Add the confidence interval
            ggplot2:: geom_ribbon(data = dataPlot[dataPlot$group == "GAM",], ggplot2::aes(x = year, ymin = infIC, ymax = supIC), fill = "#F8C457", alpha = .2) +
            
            # Add the curve associated with the GAM
            ggplot2::geom_line(data = dataPlot[dataPlot$group == "GAM",], ggplot2::aes(x = year, y = estimate, col = group), linewidth = 1.05)
          
        }
        
        # ADD YEARLY VARIATIONS & LONG-TERM TREND
        plot = plot +
          
          # Add yearlyVariations
          ggplot2::geom_point(data = dataPlot[dataPlot$group == "VAR",], ggplot2::aes(x = year, y = estimate, col = group), shape = 18, size = 2.5) + 
          
          # Add the long-term trend
          ggplot2::geom_line(data = dataPlot[dataPlot$group == "LT",], ggplot2::aes(x = year, y = estimate, group = rep, col = group), 
                             alpha = alphaLT, linewidth = lwdLT) + 
          
          # Add LT significance information
          ggplot2::annotate(geom = "text", label = signLT, color = "#429EBD", size = 8, x = xSignLT, y = ySignLT) +
          
          # Add LT growth rate
          ggplot2::annotate(geom = "text", label = grLT, color = "#429EBD", size = 8, x = xGrLT, y = yGrLT)
        
        # ADD SHORT-TERM TREND
        if(plotST){
          plot = plot +
            
            # Add the short-term trend
            ggplot2::geom_line(data = dataPlot[dataPlot$group == "ST",], ggplot2::aes(x = year, y = estimate, group = rep, col = group), 
                               alpha = alphaST, linewidth = lwdST) + 
            
            # Add LT significance information
            ggplot2::annotate(geom = "text", label = signST, color = "#FF4359", size = 8, x = xSignST, y = ySignST) +
            
            # Add LT growth rate
            ggplot2::annotate(geom = "text", label = grST, color = "#FF4359", size = 8, x = xGrST, y = yGrST)
          
          
        }
        
        # FORMAT THE PLOT : AXIS, TITLES, ...
        plot <-  plot + 
          
          # Add the reference line
          ggplot2::geom_hline(yintercept = yRef, col = "gray", linetype = "dashed", linewidth = 1.01) +
          
          # Title
          ggplot2::ggtitle(title) +
          
          # Format x and y-axis
          ggplot2::scale_x_continuous(ifelse(language == "french", "Année", "Year"), 
                                      breaks = yearValues, limits = c(dataLT_sp$minYear - 1, dataLT_sp$maxYear +1)) + 
          
          # Rename y-axis                            
          ggplot2::ylab(yName) +
          
          # Change theme
          ggplot2::theme_bw() +
          
          # Adapt y-axis range
          ggplot2::coord_cartesian(ylim =  c(0.5, 2)) +  
          
          # Format theme
          ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 1, hjust = 1), 
                         legend.text = ggplot2::element_text(size = 17),
                         axis.text = ggplot2::element_text(size = 17),
                         axis.title = ggplot2::element_text(size = 18),
                         title = ggplot2::element_text(size = 21), 
                         panel.grid.minor = ggplot2::element_blank(),
                         legend.position = "top", 
                         legend.justification = "left") +
          
          # Format type of trend colour legend
          ggplot2::scale_colour_manual(name = "", 
                                       values = colValues, 
                                       labels = colLabels, 
                                       guide = ggplot2::guide_legend(nrow = 2, 
                                                                     by.row = FALSE,
                                                                     override.aes = list(linetype = linetype,
                                                                                         color = colValues_guide,
                                                                                         shape = shape,
                                                                                         size = size,
                                                                                         alpha = alpha))) 
      } else{
        plot = ggplot2::ggplot() + ggplot2::theme_bw()
        cat("No plot for species ", sp," due to extreme values in some coefficients\n")
      }
      
    }
    
    # Save results
    if(save){
      ggplot2::ggsave(plot = plot, device = "png", filename = paste0(sp, ".png"), 
                      path = path, width = 25, height = 18, units = "cm")
    }
    
    return(plot) 
  }
  )
  
  return(plots)
}

# DATA IMPORT & FORMATTING ----
# Load data
data = read.csv("data/VF_2023/countingData.csv")

# Extract all species to analyze in the data.frame
speciesList <- c("Hedera helix", "Lolium perenne", "Rubus fruticosus",
                 "Dactylis glomerata", "Galium aparine", "Plantago lanceolata",
                 "Convolvulus arvensis", "Quercus robur", "Fraxinus excelsior",
                 "Trifolium repens", "Polygonum aviculare", "Fallopia convolvulus",
                 "Carpinus betulus", "Crataegus monogyna", "Holcus lanatus",
                 "Urtica dioica", "Arrhenatherum elatius", "Cirsium arvense",
                 "Geranium dissectum", "Daucus carota", "Geum urbanum",
                 "Lysimachia arvensis", "Lonicera periclymenum", "Chenopodium album",
                 "Poa annua",  "Prunus spinosa", "Geranium robertianum",
                 "Plantago major", "Picris hieracioides", "Anisantha sterilis",
                 "Achillea millefolium", "Trifolium pratense", "Acer pseudoplatanus",
                 "Poa trivialis", "Veronica persica", "Heracleum sphondylium",
                 "Potentilla reptans", "Medicago lupulina", "Fagus sylvatica",
                 "Pteridium aquilinum", "Cornus sanguinea", "Corylus avellana",
                 "Ranunculus repens", "Vicia sativa", "Rubia peregrina",
                 "Ligustrum vulgare", "Sonchus asper", "Brachypodium sylvaticum",
                 "Glechoma hederacea", "Alopecurus myosuroides","Acer campestre",
                 "Clematis vitalba", "Poa pratensis", "Hypochaeris radicata",
                 "Lotus corniculatus", "Viola arvensis", "Avena fatua",
                 "Geranium molle", "Castanea sativa", "Rumex acetosa",
                 "Hypericum perforatum", "Bellis perennis", "Papaver rhoeas",
                 "Taraxacum officinale", "Senecio vulgaris", "Bromus hordeaceus",
                 "Elytrigia repens", "Lapsana communis", "Convolvulus sepium",
                 "Ranunculus acris", "Veronica arvensis", "Agrostis stolonifera",
                 "Tripleurospermum inodorum","Solanum nigrum", "Sonchus oleraceus",
                 "Carex sylvatica", "Agrostis capillaris", "Fragaria vesca",
                 "Equisetum arvense", "Rubus caesius", "Cerastium fontanum",
                 "Prunella vulgaris", "Teucrium scorodonia", "Stellaria media",
                 "Galium mollugo", "Poa nemoralis", "Poterium sanguisorba",
                 "Quercus ilex", "Ulmus minor",  "Quercus petraea",
                 "Jacobaea vulgaris", "Ilex aquifolium", "Prunus avium",
                 "Anthoxanthum odoratum", "Lolium multiflorum", "Persicaria maculosa",
                 "Artemisia vulgaris", "Quercus pubescens", "Euonymus europaeus",
                 "Erigeron canadensis", "Festuca rubra", "Cirsium vulgare",
                 "Veronica chamaedrys", "Agrimonia eupatoria", "Crepis capillaris",
                 "Matricaria chamomilla", "Capsella bursa-pastoris","Polygonatum multiflorum",
                 "Rumex obtusifolius", "Rosa canina",  "Schedonorus arundinaceus",
                 "Silene latifolia", "Viola reichenbachiana", "Ervilia hirsuta",
                 "Geranium rotundifolium","Helminthotheca echioides", "Molinia caerulea",
                 "Ruscus aculeatus", "Rumex crispus", "Lactuca serriola",
                 "Myosotis arvensis", "Geranium pusillum", "Holcus mollis",
                 "Asparagus acutifolius", "Origanum vulgare", "Stellaria holostea",
                 "Carex flacca", "Alliaria petiolata", "Calluna vulgaris",
                 "Mercurialis annua", "Cytisus scoparius", "Leucanthemum vulgare",
                 "Phillyrea angustifolia","Pinus sylvestris", "Mercurialis perennis",
                 "Triticum aestivum", "Euphorbia amygdaloides","Melica uniflora",
                 "Sherardia arvensis", "Dioscorea communis", "Ranunculus bulbosus",
                 "Teucrium chamaedrys", "Hordeum murinum", "Anemone nemorosa",
                 "Anthriscus sylvestris", "Smilax aspera", "Epilobium tetragonum",
                 "Thymus vulgaris", "Buxus sempervirens", "Medicago sativa",
                 "Circaea lutetiana", "Frangula alnus", "Phleum pratense",
                 "Rubus ulmifolius", "Robinia pseudoacacia", "Ajuga reptans",
                 "Vicia sepium", "Acer platanoides", "Euphorbia cyparissias",
                 "Brachypodium retusum", "Aphyllanthes monspeliensis", "Betula pendula",
                 "Bromopsis erecta", "Sinapis arvensis", "Viburnum lantana",
                 "Brachypodium pinnatum", "Sambucus nigra", "Pilosella officinarum",
                 "Potentilla erecta", "Pastinaca sativa", "Ulex europaeus",
                 "Convallaria majalis", "Sorbus torminalis", "Viola odorata",
                 "Galeopsis tetrahit", "Hyacinthoides non-scripta", "Moehringia trinervia",
                 "Rosa arvensis", "Zea mays", "Quercus coccifera",
                 "Sorbus aucuparia")

# Add 0's to the data.frame
data = data[,c("ID","species","point", "year", "site", "julian", "r_resp", "abond_cumul", "fail")]
data = fillAbsence(data, interestVar = c("abond_cumul", "fail"), speciesList, method = "once")

# Make a dataframe with correspondency of species ID, scientific and english name
dataName = data.frame(species = c("Hedera helix","Clematis vitalba", "Rubus fruticosus", "Buxus sempervirens", "Lolium perenne"),
                      scientific_name = c("Hedera helix","Clematis vitalba", "Rubus fruticosus",
                                          "Buxus sempervirens", "Lolium perenne"),
                      english_name = c("Common ivy","Traveller's joy", "Blackberries",
                                       "Common box", "English ryegrass"))

# Trend formatting
speciesList = c("Buxus sempervirens", "Rubus fruticosus","Lolium perenne", "Clematis vitalba")

dataLongTermTrend = formatTrendEstimates(data, speciesList, repo = repo, "longTermTrend")
dataYearlyVariations = formatYearlyVariations(speciesList, repo = repo, "yearlyVariations", contr = "mean")
dataShortTermTrend = formatTrendEstimates(data[data$year >= max(data$year)-10,], speciesList, repo = repo, "shortTermTrend")

# Figures
plotLinearTrends_bis(speciesList, dataLongTermTrend, dataYearlyVariations, dataShortTermTrend,
                     language = "english", dataName = dataName, distribution = "betabinomial", plotGamm = TRUE, plotST = TRUE,
                     uncertainty = TRUE, N = 100, weight = TRUE, save = TRUE, path = result_path)
