Analytical solutions for all SFO variants

8 files changed, 128 insertions, 73 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 6c0ae40b..11559799 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -31,18 +31,27 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
   predicted_text <- character(0)
 
+  if (parent_type == "SFO") {
+    if (min_ff) {
+      targets <- c(spec[[1]]$to, if (spec[[1]]$sink) "sink" else NULL)
+      k_parent <- paste(paste0("k_", parent, "_", targets), collapse = " + ")
+    } else {
+      k_parent <- paste0("k_", parent)
+    }
+  }
+
   predicted_text[parent] <- paste0(parent_type, ".solution(t, odeini['", parent,
     if (parent_type == "SFORB") "_free", "'], ",
     switch(parent_type,
-      SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"),
-      FOMC = "alpha, beta)",
-      IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent, ")"),
-      DFOP = "k1, k2, g)",
-      SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else "", ")"),
-      HS = "k1, k2, tb)",
-      logistic = "kmax, k0, r)"
-    )
-  )
+      SFO = k_parent,
+      FOMC = "alpha, beta",
+      IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent),
+      DFOP = "k1, k2, g",
+      SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else ""),
+      HS = "k1, k2, tb",
+      logistic = "kmax, k0, r"
+    ),
+  ")")
 
   if (length(obs_vars) >= 2) {
     supported <- FALSE # except for the following cases
@@ -51,20 +60,43 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     n10 <- paste0("odeini['", parent, "']")
     n20 <- paste0("odeini['", n2, "']")
 
+    # sfo_sfo
+    if (all(spec[[1]]$sink == FALSE, length(obs_vars) == 2,
+        spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+      supported <- TRUE
+      k1 <- k_parent
+      k2 <- paste0("k_", n2, if(min_ff) "_sink" else "")
+      predicted_text[n2] <- paste0(
+        "(((", k2, "-", k1, ")*", n20, "-", k1, "*", n10, ")*exp(-", k2, "*t)+",
+        k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
+    }
+
+    # sfo_f12_sfo
     if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
         spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
       supported <- TRUE
-      k1 <- paste0("k_", n1)
+      k1 <- k_parent
       k2 <- paste0("k_", n2)
       f12 <- paste0("f_", n1, "_to_", n2)
-      if (parent_type == "SFO") {
-        predicted_text[n2] <- paste0(
-          "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
-          f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
-      }
+      predicted_text[n2] <- paste0(
+        "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
+        f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
+    }
+
+    # sfo_k120_sfo
+    if (all(use_of_ff == "min", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
+        spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+      supported <- TRUE
+      k12 <- paste0("k_", n1, "_", n2)
+      k10 <- paste0("k_", n1, "_sink")
+      k2 <- paste0("k_", n2, "_sink")
+      predicted_text[n2] <- paste0(
+        "(((", k2, "-", k12, "-", k10, ")*", n20, "-", k12, "*", n10, ")*exp(-", k2, "*t)+",
+        k12, "*", n10, "*exp(-(", k_parent, ")*t))/(", k2, "-(", k_parent, "))")
     }
   }
 
+
   if (supported) {
     deg_func <- function(observed, odeini, odeparms) {}
 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 9c3f2baa..683a6b3d 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -297,6 +297,7 @@ mkinfit <- function(mkinmod, observed,
     # This is only used for simple decline models
     if (length(obs_vars) > 1)
       stop("Decline from maximum is only implemented for models with a single observed variable")
+    observed$name <- as.character(observed$name)
 
     means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE)
     t_of_max <- means[which.max(means$value), "time"]
diff --git a/test.log b/test.log
index 74e869da..f2f9ff85 100644
--- a/test.log
+++ b/test.log
@@ -4,44 +4,33 @@ Testing mkin
 ✔ |   2       | Export dataset for reading into CAKE
 ✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1 s]
 ✔ |   4       | Calculation of FOCUS chi2 error levels [0.4 s]
-✔ |   7       | Fitting the SFORB model [3.2 s]
-✔ |   1       | Analytical solutions for coupled models [0.2 s]
+✔ |   7       | Fitting the SFORB model [3.4 s]
+✔ |   4       | Analytical solutions for coupled models [1.7 s]
 ✔ |   5       | Calculation of Akaike weights
-✔ |  10       | Confidence intervals and p-values [0.9 s]
-✔ |  14       | Error model fitting [3.7 s]
-✔ |   6       | Test fitting the decline of metabolites from their maximum [0.2 s]
+✔ |  10       | Confidence intervals and p-values [1.0 s]
+✔ |  14       | Error model fitting [4.0 s]
+✔ |   4       | Test fitting the decline of metabolites from their maximum [0.2 s]
 ✔ |   1       | Fitting the logistic model [0.2 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
 ✔ |  12       | Special cases of mkinfit calls [0.6 s]
 ✔ |   8       | mkinmod model generation and printing [0.2 s]
-✔ |   3       | Model predictions with mkinpredict [0.3 s]
-✔ |  14     2 | Evaluations according to 2015 NAFTA guidance [0.9 s]
-────────────────────────────────────────────────────────────────────────────────
-test_nafta.R:25: skip: Test data from Appendix B are correctly evaluated
-Reason: getRversion() > 4  is TRUE
-
-test_nafta.R:50: skip: Test data from Appendix D are correctly evaluated
-Reason: getRversion() > 4  is TRUE
-────────────────────────────────────────────────────────────────────────────────
-✔ |   9       | Nonlinear mixed-effects models [7.0 s]
-✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.2 s]
+✔ |   3       | Model predictions with mkinpredict [0.4 s]
+✔ |  16       | Evaluations according to 2015 NAFTA guidance [1.4 s]
+✔ |   9       | Nonlinear mixed-effects models [7.6 s]
+✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
 ✔ |   3       | Summary
-✔ |   0     1 | Plotting [0.9 s]
-────────────────────────────────────────────────────────────────────────────────
-test_plots_summary_twa.R:98: skip: Plotting mkinfit and mmkin objects is reproducible
-Reason: getRversion() > 4  is TRUE
-────────────────────────────────────────────────────────────────────────────────
+✔ |  14       | Plotting [1.4 s]
 ✔ |   4       | AIC calculation
 ✔ |   4       | Residuals extracted from mkinfit models
-✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.3 s]
+✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4 s]
 ✔ |   1       | Summaries of old mkinfit objects
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
-✔ |   9       | Hypothesis tests [7.2 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s]
+✔ |   9       | Hypothesis tests [6.5 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 32.9 s
+Duration: 36.4 s
 
-OK:       141
+OK:       158
 Failed:   0
 Warnings: 0
-Skipped:  3
+Skipped:  0
diff --git a/tests/testthat/DFOP_FOCUS_C_messages.txt b/tests/testthat/DFOP_FOCUS_C_messages.txt
index 7284dd19..6aa73e01 100644
--- a/tests/testthat/DFOP_FOCUS_C_messages.txt
+++ b/tests/testthat/DFOP_FOCUS_C_messages.txt
@@ -81,7 +81,7 @@ Sum of squared residuals at call 55: 4.364078
  85.016328  -0.776316  -4.027611   1.248897 
 Sum of squared residuals at call 56: 4.364078
  85.016327  -0.776316  -4.027611   1.248897 
-Sum of squared residuals at call 57: 4.364078
+Sum of squared residuals at call 57: 4.364077
  85.016327  -0.776316  -4.027611   1.248897 
  85.016327  -0.776316  -4.027611   1.248897 
  85.008939  -0.777792  -4.026307   1.247720 
diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R
index 86d5089f..58e328cd 100644
--- a/tests/testthat/setup_script.R
+++ b/tests/testthat/setup_script.R
@@ -39,25 +39,27 @@ fits <- mmkin(models,
   quiet = TRUE, cores = n_cores)
 
 # One metabolite
-SFO_SFO <- mkinmod(parent = list(type = "SFO", to = "m1"),
-                   m1 = list(type = "SFO"),
+SFO_SFO <- mkinmod(parent = mkinsub("SFO", to = "m1"),
+                   m1 = mkinsub("SFO"),
                    use_of_ff = "min", quiet = TRUE)
-SFO_SFO.ff <- mkinmod(parent = list(type = "SFO", to = "m1"),
-                      m1 = list(type = "SFO"),
+SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", to = "m1"),
+                      m1 = mkinsub("SFO"),
                       use_of_ff = "max", quiet = TRUE)
 SFO_SFO.ff.nosink <- mkinmod(
   parent = mkinsub("SFO", "m1", sink = FALSE),
   m1 = mkinsub("SFO"), quiet = TRUE, use_of_ff = "max")
+FOMC_SFO <- mkinmod(parent = mkinsub("FOMC", to = "m1"),
+  m1 = mkinsub("SFO"), quiet = TRUE)
 
-f_sfo_sfo_desolve <- mkinfit(SFO_SFO,
-  subset(FOCUS_2006_D, value != 0),
+# Avoid warning when fitting a dataset where zero value is removed
+FOCUS_D <- subset(FOCUS_2006_D, value != 0)
+
+f_sfo_sfo_desolve <- mkinfit(SFO_SFO, FOCUS_D,
   solution_type = "deSolve", quiet = TRUE)
-f_sfo_sfo_eigen <- mkinfit(SFO_SFO,
-  subset(FOCUS_2006_D, value != 0),
+f_sfo_sfo_eigen <- mkinfit(SFO_SFO, FOCUS_D,
   solution_type = "eigen", quiet = TRUE)
 
-f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff,
-  subset(FOCUS_2006_D, value != 0),
+f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff, FOCUS_D,
   quiet = TRUE)
 
 SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R
index 3d30e042..5972a18a 100644
--- a/tests/testthat/test_analytical.R
+++ b/tests/testthat/test_analytical.R
@@ -1,11 +1,46 @@
 context("Analytical solutions for coupled models")
 
-test_that("The analytical solution of SFO-SFO is correct", {
-  f_sfo_sfo.ff.analytical <- mkinfit(SFO_SFO.ff,
-    subset(FOCUS_2006_D, value != 0),
-    quiet = TRUE)
+test_that("The analytical solutions of SFO-SFO are correct", {
+  # No sink, no formation fractions
+  SFO_SFO_nosink <- mkinmod(
+    parent = mkinsub("SFO", to = "m1", sink = FALSE),
+    m1 = mkinsub("SFO"),
+    use_of_ff = "min", quiet = TRUE)
+  f_sfo_sfo_nosink <- mkinfit(SFO_SFO_nosink, FOCUS_D, quiet = TRUE)
+  f_sfo_sfo_nosink_deSolve <- mkinfit(SFO_SFO_nosink, FOCUS_D,
+    solution_type = "deSolve", quiet = TRUE)
+  expect_equal(
+    parms(f_sfo_sfo_nosink),
+    parms(f_sfo_sfo_nosink_deSolve)
+  )
+
+  # No sink, with formation fractions
+  SFO_SFO.ff_nosink <- mkinmod(
+    parent = mkinsub("SFO", to = "m1", sink = FALSE),
+    m1 = mkinsub("SFO"),
+    use_of_ff = "max", quiet = TRUE)
+  f_sfo_sfo_nosink <- mkinfit(SFO_SFO.ff_nosink, FOCUS_D, quiet = TRUE)
+  f_sfo_sfo_nosink_deSolve <- mkinfit(SFO_SFO.ff_nosink, FOCUS_D,
+    solution_type = "deSolve", quiet = TRUE)
+  expect_equal(
+    parms(f_sfo_sfo_nosink),
+    parms(f_sfo_sfo_nosink_deSolve)
+  )
+
+  # Without formation fraction
+  f_sfo_sfo_analytical <- mkinfit(SFO_SFO, FOCUS_D,
+    solution_type = "analytical", quiet = TRUE)
+  expect_equal(
+    parms(f_sfo_sfo_analytical),
+    parms(f_sfo_sfo_desolve)
+  )
+
+  # With formation fraction
+  f_sfo_sfo.ff_desolve <- mkinfit(SFO_SFO.ff, FOCUS_D,
+    solution_type = "deSolve", quiet = TRUE)
   expect_equal(
     parms(f_sfo_sfo.ff),
-    parms(f_sfo_sfo.ff.analytical)
+    parms(f_sfo_sfo.ff_desolve)
   )
+
 })
diff --git a/tests/testthat/test_from_max_mean.R b/tests/testthat/test_from_max_mean.R
index 002b39e2..5a66c1ef 100644
--- a/tests/testthat/test_from_max_mean.R
+++ b/tests/testthat/test_from_max_mean.R
@@ -6,15 +6,15 @@ test_that("Fitting from maximum mean value works", {
     expect_error(mkinfit(SFO_SFO, FOCUS_2006_D, from_max_mean = TRUE),
       "only implemented for models with a single observed variable"),
     "Observations with value of zero were removed")
-  
+
   # We can either explicitly create a model for m1, or subset the data
   SFO_m1 <- mkinmod(m1 = mkinsub("SFO"))
-  f.1 <- expect_warning(mkinfit(SFO_m1, FOCUS_2006_D, from_max_mean = TRUE, quiet = TRUE))
-  expect_equivalent(endpoints(f.1)$distimes["m1", ], c(170.8, 567.5), 
-                    scale = 1, tolerance = 0.1)
+  f.1 <- mkinfit(SFO_m1, FOCUS_D, from_max_mean = TRUE, quiet = TRUE)
+  expect_equivalent(endpoints(f.1)$distimes["m1", ], c(170.8, 567.5),
+    scale = 1, tolerance = 0.1)
 
-  f.2 <- expect_warning(mkinfit("SFO", subset(FOCUS_2006_D, name == "m1"), 
-                                from_max_mean = TRUE, quiet = TRUE))
-  expect_equivalent(endpoints(f.2)$distimes["m1", ], c(170.8, 567.5), 
+  f.2 <- mkinfit("SFO", subset(FOCUS_D, name == "m1"),
+    from_max_mean = TRUE, quiet = TRUE)
+  expect_equivalent(endpoints(f.2)$distimes["m1", ], c(170.8, 567.5),
                     scale = 1, tolerance = 0.1)
 })
diff --git a/tests/testthat/test_mkinfit_errors.R b/tests/testthat/test_mkinfit_errors.R
index c1e9da1d..85ee574e 100644
--- a/tests/testthat/test_mkinfit_errors.R
+++ b/tests/testthat/test_mkinfit_errors.R
@@ -4,28 +4,24 @@ test_that("mkinfit stops to prevent and/or explain user errors", {
   expect_error(mkinfit("foo", FOCUS_2006_A))
   expect_error(mkinfit(3, FOCUS_2006_A))
 
-  # We remove zero observations from FOCUS_2006_D beforehand in
-  # order to avoid another expect_warning in the code
-  FOCUS_2006_D <- subset(FOCUS_2006_D, value != 0)
-
   # We get a warning if we use transform_fractions = FALSE with formation fractions
   # and an error if any pathway to sink is turned off as well
   expect_warning(
     expect_error(
-      mkinfit(SFO_SFO.ff.nosink, FOCUS_2006_D, transform_fractions = FALSE, quiet = TRUE),
+      mkinfit(SFO_SFO.ff.nosink, FOCUS_D, transform_fractions = FALSE, quiet = TRUE),
       "turn off pathways to sink"
       ),
     "sum of formation fractions may exceed one")
 
-  expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D, transform_fractions = TRUE,
+  expect_error(mkinfit(SFO_SFO.ff, FOCUS_D, transform_fractions = TRUE,
                        parms.ini = c(f_parent_to_m1 = 0.5), fixed_parms = "f_parent_to_m1", quiet = TRUE),
    "not supported")
 
-  expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D,
+  expect_error(mkinfit(SFO_SFO.ff, FOCUS_D,
                        parms.ini = c(f_parent_to_m1 = 1.1), quiet = TRUE),
    "sum up to more than 1")
 
-  expect_error(mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "analytical"), "not implemented")
+  expect_error(mkinfit(FOMC_SFO, FOCUS_D, solution_type = "analytical"), "not implemented")
 
   expect_error(mkinfit("FOMC", FOCUS_2006_A, solution_type = "eigen"), "coefficient matrix not present")
 })