Building a model with odin

odin

• A “domain specific language”
• Originally designed for ordinary differential equations
• Subsequently extended to difference equations (enabling discrete-time stochastic models)
• Refined into second version odin2

Odin is a DSL - a domain specific language. It exists to describe a specific problem efficiently. We developed it to describe ordinary differential equations which evolve a system in continuous time. This is a great target for DSL because it’s just a collection of mathematical truths - rates exist independent of any idea of order of operation.

odin

library(odin2)
library(dust2)
library(monty)

ODE models

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)

\[\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*}\]

Things to note:

• out of order definition
• every variable has initial and deriv pair
• parameters declared with parameter(), default values can be given in brackets though not required

Compiling the model with odin2

sir_ode <- 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)
})

You can also compile code from a file by its path e.g. odin("models/sir-basic.R").

The odin call generates C code, compiles it with gcc or clang to create a shared library, and loads that plus support code into R to create a generator object gen

Running the model with dust2

sys <- dust_system_create(sir_ode, pars = list())
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)

The output has dimensions number of states x number of timepoints

dim(y)
#> [1]   3 101

Unpacking states

Output can be nicely unpacked into the different states using dust_unpack_state

• Code
• Plot

y <- dust_unpack_state(sys, y)

plot(t, y$I, type = "l", xlab = "Time", ylab = "I")

From the generator object gen we can construct a model (mod), here using no parameters – just using default parameters as defined above. We run over a set of times and output the value of the system at a number of intermediate times.

Setting the state

Upon creation all systems start with a state full of zeros.

Typically we need to set the state. We used

dust_system_set_state_initial(sys)

to set the state according to the initial conditions defined in the odin code. We can also set the state using a new state vector e.g.

dust_system_set_state(sys, c(990, 10, 0))

Setting the state

The order of states in a system can be found via

dust_unpack_index(sys)
#> $S
#> [1] 1
#> 
#> $I
#> [1] 2
#> 
#> $R
#> [1] 3

or

sys$packer_state$names()
#> [1] "S" "I" "R"

Inputting parameters (pars)

• We created our system with pars = list() - an empty list
• This means we use all default values for parameters in the odin code

• We can pass in a named list of parameters to use non-default values
• e.g. pars = list(beta = 0.35, gamma = 0.15) would use those values for beta and gamma and the default values for all other parameters
• You will need to include values for any parameters without default values!

Stochastic models

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

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

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

\[\begin{gather*} S(t + \Delta t) = S(t) - n_{SI}\\ I(t + \Delta t) = I(t) + n_{SI} - n_{IR}\\ R(t + \Delta t) = R(t) + n_{IR} \end{gather*}\]

• dt is a special parameter
• every variable has initial and update pair

…compared with ODE models

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

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

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

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)

Compiling with odin2

sir <- odin({
  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
  
  p_SI <- 1 - exp(-beta * I / N * dt)
  p_IR <- 1 - exp(-gamma * dt)
  n_SI <- Binomial(S, p_SI)
  n_IR <- Binomial(I, p_IR)
  
  N <- parameter(1000)
  I0 <- parameter(10)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
})

Functions in odin2

• we have used the exponential function exp()
• many mathematical functions are supported in odin2 - a full list can be found here

• we have used Binomial(S, p_SI) for a Binomial random draw with size = S and prob = p_SI
• many distributions are supported - a full list with their names and parameterisations can be found here

time in odin

• time is a special (continuous) variable in odin models
• In discrete-time (stochastic) models, time is updated in increments of size dt

Running a single simulation

• Code
• Plot

sys <- dust_system_create(sir, pars = list(), dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)

plot(t, y$I, type = "l", xlab = "Time", ylab = "Infected population")

Running multiple simulations

• Code
• Plot

sys <- dust_system_create(sir, pars = list(), n_particles = 50,
                                 dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)
matplot(t, t(y$I), type = "l", lty = 1, col = "#00000044",
        xlab = "Time", ylab = "Infected population")

Calculating incidence with zero_every

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI
  
  p_SI <- 1 - exp(-beta * I / N * dt)
  p_IR <- 1 - exp(-gamma * dt)
  n_SI <- Binomial(S, p_SI)
  n_IR <- Binomial(I, p_IR)
  
  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0
  
  N <- parameter(1000)
  I0 <- parameter(10)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
})

Incidence accumulates then resets

• Code
• Plot

sys <- dust_system_create(sir, pars = list(), dt = 1 / 128)
dust_system_set_state_initial(sys)
t <- seq(0, 20, by = 1 / 128)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)

