Fix bug in R/create_deg_func.R

Closes #13

12 files changed, 80 insertions, 41 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 0d924faf..6aa012b5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: mkin
 Type: Package
 Title: Kinetic Evaluation of Chemical Degradation Data
-Version: 1.2.8
-Date: 2023-11-26
+Version: 1.2.9
+Date: 2024-04-27
 Authors@R: c(
   person("Johannes", "Ranke", role = c("aut", "cre", "cph"),
     email = "johannes.ranke@jrwb.de",
@@ -36,4 +36,4 @@ VignetteBuilder: knitr
 BugReports: https://github.com/jranke/mkin/issues/
 URL: https://pkgdown.jrwb.de/mkin/
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
diff --git a/NAMESPACE b/NAMESPACE
index bcea2b1b..5b9a1c85 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -158,6 +158,7 @@ importFrom(lmtest,lrtest)
 importFrom(methods,is)
 importFrom(methods,signature)
 importFrom(nlme,intervals)
+importFrom(nlme,nlme)
 importFrom(parallel,detectCores)
 importFrom(parallel,mclapply)
 importFrom(parallel,parLapply)
diff --git a/NEWS.md b/NEWS.md
index f7c38915..d6139aae 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# mkin 1.2.9
+
+- 'R/create_deg_func.R: Make sure that no reversible reactions are specified in the case of two observed variables, as this is not supported
+
 # mkin 1.2.8
 
 - 'R/{mhmkin,status}.R': Deal with 'saem' fits that fail when updating an 'mhmkin' object
diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 5794c65c..d3e11f78 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -61,16 +61,21 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     ),
   ")")
 
