Make compiled FOMC and DFOP models work

- Many thanks @karlines for the hint in ccSolve#1 how to declare
  the time variable and for writing the ccSolve package in the first
  place
- Add tests for the compiled versions of the FOMC and DFOP models

3 files changed, 37 insertions, 6 deletions

diff --git a/R/mkinmod.R b/R/mkinmod.R
index e73c56c6..3fff0878 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -126,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

@@ -295,10 +295,9 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
       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)) {

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

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

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

+      model$compiled <- compile.ode(diffs.C, language = "C", parms = parms, declaration = "double time = *t;")

     }

   }

   # }}}

diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R
index 8eefb995..01cbaf00 100644
--- a/inst/unitTests/runit.mkinfit.R
+++ b/inst/unitTests/runit.mkinfit.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

@@ -112,6 +112,18 @@ test.FOCUS_2006_DFOP <- function()
   dev.B.DFOP <- abs(round(100 * ((median.B.DFOP - fit.B.DFOP.r)/median.B.DFOP), digits=1))

   # about 0.6% deviation for parameter f, the others are <= 0.1%

   checkIdentical(dev.B.DFOP < 1, rep(TRUE, length(dev.B.DFOP)))

+

+  # Check the compiled version of possible FOCUS_2006_B

+  if (require(ccSolve)) {

+    checkTrue(!is.null(DFOP$compiled))

+    fit.B.DFOP.compiled <- mkinfit(DFOP, FOCUS_2006_B, solution_type = "deSolve", use_compiled = TRUE, quiet=TRUE)

+

+    fit.B.DFOP.compiled.r <- as.numeric(c(fit.B.DFOP.compiled$bparms.optim, 

+                                          endpoints(fit.B.DFOP)$distimes[c("DT50", "DT90")]))

+    dev.B.DFOP.compiled <- abs(round(100 * ((median.B.DFOP - fit.B.DFOP.compiled.r)/median.B.DFOP), digits=1))

+    # about 0.6% deviation for parameter f, the others are <= 0.1%

+    checkIdentical(dev.B.DFOP.compiled < 1, rep(TRUE, length(dev.B.DFOP)))

+  }

 } # }}}

 

 # Test HS model to a relative tolerance of 1% excluding Mathematica values {{{

diff --git a/inst/unitTests/runit.mkinpredict.R b/inst/unitTests/runit.mkinpredict.R
index d30bcfab..6635ea1f 100644
--- a/inst/unitTests/runit.mkinpredict.R
+++ b/inst/unitTests/runit.mkinpredict.R
@@ -34,6 +34,26 @@ test.SFO_solution_types <- function()
   checkEquals(SFO.analytical, SFO.eigen)

 } # }}}

 

+# Check solution types for FOMC {{{

+test.FOMC_solution_types <- function()

+{

+  ot = seq(0, 100, by = 1)

+  FOMC <- mkinmod(parent = list(type = "FOMC"))

+  FOMC.analytical <- round(subset(mkinpredict(FOMC, c(alpha = 1, beta = 10),

+        c(parent = 100), ot, solution_type = "analytical"), time == 100), digits=5)

+  FOMC.deSolve <- round(subset(mkinpredict(FOMC, c(alpha = 1, beta = 10),

+        c(parent = 100), ot, solution_type = "deSolve"), use_compiled = FALSE, time == 100), digits=5)

+  checkEquals(FOMC.analytical, FOMC.deSolve)

+

+  if (require(ccSolve)) {

+    checkTrue(!is.null(FOMC$compiled))

+    FOMC.deSolve.compiled <- round(subset(mkinpredict(FOMC, c(alpha = 1, beta = 10),

+          c(parent = 100), ot, solution_type = "deSolve"), time == 100), digits=5)

+    checkEquals(FOMC.analytical, FOMC.deSolve.compiled)

+  }

+

+} # }}}

+

 # Check model specification and solution types for SFO_SFO {{{

 # Relative Tolerance is 0.01%

 # Do not use time 0, as eigenvalue based solution does not give 0 at time 0 for metabolites