

# Victor: verifier que les groupes de colonies sont toujours dans le bon ordre dans les outputs. Est-ce qu'il faudrait leur donner des noms plutôt qu'un numero?
# Victor: label doesn't work for point data



# Load libraries
library(tidyverse)
library(plyr)
library(shiny)
library(bslib) # For accordion
library(plotly)
library(rmarkdown)
library(markdown)
library(DT)
library(leaflet) # Needed for interactive map
library(mapedit) # Needed to edit interactive map
library(leaflet.extras) # For addMouseCoordinates
library(leafem)
library(gdistance) # Needed to calculate shortest path distance
library(bsicons) # Needed for accordion icons
library(spsComps) # Needed for loader
library(shinyjs) # Needed to enable / disable button
library(sf) # Needed for GIS work
library(raster) # Needed for GIS work (extent function)
library(tools) # Needed to extract extensions
library(units)
library(htmlwidgets) # Needed for onRender in leaflet map
library(reshape2) # Needed to prepare output tables

setwd(dir=ifelse(file.exists("C:/Users/TERRA"),
                 "C:/Users/TERRA/Documents/Birdynamic/Birdynamic/app", # Setwd Victor
                 "D:/Documents/JOBS/Independant/Birdynamic/app")) # Setwd Baptiste

plan(multisession)

# Load functions
source("Functions_Birdynamic_Shiny.R")

### Charger les donnees fixes
load("0.Data/clean_data_all_species.RData")
load("0.Data/spatial/transitionMatrix_CostRaster.rda")
#load("0.Data/clean_data_parcs_StNazaire.RData") #for testing
load("0.Data/clean_data_parcs.RData") #for testing
collisions <- read.csv2(file="../Test_data_Shiny/04_collisionRisk_Cormoran huppé.csv") #for testing
LimFrance <- st_read("0.Data/Limites_France/Limites_France_EEZ_extended.shp")

liste_especes = BD_list$espece_BD
liste_especes_latin = c(BD_list$Nom_latin,BD_list$Nom_latin_2[which(!is.na(BD_list$Nom_latin_2))])

## Theme colours
COL_prim = "#446e9b"
COL_sec = "#78c2ad"