-  if (length(obs_vars) >= 2) {
-    supported <- FALSE # except for the following cases
+  if (length(obs_vars) >= 2) supported <- FALSE
+  # Except for the following cases:
+
+  if (length(obs_vars) == 2) {
     n1 <- obs_vars[1]
     n2 <- obs_vars[2]
     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")) {
+    if (all(
+        spec[[1]]$sink == FALSE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k1 <- k_parent
       k2 <- paste0("k_", n2, if(min_ff) "_sink" else "")
@@ -80,8 +85,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # 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")) {
+    if (all(
+        use_of_ff == "max",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k1 <- k_parent
       k2 <- paste0("k_", n2)
@@ -92,8 +101,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # 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")) {
+    if (all(
+        use_of_ff == "min",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k12 <- paste0("k_", n1, "_", n2)
       k10 <- paste0("k_", n1, "_sink")
@@ -104,8 +117,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # dfop_f12_sfo
-    if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
-        spec[[1]]$type == "DFOP", spec[[2]]$type == "SFO")) {
+    if (all(
+        use_of_ff == "max",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "DFOP",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       f12 <- paste0("f_", n1, "_to_", n2)
       k2 <- paste0("k_", n2)
@@ -119,7 +136,6 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
   }
 
-
   if (supported) {
     deg_func <- function(observed, odeini, odeparms) {}
 
diff --git a/log/build.log b/log/build.log
index c2001830..de7e4488 100644
--- a/log/build.log
+++ b/log/build.log
@@ -7,5 +7,5 @@
 * checking for empty or unneeded directories
 Removed empty directory ‘mkin/inst/rmarkdown/templates/hierarchical_kinetics_parent’
 Removed empty directory ‘mkin/vignettes/web_only’
-* building ‘mkin_1.2.8.tar.gz’
+* building ‘mkin_1.2.9.tar.gz’
 
diff --git a/log/check.log b/log/check.log
index 5e777a87..58615282 100644
--- a/log/check.log
+++ b/log/check.log
@@ -1,6 +1,6 @@
 * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’
-* using R version 4.3.2 (2023-10-31)
-* using platform: x86_64-pc-linux-gnu (64-bit)
+* using R version 4.4.0 (2024-04-24)
+* using platform: x86_64-pc-linux-gnu
 * R was compiled by
     gcc (Debian 12.2.0-14) 12.2.0
     GNU Fortran (Debian 12.2.0-14) 12.2.0
@@ -9,9 +9,9 @@
 * using options ‘--no-tests --as-cran’
 * checking for file ‘mkin/DESCRIPTION’ ... OK
 * checking extension type ... Package
-* this is package ‘mkin’ version ‘1.2.8’
+* this is package ‘mkin’ version ‘1.2.9’
 * package encoding: UTF-8
-* checking CRAN incoming feasibility ... [2s/14s] Note_to_CRAN_maintainers
+* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking package namespace information ... OK
 * checking package dependencies ... OK
@@ -32,7 +32,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking for left-over files ... OK
 * checking index information ... OK
 * checking package subdirectories ... OK
-* checking R files for non-ASCII characters ... OK
+* checking code files for non-ASCII characters ... OK
 * checking R files for syntax errors ... OK
 * checking whether the package can be loaded ... OK
 * checking whether the package can be loaded with stated dependencies ... OK
@@ -62,11 +62,11 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking sizes of PDF files under ‘inst/doc’ ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... OK
+* checking examples ... [10s/10s] OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
-* checking package vignettes in ‘inst/doc’ ... OK
+* checking package vignettes ... OK
 * checking re-building of vignette outputs ... OK
 * checking PDF version of manual ... OK
 * checking HTML version of manual ... OK
diff --git a/log/test.log b/log/test.log
index 9dc63675..e92a7409 100644
--- a/log/test.log
+++ b/log/test.log
@@ -4,39 +4,40 @@
 ✔ |          5 | Analytical solutions for coupled models [1.5s]
 ✔ |          5 | Calculation of Akaike weights
 ✔ |          3 | Export dataset for reading into CAKE
-✔ |          6 | Use of precompiled symbols in mkinpredict [3.2s]
+✔ |          6 | Use of precompiled symbols in mkinpredict [3.1s]
 ✔ |         12 | Confidence intervals and p-values
-✔ |      1  12 | Dimethenamid data from 2018 [13.8s]
-✔ |         14 | Error model fitting [2.5s]
+✔ |      1  12 | Dimethenamid data from 2018 [13.0s]
+✔ |         14 | Error model fitting [2.6s]
 ✔ |          5 | Time step normalisation
 ✔ |          4 | Calculation of FOCUS chi2 error levels
 ✔ |         14 | Results for FOCUS D established in expertise for UBA (Ranke 2014)
 ✔ |          4 | Test fitting the decline of metabolites from their maximum
 ✔ |          1 | Fitting the logistic model
-✔ |         10 | Batch fitting and diagnosing hierarchical kinetic models [19.8s]
-✔ |      2  16 | Nonlinear mixed-effects models [149.6s]
+✔ |         10 | Batch fitting and diagnosing hierarchical kinetic models [19.7s]
+✔ |      2  16 | Nonlinear mixed-effects models [144.7s]
 ✔ |          3 | Test dataset classes mkinds and mkindsg
 ✔ |         10 | Special cases of mkinfit calls
 ✔ |          3 | mkinfit features
 ✔ |          8 | mkinmod model generation and printing
 ✔ |          3 | Model predictions with mkinpredict
-✔ |         12 | Multistart method for saem.mmkin models [23.6s]
-✔ |         16 | Evaluations according to 2015 NAFTA guidance [1.6s]
-✔ |          9 | Nonlinear mixed-effects models with nlme [4.0s]
-✔ |         15 | Plotting [4.7s]
+✔ |         12 | Multistart method for saem.mmkin models [23.3s]
+✔ |         16 | Evaluations according to 2015 NAFTA guidance [1.5s]
+✔ |          9 | Nonlinear mixed-effects models with nlme [3.7s]
+✔ |         15 | Plotting [4.5s]
 ✔ |          4 | Residuals extracted from mkinfit models
-✔ |      1  38 | saemix parent models [36.4s]
+✔ |      1  38 | saemix parent models [34.8s]
 ✔ |          2 | Complex test case from Schaefer et al. (2007) Piacenza paper
 ✔ |         11 | Processing of residue series
-✔ |         10 | Fitting the SFORB model [1.8s]
+✔ |         10 | Fitting the SFORB model [1.7s]
 ✔ |          1 | Summaries of old mkinfit objects
 ✔ |          5 | Summary
 ✔ |          4 | Results for synthetic data established in expertise for UBA (Ranke 2014)
-✔ |          9 | Hypothesis tests [3.1s]
+✔ |          9 | Hypothesis tests [2.9s]
 ✔ |          4 | Calculation of maximum time weighted average concentrations (TWAs)
+✔ |          2 | water-sediment
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 270.7 s
+Duration: 262.3 s
 
 ── Skipped tests (4) ───────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 5 minutes on my new system
@@ -47,4 +48,4 @@ Duration: 270.7 s
 • This still takes almost 2.5 minutes although we do not solve ODEs (1):
   'test_saemix_parent.R:143:3'
 
-[ FAIL 0 | WARN 0 | SKIP 4 | PASS 283 ]
+[ FAIL 0 | WARN 0 | SKIP 4 | PASS 285 ]
diff --git a/tests/testthat/print_dfop_saem_1.txt b/tests/testthat/print_dfop_saem_1.txt
index f7354320..3a1f1667 100644
--- a/tests/testthat/print_dfop_saem_1.txt
+++ b/tests/testthat/print_dfop_saem_1.txt
@@ -18,6 +18,6 @@ log_k1       -2.71 -2.94  -2.49
 log_k2       -4.14 -4.26  -4.01
 g_qlogis     -0.36 -0.54  -0.17
 a.1           0.93  0.69   1.17
-b.1           0.05  0.04   0.06
+b.1           0.05  0.04   0.05
 SD.log_k1     0.37  0.23   0.51
 SD.log_k2     0.23  0.14   0.31
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index 0618c715..ba7d1362 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -43,7 +43,7 @@ Optimised parameters:
                  est. lower upper
 parent_0        92.52 89.11  95.9
 log_k_parent    -1.66 -2.07  -1.3
-a.1              2.03  1.60   2.5
+a.1              2.03  1.61   2.5
 b.1              0.09  0.07   0.1
 SD.log_k_parent  0.51  0.22   0.8
 
@@ -57,7 +57,7 @@ SD.log_k_parent  0.5   0.2   0.8
 
 Variance model:
     est. lower upper
-a.1 2.03  1.60   2.5
+a.1 2.03  1.61   2.5
 b.1 0.09  0.07   0.1
 
 Backtransformed parameters:
diff --git a/tests/testthat/summary_saem_dfop_sfo_s.txt b/tests/testthat/summary_saem_dfop_sfo_s.txt
index 6468ff17..a19824ce 100644
--- a/tests/testthat/summary_saem_dfop_sfo_s.txt
+++ b/tests/testthat/summary_saem_dfop_sfo_s.txt
@@ -86,7 +86,7 @@ SD.g              0.21   0.06   0.4
 
 Variance model:
     est. lower upper
-a.1 0.93  0.79  1.06
+a.1 0.93  0.80  1.06
 b.1 0.05  0.05  0.06
 
 Resulting formation fractions:
diff --git a/tests/testthat/test_water-sediment.R b/tests/testthat/test_water-sediment.R
new file mode 100644
index 00000000..6d5693c9
--- /dev/null
+++ b/tests/testthat/test_water-sediment.R
@@ -0,0 +1,17 @@
+# Issue #13 on github
+water_sed_no_sed_sink <- mkinmod(
+  use_of_ff = "min",
+  water = mkinsub("SFO", "sediment"),
+  sediment = mkinsub("SFO", "water", sink = FALSE))
+
+ws_data <- FOCUS_D
+levels(ws_data$name) <- c("water", "sediment")
+
+test_that("An reversible reaction with the sink turned off in the second compartment works", {
+  # Solution method "analytical" was previously available, but erroneous
+  expect_error(
+    ws_fit_no_sed_sink <- mkinfit(water_sed_no_sed_sink, ws_data, quiet = TRUE, solution_type = "analytical"),
+    "Analytical solution not implemented")
+  ws_fit_no_sed_sink_default <- mkinfit(water_sed_no_sed_sink, ws_data, quiet = TRUE)
+  expect_equal(ws_fit_no_sed_sink_default$solution_type, "deSolve")
+})
diff --git a/vignettes/mkin.html b/vignettes/mkin.html
index 12b8671e..ea057d19 100644
--- a/vignettes/mkin.html
+++ b/vignettes/mkin.html
@@ -372,7 +372,7 @@ code > span.er { color: #a61717; background-color: #e3d2d2; }
 
 <h1 class="title toc-ignore">Short introduction to mkin</h1>
 <h4 class="author">Johannes Ranke</h4>
-<h4 class="date">Last change 18 May 2023 (rebuilt 2023-05-19)</h4>
+<h4 class="date">Last change 18 May 2023 (rebuilt 2024-04-27)</h4>
 
 
 <div id="TOC">