Low-level generation of compiled models

As it is unclear if and when ccSolve will be published on CRAN,
the generation, compilation and use of the C version of the
system of differential equations was developed for mkin, inspired and
guided by the code from the ccSolve package. Many thanks again to
Karline Soetaert for all of her work on this and other R packages.

Now all model types, including the Hockey-Stick model for the parent
compund and the IORE model for parent and/or metabolites can be compiled.

7 files changed, 104 insertions, 62 deletions

diff --git a/R/IORE.solution.R b/R/IORE.solution.R
index 9546ce56..5405be96 100644
--- a/R/IORE.solution.R
+++ b/R/IORE.solution.R
@@ -1,4 +1,4 @@
-IORE.solution <- function(t, parent.0, k.iore, N)

+IORE.solution <- function(t, parent.0, k__iore, N)

 {

-	parent = (parent.0^(1 - N) - (1 - N) * k.iore * t)^(1/(1 - N))

+	parent = (parent.0^(1 - N) - (1 - N) * k__iore * t)^(1/(1 - N))

 }

diff --git a/R/endpoints.R b/R/endpoints.R
index ae9aa3ec..001595bc 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -50,7 +50,7 @@ endpoints <- function(fit) {
       ep$distimes[obs_var, c("DT50back")] = DT50_back

     }

     if (type == "IORE") {

-      k_names = grep(paste("k.iore", obs_var, sep="_"), names(parms.all), value=TRUE)

+      k_names = grep(paste("k__iore", obs_var, sep="_"), names(parms.all), value=TRUE)

       k_tot = sum(parms.all[k_names])

       # From the NAFTA kinetics guidance, p. 5

       n = parms.all[paste("N", obs_var, sep = "_")]

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index b28235a6..6103dd1d 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -61,8 +61,8 @@ mkinerrmin <- function(fit, alpha = 0.05)
     n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
     n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), 
                                          names(parms.optim)))
-    n.k.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim)))
-    n.k.iore.optim <- n.k.iore.optim + length(grep(paste("^log_k.iore", obs_var, 
+    n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim)))
+    n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, 
 							 sep = "_"),
 						   names(parms.optim)))
 
@@ -82,7 +82,7 @@ mkinerrmin <- function(fit, alpha = 0.05)
       }
     }
 
-    n.optim <- sum(n.initials.optim, n.k.optim, n.k.iore.optim, n.N.optim, n.ff.optim)
+    n.optim <- sum(n.initials.optim, n.k.optim, n.k__iore.optim, n.N.optim, n.ff.optim)
 
     # FOMC, DFOP and HS parameters are only counted if we are looking at the
     # first variable in the model which is always the source variable
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 5118519a..7f9a55d9 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -219,7 +219,7 @@ mkinfit <- function(mkinmod, observed,
     if (length(mkinmod$spec) == 1) {

       solution_type = "analytical"

     } else {

-      if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) {

+      if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {

         solution_type = "deSolve"

       } else {

         if (is.matrix(mkinmod$coefmat)) {

@@ -319,8 +319,8 @@ mkinfit <- function(mkinmod, observed,
   if (!transform_rates) {

     index_k <- grep("^k_", names(lower))

     lower[index_k] <- 0

-    index_k.iore <- grep("^k.iore_", names(lower))

-    lower[index_k.iore] <- 0

+    index_k__iore <- grep("^k__iore_", names(lower))

+    lower[index_k__iore] <- 0

     other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))

     lower[other_rate_parms] <- 0

   }

diff --git a/R/mkinmod.R b/R/mkinmod.R
index e865fa27..a0307003 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -36,8 +36,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
   # 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

+  # Compiling the functions from the C code generated below only works if the

+  # implicit assumption about differential equations specified below

   # is satisfied

   parms <- vector()

   # }}}

@@ -99,7 +99,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
           # If sink is required, add first-order/IORE sink term

           k_compound_sink <- paste("k", box_1, "sink", sep = "_")

           if(spec[[varname]]$type == "IORE") {

-            k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_")

+            k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")

           }

           parms <- c(parms, k_compound_sink)

           decline_term <- paste(k_compound_sink, "*", box_1)

@@ -114,7 +114,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
       } else {

         k_compound <- paste("k", box_1, sep = "_")

         if(spec[[varname]]$type == "IORE") {

-          k_compound <- paste("k.iore", box_1, sep = "_")

+          k_compound <- paste("k__iore", box_1, sep = "_")

         }

         parms <- c(parms, k_compound)

         decline_term <- paste(k_compound, "*", box_1)

@@ -135,9 +135,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
       decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)

       parms <- c(parms, "k1", "k2", "g")

     } #}}}

+    HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C

     if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term

       # From p. 55 of the FOCUS kinetics report

-      decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)

+      decline_term <- paste(HS_decline, "*", box_1)

       parms <- c(parms, "k1", "k2", "tb")

     } #}}}

     # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{

@@ -272,35 +273,65 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
     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]

+  # Create a function compiled from C code if more than one observed variable and gcc is available #{{{

