From 19861806ed790dc55380c64c6a5e27ba7ecd52de Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 13 Feb 2023 06:23:02 +0100 Subject: WIP adapting to new deSolve with faster lsoda --- DESCRIPTION | 2 +- R/mkinfit.R | 13 ++--- R/mkinpredict.R | 90 +++++------------------------- vignettes/web_only/benchmarks.html | 99 +++++++++++++++++++++++++++------ vignettes/web_only/mkin_benchmarks.rda | Bin 1722 -> 1777 bytes 5 files changed, 104 insertions(+), 100 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2a3f231e..83a1a1c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -23,7 +23,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, note that no warranty is implied for correctness of results or fitness for a particular purpose. Depends: R (>= 2.15.1), -Imports: stats, graphics, methods, parallel, deSolve, R6, inline (>= 0.3.19), +Imports: stats, graphics, methods, parallel, deSolve (>= 1.35), R6, inline (>= 0.3.19), numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), saemix (>= 3.1), rlang, vctrs Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, benchmarkme, tibble, stats4, readxl diff --git a/R/mkinfit.R b/R/mkinfit.R index 0d9246dd..6cca5616 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -501,12 +501,10 @@ 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 + 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 @@ -903,9 +901,8 @@ mkinfit <- function(mkinmod, observed, fit$time <- fit_time # 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 + # 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/mkinpredict.R b/R/mkinpredict.R index 11a3b35b..f0d00319 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -24,6 +24,7 @@ #' parent compound. #' @param method.ode The solution method passed via [mkinpredict] to [ode]] in #' 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 the ode solver. @@ -34,7 +35,6 @@ #' 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. #' @return A matrix with the numeric solution in wide format @@ -56,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 @@ -172,79 +172,19 @@ 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 = 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 <- deSolve::lsoda( + y = odeini, + times = outtimes, + func = mkinmod[["symbols"]], + initfunc = "initpar", + dllname = x$dll_info[["name"]], + parms = odeparms[x$parms], # Order matters when using compiled models + atol = atol, + rtol = rtol, + maxsteps = maxsteps, + ... ) - out <- t(out_raw) colnames(out) <- c("time", mod_vars) } else { mkindiff <- function(t, state, parms) { diff --git a/vignettes/web_only/benchmarks.html b/vignettes/web_only/benchmarks.html index 451d8afe..fb799462 100644 --- a/vignettes/web_only/benchmarks.html +++ b/vignettes/web_only/benchmarks.html @@ -1592,7 +1592,7 @@ div.tocify {

Benchmark timings for mkin

Johannes Ranke

-

Last change 14 July 2022 (rebuilt 2023-01-05)

+

Last change 14 July 2022 (rebuilt 2023-02-13)

@@ -1636,11 +1636,31 @@ FOMC_SFO <- mkinmod( DFOP_SFO <- mkinmod( parent = mkinsub("FOMC", "m1"), # erroneously used FOMC twice, not fixed for consistency m1 = mkinsub("SFO")) -t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D)))[["elapsed"]] -t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D), - error_model = "tc"))[["elapsed"]] -t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D), +t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D)))[["elapsed"]] +
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D),
+    error_model = "tc"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D),
     error_model = "obs"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable

Two metabolites, synthetic data:

m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"),
                            M1 = mkinsub("SFO", "M2"),
@@ -1656,18 +1676,36 @@ SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
 
 DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
 
-t6 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a)))[["elapsed"]]
-t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))[["elapsed"]]
-
-t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
-    error_model = "tc"))[["elapsed"]]
-t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
-    error_model = "tc"))[["elapsed"]]
-
-t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
-    error_model = "obs"))[["elapsed"]]
-t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
+t6 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a)))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
+    error_model = "tc"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
+    error_model = "tc"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
+    error_model = "obs"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+
t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
     error_model = "obs"))[["elapsed"]]
+
## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable
+## Error in mkinmod[["symbols"]] : 
+##   object of type 'closure' is not subsettable

Results

@@ -1858,6 +1896,14 @@ models fitted to two datasets, i.e. eight fits for each test.

1.308 1.793 + +Linux +Ryzen 9 7950X 16-Core Processor +4.2.2 +1.3.0 +1.271 +1.787 +
@@ -2068,6 +2114,15 @@ for each test.

2.364 1.230 + +Linux +Ryzen 9 7950X 16-Core Processor +4.2.2 +1.3.0 +0.526 +0.623 +0.631 + @@ -2344,6 +2399,18 @@ dataset, i.e. one fit for each test.

0.801 1.093 + +Linux +Ryzen 9 7950X 16-Core Processor +4.2.2 +1.3.0 +0.242 +0.237 +0.238 +0.239 +0.237 +0.237 + diff --git a/vignettes/web_only/mkin_benchmarks.rda b/vignettes/web_only/mkin_benchmarks.rda index 64ac2680..93138541 100644 Binary files a/vignettes/web_only/mkin_benchmarks.rda and b/vignettes/web_only/mkin_benchmarks.rda differ -- cgit v1.2.1