From b21c601052f85e392e48d903b8903a1a392fe786 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 14 Apr 2015 19:50:57 +0200 Subject: Compile differential equation models with ccSolve package If the ccSolve package is available, and time is not in the right hand side of the equations (i.e. if only SFO and SFORB models are used), the differential equation model is compiled from auto-generated C code. Currently, one test (FOCUS 2006 D SFO_SFO) fails --- DESCRIPTION | 6 ++++-- NAMESPACE | 3 ++- R/mkinmod.R | 34 ++++++++++++++++++++++++++++++++++ R/mkinpredict.R | 22 +++++++++++++--------- man/mkinfit.Rd | 5 +++-- 5 files changed, 56 insertions(+), 14 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 03976d1e..0eaf6ebf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: mkin Type: Package Title: Routines for Fitting Kinetic Models with One or More State Variables to Chemical Degradation Data -Version: 0.9-35 +Version: 0.9-36 Date: 2015-02-20 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), @@ -13,11 +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 and a C compiler is available, differential equation + models are solved as 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, RUnit +Suggests: knitr, RUnit, ccSolve License: GPL LazyLoad: yes LazyData: yes diff --git a/NAMESPACE b/NAMESPACE index 4774475e..92f0f271 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,5 +9,6 @@ import( FME, deSolve, minpack.lm, - rootSolve + rootSolve, + ccSolve ) diff --git a/R/mkinmod.R b/R/mkinmod.R index d04f811f..37cd02b4 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -36,6 +36,9 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # 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 + # is satisfied parms <- vector() # }}} @@ -269,6 +272,37 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) model$coefmat <- m }#}}} + # Create a function compiled from C code if possible #{{{ + if (require(ccSolve)) { + diffs.C <- paste(diffs, collapse = ";\n") + diffs.C <- paste0(diffs.C, ";") + for (i in seq_along(diffs)) { + obs_var <- names(diffs)[i] + + # 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) + + # 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) + } + # Unfortunately, the models with time occuring in the diffs do not compile + if (sum(sapply(spec, function(x) x$type %in% c("SFO", "SFORB"))) == length(spec)) { + message("Compiling differential equation model from auto-generated C code...") + model$compiled <- compile.ode(diffs.C, language = "C", parms = parms) + } + } + # }}} + class(model) <- "mkinmod" return(model) } diff --git a/R/mkinpredict.R b/R/mkinpredict.R index da013d50..bed82441 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -87,17 +87,21 @@ mkinpredict <- function(mkinmod, odeparms, odeini, names(out) <- c("time", mod_vars) } if (solution_type == "deSolve") { - mkindiff <- function(t, state, parms) { + if (!is.null(mkinmod$compiled)) { + mkindiff <- mkinmod$compiled + } else { + 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]]))) + 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))) } - return(list(c(diffs))) } out <- ode( y = odeini, diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 9c82f5ff..466f587c 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -235,7 +235,8 @@ SFO_SFO <- mkinmod( parent = list(type = "SFO", to = "m1", sink = TRUE), m1 = list(type = "SFO")) # Fit the model to the FOCUS example dataset D using defaults -fit <- mkinfit(SFO_SFO, FOCUS_2006_D) +fit <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen") +fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve") summary(fit) # Use stepwise fitting, using optimised parameters from parent only fit, FOMC @@ -258,10 +259,10 @@ SFORB_SFO <- mkinmod( m1 = list(type = "SFO")) # Fit the model to the FOCUS example dataset D using defaults fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D) +fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve") # Use starting parameters from parent only SFORB fit (not really needed in this case) fit.SFORB = mkinfit(SFORB, FOCUS_2006_D) fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode) -} # Weighted fits, including IRLS SFO_SFO.ff <- mkinmod(parent = list(type = "SFO", to = "m1"), -- cgit v1.2.1