Skip to content
Snippets Groups Projects
server.R 61.1 KiB
Newer Older
Marie-Bocage's avatar
Marie-Bocage committed
server <- function(input, output, session){
Marie-Bocage's avatar
Marie-Bocage committed


  ###################################################
  ##  Fixed parameters in the server environment
  ##-------------------------------------------------
  ## Load species list
  species_data <- read.csv("./inst/ShinyApp/species_list.csv", sep = ",")

  ## Load survival and fecundities data
  data_sf <- read.csv("./inst/ShinyApp/survivals_fecundities_species.csv", sep = ",")#, encoding = "UTF-8")

  # We define theta = 1 (same as in PBR) - for simplicity, given large uncertainty of real shape of density-dependence in nature
  fixed_theta = 1

  # Coefficient of environmental variation (SD)
  ## Environnmental variance set at 8%, based on values found for birds in the literature:
  ## (Saeher & Engen 2002) : between 7% et 14 ==> average : 10%
  ## (Sther et al. 2005) : between 2.5% et 10% ==> average : 6%
  coeff_var_environ = sqrt(0.08) # SD ~28%

  # Coverage probability used for lower/upper interval input values
  CP = 0.99

  # Values of pop_growth (assumed), when the "trend" option is chosen
  growth_weak <- 1.05
  growth_average <- 1.10
  growth_strong <- 1.15

  decline_weak <- 0.97
  decline_average <- 0.94
  decline_strong <- 0.91

  pop_stable <- 1
  trend_se <- 0.05 # SE to use for pop_growth, when the "trend" option is chosen


  ##############################################
thierrychambert's avatar
thierrychambert committed
  ##--------------------------------------------

  ## Fatalities
  output$hide_fatalities <- eventReactive({
    input$button_fatalities
  },{
    if(input$button_fatalities%%2 == 1) TRUE else FALSE
  }, ignoreInit = TRUE)

  outputOptions(output, "hide_fatalities", suspendWhenHidden = FALSE)


  ## Population Size
thierrychambert's avatar
thierrychambert committed
  output$hide_pop_size <- eventReactive({
    input$button_pop_size
  },{
    if(input$button_pop_size%%2 == 1) TRUE else FALSE
  }, ignoreInit = TRUE)

  outputOptions(output, "hide_pop_size", suspendWhenHidden = FALSE)


  ## Population Growth
  output$hide_pop_growth <- eventReactive({
    input$button_pop_growth
  },{
    if(input$button_pop_growth%%2 == 1) TRUE else FALSE
  }, ignoreInit = TRUE)

  outputOptions(output, "hide_pop_growth", suspendWhenHidden = FALSE)


  ## Carrying capacity
  output$hide_carrying_cap <- eventReactive({
    input$button_carrying_cap
  },{
    if(input$button_carrying_cap%%2 == 1) TRUE else FALSE
  }, ignoreInit = TRUE)

  outputOptions(output, "hide_carrying_cap", suspendWhenHidden = FALSE)

  # Display Carrying capacity Unit Info
  output$carrying_cap_unit_info <- renderText({
    if(input$pop_size_unit == "Npair"){
      paste0("Nombre de couple")
    } else {
      paste0("Effectif total")
    }
  })


  ## Outputs / Results
  output$hide_results <- eventReactive({
    input$run
  },{
    if(input$run > 0) TRUE else FALSE
  }, ignoreInit = TRUE)

  outputOptions(output, "hide_results", suspendWhenHidden = FALSE)


  ##############################################
thierrychambert's avatar
thierrychambert committed
  ##--------------------------------------------
thierrychambert's avatar
thierrychambert committed
  observe({
    shinyjs::hide("fatalities_mean")
    shinyjs::hide("fatalities_se")
    shinyjs::hide("fatalities_lower")
    shinyjs::hide("fatalities_upper")
    shinyjs::hide("fatalities_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("fatalities_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("fatalities_run_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("farm_number_cumulated")
    shinyjs::hide("fatalities_mat_cumulated")
    shinyjs::hide("fatalities_vec_scenario")
    shinyjs::hide("pop_size_lower")
    shinyjs::hide("pop_size_upper")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_size_mean")
    shinyjs::hide("pop_size_se")
    shinyjs::hide("pop_size_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_size_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_size_run_expert")
    shinyjs::hide("pop_growth_lower")
    shinyjs::hide("pop_growth_upper")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_growth_mean")
    shinyjs::hide("pop_growth_se")
    shinyjs::hide("pop_growth_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_growth_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_growth_run_expert")
Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("pop_trend")
    shinyjs::hide("pop_trend_strength")
    shinyjs::hide("carrying_capacity_lower")
    shinyjs::hide("carrying_capacity_upper")
    shinyjs::hide("carrying_capacity_mean")
    shinyjs::hide("carrying_capacity_se")
    shinyjs::hide("carrying_cap_number_expert")
    shinyjs::hide("carrying_cap_mat_expert")
    shinyjs::hide("carrying_cap_run_expert")

Marie-Bocage's avatar
Marie-Bocage committed
    shinyjs::hide("mat_fill_vr")
    shinyjs::hide("vr_mat_number_age_classes")
Marie-Bocage's avatar
Marie-Bocage committed

    #------------
    # Show some
    #------------
    # Show inputs for fatalities part
Marie-Bocage's avatar
Marie-Bocage committed
    if(input$button_fatalities%%2 == 1){
thierrychambert's avatar
thierrychambert committed
      #shinyjs::show("fatal_constant")
Marie-Bocage's avatar
Marie-Bocage committed

      # Show inputs for single farm option (non-cumulated impacts)
      if(input$analysis_choice == "single_farm"){
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("fatalities_input_type")
        if(input$fatalities_input_type == "itvl"){
          shinyjs::show("fatalities_lower")
          shinyjs::show("fatalities_upper")
Marie-Bocage's avatar
Marie-Bocage committed
        }
        if(input$fatalities_input_type == "val"){
          shinyjs::show("fatalities_mean")
          shinyjs::show("fatalities_se")
        }
thierrychambert's avatar
thierrychambert committed
        if(input$fatalities_input_type == "eli_exp"){
          shinyjs::show("fatalities_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
          shinyjs::show("fatalities_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
          shinyjs::show("fatalities_run_expert")
      # Show inputs for cumulated impacts option
Marie-Bocage's avatar
Marie-Bocage committed
      if(input$analysis_choice == "cumulated"){
        shinyjs::hide("fatalities_input_type")
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("farm_number_cumulated")
        shinyjs::show("fatalities_mat_cumulated")
thierrychambert's avatar
thierrychambert committed
      }
Marie-Bocage's avatar
Marie-Bocage committed

      # Show inputs for multiple scenario
      if(input$analysis_choice == "multi_scenario"){
        shinyjs::hide("fatalities_input_type")
        shinyjs::show("fatalities_vec_scenario")
      }

Marie-Bocage's avatar
Marie-Bocage committed
    # Show inputs for population size part
    if(input$button_pop_size%%2 == 1){
      shinyjs::show("pop_size_input_type")
      if(input$pop_size_input_type == "itvl"){
        shinyjs::show("pop_size_lower")
        shinyjs::show("pop_size_upper")
      }
thierrychambert's avatar
thierrychambert committed
      if(input$pop_size_input_type == "val"){
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_size_mean")
        shinyjs::show("pop_size_se")
      }
thierrychambert's avatar
thierrychambert committed
      if(input$pop_size_input_type == "eli_exp"){
        shinyjs::show("pop_size_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_size_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_size_run_expert")
    # Show inputs for population trend/growth part
    if(input$button_pop_growth%%2 == 1){
      shinyjs::show("pop_growth_input_type")

      if(input$pop_growth_input_type == "itvl"){
        shinyjs::show("pop_growth_lower")
        shinyjs::show("pop_growth_upper")
      }
thierrychambert's avatar
thierrychambert committed
      if(input$pop_growth_input_type == "val"){
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_growth_mean")
        shinyjs::show("pop_growth_se")
      }
thierrychambert's avatar
thierrychambert committed
      if(input$pop_growth_input_type == "eli_exp"){
        shinyjs::show("pop_growth_number_expert")
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_growth_mat_expert")
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_growth_run_expert")
Marie-Bocage's avatar
Marie-Bocage committed
      }
thierrychambert's avatar
thierrychambert committed
      if(input$pop_growth_input_type == "trend"){
Marie-Bocage's avatar
Marie-Bocage committed
        shinyjs::show("pop_trend")
        if(input$pop_trend != "stable"){
          shinyjs::show("pop_trend_strength")
        }
Marie-Bocage's avatar
Marie-Bocage committed
      }
thierrychambert's avatar
thierrychambert committed
    }

    # Show inputs for carrying capacity part
    if(input$button_carrying_cap%%2 == 1){
      shinyjs::show("carrying_cap_input_type")
      if(input$carrying_cap_input_type == "itvl"){
        shinyjs::show("carrying_capacity_lower")
        shinyjs::show("carrying_capacity_upper")
      }
      if(input$carrying_cap_input_type == "val"){
        shinyjs::show("carrying_capacity_mean")
        shinyjs::show("carrying_capacity_se")
      }
      if(input$carrying_cap_input_type == "eli_exp"){
        shinyjs::show("carrying_cap_number_expert")
        shinyjs::show("carrying_cap_mat_expert")
        shinyjs::show("carrying_cap_run_expert")
      }
    }
Marie-Bocage's avatar
Marie-Bocage committed

Marie-Bocage's avatar
Marie-Bocage committed
    if(input$button_vital_rates%%2 == 1){
thierrychambert's avatar
thierrychambert committed
      shinyjs::show("mat_fill_vr")
      shinyjs::show("vr_mat_number_age_classes")
Marie-Bocage's avatar
Marie-Bocage committed
    }

    # Show radiobutton (output) for plot_traj graph
    if(input$run > 0){
      shinyjs::show("age_class_show")
    }

  }) # en observe show/hide
  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###
  ##############################################
thierrychambert's avatar
thierrychambert committed
  ##--------------------------------------------
  out <- reactiveValues(run = NULL, run_time = NULL, msg = NULL,
                        analysis_choice = NULL, analysis_choice_report = NULL,
                        species_choice = NULL,
                        fatalities_input_type = NULL,
                        fatalities_input_val1 = NULL,
                        fatalities_input_val2 = NULL,

                        trajectory_plot = NULL)
  rv <- reactiveValues(distAVG = NULL, dist = NULL)
  ready <- reactiveValues(fatalities = TRUE, pop_size = TRUE, pop_growth = TRUE, carrying_capacity = TRUE)
Marie-Bocage's avatar
Marie-Bocage committed
  param <- reactiveValues(N1 = NULL,
                          cumulated_impacts = FALSE,
Marie-Bocage's avatar
Marie-Bocage committed
                          fatalities_mean = NULL,
                          onset_time = NULL,
                          onset_year = NULL,
                          out_fatal = NULL,

                          pop_size_mean = NULL,
                          pop_size_se = NULL,
thierrychambert's avatar
thierrychambert committed
                          pop_size_unit = NULL,
                          pop_growth_mean_use = NULL,
Marie-Bocage's avatar
Marie-Bocage committed
                          fecundities = NULL,
                          survivals = NULL,
Marie-Bocage's avatar
Marie-Bocage committed
                          s_calibrated = NULL,
                          f_calibrated = NULL,
                          vr_calibrated = NULL,
                          carrying_capacity_mean = NULL,
                          carrying_capacity_se = NULL,


thierrychambert's avatar
thierrychambert committed
                          rMAX_species = NULL,
                          time_horizon = NULL,
                          coeff_var_environ = NULL,
                          fatal_constant = NULL,

Marie-Bocage's avatar
Marie-Bocage committed
                          fatalities_eli_result = NULL,
                          pop_size_eli_result = NULL,
                          pop_growth_eli_result = NULL,
                          carrying_cap_eli_result = NULL
  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###


  ##############################################
  ##-------------------------------------------
  # Get lambda from +/-X% growth rate
  make_lambda <- function(pop_growth)  1 + (pop_growth/100)
  #####

  #####
  ##------------------------------------------
  ## Update elicitation matrices
  ##------------------------------------------

  ###############################
  ## Cumulated Impacts Matrix
  ##-----------------------------
  observeEvent({
    input$farm_number_cumulated
  }, {
    req(input$farm_number_cumulated > 0)
    current_mat <- input$fatalities_mat_cumulated
    n_farm <- input$farm_number_cumulated
    if(n_farm > nrow(current_mat)){
      fill_mat <- c(as.vector(t(current_mat)), rep(NA,(3*(n_farm-nrow(current_mat)))))
    }else{
      fill_mat <- as.vector(t(current_mat[1:n_farm,]))
    }
    updateMatrixInput(session, inputId = "fatalities_mat_cumulated",
                      value =  matrix(fill_mat, nrow = n_farm, ncol = 3, byrow = TRUE,
                                      dimnames = list(paste("Parc", c(1:n_farm)),
                                                      c("Moyenne",
                                                        "Erreur-type",
Marie-Bocage's avatar
Marie-Bocage committed

  ########################
  ## Fatalities Matrix
  ##----------------------
  observeEvent({
    input$fatalities_number_expert
  }, {
    req(input$fatalities_number_expert > 0)
    current_mat <- input$fatalities_mat_expert
    n_experts <- input$fatalities_number_expert
    if(n_experts > nrow(current_mat)){
      fill_mat <- c(as.vector(t(current_mat)), rep(NA,(5*(n_experts-nrow(current_mat)))))
    }else{
      fill_mat <- as.vector(t(current_mat[1:n_experts,]))
    }
    updateMatrixInput(session, inputId = "fatalities_mat_expert",
                      value = matrix(fill_mat, nrow = n_experts, ncol = 5, byrow = TRUE,
                                     dimnames = list(paste0("#", 1:n_experts),
                                                     c("Poids", "Min", "Best", "Max", "% IC" ))
                                     )
                      )
  })
  #####

  ########################
  ## Pop Size Matrix
  ##----------------------
  observeEvent({
    input$pop_size_number_expert
  }, {
    req(input$pop_size_number_expert > 0)
    current_mat <- input$pop_size_mat_expert
    n_experts <- input$pop_size_number_expert
    if(n_experts > nrow(current_mat)){
      fill_mat <- c(as.vector(t(current_mat)), rep(NA,(5*(n_experts-nrow(current_mat)))))
    }else{
      fill_mat <- as.vector(t(current_mat[1:n_experts,]))
    }
    updateMatrixInput(session, inputId = "pop_size_mat_expert",
                      value = matrix(fill_mat, nrow = n_experts, ncol = 5, byrow = TRUE,
                                     dimnames = list(paste0("#", 1:n_experts),
                                                     c("Poids", "Min", "Best", "Max", "% IC" ))
                      )
    )
  })
  #####

  ########################
  ## Pop Growth Matrix
  ##----------------------
  observeEvent({
    input$pop_growth_number_expert
  }, {
    req(input$pop_growth_number_expert > 0)
    current_mat <- input$pop_growth_mat_expert
    n_experts <- input$pop_growth_number_expert
    if(n_experts > nrow(current_mat)){
      fill_mat <- c(as.vector(t(current_mat)), rep(NA,(5*(n_experts-nrow(current_mat)))))
    }else{
      fill_mat <- as.vector(t(current_mat[1:n_experts,]))
    }
    updateMatrixInput(session, inputId = "pop_growth_mat_expert",
                      value = matrix(fill_mat, nrow = n_experts, ncol = 5, byrow = TRUE,
                                     dimnames = list(paste0("#", 1:n_experts),
                                                     c("Poids", "Min", "Best", "Max", "% IC" ))
                      )
    )
  })
  #####

  ############################
  ## Carrying Capacity Matrix
  ##--------------------------
  observeEvent({
    input$carrying_cap_number_expert
  }, {
    req(input$carrying_cap_number_expert > 0)
    current_mat <- input$carrying_cap_mat_expert
    n_experts <- input$carrying_cap_number_expert
    if(n_experts > nrow(current_mat)){
      fill_mat <- c(as.vector(t(current_mat)), rep(NA,(5*(n_experts-nrow(current_mat)))))
    }else{
      fill_mat <- as.vector(t(current_mat[1:n_experts,]))
    }
    updateMatrixInput(session, inputId = "carrying_cap_mat_expert",
                      value = matrix(fill_mat, nrow = n_experts, ncol = 5, byrow = TRUE,
                                     dimnames = list(paste0("#", 1:n_experts),
                                                     c("Poids", "Min", "Best", "Max", "% IC" ))
                      )
    )
  })
  #####

  ##--------------------------------------------
  ##  Run expert elicitation
  ##--------------------------------------------
  # Function to run the elication analysis
  func_eli <- function(mat_expert){
    t_mat_expert <- t(mat_expert)
    vals <- t_mat_expert[2:4,] %>% apply(., 2, sort)
    Cp <- t_mat_expert[5,] %>% sapply(., min, 0.99) %>% sapply(., max, 0.2)
    weights <- t_mat_expert[1,]

    out <- tryCatch(
      elicitation(vals, Cp, weights), error = function(e)
        return(NULL)
        #message("Erreur : certaines valeurs dans la matrice d'experts n'ont pas de sens")
      )

    if(!is.null(out)){
      OUT <- list(out = out, mean = out$mean_smooth, SE = sqrt(out$var_smooth))
    }else{
  }

  # Function to plot the elication analysis output
  plot_expert <- function(out, show_se = TRUE, ...){
    plot_elicitation(out, ylab = "", xlab = "Valeur du paramtre", cex.lab = 1.2, yaxt = "n")
    mtext(text = "Densit de probabilit", side = 2, line = 2, cex = 1.2)

    y2 <- dgamma(x = out$mean_smooth, shape = out$shape_smooth, rate = out$rate_smooth)
    xx <- qgamma(p = c(0.01,0.99), shape = out$shape_smooth, rate = out$rate_smooth)
    clip(xx[1], xx[2], -100, y2)
    abline(v = out$mean_smooth, lwd = 3, col = "darkblue")

    mtext(text = paste("Moyenne = ", round(out$mean_smooth,2)), side = 3, line = 2.5, cex = 1.2, adj = 0)
    if(show_se) mtext(text = paste("Erreur-type = ", round(sqrt(out$var_smooth), 2)), side = 3, line = 1, cex = 1.2, adj = 0)
  }

  ########################
  ## Fatalities
  ##----------------------
  observeEvent({
    input$fatalities_run_expert
  }, {
    if( all(!is.na(input$fatalities_mat_expert)) ) {
Marie-Bocage's avatar
Marie-Bocage committed

      ## run elicitation analysis
Marie-Bocage's avatar
Marie-Bocage committed
      param$fatalities_eli_result <- func_eli(input$fatalities_mat_expert)

      ## plot distribution
      output$title_distri_plot <- renderText({ "Mortalits annuelles" })
      output$distri_plot <- renderPlot({ plot_expert(param$fatalities_eli_result$out) })

    } else {
      print("missing value")
    } # end if
  }) # end observeEvent
Marie-Bocage's avatar
Marie-Bocage committed

  ########################
  ## Population size
  ##----------------------
  observeEvent({
    input$pop_size_run_expert
  }, {
    if(all(!is.na(input$pop_size_mat_expert))) {

      ## run elicitation analysis
Marie-Bocage's avatar
Marie-Bocage committed
      param$pop_size_eli_result <- func_eli(input$pop_size_mat_expert)

      ## plot distribution
      output$title_distri_plot <- renderText({ "Taille de population" })
      output$distri_plot <- renderPlot({ plot_expert(param$pop_size_eli_result$out) })

    } else {
      print("missing value")
    } # end if
  }) # end observeEvent
Marie-Bocage's avatar
Marie-Bocage committed

  ########################
  ## Population growth
  ##----------------------
  observeEvent({
    input$pop_growth_run_expert
  },{
    if(all(!is.na(input$pop_growth_mat_expert))){

      lambda_mat_expert <- input$pop_growth_mat_expert
      lambda_mat_expert[,2:4] <- make_lambda(lambda_mat_expert[,2:4])

      ## run elicitation analysis
      param$pop_growth_eli_result <- func_eli(lambda_mat_expert)
Marie-Bocage's avatar
Marie-Bocage committed

      ## plot distribution
      output$title_distri_plot <- renderText({ "Taux de croissance de la population" })
      output$distri_plot <- renderPlot({ plot_expert(param$pop_growth_eli_result$out) })
Marie-Bocage's avatar
Marie-Bocage committed

    } else {
      print("missing value")
    } # end if
  }) # end observeEvent

  ########################
  ## Carrying capacity
  ##----------------------
  observeEvent({
    input$carrying_cap_run_expert
  },{
    if(all(!is.na(input$carrying_cap_mat_expert))) {
Marie-Bocage's avatar
Marie-Bocage committed

Marie-Bocage's avatar
Marie-Bocage committed
      param$carrying_cap_eli_result <- func_eli(input$carrying_cap_mat_expert)

thierrychambert's avatar
thierrychambert committed
      ## show output
      if(!is.na(param$carrying_cap_eli_result$out)){
        output$title_distri_plot <- renderText({ "Capacit de charge" })
        output$distri_plot <- renderPlot({ plot_expert(param$carrying_cap_eli_result$out, show_se = FALSE) })
      }else {
        output$title_distri_plot <- renderText({ "Erreur : certaines valeurs dans la matrice d'experts n'ont pas de sens" })
thierrychambert's avatar
thierrychambert committed

    } else {
thierrychambert's avatar
thierrychambert committed
      param$carrying_cap_eli_result <- NULL
      print("missing value")
      output$title_distri_plot <- renderText({ "Des valeurs sont manquantes dans la table 'experts'" })
    } # end if
  }) # end observeEvent
thierrychambert's avatar
thierrychambert committed

  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###

  #####
  ##--------------------------------------------
  ##  Display parameter distribution
  ##--------------------------------------------
thierrychambert's avatar
thierrychambert committed

  # Function to plot a gamma distribution
  plot_gamma <- function(mu, se, show_mode = TRUE, show_mean = TRUE, show_se = TRUE, ...){
thierrychambert's avatar
thierrychambert committed

    ## Define shape and scale parameter of gamma distribution
    shape = (mu/se)^2
    scale = se^2/mu
thierrychambert's avatar
thierrychambert committed

    ## Plot the curve
    par(mar = c(5, 4, 6, 2))
    curve(dgamma(x, shape=shape, scale=scale), from = max(0,mu-3*se), to = mu+3*se, lwd = 3, col = "darkblue", yaxt = "n",
          ylab = "", xlab = "Valeur du paramtre", cex.lab = 1.2)
    mtext(text = "Densit de probabilit", side = 2, line = 2, cex = 1.2)
    # show mode
    MU <- (shape-1)*scale
    y_MU <- dgamma(x = MU, shape = shape, scale = scale)
    xx <- qgamma(p = c(0.01,0.99), shape = shape, scale = scale)
    clip(xx[1], xx[2], -100, y_MU)
    abline(v = MU, lwd = 3, col = "darkblue")

    # show mean
    y_mu <- dgamma(x = mu, shape = shape, scale = scale)
    clip(xx[1], xx[2], -100, y_mu)
    abline(v = mu, lwd = 2, col = "darkblue", lty = 2)
    if(show_mode) mtext(text = paste("Mode = ", round(MU, 2)), side = 3, line = 4, cex = 1.2, adj = 0)
    if(show_mean) mtext(text = paste("Moyenne = ", round(mu, 2)), side = 3, line = 2.5, cex = 1.2, adj = 0)
    if(show_se) mtext(text = paste("Erreur-type = ", round(se, 3)), side = 3, line = 1, cex = 1.2, adj = 0)
  } # end function plot_gamma

  plot_gamma_cumulated_impacts <- function(mu, se, nparc, ...){
    ## Define shape and scale parameter of gamma distribution
    shape = (mu/se)^2
    scale = se^2/mu

    ## Define x and y lim
    xx = yy = list()
    for(j in 1:nparc){
      xx[[j]] = seq(from = max(0,mu[j]-4*se[j]), to = mu[j]+4*se[j], length.out = 1e3)
      yy[[j]] = dgamma(xx[[j]], shape=shape[j], scale=scale[j])
    }

    ylim = c(min(unlist(yy)), max(unlist(yy))*1.4)
    xlim = c(min(unlist(xx)), max(unlist(xx)))

    ## Plot
    j=1
    curve(dgamma(x, shape=shape[j], scale=scale[j]),
          from = max(0,mu[j]-4*se[j]), to = mu[j]+4*se[j], n = 1e4,
          xlim = xlim, ylim = ylim,
          lwd = 3, col = j, yaxt = "n", xaxt = "n",
          #xaxp = c(round(xlim[1]), round(xlim[2]), n = 10),
          ylab = "", xlab = "Valeur du paramtre", cex.lab = 1.2)
    axis(side = 1, at = seq(round(xlim[1]), round(xlim[2]),
                            by = max(round((round(xlim[2])-round(xlim[1]))/10),1) ))
    mtext(text = "Densit de probabilit", side = 2, line = 2, cex = 1.2)

    y1 <- dgamma(x = mu[j], shape = shape[j], scale = scale[j])
    segments(x0 = mu[j], y0 = 0, y1 = y1, lty = 2, lwd = 3, col = j)
    points(x = mu[j], y = y1, pch = 19, cex = 1.5, col = j)

    for(j in 2:nparc){
      curve(dgamma(x, shape=shape[j], scale=scale[j]),
            from = max(0,mu[j]-4*se[j]), to = mu[j]+4*se[j], n = 1e4,
            lwd = 3, col = j, yaxt = "n",
            ylab = "", xlab = "Valeur du paramtre", cex.lab = 1.2, add = TRUE)

      y1 <- dgamma(x = mu[j], shape = shape[j], scale = scale[j])
      segments(x0 = mu[j], y0 = 0, y1 = y1, lty = 2, lwd = 3, col = j)
      points(x = mu[j], y = y1, pch = 19, cex = 1.5, col = j)
    }

    legend(x = xlim[1], y = ylim[2], legend = paste("Parc", 1:nparc),
           lwd = 3, col = 1:nparc, text.col = 1:nparc, cex = 1.5,
           bty = "n", horiz = TRUE)
  } # end function plot_gamma_cumulated_impacts
  ########################
  ## Fatalities
  ##----------------------
  observeEvent({
    input$analysis_choice
    input$button_fatalities
    input$fatalities_input_type
    input$fatalities_run_expert
    input$farm_number_cumulated
    input$fatalities_mat_cumulated
  },{
    ## 1. When analysis = single farm
    if(input$analysis_choice == "single_farm"){
      # Show from input values: if button is ON and input_type is set on "value" or "itvl" (thus not "eli_exp")
      if(input$button_fatalities%%2 == 1 & input$fatalities_input_type != "eli_exp"){
        output$title_distri_plot <- renderText({ "Mortalits annuelles" })

        output$distri_plot <- renderPlot({
          req(param$fatalities_mean, param$fatalities_se > 0)
          if(input$fatalities_input_type == "itvl"){
            req(input$fatalities_lower, input$fatalities_upper)
            plot_gamma(mu = tail(param$fatalities_mean, -1), se = tail(param$fatalities_se, -1))
          }else{
            req(input$fatalities_mean, input$fatalities_se)
            plot_gamma(mu = tail(param$fatalities_mean, -1), se = tail(param$fatalities_se, -1))
          }
        })

      } else {
        # Show from elicitation expert: if button is ON and input_type is set on "expert elicitation"
        if(input$button_fatalities%%2 == 1 & input$fatalities_input_type == "eli_exp"){
          if(!is.null(param$fatalities_eli_result)){
            output$title_distri_plot <- renderText({ "Mortalits annuelles" })
            output$distri_plot <- renderPlot({ plot_expert(param$fatalities_eli_result$out) })
          } else {
            output$title_distri_plot <- NULL
            output$distri_plot <- NULL
          }
          # Hide otherwise (when button is OFF)
        }else{
          output$title_distri_plot <- NULL
          output$distri_plot <- NULL
        }
      }
    ## 2. When analysis = cumulated impacts
      if(input$analysis_choice == "cumulated"){
        output$title_distri_plot <- renderText({ "Mortalits annuelles par parc (impacts cumuls)" })
        # Plot: note we use the "NULL + delay" sequence only to avoid error message in R console
        output$distri_plot <- NULL
        delay(5,
              output$distri_plot <- renderPlot({
                req(all(!is.na(input$fatalities_mat_cumulated[,1])), all(input$fatalities_mat_cumulated[,2] > 0))
                plot_gamma_cumulated_impacts(mu = input$fatalities_mat_cumulated[,1],
                                             se = input$fatalities_mat_cumulated[,2],
                                             nparc = input$farm_number_cumulated)
              })
        )
      }else{
        ## 3. When analysis = multi_scenarios
        output$title_distri_plot <- renderText({ "Pas de graphe (pas d'incertitude dans le cas 'mulitple scnarios')" })
        output$distri_plot <- NULL
      } # end "else"
  }, ignoreInit = FALSE)
Marie-Bocage's avatar
Marie-Bocage committed

  ########################
  ## Population size
  ##----------------------
  observeEvent({
    input$button_pop_size
    input$pop_size_input_type
  },{
    # Show from input values: if button is ON and input_type is set on "value"
    if(input$button_pop_size%%2 == 1 & input$pop_size_input_type != "eli_exp"){
      output$title_distri_plot <- renderText({ "Taille initiale de la population" })

      output$distri_plot <- renderPlot({
        req(param$pop_size_mean, param$pop_size_se > 0)
        plot_gamma(mu = param$pop_size_mean, se = param$pop_size_se)
      })

    } else {
      # Show from elicitation expert: if button is ON and input_type is set on "expert elicitation"
      if(input$button_pop_size%%2 == 1 & input$pop_size_input_type == "eli_exp"){
        if(!is.null(param$pop_size_eli_result)){
          output$title_distri_plot <- renderText({ "Taille initiale de la population" })
          output$distri_plot <- renderPlot({ plot_expert(param$pop_size_eli_result$out) })
        } else {
          output$title_distri_plot <- NULL
          output$distri_plot <- NULL
        }
        # Hide otherwise (when button is OFF)
      }else{
        output$title_distri_plot <- NULL
        output$distri_plot <- NULL
      }
  }, ignoreInit = FALSE)
Marie-Bocage's avatar
Marie-Bocage committed

  ########################
  ## Population growth
  ##----------------------
  observeEvent({
    input$pop_growth_input_type
    input$button_pop_growth
  },{

    # Show from input values: if button is ON and input_type is set on "value" or "interval"
    if(input$button_pop_growth%%2 == 1 & input$pop_growth_input_type != "eli_exp" & input$pop_growth_input_type != "trend"){
      output$title_distri_plot <- renderText({ "Taux de croissance de la population" })

      output$distri_plot <- renderPlot({
        req(param$pop_growth_mean, param$pop_growth_se > 0)
        plot_gamma(mu = param$pop_growth_mean, se = param$pop_growth_se)
      })

    } else {
      # Show from elicitation expert: if button is ON and input_type is set on "expert elicitation"
      if(input$button_pop_growth%%2 == 1 & input$pop_growth_input_type == "eli_exp"){
        if(!is.null(param$pop_growth_eli_result)){
          output$title_distri_plot <- renderText({ "Taux de croissance de la population" })
          output$distri_plot <- renderPlot({ plot_expert(param$pop_growth_eli_result$out) })
        } else {
          output$title_distri_plot <- NULL
          output$distri_plot <- NULL
        }
        # Hide otherwise (when button is OFF)
      }else{
        output$title_distri_plot <- NULL
        output$distri_plot <- NULL
      }
    }
  }, ignoreInit = FALSE)
  ########################
  ## Carrying capacity
  ##----------------------
  observeEvent({
    input$carrying_cap_input_type
    input$button_carrying_cap
  },{
thierrychambert's avatar
thierrychambert committed
    # Show from input values: if button is ON and input_type is set on "value"
    if(input$button_carrying_cap%%2 == 1 & input$carrying_cap_input_type != "eli_exp"){
      output$title_distri_plot <- renderText({ "Capacit de charge" })

      output$distri_plot <- renderPlot({
        req(param$carrying_capacity_mean, param$carrying_capacity_se > 0)
thierrychambert's avatar
thierrychambert committed
        plot_gamma(mu = param$carrying_capacity_mean, se = param$carrying_capacity_se)
      })

    } else {
      # Show from elicitation expert: if button is ON and input_type is set on "expert elicitation"
      if(input$button_carrying_cap%%2 == 1 & input$carrying_cap_input_type == "eli_exp"){
        if(!is.null(param$carrying_cap_eli_result)){
          if(!is.na(param$carrying_cap_eli_result$out)){
            output$title_distri_plot <- renderText({ "Capacit de charge" })
            output$distri_plot <- renderPlot({ plot_expert(param$carrying_cap_eli_result$out) })
          }else{
            output$title_distri_plot <- renderText({ "Erreur" })
            output$distri_plot <- NULL
          }
thierrychambert's avatar
thierrychambert committed
        } else {
          output$title_distri_plot <- NULL
          output$distri_plot <- NULL
        }
thierrychambert's avatar
thierrychambert committed
        # Hide otherwise (when button is OFF)
      }else{
        output$title_distri_plot <- NULL
        output$distri_plot <- NULL
      }
    }
  }, ignoreInit = FALSE)
  #####
thierrychambert's avatar
thierrychambert committed

Marie-Bocage's avatar
Marie-Bocage committed

  #####
  ##-------------------------------------------------
  ##  Display parameter values (on the side panel)
  ##-------------------------------------------------
  #################################
Marie-Bocage's avatar
Marie-Bocage committed
  ## Fatalities
  ##-------------------------------
  ## UNIT
  output$fatalities_unit_info <- renderText({
    if(!is.null(input$fatalities_unit)){
      if(input$fatalities_unit == "h"){
        paste0("Taux de mortalit")
      } else {
        paste0("Nombre de mortalits")
      }
    }
  })

  ## Values
  output$fatalities_mean_info <- renderText({
    if(input$fatalities_unit == "h") add_perc <- "%" else add_perc <- ""
    paste0(c("Moyenne : ",
             paste0(tail(param$fatalities_mean, -1), add_perc, collapse = ", ")
    ), collapse = "")
  })


  output$fatalities_se_info <- renderText({
    if(input$fatalities_unit == "h") add_perc <- "%" else add_perc <- ""
             paste0(tail(param$fatalities_se, -1), add_perc, collapse = ", ")
Marie-Bocage's avatar
Marie-Bocage committed

  #################################
  ## Poplutation size
  ##-------------------------------
  ## UNIT
  output$pop_size_unit_info <- renderText({
    if(!is.null(param$pop_size_unit)){
      if(param$pop_size_unit == "Npair"){
        paste0("Nombre de couple")
thierrychambert's avatar
thierrychambert committed
      } else {
thierrychambert's avatar
thierrychambert committed
      }
Marie-Bocage's avatar
Marie-Bocage committed
    }
Marie-Bocage's avatar
Marie-Bocage committed
  })

  output$pop_size_mean_info <- renderText({  paste0("Moyenne : ", param$pop_size_mean) })
  output$pop_size_se_info <- renderText({  paste0("Erreur-type : ", param$pop_size_se) })
Marie-Bocage's avatar
Marie-Bocage committed

  ## Show Popsize by age (table)
  # Function to create the table
  make_mat_popsizes <- function(data_sf, species, pop_size, pop_size_unit, survivals, fecundities){
    nam <- data_sf %>%
      filter(NomEspece == species) %>%
      select(classes_age) %>%

    matrix(round(pop_vector(pop_size = pop_size, pop_size_type = pop_size_unit, s = survivals, f = fecundities)),
           nrow = 1,
           dimnames = list("Effectifs", nam)
    )
  }

  # Display the table       (Note : the "delay" piece is just there to avoid an error message - time for parameters to be "loaded in")
        output$pop_size_by_age <- renderTable({
          if(any(is.na(param$survivals)) | any(is.na(param$fecundities))){
            matrix("Valeurs de survies et/ ou de fcondits manquantes",
                   nrow = 1, dimnames = list(NULL, "Erreur"))
          }else{
            make_mat_popsizes(data_sf = data_sf, species = input$species_choice, pop_size = param$pop_size_mean,
                              pop_size_unit = input$pop_size_unit, s = param$s_calib0, f = param$f_calib0)
          } # end if
        },
        width = "500px",
        rownames = FALSE,
        digits = 0)
    )

  #################################
  ## Population growth
  ##-------------------------------
  output$pop_growth_mean_info <- renderText({  paste0("Moyenne : ", param$pop_growth_mean) })
  output$pop_growth_se_info <- renderText({  paste0("Erreur-type : ", param$pop_growth_se) })
Marie-Bocage's avatar
Marie-Bocage committed

  #################################
  ##-------------------------------
  # UNIT (like pop size)
  ## UNIT
  output$carrying_capacity_unit_info <- renderText({
    if(!is.null(param$pop_size_unit)){
      if(input$carrying_cap_input_type == "no_K"){
        "Pas de capacit de charge (K = infini)"
      }else{
        if(param$pop_size_unit == "Npair"){
          paste0("Nombre de couple")
        } else {
          paste0("Effectif total")
        }
      }
    }
  })
  ## VALUES
  output$carrying_capacity_mean_info <- renderText({
    if(input$carrying_cap_input_type == "no_K"){
      paste0("Moyenne : ", param$carrying_capacity_mean)
  output$carrying_capacity_se_info <- renderText({
    if(input$carrying_cap_input_type == "no_K"){
      NULL
    }else{
      paste0("Erreur-type : ", param$carrying_capacity_se)
Marie-Bocage's avatar
Marie-Bocage committed
    }
  })
Marie-Bocage's avatar
Marie-Bocage committed

  #################################
Marie-Bocage's avatar
Marie-Bocage committed
  ## Vital rates
  ##-------------------------------
  # Function to create the matrix
  make_mat_vr <- function(data_sf, species){
    out_mat <- data_sf %>%
      filter(NomEspece == species) %>%
      select(classes_age, survie, fecondite)
    return(out_mat)
  }

  # Update the vital rate matrix (mat_fill_vr) when changing the number of age classes
  observeEvent({
    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"))))
  # 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"))))

      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,
                                         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)
thierrychambert's avatar
thierrychambert committed
        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
thierrychambert's avatar
thierrychambert committed
          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")
Marie-Bocage's avatar
Marie-Bocage committed

  #####
  ##--------------------------------------------
  ## 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"){
    } else {
      param$cumulated_impacts = FALSE
    } # end if
  }) # end observeEvent

  #################################
  ## Fatalities
  ##-------------------------------
    # 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$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
      # Case 2 : Cumulated effects
      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
      }





  ###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###

  #################################
  ## Population size
  ##-------------------------------
    # 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)
    }
    param$pop_size_unit <- input$pop_size_unit
  })


  #################################
  ## Population growth
  ##-------------------------------
    # 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") {
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- growth_weak
          } else if(input$pop_trend_strength == "average"){
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- growth_average
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- growth_strong
          }
        } else if(input$pop_trend == "decline"){
          if(input$pop_trend_strength == "weak") {
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- decline_weak
          } else if(input$pop_trend_strength == "average"){
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- decline_average
thierrychambert's avatar
thierrychambert committed
            param$pop_growth_mean <- decline_strong
thierrychambert's avatar
thierrychambert committed
          param$pop_growth_mean <- pop_stable
thierrychambert's avatar
thierrychambert committed
        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)
      }
    }
  })



  #################################
  ## Carrying capacity
  ##------------------------------
    if(input$carrying_cap_input_type == "eli_exp"){
      if(!is.null(param$carrying_cap_eli_result)){
thierrychambert's avatar
thierrychambert committed
        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
      }
      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

        # 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)
    }
  })
  #############################################
thierrychambert's avatar
thierrychambert committed
  ## 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
    input$button_calibrate_vr
    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
thierrychambert's avatar
thierrychambert committed
    #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$fatal_constant <- input$fatalities_unit
    param$time_horizon <- input$time_horizon

thierrychambert's avatar
thierrychambert committed
    # 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
thierrychambert's avatar
thierrychambert committed

  ## 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
  }, {
Marie-Bocage's avatar
Marie-Bocage committed

    if(ready$fatalities & ready$pop_size & ready$pop_growth & ready$carrying_capacity){
      out$analysis_choice <- input$analysis_choice
      out$species_choice <- input$species_choice
thierrychambert's avatar
thierrychambert committed
      start_time <- Sys.time()

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

        out$run <- run_simul_shiny(nsim = param$nsim,
                                   cumulated_impacts = param$cumulated_impacts,
Marie-Bocage's avatar
Marie-Bocage committed

                                   fatalities_mean = param$fatalities_mean_nb,
                                   fatalities_se = param$fatalities_se_nb,
                                   pop_size_mean = param$pop_size_mean,
                                   pop_size_se = param$pop_size_se,
                                   pop_size_type = param$pop_size_unit,
Marie-Bocage's avatar
Marie-Bocage committed

Marie-Bocage's avatar
Marie-Bocage committed

                                   survivals = param$s_calibrated,
                                   fecundities = param$f_calibrated,
Marie-Bocage's avatar
Marie-Bocage committed

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

Marie-Bocage's avatar
Marie-Bocage committed

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

thierrychambert's avatar
thierrychambert committed
      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)


  }) # Close observEvent
  #####
  #####
  ##-----------------------------------------------------------------------------------
  ##                                OUTPUTS
  ##-----------------------------------------------------------------------------------
Marie-Bocage's avatar
Marie-Bocage committed

thierrychambert's avatar
thierrychambert committed
  ### 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")
thierrychambert's avatar
thierrychambert committed
    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
Marie-Bocage's avatar
Marie-Bocage committed

  # Display title
  output$title_impact_result <- renderText({
thierrychambert's avatar
thierrychambert committed
    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({
thierrychambert's avatar
thierrychambert committed
    paste("Rsultat : Probabilit d'extinction ", param$time_horizon, "ans")
Marie-Bocage's avatar
Marie-Bocage committed

  # Display impact result (table)
  output$PrExt_table <- renderTable({
    req(input$run)
    print_PrExt()
  #############################################
  ## Plot Impacts
  ##-------------------------------------------
  ## Function to plot the impact
  plot_out_impact <- function(){
thierrychambert's avatar
thierrychambert committed
    if(is.null(out$run)) {} else {
thierrychambert's avatar
thierrychambert committed

      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)

thierrychambert's avatar
thierrychambert committed
      plot_impact(N = out$run$N, onset_year = param$onset_year, percent = TRUE,
thierrychambert's avatar
thierrychambert committed
                  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({
    out$impact_plot <- plot_out_impact()
    out$impact_plot
  #############################################
  ## Plot Demographic Trajectories
  ##-------------------------------------------
  # Function to plot trajectories
  plot_out_traj <- function(){
    if(is.null(out$run)) {
    } else {
thierrychambert's avatar
thierrychambert committed

      n_scen <- dim(out$run$N)[3]
thierrychambert's avatar
thierrychambert committed
      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))}
  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({
    out$trajectory_plot <- plot_out_traj()
    out$trajectory_plot

  #############################################
  ## Save outputs for report
  ##-------------------------------------------
  # Type d'analyse
  observeEvent({
    input$run
  }, {
    if(out$analysis_choice == "single_farm") out$analysis_choice_report <- "Impacts non cumuls"
    if(out$analysis_choice == "cumulated") out$analysis_choice_report <- "Impacts cumuls"
    if(out$analysis_choice == "multi_scenario") out$analysis_choice_report <- "Multiple scnarios"
  })

  observeEvent({
    input$run
  }, {
    #out$fatalities_input_type <- input$fatalities_input_type
    if(input$fatalities_input_type == "itvl"){
      out$fatalities_input_type <- "Saisie : intervalle\n"
      out$fatalities_val1 <- paste0("Min : ", input$fatalities_lower, " ; ")
      out$fatalities_val2 <- paste0("Max : ", input$fatalities_upper)
    }
  })


  #####
  ##-----------------------------------------------------------------------------------
  ##                                REPORT
  ##-----------------------------------------------------------------------------------
  output$report <- downloadHandler(

    filename = "RapportEolpopTEST001.pdf",

    content = function(file) {
      # Copy the report file to a temporary directory before processing it, in
      # case we don't have write permissions to the current working dir (which
      # can happen when deployed).
      tempReport <- file.path(tempdir(), "report.Rmd")

      file.copy("./inst/ShinyApp/report.Rmd", tempReport, overwrite = TRUE)

      # Set up parameters to pass to Rmd document
      paramsRMD <- list(
        intro = input$intro_report,
        analysis = out$analysis_choice_report,
        species = out$species_choice,

        fatalities_input_type = out$fatalities_input_type,
        fatalities_val1 = out$fatalities_val1,
        fatalities_val2 = out$fatalities_val2,

        impact_plot = out$impact_plot,
        trajectory_plot = out$trajectory_plot
        )


      # Knit the document, passing in the `params` list, and eval it in a
      # child of the global environment (this isolates the code in the document
      # from the code in this app).
      rmarkdown::render(tempReport, output_file = file,
                        params = paramsRMD,
                        envir = new.env(parent = globalenv())
      )
    }
  ) # close downloadHandler


thierrychambert's avatar
thierrychambert committed
  ###################################################################################
thierrychambert's avatar
thierrychambert committed