diff options
| -rw-r--r-- | DESCRIPTION | 6 | ||||
| -rw-r--r-- | NAMESPACE | 3 | ||||
| -rw-r--r-- | R/mkinmod.R | 34 | ||||
| -rw-r--r-- | R/mkinpredict.R | 22 | ||||
| -rw-r--r-- | 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 @@ -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"), |