Support SFORB with formation fractions

9 files changed, 59 insertions, 51 deletions

diff --git a/NEWS.md b/NEWS.md
index 790e2e03..ba1daa87 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,8 @@
-# mkin 0.9.49.11 (unreleased)
+# mkin 0.9.50.1 (unreleased)
+
+- Support SFORB with formation fractions
+
+# mkin 0.9.49.11 (2020-04-20)
 
 - Increase a test tolerance to make it pass on all CRAN check machines
 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 1c409569..5c092612 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -378,8 +378,9 @@ mkinfit <- function(mkinmod, observed,
     if (parmname == "r") parms.ini[parmname] = 0.2
   }
   # Default values for formation fractions in case they are present
-  for (box in mod_vars) {
-    f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
+  for (obs_var in obs_vars) {
+    origin <- mkinmod$map[[obs_var]][[1]]
+    f_names <- mkinmod$parms[grep(paste0("^f_", origin), mkinmod$parms)]
     if (length(f_names) > 0) {
       # We need to differentiate between default and specified fractions
       # and set the unspecified to 1 - sum(specified)/n_unspecified
@@ -388,9 +389,9 @@ mkinfit <- function(mkinmod, observed,
       sum_f_specified = sum(parms.ini[f_specified_names])
       if (sum_f_specified > 1) {
         stop("Starting values for the formation fractions originating from ",
-             box, " sum up to more than 1.")
+             origin, " sum up to more than 1.")
       }
-      if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
+      if (mkinmod$spec[[obs_var]]$sink) n_unspecified = length(f_default_names) + 1
       else {
         n_unspecified = length(f_default_names)
       }
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 62f16e73..4587e210 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -182,7 +182,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
                 "are only supported with formation fractions (use_of_ff = 'max')")
         }
         if(spec[[varname]]$sink) {
-          # If sink is required, add first-order/IORE sink term
+          # If sink is requested, add first-order/IORE sink term
           k_compound_sink <- paste("k", box_1, "sink", sep = "_")
           if(spec[[varname]]$type == "IORE") {
             k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")
@@ -197,7 +197,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
         } else { # otherwise no decline term needed here
           decline_term = "0"
         }
-      } else {
+      } else { # Maximum use of formation fractions
         k_compound <- paste("k", box_1, sep = "_")
         if(spec[[varname]]$type == "IORE") {
           k_compound <- paste("k__iore", box_1, sep = "_")
@@ -236,26 +236,13 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
     diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
     if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms
       box_2 = map[[varname]][[2]]
-      if (use_of_ff == "min") { # Minimum use of formation fractions
-        k_free_bound <- paste("k", varname, "free", "bound", sep = "_")
-        k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
-        parms <- c(parms, k_free_bound, k_bound_free)
-        reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
-          k_bound_free, "*", box_2)
-        reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
-          k_bound_free, "*", box_2)
-      } else { # Use formation fractions also for the free compartment
-        stop("The maximum use of formation fractions is not supported for SFORB models")
-        # The problems were: Calculation of dissipation times did not work in this case
-        # and the coefficient matrix is not generated correctly by the code present
-        # in this file in this case
-        #f_free_bound <- paste("f", varname, "free", "bound", sep = "_")
-        #k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
-        #parms <- c(parms, f_free_bound, k_bound_free)
-        #reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2)
-        #reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-",
-        #  k_bound_free, "*", box_2)
-      }
+      k_free_bound <- paste("k", varname, "free", "bound", sep = "_")
+      k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
+      parms <- c(parms, k_free_bound, k_bound_free)
+      reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
+        k_bound_free, "*", box_2)
+      reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
+        k_bound_free, "*", box_2)
       diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
       diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
     } #}}}
