---
title: "1 - Running a simple simulation"
data: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{1 - Running a simple simulation}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


```{r, results="hide", message=FALSE}
library(landsepi)
```

## General presentation of the package

See `?lansdepi` for a complete description of the model, assumptions and available functions. 
See also the vignette [list of parameters](list_of_parameters.pdf) for a detailed description of all model parameters.  
Take a [quick overview](landsepi_poster.pdf) of what landsepi can do.  

## Initialisation of simulation parameters

The function `createSimulParams()` will create a directory to store simulation outputs, 
and return an object of class `LandsepiParams` that will further contain all simulation parameters. 
```{r, results="hide", message=FALSE}
simul_params <- createSimulParams(outputDir = getwd())
```

In the following, the object `simul_params` is to be updated with all simulation parameters using `set*()` functions. 
For some specific parameters already set using `set*()` functions, they may be next updated using `update*()` functions.
To help parameterisation, built-in parameters are available and may be loaded via `load*()` functions.  

To avoid dependency issues between functions, we recommend to parameterise the simulation with the following order 
(**functions in bold are crucial to run a simulation**, others are optional):  
- **createSimulParams()**  
- setSeed()  
- **setTime()**  
- setPathogen()  
- updateReproSexProb()  
- **setLandscape()**  
- setDispersal()  
- setGenes()  
- **setCultivars()**  
- allocateCultivarGenes()  
- **allocateCroptypeCultivars()**  
- **setCroptypes()**  
- **allocateLandscapeCroptypes()**  
- setInoculum()  
- updateSurvivalProb()  
- setTreatments()  
- setOutputs()  
- checkSimulParams()  
- **runSimul()**  


## Setting the seed and time parameters

A seed (for random number generator) is randomly generated when `simul_params` is initialised by 
`createSimulParams()`. If a specific seed is required, it can be set using the function `setSeed()`:
```{r}
simul_params@Seed
simul_params <- setSeed(simul_params, seed = 1)
simul_params@Seed
```

The number of cropping seasons to simulate (e.g. 6 years) and the number of time steps per 
cropping season (e.g. 120 days/year) can be set using `setTime()`:
```{r}
simul_params <- setTime(simul_params, Nyears = 6, nTSpY = 120)
simul_params@TimeParam
```


## Setting pathogen parameters

Pathogen parameters must be stored in a list of aggressiveness components defined on a susceptible host for a 
pathogen genotype not adapted to resistance.  

Buit-in parameterisation for different pathogens (e.g. rusts of cereal crops, mildew of grapevine, 
black sigatoka of banana) are available using function `loadPathogen()`:
```{r}
basic_patho_param <- loadPathogen(disease = "rust")
basic_patho_param
basic_patho_param <- loadPathogen(disease = "mildew")
basic_patho_param
basic_patho_param <- loadPathogen(disease = "sigatoka")
basic_patho_param
```

This list may be updated to change a specific parameter. For instance, to change the infection rate to 50%:
```{r}
basic_patho_param <- loadPathogen("rust")
basic_patho_param$infection_rate <- 0.5
basic_patho_param
```

Alternatively, the list may be generated manually to control all parameters: 
```{r}
basic_patho_param <- list(infection_rate = 0.4
                          , latent_period_mean = 10
                          , latent_period_var = 9
                          , propagule_prod_rate = 3.125
                          , infectious_period_mean = 24
                          , infectious_period_var = 105
                          , survival_prob = 1e-4
                          , repro_sex_prob = 0
                          , sigmoid_kappa = 5.333
                          , sigmoid_sigma = 3
                          , sigmoid_plateau = 1
                          , sex_propagule_viability_limit = 1
                          , sex_propagule_release_mean = 1
                          , clonal_propagule_gradual_release = 0)
```

Then, the list is used to fuel the object `simul_params` via the function `setPathogen()`:
```{r}
simul_params <- setPathogen(simul_params, patho_params = basic_patho_param)
simul_params@Pathogen
```


## Setting the landscape and pathogen dispersal

The landscape (i.e. boundaries of fields) must be in shapefile format. Five built-in landscapes of about 150 
fields are available using the function `loadLandscape()`:
```{r, fig.alt="The first Landscape available"}
landscape <- loadLandscape(id = 1)
length(landscape)
plot(landscape, main = "Landscape structure")
```

*See also tutorial on how to [parameterise landscape and dispersal](O3-landscape_dispersal.html) 
to use your own landscape and compute your own dispersal matrices.*

