diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 32 | ||||
| -rw-r--r-- | R/mkinmod.R | 44 | ||||
| -rw-r--r-- | R/mkinpredict.R | 27 |
3 files changed, 76 insertions, 27 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index b0df6e8c..5118519a 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2014 Johannes Ranke
+# Copyright (C) 2010-2015 Johannes Ranke
# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
# Contact: jranke@uni-bremen.de
# The summary function is an adapted and extended version of summary.modFit
@@ -26,8 +26,9 @@ mkinfit <- function(mkinmod, observed, state.ini = "auto",
fixed_parms = NULL,
fixed_initials = names(mkinmod$diffs)[-1],
- solution_type = "auto",
+ solution_type = c("auto", "analytical", "eigen", "deSolve"),
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(),
@@ -88,7 +89,7 @@ mkinfit <- function(mkinmod, observed, " not used in the model")
}
- # Warn that the sum of formation fractions may exceed one they are not
+ # Warn that the sum of formation fractions may exceed one if they are not
# fitted in the transformed way
if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) {
warning("The sum of formation fractions may exceed one if you do not use ",
@@ -111,7 +112,7 @@ mkinfit <- function(mkinmod, observed, stop("Fixing formation fractions is not supported when using the ilr ",
"transformation.")
}
- }
+ }
# Set initial parameter values, including a small increment (salt)
# to avoid linear dependencies (singular matrix) in Eigenvalue based solutions
@@ -207,8 +208,9 @@ 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
- if (!solution_type %in% c("auto", "analytical", "eigen", "deSolve"))
- stop("solution_type must be auto, analytical, eigen or de Solve")
+ # Prefer deSolve over eigen if a compiled model is present and use_compiled
+ # is not set to FALSE
+ solution_type = match.arg(solution_type)
if (solution_type == "analytical" && length(mkinmod$spec) > 1)
stop("Analytical solution not implemented for models with metabolites.")
if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
@@ -217,13 +219,17 @@ mkinfit <- function(mkinmod, observed, if (length(mkinmod$spec) == 1) {
solution_type = "analytical"
} else {
- if (is.matrix(mkinmod$coefmat)) {
- solution_type = "eigen"
- if (max(observed$value, na.rm = TRUE) < 0.1) {
- stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
- }
- } else {
+ if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) {
solution_type = "deSolve"
+ } else {
+ if (is.matrix(mkinmod$coefmat)) {
+ solution_type = "eigen"
+ if (max(observed$value, na.rm = TRUE) < 0.1) {
+ stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
+ }
+ } else {
+ solution_type = "deSolve"
+ }
}
}
}
@@ -263,6 +269,7 @@ 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, ...)
@@ -282,6 +289,7 @@ 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, ...)
diff --git a/R/mkinmod.R b/R/mkinmod.R index d04f811f..e865fa27 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2014 Johannes Ranke {{{
+# Copyright (C) 2010-2015 Johannes Ranke {{{
# Contact: jranke@uni-bremen.de
# This file is part of the R package mkin
@@ -16,7 +16,7 @@ # You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/> }}}
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
+mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE)
{
if (is.null(speclist)) spec <- list(...)
else spec <- speclist
@@ -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()
# }}}
@@ -123,8 +126,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) }
} #}}}
if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
- # From p. 53 of the FOCUS kinetics report
- decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)
+ # 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)
parms <- c(parms, "alpha", "beta")
} #}}}
if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
@@ -269,6 +272,39 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) 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]
+
+ # 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)
+ }
+ 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;")
+ }
+ }
+ # }}}
+
class(model) <- "mkinmod"
return(model)
}
diff --git a/R/mkinpredict.R b/R/mkinpredict.R index da013d50..3121f1d7 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2014 Johannes Ranke +# Copyright (C) 2010-2015 Johannes Ranke # Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG # Contact: jranke@uni-bremen.de @@ -19,6 +19,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", + use_compiled = "auto", method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, map_output = TRUE, ...) { @@ -87,23 +88,27 @@ 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) & use_compiled[1] != FALSE) { + 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, times = outtimes, func = mkindiff, - parms = odeparms, + parms = odeparms[mkinmod$parms], # Order matters when using compiled models method = method.ode, atol = atol, rtol = rtol, |