Introduction

The implementation of differential equation models directly in R (R Core Team 2018) is ideal for the development of dynamic models, because everything can be immediately tested and the whole infrastructure of R can be used. However, this adds considerable overhead, especially due to the automatic checking and type conversion of R, so that models cannot run with unlimited performance in this mode.

Fortunately, there is a way to unlock the handbreake. Both R in general and package deSolve (Soetaert, Petzoldt, and Setzer 2010a; Soetaert, Petzoldt, and Setzer 2010b) in particular allow to write time-critical parts in a so-called compiled language, e.g. Fortran or C/C++. Then, all solvers of deSolve allow direct communication between solvers and a compiled model on the machine code level, while keeping R’s facilities for data management and visualisation. Details are described in the package vignette \(\rightarrow\) compiledCode of deSolve while the short tutorial here serves as a quick introduction.

General principle

  • Implement core model (and only this) in C or Fortran,
  • Use data handling, storage and plotting facilities of R.

R starts the solver with the compiled model, both communicate directly without the overhead of calling back to R and finally the return the results tu the R level where it is accessible by the user.

The Lorenz system in pure R

library(deSolve)
library(scatterplot3d)

Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <- a * X + Y * Z
    dY <- b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c =  28)
state <- c(X = 1, Y = 1, Z = 1)
times <- seq(0, 100, by = 0.01)

out <- ode(y = state, times = times, func = Lorenz, parms = parameters)

plot(out)
scatterplot3d(out[,-1], type="l")

The Lorenz system in C

/* file lorenz.c */
#include <R.h>
static double parms[3];
#define a parms[0]
#define b parms[1]
#define c parms[2]

/* initializer  */
void initmod(void (* odeparms)(int *, double *)) {
    int N = 3;
    odeparms(&N, parms);
}

/* Derivatives */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip) {
    ydot[0] = a * y[0] + y[1] * y[2];
    ydot[1] = b * (y[1] - y[2]);
    ydot[2] = - y[0] * y[1] + c * y[1] - y[2];
}

Compile the model

  • needs compilers (RTools) installed and in the systems’ search path
system("R CMD SHLIB lorenzc.c")

Call the model from R

library(deSolve)
library(scatterplot3d)
dyn.load("lorenzc.dll")

out <- ode(state, times, func = "derivs", parms = parameters,
  dllname = "lorenzc", initfunc = "initmod")

dyn.unload("lorenzc.dll")

plot(out)
scatterplot3d(out[,-1], type="l")

Don’t forget to dyn.unload the model afterwards, otherwise the DLL will remain locked and cannot be overwritten in a next compiler run.

The Lorenz system in Fortran

c file lorenz.f
       subroutine initmod(odeparms)
         external odeparms
         double precision parms(3)
         common /myparms/parms

         call odeparms(3, parms)
         return
       end

       subroutine derivs (neq, t, y, ydot, yout, ip)
         double precision t, y, ydot, a, b, c
         integer neq, ip(*)
         dimension y(3), ydot(3), yout(*)
         common /myparms/a,b,c

         ydot(1) = a * y(1) + y(2) * y(3)
         ydot(2) = b * (y(2) - y(3))
         ydot(3) = -y(1) * y(2) + c * y(2) - y(3)

        return
       end

Compile the model

  • needs compilers (Rtools) installed and in the systems’ search path
  • this is quite standard on Linux, whereas Windows users have to download the Rtools separately.
system("R CMD SHLIB lorenzf.f")

Call the model from R

dyn.load("lorenzf.dll")

out <- ode(state, times, func = "derivs", parms = parameters,
  dllname = "lorenzf", initfunc = "initmod")

dyn.unload("lorenzf.dll")

plot(out)
scatterplot3d(out[,-1], type="l")

Don’t forget to unload the model afterwards, otherwise the DLL will remain locked and cannot be overwritten in a next compiler run.

Benchmark

Run the 3 versions (R, C, Fortran) several times in an R script and make a benchmark.

library(deSolve)
library(scatterplot3d)

## put the pure R model here

## compile the C code
system("R CMD SHLIB lorenzc.c")

## compile the fortran code
system(..............)

## load the DLLs (note .so on Linux)
dyn.load("lorenzc.dll")
dyn.load(............)

## run the R model 10 (or 100) times
system.time(
  for(i in 1:10) out   <- ode(state, times, Lorenz, parms = parameters)
)

## run the C model 10 (or 100) times
..............

## run the Fortran model 10 (or 100) times
..............

## don't forget to unload the DLLs (.so on Linux)
dyn.unload("lorenzc.dll")
dyn.unload("lorenzf.dll")

Outlook

More can be found in the package vignette compiledCode or in the paper of Soetaert, Petzoldt, and Setzer (2010b). In addition to this, compiled differential equation models are also supported by other packages like cOde and dMod (Kaschek 2018a; Kaschek 2018b) or rodeo (Kneis, Petzoldt, and Berendonk 2017; Kneis 2018).

References

Kaschek, Daniel. 2018a. cOde: Automated c Code Generation for ’deSolve’, ’bvpSolve’ and ’Sundials’. https://CRAN.R-project.org/package=cOde.

———. 2018b. dMod: Dynamic Modeling and Parameter Estimation in Ode Models. https://CRAN.R-project.org/package=dMod.

Kneis, David. 2018. Rodeo: A Code Generator for Ode-Based Models. https://CRAN.R-project.org/package=rodeo.

Kneis, David, Thomas Petzoldt, and Thomas U Berendonk. 2017. “An R-Package to Boost Fitness and Life Expectancy of Environmental Models.” Environmental Modelling & Software 96. Elsevier: 123–27. doi:10.1016/j.envsoft.2017.06.036.

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Soetaert, Karline, Thomas Petzoldt, and R. Woodrow Setzer. 2010a. “Solving Differential Equations in R.” The R Journal 2 (2): 5–15. https://journal.r-project.org/archive/2010-2/RJournal_2010-2_Soetaert~et~al.pdf.

———. 2010b. “Solving Differential Equations in R: Package deSolve.” Journal of Statistical Software 33 (9): 1–25. doi:10.18637/jss.v033.i09.

Copyright and original author: tpetzoldt, 2019-02-07