Dispersal is given by a vectorised matrix giving the probability of dispersal from any field of the 
landscape to any other field. The size of the matrix must be the square of the number of fields in the landscape. 
It is thus specific to both the pathogen and the landscape. For rusts pathogens, 
a built-in dispersal matrix is available for each landscape using the function `loadDispersalPathogen()`:
```{r}
disp_patho_clonal <- loadDispersalPathogen(id = 1)[[1]]
head(disp_patho_clonal)
length(landscape)^2 == length(disp_patho_clonal)
```

Then, the object `simul_params` is updated with the landscape and dispersal matrix via the 
functions `setLandscape()` and `setDispersalPathogen()`, respectively:
```{r}
simul_params <- setLandscape(simul_params, land = landscape)
simul_params <- setDispersalPathogen(simul_params, disp_patho_clonal)
```

## Setting croptypes, cultivars and resistance genes

Fields of the landscape are cultivated with different croptypes that can rotate through time; 
each croptype is composed of a pure cultivar or a mixture; and each cultivar may carry one or 
several resistance genes.

### Cultivars

Characteristics of each host genotype (i.e. cultivar) are summarised in a dataframe, which 
contains parameters representing the cultivar as if it was cultivated in a pure crop. 
**Note that the name of the cultivar cannot accept spaces.**  
A buit-in parameterisation for classical types of cultivars is available using `loadCultivar()` 
to generate each line of the dataframe. Type "nonCrop" allows the simulation of 
forest, fallows, etc. i.e. everything that is not planted, does not cost anything and does not 
yield anything (with regard to the subject of the study).  
**It is advised to implement a susceptible cultivar at first line of the dataframe to allow 
pathogen initial inoculation in default parameterisation.**
```{r}
cultivar1 <- loadCultivar(name = "Susceptible", type = "wheat")
cultivar2 <- loadCultivar(name = "Resistant1", type = "wheat")
cultivar3 <- loadCultivar(name = "Resistant2", type = "banana")
cultivar4 <- loadCultivar(name = "Resistant3", type = "pepper")
cultivar5 <- loadCultivar(name = "Forest", type = "nonCrop")
cultivars <- data.frame(rbind(cultivar1, cultivar2, cultivar3, cultivar4, cultivar5)
                        , stringsAsFactors = FALSE)
cultivars
```

Similarly as pathogen parameters, characteristics of the cultivars may be updated as required. 
For example, to change the growth rate of the susceptible cultivar:
```{r}
cultivars[cultivars$cultivarName == "Susceptible", "growth_rate"] <- 0.2
cultivars
```

Finally, the dataframe `cultivars` can also be generated entirely from scratch:
```{r}
cultivars_new <- data.frame(cultivarName = c("Susceptible", "Resistant"),
                            initial_density =   c(0.1, 0.2),
                            max_density =       c(2.0, 3.0),
                            growth_rate =       c(0.1, 0.2),
                            reproduction_rate = c(0.0, 0.0),
                            yield_H =           c(2.5, 2.0),
                            yield_L =           c(0.0, 0.0),
                            yield_I =           c(0.0, 0.0),
                            yield_R =           c(0.0, 0.0),
                            planting_cost =   c(225, 300),
                            market_value =      c(200, 150),
                            stringsAsFactors = FALSE)
cultivars_new
```

Note, assuming that only healthy hosts (state H) contribute to host growth, the production of healthy biomass between $t$ and $t+1$ is computed using the following logistic equation:
$$H_{t+1} = H_{t} \times  \left[1 + growth\_rate \times (1-\frac{N_{t}}{K})\right]$$
with $N_t$ the total number of individual hosts at time-step $t$ and $K$ the carrying capacity.


### Resistance genes

This part can be skipped if no resistance gene is to be simulated.
Characteristics of each plant resistance gene and each corresponding pathogenicity gene are summarised in a dataframe.  
A built-in parameterisation for classical resistance sources is available using `loadGene()` to generate each line of the 
dataframe. Type "majorGene" is for completely efficient resistance to which pathogen may adapt via a single mutation.  
Type "APR" stands for Adult Plant Resistance, i.e. a major gene which activates after a delay after planting, and type "QTL" is a 
partially efficient resistance source to which the pathogen may adapt gradually. The type "immunity" code for a completely 
efficient resistance gene that cannot be overcome, in such a way that infection is totally impossible; it is helpful 
to parameterise non-host cultivars. 
```{r}
gene1 <- loadGene(name = "MG 1", type = "majorGene")
gene2 <- loadGene(name = "Lr34", type = "APR")
gene3 <- loadGene(name = "gene 3", type = "QTL")
gene4 <- loadGene(name = "nonhost resistance", type = "immunity")
genes <- data.frame(rbind(gene1, gene2, gene3, gene4), stringsAsFactors = FALSE)
genes
```

