\documentclass[article,nojss]{jss} \DeclareGraphicsExtensions{.pdf,.eps} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts \usepackage{amsmath} \usepackage{xspace} \usepackage{verbatim} \usepackage[english]{babel} %\usepackage{mathptmx} %\usepackage{helvet} \usepackage[T1]{fontenc} \usepackage[latin1]{inputenc} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands. \newcommand{\Rmodels}{\textbf{\textsf{R models}}\xspace} \newcommand{\DLLmodels}{\textbf{\textsf{DLL models}}\xspace} \title{\proglang{R} Package \pkg{deSolve}, Writing Code in Compiled Languages} \Plaintitle{R Package deSolve, Writing Code in Compiled Languages} \Keywords{differential equation solvers, compiled code, performance, \proglang{FORTRAN}, \proglang{C}} \Plainkeywords{differential equation solvers, compiled code, performance, FORTRAN, C} \author{ Karline Soetaert\\ Royal Netherlands Institute\\ of Sea Research (NIOZ)\\ Yerseke\\ The Netherlands \And Thomas Petzoldt\\ Technische Universit\"at \\ Dresden\\ Germany \And R. Woodrow Setzer\\ National Center for\\ Computational Toxicology\\ US Environmental Protection Agency } \Plainauthor{Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer} \Abstract{This document describes how to use the \pkg{deSolve} package \citep{deSolve_jss} to solve models that are written in \proglang{FORTRAN} or \proglang{C}.} %% The address of (at least) one author should be given %% in the following format: \Address{ Karline Soetaert\\ Royal Netherlands Institute of Sea Research (NIOZ)\\ 4401 NT Yerseke, Netherlands \\ E-mail: \email{karline.soetaert@nioz.nl}\\ URL: \url{https://www.nioz.nl}\\ \\ Thomas Petzoldt\\ Institut f\"ur Hydrobiologie\\ Technische Universit\"at Dresden\\ 01062 Dresden, Germany\\ E-mail: \email{thomas.petzoldt@tu-dresden.de}\\ URL: \url{https://tu-dresden.de/Members/thomas.petzoldt/}\\ \\ R. Woodrow Setzer\\ National Center for Computational Toxicology\\ US Environmental Protection Agency } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands. %% need no \usepackage{Sweave} %\VignetteIndexEntry{R Package deSolve: Writing Code in Compiled Languages} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document \begin{document} \SweaveOpts{engine=R,eps=FALSE} <>= library("deSolve") options(prompt = "R> ") options(width=70) @ \maketitle \section{Introduction} \pkg{deSolve} \citep{deSolve_jss,deSolve}, the successor of \proglang{R} package \pkg{odesolve} \citep{Setzer01} is a package to solve ordinary differential equations (ODE), differential algebraic equations (DAE) and partial differential equations (PDE). One of the prominent features of \pkg{deSolve} is that it allows specifying the differential equations either as: \begin{itemize} \item pure \proglang{R} code \citep{Rcore}, \item functions defined in lower-level languages such as \proglang{FORTRAN}, \proglang{C}, or \proglang{C++}, which are compiled into a dynamically linked library (DLL) and loaded into \proglang{R}. \end{itemize} In what follows, these implementations will be referred to as \Rmodels and \DLLmodels respectively. Whereas \Rmodels are easy to implement, they allow simple interactive development, produce highly readible code and access to \proglang{R}s high-level procedures, \DLLmodels have the benefit of increased simulation speed. Depending on the problem, there may be a gain of up to several orders of magnitude computing time when using compiled code. Here are some rules of thumb when it is worthwhile or not to switch to \DLLmodels: \begin{itemize} \item As long as one makes use only of \proglang{R}s high-level commands, the time gain will be modest. This was demonstrated in \citet{deSolve_jss}, where a formulation of two interacting populations dispersing on a 1-dimensional or a 2-dimensional grid led to a time gain of a factor two only when using \DLLmodels. \item Generally, the more statements in the model, the higher will be the gain of using compiled code. Thus, in the same paper \citep{deSolve_jss}, a very simple, 0-D, Lotka-Volterrra type of model describing only 2 state variables was solved 50 times faster when using compiled code. \item As even \Rmodels are quite performant, the time gain induced by compiled code will often not be discernible when the model is only solved once (who can grasp the difference between a run taking 0.001 or 0.05 seconds to finish). However, if the model is to be applied multiple times, e.g. because the model is to be fitted to data, or its sensitivity is to be tested, then it may be worthwhile to implement the model in a compiled language. \end{itemize} Starting from \pkg{deSolve} version 1.4, it is now also possible to use \emph{forcing functions} in compiled code. These forcing functions are automatically updated by the integrators. See last chapter. \section{A simple ODE example} Assume the following simple ODE (which is from the \code{LSODA} source code): \begin{align*} \frac{{dy_1}}{{dt}} &= - k_1 \cdot y_1 + k_2 \cdot y_2 \cdot y_3 \\ \frac{{dy_2}}{{dt}} &= k_1 \cdot y_1 - k_2 \cdot y_2 \cdot y_3 - k_3 \cdot y_2 \cdot y_2 \\ \frac{{dy_3}}{{dt}} &= k_3 \cdot y_2 \cdot y_2 \\ \end{align*} where $y_1$, $y_2$ and $y_3$ are state variables, and $k_1$, $k_2$ and $k_3$ are parameters. We first implement and run this model in pure \proglang{R}, then show how to do this in \proglang{C} and in \proglang{FORTRAN}. \subsection{ODE model implementation in R} An ODE model implemented in \textbf{pure \proglang{R}} should be defined as: \begin{verbatim} yprime = func(t, y, parms, ...) \end{verbatim} where \code{t} is the current time point in the integration, \code{y} is the current estimate of the variables in the ODE system, and \code{parms} is a vector or list containing the parameter values. The optional dots argument (\code{\dots}) can be used to pass any other arguments to the function. The return value of \code{func} should be a list, whose first element is a vector containing the derivatives of \code{y} with respect to time, and whose next elements contain output variables that are required at each point in time. The \proglang{R} implementation of the simple ODE is given below: <>= model <- function(t, Y, parameters) { with (as.list(parameters),{ dy1 = -k1*Y[1] + k2*Y[2]*Y[3] dy3 = k3*Y[2]*Y[2] dy2 = -dy1 - dy3 list(c(dy1, dy2, dy3)) }) } @ The Jacobian ($\frac{{\partial y'}}{{\partial y}}$) associated to the above example is: <>= jac <- function (t, Y, parameters) { with (as.list(parameters),{ PD[1,1] <- -k1 PD[1,2] <- k2*Y[3] PD[1,3] <- k2*Y[2] PD[2,1] <- k1 PD[2,3] <- -PD[1,3] PD[3,2] <- k3*Y[2] PD[2,2] <- -PD[1,2] - PD[3,2] return(PD) }) } @ This model can then be run as follows: <>= parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) Y <- c(1.0, 0.0, 0.0) times <- c(0, 0.4*10^(0:11)) PD <- matrix(nrow = 3, ncol = 3, data = 0) out <- ode(Y, times, model, parms = parms, jacfunc = jac) @ \subsection{ODE model implementation in C} \label{sec:Cexamp} In order to create compiled models (.DLL = dynamic link libraries on Windows or .so = shared objects on other systems) you must have a recent version of the GNU compiler suite installed, which is quite standard for Linux. Windows users find all the required tools on \url{https://cran.r-project.org/bin/windows/Rtools/}. Getting DLLs produced by other compilers to communicate with R is much more complicated and therefore not recommended. More details can be found on \url{https://cran.r-project.org/doc/manuals/R-admin.html}. The call to the derivative and Jacobian function is more complex for compiled code compared to \proglang{R}-code, because it has to comply with the interface needed by the integrator source codes. Below is an implementation of this model in \proglang{C}: \verbatiminput{mymod.c} The implementation in \proglang{C} consists of three parts: \begin{enumerate} \item After defining the parameters in global \proglang{C}-variables, through the use of \code{\#define} statements, a function called \code{initmod} initialises the parameter values, passed from the \proglang{R}-code. This function has as its sole argument a pointer to \proglang{C}-function \code{odeparms} that fills a double array with double precision values, to copy the parameter values into the global variable. \item Function \code{derivs} then calculates the values of the derivatives. The derivative function is defined as: \begin{verbatim} void derivs (int *neq, double *t, double *y, double *ydot, double *yout, int *ip) \end{verbatim} where \code{*neq} is the number of equations, \code{*t} is the value of the independent variable, \code{*y} points to a double precision array of length \code{*neq} that contains the current value of the state variables, and \code{*ydot} points to an array that will contain the calculated derivatives. \code{*yout} points to a double precision vector whose first \code{nout} values are other output variables (different from the state variables \code{y}), and the next values are double precision values as passed by parameter \code{rpar} when calling the integrator. The key to the elements of \code{*yout} is set in \code{*ip} \code{*ip} points to an integer vector whose length is at least 3; the first element (\code{ip[0]}) contains the number of output values (which should be equal or larger than \code{nout}), its second element contains the length of \code{*yout}, and the third element contains the length of \code{*ip}; next are integer values, as passed by parameter \code{ipar} when calling the integrator.\footnote{Readers familiar with the source code of the \pkg{ODEPACK} solvers may be surprised to find the double precision vector \code{yout} and the integer vector \code{ip} at the end. Indeed none of the \pkg{ODEPACK} functions allow this, although it is standard in the \code{vode} and \code{daspk} codes. To make all integrators compatible, we have altered the \pkg{ODEPACK} \proglang{FORTRAN} codes to consistently pass these vectors.} Note that, in function \code{derivs}, we start by checking whether enough memory is allocated for the output variables (\code{if (ip[0] < 1)}), else an error is passed to \proglang{R} and the integration is stopped. \item In \proglang{C}, the call to the function that generates the Jacobian is as: \begin{verbatim} void jac(int *neq, double *t, double *y, int *ml, int *mu, double *pd, int *nrowpd, double *yout, int *ip) \end{verbatim} where \code{*ml} and \code{*mu} are the number of non-zero bands below and above the diagonal of the Jacobian respectively. These integers are only relevant if the option of a banded Jacobian is selected. \code{*nrow} contains the number of rows of the Jacobian. Only for full Jacobian matrices, is this equal to \code{*neq}. In case the Jacobian is banded, the size of \code{*nrowpd} depends on the integrator. If the method is one of \code{lsode, lsoda, vode}, then \code{*nrowpd} will be equal to \code{*mu + 2 * *ml + 1}, where the last \code{*ml} rows should be filled with $0$s. For \code{radau}, \code{*nrowpd} will be equal to \code{*mu + *ml + 1} See example \url{dynload/odeband.R} in the directory \texttt{doc/dynload}, and chapter \ref{band}. \end{enumerate} \subsection{ODE model implementation in FORTRAN} \label{sec:forexamp} Models may also be defined in \proglang{FORTRAN}. \verbatiminput{mymod.f} In \proglang{FORTRAN}, parameters may be stored in a common block (here called \code{myparms}). During the initialisation, this common block is defined to consist of a 3-valued vector (unnamed), but in the subroutines \code{derivs} and \code{jac}, the parameters are given a name (\code{k1}, ...). \subsection{Running ODE models implemented in compiled code} To run the models described above, the code in \code{mymod.f} and \code{mymod.c} must first be compiled\footnote{This requires a correctly installed GNU compiler, see above.}. This can simply be done in \proglang{R} itself, using the \code{system} command: <>= system("R CMD SHLIB mymod.f") @ for the \proglang{FORTRAN} code or <>= system("R CMD SHLIB mymod.c") @ for the \proglang{C} code. This will create file \code{mymod.dll} on windows, or \code{mymod.so} on other platforms. We load the DLL, in windows as: \begin{verbatim} dyn.load("mymod.dll") \end{verbatim} and in unix: \begin{verbatim} dyn.load("mymod.so") \end{verbatim} or, using a general statement: \begin{verbatim} dyn.load(paste("mymod", .Platform$dynlib.ext, sep = "")) \end{verbatim} The model can now be run as follows: \begin{verbatim} parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) Y <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0) times <- c(0, 0.4*10^(0:11) ) out <- ode(Y, times, func = "derivs", parms = parms, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "Sum") \end{verbatim} The integration routine (here \code{ode}) recognizes that the model is specified as a DLL due to the fact that arguments \code{func} and \code{jacfunc} are not regular \proglang{R}-functions but character strings. Thus, the integrator will check whether the function is loaded in the DLL with name \code{mymod}. Note that \code{mymod}, as specified by \code{dllname} gives the name of the shared library \emph{without extension}. This DLL should contain all the compiled function or subroutine definitions referred to in \code{func}, \code{jacfunc} and \code{initfunc}. Also, if \code{func} is specified in compiled code, then \code{jacfunc} and \code{initfunc} (if present) should also be specified in a compiled language. It is not allowed to mix \proglang{R}-functions and compiled functions. Note also that, when invoking the integrator, we have to specify the number of ordinary output variables, \code{nout}. This is because the integration routine has to allocate memory to pass these output variables back to \proglang{R}. There is no way to check for the number of output variables in a DLL automatically. If in the calling of the integration routine the number of output variables is too low, then \proglang{R} may freeze and need to be terminated! Therefore it is advised that one checks in the code whether \code{nout} has been specified correctly. In the \proglang{FORTRAN} example above, the statement \code{if (ip(1) < 1) call rexit("nout should be at least 1")} does this. Note that it is not an error (just a waste of memory) to set \code{nout} to a too large value. Finally, in order to label the output matrix, the names of the ordinary output variables have to be passed explicitly (\code{outnames}). This is not necessary for the state variables, as their names are known through their initial condition (\code{y}). \section{Alternative way of passing parameters and data in compiled code} \label{sec:parms} All of the solvers in \pkg{deSolve} take an argument \code{parms} which may be an arbitrary \proglang{R} object. In models defined in \proglang{R} code, this argument is passed unprocessed to the various functions that make up the model. It is possible, as well, to pass such R-objects to models defined in native code. The problem is that data passed to, say, \code{ode} in the argument \code{parms} is not visible by default to the routines that define the model. This is handled by a user-written initialization function, for example \code{initmod} in the \proglang{C} and \proglang{FORTRAN} examples from sections \ref{sec:Cexamp} and \ref{sec:forexamp}. However, these set only the \emph{values} of the parameters. R-objects have many attributes that may also be of interest. To have access to these, we need to do more work, and this mode of passing parameters and data is much more complex than what we saw in previous chapters. In \proglang{C}, the initialization routine is declared: \begin{verbatim} void initmod(void (* odeparms)(int *, double *)); \end{verbatim} That is, \code{initmod} has a single argument, a pointer to a function that has as arguments a pointer to an \texttt{int} and a pointer to a \texttt{double}. In \proglang{FORTRAN}, the initialization routine has a single argument, a subroutine declared to be external. The name of the initialization function is passed as an argument to the \pkg{deSolve} solver functions. In \proglang{C}, two approaches are available for making the values passed in \code{parms} visible to the model routines, while only the simpler approach is available in \proglang{FORTRAN}. The simpler requires that \code{parms} be a numeric vector. In \proglang{C}, the function passed from \pkg{deSolve} to the initialization function (called \code{odeparms} in the example) copies the values from the parameter vector to a static array declared globally in the file where the model is defined. In \proglang{FORTRAN}, the values are copied into a \code{COMMON} block. It is possible to pass more complicated structures to \proglang{C} functions. Here is an example, an initializer called \code{deltamethrin} from a model describing the pharmacokinetics of that pesticide: \begin{verbatim} #include #include #include #include "deltamethrin.h" /* initializer */ void deltamethrin(void(* odeparms)(int *, double *)) { int Nparms; DL_FUNC get_deSolve_gparms; SEXP gparms; get_deSolve_gparms = R_GetCCallable("deSolve","get_deSolve_gparms"); gparms = get_deSolve_gparms(); Nparms = LENGTH(gparms); if (Nparms != N_PARMS) { PROBLEM "Confusion over the length of parms" ERROR; } else { _RDy_deltamethrin_parms = REAL(gparms); } } \end{verbatim} In \texttt{deltamethrin.h}, the variable \code{\_RDy\_deltamethrin\_parms} and macro N\_PARMS are declared: \begin{verbatim} #define N_PARMS 63 static double *_RDy_deltamethrin_parms; \end{verbatim} The critical element of this method is the function \code{R\_GetCCallable} which returns a function (called \code{get\_deSolve\_gparms} in this implementation) that returns the parms argument as a \code{SEXP} data type. In this example, \code{parms} was just a real vector, but in principle, this method can handle arbitrarily complex objects. For more detail on handling \proglang{R} objects in native code, see \proglang{R} Development Core Team (2008). \section{deSolve integrators that support DLL models} In the most recent version of \pkg{deSolve} all integration routines can solve \DLLmodels. They are: \begin{itemize} \item all solvers of the \code{lsode} familiy: \code{lsoda}, \code{lsode}, \code{lsodar}, \code {lsodes}, \item \code{vode}, \code{zvode}, \item \code{daspk}, \item \code{radau}, \item the Runge-Kutta integration routines (including the Euler method). \end{itemize} For some of these solvers the interface is slightly different (e.g. \code{zvode, daspk}), while in others (\code{lsodar}, \code{lsodes}) different functions can be defined. How this is implemented in a compiled language is discussed next. \subsection{Complex numbers, function zvode} \code{zvode} solves ODEs that are composed of complex variables. The program below uses \code{zvode} to solve the following system of 2 ODEs: \begin{align*} \frac{{dz}}{{dt}} &= i \cdot z\\ \frac{{dw}}{{dt}} &= -i \cdot w \cdot w \cdot z\\ \end{align*} where \begin{align*} w(0) = 1/2.1 +0i\\ z(0) = 1i \end{align*} on the interval t = [0, 2 $\pi$] The example is implemented in \proglang{FORTRAN}% \footnote{this can be found in file "zvodedll.f", in the \texttt{doc/dynload} subdirectory of the package}, \code{FEX} implements the function \code{func}: \begin{verbatim} SUBROUTINE FEX (NEQ, T, Y, YDOT, RPAR, IPAR) INTEGER NEQ, IPAR(*) DOUBLE COMPLEX Y(NEQ), YDOT(NEQ), RPAR(*), CMP DOUBLE PRECISION T character(len=100) msg c the imaginary unit i CMP = DCMPLX(0.0D0,1.0D0) YDOT(1) = CMP*Y(1) YDOT(2) = -CMP*Y(2)*Y(2)*Y(1) RETURN END \end{verbatim} \code{JEX} implements the function \code{jacfunc} \begin{verbatim} SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD, RPAR, IPAR) INTEGER NEQ, ML, MU, NRPD, IPAR(*) DOUBLE COMPLEX Y(NEQ), PD(NRPD,NEQ), RPAR(*), CMP DOUBLE PRECISION T c the imaginary unit i CMP = DCMPLX(0.0D0,1.0D0) PD(2,3) = -2.0D0*CMP*Y(1)*Y(2) PD(2,1) = -CMP*Y(2)*Y(2) PD(1,1) = CMP RETURN END \end{verbatim} Assuming this code has been compiled and is in a DLL called "zvodedll.dll", this model is solved in R as follows: \begin{verbatim} dyn.load("zvodedll.dll") outF <- zvode(func = "fex", jacfunc = "jex", y = yini, parms = NULL, times = times, atol = 1e-10, rtol = 1e-10, dllname = "zvodedll", initfunc = NULL) \end{verbatim} Note that in \proglang{R} names of \proglang{FORTRAN} DLL functions (e.g. for \code{func} and \code{jacfunc}) have to be given in lowercase letters, even if they are defined upper case in \proglang{FORTRAN}. Also, there is no initialiser function here (\code{initfunc = NULL}). \subsection{DAE models, integrator daspk} \code{daspk} is one of the integrators in the package that solve DAE models. In order to be used with DASPK, DAEs are specified in implicit form: \[0 = F(t, y, y', p)\] i.e. the DAE function (passed via argument \code{res}) specifies the ``residuals'' rather than the derivatives (as for ODEs). Consequently the DAE function specification in a compiled language is also different. For code written in \proglang{C}, the calling sequence for \code{res} must be: \begin{verbatim} void myres(double *t, double *y, double *ydot, double *cj, double *delta, int *ires, double *yout, int *ip) \end{verbatim} where \code{*t} is the value of the independent variable, \code{*y} points to a double precision vector that contains the current value of the state variables, \code{*ydot} points to an array that will contain the derivatives, \code{*delta} points to a vector that will contain the calculated residuals. \code{*cj} points to a scalar, which is normally proportional to the inverse of the stepsize, while \code{*ires} points to an integer (not used). \code{*yout} points to any other output variables (different from the state variables y), followed by the double precision values as passed via argument \code{rpar}; finally \code{*ip} is an integer vector containing at least 3 elements, its first value (\code{*ip[0]}) equals the number of output variables, calculated in the function (and which should be equal to \code{nout}), its second element equals the total length of \code{*yout}, its third element equals the total length of \code{*ip}, and finally come the integer values as passed via argument \code{ipar}. For code written in \proglang{FORTRAN}, the calling sequence for \code{res} must be as in the following example: \begin{verbatim} subroutine myresf(t, y, ydot, cj, delta, ires, out, ip) integer :: ires, ip(*) integer, parameter :: neq = 3 double precision :: t, y(neq), ydot(neq), delta(neq), out(*) double precision :: K, ka, r, prod, ra, rb common /myparms/K,ka,r,prod if(ip(1) < 1) call rexit("nout should be at least 1") ra = ka* y(3) rb = ka/K *y(1) * y(2) !! residuals of rates of changes delta(3) = -ydot(3) - ra + rb + prod delta(1) = -ydot(1) + ra - rb delta(2) = -ydot(2) + ra - rb - r*y(2) out(1) = y(1) + y(2) + y(3) return end \end{verbatim} Similarly as for the ODE model discussed above, the parameters are kept in a common block which is initialised by an initialiser subroutine: \begin{verbatim} subroutine initpar(daspkparms) external daspkparms integer, parameter :: N = 4 double precision parms(N) common /myparms/parms call daspkparms(N, parms) return end \end{verbatim} See the ODE example for how to initialise parameter values in \proglang{C}. Similarly, the function that specifies the Jacobian in a DAE differs from the Jacobian when the model is an ODE. The DAE Jacobian is set with argument \code{jacres} rather than \code{jacfunc} when an ODE. For code written in \proglang{FORTRAN}, the \code{jacres} must be as: \begin{verbatim} subroutine resjacfor (t, y, dy, pd, cj, out, ipar) integer, parameter :: neq = 3 integer :: ipar(*) double precision :: K, ka, r, prod double precision :: pd(neq,neq),y(neq),dy(neq),out(*) common /myparms/K,ka,r,prod !res1 = -dD - ka*D + ka/K *A*B + prod PD(1,1) = ka/K *y(2) PD(1,2) = ka/K *y(1) PD(1,3) = -ka -cj !res2 = -dA + ka*D - ka/K *A*B PD(2,1) = -ka/K *y(2) -cj PD(2,2) = -ka/K *y(2) PD(2,3) = ka !res3 = -dB + ka*D - ka/K *A*B - r*B PD(3,1) = -ka/K *y(2) PD(3,2) = -ka/K *y(2) -r -cj PD(3,3) = ka return end \end{verbatim} \subsection{DAE models, integrator radau} Function \code{radau} solves DAEs in linearly implicit form, i.e. in the form $M y' = f(t, y, p)$. The derivative function $f$ is specified in the same way as for an ODE, i.e. \begin{verbatim} void derivs (int *neq, double *t, double *y, double *ydot, double *yout, int *ip) \end{verbatim} and \begin{verbatim} subroutine derivs (neq, t, y, ydot, out, IP) \end{verbatim} for \proglang{C} and \proglang{FORTRAN} code respectively. To show how it should be used, we implement the caraxis problem as in \citep{testset}. The implementation of this index 3 DAE, comprising 8 differential, and 2 algebraic equations in R is the last example of the \code{radau} help page. We first repeat the R implementation: <<>>= caraxisfun <- function(t, y, parms) { with(as.list(c(y, parms)), { yb <- r * sin(w * t) xb <- sqrt(L * L - yb * yb) Ll <- sqrt(xl^2 + yl^2) Lr <- sqrt((xr - xb)^2 + (yr - yb)^2) dxl <- ul; dyl <- vl; dxr <- ur; dyr <- vr dul <- (L0-Ll) * xl/Ll + 2 * lam2 * (xl-xr) + lam1*xb dvl <- (L0-Ll) * yl/Ll + 2 * lam2 * (yl-yr) + lam1*yb - k * g dur <- (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr) dvr <- (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k * g c1 <- xb * xl + yb * yl c2 <- (xl - xr)^2 + (yl - yr)^2 - L * L list(c(dxl, dyl, dxr, dyr, dul, dvl, dur, dvr, c1, c2)) }) } eps <- 0.01; M <- 10; k <- M * eps^2/2; L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 1 parameter <- c(eps = eps, M = M, k = k, L = L, L0 = L0, r = r, w = w, g = g) yini <- c(xl = 0, yl = L0, xr = L, yr = L0, ul = -L0/L, vl = 0, ur = -L0/L, vr = 0, lam1 = 0, lam2 = 0) # the mass matrix Mass <- diag(nrow = 10, 1) Mass[5,5] <- Mass[6,6] <- Mass[7,7] <- Mass[8,8] <- M * eps * eps/2 Mass[9,9] <- Mass[10,10] <- 0 Mass # index of the variables: 4 of index 1, 4 of index 2, 2 of index 3 index <- c(4, 4, 2) times <- seq(0, 3, by = 0.01) out <- radau(y = yini, mass = Mass, times = times, func = caraxisfun, parms = parameter, nind = index) @ <>= plot(out, which = 1:4, type = "l", lwd = 2) @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the caraxis model - see text for R-code} \label{fig:caraxis} \end{figure} The implementation in \proglang{FORTRAN} consists of an initialiser function and a derivative function. \begin{verbatim} c---------------------------------------------------------------- c Initialiser for parameter common block c---------------------------------------------------------------- subroutine initcaraxis(daeparms) external daeparms integer, parameter :: N = 8 double precision parms(N) common /myparms/parms call daeparms(N, parms) return end c---------------------------------------------------------------- c rate of change c---------------------------------------------------------------- subroutine caraxis(neq, t, y, ydot, out, ip) implicit none integer neq, IP(*) double precision t, y(neq), ydot(neq), out(*) double precision eps, M, k, L, L0, r, w, g common /myparms/ eps, M, k, L, L0, r, w, g double precision xl, yl, xr, yr, ul, vl, ur, vr, lam1, lam2 double precision yb, xb, Ll, Lr, dxl, dyl, dxr, dyr double precision dul, dvl, dur, dvr, c1, c2 c expand state variables xl = y(1) yl = y(2) xr = y(3) yr = y(4) ul = y(5) vl = y(6) ur = y(7) vr = y(8) lam1 = y(9) lam2 = y(10) yb = r * sin(w * t) xb = sqrt(L * L - yb * yb) Ll = sqrt(xl**2 + yl**2) Lr = sqrt((xr - xb)**2 + (yr - yb)**2) dxl = ul dyl = vl dxr = ur dyr = vr dul = (L0-Ll) * xl/Ll + 2 * lam2 * (xl-xr) + lam1*xb dvl = (L0-Ll) * yl/Ll + 2 * lam2 * (yl-yr) + lam1*yb - k*g dur = (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr) dvr = (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k*g c1 = xb * xl + yb * yl c2 = (xl - xr)**2 + (yl - yr)**2 - L * L c function values in ydot ydot(1) = dxl ydot(2) = dyl ydot(3) = dxr ydot(4) = dyr ydot(5) = dul ydot(6) = dvl ydot(7) = dur ydot(8) = dvr ydot(9) = c1 ydot(10) = c2 return end \end{verbatim} Assuming that the code is in file ``radaudae.f'', this model is compiled, loaded and solved in R as: \begin{verbatim} system("R CMD SHLIB radaudae.f") dyn.load(paste("radaudae", .Platform$dynlib.ext, sep = "")) outDLL <- radau(y = yini, mass = Mass, times = times, func = "caraxis", initfunc = "initcaraxis", parms = parameter, dllname = "radaudae", nind = index) dyn.unload(paste("radaudae", .Platform$dynlib.ext, sep = "")) \end{verbatim} \subsection{The root function from integrators lsodar and lsode} \code{lsodar} is an extended version of integrator \code{lsoda} that includes a root finding function. This function is spedified via argument \code{rootfunc}. In \code{deSolve} version 1.7, \code{lsode} has also been extended with root finding capabilities. Here is how to program such a function in a lower-level language. For code written in \proglang{C}, the calling sequence for \code{rootfunc} must be: \begin{verbatim} void myroot(int *neq, double *t, double *y, int *ng, double *gout, double *out, int *ip ) \end{verbatim} where \code{*neq} and \code{*ng} are the number of state variables and root functions respectively, \code{*t} is the value of the independent variable, \code{y} points to a double precision array that contains the current value of the state variables, and \code{gout} points to an array that will contain the values of the constraint function whose root is sought. \code{*out} and \code{*ip} are a double precision and integer vector respectively, as described in the ODE example above. For code written in \proglang{FORTRAN}, the calling sequence for \code{rootfunc} must be as in following example: \begin{verbatim} subroutine myroot(neq, t, y, ng, gout, out, ip) integer :: neq, ng, ip(*) double precision :: t, y(neq), gout(ng), out(*) gout(1) = y(1) - 1.e-4 gout(2) = y(3) - 1e-2 return end \end{verbatim} \subsection{jacvec, the Jacobian vector for integrator lsodes} Finally, in integration function \code{lsodes}, not the Jacobian \emph{matrix} is specified, but a \emph{vector}, one for each column of the Jacobian. This function is specified via argument \code{jacvec}. In \proglang{FORTRAN}, the calling sequence for \code{jacvec} is: \begin{verbatim} SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ, OUT, IP) DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*), OUT(*) INTEGER NEQ, J, IP(*) \end{verbatim} \subsection{Banded jacobians in compiled code}\label{band} In the call of the jacobian function, the number of bands below and above the diagonal (\code{ml, mu}) and the number of rows of the Jacobian matrix, \code{nrowPD} is specified, e.g. for \proglang{FORTRAN} code: \begin{verbatim} SUBROUTINE JAC (neq, T, Y, ml, mu, PD, nrowPD, RPAR, IPAR) \end{verbatim} The jacobian matrix to be returned should have dimension \code{nrowPD, neq}. In case the Jacobian is banded, the size of \code{nrowPD} depends on the integrator. If the method is one of \code{lsode, lsoda, vode}, or related, then \code{nrowPD} will be equal to \code{mu + 2 * ml + 1}, where the last ml rows should be filled with $0$s. For \code{radau}, \code{nrowpd} will be equal to \code{mu + ml + 1} Thus, it is important to write the FORTRAN or C-code in such a way that it can be used with both types of integrators - else it is likely that R will freeze if the wrong integrator is used. We implement in FORTRAN, the example of the \code{lsode} help file. The R-code reads: <<>>= ## the model, 5 state variables f1 <- function (t, y, parms) { ydot <- vector(len = 5) ydot[1] <- 0.1*y[1] -0.2*y[2] ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3] ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4] ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5] ydot[5] <- -0.3*y[4] +0.1*y[5] return(list(ydot)) } ## the Jacobian, written in banded form bandjac <- function (t, y, parms) { jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE, data = c( 0 , -0.2, -0.2, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, -0.3, -0.3, -0.3, -0.3, 0)) return(jac) } ## initial conditions and output times yini <- 1:5 times <- 1:20 ## stiff method, user-generated banded Jacobian out <- lsode(yini, times, f1, parms = 0, jactype = "bandusr", jacfunc = bandjac, bandup = 1, banddown = 1) @ In FORTRAN, the code might look like this: \begin{verbatim} c Rate of change subroutine derivsband (neq, t, y, ydot,out,IP) integer neq, IP(*) DOUBLE PRECISION T, Y(5), YDOT(5), out(*) ydot(1) = 0.1*y(1) -0.2*y(2) ydot(2) = -0.3*y(1) +0.1*y(2) -0.2*y(3) ydot(3) = -0.3*y(2) +0.1*y(3) -0.2*y(4) ydot(4) = -0.3*y(3) +0.1*y(4) -0.2*y(5) ydot(5) = -0.3*y(4) +0.1*y(5) RETURN END c The banded jacobian subroutine jacband (neq, t, y, ml, mu, pd, nrowpd, RP, IP) INTEGER neq, ml, mu, nrowpd, ip(*) DOUBLE PRECISION T, Y(5), PD(nrowpd,5), rp(*) PD(:,:) = 0.D0 PD(1,1) = 0.D0 PD(1,2) = -.02D0 PD(1,3) = -.02D0 PD(1,4) = -.02D0 PD(1,5) = -.02D0 PD(2,:) = 0.1D0 PD(3,1) = -0.3D0 PD(3,2) = -0.3D0 PD(3,3) = -0.3D0 PD(3,4) = -0.3D0 PD(3,5) = 0.D0 RETURN END \end{verbatim} Assuming that this code is in file \code{"odeband.f"}, we compile from within R and load the shared library (assuming the working directory holds the source file) with: \begin{verbatim} system("R CMD SHLIB odeband.f") dyn.load(paste("odeband", .Platform$dynlib.ext, sep = "")) \end{verbatim} To solve this problem, we write in R \begin{verbatim} out2 <- lsode(yini, times, "derivsband", parms = 0, jactype = "bandusr", jacfunc = "jacband", bandup = 1, banddown = 1, dllname = "odeband") out2 <- radau(yini, times, "derivsband", parms = 0, jactype = "bandusr", jacfunc = "jacband", bandup = 1, banddown = 1, dllname = "odeband") \end{verbatim} This will work both for the \code{lsode} family as for \code{radau}. In the first case, when entering subroutine \code{jacband}, \code{nrowpd} will have the value $5$, in the second case, it will be equal to $4$. \section{Testing functions written in compiled code} Two utilities have been included to test the function implementation in compiled code: \begin{itemize} \item \code{DLLfunc} to test the implementation of the derivative function as used in ODEs. This function returns the derivative $\frac{dy}{dt}$ and the output variables. \item \code{DLLres} to test the implementation of the residual function as used in DAEs. This function returns the residual function $\frac{dy}{dt}-f(y,t)$ and the output variables. \end{itemize} These functions serve no other purpose than to test whether the compiled code returns what it should. \subsection{DLLfunc} We test whether the ccl4 model, which is part of \code{deSolve} package, returns the proper rates of changes. (Note: see \code{example(ccl4model)} for a more comprehensive implementation) <<>>= ## Parameter values and initial conditions Parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32, 16.17, 281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671, 0.4027255, 1000, 0.02, 1.0, 3.8) yini <- c( AI=21, AAM=0, AT=0, AF=0, AL=0, CLT=0, AM=0 ) ## the rate of change DLLfunc(y = yini, dllname = "deSolve", func = "derivsccl4", initfunc = "initccl4", parms = Parms, times = 1, nout = 3, outnames = c("DOSE", "MASS", "CP") ) @ \subsection{DLLres} The deSolve package contains a FORTRAN implementation of the chemical model described above (section 4.1), where the production rate is included as a forcing function (see next section). Here we use \code{DLLres} to test it: <<>>= pars <- c(K = 1, ka = 1e6, r = 1) ## Initial conc; D is in equilibrium with A,B y <- c(A = 2, B = 3, D = 2*3/pars["K"]) ## Initial rate of change dy <- c(dA = 0, dB = 0, dD = 0) ## production increases with time prod <- matrix(nc=2,data=c(seq(0,100,by=10),seq(0.1,0.5,len=11))) DLLres(y=y,dy=dy,times=5,res="chemres", dllname="deSolve", initfunc="initparms", initforc="initforcs", parms=pars, forcings=prod, nout=2, outnames=c("CONC","Prod")) @ \section{Using forcing functions} Forcing functions in DLLs are implemented in a similar way as parameters. This means: \begin{itemize} \item They are initialised by means of an initialiser function. Its name should be passed to the solver via argument \code{initforc}. Similar as the parameter initialiser function, the function denoted by \code{initforc} has as its sole argument a pointer to the vector that contains the forcing funcion values in the compiled code. In case of \proglang{C} code, this will be a global vector; in case of \proglang{FORTRAN}, this will be a vector in a common block. The solver puts a pointer to this vector and updates the forcing functions in this memory area at each time step. Hence, within the compiled code, forcing functions can be assessed as if they are parameters (although, in contrast to the latter, their values will generally change). No need to update the values for the current time step; this has been done before entering the \code{derivs} function. \item The forcing function data series are passed to the integrator, via argument \code{forcings}; if there is only one forcing function data set, then a 2-columned matrix (time, value) will do; else the data should be passed as a list, containing (time, value) matrices with the individual forcing function data sets. Note that the data sets in this list should be \emph{in the same ordering} as the declaration of the forcings in the compiled code. \end{itemize} A number of options allow to finetune certain settings. They are in a list called \code{fcontrol} which can be supplied as argument when calling the solvers. The options are similar to the arguments from R function \code{approx}, howevers the default settings are often different. The following options can be specified: \begin{itemize} \item \code{method} specifies the interpolation method to be used. Choices are "linear" or "constant", the default is "linear", which means linear interpolation (same as \code{approx}) \item \code{rule}, an integer describing how interpolation is to take place \emph{outside} the interval [min(times), max(times)]. If \code{rule} is \code{1} then an error will be triggered and the calculation will stop if extrapolation is necessary. If it is \code{2}, the default, the value at the closest data extreme is used, a warning will be printed if \code{verbose} is TRUE. Note that the default differs from the \code{approx} default. \item \code{f}, for method=\code{"constant"} is a number between \code{0} and \code{1} inclusive, indicating a compromise between left- and right-continuous step functions. If \code{y0} and \code{y1} are the values to the left and right of the point then the value is \code{y0*(1-f)+y1*f} so that \code{f=0} is right-continuous and \code{f=1} is left-continuous. The default is to have \code{f=0}. For some data sets it may be more realistic to set \code{f=0.5}. \item \code{ties}, the handling of tied \code{times} values. Either a function with a single vector argument returning a single number result or the string "ordered". Note that the default is "ordered", hence the existence of ties will NOT be investigated; in practice this means that, if ties exist, the first value will be used; if the dataset is not ordered, then nonsense will be produced. Alternative values for \code{ties} are \code{mean}, \code{min} etc... which will average, or take the minimal value if multiple values exist at one time level. \end{itemize} The default settings of \code{fcontrol} are: \code{fcontrol=list(method="linear", rule = 2, f = 0, ties = "ordered")} Note that only ONE specification is allowed, even if there is more than one forcing function data set. (may/should change in the future). \subsection{A simple FORTRAN example} We implement the example from chapter 3 of the book \citep{Soetaert08} in FORTRAN. This model describes the oxygen consumption of a (marine) sediment in response to deposition of organic matter (the forcing function). One state variable, the organic matter content in the sediment is modeled; it changes as a function of the deposition \code{Flux} (forcing) and organic matter decay (first-order decay rate \code{k}). \[ \frac{dC}{dt}=Flux_t-k \cdot C \] with initial condition $C(t=0)=C_0$; the latter is estimated as the mean of the flux divided by the decay rate. The FORTRAN code looks like this: \begin{verbatim} c Initialiser for parameter common block subroutine scocpar(odeparms) external odeparms integer N double precision parms(2) common /myparms/parms N = 1 call odeparms(N, parms) return end c Initialiser for forcing common block subroutine scocforc(odeforcs) external odeforcs integer N double precision forcs(1) common /myforcs/forcs N = 1 call odeforcs(N, forcs) return end c Rate of change and output variables subroutine scocder (neq, t, y, ydot,out,IP) integer neq, IP(*) double precision t, y(neq), ydot(neq), out(*), k, depo common /myparms/k common /myforcs/depo if(IP(1) < 2) call rexit("nout should be at least 2") ydot(1) = -k*y(1) + depo out(1)= k*y(1) out(2)= depo return end \end{verbatim} Here the subroutine \code{scocpar} is business as usual; it initialises the parameter common block (there is only one parameter). Subroutine \code{odeforcs} does the same for the forcing function, which is also positioned in a common block, called \code{myforcs}. This common block is made available in the derivative subroutine (here called \code{scocder}), where the forcing function is named \code{depo}. At each time step, the integrator updates the value of this forcing function to the correct time point. In this way, the forcing functions can be used as if they are (time-varying) parameters. All that's left to do is to pass the forcing function data set and the name of the forcing function initialiser routine. This is how to do it in R. First the data are inputted: <<>>= Flux <- matrix(ncol=2,byrow=TRUE,data=c( 1, 0.654, 11, 0.167, 21, 0.060, 41, 0.070, 73,0.277, 83,0.186, 93,0.140,103, 0.255, 113, 0.231,123, 0.309,133,1.127,143,1.923, 153,1.091,163,1.001, 173, 1.691,183, 1.404,194,1.226,204,0.767, 214, 0.893,224,0.737, 234,0.772,244, 0.726,254,0.624,264,0.439, 274,0.168,284 ,0.280, 294,0.202,304, 0.193,315,0.286,325,0.599, 335, 1.889,345, 0.996,355,0.681,365,1.135)) head(Flux) @ and the parameter given a value (there is only one) <<>>= parms <- 0.01 @ The initial condition \code{Yini} is estimated as the annual mean of the Flux and divided by the decay rate (parameter). <<>>= meanDepo <- mean(approx(Flux[,1],Flux[,2], xout=seq(1,365,by=1))$y) Yini <- c(y=meanDepo/parms) @ After defining the output times, the model is run, using integration routine \code{ode}. The \emph{name} of the derivate function \code{"scocder"}, of the dll \code{"deSolve"}\footnote{this example is made part of the deSolve package, hence the name of the dll is "deSolve"} and of the initialiser function \code{"scocpar"} are passed, as in previous examples. In addition, the forcing function data set is also passed (\code{forcings=Flux}) as is the name of the forcing initialisation function (\code{initforc="scocforc"}). <<>>= times <- 1:365 out <- ode(y=Yini, times, func = "scocder", parms = parms, dllname = "deSolve", initforc="scocforc", forcings=Flux, initfunc = "scocpar", nout = 2, outnames = c("Mineralisation","Depo")) head(out) @ Now, the way the forcing functions are interpolated are changed: Rather than linear interpolation, constant (block, step) interpolation is used. <<>>= fcontrol <- list(method="constant") out2 <- ode(y=Yini, times, func = "scocder", parms = parms, dllname = "deSolve", initforc="scocforc", forcings=Flux, fcontrol=fcontrol, initfunc = "scocpar", nout = 2, outnames = c("Mineralisation","Depo")) @ Finally, the results are plotted: <>= par (mfrow=c(1,2)) plot(out, which = "Depo", col="red", xlab="days", ylab="mmol C/m2/ d", main="method='linear'") lines(out[,"time"], out[,"Mineralisation"], lwd=2, col="blue") legend("topleft",lwd=1:2,col=c("red","blue"), c("Flux","Mineralisation")) plot(out, which = "Depo", col="red", xlab="days", ylab="mmol C/m2/ d", main="method='constant'") lines(out2[,"time"], out2[,"Mineralisation"], lwd=2, col="blue") @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the SCOC model, implemented in compiled code, and including a forcing function - see text for R-code} \label{fig:scoc} \end{figure} \subsection{An example in C} Consider the following R-code which implements a resource-producer-consumer Lotka-Volterra type of model in R (it is a modified version of the example of function \code{ode}): <<>>= SPCmod <- function(t, x, parms, input) { with(as.list(c(parms, x)), { import <- input(t) dS <- import - b*S*P + g*C # substrate dP <- c*S*P - d*C*P # producer dC <- e*P*C - f*C # consumer res <- c(dS, dP, dC) list(res, signal = import) }) } ## The parameters parms <- c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0) ## vector of timesteps times <- seq(0, 100, by=0.1) ## external signal with several rectangle impulses signal <- as.data.frame(list(times = times, import = rep(0, length(times)))) signal$import <- ifelse((trunc(signal$times) %% 2 == 0), 0, 1) sigimp <- approxfun(signal$times, signal$import, rule = 2) ## Start values for steady state xstart <- c(S = 1, P = 1, C = 1) ## Solve model print (system.time( out <- ode(y = xstart,times = times, func = SPCmod, parms, input = sigimp) )) @ All output is printed at once: <>= plot(out) @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the Lotka-Volterra resource (S)-producer (P) - consumer (C) model with time-variable input (signal) - see text for R-code} \label{fig:lv} \end{figure} The C-code, in file \url{dynload/Forcing\_lv.c}, can be found in the packages \texttt{doc/dynload} subdirectory. It can be compiled, from within R by \begin{verbatim} system("R CMD SHLIB Forcing_lv.c") \end{verbatim} After defining the parameter and forcing vectors, and giving them comprehensible names, the parameter and forcing initialiser functions are defined (\code{parmsc} and \code{forcc} respectively). Next is the derivative function, \code{derivsc}. \begin{verbatim} #include static double parms[6]; static double forc[1]; /* A trick to keep up with the parameters and forcings */ #define b parms[0] #define c parms[1] #define d parms[2] #define e parms[3] #define f parms[4] #define g parms[5] #define import forc[0] /* initializers: */ void odec(void (* odeparms)(int *, double *)) { int N=6; odeparms(&N, parms); } void forcc(void (* odeforcs)(int *, double *)) { int N=1; odeforcs(&N, forc); } /* derivative function */ void derivsc(int *neq, double *t, double *y, double *ydot, double *yout, int*ip) { if (ip[0] <2) error("nout should be at least 2"); ydot[0] = import - b*y[0]*y[1] + g*y[2]; ydot[1] = c*y[0]*y[1] - d*y[2]*y[1]; ydot[2] = e*y[1]*y[2] - f*y[2]; yout[0] = y[0] + y[1] + y[2]; yout[1] = import; } \end{verbatim} After defining the forcing function time series, which is to be interpolated by the integration routine, and loading the DLL, the model is run: \begin{verbatim} Sigimp <- approx(signal$times, signal$import, xout=ftime,rule = 2)$y forcings <- cbind(ftime,Sigimp) dyn.load("Forcing_lv.dll") out <- ode(y=xstart, times, func = "derivsc", parms = parms, dllname = "Forcing_lv",initforc = "forcc", forcings=forcings, initfunc = "parmsc", nout = 2, outnames = c("Sum","signal"), method = rkMethod("rk34f")) dyn.unload("Forcing_lv.dll") \end{verbatim} This code executes about 30 times faster than the \proglang{R}-code. With a longer simulation time, the difference becomes more pronounced, e.g. with times till 800 days, the DLL code executes 200 times faster% \footnote{this is due to the sequential update of the forcing functions by the solvers, compared to the bisectioning approach used by approxfun}. \section{Implementing events in compiled code} An \code{event} occurs when the value of a state variable is suddenly changed, e.g. a certain amount is added, or part is removed. The integration routines cannot deal easily with such state variable changes. Typically these events occur only at specific times. In \code{deSolve}, events can be imposed by means of an input file that specifies at which time a certain state variable is altered, or via an event function. Both types of events combine with compiled code. Take the previous example, the Lotka-Volterra SPC model. Suppose that every 10 days, half of the consumer is removed. We first implement these events as a \code{data.frame} <<>>= eventdata <- data.frame(var=rep("C",10),time=seq(10,100,10),value=rep(0.5,10), method=rep("multiply",10)) eventdata @ This model is solved, and plotted as: \begin{verbatim} dyn.load("Forcing_lv.dll") out2 <- ode(y = y, times, func = "derivsc", parms = parms, dllname = "Forcing_lv", initforc="forcc", forcings = forcings, initfunc = "parmsc", nout = 2, outnames = c("Sum", "signal"), events=list(data=eventdata)) dyn.unload("Forcing_lv.dll") plot(out2, which = c("S","P","C"), type = "l") \end{verbatim} The event can also be implemented in \proglang{C} as: \begin{verbatim} void event(int *n, double *t, double *y) { y[2] = y[2]*0.5; } \end{verbatim} Here n is the length of the state variable vector \code{y}. and is then solved as: \begin{verbatim} dyn.load("Forcing_lv.dll") out3 <- ode(y = y, times, func = "derivsc", parms = parms, dllname = "Forcing_lv", initforc="forcc", forcings = forcings, initfunc = "parmsc", nout = 2, outnames = c("Sum", "signal"), events = list(func="event",time=seq(10,90,10))) dyn.unload("Forcing_lv.dll") \end{verbatim} \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} \includegraphics{comp-event} \end{center} \caption{Solution of the Lotka-Volterra resource (S)~-- producer (P)~-- consumer (C) model with time-variable input (signal) and with half of the consumer removed every 10 days - see text for R-code} \label{fig:lv2} \end{figure} \section{Delay differential equations} It is now also very simple to implement delay differential equations in compiled code and solve them with \code{dede}. In order to do so, you need to get access to the R-functions \code{lagvalue} and \code{lagderiv} that will give you the past value of the state variable or its derivative respectively. \subsection{Delays implemented in Fortran} If you use \proglang{Fortran}, then the easiest way is to link your code with a file called \url{dynload-dede/dedeUtils.c} that you will find in the packages subdirectory \texttt{doc/dynload-dede}. This file contains Fortran-callable interfaces to the delay-differential utility functions from package \pkg{deSolve}, and that are written in \proglang{C}. Its content is: \begin{verbatim} void F77_SUB(lagvalue)(double *T, int *nr, int *N, double *ytau) { static void(*fun)(double, int*, int, double*) = NULL; if (fun == NULL) fun = (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagvalue"); fun(*T, nr, *N, ytau); return; } void F77_SUB(lagderiv)(double *T, int *nr, int *N, double *ytau) { static void(*fun)(double, int*, int, double*) = NULL; if (fun == NULL) fun = (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagderiv"); fun(*T, nr, *N, ytau); return; } \end{verbatim} Here \code{T} is the time at which the value needs to be retrieved, \code{nr} is an integer that defines the number of the state variable or its derivative whose delay we want, \code{N} is the total number of state variabes and \code{ytau} will have the result. We start with an example, a Lotka-Volterra system with delay, that we will implement in \proglang{Fortran} (you will find this example in the package directory \texttt{doc/dynload-dede}, in file \texttt{dede\_lvF.f} The R-code would be: <<>>= derivs <- function(t, y, parms) { with(as.list(c(y, parms)), { if (t < tau) ytau <- c(1, 1) else ytau <- lagvalue(t - tau, c(1, 2)) dN <- f * N - g * N * P dP <- e * g * ytau[1] * ytau[2] - m * P list(c(dN, dP), tau=ytau[1], tau=ytau[2]) }) } yinit <- c(N = 1, P = 1) times <- seq(0, 500) parms <- c(f = 0.1, g = 0.2, e = 0.1, m = 0.1, tau = .2) yout <- dede(y = yinit, times = times, func = derivs, parms = parms) head(yout) @ In Fortran the code looks like this: \begin{verbatim} ! file dede_lfF.f ! Initializer for parameter common block subroutine initmod(odeparms) external odeparms double precision parms(5) common /myparms/parms call odeparms(5, parms) return end ! Derivatives and one output variable subroutine derivs(neq, t, y, ydot, yout, ip) integer neq, ip(*) double precision t, y(neq), ydot(neq), yout(*) double precision N, P, ytau(2), tlag integer nr(2) double precision f, g, e, m, tau common /myparms/f, g, e, m, tau if (ip(1) < 2) call rexit("nout should be at least 2") N = y(1) P = y(2) nr(1) = 0 nr(2) = 1 ytau(1) = 1.0 ytau(2) = 1.0 tlag = t - tau if (tlag .GT. 0.0) call lagvalue(tlag, nr, 2, ytau) ydot(1) = f * N - g * N * P ydot(2) = e * g * ytau(1) * ytau(2) - m * P yout(1) = ytau(1) yout(2) = ytau(2) return end \end{verbatim} During compilation, we need to also compile the file \code{dedeUtils.c}. Assuming that the above \proglang{Fortran} code is in file \code{dede_lvF.f}, which is found in the working directory that also contains file \code{dedeUtils.c}, the problem is compiled and run as: \begin{verbatim} system("R CMD SHLIB dede_lvF.f dedeUtils.c") dyn.load(paste("dede_lvF", .Platform$dynlib.ext, sep="")) yout3 <- dede(yinit, times = times, func = "derivs", parms = parms, dllname = "dede_lvF", initfunc = "initmod", nout = 2) \end{verbatim} \subsection{Delays implemented in C} We now give the same example in \proglang{C}-code (you will find this file as \url{dynload-dede/dede_lv.c}). \begin{verbatim} #include #include #include #include static double parms[5]; #define f parms[0] #define g parms[1] #define e parms[2] #define m parms[3] #define tau parms[4] /* Interface to dede utility functions in package deSolve */ void lagvalue(double T, int *nr, int N, double *ytau) { static void(*fun)(double, int*, int, double*) = NULL; if (fun == NULL) fun = (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagvalue"); return fun(T, nr, N, ytau); } void lagderiv(double T, int *nr, int N, double *ytau) { static void(*fun)(double, int*, int, double*) = NULL; if (fun == NULL) fun = (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagderiv"); return fun(T, nr, N, ytau); } /* Initializer */ void initmod(void (* odeparms)(int *, double *)) { int N = 5; odeparms(&N, parms); } /* Derivatives */ void derivs (int *neq, double *t, double *y, double *ydot, double *yout, int *ip) { if (ip[0] < 2) error("nout should be at least 1"); double N = y[0]; double P = y[1]; int Nout = 2; // number of returned lags ( <= n_eq !!) int nr[2] = {0, 1}; // which lags are needed? // numbering starts from zero ! double ytau[2] = {1.0, 1.0}; // array; initialize with default values ! double T = *t - tau; if (*t > tau) { lagvalue(T, nr, Nout, ytau); } ydot[0] = f * N - g * N * P; ydot[1] = e * g * ytau[0] * ytau[1] - m * P; yout[0] = ytau[0]; yout[1] = ytau[1]; } \end{verbatim} Assuming this code is in a file called \code{dede_lv.c}, which is in the working directory, this file is then compiled and run as: \begin{verbatim} system("R CMD SHLIB dede_lv.c") dyn.load(paste("dede_lv", .Platform$dynlib.ext, sep="")) yout2 <- dede(yinit, times = times, func = "derivs", parms = parms, dllname = "dede_lv", initfunc = "initmod", nout = 2) dyn.unload(paste("dede_lv", .Platform$dynlib.ext, sep="")) \end{verbatim} \section{Difference equations in compiled code} There is one special-purpose solver, triggered with \code{method = "iteration"} which can be used in cases where the new values of the state variables are estimated by the user, and need not be found by integration. This is for instance useful when the model consists of difference equations, or for 1-D models when transport is implemented by an implicit or a semi-implicit method. An example of a discrete time model, represented by a difference equation is given in the help file of solver \code{ode}. It consists of the host-parasitoid model described as from \citet[p283]{Soetaert08}. We first give the R-code, and how it is solved: \begin{verbatim} Parasite <- function (t, y, ks) { P <- y[1] H <- y[2] f <- A * P / (ks +H) Pnew <- H* (1-exp(-f)) Hnew <- H * exp(rH*(1.-H) - f) list (c(Pnew, Hnew)) } rH <- 2.82 # rate of increase A <- 100 # attack rate ks <- 15. # half-saturation density out <- ode (func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = ks, method = "iteration") \end{verbatim} Note that the function returns the updated value of the state variables rather than the rate of change (derivative). The method ``iteration'' does not perform any integration. The implementation in \proglang{FORTRAN} consists of an initialisation function to pass the parameter values (\code{initparms}) and the "update" function that returns the new values of the state variables (\code{parasite}): \begin{verbatim} subroutine initparms(odeparms) external odeparms double precision parms(3) common /myparms/parms call odeparms(3, parms) return end subroutine parasite (neq, t, y, ynew, out, iout) integer neq, iout(*) double precision t, y(neq), ynew(neq), out(*), rH, A, ks common /myparms/ rH, A, ks double precision P, H, f P = y(1) H = y(2) f = A * P / (ks + H) ynew(1) = H * (1.d0 - exp(-f)) ynew(2) = H * exp (rH * (1.d0 - H) - f) return end \end{verbatim} The model is compiled, loaded and executed in R as: \begin{verbatim} system("R CMD SHLIB difference.f") dyn.load(paste("difference", .Platform$dynlib.ext, sep = "")) require(deSolve) rH <- 2.82 # rate of increase A <- 100 # attack rate ks <- 15. # half-saturation density parms <- c(rH = rH, A = A, ks = ks) out <- ode (func = "parasite", y = c(P = 0.5, H = 0.5), times = 0:50, initfunc = "initparms", dllname = "difference", parms = parms, method = "iteration") \end{verbatim} \section{Final remark} Detailed information about communication between \proglang{C}, \proglang{FORTRAN} and \proglang{R} can be found in \citet{Rexts2009}. Notwithstanding the speed gain when using compiled code, one should not carelessly decide to always resort to this type of modelling. Because the code needs to be formally compiled and linked to \proglang{R} much of the elegance when using pure \proglang{R} models is lost. Moreover, mistakes are easily made and paid harder in compiled code: often a programming error will terminate \proglang{R}. In addition, these errors may not be simple to trace. \clearpage %% this adds References to the PDF-Index without adding an obsolete section \phantomsection \addcontentsline{toc}{section}{References} \bibliography{integration} \end{document}