Add analytical solution for DFOP-SFO

This is about twice as fast as deSolve compiled with FOCUS D

6 files changed, 65 insertions, 11 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 11559799..b29a11f7 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -16,7 +16,14 @@
 #'   benchmark(
 #'     analytical = mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
 #'     deSolve = mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
-#'     replications = 1)
+#'     replications = 2)
+#'   DFOP_SFO <- mkinmod(
+#'     parent = mkinsub("DFOP", "m1"),
+#'     m1 = mkinsub("SFO"))
+#'   benchmark(
+#'     analytical = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
+#'     deSolve = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
+#'     replications = 2)
 #' }
 create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
@@ -94,6 +101,21 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
         "(((", k2, "-", k12, "-", k10, ")*", n20, "-", k12, "*", n10, ")*exp(-", k2, "*t)+",
         k12, "*", n10, "*exp(-(", k_parent, ")*t))/(", k2, "-(", k_parent, "))")
     }
+
+    # 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")) {
+      supported <- TRUE
+      f12 <- paste0("f_", n1, "_to_", n2)
+      k2 <- paste0("k_", n2)
+      predicted_text[n2] <- paste0(
+        "((", f12, "* g - ", f12, ") * k2 * ", n10, " * exp(- k2 * t))/(k2 - ", k2, ") - ",
+        "((", f12, "* g) * k1 * ", n10, " * exp(- k1 * t))/(k1 - ", k2, ") + ",
+        "(((k1 - ", k2, ") * k2 - ", k2, "* k1 + ", k2, "^2) * ", n20, "+",
+        "((", f12, "* k1 + (", f12, "*g - ", f12, ") * ", k2, ") * k2 - ", f12, " * g * ", k2, " * k1) * ", n10, ") * ",
+        "exp( - ", k2, " * t)/((k1 - ", k2, ") * k2 - ", k2, " * k1 + ", k2, "^2)")
+    }
+
   }
 
 
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index b3cd3d30..814bf679 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -10,7 +10,7 @@ articles:
   NAFTA_examples: web_only/NAFTA_examples.html
   benchmarks: web_only/benchmarks.html
   compiled_models: web_only/compiled_models.html
-last_built: 2020-05-11T03:18Z
+last_built: 2020-05-11T13:00Z
 urls:
   reference: https://pkgdown.jrwb.de/mkin/reference
   article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/docs/reference/create_deg_func.html b/docs/reference/create_deg_func.html
index 67016be0..6d05e811 100644
--- a/docs/reference/create_deg_func.html
+++ b/docs/reference/create_deg_func.html
@@ -175,9 +175,19 @@
   <span class='fu'><a href='https://rdrr.io/pkg/rbenchmark/man/benchmark.html'>benchmark</a></span>(
     <span class='kw'>analytical</span> <span class='kw'>=</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>SFO_SFO</span>, <span class='no'>FOCUS_D</span>, <span class='kw'>solution_type</span> <span class='kw'>=</span> <span class='st'>"analytical"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>),
     <span class='kw'>deSolve</span> <span class='kw'>=</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>SFO_SFO</span>, <span class='no'>FOCUS_D</span>, <span class='kw'>solution_type</span> <span class='kw'>=</span> <span class='st'>"deSolve"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>),
