diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2023-02-13 06:23:02 +0100 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2023-02-13 06:23:02 +0100 |
| commit | 19861806ed790dc55380c64c6a5e27ba7ecd52de (patch) | |
| tree | d33be92ab644e81e7ffc0c6a6e7e22af4c48005d | |
| parent | 8d1a84ac2190538ed3bac53a303064e281595868 (diff) | |
WIP adapting to new deSolve with faster lsoda
| -rw-r--r-- | DESCRIPTION | 2 | ||||
| -rw-r--r-- | R/mkinfit.R | 13 | ||||
| -rw-r--r-- | R/mkinpredict.R | 90 | ||||
| -rw-r--r-- | vignettes/web_only/benchmarks.html | 99 | ||||
| -rw-r--r-- | 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 { <h1 class="title toc-ignore">Benchmark timings for mkin</h1> <h4 class="author">Johannes Ranke</h4> -<h4 class="date">Last change 14 July 2022 (rebuilt 2023-01-05)</h4> +<h4 class="date">Last change 14 July 2022 (rebuilt 2023-02-13)</h4> </div> @@ -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"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable +## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D), + error_model = "tc"))[["elapsed"]]</code></pre> +<pre><code>## 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</code></pre> +<pre class="r"><code>t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D), error_model = "obs"))[["elapsed"]]</code></pre> +<pre><code>## 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</code></pre> <p>Two metabolites, synthetic data:</p> <pre class="r"><code>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"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))[["elapsed"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), + error_model = "tc"))[["elapsed"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable +## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), + error_model = "tc"))[["elapsed"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable +## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), + error_model = "obs"))[["elapsed"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable +## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> +<pre class="r"><code>t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "obs"))[["elapsed"]]</code></pre> +<pre><code>## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable +## Error in mkinmod[["symbols"]] : +## object of type 'closure' is not subsettable</code></pre> </div> <div id="results" class="section level2"> <h2>Results</h2> @@ -1858,6 +1896,14 @@ models fitted to two datasets, i.e. eight fits for each test.</p> <td align="right">1.308</td> <td align="right">1.793</td> </tr> +<tr class="odd"> +<td align="left">Linux</td> +<td align="left">Ryzen 9 7950X 16-Core Processor</td> +<td align="left">4.2.2</td> +<td align="left">1.3.0</td> +<td align="right">1.271</td> +<td align="right">1.787</td> +</tr> </tbody> </table> </div> @@ -2068,6 +2114,15 @@ for each test.</p> <td align="right">2.364</td> <td align="right">1.230</td> </tr> +<tr class="odd"> +<td align="left">Linux</td> +<td align="left">Ryzen 9 7950X 16-Core Processor</td> +<td align="left">4.2.2</td> +<td align="left">1.3.0</td> +<td align="right">0.526</td> +<td align="right">0.623</td> +<td align="right">0.631</td> +</tr> </tbody> </table> </div> @@ -2344,6 +2399,18 @@ dataset, i.e. one fit for each test.</p> <td align="right">0.801</td> <td align="right">1.093</td> </tr> +<tr class="odd"> +<td align="left">Linux</td> +<td align="left">Ryzen 9 7950X 16-Core Processor</td> +<td align="left">4.2.2</td> +<td align="left">1.3.0</td> +<td align="right">0.242</td> +<td align="right">0.237</td> +<td align="right">0.238</td> +<td align="right">0.239</td> +<td align="right">0.237</td> +<td align="right">0.237</td> +</tr> </tbody> </table> </div> diff --git a/vignettes/web_only/mkin_benchmarks.rda b/vignettes/web_only/mkin_benchmarks.rda Binary files differindex 64ac2680..93138541 100644 --- a/vignettes/web_only/mkin_benchmarks.rda +++ b/vignettes/web_only/mkin_benchmarks.rda |