Similarly as pathogen parameters, characteristics of the genes may be updated as required. For example, to change the 
mutation probability of the pathogen with regard to its adaptation to "MG 1":
```{r}
genes[genes$geneName == "MG 1", "mutation_prob"] <- 1e-3
genes
```

Alternatively, the dataframe `genes` can also be generated entirely from scratch:
```{r}
genes_new <- data.frame(geneName =               c("MG1", "MG2"),
                        efficiency =             c(1.0  , 0.8  ),
                        age_of_activ_mean =      c(0.0  , 0.0  ),
                        age_of_activ_var =       c(0.0  , 0.0  ),
                        mutation_prob =          c(1E-7 , 1E-4),
                        Nlevels_aggressiveness = c(2    , 2    ),
                        adaptation_cost =        c(0.50 , 0.75 ),
                        relative_advantage =     c(0.50 , 0.75 ),
                        tradeoff_strength =      c(1.0  , 1.0  ),
                        target_trait =           c("IR" , "LAT"),
                        recombination_sd =       c(1.0,1.0),
                        stringsAsFactors = FALSE)
genes_new
```

Note that the column "efficiency" refers to the percentage of reduction of the targeted aggressiveness 
component on hosts carrying the gene. 
For example, a 80% efficient resistance against infection rate (target_trait = "IR") 
means that the infection rate of a non-adapted pathogen is reduced by 80% on a resistant host compared 
to what it is on a susceptible host. For resistances targeting the latent period, 
the percentage of reduction is applied to the inverse of the latent period in such a way that latent 
period is higher in resistant than in susceptible hosts.


### Allocating genes to cultivars

The object `simul_params` can be updated with `setGenes()` and `setCultivars()`:
```{r}
simul_params <- setGenes(simul_params, dfGenes = genes)
simul_params <- setCultivars(simul_params, dfCultivars = cultivars)
simul_params@Genes
simul_params@Cultivars
```

Then the function `allocateCultivarGenes()` allows the attribution of resistance genes to cultivars:
```{r}
simul_params <- allocateCultivarGenes(simul_params
                                      , cultivarName = "Resistant1"
                                      , listGenesNames = c("MG 1"))
simul_params <- allocateCultivarGenes(simul_params
                                      , cultivarName = "Resistant2"
                                      , listGenesNames = c("Lr34", "gene 3"))
simul_params <- allocateCultivarGenes(simul_params
                                      , cultivarName = "Resistant3"
                                      , listGenesNames = c("nonhost resistance"))
simul_params@CultivarsGenes
```

With this example of parameterisation:  
- "Susceptible" is a susceptible cultivar (initially infected by a wild-type pathogen)  
- "Resistant1" is a mono-resistant cultivar  
- "Resistant2" is a pyramided cultivar  
- "Resistant3" is a nonhost cultivar  
- "Forest" is not a crop.  
Infection is impossible on both "Resistance3" and "Forest", but the former will 
be considered for host growth, yield and planting costs, whereas the latter won't.


### Allocating cultivars to croptypes

Characteristics of each croptype (a croptype is a set of hosts cultivated in a field with specific proportions) 
are summarised in a dataframe.  
A buit-in parameterisation for classical croptypes is available 
using `loadCroptypes()` to generate the whole table filled with zeros: 
```{r}
croptypes <- loadCroptypes(simul_params, names = c("Susceptible crop"
                                                   , "Pure resistant crop"
                                                   , "Mixture"
                                                   , "Other"))
croptypes
```

Then `croptypes` is updated by `allocateCroptypeCultivars()` 
to specify the composition of every croptype with regard to cultivars (and proportion for mixtures):
```{r}
croptypes <- allocateCroptypeCultivars(croptypes
                                       , croptypeName = "Susceptible crop"
                                       , cultivarsInCroptype = "Susceptible")
croptypes <- allocateCroptypeCultivars(croptypes
                                       , croptypeName = "Pure resistant crop"
                                       , cultivarsInCroptype = "Resistant1")
croptypes <- allocateCroptypeCultivars(croptypes
                                       , croptypeName = "Mixture"
                                       , cultivarsInCroptype = c("Resistant2","Resistant3")
                                       , prop = c(0.4, 0.6))
croptypes <- allocateCroptypeCultivars(croptypes
                                       , croptypeName = "Other"
                                       , cultivarsInCroptype = "Forest")
croptypes
```

