aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-05-11 05:15:19 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-05-11 05:18:32 +0200
commit234c9059a95e104917e488a6ddd2313234a96cdc (patch)
treef6e54098f79d94578434ef727b62f7cc5d5e79b7 /R
parentd113cd79b178fdc91aecb894707ed356129dfb75 (diff)
Avoid merge() and data.frame() in cost function
also for deSolve and eigenvalue based solutions. This noticeably increases performance for these methods, see test.log and benchmark vignette.
Diffstat (limited to 'R')
-rw-r--r--R/add_err.R2
-rw-r--r--R/mkin_wide_to_long.R1
-rw-r--r--R/mkinfit.R8
-rw-r--r--R/mkinpredict.R48
-rw-r--r--R/plot.mkinfit.R1
5 files changed, 32 insertions, 28 deletions
diff --git a/R/add_err.R b/R/add_err.R
index a523e9c2..9235223f 100644
--- a/R/add_err.R
+++ b/R/add_err.R
@@ -76,6 +76,8 @@ add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"),
{
if (!is.na(seed)) set.seed(seed)
+ prediction <- as.data.frame(prediction)
+
# The output of mkinpredict is in wide format
d_long = mkin_wide_to_long(prediction, time = "time")
diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R
index bef0e408..971f5273 100644
--- a/R/mkin_wide_to_long.R
+++ b/R/mkin_wide_to_long.R
@@ -21,6 +21,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
#' @export
mkin_wide_to_long <- function(wide_data, time = "t")
{
+ wide_data <- as.data.frame(wide_data)
colnames <- names(wide_data)
if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")
vars <- subset(colnames, colnames != time)
diff --git a/R/mkinfit.R b/R/mkinfit.R
index e1089673..f6691b1b 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -592,12 +592,8 @@ mkinfit <- function(mkinmod, observed,
method.ode = method.ode,
atol = atol, rtol = rtol, ...)
- out_long <- mkin_wide_to_long(out, time = "time")
-
- out_merged <- merge(observed[c("name", "time")], out_long)
- out_merged$name <- ordered(out_merged$name, levels = obs_vars)
- out_merged <- out_merged[order(out_merged$name, out_merged$time), ]
- observed$predicted <- out_merged$value
+ observed_index <- cbind(as.character(observed$time), as.character(observed$name))
+ observed$predicted <- out[observed_index]
}
# Define standard deviation for each observation
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 75582fac..90cd45fb 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -32,12 +32,13 @@
#' is 1e-10, much lower than in \code{\link{lsoda}}.
#' @param map_output Boolean to specify if the output should list values for
#' the observed variables (default) or for all state variables (if set to
-#' FALSE).
+#' FALSE). Setting this to FALSE has no effect for analytical solutions,
+#' as these always return mapped output.
#' @param \dots Further arguments passed to the ode solver in case such a
#' solver is used.
#' @import deSolve
#' @importFrom inline getDynLib
-#' @return A data frame with the solution in wide format
+#' @return A matrix with the numeric solution in wide format
#' @author Johannes Ranke
#' @examples
#'
@@ -70,7 +71,7 @@
#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100),
#' seq(0, 20, by = 0.01))[2001,]
#'
-#' # Check compiled model versions - they are faster than the eigenvalue based solutions!
+#' # Comparison of the performance of solution types
#' SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),
#' m1 = list(type = "SFO"), use_of_ff = "max")
#' if(require(rbenchmark)) {
@@ -92,15 +93,11 @@
#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
#' solution_type = "analytical", use_compiled = FALSE)[201,])
#' }
-#' analytical = mkinpredict(SFO_SFO,
-#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
-#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-#' solution_type = "analytical", use_compiled = FALSE)[201,]
#'
#' \dontrun{
#' # Predict from a fitted model
#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE)
-#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical")
+#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "deSolve")
#' head(mkinpredict(f))
#' }
#'
@@ -127,29 +124,34 @@ mkinpredict.mkinmod <- function(x,
map_output = TRUE, ...)
{
- # Get the names of the state variables in the model
+ # Names of state variables and observed variables
mod_vars <- names(x$diffs)
+ obs_vars <- names(x$spec)
# Order the inital values for state variables if they are named
if (!is.null(names(odeini))) {
odeini <- odeini[mod_vars]
}
+ out_obs <- matrix(NA, nrow = length(outtimes), ncol = 1 + length(obs_vars),
+ dimnames = list(as.character(outtimes), c("time", obs_vars)))
+ out_obs[, "time"] <- outtimes
+
if (solution_type == "analytical") {
# This is clumsy, as we wanted fast analytical predictions for mkinfit
- out <- data.frame(time = outtimes)
- obs_vars <- names(x$spec)
pseudo_observed <-
data.frame(name = rep(obs_vars, each = length(outtimes)),
time = rep(outtimes, length(obs_vars)))
pseudo_observed$predicted <- x$deg_func(pseudo_observed, odeini, odeparms)
for (obs_var in obs_vars) {
- out[obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"]
+ out_obs[, obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"]
}
+ # We don't have solutions for unobserved state variables, the output of
+ # analytical solutions is always mapped to observed variables
+ return(out_obs)
}
if (solution_type == "eigen") {
-
evalparse <- function(string) {
eval(parse(text=string), as.list(c(odeparms, odeini)))
}
@@ -163,8 +165,8 @@ mkinpredict.mkinmod <- function(x,
}
o <- matrix(mapply(f.out, outtimes),
nrow = length(mod_vars), ncol = length(outtimes))
- out <- data.frame(outtimes, t(o))
- names(out) <- c("time", mod_vars)
+ out <- cbind(outtimes, t(o))
+ colnames(out) <- c("time", mod_vars)
}
if (solution_type == "deSolve") {
@@ -208,21 +210,23 @@ mkinpredict.mkinmod <- function(x,
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()")
- }
+ }
}
+
if (map_output) {
# Output transformation for models with unobserved compartments like SFORB
- out_mapped <- data.frame(time = out[,"time"])
+ # if not already mapped in analytical solution
for (var in names(x$map)) {
- if((length(x$map[[var]]) == 1) || solution_type == "analytical") {
- out_mapped[var] <- out[, var]
+ if((length(x$map[[var]]) == 1)) {
+ out_obs[, var] <- out[, var]
} else {
- out_mapped[var] <- rowSums(out[, x$map[[var]]])
+ out_obs[, var] <- out[, x$map[[var]][1]] + out[, x$map[[var]][2]]
}
}
- return(out_mapped)
+ return(out_obs)
} else {
- return(as.data.frame(out))
+ dimnames(out) <- list(time = as.character(outtimes), c("time", mod_vars))
+ return(out)
}
}
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index b2cb4890..48df9483 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -148,6 +148,7 @@ plot.mkinfit <- function(x, fit = x,
solution_type = solution_type, atol = fit$atol, rtol = fit$rtol,
use_compiled = FALSE)
}
+ out <- as.data.frame(out)
names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars

Contact - Imprint