Test saemix without special analytical solutions

Also increase the performance a bit (from about 210 s to about 140 s in
the case of DFOP with four chains and 300, 100 iterations).

4 files changed, 50 insertions, 28 deletions

diff --git a/NEWS.md b/NEWS.md
index 1524e1cf..4a3b352d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,8 @@
 
 - 'CAKE_export': Check for validity of the map argument, updates
 
+- 'saem()': Slightly improve speed in the case that no analytical solution for saemix is implemented, activate a test of the respective code
+
 # mkin 1.1.0 (2022-03-14)
 
 ## Mixed-effects models
diff --git a/R/saem.R b/R/saem.R
index d3b23861..36997ad7 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -239,10 +239,6 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
   degparms_start = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
 {
-  if (packageVersion("saemix") < "3.0") {
-    stop("To use the interface to saemix, you need to install a version >= 3.0\n")
-  }
-
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
@@ -444,9 +440,22 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
     solution_type = "analytical saemix"
   } else {
 
+    if (transformations == "saemix") {
+      stop("Using saemix transformations is only supported if an analytical solution is implemented for saemix")
+    }
+
     if (solution_type == "auto")
       solution_type <- object[[1]]$solution_type
 
+    # Define some variables to avoid function calls in model function
+    transparms_optim_names <- names(degparms_optim)
+    odeini_optim_names <- gsub('_0$', '', odeini_optim_parm_names)
+    diff_names <- names(mkin_model$diffs)
+    ode_transparms_optim_names <- setdiff(transparms_optim_names, odeini_optim_parm_names)
+    transform_rates <- object[[1]]$transform_rates
+    transform_fractions <- object[[1]]$transform_fractions
+
+    # Define the model function
     model_function <- function(psi, id, xidep) {
 
       uid <- unique(id)
@@ -454,22 +463,21 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
       res_list <- lapply(uid, function(i) {
 
         transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
-        names(transparms_optim) <- names(degparms_optim)
+        names(transparms_optim) <- transparms_optim_names
 
         odeini_optim <- transparms_optim[odeini_optim_parm_names]
-        names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names)
+        names(odeini_optim) <- odeini_optim_names
 
-        odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
+        odeini <- c(odeini_optim, odeini_fixed)[diff_names]
 
-        ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names)
         odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
-          transform_rates = object[[1]]$transform_rates,
-          transform_fractions = object[[1]]$transform_fractions)
+          transform_rates = transform_rates,
+          transform_fractions = transform_fractions)
         odeparms <- c(odeparms_optim, odeparms_fixed)
 
-        xidep_i <- subset(xidep, id == i)
+        xidep_i <- xidep[which(id == i), ]
 
-        if (solution_type == "analytical") {
+        if (solution_type[1] == "analytical") {
           out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
         } else {
 
diff --git a/test.log b/test.log
index 4c273c19..2fc6a43c 100644
--- a/test.log
+++ b/test.log
@@ -3,18 +3,18 @@ Loading required package: parallel
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.4s]
+✔ |         5 | Analytical solutions for coupled models [3.3s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.0s]
-✔ |     1  12 | Dimethenamid data from 2018 [32.7s]
+✔ |     1  12 | Dimethenamid data from 2018 [31.0s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ────────────────────────────────────────────────────────────────────────────────
-✔ |        14 | Error model fitting [5.0s]
+✔ |        14 | Error model fitting [4.6s]
 ✔ |         5 | Time step normalisation
-✔ |         4 | Calculation of FOCUS chi2 error levels [0.6s]
+✔ |         4 | Calculation of FOCUS chi2 error levels [0.5s]
 ✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s]
 ✔ |         4 | Test fitting the decline of metabolites from their maximum [0.3s]
 ✔ |         1 | Fitting the logistic model [0.2s]
@@ -24,28 +24,28 @@ Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |         3 | Test dataset classes mkinds and mkindsg
-✔ |        10 | Special cases of mkinfit calls [0.5s]
-✔ |         3 | mkinfit features [0.7s]
-✔ |         8 | mkinmod model generation and printing [0.3s]
-✔ |         3 | Model predictions with mkinpredict [0.4s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.6s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.6s]
-✔ |        16 | Plotting [10.4s]
+✔ |        10 | Special cases of mkinfit calls [0.4s]
+✔ |         3 | mkinfit features [0.6s]
+✔ |         8 | mkinmod model generation and printing [0.2s]
+✔ |         3 | Model predictions with mkinpredict [0.3s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.4s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [8.0s]
+✔ |        16 | Plotting [9.9s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        23 | saemix parent models [28.0s]
+✔ |        25 | saemix parent models [170.8s]
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
-✔ |         7 | Fitting the SFORB model [3.7s]
+✔ |         7 | Fitting the SFORB model [3.6s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         4 | Summary [0.1s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ |         9 | Hypothesis tests [8.1s]
+✔ |         9 | Hypothesis tests [8.0s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 112.8 s
+Duration: 251.5 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
 • Fitting with saemix takes around 10 minutes when using deSolve (1)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 222 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 224 ]
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 340db6f1..731228d9 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -92,3 +92,15 @@ test_that("Parent fits using saemix are correctly implemented", {
   # HS would likely benefit from implemenation of transformations = "saemix"
 })
 
+test_that("We can also use mkin solution methods for saem", {
+  expect_error(saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix", solution_type = "analytical"),
+    "saemix transformations is only supported if an analytical solution is implemented"
+  )
+  skip_on_cran() # This still takes almost 2.5 minutes although we do not solve ODEs
+  dfop_saemix_3 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin",
+    solution_type = "analytical")
+  distimes_dfop <- endpoints(dfop_saemix_1)$distimes
+  distimes_dfop_analytical <- endpoints(dfop_saemix_3)$distimes
+  rel_diff <- abs(distimes_dfop_analytical - distimes_dfop) / distimes_dfop
+  expect_true(all(rel_diff < 0.01))
+})