Finally the object `simul_params` is updated using `setCroptypes()`:
```{r}
simul_params <- setCroptypes(simul_params, dfCroptypes = croptypes)
simul_params@Croptypes
```

Alternatively, the dataframe `croptypes` can be generated from scratch:
```{r}
croptypes <- data.frame(croptypeID = c(0, 1, 2, 3) ## must start at 0 and match with values from landscape "croptypeID" layer
                        , croptypeName = c("Susceptible crop"
                                           , "Pure resistant crop"
                                           , "Mixture"
                                           , "Other")
                        , Susceptible = c(1,0,0  ,0)
                        , Resistant1  = c(0,1,0  ,0)
                        , Resistant2  = c(0,0,0.5,0)
                        , Resistant3  = c(0,0,0.5,0)
                        , Forest      = c(0,0,0  ,1)
                        , stringsAsFactors = FALSE)
simul_params <- setCroptypes(simul_params, croptypes)
```

### Allocating croptypes to fields of the landscape

The function `allocateLandscapeCroptypes()` manages the allocation of croptypes in time (for rotations 
of different croptypes on the same fields) and space (for mosaic of fields cultivated with different croptypes). 
See `?allocateLandscapeCroptypes` for help in parameters.  
Briefly, a rotation sequence is defined by a list. Each 
element of this list is a vector containing the indices of the croptypes that are cultivated simultaneously in the 
landscape. The "rotation_period" parameter defines the duration before switching from one element of "rotation_sequence" 
to the next. For each element of "rotation_sequence" is associated a vector of proportions of each croptype in the 
landscape (parameter "prop"). The parameter "aggreg" controls the level of spatial aggregation between croptypes. 

For example, to generate a landscape whose surface is composed of 
1/3 of forests, 1/3 of susceptible crop, and 1/3 of fields where a pure resistant crop is alternated with a mixture 
every two years:
```{r}
# croptypeIDs cultivated in each element of the rotation sequence:
rotation_sequence <- list(c(0,1,3), c(0,2,3))
rotation_period <- 2  # number of years before rotation of the landscape
prop <- list(rep(1/3, 3), rep(1/3, 3)) # proportion (in surface) of each croptype 
aggreg <-1 # level of spatial aggregation
simul_params <- allocateLandscapeCroptypes(simul_params
                                           , rotation_period = rotation_period
                                           , rotation_sequence = rotation_sequence
                                           , prop = prop
                                           , aggreg = aggreg
                                           , graphic = TRUE)
# plot(simul_params@Landscape)
```


## Setting the inoculum

To set the inoculum, the function `setInoculum()` is used. Several scenarios may be simulated and are summarized below.
For the default scenario (see below), the function is simply parameterised with the probability for individual 
hosts to be infectious (i.e. at state 'I') at the beginning of the simulation (i.e. at t=0).

For more complex scenarios (i.e. to specify the location and genetic structure of the inoculum), 
the function `setInoculum()` can be used with a 3D-array of dimensions (Nhost, Npatho, Npoly) 
indicating the initial probability to be infectious, for each cultivar, pathogen genotype and polygon, respectively.
To define this array, the functions `getMatrixGenePatho()`, `getMatrixCultivarPatho()`, `getMatrixCroptypePatho()` 
and `getMatrixPolyPatho()` acknowledge which pathogen genotypes can infect which genes, cultivars, croptypes 
and polygons respectively. Each function returns a matrix indicating if there is compatibility (value of 1, 
i.e. infection is possible) or not (value of 0, i.e. the cultivar is not present, or protected by a fully efficient 
resistance gene targeting the infection rate from the beginning of the cropping season).
```{r}
getMatrixGenePatho(simul_params)
getMatrixCultivarPatho(simul_params)
getMatrixCroptypePatho(simul_params)
getMatrixPolyPatho(simul_params)[1:10,]
```