-    <span class='kw'>replications</span> <span class='kw'>=</span> <span class='fl'>1</span>)</div><div class='output co'>#&gt; <span class='message'>Lade nötiges Paket: rbenchmark</span></div><div class='output co'>#&gt;         test replications elapsed relative user.self sys.self user.child
-#&gt; 1 analytical            1   0.198    1.000     0.198    0.000          0
-#&gt; 2    deSolve            1   0.350    1.768     0.348    0.001          0
+    <span class='kw'>replications</span> <span class='kw'>=</span> <span class='fl'>2</span>)</div><div class='output co'>#&gt; <span class='message'>Lade nötiges Paket: rbenchmark</span></div><div class='output co'>#&gt;         test replications elapsed relative user.self sys.self user.child
+#&gt; 1 analytical            2   0.395    1.000     0.393    0.002          0
+#&gt; 2    deSolve            2   0.693    1.754     0.692    0.000          0
+#&gt;   sys.child
+#&gt; 1         0
+#&gt; 2         0</div><div class='input'>  <span class='no'>DFOP_SFO</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span>(
+    <span class='kw'>parent</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"DFOP"</span>, <span class='st'>"m1"</span>),
+    <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>))</div><div class='output co'>#&gt; <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'>  <span class='fu'><a href='https://rdrr.io/pkg/rbenchmark/man/benchmark.html'>benchmark</a></span>(
+    <span class='kw'>analytical</span> <span class='kw'>=</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>DFOP_SFO</span>, <span class='no'>FOCUS_D</span>, <span class='kw'>solution_type</span> <span class='kw'>=</span> <span class='st'>"analytical"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>),
+    <span class='kw'>deSolve</span> <span class='kw'>=</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>DFOP_SFO</span>, <span class='no'>FOCUS_D</span>, <span class='kw'>solution_type</span> <span class='kw'>=</span> <span class='st'>"deSolve"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>),
+    <span class='kw'>replications</span> <span class='kw'>=</span> <span class='fl'>2</span>)</div><div class='output co'>#&gt;         test replications elapsed relative user.self sys.self user.child
+#&gt; 1 analytical            2   0.870    1.000     0.870        0          0
+#&gt; 2    deSolve            2   1.678    1.929     1.677        0          0
 #&gt;   sys.child
 #&gt; 1         0
 #&gt; 2         0</div><div class='input'># }
diff --git a/man/create_deg_func.Rd b/man/create_deg_func.Rd
index 178661ff..69058038 100644
--- a/man/create_deg_func.Rd
+++ b/man/create_deg_func.Rd
@@ -30,6 +30,13 @@ if (require(rbenchmark))
   benchmark(
     analytical = mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
     deSolve = mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
-    replications = 1)
+    replications = 2)
+  DFOP_SFO <- mkinmod(
+    parent = mkinsub("DFOP", "m1"),
+    m1 = mkinsub("SFO"))
+  benchmark(
+    analytical = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
+    deSolve = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
+    replications = 2)
 }
 }
diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R
index 58e328cd..def52697 100644
--- a/tests/testthat/setup_script.R
+++ b/tests/testthat/setup_script.R
@@ -40,16 +40,19 @@ fits <- mmkin(models,
 
 # One metabolite
 SFO_SFO <- mkinmod(parent = mkinsub("SFO", to = "m1"),
-                   m1 = mkinsub("SFO"),
-                   use_of_ff = "min", quiet = TRUE)
+  m1 = mkinsub("SFO"),
+  use_of_ff = "min", quiet = TRUE)
 SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", to = "m1"),
-                      m1 = mkinsub("SFO"),
-                      use_of_ff = "max", quiet = TRUE)
+  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)
+DFOP_SFO <- mkinmod(parent = mkinsub("DFOP", to = "m1"),
+  m1 = mkinsub("SFO"),
+  use_of_ff = "max", quiet = TRUE)
 
 # Avoid warning when fitting a dataset where zero value is removed
 FOCUS_D <- subset(FOCUS_2006_D, value != 0)
diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R
index 5972a18a..578258d3 100644
--- a/tests/testthat/test_analytical.R
+++ b/tests/testthat/test_analytical.R
@@ -1,6 +1,6 @@
 context("Analytical solutions for coupled models")
 
-test_that("The analytical solutions of SFO-SFO are correct", {
+test_that("The analytical solutions for SFO-SFO are correct", {
   # No sink, no formation fractions
   SFO_SFO_nosink <- mkinmod(
     parent = mkinsub("SFO", to = "m1", sink = FALSE),
@@ -44,3 +44,15 @@ test_that("The analytical solutions of SFO-SFO are correct", {
   )
 
 })
+
+test_that("The analytical solution for DFOP-SFO are correct", {
+  # With formation fraction
+  f_dfop_sfo_analytical <- mkinfit(DFOP_SFO, FOCUS_D,
+    solution_type = "analytical", quiet = TRUE)
+  f_dfop_sfo_desolve <- mkinfit(DFOP_SFO, FOCUS_D,
+    solution_type = "deSolve", quiet = TRUE)
+  expect_equal(
+    parms(f_dfop_sfo_analytical),
+    parms(f_dfop_sfo_desolve)
+  )
+})