dde

Models implemented in R

Models implemented in R look very similar to deSolve. Here is the Lorenz attractor implemented for dde:

lorenz_dde <- function(t, y, p) {
  sigma <- p$sigma
  R     <- p$R
  b     <- p$b
  y1 <- y[[1L]]
  y2 <- y[[2L]]
  y3 <- y[[3L]]
  c(sigma * (y2 - y1),
    R * y1 - y2 - y1 * y3,
    -b * y3 + y1 * y2)
}

The p argument is the parameters and can be any R object. Here I’ll use a list to hold the standard Lorenz attractor parameters:

p <- list(sigma = 10.0,
          R     = 28.0,
          b     =  8.0 / 3.0)

tt <- seq(0, 100, length.out = 50001)
y0 <- c(1, 1, 1)
yy <- dde::dopri(y0, tt, lorenz_dde, p)

Here is the iconic attractor

par(mar=rep(.5, 4))
plot(yy[, c(2, 4)], type = "l", lwd = 0.5, col = "#00000066",
     axes = FALSE, xlab = "", ylab = "")

The approach above is almost identical to implementing this model using deSolve:

lorenz_ds <- function(t, y, p) {
  sigma <- 10.0
  R     <- 28.0
  b     <-  8.0 / 3.0
  y1 <- y[[1L]]
  y2 <- y[[2L]]
  y3 <- y[[3L]]
  list(c(sigma * (y2 - y1),
         R * y1 - y2 - y1 * y3,
         -b * y3 + y1 * y2))
}
yy_ds <- deSolve::ode(y0, tt, lorenz_ds, p)

Dense output and interpolation

One of the nice things about the dopri solvers is that they do not need to stop the integration at the times that you request output at:

yy <- dde::dopri(y0, tt, lorenz_dde, p, return_statistics = TRUE)
attr(yy, "statistics")

##   n_eval   n_step n_accept n_reject 
##    26990     4498     4308      190

Above, the number of function evaluations (~6 per step), steps, and rejected steps is indicated (a rejected step occurs where the solver has to reduce step size multiple times to achieve the required accuracy). The number of steps here is about 1/10 the number of returned samples. This works because the solver here returns “dense output” which allows it to interpolate the solution between points that it has not visited. This is supported by many of the solvers in deSolve, too.

In contrast with deSolve, the dense output here can be collected and worked with later, though doing this requires a bit of faff.

Specify the history length; this needs to be an overestimate because once the end of the history buffer is reached it will be silently overwritten to return the last steps in history. (This is the behaviour required to support delay models without running out of memory).

yy2 <- dde::dopri(y0, range(tt), lorenz_dde, p, return_minimal = TRUE,
                  n_history = 5000, return_history = TRUE)

With these arguments yy2 is a 3 x 1 matrix, but it comes with a massive “history” matrix”:

dim(yy2)

## [1] 3 1

h <- attr(yy2, "history")
dim(h)

## [1]   17 4308

The contents of this matrix are designed to be opaque (i.e., I may change how things are represented at a future time). However, the solution can be interpolated to any number of points using this matrix:

yy2 <- dde::dopri_interpolate(h, tt)
all.equal(yy2, yy[, -1], check.attributes = FALSE)

## [1] TRUE

A model in R

This is a simple SEIR (Susceptible - Exposed - Infected - Resistant) model from epidemiology. Once exposed to the disease, an individual exists in an “Exposed” state for exactly 14 days before becoming “Infected” (you could model this with a series of compartments and get a distribution of exposed times).

seir <- function(t, y, p) {
  b <- 0.1
  N <- 1e7
  beta <- 10.0
  sigma <- 1.0 / 3.0
  delta <- 1.0 / 21.0
  t_latent <- 14.0
  Births <- N * b
  surv <- exp(-b * t_latent)

  S <- y[[1L]]
  E <- y[[2L]]
  I <- y[[3L]]
  R <- y[[4L]]

  tau <- t - t_latent
  y_lag <- dde::ylag(tau, c(1L, 3L)) # Here is ylag!
  S_lag <- y_lag[[1L]]
  I_lag <- y_lag[[2L]]

  new_inf <- beta * S * I / N
  lag_inf <- beta * S_lag * I_lag * surv / N

  c(Births  - b * S - new_inf + delta * R,
    new_inf - lag_inf - b * E,
    lag_inf - (b + sigma) * I,
    sigma * I - b * R - delta * R)
}

