[wshlavacek]

Path-agreement's BIOMD0000000338|klu_vs_dense went from PASS to EXCEPTION on x86-64 CI after #1302 (run 40565116208). That PASS was two answers agreeing by accident, and one of them was silently wrong. This PR refuses the wrong one by name and re-baselines the row.

What was silently wrong on main

GH #135: where the right-hand side (or a closed-form Jacobian) comes back non-finite, bngsim evaluates it again with every negative concentration set to 0. That is the right value for sqrt or a fractional power read a hair below 0, whose value at 0 is their own limit. It is wrong at a pole.

BIOMD0000000338 (Wajima 2009) has XIIa and K activating each other through v·S·E/(1 + E) with k = 1. So XIIa = K = 0 is an exact, unstable solution, and each species has a pole at −1. The dense LU's rounding tips XIIa below 0, and the model's own laws carry it to the pole. There the right-hand side is −inf, the clamp read XIIa and K at 0, and that held both at −1. The run returned them with no error. KLU keeps XIIa exactly 0, which is the exact solution.

run (dense, rtol 1e-6, atol 1e-9)	platform	main (d1b27d5)	this PR
jacobian="fd", t 0..0.08	arm64	XIIa = −1, K = −0.997, no error	refused, naming XIIa
jacobian="auto", t 0..0.08	arm64	XIIa = −1, K = −0.997, no error	refused, naming XIIa
jacobian="analytical", t 0..0.08	arm64	CV_CONV_FAILURE at t = 0.08	refused, naming XIIa
jacobian="fd", t 0..2 (the CI's run)	x86-64	CV_CONV_FAILURE at t = 0.72; before #1302, XIIa = −1, no error	refused, naming XIIa
jacobian="fd", t 0..2	arm64	XIIa(2) = 0.2246	the same, bit for bit
default (KLU), t 0..2	both	XIIa ≡ 0	the same, bit for bit

On arm64, BIOMD0000000810 with jacobian="fd" does the same thing at the pole of /(Interferon_gamma + k4). There main returned a cell of 3.3·10²⁶ where the default run, which agrees with RoadRunner, has 3.4. It is now refused, naming Interferon_gamma. On x86-64 the same run crosses that pole without landing on it, and returns 9.9·10¹⁵ on main and on this PR alike (under Left alone).

The fix (src/cvode_simulator.cpp)

• The clamp leaves a pole alone. Where the value first read (the right-hand side or a closed-form Jacobian, interpreted or compiled, dense or sparse) has an infinite entry, spare_at_pole looks at the state the solver's step started from. That is CVODE's accepted y_n, recovered exactly from the predicted history array as Σ_k Dky(tn, k)·(−h)^k/k!. CVODE does not report the order mid-step, so it is found by asking for one more derivative until CV_BAD_K, and that log line is dropped. A slot is left unclamped (spared) if three things hold:
  • y_n has it more than 10⁶ of its atol below 0;
  • with every other slot read at 0, its own reading is infinite;
  • it reads finite a relative 2⁻²⁰ to either side. A pole the run has landed on is infinite at that one point; an overflow such as exp(−1/S) stays infinite over a stretch (round 3).
• A run carried there is refused. If the value stays non-finite with those slots spared, and turns finite once they too are read at 0, the callback records why and returns CVODE's unrecoverable code. The run raises "CVODE integration failed … (CV_RHSFUNC_FAIL)" and the message names the concentration, its depth in atols, and the time. Because that is a CVODE failure, a failed closed-form attempt still retries on the difference quotient (GH Four Linux-only test failures, found by #169's first whole-suite run: the GH #176 FD fallback dies at a second crossing, and a steady-state fixture is exactly singular under reference LAPACK #176).
• A NaN takes main's path exactly, at any depth: sqrt, fractional powers, logarithms read below 0. So does an infinity with nothing spared. GH Extract the two forward-sensitivity jump handlers out of CvodeSimulator::run (follow-up to #109) #135's and cvodes: an analytical Jacobian still non-finite after the GH #135 clamp is used anyway, and silently freezes that species' Newton correction #867's test files are main's, unchanged.
• dq_side_point (Difference quotients that step across a jump in a rate law (#1002) #1302's side probes) and stop_would_lose don't let a refusal they record stay on the run, and the warm path clears it per run.

Extra right-hand-side evaluations

None on the ordinary path. The check runs only on a call that read an infinity at a step that started more than 10⁶ atols below 0, and only for a slot whose own reading is infinite. It then re-reads that callback's function once or twice:

• at a pole: 2 reads (both sides finite), and the run is refused;
• at an overflow: 1 read (the first side stays infinite), and the slot is clamped as on main.

Counted on arm64 with an instrumented build:

run	outcome	extra evaluations
338 dense "fd"	refused	+2
338 dense "analytical"	refused	+2
338 dense "auto"	refused	+4 (its analytical attempt, then the "fd" retry)
810 "fd"	refused	+2
round 3's exp(−1/X) extinction, auto / fd / KLU	main's answer, bit for bit	+49 / +26 / +27, against 422 / 422 / 4,957 CVODE RHS evaluations

Wall time on clamp-heavy runs that never meet an infinity is main's. Interleaved best of 7, three rounds, arm64, identical results:

• X' = −X^0.5 extinction (49,156 RHS calls, the clamp firing on most): 21.8–22.1 ms against main's 21.6–21.8 ms;
• MODEL2001200002 at its job: 9.75–9.85 ms against 9.69–9.79 ms.

The difference is the isinf scan of a clamped evaluation.

Census (arm64, main d1b27d5 against 40d7ada)

The rr ODE corpus, 1,323 models at each job's own window and tolerances, in three modes: default, jacobian="fd", and dense "fd" (dense capped at 200 species, as the path-agreement harness caps it).

	runs
answer, bit for bit	3,780
fail on both, same message	160
skipped (dense, > 200 species) or time out on both	25
answer → named refusal	2: BIOMD0000000810 "fd" and dense "fd" (main: max 3.29·10²⁶)
bits differ from the previous day's census of main	2: BIOMD0000000385 and 628, default. Main, this PR, and the PR without the 2⁻²⁰ re-read all give the same bits for both, 3 runs each, today.

• No physical state is refused. Every answer-to-refusal change is a concentration the run carried to a pole.
• SBML Test Suite (1,823 cases): every verdict and max error identical to main.
• BIOMD0000000353 and 554, the corpus's two models with an exp(−c/S) law, are bit-identical to main at their own jobs on all 8 paths tried: auto, fd, analytical, dense, dense fd, KLU, KLU fd, codegen. They are also identical without the 2⁻²⁰ re-read.

Review

Three rounds, each by an independent reviewer on a frozen build: read-only, filing nothing.

• Round 1 (4566562): a HIGH, false refusals of correct runs. The first rule judged "far below 0" against CVODE's predicted state, and the first Newton residual of every step sits exactly there. So a predictor overshooting from a step that started near 0 was refused. Example: S → P; 50·S^2.5/(0.01^2.5 + S^2.5) from S = 100 at bngsim's defaults, where main is right to 7e-9. Bounding the clamp itself also broke extinctions (X' = −k·X^a). Also a LOW: the message named the wrong species. Fixed by judging at the step's accepted start and naming only the cause (bbab179).

• Round 2 (bbab179): three HIGHs from judging depth against a tracked scale:

  • a species that spikes and empties between two output rows, or after an event's dose, was refused;
  • a run continued from an earlier run's in-tolerance undershoot was refused;
  • computing the step start on every clamp fire doubled the cost of clamp-heavy runs.

  Every one was a NaN case, and every case the change moved was an infinity at a pole. So the rule was narrowed to poles, with NaN taking main's path (df42f74).

• Round 3 (df42f74), the last round, named in advance:

  • F1, HIGH: exp(−c/S) overflows to +inf at an extinct species' in-tolerance undershoot and was taken for a pole. Fixed by the 2⁻²⁰ re-read, by the maintainer's choice (40d7ada).
  • F2, HIGH: one species' overflow got a different, NaN-only species refused. Fixed by sparing only a slot whose own reading is infinite (d815371).
  • L1, LOW: a message clause reverted to main's wording.
  • N1, same as main: under Left alone.
  • Round 3's probes on 40d7ada all match main, except the pole cases, where main parks at −1 and this PR refuses by name.

Left alone

• A step that crosses a pole without landing on it returns a trajectory where no solution exists, on main and on this PR alike: A step across a pole of the right-hand side returns a trajectory where no solution exists #1339. BIOMD0000000810 with jacobian="fd" on x86-64 is one: it returns 9.9·10¹⁵ on both builds (Rosetta), where arm64 lands on the pole and is now refused. One species, X' = 10·X/(1+X) − 1e-6 from X = 0, reaches its pole at t* = 1.512. Both builds return X(3) = −1.96e16 ("fd"), +2.33e41 ("analytical", default) and −1.6e6 to −1.9e6 (KLU) on arm64, with no error; x86-64 is the same in kind. Nothing non-finite is ever evaluated, so no clamp rule can see it.
• Domain-edge silent errors, NaN, same as main.
  • X → ; 1 with → Y; sqrt(X) from X = 1 returns X(2) = −1, Y(2) = 2/3 with no error. X goes below 0 by the model's own constant sink, and sqrt(X) has no value past t = 1.
  • Round 1's extinction X' = −k·X^0.05 from X0 = 10 at rtol 1e-6, atol 1e-9 returns X frozen at −0.125, 1.25·10⁴ rtols off.
  • Judging depth for a NaN needs a scale that round 2 showed can't be tracked reliably, so both take main's path.
• A pole that reads NaN, not ±inf, still parks, as on main. Example: 10·S·E/(1+E) − E/(1+E) split into two reactions, so that E's row is inf − inf. Round 3's N1.
• A pole within 10⁶ atols of 0 still parks, as on main, because the band widens with atol. Example: 338 at atol 1e-6, where the band reaches −1.
• BIOMD0000000338's dense runs stay chaotic in where the burst goes, because rounding decides it. That is the model's unstable equilibrium at 0, not a bngsim defect.

Tests

• python/tests/test_clamp_far_below_zero.py (new, 10 tests).
  • Fail on main (d1b27d5, arm64): 338 dense "fd", "auto" and "analytical" over 0..0.08, each either refused naming XIIa (and not K) or returning nothing 10⁶ atols below 0.
  • Guards that pass on main too:
    • 338 on KLU keeps XIIa = K = 0 exactly;
    • 404's failed closed-form attempt still retries on "fd" and returns its answer bit for bit (arm64); on x86-64 the attempt succeeds, within 3e-7 of it;
    • 810 by default is untouched, and 810 "fd" is refused naming Interferon_gamma (arm64) or returns nothing below 0 (x86-64, where it crosses the pole);
    • a Hill predictor overshoot is clamped as before (LSODA reference);
    • an extinction stopping 8.2 rtols of its start below 0 is clamped;
    • one species' overflow does not get another refused (F2).
• Fixtures: tests/data/BIOMD0000000338.xml, BIOMD0000000404.xml, BIOMD0000000810.xml (BioModels, CC0).
• parity_checks/nightly/baselines/differential.json: BIOMD0000000338|klu_vs_dense goes PASS → EXCEPTION, and its tally moves by that one row; no other row changes.
  • The record is copied verbatim from main's scheduled nightly, run 40712357505 at 9a26f5b (nightly-differential, verdicts.json), where it was that night's one alert.
  • verdicts.py diff of the new baseline against that run's verdicts reports no regressions. It reports the same with 338's note replaced by this PR's named refusal (EXCEPTION either way).
  • 810's rows don't move on x86-64, where its "fd" run crosses the pole on both builds.

Local runs

• python/tests with -n 6 on the merged branch: 15,962 passed, 435 skipped, 1 xfailed. 2 failed:
  • test_run_tests_sh: the known environment-only failure;
  • test_cache_cli's prune test: a race under xdist, which passes 3 of 3 alone.
• The clamp tests (test_rhs_nonneg_clamp, test_nonneg_clamp_scope, test_nonfinite_jacobian_after_clamp, test_nonfinite_diagnostic, test_difference_quotient_jump, test_sens_numerically_zero_base, and the new file): 110 passed.
• x86-64 under Rosetta (KLU on), 338 at the CI's settings: dense "fd", "auto" and "analytical" are refused naming XIIa, and KLU is bit for bit.
• clang-format 17.0.6 and ruff 0.14.0 through the pinned pre-commit hooks, plus ci/changelog.py check.