Finally, the function `loadInoculum()` helps build the inoculum array as follows:  
Let $\phi_{v,p,i}$ be the probability for an individual of cultivar $v$ to be infected by pathogen genotype $p$ 
in field $i$ at the beginning of the simulation. This probability is computed as follows:
$$\phi_{v, p, i} = \phi^0 \times \phi^{host}_{v} \times \phi^{patho}_{p} \times \phi^{poly}_{i} \times I^{present}_{v,i} \times I^{compatible}_{v,p}$$
with:  
$\phi^{host}$ the vector of probabilities for every host (parameter `pI0_host`),  
$\phi^{patho}$ the vector of probabilities for every pathogen genotype (parameter `pI0_patho`),  
$\phi^{poly}$ the vector of probabilities for every polygon (parameter `pI0_poly`),  
$\phi^0$ a multiplicative constant (parameter `pI0_all`).  
$I^{present}_{v,i}$ a binary variable equal to 1 if cultivar $v$ is grown in field $i$ (and 0 otherwise),  
$I^{compatible}_{v,p}$ a binary variable equal to 1 if pathogen genotype $p$ can infect cultivar $v$ at the beginning of 
the cropping season (and 0 otherwise). 


### Scenario 1 (default). Only the avirulent pathogen on the susceptible cultivar in all fields (global inoculum)
The default scenario is the presence of a pathogen genotype not adapted to any resistance gene, 
and present in all fields of the landscape where a susceptible cultivar is grown. 
In this situation, it is important that the susceptible cultivar is entered at the first line of the table cultivars.
Then, one can simply use:
```{r}
# Option 1: simply use the default parameterisation
simul_params <- setInoculum(simul_params, 5E-4)
 
# Option 2: use loadInoculum()
Npatho <- prod(simul_params@Genes$Nlevels_aggressiveness)
Nhost <- nrow(simul_params@Cultivars)
pI0 <- loadInoculum(simul_params, pI0_all=5E-4, pI0_host=c(1,rep(0, Nhost-1)), pI0_patho=c(1,rep(0, Npatho-1)))
simul_params <- setInoculum(simul_params, pI0)
```

### Scenario 2. Only the avirulent pathogen on the susceptible cultivar in only some fields (local inoculum)
To specify the location of the inoculum, the parameters `pI0_host`, `pI0_patho` and `pI0_poly` are filled with values of 1 and 0 indicating 
if the cultivar, pathogen genotype and field are inoculated or not. The probability is given by the constant `pI0_all`. 
In this example, 5 fields in the landscape are randomly chosen among those cultivated with the susceptible cultivar:
```{r, eval=FALSE}
Npatho <- prod(simul_params@Genes$Nlevels_aggressiveness)  ## Nb of pathogen genotypes
Nhost <- nrow(simul_params@Cultivars)  ## Nb of cultivars
Npoly <- nrow(simul_params@Landscape)  ## Nb of polygons in the landscape
Npoly_inoc <- 5  ## number of inoculated polygons
compatible_poly <- getMatrixPolyPatho(simul_params)[,1]  ## whether the avr pathogen can infect the polygons
id_poly <- sample(grep(1, compatible_poly), Npoly_inoc)  ## random polygon picked among compatible ones
pI0_poly <- as.numeric(1:Npoly %in% id_poly)  
pI0 <- loadInoculum(simul_params, pI0_all=5E-4, pI0_host=c(1,rep(0, Nhost-1)), pI0_patho=c(1,rep(0, Npatho-1)), 
                    pI0_poly=pI0_poly)
simul_params <- setInoculum(simul_params, pI0)
```

### Scenario 3. A diversity of pathogen genotypes in the inoculum, in all fields (global inoculum)
This scenario matches with situations where several pathogen genotypes (including those adapted to resistance) 
are initially present in the landscape at the beginning of the simulation. 
In this example, different probabilities are given to the different pathogen genotypes 
using the parameters `pI0_patho` (the same rationale can be used for different probabilities on the different hosts using 
`pI0_host`).
```{r, eval=FALSE}
## Example with 4 pathogen genotypes and 2 cultivars
pI0 <- loadInoculum(simul_params, pI0_patho=c(1E-3,1E-4,1E-4,1E-5), pI0_host=c(1,1))
simul_params <- setInoculum(simul_params, pI0)
```

