Use odeintr instead of ccSolve for compiling modelsodeintr

8 files changed, 115 insertions, 79 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..ce380243 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -26,9 +26,8 @@ mkinfit <- function(mkinmod, observed,
   state.ini = "auto", 

   fixed_parms = NULL,

   fixed_initials = names(mkinmod$diffs)[-1],

-  solution_type = c("auto", "analytical", "eigen", "deSolve"),

+  solution_type = c("auto", "analytical", "eigen", "deSolve", "odeintr"),

   method.ode = "lsoda",

-  use_compiled = "auto",

   method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"),

   maxit.modFit = "auto",

   control.modFit = list(),

@@ -208,8 +207,7 @@ mkinfit <- function(mkinmod, observed,
   # Decide if the solution of the model can be based on a simple analytical

   # formula, the spectral decomposition of the matrix (fundamental system)

   # or a numeric ode solver from the deSolve package

-  # Prefer deSolve over eigen if a compiled model is present and use_compiled

-  # is not set to FALSE

+  # Prefer odeintr if a compiled model is present

   solution_type = match.arg(solution_type)

   if (solution_type == "analytical" && length(mkinmod$spec) > 1)

      stop("Analytical solution not implemented for models with metabolites.")

@@ -219,8 +217,11 @@ mkinfit <- function(mkinmod, observed,
     if (length(mkinmod$spec) == 1) {

       solution_type = "analytical"

     } else {

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

-        solution_type = "deSolve"

+      if (!is.null(mkinmod$e$m)) {

+        solution_type = "odeintr"

+        if (!missing(atol) | !missing(rtol)) {

+          message("When fitting with solution type odeintr, atol and rtol can only be set in the call to mkinmod()")

+        }

       } else {

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

           solution_type = "eigen"

@@ -269,9 +270,8 @@ mkinfit <- function(mkinmod, observed,
     out <- mkinpredict(mkinmod, parms, 

                        odeini, outtimes, 

                        solution_type = solution_type, 

-                       use_compiled = use_compiled,

                        method.ode = method.ode,

-                       atol = atol, rtol = rtol, ...)

+                       atol.deSolve = atol, rtol.deSolve = rtol, ...)

 

     assign("out_predicted", out, inherits=TRUE)

 

@@ -289,9 +289,8 @@ mkinfit <- function(mkinmod, observed,
         out_plot <- mkinpredict(mkinmod, parms,

                                 odeini, outtimes_plot, 

                                 solution_type = solution_type, 

-                                use_compiled = use_compiled,

                                 method.ode = method.ode,

-                                atol = atol, rtol = rtol, ...)

+                                atol.deSolve = atol, rtol.deSolve = rtol, ...)

 

         plot(0, type="n", 

           xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)),

@@ -319,8 +318,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..4991be22 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -16,12 +16,17 @@
 # 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, quiet = FALSE)

+mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE,

+                    odeintr_compile = c("auto", "no", "yes"),

+                    odeintr_method = "rk5_i",

+                    odeintr_atol = 1e-8, odeintr_rtol = 1e-10)

 {

   if (is.null(speclist)) spec <- list(...)

   else spec <- speclist

   obs_vars <- names(spec)

 

+  odeintr_compile <- match.arg(odeintr_compile)

+

   # Check if any of the names of the observed variables contains any other

   for (obs_var in obs_vars) {

     if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")

@@ -99,7 +104,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 +119,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)

@@ -127,17 +132,18 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
     } #}}}

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

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

+      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

       # From p. 57 of the FOCUS kinetics report

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

+      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 +278,58 @@ 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]

+  # Try to create a function compiled from C code if more than one observed variable and g++ can be called #{{{

+  if (length(obs_vars) > 1) {

+    try(gcc_version <- system("g++ -v 2>&1", intern = TRUE))

+    if (!inherits(gcc_version, "try-error") & odeintr_compile != "no") {

 

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

+      if ((odeintr_compile == "auto" & mat == FALSE) | odeintr_compile == "yes") {

+        # HS

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

 

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

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

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

+

+        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(x[", i - 1, "], N_", state_var, ")"), 

+                              diffs.C, fixed = TRUE)

+            }

+          }

+

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

+          # First line

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

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

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

+          # Other lines

+          pattern <- paste0("\\nd_", state_var)

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

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

+

+          # Replace names of observed variables by x[i],

+          # making the implicit assumption that the observed variables only occur after "* "

+          pattern <- paste0("\\* ", state_var)

+          replacement <- paste0("* x[", i - 1, "]")

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

+

+          # Handle time, which needs to be enclosed by spaces (except for HS which is handled separately)

+          diffs.C <- gsub(" time ", " t ", diffs.C)

+        }

+        if (!quiet) message("Compiling differential equation model from auto-generated C++ code...")

+        model$e <- new.env()

+        model$e$code <- try(compile_sys("m", diffs.C, pars = parms, method = odeintr_method, 

+                                        atol = odeintr_atol, rtol = odeintr_rtol, env = model$e))

+

+        if (inherits(model$e$code, "try-error")) {

+          if (!quiet) message("Compilation failed, model solution will be slower than necessary")

+        }

+      }

     }

   }

   # }}}

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 3121f1d7..d91e87df 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -18,11 +18,12 @@
 # this program. If not, see <http://www.gnu.org/licenses/>
 
 mkinpredict <- function(mkinmod, odeparms, odeini, 
-			outtimes, solution_type = "deSolve", 
-      use_compiled = "auto",
-			method.ode = "lsoda", atol = 1e-8, rtol = 1e-10,
+			outtimes, solution_type = c("deSolve", "analytical", "eigen", "odeintr"),
+			method.ode = "lsoda", atol.deSolve = 1e-8, rtol.deSolve = 1e-10,
 			map_output = TRUE, ...) {
 
+  solution_type = match.arg(solution_type)
+
   # Get the names of the state variables in the model
   mod_vars <- names(mkinmod$diffs)
 
@@ -54,8 +55,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,30 +89,26 @@ 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) {
+     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]])))
-        }
-        return(list(c(diffs)))
+      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)))
+    }
     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,
