Custom analytical solutions for saemix

Currently SFO-SFO and DFOP-SFO. Speed increase factor about 60

6 files changed, 183 insertions, 70 deletions

diff --git a/R/saemix.R b/R/saemix.R
index a038666e..1b373078 100644
--- a/R/saemix.R
+++ b/R/saemix.R
@@ -40,7 +40,7 @@
 #' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 #' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
 #'
-#' # The returned saem.mmkin object contains an SaemixObject, we can use
+#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 #' # functions from saemix
 #' library(saemix)
 #' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
@@ -49,17 +49,26 @@
 #' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
 #' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
 #'
+#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+#'   A1 = mkinsub("SFO"))
+#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+#'   A1 = mkinsub("SFO"))
 #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
 #'   A1 = mkinsub("SFO"))
-#' f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
-#' # This takes about 4 minutes on my system
-#' f_saem <- saem(f_mmkin)
+#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
+#' # and compiled ODEs for FOMC, both are fast
+#' f_mmkin <- mmkin(list(
+#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#'   ds, quiet = TRUE)
+#' # These take about five seconds each on this system, as we use
+#' # analytical solutions written for saemix. When using the analytical
+#' # solutions written for mkin this took around four minutes
+#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
+#' f_saem_dfop_sfo <- saem(f_mmkin["SFO-SFO", ])
 #'
-#' f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
 #' # Using a single core, the following takes about 6 minutes, using 10 cores
 #' # it is slower instead of faster
-#' f_saem_des <- saem(f_mmkin_des, cores = 1)
-#' compare.saemix(list(f_saem$so, f_saem_des$so))
+#' f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
 #' }
 #' @export
 saem <- function(object, control, ...) UseMethod("saem")
