We get about 25% performance gain

with the custom lsoda call, avoiding repeated getNativeSymbolInfo calls.
It's just that we should not be calling foreign functions from different
packages, because the may change without notice. Using
getNativeSymbolInfo for "call_lsoda" avoids the CRAN note, and a similar
call could probably be used for "unlock_solver", avoiding the NOTE in
checks for cran, but we should not do this in a CRAN package.

5 files changed, 131 insertions, 33 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 693778fd..0d9246dd 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -500,6 +500,15 @@ mkinfit <- function(mkinmod, observed,
     }
   }
 
+  # Get native symbol before iterations info for speed
+  call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve")
+  if (solution_type == "deSolve" & use_compiled[1] != FALSE) {
+    mkinmod$diffs_address <- getNativeSymbolInfo("diffs",
+      PACKAGE = mkinmod$dll_info[["name"]])$address
+    mkinmod$initpar_address <- getNativeSymbolInfo("initpar",
+      PACKAGE = mkinmod$dll_info[["name"]])$address
+  }
+
   # Get the error model and the algorithm for fitting
   err_mod <- match.arg(error_model)
   error_model_algorithm = match.arg(error_model_algorithm)
@@ -610,7 +619,8 @@ mkinfit <- function(mkinmod, observed,
                          solution_type = solution_type,
                          use_compiled = use_compiled,
                          method.ode = method.ode,
-                         atol = atol, rtol = rtol, ...)
+                         atol = atol, rtol = rtol,
+                         call_lsoda = call_lsoda, ...)
 
       observed_index <- cbind(as.character(observed$time), as.character(observed$name))
       observed$predicted <- out[observed_index]
