Skip to content

differential: steady_state() against a run to equilibrium (#1311) - #1320

Open
Montana wants to merge 10 commits into
lanl:mainfrom
Montana:fix/1311-steady-vs-run
Open

Montana wants to merge 10 commits into
lanl:mainfrom
Montana:fix/1311-steady-vs-run

Conversation

@Montana

@Montana Montana commented Oct 27, 2026 •

Copy link
Copy Markdown
Contributor

Fixes #1311.

Stacked on #1321. This PR's first nightly failed Path-agreement on TIMEOUTs only: models at the edge of the 900 s cap (MODEL1006230049, 0053, 0077) lost every row once a check was added to every model. #1321 bounds each check by the model's cap and keeps a killed worker's finished rows. It is merged in here, and this PR's new checks run after every check main runs. Until #1321 merges, this PR also shows its commits.

Builds on #1291 (merged). It adds two rows to CHECKS (steady_vs_run, newton_vs_run), as #1315, #1316, #1318 and #1319 each add theirs, so whichever lands second keeps both. I'll resolve it.

What it adds

steady_state() has its own routes: the CVODE march, the KINSOL polish, the scaled criterion's run test (#1001), the stability certificate (#78). None of them was compared with the plain integrator's answer. Two checks now do it, one per method, against the same run:

check path reference
steady_vs_run steady_state() (integration) a run() to t = 1e6 (steady_state()'s own max_time) from the same initial state
newton_vs_run steady_state(method="newton") the same run

Both are passed the job's rtol/atol. Write-only accumulators (Model.is_pure_sink()) are masked from the solve and left out of the comparison, since they settle nowhere (#74).

Judged only where both have arrived, as the issue asks:

  • the solve's converged;
  • the run's last quarter of the horizon moving it by no more than 1e3·rtol on its own scale.

Otherwise the row is SKIP, saying why. The bound is max(1e3·rtol, 10·max(settle, drift)), where drift is the run's end at a tenfold tighter tolerance, the harness's usual rule.

Where the two differ, a run from the solve's state, and its spectrum, say which case it is:

the solve's state verdict
a run from it leaves it DIFF: a converged that is not
off the initial state's conserved totals (conservation_laws), which the run keeps DIFF: another system's root
an unstable root (eigvals(model.jacobian(x)), independent of the certificate) that the solve moved to DIFF
the initial state itself, an unstable root a run leaves on rounding alone SKIP, said why
another state a run rests at SKIP, both states named (the species that differs most, with both values), for the baseline to hold: the model's multistability, or a tail slower than the horizon

A refusal by name (an event, time(), a conservation law the RHS breaks: #710, #758) is UNSUPPORTED. A reference run that fails is SKIP ("no reference"): the parity suites own run failures.

 CHECKS = (
     ...
     "grid",
+    "steady_vs_run",
+    "newton_vs_run",
     "sparsity_covers_jacobian",
     "backend_as_asked",
 )

Over the corpus

Local run, --checks steady_vs_run,newton_vs_run, fd_sens's 53 .net models and 1,300 BioModels (2,706 rows):

PASS SKIP UNSUPPORTED DIFF EXCEPTION TIMEOUT
fd_sens 68 36 2 0 0 0
BioModels 1,049 829 722 0 0 0

No finding. The first full run read 2 DIFF and 6 EXCEPTION. All eight were the reference's, and the second commit reclassifies them. I reran those 5 models on the final code; the rest of the corpus is unchanged by it:

  • BIOMD0000000518 (DIFF on both methods): the initial state is itself a root (RHS ~1e-14), unstable along one mode (+5.6e-4). The solve stays on it, as the trajectory does in exact arithmetic. The run leaves it on rounding (e^558 over the horizon) and comes to rest at N2 = 15.5; a run started 1e-6 away rests at a third state, 46.4. This is steady_state()'s documented contract, not a bug.
  • BIOMD0000000609, BIOMD0000001041, MODEL1112050000, MODEL2003040001 (EXCEPTION): the run to 1e6, or its tenfold tighter twin, fails in CVODE (non-finite RHS, or no progress) far past the model's own window.

Two states, recorded as SKIP: BIOMD0000000072, BIOMD0000000546, MODEL1911120001 (Newton). On 72 Newton returns L = 0 while the run holds L at 6.02e17. L is not conserved (RL degrades), so L = 0 is the true limit, about 1e14/k away. The run settles only relative to its own size, and integration rightly reports unconverged. Hence the note's "or a tail slower than the horizon".

Cost: 1,663 s of CPU over the 1,300 BioModels (median 0.43 s per model, worst 155 s, MODEL1006230053), base run included. That is about 7 minutes on the nightly's 4 workers, against the job's 52 of 180.

Tests (python/tests/test_differential_shapes.py, appended)

New models, each called through check_steady directly rather than added to SHAPES, so they don't pay for every other check's compiles:

test expects
both methods on the cubic from x0 = 0.5 PASS
negative control: criterion="residual", #1001's defect DIFF: integration's state a run leaves (5.0e-10), Newton's unstable root 0.3 (a run from it stays, so the spectrum is what sees it)
the same wrong state, reported unconverged SKIP
a solve resting on the cubic's other stable root, 0 SKIP, "two stable states", X() at 0 against 1e-09
a stable root of A <-> B with the total off (2 for 1) DIFF, conservation law 0
the logistic's 0, which no run from it leaves DIFF, unstable root
the cubic started on 0.3, the run leaving it on rounding SKIP, initial state
the run to equilibrium raising SKIP, no reference
the Brusselator, converged planted SKIP, "has not settled"
the accumulator PASS (mask and comparison both leave it out)
x0 = 0.2 PASS (both reach 0)

The four event and time-reading shapes are declared refused by name in REFUSED.

Mutations, each caught:

  • no settle gate;
  • no eigenvalue guard;
  • no rest test;
  • no converged gate;
  • method ignored;
  • refusals not recognised;
  • the rest test run from the initial state;
  • no mask;
  • the run reading every species;
  • a 1e9·rtol bound;
  • no conservation guard;
  • no initial-root SKIP;
  • a failed reference raising.

Suites: parity_checks/tests and the shapes give 1,127 passed, 28 skipped, 1 xfailed. ruff 0.14.0 is clean. The branch merges cleanly with lanl main.

Still to do

The baseline (parity_checks/nightly/baselines/differential.json) gains these rows from this PR's own nightly run. I'll record them from its nightly-differential artifact, as #1291 did.

steady_vs_run and newton_vs_run compare each steady_state() method with
a run to t=1e6 from the same state, judged only where the solve says it
converged and the run's last quarter moves it by at most 1e3 rtol. Where
they differ, a state a run leaves, one off the initial state's conserved
totals, or an unstable root is a DIFF; another state a run rests at is
SKIP, both named, for the baseline to hold.
… no finding (lanl#1311)

The corpus run's 2 DIFFs and 6 EXCEPTIONs were the reference's. Four
models' runs to t=1e6 fail (CVODE) far past their own window: SKIP, no
reference. BIOMD0000000518 starts on an unstable root (an eigenvalue of
5.6e-4): the solve stays on it, as the trajectory does in exact
arithmetic, and the run leaves it on rounding alone.
…inished

The pool kills a model's worker at its cap (900 s on the nightly), and
every row of the model with it. MODEL1006230049, 0053 and 0077 take
716-751 s of that on main, so any check a pull request added to every
model turned all their rows TIMEOUT: lanl#1315, lanl#1318 and lanl#1319 each read
27-36 regressions, none of them a finding.

check_all now takes a deadline, 0.9 of the cap. Every run has until it
(run(timeout=)), no run or check starts past it, and a check that does
not finish is its own TIMEOUT, said why. Each verdict is written as it is
reached, so a worker the pool still kills (one compiling C past the
deadline, which no run's timeout bounds) leaves the runner every row it
finished; only the rest are TIMEOUT.
…-vs-run

steady_vs_run and newton_vs_run run after every check main runs, the
runs to equilibrium bounded by what is left of the cap. A reference run
the cap stops is this check's TIMEOUT, not a missing reference;
steady_state() takes no timeout, and a kill still leaves the other rows.
With lanl#1321, PRs stacked on this one still read 11-12 regressions: main's
own checks on MODEL1006230049 and 0053 did not fit the 810 s deadline of
a 900 s cap. Across five runs of the same checks those models took
689-894 s, so the cap flipped their rows with runner speed. 1500 s puts
the deadline at 1350, and costs only the few models that use it.
…rash

With the 1500 s cap this PR's nightly read 2 new crashes, both
BIOMD0000000385: given the time to reach its codegen compile, 4.2 MB of C
ran past bngsim's 600 s compile limit, which codegen=True reports as a
RuntimeError. Those rows were TIMEOUT at 900 s, and a timeout they stay:
one anywhere in the cause chain, or in codegen=True's refusal, which
carries the compile's timeout in its text and not as its cause. The other
five rows of the model went TIMEOUT -> PASS.
… fix/1311-steady-vs-run

# Conflicts:
#	python/tests/test_differential_shapes.py
…n EXCEPTION (lanl#1311)

Since lanl#1337, Model.jacobian() refuses a difference quotient whose step
decides an entry. check_steady asked it for the spectrum that tells an
unstable root from another stable state, and _guard turned the refusal
into EXCEPTION, a finding against steady_state(). Without a spectrum the
stability is not judged: SKIP, said why, with both states named.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

differential: steady_state() against a long run to equilibrium

2 participants