1 files changed, 40 insertions, 4 deletions
diff --git a/R/mkinmod.R b/R/mkinmod.R
index d04f811f..e865fa27 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2014 Johannes Ranke {{{
+# Copyright (C) 2010-2015 Johannes Ranke {{{
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -16,7 +16,7 @@
 # You should have received a copy of the GNU General Public License along with
 # this program. If not, see <http://www.gnu.org/licenses/> }}}
 
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
+mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
 {
   if (is.null(speclist)) spec <- list(...)
   else spec <- speclist
@@ -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()
   # }}}
 
@@ -123,8 +126,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
       }
     } #}}}
     if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
-      # From p. 53 of the FOCUS kinetics report
-      decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)
+      # From p. 53 of the FOCUS kinetics report, without the power function so it works in C
+      decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1)
       parms <- c(parms, "alpha", "beta")
     } #}}}
     if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
@@ -269,6 +272,39 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
     model$coefmat <- m
   }#}}}
 
+  # Create a function compiled from C code if more than one observed variable and ccSolve is available #{{{
+  if (length(obs_vars) > 1 & requireNamespace("ccSolve", quietly = TRUE)) {
+    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)
+    }
+    if (sum(sapply(spec, function(x) x$type %in% 
+                   c("SFO", "FOMC", "DFOP", "SFORB"))) == length(spec)) {
+      if (!quiet) message("Compiling differential equation model from auto-generated C code...")
+      model$compiled <- ccSolve::compile.ode(diffs.C, language = "C", 
+                                             parms = parms, 
+                                             declaration = "double time = *t;")
+    }
+  }
+  # }}}
+
   class(model) <- "mkinmod"
   return(model)
 }