+  if (length(obs_vars) > 1) {

+    if (Sys.which("gcc") != "") {

 

-      # 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)

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

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

+

+      # HS

+      diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)

+

+      for (i in seq_along(diffs)) {

+        state_var <- names(diffs)[i]

+

+        # IORE

+        if (state_var %in% obs_vars) {

+          if (spec[[state_var]]$type == "IORE") {

+            diffs.C <- gsub(paste0(state_var, "^N_", state_var),

+                            paste0("pow(y[", i - 1, "], N_", state_var, ")"),

+                            diffs.C, fixed = TRUE)

+          }

+        }

+

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

+        # First line

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

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

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

+        # Other lines

+        pattern <- paste0("\\nd_", state_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("\\* ", state_var)

+        replacement <- paste0("* y[", 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;")

+

+      npar <- length(parms)

+

+      initpar_code <- paste0(

+        "static double parms [", npar, "];\n",

+        paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),

+        "\n",

+        "void initpar(void (* odeparms)(int *, double *)) {\n",

+        "    int N = ", npar, ";\n",

+        "    odeparms(&N, parms);\n",

+        "}\n\n")

+

+      derivs_code <- paste0("double time = *t;\n", diffs.C)

+      derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",

+                              f = "numeric", rpar = "numeric", ipar = "integer")

+

+      model$cf <- cfunction(list(func = derivs_sig), derivs_code,

+                            otherdefs = initpar_code,

+                            convention = ".C", language = "C")

     }

   }

   # }}}

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 3121f1d7..5f11c35a 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -54,8 +54,8 @@ mkinpredict <- function(mkinmod, odeparms, odeini,
       IORE = IORE.solution(outtimes,
           evalparse(parent.name),
           ifelse(mkinmod$use_of_ff == "min", 
-	    evalparse(paste("k.iore", parent.name, "sink", sep="_")),
-	    evalparse(paste("k.iore", parent.name, sep="_"))),
+	    evalparse(paste("k__iore", parent.name, "sink", sep="_")),
+	    evalparse(paste("k__iore", parent.name, sep="_"))),
             evalparse("N_parent")),
       DFOP = DFOP.solution(outtimes,
           evalparse(parent.name),
@@ -88,10 +88,21 @@ mkinpredict <- function(mkinmod, odeparms, odeini,
     names(out) <- c("time", mod_vars)
   } 
   if (solution_type == "deSolve") {
-    if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) {
-       mkindiff <- mkinmod$compiled
-     } else {
-       mkindiff <- function(t, state, parms) {
+    if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {
+      out <- ode(
+        y = odeini,
+        times = outtimes,
+        func = "func",
+        initfunc = "initpar",
+        dllname = getDynLib(mkinmod$cf)[["name"]],
+        parms = odeparms[mkinmod$parms], # Order matters when using compiled models
+        method = method.ode,
+        atol = atol,
+        rtol = rtol,
+        ...
+      )
+    } else {
+      mkindiff <- function(t, state, parms) {
 
         time <- t
         diffs <- vector()
@@ -103,17 +114,17 @@ mkinpredict <- function(mkinmod, odeparms, odeini,
         }
         return(list(c(diffs)))
       }
+      out <- ode(
+        y = odeini,
+        times = outtimes,
+        func = mkindiff, 
+        parms = odeparms,
+        method = method.ode,
+        atol = atol,
+        rtol = rtol,
+        ...
+      )
     } 
-    out <- ode(
-      y = odeini,
-      times = outtimes,
-      func = mkindiff, 
-      parms = odeparms[mkinmod$parms], # Order matters when using compiled models
-      method = method.ode,
-      atol = atol,
-      rtol = rtol,
-      ...
-    )
     if (sum(is.na(out)) > 0) {
       stop("Differential equations were not integrated for all output times because\n",
 	   "NaN values occurred in output from ode()")
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index aa70a0b3..0946ff14 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -33,8 +33,8 @@ transform_odeparms <- function(parms, mkinmod,
 

   # Log transformation for rate constants if requested

   k <- parms[grep("^k_", names(parms))]

-  k.iore <- parms[grep("^k.iore_", names(parms))]

-  k <- c(k, k.iore)

+  k__iore <- parms[grep("^k__iore_", names(parms))]

+  k <- c(k, k__iore)

   if (length(k) > 0) {

     if(transform_rates) {

       transparms[paste0("log_", names(k))] <- log(k)

@@ -113,8 +113,8 @@ backtransform_odeparms <- function(transparms, mkinmod,
   # Exponential transformation for rate constants

   if(transform_rates) {

     trans_k <- transparms[grep("^log_k_", names(transparms))]

-    trans_k.iore <- transparms[grep("^log_k.iore_", names(transparms))]

-    trans_k = c(trans_k, trans_k.iore)

+    trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]

+    trans_k = c(trans_k, trans_k__iore)

     if (length(trans_k) > 0) {

       k_names <- gsub("^log_k", "k", names(trans_k))

       parms[k_names] <- exp(trans_k)

@@ -122,8 +122,8 @@ backtransform_odeparms <- function(transparms, mkinmod,
   } else {

     trans_k <- transparms[grep("^k_", names(transparms))]

     parms[names(trans_k)] <- trans_k

-    trans_k.iore <- transparms[grep("^k.iore_", names(transparms))]

-    parms[names(trans_k.iore)] <- trans_k.iore

+    trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]

+    parms[names(trans_k__iore)] <- trans_k__iore

   }

 

   # Do not transform exponents in IORE models