Improve calculation of DT90 in endpoints(), check and test

8 files changed, 46 insertions, 124 deletions

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 <jranke@uni-bremen.de>’
 
-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
--- a/vignettes/FOCUS_Z.pdf
+++ b/vignettes/FOCUS_Z.pdf
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]]))

-}

-```