---
title: "Getting started with odin2"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with odin2}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

`odin2` implements a high-level language for describing and implementing ordinary differential equations and difference equations in R.  It provides a "[domain specific language](https://en.wikipedia.org/wiki/Domain-specific_language)" (DSL) which _looks_ like R but is compiled directly to C++, using [`dust2`](https://mrc-ide.github.io/dust2/) to solve your system and to provide an interface to particle filters.  You can then use [`monty`](https://mrc-ide.github.io/monty/) to fit your models using MCMC.

This vignette jumps through a few of the core features of `odin2` and ways that you might use it with `dust2` and `monty`.  Other vignettes (when written!) will expand on topics covered here in more detail.

# Discrete time stochastic SIR model

A simple definition of the SIR model is:

$$\begin{align*}
\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dI}{dt} &= \beta \frac{SI}{N} - \gamma I \\
\frac{dR}{dt} &= \gamma I \\
\end{align*}$$

where $S$ is the number of susceptibles, $I$ is the number of infected and $R$ is the number recovered; the total population size $N = S + I + R$ is constant. $\beta$ is the infection rate, $\gamma$ is the recovery rate.

Discretising this model in time steps of width $dt$ gives the following update equations for each time step:

$$\begin{align*}
S_{t+1} &= S_t - n_{SI} \\
I_{t+1} &= I_t + n_{SI} - n_{IR} \\
R_{t+1} &= R_t + n_{IR}
\end{align*}$$

where

$$\begin{align*}
n_{SI} &\sim \mathrm{Binomial}(S, 1 - e^{-\beta \frac{I}{N} \cdot dt}) \\
n_{IR} &\sim \mathrm{Binomial}(I, 1 - e^{-\gamma \cdot dt})
\end{align*}$$

Here is this system, as a stochastic compartmental model:

```{r}
gen <- odin2::odin({
  p_IR <- 1 - exp(-gamma * dt)
  N <- parameter(1000)

  p_SI <- 1 - exp(-(beta * I / N * dt))
  n_SI <- Binomial(S, p_SI)
  n_IR <- Binomial(I, p_IR)

  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0

  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  I0 <- parameter(10)
})
```

This step generates C++ code for the model and compiles it; it will take a few seconds.

Once we have our system, we can pass it to [`dust2::dust_system_create`] to create and start simulating it.  Our system above has defaults for its parameters (`N`, `beta`, `gamma`, and `I0`) so we can initialise with almost no arguments:

```{r}
pars <- list(beta = 0.2, gamma = 0.1, I0 = 10, N = 1000)
sys <- dust2::dust_system_create(gen(), pars, n_particles = 10)
```

By default the system will start at time `0` and with `dt = 1`.  We can simulate 10 random epidemics starting from our initial conditions:

```{r}
dust2::dust_system_set_state_initial(sys)
time <- 0:100
y <- dust2::dust_system_simulate(sys, time)
matplot(time, t(y[2, , ]), col = "#00000055", lty = 1, type = "l",
        xlab = "Time", ylab = "Number of infecteds")
```

as this system is stochastic, each trajectory will be different.