---
title: "01. The punctual scale model"
author: "Daniele Da Re, Sophie O. Vanwambeke, Matteo Marcantonio"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: |
  %\VignetteIndexEntry{01. The punctual scale model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\number_sections: yes
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width=10, fig.height=10, fig.asp = 0.618, out.width = "95%", fig.align = "center", fig.dpi = 150, collapse = FALSE, comment = "#")
options(rmarkdown.html_vignette.check_title = FALSE)
```

This tutorial explains step-by-step the main features of **dynamAedes** package, a unified modelling framework for invasive *Aedes* mosquitoes. Users can apply the stochastic, time-discrete and spatially-explicit population dynamical model initially developed in [Da Re et al., (2021)](https://doi.org/10.1016/j.ecoinf.2020.101180) for *Aedes aegypti* and then expanded for other three species: *Ae. albopictus*, *Ae. japonicus* and *Ae. koreicus* [Da Re et al., (2022)](https://doi.org/10.1186/s13071-022-05414-4).

The model is driven by temperature, photoperiod and intra-specific larval competition and can be applied to three different spatial scales: punctual, local and regional. These spatial scales consider different degrees of spatial complexity and data availability, by accounting for both active and passive dispersal of the modelled mosquito species as well as for the heterogeneity of input temperature data.

We will describe the model applications for *Ae. albopictus* and for all spatial scales by using a simulated temperature dataset.

At the punctual scale, the model only requires a weather station temperature time series provided as a numerical matrix (in degree Celsius). For the purpose of this tutorial, we simulate a 1-year long temperature time series.

# 1. Input data
## 1.1. Simulate temperature data with seasonal trend
We first simulate a 1-year temperature time series with seasonal trend. For the time series we consider a mean value of 16°C and standard deviation of 2°C.

```{r,  message=FALSE, warning=FALSE}
#Load packages
# simulate temperatures
library(eesim)

# modelling 
library(dynamAedes)

# data manipulation and plotting
library(ggplot2)
Sys.setlocale("LC_TIME", "en_GB.UTF-8")  
```

```{r, warning=FALSE}
ndays <- 365*1 #length of the time series in days
set.seed(123)
sim_temp <- create_sims(n_reps = 1, 
  n = ndays, 
  central = 16, 
  sd = 2, 
  exposure_type = "continuous", 
  exposure_trend = "cos1", exposure_amp = -1.0, 
  average_outcome = 12,
  outcome_trend = "cos1",
  outcome_amp = 0.8, 
  rr = 1.0055)
```

A visualisation of the distribution of temperature values and temporal trend.
```{r, warning=FALSE}
hist(sim_temp[[1]]$x, 
 xlab="Temperature (°C)", 
 main="Histogram of simulated temperatures")

plot(sim_temp[[1]]$date,
 sim_temp[[1]]$x,
 main="Simulated temperatures seasonal trend", 
 xlab="Date", ylab="Temperature (°C)"
 )
```

# 2. Modelling 
## 2.1 Model settings

We define a few model parameters which need to be defined by the user.
```{r, warning=FALSE}
## Define the day of introduction (July 1st is day 1)
str <- "2000-07-01"
## Define the end-day of life cycle (August 1st is the last day)
endr <- "2000-08-01"
## Define the number of eggs to be introduced
ie <- 1000
## Define the number of model iterations
it <- 1 # The higher the number of simulations the better
## Define the number of liters for the larval density-dependent mortality
habitat_liters <- 1
#Define latitude and longitude for the diapause process
myLat <- 42
myLon <- 7
## Define the number of parallel processes (for sequential iterations set nc=1)
cl <- 1
```

Float numbers in the temperature matrix would slow the computational speed, they must be multiplied by 1000 and then transformed in integer numbers. We also transpose the matrix from long to wide format, since we conceptualized the model structure considering the rows as the spatial component (e.g., observations; here = 1) and the columns as the temporal one (e.g., variables).

```{r, warning=FALSE}
df_temp <- data.frame("Date" = sim_temp[[1]]$date, "temp" = sim_temp[[1]]$x)
w <- t(as.integer(df_temp$temp*1000)[format(as.Date(str)+1,"%j"):format(as.Date(endr)+1,"%j")])
```

## 2.2 Run the model
Running the model with the settings specified in this example takes about 2 minutes.
```{r,  message=FALSE, warning=FALSE}
simout <- dynamAedes.m(species="albopictus", 
 scale="ws",  
 jhwv=habitat_liters,  
 temps.matrix=w,  
 startd=str, 
 endd=endr,  
 n.clusters=cl, 
 iter=it,  
 intro.eggs=ie,  
 compressed.output=TRUE, 
 lat=myLat, 
 long=myLon,
 verbose=FALSE,
 seeding=TRUE)
