diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 17 | ||||
| -rw-r--r-- | R/mkinmod.R | 4 | ||||
| -rw-r--r-- | R/mkinpredict.R | 107 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 14 | ||||
| -rw-r--r-- | R/saem.R | 22 |
5 files changed, 131 insertions, 33 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 693778fd..0d9246dd 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -500,6 +500,15 @@ mkinfit <- function(mkinmod, observed, } } + # Get native symbol before iterations info for speed + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (solution_type == "deSolve" & use_compiled[1] != FALSE) { + mkinmod$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = mkinmod$dll_info[["name"]])$address + mkinmod$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = mkinmod$dll_info[["name"]])$address + } + # Get the error model and the algorithm for fitting err_mod <- match.arg(error_model) error_model_algorithm = match.arg(error_model_algorithm) @@ -610,7 +619,8 @@ mkinfit <- function(mkinmod, observed, solution_type = solution_type, use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol = atol, rtol = rtol, + call_lsoda = call_lsoda, ...) observed_index <- cbind(as.character(observed$time), as.character(observed$name)) observed$predicted <- out[observed_index] @@ -892,7 +902,10 @@ mkinfit <- function(mkinmod, observed, fit$calls <- calls fit$time <- fit_time - # We also need the model and a model name for summary and plotting + # We also need the model and a model name for summary and plotting, + # but without address info that will become invalid + mkinmod$diffs_address <- NULL + mkinmod$initpar_address <- NULL fit$mkinmod <- mkinmod fit$mkinmod$name <- mkinmod_name fit$obs_vars <- obs_vars diff --git a/R/mkinmod.R b/R/mkinmod.R index 47307ab7..215dbed6 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -45,7 +45,9 @@ #' @param dll_dir Directory where an DLL object, if generated internally by #' [inline::cfunction()], should be saved. The DLL will only be stored in a #' permanent location for use in future sessions, if 'dll_dir' and 'name' -#' are specified. +#' are specified. This is helpful if fit objects are cached e.g. by knitr, +#' as the cache remains functional across sessions if the DLL is stored in +#' a user defined location. #' @param unload If a DLL from the target location in 'dll_dir' is already #' loaded, should that be unloaded first? #' @param overwrite If a file exists at the target DLL location in 'dll_dir', diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 0dc9cf55..11a3b35b 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -19,26 +19,24 @@ #' @param solution_type The method that should be used for producing the #' predictions. This should generally be "analytical" if there is only one #' observed variable, and usually "deSolve" in the case of several observed -#' variables. The third possibility "eigen" is faster but not applicable to -#' some models e.g. using FOMC for the parent compound. +#' variables. The third possibility "eigen" is fast in comparison to uncompiled +#' ODE models, but not applicable to some models, e.g. using FOMC for the +#' parent compound. #' @param method.ode The solution method passed via [mkinpredict] to [ode]] in -#' case the solution type is "deSolve". The default "lsoda" is performant, but -#' sometimes fails to converge. +#' case the solution type is "deSolve" and we are not using compiled code. #' @param use_compiled If set to \code{FALSE}, no compiled version of the #' [mkinmod] model is used, even if is present. -#' @param atol Absolute error tolerance, passed to [ode]. Default is 1e-8, -#' lower than in [lsoda]. -#' @param rtol Absolute error tolerance, passed to [ode]. Default is 1e-10, -#' much lower than in [lsoda]. -#' @param maxsteps Maximum number of steps, passed to [ode]. +#' @param atol Absolute error tolerance, passed to the ode solver. +#' @param rtol Absolute error tolerance, passed to the ode solver. +#' @param maxsteps Maximum number of steps, passed to the ode solver. #' @param map_output Boolean to specify if the output should list values for #' 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 [ode] returns NaN values +#' @param call_lsoda The address of the compiled function "call_lsoda" #' @param \dots Further arguments passed to the ode solver in case such a #' solver is used. -#' @import deSolve #' @return A matrix with the numeric solution in wide format #' @author Johannes Ranke #' @examples @@ -116,9 +114,10 @@ mkinpredict.mkinmod <- function(x, outtimes = seq(0, 120, by = 0.1), solution_type = "deSolve", use_compiled = "auto", - method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000, + method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000L, map_output = TRUE, na_stop = TRUE, + call_lsoda = NULL, ...) { @@ -173,20 +172,80 @@ mkinpredict.mkinmod <- function(x, if (solution_type == "deSolve") { if (!is.null(x$cf) & use_compiled[1] != FALSE) { - out <- deSolve::ode( - y = odeini, - times = outtimes, - func = "diffs", - initfunc = "initpar", - dllname = if (is.null(x$dll_info)) inline::getDynLib(x$cf)[["name"]] - else x$dll_info[["name"]], - parms = odeparms[x$parms], # Order matters when using compiled models - method = method.ode, - atol = atol, - rtol = rtol, - maxsteps = maxsteps, - ... +# out <- deSolve::ode( +# y = odeini, +# times = outtimes, +# func = "diffs", +# initfunc = "initpar", +# dllname = x$dll_info[["name"]], +# parms = odeparms[x$parms], # Order matters when using compiled models +# method = method.ode, +# atol = atol, +# rtol = rtol, +# maxsteps = maxsteps, +# ... +# ) +# + # Prepare call to "call_lsoda" + # Simplified code from deSolve::lsoda() adapted to our use case + if (is.null(call_lsoda)) { + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + } + if (is.null(x$diffs_address)) { + x$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = x$dll_info[["name"]])$address + x$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = x$dll_info[["name"]])$address + } + rwork <- vector("double", 20) + rwork[] <- 0. + rwork[6] <- max(abs(diff(outtimes))) + + iwork <- vector("integer", 20) + iwork[] <- 0 + iwork[6] <- maxsteps + + n <- length(odeini) + lmat <- n^2 + 2 # from deSolve::lsoda(), for Jacobian type full, internal (2) + # hard-coded default values of lsoda(): + maxordn <- 12L + maxords <- 5L + lrn <- 20 + n * (maxordn + 1) + 3 * n # length in case non-stiff method + lrs <- 20 + n * (maxords + 1) + 3 * n + lmat # length in case stiff method + lrw <- max(lrn, lrs) # actual length: max of both + liw <- 20 + n + + on.exit(.C("unlock_solver", PACKAGE = "deSolve")) + + out_raw <- .Call(call_lsoda, + as.double(odeini), # y + as.double(outtimes), # times + x$diffs_address, # derivfunc + as.double(odeparms[x$parms]), # parms + rtol, atol, + NULL, NULL, # rho, tcrit + NULL, # jacfunc + x$initpar_address, # initfunc + NULL, # eventfunc + 0L, # verbose + 1L, # iTask + as.double(rwork), # rWork + as.integer(iwork), # iWork + 2L, # jT full Jacobian calculated internally + 0L, # nOut + as.integer(lrw), # lRw + as.integer(liw), # lIw + 1L, # Solver + NULL, # rootfunc + 0L, as.double(0), 0L, # nRoot, Rpar, Ipar + 0L, # Type + list(fmat = 0L, tmat = 0L, imat = 0L, ModelForc = NULL), # flist + list(), # elist + list(islag = 0L) # elag ) + + out <- t(out_raw) + colnames(out) <- c("time", mod_vars) } else { mkindiff <- function(t, state, parms) { diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 09cb84b8..e193e5e3 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -186,10 +186,24 @@ nlme.mmkin <- function(model, data = "auto", thisCall[["control"]] <- control } + # Provide the address of call_lsoda to the fitting function + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (model[[1]]$solution_type == "deSolve" & !is.null(model[[1]]$mkinmod$cf)) { + # The mkinmod stored in the first fit will be used by nlme + model[[1]]$mkinmod$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address + model[[1]]$mkinmod$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address + } + fit_time <- system.time(val <- do.call("nlme.formula", thisCall)) val$time <- fit_time val$mkinmod <- model[[1]]$mkinmod + # Don't return addresses that will become invalid + val$mkinmod$diffs_address <- NULL + val$mkinmod$initpar_address <- NULL + val$data <- thisCall[["data"]] val$mmkin <- model if (is.list(start)) val$mean_dp_start <- start$fixed @@ -120,12 +120,12 @@ utils::globalVariables(c("predicted", "std")) #' summary(f_saem_dfop_sfo, data = TRUE) #' #' # 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)) +#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +#' nbiter.saemix = c(200, 80)) #' -#' #saemix::compare.saemix(list( -#' # f_saem_dfop_sfo$so, -#' # f_saem_dfop_sfo_deSolve$so)) +#' #anova( +#' # f_saem_dfop_sfo, +#' # f_saem_dfop_sfo_deSolve)) #' #' # If the model supports it, we can also use eigenvalue based solutions, which #' # take a similar amount of time @@ -580,6 +580,15 @@ saemix_model <- function(object, solution_type = "auto", transform_rates <- object[[1]]$transform_rates transform_fractions <- object[[1]]$transform_fractions + # Get native symbol info for speed + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (solution_type == "deSolve" & !is.null(mkin_model$cf)) { + mkin_model$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = mkin_model$dll_info[["name"]])$address + mkin_model$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = mkin_model$dll_info[["name"]])$address + } + # Define the model function model_function <- function(psi, id, xidep) { @@ -613,7 +622,8 @@ saemix_model <- function(object, solution_type = "auto", odeparms = odeparms, odeini = odeini, solution_type = solution_type, outtimes = sort(unique(i_time)), - na_stop = FALSE + na_stop = FALSE, + call_lsoda = call_lsoda ) out_index <- cbind(as.character(i_time), as.character(i_name)) |