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.

17 files changed, 137 insertions, 89 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 59b653f4..b5045e4e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -3,7 +3,7 @@ Type: Package
 Title: Routines for Fitting Kinetic Models with One or More State
   Variables to Chemical Degradation Data
 Version: 0.9-36
-Date: 2015-05-18
+Date: 2015-06-20
 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), 
                     email = "jranke@uni-bremen.de"),
              person("Katrin", "Lindenberger", role = "ctb"),
@@ -13,14 +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 (https://github.com/karlines/ccSolve), 
-  and a C compiler (on windows: 'Rtools') are installed, differential
+  If a C compiler (on windows: 'Rtools') are installed, differential
   equation models are solved using compiled C functions.  
   Please note that no warranty is implied for correctness of results or fitness
   for a particular purpose.
-Depends: minpack.lm, rootSolve
+Depends: minpack.lm, rootSolve, inline
 Imports: FME, deSolve
-Suggests: knitr, testthat, microbenchmark, ccSolve
+Suggests: knitr, testthat, microbenchmark
 License: GPL
 LazyLoad: yes
 LazyData: yes
diff --git a/GNUmakefile b/GNUmakefile
index 9e2ee80c..7135d6ee 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -26,7 +26,11 @@ pkgfiles = NEWS \
 	   tests/* \
 	   tests/testthat* \
 	   TODO \
-	   vignettes/*
+	   vignettes/*.Rmd \
+	   vignettes/*.Rnw
+
+
+
 
 all: check clean
 
diff --git a/NAMESPACE b/NAMESPACE
index 49be428e..decb895d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,5 +9,6 @@ import(
   FME,
   deSolve,
   minpack.lm,
-  rootSolve
+  rootSolve,
+  inline
 )
diff --git a/NEWS.md b/NEWS.md
index 32ca3e45..2f58417f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 ## MAJOR CHANGES
 
-- If the ccSolve package is installed, use a version of the differential equation model compiled from C code, which is a huge performance boost for models where only the deSolve method works. Compiled models using the Hockey Stick model (HS) are not supported (yet).
+- If a compiler (gcc) is installed, use a version of the differential equation model compiled from C code, which is a huge performance boost for models where only the deSolve method works.
 - `mkinmod()`: Create a list component $compiled in the list returned by mkinmod, if a version can be compiled from autogenerated C code (see above).
 - `mkinfit()`: Set the default `solution_type` to `deSolve` when a compiled version of the model is present, except if an analytical solution is possible.
 
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

diff --git a/README.Rmd b/README.Rmd
index 7c9cf5d2..c5df2aef 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -106,10 +106,12 @@ documentation or the package vignettes referenced from the
   These have decreasing efficiency, and are automatically chosen 
   by default. 
 * As of mkin 0.9-36, model solution for models with more than one observed
-  variable is based on the
-  [`ccSolve`](https://github.com/karlines/ccSolve) package, if installed.
-  This is even faster than eigenvalue based solution, at least in the example
-  shown in the [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html)
+  variable is based on the inline package.  This is even faster than eigenvalue
+  based solution, at least in the example shown in the
+  [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html).
+  The autogeneration of C code was
+  inspired by the [`ccSolve`](https://github.com/karlines/ccSolve) package. Thanks
+  to Karline Soetaert for her work on that.
 * Model optimisation with 
   [`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html)
   internally using the `modFit` function from the `FME` package,
diff --git a/README.md b/README.md
index 98daa6f9..37af4f31 100644
--- a/README.md
+++ b/README.md
@@ -120,10 +120,12 @@ documentation or the package vignettes referenced from the
   These have decreasing efficiency, and are automatically chosen 
   by default. 
 * As of mkin 0.9-36, model solution for models with more than one observed
-  variable is based on the
-  [`ccSolve`](https://github.com/karlines/ccSolve) package, if installed.
-  This is even faster than eigenvalue based solution, at least in the example
-  shown in the [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html)
+  variable is based on the inline package.  This is even faster than eigenvalue
+  based solution, at least in the example shown in the
+  [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html).
+  The autogeneration of C code was
+  inspired by the [`ccSolve`](https://github.com/karlines/ccSolve) package. Thanks
+  to Karline Soetaert for her work on that.
 * Model optimisation with 
   [`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html)
   internally using the `modFit` function from the `FME` package,
diff --git a/man/IORE.solution.Rd b/man/IORE.solution.Rd
index 3d827bdc..bb377d3f 100644
--- a/man/IORE.solution.Rd
+++ b/man/IORE.solution.Rd
@@ -7,12 +7,12 @@
   a concentration dependent rate constant.

 }

 \usage{

-  IORE.solution(t, parent.0, k.iore, N)

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

 }

 \arguments{

   \item{t}{ Time. }

   \item{parent.0}{ Starting value for the response variable at time zero. }

-  \item{k.iore}{ Rate constant. Note that this depends on the concentration units used. }

+  \item{k__iore}{ Rate constant. Note that this depends on the concentration units used. }

   \item{N}{ Exponent describing the nonlinearity of the rate equation } 

 }

 \note{

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index b4923ddd..bf47879c 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -247,7 +247,7 @@ print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D,
 coef(fit)
 endpoints(fit)
 \dontrun{
-# deSolve is slower when ccSolve is not installed and set up
+# deSolve is slower when the gcc compiler is not available
 print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, 
                            solution_type = "deSolve")))
 coef(fit.deSolve)
diff --git a/vignettes/FOCUS_D.Rmd b/vignettes/FOCUS_D.Rmd
index 902d3d24..8ae73c16 100644
--- a/vignettes/FOCUS_D.Rmd
+++ b/vignettes/FOCUS_D.Rmd
@@ -27,7 +27,7 @@ kinetics (SFO) to one metabolite named m1, which also degrades with SFO kinetics
 

 The call to mkinmod returns a degradation model. The differential equations represented in 

 R code can be found in the character vector `$diffs` of the `mkinmod` object. If

-the `ccSolve` package is installed and functional, the differential equation model

+the gcc compiler is installed and functional, the differential equation model

 will be compiled from auto-generated C code.

 

 

diff --git a/vignettes/compiled_models.Rmd b/vignettes/compiled_models.Rmd
index bac284c5..5d83cf3a 100644
--- a/vignettes/compiled_models.Rmd
+++ b/vignettes/compiled_models.Rmd
@@ -15,16 +15,14 @@ output:
 ```{r, include = FALSE}

 library(knitr)

 opts_chunk$set(tidy = FALSE, cache = TRUE)

-if (!require("ccSolve")) 

-  message("Please install the ccSolve package for this vignette to produce sensible output")

-

 ```

 

 # Benchmark for a model that can also be solved with Eigenvalues

 

 This evaluation is taken from the example section of mkinfit. When using an mkin version

-greater than 0.9-36 and the ccSolve package is installed and functional, you will get a

-message that the model is being compiled when defining a model using mkinmod.

+equal to or greater than 0.9-36 and a compiler (gcc) is installed, you will see

+a message that the model is being compiled from autogenerated C code when

+defining a model using mkinmod.

 

 ```{r create_SFO_SFO}

 library("mkin")

@@ -83,5 +81,5 @@ smb.2["median"]/smb.2["deSolve, compiled", "median"]
 

 ```

 

-Here we get a performance benefit of more than a factor of 8 using the version

-of the differential equation model compiled from C code using the ccSolve package!

+Here we get a performance benefit of more than a factor of 10 using the version

+of the differential equation model compiled from C code using the inline package!