The model needs to know how many susceptible individuals there were 14 days ago, and how many infected there were 14 days ago. To get this from the model, we use

y_lag <- dde::ylag(tau, c(1L, 3L))

to get the values of the first and third variables (S and I) at time tau. Alternatively you can get all values with

y_lag <- dde::ylag(tau)

or get them individually

S_lag <- dde::ylag(tau, 1L)
I_lag <- dde::ylag(tau, 3L)

The ylag function can only be called from within an integration; it will throw an error if you try to call it otherwise.

What happens when we start though? If time starts at 0, then the first tau is -14 and we have no history then. dde keeps track of the initial state of the system and if a time before this is requested it returns the initial state of a variable. This is going to be reasonable for many applications but will lead to discontinuities in the derivative of your solution (and the second derivative and so on). This can make the problem hard to solve, and it may be preferable to provide your own information (see the deSolve implementation below for one possible way of implementing this).

To integrate the problem, use the dde::dopri function (which by default will use the 5th order method, which is probably the best bet for most problems). You need to provide arguments:

• n_history: number of history elements to retain. If this is too low then the integration will stop with an error and you can increase it
• return_history: set this to FALSE if you won’t want the history matrix returned; returning it costs a little time and if you don’t want to inspect it it’s better to leave it off

y0 <-  y0 <- c(1e7 - 1, 0, 1, 0)
tt <- seq(0, 365, length.out = 100)
yy <- dde::dopri(y0, tt, seir, NULL, n_history = 1000L, return_history = FALSE)
matplot(tt, yy[, 2:5], type="l")

Models implemented in R

Above, I implemented a derivative function for an SEIR model for dde as

function (t, y, p) 
{
    b <- 0.1
    N <- 1e+07
    beta <- 10
    sigma <- 1/3
    delta <- 1/21
    t_latent <- 14
    Births <- N * b
    surv <- exp(-b * t_latent)
    S <- y[[1L]]
    E <- y[[2L]]
    I <- y[[3L]]
    R <- y[[4L]]
    tau <- t - t_latent
    y_lag <- dde::ylag(tau, c(1L, 3L))
    S_lag <- y_lag[[1L]]
    I_lag <- y_lag[[2L]]
    new_inf <- beta * S * I/N
    lag_inf <- beta * S_lag * I_lag * surv/N
    c(Births - b * S - new_inf + delta * R, new_inf - lag_inf - 
        b * E, lag_inf - (b + sigma) * I, sigma * I - b * R - 
        delta * R)
}
<bytecode: 0x55bf254547d0>

The implementation using deSolve looks very similar:

seir_deSolve <- function(t, y, parms) {
  b <- 0.1
  N <- 1e7
  beta <- 10
  sigma <- 1 / 3
  delta <- 1 / 21
  t_latent <- 14.0
  I0 <- 1
  Births <- N * b
  surv <- exp(-b * t_latent)

  S <- y[[1L]]
  E <- y[[2L]]
  I <- y[[3L]]
  R <- y[[4L]]

  tau <- t - t_latent
  if (tau < 0.0) { # NOTE: assuming that t0 is always zero
    S_lag <- parms$S0
    I_lag <- parms$I0
  } else {
    y_lag <- deSolve::lagvalue(tau, c(1L, 3L))
    S_lag <- y_lag[[1L]]
    I_lag <- y_lag[[2L]]
  }

  new_inf <- beta * S * I / N
  lag_inf <- beta * S_lag * I_lag * surv / N

  list(c(Births  - b * S - new_inf + delta * R,
         new_inf - lag_inf - b * E,
         lag_inf - (b + sigma) * I,
         sigma * I - b * R - delta * R))
}

