From d113cd79b178fdc91aecb894707ed356129dfb75 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sun, 10 May 2020 21:53:00 +0200 Subject: Default to analytical for coupled models if available This revealed that transforming rates is necessary for fitting the analytical solution of the SFO-SFO model to the FOCUS D dataset. Benchmarks show that fitting coupled models with deSolve got a bit slower through the latest changes --- R/create_deg_func.R | 3 +- R/mkinfit.R | 4 +- test.log | 20 +-- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/test_FOCUS_D_UBA_expertise.R | 21 ++- vignettes/mkin_benchmarks.rda | Bin 871 -> 871 bytes vignettes/web_only/benchmarks.Rmd | 14 +- vignettes/web_only/benchmarks.html | 243 +++++++++++++--------------- 8 files changed, 150 insertions(+), 157 deletions(-) diff --git a/R/create_deg_func.R b/R/create_deg_func.R index 886c5e8b..6c0ae40b 100644 --- a/R/create_deg_func.R +++ b/R/create_deg_func.R @@ -51,7 +51,8 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { n10 <- paste0("odeini['", parent, "']") n20 <- paste0("odeini['", n2, "']") - if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2, spec[[2]]$type == "SFO")) { + if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2, + spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) { supported <- TRUE k1 <- paste0("k_", n1) k2 <- paste0("k_", n2) diff --git a/R/mkinfit.R b/R/mkinfit.R index c33afbf0..e1089673 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -464,7 +464,7 @@ mkinfit <- function(mkinmod, observed, if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat)) stop("Eigenvalue based solution not possible, coefficient matrix not present.") if (solution_type == "auto") { - if (length(mkinmod$spec) == 1) { + if (length(mkinmod$spec) == 1 || is.function(mkinmod$deg_func)) { solution_type = "analytical" } else { if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) { @@ -694,7 +694,7 @@ mkinfit <- function(mkinmod, observed, # Show parameter names if tracing is requested if(trace_parms) cat(names_optim, "\n") - # browser() + #browser() # Do the fit and take the time until the hessians are calculated fit_time <- system.time({ diff --git a/test.log b/test.log index 7788f109..2a116058 100644 --- a/test.log +++ b/test.log @@ -2,33 +2,33 @@ Loading mkin Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 s] -✔ | 4 | Calculation of FOCUS chi2 error levels [2.0 s] -✔ | 7 | Fitting the SFORB model [9.1 s] -✔ | 1 | Analytical solutions for coupled models [1.1 s] +✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [2.5 s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.4 s] +✔ | 7 | Fitting the SFORB model [9.0 s] +✔ | 1 | Analytical solutions for coupled models [0.2 s] ✔ | 5 | Calculation of Akaike weights ✔ | 10 | Confidence intervals and p-values [1.0 s] -✔ | 14 | Error model fitting [3.7 s] +✔ | 14 | Error model fitting [3.8 s] ✔ | 6 | Test fitting the decline of metabolites from their maximum [0.2 s] ✔ | 1 | Fitting the logistic model [0.2 s] ✔ | 1 | Test dataset class mkinds used in gmkin ✔ | 12 | Special cases of mkinfit calls [0.6 s] ✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.4 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.4 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.3 s] ✔ | 9 | Nonlinear mixed-effects models [12.1 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] ✔ | 3 | Summary -✔ | 14 | Plotting [3.5 s] +✔ | 14 | Plotting [3.6 s] ✔ | 4 | AIC calculation ✔ | 4 | Residuals extracted from mkinfit models ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [4.2 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.4 s] -✔ | 9 | Hypothesis tests [17.8 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.3 s] +✔ | 9 | Hypothesis tests [17.7 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 69.9 s +Duration: 66.3 s OK: 157 Failed: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index f113a668..d41f0a9f 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2020-05-09 +Date: 2020-05-10 Optimiser: IRLS [Data] diff --git a/tests/testthat/test_FOCUS_D_UBA_expertise.R b/tests/testthat/test_FOCUS_D_UBA_expertise.R index f3b12bd9..be2806de 100644 --- a/tests/testthat/test_FOCUS_D_UBA_expertise.R +++ b/tests/testthat/test_FOCUS_D_UBA_expertise.R @@ -2,19 +2,22 @@ context("Results for FOCUS D established in expertise for UBA (Ranke 2014)") # Results are from p. 40 +# Avoid warnings due to the zero value in the data for m1 at time zero +FOCUS_D <- subset(FOCUS_2006_D, value != 0) + test_that("Fits without formation fractions are correct for FOCUS D", { - fit.default <- expect_warning(mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE), "value of zero") + fit.noff <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE) - expect_equal(round(as.numeric(endpoints(fit.default)$distimes["parent", ]), 2), + expect_equal(round(as.numeric(endpoints(fit.noff)$distimes["parent", ]), 2), c(7.02, 23.33)) - expect_equal(round(as.numeric(endpoints(fit.default)$distimes["m1", ]), 1), + expect_equal(round(as.numeric(endpoints(fit.noff)$distimes["m1", ]), 1), c(131.8, 437.7)) }) test_that("Fits with formation fractions are correct for FOCUS D", { skip_on_cran() - fit.ff <- expect_warning(mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE), "value of zero") + fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_D, quiet = TRUE) expect_equivalent(round(fit.ff$bparms.optim, c(2, 4, 4, 4)), c(99.60, 0.0987, 0.0053, 0.5145)) @@ -29,9 +32,16 @@ test_that("Fits with formation fractions are correct for FOCUS D", { test_that("Fits without internal transformations are correct for FOCUS D", { skip_on_cran() + expect_warning( + # Analytical solutions (now the default for SFO_SFO with formation fractions) + # return NA at some point not internally transforming rates + expect_error(fit.ff.notrans <- mkinfit(SFO_SFO.ff, FOCUS_D, + transform_fractions = FALSE, transform_rates = FALSE, quiet = TRUE))) + expect_warning( fit.ff.notrans <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, - transform_fractions = FALSE, transform_rates = FALSE, quiet = TRUE), + transform_fractions = FALSE, transform_rates = FALSE, + quiet = TRUE, solution_type = "deSolve"), "sum of formation fractions") expect_equivalent(round(fit.ff.notrans$bparms.optim, c(2, 4, 4, 4)), @@ -40,7 +50,6 @@ test_that("Fits without internal transformations are correct for FOCUS D", { expect_equivalent(round(100 * mkinerrmin(fit.ff.notrans)$err.min, 2), c(6.40, 6.46, 4.69)) - expect_equal(round(as.numeric(endpoints(fit.ff.notrans)$distimes["parent", ]), 2), c(7.02, 23.33)) expect_equal(round(as.numeric(endpoints(fit.ff.notrans)$distimes["m1", ]), 1), diff --git a/vignettes/mkin_benchmarks.rda b/vignettes/mkin_benchmarks.rda index 1190d924..32223710 100644 Binary files a/vignettes/mkin_benchmarks.rda and b/vignettes/mkin_benchmarks.rda differ diff --git a/vignettes/web_only/benchmarks.Rmd b/vignettes/web_only/benchmarks.Rmd index 4849635f..cc28735a 100644 --- a/vignettes/web_only/benchmarks.Rmd +++ b/vignettes/web_only/benchmarks.Rmd @@ -40,9 +40,11 @@ if (mkin_version > "0.9.48.1") { ``` ```{r timings} +FOCUS_C <- FOCUS_2006_C +FOCUS_D <- subset(FOCUS_2006_D, value != 0) # Parent only -t1 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_2006_C, FOCUS_2006_D)))[["elapsed"]] -t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_2006_C, FOCUS_2006_D), error_model = "tc"))[["elapsed"]] +t1 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_C, FOCUS_D)))[["elapsed"]] +t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_C, FOCUS_D), error_model = "tc"))[["elapsed"]] # One metabolite SFO_SFO <- mkinmod( @@ -54,9 +56,9 @@ FOMC_SFO <- mkinmod( DFOP_SFO <- mkinmod( parent = mkinsub("FOMC", "m1"), m1 = mkinsub("SFO")) -t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D)))[["elapsed"]] -t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(subset(FOCUS_2006_D, value != 0)), error_model = "tc"))[["elapsed"]] -t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D), error_model = "obs"))[["elapsed"]] +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), error_model = "obs"))[["elapsed"]] # Two metabolites, synthetic data m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"), @@ -83,7 +85,7 @@ t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), error_mod t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "obs"))["elapsed"] mkin_benchmarks[system_string, paste0("t", 1:11)] <- c(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11) -mkin_benchmarks +mkin_benchmarks[, -c(1:3)] save(mkin_benchmarks, file = "~/git/mkin/vignettes/mkin_benchmarks.rda") diff --git a/vignettes/web_only/benchmarks.html b/vignettes/web_only/benchmarks.html index 054b4afc..c80a7e44 100644 --- a/vignettes/web_only/benchmarks.html +++ b/vignettes/web_only/benchmarks.html @@ -1,17 +1,17 @@ - + - + - + Benchmark timings for mkin on various systems @@ -69,8 +69,6 @@ overflow: auto; margin-left: 2%; position: fixed; border: 1px solid #ccc; -webkit-border-radius: 6px; -moz-border-radius: 6px; border-radius: 6px; } @@ -98,10 +96,15 @@ font-size: 12px; .tocify-subheader .tocify-subheader { text-indent: 30px; } - .tocify-subheader .tocify-subheader .tocify-subheader { text-indent: 40px; } +.tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader { +text-indent: 50px; +} +.tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader { +text-indent: 60px; +} .tocify .tocify-item > a, .tocify .nav-list .nav-header { margin: 0px; @@ -504,13 +507,13 @@ float: none; item.append($("", { - "text": self.text() + "html": self.html() })); } else { - item.text(self.text()); + item.html(self.html()); } @@ -1285,7 +1288,7 @@ window.initializeCodeFolding = function(show) { // create a collapsable div to wrap the code in var div = $('
'); - if (show) + if (show || $(this)[0].classList.contains('fold-show')) div.addClass('in'); var id = 'rcode-643E0F36' + currentIndex++; div.attr('id', id); @@ -1401,7 +1404,6 @@ code { } img { max-width:100%; - height: auto; } .tabbed-pane { padding-top: 12px; @@ -1463,6 +1465,7 @@ summary { border: none; display: inline-block; border-radius: 4px; + background-color: transparent; } .tabset-dropdown > .nav-tabs.nav-tabs-open > li { @@ -1475,56 +1478,12 @@ summary { } - - - - - - -