plot(t[t %% 1 == 0], y$incidence[t %% 1 == 0], type = "o", pch = 19,
     ylab = "Infection incidence", xlab = "Time")
lines(t, y$incidence, col = "red")

Time-varying inputs: using time

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI
  
  seed <- if (time == seed_time) seed_size else 0
  
  p_SI <- 1 - exp(-beta * I / N * dt)
  p_IR <- 1 - exp(-gamma * dt)
  n_SI <- min(seed + Binomial(S, p_SI), S)
  n_IR <- Binomial(I, p_IR)
  
  initial(S) <- N
  initial(I) <- 0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0
  
  N <- parameter(1000)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  seed_time <- parameter()
  seed_size <- parameter()
})

Time-varying inputs: using time

• Code
• Plot

pars <- list(seed_time = 10, seed_size = 15)
sys <- dust_system_create(sir, pars = pars, dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)

plot(t, y$I, type = "l", xlab = "Time", ylab = "Infected population")

Time-varying inputs: using interpolate

The interpolate function in odin can be used for time-varying parameters, with specification of

• the times of changepoints
• the values at those changepoints
• the type of interpolation: linear, constant or spline

Time-varying inputs: using interpolate

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI
  
  beta <- interpolate(beta_time, beta_value, "constant")
  beta_time <- parameter()
  beta_value <- parameter()
  dim(beta_time, beta_value) <- parameter(rank = 1)
  
  p_SI <- 1 - exp(-beta * I / N * dt)
  p_IR <- 1 - exp(-gamma * dt)
  n_SI <- Binomial(S, p_SI)
  n_IR <- Binomial(I, p_IR)
  
  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0
  
  N <- parameter(1000)
  I0 <- parameter(10)
  gamma <- parameter(0.1)
})

Time-varying inputs: using interpolate

• Code
• Plot

pars <- list(beta_time = c(0, 30), beta_value = c(0.2, 1))
sys <- dust_system_create(sir, pars = pars, dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)

plot(t, y$I, type = "l", xlab = "Time", ylab = "Infected population")
abline(v = 30, lty = 3)

Using arrays

sir_age <- odin({
  # Equations for transitions between compartments by age group
  update(S[]) <- S[i] - n_SI[i]
  update(I[]) <- I[i] + n_SI[i] - n_IR[i]
  update(R[]) <- R[i] + n_IR[i]
  update(incidence) <- incidence + sum(n_SI)
  
  # Individual probabilities of transition:
  p_SI[] <- 1 - exp(-lambda[i] * dt) # S to I
  p_IR <- 1 - exp(-gamma * dt) # I to R
  
  # Calculate force of infection
  
  # age-structured contact matrix: m[i, j] is mean number
  # of contacts an individual in group i has with an
  # individual in group j per time unit
  m <- parameter()
  
  # here s_ij[i, j] gives the mean number of contacts an
  # individual in group i will have with the currently
  # infectious individuals of group j
  s_ij[, ] <- m[i, j] * I[j]
  
  # lambda[i] is the total force of infection on an
  # individual in group i
  lambda[] <- beta * sum(s_ij[i, ])
  
  # Draws from binomial distributions for numbers
  # changing between compartments:
  n_SI[] <- Binomial(S[i], p_SI[i])
  n_IR[] <- Binomial(I[i], p_IR)
  
  initial(S[]) <- S0[i]
  initial(I[]) <- I0[i]
  initial(R[]) <- 0
  initial(incidence, zero_every = 1) <- 0
  
  S0 <- parameter()
  I0 <- parameter()
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  
  n_age <- parameter()
  dim(S, S0, n_SI, p_SI) <- n_age
  dim(I, I0, n_IR) <- n_age
  dim(R) <- n_age
  dim(m, s_ij) <- c(n_age, n_age)
  dim(lambda) <- n_age
})

Using arrays

pars <- list(S0 = c(990, 1000),
             I0 = c(10, 0),
             m = matrix(c(1.8, 0.4, 0.4, 1.2) / 2000, 2, 2),
             beta = 0.2,
             gamma = 0.1,
             n_age = 2)

Here m[i, j] is the mean number of contacts per time unit for an individual in group i has with an individual in group j

Using arrays

• Code
• Plot

sys <- dust_system_create(sir_age, pars = pars, dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)

matplot(t, t(y$I), type = "l", lty = 1, col = c("red", "blue"),
        xlab = "Time", ylab = "Infected population")
legend("topright", c("children", "adults"), col = c("red", "blue"), lty = 1)

Key features of arrays

