diff --git a/R/models.R b/R/models.R
index b55fb9ef9c189116f6a3921486f947e3757068d3..f732a188a417e52351693df706a2e1d859155e4e 100644
--- a/R/models.R
+++ b/R/models.R
@@ -189,23 +189,28 @@ M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params,
# Calibrate vital rates to match lam_Nt
A <- build_Leslie(s = s, f = f)
diff_rel_lam <- (lam_Nt - lambda(A))/lambda(A)
- d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
+ #d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
el <- elements_Leslie(s=s, f=f)
vr0 = el$vital_rates
- vr1 = vr0*(1+d)
+ #vr1 = vr0*(1+d)
+ vr1 = vr0*(1+diff_rel_lam)
nac = length(s)
+ s00 <- head(vr1, nac)
+ f00 <- tail(vr1, nac)
+
# Calibrate survivals
- s_Nt <- ((head(vr1, nac)) %>% sapply(., max, 0.05)) %>% sapply(., min, 0.97)
+ s_Nt <- (s00 %>% sapply(., max, 0.05)) %>% sapply(., min, 0.99)
# Calibrate fecundities
- f_Nt <- (tail(vr1, nac)) %>% sapply(., max, 0.001)
+ f00 <- f00 * (1 + (s00 - s_Nt))
+ f_Nt <- f00 %>% sapply(., max, 0.001)
f_Nt[f == 0] <- 0
## 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.05 ){
+ if( abs((lambda(build_Leslie(s = s_Nt, f = f_Nt)) - lam_Nt) / lam_Nt) > 0.05 ){
#If difference is too large : Use optimisation function for better calibration
inits <- c(f_Nt, s_Nt)
@@ -280,19 +285,24 @@ M4_WithDD_WithDemoStoch <- function(N1, s, f, h, DD_params,
# Calibrate vital rates to match lam_Nt
A <- build_Leslie(s = s, f = f)
diff_rel_lam <- (lam_Nt - lambda(A))/lambda(A)
- d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
+ #d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
el <- elements_Leslie(s=s, f=f)
vr0 = el$vital_rates
- vr1 = vr0*(1+d)
+ #vr1 = vr0*(1+d)
+ vr1 = vr0*(1+diff_rel_lam)
nac = length(s)
+ s00 <- head(vr1, nac)
+ f00 <- tail(vr1, nac)
+
# Calibrate survivals
- s_Nt <- ((head(vr1, nac)) %>% sapply(., max, 0.05)) %>% sapply(., min, 0.97)
+ s_Nt <- (s00 %>% sapply(., max, 0.05)) %>% sapply(., min, 0.99)
# Calibrate fecundities
- f_Nt <- (tail(vr1, nac)) %>% sapply(., max, 0.001)
+ f00 <- f00 * (1 + (s00 - s_Nt))
+ f_Nt <- f00 %>% sapply(., max, 0.001)
f_Nt[f == 0] <- 0
## Check if approximation is close enough to desired lambda
diff --git a/inst/ShinyApp/server.R b/inst/ShinyApp/server.R
index 4a717223a37c8c43e14ee3f9dba3ed8c08fc6424..87befd39e4fb8655154445fa2be7fe7686c3feed 100644
--- a/inst/ShinyApp/server.R
+++ b/inst/ShinyApp/server.R
@@ -1235,7 +1235,7 @@ server <- function(input, output, session){
# Case 1 : Values from expert elicitation
if(input$pop_growth_input_type == "eli_exp"){
if(!(is.null(param$pop_growth_eli_result))){
- param$pop_growth_mean <- round(min(1 + param$rMAX_species, round(param$pop_growth_eli_result$mean, 2)), 4)
+ param$pop_growth_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 {
@@ -1276,16 +1276,14 @@ server <- function(input, output, session){
if(input$pop_growth_input_type == "val"){
# Case 2 : Values directly provided as mean & SE
ready$pop_growth <- TRUE
- param$pop_growth_mean <- round(min(1 + param$rMAX_species, make_lambda(input$pop_growth_mean)), 4)
+ param$pop_growth_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(min(1 + param$rMAX_species,
- get_mu(lower = make_lambda(input$pop_growth_lower),
- upper = make_lambda(input$pop_growth_upper)))
- , 4)
+ param$pop_growth_mean <- round(get_mu(lower = make_lambda(input$pop_growth_lower),
+ upper = make_lambda(input$pop_growth_upper)), 4)
param$pop_growth_se <- round(get_sd(lower = make_lambda(input$pop_growth_lower),
upper = make_lambda(input$pop_growth_upper), coverage = CP), 5)
} # end (if3)
@@ -1355,6 +1353,8 @@ server <- function(input, output, session){
input$button_calibrate_vr
},{
+ print(param$pop_growth_mean)
+
vr_calib0 <- calibrate_params(
inits = init_calib(s = param$survivals, f = param$fecundities, lam0 = param$pop_growth_mean),
f = param$fecundities, s = param$survivals, lam0 = param$pop_growth_mean
@@ -1372,12 +1372,14 @@ server <- function(input, output, session){
rMAX_species <- rMAX_spp(surv = tail(param$survivals,1), afr = min(which(param$fecundities != 0)))
param$rMAX_species <- rMAX_species
+ param$pop_growth_mean_use <- round(min(1 + rMAX_species, param$pop_growth_mean), 4)
+
param$theta <- fixed_theta
#param$theta <- theta_spp(rMAX_species)
param$vr_calibrated <- calibrate_params(
- inits = init_calib(s = param$survivals, f = param$fecundities, lam0 = param$pop_growth_mean),
- f = param$fecundities, s = param$survivals, lam0 = param$pop_growth_mean
+ 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))
@@ -1428,6 +1430,8 @@ server <- function(input, output, session){
input$run
}, {
+ print(param$pop_growth_mean_use)
+
if(ready$fatalities & ready$pop_size & ready$pop_growth & ready$carrying_capacity){
out$analysis_choice <- input$analysis_choice
@@ -1446,7 +1450,7 @@ server <- function(input, output, session){
pop_size_se = param$pop_size_se,
pop_size_type = param$pop_size_unit,
- pop_growth_mean = param$pop_growth_mean,
+ pop_growth_mean = param$pop_growth_mean_use,
pop_growth_se = param$pop_growth_se,
survivals = param$s_calibrated,
diff --git a/inst/ShinyApp/ui.R b/inst/ShinyApp/ui.R
index 02f17cbc899a6680cd104a5e06a01aa7bfc79067..b0bd9f1e1286fdfa3a24ef5b60290c057d68bf2f 100644
--- a/inst/ShinyApp/ui.R
+++ b/inst/ShinyApp/ui.R
@@ -512,18 +512,18 @@ rm(list = ls(all.names = TRUE))
numericInput(inputId = "pop_growth_lower",
label = HTML("Borne inf�rieure<br>(taux d'accroissement en %)"),
value = -6,
- min = -100, max = 100, step = 1),
+ min = -100, max = Inf, step = 1),
numericInput(inputId = "pop_growth_upper",
label = HTML("Borne sup�rieure<br>(taux d'accroissement en %)"),
value = -6,
- min = -100, max = 100, step = 1),
+ min = -100, max = Inf, step = 1),
## Input values: mean and se
numericInput(inputId = "pop_growth_mean",
label = "Moyenne (taux d'accroissement en %)",
value = -1.5,
- min = -100, max = 100, step = 1),
+ min = -100, max = Inf, step = 1),
numericInput(inputId = "pop_growth_se",
label = "Erreur-type (aussi en %)",
diff --git a/run_analysis.R b/run_analysis.R
index b90cdc4d295dd5537973716f581a170208d33c15..4a83a3fbb1770ea245ae56c82cad05e091e88b20 100644
--- a/run_analysis.R
+++ b/run_analysis.R
@@ -23,15 +23,17 @@ fatalities_mean = c(0, 3) #c(0, 5, 3, 4, 2, 1, 4, 2, 2, 3)
fatalities_se = c(0, 0.582) # c(0, rep(0.5,9))
length(fatalities_mean)
-survivals <- c(0.65, 0.75, 0.85, 0.94)
-fecundities <- c(0, 0, 0.05, 0.40)
-
+#survivals <- c(0.65, 0.75, 0.85, 0.94)
+#fecundities <- c(0, 0, 0.05, 0.40)
#survivals <- c(0.47, 0.67, 0.67)
#fecundities <- c(0, 0.30, 1.16)
#survivals <- c(0.25, 0.30)
#fecundities <- c(0, 19.8)
-pop_growth_mean = 0.94
+survivals <- c(0.3, 0.65)
+fecundities <- c(0, 4.5)
+
+pop_growth_mean = 1.98
# lambda( build_Leslie(s = survivals, f = fecundities) )
pop_growth_se = 0
@@ -70,8 +72,8 @@ theta = 1
##
rMAX_use <- infer_rMAX(K = K, theta = theta,
- pop_size_current = sum(N000), pop_growth_current = pop_growth_mean,
- rMAX_theoretical = rMAX_species)
+ pop_size_current = sum(N000), pop_growth_current = pop_growth_mean,
+ rMAX_theoretical = rMAX_species)
rMAX_use
rMAX_species
@@ -93,43 +95,41 @@ vr_calibrated <- calibrate_params(inits = inits, f = fecundities, s = survivals,
s_calibrated <- head(vr_calibrated, length(survivals))
f_calibrated <- tail(vr_calibrated, length(fecundities))
-lambda( build_Leslie(s = s_calibrated, f = f_calibrated) )
s_calibrated
f_calibrated
+lambda( build_Leslie(s = s_calibrated, f = f_calibrated) )
-
-length(survivals)
##==============================================================================
## Analyses (simulations) ==
##==============================================================================
time <- system.time(
-run0 <- run_simul(nsim = nsim,
- cumulated_impacts = cumulated_impacts,
+ run0 <- run_simul(nsim = nsim,
+ cumulated_impacts = cumulated_impacts,
- fatalities_mean = fatalities_mean,
- fatalities_se = fatalities_se,
- onset_time = onset_time,
+ fatalities_mean = fatalities_mean,
+ fatalities_se = fatalities_se,
+ onset_time = onset_time,
- pop_size_mean = pop_size_mean,
- pop_size_se = pop_size_se,
- pop_size_type = pop_size_type,
+ pop_size_mean = pop_size_mean,
+ pop_size_se = pop_size_se,
+ pop_size_type = pop_size_type,
- pop_growth_mean = pop_growth_mean,
- pop_growth_se = pop_growth_se,
+ pop_growth_mean = pop_growth_mean,
+ pop_growth_se = pop_growth_se,
- survivals = s_calibrated,
- fecundities = f_calibrated,
+ survivals = s_calibrated,
+ fecundities = f_calibrated,
- carrying_capacity_mean = carrying_capacity_mean,
- carrying_capacity_se = carrying_capacity_se,
+ carrying_capacity_mean = carrying_capacity_mean,
+ carrying_capacity_se = carrying_capacity_se,
- theta = theta,
- rMAX_species = rMAX_species,
+ theta = theta,
+ rMAX_species = rMAX_species,
- model_demo = NULL,
- time_horizon = time_horizon,
- coeff_var_environ = coeff_var_environ,
- fatal_constant = fatal_constant)
+ model_demo = NULL,
+ time_horizon = time_horizon,
+ coeff_var_environ = coeff_var_environ,
+ fatal_constant = fatal_constant)
)
#####################################################
@@ -181,7 +181,7 @@ sum(NN[,1,1])
x11()
plot_traj(N, age_class_use = "pairs", fecundities = fecundities,
- Legend = paste("sc", 1:length(fatalities_mean)), ylim = c(0, NA))
+ Legend = paste("sc", 1:length(fatalities_mean)), ylim = c(0, NA))
plot_traj(N, age_class_use = "NotJuv0", fecundities = fecundities,
Legend = paste("sc", 1:length(fatalities_mean)), ylim = c(0, NA))