Compile differential equation models with ccSolve package

If the ccSolve package is available, and time is not in the right hand
side of the equations (i.e. if only SFO and SFORB models are used),
the differential equation model is compiled from auto-generated C code.

Currently, one test (FOCUS 2006 D SFO_SFO) fails

5 files changed, 56 insertions, 14 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 03976d1e..0eaf6ebf 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: mkin
 Type: Package
 Title: Routines for Fitting Kinetic Models with One or More State
   Variables to Chemical Degradation Data
-Version: 0.9-35
+Version: 0.9-36
 Date: 2015-02-20
 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), 
                     email = "jranke@uni-bremen.de"),
@@ -13,11 +13,13 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006).
   Includes a function for conveniently defining differential equation models,
   model solution based on eigenvalues if possible or using numerical solvers
   and a choice of the optimisation methods made available by the FME package.
+  If the ccSolve package and a C compiler is available, differential equation
+  models are solved as compiled C functions.
   Please note that no warranty is implied for correctness of results or fitness
   for a particular purpose.
 Depends: minpack.lm, rootSolve
 Imports: FME, deSolve
-Suggests: knitr, RUnit
+Suggests: knitr, RUnit, ccSolve
 License: GPL
 LazyLoad: yes
 LazyData: yes
diff --git a/NAMESPACE b/NAMESPACE
index 4774475e..92f0f271 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,5 +9,6 @@ import(
   FME,
   deSolve,
   minpack.lm,
-  rootSolve
+  rootSolve,
+  ccSolve
 )
diff --git a/R/mkinmod.R b/R/mkinmod.R
index d04f811f..37cd02b4 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -36,6 +36,9 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
   # differential equations (if supported), parameter names and a mapping from

   # model variables to observed variables. If possible, a matrix representation

   # of the differential equations is included

+  # Compiling the functions from the C code generated below using the ccSolve package

+  # only works if the implicit assumption about differential equations specified below

+  # is satisfied

   parms <- vector()

   # }}}

 

@@ -269,6 +272,37 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
     model$coefmat <- m

   }#}}}

 

+  # Create a function compiled from C code if possible #{{{

+  if (require(ccSolve)) {

+    diffs.C <- paste(diffs, collapse = ";\n")

+    diffs.C <- paste0(diffs.C, ";")

+    for (i in seq_along(diffs)) {

+      obs_var <- names(diffs)[i]

+

+      # Replace d_... terms by f[i-1]

+      # First line

+      pattern <- paste0("^d_", obs_var)

+      replacement <- paste0("\nf[", i - 1, "]")

+      diffs.C <- gsub(pattern, replacement, diffs.C)

+      # Other lines

+      pattern <- paste0("\\nd_", obs_var)

+      replacement <- paste0("\nf[", i - 1, "]")

+      diffs.C <- gsub(pattern, replacement, diffs.C)

+

+      # Replace names of observed variables by y[i],

+      # making the implicit assumption that the observed variables only occur after "* "

+      pattern <- paste0("\\* ", obs_var)

+      replacement <- paste0("* y[", i - 1, "]")

+      diffs.C <- gsub(pattern, replacement, diffs.C)

+    }

+    # Unfortunately, the models with time occuring in the diffs do not compile

+    if (sum(sapply(spec, function(x) x$type %in% c("SFO", "SFORB"))) == length(spec)) {

+      message("Compiling differential equation model from auto-generated C code...")

+      model$compiled <- compile.ode(diffs.C, language = "C", parms = parms)

+    }

+  }

+  # }}}

+

   class(model) <- "mkinmod"

   return(model)

 }

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index da013d50..bed82441 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -87,17 +87,21 @@ mkinpredict <- function(mkinmod, odeparms, odeini,
     names(out) <- c("time", mod_vars)
   } 
   if (solution_type == "deSolve") {
-    mkindiff <- function(t, state, parms) {
+    if (!is.null(mkinmod$compiled)) {
+       mkindiff <- mkinmod$compiled
+     } else {
+       mkindiff <- function(t, state, parms) {
 
-      time <- t
-      diffs <- vector()
-      for (box in names(mkinmod$diffs))
-      {
-        diffname <- paste("d", box, sep="_")      
-        diffs[diffname] <- with(as.list(c(time, state, parms)),
-          eval(parse(text=mkinmod$diffs[[box]])))
+        time <- t
+        diffs <- vector()
+        for (box in names(mkinmod$diffs))
+        {
+          diffname <- paste("d", box, sep="_")      
+          diffs[diffname] <- with(as.list(c(time, state, parms)),
+            eval(parse(text=mkinmod$diffs[[box]])))
+        }
+        return(list(c(diffs)))
       }
-      return(list(c(diffs)))
     } 
     out <- ode(
       y = odeini,
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 9c82f5ff..466f587c 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -235,7 +235,8 @@ SFO_SFO <- mkinmod(
   parent = list(type = "SFO", to = "m1", sink = TRUE),
   m1 = list(type = "SFO"))
 # Fit the model to the FOCUS example dataset D using defaults
-fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
+fit <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen")
+fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve")
 summary(fit)
 
 # Use stepwise fitting, using optimised parameters from parent only fit, FOMC
@@ -258,10 +259,10 @@ SFORB_SFO <- mkinmod(
   m1 = list(type = "SFO"))
 # Fit the model to the FOCUS example dataset D using defaults
 fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D)
+fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve")
 # Use starting parameters from parent only SFORB fit (not really needed in this case)
 fit.SFORB = mkinfit(SFORB, FOCUS_2006_D)
 fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode)
-}
 
 # Weighted fits, including IRLS
 SFO_SFO.ff <- mkinmod(parent = list(type = "SFO", to = "m1"),