+      atol = atol.deSolve,
+      rtol = rtol.deSolve,
       ...
     )
     if (sum(is.na(out)) > 0) {
@@ -119,6 +116,17 @@ mkinpredict <- function(mkinmod, odeparms, odeini,
 	   "NaN values occurred in output from ode()")
       }
   }
+  if (solution_type == "odeintr") {
+    if (is.null(mkinmod$e$m)) stop("Method odeintr can not be used as no compiled version of the model is available")
+    odeparms_argstring = ""
+    for (parname in names(odeparms)) {
+      odeparms_argstring = paste0(odeparms_argstring, parname, " = ", odeparms[parname], ", ")
+    }
+    odeparms_argstring = gsub(", $", "", odeparms_argstring)
+    with(as.list(odeparms_argstring), eval(parse(text=paste0("mkinmod$e$m_set_params(", odeparms_argstring, ")"))))
+    out <- mkinmod$e$m_at(init = odeini, times = outtimes)
+    names(out) <- c("time", names(mkinmod$diffs))
+  }
   if (map_output) {
     # Output transformation for models with unobserved compartments like SFORB
     out_mapped <- data.frame(time = out[,"time"])
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index 31746fb8..1523be7d 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -55,7 +55,7 @@ plot.mkinfit <- function(x, fit = x,
   odeparms <- parms.all[odenames]
 
   out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, 
-          solution_type = solution_type, atol = fit$atol, rtol = fit$rtol)
+          solution_type = solution_type, atol.deSolve = fit$atol, rtol.deSolve = fit$rtol)
 
   # Set up the plot if not to be added to an existing plot
   if (add == FALSE) {
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