@@ -83,7 +92,12 @@ saem.mmkin <- function(object,
   }
   fit_time <- system.time({
     f_saemix <- saemix::saemix(m_saemix, d_saemix, control)
-    f_saemix <- try(saemix::saemix.predict(f_saemix), silent = TRUE)
+    f_pred <- try(saemix::saemix.predict(f_saemix), silent = TRUE)
+    if (!inherits(f_pred, "try-error")) {
+      f_saemix <- f_pred
+    } else {
+      warning("Creating predictions from the saemix model failed")
+    }
   })
   if (suppressPlot) {
     grDevices::dev.off()
@@ -145,37 +159,43 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
 
   model_function <- FALSE
 
-  if (length(mkin_model$spec) == 1 & mkin_model$use_of_ff == "max") {
-    parent_type <- mkin_model$spec[[1]]$type
-    if (length(odeini_fixed) == 1) {
-      if (parent_type == "SFO") {
-        stop("saemix needs at least two parameters to work on.")
-      }
-      if (parent_type == "FOMC") {
-        model_function <- function(psi, id, xidep) {
-          odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
+  # Model functions with analytical solutions
+  # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here
+  # In general, we need to consider exactly how the parameters in mkinfit were specified,
+  # as the parameters are currently mapped by position in these solutions
+  sinks <- sapply(mkin_model$spec, function(x) x$sink)
+  if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) {
+    # Parent only
+    if (length(mkin_model$spec) == 1) {
+      parent_type <- mkin_model$spec[[1]]$type
+      if (length(odeini_fixed) == 1) {
+        if (parent_type == "SFO") {
+          stop("saemix needs at least two parameters to work on.")
         }
-      }
-      if (parent_type == "DFOP") {
-        model_function <- function(psi, id, xidep) {
-          g <- plogis(psi[id, 3])
-          t = xidep[, "time"]
-          odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
-            (1 - g) * exp(- exp(psi[id, 2]) * t))
+        if (parent_type == "FOMC") {
+          model_function <- function(psi, id, xidep) {
+            odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
+          }
         }
-      }
-      if (parent_type == "HS") {
-        model_function <- function(psi, id, xidep) {
-          tb <- exp(psi[id, 3])
-          t = xidep[, "time"]
-          k1 = exp(psi[id, 1])
-          odeini_fixed * ifelse(t <= tb,
-            exp(- k1 * t),
-            exp(- k1 * t) * exp(- exp(psi[id, 2]) * (t - tb)))
+        if (parent_type == "DFOP") {
+          model_function <- function(psi, id, xidep) {
+            g <- plogis(psi[id, 3])
+            t <- xidep[, "time"]
+            odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
+              (1 - g) * exp(- exp(psi[id, 2]) * t))
+          }
         }
-      }
-    } else {
-      if (length(odeparms_fixed) == 0) {
+        if (parent_type == "HS") {
+          model_function <- function(psi, id, xidep) {
+            tb <- exp(psi[id, 3])
+            t <- xidep[, "time"]
+            k1 = exp(psi[id, 1])
+            odeini_fixed * ifelse(t <= tb,
+              exp(- k1 * t),
+              exp(- k1 * t) * exp(- exp(psi[id, 2]) * (t - tb)))
+          }
+        }
+      } else {
         if (parent_type == "SFO") {
           model_function <- function(psi, id, xidep) {
             psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
@@ -189,7 +209,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
         if (parent_type == "DFOP") {
           model_function <- function(psi, id, xidep) {
             g <- plogis(psi[id, 4])
-            t = xidep[, "time"]
+            t <- xidep[, "time"]
             psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
               (1 - g) * exp(- exp(psi[id, 3]) * t))
           }
@@ -197,7 +217,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
         if (parent_type == "HS") {
           model_function <- function(psi, id, xidep) {
             tb <- exp(psi[id, 4])
-            t = xidep[, "time"]
+            t <- xidep[, "time"]
             k1 = exp(psi[id, 2])
             psi[id, 1] * ifelse(t <= tb,
               exp(- k1 * t),
@@ -206,9 +226,57 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
         }
       }
     }
+
+    # Parent with one metabolite
+    # Parameter names used in the model functions are as in
+    # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
+    if (length(mkin_model$spec) == 2) {
+      types <- unname(sapply(mkin_model$spec, function(x) x$type))
+      # Initial value for the metabolite (n20) must be fixed
+      if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
+        n20 <- odeini_fixed
+        parent_name <- names(mkin_model$spec)[1]
+        if (identical(types, c("SFO", "SFO"))) {
+          model_function <- function(psi, id, xidep) {
+            t <- xidep[, "time"]
+            n10 <- psi[id, 1]
+            k1 <- exp(psi[id, 2])
+            k2 <- exp(psi[id, 3])
+            f12 <- plogis(psi[id, 4])
+            ifelse(xidep[, "name"] == parent_name,
+              n10 * exp(- k1 * t),
+              (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+                (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+            )
+          }
+        }
+        if (identical(types, c("DFOP", "SFO"))) {
+          model_function <- function(psi, id, xidep) {
+            t <- xidep[, "time"]
+            n10 <- psi[id, 1]
+            k2 <- exp(psi[id, 2])
+            f12 <- plogis(psi[id, 3])
+            l1 <- exp(psi[id, 4])
+            l2 <- exp(psi[id, 5])
+            g <- plogis(psi[id, 6])
+            ifelse(xidep[, "name"] == parent_name,
+              n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+              ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+                (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+                ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+                    ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+                      f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+                ((l1 - k2) * l2 - k2 * l1 + k2^2)
+            )
+          }
+        }
+      }
+    }
   }
 
-  if (!is.function(model_function)) {
+  if (is.function(model_function)) {
+    solution_type = "analytical saemix"
+  } else {
     model_function <- function(psi, id, xidep) {
 
       uid <- unique(id)
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index f7110dd0..4c9776a6 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -96,7 +96,7 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   colnames(confint_ranef)[1] <- "est."
 
   # Error model
-  enames <- object$so@results@name.sigma
+  enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1")
   confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")])
   colnames(confint_errmod)[1] <- "est."
 
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index d9720761..8d117fb1 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -10,7 +10,7 @@ articles:
   web_only/NAFTA_examples: NAFTA_examples.html
   web_only/benchmarks: benchmarks.html
   web_only/compiled_models: compiled_models.html
-last_built: 2020-11-08T02:20Z
+last_built: 2020-11-09T06:03Z
 urls:
   reference: https://pkgdown.jrwb.de/mkin/reference
   article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png
index cfc5bc2b..17a35806 100644
--- a/docs/dev/reference/Rplot001.png
+++ b/docs/dev/reference/Rplot001.png
diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html
index 06fcfaa7..25608fc8 100644
--- a/docs/dev/reference/saem.html
+++ b/docs/dev/reference/saem.html
@@ -229,28 +229,28 @@ using <a href='mmkin.html'>mmkin</a>.</p>
   state.ini <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fl'>100</span><span class='op'>)</span>, fixed_initials <span class='op'>=</span> <span class='st'>"parent"</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
 <span class='va'>f_saem_p0_fixed</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent_p0_fixed</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:44:42 2020"
+#&gt; [1] "Mon Nov  9 07:04:09 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:44:43 2020"</div><div class='input'>
+#&gt; [1] "Mon Nov  9 07:04:11 2020"</div><div class='input'>
 <span class='va'>f_mmkin_parent</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='st'>"FOMC"</span>, <span class='st'>"DFOP"</span><span class='op'>)</span>, <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
 <span class='va'>f_saem_sfo</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent</span><span class='op'>[</span><span class='st'>"SFO"</span>, <span class='op'>]</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:44:45 2020"
+#&gt; [1] "Mon Nov  9 07:04:12 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:44:46 2020"</div><div class='input'><span class='va'>f_saem_fomc</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent</span><span class='op'>[</span><span class='st'>"FOMC"</span>, <span class='op'>]</span><span class='op'>)</span>
+#&gt; [1] "Mon Nov  9 07:04:13 2020"</div><div class='input'><span class='va'>f_saem_fomc</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent</span><span class='op'>[</span><span class='st'>"FOMC"</span>, <span class='op'>]</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:44:47 2020"
+#&gt; [1] "Mon Nov  9 07:04:14 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:44:49 2020"</div><div class='input'><span class='va'>f_saem_dfop</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent</span><span class='op'>[</span><span class='st'>"DFOP"</span>, <span class='op'>]</span><span class='op'>)</span>
+#&gt; [1] "Mon Nov  9 07:04:16 2020"</div><div class='input'><span class='va'>f_saem_dfop</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent</span><span class='op'>[</span><span class='st'>"DFOP"</span>, <span class='op'>]</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:44:49 2020"
+#&gt; [1] "Mon Nov  9 07:04:16 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:44:52 2020"</div><div class='input'>
-<span class='co'># The returned saem.mmkin object contains an SaemixObject, we can use</span>
+#&gt; [1] "Mon Nov  9 07:04:19 2020"</div><div class='input'>
+<span class='co'># The returned saem.mmkin object contains an SaemixObject, therefore we can use</span>
 <span class='co'># functions from saemix</span>
 <span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'>saemix</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; <span class='message'>Package saemix, version 3.1.9000</span>
@@ -262,33 +262,69 @@ using <a href='mmkin.html'>mmkin</a>.</p>
 <span class='va'>f_mmkin_parent_tc</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span><span class='op'>(</span><span class='va'>f_mmkin_parent</span>, error_model <span class='op'>=</span> <span class='st'>"tc"</span><span class='op'>)</span>
 <span class='va'>f_saem_fomc_tc</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_parent_tc</span><span class='op'>[</span><span class='st'>"FOMC"</span>, <span class='op'>]</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:44:54 2020"
+#&gt; [1] "Mon Nov  9 07:04:21 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:44:59 2020"</div><div class='input'><span class='fu'><a href='https://rdrr.io/pkg/saemix/man/compare.saemix.html'>compare.saemix</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='va'>f_saem_fomc</span><span class='op'>$</span><span class='va'>so</span>, <span class='va'>f_saem_fomc_tc</span><span class='op'>$</span><span class='va'>so</span><span class='op'>)</span><span class='op'>)</span>
+#&gt; [1] "Mon Nov  9 07:04:26 2020"</div><div class='input'><span class='fu'><a href='https://rdrr.io/pkg/saemix/man/compare.saemix.html'>compare.saemix</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='va'>f_saem_fomc</span><span class='op'>$</span><span class='va'>so</span>, <span class='va'>f_saem_fomc_tc</span><span class='op'>$</span><span class='va'>so</span><span class='op'>)</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Likelihoods computed by importance sampling </div><div class='output co'>#&gt;        AIC      BIC
 #&gt; 1 467.7644 465.0305
 #&gt; 2 469.4862 466.3617</div><div class='input'>
-<span class='va'>dfop_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"DFOP"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
+<span class='va'>sfo_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
   A1 <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span><span class='op'>)</span>
-</div><div class='output co'>#&gt; <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='va'>f_mmkin</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>dfop_sfo</span><span class='op'>)</span>, <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span>, solution_type <span class='op'>=</span> <span class='st'>"analytical"</span><span class='op'>)</span>
-<span class='co'># This takes about 4 minutes on my system</span>
-<span class='va'>f_saem</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='va'>fomc_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"FOMC"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
+  A1 <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='va'>dfop_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"DFOP"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
+  A1 <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='co'># The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,</span>
+<span class='co'># and compiled ODEs for FOMC, both are fast</span>
+<span class='va'>f_mmkin</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>
+    <span class='st'>"SFO-SFO"</span> <span class='op'>=</span> <span class='va'>sfo_sfo</span>, <span class='st'>"FOMC-SFO"</span> <span class='op'>=</span> <span class='va'>fomc_sfo</span>, <span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>dfop_sfo</span><span class='op'>)</span>,
+  <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
+<span class='co'># These take about five seconds each on this system, as we use</span>
+<span class='co'># analytical solutions written for saemix. When using the analytical</span>
+<span class='co'># solutions written for mkin this took around four minutes</span>
+<span class='va'>f_saem_sfo_sfo</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin</span><span class='op'>[</span><span class='st'>"SFO-SFO"</span>, <span class='op'>]</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; Running main SAEM algorithm
+#&gt; [1] "Mon Nov  9 07:04:28 2020"
+#&gt; ....
+#&gt;     Minimisation finished
+#&gt; [1] "Mon Nov  9 07:04:33 2020"</div><div class='input'><span class='va'>f_saem_dfop_sfo</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin</span><span class='op'>[</span><span class='st'>"SFO-SFO"</span>, <span class='op'>]</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:45:00 2020"
+#&gt; [1] "Mon Nov  9 07:04:33 2020"
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:49:01 2020"</div><div class='output co'>#&gt; <span class='error'>Error in exp(trans_k): non-numeric argument to mathematical function</span></div><div class='output co'>#&gt; <span class='message'>Timing stopped at: 317.7 169.9 258.9</span></div><div class='input'>
-<span class='va'>f_mmkin_des</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>dfop_sfo</span><span class='op'>)</span>, <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span>, solution_type <span class='op'>=</span> <span class='st'>"deSolve"</span><span class='op'>)</span>
+#&gt; [1] "Mon Nov  9 07:04:39 2020"</div><div class='input'>
 <span class='co'># Using a single core, the following takes about 6 minutes, using 10 cores</span>
 <span class='co'># it is slower instead of faster</span>
-<span class='va'>f_saem_des</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin_des</span>, cores <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span>
+<span class='va'>f_saem_fomc</span> <span class='op'>&lt;-</span> <span class='fu'>saem</span><span class='op'>(</span><span class='va'>f_mmkin</span><span class='op'>[</span><span class='st'>"FOMC-SFO"</span>, <span class='op'>]</span>, cores <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; Running main SAEM algorithm
-#&gt; [1] "Sun Nov  8 02:49:20 2020"
+#&gt; [1] "Mon Nov  9 07:04:39 2020"
+#&gt; DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
+#&gt;       taken on this call before reaching TOUT     
+#&gt; In above message, I1 = 5000
+#&gt;  
+#&gt; In above message, R1 = 0.00156238
+#&gt;  
+#&gt; DLSODA-  At T (=R1) and step size H (=R2), the    
+#&gt;       corrector convergence failed repeatedly     
+#&gt;       or with ABS(H) = HMIN   
+#&gt; In above message, R1 = 0, R2 = 1.1373e-10
+#&gt;  
+#&gt; DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
+#&gt;       taken on this call before reaching TOUT     
+#&gt; In above message, I1 = 5000
+#&gt;  
+#&gt; In above message, R1 = 2.24752e-06
+#&gt;  
+#&gt; DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
+#&gt;       taken on this call before reaching TOUT     
+#&gt; In above message, I1 = 5000
+#&gt;  
+#&gt; In above message, R1 = 0.000585935
+#&gt;  
 #&gt; ....
 #&gt;     Minimisation finished
-#&gt; [1] "Sun Nov  8 02:57:30 2020"</div><div class='output co'>#&gt; <span class='error'>Error in exp(trans_k): non-numeric argument to mathematical function</span></div><div class='output co'>#&gt; <span class='message'>Timing stopped at: 591.7 205.6 521.4</span></div><div class='input'><span class='fu'><a href='https://rdrr.io/pkg/saemix/man/compare.saemix.html'>compare.saemix</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='va'>f_saem</span><span class='op'>$</span><span class='va'>so</span>, <span class='va'>f_saem_des</span><span class='op'>$</span><span class='va'>so</span><span class='op'>)</span><span class='op'>)</span>
-</div><div class='output co'>#&gt; <span class='error'>Error in compare.saemix(list(f_saem$so, f_saem_des$so)): object 'f_saem' not found</span></div><div class='input'><span class='co'># }</span>
+#&gt; [1] "Mon Nov  9 07:11:24 2020"</div><div class='output co'>#&gt; <span class='warning'>Warning: Creating predictions from the saemix model failed</span></div><div class='input'><span class='co'># }</span>
 </div></pre>
   </div>
   <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">
diff --git a/man/saem.Rd b/man/saem.Rd
index 96b8b55a..b2daf419 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -77,7 +77,7 @@ f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
 f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
 
-# The returned saem.mmkin object contains an SaemixObject, we can use
+# The returned saem.mmkin object contains an SaemixObject, therefore we can use
 # functions from saemix
 library(saemix)
 compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
@@ -86,17 +86,26 @@ f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
 f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
 compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
 
+sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+  A1 = mkinsub("SFO"))
+fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+  A1 = mkinsub("SFO"))
 dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
   A1 = mkinsub("SFO"))
-f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
-# This takes about 4 minutes on my system
-f_saem <- saem(f_mmkin)
+# The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
+# and compiled ODEs for FOMC, both are fast
+f_mmkin <- mmkin(list(
+    "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+  ds, quiet = TRUE)
+# These take about five seconds each on this system, as we use
+# analytical solutions written for saemix. When using the analytical
+# solutions written for mkin this took around four minutes
+f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
+f_saem_dfop_sfo <- saem(f_mmkin["SFO-SFO", ])
 
-f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
 # Using a single core, the following takes about 6 minutes, using 10 cores
 # it is slower instead of faster
-f_saem_des <- saem(f_mmkin_des, cores = 1)
-compare.saemix(list(f_saem$so, f_saem_des$so))
+f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
 }
 }
 \seealso{