Another overhaul of analytical solutions

Still in preparation for analytical solutions of coupled models

7 files changed, 74 insertions, 116 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
new file mode 100644
index 00000000..40efb3a3
--- /dev/null
+++ b/R/create_deg_func.R
@@ -0,0 +1,54 @@
+#' Create degradation functions for known analytical solutions
+#'
+#' @param spec List of model specifications as contained in mkinmod objects
+#' @param use_of_ff Minimum or maximum use of formation fractions
+#' @return Degradation function to be attached to mkinmod objects
+#' @examples
+#'
+#' SFO_SFO <- mkinmod(
+#'   parent = mkinsub("SFO", "m1"),
+#'   m1 = mkinsub("SFO"))
+#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
+
+create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
+
+  use_of_ff <- match.arg(use_of_ff)
+
+  min_ff <- use_of_ff == "min"
+
+  obs_vars <- names(spec)
+
+  n <- character(0)
+
+  parent <- obs_vars[1]
+
+  n[1] <- paste0(parent, " = ", spec[[1]]$type, ".solution(outtimes, odeini['", parent, 
+    if (spec[[1]]$type == "SFORB") "_free", "'], ",
+    switch(spec[[1]]$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, ")"),
+      DFOP = "k1, k2, g)",
+      SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else "", ")"),
+      HS = "k1, k2, tb)",
+      logistic = "kmax, k0, r)"
+    )
+  )
+
+  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"
+  )
+
+  deg_func <- function(odeini, odeparms, outtimes) {}
+
+  body(deg_func) <- parse(text = f_body)
+
+  return(deg_func)
+}
diff --git a/R/logistic.solution.R b/R/logistic.solution.R
deleted file mode 100644
index d9db13d7..00000000
--- a/R/logistic.solution.R
+++ /dev/null
@@ -1,59 +0,0 @@
-#' Logistic kinetics
-#' 
-#' Function describing exponential decline from a defined starting value, with
-#' an increasing rate constant, supposedly caused by microbial growth
-#' 
-#' @param t Time.
-#' @param parent.0 Starting value for the response variable at time zero.
-#' @param kmax Maximum rate constant.
-#' @param k0 Minumum rate constant effective at time zero.
-#' @param r Growth rate of the increase in the rate constant.
-#' @return The value of the response variable at time \code{t}.
-#' @note The solution of the logistic model reduces to the
-#'   \code{\link{SFO.solution}} if \code{k0} is equal to \code{kmax}.
-#' @references FOCUS (2014) \dQuote{Generic guidance for Estimating Persistence
-#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
-#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
-#'   Version 1.1, 18 December 2014
-#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
-#' @examples
-#' 
-#'   # Reproduce the plot on page 57 of FOCUS (2014)
-#'   plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.2),
-#'        from = 0, to = 100, ylim = c(0, 100),
-#'        xlab = "Time", ylab = "Residue")
-#'   plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.4),
-#'        from = 0, to = 100, add = TRUE, lty = 2, col = 2)
-#'   plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.8),
-#'        from = 0, to = 100, add = TRUE, lty = 3, col = 3)
-#'   plot(function(x) logistic.solution(x, 100, 0.08, 0.001, 0.2),
-#'        from = 0, to = 100, add = TRUE, lty = 4, col = 4)
-#'   plot(function(x) logistic.solution(x, 100, 0.08, 0.08, 0.2),
-#'        from = 0, to = 100, add = TRUE, lty = 5, col = 5)
-#'   legend("topright", inset = 0.05,
-#'          legend = paste0("k0 = ", c(0.0001, 0.0001, 0.0001, 0.001, 0.08),
-#'                          ", r = ", c(0.2, 0.4, 0.8, 0.2, 0.2)),
-#'          lty = 1:5, col = 1:5)
-#' 
-#'   # Fit with synthetic data
-#'   logistic <- mkinmod(parent = mkinsub("logistic"))
-#' 
-#'   sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-#'   parms_logistic <- c(kmax = 0.08, k0 = 0.0001, r = 0.2)
-#'   d_logistic <- mkinpredict(logistic,
-#'     parms_logistic, c(parent = 100),
-#'     sampling_times)
-#'   d_2_1 <- add_err(d_logistic,
-#'     sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
-#'     n = 1, reps = 2, digits = 5, LOD = 0.1, seed = 123456)[[1]]
-#' 
-#'   m <- mkinfit("logistic", d_2_1, quiet = TRUE)
-#'   plot_sep(m)
-#'   summary(m)$bpar
-#'   endpoints(m)$distimes
-#' 
-#' @export
-logistic.solution <- function(t, parent.0, kmax, k0, r)
-{
-	parent = parent.0 * (kmax / (kmax - k0 + k0 * exp (r * t))) ^(kmax/r)
-}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 2b7e71cb..54dd75c2 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -269,7 +269,7 @@ mkinfit <- function(mkinmod, observed,
     if (mkinmod[[1]] %in% parent_models_available) {
       speclist <- list(list(type = mkinmod, sink = TRUE))
       names(speclist) <- presumed_parent_name
-      mkinmod <- mkinmod(speclist = speclist)
+      mkinmod <- mkinmod(speclist = speclist, use_of_ff = "min")
     } else {
       stop("Argument mkinmod must be of class mkinmod or a string containing one of\n  ",
            paste(parent_models_available, collapse = ", "))
diff --git a/R/mkinmod.R b/R/mkinmod.R
index f52baa4f..21551861 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -101,7 +101,7 @@
 #' }
 #'
 #' @export mkinmod
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE)
+mkinmod <- function(..., use_of_ff = "max", speclist = NULL, quiet = FALSE, verbose = FALSE)
 {
   if (is.null(speclist)) spec <- list(...)
   else spec <- speclist
@@ -421,45 +421,8 @@ 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)
-    }
-  }
-  # }}}
+  # Attach a degradation function if an analytical solution is available
+  model$deg_func <- create_deg_func(spec, use_of_ff)
 
   class(model) <- "mkinmod"
   return(model)
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 0f8e83bb..df51dbe3 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -43,36 +43,36 @@
 #'
 #' SFO <- mkinmod(degradinol = mkinsub("SFO"))
 #' # Compare solution types
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
 #'       solution_type = "analytical")
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
 #'       solution_type = "deSolve")
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 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,
+#' mkinpredict(SFO, c(k_degradinol = 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,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
 #'       solution_type = "analytical")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
 #'       method = "lsoda")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
 #'       method = "ode45")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 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),
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100),
 #'       seq(0, 20, by = 0.1))[201,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