### Scenario 4. A diversity of pathogen genotypes in the inoculum, in some fields only (local inoculum)
Here, the example is similar as the previous one but only 5 fields are inoculated:
```{r, eval=FALSE}
Npoly <- nrow(simul_params@Landscape)
Npoly_inoc <- 5  ## number of inoculated polygons 
id_poly <- sample(1:Npoly, Npoly_inoc)  ## random polygon 
pI0_poly <- as.numeric(1:Npoly %in% id_poly) 
pI0 <- loadInoculum(simul_params, pI0_patho=c(1E-3,1E-4,1E-4,1E-5), pI0_host=c(1,1), pI0_poly=pI0_poly)
simul_params <- setInoculum(simul_params, pI0)
```

### Scenario 5. Custom 3D array to define the inoculum
At last, if one wants to run a simulation with a custom inoculum on the whole landscape, it can simply be entered in 
the function `setInoculum()` as an array of dimension (Nhost, Npatho, Npoly). 
Note that in this situation, host individuals may be infected regardless the resistance 
gene they carry. 
```{r, eval=FALSE}
## example with 2 cultivars, 4 pathogen genotypes and 5 fields
Nhost=2
Npatho=4
Npoly=5
pI0 <- array(data = 1:40 / 100, dim = c(Nhost, Npatho, Npoly))
simul_params <- setInoculum(simul_params, pI0)
```

To generate an array that accounts for the fact that (i) the cultivars are not grown in all polygons, and 
(ii) cultivars may carry a resistance gene that prevent initial infection by some pathogen genotypes, one can 
use the function `loadInoculum()` as follows:
```{r, eval=FALSE}
corrected_pI0 <- loadInoculum(simul_params, pI0_mat=pI0)
simul_params <- setInoculum(simul_params, corrected_pI0)
```

### Visualization

The 3D-array inoculum can be vizualised using the function `inoculumToMatrix()`:
```{r}
inoculumToMatrix(simul_params)[,,1:5]
```


## Updating the survival probability during the off-season

At the end of each cropping season, pathogens experience a bottleneck representing the
off-season. The probability for a pathogen propagule to survive the off-season depends on the capacity 
of the "green bridge" to host the pathogen. This green bridge can, for example, be a wild reservoir or 
volunteer plants remaining in the field (e.g. owing to seedlings or incomplete harvest). It is assumed 
that the probability of survival is the same every year and in every polygon, but this assumption can 
be relaxed by updating the probability with a matrix indicating a probability for every croptype and 
every year. This simulates different management strategies or agronomic practices during the off-season. 
The function `updateSurvivalProb()` creates this matrix as follows:  
  
Let $\lambda_{y,c}$ be the probability for a pathogen propagule to survive the off-season between year 
$y$ and year $y+1$ in a polygon cultivated with croptype $c$. Unless the matrix is directly entered 
via the parameter `mat`, it is computed by $\lambda_{y,c} = \lambda_{y} \times \lambda_{c}$ with:  
$\lambda_{y}$ the vector of probabilities for every year (parameter `mat_year`),  
$\lambda_{c}$ the vector of probabilities for every croptype (parameter `mat_croptype`).  
```{r}
Ncroptypes <- nrow(simul_params@Croptypes)
Nyears <- simul_params@TimeParam$Nyears
## Same probability in every croptype:
simul_params <- updateSurvivalProb(simul_params, mat_year=1:Nyears/100)
simul_params@Pathogen
## Same probability every year:
simul_params <- updateSurvivalProb(simul_params, mat_croptype=1:Ncroptypes/10)
simul_params@Pathogen
## specific probability for different croptypes and years:
simul_params <- updateSurvivalProb(simul_params, mat_year=1:Nyears/100, mat_croptype=1:Ncroptypes/10)
simul_params@Pathogen
## One probability per year and per croptype:
simul_params <- updateSurvivalProb(simul_params, mat=matrix(runif(Nyears*Ncroptypes), ncol=Ncroptypes))
simul_params@Pathogen
```


### Visualization

The matrix of survival probabilities for every year and croptype, as well as for every polygon and year, 
can be vizualised using the function `survivalProbToMatrix()`:
```{r}
survivalProbToMatrix(simul_params)
```


## Setting chemical treatments

This part can be skipped if no chemical treatment is to be simulated.

Cultivars may be treated with chemicals which reduce the pathogen infection rate. 
Treatment efficiency is maximum (i.e. equal to the parameter `treatment_efficiency`)
at the time of treatment application (noted $t^*$); then it decreases with time (i.e. natural pesticide degradation) 
and host growth (i.e. new biomass is not protected by treatments): 

