<- odin2::odin("models/sir.R")
sir #> ✔ 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
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
<- parameter(1e6)
N <- parameter(1)
I0 <- parameter(4)
beta <- parameter(2) sigma
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 viaderiv()
2.2 Compiling models
Suppose we have this model stored in a file models/sir.R
, we can compile this model with odin()
:
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:
<- list()
pars <- dust2::dust_system_create(sir(), pars, n_particles = 1)
sys
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:
::dust_system_state(sys)
dust2#> [1] 0 0 0
There are two ways of setting state:
- Provide a new state vector and set it with
dust_system_set_state
- Use the initial condition from the model itself (the expressions with
initial()
on the left hand side)
Here, we’ll use the latter
::dust_system_set_state_initial(sys)
dust2::dust_system_state(sys)
dust2#> [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
:
<- dust2::dust_system_state(sys)
s ::dust_unpack_state(sys, s)
dust2#> $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:
dust_system_run_to_time
which runs the system up to some point in time, but returns nothingdust_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!
<- seq(0, 150, by = 0.25)
t <- dust2::dust_system_simulate(sys, t)
y 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
::dust_system_simulate(sys, t)
dust2#> 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:
::dust_system_time(sys)
dust2#> [1] 150
You can reset the time back to zero (or any other time) using dust_system_set_time
:
::dust_system_set_time(sys, 0) dust2
this does not reset the state, however:
::dust_system_state(sys)
dust2#> [1] 2.031875e+05 6.508638e-69 7.968125e+05
We can set new parameters, perhaps, and then update the state:
::dust_system_update_pars(sys, list(I0 = 5))
dust2::dust_system_set_state_initial(sys)
dust2::dust_system_state(sys)
dust2#> [1] 999995 5 0