Make predefined symbols safer

We still need to create a parallel processing cluster _after_ creating a
compiled model that is moved to a user defined location, at least I did
not find another way to make it work. This is not a problem with
parallel processing without a cluster, which is not available on
Windows.

9 files changed, 113 insertions, 28 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index b97bc7e2..c851fddb 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -501,10 +501,15 @@ mkinfit <- function(mkinmod, observed,
   }
 
   # Get native symbol before iterations info for speed
+  use_symbols = FALSE
   if (solution_type == "deSolve" & use_compiled[1] != FALSE) {
-    mkinmod[["symbols"]] <- deSolve::checkDLL(dllname = mkinmod$dll_info[["name"]],
-      func = "diffs", initfunc = "initpar",
-      jacfunc = NULL, nout = 0, outnames = NULL)
+    mkinmod[["symbols"]] <- try(
+      deSolve::checkDLL(dllname = mkinmod$dll_info[["name"]],
+        func = "diffs", initfunc = "initpar",
+        jacfunc = NULL, nout = 0, outnames = NULL))
+    if (!inherits(mkinmod[["symbols"]], "try-error")) {
+      use_symbols = TRUE
+    }
   }
 
   # Get the error model and the algorithm for fitting
@@ -616,8 +621,9 @@ mkinfit <- function(mkinmod, observed,
                          odeini, outtimes,
                          solution_type = solution_type,
                          use_compiled = use_compiled,
+                         use_symbols = use_symbols,
                          method.ode = method.ode,
-                         atol = atol, rtol = rtol, 
+                         atol = atol, rtol = rtol,
                          ...)
 
       observed_index <- cbind(as.character(observed$time), as.character(observed$name))
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 29542b71..2f930adb 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -464,11 +464,11 @@ mkinmod <- function(..., use_of_ff = "max", name = NULL,
       silent = TRUE)
 
     if (!inherits(model$cf, "try-error")) {
+      if (!quiet) message("Temporary DLL for differentials generated and loaded")
       if (!is.null(dll_dir)) {
         model$dll_info <- inline::moveDLL(model$cf, name, dll_dir,
           unload = unload, overwrite = overwrite, verbose = !quiet)
       }
-      if (!quiet) message("Temporary DLL for differentials generated and loaded")
       model$dll_info <- inline::getDynLib(model$cf)
     }
   }
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 8fa41217..60456fb2 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -27,8 +27,8 @@
 #' When using compiled code, only lsoda is supported.
 #' @param use_compiled If set to \code{FALSE}, no compiled version of the
 #' [mkinmod] model is used, even if is present.
-#' @param use_symbols If set to \code{TRUE}, symbol info present in the 
-#' [mkinmod] object is used if available for accessing compiled code
+#' @param use_symbols If set to \code{TRUE} (default), symbol info present in
+#' the [mkinmod] object is used if available for accessing compiled code
 #' @param atol Absolute error tolerance, passed to the ode solver.
 #' @param rtol Absolute error tolerance, passed to the ode solver.
 #' @param maxsteps Maximum number of steps, passed to the ode solver.
