From 1576c4295989918e27ca36434524b5784d869968 Mon Sep 17 00:00:00 2001
From: Maxime Jaunatre <maxime.jaunatre@yahoo.fr>
Date: Thu, 18 Aug 2022 17:59:48 +0200
Subject: [PATCH 1/3] Expl work

I test if I can retrieve binded env values, and reconstruct function
with it in plain text inside.
And this is working just nicely now. I'll build this tomorrow !
---
 dev/dev_way_for_recru_mod.R | 108 ++++++++++++++++++++++++++++++++++++
 1 file changed, 108 insertions(+)
 create mode 100644 dev/dev_way_for_recru_mod.R

diff --git a/dev/dev_way_for_recru_mod.R b/dev/dev_way_for_recru_mod.R
new file mode 100644
index 0000000..2b64df6
--- /dev/null
+++ b/dev/dev_way_for_recru_mod.R
@@ -0,0 +1,108 @@
+# Issue 1 : recrutment function comes with env ####
+
+#' In recruitmenent function, the function is written and get values from
+#' a binded environment coming from function definition.
+
+species <- "Yggdrasil"
+climatic <- 1
+path <- here()
+replicat <- 42
+
+fIPM <- here(path, "output", species, paste0("IPM_Clim_", climatic, ".Rds"))
+raw_IPM <- readRDS(assertFileExists(fIPM)) # NOTE 10" to load...
+assertNumber(replicat, lower = 1, upper = length(raw_IPM))
+raw_IPM <- raw_IPM[[replicat]]
+
+foo <- function(x) {
+    browser()
+    raw_IPM$RecFun(x)
+}
+
+foo(1)
+
+environment(raw_IPM$RecFun)$K_i
+
+# function (BATOTSP)
+# {
+#     as.numeric(K_i + K_BA * BATOTSP + K_logBA * log(BATOTSP))
+# }
+# <environment: 0x5616bc88bc18> # WHY IS THERE A BINDED env HERE ???
+
+raw_IPM$rec$params_m
+raw_IPM$list_m
+
+# Browse[3]> K_i
+# intercept
+# -0.7360761
+# Browse[3]> K_BA
+# BATOTSP
+# -0.0179456
+# Browse[3]> K_logBA
+# logBATOTSP
+# -0.2630412
+
+# Browse[3]> ls(envir = where("K_i"))
+# [1] "IsFunc"               "K_BA"
+# [3] "K_i"                  "K_i_inter"
+# [5] "K_logBA"              "L"
+# [7] "L1"                   "L2"
+# [9] "list_covs"            "mu_rec"
+# [11] "params_i"             "params_i_inter"
+# [13] "params_i_no"          "params_logsize_inter"
+# [15] "params_rec"           "params_size_inter"
+# [17] "r_rec"
+
+# Finding a way to build function from reg ####
+
+#' Is it possible to write numeric values as is in a function ??
+#' maybe with expr
+
+x <- 2
+library(rlang)
+foo <- expr(x + 2)
+foo
+eval_bare(foo, x = 3)
+
+x <- 1
+y <- 2
+z <- call2("<-", expr(x), y + 1)
+z
+eval(z)
+x
+y
+
+call2("{",
+      call2(
+          "+", expr(x), {1}
+      )
+)
+
+z <- 3
+tot <- function(x, y){x + z * y}
+tot
+tot(1, y = 0.5) # 2.5
+formals(tot)
+names(formals(tot))[2] <- "Y" # Change argument name !
+formals(tot)
+body(tot)
+body(tot) <- call2("{",
+                   call2( "+", expr(x), call2("*", z, expr(Y))
+                   ) )
+tot
+tot(1, Y = 0.5) # 2.5
+# tot(1, y = 0.5)
+
+test <- call2(
+    "<-", expr(res), expr(function(x){1})
+)
+eval(test)
+res
+res()
+
+function(x) {1}
+
+
+lobstr::ast(foo <- function(x){x +1})
+
+# it works !
+
-- 
GitLab


From f00c55a5f262f28262874ee52bc39b797ffa4e64 Mon Sep 17 00:00:00 2001
From: Maxime Jaunatre <maxime.jaunatre@yahoo.fr>
Date: Fri, 19 Aug 2022 18:25:40 +0200
Subject: [PATCH 2/3] Almost done, missing full usage

I can recreate the function way more easily but the new function has a
binded environment and I can't manage to get rid of it.

But it's way better in the end so I'll use this solution !
---
 .Rbuildignore                    |   3 +-
 R/class_species.R                |   2 +-
 R/load_oldIPM_rec.R              | 130 +++++++++++++++++++++++++++++++
 dev/dev_way_for_recru_mod.R      |  37 ++++++++-
 dev/tryhard.R                    |   2 +
 tests/testthat/test-oldIPM_rec.R |  67 ++++++++++++++++
 6 files changed, 238 insertions(+), 3 deletions(-)
 create mode 100644 R/load_oldIPM_rec.R
 create mode 100644 tests/testthat/test-oldIPM_rec.R

diff --git a/.Rbuildignore b/.Rbuildignore
index 0614dd7..4c87069 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -2,10 +2,11 @@
 ^\.Rproj\.user$
 ^LICENSE\.md$
 
-dev/
+^dev$
 ^data-raw$
 
 R/R_0_sp.R
 R/Sim_DemStock.R
 R/old_Sim_NonDem.R
 ^output$