The differences are that:

• deSolve requires that the derivatives are returned as a list, whereas dde uses a numeric vector (see below for details about this)
• deSolve requires that you provide the initial values for the lagged values (and we also need to know what the initial time is too, but I’m assuming that as zero)
• The appropriate function for pulling previous values from the history buffer is deSolve::lagvalue (for dde it is dde::ylag)

Aside from this the code is essentially identical.

To run the model with deSolve, use deSolve::dede which automatically sets up a history buffer of 10000 elements (the mxhist element of the control list alters this).

y0 <-  y0 <- c(1e7 - 1, 0, 1, 0)
tt <- seq(0, 365, length.out = 100)
initial <- list(S0 = y0[[1]], I0 = y0[[3]])
yy_ds <- deSolve::dede(y0, tt, seir_deSolve, initial)

This produces output that the same as dde:

yy_dde <- dde::dopri(y0, tt, seir, NULL, n_history = 1000L,
                     return_history = FALSE)
op <- par(mfrow=c(1, 2), mar=c(4, .5, 1.4, .5), oma=c(0, 2, 0, 0))
matplot(tt, yy_dde[, -1], type="l", main = "dde")
matplot(tt, yy_ds[, -1], type="l", main = "deSolve", yaxt="n")

The performance of both packages is fairly similar, taking a few tens of milliseconds to run on my machines

tR <- microbenchmark::microbenchmark(times = 30,
  deSolve = deSolve::dede(y0, tt, seir_deSolve, initial),
  dde = dde::dopri(y0, tt, seir, NULL, n_history = 1000L,
                   return_history = FALSE))
tR

## Unit: milliseconds
##     expr       min        lq      mean   median        uq       max neval
##  deSolve 16.164625 16.392670 17.315180 17.25864 17.593618 22.478953    30
##      dde  6.806705  6.872548  7.597816  7.08565  7.927574  9.941387    30

Models implemented in C

The compiled code interface for deSolve has greatly influenced dde and models implemented in either framework will be similar. Eventually dde may support a fully deSolve compatible interface but for now there are a few differences.

#include <R.h>
#include <R_ext/Rdynload.h>

void lagvalue(double tau, int *nr, int N, double *ytau);

// The parameters are going to be arranged:
//
//   t0
//   S0, I0
//   (b, N, beta, sigma, delta, t_latent)
//
// See below for why t0, S0 and I0 are stored
static double parms[3];

// The standard deSolve initialisation function
void seir_initmod(void (* odeparms)(int *, double *)) {
  int N = 3;
  odeparms(&N, parms);
}

// The RHS
void seir_deSolve(int *n, double *t, double *y, double *dydt,
                  double *yout, int *ip) {
  // again, hard-coded parameters for now; will change this shortly
  // once I get the same working with the dde impementation.
  double b = 0.1, N = 1e7, beta = 10.0, sigma = 1.0 / 3.0,
    delta = 1.0 / 21.0, t_latent = 14.0;
  double Births = N * b, surv = exp(-b * t_latent);

  // Because of the way that deSolve implements delays we need to
  // store the initial time and values in the parameters vector; if
  // the requested time is earlier than the time we started at then
  // the initial values need to be used, which we also store in the
  // parameters.
  double t0 = parms[0];
  const double tau = *t - t_latent;
  static int idx[2] = {0, 2};
  double S_lag, I_lag;
  if (tau <= t0) {
    S_lag = parms[1];
    I_lag = parms[2];
  } else {
    double ylag[2];
    lagvalue(tau, idx, 2, ylag);
    S_lag = ylag[0];
    I_lag = ylag[1];
  }

  const double S = y[0], E = y[1], I = y[2], R = y[3];
  const double new_inf = beta * S * I / N;
  const double lag_inf = beta * S_lag * I_lag * surv / N;

  dydt[0] = Births  - b * S - new_inf + delta * R;
  dydt[1] = new_inf - lag_inf - b * E;
  dydt[2] = lag_inf - (b + sigma) * I;
  dydt[3] = sigma * I - b * R - delta * R;
}

