odin.qmd

# Getting started with odin {#sec-getting-started}

```{r}
#| include: false
source("common.R")
```

`odin` is a "[Domain Specific Language](https://en.wikipedia.org/wiki/Domain-specific_language)" (DSL), that is a mini-language that is specifically tailored to a domain.  By targeting a narrow domain we can have a language that expressively captures a problem and gives an efficient solution, as opposed to a general purpose language where you can write code that expresses *anything*, but your solution might either be quite verbose or have poor performance.  You might be familiar with other DSLs such as the `dplyr` or `ggplot2` syntax, or SQL if you have used databases.

```{r}
library(odin2)
```

## A simple example {#sec-odin-sir}

Here is a small system of differential equations for an "SIR" (Susceptible-Infected-Recovered) model:

\begin{gather*}
\frac{dS}{dt} = -\beta S \frac{I}{N}\\
\frac{dI}{dt} = \beta S \frac{I}{N} - \gamma I\\
\frac{dR}{dt} = \gamma I
\end{gather*}

And here is an implementation of these equations in `odin`:

```{r}
sir <- odin({
  deriv(S) <- -beta * S * I / N
  deriv(I) <- beta * S * I / N - gamma * I
  deriv(R) <- gamma * I

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

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

There are a few things to note here:

* equations are defined out of order; we just will things into existence and trust `odin` to put them in the right order for us (much more on this later)
* our system has "parameters"; these are things we'll be able to change easily in it after creation
* every variable (as in, model variable --- those that make up the state) has a pair of equations; an initial condition via `initial()` and an equation for the derivative with respect to time via `deriv()`

::: {.callout-note}
Above you will see us struggle a bit with terminology, particularly between "system" and "model".  There are going to be two places where "model" can be usefully used - a model could be the **generative model** (here an SIR model) or it could be the **statistical model**, more the focus on the second half of the book.  In order to make things slightly less ambiguous we will often refer to the generative model as a "system", and this is reflected in the calls to functions in dust.
:::

The `odin()` function will compile (more strictly, [transpile](https://en.wikipedia.org/wiki/Source-to-source_compiler)) the odin code to C++ and then compile this into machine code and load it into your session.  The resulting object is a `dust_sytem_generator` object:

```{r}
sir
```

As its name suggests, we can use this object to generate different versions of our system with different configurations. We pass this object to other functions from `dust2` to create, run and examine our simulations.

## Running systems using `dust2`

```{r}
library(dust2)
```

We create a system by using `dust_system_create` and passing in a list of parameters:

```{r}
pars <- list()
sys <- dust_system_create(sir, pars)
sys
```

To interact with a system, you will use functions from `dust2`, which we'll show below.

All systems start off with a state of all zeros:

```{r}
dust_system_state(sys)
```

There are two ways of setting state:

1. Provide a new state vector and set it with `dust_system_set_state`
2. Use the initial condition from the model itself (the expressions with `initial()` on the left hand side)

Here, we'll use the latter, and use the initial condition defined in the odin code:

```{r}
dust_system_set_state_initial(sys)
dust_system_state(sys)
```

In general, the order of the state is arbitrary (though in practice it is fairly stable).  To turn this state vector into something more interpretable you can use `dust_unpack_state`:

```{r}
s <- dust_system_state(sys)
dust_unpack_state(sys, s)
```

This gives a named list of numbers, which means we can work with the output of the model in a more reasonable way.

Next, we want to run the system.  There are two main functions for doing this:

1. `dust_system_run_to_time` which runs the system up to some point in time, but returns nothing
2. `dust_system_simulate` which runs the system over a series of points in time and returns the state at these points in time

Here we'll do the latter as it will produce something we can look at!

```{r}
t <- seq(0, 150, by = 0.25)
y <- dust_system_simulate(sys, t)
dim(y)
```

This produces a `r nrow(y)` x `r ncol(y)` matrix.  This is typical of how time-series data are produced from odin2/dust2 and may be a surprise:

* state will always be on the first dimension
* time will always be on the last dimension

This may take some getting used to at first, but practically some degree of manipulation will be required in any case if you want to plot things.

We'll explore some other ways of plotting later, but here's the number of individuals in the infectious class over time, as the epidemic proceeds.

```{r}
plot(t, dust_unpack_state(sys, y)$I, type = "l",
     xlab = "Time", ylab = "Infected population")
```

Once a system has been run to a time, you'll need to reset it to run it again.  For example, this won't work

```{r}
#| error: true
dust_system_simulate(sys, t)
```

The error here is trying to tell us that the first time in our vector `t`, which is 0, is smaller than the current time in the system, which is `r max(t)`.  We can query the system to get its current time, too:

```{r}
dust_system_time(sys)
```

You can reset the time back to zero (or any other time) using `dust_system_set_time`:

```{r}
dust_system_set_time(sys, 0)
```

this does not reset the state, however:

```{r}
dust_system_state(sys)
```

We can set new parameters, perhaps, and then update the state:

```{r}
dust_system_update_pars(sys, list(I0 = 5))
dust_system_set_state_initial(sys)
dust_system_state(sys)
```

## Going further

These are the lowest-level functions you would typically use, and we expect that you would combine these together yourself to perform specific tasks.  For example suppose you wanted to explore how `beta` affected the epidemic trajectory over this set of times, we might write:

```{r}
run_with_beta <- function(beta, t) {
  sys <- dust_system_create(sir, list(beta = beta))
  dust_system_set_state_initial(sys)
  idx <- dust_unpack_index(sys)$I
  drop(dust_system_simulate(sys, t, index_state = idx))
}
```

Here, we've used `dust_unpack_index` to get the index of the `I` compartment and passed that in as the index to `dust_system_simulate`.

We could then run this over a series of `beta` values:

```{r}
beta <- seq(0.1, 0.3, by = 0.02)
y <- vapply(beta, run_with_beta, t, FUN.VALUE = numeric(length(t)))
matplot(t, y, type = "l", lty = 1, col = "steelblue4",
        xlab = "Time", ylab = "Infected population")
```

We could modify this to take a single system object and reset it each time (which might be useful if our model was slow to initialise).  Alternatively we could simulate the whole grid of beta values at once:

```{r}
pars <- lapply(beta, function(b) list(beta = b))
sys <- dust_system_create(sir, pars, n_groups = length(pars))
dust_system_set_state_initial(sys)
idx <- dust_unpack_index(sys)$I
y <- dust_system_simulate(sys, t, index_state = idx)
matplot(t, t(y[1, , ]), type = "l", lty = 1, col = "steelblue4",
        xlab = "Time", ylab = "Infected population")
```

Or perhaps you would want to combine these approaches and start the simulation from a range of stochastic starting points.  Given the unbounded range of things people want to do with dynamical models, the hope is that you can use the functions in `dust2` in order to build workflows that match your needs, and that these functions should be well documented and efficient.