From 438a889c37ffdf8f0c6585092da6abdb63b4575e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 30 Jun 2015 14:17:13 +0200 Subject: Improve calculation of DT90 in endpoints(), check and test --- DESCRIPTION | 4 +- GNUmakefile | 4 ++ NEWS.md | 12 ++++- R/endpoints.R | 38 +++++++++------- check.log | 10 +++-- test.log | 1 + vignettes/FOCUS_Z.pdf | Bin 225060 -> 219834 bytes vignettes/compiled_models.Rmd | 101 ------------------------------------------ 8 files changed, 46 insertions(+), 124 deletions(-) delete mode 100644 vignettes/compiled_models.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index a8cfb15b..85846fa0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,7 @@ Type: Package Title: Routines for Fitting Kinetic Models with One or More State Variables to Chemical Degradation Data Version: 0.9-40 -Date: 2015-06-29 +Date: 2015-06-30 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), @@ -18,7 +18,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006). Please note that no warranty is implied for correctness of results or fitness for a particular purpose. Depends: minpack.lm, rootSolve, inline, parallel -Imports: FME, deSolve +Imports: stats, graphics, FME, deSolve Suggests: knitr, testthat, microbenchmark License: GPL LazyLoad: yes diff --git a/GNUmakefile b/GNUmakefile index 80b88eb4..b6bec285 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -61,8 +61,12 @@ quickcheck: build-no-vignettes clean: $(RM) -r $(PKGNAME).Rcheck/ + $(RM) -r vignettes/*.bbl + $(RM) -r vignettes/*.blg $(RM) -r vignettes/*.fls $(RM) -r vignettes/*.fdb_latexmk + $(RM) -r vignettes/cache + $(RM) -r vignettes/files $(RM) -r vignettes/*_cache $(RM) -r vignettes/*_files $(RM) -r vignettes/*-concordance.tex diff --git a/NEWS.md b/NEWS.md index 3a8900ac..c2dec7ad 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,16 @@ # NEWS for package 'mkin' -## CHANGES in mkin VERSION 0.9-39 +## CHANGES in mkin VERSION 0.9-40 + +### BUG FIXES + +- `endpoints()`: For DFOP and SFORB models, where `optimize()` is used, make use of the fact that the DT50 must be between DT50_k1 and DT50_k2 (DFOP) or DT50_b1 and DT50_b2 (SFORB), as `optimize()` sometimes did not find the minimum. Likewise for finding DT90 values. Also fit on the log scale to make the function more efficient. + +### INTERNAL CHANGES + +- `DESCRIPTION`, `NAMESPACE`, `R/*.R`: Import stats and graphics package, qualify function calls for non-base packages installed with R to avoid NOTES made by R CMD check --as-cran with upcoming R versions. + +## CHANGES in mkin VERSION 0.9-39 (2015-06-26) ### MAJOR CHANGES diff --git a/R/endpoints.R b/R/endpoints.R index 001595bc..a930acb7 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -75,17 +75,20 @@ endpoints <- function(fit) { k1 = parms.all["k1"] k2 = parms.all["k2"] g = parms.all["g"] - f <- function(t, x) { + f <- function(log_t, x) { + t <- exp(log_t) fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) (fraction - (1 - x/100))^2 } - DTmax <- 1000 - DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum - DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o) - DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum - DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o) DT50_k1 = log(2)/k1 DT50_k2 = log(2)/k2 + DT90_k1 = log(10)/k1 + DT90_k2 = log(10)/k2 + DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum)) + DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum)) + if (inherits(DT50, "try-error")) DT50 = NA + if (inherits(DT90, "try-error")) DT90 = NA + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 } @@ -119,26 +122,29 @@ endpoints <- function(fit) { b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + DT50_b1 = log(2)/b1 + DT50_b2 = log(2)/b2 + DT90_b1 = log(10)/b1 + DT90_b2 = log(10)/b2 + + SFORB_fraction = function(t) { ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t) } - f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2 - max_DT <- 1000 - DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum - if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o - f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2 - DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum - if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o + f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2 + DT50 <- try(exp(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum)) + f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2 + DT90 <- try(exp(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum)) + + if (inherits(DT50, "try-error")) DT50 = NA + if (inherits(DT90, "try-error")) DT90 = NA for (k_out_name in k_out_names) { ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output } - DT50_b1 = log(2)/b1 - DT50_b2 = log(2)/b2 - # Return the eigenvalues for comparison with DFOP rate constants ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1 ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2 diff --git a/check.log b/check.log index ef7519fd..385ee730 100644 --- a/check.log +++ b/check.log @@ -10,7 +10,7 @@ * checking CRAN incoming feasibility ... NOTE Maintainer: ‘Johannes Ranke ’ -Days since last update: 3 +Days since last update: 4 * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK @@ -56,7 +56,10 @@ Days since last update: 3 * checking data for ASCII and uncompressed saves ... OK * checking sizes of PDF files under ‘inst/doc’ ... OK * checking installed files from ‘inst/doc’ ... OK -* checking files in ‘vignettes’ ... OK +* checking files in ‘vignettes’ ... NOTE +The following files look like leftovers/mistakes: + ‘FOCUS_Z.bbl’, ‘FOCUS_Z.blg’ +Please remove them from your package. * checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED @@ -65,7 +68,6 @@ Days since last update: 3 * checking running R code from vignettes ... ‘FOCUS_D.Rmd’ using ‘UTF-8’ ... OK ‘FOCUS_L.Rmd’ using ‘UTF-8’ ... OK - ‘compiled_models.Rmd’ using ‘UTF-8’ ... OK ‘FOCUS_Z.Rnw’ using ‘UTF-8’ ... OK ‘mkin.Rnw’ using ‘UTF-8’ ... OK OK @@ -73,7 +75,7 @@ Days since last update: 3 * checking PDF version of manual ... OK * DONE -Status: 1 NOTE +Status: 2 NOTEs See ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ for details. diff --git a/test.log b/test.log index 1449fdb5..dfabfdcd 100644 --- a/test.log +++ b/test.log @@ -11,3 +11,4 @@ Fitting of parent only models : ..................... Complex test case from Schaefer et al. (2007) Piacenza paper : .. Results for synthetic data established in expertise for UBA (Ranke 2014) : .... +Way to go! diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf index 8b8f7416..3fb4f69c 100644 Binary files a/vignettes/FOCUS_Z.pdf and b/vignettes/FOCUS_Z.pdf differ diff --git a/vignettes/compiled_models.Rmd b/vignettes/compiled_models.Rmd deleted file mode 100644 index 8dc74692..00000000 --- a/vignettes/compiled_models.Rmd +++ /dev/null @@ -1,101 +0,0 @@ ---- -title: "Performance benefit by using compiled model definitions in mkin" -author: "Johannes Ranke" -date: "`r Sys.Date()`" -output: - html_document: - css: mkin_vignettes.css - toc: true - mathjax: null - theme: united -vignette: > - %\VignetteIndexEntry{Performance benefit by using compiled model definitions in mkin} - %\VignetteEngine{knitr::rmarkdown} - \usepackage[utf8]{inputenc} ---- - -```{r, include = FALSE} -library(knitr) -opts_chunk$set(tidy = FALSE, cache = TRUE) -``` - -## Benchmark for a model that can also be solved with Eigenvalues - -This evaluation is taken from the example section of mkinfit. When using an mkin version -equal to or greater than 0.9-36 and a C compiler (gcc) is available, you will see -a message that the model is being compiled from autogenerated C code when -defining a model using mkinmod. The `mkinmod()` function checks for presence of -the gcc compiler using - -```{r check_gcc} -Sys.which("gcc") -``` -First, we build a simple degradation model for a parent compound with one metabolite. - -```{r create_SFO_SFO} -library("mkin") -SFO_SFO <- mkinmod( - parent = mkinsub("SFO", "m1"), - m1 = mkinsub("SFO")) -``` - -We can compare the performance of the Eigenvalue based solution against the -compiled version and the R implementation of the differential equations using -the microbenchmark package. - - -```{r benchmark_SFO_SFO} -library("microbenchmark") -mb.1 <- microbenchmark( - mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve", use_compiled = FALSE, - quiet = TRUE), - mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen", quiet = TRUE), - mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve", quiet = TRUE), - times = 3, control = list(warmup = 1)) -smb.1 <- summary(mb.1)[-1] -rownames(smb.1) <- c("deSolve, not compiled", "Eigenvalue based", "deSolve, compiled") -print(smb.1) -``` - -We see that using the compiled model is by a factor of -`r round(smb.1["deSolve, not compiled", "median"]/smb.1["deSolve, compiled", "median"], 1)` -faster than using the R version with the default ode solver, and it is even -faster than the Eigenvalue based solution implemented in R which does not need -iterative solution of the ODEs: - -```{r} -smb.1["median"]/smb.1["deSolve, compiled", "median"] -``` - -## Benchmark for a model that can not be solved with Eigenvalues - -This evaluation is also taken from the example section of mkinfit. - -```{r benchmark_FOMC_SFO} -FOMC_SFO <- mkinmod( - parent = mkinsub("FOMC", "m1"), - m1 = mkinsub( "SFO")) - -mb.2 <- microbenchmark( - mkinfit(FOMC_SFO, FOCUS_2006_D, use_compiled = FALSE, quiet = TRUE), - mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE), - times = 3, control = list(warmup = 1)) -smb.2 <- summary(mb.2)[-1] -rownames(smb.2) <- c("deSolve, not compiled", "deSolve, compiled") -print(smb.2) -smb.2["median"]/smb.2["deSolve, compiled", "median"] -``` - -Here we get a performance benefit of a factor of -`r round(smb.2["deSolve, not compiled", "median"]/smb.2["deSolve, compiled", "median"], 1)` -using the version of the differential equation model compiled from C code using -the inline package! - -This vignette was built with mkin `r packageVersion("mkin")` on - -```{r sessionInfo, echo = FALSE} -cat(capture.output(sessionInfo())[1:3], sep = "\n") -if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) { - cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]])) -} -``` -- cgit v1.2.1