[elyfmiller]

Part of #907: cases 5 and 3. Case 1 was #1097. Cases 2 and 4 are left, so this doesn't close the issue.

Problem

A plain run restarts CVODE at a state-switch root only where the crossing is resolved: the solver's interpolant must go past the surface by more than the band the tolerance leaves the residual, times 10·√N (#903). The band was

δg = Σ_j |∂g/∂x_j|·(rtol·|x_j| + atol_j)

with atol_j read from the run's static tolerance, and the species summed one by one. Two of the issue's cases come from that. The review that found them kept its probes locally, so they were rebuilt here from the issue's descriptions: a window that opens where a residual crosses and that a clock (T′ = 1) closes 0.5 later, r() = if((Obs > thr) && (Tobs < 1), k, 0), with exact Y = k·0.5 = 5.

• Case 5, a tracking atol. Under TrackingAtol(decades=6, ceiling=1e-8), CVODE holds S, about 1e-6, to 2e-14. The band read the 1e-8 ceiling, though, so a window 5e-9 wide in S read as unresolved.
• Case 3, an observable over many species. Obs sums 1,000 species of 1e-3 each and moves 0.1% of itself, 1e-3, across the window. Its band was rtol·Obs + 1000·atol ≈ 1e-5, and 3.2e-3 with the factor, past the whole excursion.

On main (aedadd8), rtol = atol = 1e-8, five phases of the window:

	Y/exact, phases 5.0, 5.37, 6.13, 7.71, 9.19
case 3, Obs over 1 species (control)	1 1 1 1 1
case 3, Obs over 1,000 species	1 0 0 0 1
case 5, tracking atol	1 0 0 0 1
case 5, scalar atol 1e-8 (control: the window is inside the tolerance)	1 0 0 0 1

Fix

Each species' term in the band is now what CVODE's own error-weight rule holds it to at the state read, and the terms add in quadrature:

δg = sqrt(Σ_j (|∂g/∂x_j|·u_j)²),   u_j = 1/ewt_j(x) = rtol·|x_j| + atol_j

• tolerance_unit (include/bngsim/atol_vector.hpp) gives u_j:

  • rtol·|x| + atol_j, the same expression as before, under a scalar or vector atol;
  • the denominator of fill_tracking_ewt (rtol·|x| + clamp(rtol·|x|, lo, hi)) under a tracking atol.

  It is evaluated at the state of each reading, as the old band was. Under tracking it is also floored at the weights CVODE tested the step with (CVodeGetErrWeights). A tracking atol moves with the state, so where a species falls within a step, the rule at the reading's state alone could hold it tighter than CVODE did. Without tracking the atol floors both, and the rule alone is main's band.

• Quadrature is what the WRMS error test bounds: with ‖e‖_WRMS ≤ 1, |Σ_j ∂g/∂x_j·e_j| ≤ √N·sqrt(Σ_j (∂g/∂x_j·u_j)²). The old sum is the looser bound, which made an observable's band grow with its number of species. state_switch_band scales by the largest term, so no square under- or overflows, and a residual of one species reads |∂g/∂x|·u to the bit. A term that isn't a finite number gives the sum, as before.

• Used by:

  • the restart test (flow_carries_state_switch): its readings, the root's band where the halvings stop, and each species' reach;
  • the A .net rate law that slides along its state switch runs up to 12x slower since #903 (60 s timeout where 6.7 s sufficed) #908 hold. It stands for half the test's band, so its band at the root, its box, and shown() in the root function read the same way.

  The Forward sensitivity is silently wrong once a rate law makes the state slide along its switching surface (dS/damp = 1 against 0) #926 slide detector keeps its own band. It runs only in sensitivity runs, where the restart test and the hold do not.

• Who else sees it: any residual with two or more nonzero terms (A < B, A + B > c) now gets a bar up to √n narrower, so it can restart where main stepped on. A .net rate law that slides along its state switch runs up to 12x slower since #903 (60 s timeout where 6.7 s sufficed) #908's holds engage a little less on such residuals: the gate |g_root| < band is harder to meet, and the band is narrower. A residual with one nonzero term reads the same bar to the bit.

