aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/mkinfit.R17
-rw-r--r--R/mkinmod.R4
-rw-r--r--R/mkinpredict.R107
-rw-r--r--R/nlme.mmkin.R14
-rw-r--r--R/saem.R22
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 d8740aed..b1fb57cb 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
diff --git a/R/saem.R b/R/saem.R
index 696ea0ee..5b8021de 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -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))

Contact - Imprint