Key warm-start caches per condition in Pequil/Pimpl

Replace the single warm-start cache in Pequil/Pimpl with a per-condition
registry. A condition-less transformation composed across conditions now
warm-starts each condition from its own previous root instead of from
whichever condition was solved last, which could cross basins of
attraction in multistable models. parfn forwards the active condition to
the warm-starting p2p closure.

Flush warm-start caches at the start of every mstrust fit so each
multistart fit rediscovers its own steady-state roots rather than
inheriting a stale basin. This covers both the sequential path, where all
fits share one objfun closure, and forked workers, which inherit the
parent cache via copy-on-write.

Make lpSolve an unconditional dependency since the sink-cluster LP now
always runs; move MASS, pander, and lpSolve from Suggests to Imports.

Update tests for the registry structure and add a per-condition warm-start
test.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: simonbeyer1

DESCRIPTION

@@ -34,13 +34,14 @@ Imports:
     cowplot,
     digest,
     nleqslv,
-    scales
+    scales,
+    MASS,
+    pander,
+    lpSolve
 Remotes:
     simonbeyer1/CppODE,
     dkaschek/cOde
 Suggests:
-    MASS,
-    pander,
     knitr,
     rmarkdown,
     testthat (>= 3.0.0),
@@ -49,8 +50,7 @@ Suggests:
     deSolve,
     JuliaCall,
     openxlsx,
-    Polychrome,
-    lpSolve
+    Polychrome
 VignetteBuilder: knitr
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)

R/classes.R