```

# 3. Analyse the results
A first summary of simulations can be obtained with:
```{r, warning=FALSE}
summary(simout)
```

The *simout* object is a S4 object where simulation outputs and related details are saved in different slot:

For example, the number of model iterations is saved in: 
```{r, warning=FALSE, message=FALSE, results='hide'}
simout@n_iterations
```

The simulation output is stored in:
```{r, warning=FALSE, message=FALSE, results='hide'}
simout@simulation
```

Which is a list where the the **first** level stores simulation of different iteration, while the **second**corresponds to the simulated days in the corresponding iteration.
If we inspect the first iteration, we observe that the model has computed `length(simout[[1]])` days, since we have started the simulation on the 1st of July and ended on the 1st of August.

```{r, warning=FALSE}
length(simout@simulation[[1]])
```

The **third** level corresponds to the quantity of individuals for each stage (rows) in each day for each iteration.
If we inspect the 1st and the 15th day within the first iteration, we obtain a matrix having:
```{r, warning=FALSE}
simout@simulation[[1]][[1]]
simout@simulation[[1]][[15]]
```
We can use the auxiliary functions of the package to analyse the results. 

## 3.1 Derive probability of a successfull introduction at the end of the simulated period
First, we can retrieve the "probability of successful introduction", computed as the proportion of model iterations that resulted in a viable mosquito population at a given date. In this case the results is 1, since we have only one iteration and the population is still viable at the end of the simulation.

```{r, message=FALSE, warning=FALSE}
psi(input_sim = simout, eval_date = 30)
```

## 3.2 Derive abundance 95% CI for each life stage and in each day
We now compute the interquantile range abundance for all the stages of the simulated population using the function *adci*.
```{r, message=FALSE, warning=FALSE}

# Retrieve the maximum number of simulated days
dd <- max(simout)

# Compute the inter-quartile of abundances along the iterations
breaks=c(0.25,0.50,0.75)
ed=1:dd

outdf <- rbind(
  adci(simout, eval_date=ed, breaks=breaks, stage="Eggs"),
  adci(simout, eval_date=ed, breaks=breaks, stage="Juvenile"),
  adci(simout, eval_date=ed, breaks=breaks, stage="Adults"),
  adci(simout, eval_date=ed, breaks=breaks, stage="Dia"))
```

Then we can look at the time series of the population dynamics stage by stage.
```{r, message=FALSE, warning=FALSE}
outdf$stage <- factor(outdf$stage, levels= c('Egg', 'DiapauseEgg', 'Juvenile', 'Adult'))
outdf$Date <- rep(seq.Date(as.Date(str), as.Date(endr) - 2, by="day"), 4)

ggplot(outdf, aes(x=Date, y=X50., group=factor(stage), col=factor(stage))) +
ggtitle("Ae. albopictus Interquantile range abundance") +
geom_ribbon(aes(ymin=X25., ymax=X75., fill=factor(stage)), 
  col="white", 
  alpha=0.2, 
  outline.type="full") +
geom_line(linewidth=0.8) +
labs(x="Date", y="Interquantile range abundance", col="Stage", fill="Stage") +
facet_wrap(~stage, scales = "free_y") +
theme_light() +
theme(legend.position="bottom",  
  text = element_text(size=16), 
  strip.text = element_text(face = "italic"))

```