@@ -174,7 +174,7 @@ mkinpredict.mkinmod <- function(x,
   if (solution_type == "deSolve") {
     if (!is.null(x$cf) & use_compiled[1] != FALSE) {
 
-      if (!is.null(x$symbols) & use_symbols[1] == TRUE) {
+      if (!is.null(x$symbols) & use_symbols) {
         lsoda_func <- x$symbols
       } else {
         lsoda_func <- "diffs"
diff --git a/R/saem.R b/R/saem.R
index 2fa770bb..83de97b0 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -583,11 +583,15 @@ saemix_model <- function(object, solution_type = "auto",
     transform_fractions <- object[[1]]$transform_fractions
 
     # Get native symbol info for speed
+    use_symbols = FALSE
     if (solution_type == "deSolve" & !is.null(mkin_model$cf)) {
-      mkin_model$symbols <- deSolve::checkDLL(
+      mkin_model$symbols <- try(deSolve::checkDLL(
         dllname = mkin_model$dll_info[["name"]],
         func = "diffs", initfunc = "initpar",
-        jacfunc = NULL, nout = 0, outnames = NULL)
+        jacfunc = NULL, nout = 0, outnames = NULL))
+      if (!inherits(mkin_model$symbols, "try-error")) {
+        use_symbols = TRUE
+      }
     }
 
     # Define the model function
diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
index 38a6bd20..cb328308 100644
--- a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
+++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
@@ -21,11 +21,18 @@ library(readxl)
 ```{r n_cores, cache = FALSE}
 n_cores <- detectCores()
 
-if (Sys.info()["sysname"] == "Windows") {
-  cl <- makePSOCKcluster(n_cores)
-} else {
-  cl <- makeForkCluster(n_cores)
+# We need to start a new cluster after defining a compiled model that is
+# saved as a DLL to the user directory, therefore we define a function
+# This is used again after defining the pathway model
+start_cluster <- function(n_cores) {
+  if (Sys.info()["sysname"] == "Windows") {
+    ret <- makePSOCKcluster(n_cores)
+  } else {
+    ret <- makeForkCluster(n_cores)
+  }
+  return(ret)
 }
+cl <- start_cluster(n_cores)
 ```
 
 \clearpage
@@ -179,6 +186,10 @@ illparms(parent_best_pH_2)
 parms(parent_best_pH_2, ci = TRUE) |> kable(digits = 3)
 ```
 
+```{r}
+stopCluster(cl)
+```
+
 \clearpage
 
 # Pathway fits
@@ -203,6 +214,7 @@ Separate evaluations of all datasets are performed with constant variance
 and using two-component error.
 
 ```{r path-1-sep, dependson = c("path-1-degmod", "ds")}
+cl <- start_cluster(n_cores)
 sforb_sep_const <- mmkin(list(sforb_path = m_sforb_sfo2), ds,
   cluster = cl, quiet = TRUE)
 sforb_sep_tc <- update(sforb_sep_const, error_model = "tc")
@@ -268,6 +280,10 @@ plot(path_1_refined)
 parms(path_1_refined, ci = TRUE) |> kable(digits = 3)
 ```
 
+```{r}
+stopCluster(cl)
+```
+
 \clearpage
 
 # Appendix
diff --git a/log/test.log b/log/test.log
index c89cf02e..fa808f22 100644
--- a/log/test.log
+++ b/log/test.log
@@ -4,20 +4,21 @@
 ✔ |         5 | Analytical solutions for coupled models [1.6s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
+✔ |         7 | Use of precompiled symbols in mkinpredict [3.3s]
 ✔ |        12 | Confidence intervals and p-values [0.4s]
-✔ |     1  12 | Dimethenamid data from 2018 [12.2s]
+✔ |     1  12 | Dimethenamid data from 2018 [13.4s]
 ────────────────────────────────────────────────────────────────────────────────
-Skip ('test_dmta.R:99'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
+Skip ('test_dmta.R:88'): 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 [2.4s]
+✔ |        14 | Error model fitting [2.5s]
 ✔ |         5 | Time step normalisation
 ✔ |         4 | Calculation of FOCUS chi2 error levels [0.3s]
-✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.4s]
-✔ |         4 | Test fitting the decline of metabolites from their maximum [0.2s]
+✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.5s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.3s]
 ✔ |         1 | Fitting the logistic model [0.1s]
-✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [19.5s]
-✔ |     1  11 | Nonlinear mixed-effects models [5.8s]
+✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [19.8s]
+✔ |     1  11 | Nonlinear mixed-effects models [6.1s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
@@ -27,12 +28,12 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | mkinfit features [0.5s]
 ✔ |         8 | mkinmod model generation and printing
 ✔ |         3 | Model predictions with mkinpredict [0.1s]
-✔ |        12 | Multistart method for saem.mmkin models [22.0s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.6s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [3.9s]
-✔ |        15 | Plotting [5.1s]
+✔ |        12 | Multistart method for saem.mmkin models [21.7s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.5s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [3.8s]
+✔ |        15 | Plotting [4.7s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [31.4s]
+✔ |     1  36 | saemix parent models [32.0s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
 Reason: This still takes almost 2.5 minutes although we do not solve ODEs
@@ -47,11 +48,11 @@ Reason: This still takes almost 2.5 minutes although we do not solve ODEs
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [0.7s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 115.0 s
+Duration: 120.2 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)
 • This still takes almost 2.5 minutes although we do not solve ODEs (1)
 
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 277 ]
diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd
index ad749fa4..792d0e47 100644
--- a/man/mkinpredict.Rd
+++ b/man/mkinpredict.Rd
@@ -15,6 +15,7 @@ mkinpredict(x, odeparms, odeini, outtimes, ...)
   outtimes = seq(0, 120, by = 0.1),
   solution_type = "deSolve",
   use_compiled = "auto",
+  use_symbols = FALSE,
   method.ode = "lsoda",
   atol = 1e-08,
   rtol = 1e-10,
@@ -67,6 +68,9 @@ parent compound.}
 \item{use_compiled}{If set to \code{FALSE}, no compiled version of the
 \link{mkinmod} model is used, even if is present.}
 
+\item{use_symbols}{If set to \code{TRUE} (default), symbol info present in
+the \link{mkinmod} object is used if available for accessing compiled code}
+
 \item{method.ode}{The solution method passed via \link{mkinpredict} to \link{ode}] in
 case the solution type is "deSolve" and we are not using compiled code.
 When using compiled code, only lsoda is supported.}
diff --git a/tests/testthat/print_mmkin_sfo_sfo_dmta.txt b/tests/testthat/print_mmkin_sfo_sfo_dmta.txt
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/tests/testthat/print_mmkin_sfo_sfo_dmta.txt
diff --git a/tests/testthat/test_compiled_symbols.R b/tests/testthat/test_compiled_symbols.R
new file mode 100644
index 00000000..e5f151eb
--- /dev/null
+++ b/tests/testthat/test_compiled_symbols.R
@@ -0,0 +1,54 @@
+context("Use of precompiled symbols in mkinpredict")
+
+test_that("We can safely use compiled code", {
+
+  # Generate temporary DLL
+  sfo_sfo_tmp <- mkinmod(DMTA = mkinsub("SFO", to = "M23"),
+    M23 = mkinsub("SFO"))
+
+  # Generate temporary DLL and move to user specified location
+  if (!dir.exists("test_dlls")) dir.create("test_dlls")
+  sfo_sfo_dll <- mkinmod(DMTA = mkinsub("SFO", to = "M23"),
+    M23 = mkinsub("SFO"),
+    dll_dir = "test_dlls",
+    name = "sfo_sfo",
+    unload = TRUE, overwrite = TRUE
+  )
+
+  if (Sys.info()["sysname"] != "Windows") {
+    # mclapply using forks
+    expect_known_output(
+      mmkin(list(sfo_sfo_dll), dmta_ds, cores = n_cores, quiet = TRUE),
+      "print_mmkin_sfo_sfo_dmta.txt"
+    )
+
+    # cluster describing itself as socket cluster
+    cl_fork <- parallel::makeForkCluster(n_cores)
+    expect_known_output(
+      mmkin(list(sfo_sfo_tmp), dmta_ds, cluster = cl_fork, quiet = TRUE),
+      "print_mmkin_sfo_sfo_dmta.txt"
+    )
+    expect_known_output(
+      mmkin(list(sfo_sfo_dll), dmta_ds, cluster = cl_fork, quiet = TRUE),
+      "print_mmkin_sfo_sfo_dmta.txt"
+    )
+    parallel::stopCluster(cl_fork)
+  }
+
+  # PSOCK cluster
+  cl_psock <- parallel::makePSOCKcluster(n_cores)
+  expect_known_output(
+    mmkin(list(sfo_sfo_tmp), dmta_ds, cluster = cl_psock, quiet = TRUE),
+    "print_mmkin_sfo_sfo_dmta.txt"
+  )
+  expect_known_output(
+    mmkin(list(sfo_sfo_dll), dmta_ds, cluster = cl_psock, quiet = TRUE),
+    "print_mmkin_sfo_sfo_dmta.txt"
+  )
+  parallel::stopCluster(cl_psock)
+
+  # Clean up
+  expect_true(file.remove("test_dlls/sfo_sfo.so"))
+  expect_true(file.remove("test_dlls"))
+})
+