@@ -296,6 +296,11 @@ parfn <- function(p2p, parameters = NULL, condition = NULL) {
   mappings[[1]] <- p2p
   names(mappings) <- condition
 
+  # p2p of warm-starting transformations (Pequil / Pimpl) accepts a `condition`
+  # argument so it can keep one warm-start cache per condition. Forward it only
+  # when present; Pexpl and other p2p's keep their original signature untouched.
+  p2p_has_cond <- "condition" %in% names(formals(p2p))
+
   outfn <- function(..., fixed = NULL, deriv = TRUE, deriv2 = FALSE, conditions = condition, env = NULL) {
 
 
@@ -310,8 +315,18 @@ parfn <- function(p2p, parameters = NULL, condition = NULL) {
 
     if (is.null(overlap)) conditions <- union(condition, conditions)
 
+    # Per-condition warm-start key: prefer this parfn's own condition, else the
+    # single condition it is currently being evaluated under (set by prodfn when
+    # a condition-less parfn is composed with a condition-specifying one).
+    cond_key <- if (!is.null(condition)) condition[[1]]
+                else if (length(conditions) == 1L) conditions
+                else NULL
+
     if (is.null(overlap) | length(overlap) > 0)
-      result <- p2p(pars = pars, fixed = fixed, deriv = deriv, deriv2 = deriv2)
+      result <- if (p2p_has_cond)
+        p2p(pars = pars, fixed = fixed, deriv = deriv, deriv2 = deriv2, condition = cond_key)
+      else
+        p2p(pars = pars, fixed = fixed, deriv = deriv, deriv2 = deriv2)
     else
       result <- NULL
 

R/parameters.R

@@ -63,6 +63,37 @@ P <- function(trafo = NULL, parameters = NULL, condition = NULL,
 }
 
 
+## Per-condition warm-start registry shared by Pimpl / Pequil / .Pequil_totals.
+## Each of those builds ONE p2p closure even when condition-less, and the common
+## usage pattern is a condition-less `Pequil` (or `Pimpl`) composed via `*` with a
+## condition-specifying `Pexpl`. The prodfn loop then calls that single p2p once
+## per condition. With a single cache env this means every condition warm-starts
+## from whichever condition was solved last (order-dependent, can cross basins of
+## attraction). The registry keeps one cache env PER condition key instead, so a
+## condition is always warm-started from its OWN previous root. `parfn()` forwards
+## the active condition to p2p (see the `condition` argument there); a NULL/empty
+## key (parfn called outside any condition context) falls back to a shared slot.
+.warmstart_registry <- function() {
+  caches <- new.env(parent = emptyenv())
+  get_cache <- function(key) {
+    k <- if (is.null(key) || !nzchar(key)) "__default__" else key
+    cc <- get0(k, envir = caches, inherits = FALSE)
+    if (is.null(cc)) { cc <- new.env(parent = emptyenv()); assign(k, cc, envir = caches) }
+    cc
+  }
+  reset <- function() {
+    nms <- ls(caches, all.names = TRUE)
+    for (k in nms) {
+      cc <- get(k, envir = caches, inherits = FALSE)
+      inner <- ls(cc, all.names = TRUE)
+      if (length(inner)) rm(list = inner, envir = cc)
+    }
+    invisible(nms)
+  }
+  list(get = get_cache, reset = reset, caches = caches)
+}
+
+
 #' Parameter transformation (explicit, algebraic)
 #'
 #' Builds `p_inner = f(p_outer)` from symbolic expressions via
@@ -338,7 +369,7 @@ Pexpl <- function(trafo, parameters = NULL, attach.input = FALSE, condition = NU
 #'     vanishes at SS; if a feeder's rate involves exactly one state, that
 #'     state must be zero.
 #'   \item *Sink cluster (LP)*: subsets whose combined mass leaks
-#'     monotonically (mass-balance LP via `lpSolve::lp`, only when installed).
+#'     monotonically (mass-balance LP via `lpSolve::lp`).
 #' }
 #' For each zero-state the column is dropped, the state symbol is substituted
 #' by `"0"` in remaining rates, structurally-zero reactions are removed, and
@@ -403,7 +434,6 @@ Pexpl <- function(trafo, parameters = NULL, attach.input = FALSE, condition = NU
   }
 
   .sink_cluster <- function(M, eps = 1e-8, Mbig = 1e4) {
-    if (!requireNamespace("lpSolve", quietly = TRUE)) return(integer(0))
     nF <- nrow(M); nS <- ncol(M)
     if (nF == 0L || nS == 0L) return(integer(0))
     c_obj <- colSums(M)
@@ -805,8 +835,7 @@ Pimpl <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
     array(c(H4), dim(H4)[2:4], dimnames = dimnames(H4)[2:4])
   }
 
-  cache <- new.env(parent = emptyenv())
-  cache$guess <- NULL
+  reg <- .warmstart_registry()
 
   ## Biological SS magnitudes span ~10 orders when rate constants span 2-3:
   ## default log-uniform sampling over [1e-5, 1e5] when positive = TRUE.
@@ -962,12 +991,13 @@ Pimpl <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   }
 
 
-  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE) {
+  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE, condition = NULL) {
     if (deriv2 && !emit_d2)
       stop("Pimpl was built with deriv2 = FALSE; rebuild with deriv2 = TRUE.", call. = FALSE)
     if (!emit_d1) deriv <- FALSE
     if (deriv2 && !deriv) deriv <- TRUE
 
+    cache <- reg$get(condition)
     p   <- pars
     dP  <- attr(p, "deriv")
     dP2 <- if (deriv2) attr(p, "deriv2") else NULL
@@ -1086,9 +1116,9 @@ Pimpl <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   attr(p2p, "modelname")   <- modelname
   attr(p2p, "compileInfo") <- collectCompileInfo(PEval$func, PEval$jac, PEval$hess)
   attr(p2p, "resetWarmStart") <- local({
-    cache_ref <- cache; mn <- modelname; cond <- condition
+    reg_ref <- reg; mn <- modelname; cond <- condition
     function() {
-      cache_ref$guess <- NULL
+      reg_ref$reset()
       paste0("Pimpl(", mn, if (!is.null(cond)) paste0(":", cond) else "", ")")
     }
   })
@@ -1198,21 +1228,21 @@ Pimpl <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   controls <- c(list(keep.root = keep.root, attach.input = attach.input,
                      start.time = start.time, end.time = end.time), ode_ctrl)
 
-  cache <- new.env(parent = emptyenv())
-  cache$yini <- cache$last_hash <- cache$last_result <- NULL
+  reg <- .warmstart_registry()
 
   default_sens <- matrix(0, n_dep, length(all_sens), dimnames = list(dependent, all_sens))
   if (length(pivots)) default_sens[cbind(pivots, pivots)] <- 1
   default_sens2 <- if (emit_d2)
     array(0, c(n_dep, length(all_sens), length(all_sens)),
           dimnames = list(dependent, all_sens, all_sens)) else NULL
 
-  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE) {
+  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE, condition = NULL) {
     if (deriv2 && !emit_d2)
       stop("Pequil(deriv2 = TRUE) requires the model to be built with deriv2 = TRUE.",
            call. = FALSE)
     if (!emit_d1) deriv <- FALSE
     if (deriv2 && !deriv) deriv <- TRUE
+    cache <- reg$get(condition)
     p   <- pars
     dP  <- attr(p, "deriv")
     dP2 <- if (deriv2) attr(p, "deriv2") else NULL
@@ -1352,9 +1382,9 @@ Pimpl <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   attr(p2p, "modelname")   <- modelname
   attr(p2p, "compileInfo") <- collectCompileInfo(model, model_s, model_s2)
   attr(p2p, "resetWarmStart") <- local({
-    cache_ref <- cache; mn <- modelname; cond <- condition
+    reg_ref <- reg; mn <- modelname; cond <- condition
     function() {
-      cache_ref$yini <- cache_ref$last_hash <- cache_ref$last_result <- NULL
+      reg_ref$reset()
       paste0("Pequil(", mn, if (!is.null(cond)) paste0(":", cond) else "", ")")
     }
   })
@@ -1465,9 +1495,7 @@ Pequil <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
                               maxprogress = 100L, hini = 0, roottol = 1e-6, maxroot = 1L),
                          controlsODE)
 
-  cache <- new.env(parent = emptyenv())
-  cache$yini <- cache$sensini <- cache$sens2ini <-
-    cache$last_hash <- cache$last_result <- NULL
+  reg <- .warmstart_registry()
 
   default_sens <- matrix(0, n_dep, length(all_sens),
                          dimnames = list(dependent, all_sens))
@@ -1480,13 +1508,14 @@ Pequil <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   controls <- c(list(keep.root = keep.root, attach.input = attach.input,
                      start.time = start.time, end.time = end.time), ode_ctrl)
 
-  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE) {
+  p2p <- function(pars, fixed = NULL, deriv = TRUE, deriv2 = FALSE, condition = NULL) {
     if (deriv2 && !emit_d2)
       stop("Pequil(deriv2 = TRUE) requires the model to be built with deriv2 = TRUE.",
            call. = FALSE)
     if (!emit_d1) deriv <- FALSE
     if (deriv2 && !deriv) deriv <- TRUE
 
+    cache <- reg$get(condition)
     p   <- pars
     dP  <- attr(p, "deriv")
     dP2 <- if (deriv2) attr(p, "deriv2") else NULL
@@ -1658,10 +1687,9 @@ Pequil <- function(trafo, parameters = NULL, forcings = NULL, condition = NULL,
   attr(p2p, "modelname")   <- modelname
   attr(p2p, "compileInfo") <- collectCompileInfo(model, model_s, model_s2)
   attr(p2p, "resetWarmStart") <- local({
-    cache_ref <- cache; mn <- modelname; cond <- condition
+    reg_ref <- reg; mn <- modelname; cond <- condition
     function() {
-      cache_ref$yini <- cache_ref$sensini <- cache_ref$sens2ini <-
-        cache_ref$last_hash <- cache_ref$last_result <- NULL
+      reg_ref$reset()
       paste0("Pequil(", mn, if (!is.null(cond)) paste0(":", cond) else "", ")")
     }
   })

R/statistics.R

@@ -954,6 +954,19 @@ mstrust <- function(objfun, center, studyname, rinit = .1, rmax = 10, fits = 20,
       argstrust[["traceFile"]] <- file.path(resultFolder, paste0(formatC(i, digits = digits, flag = "0"), "_", argslist[["traceFile"]]))
     }
 
+    # Each fit must rediscover its own condition-specific steady-state roots.
+    # A stale Pequil/Pimpl warm-start root would seed this fit's equilibration
+    # and short-circuit the parfn-internal multistart, silently pinning the
+    # wrong basin in multistable models. This body runs once per fit in BOTH
+    # backends, so flushing here covers both leak paths:
+    #   - cores = 1 (%do%): all fits share one objfun closure; clear the root
+    #     left by the previous fit.
+    #   - cores > 1 (%dopar%): the body runs inside the fork, so this clears the
+    #     fork's copy-on-write inheritance of whatever the parent cache held at
+    #     fork time (non-empty e.g. if objfun was evaluated before mstrust),
+    #     which would otherwise be the SAME stale root in every fork.
+    try(resetWarmStarts(objfun, verbose = FALSE), silent = TRUE)
+
     # Retry loop: a try-error or fit$error triggers re-sampling parinit and
     # re-running optmethod, up to nTries times. parframe-supplied centers
     # are skipped (rows are taken as given). Warm-start caches in

tests/testthat/test-Pequil-Pimpl.R

@@ -454,7 +454,9 @@ test_that("resetWarmStarts clears Pequil's cache by name", {
   pars <- c(k_in = 1.5, k_out = 0.3, A = 0.1)
   pf(pars)
 
-  reset_env <- environment(attr(pf, "resetWarmStart"))$cache_ref
+  # Warm-start caches are now kept per condition in a registry; a condition-less
+  # call lands in the "__default__" slot.
+  reset_env <- environment(attr(pf, "resetWarmStart"))$reg_ref$caches[["__default__"]]
   expect_false(is.null(reset_env$yini))
 
   labels <- resetWarmStarts(pf, verbose = FALSE)
@@ -475,7 +477,7 @@ test_that("resetWarmStarts clears Pimpl's cache by name", {
               compile = TRUE, verbose = FALSE)
   pf(c(a = 1, x = 0.5))
 
-  reset_env <- environment(attr(pf, "resetWarmStart"))$cache_ref
+  reset_env <- environment(attr(pf, "resetWarmStart"))$reg_ref$caches[["__default__"]]
   expect_false(is.null(reset_env$guess))
 
   labels <- resetWarmStarts(pf, verbose = FALSE)
@@ -501,8 +503,8 @@ test_that("resetWarmStarts walks into composed functions", {
   p1(c(k1 = 1, kdg1 = 0.5, A = 0.1))
   p2(c(k2 = 2, kdg2 = 0.4, B = 0.2))
 
-  ref1 <- environment(attr(p1, "resetWarmStart"))$cache_ref
-  ref2 <- environment(attr(p2, "resetWarmStart"))$cache_ref
+  ref1 <- environment(attr(p1, "resetWarmStart"))$reg_ref$caches[["__default__"]]
+  ref2 <- environment(attr(p2, "resetWarmStart"))$reg_ref$caches[["__default__"]]
   expect_false(is.null(ref1$yini))
   expect_false(is.null(ref2$yini))
 
@@ -515,6 +517,40 @@ test_that("resetWarmStarts walks into composed functions", {
 })
 
 
+test_that("a condition-less Pequil keeps an independent warm start per condition", {
+  skip_if_no_compile()
+  oldwd <- setwd(.dmod_fx_workdir()); on.exit(setwd(oldwd), add = TRUE)
+
+  stamp <- as.integer(Sys.time())
+
+  # Condition-less equilibration: SS A* = k_in / k_out.
+  g_eq <- Pequil(c(A = "k_in - k_out * A"),
+                 parameters = c("k_in", "k_out"),
+                 modelname = paste0("test_percond_eq_", stamp),
+                 compile = TRUE, verbose = FALSE, attach.input = TRUE)
+
+  # Two conditions with different k_in -> different steady states (3 and 6).
+  trafos <- list(
+    C1 = c(k_in = "s",     k_out = "1", A = "1"),
+    C2 = c(k_in = "2 * s", k_out = "1", A = "1"))
+  px <- P(trafos, modelname = paste0("test_percond_px_", stamp),
+          compile = TRUE, verbose = FALSE)
+
+  pf  <- g_eq * px
+  out <- pf(c(s = 3), deriv = FALSE)
+
+  expect_equal(as.numeric(out$C1["A"]), 3, tolerance = 1e-4)
+  expect_equal(as.numeric(out$C2["A"]), 6, tolerance = 1e-4)
+
+  # The single Pequil closure must now hold one warm-start slot per condition,
+  # each cached at its own root -- not a single shared/overwritten cache.
+  caches <- environment(attr(g_eq, "resetWarmStart"))$reg_ref$caches
+  expect_setequal(ls(caches), c("C1", "C2"))
+  expect_equal(as.numeric(caches[["C1"]]$yini["A"]), 3, tolerance = 1e-4)
+  expect_equal(as.numeric(caches[["C2"]]$yini["A"]), 6, tolerance = 1e-4)
+})
+
+
 test_that("resetWarmStarts is idempotent and silent when nothing to reset", {
   f <- function(x) x + 1
   labels <- resetWarmStarts(f, verbose = FALSE)

tests/testthat/test-normL2.R

@@ -350,8 +350,10 @@ test_that("normL2 rejects unknown opt.BLOQ values", {
   bench <- fx_decay_compiled()
   data  <- fx_decay_data_bloq(sigma = 0.05, lloq = 0.1,
                               times = seq(0, 10, by = 1))
+  # match.arg() rejects the unknown value; assert on the listed choices rather
+  # than the "should be one of" prefix, which match.arg translates per locale.
   expect_error(normL2(data, bench$prd_id, opt.BLOQ = "M2"),
-               "should be one of")
+               "M4BEAL")
 })