diff --git a/DESCRIPTION b/DESCRIPTION
index b6e6e757e51dbc804f99e2957705fefc2d3fd6fc..e9b1450ef85a6858b4ae4c8fb9c698daa901847c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,5 +23,6 @@ Imports:
     popbio,
     RColorBrewer,
     SHELF,
+    shiny,
     stats,
     utils
diff --git a/NAMESPACE b/NAMESPACE
index 489003217c62862a66bfc971be54985f1a937e29..d0ab66e17e98e1f8739a88485f243a3f96ab7b36 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,6 +9,7 @@ export(calibrate_params)
 export(elicitation)
 export(get_metrics)
 export(infer_DD)
+export(infer_rMAX)
 export(init_calib)
 export(make_transparent)
 export(match_lam_delta)
@@ -18,7 +19,9 @@ export(plot_traj)
 export(pop_project)
 export(pop_project_cumulated_impacts)
 export(pop_vector)
+export(rMAX_spp)
 export(run_simul)
+export(run_simul_shiny)
 export(sample_elicitation)
 export(sample_gamma)
 import(RColorBrewer)
@@ -33,3 +36,4 @@ import(stats, except = c(filter, lag))
 
 import(utils)
 importFrom(dplyr,filter)
+importFrom(shiny,incProgress)
diff --git a/R/calibrate_params.R b/R/calibrate_params.R
index 559b1550af9ab9db48f9fab91274687b7a8e0845..ea57ed8c02586d12dd3e08b0ca7c51aee81ea1a4 100644
--- a/R/calibrate_params.R
+++ b/R/calibrate_params.R
@@ -40,9 +40,9 @@ calibrate_params <- function(inits = NULL, s, f, lam0){
   # Set parameter boundaries for the optimization
   if(lam0 - lam00 < 0){
     lower = c(rep(0, length(fu)), apply(cbind((s*0.5), 0.05), 1, max))
-    upper = c(fu, s)
+    upper = c(fu, s)*1.01
   }else{
-    lower = c(fu, s)
+    lower = c(fu, s)*0.99
     upper = c(rep(10, length(fu)), apply(cbind((s*1.25), 0.98), 1, min))
   }
 
@@ -53,6 +53,7 @@ calibrate_params <- function(inits = NULL, s, f, lam0){
   opt <- stats::optim(par = inits, fn = uti_fun, fo=fo, nac=nac, lam0=lam0,
                       lower = lower, upper = upper, method="L-BFGS-B")
 
+
   # Return output : New fecundity vector
   params <- c(tail(opt$par, nac), fo, head(opt$par, -nac))
   names(params) <- c(paste0("s", 1:nac), paste0("f", 1:nac))
diff --git a/R/infer_DD.R b/R/infer_DD.R
index 46a0176632a626e1100cb301ba371adb474b2413..fa9cba546a3ba93698f292c9a38d22b5fc674ec1 100644
--- a/R/infer_DD.R
+++ b/R/infer_DD.R
@@ -23,8 +23,9 @@
 infer_DD <- function(rMAX = NULL, K = NULL, theta = 1, pop_size_current, pop_growth_current){
 
   if(!is.null(rMAX)){
-    # Infer K
-    #rMAX = lambda_MAX - 1
+
+    ## Infer K
+        #rMAX = lambda_MAX - 1
     r_a = pop_growth_current - 1
     N_a = pop_size_current
 
@@ -33,12 +34,13 @@ infer_DD <- function(rMAX = NULL, K = NULL, theta = 1, pop_size_current, pop_gro
   }else{
 
     if(!is.null(K)){
-      # Infer rMAX
+
+      ## Infer rMAX
       N_a = pop_size_current
       r_a = pop_growth_current - 1
 
       rMAX <- r_a/((1-(N_a/K))^theta)
-      #lambda_MAX = rMAX + 1
+          #lambda_MAX = rMAX + 1
 
     }else{
 
diff --git a/R/infer_rMAX.R b/R/infer_rMAX.R
new file mode 100644
index 0000000000000000000000000000000000000000..b252ced908e991943c2871856cc854205f73c277
--- /dev/null
+++ b/R/infer_rMAX.R
@@ -0,0 +1,39 @@
+#' Infer the rMAX
+#'
+#' @description
+#' Calculates the rMAX of a density-dependence relationship, based on K and current population status (size and trend)
+#'
+#' @param K a strictly positive number. Carrying capacity (= maximum size that the population can ever reach)
+#' @param theta a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+#' The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+#' Note lambda = r + 1
+#' @param pop_size_current a strictly positive number. Current population size.
+#' @param pop_growth_current a strictly positive number. Current population growth rate.
+#' @param rMAX_theoretical the maximum value for rMAX, usually calculated from the function rMAX_spp. See ?rMAX_spp for details.
+#'
+#' @return the r_MAX value.
+#' @export
+#'
+#' @examples
+#'## Infer rMAX
+#'infer_rMAX(rMAX = NULL, K = 2000, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+#'
+#'
+infer_rMAX <- function(K, theta = 1, pop_size_current, pop_growth_current, rMAX_theoretical = Inf){
+
+      ## Infer rMAX
+      N_a = pop_size_current
+      r_a = pop_growth_current - 1
+
+      rMAX <- r_a/((1-(N_a/K))^theta)
+
+
+      ## Take the minimum between this value and the theoretical species rMAX
+      rMAX <- min(rMAX, rMAX_theoretical)
+
+  return(rMAX)
+
+}  # END FUNCTION
+################################################################################
+
+
diff --git a/R/models.R b/R/models.R
index 2f9e495f5ed709df33d56ddc6eb01056a89c0555..744b7c3cf370c74aae528c663ad8b1895e7137f9 100644
--- a/R/models.R
+++ b/R/models.R
@@ -26,6 +26,8 @@
 #'
 M1_noDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
+  ## M1_noDD_noDemoStoch
+
   # Build the LESLIE matrix
   A <- build_Leslie(s = s, f = f)
 
@@ -68,6 +70,8 @@ M1_noDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 #'
 M2_noDD_WithDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
+  ## M2_noDD_WithDemoStoch
+
   # Number of age classes
   nac = length(s)
 
@@ -126,6 +130,8 @@ M2_noDD_WithDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 #'
 M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params){
 
+  ## M3_WithDD_noDemoStoch
+
   # Extract DD parameters from list
   rMAX = DD_params$rMAX
   K = DD_params$K
@@ -143,6 +149,19 @@ M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params){
   s_Nt <- head(vr_Nt, length(s)) %>% sapply(min, 0.999)
   f_Nt <- tail(vr_Nt, length(f))
 
+
+  ## Check if approximation is close enough to desired lambda
+  if( abs((lambda(build_Leslie(s = s_Nt, f = f_Nt)) - lam_Nt) / lam_Nt) > 0.005 ){
+
+    # If difference is too large : Use optimisation function for better calibration
+    inits <- c(tail(vr_Nt, length(f)), head(vr_Nt, length(s)) %>% sapply(min, 0.999))
+    inits <- inits[inits != 0]
+    vr_calib <- calibrate_params(inits = inits, f = f_Nt, s = s_Nt, lam0 = lam_Nt)
+    s_Nt <- head(vr_calib, length(s_Nt))
+    f_Nt <- tail(vr_calib, length(f_Nt))
+
+  } # if
+
   # Build the LESLIE matrix
   A_Nt <- build_Leslie(s = s_Nt, f = f_Nt)
 
@@ -186,6 +205,8 @@ M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params){
 #'
 M4_WithDD_WithDemoStoch <- function(N1, s, f, h, DD_params){
 
+  ## M4_WithDD_WithDemoStoch
+
   # Extract DD parameters from list
   rMAX = DD_params$rMAX
   K = DD_params$K
@@ -203,6 +224,20 @@ M4_WithDD_WithDemoStoch <- function(N1, s, f, h, DD_params){
   s_Nt <- head(vr_Nt, length(s)) %>% sapply(min, 0.999)
   f_Nt <- tail(vr_Nt, length(f))
 
+
+  ## Check if approximation is close enough to desired lambda
+  if( abs((lambda(build_Leslie(s = s_Nt, f = f_Nt)) - lam_Nt) / lam_Nt) > 0.005 ){
+
+    # If difference is too large : Use optimisation function for better calibration
+    inits <- c(tail(vr_Nt, length(f)), head(vr_Nt, length(s)) %>% sapply(min, 0.999))
+    inits <- inits[inits != 0]
+    vr_calib <- calibrate_params(inits = inits, f = f_Nt, s = s_Nt, lam0 = lam_Nt)
+    s_Nt <- head(vr_calib, length(s_Nt))
+    f_Nt <- tail(vr_calib, length(f_Nt))
+
+  } # if
+
+
   # Number of age classes
   nac = length(s)
 
diff --git a/R/pop_project.R b/R/pop_project.R
index f731a1d0662550e3be677abf78f2eaebf3d0e71d..446e9a5a0c364d8ad18ca748b9f7507c73fd8728 100644
--- a/R/pop_project.R
+++ b/R/pop_project.R
@@ -55,6 +55,7 @@ pop_project <- function(fatalities,
   # Number of fatalities scenario
   nsc <- length(fatalities)
 
+
   # Initiate Pop Size (output) Array
   N <- array(NA, dim = c(nac, nyr, nsc), dimnames = list(paste0("age", 1:nac),
                                                          paste0("year", 1:nyr),
@@ -71,13 +72,14 @@ pop_project <- function(fatalities,
 
     # Fatalities : constant number (M) or constant rate (h)
     if(fatal_constant == "M"){
-      h <- M/apply(N[,t-1,], 2, sum)
+      h <- sapply(M/apply(N[,t-1,], 2, sum), min, 1)
     } else {
-      h <- M/apply(N[,1,], 2, sum)
+      h <- sapply(M/apply(N[,1,], 2, sum), min, 1)
     }
 
     # Sample a seed for RNG
-    seed <- runif(1, 0, 1e6)
+    seed <- ((((Sys.time() %>% as.numeric) %% 1e5) * 1e5) %% 1e5) %>% round
+      #runif(1, 0, 1e6)
 
     ## Projection : apply the LESLIE matrix calculation forward
     # Scenario 0
diff --git a/R/pop_project_cumulated_impacts.R b/R/pop_project_cumulated_impacts.R
index b8e042f0d43931dfc5058e85a58d11e9cfdd138a..0dcf41dad707d80280ba28e9788089b920875fc0 100644
--- a/R/pop_project_cumulated_impacts.R
+++ b/R/pop_project_cumulated_impacts.R
@@ -83,9 +83,9 @@ pop_project_cumulated_impacts <- function(fatalities,
 
     # Fatalities : constant number (M) or constant rate (h)
     if(fatal_constant == "M"){
-      h <- Mc[,t-1]/apply(N[,t-1,], 2, sum)
+      h <- sapply(Mc[,t-1]/apply(N[,t-1,], 2, sum), min, 1)
     } else {
-      h <- Mc[,t-1]/apply(N[,1,], 2, sum)
+      h <- sapply(Mc[,t-1]/apply(N[,1,], 2, sum), min, 1)
     }
 
     # Sample a seed for RNG
diff --git a/R/rMAX_spp.R b/R/rMAX_spp.R
new file mode 100644
index 0000000000000000000000000000000000000000..4dfbabcd7669d4fe7086a47b3eb433409907fb59
--- /dev/null
+++ b/R/rMAX_spp.R
@@ -0,0 +1,24 @@
+#' Rmax of a species
+#'
+#' @description
+#' Calculate the theoretical Rmax of a species using the equation of Niel and Lebreton 2005.
+#'
+#' References.
+#'
+#' Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+#' populations from incomplete data. Conservation Biology 19:826–835.
+#'
+#' @param surv Adult survival. Single value between 0 and 1.
+#' @param afr Age at first reproduction. Single integer value.
+#'
+#' @return the theoretical Rmax of the species
+#' @export
+#'
+#' @examples
+#' rMAX_spp(surv = 0.6, afr = 2)
+#'
+rMAX_spp <- function(surv, afr){
+  ( ((surv*afr - surv + afr + 1) + sqrt((surv - surv*afr - afr - 1)^2 - (4*surv*(afr)^2))) / (2*afr) ) - 1
+}
+
+
diff --git a/R/run_simul.R b/R/run_simul.R
index 9fd71801b6e577532ceebdeb1965b61f6eed696e..19b15e7488fbffa55b8527ece0eef1c129d16c88 100644
--- a/R/run_simul.R
+++ b/R/run_simul.R
@@ -16,7 +16,22 @@
 #' @param pop_growth_se Standard Error for population growth rate (= uncertainty around the value provided).
 #' @param survivals a vector. Average survival probabilities for each age class.
 #' @param fecundities a vector of fecundity values for each age class.
-#' @param DD_params NULL or a list. Density-dependence parameters (rMAX, K, theta). Only used in DD models M3 and M4.
+#'
+#' @param carrying_capacity a strictly positive number.
+#' Carrying capacity (= maximum size that the population can reach). Here, the unit is the same as pop_size_type.
+#' It can thus be expressed as the total population or the number of pair.
+#' @param theta a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+#' The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+#' Note lambda = r + 1
+#'
+#' @param rMAX_species the maximum value of rMAX for the species under consideration,
+#' usually calculated using the Niel & Lebreton (2005) equation.
+#' It can be calculated using the function rMAX_spp. See ?rMAX_spp for details.
+#'
+#' References :
+#' Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+#' populations from incomplete data. Conservation Biology 19:826–835.
+#'
 #' @param model_demo is NULL, by default, because the model choice will be made inside each iteration (simulation),
 #' base on the values of N0 and lam0 that are drawn.
 #' But it can be forced by setting the value, which must then be an R object corresponding to the demographic model to be used.
@@ -30,7 +45,6 @@
 #' each simulation iteration (dim 4)
 #' @export
 #'
-#'
 #' @examples
 #' fatalities_mean = c(0, 5, 10, 15, 20)
 #' fatalities_se = fatalities_mean*0.05
@@ -51,13 +65,17 @@
 #' coeff_var_environ = 0.10
 #' fatal_constant = "h"
 #'
-#' DD_params <- list(rMAX = NULL, K = 1200, theta = 1)
+#' carrying_capacity = 1200
+#' theta = 1
+#' rMAX_species <- 0.15
 #'
 #' run_simul(nsim = 10, cumuated_impacts = FALSE,
 #'            fatalities_mean, fatalities_se, onset_time = NULL,
 #'            pop_size_mean, pop_size_se, pop_size_type,
 #'            pop_growth_mean, pop_growth_se,
-#'            survivals, fecundities, DD_params = DD_params,
+#'            survivals, fecundities,
+#'            carrying_capacity, theta,
+#'            rMAX_species,
 #'            model_demo = NULL, time_horzion, coeff_var_environ, fatal_constant)
 #'
 #'
@@ -65,9 +83,22 @@ run_simul <- function(nsim, cumuated_impacts,
                       fatalities_mean, fatalities_se, onset_time,
                       pop_size_mean, pop_size_se, pop_size_type,
                       pop_growth_mean, pop_growth_se,
-                      survivals, fecundities, DD_params,
+                      survivals, fecundities,
+                      carrying_capacity, theta = 1, rMAX_species,
                       model_demo = NULL, time_horzion, coeff_var_environ, fatal_constant){
 
+
+  # Create object to store DD parameters
+  DD_params <- list()
+
+  # Define K
+  K <- sum(pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+
+  # Fill the list of DD parameters
+  DD_params$K <- NULL
+  DD_params$theta <- theta
+  DD_params$rMAX <- rMAX_species
+
   # Coefficient of variation for environment stochasticity
   cv_env <- coeff_var_environ
 
@@ -85,10 +116,15 @@ run_simul <- function(nsim, cumuated_impacts,
                                                                paste0("year", 1:nyr),
                                                                paste0("sc", (1:nsc)-1)
   ))
-  # object to store values of population growth drawn at each iteration
+  # Object to store values of population growth drawn at each iteration
   lam_it <- rep(NA, nsim)
 
+  # Time
   time_run <- system.time(
+
+    ##--------------------------------------------
+    # Start Loops over simulations              --
+    ##--------------------------------------------
     for(sim in 1:nsim){
 
       ## PARAMETER UNCERTAINTY : draw values for each input
@@ -103,6 +139,12 @@ run_simul <- function(nsim, cumuated_impacts,
         round %>%
         pop_vector(pop_size_type = pop_size_type, s = survivals, f = fecundities)
 
+      # Define K
+      K <- sum(pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+      if(K < sum(N0)) K <- round(sum(N0)*1.05)
+      DD_params$K <- K
+
+
       if(pop_growth_se > 0){
 
         # 3. Population Growth Rate
@@ -125,8 +167,9 @@ run_simul <- function(nsim, cumuated_impacts,
 
       } # End if/else
 
+model_demo = NULL
 
-      # Choose the model demographique to use (it choice was not forced)
+            # Choose the model demographique to use (if choice was not forced)
       if(is.null(model_demo)){
 
         ## Define the complete model by default
@@ -143,9 +186,9 @@ run_simul <- function(nsim, cumuated_impacts,
         if(lam_it[sim] > 1){
 
           # Extract rMAX
-          DD_params$rMAX <-
-            infer_DD(K = DD_params$K, theta = DD_params$theta,
-                     pop_size_current = sum(N0), pop_growth_current = lam_it[sim])$rMAX
+          DD_params$rMAX <- infer_rMAX(K = K, theta = theta,
+                     pop_size_current = sum(N0), pop_growth_current = lam_it[sim],
+                     rMAX_theoretical = rMAX_species)
 
           # ... and initially LARGE population
           if(sum(N0) > 500) model_demo <- M3_WithDD_noDemoStoch
@@ -173,7 +216,8 @@ run_simul <- function(nsim, cumuated_impacts,
                                s = s, f = f, DD_params = DD_params,
                                model_demo = model_demo, time_horzion = time_horzion,
                                coeff_var_environ = coeff_var_environ, fatal_constant = fatal_constant)
-    } # sim
+    } # sim ##-----------------------------------------------------------------------------------------
+
   ) # system.time
 
   return(list(time_run = time_run, N = N, lambdas = lam_it))
diff --git a/R/run_simul_shiny.R b/R/run_simul_shiny.R
new file mode 100644
index 0000000000000000000000000000000000000000..35ae8553ec4627c05ed85fcb98c394c7664703fe
--- /dev/null
+++ b/R/run_simul_shiny.R
@@ -0,0 +1,203 @@
+##==============================================================================
+##                     Function to run simulations                            ==
+##==============================================================================
+#' Run multiple population projections (simulations).
+#' Exactly the same function as run_simul, except that it includes a line of code to display
+#' the simulation progress bar in Shiny (incProgress function)
+#'
+#' @param nsim number of simulation
+#' @param cumuated_impacts Logical. If TRUE, we used the projection model for cumulated impacts.
+#' @param fatalities_mean a vector (numeric). Average number of fatalities, for each scenario.
+#' @param fatalities_se a vector (numeric). Standard Error for the number of fatalities, for each scenario (= uncertainties around the values provided).
+#' @param onset_time a vector (numeric). The times at which each wind farm fatality starts applying.
+#' @param pop_size_mean a single number. Average population size (either total population or number of pairs - see Ntype below).
+#' @param pop_size_se Standard Error for population size (= uncertainty around the value provided).
+#' @param pop_size_type character value indicating if the provided value pop_size correpsonds to Total Population Size ("Ntotal")
+#' or the Number of Pairs ("Npair"). A stable age distribution is used to infer the size of each age class.
+#' @param pop_growth_mean a number. Average population growth rate (lambda).
+#' @param pop_growth_se Standard Error for population growth rate (= uncertainty around the value provided).
+#' @param survivals a vector. Average survival probabilities for each age class.
+#' @param fecundities a vector of fecundity values for each age class.
+#'
+#' @param carrying_capacity a strictly positive number.
+#' Carrying capacity (= maximum size that the population can reach). Here, the unit is the same as pop_size_type.
+#' It can thus be expressed as the total population or the number of pair.
+#' @param theta a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+#' The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+#' Note lambda = r + 1
+#'
+#' @param rMAX_species the maximum value of rMAX for the species under consideration,
+#' usually calculated using the Niel & Lebreton (2005) equation.
+#' It can be calculated using the function rMAX_spp. See ?rMAX_spp for details.
+#'
+#' References :
+#' Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+#' populations from incomplete data. Conservation Biology 19:826–835.
+#'
+#' @param model_demo is NULL, by default, because the model choice will be made inside each iteration (simulation),
+#' base on the values of N0 and lam0 that are drawn.
+#' But it can be forced by setting the value, which must then be an R object corresponding to the demographic model to be used.
+#' The 4 possible models currently are: M1_noDD_noDemoStoch, M2_noDD_WithDemoStoch, M3_WithDD_noDemoStoch, M4_WithDD_WithDemoStoch,
+#' @param time_horzion a number. The number of years (time horizon) over which to project the population dynamics.
+#' @param coeff_var_environ a number. The coefficient of variation to model environment stochasticity.
+#' @param fatal_constant text (character). Either "h" or "M". Using "h" sets the fatality RATE as the constant value across years.
+#' Using "M" sets the NUMBER of fatalities as the constant value across years.
+#'
+#' @return a 4D array containing the size of each age class (dim 1), for each year (dim 2), each scenario (dim 3), and
+#' each simulation iteration (dim 4)
+#' @export
+#'
+#'
+#' @importFrom shiny incProgress
+#'
+
+run_simul_shiny <- function(nsim, cumuated_impacts,
+                      fatalities_mean, fatalities_se, onset_time,
+                      pop_size_mean, pop_size_se, pop_size_type,
+                      pop_growth_mean, pop_growth_se,
+                      survivals, fecundities,
+                      carrying_capacity, theta = 1, rMAX_species,
+                      model_demo = NULL, time_horzion, coeff_var_environ, fatal_constant){
+
+
+  # Create object to store DD parameters
+  DD_params <- list()
+
+  # Define K
+  K <- sum(pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+
+  # Fill the list of DD parameters
+  DD_params$K <- NULL
+  DD_params$theta <- theta
+  DD_params$rMAX <- rMAX_species
+
+  # Coefficient of variation for environment stochasticity
+  cv_env <- coeff_var_environ
+
+  # Number of years
+  nyr <- time_horzion
+
+  # Number of age classes
+  nac <- length(survivals)
+
+  # Number of fatalities scenario (+1 because we include a base scenario of NO fatality)
+  nsc <- length(fatalities_mean)
+
+  # Initiate Pop Size (output) Array
+  N <- array(NA, dim = c(nac, nyr, nsc, nsim), dimnames = list(paste0("age", 1:nac),
+                                                               paste0("year", 1:nyr),
+                                                               paste0("sc", (1:nsc)-1)
+  ))
+  # Object to store values of population growth drawn at each iteration
+  lam_it <- rep(NA, nsim)
+
+  # Time
+  time_run <- system.time(
+
+    ##--------------------------------------------
+    # Start Loops over simulations              --
+    ##--------------------------------------------
+    for(sim in 1:nsim){
+
+
+      # Increments to be shown in Shiny's Progess bar
+      shiny::incProgress(1/nsim, detail = paste("simulation", sim))
+      Sys.sleep(0.001)
+
+      ## PARAMETER UNCERTAINTY : draw values for each input
+      # 1. Nomber of fatalities
+      M <- NA
+      for(j in 1:nsc){
+        M[j] <- sample_gamma(1, mu = fatalities_mean[j], sd = fatalities_se[j])
+      }
+
+      # 2. Population size : draw and distribute by age class
+      N0 <- sample_gamma(1, mu = pop_size_mean, sd =  pop_size_se) %>%
+        round %>%
+        pop_vector(pop_size_type = pop_size_type, s = survivals, f = fecundities)
+
+      # Define K
+      K <- sum(pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+      if(K < sum(N0)) K <- round(sum(N0)*1.05)
+      DD_params$K <- K
+
+
+      if(pop_growth_se > 0){
+
+        # 3. Population Growth Rate
+        lam0 <- sample_gamma(1, mu = pop_growth_mean, sd = pop_growth_se)
+
+        # 4. Calibrate vital rates to match the the desired lambda
+        inits <- init_calib(s = survivals, f = fecundities, lam0 = lam0)
+
+        vr <- calibrate_params(inits = inits, f = fecundities, s = survivals, lam0 = lam0)
+        s <- head(vr, length(survivals))
+        f <- tail(vr, length(fecundities))
+        lam_it[sim] <- lambda(build_Leslie(s,f))
+
+      }else{
+
+        # No parameter uncertainty on population growth
+        s <- survivals
+        f <- fecundities
+        lam_it[sim] <- lambda(build_Leslie(s,f))
+
+      } # End if/else
+
+      model_demo = NULL
+
+      # Choose the model demographique to use (if choice was not forced)
+      if(is.null(model_demo)){
+
+        ## Define the complete model by default
+        model_demo <- M4_WithDD_WithDemoStoch
+
+        # DECLINING (or stable), but initially LARGE population
+        if(lam_it[sim] <= 1 & sum(N0) > 3000) model_demo <- M1_noDD_noDemoStoch
+
+        # DECLINING  (or stable), and initially SMALL population
+        if(lam_it[sim] <= 1 & sum(N0) <= 3000) model_demo <- M2_noDD_WithDemoStoch
+
+
+        # GROWING population...
+        if(lam_it[sim] > 1){
+
+          # Extract rMAX
+          DD_params$rMAX <- infer_rMAX(K = K, theta = theta,
+                                       pop_size_current = sum(N0), pop_growth_current = lam_it[sim],
+                                       rMAX_theoretical = rMAX_species)
+
+          # ... and initially LARGE population
+          if(sum(N0) > 500) model_demo <- M3_WithDD_noDemoStoch
+
+
+          # ... but initially SMALL population
+          if(sum(N0) <= 500) model_demo <- M4_WithDD_WithDemoStoch
+
+        } # if lam > 1
+
+
+      } # end if "is.null"
+
+
+      #
+      if(cumuated_impacts){
+        fun_project <- pop_project_cumulated_impacts
+      }else{
+        fun_project <- pop_project
+        onset_time = NULL
+      } # end if
+
+      # Project population trajectory
+      N[,,,sim] <- fun_project(fatalities = M, onset_time = onset_time, intial_pop_vector = N0,
+                               s = s, f = f, DD_params = DD_params,
+                               model_demo = model_demo, time_horzion = time_horzion,
+                               coeff_var_environ = coeff_var_environ, fatal_constant = fatal_constant)
+    } # sim ##-----------------------------------------------------------------------------------------
+
+  ) # system.time
+
+  return(list(time_run = time_run, N = N, lambdas = lam_it))
+
+} # End of function
+################################################################################
diff --git a/devtools_history.R b/devtools_history.R
index dfb98880e623b6b4f3bb0ba374b7115756c545a0..f4ee12f6c8b61a34f48764cf737e55b7c27d53f1 100644
--- a/devtools_history.R
+++ b/devtools_history.R
@@ -58,3 +58,6 @@ usethis::use_package("dplyr")
 usethis::use_package("ggplot2")
 
 
+# Add package in the DESCRIPTION file: Imports
+usethis::use_package("shiny")
+
diff --git a/inst/ShinyApp/expert_plot_piece.R b/inst/ShinyApp/expert_plot_piece.R
deleted file mode 100644
index 28d930f6e87fd312da953c353268b73fbb303281..0000000000000000000000000000000000000000
--- a/inst/ShinyApp/expert_plot_piece.R
+++ /dev/null
@@ -1,30 +0,0 @@
-
-
-# Expert stuff
-
-values <- reactiveValues(mat = NULL, t_mat = NULL, vals = NULL, Cp = NULL, weights = NULL, out = NULL)
-observe({
-  values$t_mat <- t(input$expert)
-})
-
-observe({
-  values$vals = values$t_mat[3:5,]
-  values$Cp = values$t_mat[6,]
-  values$weights = values$t_mat[2,]
-})
-
-func_out <- eventReactive({input$run_expert}, {if(all(is.na(values$t_mat))){print(values$vals)}
-  else {values$mat <- exp_estim_gamma(values$vals, values$Cp, values$weights)
-  SE <- sqrt(values$mat$var_smooth)
-  return(list(values$mat$mean_smooth, SE))}})
-
-output$Mean <- renderPrint({
-  func_out()
-})
-
-func_plot <- function(){if(is.null(values$mat)) {} else {plot_exp_estim(values$mat)}}
-output$plot <- renderPlot({
-  func_plot()
-})
-
-
diff --git a/inst/ShinyApp/f_output.R b/inst/ShinyApp/f_output.R
deleted file mode 100644
index 460821333bfe041b39bf86d2c233c9bd71663b57..0000000000000000000000000000000000000000
--- a/inst/ShinyApp/f_output.R
+++ /dev/null
@@ -1,37 +0,0 @@
-##==============================================================================
-##                           Plot trajectories                                ==
-##==============================================================================
-
-plot_it <- function(N){
-  TH <- dim(N)[2]
-
-  # Average trend
-  N_avg <- apply(N, c(1,2,3), mean)
-
-  # Color palette
-  col_sc0 <- colorRampPalette(colors = c("black", "grey"))(nsim) %>% make_transparent(., percent = 85)
-  col_h <- colorRampPalette(colors = c("red", "orange"))(nsim) %>% make_transparent(., percent = 85)
-
-  # Initiate plot
-  plot(x = 1:TH, y = colSums(N_avg[,,"sc0"]), 'l', col=1, lwd=3, ylim = c(0,max(colSums(N))),
-       xlab = "Annee", ylab = "Taille de population (totale)", main=NULL)
-
-  for(sim in 1:nsim) points(x = 1:TH, y = colSums(N[,,"sc0",sim]), 'l', col=col_sc0[sim])
-
-  for(sim in 1:nsim) points(x = 1:TH, y = colSums(N[,,"sc1",sim]), 'l', col=col_h[sim])
-
-  points(x = 1:TH, y = colSums(N_avg[,,"sc0"]), 'l', col=1, lwd=3)
-  points(x = 1:TH, y = colSums(N_avg[,,"sc1"]), 'l', col="red", lwd=3)
-} # End function
-
-
-
-##==============================================================================
-##                              Print impact                                  ==
-##==============================================================================
-print_it <- function(impact, lci, uci){
-  #print(
-    paste0("Impact sur la taille de population : ", round(impact, 2)*100, "%",
-         " [", round(lci, 2)*100, "% ; ", round(uci, 2)*100, "%]")
-  #)
-} # End function
diff --git a/inst/ShinyApp/param_fixes.R b/inst/ShinyApp/param_fixes.R
deleted file mode 100644
index 2fb0f238105d0fa4686522c11476bbd796f44311..0000000000000000000000000000000000000000
--- a/inst/ShinyApp/param_fixes.R
+++ /dev/null
@@ -1,48 +0,0 @@
-library(popbio)
-library(magrittr)
-library(RColorBrewer)
-
-##==============================================================================
-##                               Fixed inputs                                 ==
-##==============================================================================
-# Nb simulations
-nsim = 50
-
-# Time Horizon
-TH = 30
-
-# Bird fatalities for sc0
-M0 = 0
-M0_se = 0
-
-#M1 = 5
-#M1_se = 1
-
-# Population size
-#N00_mu = 200
-N00_se = 0
-
-# Population trend
-#lam0 = 1.01
-lam0_se = 0
-
-# Environmental stochasticity
-cv_env = 0
-
-# Number of age classes
-nac = 4
-
-# Survival of each age class
-s_input <- c(0.5, 0.7, 0.8, 0.95)
-# sm <- matrix(0, nrow=nac, ncol=1, dimnames = list(paste0("age", 1:nac), "Survival"))
-
-# Fecundity of each age class
-f_input <- c(0, 0, 0.05, 0.5)
-
-# Model
- #M2_noDD_WithDemoStoch
-model_demo = M1_noDD_noDemoStoch
-
-## Build the Leslie matrix
-#A <- build_Leslie(s = s, f = f)
-#round(lambda(A), 2)
diff --git a/inst/ShinyApp/print_out_piece.R b/inst/ShinyApp/print_out_piece.R
deleted file mode 100644
index a425d9dc323308961fd6c43c3cf375e944f2628f..0000000000000000000000000000000000000000
--- a/inst/ShinyApp/print_out_piece.R
+++ /dev/null
@@ -1,9 +0,0 @@
-print_out <- reactive({
-  if(is.null(out$N)){
-    "Waiting"
-  } else {
-    print_it(impact = get_metrics(N = out$N)[30,"avg","sc1"],
-             lci = get_metrics(N = out$N)[30,"lci","sc1"],
-             uci = get_metrics(N = out$N)[30,"uci","sc1"])
-  }
-}) # end reactive
diff --git a/inst/ShinyApp/runApp.R b/inst/ShinyApp/runApp.R
new file mode 100644
index 0000000000000000000000000000000000000000..cafa2b5e825a8220c8670225e8679b74e76d94a2
--- /dev/null
+++ b/inst/ShinyApp/runApp.R
@@ -0,0 +1,3 @@
+source("./inst/ShinyApp/ui.R")
+source("./inst/ShinyApp/server.R")
+shinyApp(ui = ui, server = server)
diff --git a/inst/ShinyApp/run_piece.R b/inst/ShinyApp/run_piece.R
deleted file mode 100644
index aad5607fd9d34ef0cb76596fca18fb65d6651acd..0000000000000000000000000000000000000000
--- a/inst/ShinyApp/run_piece.R
+++ /dev/null
@@ -1,31 +0,0 @@
-out$N <- run_simul(nsim = nsim,
-                   fatalities_mean = c(M0, input$M1),
-                   fatalities_se = c(M0_se, input$M1_se),
-                   pop_size_mean = input$N00_mu,
-                   pop_size_se = N00_se,
-                   N_type = input$N_type,
-                   pop_growth_mean = input$lam0,
-                   pop_growth_se = lam0_se,
-                   survivals_mean = s_input,
-                   fecundities_mean = f_input,
-                   model_demo = model_demo,
-                   time_horzion = TH,
-                   coeff_var_environ = cv_env,
-                   fatal_constant = input$mort_cons)
-
-
-
-out <- run_simul(nsim = nsim,
-                 fatalities_mean = c(M0, input$M1),
-                 fatalities_se = c(M0_se, M1_se),
-                 pop_size_mean = N00_mu,
-                 pop_size_se = N00_se,
-                 N_type = "Npair",
-                 pop_growth_mean = lam0,
-                 pop_growth_se = lam0_se,
-                 survivals_mean = s_input,
-                 fecundities_mean = f_input,
-                 model_demo = model_demo,
-                 time_horzion = TH,
-                 coeff_var_environ = cv_env,
-                 fatal_constant = "h")
diff --git a/junk.R b/junk.R
index ca62fdbc1808ebfd2a067c183f76b019414f5c40..25bc685eca65c156690a4b24d7e77a0e7c069a66 100644
--- a/junk.R
+++ b/junk.R
@@ -1,233 +1,284 @@
+rm(list = ls(all.names = TRUE))
+graphics.off()
+
+
+## Libraries
 library(eolpop)
 
-fatalities_mean = c(0, 5, 10, 15, 20)
-fatalities_se = fatalities_mean*0.05
+## Inputs
+nsim = 10
+
+fatalities_mean = c(0, 10, 5, 8)
+fatalities_se = c(0, 0.05, 0.05, 0.05)
 
 pop_size_mean = 200
-pop_size_se = 30
-pop_size_type = "Npair"
+pop_size_se = 25
 
 pop_growth_mean = 1
-pop_growth_se = 0.03
+pop_growth_se = 0
 
 survivals <- c(0.5, 0.7, 0.8, 0.95)
-
 fecundities <- c(0, 0, 0.05, 0.55)
 
-DD_params <- list(rMAX = NULL, K = 1200, theta = 1)
-
-time_horzion = 30
+model_demo = NULL # M2_noDD_WithDemoStoch #M1_noDD_noDemoStoch #M4_WithDD_WithDemoStoch #M3_WithDD_noDemoStoch #
+time_horzion = 50
 coeff_var_environ = 0.10
-fatal_constant = "h"
-
-run_test <- run_simul(nsim = 10, cumuated_impacts = FALSE,
-                      fatalities_mean, fatalities_se, onset_time = NULL,
-                      pop_size_mean, pop_size_se, pop_size_type,
-                      pop_growth_mean, pop_growth_se,
-                      survivals, fecundities, DD_params = DD_params,
-                      model_demo = NULL, time_horzion, coeff_var_environ, fatal_constant)
-
+fatal_constant = "M"
+pop_size_type = "Npair"
 
+cumuated_impacts = TRUE
 
-run_test
-run0 <- run_test
+onset_year = c(2010, 2013, 2016)
+onset_time = onset_year - min(onset_year) + 1
+onset_time = c(min(onset_time), onset_time)
 
+# Pop size total
+sum(pop_vector(pop_size = pop_size_mean, pop_size_type = pop_size_type, s = survivals, f = fecundities))
 
 
+# Define K
+carrying_capacity = 500
+theta = 1
+K = pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities) %>% sum
+K
 
+# Define theoretical rMAX for the species
+rMAX_species <- rMAX_spp(surv = tail(survivals,1), afr = min(which(fecundities != 0)))
+rMAX_species
 
+##  Avoid unrealistic scenarios
+pop_growth_mean <- min(1 + rMAX_species, pop_growth_mean)
+pop_growth_mean
 
 
+##--------------------------------------------
+##  Calibration                             --
+##--------------------------------------------
+# Calibrate vital rates to match the the desired lambda
+inits <- init_calib(s = survivals, f = fecundities, lam0 = pop_growth_mean)
+vr_calibrated <- calibrate_params(inits = inits, f = fecundities, s = survivals, lam0 = pop_growth_mean)
+s_calibrated <- head(vr_calibrated, length(survivals))
+f_calibrated <- tail(vr_calibrated, length(fecundities))
 
+build_Leslie(s = s_calibrated, f = f_calibrated) %>% lambda
 
+survivals = s_calibrated
+fecundities = f_calibrated
 
 
 
 
 
 
+#################################################################################################################
 
+##--------------------------------------------
+##  run simul                               --
+##--------------------------------------------
 
+# Create object to store DD parameters
+DD_params <- list()
 
+# Fill the list of DD parameters
+DD_params$K <- NULL
+DD_params$theta <- theta
+DD_params$rMAX <- rMAX_species
 
+# Coefficient of variation for environment stochasticity
+cv_env <- coeff_var_environ
 
+# Number of years
+nyr <- time_horzion
 
+# Number of age classes
+nac <- length(survivals)
 
+# Number of fatalities scenario (+1 because we include a base scenario of NO fatality)
+nsc <- length(fatalities_mean)
 
+# Initiate Pop Size (output) Array
+N <- array(NA, dim = c(nac, nyr, nsc, nsim), dimnames = list(paste0("age", 1:nac),
+                                                             paste0("year", 1:nyr),
+                                                             paste0("sc", (1:nsc)-1)
+))
+# object to store values of population growth drawn at each iteration
+lam_it <- rep(NA, nsim)
+sim=1
 
 
-library(dplyr)
-library(ggplot2)
+##########################
 
-# Get metrics and dimensions
-out <- get_metrics(N)$scenario$impact_sc
-TH <- dim(N)[2]
-nsc <- dim(N)[3]
+    ## PARAMETER UNCERTAINTY : draw values for each input
+    # 1. Nomber of fatalities
+    M <- NA
+    for(j in 1:nsc){
+        M[j] <- sample_gamma(1, mu = fatalities_mean[j], sd = fatalities_se[j])
+    }
 
-# Build dataframe
-df <- as.data.frame(cbind(year = 1:nrow(out), out[,,1], scenario = 1))
-for(j in 2:dim(out)[3]) df <- rbind(df, cbind(year = 1:nrow(out), out[,,j], scenario = j))
+    # 2. Population size : draw and distribute by age class
+    N0 <- sample_gamma(1, mu = pop_size_mean, sd =  pop_size_se) %>%
+      round %>%
+      pop_vector(pop_size_type = pop_size_type, s = survivals, f = fecundities)
 
-## Define Graphic Parameters
-size = 1.5
+    # Define K
+    K <- sum(pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+    if(K < sum(N0)) K <- round(sum(N0)*1.05)
+    DD_params$K <- K
 
-# Plot lines
-p <-
-  ggplot(data = df, aes(x = year, y = avg)) +
-  geom_line(data = filter(df, scenario > 1), size = size, aes(colour = factor(scenario))) +
-  geom_line(data = filter(df, scenario == 1), size = size, colour = "black")
 
-# Plot CIs
-p <- p + geom_ribbon(data = filter(df, scenario > 1),
-                     aes(ymin = uci, ymax = lci, fill = factor(scenario)), linetype = 0, alpha = 0.100)
 
-# Add legend
-Legend <- "parc"
-nsc <- max(df$scenario)-1
-p <- p + labs(x = "Annee", y = "Impact relatif",
-              col = "Scenario", fill = "Scenario") +
-  scale_color_hue(labels = paste(Legend, 1:nsc), aesthetics = c("colour", "fill"))
 
 
+    if(pop_growth_se > 0){
 
+      # 3. Population Growth Rate
+      lam0 <- sample_gamma(1, mu = pop_growth_mean, sd = pop_growth_se)
 
+      # 4. Calibrate vital rates to match the the desired lambda
+      inits <- init_calib(s = survivals, f = fecundities, lam0 = lam0)
 
+      vr <- calibrate_params(inits = inits, f = fecundities, s = survivals, lam0 = lam0)
+      s <- head(vr, length(survivals))
+      f <- tail(vr, length(fecundities))
+      lam_it[sim] <- lambda(build_Leslie(s,f))
 
+    }else{
 
+      # No parameter uncertainty on population growth
+      s <- survivals
+      f <- fecundities
+      lam_it[sim] <- lambda(build_Leslie(s,f))
 
+    } # End if/else
 
 
+    lam0 > 1+rMAX_species
+    pop_growth_mean
 
 
+    model_demo = NULL
 
+    # Choose the model demographique to use (if choice was not forced)
+    if(is.null(model_demo)){
 
+      ## Define the complete model by default
+      model_demo <- M4_WithDD_WithDemoStoch
 
+      # DECLINING (or stable), but initially LARGE population
+      if(lam_it[sim] <= 1 & sum(N0) > 3000) model_demo <- M1_noDD_noDemoStoch
 
+      # DECLINING  (or stable), and initially SMALL population
+      if(lam_it[sim] <= 1 & sum(N0) <= 3000) model_demo <- M2_noDD_WithDemoStoch
 
 
-out <- get_metrics(N)$scenario$impact_sc
-TH <- dim(N)[2]
-nsc <- dim(N)[3]
+      # GROWING population...
+      if(lam_it[sim] > 1){
 
-library(dplyr)
+        # Extract rMAX
+        DD_params$rMAX <- infer_rMAX(K = K, theta = theta,
+                                     pop_size_current = sum(N0), pop_growth_current = lam_it[sim],
+                                     rMAX_theoretical = rMAX_species)
 
-df <- as.data.frame(cbind(year = 1:nrow(out), out[,,1], scenario = 1))
-for(j in 2:dim(out)[3]) df <- rbind(df, cbind(year = 1:nrow(out), out[,,j], scenario = j))
-dim(df)
-class(df)
-names(df)
+        # ... and initially LARGE population
+        if(sum(N0) > 500) model_demo <- M3_WithDD_noDemoStoch
 
-filter(df, scenario == 1)
-filter(df, scenario > 1)
 
+        # ... but initially SMALL population
+        if(sum(N0) <= 500) model_demo <- M4_WithDD_WithDemoStoch
 
-library(ggplot2)
-library(plotly)
+      } # if lam > 1
 
-## pars
-size = 1.5
 
-# Plot lines
-p <-
-  ggplot(data = df, aes(x = year, y = avg)) +
-  geom_line(data = filter(df, scenario > 1), size = size, aes(colour = factor(scenario))) +
-  geom_line(data = filter(df, scenario == 1), size = size, colour = "black")
+    } # end if "is.null"
 
-# Plot CIs
-p <- p + geom_ribbon(data = filter(df, scenario > 1),
-                     aes(ymin = uci, ymax = lci, fill = factor(scenario)), linetype = 0, alpha = 0.100)
 
-# Add legend
-Legend <- "parc"
-nsc <- max(df$scenario)-1
-p <- p + labs(x = "Annee", y = "Impact relatif",
-              col = "Scenario", fill = "Scenario") +
-  scale_color_hue(labels = paste(Legend, 1:nsc), aesthetics = c("colour", "fill"))
-p
 
-# Plot lines
-p <-
-  ggplot(data = df, aes(x = year, y = avg)) +
-  geom_line(data = filter(df, scenario > 1), size = size, aes(colour = factor(scenario))) +
-  geom_line(data = filter(df, scenario == 1), size = size, colour = "black")
 
-# Plot CIs
-p <- p + geom_ribbon(data = filter(df, scenario > 1),
-                     aes(ymin = uci, ymax = lci, fill = factor(scenario)), linetype = 0, alpha = 0.100)
+    fatalities = M
+#######################################################################
+    ##--------------------------------------------
+    ##  Pop project cumulated                               --
+    ##--------------------------------------------
 
-# Add legend
-Legend <- "parc"
-nsc <- max(df$scenario)-1
-p <- p + labs(x = "Annee", y = "Impact relatif",
-              col = "Scenario", fill = "Scenario") +
-  scale_color_hue(labels = paste(Legend, 1:nsc)) +
-  scale_fill_hue(labels = paste(Legend, 1:nsc))
 
-p
-ggplotly(p)
 
-rm(fig)
-fig <- ggplotly(p)
-fig
+    ## Fatalities taking onset time into account
+    # Initiate matrix
+    Mi <- matrix(fatalities, nrow = length(fatalities), ncol = nyr)
 
+    # Fatalities from each wind farm
+    for(j in 2:nrow(Mi)){
+      if(onset_time[j] > 1) Mi[j,1:(onset_time[j]-1)] <- 0
+    } # j
 
+    # Cumulated Fatalities
+    Mc <- Mi
+    for(j in 2:nrow(Mc)) Mc[j,] <- apply(Mc[(j-1):j,], 2, sum)
 
+    # Initiate Pop Size (output) Array
+    N <- array(0, dim = c(nac, nyr, nsc), dimnames = list(paste0("age", 1:nac),
+                                                           paste0("year", 1:nyr),
+                                                           paste0("scenario", 1:nsc)
+    ))
+    N[,1,] <-  N0
 
 
+    ## Loops over time (years)
+    for(t in 2:nyr){
 
-#### Plotly
-p <-
-  ggplot(data = filter(df, scenario > 1), aes(x = year, y = avg)) +
-  geom_line(size = size, aes(colour = factor(scenario)))
+      # Environmental Stochasticity
+      ss = 1 - rnorm(nac, mean = qlogis(1-s), sd = cv_env/(s)) %>% plogis        ## sample thru the mortality rate : 1 - s
+      ff = rnorm(nac, mean = log(f), sd = cv_env) %>% exp
 
+      # Fatalities : constant number (M) or constant rate (h)
+      if(fatal_constant == "M"){
+        h <- Mc[,t-1]/apply(N[,t-1,], 2, sum)
+      } else {
+        h <- Mc[,t-1]/apply(N[,1,], 2, sum)
+      }
 
-# Plot CIs
-p <- p + geom_ribbon(aes(ymin = uci, ymax = lci, fill = factor(scenario)),
-                     linetype = 0, alpha = 0.100)
+      # Sample a seed for RNG
+      seed <- runif(1, 0, 1e6)
 
-# Add legend
-Legend <- "parc"
-nsc <- max(df$scenario)-1
-p <- p + labs(x = "Annee", y = "Impact relatif",
-              col = "Scenario", fill = "Scenario") +
-  scale_color_hue(labels = paste(Legend, 1:nsc), aesthetics = c("colour", "fill"))
-p
+      ## Projection : apply the LESLIE matrix calculation forward
+      # Scenario 0
+      for(j in 1:nsc){
+        set.seed(seed)
+        N[,t,j] <- model_demo(N1 = N[,t-1,j], s = ss, f = ff, h = h[j], DD_params = DD_params)
+      } # j
 
-rm(fig)
-fig <- ggplotly(p)
-fig
+    } # t
 
 
-fig <- ggplotly(p)
-fig
 
 
 
 
+    t = 31
+    j = 4
 
-#### Plotly
 
-p <-
-  ggplot(data = df, aes(x = year, y = avg)) +
-  geom_line(size = size, aes(colour = factor(scenario)))
-p
 
-fig <- ggplotly(p, layerData = 1)
-fig
+    ## M2_noDD_WithDemoStoch
 
+    t = t+1
 
-##################
+    h = sapply(Mc[,t-1]/apply(N[,t-1,], 2, sum), min, 1)
+    h
 
+    #N1 = N[,t,j]
+    #N1 = N2
+    N1 = c(2,1,1,3)
 
-p <- ggplot(data = df, aes(x = year, y = avg, colour = factor(scenario)))
-p <- p + geom_line(size = 2)
-p
+    # Number of age classes
+    nac = length(s)
 
-p <- p + geom_ribbon(aes(ymin = uci, ymax = lci, fill = factor(scenario)), linetype = 0, alpha = 0.075)
-p
+    # Survivors (to "natural mortality" (s) and Wind Turbine Fatalities (1-h))
+    S <- rbinom(nac, N1, (1-h)*s)
+    N2 <- c(rep(0, nac-1), tail(S,1)) + c(0, head(S,-1))
 
+    # Births
+    N2[1] <- sum(rpois(nac, f*N2))
 
-#########
-p <- ggplot(data = df, aes(x = year, y = avg, colour = factor(scenario)))
-p <- p + geom_line(size = 2)
-p
+    N2
+    N[,t+1,j] = N2
diff --git a/man/infer_rMAX.Rd b/man/infer_rMAX.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2fef41001c2777e0fa344ebed79a0e91e837579d
--- /dev/null
+++ b/man/infer_rMAX.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/infer_rMAX.R
+\name{infer_rMAX}
+\alias{infer_rMAX}
+\title{Infer the rMAX}
+\usage{
+infer_rMAX(
+  K,
+  theta = 1,
+  pop_size_current,
+  pop_growth_current,
+  rMAX_theoretical = Inf
+)
+}
+\arguments{
+\item{K}{a strictly positive number. Carrying capacity (= maximum size that the population can ever reach)}
+
+\item{theta}{a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+Note lambda = r + 1}
+
+\item{pop_size_current}{a strictly positive number. Current population size.}
+
+\item{pop_growth_current}{a strictly positive number. Current population growth rate.}
+
+\item{rMAX_theoretical}{the maximum value for rMAX, usually calculated from the function rMAX_spp. See ?rMAX_spp for details.}
+}
+\value{
+the r_MAX value.
+}
+\description{
+Calculates the rMAX of a density-dependence relationship, based on K and current population status (size and trend)
+}
+\examples{
+## Infer rMAX
+infer_rMAX(rMAX = NULL, K = 2000, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+
+
+}
diff --git a/man/rMAX_spp.Rd b/man/rMAX_spp.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c0ee1fec7fc101ee2556d766f670e9191b7a9a72
--- /dev/null
+++ b/man/rMAX_spp.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rMAX_spp.R
+\name{rMAX_spp}
+\alias{rMAX_spp}
+\title{Rmax of a species}
+\usage{
+rMAX_spp(surv, afr)
+}
+\arguments{
+\item{surv}{Adult survival. Single value between 0 and 1.}
+
+\item{afr}{Age at first reproduction. Single integer value.}
+}
+\value{
+the theoretical Rmax of the species
+}
+\description{
+Calculate the theoretical Rmax of a species using the equation of Niel and Lebreton 2005.
+
+References.
+
+Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+populations from incomplete data. Conservation Biology 19:826–835.
+}
+\examples{
+rMAX_spp(surv = 0.6, afr = 2)
+
+}
diff --git a/man/run_simul.Rd b/man/run_simul.Rd
index 94c03ee947ab85864b14d7326a777ae765aad09b..e320bff35f059ea974730313f5835b2c38f918ee 100644
--- a/man/run_simul.Rd
+++ b/man/run_simul.Rd
@@ -17,7 +17,9 @@ run_simul(
   pop_growth_se,
   survivals,
   fecundities,
-  DD_params,
+  carrying_capacity,
+  theta = 1,
+  rMAX_species,
   model_demo = NULL,
   time_horzion,
   coeff_var_environ,
@@ -50,7 +52,21 @@ or the Number of Pairs ("Npair"). A stable age distribution is used to infer the
 
 \item{fecundities}{a vector of fecundity values for each age class.}
 
-\item{DD_params}{NULL or a list. Density-dependence parameters (rMAX, K, theta). Only used in DD models M3 and M4.}
+\item{carrying_capacity}{a strictly positive number.
+Carrying capacity (= maximum size that the population can reach). Here, the unit is the same as pop_size_type.
+It can thus be expressed as the total population or the number of pair.}
+
+\item{theta}{a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+Note lambda = r + 1}
+
+\item{rMAX_species}{the maximum value of rMAX for the species under consideration,
+usually calculated using the Niel & Lebreton (2005) equation.
+It can be calculated using the function rMAX_spp. See ?rMAX_spp for details.
+
+References :
+Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+populations from incomplete data. Conservation Biology 19:826–835.}
 
 \item{model_demo}{is NULL, by default, because the model choice will be made inside each iteration (simulation),
 base on the values of N0 and lam0 that are drawn.
@@ -91,13 +107,17 @@ time_horzion = 30
 coeff_var_environ = 0.10
 fatal_constant = "h"
 
-DD_params <- list(rMAX = NULL, K = 1200, theta = 1)
+carrying_capacity = 1200
+theta = 1
+rMAX_species <- 0.15
 
 run_simul(nsim = 10, cumuated_impacts = FALSE,
            fatalities_mean, fatalities_se, onset_time = NULL,
            pop_size_mean, pop_size_se, pop_size_type,
            pop_growth_mean, pop_growth_se,
-           survivals, fecundities, DD_params = DD_params,
+           survivals, fecundities,
+           carrying_capacity, theta,
+           rMAX_species,
            model_demo = NULL, time_horzion, coeff_var_environ, fatal_constant)
 
 
diff --git a/man/run_simul_shiny.Rd b/man/run_simul_shiny.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1e8bc3385736e6988fba29bdee16c19a6b9bf4af
--- /dev/null
+++ b/man/run_simul_shiny.Rd
@@ -0,0 +1,93 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/run_simul_shiny.R
+\name{run_simul_shiny}
+\alias{run_simul_shiny}
+\title{Run multiple population projections (simulations).
+Exactly the same function as run_simul, except that it includes a line of code to display
+the simulation progress bar in Shiny (incProgress function)}
+\usage{
+run_simul_shiny(
+  nsim,
+  cumuated_impacts,
+  fatalities_mean,
+  fatalities_se,
+  onset_time,
+  pop_size_mean,
+  pop_size_se,
+  pop_size_type,
+  pop_growth_mean,
+  pop_growth_se,
+  survivals,
+  fecundities,
+  carrying_capacity,
+  theta = 1,
+  rMAX_species,
+  model_demo = NULL,
+  time_horzion,
+  coeff_var_environ,
+  fatal_constant
+)
+}
+\arguments{
+\item{nsim}{number of simulation}
+
+\item{cumuated_impacts}{Logical. If TRUE, we used the projection model for cumulated impacts.}
+
+\item{fatalities_mean}{a vector (numeric). Average number of fatalities, for each scenario.}
+
+\item{fatalities_se}{a vector (numeric). Standard Error for the number of fatalities, for each scenario (= uncertainties around the values provided).}
+
+\item{onset_time}{a vector (numeric). The times at which each wind farm fatality starts applying.}
+
+\item{pop_size_mean}{a single number. Average population size (either total population or number of pairs - see Ntype below).}
+
+\item{pop_size_se}{Standard Error for population size (= uncertainty around the value provided).}
+
+\item{pop_size_type}{character value indicating if the provided value pop_size correpsonds to Total Population Size ("Ntotal")
+or the Number of Pairs ("Npair"). A stable age distribution is used to infer the size of each age class.}
+
+\item{pop_growth_mean}{a number. Average population growth rate (lambda).}
+
+\item{pop_growth_se}{Standard Error for population growth rate (= uncertainty around the value provided).}
+
+\item{survivals}{a vector. Average survival probabilities for each age class.}
+
+\item{fecundities}{a vector of fecundity values for each age class.}
+
+\item{carrying_capacity}{a strictly positive number.
+Carrying capacity (= maximum size that the population can reach). Here, the unit is the same as pop_size_type.
+It can thus be expressed as the total population or the number of pair.}
+
+\item{theta}{a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+The relationship is defined as : r <- rMAX*(1-(N/K)^theta)
+Note lambda = r + 1}
+
+\item{rMAX_species}{the maximum value of rMAX for the species under consideration,
+usually calculated using the Niel & Lebreton (2005) equation.
+It can be calculated using the function rMAX_spp. See ?rMAX_spp for details.
+
+References :
+Niel, C., and J. Lebreton. 2005. Using demographic invariants to detect overharvested bird
+populations from incomplete data. Conservation Biology 19:826–835.}
+
+\item{model_demo}{is NULL, by default, because the model choice will be made inside each iteration (simulation),
+base on the values of N0 and lam0 that are drawn.
+But it can be forced by setting the value, which must then be an R object corresponding to the demographic model to be used.
+The 4 possible models currently are: M1_noDD_noDemoStoch, M2_noDD_WithDemoStoch, M3_WithDD_noDemoStoch, M4_WithDD_WithDemoStoch,}
+
+\item{time_horzion}{a number. The number of years (time horizon) over which to project the population dynamics.}
+
+\item{coeff_var_environ}{a number. The coefficient of variation to model environment stochasticity.}
+
+\item{fatal_constant}{text (character). Either "h" or "M". Using "h" sets the fatality RATE as the constant value across years.
+Using "M" sets the NUMBER of fatalities as the constant value across years.}
+}
+\value{
+a 4D array containing the size of each age class (dim 1), for each year (dim 2), each scenario (dim 3), and
+each simulation iteration (dim 4)
+}
+\description{
+Run multiple population projections (simulations).
+Exactly the same function as run_simul, except that it includes a line of code to display
+the simulation progress bar in Shiny (incProgress function)
+}
diff --git a/run_analysis.R b/run_analysis.R
index a615eae305740416f58941e2a659517529e8c1eb..ee7136dc5352daba1254f90b18396a0d400658c8 100644
--- a/run_analysis.R
+++ b/run_analysis.R
@@ -6,35 +6,53 @@ graphics.off()
 library(eolpop)
 
 ## Inputs
-nsim = 100
+nsim = 10
 
-fatalities_mean = c(0, 2, 3, 5, 2)
-fatalities_se = fatalities_mean*0.05
+fatalities_mean = c(0, 10, 5, 8)
+fatalities_se = c(0, 0.05, 0.05, 0.05)
 
 pop_size_mean = 200
-pop_size_se = 30
+pop_size_se = 25
 
-pop_growth_mean = 1.1
+pop_growth_mean = 1
 pop_growth_se = 0
 
 survivals <- c(0.5, 0.7, 0.8, 0.95)
 fecundities <- c(0, 0, 0.05, 0.55)
 
 model_demo = NULL # M2_noDD_WithDemoStoch #M1_noDD_noDemoStoch #M4_WithDD_WithDemoStoch #M3_WithDD_noDemoStoch #
-time_horzion = 30
+time_horzion = 50
 coeff_var_environ = 0.10
-fatal_constant = "h"
-pop_size_type = "Ntotal"
+fatal_constant = "M"
+pop_size_type = "Npair"
 
-cumuated_impacts = TRUE
+cumulated_impacts = TRUE
 
-onset_year = 2000 + c(1, 3, 7, 15, 20)
+onset_year = c(2010, 2013, 2016)
 onset_time = onset_year - min(onset_year) + 1
+onset_time = c(min(onset_time), onset_time)
+
+# Pop size total
+sum(pop_vector(pop_size = pop_size_mean, pop_size_type = pop_size_type, s = survivals, f = fecundities))
+
+
+# Define K
+carrying_capacity = 500
+theta = 1
+K = pop_vector(pop_size = carrying_capacity, pop_size_type = pop_size_type, s = survivals, f = fecundities) %>% sum
+K
+
+# Define theoretical rMAX for the species
+rMAX_species <- rMAX_spp(surv = tail(survivals,1), afr = min(which(fecundities != 0)))
+rMAX_species
+
+##  Avoid unrealistic scenarios
+pop_growth_mean <- min(1 + rMAX_species, pop_growth_mean)
+pop_growth_mean
 
-DD_params <- list(rMAX = NULL, K = 1200, theta = 1)
 
 ##--------------------------------------------
-##  Calibration : FYI, for table dsiply     --
+##  Calibration                             --
 ##--------------------------------------------
 # Calibrate vital rates to match the the desired lambda
 inits <- init_calib(s = survivals, f = fecundities, lam0 = pop_growth_mean)
@@ -42,41 +60,40 @@ vr_calibrated <- calibrate_params(inits = inits, f = fecundities, s = survivals,
 s_calibrated <- head(vr_calibrated, length(survivals))
 f_calibrated <- tail(vr_calibrated, length(fecundities))
 
+build_Leslie(s = s_calibrated, f = f_calibrated) %>% lambda
+
+
+
+
 ##==============================================================================
 ##                         Analyses (simulations)                             ==
 ##==============================================================================
-run0 <- run_simul(nsim, cumuated_impacts,
+run0 <- run_simul(nsim, cumulated_impacts,
                   fatalities_mean, fatalities_se, onset_time,
                   pop_size_mean, pop_size_se, pop_size_type,
                   pop_growth_mean, pop_growth_se,
                   survivals = s_calibrated, fecundities = f_calibrated,
-                  DD_params = DD_params,
+                  carrying_capacity = carrying_capacity, theta = theta,
+                  rMAX_species = rMAX_species,
                   model_demo, time_horzion, coeff_var_environ, fatal_constant)
 
 
-# save(run0, file = "./data/run0.rda")
-names(run0)
-run0$time_run
-# saved time (ratio): 493/12
 
-N <- run0$N ; dim(N)
-out <- get_metrics(N, cumuated_impacts = cumuated_impacts)
-names(out)
 
-out$warning
-names(out$scenario)
+names(run0)
+N <- run0$N ; dim(N)
+plot_traj(N, xlab = "Annee", ylab = "Taille de population (totale)")
+abline(h = K)
 
-fatalities_mean
-out$indiv_farm$impact[time_horzion,,]
+colSums(N[,,,]) %>% max
 
-out$scenario$impact[time_horzion,,]
-out$scenario$Pext
-out$scenario$DR_Pext
+plot_impact(N, onset_year = onset_year , xlab = "Annee", ylab = "Impact relatif")
 
-plot_traj(N, xlab = "Annee", ylab = "Taille de population (totale)")
 
-p <- plot_impact(N, onset_year = onset_year , xlab = "Annee", ylab = "Impact relatif")
-p
+N <- run0$N
+output <- get_metrics(N, cumuated_impacts = cumulated_impacts)
+output$scenario$Pext
 
+#plot_impact(N = N, xlab = "year", ylab = "pop size")
 #source("draws_histog.R")
 #draws_histog(draws = run0$lambdas, mu = pop_growth_mean, se = pop_growth_se)