9  Debugging

There is functionality in odin2 to help with debugging, including

library(odin2)
library(dust2)

9.1 Compilation error messages

Compilation error messages in odin2 are designed to help you identify why compilation of your model code with odin() has failed.

Let’s revisit the simple stochastic SIR model from Section 3.2:

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

We can see that this model successfully compiled. Now let’s suppose we omitted the initial(R) <- 0 expression:

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")`

Here the error message is telling us that we need an initial() call for R as we have an update() call for it. All error messages include a prompt to run odin2::odin_error_explain("Exxxx") for more information where Exxxx will be the error code for the specific type of error (the full list of error codes can be found here). In RStudio, you can click on the link in the console to show information about the error message; alternatively, running in this case gives:

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()`.
#>

9.2 Debugging tools

Revisiting the the simple stochastic SIR model from Section 3.2 again, let’s suppose that instead of

update(I) <- I + n_SI - n_IR

we write

update(I) <- I + n_SI

In this case we have a model that successfully compiles but it contains 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

If we were to create a dust system from this and 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] 8
#> 
#> $I
#> [1] 992
#> 
#> $R
#> [1] 2596
#> 
#> $incidence
#> [1] 1

then it is clear that there is a bug in our model as the number of individuals in the R compartment has exceeded N.

When you believe there is a bug in a successfully compiled model, you can use special functionality provided by odin2to help you debug:

  • 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.
  • You can examine the code generated from your model with odin_show().

9.2.1 Using print()

In the code below, we can print out the value of R as we go.

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("R: {R}")
})

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)
dust_system_run_to_time(sys, 10)
#> [0.000000] R: 0.000000
#> [1.000000] R: 0.000000
#> [2.000000] R: 1.000000
#> [3.000000] R: 1.000000
#> [4.000000] R: 2.000000
#> [5.000000] R: 3.000000
#> [6.000000] R: 7.000000
#> [7.000000] R: 10.000000
#> [8.000000] R: 13.000000
#> [9.000000] R: 18.000000

Here we see the value of R for all of our 10 time points. We can also print expressions e.g. S + I + R. We can also be more selective, choosing when we would like to print using the when argument to print() like this:

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 + R > N: {S + I + R}", when = S + I + R > N)
})

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)
dust_system_run_to_time(sys, 10)
#> [3.000000] S + I + R > N: 1001.000000
#> [4.000000] S + I + R > N: 1003.000000
#> [5.000000] S + I + R > N: 1005.000000
#> [6.000000] S + I + R > N: 1008.000000
#> [7.000000] S + I + R > N: 1016.000000
#> [8.000000] S + I + R > N: 1018.000000
#> [9.000000] S + I + R > N: 1021.000000

The when function can take more complicated logical expressions, including brackets, and || and && for or and and respectively. While at present we can only have one print() line within your model, you can print multiple variables together like this:

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

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: 12.000000
#> [7.000000] S + I: 1000.000000, R: 16.000000
#> [8.000000] S + I: 1000.000000, R: 19.000000
#> [9.000000] S + I: 1000.000000, R: 26.000000

9.2.1.1 Formatting

We use glue to drive the formatting, which may feel familiar to you already. Above, we surrounded variables with {curly braces} to print their value instead of their name, where all other text was printed verbatim.

You can control the precision of numbers being printed in the form below:

m <- odin({
  update(x) <- Normal(x, 1)
  initial(x) <- 1
  print("x: {x; .2f}")
})
sys <- dust_system_create(m)
dust_system_run_to_time(sys, 10)
#> [0.000000] x: 0.00
#> [1.000000] x: -0.33
#> [2.000000] x: 0.13
#> [3.000000] x: 0.86
#> [4.000000] x: 3.13
#> [5.000000] x: 1.48
#> [6.000000] x: 2.32
#> [7.000000] x: 2.62
#> [8.000000] x: 3.11
#> [9.000000] x: 2.78

Here, the .2f gives us two decimal places; you may recognise that format if you’ve used sprintf(). By default, variables are printed as floats (f), and you may sometimes find g useful for particularly small or large values.

Note that this is an experimental; see the debugging vignette for current limitations.

9.2.2 Using the interactive debugger, browser()

You can use the browser() function to ask odin to drop you into a debugger which will appear very similar to doing so in R, but has a couple of arguments that are specific to odin.

  • phase: the system phase where the debugger should be inserted; this will typically be update or deriv.
  • when: optionally a condition that should be satisfied for the debugger to be triggered. You will typically want to set this or it will be called at every step.

For example, for our SIR model with a bug we might write:

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

The location of the call to browser() does not matter; it will be activated at the end of the phase. The condition here might be something we cook up to look at what happens as the number of individuals in the infected class starts tailing off at the end of the simulation.

We create and initialise the system as normal:

sys <- dust_system_create(sir)
dust_system_set_state_initial(sys)

However, when you run the system you will pause part way through evaluation:

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]>

Now you can explore things that odin and dust know about, in a similar way that you would in R.

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

and you can inspect values or perform calculations:

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

You can change values for testing just within this browser environment, but note that the changes will not be injected back into your model.

If you press c or n, then the system will proceed to the next step and drop you back into the debugger, when it reaches the browser command again, just like in R. So if your model had a variable called n, you would have to type message(n), or use sprintf to see its value. You can exit with Q which will return you to the console with an error.

To cancel all future browser calls and run to the end of your simulation:

dust_browser_continue()

and then the next time you press c to continue, your simulation will run to completion.

9.2.3 Showing the generated code

Sometimes just looking at the generated code can be helpful. You can do this with odin_show:

odin_show({
  initial(x) <- 2
  update(x) <- Normal(x, 3)
})
#> 
#> ── odin code: ──────────────────────────────────────────────────────────────────
#> #include <dust2/common.hpp>
#> // [[dust2::class(odin_system)]]
#> // [[dust2::time_type(discrete)]]
#> class odin_system {
#> public:
#>   odin_system() = delete;
#>   using real_type = double;
#>   using rng_state_type = monty::random::generator<real_type>;
#>   struct shared_state {
#>     struct odin_internals_type {
#>       struct {
#>         dust2::packing state;
#>       } packing;
#>       struct {
#>         std::array<size_t, 1> state;
#>       } offset;
#>     } odin;
#>   };
#>   struct internal_state {};
#>   using data_type = dust2::no_data;
#>   static dust2::packing packing_state(const shared_state& shared) {
#>     return shared.odin.packing.state;
#>   }
#>   static shared_state build_shared(cpp11::list parameters) {
#>     shared_state::odin_internals_type odin;
#>     odin.packing.state = dust2::packing{
#>       {"x", {}}
#>     };
#>     odin.packing.state.copy_offset(odin.offset.state.begin());
#>     return shared_state{odin};
#>   }
#>   static internal_state build_internal(const shared_state& shared) {
#>     return internal_state{};
#>   }
#>   static void update_shared(cpp11::list parameters, shared_state& shared) {
#>   }
#>   static void update_internal(const shared_state& shared, internal_state& internal) {
#>   }
#>   static void initial(real_type time, const shared_state& shared, internal_state& internal, rng_state_type& rng_state, real_type* state) {
#>     state[0] = 2;
#>   }
#>   static void update(real_type time, real_type dt, const real_type* state, const shared_state& shared, internal_state& internal, rng_state_type& rng_state, real_type* state_next) {
#>     const auto x = state[0];
#>     state_next[0] = monty::random::normal<real_type>(rng_state, x, 3);
#>   }
#> };

Above may seem quite a bit of code (for only two lines of odin code!), but if you look for the initial function and the update function at the bottom you will see the code representing your two model lines. To simplify things a little, we can show just a particular method (usually update or rhs), with the what argument, for example:

odin_show({
  initial(x) <- 2
  update(x) <- Normal(x, 3)
}, what = "update")
#> 
#> ── odin code (update): ─────────────────────────────────────────────────────────
#> static void update(real_type time, real_type dt, const real_type* state, const shared_state& shared, internal_state& internal, rng_state_type& rng_state, real_type* state_next) {
#>   const auto x = state[0];
#>   state_next[0] = monty::random::normal<real_type>(rng_state, x, 3);
#> }