From cda47972e2b6a9610e3118dcd2270d7a1c76de3d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 7 Nov 2020 11:54:13 +0100 Subject: Make deSolve predictions within saemix robust Also, exclude the saemix function when loading saemix in the example code, to prevent overriding our generic --- R/mkinpredict.R | 13 +++++++++++-- R/saemix.R | 36 +++++++++++++++++++++++------------- man/mkinpredict.Rd | 3 +++ man/saemix.Rd | 50 ++++++++++++++++++++++++++------------------------ 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)) } } -- cgit v1.2.1