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

2 files changed, 47 insertions, 9 deletions

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,