2  Getting started with odin

odin is a “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 expresively 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.

2.1 A simple example

Here is a simple “SIR” (Susceptible-Infected-Recovered) model, implemented as a set of ordinary differential equations:

deriv(S) <- -beta * S * I / N
deriv(I) <- beta * S * I / N - sigma * I
deriv(R) <- sigma * I

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

N <- parameter(1e6)
I0 <- parameter(1)
beta <- parameter(4)
sigma <- parameter(2)

There are a few things to note here:

  • equations are defined out of order; we just will things into existance 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()

2.2 Compiling models

Suppose we have this model stored in a file models/sir.R, we can compile this model with odin():

sir <- odin2::odin("models/sir.R")
#> ✔ Wrote 'DESCRIPTION'
#> ✔ Wrote 'NAMESPACE'
#> ✔ Wrote 'R/dust.R'
#> ✔ Wrote 'src/dust.cpp'
#> ✔ Wrote 'src/Makevars'
#> ℹ 13 functions decorated with [[cpp11::register]]
#> ✔ generated file 'cpp11.R'
#> ✔ generated file 'cpp11.cpp'
#> ℹ Re-compiling odina3f6ff4b
#> ── R CMD INSTALL ───────────────────────────────────────────────────────────────
#> * installing *source* package ‘odina3f6ff4b’ ...
#> ** using staged installation
#> ** libs
#> using C++ compiler: ‘g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> g++ -std=gnu++17 -I"/opt/R/4.4.1/lib/R/include" -DNDEBUG  -I'/home/runner/work/_temp/Library/cpp11/include' -I'/home/runner/work/_temp/Library/dust2/include' -I'/home/runner/work/_temp/Library/monty/include' -I/usr/local/include   -DHAVE_INLINE -fopenmp  -fpic  -g -O2  -Wall -pedantic  -c cpp11.cpp -o cpp11.o
#> g++ -std=gnu++17 -I"/opt/R/4.4.1/lib/R/include" -DNDEBUG  -I'/home/runner/work/_temp/Library/cpp11/include' -I'/home/runner/work/_temp/Library/dust2/include' -I'/home/runner/work/_temp/Library/monty/include' -I/usr/local/include   -DHAVE_INLINE -fopenmp  -fpic  -g -O2  -Wall -pedantic  -c dust.cpp -o dust.o
#> g++ -std=gnu++17 -shared -L/opt/R/4.4.1/lib/R/lib -L/usr/local/lib -o odina3f6ff4b.so cpp11.o dust.o -fopenmp -L/opt/R/4.4.1/lib/R/lib -lR
#> installing to /tmp/RtmpDpNQIR/devtools_install_18616fbaa42d/00LOCK-dust_cb2273f/00new/odina3f6ff4b/libs
#> ** checking absolute paths in shared objects and dynamic libraries
#> * DONE (odina3f6ff4b)
#> ℹ Loading odina3f6ff4b

The result of this is a function which creates a dust_system_generator object, which we can use with functions from dust2 to simulate from this model.

sir()
#> 
#> ── <dust_system_generator: odin> ───────────────────────────────────────────────
#> ℹ Use 'dust2::dust_system_create()' to create a system with this generator
#> ℹ This system runs in continuous time

2.3 Running systems using dust2

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

pars <- list()
sys <- dust2::dust_system_create(sir(), pars, n_particles = 1)
sys
#> 
#> ── <dust_system: odin> ─────────────────────────────────────────────────────────
#> ℹ single particle with 3 states
#> ℹ This system runs in continuous time

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:

dust2::dust_system_state(sys)
#> [1] 0 0 0

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

dust2::dust_system_set_state_initial(sys)
dust2::dust_system_state(sys)
#> [1] 999999      1      0

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:

s <- dust2::dust_system_state(sys)
dust2::dust_unpack_state(sys, s)
#> $S
#> [1] 999999
#> 
#> $I
#> [1] 1
#> 
#> $R
#> [1] 0

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!

t <- seq(0, 150, by = 0.25)
y <- dust2::dust_system_simulate(sys, t)
dim(y)
#> [1]   3 601

This produces a 3 x 601 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.

plot(t, dust2::dust_unpack_state(sys, y)$I, type = "l",
     xlab = "Time", ylab = "Total infectious")

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

dust2::dust_system_simulate(sys, t)
#> Error: Expected 'time[1]' (0) to be larger than the previous value (150)

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 150. We can query the system to get its current time, too:

dust2::dust_system_time(sys)
#> [1] 150

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

dust2::dust_system_set_time(sys, 0)

this does not reset the state, however:

dust2::dust_system_state(sys)
#> [1] 2.031875e+05 6.508638e-69 7.968125e+05

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

dust2::dust_system_update_pars(sys, list(I0 = 5))
dust2::dust_system_set_state_initial(sys)
dust2::dust_system_state(sys)
#> [1] 999995      5      0