• All arrays (whether state variable or parameter) need a dim equation
• No use of i, j etc on LHS - indexing on the LHS is implicit
• Support for up to 8 dimensions, with index variables i, j, k, l, i5, i6, i7, i8
• Functions for reducing arrays such as sum, prod, min, max - can be applied over entire array or slices

Arrays: model with age and vaccine

sir_age_vax <- odin({
  # Equations for transitions between compartments by age group
  update(S[, ]) <- new_S[i, j]
  update(I[, ]) <- I[i, j] + n_SI[i, j] - n_IR[i, j]
  update(R[, ]) <- R[i, j] + n_IR[i, j]
  update(incidence) <- incidence + sum(n_SI)
  
  # Individual probabilities of transition:
  p_SI[, ] <- 1 - exp(-rel_susceptibility[j] * lambda[i] * dt) # S to I
  p_IR <- 1 - exp(-gamma * dt) # I to R
  p_vax[, ] <- 1 - exp(-eta[i, j] * dt)
  
  # Force of infection
  m <- parameter() # age-structured contact matrix
  s_ij[, ] <- m[i, j] * sum(I[j, ])
  lambda[] <- beta * sum(s_ij[i, ])
  
  # Draws from binomial distributions for numbers changing between
  # compartments:
  n_SI[, ] <- Binomial(S[i, j], p_SI[i, j])
  n_IR[, ] <- Binomial(I[i, j], p_IR)
  
  # Nested binomial draw for vaccination in S
  # Assume you cannot move vaccine class and get infected in same step
  n_S_vax[, ] <- Binomial(S[i, j] - n_SI[i, j], p_vax[i, j])
  new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, n_vax]
  new_S[, 2:n_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1]
  
  initial(S[, ]) <- S0[i, j]
  initial(I[, ]) <- I0[i, j]
  initial(R[, ]) <- 0
  initial(incidence, zero_every = 1) <- 0
  
  # User defined parameters - default in parentheses:
  S0 <- parameter()
  I0 <- parameter()
  beta <- parameter(0.0165)
  gamma <- parameter(0.1)
  eta <- parameter()
  rel_susceptibility <- parameter()
  
  # Dimensions of arrays
  n_age <- parameter()
  n_vax <- parameter()
  dim(S, S0, n_SI, p_SI) <- c(n_age, n_vax)
  dim(I, I0, n_IR) <- c(n_age, n_vax)
  dim(R) <- c(n_age, n_vax)
  dim(m, s_ij) <- c(n_age, n_age)
  dim(lambda) <- n_age
  dim(eta) <- c(n_age, n_vax)
  dim(rel_susceptibility) <- c(n_vax)
  dim(p_vax, n_S_vax, new_S) <- c(n_age, n_vax)
})

Arrays: boundary conditions

Multiline equations can be used to deal with boundary conditions, e.g. we have

new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, n_vax]
new_S[, 2:n_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1]

which we could also write as

new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j]
new_S[, 1] <- new_S[i, j] + n_S_vax[i, n_vax]
new_S[, 2:n_vax] <- new_S[i, j] + n_S_vax[i, j - 1]

or another way of writing this would be to use if ... else

new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] +
    (if (j == 1) n_S_vax[i, n_vax] else n_S_vax[i, j - 1])

Dimensions

Dimensions can be hard-coded:

dim(S) <- 2

or can be soft-coded via parameters:

dim(S) <- n_age
n_age <- parameter()

Dimensions

Array parameters can have dimensions explicitly declared

m <- parameter()
dim(m) <- c(n_age, n_age)
n_age <- parameter()

or detected when the parameters are passed into the system - but a dim equation is needed where you must explicitly state the rank:

dim(m) <- parameter(rank = 2)

Dimensions

Dimensions of one array can be declared as the same as the dimensions of another:

dim(m) <- c(n_age, n_age)
dim(s_ij) <- dim(m)

or dimensions of arrays that have the same dimensions can be declared collectively:

dim(m, s_ij) <- c(n_age, n_age)

Compilation error messages

Let’s revisit the stochastic SIR model, which compiles successfully:

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI

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

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0

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

sir
#> 
#> ── <dust_system_generator: odin_system> ────────────────────────────────────────
#> ℹ This system runs in discrete time with a default dt of 1
#> ℹ This system has 4 parameters
#> → 'N', 'I0', 'beta', and 'gamma'
#> ℹ Use dust2::dust_system_create() (`?dust2::dust_system_create()`) to create a system with this generator
#> ℹ Use coef() (`?stats::coef()`) to get more information on parameters

Compilation error messages

