Change implementation of analytical solutions

Preparing for symbolic solutions for more than one compound

3 files changed, 104 insertions, 94 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 5c092612..2b7e71cb 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -185,13 +185,13 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
 #' # Fit the model to the FOCUS example dataset D using defaults
 #' print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D,
 #'                            solution_type = "eigen", quiet = TRUE)))
-#' coef(fit)
+#' parms(fit)
 #' endpoints(fit)
 #' \dontrun{
 #' # deSolve is slower when no C compiler (gcc) was available during model generation
 #' print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D,
 #'                            solution_type = "deSolve")))
-#' coef(fit.deSolve)
+#' parms(fit.deSolve)
 #' endpoints(fit.deSolve)
 #' }
 #'
@@ -926,6 +926,6 @@ mkinfit <- function(mkinmod, observed,
   fit$version <- as.character(utils::packageVersion("mkin"))
   fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")
 
-  class(fit) <- c("mkinfit", "modFit")
+  class(fit) <- c("mkinfit")
   return(fit)
 }
diff --git a/R/mkinmod.R b/R/mkinmod.R
index ca1402fd..099e1155 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -421,6 +421,46 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
   }
   # }}}
 
+  # Attach a degradation function if an analytical solution is available {{{
+  parent_type = spec[[1]]$type
+  parent_name = names(spec)[[1]]
+  if (length(spec) == 1) {
+    odeparm_map <- switch(parent_type,
+      SFO = c(
+        k = if(use_of_ff == "min") paste("k", parent_name, "sink", sep = "_")
+          else  paste("k", parent_name, sep = "_")),
+      FOMC = c(alpha = "alpha", beta = "beta"),
+      IORE = c(
+        k__iore = if(use_of_ff == "min") paste("k__iore", parent_name, "sink", sep = "_")
+          else paste("k__iore", parent_name, sep = "_"),
+        N = paste("N", parent_name, sep = "_")),
+      DFOP = c(k1 = "k1", k2 = "k2", g = "g"),
+      HS = c(k1 = "k1", k2 = "k2", tb = "tb"),
+      SFORB = c(
+        k_12 = paste("k", parent_name, "free_bound", sep = "_"),
+        k_21 = paste("k", parent_name, "bound_free", sep = "_"),
+        k_1output = paste("k", parent_name, "free_sink", sep = "_")),
+      logistic = c(kmax = "kmax", k0 = "k0", r = "r")
+    )
+    odeparm_rev_map <- names(odeparm_map)
+    names(odeparm_rev_map) <- odeparm_map
+
+    model$deg_func <- function(odeini, odeparms, outtimes) {
+      parent_func <- getFromNamespace(paste0(parent_type, ".solution"), "mkin")
+      odeparm_list <- as.list(odeparms)
+      names(odeparm_list) <- odeparm_rev_map[names(odeparm_list)]
+
+      values <- do.call(parent_func,
+        args = c(
+          list(t = outtimes, parent.0 = odeini[1]),
+          odeparm_list))
+      out <- data.frame(outtimes, values)
+      names(out) <- c("time", parent_name)
+      return(out)
+    }
+  }
+  # }}}
+
   class(model) <- "mkinmod"
   return(model)
 }
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 16ee7903..0f8e83bb 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -41,57 +41,63 @@
 #' @author Johannes Ranke
 #' @examples
 #'
-#'   SFO <- mkinmod(degradinol = mkinsub("SFO"))
-#'   # Compare solution types
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         solution_type = "analytical")
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         solution_type = "deSolve")
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         solution_type = "deSolve", use_compiled = FALSE)
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         solution_type = "eigen")
+#' SFO <- mkinmod(degradinol = mkinsub("SFO"))
+#' # Compare solution types
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       solution_type = "analytical")
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       solution_type = "deSolve")
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       solution_type = "deSolve", use_compiled = FALSE)
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       solution_type = "eigen")
 #'
+#' # Compare integration methods to analytical solution
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       solution_type = "analytical")[21,]
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       method = "lsoda")[21,]
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       method = "ode45")[21,]
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#'       method = "rk4")[21,]
+#' # rk4 is not as precise here
 #'
-#'   # Compare integration methods to analytical solution
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         solution_type = "analytical")[21,]
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         method = "lsoda")[21,]
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         method = "ode45")[21,]
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
-#'         method = "rk4")[21,]
-#'  # rk4 is not as precise here
+#' # The number of output times used to make a lot of difference until the
+#' # default for atol was adjusted
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
+#'       seq(0, 20, by = 0.1))[201,]
+#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
+#'       seq(0, 20, by = 0.01))[2001,]
 #'