$$IR(t) = IR_{max} \times  \left(1 - \frac{treatment\_efficiency}{1+e^{4.0- 8.5 \times C(t)}}\right)$$
$C(t) = C_1 \times C_2$ is the treatment concentration at $t$, which depends on:

* the timelag $\Delta t = t - t^*$ passed since the time of treatment application $t^*$ 
(i.e. natural pesticide degradation via the parameter `treatment_degradation_rate`)

$$C_1 = e^{-treatment\_degradation\_rate \times \Delta t}$$

* the newly produced host biomass, which is not protected by treatments ($N(t^*)$ and $N(t)$ are host biomass at the 
time of treatment application and at time $t$, respectively)

$$C_2 = min(1.0, N(t^*)/N(t))$$
Cultivars to be treated with chemicals and dates of (possible) applications are defined with parameters `treatment_cultivars` and `treatment_timesteps`, and are the same for all polygons cultivated with the cultivars to be treated. 
However, the chemicals are applied in a polygon only if the disease severity (i.e. $I(t^*) / N(t^*)$) at the application date is higher than a given threshold, defined by `treatment_application_threshold`.
The cost of a single treatment application (monetary units/ha) is defined by `treatment_cost` and will impact the economic outputs.
Treatment parameters can be loaded via the function `loadTreatment` and next set via the function `setTreatment`:

```{r}
## Via function loadTreatment:
treatment <- loadTreatment(disease="mildew")
## Direct implementation:
treatment <- list(treatment_degradation_rate = 0.1,
                  treatment_efficiency = 0.8,
                  treatment_timesteps =  seq(1,120,14) ,
                  treatment_cultivars  = c(0),
                  treatment_cost = 0,
                  treatment_application_threshold = c(0.0))
## Set the parameters
simul_params <- setTreatment(simul_params, treatment)
```

## Choosing output variables

Several epidemiological, evolutionary and economic outputs can be generated by landsepi and represented in 
text files (`writeTXT=TRUE`), graphics (`graphic=TRUE`) and a video (`videoMP4=TRUE`).  
See equation below, as well as `?epid_output` and `?evol_output` for details on the different output variables.  

>Possible epidemiological outputs include:  
- **"audpc"** : Area Under Disease Progress Curve (average number of diseased hosts per square meter, $AUDPC_{v,y}$)  
- **"audpc_rel"** : relative Area Under Disease Progress Curve (average proportion of diseased hosts 
relative to the total number of existing hosts, $AUDPC_{v,y}^r$)  
- **"gla"** : Green Leaf Area (average number of healthy hosts per square meter, $GLA_{v,y}$)  
- **"gla_rel"** : relative Green Leaf Area  (average proportion of healthy hosts relative to the 
total number of existing hosts, $GLA_{v,y}^r$)  
- **"eco_yield"** : total crop yield (in weight or volume units per ha, $Yield_{v,y}$)  
- **"eco_cost"** : operational crop costs (in monetary units per ha, $Operational\_cost_{v,y}$)   
- **"eco_product"** : total crop products (in monetary units per ha, $Product_{v,y}$)   
- **"eco_margin"** : margin (products - operational costs, in monetary units per ha, $Margin_{v,y}$)  
- **"contrib"**: contribution of pathogen genotypes to LIR dynamics ($Contrib_{p,v,y}$) 
- **"HLIR_dynamics", "H_dynamics", "L_dynamics", "IR_dynamics", "HLI_dynamics"**, etc.: Epidemic dynamics 
related to the specified sanitary status (H, L, I or R and all their combinations). Graphics only, 
works only if graphic=TRUE.  
- **"all"** : compute all these outputs (default)  
- **""** : none of these outputs will be generated.  

>Possible evolutionary outputs, based on the computation of genotype frequencies ($freq(I)_{p,t}$) include:  
- **"evol_patho"**: evolution of pathogen genotypes  
- **"evol_aggr"**: evolution of pathogen aggressiveness (i.e. phenotype)  
- **"durability"**: durability of resistance genes   
- **"all"**: compute all these outputs (default)  
- **""**: none of these outputs will be generated.  

### Equations

In the following equations, $H_{i,v,t}$, $L_{i,v,p,t}$, $I_{i,v,p,t}$ and $R_{i,v,p,t}$ respectively denote the number of healthy, latent, infectious and removed host individuals, respectively, in field $i$ ($i=1,…,J$), for cultivar $v$ ($v=1,…,V$) at time step $t$ ($t=1,…,T \times Y$ with $Y$ the total number of cropping seasons and $T$ the number of time-steps per season).

