From b8b60ef92c605e862294fd29c51e20e31e0ded81 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Dec 2020 21:29:18 +0100 Subject: Make saem using mkinpredict work again I threw out mclapply as it did not play well with the linear algebra routines used in the saemix code. Most of the change is actually indentation in the code creating the model function. But there is an important fix in mkinpredict which I had broken. --- R/mkinpredict.R | 4 ++- R/saem.R | 97 ++++++++++++++++++++++++++++++--------------------------- man/saem.Rd | 30 ++++++++++-------- test.log | 14 ++++----- 4 files changed, 78 insertions(+), 67 deletions(-) diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 4b618490..3402a7ba 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -134,6 +134,8 @@ mkinpredict.mkinmod <- function(x, dimnames = list(as.character(outtimes), c("time", obs_vars))) out_obs[, "time"] <- outtimes + n_out_na <- 0 # to check if we get NA values with deSolve + if (solution_type == "analytical") { # This is clumsy, as we wanted fast analytical predictions for mkinfit, # which bypasses mkinpredict in the case of analytical solutions @@ -222,7 +224,7 @@ mkinpredict.mkinmod <- function(x, if (map_output) { # Output transformation for models with unobserved compartments like SFORB # if not already mapped in analytical solution - if (!na_stop) { + if (n_out_na > 0 & !na_stop) { available <- c(TRUE, rep(FALSE, length(outtimes) - 1)) } else { available <- rep(TRUE, length(outtimes)) diff --git a/R/saem.R b/R/saem.R index fda5569e..a21f9990 100644 --- a/R/saem.R +++ b/R/saem.R @@ -26,10 +26,6 @@ utils::globalVariables(c("predicted", "std")) #' automatic choice is not desired #' @param quiet Should we suppress the messages saemix prints at the beginning #' and the end of the optimisation process? -#' @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 excessive forking, apparently depending on the BLAS version -#' used. #' @param suppressPlot Should we suppress any plotting that is done #' by the saemix function? #' @param control Passed to [saemix::saemix] @@ -43,7 +39,7 @@ utils::globalVariables(c("predicted", "std")) #' 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, +#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, #' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) #' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed) #' @@ -76,9 +72,10 @@ utils::globalVariables(c("predicted", "std")) #' 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 +#' # saem fits of SFO-SFO and DFOP-SFO to these data 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["DFOP-SFO", ]) #' # We can use print, plot and summary methods to check the results @@ -86,10 +83,18 @@ utils::globalVariables(c("predicted", "std")) #' plot(f_saem_dfop_sfo) #' summary(f_saem_dfop_sfo, data = TRUE) #' -#' # Using a single core, the following takes about 6 minutes as we do not have an -#' # analytical solution. Using 10 cores it is slower instead of faster -#' f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1) -#' plot(f_saem_fomc) +#' # The following takes about 6 minutes +#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) +#' +#' #saemix::compare.saemix(list( +#' # f_saem_dfop_sfo$so, +#' # f_saem_dfop_sfo_deSolve$so)) +#' +#' # If the model supports it, we can also use eigenvalue based solutions, which +#' # take a similar amount of time +#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen", +#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) #' } #' @export saem <- function(object, ...) UseMethod("saem") @@ -115,14 +120,13 @@ saem.mmkin <- function(object, # that we discard afterwards tmp <- tempfile() grDevices::png(tmp) + on.exit(grDevices::dev.off()) + on.exit(unlink(tmp), add = TRUE) } fit_time <- system.time({ utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) }) - if (suppressPlot) { - grDevices::dev.off() - unlink(tmp) - } + transparms_optim <- f_saemix@results@fixed.effects names(transparms_optim) <- f_saemix@results@name.fixed @@ -417,44 +421,45 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk uid <- unique(id) - res_list <- parallel::mclapply(uid, function(i) { - transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict - names(transparms_optim) <- names(degparms_optim) + res_list <- lapply(uid, function(i) { - odeini_optim <- transparms_optim[odeini_optim_parm_names] - names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names) + transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict + names(transparms_optim) <- names(degparms_optim) - odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)] + odeini_optim <- transparms_optim[odeini_optim_parm_names] + names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names) - ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names) - odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model, - transform_rates = object[[1]]$transform_rates, - transform_fractions = object[[1]]$transform_fractions) - odeparms <- c(odeparms_optim, odeparms_fixed) + odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)] - xidep_i <- subset(xidep, id == i) + ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names) + odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model, + transform_rates = object[[1]]$transform_rates, + transform_fractions = object[[1]]$transform_fractions) + odeparms <- c(odeparms_optim, odeparms_fixed) - if (solution_type == "analytical") { - out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms) - } else { + xidep_i <- subset(xidep, id == i) - i_time <- xidep_i$time - i_name <- xidep_i$name + if (solution_type == "analytical") { + out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms) + } else { - out_wide <- mkinpredict(mkin_model, - odeparms = odeparms, odeini = odeini, - solution_type = solution_type, - outtimes = sort(unique(i_time)), - na_stop = FALSE - ) + i_time <- xidep_i$time + i_name <- xidep_i$name - out_index <- cbind(as.character(i_time), as.character(i_name)) - out_values <- out_wide[out_index] - } - return(out_values) - }, mc.cores = cores) - res <- unlist(res_list) - return(res) + out_wide <- mkinpredict(mkin_model, + odeparms = odeparms, odeini = odeini, + solution_type = solution_type, + 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] + } + return(out_values) + }) + res <- unlist(res_list) + return(res) } } diff --git a/man/saem.Rd b/man/saem.Rd index d8d6ea0f..1fdfb2c4 100644 --- a/man/saem.Rd +++ b/man/saem.Rd @@ -53,11 +53,6 @@ automatic choice is not desired} \item{control}{Passed to \link[saemix:saemix]{saemix::saemix}} -\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 excessive forking, apparently depending on the BLAS version -used.} - \item{verbose}{Should we print information about created objects of type \link[saemix:SaemixModel-class]{saemix::SaemixModel} and \link[saemix:SaemixData-class]{saemix::SaemixData}?} @@ -98,7 +93,7 @@ using \link{mmkin}. 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, +f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed) @@ -131,9 +126,10 @@ dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), 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 +# saem fits of SFO-SFO and DFOP-SFO to these data 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["DFOP-SFO", ]) # We can use print, plot and summary methods to check the results @@ -141,10 +137,18 @@ print(f_saem_dfop_sfo) plot(f_saem_dfop_sfo) summary(f_saem_dfop_sfo, data = TRUE) -# Using a single core, the following takes about 6 minutes as we do not have an -# analytical solution. Using 10 cores it is slower instead of faster -f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1) -plot(f_saem_fomc) +# The following takes about 6 minutes +#f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) + +#saemix::compare.saemix(list( +# f_saem_dfop_sfo$so, +# f_saem_dfop_sfo_deSolve$so)) + +# If the model supports it, we can also use eigenvalue based solutions, which +# take a similar amount of time +#f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen", +# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) } } \seealso{ diff --git a/test.log b/test.log index 6105128f..13fa7e70 100644 --- a/test.log +++ b/test.log @@ -6,32 +6,32 @@ Testing mkin ✔ | 2 | Export dataset for reading into CAKE ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.0 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s] -✔ | 7 | Fitting the SFORB model [3.5 s] -✔ | 5 | Analytical solutions for coupled models [3.2 s] +✔ | 7 | Fitting the SFORB model [3.7 s] +✔ | 5 | Analytical solutions for coupled models [3.3 s] ✔ | 5 | Calculation of Akaike weights ✔ | 14 | Confidence intervals and p-values [1.2 s] ✔ | 14 | Error model fitting [4.7 s] ✔ | 3 | Time step normalisation ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] ✔ | 1 | Fitting the logistic model [0.2 s] -✔ | 31 | Nonlinear mixed effects models [27.5 s] +✔ | 31 | Nonlinear mixed effects models [27.4 s] ✔ | 2 | Test dataset classes mkinds and mkindsg ✔ | 1 | mkinfit features [0.3 s] ✔ | 12 | Special cases of mkinfit calls [0.8 s] ✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.4 s] ✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.6 s] -✔ | 9 | Nonlinear mixed-effects models [7.8 s] +✔ | 9 | Nonlinear mixed-effects models [7.9 s] ✔ | 16 | Plotting [1.9 s] ✔ | 4 | Residuals extracted from mkinfit models ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] ✔ | 4 | Summary [0.1 s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] -✔ | 9 | Hypothesis tests [7.3 s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] +✔ | 9 | Hypothesis tests [7.5 s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 69.3 s +Duration: 69.5 s [ FAIL 0 | WARN 0 | SKIP 0 | PASS 204 ] -- cgit v1.2.1