diff options
| -rw-r--r-- | DESCRIPTION | 9 | ||||
| -rw-r--r-- | GNUmakefile | 6 | ||||
| -rw-r--r-- | NAMESPACE | 3 | ||||
| -rw-r--r-- | NEWS.md | 2 | ||||
| -rw-r--r-- | R/IORE.solution.R | 4 | ||||
| -rw-r--r-- | R/endpoints.R | 2 | ||||
| -rw-r--r-- | R/mkinerrmin.R | 6 | ||||
| -rw-r--r-- | R/mkinfit.R | 6 | ||||
| -rw-r--r-- | R/mkinmod.R | 93 | ||||
| -rw-r--r-- | R/mkinpredict.R | 43 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 12 | ||||
| -rw-r--r-- | README.Rmd | 10 | ||||
| -rw-r--r-- | README.md | 10 | ||||
| -rw-r--r-- | man/IORE.solution.Rd | 4 | ||||
| -rw-r--r-- | man/mkinfit.Rd | 2 | ||||
| -rw-r--r-- | vignettes/FOCUS_D.Rmd | 2 | ||||
| -rw-r--r-- | vignettes/compiled_models.Rmd | 12 |
17 files changed, 137 insertions, 89 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 59b653f4..b5045e4e 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-36 -Date: 2015-05-18 +Date: 2015-06-20 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), @@ -13,14 +13,13 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006). Includes a function for conveniently defining differential equation models, model solution based on eigenvalues if possible or using numerical solvers and a choice of the optimisation methods made available by the 'FME' package. - If the 'ccSolve' package (https://github.com/karlines/ccSolve), - and a C compiler (on windows: 'Rtools') are installed, differential + If a C compiler (on windows: 'Rtools') are installed, differential equation models are solved using compiled C functions. Please note that no warranty is implied for correctness of results or fitness for a particular purpose. -Depends: minpack.lm, rootSolve +Depends: minpack.lm, rootSolve, inline Imports: FME, deSolve -Suggests: knitr, testthat, microbenchmark, ccSolve +Suggests: knitr, testthat, microbenchmark License: GPL LazyLoad: yes LazyData: yes diff --git a/GNUmakefile b/GNUmakefile index 9e2ee80c..7135d6ee 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -26,7 +26,11 @@ pkgfiles = NEWS \ tests/* \ tests/testthat* \ TODO \ - vignettes/* + vignettes/*.Rmd \ + vignettes/*.Rnw + + + all: check clean @@ -9,5 +9,6 @@ import( FME, deSolve, minpack.lm, - rootSolve + rootSolve, + inline ) @@ -2,7 +2,7 @@ ## MAJOR CHANGES -- If the ccSolve package is installed, use a version of the differential equation model compiled from C code, which is a huge performance boost for models where only the deSolve method works. Compiled models using the Hockey Stick model (HS) are not supported (yet). +- If a compiler (gcc) is installed, use a version of the differential equation model compiled from C code, which is a huge performance boost for models where only the deSolve method works. - `mkinmod()`: Create a list component $compiled in the list returned by mkinmod, if a version can be compiled from autogenerated C code (see above). - `mkinfit()`: Set the default `solution_type` to `deSolve` when a compiled version of the model is present, except if an analytical solution is possible. diff --git a/R/IORE.solution.R b/R/IORE.solution.R index 9546ce56..5405be96 100644 --- a/R/IORE.solution.R +++ b/R/IORE.solution.R @@ -1,4 +1,4 @@ -IORE.solution <- function(t, parent.0, k.iore, N)
+IORE.solution <- function(t, parent.0, k__iore, N)
{
- parent = (parent.0^(1 - N) - (1 - N) * k.iore * t)^(1/(1 - N))
+ parent = (parent.0^(1 - N) - (1 - N) * k__iore * t)^(1/(1 - N))
}
diff --git a/R/endpoints.R b/R/endpoints.R index ae9aa3ec..001595bc 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -50,7 +50,7 @@ endpoints <- function(fit) { ep$distimes[obs_var, c("DT50back")] = DT50_back
}
if (type == "IORE") {
- k_names = grep(paste("k.iore", obs_var, sep="_"), names(parms.all), value=TRUE)
+ k_names = grep(paste("k__iore", obs_var, sep="_"), names(parms.all), value=TRUE)
k_tot = sum(parms.all[k_names])
# From the NAFTA kinetics guidance, p. 5
n = parms.all[paste("N", obs_var, sep = "_")]
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index b28235a6..6103dd1d 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -61,8 +61,8 @@ mkinerrmin <- function(fit, alpha = 0.05) n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim))) n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), names(parms.optim))) - n.k.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim))) - n.k.iore.optim <- n.k.iore.optim + length(grep(paste("^log_k.iore", obs_var, + n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim))) + n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, sep = "_"), names(parms.optim))) @@ -82,7 +82,7 @@ mkinerrmin <- function(fit, alpha = 0.05) } } - n.optim <- sum(n.initials.optim, n.k.optim, n.k.iore.optim, n.N.optim, n.ff.optim) + n.optim <- sum(n.initials.optim, n.k.optim, n.k__iore.optim, n.N.optim, n.ff.optim) # FOMC, DFOP and HS parameters are only counted if we are looking at the # first variable in the model which is always the source variable diff --git a/R/mkinfit.R b/R/mkinfit.R index 5118519a..7f9a55d9 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -219,7 +219,7 @@ mkinfit <- function(mkinmod, observed, if (length(mkinmod$spec) == 1) {
solution_type = "analytical"
} else {
- if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) {
+ if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {
solution_type = "deSolve"
} else {
if (is.matrix(mkinmod$coefmat)) {
@@ -319,8 +319,8 @@ mkinfit <- function(mkinmod, observed, if (!transform_rates) {
index_k <- grep("^k_", names(lower))
lower[index_k] <- 0
- index_k.iore <- grep("^k.iore_", names(lower))
- lower[index_k.iore] <- 0
+ index_k__iore <- grep("^k__iore_", names(lower))
+ lower[index_k__iore] <- 0
other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))
lower[other_rate_parms] <- 0
}
diff --git a/R/mkinmod.R b/R/mkinmod.R index e865fa27..a0307003 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -36,8 +36,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) # differential equations (if supported), parameter names and a mapping from
# model variables to observed variables. If possible, a matrix representation
# of the differential equations is included
- # Compiling the functions from the C code generated below using the ccSolve package
- # only works if the implicit assumption about differential equations specified below
+ # Compiling the functions from the C code generated below only works if the
+ # implicit assumption about differential equations specified below
# is satisfied
parms <- vector()
# }}}
@@ -99,7 +99,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) # If sink is required, add first-order/IORE sink term
k_compound_sink <- paste("k", box_1, "sink", sep = "_")
if(spec[[varname]]$type == "IORE") {
- k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_")
+ k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")
}
parms <- c(parms, k_compound_sink)
decline_term <- paste(k_compound_sink, "*", box_1)
@@ -114,7 +114,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) } else {
k_compound <- paste("k", box_1, sep = "_")
if(spec[[varname]]$type == "IORE") {
- k_compound <- paste("k.iore", box_1, sep = "_")
+ k_compound <- paste("k__iore", box_1, sep = "_")
}
parms <- c(parms, k_compound)
decline_term <- paste(k_compound, "*", box_1)
@@ -135,9 +135,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)
parms <- c(parms, "k1", "k2", "g")
} #}}}
+ HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C
if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
# From p. 55 of the FOCUS kinetics report
- decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)
+ decline_term <- paste(HS_decline, "*", box_1)
parms <- c(parms, "k1", "k2", "tb")
} #}}}
# Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
@@ -272,35 +273,65 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) model$coefmat <- m
}#}}}
- # Create a function compiled from C code if more than one observed variable and ccSolve is available #{{{
- if (length(obs_vars) > 1 & requireNamespace("ccSolve", quietly = TRUE)) {
- diffs.C <- paste(diffs, collapse = ";\n")
- diffs.C <- paste0(diffs.C, ";")
- for (i in seq_along(diffs)) {
- obs_var <- names(diffs)[i]
+ # Create a function compiled from C code if more than one observed variable and gcc is available #{{{
+ if (length(obs_vars) > 1) {
+ if (Sys.which("gcc") != "") {
- # Replace d_... terms by f[i-1]
- # First line
- pattern <- paste0("^d_", obs_var)
- replacement <- paste0("\nf[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
- # Other lines
- pattern <- paste0("\\nd_", obs_var)
- replacement <- paste0("\nf[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
+ diffs.C <- paste(diffs, collapse = ";\n")
+ diffs.C <- paste0(diffs.C, ";")
+
+ # HS
+ diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)
+
+ for (i in seq_along(diffs)) {
+ state_var <- names(diffs)[i]
+
+ # IORE
+ if (state_var %in% obs_vars) {
+ if (spec[[state_var]]$type == "IORE") {
+ diffs.C <- gsub(paste0(state_var, "^N_", state_var),
+ paste0("pow(y[", i - 1, "], N_", state_var, ")"),
+ diffs.C, fixed = TRUE)
+ }
+ }
+
+ # Replace d_... terms by f[i-1]
+ # First line
+ pattern <- paste0("^d_", state_var)
+ replacement <- paste0("\nf[", i - 1, "]")
+ diffs.C <- gsub(pattern, replacement, diffs.C)
+ # Other lines
+ pattern <- paste0("\\nd_", state_var)
+ replacement <- paste0("\nf[", i - 1, "]")
+ diffs.C <- gsub(pattern, replacement, diffs.C)
+
+ # Replace names of observed variables by y[i],
+ # making the implicit assumption that the observed variables only occur after "* "
+ pattern <- paste0("\\* ", state_var)
+ replacement <- paste0("* y[", i - 1, "]")
+ diffs.C <- gsub(pattern, replacement, diffs.C)
+ }
- # Replace names of observed variables by y[i],
- # making the implicit assumption that the observed variables only occur after "* "
- pattern <- paste0("\\* ", obs_var)
- replacement <- paste0("* y[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
- }
- if (sum(sapply(spec, function(x) x$type %in%
- c("SFO", "FOMC", "DFOP", "SFORB"))) == length(spec)) {
if (!quiet) message("Compiling differential equation model from auto-generated C code...")
- model$compiled <- ccSolve::compile.ode(diffs.C, language = "C",
- parms = parms,
- declaration = "double time = *t;")
+
+ npar <- length(parms)
+
+ initpar_code <- paste0(
+ "static double parms [", npar, "];\n",
+ paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),
+ "\n",
+ "void initpar(void (* odeparms)(int *, double *)) {\n",
+ " int N = ", npar, ";\n",
+ " odeparms(&N, parms);\n",
+ "}\n\n")
+
+ derivs_code <- paste0("double time = *t;\n", diffs.C)
+ derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",
+ f = "numeric", rpar = "numeric", ipar = "integer")
+
+ model$cf <- cfunction(list(func = derivs_sig), derivs_code,
+ otherdefs = initpar_code,
+ convention = ".C", language = "C")
}
}
# }}}
diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 3121f1d7..5f11c35a 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -54,8 +54,8 @@ mkinpredict <- function(mkinmod, odeparms, odeini, IORE = IORE.solution(outtimes, evalparse(parent.name), ifelse(mkinmod$use_of_ff == "min", - evalparse(paste("k.iore", parent.name, "sink", sep="_")), - evalparse(paste("k.iore", parent.name, sep="_"))), + evalparse(paste("k__iore", parent.name, "sink", sep="_")), + evalparse(paste("k__iore", parent.name, sep="_"))), evalparse("N_parent")), DFOP = DFOP.solution(outtimes, evalparse(parent.name), @@ -88,10 +88,21 @@ mkinpredict <- function(mkinmod, odeparms, odeini, names(out) <- c("time", mod_vars) } if (solution_type == "deSolve") { - if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) { - mkindiff <- mkinmod$compiled - } else { - mkindiff <- function(t, state, parms) { + if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) { + out <- ode( + y = odeini, + times = outtimes, + func = "func", + initfunc = "initpar", + dllname = getDynLib(mkinmod$cf)[["name"]], + parms = odeparms[mkinmod$parms], # Order matters when using compiled models + method = method.ode, + atol = atol, + rtol = rtol, + ... + ) + } else { + mkindiff <- function(t, state, parms) { time <- t diffs <- vector() @@ -103,17 +114,17 @@ mkinpredict <- function(mkinmod, odeparms, odeini, } return(list(c(diffs))) } + out <- ode( + y = odeini, + times = outtimes, + func = mkindiff, + parms = odeparms, + method = method.ode, + atol = atol, + rtol = rtol, + ... + ) } - out <- ode( - y = odeini, - times = outtimes, - func = mkindiff, - parms = odeparms[mkinmod$parms], # Order matters when using compiled models - method = method.ode, - atol = atol, - rtol = rtol, - ... - ) if (sum(is.na(out)) > 0) { stop("Differential equations were not integrated for all output times because\n", "NaN values occurred in output from ode()") diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index aa70a0b3..0946ff14 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -33,8 +33,8 @@ transform_odeparms <- function(parms, mkinmod, # Log transformation for rate constants if requested
k <- parms[grep("^k_", names(parms))]
- k.iore <- parms[grep("^k.iore_", names(parms))]
- k <- c(k, k.iore)
+ k__iore <- parms[grep("^k__iore_", names(parms))]
+ k <- c(k, k__iore)
if (length(k) > 0) {
if(transform_rates) {
transparms[paste0("log_", names(k))] <- log(k)
@@ -113,8 +113,8 @@ backtransform_odeparms <- function(transparms, mkinmod, # Exponential transformation for rate constants
if(transform_rates) {
trans_k <- transparms[grep("^log_k_", names(transparms))]
- trans_k.iore <- transparms[grep("^log_k.iore_", names(transparms))]
- trans_k = c(trans_k, trans_k.iore)
+ trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
+ trans_k = c(trans_k, trans_k__iore)
if (length(trans_k) > 0) {
k_names <- gsub("^log_k", "k", names(trans_k))
parms[k_names] <- exp(trans_k)
@@ -122,8 +122,8 @@ backtransform_odeparms <- function(transparms, mkinmod, } else {
trans_k <- transparms[grep("^k_", names(transparms))]
parms[names(trans_k)] <- trans_k
- trans_k.iore <- transparms[grep("^k.iore_", names(transparms))]
- parms[names(trans_k.iore)] <- trans_k.iore
+ trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
+ parms[names(trans_k__iore)] <- trans_k__iore
}
# Do not transform exponents in IORE models
@@ -106,10 +106,12 @@ documentation or the package vignettes referenced from the These have decreasing efficiency, and are automatically chosen by default. * As of mkin 0.9-36, model solution for models with more than one observed - variable is based on the - [`ccSolve`](https://github.com/karlines/ccSolve) package, if installed. - This is even faster than eigenvalue based solution, at least in the example - shown in the [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html) + variable is based on the inline package. This is even faster than eigenvalue + based solution, at least in the example shown in the + [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html). + The autogeneration of C code was + inspired by the [`ccSolve`](https://github.com/karlines/ccSolve) package. Thanks + to Karline Soetaert for her work on that. * Model optimisation with [`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html) internally using the `modFit` function from the `FME` package, @@ -120,10 +120,12 @@ documentation or the package vignettes referenced from the These have decreasing efficiency, and are automatically chosen by default. * As of mkin 0.9-36, model solution for models with more than one observed - variable is based on the - [`ccSolve`](https://github.com/karlines/ccSolve) package, if installed. - This is even faster than eigenvalue based solution, at least in the example - shown in the [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html) + variable is based on the inline package. This is even faster than eigenvalue + based solution, at least in the example shown in the + [vignette `compiled_models`](http://rawgit.com/jranke/mkin/master/vignettes/compiled_models.html). + The autogeneration of C code was + inspired by the [`ccSolve`](https://github.com/karlines/ccSolve) package. Thanks + to Karline Soetaert for her work on that. * Model optimisation with [`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html) internally using the `modFit` function from the `FME` package, diff --git a/man/IORE.solution.Rd b/man/IORE.solution.Rd index 3d827bdc..bb377d3f 100644 --- a/man/IORE.solution.Rd +++ b/man/IORE.solution.Rd @@ -7,12 +7,12 @@ a concentration dependent rate constant.
}
\usage{
- IORE.solution(t, parent.0, k.iore, N)
+ IORE.solution(t, parent.0, k__iore, N)
}
\arguments{
\item{t}{ Time. }
\item{parent.0}{ Starting value for the response variable at time zero. }
- \item{k.iore}{ Rate constant. Note that this depends on the concentration units used. }
+ \item{k__iore}{ Rate constant. Note that this depends on the concentration units used. }
\item{N}{ Exponent describing the nonlinearity of the rate equation }
}
\note{
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index b4923ddd..bf47879c 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -247,7 +247,7 @@ print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, coef(fit) endpoints(fit) \dontrun{ -# deSolve is slower when ccSolve is not installed and set up +# deSolve is slower when the gcc compiler is not available print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve"))) coef(fit.deSolve) diff --git a/vignettes/FOCUS_D.Rmd b/vignettes/FOCUS_D.Rmd index 902d3d24..8ae73c16 100644 --- a/vignettes/FOCUS_D.Rmd +++ b/vignettes/FOCUS_D.Rmd @@ -27,7 +27,7 @@ kinetics (SFO) to one metabolite named m1, which also degrades with SFO kinetics The call to mkinmod returns a degradation model. The differential equations represented in
R code can be found in the character vector `$diffs` of the `mkinmod` object. If
-the `ccSolve` package is installed and functional, the differential equation model
+the gcc compiler is installed and functional, the differential equation model
will be compiled from auto-generated C code.
diff --git a/vignettes/compiled_models.Rmd b/vignettes/compiled_models.Rmd index bac284c5..5d83cf3a 100644 --- a/vignettes/compiled_models.Rmd +++ b/vignettes/compiled_models.Rmd @@ -15,16 +15,14 @@ output: ```{r, include = FALSE}
library(knitr)
opts_chunk$set(tidy = FALSE, cache = TRUE)
-if (!require("ccSolve"))
- message("Please install the ccSolve package for this vignette to produce sensible output")
-
```
# 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
-greater than 0.9-36 and the ccSolve package is installed and functional, you will get a
-message that the model is being compiled when defining a model using mkinmod.
+equal to or greater than 0.9-36 and a compiler (gcc) is installed, you will see
+a message that the model is being compiled from autogenerated C code when
+defining a model using mkinmod.
```{r create_SFO_SFO}
library("mkin")
@@ -83,5 +81,5 @@ smb.2["median"]/smb.2["deSolve, compiled", "median"] ```
-Here we get a performance benefit of more than a factor of 8 using the version
-of the differential equation model compiled from C code using the ccSolve package!
+Here we get a performance benefit of more than a factor of 10 using the version
+of the differential equation model compiled from C code using the inline package!
|