• Not changed: √N over all species (case 2) and the ×10 margin (case 4). Both change what 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's declined crossings rest on, and they are left for their own change, as cvode: sample the restart test's step halving toward the root as well (#907) #1097 said.

Results

The same table, this branch:

	Y/exact
case 3, Obs over 1 species	1 1 1 1 1
case 3, Obs over 1,000 species	1 1 1 1 1
case 5, tracking atol	1 1 1 1 1
case 5, scalar atol 1e-8	1 0 0 0 1, unchanged: a 5e-9 window inside atol 1e-8 is one the tolerance cannot decide

GH #176's fixture and its five twins step exactly as main. That is the split, ramp, latch, damped and moving variants of _ltype_variant, at rtol 1e-6, 1e-8 and 1e-10 (atol = rtol, t ∈ [0, 150], 301 points), with the FD and the automatic Jacobian. All 36 runs are identical to main in n_steps, n_rhs_evals, n_jac_evals and n_root_evals, and in the observables' bytes. The same holds with a vector atol of the same values: 36 of 36. "split" reads two species, so its band goes through the quadrature.

twin	FD steps at rtol 1e-6 / 1e-8 / 1e-10 (main = this branch)
fixture	fails on FD at 1e-6 (both) / 588 / 1574
split	669 / 1065 / 1688
ramp	fails on FD at 1e-6 (both) / 557 / 1314
latch	fails on FD at 1e-6 (both) / 1630 / 2324
damped	616 / 976 / 2548
moving	281 / 526 / 1112

Tests

python/tests/test_net_state_switch_stepping.py gains 30 tests, through a _clock_closed_window builder:

test	main	here
test_a_window_on_an_observable_over_many_species_is_integrated, n = 1000, 5 phases × t_end 20 and 50	6 fail	pass
the same, n = 1 (control)	pass	pass
test_a_window_under_a_tracking_atol_is_integrated, 5 phases × t_end 20 and 50	6 fail	pass

The 12 that fail on main are the lost phases above, at both run lengths. Whether main fails one depends on where the step that spans the window ends, which is the step-geometry caveat the issue raises about case 1's test. Hence the two run lengths.

Test audit

Checked against the test-writing rules Bill sent: each expected value's oracle, each tolerance's reason, and each test shown to fail against broken code. Commit cd3d777.

• Oracle: Y = k·0.5 = 5. The model fixes the window's edges: Obs reaches thr at t_open, and T reaches 1 at t_open + 0.5. A lost window gives Y = 0. The derivation now sits beside the value.
• Tolerance: rel 1e-6, now with its reason. CVODE holds Y to rtol 1e-8, and the 30 runs are at most 3.1e-8 off (measured), so 1e-6 leaves 30×.
• Mutations, each rebuilt and run:

mutation	caught by
the band's terms added as a plain sum (ℓ1)	the 6 case-3 runs over 1,000 species that fail on main
the unit ignores the tracking atol	the 6 case-5 runs that fail on main
the floor at CVODE's error weights dropped	nothing: 128 of 128 pass (see Not covered)

Other checks:

• Python suite (macOS arm64, 10 workers): 15,902 passed, 464 skipped, 1 xfailed. 1 failed: test_run_tests_sh.py::test_the_script_itself_runs_walk_up_dependent_tests_green, which fails the same way on main in this environment, because the venv has no jax, diffrax or roadrunner and the stand-in run reports their modules as skips.
• Formatting: clang-format 17.0.6 and ruff 0.14.0 are clean. The changelog fragment, changelog.d/907.fixed.2.md, sits beside cvode: sample the restart test's step halving toward the root as well (#907) #1097's 907.fixed.md.
• perf_ab: benchmarks/perf_ab.py --quick, main aedadd8 against this branch (21 models up to 150 species, ode and codegen). No significant slowdown, and the solver's work counters are the same on all 42 rows. The geometric-mean ratios are 0.982 warm and 1.002 cold for ode, and 0.980 warm and 0.995 cold for codegen. Two of the 21 networks root a state switch in a rate law, oscillatory_system (1) and gene_expression (3), and both take main's steps exactly. None of the seven larger models in the full suite has one.

Review

There was no second session to review this, so a separate review agent read the diff, read-only. It found no correctness bug. It confirmed:

• bit-identity for one-term residuals, with the same expression and contraction, and v/v = 1;
• that the ℓ2 times √N bar is the exact bound from CVODE's WRMS test;
• that the hold still stays under two thirds of the test's band;
• that the tracking pointer's lifetime, the run paths and the slide detector are kept apart.

Its low findings are all fixed here:

• the tracking unit is floored at the weights the step was tested with;
• the changelog's and comment's arithmetic is made explicit, with band and bar told apart;
• the changelog says who else sees the change;
• a second run length widens the tests;
• a null guard is added on the hold's tracking pointer;
• the helper's comment is moved and rewrapped.

Not covered

• Cases 2 and 4 of A state window a restart would catch is still stepped over when #903's restart test calls the crossing unresolved (16-point sampling, 10·√N band) #907.
• The Forward sensitivity is silently wrong once a rate law makes the state slide along its switching surface (dS/damp = 1 against 0) #926 slide detector's band, which still reads the static atol and the sum.
• A hold under a tracking atol, and a slide over many species: no test here exercises them.
• The floor at CVODE's error weights under a tracking atol. Removing it fails no test.
  • Why it's hard to test: the floor moves the restart bar by at most how far a species' tracked atol changes within one step, well inside the ×10 margin.
  • What a test would need: a model tuned to sit in that sliver, which only this change's own reasoning could say the right answer for.
• Linux and Windows, and the SBML semantic suite gate: CI runs them.