-#'   # The number of output times used to make a lot of difference until the
-#'   # default for atol was adjusted
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
-#'         seq(0, 20, by = 0.1))[201,]
-#'   mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
-#'         seq(0, 20, by = 0.01))[2001,]
+#' # Check compiled model versions - they are faster than the eigenvalue based solutions!
+#' SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),
+#'                   m1 = list(type = "SFO"))
+#' if(require(rbenchmark)) {
+#'   benchmark(
+#'     eigen = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+#'       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+#'       solution_type = "eigen")[201,],
+#'     deSolve_compiled = mkinpredict(SFO_SFO, 
+#'       c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+#'       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+#'       solution_type = "deSolve")[201,],
+#'     deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+#'      c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+#'      solution_type = "deSolve", use_compiled = FALSE)[201,],
+#'     replications = 10)
+#' }
 #'
-#'   # Check compiled model versions - they are faster than the eigenvalue based solutions!
-#'   SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),
-#'                     m1 = list(type = "SFO"))
-#'   system.time(
-#'     print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
-#'                 c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-#'                 solution_type = "eigen")[201,]))
-#'   system.time(
-#'     print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
-#'                 c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-#'                 solution_type = "deSolve")[201,]))
-#'   system.time(
-#'     print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
-#'                 c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-#'                 solution_type = "deSolve", use_compiled = FALSE)[201,]))
+#' # Since mkin 0.9.49.11 we also have analytical solutions for some models, including SFO-SFO
+#' #    deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+#' #     c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+#' #     solution_type = "analytical", use_compiled = FALSE)[201,],
 #'
-#'   \dontrun{
-#'     # Predict from a fitted model
-#'     f <- mkinfit(SFO_SFO, FOCUS_2006_C)
-#'     head(mkinpredict(f))
-#'   }
+#' \dontrun{
+#'   # Predict from a fitted model
+#'   f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE)
+#'   head(mkinpredict(f))
+#' }
 #'
 #' @export
 mkinpredict <- function(x, odeparms, odeini,
@@ -124,54 +130,17 @@ mkinpredict.mkinmod <- function(x,
     odeini <- odeini[mod_vars]
   }
 
-  # Create function for evaluation of expressions with ode parameters and initial values
-  evalparse <- function(string)
-  {
-    eval(parse(text=string), as.list(c(odeparms, odeini)))
-  }
-
-  # Create a function calculating the differentials specified by the model
-  # if necessary
   if (solution_type == "analytical") {
-    parent.type = names(x$map[[1]])[1]
-    parent.name = names(x$diffs)[[1]]
-    o <- switch(parent.type,
-      SFO = SFO.solution(outtimes,
-        evalparse(parent.name),
-          ifelse(x$use_of_ff == "min",
-            evalparse(paste("k", parent.name, "sink", sep="_")),
-            evalparse(paste("k", parent.name, sep="_")))),
-      FOMC = FOMC.solution(outtimes,
-        evalparse(parent.name),
-        evalparse("alpha"), evalparse("beta")),
-      IORE = IORE.solution(outtimes,
-        evalparse(parent.name),
-        ifelse(x$use_of_ff == "min",
-          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),
-        evalparse("k1"), evalparse("k2"),
-        evalparse("g")),
-      HS = HS.solution(outtimes,
-        evalparse(parent.name),
-        evalparse("k1"), evalparse("k2"),
-        evalparse("tb")),
-      SFORB = SFORB.solution(outtimes,
-        evalparse(parent.name),
-        evalparse(paste("k", parent.name, "bound", sep="_")),
-        evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
-        evalparse(paste("k", parent.name, "sink", sep="_"))),
-      logistic = logistic.solution(outtimes,
-        evalparse(parent.name),
-        evalparse("kmax"), evalparse("k0"),
-        evalparse("r"))
-    )
-    out <- data.frame(outtimes, o)
-    names(out) <- c("time", sub("_free", "", parent.name))
+    out <- x$deg_func(odeini = odeini,
+      odeparms = odeparms, outtimes = outtimes)
   }
+
   if (solution_type == "eigen") {
+
+    evalparse <- function(string) {
+      eval(parse(text=string), as.list(c(odeparms, odeini)))
+    }
+
     coefmat.num <- matrix(sapply(as.vector(x$coefmat), evalparse),
       nrow = length(mod_vars))
     e <- eigen(coefmat.num)
@@ -184,6 +153,7 @@ mkinpredict.mkinmod <- function(x,
     out <- data.frame(outtimes, t(o))
     names(out) <- c("time", mod_vars)
   }
+
   if (solution_type == "deSolve") {
     if (!is.null(x$cf) & use_compiled[1] != FALSE) {
       out <- ode(