From ec574cff822a1238138c0aa69b3d1459bdc3dfa8 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 19 Jun 2015 17:46:11 +0200 Subject: Use odeintr instead of ccSolve for compiling models --- DESCRIPTION | 11 ++-- GNUmakefile | 3 +- NAMESPACE | 3 +- NEWS.md | 6 +- R/IORE.solution.R | 4 +- R/endpoints.R | 2 +- R/mkinerrmin.R | 6 +- R/mkinfit.R | 23 ++++---- R/mkinmod.R | 95 +++++++++++++++++++----------- R/mkinpredict.R | 50 +++++++++------- R/plot.mkinfit.R | 2 +- R/transform_odeparms.R | 12 ++-- README.md | 2 +- man/IORE.solution.Rd | 4 +- man/mkinfit.Rd | 24 ++------ man/mkinmod.Rd | 34 ++++++++++- man/mkinpredict.Rd | 34 ++++++----- tests/testthat/test_mkinpredict_SFO_SFO.R | 27 ++++++--- vignettes/FOCUS_D.Rmd | 2 +- vignettes/FOCUS_D.html | 16 ++--- vignettes/FOCUS_Z.pdf | Bin 215014 -> 223035 bytes vignettes/compiled_models.Rmd | 35 +++++------ vignettes/compiled_models.html | 63 +++++++++----------- 23 files changed, 256 insertions(+), 202 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 59b653f4..e3f0671d 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-19 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 - equation models are solved using compiled C functions. + If a C++ compiler (on windows: 'Rtools') is 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 -Imports: FME, deSolve -Suggests: knitr, testthat, microbenchmark, ccSolve +Imports: FME, odeintr, deSolve +Suggests: knitr, testthat, microbenchmark License: GPL LazyLoad: yes LazyData: yes diff --git a/GNUmakefile b/GNUmakefile index 9e2ee80c..97557f04 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -26,7 +26,8 @@ pkgfiles = NEWS \ tests/* \ tests/testthat* \ TODO \ - vignettes/* + vignettes/*.Rmd \ + vignettes/*.Rnw all: check clean diff --git a/NAMESPACE b/NAMESPACE index 49be428e..1f769331 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,5 +9,6 @@ import( FME, deSolve, minpack.lm, - rootSolve + rootSolve, + odeintr ) diff --git a/NEWS.md b/NEWS.md index 32ca3e45..d4b74f48 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,13 +2,13 @@ ## 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). -- `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. +- `mkinmod()`: If the C++ compiler g++ is installed and eigenvalue based solution is not possible, produce 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. The integration method and its tolerances are set in the call to `mkinmod()` in this case. The compiled model is stored in an environment component $e in the list returned by mkinmod. +- `mkinfit()`: Set the default `solution_type` to `odeintr` when a compiled version of the model is present, except if an analytical solution is possible. ## MINOR CHANGES - Added a simple showcase vignette with an evaluation of FOCUS example dataset D +- Change the naming of the IORE rate constante from k.iore to k__iore to make int compatible with C++ # CHANGES in mkin VERSION 0.9-35 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..ce380243 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -26,9 +26,8 @@ mkinfit <- function(mkinmod, observed, state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - solution_type = c("auto", "analytical", "eigen", "deSolve"), + solution_type = c("auto", "analytical", "eigen", "deSolve", "odeintr"), method.ode = "lsoda", - use_compiled = "auto", method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), @@ -208,8 +207,7 @@ mkinfit <- function(mkinmod, observed, # Decide if the solution of the model can be based on a simple analytical # formula, the spectral decomposition of the matrix (fundamental system) # or a numeric ode solver from the deSolve package - # Prefer deSolve over eigen if a compiled model is present and use_compiled - # is not set to FALSE + # Prefer odeintr if a compiled model is present solution_type = match.arg(solution_type) if (solution_type == "analytical" && length(mkinmod$spec) > 1) stop("Analytical solution not implemented for models with metabolites.") @@ -219,8 +217,11 @@ mkinfit <- function(mkinmod, observed, if (length(mkinmod$spec) == 1) { solution_type = "analytical" } else { - if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) { - solution_type = "deSolve" + if (!is.null(mkinmod$e$m)) { + solution_type = "odeintr" + if (!missing(atol) | !missing(rtol)) { + message("When fitting with solution type odeintr, atol and rtol can only be set in the call to mkinmod()") + } } else { if (is.matrix(mkinmod$coefmat)) { solution_type = "eigen" @@ -269,9 +270,8 @@ mkinfit <- function(mkinmod, observed, out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution_type, - use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol.deSolve = atol, rtol.deSolve = rtol, ...) assign("out_predicted", out, inherits=TRUE) @@ -289,9 +289,8 @@ mkinfit <- function(mkinmod, observed, out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot, solution_type = solution_type, - use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol.deSolve = atol, rtol.deSolve = rtol, ...) plot(0, type="n", xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)), @@ -319,8 +318,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..4991be22 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -16,12 +16,17 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see }}} -mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) +mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, + odeintr_compile = c("auto", "no", "yes"), + odeintr_method = "rk5_i", + odeintr_atol = 1e-8, odeintr_rtol = 1e-10) { if (is.null(speclist)) spec <- list(...) else spec <- speclist obs_vars <- names(spec) + odeintr_compile <- match.arg(odeintr_compile) + # Check if any of the names of the observed variables contains any other for (obs_var in obs_vars) { if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other") @@ -99,7 +104,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 +119,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) @@ -127,17 +132,18 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) } #}}} if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term # From p. 53 of the FOCUS kinetics report, without the power function so it works in C - decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1) + decline_term <- paste("(alpha/beta) * 1/(( time /beta) + 1) *", box_1) parms <- c(parms, "alpha", "beta") } #}}} if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term # From p. 57 of the FOCUS kinetics report - 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) + 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 +278,58 @@ 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] + # Try to create a function compiled from C code if more than one observed variable and g++ can be called #{{{ + if (length(obs_vars) > 1) { + try(gcc_version <- system("g++ -v 2>&1", intern = TRUE)) + if (!inherits(gcc_version, "try-error") & odeintr_compile != "no") { - # 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) + if ((odeintr_compile == "auto" & mat == FALSE) | odeintr_compile == "yes") { + # HS + diffs <- gsub(HS_decline, "(t <= tb ? k1 : k2)", diffs, fixed = TRUE) - # 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;") + diffs.C <- paste(diffs, collapse = ";\n") + diffs.C <- paste0(diffs.C, ";") + + 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(x[", i - 1, "], N_", state_var, ")"), + diffs.C, fixed = TRUE) + } + } + + # Replace d_... terms by dxdt[i-1] + # First line + pattern <- paste0("^d_", state_var) + replacement <- paste0("\ndxdt[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + # Other lines + pattern <- paste0("\\nd_", state_var) + replacement <- paste0("\ndxdt[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + + # Replace names of observed variables by x[i], + # making the implicit assumption that the observed variables only occur after "* " + pattern <- paste0("\\* ", state_var) + replacement <- paste0("* x[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + + # Handle time, which needs to be enclosed by spaces (except for HS which is handled separately) + diffs.C <- gsub(" time ", " t ", diffs.C) + } + if (!quiet) message("Compiling differential equation model from auto-generated C++ code...") + model$e <- new.env() + model$e$code <- try(compile_sys("m", diffs.C, pars = parms, method = odeintr_method, + atol = odeintr_atol, rtol = odeintr_rtol, env = model$e)) + + if (inherits(model$e$code, "try-error")) { + if (!quiet) message("Compilation failed, model solution will be slower than necessary") + } + } } } # }}} diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 3121f1d7..d91e87df 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -18,11 +18,12 @@ # this program. If not, see mkinpredict <- function(mkinmod, odeparms, odeini, - outtimes, solution_type = "deSolve", - use_compiled = "auto", - method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, + outtimes, solution_type = c("deSolve", "analytical", "eigen", "odeintr"), + method.ode = "lsoda", atol.deSolve = 1e-8, rtol.deSolve = 1e-10, map_output = TRUE, ...) { + solution_type = match.arg(solution_type) + # Get the names of the state variables in the model mod_vars <- names(mkinmod$diffs) @@ -54,8 +55,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,30 +89,26 @@ 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) { + mkindiff <- function(t, state, parms) { - time <- t - diffs <- vector() - for (box in names(mkinmod$diffs)) - { - diffname <- paste("d", box, sep="_") - diffs[diffname] <- with(as.list(c(time, state, parms)), - eval(parse(text=mkinmod$diffs[[box]]))) - } - return(list(c(diffs))) + time <- t + diffs <- vector() + for (box in names(mkinmod$diffs)) + { + diffname <- paste("d", box, sep="_") + diffs[diffname] <- with(as.list(c(time, state, parms)), + eval(parse(text=mkinmod$diffs[[box]]))) } - } + return(list(c(diffs))) + } 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, + atol = atol.deSolve, + rtol = rtol.deSolve, ... ) if (sum(is.na(out)) > 0) { @@ -119,6 +116,17 @@ mkinpredict <- function(mkinmod, odeparms, odeini, "NaN values occurred in output from ode()") } } + if (solution_type == "odeintr") { + if (is.null(mkinmod$e$m)) stop("Method odeintr can not be used as no compiled version of the model is available") + odeparms_argstring = "" + for (parname in names(odeparms)) { + odeparms_argstring = paste0(odeparms_argstring, parname, " = ", odeparms[parname], ", ") + } + odeparms_argstring = gsub(", $", "", odeparms_argstring) + with(as.list(odeparms_argstring), eval(parse(text=paste0("mkinmod$e$m_set_params(", odeparms_argstring, ")")))) + out <- mkinmod$e$m_at(init = odeini, times = outtimes) + names(out) <- c("time", names(mkinmod$diffs)) + } if (map_output) { # Output transformation for models with unobserved compartments like SFORB out_mapped <- data.frame(time = out[,"time"]) diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 31746fb8..1523be7d 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -55,7 +55,7 @@ plot.mkinfit <- function(x, fit = x, odeparms <- parms.all[odenames] out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, - solution_type = solution_type, atol = fit$atol, rtol = fit$rtol) + solution_type = solution_type, atol.deSolve = fit$atol, rtol.deSolve = fit$rtol) # Set up the plot if not to be added to an existing plot if (add == FALSE) { 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 diff --git a/README.md b/README.md index 98daa6f9..d60dd9b1 100644 --- a/README.md +++ b/README.md @@ -121,7 +121,7 @@ documentation or the package vignettes referenced from the 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. + [`odeintr`](http://cran.r-project.org/package=odeintr) 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) * Model optimisation with 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..f3ce6a08 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -21,9 +21,8 @@ mkinfit(mkinmod, observed, parms.ini = "auto", state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - solution_type = c("auto", "analytical", "eigen", "deSolve"), + solution_type = c("auto", "analytical", "eigen", "deSolve", "odeintr"), method.ode = "lsoda", - use_compiled = "auto", method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), @@ -103,11 +102,6 @@ mkinfit(mkinmod, observed, \code{\link{ode}} in case the solution type is "deSolve". The default "lsoda" is performant, but sometimes fails to converge. } - \item{use_compiled}{ - If set to \code{FALSE}, no compiled version of the \code{\link{mkinmod}} - model is used, in the calls to \code{\link{mkinpredict}} even if - a compiled verion is present. - } \item{method.modFit}{ The optimisation method passed to \code{\link{modFit}}. @@ -175,12 +169,13 @@ mkinfit(mkinmod, observed, according to the number of observations. } \item{atol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-8, - lower than in \code{\link{lsoda}}. + Absolute error tolerance, passed to \code{\link{ode}} when the solution + type is "deSolve". Default is 1e-8, lower than in \code{\link{lsoda}}. } \item{rtol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-10, - much lower than in \code{\link{lsoda}}. + Absolute error tolerance, passed to \code{\link{ode}} when the solution + type is "deSolve". Default is 1e-10, much lower than in + \code{\link{lsoda}}. } \item{n.outtimes}{ The length of the dataseries that is produced by the model prediction @@ -246,13 +241,6 @@ print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen", quiet = TRUE))) coef(fit) endpoints(fit) -\dontrun{ -# deSolve is slower when ccSolve is not installed and set up -print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, - solution_type = "deSolve"))) -coef(fit.deSolve) -endpoints(fit.deSolve) -} # Use stepwise fitting, using optimised parameters from parent only fit, FOMC \dontrun{ diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd index f6a8fe5d..c5d59403 100644 --- a/man/mkinmod.Rd +++ b/man/mkinmod.Rd @@ -11,9 +11,16 @@ For the definition of model types and their parameters, the equations given in the FOCUS and NAFTA guidance documents are used. + + When a C++ compiler is available (currently checking for g++), the model is + compiled using the \code{\link{odeintr}} package, so that model fitting will + be faster. } \usage{ -mkinmod(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) +mkinmod(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, + odeintr_compile = c("auto", "no", "yes"), + odeintr_method = "rk5_i", + odeintr_atol = 1e-8, odeintr_rtol = 1e-10) } \arguments{ \item{...}{ @@ -45,6 +52,24 @@ mkinmod(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) \item{quiet}{ Should messages be suppressed? } + \item{odeintr_compile}{ + When "no", the model is not compiled using \code{\link{compile_sys}}. + When "auto", it is only compiled when an eigenvalue based solution + is not possible, as compiling the model takes some time and we + want to avoid this in case it does not bring a large performance benefit. + } + \item{odeintr_method}{ + The integration method passed to \code{\link{compile_sys}}. The default "rk5_i" + is an interpolating adaptive method. + } + \item{odeintr_atol}{ + Absolute error tolerance, only used when the model is compiled using the odeintr + package. Default is 1e-8, lower than in \code{\link{compile_sys}}. + } + \item{odeintr_rtol}{ + Relative error tolerance, only used when the model is compiled using the odeintr + package. Default is 1e-10, lower than in \code{\link{compile_sys}}. + } } \value{ A list of class \code{mkinmod} for use with \code{\link{mkinfit}}, containing @@ -56,6 +81,11 @@ mkinmod(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) \item{use_of_ff}{ The content of \code{use_of_ff} is passed on in this list component. } \item{coefmat}{ The coefficient matrix, if the system of differential equations can be represented by one. } + \item{e}{An environment where the functions generated by \{code\{link{compile_sys}} + are installed when the model is compiled. The name used for the model is "m", so + for example the function for setting parameters is e$m_set_params(). See + the help page for \code{\link{compile_sys}} for the available functions + and their arguments. } \note{ The IORE submodel is not well tested (yet). When using this model for metabolites, @@ -92,12 +122,12 @@ SFO_SFO <- mkinmod( # If we have several parallel metabolites # (compare tests/testthat/test_synthetic_data_for_UBA_2014.R) +\dontrun{ m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), M1 = mkinsub("SFO"), M2 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) -\dontrun{ fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par, synthetic_data_for_UBA_2014[[12]]$data, quiet = TRUE) diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd index c6aee75f..745f1922 100644 --- a/man/mkinpredict.Rd +++ b/man/mkinpredict.Rd @@ -9,9 +9,10 @@ kinetic parameters and initial values for the state variables. } \usage{ - mkinpredict(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", - use_compiled = "auto", method.ode = "lsoda", atol = 1e-08, rtol = 1e-10, - map_output = TRUE, ...) + mkinpredict(mkinmod, odeparms, odeini, outtimes, + solution_type = c("deSolve", "analytical", "eigen", "odeintr"), + method.ode = "lsoda", atol.deSolve = 1e-08, rtol.deSolve = 1e-10, + map_output = TRUE, ...) } \arguments{ \item{mkinmod}{ @@ -35,22 +36,21 @@ generally be "analytical" if there is only one observed variable, and usually "deSolve" in the case of several observed variables. The third possibility "eigen" is faster but not applicable to some models e.g. - using FOMC for the parent compound. + using FOMC for the parent compound. If a compiler is installed and functional, + method "odeintr" is also available, which is faster for complex models. + The default is "deSolve", as it works for all models and does not depend + on a compiler to be present. } \item{method.ode}{ The solution method passed via \code{\link{mkinpredict}} to \code{\link{ode}} in case the solution type is "deSolve". The default "lsoda" is performant, but sometimes fails to converge. } - \item{use_compiled}{ - If set to \code{FALSE}, no compiled version of the \code{\link{mkinmod}} - model is used, even if is present. - } - \item{atol}{ + \item{atol.deSolve}{ Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-8, lower than in \code{\link{lsoda}}. } - \item{rtol}{ + \item{rtol.deSolve}{ Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-10, much lower than in \code{\link{lsoda}}. } @@ -75,8 +75,6 @@ solution_type = "analytical") mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, solution_type = "deSolve") - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - solution_type = "deSolve", use_compiled = FALSE) mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, solution_type = "eigen") @@ -99,20 +97,26 @@ mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), seq(0, 20, by = 0.01))[2001,] - # Check compiled model versions - they are faster than the eigenvalue based solutions! SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), m1 = list(type = "SFO")) + \dontrun{ + # Check compiled model versions - they are faster than the eigenvalue based solutions! + SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), + m1 = list(type = "SFO"), odeintr_compile = "yes") + } system.time( print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), solution_type = "eigen")[201,])) + \dontrun{ system.time( print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "deSolve")[201,])) + solution_type = "odeintr")[201,])) system.time( print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "deSolve", use_compiled = FALSE)[201,])) + solution_type = "deSolve")[201,])) + } } \keyword{ manip } diff --git a/tests/testthat/test_mkinpredict_SFO_SFO.R b/tests/testthat/test_mkinpredict_SFO_SFO.R index b7004f64..b6e5f9e2 100644 --- a/tests/testthat/test_mkinpredict_SFO_SFO.R +++ b/tests/testthat/test_mkinpredict_SFO_SFO.R @@ -28,19 +28,22 @@ test_that("Variants of model predictions for SFO_SFO model give equivalent resul # and relative tolerance is thus not met tol = 0.01 SFO_SFO.1 <- mkinmod(parent = list(type = "SFO", to = "m1"), - m1 = list(type = "SFO"), use_of_ff = "min", quiet = TRUE) + m1 = list(type = "SFO"), use_of_ff = "min", quiet = TRUE, odeintr_compile = "yes") SFO_SFO.2 <- mkinmod(parent = list(type = "SFO", to = "m1"), m1 = list(type = "SFO"), use_of_ff = "max", quiet = TRUE) ot = seq(0, 100, by = 1) - r.1.e <- subset(mkinpredict(SFO_SFO.1, - c(k_parent_m1 = 0.1, k_parent_sink = 0.1, k_m1_sink = 0.1), - c(parent = 100, m1 = 0), ot, solution_type = "eigen"), - time %in% c(1, 10, 50, 100)) - r.1.d <- subset(mkinpredict(SFO_SFO.1, - c(k_parent_m1 = 0.1, k_parent_sink = 0.1, k_m1_sink = 0.1), - c(parent = 100, m1 = 0), ot, solution_type = "deSolve"), - time %in% c(1, 10, 50, 100)) + test_parms = c(k_parent_m1 = 0.1, k_parent_sink = 0.1, k_m1_sink = 0.1) + test_ini = c(parent = 100, m1 = 0) + r.1.e <- subset(mkinpredict(SFO_SFO.1, test_parms, test_ini, ot, + solution_type = "eigen"), + time %in% c(1, 10, 50, 100)) + r.1.d <- subset(mkinpredict(SFO_SFO.1, test_parms, test_ini, ot, + solution_type = "deSolve"), + time %in% c(1, 10, 50, 100)) + r.1.o <- subset(mkinpredict(SFO_SFO.1, test_parms, test_ini, ot, + solution_type = "odeintr"), + time %in% c(1, 10, 50, 100)) r.2.e <- subset(mkinpredict(SFO_SFO.2, c(k_parent = 0.2, f_parent_to_m1 = 0.5, k_m1 = 0.1), c(parent = 100, m1 = 0), ot, solution_type = "eigen"), @@ -55,6 +58,12 @@ test_that("Variants of model predictions for SFO_SFO model give equivalent resul dev.1.e_d.percent = ifelse(is.na(dev.1.e_d.percent), 0, dev.1.e_d.percent) expect_equivalent(dev.1.e_d.percent < tol, rep(TRUE, length(dev.1.e_d.percent))) + # Compare eigen and odeintr for minimum use of formation fractions + dev.1.e_o.percent = 100 * (r.1.e[-1] - r.1.o[-1])/r.1.e[-1] + dev.1.e_o.percent = as.numeric(unlist((dev.1.e_o.percent))) + dev.1.e_o.percent = ifelse(is.na(dev.1.e_o.percent), 0, dev.1.e_o.percent) + expect_equivalent(dev.1.e_o.percent < tol, rep(TRUE, length(dev.1.e_o.percent))) + # Compare eigen and deSolve for maximum use of formation fractions dev.2.e_d.percent = 100 * (r.1.e[-1] - r.1.d[-1])/r.1.e[-1] dev.2.e_d.percent = as.numeric(unlist((dev.2.e_d.percent))) diff --git a/vignettes/FOCUS_D.Rmd b/vignettes/FOCUS_D.Rmd index 902d3d24..f3dd6661 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 +a compiler (g++) is installed and functional, the differential equation model will be compiled from auto-generated C code. diff --git a/vignettes/FOCUS_D.html b/vignettes/FOCUS_D.html index 6573cc7a..9522e881 100644 --- a/vignettes/FOCUS_D.html +++ b/vignettes/FOCUS_D.html @@ -276,13 +276,13 @@ 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 +a compiler (g++) is installed and functional, the differential equation model will be compiled from auto-generated C code.

SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
-
## Compiling differential equation model from auto-generated C code...
+
## Compiling differential equation model from auto-generated C++ code...
 

 
 
print(SFO_SFO$diffs)
@@ -312,7 +312,7 @@ using the plot method for mkinfit objects.

 
mkinparplot(fit)
 

 
-
plot of chunk unnamed-chunk-6 

+
plot of chunk unnamed-chunk-6 

 
 
A comprehensive report of the results is obtained using the summary method for mkinfit
 objects.

@@ -322,16 +322,16 @@ objects.

 
 
## mkin version:    0.9.36 
 ## R version:       3.2.0 
-## Date of fit:     Fri Jun  5 14:20:31 2015 
-## Date of summary: Fri Jun  5 14:20:31 2015 
+## Date of fit:     Fri Jun 19 16:21:21 2015 
+## Date of summary: Fri Jun 19 16:21:21 2015 
 ## 
 ## Equations:
 ## d_parent = - k_parent_sink * parent - k_parent_m1 * parent
 ## d_m1 = + k_parent_m1 * parent - k_m1_sink * m1
 ## 
-## Model predictions using solution type deSolve 
+## Model predictions using solution type odeintr 
 ## 
-## Fitted with method Port using 153 model solutions performed in 0.621 s
+## Fitted with method Port using 153 model solutions performed in 0.562 s
 ## 
 ## Weighting: none
 ## 
@@ -370,7 +370,7 @@ objects.

 ## parent_0           1.00000            0.6075        -0.06625       -0.1701
 ## log_k_parent_sink  0.60752            1.0000        -0.08740       -0.6253
 ## log_k_parent_m1   -0.06625           -0.0874         1.00000        0.4716
-## log_k_m1_sink     -0.17006           -0.6253         0.47163        1.0000
+## log_k_m1_sink     -0.17006           -0.6253         0.47164        1.0000
 ## 
 ## Residual standard error: 3.211 on 36 degrees of freedom
 ## 
diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf
index 3174a23a..e2a4baa9 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
index bac284c5..b6d54710 100644
--- a/vignettes/compiled_models.Rmd
+++ b/vignettes/compiled_models.Rmd
@@ -15,22 +15,20 @@ 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.
+greater or equal than 0.9-36 and the C++ compiler g++ is installed and functional (on Windows,
+install Rtools), you will get a message that the model is being compiled when
+defining a model using mkinmod.
 
 ```{r create_SFO_SFO}
 library("mkin")
 SFO_SFO <- mkinmod(
   parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "SFO"))
+  m1 = list(type = "SFO"), odeintr_compile = "yes")
 ```
 
 We can compare the performance of the Eigenvalue based solution against the
@@ -39,28 +37,23 @@ the microbenchmark package.
 
 
 ```{r benchmark_SFO_SFO, echo=-(1:2)}
-# Redefining the model, in order not to confuse the knitr cache which leads to segfaults
-suppressMessages(SFO_SFO <- mkinmod(
-  parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "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),
+  mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen", quiet = TRUE),
+  mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "odeintr", 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")
+rownames(smb.1) <- c("deSolve, not compiled", "Eigenvalue based", "odeintr, compiled")
 print(smb.1)
 ```
 
-We see that using the compiled model is almost a factor of 8 faster than using the R version
+We see that using the compiled model is more than a factor of 7 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"]
+smb.1["median"]/smb.1["odeintr, compiled", "median"]
 ```
 
 # Benchmark for a model that can not be solved with Eigenvalues
@@ -73,15 +66,15 @@ FOMC_SFO <- mkinmod(
   m1 = list(type = "SFO"))
 
 mb.2 <- microbenchmark(
-  mkinfit(FOMC_SFO, FOCUS_2006_D, use_compiled = FALSE, quiet = TRUE),
-  mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE),
+  mkinfit(FOMC_SFO, FOCUS_2006_D, solution_type = "deSolve", quiet = TRUE),
+  mkinfit(FOMC_SFO, FOCUS_2006_D, solution_type = "odeintr", quiet = TRUE),
   times = 3, control = list(warmup = 1))
 smb.2 <- summary(mb.2)[-1]
-rownames(smb.2) <- c("deSolve, not compiled", "deSolve, compiled")
+rownames(smb.2) <- c("deSolve, not compiled", "odeintr, compiled")
 print(smb.2)
-smb.2["median"]/smb.2["deSolve, compiled", "median"]
+smb.2["median"]/smb.2["odeintr, 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!
+of the differential equation model compiled from C++ code using the odeintr package!
diff --git a/vignettes/compiled_models.html b/vignettes/compiled_models.html
index 2f2a6edb..efdbe20d 100644
--- a/vignettes/compiled_models.html
+++ b/vignettes/compiled_models.html
@@ -77,37 +77,30 @@ img {
 -->
 

 
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.

+
This evaluation is taken from the example section of mkinfit. When using an mkin version greater or equal than 0.9-36 and the C++ compiler g++ is installed and functional (on Windows, install Rtools), you will get a message that the model is being compiled when defining a model using mkinmod.

 
library("mkin")
 SFO_SFO <- mkinmod(
   parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "SFO"))

-
## Compiling differential equation model from auto-generated C code...

+  m1 = list(type = "SFO"), odeintr_compile = "yes")

+
## Compiling differential equation model from auto-generated C++ code...

 
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.

-
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")
+
smb.1 <- summary(mb.1)[-1]
+rownames(smb.1) <- c("deSolve, not compiled", "Eigenvalue based", "odeintr, compiled")
 print(smb.1)

 
##                             min        lq      mean    median        uq
-## deSolve, not compiled 6192.0125 6195.3470 6211.0309 6198.6816 6220.5401
-## Eigenvalue based       956.7604 1008.7224 1026.2572 1060.6844 1061.0055
-## deSolve, compiled      869.6880  871.9315  883.4929  874.1751  890.3953
+## deSolve, not compiled 5254.1030 5261.3501 5277.1074 5268.5973 5288.6096
+## Eigenvalue based       897.1575  921.6935  930.9546  946.2296  947.8531
+## odeintr, compiled      693.6001  709.0719  719.5530  724.5438  732.5295
 ##                             max neval
-## deSolve, not compiled 6242.3986     3
-## Eigenvalue based      1061.3266     3
-## deSolve, compiled      906.6155     3

-
We see that using the compiled model is almost a factor of 8 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:

-
smb.1["median"]/smb.1["deSolve, compiled", "median"]

+## deSolve, not compiled 5308.6218     3
+## Eigenvalue based       949.4766     3
+## odeintr, compiled      740.5151     3

+
We see that using the compiled model is more than a factor of 7 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:

+
smb.1["median"]/smb.1["odeintr, compiled", "median"]

 
##                         median
-## deSolve, not compiled 7.120877
-## Eigenvalue based      1.205328
-## deSolve, compiled     1.000000

+## deSolve, not compiled 7.290796
+## Eigenvalue based      1.370242
+## odeintr, compiled     1.000000

 
 

 
Benchmark for a model that can not be solved with Eigenvalues

@@ -115,25 +108,25 @@ print(smb.1)
FOMC_SFO <- mkinmod(
   parent = list(type = "FOMC", to = "m1", sink = TRUE),
   m1 = list(type = "SFO"))
-
## Compiling differential equation model from auto-generated C code...
+
## Compiling differential equation model from auto-generated C++ code...
mb.2 <- microbenchmark(
-  mkinfit(FOMC_SFO, FOCUS_2006_D, use_compiled = FALSE, quiet = TRUE),
-  mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE),
+  mkinfit(FOMC_SFO, FOCUS_2006_D, solution_type = "deSolve", quiet = TRUE),
+  mkinfit(FOMC_SFO, FOCUS_2006_D, solution_type = "odeintr", quiet = TRUE),
   times = 3, control = list(warmup = 1))
 smb.2 <- summary(mb.2)[-1]
-rownames(smb.2) <- c("deSolve, not compiled", "deSolve, compiled")
+rownames(smb.2) <- c("deSolve, not compiled", "odeintr, compiled")
 print(smb.2)
##                             min        lq      mean    median        uq
-## deSolve, not compiled 13.297283 13.427702 13.481155 13.558121 13.573092
-## deSolve, compiled      1.486926  1.526887  1.546851  1.566848  1.576813
+## deSolve, not compiled 11.243675 11.324875 11.382415 11.406074 11.451785
+## odeintr, compiled      1.207114  1.209908  1.239989  1.212703  1.256426
 ##                             max neval
-## deSolve, not compiled 13.588063     3
-## deSolve, compiled      1.586778     3
-
smb.2["median"]/smb.2["deSolve, compiled", "median"]
+## deSolve, not compiled 11.497496 3 +## odeintr, compiled 1.300149 3 +
smb.2["median"]/smb.2["odeintr, compiled", "median"]
##                         median
-## deSolve, not compiled 8.653119
-## deSolve, compiled     1.000000
-

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!

+## deSolve, not compiled 9.405494 +## odeintr, compiled 1.000000 +

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 odeintr package!

-- cgit v1.2.1