$AUDPC_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} \{I_{i,v,p,t}+R_{i,v,p,t}\}}{T \times \sum_{i=1}^{J} A_i}$

$AUDPC^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} \{I_{i,v,p,t}+R_{i,v,p,t}\}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} \{L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}\}\right)}$

$GLA_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{T \times \sum_{i=1}^{J} A_i}$

$GLA^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} \{L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}\}\right)}$

$Yield_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left( y_{H,v} \times H_{i,v,t} + \sum_{p=1}^{P} \{y_{L,v} \times L_{i,v,p,t} + y_{I,v} \times I_{i,v,p,t} + y_{R,v} \times R_{i,v,p,t}\} \right)} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H^*_{i,v,t}}$
with $y_{H,v}$, $y_{L,v}$, $y_{I,v}$ and $y_{R,v}$ the theoretical yield of cultivar v in pure crop (in weight or volume unit/ha/season) associated with the sanitary statuses ‘H’, ‘L’, ‘I’ and ‘R’, respectively. $H^*_{i,v,t}$ is the number of healthy hosts in a pure crop and in absence of disease. 

$Operational\_cost_{v,y}=planting\_cost_v \times \frac{\sum_{i \in \Omega_{v,y}}(area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i} + treatment\_cost \times \frac{\sum_{i \in \Omega_{v,y}}(TFI_{i,v,y} \times area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i}$
with $\Omega_{v,y}$ the set of polygon indices where cultivar $v$ is cultivated at year $y$, and $prop_{v,i}$ the proportion of cultivar $v$ in polygon $i$ (for mixtures). $TFI$ stands for the Treatment Frequency Indicator (number of treatment applications per ha).

$Product_{v,y}=yield_{v,y} \times market\_value$

$Margin_{v,y} = Product_{v,y} - Operational\_Cost_{v,y}$

$Contrib_{p,v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \{L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}\}} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} \{L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}\}}$

$freq(I)_{p,t} = \frac{\sum_{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}} {\sum_{p=1}^{P} \sum_{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}}$

With respect to evolutionary outputs, for each pathogen genotype (`evol_patho`) or phenotype (`evol_aggr`, note that different pathogen genotypes may lead to the same phenotype on a resistant host, i.e. level of aggressiveness), several computations are performed:   
- appearance: time to first appearance (as propagule);  
- R_infection: time to first true infection of a resistant host;  
- R_invasion: time to invasion, when the number of infections of resistant hosts reaches `thres_breakdown`, above which the genotype or phenotype is unlikely to go extinct.  
The value `Nyears + 1 time step` is used if the genotype or phenotype never appeared/infected/invaded.  
Durability is defined as the time to invasion of completely adapted pathogen individuals.



### Parameterisation

A list of outputs can be generated using `loadOutputs()`:
```{r}
outputlist <- loadOutputs(epid_outputs = "all", evol_outputs = "all", disease="rust")
outputlist
```
Among the elements of `outputList`, "audpc100S" is the audpc in a fully susceptible landscape (used as reference value for graphics).
If necessary, the function `compute_audpc100S` helps compute this value in a single 1-km^2 field:
```{r}
audpc100S <- compute_audpc100S("rust", "wheat", area=1E6)
audpc100S <- compute_audpc100S("mildew", "grapevine", area=1E6)
audpc100S <- compute_audpc100S("sigatoka", "banana", area=1E6, nTSpY=182)
```


Then `simul_params` can be updated via `setOutputs()`:
```{r}
simul_params <- setOutputs(simul_params, outputlist)
```

*See also tutorial on how to [run a numerical experimental design](O2-run_exp_design.html) 
to compute your own output variables and to run several simulations within an experimental design*

## Running the simulation

The functions `checkSimulParams()` and `saveDeploymentStrategy()` check simulation parameters and 
save the object `simul_params` (which contains all parameters associated with the deployment strategy) 
into a GPKG file, respectively. 
```{r, eval=FALSE}
checkSimulParams(simul_params)
simul_params <- saveDeploymentStrategy(simul_params)
```

Then, the function `runSimul()` launches the simulation. Use `?runSimul` to get all available options.
```{r, eval=FALSE}
runSimul(simul_params, graphic = TRUE, videoMP4 = FALSE)
```

```{r, include=FALSE}
system(paste("rm -rf ", simul_params@OutputDir))
```