### UI ----------
ui <- page_fillable(
  useShinyjs(), # Needed to enable / disable buttons
  
  theme = bs_theme(version = 5, bootswatch = "sandstone", primary=COL_prim, secondary=COL_sec),
  
  
  tags$style(paste0("
     .accordion .accordion-button {
    color: white;
    background-color: ", COL_prim, ";
     }")),
  
  
  accordion_filters <- accordion(open=F, # FALSE to have no accordion opened at start; NULL to have the first one
                                 multiple=F,
                                 
                                 accordion_panel("Comment utiliser l'application ?", icon = bsicons::bs_icon("question-circle"),
                                                 HTML("<h3>Méthode</h3>"), 
                                                 TEXT_intro,
                                                 HTML("<br><br><h3>Cheatsheet</h3>"), 
                                                 includeMarkdown("www/MD_cheatsheet.md"),
                                                 HTML("<br><br><h3>Tutoriel video</h3>"), 
                                                 HTML('<iframe width="560" height="315" src="https://www.youtube.com/embed/T1-k7VYwsHg" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>')
                                 ),
                                 
                                 accordion_panel("Données de mortalité", icon = bsicons::bs_icon("table"),
                                                 htmlOutput("Warning_mortality"),
                                                 fluidRow(
                                                   column(4, align="center",
                                                          card(fileInput("Upload_Mortality", "Chargez vos données de mortalité", accept=c('.csv'), multiple=FALSE, buttonLabel="Sélectionner", placeholder="Aucun fichier", width="90%")),
                                                          
                                                   ),
                                                   column(8,align="center",
                                                          DT::dataTableOutput("Mortality_Data")
                                                   )
                                                 )
                                 ),
                                 
                                 accordion_panel("Coordonnées du parc éolien", icon = bsicons::bs_icon("map"),
                                                 htmlOutput("Warning_parc"),
                                                 fluidRow(
                                                   column(4, align="center",
                                                          card(fileInput("Upload_Shp", HTML("<b>Charger le shapefile du parc</b>"), buttonLabel="Sélectionner", accept=c('.shp','.dbf','.sbn','.sbx','.shx',".prj"), multiple=TRUE)),
                                                          card(HTML("<b><br>Ajouter un parc manuellement</b>"),
                                                               numericInput(inputId="NewPtsLon", label="Longitude", value=NULL),
                                                               numericInput(inputId="NewPtsLat", label="Latitude", value=NULL),
                                                               textInput(inputId="NewPtsName", label="Nom du parc", value=""),
                                                               actionButton(inputId="NewPtsButt", label="Ajouter"))
                                                   ),
                                                   column(8, align="center",
                                                          editModUI("map", height=600)
                                                   )
                                                 )
                                 ),
                                 
                                 accordion_panel("Analyse", icon = bsicons::bs_icon("graph-down"),
                                                 htmlOutput("Warning_general"),
                                                 fluidRow(
                                                   numericInput("iteration_nbr","Nombre d'itérations", 
                                                                value = 1000, 
                                                                min = 100, 
                                                                max = 10000),
                                                   selectInput("selected_species", 
                                                               "Choisir les espèces à analiser:", 
                                                               NULL,
                                                               multiple=T),
                                                   fileInput("Upload_distances", 
                                                             "Fichier de distance issu d'une simulation précédente (facultatif)", 
                                                             accept=c('.csv'), multiple=FALSE, 
                                                             buttonLabel="Sélectionner", 
                                                             placeholder="Aucun fichier", width="30%"),
                                                   #this is to make testing easier, to be removed in final app
                                                   checkboxInput("do_euclidian","Use euclidian? (for testing)",T),
                                                 ),
                                                 htmlOutput("Warning_distances"),
                                                 actionButton('RunButt', 'Analyser', icon=icon("arrow-right", lib="font-awesome")),
                                                 conditionalPanel(condition = "input.RunButt >= 1", # Show loader and results only if we clicked on calculate
                                                                  fluidRow(
                                                                    DTOutput("TableNational")
                                                                  ),
                                                                  fluidRow(
                                                                    column(6, 
                                                                           fluidRow(plotOutput("FigureSummary1", width="90%", height="300px")),
                                                                           fluidRow(plotOutput("FigureSummary2", width="90%", height="300px"))),
                                                                    column(6, leafletOutput("MapSummary", width="90%", height="600px"))
                                                                  ),
                                                                  HTML("<br><br><h3>Télécharger les résultats</h3><br>"),
                                                                  downloadButton("downloadOutputs", "Télécharger")
                                                 )
                                 )
  ),
  
  includeMarkdown("www/MD_include_logos.md")
)




### SERVER -----
server <- function(input, output, session) {
  
  #increase file upload limit to 50MB
  options(shiny.maxRequestSize=50*1024^2)
  
  ##### INITALISATION -----
  
  ### Create reactive values
  LoadedMortality <- reactiveVal(data.frame()) # Mortality csv loaded by user
  LoadedDistances <- reactiveVal(NULL) # csv of distances between parcs and colonies uploaded by user from previous simulation
  FormattedMortality <- reactiveVal(data.frame()) # Mortality table once formated
  FOLDER_name <- reactiveVal(paste0("test_", as.character(Sys.Date()))) # Name of the output files
  Warning_mortality <- reactiveVal("<p style='color:red'>Les données de mortalité doivent être chargées pour lancer l'analyse</p>") # Warning if problem in mortality data
  Warning_parc <- reactiveVal("<p style='color:red'>Les coordonnées d'un parc éolien doivent être fournies pour lancer l'analyse</p>") # Warning if problem in parc data
  Warning_general <- reactiveVal("<p style='color:red'>Données de parc et de mortalité manquantes pour lancer l'analyse</p>")
  Warning_distances <- reactiveVal("") # non-blocking warning if problem in distances table
  Stored_results <- reactiveValues('Res_list'=list()) # Plots and tables stored (I use Res_list inside to be able to reset to empty)
  Shapefile_stored <- reactiveValues('Shp_list'=list(Upload_shp=NULL, Drawn_shp=NULL, Longlat_pts=NULL, Merged_shp=NULL)) # Plots and tables stored (I use Shp_list inside to be able to reset to empty)
  Last_shp <- reactiveVal(NULL) # Becomes "Point" if the last parc edit is a drawn point, and "Upload" if it's an uploaded shapefile; used to choose whether the interactive map is rescaled or not (only if uploaded)
  Parc_names_mortality <- reactiveVal(NULL) #stores the names of the parcs read from the collision table
  Species_names_mortality <- reactiveVal(NULL) #stores the names of the species read from the collision table
  
  ### Create interactive map
  edits <- Bdy_CreateMap()
  
  
  ##### EVENTS -----
  
  ### Enable / disable RunButt button depending if mortality and parc data are ready and valid
  observe({
    if(Warning_general()==""){ # Will be TRUE if at least one of the two data sources is invalid
      shinyjs::enable("RunButt") # Need to integrate useShinyjs in UI
    }else {
      shinyjs::disable("RunButt")
    }
  })
  
  ### Load mortality ----
  observeEvent(input$Upload_Mortality, {
    Bdy_log("START - Mortality data uploaded")
    
    ### Load data (if only 1 column, meaning we have a problem with separator, it tries other separators)
    Mortality_raw <- read.csv2(input$Upload_Mortality$datapath, header = TRUE, check.names = FALSE)
    if(ncol(Mortality_raw)==1){Mortality_raw = read.csv2(input$Upload_Mortality$datapath, header = TRUE, sep=",", check.names = FALSE)}
    if(ncol(Mortality_raw)==1){Mortality_raw = read.csv2(input$Upload_Mortality$datapath, header = TRUE, sep=";", check.names = FALSE)}
    if(ncol(Mortality_raw)==1){Mortality_raw = read.csv2(input$Upload_Mortality$datapath, header = TRUE, sep="\t", check.names = FALSE)}
    
    ### Update data warning and format the table if it is correct
    checkMorta <- Bdy_CheckMortality(Mortality_raw,liste_especes_latin)
    Mortality_ready <- checkMorta$table
    Warning_mortality(checkMorta$warning)
    
    Warning_parc(Bdy_CheckShapefile(Shapefile_stored$Shp_list$Merged_shp, unique(Mortality_ready$Parc)))
    
    #if the warning is green then the appropriate columns are present
    if(grepl("color:green", Warning_mortality())){
      
      #adding new column with french names 
      Mortality_ready$Espece <- BD_list$espece_BD[match(Mortality_ready$espece_latin,liste_especes_latin)]
      
      Parc_names_mortality(unique(Mortality_ready$parc)) 
      Species_names_mortality(unique(Mortality_ready$Espece))
      updateSelectInput(session,inputId ="selected_species",choices=Species_names_mortality(),selected=Species_names_mortality())
    }
    
    ### Update LoadedMortality interactive value
    LoadedMortality(Mortality_raw)
    FormattedMortality(Mortality_ready)
    
    Bdy_log("END - Mortality data uploaded")
    
  })
  
  ### Load parc data (4 observeEvents) ----
  # 1st: when new shapefile uploaded
  observeEvent(input$Upload_Shp, {
    Bdy_log("START - Upload parc shapefile")
    
    Upload_shp <- Bdy_UploadShapefile(input$Upload_Shp)
    req(nrow(Upload_shp)>0)
    
    if(! "NAME" %in% names(Upload_shp)){Upload_shp$NAME <- ""}
    
    Shapefile_stored$Shp_list$Upload_shp <- rbind(
      Shapefile_stored$Shp_list$Upload_shp[,c("leaflet_id", "NAME", "geometry")], # Previously uploaded shapefile (to allow uploading 2 parcs)
      Upload_shp[,c("leaflet_id", "NAME", "geometry")] # New upload
    )
    Last_shp("Upload")
    
    Bdy_log("END - Upload parc shapefile")
  })
  
  #2nd: when user draw points on the interactive map
  observeEvent(edits()$finished, {
    req(nrow(edits()$finished)>0)
    req((F %in% (edits()$finished$`_leaflet_id` %in% Shapefile_stored$Shp_list$Drawn_shp$leaflet_id)) | (T %in% (edits()$deleted$`_leaflet_id` %in% Shapefile_stored$Shp_list$Drawn_shp$leaflet_id)))
    
    Bdy_log("START - Draw parc points")
    
    Shapefile_stored$Shp_list$Drawn_shp <- edits()$finished %>% mutate(NAME="")
    Shapefile_stored$Shp_list$Drawn_shp$leaflet_id <- Shapefile_stored$Shp_list$Drawn_shp$`_leaflet_id`
    
    Last_shp("Point")
    Bdy_log("END - Draw parc points")
    
  })
  observeEvent(edits()$deleted, {Last_shp("Deleted")})
  
  # 3rd: when a point is added with Lon/Lat
  observeEvent(input$NewPtsButt, {
    req(input$NewPtsLon) ; req(input$NewPtsLat) ; req(input$NewPtsName)
    
    Bdy_log("START - Add points from Lon / Lat")
    
    Num_pt <- ifelse(is.null(nrow(Shapefile_stored$Shp_list$Longlat_pts)), 1, (nrow(Shapefile_stored$Shp_list$Longlat_pts)+1))
    NewPt <- st_as_sf(data.frame(Lon=input$NewPtsLon, Lat=input$NewPtsLat, leaflet_id=paste0("LonLat_", Num_pt), NAME=input$NewPtsName), coords = c("Lon","Lat"), crs=st_crs(4326))
    
    Shapefile_stored$Shp_list$Longlat_pts <- rbind(Shapefile_stored$Shp_list$Longlat_pts, NewPt[,c("leaflet_id", "NAME", "geometry")])
    Last_shp("LonLat")
    
    Bdy_log("END - Add points from Lon / Lat")
    
  })
  
  # 4th: each time Event 1 / 2 / 3 occurs, merge results and update leaflet map
  observeEvent(c(Shapefile_stored$Shp_list, edits()$deleted), {
    req(Last_shp() != "Merged")
    Bdy_log("START - Merge parc shapefile")
    
    ### Merge uploaded shapefile and drawn points
    Merged_shp <- rbind(Shapefile_stored$Shp_list$Upload_shp[,c("leaflet_id", "NAME", "geometry")], 
                        Shapefile_stored$Shp_list$Drawn_shp[,c("leaflet_id", "NAME", "geometry")]) %>% rbind(
                          .,
                          Shapefile_stored$Shp_list$Longlat_pts[,c("leaflet_id", "NAME", "geometry")]
                        )
    
    # Take edited names
    Old_Merged <- Shapefile_stored$Shp_list$Merged_shp
    Merged_shp$NAME[Merged_shp$leaflet_id %in% Old_Merged$leaflet_id] <- Old_Merged$NAME[match(Merged_shp$leaflet_id[Merged_shp$leaflet_id %in% Old_Merged$leaflet_id], Old_Merged$leaflet_id)]
    req(is.null(nrow(Merged_shp))==F)
    
    # Subset points or polygons removed (based on their leaflet_id or geometry)
    if(is.null(nrow(edits()$deleted))==F){
      Merged_shp <- Merged_shp %>%
        subset(., ! .$leaflet_id %in% edits()$deleted$leaflet_id) %>%
        subset(., ! .$leaflet_id %in% st_filter(st_make_valid(Merged_shp), st_make_valid(edits()$deleted), .predicate=st_intersects)$leaflet_id) %>%
        subset(., ! .$leaflet_id %in% st_filter(st_buffer(st_make_valid(Merged_shp[st_geometry_type(Merged_shp)=="POINT",]), 1), st_make_valid(edits()$deleted), .predicate=st_intersects)$leaflet_id) %>% # For points I need a buffer to st_intersects the points in "deleted"
        st_make_valid(.)
    }
    
    ### TEMPO plot to see the data considered for analyses
    Stored_results$Res_list$Figure1 <- ggplot(Merged_shp)+geom_sf()
    
    ### Check data is ok
    # Check all in allowed zone
    Merged_shp$CoverLimit <- lengths(st_covered_by(Merged_shp, st_transform(LimFrance, st_crs(Merged_shp))))
    Warning_parc(Bdy_CheckShapefile(Merged_shp, Parc_names_mortality(), Merged_shp$CoverLimit))
    
    ### Add popup
    if(is.null(nrow(Merged_shp))==F){if(nrow(Merged_shp)>0){
      Merged_shp$Color <- ifelse(is.na(Merged_shp$NAME) | (Merged_shp$NAME == "") | (! Merged_shp$NAME %in% Parc_names_mortality()), "#CD34B5", "#409808") %>% ifelse(Merged_shp$CoverLimit>0, ., "#CF080D")
      Merged_shp$Popup <- Bdy_MapPopup(Merged_shp)
      #Merged_shp$CoverLimit<-NULL
    }}
    
    ### Store
    Shapefile_stored$Shp_list$Merged_shp <- Merged_shp
    
    ### Add to the interactive map
    if(nrow(Merged_shp)>0){
      
      # Get extent to map
      EXT1 <- 3*extent(Merged_shp)
      EXT2 <- c(NA, NA, NA, NA)
      if(is.null(nrow(edits()$finished))==F){EXT2 <- extent(edits()$finished)}
      EXT_shp <- c(min(EXT1[1], EXT2[1], na.rm=T), max(EXT1[2], EXT2[2], na.rm=T), min(EXT1[3], EXT2[3], na.rm=T), max(EXT1[4], EXT2[4], na.rm=T))
      
      # If some features are deleted, create a new map (there are some bugs in leaflet functions to remove polygons one by one in such cases)
      if(is.null(nrow(edits()$deleted))==F){
        edits <- Bdy_CreateMap(EXT_shp[c(1,3,2,4)])
      }
      
      # Split polygons and points
      Poly_toadd <- Merged_shp[st_geometry_type(Merged_shp) %in% c("POLYGON", "MULTIPOLYGON"),]
      Pts_toadd <- Merged_shp[st_geometry_type(Merged_shp)  %in% c("POINT", "MULTIPOINT"),]
      
      leafletProxy("map-map") %>% 
        clearGroup("Parc") %>% # Needed to avoid overwritting the same polygons multiple times
        addPolygons(data = Poly_toadd, group="Parc", color=Poly_toadd$Color, fillOpacity=1, popup=Poly_toadd$Popup, popupOptions = popupOptions(autoPan=T, closeOnClick=T)) %>%
        addCircleMarkers(data = Pts_toadd, group="Parc", color=Pts_toadd$Color, fillOpacity=1, popup=Pts_toadd$Popup, popupOptions = popupOptions(autoPan=T, closeOnClick=T)) %>%
        addLabelOnlyMarkers(data=st_centroid(Merged_shp), group="Parc", label=Merged_shp$NAME, labelOptions=labelOptions(permanent=T, textOnly=T, direction="top"))
      
      # Rescale map if shapefile uploaded or lon/lat point added (not if I just drew a point)
      if(Last_shp() %in% c("Upload", "LonLat")){
        leafletProxy("map-map") %>% 
          fitBounds(lng1 = EXT_shp[1], lng2 = EXT_shp[2], lat1 = EXT_shp[3], lat2 = EXT_shp[4])
      }
      Last_shp("Merged")
    } #else{edits <- Bdy_CreateMap()}
    
    Bdy_log("END - Merge parc shapefile")
    
  })
  
  ### Edit name
  observeEvent(c(input$text_onchange, input$Upload_Mortality), { # Includes Upload_mortality to update warning when mortality file is changed
    req(is.null(nrow(Shapefile_stored$Shp_list$Merged_shp))==F)
    
    Bdy_log("START - Edit names")
    
    # Edit stored shapefile
    Merged_shp <- Shapefile_stored$Shp_list$Merged_shp
    Merged_shp$NAME[Merged_shp$leaflet_id==sub("EditName_", "", input$text_onchange$id)] <- input$text_onchange$value
    Merged_shp$Color <- ifelse(is.na(Merged_shp$NAME) | (Merged_shp$NAME == "") | (! Merged_shp$NAME %in% Parc_names_mortality()), "#CD34B5", "#409808") %>% ifelse(Merged_shp$CoverLimit>0, ., "#CF080D")
    Merged_shp$Popup <- Bdy_MapPopup(Merged_shp)
    Shapefile_stored$Shp_list$Merged_shp <- Merged_shp
    Warning_parc(Bdy_CheckShapefile(Merged_shp, Parc_names_mortality()))
    
    # Edit leaflet
    Poly_toadd <- Merged_shp[st_geometry_type(Merged_shp) %in% c("POLYGON", "MULTIPOLYGON"),]
    Pts_toadd <- Merged_shp[st_geometry_type(Merged_shp)  %in% c("POINT", "MULTIPOINT"),]
    
    leafletProxy("map-map") %>% 
      clearGroup("Parc") %>% # Needed to avoid overwritting the same polygons multiple times
      addPolygons(data = Poly_toadd, group="Parc", color=Poly_toadd$Color, fillOpacity=1, popup=Poly_toadd$Popup, popupOptions = popupOptions(autoPan=T, closeOnClick=T)) %>%
      addCircleMarkers(data = Pts_toadd, group="Parc", color=Pts_toadd$Color, fillOpacity=1, popup=Pts_toadd$Popup, popupOptions = popupOptions(autoPan=T, closeOnClick=T)) %>%
      addLabelOnlyMarkers(data=st_centroid(Merged_shp), group="Parc", label=Merged_shp$NAME, labelOptions=labelOptions(permanent=T, textOnly=T, direction="top"))
    
    Bdy_log("END - Edit names")
  })
  
  
  
  ### Update Warning_general when warnings from spatial data or mortality data are updated
  observeEvent(c(Warning_mortality(), Warning_parc(),input$selected_species), {
    Bdy_log("START - Update Warning_general")
    
    Warning_general(
      
      if(grepl("color:red", Warning_mortality())){
        if(grepl("color:red", Warning_parc())){
          "<p style='color:red'>Données de parc et de mortalité manquantes pour lancer l'analyse</p>"
        }else{
          "<p style='color:red'>Données de mortalité manquantes pour lancer l'analyse</p>"
        }
      }else if(grepl("color:red", Warning_parc())){
        "<p style='color:red'>Données de parc manquantes pour lancer l'analyse</p>"
      }else if(is.null(input$selected_species)){
        "<p style='color:red'>Au moins une espèce doit être sélectionnée pour lancer l'analyse</p>"
      }else{
        ""
      }
      
    )
    
    #checking again the loaded distance table
    Warning_distances(Bdy_checkDistances(LoadedDistances(),colonies00_L93$code_colonie,Parc_names_mortality()))
    
    Bdy_log("END - Update Warning_general")
  })
  
  ### load and check shortest path distance from previous simulation
  observeEvent(input$Upload_distances,{
    Distances_raw <- read.csv(input$Upload_distances$datapath, header = TRUE, check.names = FALSE,row.names = 1)
    Warning_distances(Bdy_checkDistances(Distances_raw,colonies00_L93$code_colonie,Parc_names_mortality()))
    LoadedDistances(Distances_raw[])
  })
  
  ### Run analyses ----
  observeEvent(input$RunButt, {
    Bdy_log("START - Run analyses")
    
    ### Start loader
    loader_analyses <- addLoader$new("RunButt", color = COL_prim, method = "replace", height = "30rem", opacity=0.4) ; loader_analyses$show()
    
    #checking if we use the loaded distance table or not
    if(grepl("color:purple",Warning_distances())){
      LoadedDistances(NULL)
    }else{
      LoadedDistances <- LoadedDistances()[,which(colnames(LoadedDistances()) %in% Parc_names_mortality())]
    }
    
    ### Transform parc to Lambert projection
    parcs_L93 <- st_transform(Shapefile_stored$Shp_list$Merged_shp, st_crs(2154))
    
    moeurs_terrestres <- c("Goéland marin", "Goéland brun", "Goéland argenté", "Grand cormoran", "Mouette mélanocéphale", 
                           "Mouette rieuse", "Goéland cendré", "Sterne caugek", "Sterne pierregarin", "Guifette noire", 
                           "Goéland leucophée", "Mouette pygmée")
    #species present in less than 10 colonies, mostly for testing
    espèces_rare <- c("Fou de Bassan","Goéland cendré","Guifette noire","Guillemot de Troïl",
                      "Macareux moine","Pingouin torda","Puffin des Anglais")
    
    #sel_sp = unique(collisions$espece) #for testing
    sel_sp <- input$selected_species
    sel_sp_marine <- sel_sp[which(!sel_sp %in% moeurs_terrestres)]
    ## combining the colonies table for each bird with marine habits into one 
    ## useful to avoid recalculating distances
    first = T
    #for(sp in input$selected_species){
    for(sp in sel_sp_marine){ 
      load(paste0("0.Data/processed_data_per_species/clean_data_",sp,".RData"))
      if(first){
        colonies_combined = colonies
        first = F
      }else{
        colonies_combined = rbind(colonies_combined,colonies)
      }
      rm(colonies,group_counts,ppa_yr)
    }
    
    #If at least one species with moeurs marines
    if(length(sel_sp_marine) != 0){
      #removing duplicated colonies (if both species are present in the same)
      colonies_combined <- colonies_combined[which(!duplicated(colonies_combined$code_colonie)),]
      
      #progress bar for distance calculations
      progressDistances <- shiny::Progress$new()
      progressDistances$set(message = "Calcul des distances entre colonies et parc éoliens", value = 0)
      
      #calculating shortest path distances for theses species
      maxRange = max(foraging_range_tbl$Max_km[match(sel_sp_marine,foraging_range_tbl$Espece)],na.rm=T)
      shortest_path_distances <- Bdy_GetColoParcDistances(colonies=colonies_combined,
                                                          parcs = parcs_L93,
                                                          costMatrix=trCost,
                                                          #loadedDistances=LoadedDistances(),
                                                          maxForagingRange=maxRange,
                                                          euclidian=input$do_euclidian, #for testing
                                                          progress=progressDistances)
      Sys.sleep(1)
      progressDistances$close()
    }
    
    model_output <- list()
    
    #progress bar for going through species
    progressSpecies <- shiny::Progress$new()
    progressSpecies$set(message = "Analyses pour l'espèce:", value = 0)
    
    out <- list() ; length(out) <- length(sel_sp)
    
    for(sp in sel_sp)
    {
      out[[sp]] <- future({
      
      print(paste0("Starting analysis for ",sp))
      
      #sp = BD_list$espece_BD[1] #for testing
      run_noImpact_model = F
      
      progressSpecies$inc(0,detail=paste0(sp," (", match(sp,sel_sp),"/",length(sel_sp),")"))
      
      ### loading species data
      load(paste0("0.Data/processed_data_per_species/clean_data_",sp,".RData"))
      
      ### Run analyses
      
      ## calculate shortest path or euclidian distance between wind parcs and colonies
      print("Calculating distances between parcs and colonies:")
      
      # Distances colonies - parcs : on utilise la distance par la mer pour la plupart des espèces ; 
      # mais pour les espèces aux moeurs terrestres on utilise la distance euclidienne
      if(sp %in% moeurs_terrestres){
        
        colo_parc_distances <- Bdy_GetColoParcDistances(colonies=colonies,
                                                        parcs = parcs_L93,
                                                        costMatrix=trCost,
                                                        euclidian=T)
      }else{
        selec <- which(row.names(shortest_path_distances) %in% colonies$code_colonie)
        colo_parc_distances <- shortest_path_distances[selec,,drop=F]
      }
      
      ## Apportionning
      print("Apportionning")
      # Load sea area
      #(Attention : info non disponible pour la Guifette noire)
      if(sp != "Guifette noire") load(paste0("0.Data/sea_area/sea_area_", sp, ".rda"))
      
      apportionning <- Bdy_apportionning(species = sp,
                                         max_foraging_range_km = foraging_range_tbl$Max_km[which(foraging_range_tbl$Espece==sp)], #B
                                         colonies = colonies,
                                         sea_area=sea_area,
                                         tbl_dist=colo_parc_distances)  
      
      ## extract collision risks from mortality data
      print("Processing mortality")
      
      #processing the mortality data for the species of interest
      #collisionRisks <- read.csv2("../Test_data_Shiny/04_CollisionRisk_Cormoran huppé.csv") #for testing
      morta_iter_group <- Bdy_processing_mortality(
        collision = FormattedMortality()[which(FormattedMortality()$Espece==sp),],
        #collision = collisionRisks, #for testing
        season=seasons[which(seasons$Espèce==sp),][-1],
        n_iteration=input$iteration_nbr,
        parcNames = Parc_names_mortality(),
        #parcNames = parcs_L93$NAME, # for testing
        RW_group = apportionning$RW_group
      )
      
      ## running model with no impact
      print("Running or loading null model")
      
      #progress bar for model with impacts
      progressNullModel <- shiny::Progress$new()
      progressNullModel$set(message = "Chargement du modèle nul", value = 0.3)
      
      if(run_noImpact_model){
        ## Making the model without collisions
        no_impact_output <- Bdy_model_no_impact(count_data=group_counts,
                                                PI=ppa_yr,
                                                survival=vital_rates[vital_rates$Espèce == sp, "Survie"],
                                                fecundity=vital_rates[vital_rates$Espèce == sp, "Fecund"],
                                                propRepro=prop_repro[vital_rates$Espèce == sp, "PropRepro"],
                                                modelFile="0.Data/model_01.txt",nimble=F,
                                                ni=input$iteration_nbr,)
      }else{
        #load("0.Data/model_output_no_impact/model_no_impact_Fou de Bassan.RData") #for testing
        load(paste0("0.Data/model_output_no_impact/model_no_impact_",sp,"_light.RData"))
      }
      progressNullModel$inc(0.7)
      Sys.sleep(1)
      progressNullModel$close()
      
      
      ## Runnning the model with collisions
      print("Running the model with collisions")
      
      #progress bar for model with impacts
      progressImpactModel <- shiny::Progress$new()
      progressImpactModel$set(message = "Modèle avec impact", value = 0)
      
      with_impact_output <- Bdy_model_with_impact(n_group=nlevels(colonies$group),
                                                  ny_data=ncol(group_counts),
                                                  posterior=no_impact_output,
                                                  mortality=morta_iter_group,
                                                  ni=input$iteration_nbr,
                                                  progress=progressImpactModel)
      progressImpactModel$close()
      
      
      progressSpecies$inc(1/length(sel_sp),detail=paste0(sp," (", match(sp,sel_sp),"/",length(sel_sp),")"))
      
      model_output[[sp]] <- list(no_impact=with_impact_output$sc0, with_impact=with_impact_output$sc1, colonies=colonies, mortality=morta_iter_group)
      
    }) #future
    } #end species loop
    
    Sys.sleep(1)
    progressSpecies$close()
    
    print("Preparing figures and tables")
    
    progressResults <- shiny::Progress$new()
    progressResults$set(message = "Préparation des résultats", value = 0.1)
    
    ### Make sure to empty all results before storing figures (in case they were drawn previously)
    Stored_results$Res_list <- NULL 
    
    ### Prepare tables (Subpop and National)
    Raw_ResTables <- Bdy_RawResTables(model_output)
    Stored_results$Res_list$Raw_ResTables <- Raw_ResTables
    progressResults$set(0.2)
    
    Pretty_table_national <- Bdy_PrettyResultTable(Raw_ResTables, type="national")
    Stored_results$Res_list$DT_table_national <- Bdy_InteractiveResultTable(Pretty_table_national, type="national")
    progressResults$set(0.4)
    
    ### Prepare figures
    Summary_Figures<- Bdy_SummaryFigure(Raw_ResTables$Tableau_National)
    Stored_results$Res_list$Summary_Figure1 <- Summary_Figures$G_tendances
    Stored_results$Res_list$Summary_Figure2 <- Summary_Figures$G_extinction
    progressResults$set(0.6)
    
    Stored_results$Res_list$Summary_Map <- Bdy_SummaryMap(model_output, Raw_ResTables, Shapefile_stored$Shp_list$Merged_shp)
    Stored_results$Res_list$Text1 <- "Resultat 1" %>% HTML()
    progressResults$set(0.8)
    
    ### Prepare outputs
    Bdy_PrepareOutputs(FOLDER_name=FOLDER_name(),
                       mortality=FormattedMortality(),
                       parcs=parcs_L93,
                       distances=shortest_path_distances,
                       loadedDistances=LoadedDistances())
    progressResults$inc(1) 
    
    ### Close loader
    loader_analyses$hide()
    
    Sys.sleep(1)
    progressResults$close()
    
    Bdy_log("END - Run analyses")
  })
  
  
  ### Click + to get more results on a species
  observeEvent(input$select_button, {
    
    # Get species to detail
    SP <- input$select_button %>% sub("button_", "", .) %>% gsub("_", " ", .)
    Bdy_log(paste0("Detail button clicked for ", SP))
    
    # Prepare table
    Tab_SP <- Stored_results$Res_list$Raw_ResTables %>%
      Bdy_PrettyResultTable(., type="subpop") %>%
      subset(., Espece==SP) %>% 
      Bdy_InteractiveResultTable(., type="subpop")
    
    # Create popup
    showNotification(ui=renderDT({Tab_SP}, width="100%", height="auto"), type="default", duration=NULL)
  })
  
  
  ##### OUTPUTS -----
  output$Warning_mortality <- renderText({Warning_mortality() %>% paste0("<center><i>", ., "</center></i>") %>% HTML()})
  output$Warning_parc <- renderText({Warning_parc() %>% paste0("<center><i>", ., "</center></i>") %>% HTML()})
  output$Warning_general <- renderText({Warning_general() %>% paste0("<center><i>", ., "</center></i><br>") %>% HTML()})
  output$Warning_distances <- renderText({Warning_distances() %>% HTML()})
  
  output$Mortality_Data <- DT::renderDataTable({
    req(nrow(LoadedMortality())>0)
    #req("espece_latin" %in% names(LoadedMortality()) & "parc" %in% names(LoadedMortality()) & "espece" %in% names(LoadedMortality()) & "estimate" %in% names(LoadedMortality()) & "month" %in% names(LoadedMortality()))
    #DT::datatable(LoadedMortality() %>% mutate(espece_latin=as.factor(espece_latin), parc=as.factor(Parc), Espece=as.factor(Espece), coefficient=as.numeric(Estimate), mois=as.factor(Month)),
    DT::datatable(LoadedMortality(),
                  rownames=F,
                  filter="top",
                  options=list( headerCallback = DT::JS(# Needed to reduce column name font
                    "function(thead) {",
                    "  $(thead).css('font-size', '0.75em');",
                    "}"))
    ) %>% DT::formatStyle(columns = c(1:ncol(LoadedMortality())), fontSize = '75%')
  })
  
  output$TableNational <- renderDT({Stored_results$Res_list$DT_table_national})
  output$FigureSummary1 <- renderPlot({Stored_results$Res_list$Summary_Figure1})
  output$FigureSummary2 <- renderPlot({Stored_results$Res_list$Summary_Figure2})
  output$MapSummary <- renderLeaflet({Stored_results$Res_list$Summary_Map})
  
  output$downloadOutputs <- downloadHandler(
    filename = function() {paste0(FOLDER_name(), ".zip")},
    content = function(file) {
      file.copy(paste0("Outputs/", FOLDER_name(), ".zip"), file, overwrite = TRUE)
    }
  )
  
}



### RUN -------
shinyApp(ui = ui, server = server)




