Better support metabolites in `sawtooth`

Also, add more example code to make it easier to use
`one_box`, `sawtooth` and `pfm_degradation`.

However, for generating sawtooth curves with specified
output times, overlaying the output of several calls to
`mkinpredict` is probably still easier.

2 files changed, 64 insertions, 4 deletions

diff --git a/R/pfm_degradation.R b/R/pfm_degradation.R
index 832a797..63a6e19 100644
--- a/R/pfm_degradation.R
+++ b/R/pfm_degradation.R
@@ -9,21 +9,48 @@
 #' @param years For how many years should the degradation be predicted?
 #' @param step_days What step size in days should the output have?
 #' @param times The output times
+#' @return A data frame containing the output times and the concentrations
+#'   assuming initial concentrations of 1 for the parent and zero for
+#'   metabolites, if any.
 #' @export
 #' @author Johannes Ranke
 #' @examples
-#' head(pfm_degradation("SFO", DT50 = 10))
+#' # Simple example of an SFO decline curve
+#' sfo_out <- pfm_degradation("SFO", DT50 = 10)
+#' head(sfo_out)
+#'
+#' # Fictive example with a metabolite where we first generate an SFO-SFO model
+#' sfo_sfo <- mkinmod(
+#'   parent = mkinsub("SFO", to = "metabolite"),
+#'   metabolite = mkinsub("SFO"))
+#'
+#' sfo_sfo_out <- pfm_degradation(sfo_sfo,
+#'   parms = c(k_parent = 0.1, f_parent_to_metabolite = 0.5, k_metabolite = 0.02))
+#'
+#' plot(
+#'   sfo_sfo_out[, "time"],
+#'   sfo_sfo_out[, "parent"], type = "l",
+#'   xlab = "Time", ylab = "Relative concentration",
+#'   xlim = c(0, 100))
+#' lines(
+#'   sfo_sfo_out[, "time"],
+#'   sfo_sfo_out[, "metabolite"], lty = 2)
+#'
 pfm_degradation <- function(model = "SFO",
   DT50 = 1000, parms = c(k_parent = log(2)/DT50),
   years = 1, step_days = 1,
   times = seq(0, years * 365, by = step_days))
 {
-  if (model %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) {
-    model <- mkinmod(parent = list(type = model))
+  if (!inherits(model, "mkinmod")) {
+    if (model[1] %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) {
+      model <- mkinmod(parent = list(type = model))
+    } else {
+      stop("Please specify the model with a suitable name or an mkinmod object")
+    }
   }
   initial_state = c(1, rep(0, length(model$diffs) - 1))
   names(initial_state) <- names(model$diffs)
-  time_course <- mkinpredict(model, odeparms = parms, 
+  time_course <- mkinpredict(model, odeparms = parms,
     odeini = initial_state,
     outtimes = times,
     solution_type = ifelse(length(model$spec) == 1, "analytical", "deSolve"))
diff --git a/R/twa.R b/R/twa.R
index 0506334..8ca490d 100644
--- a/R/twa.R
+++ b/R/twa.R
@@ -82,8 +82,41 @@ one_box.character <- function(x, ini = 1, parms, ...,
 }
 
 #' @rdname one_box
+#' @param odeparms Will be passed to [mkinpredict]
+#' @param solution_type Will be passed to [mkinpredict]
 #' @importFrom mkin mkinpredict
 #' @export
+#' @examples
+#' library(mkin)
+#' SFO_SFO <- mkinmod(
+#'   parent = mkinsub("SFO", to = "m1"),
+#'   m1 = mkinsub("SFO"))
+#' c_0 = c(parent = 100, m1 = 0)
+#' deg_parms = c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01)
+#' sfo_sfo_box <- one_box(SFO_SFO, odeparms = deg_parms, ini = c_0)
+#' head(sfo_sfo_box)
+one_box.mkinmod <- function(x,
+  ini = c(1, rep(0, length(x$diffs) - 1)),
+  odeparms,
+  solution_type = "deSolve",
+  ...,
+  t_end = 100, res = 0.01)
+{
+
+  # Set names for initial amounts to be compatible with the model
+  if (is.null(names(ini))) names(ini) <- names(x$diffs)
+
+  t_out = seq(0, t_end, by = res)
+
+  tmp <- mkinpredict(x, odeparms = odeparms, odeini = ini,
+    outtimes = t_out, solution_type = solution_type)[, -1, drop = FALSE]
+  result <- ts(tmp, 0, t_end, frequency = 1/res)
+  class(result) <- c("one_box", "ts")
+  return(result)
+}
+
+#' @rdname one_box
+#' @export
 one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) {
   fit <- x
   if (ini[1] == "model") {