Merge branch 'compile_odes'

5 files changed, 57 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..4b0ad1a7 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"),
@@ -275,6 +276,7 @@ summary(f.w.mean)
 f.w.mean.irls <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, weight = "mean",
                          reweight.method = "obs")
 summary(f.w.mean.irls)
+}
 
 \dontrun{
 # Manual weighting