+#' mkinpredict(SFO, c(k_degradinol = 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"))
+#'                   m1 = list(type = "SFO"), use_of_ff = "min")
 #' if(require(rbenchmark)) {
 #'   benchmark(
 #'     eigen = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
diff --git a/R/nlme.R b/R/nlme.R
index 60cd2ca1..fa69db3c 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -13,15 +13,15 @@
 #' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
 #' m_SFO <- mkinmod(parent = mkinsub("SFO"))
 #' d_SFO_1 <- mkinpredict(m_SFO,
-#'   c(k_parent_sink = 0.1),
+#'   c(k_parent = 0.1),
 #'   c(parent = 98), sampling_times)
 #' d_SFO_1_long <- mkin_wide_to_long(d_SFO_1, time = "time")
 #' d_SFO_2 <- mkinpredict(m_SFO,
-#'   c(k_parent_sink = 0.05),
+#'   c(k_parent = 0.05),
 #'   c(parent = 102), sampling_times)
 #' d_SFO_2_long <- mkin_wide_to_long(d_SFO_2, time = "time")
 #' d_SFO_3 <- mkinpredict(m_SFO,
-#'   c(k_parent_sink = 0.02),
+#'   c(k_parent = 0.02),
 #'   c(parent = 103), sampling_times)
 #' d_SFO_3_long <- mkin_wide_to_long(d_SFO_3, time = "time")
 #'
diff --git a/R/parent_solutions.R b/R/parent_solutions.R
index c33d6d13..e02bcda7 100644
--- a/R/parent_solutions.R
+++ b/R/parent_solutions.R
@@ -136,7 +136,7 @@ DFOP.solution <- function(t, parent_0, k1, k2, g)
 #' between them.
 #' 
 #' @family parent solutions
-#' @inherit HS.solution 
+#' @inherit DFOP.solution 
 #' @param tb Break point. Before this time, exponential decline according to
 #'   \code{k1} is calculated, after this time, exponential decline proceeds
 #'   according to \code{k2}.
@@ -161,7 +161,7 @@ HS.solution <- function(t, parent_0, k1, k2, tb)
 #' and no substance in the bound fraction.
 #' 
 #' @family parent solutions
-#' @inherit HS.solution 
+#' @inherit SFO.solution 
 #' @param k_12 Kinetic constant describing transfer from free to bound.
 #' @param k_21 Kinetic constant describing transfer from bound to free.
 #' @param k_1output Kinetic constant describing degradation of the free