// This is the interface to deSolve's lag functions.  Note that unlike
// dde you are responsible for checking for underflows and providing
// values for underflowed times.
void lagvalue(double tau, int *nr, int N, double *ytau) {
  typedef void lagvalue_t(double, int *, int, double *);
  static lagvalue_t *fun = NULL;
  if (fun == NULL) {
    fun = (lagvalue_t*) R_GetCCallable("deSolve", "lagvalue");
  }
  fun(tau, nr, N, ytau);
}

This looks very similar to the dde version above but:

• parms (or whatever parameters are called) are handled as a global variable that is updated via a model initialisation function, whereas in dde they’re passed in as a void pointer
• We need to keep track of the initial state of the system via passing in t0 and initial conditions for S and I * There is an argument double *yout for additional output variables (of length *ip; in dde these are handled via a separate function.
• The lagvalue function must be explicitly defined (which requires loading some R-related headers (in dde this is achieved by including <dde/dde.h> and <dde/dde.c>.

Apart from these details, the model definition should appear very similar.

initial <- c(0.0, y0[[1]], y0[[3]])

zz_ds <- deSolve::dede(y0, tt, "seir_deSolve", initial,
                       initfunc = "seir_initmod", dllname = "dde_seir_ds")
zz_dde <- dde::dopri(y0, tt, "seir", numeric(), dllname = "dde_seir",
                     n_history = 1000L, return_history = FALSE)

Check that outputs of these models are the same as the R version above:

all.equal(zz_ds, yy_ds, check.attributes = FALSE)

## [1] TRUE

all.equal(zz_dde, yy_dde, check.attributes = FALSE)

## [1] TRUE

Here, the timings are even closer and have dropped from on the order of 20 milliseconds to 0.5 milliseconds; so we’re getting a ~40x speed up from using compiled code.

tC <- microbenchmark::microbenchmark(
  deSolve = deSolve::dede(y0, tt, "seir_deSolve", initial,
                          initfunc = "seir_initmod", dllname = "dde_seir_ds"),
  dde = dde::dopri(y0, tt, "seir", numeric(), dllname = "dde_seir",
                   n_history = 1000L, return_history = FALSE))
tC

## Unit: microseconds
##     expr     min       lq     mean  median      uq      max neval
##  deSolve 460.518 478.5015 526.5187 500.508 518.161 2922.566   100
##      dde 939.231 950.6820 970.8470 961.642 982.481 1065.436   100

The difference in speed will tend to increase as the models become larger (in terms of numbers of equations and parameters). On the other hand, constructing large models in C can be a hassle (but see odin for a possible solution).

You can extract a little more performance by tweaking options to dde::dopri; in particular, adding return_minimal = TRUE will avoid transposing the output, binding the times on, and (if given) avoiding binding output variables. These costs may be nontrivial with bigger models, though the cost of running a larger model will likely be larger still. Previous version of R suffered from a large cost of looking up the address of the compiled function (Windows may still take longer to do this than macOS/Linux). In that case, use getNativeSymbolInfo("seir") and pass that through to dopri as the func argument.

ptr <- getNativeSymbolInfo("seir")
tC2 <- microbenchmark::microbenchmark(
  deSolve = deSolve::dede(y0, tt, "seir_deSolve", initial,
                          initfunc = "seir_initmod", dllname = "dde_seir_ds"),
  dde = dde::dopri(y0, tt, "seir", numeric(), dllname = "dde_seir",
                   n_history = 1000L, return_history = FALSE),
  dde2 = dde::dopri(y0, tt, ptr, numeric(), n_history = 1000L,
                    return_history = FALSE, return_minimal = TRUE))
tC2

## Unit: microseconds
##     expr     min       lq     mean   median       uq      max neval
##  deSolve 459.066 484.5175 518.2269 501.5400 515.4155 2039.531   100
##      dde 939.761 947.7815 969.2676 958.9820 974.3955 1051.289   100
##     dde2 913.412 920.9315 934.8000 928.0845 934.8820 1011.285   100