feat: cppDE 0.11.0, ptc() for steady states, equilibrate relative to each state

- ptc() solves f(x, p) = 0 for a cppFUN() model by pseudo-transient
  continuation: least-squares steps (I/dt - J) d = f with the rows
  C d = total - C x, dt controlled by the local error of the implicit
  Euler step. It follows the flow to a stable steady state and turns into
  Newton near it; flow = FALSE solves plain equations. Positive variables
  are updated multiplicatively and projected onto C x = total. The loop
  runs in C++ on the model's forward entry and LAPACK.
- rootfunc = "equilibrate" stops once |dx/dt| <= roottol |x| + abstol for
  every state, in both backends. The absolute test stopped small states
  too early and large ones too late.

Example in inst/examples/ptc.R, tests in test-ptc.R and
test-ode-equilibrate.R.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: simonbeyer1

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: cppDE
 Title: 'C++' Code Generation for ODEs with Sensitivity Propagation
-Version: 0.10.5
+Version: 0.11.0
 Authors@R: c(
     person("Simon", "Beyer", , "simon.beyer@fdm.uni-freiburg.de", role = c("aut", "cre")),
     person("Lawrence Livermore National Security", role = "cph",

NAMESPACE

@@ -12,6 +12,7 @@ export(diagnostics)
 export(forkGuard)
 export(install_libs)
 export(prepareBatch)
+export(ptc)
 export(solveBatch)
 export(solveODE)
 export(solveODEBatch)

NEWS.md

@@ -1,3 +1,14 @@
+# cppDE 0.11.0
+
+* New `ptc()`: steady state of a `cppFUN()` model by pseudo-transient
+  continuation, with linear constraints `C x = total` and positive variables.
+  Implicit Euler steps under local error control follow the flow to a stable
+  steady state and become Newton steps near it; `flow = FALSE` solves plain
+  equations.
+* `rootfunc = "equilibrate"` stops once every `|dx/dt| <= roottol |x| + abstol`,
+  in both backends. The absolute `|dx/dt| < roottol` stopped small states too
+  early and large ones too late.
+
 # cppDE 0.10.5
 
 * **Bug fix.** `derivMode = "reverse"` and `"forward-reverse"` return the

R/cppODE.R

@@ -11,8 +11,8 @@
 #' variable), `value` (expression to apply), `method` (`"replace"`,
 #' `"add"`, or `"multiply"`), and exactly one of `time` (event time) or
 #' `root` (root expression). `rootfunc` terminates integration:
-#' `"equilibrate"` stops at steady state; a character vector of
-#' expressions stops at the first zero crossing.
+#' `"equilibrate"` stops at steady state (see `roottol` in [solveODE()]); a
+#' character vector of expressions stops at the first zero crossing.
 #'
 #' A time event fires when its time lies between the first and the last
 #' integration time, the first included and the last not; with

R/ptc.R

@@ -0,0 +1,74 @@
+#' Steady state by pseudo-transient continuation
+#'
+#' Solves `f(x, p) = 0` for variables of a [cppFUN()] model, one output per
+#' solved variable, optionally under linear constraints `C x = total`.
+#'
+#' Each step solves `(I/dt - J) d = f` together with the constraint rows in the
+#' least-squares sense. Small `dt` gives implicit Euler steps along
+#' `dx/dt = f`, large `dt` Newton steps; `dt` follows the local error of the
+#' pseudo-time step, so the iteration follows the flow to a stable steady
+#' state and converges quadratically near it. With `flow = FALSE` the steps
+#' follow the Newton flow and any regular root of plain equations is returned.
+#'
+#' @param fun A compiled [cppFUN()] object with forward derivatives.
+#' @param x Named numeric, values of all variables of `fun`. The solved
+#'   variables start there, the others keep their values.
+#' @param parms Named numeric, values of the parameters of `fun`.
+#' @param solve Names of the variables to solve for. Default all.
+#' @param rows Names of the outputs of `fun` paired with `solve`, in the same
+#'   order. Default `solve`.
+#' @param C,total Constraint matrix with columns named by `solve`, and its
+#'   right-hand side: `C x = total`. Default none.
+#' @param flow `TRUE` if the outputs are the right-hand side of an ODE,
+#'   `FALSE` for plain equations.
+#' @param positive Keep the solved variables positive: the step is taken
+#'   multiplicatively and the constraints are restored by `x exp(C' lambda)`.
+#' @param controls Named list: `rtol` (default `1e-10`) and `atol` (`1e-14`),
+#'   every row within `rtol` of its turnover `sum_k |df_i/dx_k| |x_k|` plus
+#'   `atol` times its diagonal; `flowTol` (`1`), relative local error of a
+#'   pseudo-time step; `maxit` (`400`); `dtInit` (`1e-2`), the first
+#'   pseudo-time step, relative to the fastest rate if `flow = TRUE`; `zmax`
+#'   (`5`), the largest log change of a variable per step.
+#' @return A list with the named solution `x` (all variables), `converged`,
+#'   `iterations` and `message`.
+#' @references Kelley CT, Keyes DE (1998). Convergence analysis of
+#'   pseudo-transient continuation. SIAM J Numer Anal 35(2):508-523.
+#' @example inst/examples/ptc.R
+#' @export
+ptc <- function(fun, x, parms, solve = names(x), rows = solve, C = NULL, total = NULL,
+                flow = TRUE, positive = TRUE, controls = list()) {
+  vars <- attr(fun, "variables"); pars <- attr(fun, "parameters")
+  outs <- names(attr(fun, "equations")); mn <- attr(fun, "modelname")
+  if (is.null(vars) || is.null(mn)) stop("'fun' must be a cppFUN() object")
+  if (!"forward" %in% attr(fun, "derivMode")) stop("'fun' needs forward derivatives")
+  if (!all(vars %in% names(x))) stop("'x' misses ", toString(setdiff(vars, names(x))))
+  if (!all(pars %in% names(parms))) stop("'parms' misses ", toString(setdiff(pars, names(parms))))
+  if (!all(solve %in% vars)) stop("unknown variables in 'solve': ", toString(setdiff(solve, vars)))
+  if (length(rows) != length(solve) || !all(rows %in% outs))
+    stop("'rows' must name one output of 'fun' per solved variable")
+  if (is.null(C)) {
+    C <- matrix(0, 0, length(solve)); total <- numeric(0)
+  } else {
+    if (!all(solve %in% colnames(C))) stop("'C' needs a column per solved variable")
+    C <- C[, solve, drop = FALSE]
+    if (length(total) != nrow(C)) stop("'total' needs one value per row of 'C'")
+  }
+  def <- list(rtol = 1e-10, atol = 1e-14, flowTol = 1, maxit = 400L, dtInit = 1e-2, zmax = 5)
+  bad <- setdiff(names(controls), names(def))
+  if (length(bad)) stop("unknown controls: ", toString(bad))
+  ctrl <- utils::modifyList(def, as.list(controls))
+  ctrl$flow <- as.numeric(isTRUE(flow)); ctrl$positive <- as.numeric(isTRUE(positive))
+  ctrl <- lapply(ctrl, as.double)
+  x0 <- as.double(x[vars])
+  if (isTRUE(positive) && any(x0[match(solve, vars)] <= 0))
+    stop("with positive = TRUE the solved variables must start positive")
+
+  sym <- .nativeSym(paste0(mn, "_eval_ad"))
+  if (is.null(sym)) stop("compile '", mn, "' first")
+  addr <- getNativeSymbolInfo(sym$name, PACKAGE = sym$dll)$address
+  res <- .Call(C_cppde_ptc, addr, length(outs), x0, as.double(parms[pars]),
+               as.integer(match(solve, vars) - 1L), as.integer(match(rows, outs) - 1L),
+               matrix(as.double(C), nrow(C)), as.double(total), ctrl)
+  xo <- setNames(x0, vars); xo[solve] <- res$x
+  list(x = xo, converged = res$converged, iterations = res$iterations, message = res$message)
+}

R/solveODE.R

@@ -659,7 +659,9 @@
 #' @param roottol Tolerance for root finding in root-triggered events.
 #'   Default `1e-6`. An event `value` that depends on `time` evaluates the
 #'   firing time directly, so its derivatives inherit this tolerance rather
-#'   than `reltol`.
+#'   than `reltol`. With `rootfunc = "equilibrate"` the integration stops once
+#'   every `|dx/dt| <= roottol |x| + abstol`; the native backend includes the
+#'   sensitivities.
 #' @param maxroot Maximum number of triggers per root event. Default `1`.
 #' @param onFailure How to react when the solver returns a non-zero
 #'   return code. One of `"stop"` (default; raise an error with the solver

inst/examples/example_equilibrate.R

@@ -14,7 +14,8 @@ rhs <- c(
   pA = "k1*A*R-k2*pA"
 )
 
-# Equilibrate - stoppt wenn alle Ableitungen (inkl. Sensitivitaeten) < roottol
+# Equilibrate: stops once every |dx/dt| <= roottol |x| + abstol,
+# sensitivities included
 model <- cppODE(
   rhs = rhs,
   rootfunc = "equilibrate",

inst/examples/ptc.R

@@ -0,0 +1,20 @@
+\donttest{
+if (codegenAvailable()) {
+  ## Enzyme E converts S to P, P returns to S. E + ES and S + ES + P are
+  ## conserved, so the steady state is fixed by the two totals.
+  f <- c(E  = "-kon * E * S + (koff + kcat) * ES",
+         S  = "-kon * E * S + koff * ES + kback * P",
+         ES = " kon * E * S - (koff + kcat) * ES",
+         P  = " kcat * ES - kback * P")
+  enz <- cppFUN(f, variables = names(f), parameters = c("kon", "koff", "kcat", "kback"),
+                modelname = "ptc_enzyme", outdir = tempdir(), compile = TRUE)
+  C <- rbind(Etot = c(E = 1, S = 0, ES = 1, P = 0),
+             Stot = c(E = 0, S = 1, ES = 1, P = 1))
+  ss <- ptc(enz, x = c(E = 1, S = 1, ES = 1, P = 1),
+            parms = c(kon = 1, koff = 0.5, kcat = 0.3, kback = 0.1),
+            C = C, total = c(1, 10))
+  ss$x
+  enz$func(E = ss$x[["E"]], S = ss$x[["S"]], ES = ss$x[["ES"]], P = ss$x[["P"]],
+           kon = 1, koff = 0.5, kcat = 0.3, kback = 0.1)
+}
+}

inst/include/cppde/cppde_dual_math.hpp

@@ -647,6 +647,33 @@ inline double max_abs_all_levels_vec(const State& v) {
   return m;
 }
 
+// steady_all_levels(f, x, tol, atol): every |f| <= tol |x| + atol over a value
+// and all of its tangent slots, with f the time derivative of x. The steady-state
+// check of rootfunc = "equilibrate".
+template<class T>
+inline std::enable_if_t<std::is_arithmetic_v<T>, bool>
+steady_all_levels(const T& f, const T& x, double tol, double atol) {
+  return std::abs(static_cast<double>(f)) <= tol * std::abs(static_cast<double>(x)) + atol;
+}
+
+template<class T, unsigned N>
+inline bool steady_all_levels(const dual<T, N>& f, const dual<T, N>& x, double tol, double atol) {
+  auto& fm = const_cast<dual<T, N>&>(f);
+  auto& xm = const_cast<dual<T, N>&>(x);
+  if (!steady_all_levels(fm.x(), xm.x(), tol, atol)) return false;
+  unsigned nd = fm.size(), nx = xm.size();
+  for (unsigned i = 0; i < nd; ++i)
+    if (!steady_all_levels(fm.d(i), i < nx ? xm.d(i) : T(0), tol, atol)) return false;
+  return true;
+}
+
+template<class State>
+inline bool steady_all_levels_vec(const State& f, const State& x, double tol, double atol) {
+  for (std::size_t i = 0; i < f.size(); ++i)
+    if (!steady_all_levels(f[i], x[i], tol, atol)) return false;
+  return true;
+}
+
 } // namespace cppde
 
 #endif // CPPDE_DUAL_MATH_HPP

inst/include/cppde/cppde_events.hpp

@@ -108,12 +108,14 @@ inline bool direction_matches(double last_val, int direction) {
 // Steady-state termination helper (threshold check, no root-finding)
 // ============================================================================
 
+// Steady state: every |dx/dt| <= tol |x| + atol, on every AD level.
 template<class System, class State, class Time>
-std::function<bool(const State&, const Time&)> make_steady_state_termination(System& sys, double tol) {
-  return [&sys, tol](const State& x, const Time& t) -> bool {
+std::function<bool(const State&, const Time&)> make_steady_state_termination(System& sys, double tol,
+                                                                              double atol) {
+  return [&sys, tol, atol](const State& x, const Time& t) -> bool {
     State dxdt(x.size());
     sys(x, dxdt, t);
-    return cppde::max_abs_all_levels_vec(dxdt) < tol;
+    return cppde::steady_all_levels_vec(dxdt, x, tol, atol);
   };
 }