---
title: "Basic functions and examples"
package: matreex
output: 
    github_document:
    rmarkdown::html_vignette:
vignette: >
  %\VignetteIndexEntry{Basic functions and examples}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
bibliography: references.bib 
link-citations: true
---

```{r, include = FALSE}
# load_all()
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(cli.progress_show_after = 600) 
# 10 minutes before showing a cli progress bar
```

This vignette illustrates the basic functions used to run simulations projecting the population of size-structured trees with `{matreex}` package. The package can model the dynamic of monospecific or plurispecific tree communities. Different modules allow also to simulate harvesting and disturbances. To project forward the dynamic of size-structured tree populations, this package relies on **integrated projection models** (later named IPM).


# Rapid description of the Integral Projection Model

Details on the fitting and the integration of IPM model can be found in @ellner2016.
Briefly, an IPM predicts the size distribution, $n(z', t+1)$, of a population at time t+1 from its size distribution at time t, $n(z, t)$, with $z$ the size at t and $z'$ the size at $t + 1$,
based on the following equation (@ellner2016): 

$$
n(z', t+1) = \int_{L}^{U} K(z', z) n(z, t) dz 
$$

with $L$ and $U$ being, respectively, the lower and upper observed sizes for integration of the kernel $K$.

The kernel $K(z' ,z)$ can be split into the survival and growth kernel $P(z' ,z)$ and the fecundity kernel $F(z', z)$, as follows : $K(z' ,z) = P(z' ,z) + F(z' ,z)$ .

The fecundity kernel $F(z', z)$ gives the size distribution of newly recruited trees at time $t+1$ as a function of the size distribution at time t. The survival and growth kernel $P(z', z)$ is defined as $P(z', z) = s(z) \times G(z', z)$, $s$ being the survival function and $G$ the growth function. 
The kernel $K(z' ,z)$, thus, integrate the three key vital rates functions: growth, survival, and recruitment. The kernel $P$ is numerically approximated with a big iteration matrix and the continuous size distribution $n$ is approximated by a big state vector. The dimension and the width of the size class are selected to ensure a good numerical integration of the kernel $P$. Details on the numerical integration are given below and in @kunstler2021 and @guyennon2023 . 
<!--The figure 1 presents the structure and approach of the IPM. -->
Note, that this package do not cover the statistical fitting of the vital rates functions.


# Simulations input

## Define a species

Before simulating forest, we need first to define tree species in R. Species regroups basic vital rates functions: growth, recruitment and survival. In this package, we provide the fitted vital rates functions used in @kunstler2021 and @guyennon2023 . These functions depend on tree size, local competition based on the sum of basal area of competitors, and two climatic variables. The climatic variables are the sum of the sum of growing degree days *sgdd*, the water aridity index *wai* (see @kunstler2021 for more details). @kunstler2021 and @guyennon2023 used a resampling procedure to estimate the vital rates functions resulting in 100 resampled estimate of the parameters of each functions, here to simplify the simulations we provide only the averaged parameters over those 100 resamples. It is however, in theory, possible to run simulations with other vital rates functions if they are provided as glm objects (please contact the authors to test new variables).

```{r hide library, echo = FALSE}
# I do this before to show data inline
library(matreex)
data("fit_Picea_abies")
```

To build an IPM for a species we start from fitted function @kunstler2021 and @guyennon2023 . To do the numerical integration of $P$, we need to define the mesh dimension, here $700$, the lower $L$ and upper limits $U$, here respectively 90mm and `get_maxdbh(fit_Picea_abies) * 1.1` = `r get_maxdbh(fit_Picea_abies) * 1.1`mm, and range of value of competition index BA (here from 0 to 60 $m^2/ha$). These are the values used in @kunstler2021 and @guyennon2023 that were optimized to provide good numerical integration.

**Please keep in mind this computation is intensive and may take few minutes!**

```{r make_IPM_Picea}
library(matreex)
library(dplyr)
library(ggplot2)

# Load fitted model for a species
# fit_species # list of all species in dataset
data("fit_Picea_abies")

# Load associated climate
data("climate_species")
climate <- subset(climate_species, N == 2 & sp == "Picea_abies", select = -sp)
# N here is a climate defined in Kunstler et al 2021. 
# N == 2 is the optimum climate for the species. 
# see ?climate_species for more info.
climate

Picea_ipm <- make_IPM(
    species = "Picea_abies", 
    climate = climate, 
    fit = fit_Picea_abies,
    clim_lab = "optimum clim", 
    mesh = c(m = 700, L = 90, U = get_maxdbh(fit_Picea_abies) * 1.1),
    BA = 0:60, # Default values are 0:200, smaller values speed up this vignette.
    verbose = TRUE
)
```

```{r hide warn, echo=FALSE, }
options(W_matreex_edist = FALSE)
```

Once the IPM is integrated on a BA range, we can use it to build a species object. In R, a species is a list object that is constructed with the `species()` function. In addition to the IPM kernel $P$, this list require few more functions to work during simulations :

* `init_pop` : Function to draw the initial size distribution. The default one draw distribution for a basal area (later named BA) of 1 $m^2/ha$ with random functions (see the help of the function for more details). The package provide other functions to draw random distribution at a selected BA (`def_initBA()`) or a given distribution (`def_init_k()`).

* `harvest_fun` : Function that cut tree given the size distribution. The default function cut 0.6% per year of the trees regardless of their size. Further functions allow to harvest according to Uneven and Even rules (see Harvesting Vignette).

* `disturb_fun` : Function that return tree mortality after a disturbance. The default is no disturbance.

A species also comes with few parameters, but they are only used for harvest and disturbance modules.

For this example, we will just modify the initial size distribution to start at a basal area of 30 $m^2/ha$.

```{r species}
Picea_sp <- species(IPM = Picea_ipm, init_pop = def_initBA(30))
```

## Define a forest

Once each species objects are built, we can assemble them in a `forest` object. This function also require additional parameters, but they are only used for harvest and disturbance modules.

```{r forest}
Picea_for <- forest(species = list(Picea = Picea_sp))
```

# Running simulations

A simulations project forward a size-structured population from its initial state with the matrix kernel $P$ and the recruitment function. This function requires the length of the simulation and the time limit to search for an equilibrium. A simulation is run until the time limit is reached and can continue further if an equilibrium is not reached. Another parameter is the simulated surface `SurfEch`, it's define the surface of the studied forest. This parameter is mainly here for historical purpose as models were fitted on $300m^2$ plot, and output is scaled to one hectare.

```{r sp1sim}
set.seed(42) # The seed is here for initial population random functions.
Picea_sim <- sim_deter_forest(
    Picea_for, 
    tlim = 200, 
    equil_time = 300, equil_dist = 50, equil_diff = 1,
    SurfEch = 0.03,
    verbose = TRUE
)
```

In the code above, we simulate for 200 years (`tlim`) years and past this time we continue the simulation till it reach an equilibrium up to 300 years (see figure A). After `tlim`, the  simulation will continue until the population reach an equilibrium up to a maximum of 300 year (`equil_time`) . 


The criteria for reaching an equilibrium (in green) is based on computing the range of variation of BA for a moving window of length `equil_dist` since the current step (in grey). The equilibrium is reached if the range of variation within this moving window is less than `equil_diff`. The equilibrium is reached at the first timestep for which the basal area range is lower than `equil_diff`. The steps between `tlim` and `t_equil` are not recorded. The search for the equilibrium start at `tlim` over the last `equil_dist` steps (see figure B). This is why `equil_dist` must not be higher than `tlim`. 
If we want to register the full dynamic, we can set `tlim = equil_time` (see figure C). **The equilibrium is always the last size distribution** (shown in green in figure). Note that in this case, the final distribution will be returned in the result twice. The `equil_dist` and `equil_diff` parameters are not important in this case.

```{r timeline_explain, engine='tikz', fig.ext = 'png', engine.opts = list(template = "tex/tikz.tex"), echo = FALSE, eval=TRUE}
\begin{tikzpicture}

    \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{A}};
    \draw[ultra thick] (0,0) -- (12,0);
    \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
    \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
    \draw[ultra thick] (6,0) node[below=3pt, thick] {$t_{lim}$};
    \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{equil\_time}$};
    \draw[ultra thick] (10,0) node[below=3pt, thick] {$t_{equil}$} {};
    \foreach \x in {0, 2,4,6,8,10, 12}
    \draw (\x cm,3pt) -- (\x cm,-3pt);

    \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(6,0) 
       node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation};
   %  \draw[green, line width=40pt](0,1) -- +(6,0);
    \fill[blue!10] (0,0.5) rectangle ++(6,2.5);
    \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (6,3.2)  -- +(3.9,0) 
         node [black,midway,above=4pt, font=\scriptsize] {Discarded steps};
    \fill[orange!10] (6,0.5) rectangle ++ (3.9,2.5);

      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (7,1.8)  -- +(3,0) 
         node [black,midway,above=4pt, font=\scriptsize] {equil\_dist};
      \draw[lightgray, line width=10pt](7,1.6) -- +(2.9,0);
      \draw [thick ,decorate,decoration={brace,amplitude=5pt, mirror}] (10,1.4)  -- +(0,0.4) 
      node [black,midway,right=4pt, font=\scriptsize] {equil\_diff};
      \draw[green, line width=10pt](9.9,1.6) -- +(.1,0);

    \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
      xtick=\empty,ytick=\empty, width=11.5cm,height=5.5cm,
      xmin=0, xmax=4, ymin=0, ymax=3, ylabel=$BA$, ylabel near ticks
      ]
  \addplot[mark=*,mark size=1pt,color=blue!70!black,samples=200] {1.2+exp(-x)*sin(5*x)};
  \end{axis}
  
  \node[font=\scriptsize, above] at (11.5, 3.2) {Recorded equilibrium};
  \draw[green, ->] plot [smooth, tension=1] coordinates {(10, 1.8) (10.2,2.2) (11.1, 2.7) (11.5, 3.2)};
  

   \begin{scope}[yshift=-5.5cm]

      \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{B}};
      \draw[ultra thick] (0,0) -- (12,0);
      \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
      \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
      \draw[ultra thick] (8,0) node[below=3pt, thick] {$t_{lim}$};
      \draw[ultra thick] (10,0) node[below=3pt, thick] {$t_{equil}$} {};
      \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{equil\_time}$};
      \foreach \x in {0, 2,4,6,8,10, 12}
      \draw (\x cm,3pt) -- (\x cm,-3pt);
  
      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(8,0) 
         node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation};
     %  \draw[green, line width=40pt](0,1) -- +(6,0);
      \fill[blue!10] (0,0.5) rectangle ++(8,2.5);
      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (8,3.2)  -- +(1.9,0) 
         node [black,midway,above=4pt, font=\scriptsize] {Discarded steps};
      \fill[orange!10] (8,0.5) rectangle ++ (1.9,2.5);
  
        \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (7,1.8)  -- +(3,0) 
           node [black,midway,above=4pt, font=\scriptsize] {equil\_dist};
        \draw[lightgray, line width=10pt](7,1.6) -- +(2.9,0);
        \draw [thick ,decorate,decoration={brace,amplitude=5pt, mirror}] (10,1.4)  -- +(0,0.4) 
        node [black,midway,right=4pt, font=\scriptsize] {equil\_diff};
        \draw[green, line width=10pt](9.9,1.6) -- +(.1,0);
  
      \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
        xtick=\empty,ytick=\empty, width=11.5cm,height=5.5cm,
        xmin=0, xmax=4, ymin=0, ymax=3, ylabel=$BA$, ylabel near ticks
        ]
    \addplot[mark=*,mark size=1pt,color=blue!70!black,samples=200] {1.2+exp(-x)*sin(5*x)};
    \end{axis}

    \node[font=\scriptsize, above] at (11.5, 3.2) {Recorded equilibrium};
    \draw[green, ->] plot [smooth, tension=1] coordinates {(10, 1.8) (10.2,2.2) (11.1, 2.7) (11.5, 3.2)};
    

     \end{scope}

     \begin{scope}[yshift=-11cm]

      \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{C}};
      \draw[ultra thick] (0,0) -- (12,0);
      \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
      \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
      \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{lim}$};
      \draw[ultra thick] (12,0) node[below=16pt, thick] {$t_{equil\_time}$};
      \foreach \x in {0, 2,4,6,8,10, 12}
      \draw (\x cm,3pt) -- (\x cm,-3pt);
   
      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(12,0) 
      node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation (include equilibrium)};
   \fill[blue!10] (0,0.5) rectangle ++(12,2.5);
   \draw [gray, thick ,decorate,decoration={brace,amplitude=5pt, mirror}] (12,1.4)  -- +(0,0.4) 
   node [gray,midway,right=4pt, font=\scriptsize] {equil\_diff};


   \draw[green, line width=10pt](11.9,1.6) -- +(.1,0);

   \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
     xtick=\empty,ytick=\empty, width=13.5cm,height=5.5cm,
     xmin=0, xmax=5, ymin=0, ymax=3, ylabel=$BA$, ylabel near ticks
     ]
      \addplot[mark=*,mark size=1pt,color=blue!70!black,samples=200] {1.2+exp(-x)*sin(5*x)};
   \end{axis}

        \begin{scope}
         \clip(9.5, 3) rectangle (12, 0);
         \draw [gray, decorate,decoration={brace,amplitude=12pt}] (9,1.8) -- +(3,0)
            node [gray,midway,above=8pt, font=\scriptsize] {equil\_dist};
      \end{scope}
      \draw[gray, dotted] ([yshift=6pt]9, 1.8) -- ([yshift=6pt]9.5,1.8);

      \node[font=\scriptsize, above] at (13.5, 3.2) {Recorded equilibrium};
      \draw[green, ->] plot [smooth, tension=1] coordinates {(12, 1.1) (12.2,1.8) (13.1, 2.5) (13.5, 3.2)};
      
     \end{scope}

\end{tikzpicture}
```

Keep in mind that despite high `tlim` and `equil_time` values, the equilibrium may not be reached at the end of the simulation. There is currently no way in the algorithm to report this to the user. The best way to detect "false equilibrium" is when the last step takes place at `t == equil_time` and by plotting the basal area along time. These cases are illustrated in figure D and E.

Also, this equilibrium is only computed on total basal area, and the distribution can change. We welcome any suggestions you may have regarding the equilibrium definition.

```{r timeline_explain_no_eq, engine='tikz', fig.ext = 'png', engine.opts = list(template = "tex/tikz.tex"), echo = FALSE, eval=TRUE}
\begin{tikzpicture}

    \begin{scope}

      \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{D}};
      \draw[ultra thick] (0,0) -- (12,0);
      \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
      \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
      \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{lim}$};
      \draw[ultra thick] (12,0) node[below=16pt, thick] {$t_{equil\_time}$};
      \foreach \x in {0, 2,4,6,8,10, 12}
      \draw (\x cm,3pt) -- (\x cm,-3pt);
   
      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(12,0) 
      node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation (include "equilibrium")};
   \fill[blue!10] (0,0.5) rectangle ++(12,2.5);
   \draw[green, line width=10pt](11.9,1.1) -- +(.1,0);

   \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
     xtick=\empty,ytick=\empty, width=13.5cm,height=5.5cm,
     xmin=0, xmax=5, ymin=0, ymax=5, ylabel=$BA$, ylabel near ticks
     ]
      \addplot[mark=*,mark size=1pt,color=blue!70!black,samples=200] {2+sin(2*x)};
   \end{axis}

   \node[font=\scriptsize, above] at (13.5, 3.2) {Recorded "equilibrium"};
   \draw[green, ->] plot [smooth, tension=1] coordinates {(12, 1.1) (12.2,1.8) (13.1, 2.5) (13.5, 3.2)};
      
     \end{scope}

     \begin{scope}[yshift=-5.5cm]

      \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{E}};
      \draw[ultra thick] (0,0) -- (12,0);
      \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
      \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
      \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{lim}$};
      \draw[ultra thick] (12,0) node[below=16pt, thick] {$t_{equil\_time}$};
      \foreach \x in {0, 2,4,6,8,10, 12}
      \draw (\x cm,3pt) -- (\x cm,-3pt);
   
      \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(6,0) 
      node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation };
   \fill[blue!10] (0,0.5) rectangle ++(6,2.5);
   \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (6,3.2)  -- +(5.9,0) 
     node [black,midway,above=4pt, font=\scriptsize] {Discarded steps};
   \fill[orange!10] (6,0.5) rectangle ++ (5.9,2.5);

   \draw[lightgray, line width=10pt](9,1.1) -- +(2.9,0);
   \draw[green, line width=10pt](11.9,1.1) -- +(.1,0);

     \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
      xtick=\empty,ytick=\empty, width=13.5cm,height=5.5cm,
      xmin=0, xmax=5, ymin=0, ymax=5, ylabel=$BA$, ylabel near ticks
      ]
       \addplot[mark=*,mark size=1pt,color=blue!70!black,samples=200] {2+sin(2*x)};
    \end{axis}

    \node[font=\scriptsize, above] at (13.5, 3.2) {Recorded "equilibrium"};
    \draw[green, ->] plot [smooth, tension=1] coordinates {(12, 1.1) (12.2,1.8) (13.1, 2.5) (13.5, 3.2)};
            
     \end{scope}

\end{tikzpicture}
```

The output of a simulation is a data.frame in long format (according to tidyverse style). This is very helpful to filter the output and plot it with `{ggplot2}`. Variables exported are the basal area per species `BAsp`, `n` and `h`the number of alive and harvested individuals per mesh, and `N` and `H` the total number of alive and harvested individuals in the forest (per hectar).

```{r sp1plot}
Picea_sim  %>%
    dplyr::filter(var == "BAsp", ! equil) %>%
    ggplot(aes(x = time, y = value)) +
    geom_line(linewidth = .4) + ylab("BA")
```

If size distributions needs to be extracted, it can be easily done with `{dplyr}` functions. The equilibrium step is associated with a logical variable to extract it.

```{r sp1head}
head(Picea_sim)

# get the maximum time
max_t <- max(Picea_sim$time)

# Filter example to extract the size distribution
Picea_sim %>% 
    dplyr::filter(grepl("n", var), time == max_t) %>% 
    dplyr::select(size, value)
```

# Customizing the simulations

The above simulation is one of the simplest we can produce with this package. This chapter will describe some basic customization we can add before running a simulation. 

## Initialisation step

By default, the initialization of the population run random process to draw a size distribution for each species. We already show a function (`def_initBA()`) that scale this distribution to a given basal area. However, for a basal area value, multiple distribution are possible. To control the exact distribution at start, we use `def_init_k()`. This choice of starting distribution can be used to reproduce simulations, starting from an equilibrium or a post disturbance state.

Here is an example where we start from $t = 150$ of the previous simulation. This will illustrate that despite the simulation said it reached equilibrium at time $t =$ `r (prev_equil <- max(Picea_sim$time))`, our parameters have introduced failed to identify the true equilibrium. The previous equilibrium is highlighted in blue rectangle.

```{r sp1initk}
# extract distribution
distrib_t150 <- Picea_sim %>% 
    dplyr::filter(var == "n", time == 150) %>%
    dplyr::pull(value)

Picea_sp$init_pop <- def_init_k(distrib_t150)

Picea_sim_k <- sim_deter_forest(
    forest(species = list(Picea = Picea_sp)), 
    tlim = 200, 
    equil_time = 300, equil_dist = 50,
    SurfEch = 0.03,
    verbose = TRUE
)

Picea_sim_k  %>%
    dplyr::filter(var == "BAsp", ! equil) %>%
    # below, we keep the time reference of the previous simulation 
    # to simplify the understanding of the full document.
    dplyr::mutate(time = time + 150) %>% 
    ggplot(aes(x = time, y = value)) +
    geom_line(linewidth = .4) + ylab("BA") +
    geom_rect(mapping = aes(xmin = max(Picea_sim$time) - 50, xmax = max(Picea_sim$time), 
                                     ymin = max(value-1), ymax = max(value)),
                       alpha = 0.002, fill = "blue") +
    geom_text(aes(label = "False equilibrium", 
                  x = max(Picea_sim$time) - 25, y = max(value) - 3), size = 4) 
```


## Recruitment delay

New trees are recruited at $L$ (90mm). Trees takes however several years to grow from seed to the minimum size $L$. To represents the time lag for a tree to recruit up to $L$, we can modify a species by adding a delay for recruitment of new individuals. By default, the recruitment is a given number of new individuals. This number is split in half and adds to the first two class of size distribution. Adding delay expand the IPM with `n_delay` age based classes to represent the year its takes (here 5 years) to grow up to $L$. The new recruit will age from one age class to another until they enter the size-based IPM.

A default delay is used by `build_IPM()` for each species. These values are computed from regressions.

```{r, Forest_delay}
n_delay <- 5
Picea_sp_d5 <- delay(Picea_sp, n_delay)
Picea_sp_d5$info["species"] <- "Picea_delayed" # We rename the species for easier plot.
Picea_sp_d5$init_pop <- def_initBA(30)
```

The simulation is run in the same way as the delay is only defined at the IPM level.

```{r, delay_sim}
set.seed(42)
Picea_sim_d5 <- sim_deter_forest(
    forest(species = list(Picea = Picea_sp_d5)),
    tlim = 200, 
    equil_time = 200, equil_dist = 50,
    SurfEch = 0.03,
    verbose = TRUE
)
```

Equilibrium BA should be really close with or without delay ($\Delta_{BA} < 1$). N is expected to increase with delay since delayed mesh cell with seeds are counted in.

```{r, delay_plot}
Picea_sim_d5 %>%
    rbind(Picea_sim) %>%
    dplyr::filter(var %in% c("BAsp", "N"), !equil) %>%
    ggplot(aes(x = time, y = value, color = species)) +
    geom_line(linetype = "dashed", linewidth = .3) +
    geom_point(size = .7) + ylab("BA") +
    facet_wrap(~ var, scales = "free_y") +
    NULL
```

<!-- Despite a really close BA, the size distribution is different at equilibrium. -->

<!-- *Note : Values below the redline does not count in BA computation.* -->

```{r, delay_dist_plot, eval = FALSE, echo = FALSE}

base_d <- as.numeric(Picea_sp$IPM$info["delay"])
Picea_sim_d5 %>%
    dplyr::mutate(mesh = mesh - n_delay) %>%
    rbind(Picea_sim) %>%
    dplyr::filter(var == "n", time == 200) %>%
    ggplot(aes(x = mesh, y = value, fill = species, group = desc(species))) +
    geom_area(color = "black", alpha = .3, size = .4, position = "identity") +
    geom_vline(xintercept = base_d, color = "red") +
    annotate("text", x = base_d / 2, y = .5, label = "Minimal mesh in BA", size = 3,
             color = "red", angle = 90, ) +
    NULL
```


## Multiple species

Multi-specific simulations are performed like the simulations previously illustrated. The only difference is in the construction of the forest object. This explain why the `species` argument for `forest()` function require a list for input. 

We need to modelise a second species. Be careful to select the same climate as the first species.

```{r second_species}
data("fit_Abies_alba")

Abies_ipm <- make_IPM(
    species = "Abies_alba", 
    climate = climate, # this variable is defined at the top of the doc.
    fit = fit_Abies_alba,
    clim_lab = "optimum clim",
    mesh = c(m = 700, L = 90, U = get_maxdbh(fit_Abies_alba) * 1.1),
    BA = 0:60, # Default values are 0:200, smaller values speed up this vignette.
    verbose = TRUE
)
Abies_sp <- species(IPM = Abies_ipm, init_pop = def_initBA(35))
```

```{r nsp_forest}
# We edit back the init_fun for Picea
Picea_sp$init_pop <- def_initBA(15)
Picea_Abies_for <- forest(species = list(Picea = Picea_sp, Abies = Abies_sp))

set.seed(42)
Picea_Abies_sim <- sim_deter_forest(
    Picea_Abies_for, 
    tlim = 500, 
    equil_time = 500, equil_dist = 50,
    SurfEch = 0.03,
    verbose = TRUE
)

Picea_Abies_sim  %>%
    dplyr::filter(var == "BAsp", ! equil) %>%
    ggplot(aes(x = time, y = value, color = species)) +
    geom_line(linewidth = .4) + ylab("BA") +
    stat_summary(fun = "sum",  aes(col="Total"),
                 geom ='line', linetype = "dashed", linewidth = .3)
```

An important point to note for multiple species is that the equilibrium is defined at the forest level, that is the sum of species basal area.

```{r nsp_timeline, engine='tikz', fig.ext = 'png', engine.opts = list(template = "tex/tikz.tex"), echo = FALSE, eval=TRUE}
\begin{tikzpicture}

   \draw[ultra thick] (0,3) node[above= 25pt] {\textbf{A}};
   \draw[ultra thick] (0,0) -- (12,0);
   \draw[ultra thick] (0,0) node[below=3pt,thick] {$t_0$};
   \draw[ultra thick] (2,0) node[below=3pt,thick] {$t_1$};
   \draw[ultra thick] (6,0) node[below=3pt, thick] {$t_{lim}$};
   \draw[ultra thick] (12,0) node[below=3pt, thick] {$t_{equil\_time}$};
   \draw[ultra thick] (10,0) node[below=3pt, thick] {$t_{equil}$} {};
   \foreach \x in {0, 2,4,6,8,10, 12}
   \draw (\x cm,3pt) -- (\x cm,-3pt);

   \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (0,3.2)  -- +(6,0) 
      node [black,midway,above=4pt, font=\scriptsize] {Mandatory recorded simulation};
  %  \draw[green, line width=40pt](0,1) -- +(6,0);
   \fill[blue!10] (0,0.5) rectangle ++(6,2.5);
   \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (6,3.2)  -- +(3.9,0) 
        node [black,midway,above=4pt, font=\scriptsize] {Discarded steps};
   \fill[orange!10] (6,0.5) rectangle ++ (3.9,2.5);

     \draw [thick ,decorate,decoration={brace,amplitude=5pt}] (7,1.9)  -- +(2.9,0) 
        node [black,midway,above=4pt, font=\scriptsize] {equil\_dist};
     \draw[lightgray, line width=8pt](7,1.7) -- +(2.9,0);
     \draw [thick ,decorate,decoration={brace,amplitude=5pt, mirror}] (9.9,1.5)  -- +(0,0.4) 
     node [black,midway,right=4pt, font=\scriptsize] {equil\_diff};
     \draw[green, line width=12pt](9.9,0.85) -- +(.1,0);

   \begin{axis}[trig format plots=rad,axis lines=left, axis line style={draw=none},
     xtick=\empty,ytick=\empty, width=11.5cm,height=5.5cm,
     xmin=0, xmax=4, ymin=0, ymax=6, ylabel=$BA$, ylabel near ticks
     ]
 \addplot[mark=*,mark size=1pt,color=cyan!70!black,samples=200] {1.4+exp(-x)*sin(10*x)};
 \addplot[mark=*,mark size=1pt,color=teal!70!black,samples=200] {1.2+exp(-x)*sin(2*x)};
 \addplot[mark=*,mark size=1pt,color=gray!70!black,samples=200] {1.4+exp(-x)*sin(10*x)+ 1.2+exp(-x)*sin(2*x)};
 \end{axis}
 
 \node[font=\scriptsize, above] at (11.5, 3.2) {Recorded equilibrium};
 \draw[green, ->] plot [smooth, tension=1] coordinates {(10, 0.95) (10.2,1.2) (11.8, 1.7) (11.5, 3.2)};
 
 \node[font=\scriptsize, right, cyan!70!black] at (10, 0.97) {species A};
 \node[font=\scriptsize, right, teal!70!black] at (10, 0.78) {species B};
 \node[font=\scriptsize, align= center, above, gray!70!black] at (5, 1.8) {Species sum \\(not recorded)};


\end{tikzpicture}
```

# References

<div id="refs"></div>