Make deSolve predictions within saemix robust

Also, exclude the saemix function when loading saemix in the example
code, to prevent overriding our generic

4 files changed, 63 insertions, 39 deletions

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 350ee56a..a6e7ca1c 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -34,6 +34,7 @@
 #'   the observed variables (default) or for all state variables (if set to
 #'   FALSE). Setting this to FALSE has no effect for analytical solutions,
 #'   as these always return mapped output.
+#' @param na_stop Should it be an error if deSolve::ode returns NaN values
 #' @param \dots Further arguments passed to the ode solver in case such a
 #'   solver is used.
 #' @import deSolve
@@ -121,6 +122,7 @@ mkinpredict.mkinmod <- function(x,
   solution_type = "deSolve",
   use_compiled = "auto",
   method.ode = "lsoda", atol = 1e-8, rtol = 1e-10,
+  na_stop = TRUE,
   map_output = TRUE, ...)
 {
 
@@ -208,9 +210,16 @@ mkinpredict.mkinmod <- function(x,
         ...
       )
     }
-    if (sum(is.na(out)) > 0) {
+    n_out_na <- sum(is.na(out))
+    if (n_out_na > 0 & na_stop) {
+      cat("odeini:\n")
+      print(odeini)
+      cat("odeparms:\n")
+      print(odeparms)
+      cat("out:\n")
+      print(out)
       stop("Differential equations were not integrated for all output times because\n",
-     "NaN values occurred in output from ode()")
+        n_out_na, " NaN values occurred in output from ode()")
     }
   }
 
diff --git a/R/saemix.R b/R/saemix.R
index 7a225601..280490a0 100644
--- a/R/saemix.R
+++ b/R/saemix.R
@@ -18,13 +18,18 @@
 #' @param object An [mmkin] row object containing several fits of the same
 #'   [mkinmod] model to different datasets
 #' @param verbose Should we print information about created objects?
+#' @param cores The number of cores to be used for multicore processing using
+#'   [parallel::mclapply()]. Using more than 1 core is experimental and may
+#'   lead to uncontrolled forking, apparently depending on the BLAS version
+#'   used.
 #' @param \dots Further parameters passed to [saemix::saemixData]
 #'   and [saemix::saemixModel].
 #' @return An [saemix::SaemixObject].
 #' @examples
 #' \dontrun{
-#' # We do not load the saemix package, as this would override our saemix
-#' # generic
+#' # We can load saemix, but should exclude the saemix function
+#' # as it would mask our generic version of it
+#' library(saemix, exclude = "saemix")
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")]))
 #' names(ds) <- paste("Dataset", 6:10)
@@ -37,18 +42,25 @@
 #' f_saemix_fomc <- saemix(f_mmkin_parent["FOMC", ])
 #' f_saemix_dfop <- saemix(f_mmkin_parent["DFOP", ])
 #'
-#' # We can use functions from the saemix package by prepending saemix::
-#' saemix::compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
+#' # As this returns an SaemixObject, we can use functions from saemix
+#' compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
 #'
 #' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
 #' f_saemix_fomc_tc <- saemix(f_mmkin_parent_tc["FOMC", ])
-#' saemix::compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
+#' compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
 #'
 #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
 #'   A1 = mkinsub("SFO"))
-#' f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
+#' f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
+#' # This takes about 4 minutes on my system
 #' f_saemix <- saemix(f_mmkin)
 #'
+#' # Using a single core, it takes about 6 minutes, using 10 cores it is slower
+#' # instead of faster
+#' f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
+#' f_saemix_des <- saemix(f_mmkin_des, cores = 1)
+#' compare.saemix(list(f_saemix, f_saemix_des))
+#'
 #' }
 #' @export
 saemix <- function(model, data, control, ...) UseMethod("saemix")
@@ -58,9 +70,10 @@ saemix <- function(model, data, control, ...) UseMethod("saemix")
 saemix.mmkin <- function(model, data,
   control = list(displayProgress = FALSE, print = FALSE,
     save = FALSE, save.graphs = FALSE),
+  cores = 1,
   verbose = FALSE, suppressPlot = TRUE, ...)
 {
-  m_saemix <- saemix_model(model, verbose = verbose)
+  m_saemix <- saemix_model(model, cores = cores, verbose = verbose)
   d_saemix <- saemix_data(model, verbose = verbose)
   if (suppressPlot) {
     # We suppress the log-likelihood curve that saemix currently
@@ -74,15 +87,10 @@ saemix.mmkin <- function(model, data,
     dev.off()
     unlink(tmp)
   }
-  class(result) <- c("saemix.mmkin", "saemix")
   return(result)
 }
 
 #' @rdname saemix
-#' @param cores The number of cores to be used for multicore processing using
-#'   [parallel::mclapply()]. Using more than 1 core is experimental and may
-#'   lead to uncontrolled forking, apparently depending on the BLAS version
-#'   used.
 #' @return An [saemix::SaemixModel] object.
 #' @export
 saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
@@ -203,7 +211,9 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
             out_wide <- mkinpredict(mkin_model,
               odeparms = odeparms, odeini = odeini,
               solution_type = solution_type,
-              outtimes = sort(unique(i_time)))
+              outtimes = sort(unique(i_time)),
+              na_stop = FALSE
+            )
 
             out_index <- cbind(as.character(i_time), as.character(i_name))
             out_values <- out_wide[out_index]
diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd
index 17916356..f2799bf4 100644
--- a/man/mkinpredict.Rd
+++ b/man/mkinpredict.Rd
@@ -30,6 +30,7 @@ mkinpredict(
   method.ode = "lsoda",
   atol = 1e-08,
   rtol = 1e-10,
+  na_stop = TRUE,
   map_output = TRUE,
   ...
 )
@@ -90,6 +91,8 @@ as these always return mapped output.}
 
 \item{\dots}{Further arguments passed to the ode solver in case such a
 solver is used.}
+
+\item{na_stop}{Should it be an error if deSolve::ode returns NaN values}
 }
 \value{
 A matrix with the numeric solution in wide format
diff --git a/man/saemix.Rd b/man/saemix.Rd
index 1959817a..a664b0cc 100644
--- a/man/saemix.Rd
+++ b/man/saemix.Rd
@@ -14,7 +14,9 @@ saemix(model, data, control, ...)
   data,
   control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs =
     FALSE),
+  cores = 1,
   verbose = FALSE,
+  suppressPlot = TRUE,
   ...
 )
 
@@ -30,15 +32,15 @@ internally from the \link{mmkin} object.}
 \item{\dots}{Further parameters passed to \link[saemix:saemixData]{saemix::saemixData}
 and \link[saemix:saemixModel]{saemix::saemixModel}.}
 
-\item{verbose}{Should we print information about created objects?}
-
-\item{object}{An \link{mmkin} row object containing several fits of the same
-\link{mkinmod} model to different datasets}
-
 \item{cores}{The number of cores to be used for multicore processing using
 \code{\link[parallel:mclapply]{parallel::mclapply()}}. Using more than 1 core is experimental and may
 lead to uncontrolled forking, apparently depending on the BLAS version
 used.}
+
+\item{verbose}{Should we print information about created objects?}
+
+\item{object}{An \link{mmkin} row object containing several fits of the same
+\link{mkinmod} model to different datasets}
 }
 \value{
 An \link[saemix:SaemixObject-class]{saemix::SaemixObject}.
@@ -63,39 +65,39 @@ Starting values for the fixed effects (population mean parameters, argument psi0
 }
 \examples{
 \dontrun{
-library(saemix)
+# We can load saemix, but should exclude the saemix function
+# as it would mask our generic version of it
+library(saemix, exclude = "saemix")
 ds <- lapply(experimental_data_for_UBA_2019[6:10],
  function(x) subset(x$data[c("name", "time", "value")]))
 names(ds) <- paste("Dataset", 6:10)
 f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
   state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
 f_saemix_p0_fixed <- saemix(f_mmkin_parent_p0_fixed)
-m_saemix_p0_fixed <- saemix_model(f_mmkin_parent_p0_fixed["FOMC", ])
-d_saemix_parent <- saemix_data(f_mmkin_parent_p0_fixed)
-saemix_options <- list(seed = 123456, displayProgress = FALSE,
-  save = FALSE, save.graphs = FALSE, nbiter.saemix = c(200, 80))
-f_saemix_p0_fixed <- saemix(m_saemix_p0_fixed, d_saemix_parent, saemix_options)
 
 f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
-m_saemix_sfo <- saemix_model(f_mmkin_parent["SFO", ])
-m_saemix_fomc <- saemix_model(f_mmkin_parent["FOMC", ])
-m_saemix_dfop <- saemix_model(f_mmkin_parent["DFOP", ])
-d_saemix_parent <- saemix_data(f_mmkin_parent["SFO", ])
-f_saemix_sfo <- saemix(m_saemix_sfo, d_saemix_parent, saemix_options)
-f_saemix_fomc <- saemix(m_saemix_fomc, d_saemix_parent, saemix_options)
-f_saemix_dfop <- saemix(m_saemix_dfop, d_saemix_parent, saemix_options)
+f_saemix_sfo <- saemix(f_mmkin_parent["SFO", ])
+f_saemix_fomc <- saemix(f_mmkin_parent["FOMC", ])
+f_saemix_dfop <- saemix(f_mmkin_parent["DFOP", ])
+
+# As this returns an SaemixObject, we can use functions from saemix
 compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
+
 f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
-m_saemix_fomc_tc <- saemix_model(f_mmkin_parent_tc["FOMC", ])
-f_saemix_fomc_tc <- saemix(m_saemix_fomc_tc, d_saemix_parent, saemix_options)
+f_saemix_fomc_tc <- saemix(f_mmkin_parent_tc["FOMC", ])
 compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
 
 dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
   A1 = mkinsub("SFO"))
-f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
-m_saemix <- saemix_model(f_mmkin)
-d_saemix <- saemix_data(f_mmkin)
-f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
+f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
+# This takes about 4 minutes on my system
+f_saemix <- saemix(f_mmkin)
+
+# Using a single core, it takes about 6 minutes, using 10 cores it is slower
+# instead of faster
+f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
+f_saemix_des <- saemix(f_mmkin_des, cores = 1)
+compare.saemix(list(f_saemix, f_saemix_des))
 
 }
 }