server.R

    input$vr_mat_number_age_classes
  }, {
    req(input$vr_mat_number_age_classes)
    number_age_class <- input$vr_mat_number_age_classes
    updateMatrixInput(session, inputId = "mat_fill_vr",
                        value = matrix(data = NA,
                                       nrow = number_age_class,
                                       ncol = 2,
                                       dimnames = list(c(paste("Age", (1:number_age_class)-1)), c("Survie", "Fcondit"))))
  }) # end observeEvent



  # Update the vital rate matrix (mat_fill_vr) when changing species in the list
  observeEvent({
    input$species_choice
  }, {

    if(input$species_choice == "Espce gnrique") {

      number_age_class <- input$vr_mat_number_age_classes
      updateMatrixInput(session, inputId = "mat_fill_vr",
                        value = matrix(data = NA,
                                       nrow = number_age_class,
                                       ncol = 2,
                                       dimnames = list(c(paste("Age", (1:number_age_class)-1)), c("Survie", "Fcondit"))))
    } else {

      tab_species <- make_mat_vr(data_sf = data_sf, species = input$species_choice)

      if(all(is.na(tab_species))) {
        number_age_class <- input$vr_mat_number_age_classes
        updateMatrixInput(session, inputId = "mat_fill_vr",
                          value = matrix(data = NA,
                                         nrow = number_age_class,
                                         ncol = 2,
                                         dimnames = list(c(paste("Age", (1:number_age_class)-1)), c("Survie", "Fcondit"))))

      } else {
        number_age_class <- nrow(tab_species)
        ages <- tab_species$classes_age
        survivals <- tab_species$survie
        fecundities <- tab_species$fecondite

        updateMatrixInput(session, inputId = "mat_fill_vr",
                          value = matrix(data = c(survivals, fecundities),
                                         nrow = number_age_class,
                                         ncol = 2,
                                         dimnames = list(ages, c("Survie", "Fcondit"))))
      } # end if 2
    } # end if 1
  }) # end observeEvent species_choice


  # Display vital rates output table
  delay(ms = 300,
        output$vital_rates_info <- renderTable({

          #input$mat_fill_vr

          tab_species <- make_mat_vr(data_sf = data_sf, species = input$species_choice)
          ages <- tab_species$classes_age
          matrix(data = c(param$s_calib0, param$f_calib0),
                  nrow = length(param$s_calib0),
                  ncol = 2,
                  dimnames = list(ages, c("Survie", "Fcondit"))
                 )
        }, rownames = TRUE)
  )




  # Display intrinsic lambda (based solely on Leslie matrix)
  delay(ms = 300,
        output$lambda0_info <- renderText({
          req(all(!is.na(input$mat_fill_vr)))
          lam <- lambda(build_Leslie(s = param$s_calib0, f = param$f_calib0))
          taux <- round(lam-1,4)*100
          if(taux < 0) Text <- "Dclin : " else Text <- "Croissance : "
          if(taux == 0) Text <- "Population stable : "
          paste0(Text, taux, "% par an")
        })
  )


  #####

  #################################
  ## Dispersal
  ##-------------------------------
  observeEvent({
    input$species_choice
  }, {
    distAVG <- species_data %>%
      filter(NomEspece == input$species_choice) %>%
      select(DistDispMoyKM)

    rv$distAVG <- round(distAVG, 1)

    rv$dist <- c(round(-log(0.03)*distAVG, 1),
                 round(-log(0.05)*distAVG, 1),
                 round(-log(0.10)*distAVG, 1))
  })

  output$dispersal_mean_info <- renderText({
    paste0("Distance moyenne de dispersion : ", rv$distAVG, " km")
    })

  output$dispersal_d03p_info <- renderText({
    paste0("Seuil de distance quiv. 3% de dispersion : ", rv$dist[1], " km")
  })

  output$dispersal_d05p_info <- renderText({
    paste0("Seuil de distance quiv. 5% de dispersion : ", rv$dist[2], " km")
  })

  output$dispersal_d10p_info <- renderText({
    paste0("Seuil de distance quiv. 10% de dispersion : ", rv$dist[3], " km")
  })

  #####



  #####
  ##--------------------------------------------
  ## Select parameter values for simulations
  ##--------------------------------------------
  # Functions to calculate mean and SD from lower & upper values
  get_mu <- function(lower, upper) (lower + upper)/2
  get_sd <- function(lower, upper, coverage) ((abs(upper - lower)/2))/qnorm(1-((1-coverage)/2))

  #################################
  ## Cumulated impacts or not ?
  ##-------------------------------
  observeEvent({
    input$run
  }, {
    if(input$analysis_choice == "cumulated"){
      param$cumulated_impacts = TRUE
    } else {
      param$cumulated_impacts = FALSE
    } # end if
  }) # end observeEvent

  #################################
  ## Fatalities
  ##-------------------------------
  observe({
    # Case 1 : Not cumulated effects (if1)
    if(input$analysis_choice == "single_farm"){

      # Case 1.1 : Values from expert elicitation (if2)
      if(input$fatalities_input_type == "eli_exp"){
        if(!(is.null(param$fatalities_eli_result))){
          param$fatalities_mean <- c(0, round(param$fatalities_eli_result$mean, 2))
          param$onset_time <- NULL
          param$fatalities_se <- c(0, round(param$fatalities_eli_result$SE, 3))
          ready$fatalities <- TRUE
        } else {
          ready$fatalities <- FALSE
        }

      } else {

        if(input$fatalities_input_type == "val"){
          # Case 1.2 : Values directly provided as mean & SE
          param$fatalities_mean <- c(0, input$fatalities_mean)
          param$onset_time <- NULL
          param$fatalities_se <- c(0, input$fatalities_se)
          ready$fatalities <- TRUE

        }else{
          # Case 1.3 : Values directly provided as lower/upper interval
          param$fatalities_mean <- c(0, round(get_mu(lower = input$fatalities_lower, upper = input$fatalities_upper), 2))
          param$onset_time <- NULL
          param$fatalities_se <- c(0, round(get_sd(lower = input$fatalities_lower, upper = input$fatalities_upper, coverage = CP), 3))
          ready$fatalities <- TRUE
        } # end (if3)

      } # end (if2)

      # Case 2 : Cumulated effects
    } else {
      if(input$analysis_choice == "cumulated"){
        ready$fatalities <- TRUE
        param$fatalities_mean <- c(0, input$fatalities_mat_cumulated[,1])
        param$fatalities_se <- c(0, input$fatalities_mat_cumulated[,2])
        param$onset_year <- c(min(input$fatalities_mat_cumulated[,3]), input$fatalities_mat_cumulated[,3])
        param$onset_time <- param$onset_year - min(param$onset_year) + 1

        # Case 3 : Scenarios
      }else{
        req(input$fatalities_vec_scenario)
        vec01 <- as.numeric(unlist(strsplit(input$fatalities_vec_scenario, " ")))
        param$fatalities_mean <- c(0, vec01)
        param$fatalities_se <- rep(0, length(vec01)+1)
        param$onset_time <- NULL
        ready$fatalities <- TRUE
      }





    } # end (if1)

  }) # end observe
  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###

  #################################
  ## Population size
  ##-------------------------------
  observe({
    # Case 1 : Values from expert elicitation
    if(input$pop_size_input_type == "eli_exp"){
      if(!(is.null(param$pop_size_eli_result))){
        param$pop_size_mean <- round(param$pop_size_eli_result$mean, 1)
        param$pop_size_se <- round(param$pop_size_eli_result$SE, 1)
        ready$pop_size <- TRUE
      } else {
        ready$pop_size <- FALSE
      }

    } else {

      if(input$pop_size_input_type == "val"){
        # Case 2 : Values directly provided as mean & SE
        ready$pop_size <- TRUE
        param$pop_size_mean <- round(input$pop_size_mean, 1)
        param$pop_size_se <- round(input$pop_size_se, 1)

      }else{
        # Case 3 : Values directly provided as lower/upper interval
        ready$pop_size <- TRUE
        param$pop_size_mean <- round(get_mu(lower = input$pop_size_lower, upper = input$pop_size_upper), 1)
        param$pop_size_se <- round(get_sd(lower = input$pop_size_lower, upper = input$pop_size_upper, coverage = CP), 1)
      } # end (if3)

    }
    param$pop_size_unit <- input$pop_size_unit
  })


  #################################
  ## Population growth
  ##-------------------------------
  observe({
    # Case 1 : Values from expert elicitation
    if(input$pop_growth_input_type == "eli_exp"){
      if(!(is.null(param$pop_growth_eli_result))){
        param$pop_growth_mean <- round(param$pop_growth_eli_result$mean, 4)
        param$pop_growth_se <- round(param$pop_growth_eli_result$SE, 5)
        ready$pop_growth <- TRUE
      } else {
        ready$pop_growth <- FALSE
      }

    } else {

      # Case 2 : Trend information
      if(input$pop_growth_input_type == "trend"){
        ready$pop_growth <- TRUE

        if(input$pop_trend == "growth") {
          if(input$pop_trend_strength == "weak") {
            param$pop_growth_mean <- growth_weak
          } else if(input$pop_trend_strength == "average"){
            param$pop_growth_mean <- growth_average
          } else {
            param$pop_growth_mean <- growth_strong
          }
        } else if(input$pop_trend == "decline"){
          if(input$pop_trend_strength == "weak") {
            param$pop_growth_mean <- decline_weak
          } else if(input$pop_trend_strength == "average"){
            param$pop_growth_mean <- decline_average
          } else {
            param$pop_growth_mean <- decline_strong
          }
        } else {
          param$pop_growth_mean <- pop_stable
        }
        param$pop_growth_se <- trend_se


        # Case 3 : Values directly provided (i.e., not from expert elicitation)
      } else {

        if(input$pop_growth_input_type == "val"){
          # Case 2 : Values directly provided as mean & SE
          ready$pop_growth <- TRUE
          param$pop_growth_mean <- round(make_lambda(input$pop_growth_mean), 4)
          param$pop_growth_se <- round(input$pop_growth_se/100, 5)

        }else{
          # Case 3 : Values directly provided as lower/upper interval
          ready$pop_growth <- TRUE
          param$pop_growth_mean <- round(get_mu(lower = make_lambda(input$pop_growth_lower),
                                                          upper = make_lambda(input$pop_growth_upper)), 4)
          param$pop_growth_se <- round(get_sd(lower = make_lambda(input$pop_growth_lower),
                                              upper = make_lambda(input$pop_growth_upper), coverage = CP), 5)
        } # end (if3)

      }
    }
  })



  #################################
  ## Carrying capacity
  ##------------------------------
  observe({
    if(input$carrying_cap_input_type == "eli_exp"){
      if(!is.null(param$carrying_cap_eli_result)){
        param$carrying_capacity_mean <- round(param$carrying_cap_eli_result$mean)
        param$carrying_capacity_se <- round(param$carrying_cap_eli_result$SE, 1)
        ready$carrying_capacity <- TRUE
      } else {
        ready$carrying_capacity <- FALSE
      }

    } else {

      if(input$carrying_cap_input_type == "no_K"){
        ready$carrying_capacity <- TRUE
        param$carrying_capacity_mean <- max(param$pop_size_mean*100, 1e30) # use a very large K
        param$carrying_capacity_se <- 0

      }else{
        # values: mean and se
        if(input$carrying_cap_input_type == "val"){
          ready$carrying_capacity <- TRUE
          param$carrying_capacity_mean <- input$carrying_capacity_mean
          param$carrying_capacity_se <- input$carrying_capacity_se

        }else{
          # lower/upper interval
          ready$carrying_capacity <- TRUE
          param$carrying_capacity_mean <- round(get_mu(lower = input$carrying_capacity_lower, upper = input$carrying_capacity_upper), 0)
          param$carrying_capacity_se <- round(get_sd(lower = input$carrying_capacity_lower, upper = input$carrying_capacity_upper, coverage = CP), 1)
        } # end (if3)

      }
    }
  })
  #############################################
  ## Survivals, fecundities
  ##-------------------------------------------
  observe({
    param$survivals <- input$mat_fill_vr[,1]
    param$fecundities <- input$mat_fill_vr[,2]

    # for now, until calibration is really done
    param$s_calib0 <- param$survivals
    param$f_calib0 <- param$fecundities

  }) # end observeEvent
  #####

  #############################################
  ## Calibration of survivals & fecundities
  ##-------------------------------------------
  ## Calibration 1 : just for information display
  observeEvent({
    input$button_calibrate_vr
  },{

    print(param$pop_growth_mean)

    vr_calib0 <- calibrate_params(
      inits = init_calib(s = param$survivals, f = param$fecundities, lam0 = param$pop_growth_mean),
      f = param$fecundities, s = param$survivals, lam0 = param$pop_growth_mean
    )
    param$s_calib0 <- head(vr_calib0, length(param$survivals))
    param$f_calib0 <- tail(vr_calib0, length(param$fecundities))
  })

  ## Calibration 2 : for simulation run
  observeEvent({
    input$run
  },{

    # We also define rMAX and theta here
    rMAX_species <- rMAX_spp(surv = tail(param$survivals,1), afr = min(which(param$fecundities != 0)))
    param$rMAX_species <- rMAX_species

    param$pop_growth_mean_use <- round(min(1 + rMAX_species, param$pop_growth_mean), 4)

    param$theta <- fixed_theta
    #param$theta <- theta_spp(rMAX_species)

    param$vr_calibrated <- calibrate_params(
      inits = init_calib(s = param$survivals, f = param$fecundities, lam0 = param$pop_growth_mean_use),
      f = param$fecundities, s = param$survivals, lam0 = param$pop_growth_mean_use
    )
    param$s_calibrated <- head(param$vr_calibrated, length(param$survivals))
    param$f_calibrated <- tail(param$vr_calibrated, length(param$fecundities))
  })
  #####

  ############################################################
  ## Convert Fatalities as numbers (not rates)
  ##----------------------------------------------------------
  # Make sure fatalities are expressed as "number" (not rate) for the run_simul function
  se_prod2 <- function(mu1, se1, mu2, se2) sqrt((se1^2 * se2^2) + (se1^2 * mu2^2) + (mu1^2 * se2^2))

  observeEvent({
    input$run
  },{
    if(input$fatalities_unit == "h"){
      pop_size_tot <- sum(pop_vector(pop_size = param$pop_size_mean, pop_size_type = param$pop_size_type, s = param$s_calibrated, f = param$f_calibrated)[-1])
      param$fatalities_mean_nb <- (param$fatalities_mean/100) * pop_size_tot
      param$fatalities_se_nb <- se_prod2(mu1 = param$fatalities_mean/100,
                                         se1 = param$fatalities_se/100,
                                         mu2 = pop_size_tot,
                                         se2 = (pop_size_tot/param$pop_size_mean) * param$pop_size_se)
    }else{
      param$fatalities_mean_nb <- param$fatalities_mean
      param$fatalities_se_nb <- param$fatalities_se
    }
  })

  ############################################################
  ## Observe parameter values to be used in simulations run
  ##----------------------------------------------------------
  observeEvent({
    input$run
  }, {
    param$nsim <- input$nsim
    param$fatal_constant <- input$fatalities_unit
    param$time_horizon <- input$time_horizon

    # This condition is used to avoid wild population swings in fast-paced species
    if(max(param$fecundities) < 1.5){
      param$coeff_var_environ <- coeff_var_environ
    }else{
      param$coeff_var_environ <- 0
    }

  }) # Close observEvent


  observe ({
    param # to ensure up-to-date values are run
  }) # end observe
  #####


  ## Function to translate time units in french
  units_time_french <- function(u){
    if(u == "secs")  u_fr <- "secondes"
    if(u == "mins")  u_fr <- "minutes"
    if(u == "hours") u_fr <- "heures"
    if(u == "days")  u_fr <- "jours"
    if(u == "weeks") u_fr <- "semaines"
    return(u_fr)
  }

  #####
  ##-----------------------------------------------------------------------------------
  ##                                RUN SIMULATIONS
  ##-----------------------------------------------------------------------------------
  observeEvent({
    input$run
  }, {

    print(param$pop_growth_mean_use)

    if(ready$fatalities & ready$pop_size & ready$pop_growth & ready$carrying_capacity){

      out$analysis_choice <- input$analysis_choice
      start_time <- Sys.time()

      withProgress(message = 'Simulation progress', value = 0, {

        out$run <- run_simul_shiny(nsim = param$nsim,
                                   cumulated_impacts = param$cumulated_impacts,

                                   fatalities_mean = param$fatalities_mean_nb,
                                   fatalities_se = param$fatalities_se_nb,
                                   onset_time = param$onset_time,

                                   pop_size_mean = param$pop_size_mean,
                                   pop_size_se = param$pop_size_se,
                                   pop_size_type = param$pop_size_unit,

                                   pop_growth_mean = param$pop_growth_mean_use,
                                   pop_growth_se = param$pop_growth_se,

                                   survivals = param$s_calibrated,
                                   fecundities = param$f_calibrated,

                                   carrying_capacity_mean = param$carrying_capacity_mean,
                                   carrying_capacity_se = param$carrying_capacity_se,

                                   theta = param$theta,
                                   rMAX_species = param$rMAX_species,

                                   model_demo = NULL,
                                   time_horizon = param$time_horizon,
                                   coeff_var_environ = param$coeff_var_environ,
                                   fatal_constant = param$fatal_constant)
      }) # Close withProgress

      end_time <- Sys.time()
      duration <- end_time - start_time
      out$run_time <- paste(round(as.numeric(duration), 2), units_time_french(units(duration)))
      print(out$run_time)


    }else{
      out$run <- NULL
      out$msg <- "error_not_ready"
    }
  }) # Close observEvent
  #####



  #####
  ##-----------------------------------------------------------------------------------
  ##                                OUTPUTS
  ##-----------------------------------------------------------------------------------

  ### Run time
  output$run_time <- renderText({
    req(input$run > 0)
    paste("Temps de calcul (simulations) :", out$run_time)
  })


  #######################################################################
  ## Impact : individual farms (for "cumulated impact" analysis only)
  ##---------------------------------------------------------------------
  print_indiv_impact <- function(){
    req(out$run)
    res <- get_metrics(N = out$run$N, cumulated_impacts = TRUE)
    n_farm <- (dim(res$indiv_farm$impact)[3]-1)
    fil <- paste0(round(t(res$indiv_farm$impact[param$time_horizon, -2, -1]),2)*100, "%")
    matrix(fil,
           nrow = n_farm,
           dimnames = list(paste("Parc",1:n_farm), c("Impact", "IC (min)", "IC (max)"))
    )
  } # end function print_impact

  # Display title
  output$title_indiv_impact_result <- renderText({
    req(input$run > 0, out$analysis_choice == "cumulated")
    paste("Rsultat : Impact de chaque parc olien, estim au bout de" , param$time_horizon, "ans")
  })

  # Display impact result (table)
  output$indiv_impact_table <- renderTable({
    req(input$run & out$analysis_choice == "cumulated")
    print_indiv_impact()
  }, rownames = TRUE)


  ##################################################
  ## Impact : GLOBAL (for all types of analysis)
  ##------------------------------------------------
  print_impact <- function(){
    req(out$run)
    res <- get_metrics(N = out$run$N, cumulated_impacts = FALSE)
    n_scen <- (dim(res$scenario$impact)[3]-1)

    RowNam <- NULL
    if(out$analysis_choice == "single_farm") RowNam <- c("Parc 1")
    if(out$analysis_choice == "cumulated") RowNam <- c("Parc 1", paste("... + Parc", (2:n_scen)))
    if(out$analysis_choice == "multi_scenario") RowNam <- paste("Scenario", (1:n_scen))

    fil <- paste0(round(t(res$scenario$impact[param$time_horizon, -2, -1]),2)*100, "%")
    matrix(fil,
           nrow = n_scen,
           dimnames = list(RowNam, c("Impact", "IC (min)", "IC (max)"))
    )
  } # end function print_impact

  # Display title
  output$title_impact_result <- renderText({
    req(input$run)
    paste("Rsultat : Impact global estim au bout de" , param$time_horizon, "ans")
  })

  # Display impact result (table)
  output$impact_table <- renderTable({
    req(input$run)
    print_impact()
  }, rownames = TRUE)


  #############################################
  ## Probability of extinction
  ##-------------------------------------------
  print_PrExt <- function(){
    req(out$run)
    res <- get_metrics(N = out$run$N, cumulated_impacts = FALSE)
    n_scen <- dim(res$scenario$impact)[3]

    RowNam <- NULL
    if(out$analysis_choice == "single_farm") RowNam <- c("Sans parc", "Avec parc")
    if(out$analysis_choice == "cumulated") RowNam <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
    if(out$analysis_choice == "multi_scenario") RowNam <- paste("Scenario", (1:n_scen)-1)

    fil <- paste0(round(t(res$scenario$Pext),2)*100, "%")
    matrix(fil,
           nrow = n_scen,
           dimnames = list(RowNam, c("Probabilit d'extinction"))
    )
  } # end function print_PrExt

  # Display title
  output$title_PrExt_result <- renderText({
    req(input$run)
    paste("Rsultat : Probabilit d'extinction ", param$time_horizon, "ans")
  })

  # Display impact result (table)
  output$PrExt_table <- renderTable({
    req(input$run)
    print_PrExt()
  }, rownames = TRUE)


  #############################################
  ## Plot Impacts
  ##-------------------------------------------
  ## Function to plot the impact
  plot_out_impact <- function(){
    if(is.null(out$run)) {} else {

      n_scen <- dim(out$run$N)[3]
      Legend <- NULL
      if(out$analysis_choice == "single_farm") Legend <- c("Sans parc", "Avec parc")
      if(out$analysis_choice == "cumulated") Legend <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
      if(out$analysis_choice == "multi_scenario") Legend <- paste("Scenario", (1:n_scen)-1)

      plot_impact(N = out$run$N, onset_year = param$onset_year, percent = TRUE,
                  xlab = "\nAnne", ylab = "Impact relatif (%)\n", Legend = Legend)
      }
  }

  output$title_impact_plot <- renderText({
    if(input$run > 0){
      "Rsultat : Impact relatif au cours du temps"
    }
  })

  output$impact_plot <- renderPlot({
    plot_out_impact()
  })

  #############################################
  ## Plot Demographic Trajectories
  ##-------------------------------------------
  # Function to plot trajectories
  plot_out_traj <- function(){
    if(is.null(out$run)) {
    } else {

      n_scen <- dim(out$run$N)[3]

      # Define Legend
      Legend <- NULL
      if(out$analysis_choice == "single_farm") Legend <- c("Sans parc", "Avec parc")
      if(out$analysis_choice == "cumulated") Legend <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
      if(out$analysis_choice == "multi_scenario") Legend <- paste("Scenario", (1:n_scen)-1)

      # Plot population trajectories
      plot_traj(N = out$run$N, age_class_use = input$age_class_show, fecundities = param$f_calibrated, onset_year = param$onset_year,
                xlab = "\nAnne", ylab = "Taille de population\n", Legend = Legend, ylim = c(0, NA))}
  } # End function

  output$title_traj_plot <- renderText({
    if(input$run > 0){
      "Graphique : Trajectoires dmographiques"
    }
  })


  output$warning_traj_plot <- renderText({
    if(input$run > 0){
      "Attention : Il s'agit de prdictions en l'tat actuel des connaissances.
      Personne ne peut prdire comment les facteurs d'influence dmographique (environnement, etc.)
      vont voluer dans le futur. Donc personne ne peut prdire de faon exacte ce que sera la taille
      de population dans plusieurs annes.\n
      Ce graphe est simplementfourni  titre informatif. Attention  ne pas le sur-interprter.\n
      L'impact des collisions doit tre valu  partir des valeurs et du graphe ci-dessus, qui fournissent
      une estimation plus robuste (c--d. moins sensibles aux postulats et incertitudes) de cet impact."
    }
  })


  output$traj_plot <- renderPlot({
    plot_out_traj()
  })
  #####


  ###################################################################################
} # End server