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.

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