+^outputSim$
diff --git a/R/class_species.R b/R/class_species.R
index aff11a4..dd211ea 100644
--- a/R/class_species.R
+++ b/R/class_species.R
@@ -117,7 +117,7 @@ species <- function(IPM, init_pop, harvest_fun, recruit_fun){
 #' @export
 old_ipm2species <- function(species, climatic = 1, path = here(), replicat = 42,
                             harvest = def_harv, init_pop = def_init){
-
+    species <- "Yggdrasil"
     assertCharacter(species, len = 1)
     assertCharacter(path, len = 1)
     assertCount(climatic)
diff --git a/R/load_oldIPM_rec.R b/R/load_oldIPM_rec.R
new file mode 100644
index 0000000..77766f9
--- /dev/null
+++ b/R/load_oldIPM_rec.R
@@ -0,0 +1,130 @@
+#' Create call from old name
+#'
+#' Get variable name and value in fitted derivated object and return a call
+#'
+#' This function is used to load old IPM and modify the recruitment function
+#' to use more arguments (BAtot and BAsp) in case of multiple species
+#'
+#' @param x Name of the argument. Can be preceded with log
+#' (log variable is computed) or end with 2 (square)chr.
+#' @param df2 Dataframe with x in var1 column and a column named K with numeric
+#' values. Data.frame
+#'
+#' @examples
+#' df <- data.frame(var1 = "BATOTSP", var2 = NA_character_, params = -0.0179,
+#' value.x = 1, value.y = 1, K = -0.0179)
+#' multi(x = "BATOTSP", df)
+#' #> BATOTSP_in <- -0.0179 * BATOTSP
+#'
+#' # Square variable
+#' df <- data.frame(var1 = "BATOTSP2", var2 = NA_character_, params = -0.0179,
+#' value.x = 1, value.y = 1, K = -0.0179)
+#' multi(x = "BATOTSP2", df)
+#' #> BATOTSP2_in <- -0.0179 * BATOTSP ^ 2
+#'
+#' @importFrom rlang ensym call2
+#' @import checkmate
+#' @noRd
+multi <- function(x, df){
+
+    assertCharacter(x, any.missing = FALSE, len = 1)
+
+    if(grepl("log", x)){
+        evar <- sub("^log", "", x)
+        evar <- call2("log", ensym(evar))
+    } else if(grepl("2$", x)){
+        evar <- sub("2$", "", x)
+        evar <- call2("^", ensym(evar), 2)
+    } else {
+        evar <- ensym(x)
+    }
+
+    tmp <- call2("*", df[df$var1 == x, "K"], evar)
+    new_name <- paste0(x, "_in")
+    res <- call2("<-", ensym(new_name), tmp)
+
+    return(res)
+}
+
+
+#' Format fitted derivated object to table
+#'
+#' Format values to be used in other functions
+#'
+#' @param params Estimated parameters for the fit of the model.
+#' @param list_covs Climatic covariates values.
+#'
+#' @importFrom dplyr left_join mutate
+#' @importFrom tidyr pivot_longer everything separate
+#' @importFrom tibble rownames_to_column
+#'
+#' @noRd
+format_fit <- function(params, list_covs){
+
+    invar <- names(params)[!names(params) %in% names(list_covs)]
+
+    lc <- pivot_longer(list_covs, cols = everything())
+    lc <- rbind(lc, data.frame(name= c(invar, NA),  value=1))
+    p <- as.data.frame(params) %>% rownames_to_column(var = "var") %>%
+        separate(var, c("var1", "var2"), sep = "\\:", fill = "right")
+    res <- left_join(p, lc, by=c('var1'='name')) %>%
+        left_join(lc, by=c('var2'='name')) %>%
+        mutate(K = params * value.x * value.y)
+
+    return(res)
+}
+
+
+
+#' Export recruitment function from estimated parameters.
+#'
+#' @param params Estimated parameters for the fit of the model.
+#' @param list_covs Climatic covariates values.
+#'
+#' @importFrom purrr map
+#' @importFrom dplyr filter
+#' @importFrom rlang expr call2
+#'
+#' @return
+#' Function with 4 parameters : BATOTSP, BATOTNonSP, mesh and SurfEch
+#'
+#' @noRd
+exp_recFun <- function(params, list_covs){
+
+    df2 <- format_fit(params, list_covs)
+
+    invar <- names(params)[!names(params) %in% names(list_covs)]
+    invar <- invar[! grepl("ntercept", invar)]
+    inter <- sum(filter(df2, ! var1 %in% invar | var2 %in% invar)$K)
+
+
+    exp_invar <- map(invar, multi, df2)
+    add_invar <- map(c(list(expr(intercept <- 1)),exp_invar),
+                     ~ call2("<-", expr(res), call2("+", expr(res), .x[[2]] )))
+
+    final_res <- list(
+        expr(mesh <- length(mesh)),
+        expr(distrib <- c(rep(1/2, 2), numeric(mesh - 2)) ),
+        expr(final <- exp(res) * SurfEch / 0.03 * distrib ),
+        expr(return(final))
+    )
+    calls <- c(exp_invar, add_invar, final_res)
+
+    empty <- function(BATOTSP, BATOTNonSP, mesh, SurfEch = 0.003){}
+
+    body(empty)[[2]] <- call2("<-", expr(intercept), inter)
+    body(empty)[[3]] <- expr(res <- 0)
+    for(i in seq_along(calls)){
+        body(empty)[[i + 3]] <- calls[[i]]
+    }
+
+    # this is messy and is to remove binded env.
+    env_unbind(env = environment(empty), c("i", "calls", "final_res",
+                                           "add_invar", "exp_invar",
+                                           "inter", "invar", "df2"),
+               inherit = FALSE)
+    # empty
+    return(empty)
+}
+
+
diff --git a/dev/dev_way_for_recru_mod.R b/dev/dev_way_for_recru_mod.R
index 2b64df6..b0a5e64 100644
--- a/dev/dev_way_for_recru_mod.R
+++ b/dev/dev_way_for_recru_mod.R
@@ -2,6 +2,8 @@
 
 #' In recruitmenent function, the function is written and get values from
 #' a binded environment coming from function definition.
+load_all()
+library(rlang)
 
 species <- "Yggdrasil"
 climatic <- 1
@@ -13,6 +15,24 @@ raw_IPM <- readRDS(assertFileExists(fIPM)) # NOTE 10" to load...
 assertNumber(replicat, lower = 1, upper = length(raw_IPM))
 raw_IPM <- raw_IPM[[replicat]]
 
+
+func <- raw_IPM$RecFun
+raw_IPM$RecFun
+raw_IPM$rec$formula
+raw_IPM$rec$params_m
+raw_IPM$list_m
+
+str(lobstr::ast(!!eval(parse(text = raw_IPM$rec$formula))))
+
+e <- body(func)
+ts <- e[[2]]
+lobstr::ast(!!ts)
+lobstr::sxp(ts)
+rlang::parse_expr(!!ts)
+parse(text  = parse(text = ts)[[2]])[[1]]
+
+accumulate(ts, ~ parse(text = .x), .dir = "backward")
+
 foo <- function(x) {
     browser()
     raw_IPM$RecFun(x)
@@ -78,12 +98,18 @@ call2("{",
 )
 
 z <- 3
-tot <- function(x, y){x + z * y}
+tot <- function(x, y){
+    tmp <- z * y
+    tmp <- x +  tmp
+    tmp
+    }
 tot
 tot(1, y = 0.5) # 2.5
 formals(tot)
 names(formals(tot))[2] <- "Y" # Change argument name !
 formals(tot)
+body(tot)[[4]] <- call2("<-", expr(tmp), call2("+", 1, expr(tmp)))
+body(tot)[[5]] <- call2("return", expr(tmp))
 body(tot)
 body(tot) <- call2("{",
                    call2( "+", expr(x), call2("*", z, expr(Y))
@@ -106,3 +132,12 @@ lobstr::ast(foo <- function(x){x +1})
 
 # it works !
 
+# Final test ####
+
+list_covs <- raw_IPM$list_m
+params <- raw_IPM$rec$params_m
+
+foo <- exp_recFun(params = raw_IPM$rec$params_m, list_covs = raw_IPM$list_m)
+foo
+foo(1, 2, 1:5, 0.03)
+
diff --git a/dev/tryhard.R b/dev/tryhard.R
index a593e8f..4aef777 100644
--- a/dev/tryhard.R
+++ b/dev/tryhard.R
@@ -8,6 +8,8 @@ o_Yggdrasil$info["species"] <- "Pinea"
 Forest <- forest(list(Yggdrasil, o_Yggdrasil))
 Forest <- forest(list(Yggdrasil))
 
+Forest$species$Yggdrasil$recruit_fun
+
 set.seed(42)
 res <- sim_deter_forest(Forest, tlim = 10, equil_time = 1e3,
                  correction = "cut",
diff --git a/tests/testthat/test-oldIPM_rec.R b/tests/testthat/test-oldIPM_rec.R
new file mode 100644
index 0000000..a20ca82
--- /dev/null
+++ b/tests/testthat/test-oldIPM_rec.R
@@ -0,0 +1,67 @@
+test_that("multi works", {
+    df <- data.frame(var1 = "BATOTSP", var2 = NA_character_, params = -0.0179,
+                     value.x = 1, value.y = 1, K = -0.0179)
+    expect_identical(
+        multi(x = "BATOTSP", df = df),
+        call2("<-", expr(BATOTSP_in), call2("*", -0.0179, expr(BATOTSP)))
+    )
+
+    df <- data.frame(var1 = "BATOTSP2", var2 = NA_character_, params = -0.0179,
+                     value.x = 1, value.y = 1, K = -0.0179)
+    expect_identical(
+        multi(x = "BATOTSP2", df = df),
+        call2("<-", expr(BATOTSP2_in), call2("*", -0.0179, expr(BATOTSP ^ 2)))
+    )
+
+    df <- data.frame(var1 = "logBATOTSP", var2 = NA_character_, params = -0.0179,
+                     value.x = 1, value.y = 1, K = -0.0179)
+    expect_identical(
+        multi(x = "logBATOTSP", df = df),
+        call2("<-", expr(logBATOTSP_in), call2("*", -0.0179, expr(log(BATOTSP))))
+    )
+})
+
+
+test_that("format_fit works", {
+    params <- c(intercept = -0.864, BATOTSP = -0.018, sgddb = 286.813,
+                wai = -0.057, wai2 = 0.288
+    )
+
+    list_covs <- data.frame(wai = -0.187, sgddb = 0, waib = 1.23, wai2 = 0.34)
+
+    exp <- data.frame(
+        var1 = c("intercept", "BATOTSP", "sgddb", "wai","wai2"),
+        var2 = NA_character_,
+        params = c(-0.864, -0.018, 286.813,             -0.057, 0.288),
+        value.x = c(1, 1, 0, -0.187, 0.34),
+        value.y = c(1, 1, 1, 1, 1),
+        K = c(-0.864, -0.018, 0, 0.010659, 0.097920))
+
+    expect_identical( format_fit(params, list_covs), exp)
+})
+
+load_all()
+test_that("format_fit works", {
+    params <- c(intercept = -0.864, BATOTSP = -0.018, sgddb = 286.813,
+                wai = -0.057, wai2 = 0.288
+    )
+
+    list_covs <- data.frame(wai = -0.187, sgddb = 0, waib = 1.23, wai2 = 0.34)
+
+    empty <- function (BATOTSP, BATOTNonSP, mesh, SurfEch = 0.003) {
+        intercept <- -0.755421
+        res <- 0
+        BATOTSP_in <- -0.018 * BATOTSP
+        res <- res + intercept
+        res <- res + BATOTSP_in
+        mesh <- length(mesh)
+        distrib <- c(rep(1/2, 2), numeric(mesh - 2))
+        final <- exp(res) * SurfEch/0.03 * distrib
+        return(final)
+    }
+    body(empty)[[2]][[3]] <- -0.755421
+    body(empty)[[4]][[3]][[2]] <- -0.018
+
+
+    expect_identical(exp_recFun(params, list_covs),  empty)
+})
-- 
GitLab


From 49bb79bc89b65a1fb40c0f825fdc2ffdc7be250d Mon Sep 17 00:00:00 2001
From: Maxime Jaunatre <maxime.jaunatre@yahoo.fr>
Date: Mon, 22 Aug 2022 13:39:41 +0200
Subject: [PATCH 3/3] Almost working but not yet

Integrated modification in Sim_Deter, tested all recrutment function
eport.

I tested in dev/tryhard.R and it's almost here but this is still
different.

After discussion, this difference is expected since now some intra sp
competition goes to extra sp competition during recruitment.
As the fit give more impact to intra sp competition, higher BA is now
expected. We can also observe that this competition impact the
equilibrium state with more low mesh individuals.
---
 NAMESPACE                           |   9 ++
 R/Sim_Deter.R                       |  25 +++--
 R/class_species.R                   |   8 +-
 R/load_oldIPM_rec.R                 |  17 +++-
 dev/dev_issue1.R                    |  58 +++++++++++
 dev/dev_way_for_recru_mod.R         | 143 ----------------------------
 dev/run_all_test.R                  |   4 +-
 dev/tryhard.R                       |  43 ++++++---
 tests/testthat/test-class_species.R |   7 +-
 tests/testthat/test-gen_meth.R      |  25 ++---
 tests/testthat/test-oldIPM_rec.R    |   6 +-
 vignettes/multisp_deter_sim.Rmd     |  33 +++++--
 12 files changed, 185 insertions(+), 193 deletions(-)
 create mode 100644 dev/dev_issue1.R
 delete mode 100644 dev/dev_way_for_recru_mod.R

diff --git a/NAMESPACE b/NAMESPACE
index cbf9ff6..b332a5e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -41,12 +41,21 @@ import(checkmate)
 import(here)
 import(purrr)
 importFrom(dplyr,between)
+importFrom(dplyr,filter)
+importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,relocate)
 importFrom(purrr,"%>%")
+importFrom(purrr,map)
 importFrom(purrr,map_chr)
 importFrom(rlang,.data)
+importFrom(rlang,call2)
+importFrom(rlang,ensym)
+importFrom(rlang,env_unbind)
+importFrom(rlang,expr)
 importFrom(stats,rbinom)
 importFrom(stats,runif)
 importFrom(tibble,rownames_to_column)
+importFrom(tidyr,everything)
 importFrom(tidyr,pivot_longer)
+importFrom(tidyr,separate)
diff --git a/R/Sim_Deter.R b/R/Sim_Deter.R
index 6b48380..7feef57 100644
--- a/R/Sim_Deter.R
+++ b/R/Sim_Deter.R
@@ -155,6 +155,7 @@ sim_deter_forest <- function(Forest,
         ncol = nsp, nrow = tlim + 2, dimnames = list(NULL, names(Forest$species))
     ))
     sim_BA <- rep(NA_real_, equil_time)
+    sim_BAnonSp <- rep(NA_real_, nsp)
 
     ## Initiate pop ####
     X <- map2(map(Forest$species, `[[`, "init_pop"),
@@ -167,6 +168,8 @@ sim_deter_forest <- function(Forest,
     # save first pop
     sim_BAsp[1, ] <- map2_dbl(X, ct, ~ .x %*% .y )
     sim_BA[1] <- sum(sim_BAsp[1,])
+    sim_BAnonSp <- map2_dbl( - sim_BAsp[1, ,drop = FALSE], sim_BA[1],  `+`)
+
     # tmp <- map2(X, sim_BAsp[1, ,drop = TRUE], ~ c(.x, .y, sum(.x)))
     tmp <- imap(X, function(x, .y, ba, harv){
         c(x, ba[[.y]], sum(x), harv[[.y]], sum(harv[[.y]]) )
@@ -202,11 +205,6 @@ sim_deter_forest <- function(Forest,
             x} )
     }
 
-    # QUESTION change this ?
-    b  <- map(map(Forest$species, ~ length(get_mesh(.x))),
-               ~ c(rep(1 / 2, 2), numeric(.x - 2))
-    )
-
     if (verbose) {
         message("Starting while loop. Maximum t = ", equil_time)
     }
@@ -223,14 +221,21 @@ sim_deter_forest <- function(Forest,
         X <- map2(sim_ipm, X, ~ drop( .x %*% .y ) ) # Growth
         Harv <- map2(map(Forest$species, `[[`, "harvest_fun"), X, exec) # Harvest
         X <- map2(X, Harv, `-`)
-        recrues <- map2(map(Forest$species, `[[`, "recruit_fun"),
-                        # QUESTION is this BA below ??
-                        map2(ct, X, `%*%`), exec)
-        recrues <- map2(recrues, b, ~ exp(.x) * SurfEch / 0.03 * .y)
+
+        recrues <- imap(
+            map(Forest$species, `[[`, "recruit_fun"),
+            function(x, .y, basp, banonsp, mesh, SurfEch){
+                exec(x, basp[[.y]], banonsp[.y], mesh[[.y]], SurfEch)
+            }, basp = sim_BAsp[t-1,,drop = FALSE], banonsp = sim_BAnonSp,
+            mesh = map(Forest$species, get_mesh), SurfEch = SurfEch )
+
         X <- map2(X, recrues, `+`) # Recruitment
         # compute new BA for selecting the right IPM and save values
-        sim_BAsp[t, ] <- map2_dbl(X, ct, ~ .x %*% .y )
+        sim_BAsp[t, ] <- map2_dbl(X, ct, `%*%`)
         sim_BA[t] <- sum(sim_BAsp[t,])
+        sim_BAnonSp <- map2_dbl( - sim_BAsp[t, ,drop = FALSE], sim_BA[t],  `+`)
+
+        # Update X
         if (t <= tlim) {
             # tmp <- map2(X, sim_BAsp[t, ,drop = TRUE], ~ c(.x, .y, sum(.x)))
             tmp <- imap(X, function(x, .y, ba, harv){
diff --git a/R/class_species.R b/R/class_species.R
index dd211ea..8a593ca 100644
--- a/R/class_species.R
+++ b/R/class_species.R
@@ -55,7 +55,8 @@ validate_species <- function(x){
     assertFunction(values$init_pop, args = c("mesh", "SurfEch"))
     assertFunction(values$harvest_fun, args = c("x"))
     # TODO : check that X return >= 0 values of same length
-    assertFunction(values$recruit_fun, args = c("BATOTSP")) # FIXME args must be changed.
+    assertFunction(values$recruit_fun, args = c("BATOTSP", "BATOTNonSP",
+                                                "mesh", "SurfEch"))
     # check infos ####
     assertCharacter(values$info, any.missing = FALSE)
     if(any(names(values$info) != c("species", "climatic"))){
@@ -117,7 +118,7 @@ species <- function(IPM, init_pop, harvest_fun, recruit_fun){
 #' @export
 old_ipm2species <- function(species, climatic = 1, path = here(), replicat = 42,
                             harvest = def_harv, init_pop = def_init){
-    species <- "Yggdrasil"
+
     assertCharacter(species, len = 1)
     assertCharacter(path, len = 1)
     assertCount(climatic)
@@ -135,7 +136,8 @@ old_ipm2species <- function(species, climatic = 1, path = here(), replicat = 42,
 
     res <- species(
         IPM = res_ipm, init_pop = init_pop, harvest_fun = harvest,
-        recruit_fun = raw_IPM$RecFun
+        recruit_fun = exp_recFun(params = raw_IPM$rec$params_m,
+                                 list_covs = raw_IPM$list_m)
     )
 
     return(res)
diff --git a/R/load_oldIPM_rec.R b/R/load_oldIPM_rec.R
index 77766f9..38c15f2 100644
--- a/R/load_oldIPM_rec.R
+++ b/R/load_oldIPM_rec.R
@@ -78,16 +78,31 @@ format_fit <- function(params, list_covs){
 
 #' Export recruitment function from estimated parameters.
 #'
+#' Rebuild the function to use BASp and BAnonSp for a species.
+#'
 #' @param params Estimated parameters for the fit of the model.
 #' @param list_covs Climatic covariates values.
 #'
 #' @importFrom purrr map
 #' @importFrom dplyr filter
-#' @importFrom rlang expr call2
+#' @importFrom rlang expr call2 env_unbind
+#'
+#' @details
+#' Each function has an environment binded with params and list_covs.
+#' I can't remove it and it may be usefull later after all.
 #'
 #' @return
 #' Function with 4 parameters : BATOTSP, BATOTNonSP, mesh and SurfEch
 #'
+#' @examples
+#' params <- c(intercept = -0.864, BATOTSP = -0.018, sgddb = 286.813,
+#' wai = -0.057, wai2 = 0.288 )
+#' list_covs <- data.frame(wai = -0.187, sgddb = 0, waib = 1.23, wai2 = 0.34)
+#'
+#' foo <- exp_recFun(params, list_covs)
+#' foo
+#' foo(1, 2, 1:5, 0.03)
+#'
 #' @noRd
 exp_recFun <- function(params, list_covs){
 
diff --git a/dev/dev_issue1.R b/dev/dev_issue1.R
new file mode 100644
index 0000000..e8a46cc
--- /dev/null
+++ b/dev/dev_issue1.R
@@ -0,0 +1,58 @@
+# Working on multiple species systems ####
+
+document()
+load_all()
+spe <- "Yggdrasil"
+# Yggdrasil <- old_ipm2species(spe, path = "tests/testthat/testdata/", replicat = 1)
+Yggdrasil <- old_ipm2species(spe, path = here(), replicat = 1)
+Ents <- Yggdrasil
+Ents$info["species"] <- "Ents"
+
+Forest <- forest(list(Yggdrasil))
+
+# Forest$species$Yggdrasil$recruit_fun
+
+set.seed(42)
+res <- sim_deter_forest(Forest, tlim = 600, equil_time = 1e3,
+                        correction = "cut",
+                        verbose = TRUE)
+
+times <- as.numeric(sub("t", "", colnames(res)))
+plot(times, res[grepl("BA",rownames(res)),], ylab = "Total BA", xlab = "time",
+     cex = 0.1, type = "b", ylim = c(0, 100))
+text(700, res[grepl("BA",rownames(res)),ncol(res)],
+     round(res[grepl("BA",rownames(res)),ncol(res)], 2), cex = .7, pos = 3)
+
+set.seed(42)
+Forest2 <- forest(list(Yggdrasil, Ents))
+res2 <- sim_deter_forest(Forest2, tlim = 600, equil_time = 1e3,
+                         correction = "cut",
+                         verbose = TRUE)
+
+
+value <- colSums(res2[grepl("BA",rownames(res2)),])
+points(times, value, type = "b", cex = 0.1, col = "darkred")
+text(700, value[ncol(res2)], round(value[ncol(res2)], 2), cex = .7, pos = 3)
+points(times, res2[grepl("Ygg.*BA",rownames(res2)),],
+       type = "b", cex = 0.1, col = "darkblue")
+points(times, res2[grepl("Ent.*BA",rownames(res2)),],
+       type = "b", cex = 0.1, col = "darkgreen")
+legend("top", c("Yddgrasil solo", "Yggdrasil + Ents", "Yggdrasil", "Ents"),
+       fill = c("black", "darkred", "darkblue", "darkgreen"))
+text(700, res2[grepl("BA",rownames(res2)),ncol(res2)],
+     round(res2[grepl("BA",rownames(res2)),ncol(res2)], 2), cex = .7, pos = 3)
+
+#' La difference de BA s'expliquerais par une plus faible competition sur
+#' les cellules de mesh petits car la competition extrasp s'applique
+#' lors du recrutement.
+
+eq_1 <- res[grepl("m",rownames(res)),ncol(res)]
+eq_2_Ygg <- res2[grepl("Ygg.*m",rownames(res2)),ncol(res2)]
+eq_2_Ent <- res2[grepl("Ent.*m",rownames(res2)),ncol(res2)]
+eq_2 <- eq_2_Ent + eq_2_Ygg
+
+plot(eq_2- - eq_1, type = "b", cex = .2,
+     xlab = "mesh cell", ylab = "Net difference",
+     main = paste("Difference in equilibrium distribution between",
+                  "\nintra-specific and hetero-specific competition"))
+
diff --git a/dev/dev_way_for_recru_mod.R b/dev/dev_way_for_recru_mod.R
deleted file mode 100644
index b0a5e64..0000000
--- a/dev/dev_way_for_recru_mod.R
+++ /dev/null
@@ -1,143 +0,0 @@
-# Issue 1 : recrutment function comes with env ####
-
-#' In recruitmenent function, the function is written and get values from
-#' a binded environment coming from function definition.
-load_all()
-library(rlang)
-
-species <- "Yggdrasil"
-climatic <- 1
-path <- here()
-replicat <- 42
-
-fIPM <- here(path, "output", species, paste0("IPM_Clim_", climatic, ".Rds"))
-raw_IPM <- readRDS(assertFileExists(fIPM)) # NOTE 10" to load...
-assertNumber(replicat, lower = 1, upper = length(raw_IPM))
-raw_IPM <- raw_IPM[[replicat]]
-
-
-func <- raw_IPM$RecFun
-raw_IPM$RecFun
-raw_IPM$rec$formula
-raw_IPM$rec$params_m
-raw_IPM$list_m
-
-str(lobstr::ast(!!eval(parse(text = raw_IPM$rec$formula))))
-
-e <- body(func)
-ts <- e[[2]]
-lobstr::ast(!!ts)
-lobstr::sxp(ts)
-rlang::parse_expr(!!ts)
-parse(text  = parse(text = ts)[[2]])[[1]]
-
-accumulate(ts, ~ parse(text = .x), .dir = "backward")
-
-foo <- function(x) {
-    browser()
-    raw_IPM$RecFun(x)
-}
-
-foo(1)
-
-environment(raw_IPM$RecFun)$K_i
-
-# function (BATOTSP)
-# {
-#     as.numeric(K_i + K_BA * BATOTSP + K_logBA * log(BATOTSP))
-# }
-# <environment: 0x5616bc88bc18> # WHY IS THERE A BINDED env HERE ???
-
-raw_IPM$rec$params_m
-raw_IPM$list_m
-
-# Browse[3]> K_i
-# intercept
-# -0.7360761
-# Browse[3]> K_BA
-# BATOTSP
-# -0.0179456
-# Browse[3]> K_logBA
-# logBATOTSP
-# -0.2630412
-
-# Browse[3]> ls(envir = where("K_i"))
-# [1] "IsFunc"               "K_BA"
-# [3] "K_i"                  "K_i_inter"
-# [5] "K_logBA"              "L"
-# [7] "L1"                   "L2"
-# [9] "list_covs"            "mu_rec"
-# [11] "params_i"             "params_i_inter"
-# [13] "params_i_no"          "params_logsize_inter"
-# [15] "params_rec"           "params_size_inter"
-# [17] "r_rec"
-
-# Finding a way to build function from reg ####
-
-#' Is it possible to write numeric values as is in a function ??
-#' maybe with expr
-
-x <- 2
-library(rlang)
-foo <- expr(x + 2)
-foo
-eval_bare(foo, x = 3)
-
-x <- 1
-y <- 2
-z <- call2("<-", expr(x), y + 1)
-z
-eval(z)
-x
-y
-
-call2("{",
-      call2(
-          "+", expr(x), {1}
-      )
-)
-
-z <- 3
-tot <- function(x, y){
-    tmp <- z * y
-    tmp <- x +  tmp
-    tmp
-    }
-tot
-tot(1, y = 0.5) # 2.5
-formals(tot)
-names(formals(tot))[2] <- "Y" # Change argument name !
-formals(tot)
-body(tot)[[4]] <- call2("<-", expr(tmp), call2("+", 1, expr(tmp)))
-body(tot)[[5]] <- call2("return", expr(tmp))
-body(tot)
-body(tot) <- call2("{",
-                   call2( "+", expr(x), call2("*", z, expr(Y))
-                   ) )
-tot
-tot(1, Y = 0.5) # 2.5
-# tot(1, y = 0.5)
-
-test <- call2(
-    "<-", expr(res), expr(function(x){1})
-)
-eval(test)
-res
-res()
-
-function(x) {1}
-
-
-lobstr::ast(foo <- function(x){x +1})
-
-# it works !
-
-# Final test ####
-
-list_covs <- raw_IPM$list_m
-params <- raw_IPM$rec$params_m
-
-foo <- exp_recFun(params = raw_IPM$rec$params_m, list_covs = raw_IPM$list_m)
-foo
-foo(1, 2, 1:5, 0.03)
-
diff --git a/dev/run_all_test.R b/dev/run_all_test.R
index c2ed075..a2162db 100644
--- a/dev/run_all_test.R
+++ b/dev/run_all_test.R
@@ -8,7 +8,8 @@ devtools::check()
 # library(covr, testthat)
 # get a shiny to find which line is not yet tested. Very helpfull
 report(x = package_coverage(line_exclusions = list("R/Sim_NonDem.R")))
-report(x = package_coverage())
+x <- package_coverage()
+report(x)
 # covr::package_coverage()
 
 # rerun all coverage for dataframe analysis.
@@ -24,6 +25,7 @@ x <- file_coverage('R/class_ipm.R', 'tests/testthat/test-class_ipm.R')
 x <- file_coverage('R/class_species.R', 'tests/testthat/test-class_species.R')
 x <- file_coverage('R/class_forest.R', 'tests/testthat/test-class_forest.R')
 x <- file_coverage('R/generic_methods.R', 'tests/testthat/test-gen_meth.R')
+x <- file_coverage('R/load_oldIPM_rec.R', 'tests/testthat/test-oldIPM_rec.R')
 report(x)
 # zero_coverage() shows only uncovered lines.
 # If run within RStudio, `zero_coverage()` will open a marker pane with the
diff --git a/dev/tryhard.R b/dev/tryhard.R
index 4aef777..6f3b63b 100644
--- a/dev/tryhard.R
+++ b/dev/tryhard.R
@@ -1,24 +1,43 @@
 document()
 load_all()
 spe <- "Yggdrasil"
-Yggdrasil <- old_ipm2species(spe, path = "tests/testthat/testdata/", replicat = 1)
-o_Yggdrasil <- Yggdrasil
-o_Yggdrasil$info["species"] <- "Pinea"
+# Yggdrasil <- old_ipm2species(spe, path = "tests/testthat/testdata/", replicat = 1)
+Yggdrasil <- old_ipm2species(spe, path = here(), replicat = 1)
+Ents <- Yggdrasil
+Ents$info["species"] <- "Ents"
 
-Forest <- forest(list(Yggdrasil, o_Yggdrasil))
 Forest <- forest(list(Yggdrasil))
-
-Forest$species$Yggdrasil$recruit_fun
+# Forest$species$Yggdrasil$recruit_fun
 
 set.seed(42)
-res <- sim_deter_forest(Forest, tlim = 10, equil_time = 1e3,
+res <- sim_deter_forest(Forest, tlim = 600, equil_time = 1e3,
                  correction = "cut",
                  verbose = TRUE)
 
-res[c(1:2, 700:702, 703,1402:1403, 1404:1405, 2103:2106, 2805:2806),
-         c(1:3,30:31)]
-res[c(1:2, 700:702, 703,1402:1403), c(1:3,30:31)]
+times <- as.numeric(sub("t", "", colnames(res)))
+plot(times, res[grepl("BA",rownames(res)),], ylab = "Total BA", xlab = "time",
+     cex = 0.1, type = "b", ylim = c(0, 100))
+# plot(res[grepl("N",rownames(res)),])
+text(700, res[grepl("BA",rownames(res)),ncol(res)],
+     round(res[grepl("BA",rownames(res)),ncol(res)], 2), cex = .7, pos = 3)
+
+set.seed(42)
+Forest2 <- forest(list(Yggdrasil, Ents))
+res2 <- sim_deter_forest(Forest2, tlim = 600, equil_time = 1e3,
+                        correction = "cut",
+                        verbose = TRUE)
+
+
+value <- colSums(res2[grepl("BA",rownames(res2)),])
+points(times, value, type = "b", cex = 0.1, col = "darkred")
+text(700, value[ncol(res2)], round(value[ncol(res2)], 2), cex = .7, pos = 3)
+points(times, res2[grepl("Ygg.*BA",rownames(res2)),],
+       type = "b", cex = 0.1, col = "darkblue")
+points(times, res2[grepl("Ent.*BA",rownames(res2)),],
+       type = "b", cex = 0.1, col = "darkgreen")
+legend("top", c("Yddgrasil solo", "Yggdrasil + Ents", "Yggdrasil", "Ents"),
+       fill = c("black", "darkred", "darkblue", "darkgreen"))
+text(700, res2[grepl("BA",rownames(res2)),ncol(res2)],
+     round(res2[grepl("BA",rownames(res2)),ncol(res2)], 2), cex = .7, pos = 3)
 
-plot(res[grepl("BA",rownames(res)),])
-plot(res[grepl("N",rownames(res)),])
 
diff --git a/tests/testthat/test-class_species.R b/tests/testthat/test-class_species.R
index cbbc9fd..33529a6 100644
--- a/tests/testthat/test-class_species.R
+++ b/tests/testthat/test-class_species.R
@@ -27,7 +27,9 @@ test_that("validate_species works", {
     raw_IPM <- raw_IPM[[1]]
     IPM <- old_ipm2ipm("Yggdrasil", climatic = 1, path = path, replicat = 1)
 
-    x <- new_species(IPM, def_init, def_harv, raw_IPM$RecFun)
+    x <- new_species(IPM, def_init, def_harv,
+                     exp_recFun(params = raw_IPM$rec$params_m,
+                                list_covs = raw_IPM$list_m))
 
 
     expect_identical(x, validate_species(x))
@@ -59,7 +61,8 @@ test_that("old_ipm2species works", {
         old_ipm2species("Yggdrasil", climatic = 1, path = path, replicat = 1),
         new_species(
             old_ipm2ipm("Yggdrasil", climatic = 1, path = path, replicat = 1),
-            def_init, def_harv, raw_IPM$RecFun)
+            def_init, def_harv, exp_recFun(params = raw_IPM$rec$params_m,
+                                           list_covs = raw_IPM$list_m))
     )
 
 })
diff --git a/tests/testthat/test-gen_meth.R b/tests/testthat/test-gen_meth.R
index fd789f0..ff71df8 100644
--- a/tests/testthat/test-gen_meth.R
+++ b/tests/testthat/test-gen_meth.R
@@ -110,15 +110,18 @@ test_that("delay forest works", {
     species = "darwin", climatic = 1, compress = FALSE)
     validate_ipm(IPM)
 
-    harv <- function(BATOTSP){NULL}
+    rec <- function(BATOTSP, BATOTNonSP, mesh, SurfEch){NULL}
 
     sp <- new_species(IPM = IPM, init_pop = def_init,
-                     harvest_fun = def_harv, recruit_fun = harv)
+                     harvest_fun = def_harv, recruit_fun = rec)
 
     x <- new_forest(list(darwin = sp))
     exp <- new_forest(list(darwin = delay(sp, 2)))
     # HACK ne pas nommer les elements ici fout la merde.
 
+    # validate_forest(x)
+    # validate_forest(exp)
+
     expect_identical( delay(x, 2), exp )
     expect_identical( delay(x, 0), x )
 
@@ -167,18 +170,18 @@ test_that("delay species works", {
     species = "darwin", climatic = 1, compress = FALSE)
     validate_ipm(IPM)
 
-    harv <- function(BATOTSP){NULL}
+    rec <- function(BATOTSP, BATOTNonSP, mesh, SurfEch){NULL}
 
-    x <- new_species(IPM = IPM, init_pop = def_init,
-                     harvest_fun = def_harv, recruit_fun = harv)
+    sp <- new_species(IPM = IPM, init_pop = def_init,
+                      harvest_fun = def_harv, recruit_fun = rec)
 
     exp <- new_species(correction(IPM, "cut"), init_pop = def_init,
-                       harvest_fun = def_harv, recruit_fun = harv)
-    # validate_species(x)
+                       harvest_fun = def_harv, recruit_fun = rec)
+    # validate_species(sp)
     # validate_species(exp)
 
-    expect_identical( correction(x, "cut"), exp )
-    expect_identical( correction(x), x )
+    expect_identical( correction(sp, "cut"), exp )
+    expect_identical( correction(sp), sp )
 
 })
 
@@ -193,10 +196,10 @@ test_that("delay forest works", {
     species = "darwin", climatic = 1, compress = FALSE)
     validate_ipm(IPM)
 
-    harv <- function(BATOTSP){NULL}
+    rec <- function(BATOTSP, BATOTNonSP, mesh, SurfEch){NULL}
 
     sp <- new_species(IPM = IPM, init_pop = def_init,
-                      harvest_fun = def_harv, recruit_fun = harv)
+                      harvest_fun = def_harv, recruit_fun = rec)
 
     x <- new_forest(list(darwin = sp))
     exp <- new_forest(list(darwin = correction(sp, "cut")))
diff --git a/tests/testthat/test-oldIPM_rec.R b/tests/testthat/test-oldIPM_rec.R
index a20ca82..07e241e 100644
--- a/tests/testthat/test-oldIPM_rec.R
+++ b/tests/testthat/test-oldIPM_rec.R
@@ -40,7 +40,7 @@ test_that("format_fit works", {
     expect_identical( format_fit(params, list_covs), exp)
 })
 
-load_all()
+
 test_that("format_fit works", {
     params <- c(intercept = -0.864, BATOTSP = -0.018, sgddb = 286.813,
                 wai = -0.057, wai2 = 0.288
@@ -62,6 +62,6 @@ test_that("format_fit works", {
     body(empty)[[2]][[3]] <- -0.755421
     body(empty)[[4]][[3]][[2]] <- -0.018
 
-
-    expect_identical(exp_recFun(params, list_covs),  empty)
+    expect_identical(formals(exp_recFun(params, list_covs)),  formals(empty))
+    expect_identical(body(exp_recFun(params, list_covs)),  body(empty))
 })
diff --git a/vignettes/multisp_deter_sim.Rmd b/vignettes/multisp_deter_sim.Rmd
index da11ca2..ada2211 100644
--- a/vignettes/multisp_deter_sim.Rmd
+++ b/vignettes/multisp_deter_sim.Rmd
@@ -15,7 +15,8 @@ knitr::opts_chunk$set(
 ```
 
 ```{r setup}
-library(treeforce)
+# library(treeforce)
+devtools::load_all()
 ```
 
 # Introduction
@@ -35,7 +36,9 @@ Functions are defined to load the previous IPM saved as *.Rds*. The files must l
 
 Once the files are in place, the code below allow to create the species of interest. One last argument is *replicat* which select one of the 100 simulations of previous work. **A mean IPM should replace this later**.
 
-Default functions of initialization and harvest are used. You can write your own functions as long as they use the same arguments. Try to write function that don't initialize negative populations
+Default functions of initialization and harvest are used. You can write your own functions as long as they use the same arguments. Try to write function that don't initialize negative populations.
+
+Example data set is a cropped IPM with a mesh of length 30 with BA from 1 to 10.
 
 ```{r, loading species Yggdrasil}
 spe <- "Yggdrasil" # Example dataset given with the package
@@ -55,10 +58,13 @@ A forest is just a list of the species. **Later it will also require harvest rul
 Forest <- forest(list(Yggdrasil))
 
 set.seed(42)
-res <- sim_deter_forest(Forest, tlim = 100, equil_time = 1e3,
+res <- sim_deter_forest(Forest, tlim = 30, equil_time = 1e3,
                  correction = "cut", equil_dist = 50,
                  verbose = TRUE)
+```
 
+
+```{r, single_print}
 m <- length(Forest$species$Yggdrasil$IPM$mesh)
 res[c(1:2, m:(m+3), (2 *m + 2):(2 *m + 3)), c(1:3,30:31)]
 times <- as.numeric(sub("t", "", colnames(res)))
@@ -74,21 +80,34 @@ We expect the same equilibrium result (minus the initial population random effec
 ```{r, two sp sim}
 Ents <- Yggdrasil
 Ents$info["species"] <- "Ents"
-Forest2 <- forest(list(Yggdrasil, Ents))
+Roger <- Yggdrasil
 
 set.seed(42)
-res <- sim_deter_forest(Forest2, tlim = 100, equil_time = 1e3,
+res <- sim_deter_forest(Forest2, tlim = 30, equil_time = 1e3,
                  correction = "cut", equil_dist = 50,
                  verbose = TRUE)
+```
 
+
+```{r, two_print, echo = FALSE}
 m2 <- m *2 + 3
 res[c(1:2, m:(m+3), (2 *m + 2):m2, 
       (m2 +1):(m2 +2), (m2 + m):(m2 + m +3), (m2 + 2 *m + 2):(2*m2)),
          c(1:2,30:31)]
 
 times <- as.numeric(sub("t", "", colnames(res)))
+value <- colSums(res[grepl("BA",rownames(res)),])
+plot(times, value, type = "b", ylab = "Total BA", xlab = "time", 
+     cex = 0.1, ylim = c(1, max(value) * 1.05))
+points(times, res[grepl("Ygg.*BA",rownames(res)),], 
+     type = "b", ylab = "Total BA", xlab = "time", cex = 0.1, col = "darkblue")
+points(times, res[grepl("Ent.*BA",rownames(res)),], 
+     type = "b", ylab = "Total BA", xlab = "time", cex = 0.1, col = "darkgreen")
+```
+
+Run all previous code
+
+```{r}
 
-plot(times, colSums(res[grepl("BA",rownames(res)),]), 
-     type = "b", ylab = "Total BA", xlab = "time", cex = 0.1)
 ```
 
-- 
GitLab