Movatterモバイル変換

odin

odin implements a high-level language for describing andimplementing ordinary differential equations in R. It provides a “domainspecific language” (DSL) whichlooks like R but is compileddirectly to C. The actual solution of the differential equations is donewith thedeSolvepackage, giving access to the excellent Livermore solvers(lsoda,lsode, etc), or withdde foruse with delay differential equations.

• The DSL isdeclarative reflecting the mathematical natureof the equations (typically ordinary differential equations are simplemathematical relationships, so the order should not matter).
• It includes support for equations that involve vectors, matrices andhigher dimensional arrays (up to 8!), including a high-level arrayindexing notation that removes the need for explicit looping.
• Delay differential equations are supported, including when thedelayed quantities are arbitrarily complicated expressions ofvariables.
• Interpolation functions can be used to include time-varyingquantities into the model (piecewise constant, linear and splineinterpolation is supported, usingcinterpolate.
• The equations are analysed before compilation so that parts that donot depend on time are not included in the final derivativecalculations.
• Supports user-supplied parameters for any part of the system.
• Supports a large number of mathematical functions (see thefunctionsvignette) for a complete list.

In addition, the same machinery can be used to generate discrete-timemodels that proceed over a set of steps (rather than through continuoustime). These may be stochastic and make use of any of R’s random numberfunctions.

odin works using code generation; the nice thing aboutthis approach is that it never gets bored. So if the generated code haslots of tedious repetitive bits, they’re at least likely to be correct(compared with implementing yourself).

Background

The “deSolve” package for R is the de-facto way of solvingdifferential equations in R; it provides excellent solvers and hasremained stable for over a decade. However, users must implementequations in R and suffer a large speed cost, or implement theirequations in C which is (depending on the complexity of the system)either routine and a bit boring, or complicated and error prone. Thistranslation can be especially complicated with delay differentialequations, or with models where the variables are more naturally storedas variable sized arrays.

Apparently not many people know thatdeSolve can usetarget functions written in C rather than just in R. This is describedin detail in the excellent “compiledCode” vignette(vignette("compiledCode") oronline.

While thedeSolve authors are bearish on the benefits ofthis, I have often seen performance improvements of over 100x. Where anODE is being used in application where it is called repeatedly (e.g., anoptimisation or MCMC) the cost of rewriting the system pays itselfback.

For simple systems the rewriting is essentially mechanical. Thelorenz attractor could be implemented in R as:

lorenz<-function(t, y, parms) {  sigma<- parms[1]  R<- parms[2]  b<- parms[3]  y1<- y[1]  y2<- y[2]  y3<- y[3]list(c(sigma* (y2- y1),         R* y1- y2- y1* y3,-b* y3+ y1* y2))}

and in C as

void initmod(void(* odeparms)(int*,double*)){int N=3;  odeparms(&N, parms);}void lorenz(int*n,double*t,double*y,double*dydt,double*yout,int*ip){double sigma= parms[0];double R= parms[1];double b= parms[2];double y1= y[0];double y2= y[2];double y3= y[3];  dydt[0]= sigma*(y2- y1);  dydt[1]= R* y1- y2- y1* y3;  dydt[2]=-b* y3+ y1* y2;}

The connection between the two languages should be fairly obvious. Assystems get more complicated much of the difficulty of writing thesystems in C becomes the tedium of book keeping as parameters and statevectors are unpacked, rather than any deep programming challenges.Modifying large systems is a particular challenge as technical debt canaccrue quickly.

The core job ofodin is to simplify this transition sothat models can be both developed and solved rapidly.

Example

The Lorenz attractor above can be implemented as:

lorenz<- odin::odin({## Derivativesderiv(y1)<- sigma* (y2- y1)deriv(y2)<- R* y1- y2- y1* y3deriv(y3)<--b* y3+ y1* y2## Initial conditionsinitial(y1)<-10.0initial(y2)<-1.0initial(y3)<-1.0## parameters  sigma<-10.0  R<-28.0  b<-8.0/3.0})

The connection to the R and C versions in the section above should befairly clear. The code above is never actually evaluated though; insteadit is parsed and used to build up C code for the model.

Note that this includes initial conditions; all odin models includespecifications for initial conditions because the ordering of thevariables is arbitrary and may be re-ordered.

This generates an object that can be used to integrate the set ofdifferential equations, by default starting at the initial conditionsspecified above (though custom initial conditions can be given). Theequations are translated into C, compiled, loaded, and bundled into anobject.lorenz here is a function that generates aninstance of the model.

mod<-lorenz()t<-seq(0,100,length.out =50000)y<- mod$run(t)

For more complicated examples, check out anagestructured SIR model, and for more details see themain packagevignette

Limitations

Writing this has given me a much greater appreciation of thedifficulties of writing compiler error messages.

This does not attempt togenerally translate R into C(though very simple expressions are handled) but only a small subsetthat follows the stereotyped way that R+C ODE models tend to be written.It tries to do things like minimise the number of memory allocationswhile preventing leaks. The generated code is designed to bestraightforward to read, leaving any really funky optimisation to thecompiler.

Because this relies on code generation, and the approach is partlytextual, some oddities will appear in the generated code (things liken + 0). Over time I’ll remove the most egregious of these.It’s probable that there will be some unused variables, and unusedelements in the parameters struct.

Prior work

ODEs seem particularly suitable for code generation, perhaps becauseof the relative simplicity of the code. As such, there is a lot of priorwork in this area. Many of these tools are heavily tailored to suit aparticular domain.

In R:

• RxODE -focussed on pharmacokinetic models, but suitable in the same domain asmany odin models. Does not include support for delay equations,automatic arrays or discrete/stochastic systems and uses it’s ownsolvers rather than interfacing with existing ones. Notably it also usesR as the host language for the DSL rather than requiring the user towrite code in strings or in a custom language.
• rodeofocussed on biochemical reactions based around thePetersenmatrix. Creates code for use withdeSolve
• cOdecreates code for use withdeSolveandbvpSolve.Models are entered as vector of strings which resembles C or R code.Automatic generation of Jacobian matrices is supported.
• mrgsolveis focussed on models in quantitative pharmacology and systems biology.It bundles its own solvers, and uses it’s ownPKMODELlanguage (example).

In other languages:

• Paraisoin Haskell
• VFGENin C++
• pygomin Python for solving compartmental models

Installation

Install odin from CRAN with

install.packages("odin")

Alternatively, you can install a potentially more recent version ofodin from themrc-ide dratrepository

# install.packages("drat") # -- if you don't have drat installeddrat:::add("mrc-ide")install.packages("odin")

You will need a compiler to install dependencies for the package, andto build any models with odin. Windows users should installRtools. Seethe relevant section inR-adminfor advice. Be sure to select the “edit PATH” checkbox duringinstallation or the tools will not be found.

The functionodin::can_compile() will check if it isable to compile things, but by the time you install the package thatwill probably have been satisfied.

The development version of the package can be installed directly fromgithub if you prefer with:

devtools::install_github("mrc-ide/odin",upgrade =FALSE)

License

MIT © Imperial College of Science, Technology and Medicine