server.R

        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
  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###


  # 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$survivals, f = param$fecundities)[-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
    }
  })


  #################################
  ## 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)
        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 <- input$pop_size_mean
        param$pop_size_se <- input$pop_size_se

      }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), 2)
        param$pop_size_se <- round(get_sd(lower = input$pop_size_lower, upper = input$pop_size_upper, coverage = CP), 3)
      } # 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(min(1 + param$rMAX_species, round(param$pop_growth_eli_result$mean, 2)), 2)
        param$pop_growth_se <- round(param$pop_growth_eli_result$SE, 2)
        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 <- 1.01
          } else if(input$pop_trend_strength == "average"){
            param$pop_growth_mean <- 1.03
          } else {
            param$pop_growth_mean <- 1.06
          }
        } else if(input$pop_trend == "decline"){
          if(input$pop_trend_strength == "weak") {
            param$pop_growth_mean <- 0.99
          } else if(input$pop_trend_strength == "average"){
            param$pop_growth_mean <- 0.97
          } else {
            param$pop_growth_mean <- 0.94
          }
        } else {
          param$pop_growth_mean <- 1
        }
        param$pop_growth_se <- 0


        # 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(min(1 + param$rMAX_species, make_lambda(input$pop_growth_mean)), 3)
          param$pop_growth_se <- input$pop_growth_se/100

        }else{
          # Case 3 : Values directly provided as lower/upper interval
          ready$pop_growth <- TRUE
          param$pop_growth_mean <- round(min(1 + param$rMAX_species,
                                             round(get_mu(lower = make_lambda(input$pop_growth_lower),
                                                          upper = make_lambda(input$pop_growth_upper)), 2)
                                             ), 3)
          param$pop_growth_se <- round(get_sd(lower = make_lambda(input$pop_growth_lower),
                                              upper = make_lambda(input$pop_growth_upper), coverage = CP), 3)
        } # 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 == "unknown"){
        ready$carrying_capacity <- TRUE
        param$carrying_capacity_mean <- max(param$pop_size_mean*100, 1e8) # use a very large K
        param$carrying_capacity_se <- 0

      }else{
        # values: mean and se
        if(input$pop_size_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), 2)
          param$carrying_capacity_se <- round(get_sd(lower = input$carrying_capacity_lower, upper = input$carrying_capacity_upper, coverage = CP), 3)
        } # end (if3)

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

  #############################################
  ## Calibration of survivals & fecundities
  ##-------------------------------------------
  observeEvent({
    input$run
  },{

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

    param$vr_calibrated <- 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_calibrated <- head(param$vr_calibrated, length(param$survivals))
    param$f_calibrated <- tail(param$vr_calibrated, length(param$fecundities))
  })
  #####

  ############################################################
  ## 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

  }) # Close observEvent


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

    # fixed in global environment (for now)
    param$theta <- theta
    param$coeff_var_environ <- coeff_var_environ

  }) # end observe
  #####

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

    if(ready$fatalities & ready$pop_size & ready$pop_growth & ready$carrying_capacity){
      out$analysis_choice <- input$analysis_choice
      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,
                                   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

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



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

  #######################################################################
  ## 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")
    "Rsultat : Impact de chaque parc olien, estim au bout de 30 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)
    "Rsultat : Impact global estim au bout de 30 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)
    "Rsultat : Probabilit d'extinction  30 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]
      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_traj(N = out$run$N, onset_year = param$onset_year,
                xlab = "\nAnne", ylab = "Taille de population\n", Legend = Legend)}
  }

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

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


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