@@ -892,7 +902,10 @@ mkinfit <- function(mkinmod, observed,
   fit$calls <- calls
   fit$time <- fit_time
 
-  # We also need the model and a model name for summary and plotting
+  # We also need the model and a model name for summary and plotting,
+  # but without address info that will become invalid
+  mkinmod$diffs_address <- NULL
+  mkinmod$initpar_address <- NULL
   fit$mkinmod <- mkinmod
   fit$mkinmod$name <- mkinmod_name
   fit$obs_vars <- obs_vars
diff --git a/R/mkinmod.R b/R/mkinmod.R
index d8740aed..b1fb57cb 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -45,7 +45,9 @@
 #' @param dll_dir Directory where an DLL object, if generated internally by
 #'   [inline::cfunction()], should be saved.  The DLL will only be stored in a
 #'   permanent location for use in future sessions, if 'dll_dir' and 'name'
-#'   are specified.
+#'   are specified. This is helpful if fit objects are cached e.g. by knitr,
+#'   as the cache remains functional across sessions if the DLL is stored in
+#'   a user defined location.
 #' @param unload If a DLL from the target location in 'dll_dir' is already
 #'   loaded, should that be unloaded first?
 #' @param overwrite If a file exists at the target DLL location in 'dll_dir',
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 0dc9cf55..11a3b35b 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -19,26 +19,24 @@
 #' @param solution_type The method that should be used for producing the
 #' predictions. This should generally be "analytical" if there is only one
 #' observed variable, and usually "deSolve" in the case of several observed
-#' variables. The third possibility "eigen" is faster but not applicable to
-#' some models e.g.  using FOMC for the parent compound.
+#' variables. The third possibility "eigen" is fast in comparison to uncompiled
+#' ODE models, but not applicable to some models, e.g. using FOMC for the
+#' parent compound.
 #' @param method.ode The solution method passed via [mkinpredict] to [ode]] in
-#' case the solution type is "deSolve". The default "lsoda" is performant, but
-#' sometimes fails to converge.
+#' case the solution type is "deSolve" and we are not using compiled code.
 #' @param use_compiled If set to \code{FALSE}, no compiled version of the
 #' [mkinmod] model is used, even if is present.
-#' @param atol Absolute error tolerance, passed to [ode]. Default is 1e-8,
-#' lower than in [lsoda].
-#' @param rtol Absolute error tolerance, passed to [ode]. Default is 1e-10,
-#' much lower than in [lsoda].
-#' @param maxsteps Maximum number of steps, passed to [ode].
+#' @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.
 #' @param map_output Boolean to specify if the output should list values for
 #'   the observed variables (default) or for all state variables (if set to
 #'   FALSE). Setting this to FALSE has no effect for analytical solutions,
 #'   as these always return mapped output.
 #' @param na_stop Should it be an error if [ode] returns NaN values
+#' @param call_lsoda The address of the compiled function "call_lsoda"
 #' @param \dots Further arguments passed to the ode solver in case such a
 #'   solver is used.
-#' @import deSolve
 #' @return A matrix with the numeric solution in wide format
 #' @author Johannes Ranke
 #' @examples
@@ -116,9 +114,10 @@ mkinpredict.mkinmod <- function(x,
   outtimes = seq(0, 120, by = 0.1),
   solution_type = "deSolve",
   use_compiled = "auto",
-  method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000,
+  method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000L,
   map_output = TRUE,
   na_stop = TRUE,
+  call_lsoda = NULL,
   ...)
 {
 
@@ -173,20 +172,80 @@ mkinpredict.mkinmod <- function(x,
   if (solution_type == "deSolve") {
     if (!is.null(x$cf) & use_compiled[1] != FALSE) {
 
-      out <- deSolve::ode(
-        y = odeini,
-        times = outtimes,
-        func = "diffs",
-        initfunc = "initpar",
-        dllname = if (is.null(x$dll_info)) inline::getDynLib(x$cf)[["name"]]
-          else x$dll_info[["name"]],
-        parms = odeparms[x$parms], # Order matters when using compiled models
-        method = method.ode,
-        atol = atol,
-        rtol = rtol,
-        maxsteps = maxsteps,
-        ...
+#       out <- deSolve::ode(
+#         y = odeini,
+#         times = outtimes,
+#         func = "diffs",
+#         initfunc = "initpar",
+#         dllname = x$dll_info[["name"]],
+#         parms = odeparms[x$parms], # Order matters when using compiled models
+#         method = method.ode,
+#         atol = atol,
+#         rtol = rtol,
+#         maxsteps = maxsteps,
+#         ...
+#       )
+#
+      # Prepare call to "call_lsoda"
+      # Simplified code from deSolve::lsoda() adapted to our use case
+      if (is.null(call_lsoda)) {
+        call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve")
+      }
+      if (is.null(x$diffs_address)) {
+        x$diffs_address <- getNativeSymbolInfo("diffs",
+          PACKAGE = x$dll_info[["name"]])$address
+        x$initpar_address <- getNativeSymbolInfo("initpar",
+          PACKAGE = x$dll_info[["name"]])$address
+      }
+      rwork <- vector("double", 20)
+      rwork[] <- 0.
+      rwork[6] <- max(abs(diff(outtimes)))
+
+      iwork <- vector("integer", 20)
+      iwork[] <- 0
+      iwork[6] <- maxsteps
+
+      n <- length(odeini)
+      lmat <- n^2 + 2     # from deSolve::lsoda(), for Jacobian type full, internal (2)
+      # hard-coded default values of lsoda():
+      maxordn <- 12L
+      maxords <- 5L
+      lrn <- 20 + n * (maxordn + 1) + 3 * n  # length in case non-stiff method
+      lrs <- 20 + n * (maxords + 1) + 3 * n + lmat        # length in case stiff method
+      lrw <- max(lrn, lrs)                       # actual length: max of both
+      liw <- 20 + n
+
+      on.exit(.C("unlock_solver", PACKAGE = "deSolve"))
+
+      out_raw <- .Call(call_lsoda,
+        as.double(odeini),   # y
+        as.double(outtimes), # times
+        x$diffs_address,     # derivfunc
+        as.double(odeparms[x$parms]), # parms
+        rtol, atol,
+        NULL, NULL, # rho, tcrit
+        NULL, # jacfunc
+        x$initpar_address, # initfunc
+        NULL, # eventfunc
+        0L,   # verbose
+        1L,   # iTask
+        as.double(rwork), # rWork
+        as.integer(iwork), # iWork
+        2L,     # jT full Jacobian calculated internally
+        0L,    # nOut
+        as.integer(lrw), # lRw
+        as.integer(liw), # lIw
+        1L, # Solver
+        NULL, # rootfunc
+        0L, as.double(0), 0L, # nRoot, Rpar, Ipar
+        0L, # Type
+        list(fmat = 0L, tmat = 0L, imat = 0L, ModelForc = NULL), # flist
+        list(), # elist
+        list(islag = 0L) # elag
       )
+
+      out <- t(out_raw)
+      colnames(out) <- c("time", mod_vars)
     } else {
       mkindiff <- function(t, state, parms) {
 
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index 09cb84b8..e193e5e3 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -186,10 +186,24 @@ nlme.mmkin <- function(model, data = "auto",
     thisCall[["control"]] <- control
   }
 
+  # Provide the address of call_lsoda to the fitting function
+  call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve")
+  if (model[[1]]$solution_type == "deSolve" & !is.null(model[[1]]$mkinmod$cf)) {
+    # The mkinmod stored in the first fit will be used by nlme
+    model[[1]]$mkinmod$diffs_address <- getNativeSymbolInfo("diffs",
+      PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address
+    model[[1]]$mkinmod$initpar_address <- getNativeSymbolInfo("initpar",
+      PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address
+  }
+
   fit_time <- system.time(val <- do.call("nlme.formula", thisCall))
   val$time <- fit_time
 
   val$mkinmod <- model[[1]]$mkinmod
+  # Don't return addresses that will become invalid
+  val$mkinmod$diffs_address <- NULL
+  val$mkinmod$initpar_address <- NULL
+
   val$data <- thisCall[["data"]]
   val$mmkin <- model
   if (is.list(start)) val$mean_dp_start <- start$fixed
diff --git a/R/saem.R b/R/saem.R
index 696ea0ee..5b8021de 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -120,12 +120,12 @@ utils::globalVariables(c("predicted", "std"))
 #' summary(f_saem_dfop_sfo, data = TRUE)
 #'
 #' # The following takes about 6 minutes
-#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
-#' #  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+#'   nbiter.saemix = c(200, 80))
 #'
-#' #saemix::compare.saemix(list(
-#' #  f_saem_dfop_sfo$so,
-#' #  f_saem_dfop_sfo_deSolve$so))
+#' #anova(
+#' #  f_saem_dfop_sfo,
+#' #  f_saem_dfop_sfo_deSolve))
 #'
 #' # If the model supports it, we can also use eigenvalue based solutions, which
 #' # take a similar amount of time
@@ -580,6 +580,15 @@ saemix_model <- function(object, solution_type = "auto",
     transform_rates <- object[[1]]$transform_rates
     transform_fractions <- object[[1]]$transform_fractions
 
+    # Get native symbol info for speed
+    call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve")
+    if (solution_type == "deSolve" & !is.null(mkin_model$cf)) {
+      mkin_model$diffs_address <- getNativeSymbolInfo("diffs",
+        PACKAGE = mkin_model$dll_info[["name"]])$address
+      mkin_model$initpar_address <- getNativeSymbolInfo("initpar",
+        PACKAGE = mkin_model$dll_info[["name"]])$address
+    }
+
     # Define the model function
     model_function <- function(psi, id, xidep) {
 
@@ -613,7 +622,8 @@ saemix_model <- function(object, solution_type = "auto",
             odeparms = odeparms, odeini = odeini,
             solution_type = solution_type,
             outtimes = sort(unique(i_time)),
-            na_stop = FALSE
+            na_stop = FALSE,
+            call_lsoda = call_lsoda
           )
 
           out_index <- cbind(as.character(i_time), as.character(i_name))