Now let’s remove the initial(R) equation

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI

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

  initial(S) <- N - I0
  initial(I) <- I0
  initial(incidence, zero_every = 1) <- 0

  N <- parameter(1000)
  I0 <- parameter(10)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
})
#> Error in `odin()`:
#> ! Variables used in 'update()' do not have 'initial()' calls
#> ✖ Did not find 'initial()' calls for 'R'
#> → Context:
#> update(R) <- R + n_IR
#> ℹ For more information, run `odin2::odin_error_explain("E2003")`

Compilation error messages

We can run odin_error_explain() with the corresponding error code as suggested (clickable in RStudio) for more information

odin_error_explain("E2003")
#> 
#> ── E2003 ───────────────────────────────────────────────────────────────────────
#> Variables are missing initial conditions.
#> 
#> All variables used in `deriv()` or `update()` require a corresponding entry in
#> `initial()` to set their initial conditions.  The error will highlight all
#> lines that have a `deriv()` or `update()` call that lacks a call to
#> `initial()`.  This can sometimes be because you have a spelling error in your
#> call to `initial()`.
#> 

Debugging

Here’s a model that compiles but has a bug

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI

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

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0

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

sir
#> ── <dust_system_generator: odin_system> ────────────────────────────────────────
#> ℹ This system runs in discrete time with a default dt of 1
#> ℹ This system has 4 parameters
#> → 'N', 'I0', 'beta', and 'gamma'
#> ℹ Use dust2::dust_system_create() (`?dust2::dust_system_create()`) to create a system with this generator
#> ℹ Use coef() (`?stats::coef()`) to get more information on parameters

Debugging

Let’s run it

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)
t <- 50
res <- dust_system_simulate(sys, t)
dust_unpack_state(sys, res)
#> $S
#> [1] 7
#> 
#> $I
#> [1] 993
#> 
#> $R
#> [1] 2469
#> 
#> $incidence
#> [1] 2

The number of individuals in the R compartment has exceeded N!

Debugging

You can stare at the code and hope to spot what’s wrong, but that can be like trying to find a needle in a haystack, especially for larger models.

So there are tools to help you:

• You can print() the values of variables in your model as it runs
• You can use the browser() function in the same way as you would in R.

Debugging: using print()

With print() we can tell it what to print and (optionally) conditions for when to print.

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI

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

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0

  N <- parameter(1000)
  I0 <- parameter(10)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  
  print("S + I: {S + I}, R: {R}", when = S + I + R > N && time > 5)
})

Debugging: using print()

Then when we run it prints what we wanted with the times in square brackets:

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)
dust_system_run_to_time(sys, 10)
#> [6.000000] S + I: 1000.000000, R: 15.000000
#> [7.000000] S + I: 1000.000000, R: 17.000000
#> [8.000000] S + I: 1000.000000, R: 21.000000
#> [9.000000] S + I: 1000.000000, R: 26.000000

This may help us spot that we’re not removing individuals from I when they move to R!

Debugging: using browser()

With the browser() we need to declare the phase for inserting (deriv for ODE models, update for stochastic models) and optionally when to trigger it.

sir <- odin({
  update(S) <- S - n_SI
  update(I) <- I + n_SI
  update(R) <- R + n_IR
  update(incidence) <- incidence + n_SI

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

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0
  initial(incidence, zero_every = 1) <- 0

  N <- parameter(1000)
  I0 <- parameter(10)
  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  
  browser(phase = "update", when = S + I + R > N)
})

Debugging: using browser()

When we try to run it will pause during evaluation:

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)
dust_system_run_to_time(sys, 50)
ℹ dust debug ('update'; time = 2): see dust_browser for help
Called from: eval(substitute(expr), data, enclos = parent.frame())
Browse[1]> 

You can list all the variables

Browse[1]> ls()
 [1] "beta"      "gamma"     "I"         "I0"        "incidence" "N"         "n_IR"      "n_SI"      "p_IR"     
[10] "p_SI"      "R"         "S"         "time"  

Debugging: using browser()

You can inspect values and perform calculations

Browse[1]> N
[1] 1000
Browse[1]> I
[1] 18
Browse[1]> S
[1] 982
Browse[1]> I / N
[1] 0.018

Debugging: using browser()

• Enter c or n, to continue until the next time the browser() is triggered
  
• If your model has a variable called n, you would need to type message(n) to see its value
• Enter Q to immediately exit to the console with an error.
• Run dust_browser_continue() to cancel all future browser() calls and continue running the model

Next steps

• Get started fitting odin models with monty
• Learn how to use delay() for delay differential equations
• Understand the order of events for discrete-time stochastic models
• For larger odin models, learn about packaging them