diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 14 | ||||
| -rw-r--r-- | R/mkinmod.R | 4 | ||||
| -rw-r--r-- | R/mkinpredict.R | 37 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 14 | ||||
| -rw-r--r-- | R/saem.R | 22 |
5 files changed, 63 insertions, 28 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 693778fd..6cca5616 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -500,6 +500,13 @@ mkinfit <- function(mkinmod, observed, } } + # Get native symbol before iterations info for speed + if (solution_type == "deSolve" & use_compiled[1] != FALSE) { + mkinmod[["symbols"]] <- deSolve::checkDLL(dllname = mkinmod$dll_info[["name"]], + func = "diffs", initfunc = "initpar", + jacfunc = NULL, nout = 0, outnames = NULL) + } + # 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 +617,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 +900,9 @@ 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 symbols because they could become invalid + fit$symbols <- 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..f0d00319 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -19,18 +19,17 @@ #' @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. +#' When using compiled code, only lsoda is supported. #' @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, @@ -38,7 +37,6 @@ #' @param na_stop Should it be an error if [ode] returns NaN values #' @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 @@ -58,11 +56,11 @@ #' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, #' solution_type = "analytical")[21,] #' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, -#' method = "lsoda")[21,] +#' method = "lsoda", use_compiled = FALSE)[21,] #' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, -#' method = "ode45")[21,] +#' method = "ode45", use_compiled = FALSE)[21,] #' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, -#' method = "rk4")[21,] +#' method = "rk4", use_compiled = FALSE)[21,] #' # rk4 is not as precise here #' #' # The number of output times used to make a lot of difference until the @@ -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,20 @@ mkinpredict.mkinmod <- function(x, if (solution_type == "deSolve") { if (!is.null(x$cf) & use_compiled[1] != FALSE) { - out <- deSolve::ode( + out <- deSolve::lsoda( y = odeini, times = outtimes, - func = "diffs", + func = mkinmod[["symbols"]], initfunc = "initpar", - dllname = if (is.null(x$dll_info)) inline::getDynLib(x$cf)[["name"]] - else x$dll_info[["name"]], + dllname = x$dll_info[["name"]], parms = odeparms[x$parms], # Order matters when using compiled models - method = method.ode, atol = atol, rtol = rtol, maxsteps = maxsteps, ... ) + + 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)) |