server.R

  )


  #####

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


  # 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, 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(min(1 + param$rMAX_species, round(param$pop_growth_eli_result$mean, 2)), 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(min(1 + param$rMAX_species, 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(min(1 + param$rMAX_species,
                                             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)){
        #if(param$carrying_cap_eli_result != "missing value"){
          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$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]
  }) # end observeEvent
  #####

  #############################################
  ## Calibration of survivals & fecundities
  ##-------------------------------------------
  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$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),
      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$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")
    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]
      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