@@ -286,7 +273,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
             k_from_to, "*", origin_box)
         } else {
           # Do not introduce a formation fraction if this is the only target
-          if (spec[[origin_box]]$sink == FALSE && n_targets == 1) {
+          if (spec[[varname]]$sink == FALSE && n_targets == 1) {
             diffs[[target_box]] <- paste(diffs[[target_box]], "+",
                                          decline_term)
           } else {
@@ -302,7 +289,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
 
   model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff)
 
-  # Create coefficient matrix if appropriate#{{{
+  # Create coefficient matrix if possible #{{{
   if (mat) {
     boxes <- names(diffs)
     n <- length(boxes)
@@ -321,12 +308,12 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
 
           } else {          # off-diagonal elements
             k.candidate = paste("k", from, to, sep = "_")
-      if (sub("_free$", "", from) == sub("_bound$", "", to)) {
+            if (sub("_free$", "", from) == sub("_bound$", "", to)) {
               k.candidate = paste("k", sub("_free$", "_free_bound", from), sep = "_")
-      }
-      if (sub("_bound$", "", from) == sub("_free$", "", to)) {
+            }
+            if (sub("_bound$", "", from) == sub("_free$", "", to)) {
               k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep = "_")
-      }
+            }
             k.effective = intersect(model$parms, k.candidate)
             m[to, from] = ifelse(length(k.effective) > 0,
                 k.effective, "0")
@@ -350,7 +337,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
           } else {          # off-diagonal elements
             f.candidate = paste("f", from, "to", to, sep = "_")
             k.candidate = paste("k", from, to, sep = "_")
-            # SFORB with maximum use of formation fractions not implemented, see above
+            k.candidate = sub("free.*bound", "free_bound", k.candidate)
+            k.candidate = sub("bound.*free", "bound_free", k.candidate)
             m[to, from] = ifelse(f.candidate %in% model$parms,
               paste(f.candidate, " * k_", from, sep = ""),
               ifelse(k.candidate %in% model$parms, k.candidate, "0"))
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index d2b851b6..020917b9 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -60,8 +60,8 @@ A list of class \code{mkinmod} for use with \code{\link{mkinfit}},
     The coefficient matrix, if the system of differential equations can be
     represented by one.
   }
-  \item{ll}{
-    The likelihood function, taking the parameter vector as the first argument.
+  \item{cf}{
+    If generated, the compiled function as returned by cfunction.
   }
 }
 \description{
diff --git a/test.log b/test.log
index 684f6729..a7879c33 100644
--- a/test.log
+++ b/test.log
@@ -3,16 +3,16 @@ Testing mkin
 ✔ |  OK F W S | Context
 ✔ |   2       | Export dataset for reading into CAKE
 ✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.5 s]
-✔ |   4       | Calculation of FOCUS chi2 error levels [2.2 s]
-✔ |   4       | Fitting the SFORB model [1.7 s]
+✔ |   4       | Calculation of FOCUS chi2 error levels [2.1 s]
+✔ |   6       | Fitting the SFORB model [8.8 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |  10       | Confidence intervals and p-values [9.4 s]
-✔ |  14       | Error model fitting [38.0 s]
+✔ |  14       | Error model fitting [37.2 s]
 ✔ |   6       | Test fitting the decline of metabolites from their maximum [0.8 s]
 ✔ |   1       | Fitting the logistic model [0.9 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
 ✔ |  12       | Special cases of mkinfit calls [2.3 s]
-✔ |   9       | mkinmod model generation and printing [0.2 s]
+✔ |   8       | mkinmod model generation and printing [0.2 s]
 ✔ |   3       | Model predictions with mkinpredict [0.3 s]
 ✔ |  16       | Evaluations according to 2015 NAFTA guidance [4.1 s]
 ✔ |   9       | Nonlinear mixed-effects models [11.9 s]
@@ -21,15 +21,15 @@ Testing mkin
 ✔ |  14       | Plotting [4.9 s]
 ✔ |   4       | AIC calculation
 ✔ |   4       | Residuals extracted from mkinfit models
-✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [5.3 s]
+✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [5.4 s]
 ✔ |   1       | Summaries of old mkinfit objects
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.1 s]
-✔ |   9       | Hypothesis tests [36.8 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.9 s]
+✔ |   9       | Hypothesis tests [36.7 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 131.8 s
+Duration: 137.8 s
 
-OK:       154
+OK:       155
 Failed:   0
 Warnings: 0
 Skipped:  0
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index 81a4ed74..8faa19db 100644
--- a/tests/testthat/FOCUS_2006_D.csf
+++ b/tests/testthat/FOCUS_2006_D.csf
@@ -5,7 +5,7 @@ Description:
 MeasurementUnits: % AR
 TimeUnits: days
 Comments: Created using mkin::CAKE_export
-Date: 2020-04-20
+Date: 2020-04-22
 Optimiser: IRLS
 
 [Data]
diff --git a/tests/testthat/test_SFORB.R b/tests/testthat/test_SFORB.R
index b168a4ee..ad9881a8 100644
--- a/tests/testthat/test_SFORB.R
+++ b/tests/testthat/test_SFORB.R
@@ -9,4 +9,22 @@ test_that("Fitting the SFORB model is equivalent to fitting DFOP", {
   expect_match(s_sforb, "Estimated Eigenvalues of SFORB model\\(s\\):")
   expect_match(s_sforb, "parent_b1 parent_b2")
   expect_match(s_sforb, "0.45956 *0.01785")
+
+  DFOP_SFO <- mkinmod(parent = mkinsub("DFOP", "M1"),
+    M1 = mkinsub("SFO"),
+    use_of_ff = "max", quiet = TRUE)
+  SFORB_SFO <- mkinmod(parent = mkinsub("SFORB", "M1"),
+    M1 = mkinsub("SFO"),
+    use_of_ff = "max", quiet = TRUE)
+
+  SFORB_SFO$coefmat
+
+  f_dfop_sfo <- mkinfit(DFOP_SFO, DFOP_par_c, quiet = TRUE)
+  f_sforb_sfo <- mkinfit(SFORB_SFO, DFOP_par_c, quiet = TRUE)
+  f_sforb_sfo_eigen <- mkinfit(SFORB_SFO, DFOP_par_c, solution_type = "eigen", quiet = TRUE)
+
+  expect_equivalent(endpoints(f_sforb_sfo)$distimes, endpoints(f_dfop_sfo)$distimes,
+    tolerance = 1e-6)
+  expect_equivalent(endpoints(f_sforb_sfo_eigen)$distimes, endpoints(f_dfop_sfo)$distimes,
+    tolerance = 1e-6)
 })
diff --git a/tests/testthat/test_mkinfit_errors.R b/tests/testthat/test_mkinfit_errors.R
index 7987d291..940e0b2e 100644
--- a/tests/testthat/test_mkinfit_errors.R
+++ b/tests/testthat/test_mkinfit_errors.R
@@ -7,6 +7,7 @@ test_that("mkinfit stops to prevent and/or explain user errors", {
   # 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(
@@ -14,7 +15,7 @@ test_that("mkinfit stops to prevent and/or explain user errors", {
       mkinfit(SFO_SFO.ff.nosink, FOCUS_2006_D, transform_fractions = FALSE, quiet = TRUE),
       "turn off pathways to sink"
       ),
-    "sum of formation fractions")
+    "sum of formation fractions may exceed one")
 
   expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D, transform_fractions = TRUE,
                        parms.ini = c(f_parent_to_m1 = 0.5), fixed_parms = "f_parent_to_m1", quiet = TRUE),
diff --git a/tests/testthat/test_mkinmod.R b/tests/testthat/test_mkinmod.R
index 1bdf094b..3178186a 100644
--- a/tests/testthat/test_mkinmod.R
+++ b/tests/testthat/test_mkinmod.R
@@ -20,10 +20,6 @@ test_that("mkinmod stops to prevent and/or explain user errors", {
   expect_error(mkinmod(parent = mkinsub("IORE", "m1"),
                        m1 = mkinsub("SFO"), use_of_ff = "min"),
                "only supported with formation fractions")
-
-  expect_error(mkinmod(parent = mkinsub("SFORB", "m1"),
-                       m1 = mkinsub("SFO"), use_of_ff = "max"),
-               "not supported")
 })
 
 test_that("Printing mkinmod models is reproducible", {