Data: nightly parity reference archive (reference-data-v3)

Reference-engine outputs for bngsim's nightly parity check (issue #805), recorded by https://github.com/lanl/bngsim/actions/runs/40899940910 at e425546.

Not a bngsim release. One asset per suite; each holds a README, entries.json (every entry's model, source and license) and, where any entry is not CC0, bngsim's NOTICE. The nightly replays from these once parity_checks/nightly/reference_archive.json pins them; parity_checks/nightly/README.md has the procedure.

Published by hand from that run's assets. The run's publish job was refused (HTTP 403) when it tried to create the tag at the recorded commit, so this tag points at main's head when it was published. The recorded commit is the one named above, and the manifest pins it.

Data: nightly parity reference archive (reference-data-v2)

Reference-engine outputs for bngsim's nightly parity check (issue #805), recorded by https://github.com/lanl/bngsim/actions/runs/40404704911 at 104078d.

Not a bngsim release. One asset per suite; each holds a README, entries.json (every entry's model, source and license) and, where any entry is not CC0, bngsim's NOTICE. The nightly replays from these once parity_checks/nightly/reference_archive.json pins them; parity_checks/nightly/README.md has the procedure.

Data: nightly parity reference archive (reference-data-v1)

Reference-engine outputs for bngsim's nightly parity check (issue #805), recorded by https://github.com/lanl/bngsim/actions/runs/40125305877 at 2b5a602.

Not a bngsim release. One asset per suite; each holds a README, entries.json (every entry's model, source and license) and, where any entry is not CC0, bngsim's NOTICE. The nightly replays from these once parity_checks/nightly/reference_archive.json pins them; parity_checks/nightly/README.md has the procedure.

Data: BioModels SBML for rr_parity (331 models, snapshot 2026-07-03)

Data, not a bngsim release. No package, nothing on PyPI.

331 SBML files that only BioModels serves: the biomodels_dir jobs of parity_checks/rr_parity/ode_jobs.json. Each is the version pinned by benchmarks/suites/biomodels/biomodels_provenance.json (snapshot 2026-07-03), with BioModels' bytes unmodified.

Why it is here. BioModels sits behind Amazon CloudFront, which refuses every request from a GitHub-hosted runner with its own "Request blocked" page (#807). The nightly parity workflow reads this copy instead, so it can run all 1,323 BioModels models rather than 1,006.

License. BioModels dedicates its whole dataset to the public domain under CC0 1.0: https://www.biomodels.org/termsofuse

Checking it. biomodels-sbml-2026-07-03.tar.gz (11.2 MB), sha256 09e4337984bb03737f6120380d8b4cb40b7d23f3206591bbab2bfcd641a1c407.

parity_checks/rr_parity/fetch_biomodels_mirror.py downloads it, checks that sha256, then checks every file against its record in the provenance manifest. --build rebuilds the identical archive from a corpus fetched with benchmarks/suites/biomodels/fetch.py.

bngsim 0.16.0

pip install bngsim==0.16.0

Full changelog for 0.16.0

---

Added

• Per-reaction firing counts and propensity integrals from an exact SSA
  run (issue #616). Simulator(model, method="ssa", reaction_stats=True)
  records, at every output time, N_r(t) and ∫₀ᵗ |a_r| ds for each
  reaction — the two accumulators the likelihood-ratio (Girsanov) parameter
  gradient of an ensemble is built from — on Result.reaction_firing_counts,
  Result.reaction_propensity_integrals and Result.reaction_labels
  ((n_times, n_reactions); empty unless asked for). The blocks survive
  run_replicates (sequential and parallel; stacked to
  (n_sims, n_times, n_reactions) under squeeze=True), save/load and
  to_dataset. The new bngsim.girsanov module turns them into scores for
  named mass-action rate constants (scores, per parameter or per log₁₀)
  and into the ensemble gradient with its standard error (gradient, the
  sample covariance of observable and score). Off by default, and enabling
  it leaves every trajectory byte-identical: the integral is banked at the
  moment a propensity is rewritten, so the incremental loop pays
  O(affected) per step for it, the same as for the propensities. Exact SSA
  only — PSA, ODE and network-free refuse the flag — and scores refuses a
  model with a functional or Michaelis–Menten rate law, whose score is not a
  function of the two accumulators.

• Model.analytical_jacobian_status: why a model runs on the
  finite-difference Jacobian, in words (issue #506). "complete",
  "pending" (functional rate laws not derived yet), or "declined: <reason>"
  — the rate law and what stopped its derivation, the derivation budget, the
  C++ attach gate, or the environment switch — so a harness can record the
  reason instead of only the boolean analytical_jacobian_complete. Carried
  by clone().

• Result.state / Result.state_names: the full integrator state
  trajectory (issue #508). Result.species is the reported block. An SBML
  event that assigns a parameter or a compartment promotes its target to a
  state entry the species projection hides (GH #71), and a species column in
  an event-resized or rate-ruled compartment, or one an assignment rule
  targets, is remapped to its reported value after the solve (GH #85, #131) —
  so a row of species is not a state Model.rhs / Model.jacobian accept:
  too short, or not the point the integrator was at. state is every entry
  as the integrator held it, (n_times, n_state) in Model.species_names
  order, and state_names names its columns; on a .net model, and on an
  SBML model with nothing promoted or remapped, it is the same array as
  species. Solve results only: a Result loaded from disk or stacked from
  a batch raises ValueError for it.

• The model's RHS, Jacobian, stoichiometry matrix and propensity vector are
  public NumPy evaluators, and a steady-state result carries the spectrum its
  stability verdict was read off (issue #523). Model.rhs(y, t),
  Model.jacobian(y, t, sparse=), Model.propensities(y, t) and
  Model.stoichiometry_matrix(sparse=) evaluate at any state the caller
  passes, under the live parameters, without touching the stored
  concentrations. Each is the evaluator the corresponding engine runs, not a
  re-derivation: rhs is the CVODE callback's interpreted RHS with its
  observable-and-function refresh run at y; jacobian takes the closed form
  when it covers every reaction and the steady-state solver's own difference
  quotient otherwise — the same rule, lifted into bngsim/fd_jacobian.hpp so
  the two cannot diverge — and reports which through .source
  ("analytical" / "finite-difference", the spellings
  solver_jacobian_source already used); a partial analytical Jacobian is
  never returned as if whole. propensities is the SSA loop's propensity pass,
  refresh included, in the SSA volume convention (falling factorial for a
  repeated reactant, times the compartment volume), and says so.
  stoichiometry_matrix is the net-coefficient matrix conservation-law
  detection row-reduces, 0-based, with a zero row for a $-fixed species. The
  dense Jacobian is a bngsim.JacobianMatrix (an ndarray with source);
  sparse=True on either matrix returns a scipy CSC array over the model's
  structural pattern. The C++ names are bound one-to-one on the core
  (compute_derivs, fill_dense_analytical_jacobian, fill_dense_fd_jacobian,
  fill_sparse_analytical_jacobian, compute_propensity,
  compute_propensities, stoichiometry, jacobian_sparsity,
  species_is_fixed).

  SteadyStateResult.eigenvalues is the spectrum of the Jacobian restricted to
  the species the Newton polish solved for — the matrix the issue #78
  certificate already took and discarded after reading one bit — sorted by
  descending real part with a conjugate pair adjacent, +imag first. It is
  filled whenever the certificate computed a spectrum and empty when it did
  not: an integration result, or a root it declined on above 512 unknowns,
  consistent with "undetermined" there. Together these are the provider the
  continuation (#501), multistart (#524) and generator (#504) work assumes.
  New user-guide page: Model evaluators.

Changed

• copy.copy on a Model or a Simulator now refuses instead of handing
  back an alias (issue #643). copy.copy(model) returned a Model sharing
  the original's engine handle, which is textbook shallow-copy behaviour —
  copy.copy is documented to insert references to the objects found in the
  original. What made it a trap here is that all of a Model's state lives
  behind that one handle, so there was nothing left for the copy to own: it was
  an alias with a different id(). A set_param or set_concentration through
  the "copy" landed on the original, which then integrated to a different answer
  with nothing raised. A Simulator shares its Model outright and aliased by
  the same route. Neither of the two things a caller reaches for next corrected
  it: copy.deepcopy raises (it cannot pickle the handle), and clone(), which
  does give an independent model, is invisible from the copy path.

  Both now raise a TypeError naming what the copy would have shared and the
  remedy — clone() for a Model, Simulator(model.clone(), method=...) for a
  Simulator. Breaking, and deliberately so: it removes an operation that
  ran before, which is the point, because it ran, looked right, and quietly made
  the original wrong. Returning clone() from __copy__ was considered and
  rejected — it would trade a loud failure for a silently different one for any
  caller who copies and then mutates.

  copy.deepcopy does not route through __copy__ — it falls to __reduce_ex__,
  so it still reports the issue #636 pickle refusal. That message now ends with
  the in-process half for a Simulator too ("Within this process, build a second
  Simulator over model.clone()"), which Model's has carried since #636 and
  which is the only advice that fits a deepcopy caller, who is not crossing a
  process boundary at all. It is a separate sentence, guarded by "within this
  process", so it cannot be read as a way to make a Simulator picklable.

  Result keeps its shallow copy. Its API is read-only — no setter, no mutating
  method — so no call on the copy can write into the original; a test pins that
  premise, and fails if a setter is ever added. Its arrays are shared, so an
  in-place copy.copy(result).species[0] = ... does reach the original — but so
  does the same write through result itself, since every access returns the
  stored array rather than a fresh one. That is documented on Result rather
  than fixed, being numpy's aliasing and not the copy's.

• Pickling a Result, Model or Simulator now says what to do instead of
  naming an internal C++ type (issue #636). Each holds a handle into the
  compiled extension, so none of them can cross a process boundary — a design
  boundary, not a defect. What a caller got was
  TypeError: cannot pickle 'bngsim._bngsim_core.ResultCore' object, which
  names a type they have never seen, at the moment they are parallelising a fit
  across a process pool and need to know that the way through is to convert
  first. The refusal now names itself, the reason, and the conversions that
  carry what a worker actually needs: as_roadrunner(), to_xarray(), the
  plain .time / .species / .observables arrays, or save(path) here and
  Result.load(path) on the other side; for a Model, rebuilding it there
  from the same source, with clone() for an independent copy in this process.
  bngsim's own parallelism is threads, so nothing in the library ever reached
  this message — a consumer reaches it on the first return value.

  copy.deepcopy gets the same sentence, having previously failed with the same
  internal type name. copy.copy was deliberately left as it was: it falls back
  to __reduce_ex__, so the refusal would otherwise have taken a shallow copy
  with it as a side effect of a change about an error message. What a shallow
  copy of these three should mean was then settled on its own terms — see
  the issue #643 entry above. Every remedy the messages name is checked by a test to
  exist and to survive a round trip, because an error message that recommends a
  method is a claim about the API.

• SIR_v4 is re-sourced from wshlavacek/BNGL-Models after its model fix, so
  the corpus model keeps its analytical Jacobian and analytic sensitivity RHS
  (issue #543). Its four year-selection if() chains fell back to 0 after its
  last modelled year, so its rate law divided by zero past day 1461
  (wshlavacek/BNGL-Models#55 ends them with the 2024 values), and each season
  opened on t>=t_start(), the one instant whe...

bngsim 0.15.1

pip install bngsim==0.15.1

Full changelog for 0.15.1

---

Changed

• Re-vendored RuleMonkey 3.10.0 (3d97f66) → 3.10.1 (87628b3). A patch
  release of two performance fixes. Nothing a model does under method="rm"
  changes: both upstream entries record all 164 corpus, feature-coverage and
  BNG2-oracle models producing byte-identical output, RuleMonkey's include/ is
  byte-identical to 3.10.0, and no Tier-0 refusal, warning or diagnostic is added
  or removed. Only cost changes, and it changes by orders of magnitude.

  A rate law that reads a Species observable no longer full-walks the pool
  after every SSA event. A Species observable is tracked incrementally, but
  the tracked value only settled when its dirty complexes were flushed, and that
  flush fired at sample time — far too late for a rate law reading the value on
  the very next propensity recompute. So RuleMonkey excluded rate-dependent
  Species observables from the tracker outright and walked every complex holding
  a molecule of the pattern's seed type, from scratch, after every event. That
  made such a model quadratic in its own seed population; the flush now also runs
  after each event, for the observables a rate law actually reads.

  seed	before	after
  10 000	88.9 us/event	5.3 us/event
  100 000	559.0 us/event	3.3 us/event
  300 000	did not finish	3.6 us/event

  The charge did not need the observable to move: a rate over a Species
  observable of an inert pool no rule touches paid a flat ~650 us/event for a
  value that could not have changed.

  The same fix then had to reach multi-pattern Species observables, which
  the first gate admitted only one pattern at a time — Species has taken any
  number of patterns all along, so those kept falling through to the full walk
  (3 502 us/event at 10 000 seed, against 14.3 after). Lifting it was a semantics
  change rather than a gate change: the walk applies each pattern's quantity
  relation to that pattern's own complex-wide count, so a complex matched by two
  patterns counts twice, and the tracker's one-tally-per-complex tables had to be
  split per pattern first. An n-pattern observable now holds n times as much;
  the peak-RSS sweep across all 164 models moves by less than its own
  run-to-run spread.

  Unlike the 3.9.0 → 3.10.0 refresh, this one needed no two-repo ordering.
  RuleMonkey's standalone third_party/bngsim_expr copy already matched this
  tree on all four EXPRTK_SYNC_FILES paths, so the guard passed as-is and the
  vendored surface moved by four files: CHANGELOG.md, CMakeLists.txt,
  cpp/rulemonkey/engine.cpp and docs/model_semantics.md.

• MODEL2002070001's INVALID_REFERENCE premise is relaxed: a reference that
  answers NaN is not a comparison basis, whatever bngsim did (#485). The
  0.15.0 sweep flagged the entry stale — natural=DIFF, rr_finite=False — and it
  was half right. The entry claimed both engines were broken; bngsim's half of
  that stopped being true in 0.14.0, when the #353 unit-volume substitution for a
  non-finite compartment size let it integrate the model. RoadRunner's half is
  unchanged. So the natural verdict resurfaced as a DIFF with max_rel_err
inf, 606/707 cells failing — RoadRunner's NaN scored against bngsim as a
  divergence, the one thing it is not.

  The model is well posed and the NaN is not in its dynamics. Both compartments
  declare size="NaN", but every species is hasOnlySubstanceUnits="true", the
  dynamics are 7 rate rules with no reaction, and no <math> block in the file
  names mucosa or lumen — nothing reads a size, so the IVP is well posed in
  amounts and bngsim's substitution warning describes a value nothing reads.
  RoadRunner carries species as concentrations, so x1..x6 — the six species in
  those compartments — are NaN from t = 0, while x7, the one rate-rule
  variable that is a parameter rather than a species, stays finite and agrees
  with bngsim cell for cell. The NaN reaches RoadRunner's state through the
  compartment, not through the dynamics.

  At the sweep tolerance (rtol 1e-9 / atol 1e-12, t = 0..100, 101 points):
  bngsim 707/707 cells finite; RoadRunner 2.9.2 runs without raising and returns
  101/707; AMICI 1.0.1 confirms bngsim's trajectory independently, 0/707 cells
  failing at max_rel_err 0.0. Two engines produce a trajectory here and the
  third produces NaN, which is all the override needs to claim: RoadRunner's
  output is not a usable comparison basis.

  _apply_invalid_reference now keys the premise on the reference alone — RR ran
  and emitted non-finite output — instead of also requiring that bngsim failed.
  What that requirement was really doing was guarding against burying a bngsim
  defect, and a model bngsim can now integrate needs a guard that still works:
  _reference_nonfinite_covers re-runs the same verdict over just the columns
  the reference kept finite, and the override holds on a DIFF only if that
  passes. A real bngsim divergence beside the NaN leaves a failing cell in a
  finite reference column, and the row stays a DIFF. Restricting the columns
  can only tighten the re-run (scale is the peak over the retained block), so
  the claim it licenses is never a loosening.

  Two shapes the relaxation deliberately does not cover. rr_finite is None
  — RoadRunner raised, or a segfaulted child left no per-engine status — stays
  stale on anything but a natural BAD_TEST: "the reference raised" and "we
  never learned" are indistinguishable there, and REFERENCE_FAILED is already
  the right bucket for the first. A natural PASS stays stale too: if the
  compared columns agreed, the reference was usable for the comparison whatever
  it did elsewhere.

  The row lands on BAD_TEST, not PASS. PASS is the #482/#483 treatment for
  a divergence attributable to the reference engine, and it would score
  RoadRunner as having validated a model it answered NaN on. BAD_TEST says
  what happened: the model's own SBML put a NaN in the reference's state, and
  no comparison exists. Both are non-scoring, so neither costs bngsim a model.

  Suite-side only — no bngsim behavior changes. overrides.py rewrites the entry
  (now citing #485, the AMICI confirmation, and the x7 control),
  ode_jobs.json re-bakes the one job record, and test_outcome_classify.py
  pins the relaxed premise and the new guard. This was the last DIFF in the ODE
  sweep: 1237 PASS / 27 REFERENCE_FAILED / 54 BAD_TEST / 5 TIMEOUT over
  1323 jobs, no DIFF and no EXCEPTION row, and no override flagged stale.

• BIOMD0000000627's rr_parity divergence is allow-listed as a KNOWN_ARTIFACT:
  RoadRunner steps over a smooth pre-stimulus ramp, and bngsim is the correct
  engine (#482). Reported as a 0.13.0 regression on the premise that the
  pre-switch trajectory is a steady state, so every sample before t = 200 must
  sit on a flat baseline. It is not, and they must not.

  Jolivet2015 gates its stimulus on an assignment rule,
  is_stimulated = piecewise(0, (time<=200) or (time>=t_0+t_end), 1), but the
  blood-flow multiplier alongside it is ungated: f_CBF_dyn is a logistic
  centred at t_0 + t_1 - 3 = 199, so perfusion (F_in = F_0*f_CBF_dyn) runs
  1.0000004 → 1.4157864 over t = 196..200 — the response leads the stimulus by
  one second, by construction. species_24 (CO2) is washed out along that ramp,
  so the correct trajectory is already falling well before the switch.

  bngsim resolves the ramp because the #305 crossing stop lands the step on
  t = 200, which pulls the step size down on the approach. The ramp itself is
  smooth: it brackets no #72 root, and the model has no floor/ceiling/
  modulo, so _periodic_disc_max_step is None and #274 never enters. Remove
  the #305 stop and bngsim reproduces RoadRunner's answer bit for bit — which is
  exactly what bngsim ≤ 0.12.2 did, and why the fix reads as a regression.

  RoadRunner at the sweep grid steps clean over the ramp and reports the flat
  baseline 2.2084924 through t = 200 (worst species_24 at t = 199.6:
  bngsim 1.9434505 vs RR 2.2084924, reld 0.136). Its answer gets worse as tol
  tightens — at rtol/atol 1e-11/1e-14 and 1e-12/1e-16 it stays flat through
  t = 200.4 as well, then missing part of the stimulus response too. Tol-stable
  wrongness is the signature of a feature never sampled, and a denser output grid
  cannot help, because CV_NORMAL interpolates output points. TOL_OVERRIDES is
  therefore not the remedy.

  Three independent converged references agree with bngsim to 7 digits at
  t = 197.6/198.8/199.6/200.0/200.4
  (2.2081597/2.1457629/1.9434505/1.9173758/1.9130018): RoadRunner with
  maximum_time_step 0.05; RoadRunner with maximum_time_step 0.01, agreeing
  with the first to 7 digits and so converged; and RoadRunner with no step
  bound at all, simply restarted at t = 190 so its own step control has to
  resolve the window.

  Suite-side only — no bngsim behavior changes. parity_checks/rr_parity/ overrides.py gains the entry and ode_jobs.json re-bakes the one job record.

• Three rr_parity override keys that could never match are retired, and a
  test now says so (#482 follow-up). MODEL0913003363:ode (TOL_OVERRIDES),
  MODEL1006230083:ode and MODEL6963432821:ode (both NO_ORACLE_ADJUDICATED)
  name models the ODE manifest does not build: all three are drop,no_model in
  the stage-1 corpus filter, their fetched files zero bytes (sha256 = the empty
  string's), one of 305 corpus entries in that state. They have been stale since
  the initial public release.

  A stale key is silent rot. It reads as a live disposition while covering
  nothing, and it is the shape a typo'd or renamed model_id takes — an override
  ...

bngsim 0.15.0

pip install bngsim==0.15.0

Full changelog for 0.15.0

---

Added

• A model that calls mratio() in a rate law now gets an analytic forward
  sensitivity instead of CVODES' difference quotient (issue #457). The
  differentiation layer had never heard of mratio, so any derivative through
  it came back unevaluated and the whole model fell back to the difference
  quotient. That answer is correct but slower, and the derivative it was falling
  back from is available in closed form.

  mratio(a,b,z) is M(a+1,b+1,z) over M(a,b,z), and Kummer's identity is
  dM/dz = (a/b)*M(a+1,b+1,z). Through the quotient rule, writing R for
  mratio,

  dR/dz = R(a,b,z) * [ (a+1)/(b+1)*R(a+1,b+1,z) - (a/b)*R(a,b,z) ]
  

  so the derivative in the third argument is two more mratio calls. No new
  special function and no new numerics. This is the derivative worth having:
  BNG builds a and b from molecule counts and puts the rate constant in z,
  as z = -1/Keq, so fitting a rate constant goes through exactly this. The
  first two arguments have no comparable closed form and get none. A model
  differentiating through one of them declines and keeps the difference
  quotient, as it did before.

  The second call is the awkward part. Issue #453 gave mratio a region it
  trusts and made it refuse outside it, and (a+1, b+1, z) is not automatically
  inside: 97.6% of the arguments mratio answers for have the shifted call
  answered too, and the whole of the rest has a above -1, where a+1 turns
  positive. Every model BNG generates has a at or below -1. Refusing the rest
  would turn a model that runs today into one that fails, since the compiled
  path has no retry, so the emitted helper falls back to a second exact
  expression for the same derivative. Kummer's equation
  z*M'' + (b-z)*M' - a*M = 0 gives the contiguous relation

  (a+1)/(b+1)*R(a+1,b+1,z)*R(a,b,z) = ( b - (b-z)*R(a,b,z) ) / z
  

  so the whole derivative follows from the one call the value already makes.
  That subtraction cancels, which is why it is the second choice and not the
  first: against a 60 digit reference over 1293 arguments the shifted-call form
  is right to 9e-16 at the median and this one only to 5e-13, and for a small
  |z| it loses everything. A small |z| is where mratio trusts the shifted
  call unconditionally, so the fallback only ever runs at a large |z|, where
  its worst measured error is 1.7e-8.

  The interpreted engine has mratio but nothing for its derivative, so the
  interpreted Jacobian declines it and uses finite differences, which is what it
  did before. Only the generated C has the helper.

• mratio() answers for 40% more of the arguments it used to refuse (issue
  #456). Issue #453 made mratio(a, b, z) refuse the arguments its continued
  fraction cannot be trusted with, which removed a class of silently wrong
  answers. The cost was that it also refused arguments the fraction would have
  got right, since in the uncertain region there is no way to tell which ones
  those are. This narrows that region with a second method that can vouch for
  its own answer.

  For a large |z| the ratio has an asymptotic expansion in which the Gamma
  factors cancel between numerator and denominator. That cancellation is the
  point: the individual Kummer functions overflow a double long before their
  ratio does, which is the whole reason mratio exists.

  z to -infinity:  R(a,b,z) = (b/(-z)) * 2F0(a+1, a-b+1; ; -1/z)
                                       / 2F0(a,   a-b+1; ; -1/z)
  z to +infinity:  R(a,b,z) = (b/a)    * 2F0(b-a, -a;    ;  1/z)
                                       / 2F0(b-a, 1-a;   ;  1/z)
  

  2F0 is a divergent asymptotic series, so each is summed to its smallest
  term, which is where such a series comes closest to the function it stands
  for, and that term estimates how close it got. Certifying its own accuracy is
  what makes the route safe to add at all: it is consulted only for arguments
  that were going to be refused, and only its own estimate lets it answer, so it
  cannot change a value mratio already returns and cannot turn a refusal into
  a wrong number. Where the estimate is not small enough the refusal stands.

  Two things it gets wrong if they are not guarded, both found by breaking it.
  An expansion whose leading coefficient is 1/Gamma(b-a) has no leading term
  at all when b - a is a non-positive integer, so what is left is the
  exponentially small piece the expansion drops; unguarded it returned 5e-06 for
  M(a,a,z) over itself, which is exactly 1. And a series that stopped early
  because a term hit exactly zero has demonstrated no decay, so its smallest
  term says nothing; unguarded it returned 14.33 where the ratio is 0.9998.
  Neither needs a test of its own in the end. lgamma is +inf at exactly the
  first set of arguments, so estimating the dropped branch in logs refuses them
  by arithmetic, and taking the estimate to be the last term actually added
  rather than the zero that ended the series covers the second.

  Measured over 5955 argument triples against a 40 to 60 digit reference, the
  refused fraction of the grid falls from 46.4% to 28.5% and the worst error
  among the newly returned values is 2.7e-11. The estimate itself is honest to
  about a factor of ten, so the 1e-10 it is held to is worth roughly 1e-9 in the
  worst case.

  As with the fraction, src/expression.cpp is the single source of truth and
  the generated C carries a copy held to it. The two agree bit for bit at every
  one of the 5955 triples.

• A dosing schedule that went through libSBML now gets the same gradient as one
  written by hand (issue #465). libSBML does not emit rem(a, b). It expands it
  into a sign test over two remainders,

  if(sign(a) != sign(b), a - b*ceil(a/b), a - b*floor(a/b))
  

  so a model writing rem(time(), P) >= d — the same repeating schedule issue #436
  compensated — arrived with the schedule hidden behind an if() inside the
  condition, and was refused for having two step functions where the recogniser
  wanted one. Five corpus models write it. The gradient a model got therefore
  depended on which tool had written its SBML, which is not a distinction anyone
  fitting it can act on.

  The schedule is unchanged; what is new is getting to it. An if() inside a
  condition arrives as a Piecewise, and it is collapsed to the branch the run
  takes before the step functions are counted. Which branch that is cannot be read
  off the text: for a positive P the floor remainder runs over [0, P) and the
  ceil one over (-P, 0], so the two differ by P everywhere except at the
  multiples of it. So the guard is resolved against the run's own clock domain and
  the model's parameter point — sign(clock) is 1 for the whole of a run, what is
  left must be clock-free, and it must evaluate at the parameter point. A negative
  period selects the other branch and is read just as well.

  The route not taken is the interesting part. The obvious fix is to propose
  both branches and keep whichever _schedule_matches_residual accepts, since that
  check exists to catch a schedule the model does not follow. It is unsafe, and the
  failure is silent. The ceil branch of a positive-period schedule reads as a
  negative period, so the check's probe points land at negative clock values —
  outside any run window — and there the guard really does select that branch and
  the residual really does hold one sign for the whole period. It is accepted as a
  schedule that never turns over, and the model is admitted with every one of its
  real crossings uncompensated: exactly the quiet failure that check was built to
  prevent. That is pinned as a test rather than left as a comment.

  A guard that still reads the clock after the substitution is declined, because it
  changes which branch is live partway through the run and no single period, offset
  and duty describes the window — MODEL1006230027's shape, which stays refused
  along with MODEL1708310001's remainder of a remainder and Sturrock2015's
  sinusoid.

  One thing had to move besides the recogniser. _iter_condition_atoms descends
  into an if()'s condition as readily as into its branches — rightly, since a
  nested if() in a branch is a real discontinuity — so the guard
  sign(time()) != sign(P) also arrived as an atom in its own right, was neither a
  clock threshold nor a comparison over state, and refused the model even with the
  schedule around it fully compensated. It is not a crossing: sign(clock) is 1
  for the whole of a run, so it picks its branch before the first step and holds
  it. That is the same ground condition_cannot_cross already applies to a
  comparison over run-constants, reached by reading a value rather than a name —
  which is why it could not go in that function, documented as structural and
  never numeric, and is asked from clock_crossing_compensated instead, which is
  value-dependent already and is what switch_gate_cache_digest carries. It is
  narrow on purpose: only an atom that reads the clock and stops reading it once
  sign(clock) is resolved is claimed, so time() >= sigma is left to the
  recognisers exactly as before.

  The gradient is checked against a central difference of two trajectories on every
  column, through the same h-ladder oracle issue #436 uses, and the same model
  written the two ways places identical edges with identical ∂t*/∂p. Forcing the
  collapse to the wrong branch fails seven of the new tests, the P and d
  columns by an order-one fraction.

• A rate law that switches on a repeating schedule now has a forward-sensitivity
  gradient (issue #436). if(time() - 24*floor(time()/24) >= 7, on, off) is "on
  for the last 17 hours of every 24 hour day", and with a start ...

bngsim 0.14.0

pip install bngsim==0.14.0

Full changelog for 0.14.0

---

Added

• A rate-law switch condition that is a single power of the clock is now
  compensated, so its forward-sensitivity run is no longer refused (issue #418).
  if(time()*time() >= thresh, ...) — the shape issue #414 refused, since neither
  issue #48's affine solver nor issue #150's state root brackets it — has a
  crossing in closed form: c·clock^n is strictly monotonic on clock ≥ 0, so it
  crosses exactly once, at clock = (thresh/c)^(1/n). _clock_monomial_threshold
  solves that (time = sqrt(thresh) for the quadratic), and because the result is
  the clock value at the crossing — the same (clock, threshold_expr) contract
  _clock_affine_threshold (issue #355) already returns — the whole issue #48
  machinery jumps it unchanged: ∂t*/∂thresh = 1/(2·sqrt(thresh)), the crossing
  registered as a stop time, the in-branch ∂f/∂p a clean Piecewise zero. The
  detector and the issue #68 gate share the one recognizer, so a model that was
  declined for this shape is now admitted on both. This is the first step of the
  machinery half of issue #414; a clock threshold that is not a bare power
  ((time-5)^2, two crossings; time^2 + time, mixed) has no single crossing to
  name from the text alone and stays refused, and the difference-quotient
  reproducer that #415 exercises under the MIR JIT now runs the analytic path with
  its self-multiply in the sensitivity RHS as well as the state RHS.

• A codegen artifact says which codegen key built it, so a cache sweep can tell
  live from orphaned (issue #363). Artifacts are now named
  rhs_<key>_<hash><suffix> — e.g. rhs_28+317a5b34d5dc9959_9e1f… — where <key>
  is _CODEGEN_CACHE_KEY, the _CODEGEN_VERSION constant plus a digest of the
  emitters' own source (#51). The key is still mixed into <hash>, so nothing
  about invalidation changes; carrying it beside the hash is what makes the dead
  corpus countable.

  That matters because #51's key is deliberately conservative: any edit to
  _codegen.py / _jacobian.py / _saturable_jacobian.py / _switch_sensitivity.py
  — a comment included — orphans every artifact on the machine at once. So on a
  machine that tracks bngsim development, "everything from before my last emitter
  edit is garbage" is the common case, and until now it was unexpressible: two
  opaque hex strings, with time as the only signal a filename carried.

  • bngsim-cache info reports live and orphaned as two numbers, plus a
    per-key table — the audit of a shared or pre-warmed artifact directory, which
    says which bngsim's artifacts are in it and how much each holds.
  • bngsim-cache prune --orphaned sweeps exactly the artifacts this install can
    never load again. Much better targeted than --older-than, which keeps orphans
    that happen to be recent and evicts live artifacts that are merely idle. The
    --min-age floor still applies: a fresh orphan is a sibling process's compile
    under its own key, about to be dlopened.
  • --keep-key KEY (repeatable) spares another install's artifacts, because a venv
    per project is ordinary and each has its own key — "orphaned" from one venv's
    point of view is live from another's, and sweeping on that would make a shared
    cache thrash. prune_codegen_cache(orphaned=True, keep_keys=[…]) is the API
    form; passing keep_keys without orphaned raises rather than reading as a
    protection the age and size bounds do not honor.
  • CacheEntry.codegen_key, CacheInfo.live / .orphaned / .by_key, and
    bngsim.cache.artifact_key() expose the same thing in process.

  This is a one-time full invalidation. Every artifact already on disk is under
  the old scheme and unreachable under the new one, so every model recompiles once.
  info counts those pre-#363 names as orphaned (nothing will ever look one up
  again) and lists them under -; prune --orphaned --keep-key - spares them, for
  a directory shared with an install too old to write a key. No _CODEGEN_VERSION
  bump: this edits _codegen.py, so #51's source digest invalidates every cache on
  every machine anyway — which is the argument for landing the rename together with
  it, since the cost is paid once either way.

• The codegen artifact cache is inspectable and prunable: bngsim-cache
  (issue #205). ~/.cache/bngsim/codegen grew without bound and there was no
  supported way to look at it — 2.0 GB across 14,377 entries after six weeks of
  ordinary work on one developer machine, with rm -rf on a path you had to know
  by heart as the only remedy. New console script bngsim-cache (also
  python -m bngsim.cache, which needs no reinstall) with four verbs:

  • info — path, entry count, total size, build/last-used dates, and a
    breakdown by artifact kind (model RHS, SSA propensity, the src_ fallback
    key of #174, plus the leaked partials below).
  • clean — remove only the debris of interrupted compiles: bngsim_shard_*
    scratch directories, the stray rhs_<key>_<hash>.<pid>_<n>.c beside them, and
    temp/sidecar libraries. Nothing else cleans these up, and no compiled
    artifact is touched, so no cache hit is lost.
  • prune --older-than 30d / --max-size 2G — evict least-recently-used
    artifacts to fit an age and/or size bound. The partial sweep runs first, so
    the size cap is a bound on the whole directory.
  • clear — everything, behind a confirmation prompt (--yes for scripts;
    required, not assumed, when stdin is not a terminal).

  --dry-run on every mutating verb, --json on info, and -C/--cache-dir to
  point any of them at a directory other than the configured cache — e.g. to
  audit a pre-warmed artifact directory from a login node.
  bngsim.codegen_cache_info(), clean_codegen_cache(),
  prune_codegen_cache() and clear_codegen_cache() are the same four in
  process, for a notebook or a fitting harness.

  Nothing prunes automatically, deliberately. No size cap at compile_rhs
  time, no sweep on import or on Simulator construction. A library that deletes
  files as a side effect of being used is a surprise, and the failure mode —
  evicting the artifact another process is about to dlopen — is the exact class
  of bug the cache exists to avoid.

  Two safety properties hold across every verb. Only files bngsim wrote are
  ever removed: BNGSIM_CODEGEN_CACHE_DIR is a user-supplied path that people
  point at shared scratch, so anything unrecognized is classified foreign,
  reported, and left alone — by clear as much as by clean. And nothing
  recently touched is removed: a compile in flight writes its scratch files
  into this very directory, so every verb holds off on entries used or written
  within --min-age (default 1h, over the 600 s default BNGSIM_CODEGEN_TIMEOUT) and
  on POSIX additionally holds a partial whose compile is still running.

  LRU order is max(atime, mtime), which was measured rather than assumed: a
  plain read() on macOS APFS leaves atime untouched while dlopen — the only
  way this cache is ever used — advances it, and a noatime mount never moves it
  at all. Taking the newer of the two gives true recency wherever the filesystem
  records it and degrades to build order where it does not, and info reports
  which one you are getting.

  Not addressed here, and filed rather than guessed at: an entry carried no
  record of the codegen key that built it, so nothing could report how much of the
  cache was orphaned (#363, above, which put the key in the filename); and MSVC
  leaked a .lib/.exp pair per successful Windows compile (#362, whose fix rode
  along with an invalidation this release was paying anyway — clean collects any
  debris left from before it).

• BNGL models load, behind a bngl extra (issue #162).
  Model.from_bngl("m.bngl") and Model.load("m.bngl") now work;
  _LOAD_DISPATCH gained .bngl and the "expand it yourself first" refusal is
  gone. bngsim still has no BNGL parser and does not want one — BNGL describes
  rules, and turning them into a network is BNG2.pl generate_network's job —
  so the loader shells out and reads the emitted .net through
  Model.from_net. pip install 'bngsim[bngl]' supplies BNG2.pl via
  PyBioNetGen; $BNG2_PL / $BNGPATH / a BNG2.pl on PATH all outrank it,
  and perl has to be there too (so BNGL loading reports itself unavailable on
  stock Windows rather than failing at load time).

  The file's experiment is not executed — and only the experiment is
  stripped. A .bngl in the wild ends in simulate({...}) or
  parameter_scan({...}), and BNG2.pl runs whatever it is handed, so loading a
  model would have meant running the author's whole experiment to obtain a
  network generate_network alone produces in seconds. But dropping every
  action is the opposite mistake, and the first cut here made it: an action-scope
  setOption("NumberPerQuantityUnit",6.0221e23) above the model block
  configures the generation, and without it BNG2.pl emitted the same topology
  with every bimolecular rate constant off by that factor — 1e12 where the
  reference network says 1.66e-12, caught by diffing
  benchmarks/models/bngl/ode/catalysis.bngl against a direct BNG2.pl run.
  What is dropped is now exactly the verbs that would run the experiment, write
  an artifact (a stray writeNET would land on the network being read), or stop
  BNG2.pl early — plus everything after the source's own generate_network,
  which is where the protocol begins. That call's max_iter / max_agg /
  max_stoich are forwarded: they are what make an unbounded rule set finite, so
  regenerating with BNG2.pl's defaults would silently give a different model (or
  never return). from_bngl(..., protocol=True) r...

bngsim 0.13.0

pip install bngsim==0.13.0

Full changelog for 0.13.0

---

Added

• bngsim.HAS_LAPACK_DENSE, and a Linux CI leg that actually runs the BLAS
  dense solver (issue #269). The optional dgetrf dense factor (GH #84) had
  no supported capability probe and, after #265, no Linux coverage at all.

  The probe first, because the CI guard needed it: the flag existed only as
  bngsim._bngsim_core.HAS_LAPACK_DENSE, so the one test file that gated on it
  reached into the private extension module with a getattr default. It is now
  bngsim.HAS_LAPACK_DENSE, in __all__, mirrored by
  capabilities()["features"]["lapack_dense"] with a missing entry that names
  a concrete rebuild path. False costs speed on large dense Jacobians and
  nothing else — the built-in LU is the default on every platform and the BLAS
  path is opt-in via BNGSIM_LAPACK_DENSE=1 — so the flag answers "will that
  variable do anything here?". Of the published wheels only the macOS ones say
  yes; the manylinux and Windows legs resolve no LAPACK. The population this
  backend actually reaches is source installs on a LAPACK-equipped Linux box.

  That population had never run a CI job. python-tests.yml built the backend
  on Linux only as a side effect of libsuitesparse-dev, which #265 removed on
  purpose (it is what makes that leg the bare host an sdist install lands on),
  and native-tests.yml builds on a stock ubuntu-latest where
  find_package(LAPACK QUIET) finds nothing. Two new legs, neither touching the
  bare-host proof:

  • native-tests.yml gains a blas: liblapack matrix entry. It configures
    -DBNGSIM_ENABLE_KLU=OFF, so it never reaches the SuiteSparse autobuild and
    cannot weaken #178 by construction. The stock entry stays — that is the
    configuration the manylinux and Windows wheels ship.
  • python-tests.yml gains a third leg, ubuntu-latest + liblapack-dev and
    nothing else. Still no SuiteSparse, so it still exercises
    BNGSIM_KLU_AUTOBUILD. The bare leg is untouched.

  Both assert the backend they expect, in both directions: a leg that quietly
  loses its BLAS fails, and so does a bare leg that quietly gains one.

• An absolute tolerance that follows the trajectory: atol="tracking" /
  bngsim.TrackingAtol (issue #213). #196 gave the state axis a per-species
  atol and removed the cross-species compromise. It left the within-species,
  over-time half untouched: whatever number species i gets, it keeps for the
  whole run, so a species that starts at order one and decays to something tiny
  outgrows its own tolerance partway through and stops being error-controlled
  from there on. CVODE's construct for that is CVodeWFtolerances, an
  error-weight function evaluated at the state actually being integrated, and
  this is that third mode. The rule:

  atol_i(y) = clamp(rtol*|y_i|, ceiling_i * 10**-decades, ceiling_i)
  

  ceiling is the #196 vector ("auto" by default), so the mode is a strict
  extension of it: at decades=0 it reduces to that vector exactly, and the
  clamp's upper end means tracking is never looser than the vector it was
  built from — only tighter, by at most decades.

  New tests/data/deep_decay.net is the reproducer, and it is built so the
  #196 vector is provably a no-op on it: both live species start at exactly
  1.0, so atol="auto" derives 1e-8 for every species — the number the scalar
  it replaces would have used. One then decays sixteen decades. Worst relative
  error in that species against the analytical exp(-t):

  mode	rel. err.	steps
  scalar atol=1e-8	1.8e+04	162
  atol="auto" (#196 vector)	1.8e+04	162
  TrackingAtol(decades=3)	3.6e+01	275
  TrackingAtol(decades=6)	9.5e-02	366
  TrackingAtol() (12)	2.6e-06	605
  TrackingAtol(decades=20)	9.4e-06	477

  The default depth of 12 is where the accuracy stops improving; past it the
  limit is roundoff rather than the tolerance, which is why it gets slightly
  worse at 16 and 20. The cost is real and is not hidden: roughly 4x the
  steps on this model.

  Opt-in and orthogonal to #196 — a scalar atol stays on CVodeSStolerances
  and a vector on CVodeSVtolerances, unchanged. Accepted by Simulator.run,
  run_batch, parameter_scan, compute_all_sensitivities, steady_state,
  steady_state_batch and set_tolerances, and frozen once per batch/scan the
  same way "auto" is. hasattr(bngsim, "TrackingAtol") is the capability
  probe.

  The sensitivity axis is re-decided rather than left to inherit: atolS
  keeps reading the ceiling, so turning tracking on leaves the sensitivity
  tolerances exactly where the same vector would have put them. Re-deriving
  atolS from the live state would make that base tighten mid-run, which is
  the step-controller hazard the issue #183 high-water mark exists to avoid.
  Measured on the same fixture, that costs nothing: the sensitivity column's
  worst relative error against -t exp(-t) goes from 2.4e+02 to 1.8e-06 under
  tracking anyway, because the state axis is what drives the step size.

  What it costs on real models, swept over the first 400 rr_parity SBML
  models (391 of which integrate at the default tolerance over t = 0..100).
  37 of those 391 have a species that falls six or more decades below its own
  initial value, so the mechanism is not exotic:

  arm	integrates	lost	gained	steps vs default
  default (control)	391	—	—	1.00
  decades=3	393	0	2	1.18
  decades=6	392	1	2	1.24
  decades=12	387	6	2	1.58

  Six models that integrate at the default tolerance do not at the default
  depth — a species whose value is a difference of large fluxes cannot be
  resolved below its own roundoff, and the step size collapses instead. Those
  failures are loud, not silent, and a solver failure under tracking now names
  the depth and suggests lowering it, because CVODE's own report ("made no
  progress", flag=-3) mentions nothing about the tolerance mode that caused
  it. Two models go the other way and integrate only under tracking.

  The control arm is the important row: the same 391 models, run with the
  default scalar atol before and after this change, produce 391/391
  byte-identical trajectories with identical step counts. Nothing moves
  unless it is asked to.

  On Bruno_JExpBot2016, which #213 cites as the measured case, the
  premise does not reproduce. The issue reports the objective moving 1.3e-05
  between atol=1e-8 and "a converged tail near 1e-16". There is no tail:
  reconstructing the shipped job's six conditions and sweeping both axes, the
  worst deviation of any observable from a rtol=1e-12, atol=1e-16 reference
  is 7.3e-08 at atol=1e-8 and 2.4e-08 at atol=1e-16 — both at the rtol=1e-8
  truncation floor, which is where the sequence stops. Drop rtol to 1e-10 and
  the whole atol axis below 1e-10 collapses onto 5e-10. The reason is
  structural: every live species in that model stays between about 0.2 and 6.4
  for the whole run, so none of them ever goes under its own atol. Tracking
  moves it by 7.9e-08, i.e. into the same band — which is the right answer for a
  model that has nothing for it to fix. The mechanism is real; that model is not
  an instance of it, and tests/data/deep_decay.net is.

• The per-species atol derivation a caller needs is now public (issue
  #212). #196 shipped the capability and exported the wrong half of it.
  Simulator.auto_atol() — which derives from the model's live state — was
  public; the stateless derive_atol(state, rtol, floor=...), which derives
  from a state you hand it, was in the private bngsim._atol, along with the
  AUTO token. That is backwards for the primary consumer. A parameter fit has
  to derive from the model's nominal state and hold the vector constant: a
  tolerance that moved with a fitted initial condition would put a step in the
  objective everywhere the derivation crossed a rounding boundary, and it fails
  invisibly — the objective still looks correct, the finite-difference gradient
  check still passes, and only the search behaves oddly. So the function a
  fitting frontend needs was the one it could not reach, and the one it could
  reach was the one that would break it. bngsim.derive_atol,
  bngsim.normalize_atol_vector and bngsim.AUTO are now exported from the
  package namespace and listed in __all__.

  hasattr(bngsim, "AUTO") is consequently a working capability probe for the
  whole feature. It previously returned False on a build that has the
  capability, which is the worst failure mode a probe has — it silently routes a
  capable install down the scalar fallback and the result still looks right. The
  version string is not a substitute: the checkout that first carried #196 still
  declared 0.12.2, the version of the wheel 25 commits behind it.

  normalize_atol_vector went public rather than staying private on a narrower
  argument: every atol= entry point already runs a vector through it, so
  calling it yourself changes only when the error arrives — and for a caller
  that assembles its own vector (per-species clamping, a table, a nominal
  state), taking the length-and-position check once at setup beats taking it at
  the first run() of a fit. is_scalar_atol and the AtolLike alias stay
  internal.

  No behavior changed and nothing was renamed; these are re-exports of the same
  objects. The docs gained the distinction that motivated all of it: which state
  each derivation reads, and the recipe for holding one constant across a fit
  (docs/user-guide/solvers.md, "Which st...

bngsim 0.12.2

pip install bngsim==0.12.2

Full changelog for 0.12.2

---

Added

• Model.effective_ic_sensitivity() — the ∂x(0)/∂θ the parameter axis
  already carries (issue #155). output_sensitivities(axis="parameter") is a
  total derivative: it carries the right-hand-side path and the
  initial-condition seeding. axis="ic" is the companion basis ∂y(t)/∂x_k(0)
  with the initial value held independent. The two are therefore not
  orthogonal —

  d_param[θ] = (right-hand-side path) + Σ_k (∂x_k(0)/∂θ)·d_ic[x_k]
  

  — so a consumer that routes a fitted parameter to every native column it
  reaches and sums them double-counts any seeded initial condition. That was
  always the contract; nothing about the numbers changes here. What changes is
  that it is now stated, and answerable in code rather than by measurement.

  The reader is paired with the existing writer declare_ic_sensitivity, and
  answers from model structure alone — parameter graph, live initial
  conditions, declarations — with no integration, so a fitting frontend builds
  its gradient routing once at setup instead of burning a throwaway run or
  re-deriving per evaluation. It reports the effective matrix: after the issue
  #113 retirement of a species moved off its declared initial condition, and
  after the issue #111 declaration overlay. Keys are the ids
  sensitivity_params= accepts — a compound <initialAssignment> lowered by
  issue #147 reports the original symbols, never the synthetic
  _ic_<species> carrier.

  Present-and-zero is not absent. A present entry valued 0.0 means
  "seeded, coefficient zero at this state" — ∂(a*R0)/∂R0 is a, which
  vanishes at a = 0 without the seeding path ceasing to exist — while an
  absent entry means there is no seeding path at all and the caller's own
  initial-condition term is the missing piece. Only the second is a signal to
  add an ic-axis term. Making that distinction visible required separating
  structural reachability from numeric value in the seed derivation: the
  shared derived-parameter DAG walk drops numeric zeros (it is also the
  switch-time scan's "not a parameter threshold" signal), so
  compute_ic_param_sens_seed now recovers presence from the token closure and
  emits a zero-valued row. The solver still receives no zero rows — they seed
  nothing — so no gradient changes.

  Result.ic_sensitivity_seed carries the same matrix as the per-run record,
  for the cases the model-level reader structurally cannot answer: a batch or
  scan over a nonlinear derived initial condition (Rtot = R0*scale) has
  point-dependent coefficients, so each row stamps its own, squeeze keeps the
  matrix only if every row agreed, and compute_all_sensitivities takes the
  union over parameter chunks. It reports None — not recorded, distinct from
  {} — for a carry_sensitivities=True phase, where the engine seeds from the
  prior phase's dx/dθ and discards these rows entirely (#210/#81); reporting
  the computed-but-unused rows there would have been the same silent-wrong-answer
  class this whole property exists to end.

  Both readers and the solver seeding share one derivation
  (Model._ic_sensitivity_triples), so a reported matrix cannot drift from the
  seed it describes.

  Also: a capabilities()["features"]["effective_ic_sensitivity"] flag to gate
  on instead of a version string (a build without the reader cannot say what the
  seed carries, and every answer a consumer could guess is silently wrong, so
  refusing a gradient fit is the honest behaviour); the contract written into the
  output_sensitivities docstring, docs/reference/api.md and the PyBNF
  integration guide; and a note that SteadyStateResult raises no such question
  (its parameter axis is the implicit-function derivative J·∂x*/∂p = −∂f/∂p,
  with no seeding term, and axis="ic" raises).

Fixed

• Forward sensitivity refused two state-switch crossings at one instant, but
  on the corpus they are always one crossing written twice (issue #153). Issue
  #150 registers one CVODE root per switch condition and deduplicates them by
  the residual's text, which merges X<1 with X<=1 and nothing else. Two
  spellings of ONE crossing therefore arrive as two roots that fire together,
  and the batch was refused outright — each jump reads f on the two branches
  of its own condition, and one step across a shared crossing cannot separate
  them. The reasoning is right for genuinely independent switches; it is just
  not what reaches it. Both models that hit the refusal write one crossing
  twice: sp_fourier_synthesizer roots on ds1 and on 3·ds1 − 12·s1²·ds1
  (five residuals in all, every one a multiple of Cos1 − amp_offset, all
  crossing at t = π/2), and ml_hopfield on dS1/dt and dS3/dt, which are
  identically equal along the trajectory because its own weight matrix leaves
  S1 ≡ S3 invariant. One is visible in the text and one is not, which is why
  the decision belongs at the crossing rather than in the dedup.

  The batch is now carried into the jump together and the two halves of the
  saltation term are checked separately, at the crossing:

  • f⁻ − f⁺ has to carry every branch change, which the one flow probe
    does exactly when it crosses every residual in the batch — so the δt
    ladder now grows until they all flip together, and a batch that never
    does is several conditions meeting by coincidence rather than one crossing.
  • dt*/dθ has to be one vector, and flipping together does not
    establish that. A common factor does: h cancels out of the
    implicit-function ratio, scaling numerator and denominator alike where the
    residual vanishes. An equality that holds only along the trajectory does
    not — ml_hopfield's two residuals have non-parallel gradients (cosine
    −0.30 at the crossing) and their dt*/dθ are permutations of each other
    under the model's W12 ↔ W23 symmetry, so a perturbation that breaks it
    splits the crossing in two. The vector is therefore formed from each
    residual in turn and compared, and a batch that disagrees is refused with
    the numbers.

  The second is the criterion, and it is deliberately weaker than "one
  surface": what the jump needs is a single t*(θ) to shift the flow along,
  which two independent crossings also have when the requested columns move
  them together. A fixture of two unrelated species decaying through a shared
  threshold merges exactly (its ∂W/∂c is the closed-form −8 to eight digits)
  for the threshold column and is refused for the column that splits the pair.

  Order matters and is the cheap way round: the branch gap is measured before
  any dt*/dθ, and both corpus models are the BNGL signed-rate idiom
  (if(r>0, r, 0) against if(r<0, −r, 0), continuous where r = 0), so they
  return with no jump at all — and never reach the second check. Their columns
  come from the in-branch sensitivity RHS and match a finite difference of the
  model's own trajectory.

  Corpus A/B on the netcond arm — the 80 of 585 .net models whose text
  carries an if(), under forward sensitivities, against the issue #150 binary:
  76 identical, 2 newly answered, 1 wall-clock capped, 1 identical refusal.
  Nothing moved: a single-switch crossing takes the same path it did, and the
  two newly answered rows are exactly ml_hopfield and sp_fourier_synthesizer.
  The remaining refusal is ml_q_learning, whose residual is not
  differentiable across its own surface (issue #154).

• Forward sensitivity through a state-dependent switch in a rate law missed
  the saltation jump at its crossing (issue #150). A condition that reads the
  state — piecewise(0, Virus < 1, Virus*rho_V) — flips a branch of f at a
  crossing whose time t*(θ) moves with every parameter through the
  trajectory. The in-branch derivative was never the problem: sympy.diff of
  the Piecewise is right on both sides and carries no boundary delta. What is
  discontinuous is ∂x/∂θ itself, by the saltation term

  s(t*⁺) = s(t*⁻) + (f⁻ − f⁺)·dt*/dθ
  

  and neither way bngsim could produce a sensitivity RHS carried it — the
  analytic one differentiates each branch where it is live, CVODES' internal
  difference quotient integrates the variational equation straight across. GH
  #68 declined the analytic path here and issue #146 corrected the warning to
  say the fallback was wrong too; this is the fix.

  Three pieces, all of them the rate-law twin of something issue #144 already
  built for a state-dependent event trigger: split the condition into a
  residual lhs − rhs (NetworkModel::state_switch, cached per condition and
  re-derived by clone() — an expression id means something else in another
  evaluator); register that residual as a CVODE root, so the crossing is
  located rather than chased; and at the root, read f on both branches at
  x(t*), form dt*/dθ by the implicit function theorem, and jump. The
  ∂t*/∂θ solve is now one function (Impl::residual_dtstar) reached from both
  the event trigger and the rate law, transversality floor included.

  The branch selection is the one genuinely new mechanic. A clock switch picks
  its branch by nudging the clock across its threshold; a state switch has no
  such variable, so the whole state is nudged along the flow, x^± = x(t*) ± δt·f, which for a unit-rate counter reduces to the clock nudge exactly. δt
  is measured, not assumed: CVODE locates a root only to ~100·ε·(|t|+|h|), so
  the size starts just past that and grows until the residual actually changes
  sign across the pair. Failing to flip inside the ladder is refused rather than
  answered with a jump of zero — "we could not tell the branches apart" and "the
  branches agree" are otherwise indistinguishable. The ver...