Analytical SFO_SFO about as fast as deSolve compiled

2 files changed, 52 insertions, 19 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 40efb3a3..fb419a98 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -8,8 +8,16 @@
 #' SFO_SFO <- mkinmod(
 #'   parent = mkinsub("SFO", "m1"),
 #'   m1 = mkinsub("SFO"))
-#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
-
+#' FOCUS_D <- subset(FOCUS_2006_D, value != 0) # to avoid warnings
+#' fit_1 <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE)
+#' fit_2 <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE)
+#' \dontrun{
+#' if (require(rbenchmark))
+#'   benchmark(
+#'     analytical = mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
+#'     deSolve = mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
+#'     replications = 1)
+#' }
 create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
   use_of_ff <- match.arg(use_of_ff)
@@ -21,10 +29,13 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
   n <- character(0)
 
   parent <- obs_vars[1]
+  parent_type <- spec[[1]]$type
+
+  supported <- TRUE # This may be modified below
 
-  n[1] <- paste0(parent, " = ", spec[[1]]$type, ".solution(outtimes, odeini['", parent, 
-    if (spec[[1]]$type == "SFORB") "_free", "'], ",
-    switch(spec[[1]]$type,
+  n[1] <- paste0(parent, " = ", parent_type, ".solution(outtimes, odeini['", parent,
+    if (parent_type == "SFORB") "_free", "'], ",
+    switch(parent_type,
       SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"),
       FOMC = "alpha, beta)",
       IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent, ")"),
@@ -35,20 +46,42 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     )
   )
 
-  all_n <- paste(n, collapse = ",\n")
+  if (length(obs_vars) >= 2) {
+    supported <- FALSE # except for the following cases
+    n1 <- obs_vars[1]
+    n2 <- obs_vars[2]
+    n10 <- paste0("odeini['", parent, "']")
+    n20 <- paste0("odeini['", n2, "']")
 
-  f_body <- paste0("{\n",
-    "out <- with(as.list(odeparms), {\n",
-    "data.frame(\n",
-      "time = outtimes,\n",
-      all_n, "\n",
-    ")})\n",
-    "return(out)\n}\n"
-  )
+    if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2, spec[[2]]$type == "SFO")) {
+      supported <- TRUE
+      k1 <- paste0("k_", n1)
+      k2 <- paste0("k_", n2)
+      f12 <- paste0("f_", n1, "_to_", n2)
+      if (parent_type == "SFO") {
+        n[2] <- paste0(n2, " = (((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*outtimes)+",
+          f12, "*", k1, "*", n10, "*exp(-", k1, "*outtimes))/(", k2, "-", k1, ")")
+      }
+    }
+  }
 
-  deg_func <- function(odeini, odeparms, outtimes) {}
+  if (supported) {
+    all_n <- paste(n, collapse = ",\n")
+
+    f_body <- paste0("{\n",
+      "out <- with(as.list(odeparms), {\n",
+      "data.frame(\n",
+        "time = outtimes,\n",
+        all_n, ")\n",
+      "})\n",
+      "return(out)\n}\n"
+    )
 
-  body(deg_func) <- parse(text = f_body)
+    deg_func <- function(odeini, odeparms, outtimes) {}
 
-  return(deg_func)
+    body(deg_func) <- parse(text = f_body)
+    return(deg_func)
+  } else {
+    return(NULL)
+  }
 }
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 54dd75c2..1ddb6f36 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -459,8 +459,8 @@ mkinfit <- function(mkinmod, observed,
   # Prefer deSolve over eigen if a compiled model is present and use_compiled
   # is not set to FALSE
   solution_type = match.arg(solution_type)
-  if (solution_type == "analytical" && length(mkinmod$spec) > 1)
-     stop("Analytical solution not implemented for models with metabolites.")
+  if (solution_type == "analytical" && !is.function(mkinmod$deg_func))
+     stop("Analytical solution not implemented for this model.")
   if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
      stop("Eigenvalue based solution not possible, coefficient matrix not present.")
   if (solution_type == "auto") {