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 where the pulse's derivative with
  respect to the onset is infinite for its a_2023 = 1.5
  (wshlavacek/BNGL-Models#56 opens them on t>t_start(), the model-side stand-in
  issue #545 describes). vendor_corpus.py takes the model from a dedicated
  bngl_fixes pin at 9ab070d, keyed in a new UPSTREAM_FIXES table the way
  CURATED_SIX is, so every other bngl_models model stays at 81c90d8: moving
  that pin would bring 181 upstream commits with it. manifest.json,
  manifest.csv and jobs.json change only in SIR_v4's record, and
  parity_checks/tests/test_upstream_fixes.py locks the re-source. On the
  regenerated network the analytical Jacobian is complete; a 19-parameter
  forward sensitivity over [0, 1460] runs on the analytic RHS at rtol / atol
  1e-8 / 1e-8, 1e-10 / 1e-12 and 1e-12 / 1e-14, with its d_2023 column within
  1.8e-5 of a central difference of the trajectory; and the bng_parity job still
  passes against BNG2.pl, at max_rel_err 0. The golden record for SIR_v4 is
  stale until the next golden regeneration, as the CURATED_SIX records were.
  gen_networks.py skips a model whose .net is already cached, so a checkout
  with an ode_fullnet cache regenerates this one by removing
  nets/original__bngl_models__my_models__ode__SIR_v4.bngl.net and running
  gen_networks.py --models SIR_v4.

• The TotalRate symmetry carry is gone: RuleWorld/nfsim#92 merged and the
  vendored tree now takes the fix from upstream (issue #428). The carry queue
  is 13 topics, down from 14, and third_party/nfsim is rebuilt on upstream
  master at 1a54f387, the merge commit of #92 — deliberately not the
  master tip, since the 160 commits after it are upstream's performance
  campaign and a separate refresh (issue #519). The carry that left is the one bngsim
  contributed upstream after #89 brought the symmetry factor to every rate law
  (issue #426): a TotalRate rule states its whole propensity and has no match
  count to correct, so the factor is skipped there. Upstream's commit is the
  same guard; only the comment above it was reworded in review, which is the
  one line-level change in the vendored tree. git am of the remaining 13
  patches replays onto the new base without a conflict, and
  python/tests/test_nfsim_symmetry_factor.py with its fixture pair
  tests/data/nfsim/symmetry_factor_total_rate.{bngl,xml} passes unchanged
  against the carry-free tree, which is what shows the upstream version behaves
  the same.

Fixed

• A pickle state bngsim wrote but this class cannot read is refused by name,
  instead of being handed to numpy as if it were numpy's own (issue #646).
  _ndarray_pickle.unwrap_state was asked two questions at once — did bngsim
  write this state, and did this class write it — and answered both with one
  equality test against the caller's tag. A tag it did not recognise therefore
  came back as "not ours at all", which both callers read as a stream older
  than the tagging and passed down to numpy, where it could only surface as
  TypeError: __setstate__() argument 1, item 0 must be tuple, not str — for a
  stream bngsim itself wrote and could have named.

  A state is now recognised as bngsim's by its bngsim. tag namespace, and a
  tagged state whose tag this class does not read is refused naming both tags.
  Nothing reaches it through ordinary pickling, since pickle pairs a state with
  the class that reduced it; it is the diagnostic path the version check already
  builds, which covered a stream from a newer build and not one whose tag moved.
  The pre-tag passthrough is unchanged, and the namespace check is what keeps a
  future 4-tuple from numpy out of the new branch.

  Two more from @Montana's review of #641, in the same function. The version
  refusal now says which direction the mismatch runs — it told every caller the
  stream "was written by a newer bngsim; upgrade to read it", which is the wrong
  way to send someone holding a version this build has dropped. And
  unwrap_state returns an explicit tagged flag rather than signalling a
  pre-tag stream with a None payload: neither class writes a None payload
  today, so that was a trap set for a third class joining the wrapper, whose
  legitimate None would have been read as an untagged stream and silently
  given the unlabeled default.

• JacobianMatrix.source survives a pickle round trip, so a Jacobian
  computed in a worker process or read back from a disk cache still says how
  it was built (issue #635). numpy's reduction carries the array and nothing
  a subclass added, and on unpickling __array_finalize__ is called with no
  parent, where it returned before setting anything — so source was not
  merely wrong after a round trip but missing, and reading it raised
  AttributeError. multiprocessing, concurrent.futures and
  joblib.Memory all move or cache results through pickle, and source is
  what tells a consumer whether to trust the matrix as closed-form or treat it
  as a difference quotient.

  The fix is the one issue #629 applied to NamedArray: __reduce__ wraps
  numpy's own state with source under a tag and version, and __setstate__
  unwraps it. __array_finalize__ now sets source = "" when there is no
  parent, so an array built from nothing has the attribute. A stream written
  by an older build still loads, with source == ""; a state from a future
  format version is refused by name. Pinned at every pickle protocol, for both
  labels, and on views, transposes and Fortran-order arrays.

  The wire format itself is pinned too, with the tag and version spelled out
  rather than imported from the module they describe: a test that read them
  from _evaluators would follow a rename or a version bump instead of
  catching one. That was not hypothetical — before the pin, renaming
  _PICKLE_TAG left every test in test_jacobian_matrix_pickle.py and
  test_named_array_colnames.py green while a stream pickled beforehand
  stopped loading, and not cleanly: an unrecognized tag falls through to the
  legacy-stream path and surfaces as numpy's __setstate__() argument 1, item 0 must be tuple, not str. NamedArray has had that guard since #629, which
  is why the same rename on it was caught at once; JacobianMatrix now has its
  own copy, and the version-refusal test stops importing the constants as well.
  Reported by Michael Mendy (@Montana).

• A NamedArray keeps its column names through a pickle round trip, so an
  as_roadrunner table collected from a worker process arrives labeled
  (issue #629). colnames is a subclass attribute, and numpy's reduction
  carries the array and nothing a subclass added, so every round trip through
  pickle — which is how multiprocessing, joblib and the rest move an object
  between processes — dropped it. A fit that farms simulations out to workers
  got its tables back nameless and every arr["[X]"] in the parent raised,
  while the same code against libroadrunner worked: RoadRunner implements the
  pickle hooks, so this was a divergence from the compatibility target on the
  workflow as_roadrunner exists for. (Before issue #561 it was worse: the
  unpickled array had no colnames attribute at all and reading it raised
  AttributeError.)

  __reduce__ now wraps numpy's own state with the names, and __setstate__
  unwraps it, so the array's buffer handling stays numpy's. The names ride
  along at every pickle protocol, each of which is pinned by a test: numpy
  picks the reduction per protocol, and the no-state _frombuffer form it uses
  for a base ndarray at protocol 5 could not carry a name if a subclass ever
  took that path. A stream written by an older build still loads, and arrives
  unlabeled exactly as that build would have left it; a state from a future
  format version is refused by name rather than misread.

  What a round trip preserves, it does not invent: a sliced array comes back
  with the names the slice kept, and a product or a transpose comes back with
  none, per issue #561. copy() and copy.deepcopy() are unchanged and still
  drop the names — they go through numpy's own copy hooks, and libroadrunner
  drops them there too.

• A NamedArray column name never points at a column it does not name, so a
  sliced as_roadrunner array stops answering a lookup with its neighbour's
  species (issue #561). colnames labels columns by position, and
  __array_finalize__ copied the parent's list onto every array numpy derived
  from a labeled one. A plain column slice — arr[:, 1:], on the "all numpy
  operations work" path Result.as_roadrunner advertises — therefore came back
  two columns wide still carrying three names, with the __new__ invariant
  len(colnames) == shape[1] violated on a live object. Lookup resolved the
  name against the stale list and indexed the sliced data with it:
  sub["[X]"] returned the [Y] trajectory, sub["time"] returned [X], and
  only the last name ran off the end and raised. A PyBNF-style fit taking the
  time column off the front could feed the wrong species into a residual with
  nothing to read as an error.

  Names are now attached only where the columns behind them are known.
  __getitem__ computes which columns a key kept and relabels the result, so
  arr[:, 1:], arr[:, ::-1], arr[:, [2, 0]], a boolean column mask, a row
  slice and a chained slice all come back correctly labeled — more than the
  reference implementation does, since libroadrunner returns colnames == []
  for every one of them. Anything else numpy hands back — a product, a
  transpose, np.sort(arr, axis=1), a copy, an unpickled array — carries no
  names, because this hook is told nothing about the transform that produced
  the array and a same-shape result is as likely to have had its columns
  permuted as left alone. Lookup on such an array raises KeyError naming the
  cause and the way to carry names across, instead of resolving against a
  column that moved; libroadrunner drops the labels on those too. The numbers
  are untouched in every case — only the labels changed.

• A cyclic function graph is solved instead of swept, so the RHS stops being a
  function of how many times it has been called (issue #621). Functions that
  read each other are a simultaneous system. evaluate_functions() walked them
  once per RHS evaluation in whatever order the sort's cycle fallback left —
  a Gauss-Seidel sweep with no convergence check and no fixed number of steps —
  so the value a function held depended on the call history. Model.rhs(y, t)
  returned a different number on every call at the same state, and
  compute_derivs was not a function of (t, y) at all. Below a loop gain of 1
  it crept toward the answer without arriving (-1.0, -1.75, -1.9375, … for a
  system whose solution is -2); above 1 it ran away.

  build() now groups the functions into strongly-connected components
  (Tarjan), and evaluate_functions() walks the groups in dependency order. A
  group of one — every function in every corpus model but one — is a single
  evaluation, exactly the walk it has always been, so nothing acyclic changes or
  pays. A group of two or more is solved by Newton on the residual
  g(x) = F(x) − x, seeded from a constant so the answer is a function of the
  state alone rather than of the previous call's answer.

  Newton rather than iteration is forced by the one model that needs it.
  MODEL1006230117's group is linear in its unknowns (R reads LR, LRG,
  RG; each of those reads only R), so the system has one exact solution —
  and a linear system is precisely where Newton lands on it whatever the gain,
  while fixed-point iteration diverges for any gain above 1. That model's gain
  is ~36: it reached inf before t = 1e-13 and could not be simulated at
  all. It now integrates, and its R resolves to the closed-form
  R_Total/(1 + LR/R + LRG/R + RG/R).

  The compiled path cannot express the solve, so codegen declines a cyclic
  function graph and the caller falls back to the interpreted engine. It had
  been emitting func[i] = … in an order that cannot exist — 12 use-before-def
  reads of an uninitialised func[] on that model, where the compiled run died
  on a non-finite RHS at t = 0 while the interpreted one was correct. Two
  engines disagreeing in silence is worse than either being wrong alone.

  Acyclic behaviour is unchanged and cross-checked: on function-heavy models
  (8–37 assignment rules) the newly-ordered interpreted engine and the
  untouched compiled emitter agree to 8e-12 or better, and all 2,774 .net
  files and 1,291 vendored SBML documents still load.

• A parameter defined through a reference cycle, or in terms of itself, is
  refused at build time instead of being assigned a number (issue #617).
  x = y+1 with y = x+1 denotes no value and s = s*2 denotes none either,
  but both used to build. What came out was worse than an error: with no
  topological order to evaluate the cycle in, the fallback order ran and the
  values were manufactured from the front end's seeds — the same model loaded
  x = 2.0, 1.0 or 101.0 on three different seed pairs — and each later
  re-derivation relaxed the cycle one more step, so writing an unrelated
  parameter moved x by a fixed amount every time. A rate constant that depends
  on how many writes have happened is not a model, and it reached the solver:
  the run completed and returned an ordinary-looking trajectory. A
  self-reference was quieter still, references_model_symbol skipping the
  parameter's own name so the row was demoted to a plain Constant holding
  seed*2 with the expression dropped.

  This is the residue issue #568 left. That fix made an acyclic chain converge
  in one pass; a cycle has no topological order, so the drift it removed
  everywhere else survived in the one place the sort could not help. The
  refusal names the cycle in reading order — x -> y -> x, each reading the
  next — rather than just listing the parameters involved, and it names the
  cycle rather than a parameter merely downstream of one, which Kahn's
  algorithm also leaves unplaced and which is perfectly well defined.

  It matches both reference implementations, which is what makes it a fix
  rather than a preference: BNG2.pl refuses every parameter dependency cycle,
  solvable or not (ABORT: Parameter y has a dependency cycle y->x->y, and
  ABORT: Parameter s is defined recursively), and run_network refuses the
  same .net. A solvable cycle (a = 10-b with b = a, i.e. a = 5) is
  refused too: finding it needs a simultaneous solve, which one evaluation pass
  is not. Nothing in the tree is affected — all 2,774 .net files and all 1,291
  vendored BioModels SBML documents that loaded before still load, with the same
  32 unrelated refusals as before.

  Two things are deliberately exempt, because neither can denote nothing.
  A numeric literal is stepped over whole, so the e of 1e-5 is not read
  as a symbol — without that, a parameter named e or E (the enzyme in every
  Michaelis-Menten model; e appears in four .net files here) whose
  expression carried scientific notation was refused as defining itself. And a
  function-bound slot is not a definition at all: evaluate_functions()
  overwrites it every step, so its expression is only the seed it holds until
  the first pass, and refusing it would refuse the very shape the function sort
  admits — the .net spelling of an SBML <assignmentRule>, which emits both a
  # ConstantExpression row and a same-named function.

  A cyclic function graph deliberately still builds. SBML forbids one, but
  MODEL1006230117 in the vendored corpus has mutually recursive assignment
  rules — R = R_Total-LR-LRG-RG and LRG = (L_iso*R*Gs)/(K_H*K_C) — that are
  linear in their two unknowns and so a simultaneous system with a solution,
  which the per-RHS pass relaxes toward rather than invents. Refusing it would
  drop a real model whose equations mean something, and solving it needs
  machinery the sort is not; a parameter cycle has neither excuse. Until that
  linear case is actually solved, a function cycle keeps the path-dependent RHS
  GH #76 describes, now narrowed to cycles alone.

• A derived parameter that reads another derived parameter no longer keeps a
  stale value — and a stale rate constant (issue #568). Derived (expression-
  valued) parameters were re-evaluated in a single pass in declaration order,
  which silently assumes declaration order is dependency order. Nothing enforced
  that: ModelBuilder::add_parameter appends in call order, the .net reader
  hands rows over in file order, and build() defines every parameter as an
  ExprTk variable before compiling any expression — so a chain declared
  bottom-last built and ran with no complaint, one link at a time.

  With bb = a*3 declared before a = base*2, set_param("base", 5.0) moved
  a to 10.0 and left bb at its load-time 6.0 instead of 30.0, because the
  pass read a before reaching a's own row; a second, value-identical
  set_param("base", 5.0) then moved bb to 30.0, which is the same statement
  from the other side — the write was one relaxation step, not a solve. Load
  time had the same hole: on tests/data/derived_param_reverse_order.net
  (k = k2*2, k2 = k3*1.0, k3 = kbase*3, written top-down) the model built
  with k = 0.0 — the reader's seed for a row whose value column is an
  expression — so the one reaction had a zero rate constant, the species did not
  decay at all, and the run reported success. The finite-difference sensitivity
  probes carried the same defect: ∂S/∂kbase on that file came back
  -30/-60/-90/-120 against a closed form of -22.2/-32.9/-36.6/-36.1. Because
  a derived parameter is routinely a rate law, every one of these is a plausible
  wrong number rather than an error. The same file with its four rows written
  bottom-up was correct throughout.

  ModelBuilder::build() now sorts the derived parameters into dependency order
  once (SharedModelData::derived_param_order) and evaluates them in it, and
  the new NetworkModel::refresh_derived_params() is the single pass over that
  order. The loop it replaces had been hand-copied at seven sites — set_param,
  the interpreted and the codegen CVODES sensitivity RHS sync,
  restore_nominal_params, rederive_expression_params (reached from both the
  state-trigger and the event-sensitivity perturbation syncs), the
  switch-isolation rebuild, and the steady-state FD probe — each with its own
  spelling of the same three lines and its own local for the evaluator. That is
  what made this a seven-way bug rather than a one-way one, and why the fix is
  one function rather than seven edits: the ordering rule now has exactly one
  place to be wrong. The seven copies are nine call sites of the one method (two
  of them were helpers with two callers each), and
  CvodeSimulator::Impl::rederive_expression_params is gone, having become pure
  delegation. One call converges a chain of any depth, so the order the rows
  were declared in no longer changes any answer. This is GH #76's fix applied one field over — that
  issue was the identical single-declaration-order-pass defect for functions,
  and both now share one Kahn sort, which seeds ready nodes in ascending index
  so the output is deterministic and, among parameters that become ready
  together, in declaration order. A reference cycle (x = y+1 with y = x+1)
  is unsatisfiable however it is ordered, and is refused at build time — see
  the issue #617 entry above, which landed after this one and supersedes what
  it originally said here about such a model still building.

  No shipped model changes, and that is measured rather than assumed: of the
  2,774 .net files in the tree, 1,285 have at least one derived parameter and
  exactly zero declare one before a derived parameter it reads; of the 1,291
  vendored BioModels SBML documents that load, 382 do and again exactly zero are
  out of order. The only file the sort reorders is the new regression fixture.
  The SBML figure is not luck — the loader had been hand-sorting the parameters
  it lifts precisely to stay inside the old one-pass rule (its comments named
  the failure mode and BIOMD0000000569 as the case). That sort stays, but it is
  no longer what makes the numbers right.

• An SBML document that declares a required="true" package bngsim does not
  interpret is refused by name, instead of loading as its bare core layer
  (issue #592). SBML Level 3 multi is the package rule-based models are
  encoded in — species types, feature values, binding sites, reactions whose
  meaning comes from component maps. Model.from_sbml accepted a multi
  document, dropped the entire multi layer, and built a model out of whatever
  the core layer happened to contain: no exception, no warning, and a
  simulation of something that is not the model on disk. The loader consulted
  exactly one package (doc.getPlugin("comp"), for the GH #230 flattening) and
  read past every other one as if it were not there.

  The guard keys on required, not on a list of package names, which is what
  makes it general and correct: SBML defines required="true" to mean the
  package changes the model's mathematical meaning, so the same check covers
  qual, spatial and any future package, while the presentation-only ones
  (layout, render, fbc) declare required="false" and pass through
  untouched. comp and distrib are the handled list —
  bngsim._sbml_unsupported.HANDLED_PACKAGES, with the reason each is on it —
  and a package this libSBML has no extension for is checked too, from the
  document's unknown-package list. The refusal is a ModelError naming each
  package and its namespace URI, with the same
  BNGSIM_ALLOW_UNSUPPORTED_CONSTRUCTS=1 opt-out the delay / AlgebraicRule
  refusals have; under the opt-out it warns per package and loads the core
  layer.

  Nothing that loaded before stops loading: across the 14,885 XML documents in
  the tree and the full 1,824-case SBML semantic suite, the only ones the guard
  refuses are three cached BioModels qual documents (BIOMD0000000562 / 592 /
  593), each of whose core layer is empty — 0 species and 0 reactions against
  10–28 qualitative species and their transitions — and each of which bngsim
  had been loading as a model with zero species. The suite declares only
  comp:required="true" and
  fbc:required="false", so its grading is unchanged. Two libSBML behaviors
  would have made this refuse whole corpora and are handled explicitly: on a
  Level 2 document it reports its annotation-based layout / render plugins
  as required (there is no required attribute before Level 3), and on any
  L3V2 document it reports the l3v2extendedmath plugin — which carries the
  core namespace — as required.

• clear_codegen_cache() no longer deletes a running sharded compile's
  scratch directory (issue #557). Every sweep in bngsim.cache holds back an
  entry whose compile is still alive, by reading the PID out of the name and
  probing it with os.kill(pid, 0) — and clear documents that its zero
  min_age floor switches everything off except that check, "because emptying
  a cache is a reason to reclaim disk and not a reason to break a build that is
  in progress". For the sharded path the promise could never fire:
  _compile_sharded created its scratch directory with
  tempfile.mkdtemp(prefix="bngsim_shard_"), whose name carries no PID, so the
  check answered "not running" for every shard directory while is_removable
  said yes. clear_codegen_cache(), bngsim-cache clear and
  bngsim-cache clean --min-age 0 deleted the driver.c / unit_*.c / *.o
  out from under a live cc -c pool, and the concurrent compile died with
  "Codegen shard compilation failed" or "Codegen shard link failed". The sharded
  path is the many-core/HPC case of GH #160, which is exactly where a sibling
  worker reclaiming disk is most likely.

  The directory is now named bngsim_shard_p<pid>_<random>, the same way a
  serial partial carries .<pid>_<n>, and the liveness check reads either shape.
  The prefix is a _codegen constant bngsim.cache imports rather than a
  literal spelled in both: the two disagreeing with nothing to notice is the
  whole of this bug. A shard directory from a bngsim that predates this carries
  no PID and so is still collected, which is right for what it is — debris of a
  process that is gone.

• A manifest location can no longer name a file outside the extracted
  archive (issue #562). read_omex() refuses a zip member whose path escapes
  the extraction root, but the manifest half of the same guard was
  entry.location.lstrip("./") — which strips a leading character set
  {'.', '/'} and leaves interior .. segments in place. A <content location="./a/../../../../etc/passwd"> therefore resolved to
  root/a/../../../../etc/passwd, and load_model() / load_protocol() /
  load_full_protocol() opened that host file and fed it to Model.from_sbml,
  Model.from_net or the SED-ML reader. The only brake was the is_file()
  check, which stats the un-normalized path and so passes as long as the first
  component exists in the archive — which any archive guarantees by shipping
  that directory. Reading an untrusted .omex is this module's documented
  workflow (BioModels distributes archives this way), and the attacker controls
  both the manifest and the layout.

  Locations now go through one normalizer that drops . segments, resolves
  .. against what it has already consumed, and reports an escape.
  _parse_manifest refuses the archive outright, naming the location — one row
  pointing at a host file condemns the manifest the way one escaping member
  condemns the zip — and OmexArchive.path_of keeps its own containment check,
  resolved, since OmexArchive and OmexEntry are public and an archive can be
  assembled entry by entry without read_omex.

  Two things fixed in passing, both from the same lstrip: a location under a
  hidden directory lost the dot (./.hidden/model.xml was written into the
  zip as hidden/model.xml, on both the write and the read side, so it
  round-tripped under a name that was not the one in the manifest); and
  write_omex now refuses an escaping or empty *_location, so the writer
  cannot mint a manifest the reader must refuse. A ConversionError raised
  while parsing the manifest also cleans up the temporary extraction directory
  now, as every other read_omex refusal already did.

• build_model_from_parsed builds the eight .net files whose table
  functions and MM rate laws used to reach ExprTk unregistered, and names the
  two it will not build (issue #597). Model.from_net ran two post-parse
  steps that lived only in net_file_loader.cpp, so ten .net files in
  tests/data loaded under one documented loader and raised under the other —
  bngsim._net_reader, the reader bngsim.__all__ exports and
  docs/user-guide/loading-models.md documents. The failure was
  ModelBuilder: failed to compile function '__net_reader_func_0': Hill — ... ERR239 - Undefined symbol: 'Hill', which named neither the file, the
  reaction, nor the loader that does read it. The same two-loaders-disagree
  shape as issue #554, in its loud form.

  Two causes, and the reader had been losing data on both: it read only the
  rate column, so the operands after a legacy token — Sat k3 K4, MM kcat Km
  — were discarded and the bare token was left behind to be wrapped as if it
  were an expression.

  • tfun(...) table functions. The loader reads the spec out of the
    functions line and registers the table before build() compiles the
    expression that calls it, and points relative .tfun paths at the .net
    file's own directory — where BNG writes the table. That step is now
    bngsim::net_function_tables(), called by the loader and bound for the
    reader, rather than a second reading of tfun(...) syntax in Python to
    drift from the first (ModelBuilder.set_net_file_dir,
    add_table_function_spec and add_inline_table_function_spec are bound
    too). Whole-body and embedded forms both, so wrap_single.net's
    (tfun('drive.tfun',time)+5)/k_scale keeps its wrapper arithmetic.
  • MM kcat Km needed no rewrite at all, only parsing. ModelBuilder has
    the Michaelis-Menten rate law and takes its two parameters as one
    "kcat,Km" string; the reader now hands it that, and a short MM line
    falls through to the elementary branch exactly as it does in C++.

  All eight build a model identical to Model.from_net's — same
  codegen_data(), same table-function names, same trajectory — because loading
  was never the bar (#554 was two loaders that both loaded and disagreed by a
  factor).

  The deprecated Sat and Hill tokens are refused by name instead.
  Their rewrite is ~150 lines of C++ that also allocates non-colliding function
  and observable names and emits the deprecation warning Model.from_net raises
  as a UserWarning; reimplementing it in the reader would be the second dialect
  #554 was about. build_model_from_parsed now raises ValueError naming the
  reaction, the token and its operands, Model.from_net as the loader that
  rewrites them, and the explicit rate law to write by hand. parse_net_file
  still reports them faithfully — type="legacy" with the operands in
  legacy_constants, so the interchange dict stays inspectable — and its
  reaction dicts gain legacy_constants, the dict gains net_file_dir, and
  type gains "mm" and "legacy".

• A .net functions line that is unindexed, or written name() = expression,
  is refused instead of dropped or read with = as the function name (issue
  #606). The functions block was the last of the five .net block parsers
  still discarding or mis-reading a line it did not recognise, and it did both.
  rate_fn() (kcat*Atot)/(10+Atot) has two fields, so it fell through the
  three-field gate and vanished; rate_fn() = (kcat*Atot)/(10+Atot) — what
  writeFile({format=>"net"}) writes — has three, so it passed the gate and the
  name was read from the fixed second field, which is the =. bngsim invented a
  function literally called =, carrying the right expression under the wrong
  name. Either way the C++ loader then failed against the wrong block
  (reaction 0 (Elementary) references unknown parameter 'rate_fn' — the
  reaction was fine), and parse_net_file did not fail at all: it handed the
  caller the wrong functions list with no error and no warning.

  run_network refuses both shapes (ERROR: Found invalid line "rate_fn()" while reading functions block), so there is nothing new to support here:
  the functions block is the one block where BNG2.pl's two writers differ in the
  separator and not just in the leading index, and writeFile's functions
  block is not a network any BNG reader loads. A functions line must now be
  <index> <name>() <expression> — the index is mandatory here, as in the
  reactions block (issue #603) and unlike species, parameters and groups (issue
  #601); the name field must be an identifier, so = is refused as a name; and
  an expression that starts with = is refused too, since that is the same
  writeFile shape with an index added by hand, which was silently wrong in the
  Python reader and an unattributed ExprTk ERR248 in the C++ one. Both
  documented loaders raise the same message for the same file (issue #554). The
  corpus is unaffected: all 13,690 functions lines across the 2,773 .net files
  in the tree are generate_network's indexed, =-free shape.

• An expression-valued .net species initial concentration is evaluated
  instead of refused (issue #600). parse_species() read the IC column as a
  number or a bare parameter name, so 1 A() 2*A0 — which run_network
  evaluates to 200 — was rejected. It is now lifted into a synthetic
  _InitialConc<N> parameter and the species points at that, which is exactly
  what BNG2.pl's own generate_network does with a BNGL seed-species
  expression: it writes _InitialConc1 2*A0 into the parameters block and the
  name into the species column. Reusing that shape means the lift rides
  machinery that already exists — the parameters block compiles expressions
  through ExprTk, build() re-resolves the species IC from the evaluated value
  (issue #79) and seeds forward sensitivities from the reference, and since
  issue #602 it refuses an expression naming a symbol the model never declares.

  A token that is a name and nothing else is still refused with the message
  that names the species (issue #571), rather than being lifted: B0_typo is a
  reference that should have resolved, not arithmetic to evaluate, and the
  clearer error is worth keeping for by far the likelier mistake. A typo inside
  an expression — 2*B0_typo — cannot be told apart by that rule, so it reaches
  the parameter compile and is refused there. Both documented .net loaders do
  this identically, down to the synthetic parameter's name, so parse_net_file
  and Model.from_net still report the same numbers (issue #554).

• A .net parameter whose expression names an unknown symbol now refuses the
  model instead of silently keeping the numeric prefix of the malformed text
  (issue #602). ModelBuilder::build() compiled expression-valued parameters
  inside a try whose catch swallowed every failure with continue, annotated
  "May reference functions — handle below". There is no below: nothing
  re-compiles a parameter skipped there, and there is nothing to re-compile for,
  since every function already owns a parameter slot and all parameters are
  registered as evaluator variables before the loop — a function name is a
  symbol the compile can already see. What the swallow actually caught was an
  expression naming something that exists nowhere, and it left the parameter
  holding whatever the front end's partial std::stod had scraped off the front.

  So the value was not merely zero, which is the shape that made this dangerous:
  k1 = 2*kbse (a typo for 2*kbase) came back as 2.0, and the model
  integrated to completion at twice the intended rate with no error and no
  warning — a plausible-looking rate constant rather than an obvious zero, on a
  parameter that feeds rate laws and not just initial conditions. The two
  documented .net loaders also disagreed about it, 0.0 from
  bngsim._net_reader against 2.0 from the C++ loader, which is the split
  #554 exists to prevent. run_network refuses the same file outright ("Could
  not find parameter kbse. Exiting."), and a function body that will not
  compile already threw a hundred lines further down; only the parameter was
  silent.

  A failed parameter compile now throws, naming the parameter and its
  expression, in the same shape as the function-compile error. The Python reader
  evaluates through a parameters-only ModelBuilder.build(), so it inherits the
  refusal and the two loaders agree again. parse_net_file's
  warn-and-report-0.0 path is removed with the defect it was compensating for —
  including the extra probe build it ran per suspect parameter — and the
  engine-free fallback already refused rather than substituting, so it is
  unchanged.

• A .net species or parameter line without a leading index is read instead of
  silently dropped (issue #600). parse_species() and parse_parameters()
  required three whitespace fields, which made the leading index mandatory.
  BNG2.pl treats it as optional — Perl2/SpeciesList.pm strips one with
  s/^\s*\d+\s+//, kept with the comment "Can't deprecate this because indices
  used in NET files" — and BNG2.pl writes both shapes into a file named .net:
  generate_network emits 1 A() A0, while writeFile({format=>"net"}) emits a
  bare A() A0, and likewise A0 100 in the parameters block. run_network
  reads either. bngsim skipped the unindexed form line by line, with no error and
  no warning. A discarded species block left a model that either failed against
  the wrong block (reaction 0 has reactant species index 1 out of range [1, 0])
  or, with no reactions, loaded clean as a zero-species model; a discarded
  parameters block left every species IC and rate law naming a parameter unable
  to resolve against declarations that had simply been thrown away.

  Both blocks now take the index when it is present and read the line without it
  when it is not, in both documented .net loaders (#554), and a line that
  cannot be parsed raises rather than disappearing. A
  writeFile({format=>"net"}) network loads to the same initial state
  run_network reports for it. The value column is also joined from the whole
  remainder of the line, so a spaced 2 * A0 is no longer truncated to its first
  token and silently seeded as 2.

  The groups and reactions parsers carried the same assumption, and an audit
  of both found it broken in opposite directions. run_network reads an
  unindexed groups block, but bngsim took the observable name from a fixed second
  field, so Atot 1 became an observable named 1 with no entries and the load
  died on a duplicate-symbol error that blamed the file for names it did not
  contain — refusing a .net BNG2.pl accepts, with a diagnostic pointing at the
  wrong thing. The groups index is now optional too, and the entries column stays
  optional so an observable whose pattern matches no species keeps loading with
  none (BNG2.pl writes a bare <index> <name> for those; 55 such lines exist in
  the parity corpus). A reaction line's index, by contrast, is genuinely
  mandatory — run_network refuses a reactions block written without it — but a
  short line was skipped rather than refused, so the loader built a model with no
  reactions at all that loaded clean and integrated a flat trajectory while
  reporting success. It now raises and names the line.

• A .net species whose initial concentration names an unknown parameter now
  fails loudly instead of loading as a silent 0.0 (issue #571).
  parse_species() reads the IC column as a number or a parameter name; on a
  lookup miss it took the default s.concentration = 0.0 and loaded anyway. A
  typo, a parameter declared after the species block, or an arithmetic IC such
  as 2*A0 (std::stod stops at the 2, so the whole token is looked up as a
  parameter name and misses) each turned a real initial condition into zero. The
  species then started at 0, every observable and downstream flux built on it
  was wrong, and the run reported success — while the same loader already failed
  loudly for an unknown rate-law parameter and an out-of-range observable index.
  The loader now raises species '<name>' initial concentration '<token>' is neither a number nor a declared parameter, surfaced as ModelError through
  Model.from_net and as ValueError through the core loader. The Python .net
  reader (bngsim._net_reader._parse_species), which #554 had deliberately made
  agree with the C++ loader on the 0.0 default, is changed in step so the two
  still agree — now on the error rather than on the wrong number. BNGL
  expression-valued seed species stay supported: BNG2.pl lifts each expression
  into a synthetic _InitialConcN parameter and writes that parameter's name
  — a single token — into the species block, so it arrives as a declared
  reference and still resolves through the expression evaluator. Only raw
  arithmetic written straight into the IC column — a form no producer emits — is
  rejected rather than given a second expression path.

• The interactive clock moves with the model's state, so a leg after a stop
  condition stops mislabelling its own time axis (issue #553). run_until
  assigned self._current_time = t after run() returned. Stop conditions are
  evaluated against a completed result, so run() raises StopConditionMet
  once the backend has integrated the whole interval and written the
  end-of-span concentrations back into the model — and that assignment was
  skipped. The simulator was left holding the t=100 state with its clock still
  reading 0.0, and the next run_until(t=50) then built t_span=(0.0, 50.0)
  over that advanced state and returned a Result whose time axis claimed to
  start at 0 while species[0] held the t=100 values. A silently wrong
  trajectory from two documented features used together, with no error or
  warning. StopConditionMet was the reliable way in; any post-solve raise did
  it, including SsaBoundaryWarning under -W error.

  run() now advances _current_time itself, the moment the backend writes the
  state back and before either post-solve step that can raise. The clock is no
  longer maintained by the caller, so there is no assignment left for an
  exception to skip, and current_time means what its name says on every path:
  the time the model's stored state is at. It is read off result.time[-1]
  rather than t_span[1], which is the same number for an ordinary run and the
  right one when sample_times or steady_state=True picks the final output
  time. parameter_scan rewinds the model when it finishes, so it now puts the
  clock back with it. The original repro's second leg reports
  t_span=(100.0, 150.0) over the t=100 state; its run_until(t=50) raises
  ValueError: Target time (50) must be > current time (100.0) instead of
  relabelling.

  Rolling the model back to the trigger point — the issue's other suggestion —
  is not implementable from what StopConditionMet carries. _truncate_result
  builds the partial Result with no core, so Result.state refuses on it, and
  its species block is the reported projection: on a model whose SBML event
  assigns a parameter, the promoted entry is missing and the row is narrower
  than the model's state vector (verified: species is (n, 1) against a state
  of width 2). A rollback would write a short or wrong state. The semantics are
  documented instead — on StopConditionMet, Simulator.current_time,
  run_until, and in docs/user-guide/simulation.md: a stop condition truncates
  the reported result and does not rewind the model, the clock names the end of
  the integrated span, and snapshot() / restore() is the way to continue from
  before the stop.

• parse_net_file evaluates a .net parameter expression with the engine's
  own evaluator, so the reader and Model.from_net seed a model from the same
  numbers (issue #554). The reader evaluated each expression with Python's
  eval. BNGL spells exponentiation ^, which Python reads as bitwise XOR, so
  A0 10^2 came back as 8 — and rad_cell^3 on float operands was a
  TypeError that a bare except Exception turned into 0.0, as were BNGL's
  if(c,t,f) and &&, and a parameter named lambda. The wrong number did not
  stay in the parameter slot: _parse_species resolves a species initial
  condition written as a parameter name through it, so build_model_from_parsed
  seeded the state from it while get_param reported the value the engine
  computes — the model integrating from one initial state and reporting another,
  with no warning, through two functions bngsim.__all__ exports and documents
  as the way to load a .net for inspection. Eleven of the 133 .net files in
  this tree were affected, eight of them in the initial state:
  energy_example1.net started every species at 0 instead of 2522.5 and up,
  prion_aggregation.net at 0 instead of 7500.

  Evaluation now runs through a parameters-only ModelBuilder, whose build()
  compiles and evaluates every expression exactly as the C++ .net loader's
  does. There is no second BNGL dialect to keep in step: the one evaluator
  answers both readers, and across all 133 files every parameter value and every
  initial state is now identical to Model.from_net's. The extra build is a
  parameters-only one and does not show up in the parse: 133 files still parse in
  339 ms, as they did before.

  A species initial condition that names a parameter is handed to the builder as
  that reference — add_species_param_ref, the call net_file_loader.cpp
  makes — instead of as a number resolved beforehand, so build() re-resolves it
  from the compiled parameter. That also records the (species, parameter) pair
  forward sensitivity seeds ∂y(0)/∂p from, which a model built this way never
  carried. parse_net_file reports the pairs under a new species_ic_params
  key; build_model_from_parsed treats it as optional, so a hand-built dict
  still builds.

  An expression the evaluator cannot compile — one naming a symbol the model
  never declares, or malformed — is still reported as 0.0, because
  Model.from_net loads the same file with the same zeros and the two readers
  should not disagree in the other direction. It is no longer silent: the reader
  names the parameter and its expression in a UserWarning, having re-tested
  each zero-valued expression against a NaN seed to tell "compiled, and the
  answer is zero" from "never compiled". Only the root cause is named, not the
  parameters downstream of it. No .net file in this tree trips it.

  parse_net_file is documented as working with no compiled extension present —
  its dict is the interchange format for handing a .net model to scipy,
  gillespy2 or a hand-written RHS — and still does. That path cannot reach the
  engine's evaluator, so it evaluates ordinary arithmetic, reads ^ as
  exponentiation (the one operator Python spells the same and means differently),
  and raises on BNGL syntax it cannot do rather than substituting a number: a
  ValueError naming the parameter, its expression and the reason. Refusing is
  the change there — a file with if() in a parameter used to parse to a
  plausible-looking 0.0.

• run_tests.sh now stands the repo up instead of hiding it, so the suite
  passes there as it does in place (issue #590). The script relocates the tests
  so the source tree's python/bngsim/ cannot shadow the installed package it
  exists to exercise. It did that by copying them into a bare temp directory —
  which hid far more than python/bngsim/. Sixty-nine modules resolve
  tests/data, benchmarks/ or scripts/ by walking up from __file__, and
  from a bare temp directory that walk-up leaves the repo entirely: 64 failures,
  181 errors, and 359 skips reported as "not in this checkout (installed
  package)" while the source tree sat right there. None was a defect in bngsim.

  The temp directory is now a stand-in for the repo: python/tests/ copied to
  <tmp>/python/tests/, and every other top-level entry symlinked in. A walk-up
  then lands on something that is the repo, while <tmp>/python/ still holds
  tests and nothing else — so the one directory the script exists to hide is
  still hidden, and nothing else is. Linking costs nothing however large the tree
  is, and pyproject.toml comes along with it, so pytest finds the project
  config there too: [tool.pytest.ini_options] applies and the run is configured
  exactly as an in-place one rather than silently diverging.

  Four modules needed fixing, each a defect the old layout had masked rather than
  a concession to the new one:

  • test_core_stub_covers_bindings.py guarded on src/_bngsim_core.cpp but
    read python/bngsim/_bngsim_core.pyi. Two different files, one guard; a
    tree with one and not the other passed the guard and died on a
    FileNotFoundError. Its siblings already skip on the stub.
  • test_rebuild_editable_source_digest.py, the same shape: guarded on
    scripts/rebuild_editable.py, loads python/bngsim/_build_provenance.py.
  • test_primary_param_names.py rglob-ed the whole repo root for .net
    fixtures, which also swept in .pytest_cache — so "every .net in the tree"
    was partly an assertion about the developer's own cache. It names the three
    roots that hold fixtures now, which is also what makes it work through a
    symlink: Path.rglob will not descend into one, though it follows one at
    the root of the walk.
  • test_sbml_reversible_split.py counted parents[3] to find a sibling
    PyBNF checkout. Its own docstring records that count having been off by one
    before, which made the skip unconditional and the test had never run
    anywhere; relocation would have done it again.

  And the script no longer names its working directory TMPDIR. Where that
  variable is already exported — launchd exports it on macOS, and many Linux CI
  images set it — assigning to it redirects every child's temp files into the
  stand-in, because Python's tempfile reads it. pytest's tmp_path then lands
  inside a tree whose root carries bngsim's pyproject.toml, so a walk-up out of
  any temp directory "finds" a source root that is a link farm, and the per-run
  artifact cache lands inside rootdir, which test_artifact_cache_isolation.py
  asserts against. The script had done this since it was written; it surfaced
  only once the stand-in became the rootdir.

  The link farm needs real symlinks, which Windows grants only under Developer
  Mode or an elevated token — ln -s under Git Bash/MSYS copies instead. A
  silent deep copy would duplicate build/, third_party/, .git and .venv,
  so the script checks one entry and exits non-zero naming the fix. The full
  suite runs on ubuntu and macOS; test_run_tests_sh.py skips on Windows,
  because what it tests is a bash script that builds a symlink farm.

  One thing the stand-in must not do is make python/bngsim/ unreadable
  (issue #594). What shadows the installed package is that directory being on
  sys.path; its bytes being readable shadows nothing. Eleven tests read it as
  files — the committed _bngsim_core.pyi, and _build_provenance.py, which is
  loaded by path and asserts that doing so imports no bngsim — and they resolved
  it by walking up from __file__, so under the stand-in they skipped themselves
  as "not in this checkout" while the script was running from inside the
  checkout. test_core_stub_covers_bindings.py,
  test_rebuild_editable_source_digest.py and
  test_rebuild_editable_feature_options.py take the source root from the rig
  now, as test_version_consistency.py already did — which is why that one read
  the same stub and never skipped. Their .exists() guards stay: those are right
  against a real wheel checkout, which is what their reasons describe.

  Measured, BNGPATH unset both ways: 5981 passed, 62 skipped, 0 failed,
  exit 0 under the script, and the same 5981 / 62 / 0 in place — every test the
  suite has, executed either way. Before #590: 5371 passed, 64 failed, 181
  errors, 421 skipped, exit 1.

• run_tests.sh runs the suite it collects, instead of aborting on six import
  errors having executed nothing (issue #578). The script copies the tests to a
  temp directory so the source tree's python/bngsim/ cannot shadow the
  installed package it exists to exercise — but the copy was test_*.py plus
  conftest.py, which is not what python/tests/ holds. Left behind:
  __init__.py, without which test_sbml_psa.py's
  from .test_ssa_psa_volume import … has no parent package to resolve against;
  and the sibling helpers _ci_workflow.py and _extrande_reference.py, which
  four more modules import at module scope. All six raised at collection time,
  and pytest aborts the run on a collection error — so the script exited 2 having
  executed none of the 5775 tests it had just counted.

  It now copies every .py, into a directory named tests. The name is
  load-bearing: with __init__.py present pytest derives the package name from
  the directory and puts its parent on sys.path, and a mktemp name like
  tmp.AbC123 is not a legal identifier. tests reproduces the source layout,
  and the parent it adds is the temp root rather than the repo's python/ —
  which is the shadowing the whole exercise exists for.

  Fixing only the copy would have left the worse half in place, 168 tests' worth.
  Five more modules did not raise: they resolved a repo path by walking up from
  __file__, which under the copy lands outside the repo, so their globs came
  back empty. Two skipped themselves whole (test_ssa_baseline_regression.py,
  test_ssa_biomodels_parity_zgate.py — 24 tests) and three quietly collected
  fewer parametrized cases: test_benchmark_runner_help.py 127 down to 11,
  test_ci_run_list_coverage.py 21 to 6, test_codegen_structural_key.py 37 to
  24. A green run that had dropped a fifth of what it should have checked, saying
  so nowhere. The script now also exports BNGSIM_SOURCE_ROOT, and those modules
  resolve the repo through the new python/tests/_source_root.py:
  $BNGSIM_TEST_DATA's parents, then $BNGSIM_SOURCE_ROOT, then the walk-up,
  each candidate confirmed by finding bngsim's own pyproject.toml there so a
  wrong guess cannot be returned as a right one — the search four other modules
  already carried inline. A copied run now collects exactly what an in-place run
  collects, 6029 either way, and test_run_tests_sh.py asserts that equality
  rather than merely that nothing errored.

  That still left the run red — 64 failures, 181 errors and 359 spurious skips,
  every one of them a module resolving tests/data, benchmarks/ or scripts/
  by a walk-up the copy invalidated. Closed by the entry above (issue #590),
  which stopped hiding those directories rather than teaching 69 modules to find
  them.

• A golden regeneration whose variable blew up now fails the gate instead of
  scoring a perfect match (issue #572). _core.fingerprint is what a consumer
  regenerating our golden corpus checks itself against: checksum must match
  byte-for-byte on the same platform, and fingerprint_max_rel is the
  cross-platform fallback when a different BLAS or rounding legitimately moves
  the last digits. Both forgave a blow-up, in the one direction a gate must
  never fail:

  • fingerprint_max_rel returned 0.0 — perfect agreement — for a variable
    that went NaN. It scored each stat as abs(x - y) / max(abs(y), abs_floor)
    and kept the running worst with max(worst, …), both of them Python's
    builtin max, which keeps its first argument when the comparison is False.
    A NaN stat makes nan > 1e-9 False, so the denominator is NaN, the ratio is
    NaN, and max(0.0, nan) is 0.0. No argument order helped, and the
    structural guard could not catch it either: a blown-up column has the same
    variable set and the same n_time as a healthy one — only the values went
    NaN. Since the fingerprint is consulted only after a checksum mismatch,
    which a blow-up always produces, this was the entire cross-platform verdict.
    A stat that is non-finite on exactly one side, or a different non-finite on
    each side (+inf vs -inf, inf vs NaN), now returns inf — above any
    tolerance a caller can pick, the way differ.deterministic_verdict makes a
    one-side-non-finite cell an unconditional hard fail. The same blow-up on both
    sides is agreement and still scores 0.0.
  • The checksum could not tell the three blow-up modes apart. _round_sig
    canonicalized every non-finite cell to NaN, so
    checksum(+inf) == checksum(-inf) == checksum(NaN). A run that blew up one
    way byte-matched a golden that blew up another and was certified an exact
    reproduction by the strongest check, which ends the comparison on a match.
    Each mode now keeps its own value, canonical rather than the input's own
    bits (a NaN payload is not portable). Finite data is untouched: every
    committed golden reproduces its recorded checksum unchanged, and a test
    re-hashes the committed trajectories to keep it that way.

• Forward sensitivity to the onset of a logistic in time runs, instead of
  handing back NaN at the first sensitivity RHS call (issue #549).
  BIOMD0000000554 and BIOMD0000000627 write their schedules as
  1/(1 + exp(-k·(t - (onset + …)))), and ∂/∂onset of one comes out
  A·exp(x)/(B·exp(x) + 1)^2 — bounded by A/4B² for every real x, and
  inf/inf = NaN the moment x passes ln(DBL_MAX) ≈ 709.8, which 4.59186·199
  at t = 0 does without trying. Both models failed at t = 0 for a term whose
  true value is 1e-397, while their trajectories matched the reference, since
  1/(1 + inf) is a clean 0 and only the differentiated form fails. The GH #388
  rewrite that divides such a ratio through by its exp exists for this exact
  shape — it names BIOMD0000000554 — and had never fired on either model, for two
  independent reasons, each sufficient on its own:

  • A constant shift in the exponent folds into a coefficient.
    exp(-4.59186·(t - (t_0 + t_1 - 3))) multiplies out to an exponent with a
    number in it, and sympy folds that number into a factor out front, so the
    denominator summand is 960856.16·exp(u) — a product, where the match wanted a
    bare exp. A summand may now carry a numeric coefficient c: it divides
    through as f does against rest/c, with the quotient scaled by c^-m, which
    is exact for any nonzero c and moves no pole, the denominator vanishing where
    it always did. A symbolic coefficient is still refused — it is a zero the
    original expression is perfectly happy at.
  • sp.Float(2.0) == 2 is False. The helper deciding whether one exponent
    is a whole multiple of another compared a sympy number against a Python int,
    so it answered "not an integer" for every Float it was ever handed, and had
    since GH #393 introduced it. libSBML writes its coefficients 1.0, so an
    SBML-sourced exponent is -1.0·sr_GLY·(t - to - to_GLY) and the ratio that
    decides the match is Float(1.0): the rewrite declined every SBML-sourced
    sigmoid there is. Compared as doubles now, which is what the helper's docstring
    always claimed; 3/2 and 2.5 are refused exactly as before.

  Both misses were also in the native saturable differentiator
  (bngsim._saturable_jacobian), which emits the analytical Jacobian and the
  J·yS half of the analytic sensitivity RHS without sympy and so cannot inherit
  the fix — narrower there, since it parses text and folds nothing, so it takes a
  law that writes the coefficient itself (1 + 5·exp(w)), but the same NaN one
  term later. Fixed alongside, as GH #336 and GH #388 both had to. An expression
  whose summand carries no coefficient emits the text it always did, character for
  character, through either emitter. The two models now return finite sensitivity
  tensors that agree with a central difference of the trajectory to the
  difference's own accuracy at every output point except the one landing exactly on
  a discontinuity the parameter moves, which is issue #545's territory and
  unchanged here.

• A forward-sensitivity failure no longer sends the reader after a pulse that is
  not there (issue #549). The note issue #547 appends whenever the non-finite
  half is ∂f/∂p described a power s^(a-1) rising from an onset, and the comoving
  column of issue #545 — the only shape it knew, asserted as the diagnosis. An
  overflow is the other way a finite rate law hands back a non-finite derivative,
  and a reader who had hit one was sent looking for a threshold the run had not
  reached. The message now names both routes and says which is the likelier when
  the column goes non-finite at the very first call. It is said only where a value
  actually went non-finite: an overflow leaves a NaN, and so always leaves a
  witness, which is precisely what the "step gave out just after a restart" hint is
  the absence of — that message is unchanged.

• Forward sensitivity to the onset of a pulse that rises as s^(a-1) runs, and
  is accurate, for 1 < a < 2 (issue #545). Such a pulse is continuous, but its
  ∂f/∂on goes as s^(a-2): unbounded past the onset, with most of its integral
  within a few ulp of it when a is near 1. No polynomial step resolves that
  forcing, so the onset column failed at the onset, stalled just past it, or
  finished wrong — 2.2% on SIR_v4 at a_2023 = 1.2, with no warning. Past a
  crossing that its parameter moves at c = ∂t*/∂p, the column is now integrated
  as V = S + c·f, which obeys V' = J·V + β with β the derivative of f along
  the parameter and every clock together: the shift leaves the vanishing base
  unchanged, so the singular terms cancel in β, and V is as smooth as the
  state. The generator emits β as bngsim_codegen_sens_rhs_comoving for each
  parameter whose crossing makes a power singular — a clock-dependent base that
  can reach 0, raised to an exponent between 0 and 1 or to one that is not a
  number — with c read off that base (_CODEGEN_VERSION 30). Every other model
  emits the source it did, including one whose onset is a logistic or a Gaussian.
  The solver switches a column to V at a #48 crossing or a #150 state switch
  whose ∂t*/∂p matches, back to S at the next restart, and reports
  S = V − c·f at every output, in the error floor, and in the carry-over seed. A
  run with events keeps plain columns, and BNGSIM_SENS_COMOVING=0 restores them
  for an A/B. Separately, a threshold written through clock guards — SIR_v4's
  t > if(t<365, d_2021, if(t<730, 365+d_2022, …)) — is now a clock crossing per
  branch, stopped on exactly. It used to be a state switch that CVODE roots on, and
  CVODE tests for a root only on a step it accepts: where every step spanning the
  onset failed the error test, the run wedged one ulp short of a root it never
  reached, with the comoving column or without it. On the issue's pulse —
  time() and counter clocks, >= and > windows, two pulse heights, a from
  1.05 to 1.9, four tolerances — 134 of 160 cells failed and none do, every column
  within the tolerance the smooth pulse gets of the closed form. On SIR_v4's
  d_2023 column against a Richardson central difference (on rows more than two
  days from the onset and the end, where the difference itself is not accurate),
  a_2023 = 1.2 goes from 2.2% at rtol/atol 1e-8/1e-8 and a stall at every tighter
  setting to 1.9e-4, 1.9e-6 and 2.0e-6 at 1e-8/1e-8, 1e-10/1e-12 and 1e-12/1e-14,
  the last two the reference's own agreement; the vendored a_2023 = 1.5 goes from
  a 2e-6 plateau to 7.7e-9 at 1e-12/1e-14.

• A forward-sensitivity run that fails on a derivative says so, instead of
  blaming a species or a rate jump (issue #545). A pulse such as
  k0 + k1*s^(a-1)*(1-s) over a window opening at an onset is finite and
  continuous, but for 1 < a < 2 its derivative with respect to the onset goes
  as s^(a-2): infinite at the onset, and unbounded just past it. When the
  sensitivity RHS came back non-finite at the onset, the message advised
  constraining a species, though every species and every rate law was finite.
  When the step gave out just past a strict onset, it blamed an if() rate jump,
  or, as CV_ERR_FAILURE, gave no reason at all. The non-finite witness now
  records which sensitivity column went non-finite, in which rows, and — from the
  #177 term scale — whether ∂f/∂p itself is the non-finite half; where every
  rate law and species is finite, the message names that column and describes the
  shape and its remedies. With the analytic sensitivity RHS in use, a stall, or an
  error-test or convergence failure, at the time the run last restarted (a
  crossing it stopped on, an event, a root) now says the step gave out where the
  run had just restarted, unless such a witness already explains it. The domain
  advice is unchanged wherever a rate law or a species really is non-finite or
  negative. The sensitivities guide describes the shape under "A rate that rises
  from zero at the crossing". The sensitivity RHS keeps evaluating the true
  derivative at the onset itself: a finite value there would remove the failure
  at the onset but not the unbounded forcing just past it, and on some models
  would turn a run that fails into one that finishes with a wrong column and no
  warning. The comoving column in the entry above integrates past the onset
  instead, so these messages are now for the runs it does not reach.

• SIR_v5 keeps its analytical Jacobian and its analytic sensitivity RHS: a
  power of an if() whose other branch is 0 differentiates to the step it is
  (issue #541). sympy distributes a power over the branches of a Piecewise,
  so the seasonal pulse if(c, u, 0)^(a()-1) holds 0^(a()-1) off-season, and
  sympy's power rule, written for a nonzero base, turned every derivative
  through that node into nan, although 0^e is a step in e and constant on
  each side. With a() an if() chain on the day counter, ∂/∂t declined the
  Jacobian and ∂/∂a_2021 the sensitivity RHS. Such a power now differentiates
  to 0 wherever its exponent's values are run-constants; one whose exponent
  reads the state in a value is a step nothing locates, and still declines.
  Three defects behind it are fixed as well. The zero-base logarithm guard (#310)
  builds Eq(base, 0) & (exp > 0), which sympy folds into an ITE when the
  base or the exponent is a Piecewise, and the emitters normalized booleans
  only before the guard ran. So a Hill exponent selected by an if() under a
  logarithm (k*A^if(time()<5,n1,n2)*log(A+2)) declined the Jacobian over an
  ExprTk compile error, and a forward-sensitivity run raised on a C compile
  error — BIOMD0000000276 and BIOMD0000000277 among them, which now run on the
  analytic RHS and agree with CVODES' difference quotient to 2e-9. And GH #96's
  removable-denominator rewrite cancelled (u·r)^n/u only for a bare-symbol
  u, so at a season onset (u = t − t_start) the derivative was 0/0: a run
  whose counter started on an onset stopped at the first sensitivity RHS call
  (CV_FIRST_SRHSFUNC_ERR), where the difference quotient ran. The rewrite now
  accepts any factor of the base. Last, a rewrite can fold a non-finite atom of
  its own: the guard's Eq(base, 0) over a season window whose bounds are if()
  chains ending in 0 folds to zoo*t on the last branch, and the emitters looked
  for non-finite atoms only before rewriting, so that build failed to compile.
  They now look again afterwards and decline, and the sensitivity RHS names the
  fold. On SIR_v5 the analytic columns match central differences of the
  trajectory to 1.6e-5, and a sensitivity run at rtol 1e-11 takes 0.35 s where
  the difference quotient took 409 s. SIR_v4, the issue's
  other model, divides by its shape parameters after its last modelled year
  (if(t<=1461, a_2024, 0)), so its rate law is non-finite there before
  anything is differentiated. It still declines — sympy's unsigned zoo does
  not follow IEEE arithmetic, so a derivative emitted as non-finite could stop a
  solve that finite differences complete — but both reasons now name the
  branch (… which the rate law itself already has where (t >= 365) and … and (t > 1461)) instead of blaming the derivative, and the sensitivity RHS no
  longer calls a non-finite atom a "non-differentiable or unsupported
  function". Ending the model's four year-selection chains with their 2024
  values restores both derivatives; that model-side fix is issue #543, whose
  re-source is under Changed. Across
  the 1,909-model rr_parity + ode_fullnet corpus with the
  derivation budgets off, 64 models reach a changed path, and 40 sampled others
  hand C++ and the C emitters byte-identical terms and source. No attach or
  sensitivity verdict was lost: SIR_v5 gained both, and BIOMD0000000429 and
  BIOMD0000000568 gained the analytic sensitivity RHS, which matches the
  difference quotient to 8e-6. The widened rewrite changes the emitted text of
  62 models; over 1,694 changed derivatives at 9,782 random points no value
  turned non-finite, and where the two forms differ by more than 1e-9 relative
  the derivative nearly cancels, and at 50 digits both equal the unrewritten
  derivative. MODEL1006230090, whose unbudgeted derivation runs for minutes, was
  not compared, and BIOMD0000000385 is not in the random-point count for the
  same reason. Where the widened rewrite applies it costs about 0.2 ms per
  derivative (BIOMD0000000628, measured in one process against main's); where it
  does not, nothing measurable (BIOMD0000000469's 1,559 derivatives take 0.04 s
  either way). The four affected models whose unbudgeted attach takes more than
  20 s decline at the default derivation budget on main and here alike.

• The SSA/PSA rejection of a delayed event points at the issue that tracks
  delay support, not at a planning document that is no longer in the tree
  (issue #526). The message named the event and the remedy, then sent the
  reader to "Phase 5c of the SBML SSA Support Plan". It now names issue #526,
  and still says to use method='ode', which runs delayed events, or remove
  the delay. Delay support under SSA/PSA itself is still open in #526.

• A rate law shared by several reactions is derived once per model, so the
  linlog BioModels BIOMD0000000469–473 keep their analytical Jacobian on any
  machine (issue #532). What a Functional rate law derives to depends on its
  text and on inputs fixed for the whole model — the function map, observable
  groups, species metadata, constants and writable volumes — and never on which
  reaction carries it: each reaction's stoichiometry is applied afterwards, by
  the C++ scatter at attach and by the C scatter in the compiled Jacobian. Yet
  attach_functional_jacobian derived it again for every reaction. Sharing is
  common: a .net function drives every reaction its rule generated, and the
  SBML loader emits a reaction with non-integer stoichiometry as one reaction
  per species, each carrying the whole kinetic law. The five Smallbone 2013
  linlog models have a biomass reaction with 65 reactants and 3 products; it
  becomes 68 per-species reactions sharing one 5,413-character law, and deriving
  that law 68 times was 92% of a 29 s derivation. Since #531 the models ran
  into the 20 s budget and attached or declined by machine load. The attach now
  derives each distinct rate law once per path and hands the terms to every
  reaction that carries it; a decline still ends the attach with the same
  reason. The compiled Jacobian's reconstruction
  (_codegen._functional_jacobian_groups, shared with the compiled sensitivity
  RHS) does the same, which matters because it only runs once the attach has
  succeeded: on BIOMD0000000469 it re-derived for 30 s. Both caches live for
  one derivation of one model. Measured one model at a time under the default
  budget, BIOMD0000000469 and 470 now attach in 5 s and 3 s where they declined
  at the budget, 471–473 attach in 2–4 s instead of 19–22 s, and their compiled
  Jacobians are emitted in under 3 s instead of 17–20 s. Across the corpus, 94
  models have reactions that share a rate law; for the 93 compared with the
  budget off, the terms handed to the C++ attach and the emitted C are
  identical to what main produced, and no attach verdict changed. The 94th
  shares only a constant zero rate and was not compared, since its other rate
  laws derive for minutes with no budget. mt_music_sequencer, the other model
  #532 names, spends its budget on a few genuinely expensive rate laws that no
  two reactions share, so it still declines at the budget.

• The stale-binary guard no longer calls a current extension STALE after a
  checkout, stash or rebase rewrites the C++ with identical bytes (issue
  #528). The guard compares the loaded extension's mtime with the newest
  C++/CMake source, so the routine end of a C++ pull request —
  git checkout main && git pull after a squash merge, where checkout writes
  the pre-merge bytes and pull writes the merged ones back — printed The binary does NOT reflect current C++, failed the pytest preflight, and
  demanded a rebuild of a binary that was current. scripts/rebuild_editable.py
  now digests the sources the guard watches (third_party/ included) under
  its rebuild lock before the configure and again after the install, and when
  the two agree leaves a record beside the installed extension: that digest
  and the SHA-256 of the extension file. When the mtime comparison says stale,
  the guard reads the record, and if it names the loaded extension's exact
  bytes and the sources as they are now, the verdict is fresh and the banner
  says fresh (digest). The mtime path is unchanged; the hashing, about 25 ms,
  runs only where the guard would otherwise demand a rebuild. Sources that
  change during a rebuild leave no record, and an extension built or copied
  any other way has none that names its bytes, so both keep the mtime verdict.
  The digest is not baked into the extension at configure time the way
  __build_commit__ is: cmake --build recompiles an edited source without
  reconfiguring, so a configure-time digest could vouch for a binary compiled
  from sources it does not describe.

• The C++ attach gate's verdict reaches Model.analytical_jacobian_status
  as text — the mismatching or non-finite entry, the probe state and both
  values — instead of a pointer at BNGSIM_JAC_DEBUG=1 (issue #534).
  NetworkModel::set_functional_jacobian returns a bool, and when it declined
  the reason existed only as a stderr line under BNGSIM_JAC_DEBUG=1, so the
  status read the same generic sentence for every refusal: 17 corpus models,
  in two classes invisible from the status — eleven whose analytical entry is
  non-finite at the seed while the RHS is finite (∂(k·sqrt(A))/∂A at
  A = 0), and six with a mismatch, five of them with their seed on a
  switching surface of their own if(c > 0) / if(c < 0) rate-law pairs,
  where the analytical entry is 0 at probe 0 and the finite difference is the
  smooth slope. The gate now records its verdict on the model — the kind of
  refusal, the reaction or the entry, the probe, the analytical and
  finite-difference values — readable as
  NetworkModel.last_functional_jacobian_decline() (a dict, or None after
  an attach that succeeded; carried by clone()), and the attach driver folds
  it into the status with species names: declined: the finite-difference self-check found a trustworthy mismatch at probe 0 (the seed state) in ∂f[LX()]/∂[LX()]: analytical 0, finite-difference -10, or declined: ∂f[B()]/∂[A()] is non-finite (inf) at probe 0 (the seed state) while the RHS is finite. Every other refusal is named too — a term that did not
  compile, with the compiler's message and the expression; an entry with no
  slot in the sparsity pattern; a reaction or observable index the model does
  not have — and the sparse self-check's non-finite line now carries the
  column as well as the row. The stderr line under BNGSIM_JAC_DEBUG=1
  stays.

• The analytical-Jacobian self-check tests finite-difference convergence
  across three step sizes, so a double-precision rounding plateau is no
  longer mistaken for a converged slope (issue #533). The gate judged an
  entry once its central differences at 1e-5·|y| and 5e-7·|y| agreed.
  BIOMD0000000393's rate laws carry the quadratic root
  (N/2)·(sqrt(1 + 4x/N) − 1) with 4x/N ≈ 1e-6, which doubles resolve to
  only 8e-10 relative; its rounding structure is the same at both steps, so
  the two differences agreed to 1e-6 while both sat 3% off the true slope on
  an entry whose value is a small residual of ten O(1) contributions, and
  the correct Jacobian was rejected at the third probe state. Rounding error
  in a central difference grows as the step shrinks while truncation error
  shrinks, so a third, coarser step at 2e-4·|y| moves off the plateau: the
  gate now requires all three differences to agree pairwise within the same
  5e-3 band before it judges an entry, in both the dense and the sparse
  path, and skips the entry as an FD artifact otherwise. A smooth entry
  converges across all three (truncation at 2e-4·|y| is ~1e-7 relative for
  O(1) curvature), so a genuinely wrong derivative is still caught there.
  BIOMD0000000393 attaches, with zero mismatches beyond 1% over ~3,300
  converged entries at each of its four probe states, and the extracted
  fixture tests/data/jac_selfcheck_rounding_plateau.net pins the
  mechanism.

• max, min and abs over a differentiation variable are
  differentiated as the piecewise functions they are, so the model keeps its
  analytical Jacobian (issue #507). sympy differentiates Max to
  Heaviside and Abs to re()/im(), which no emitter spells, so
  differentiate_rate_law declined the rate law and the whole model — 28 in
  the corpus, a Lorenz attractor among them — ran on finite differences. Now
  abs differentiates to sign(a)·a' through a twin of Abs carrying that
  derivative, and the Heaviside factors of a max/min derivative are
  spelled if(u >= 0, 1, 0) afterwards — both forms the emitters already
  print, and neither copies the argument the way a Piecewise rewrite of the
  value would (nested abs grew exponentially that way; BIOMD0000000385's
  derivatives are three times smaller and derive 2.6× faster than under the
  rewrite). Only an abs over a differentiation variable is twinned
  (abs(k) over a parameter still prints as abs(k)), and Piecewise
  conditions are left alone. At a tie the derivative is a one-sided
  convention; that is a switching surface, where the self-check judges
  nothing since #530. ceiling, floor and sign stay pre-declined. Of the
  28, 14 now attach and match finite differences to ~1e-10 at random states,
  BIOMD0000000385 among them (it used to spend 138 s discovering a decline).
  Five 472-reaction models (BIOMD0000000469–473) now derive to the edge of
  the 20 s budget — they attach on an idle machine and not under load —
  where they used to decline instantly; raise BNGSIM_JAC_DERIV_BUDGET_S to
  keep them. The rest stop at the C++ gate or elsewhere for reasons of their
  own, which the new analytical_jacobian_status names: five sit with their
  seed exactly on a switching surface of their own if(c > 0) / if(c < 0)
  rate-law pairs (ph_lorenz_attractor is one), where the Piecewise derivative
  has a point defect the gate is right to refuse.

• Every path that declines the analytical Jacobian says why (issue #506).
  Four decline paths — a function that could not be inlined, a rate law that
  could not be parsed, a symbol that survived inlining, a derivative no
  emitter prints — and the C++ attach gate went silent at every log level, so
  prepare_analytical_jacobian() returned False with no way to learn the
  reason short of monkeypatching the differentiator. Each site now logs one
  INFO line on the bngsim logger naming the rate law and the reason (the
  surviving symbols, the offending node of the derivative, the failed parse,
  the C++ gate's verdict), the former DEBUG line for the pre-declined
  functions is at INFO with the rest, and the reason is readable as
  Model.analytical_jacobian_status.

• The analytical-Jacobian self-check no longer rejects a correct Jacobian
  because the initial state sits on a switching surface of a piecewise rate
  law (issue #511). NetworkModel::set_functional_jacobian validates the
  freshly attached closed-form Jacobian against reliability-gated central
  finite differences, first at the initial state, and drops the whole
  Jacobian on the first trustworthy mismatch. A rate law such as
  if(V > 0, V*k, 0) switches on a quantity that is exactly 0 at the seed
  state whenever the species it compares are still at their seed values, so
  the initial probe sat on the switching surface: the forward probe took
  one branch and the backward probe the other, the central difference was the
  mean of the two one-sided slopes — Richardson-converged, since both step
  sizes straddle the same surface — and the analytical value, the slope of
  the branch the point itself is on, was flagged as a mismatch. Five corpus
  models (a wave equation, a Schrödinger discretization,
  BIOMD0000000075/161/613) ran their whole integration on the
  finite-difference Jacobian that way, visibly only under
  BNGSIM_JAC_DEBUG=1. The per-entry verdict now also requires the two
  one-sided differences at the smaller step to agree — they come free from
  the unperturbed RHS the gate already evaluates — and treats an entry where
  they do not as a point with no derivative to check, leaving it to the
  spread-out probe states, which do not sit on the seed's surfaces. A
  Jacobian that is wrong on either branch is still caught there. The eleven
  corpus models whose analytical entry is genuinely non-finite at the seed
  (∂(k·sqrt(A))/∂A at A = 0) stay on finite differences, as before. The
  dense self-check now prints the first mismatching or non-finite entry under
  BNGSIM_JAC_DEBUG=1, as the sparse one already did, and an attach the C++
  gate declines is logged at INFO instead of silently (the self-check path
  of issue #506).

• The _eval_rhs / _dense_analytical_jacobian test hooks refuse a state
  of the wrong length (issue #508). They copied whatever list they were
  given and read n_species entries from it, so a shorter vector — a
  Result.species row on a model with event-promoted state — was read past
  its end: undefined behaviour that returned NaN on one call and a number on
  the next. Both now raise ValueError naming the two lengths, as the public
  evaluators of #527 already did.

• parity_checks: the two Jacobian characterization harnesses run on the
  parity gate's own settings, from the full integrator state, tolerant of a
  degenerate sample, and from the Jacobian the model actually carries (issues
  #508, #509, #510, #512). Four defects in
  bng_parity/jacobian_characterization.py and
  rr_parity/jacobian_characterization_sbml.py, all found re-measuring the
  paper's Figure 2. States (#508): the SBML harness fed Result.species
  rows to a Jacobian hook sized by the full core state, reading past the
  buffer on the 219 corpus models with an event-promoted parameter or
  compartment; both harnesses now sample Result.state. Settings (#509):
  the BNGL harness silently integrated every model on t_end = 100 when the
  gate's runs/report_ode.json was absent — the state of every fresh
  checkout — and the SBML harness ran at the SED-ML tolerances the gate
  overrides. The BNGL run now refuses to start without the gate report (or an
  explicit --horizons), a model the gate has no row for gets its density
  only, and the SBML harness applies the gate's rule — the shared
  DEFAULT_RTOL / DEFAULT_ATOL unless the job carries a per-model tol
  override, integrating from the SED-ML initial_time as the gate does. Every
  row records the run block it was integrated on (horizon, tolerances, their
  sources) and the gate's outcome for the model. Samples (#510): one
  output-grid state with a non-finite RHS or Jacobian — an interpolant
  slightly negative where a rate law takes a root or a log — raised
  LinAlgError out of the eigenvalue sweep and voided the model's stiffness
  (seven models that pass parity). A negative entry smaller than atol is
  clamped to zero before the evaluation, a sample still non-finite is skipped
  and counted (n_time_points_skipped, n_time_points_clamped), and only a
  sweep with no surviving sample fails the row. Source (#512): both
  harnesses wrote analytical_jacobian_complete and then characterized from
  the private hook's partial assembly regardless, under jacobian_method: "native_analytical" — a Lorenz attractor came out at density 0 and was
  filed as degenerate. They now evaluate through Model.jacobian (the
  closed form when complete, the difference quotient of Model.rhs
  otherwise) and record which as jacobian_method
  (native_analytical / finite_difference); --analyze files a row from
  an older report that carries the partial assembly as
  incomplete_jacobian, excluded from every census, and no longer calls a
  row without a stiffness sweep degenerate.

• sbml_to_bngl keeps a parameter-valued species initial condition
  symbolic (#521). The .bngl writer had the defect #517 fixed in
  write_net: every species' initial value was written as a number, so a
  species whose initial condition is a parameter — an SBML
  initialAssignment naming one, or a compound parameter-only initial
  condition lowered to a derived _ic_* parameter — lost that reference, and
  a forward sensitivity with respect to that parameter came back identically
  zero from Model.from_bngl — finite, no warning — while Model.from_sbml
  and, since #517, Model.from_net gave the right answer. The species line
  now names the parameter; BNG2.pl carries the name into its .net and the
  reader records the reference. An amount-valued species in a compartment
  whose volume is not 1, which the flat .net has to leave as a literal, is
  written as name / V here: BNG2.pl lifts the expression into a derived
  _InitialConc* parameter and the reload reaches the model parameter
  through it.

• A species assignment rule survives the .net round trip, and the default
  gate can tell when it does not (#515). Model.from_sbml reports an
  assignment-rule species at its rule's live value through a report map the
  simulator applies at output time; the flat .net never carried that map, so
  the reloaded model reported the species frozen at its initial value while its
  dynamics were already right — every rate law bound the rule's observable or
  function, not the clamped species. Model.from_net now rebuilds the map from
  the file's own structure: a fixed species that shares its name with a
  non-identity observable (a linear rule) or with a function (any other rule)
  reports through it, exactly as the source model does. The sbml_to_net
  L2 gate compares the two models' reported values at the probe states
  instead of asking only whether the source rule varies, and probes a
  non-uniform state as well: the old uniform perturbation left a difference
  rule whose summands start equal (ppERKc = x1 - x2) at zero on both sides,
  which is how BIOMD0000000251 and BIOMD0000000269 were certified faithful
  while their reloaded .net reported a species that never moved. A varying
  assignment-rule species is therefore no longer refused under strict (GH #18):
  it round-trips.

• sbml_to_net keeps a parameter-valued species initial condition
  symbolic (#514). The sibling of #496 at the species block: write_net
  wrote every species' initial value as a number, so a species whose initial
  condition is a parameter — an SBML initialAssignment naming one, or a
  compound parameter-only initial condition lowered to a derived _ic_*
  parameter — lost that reference on the round trip. The trajectories of the
  two loaders were identical, but the reloaded model no longer seeded
  ∂y(0)/∂p, so a forward sensitivity with respect to that parameter came
  back identically zero from Model.from_net while Model.from_sbml gave the
  right answer. The species line now names the parameter, which the reader
  has accepted all along (issue #79); an amount-valued species in a
  compartment whose volume is not 1 keeps its literal, since the flat .net
  carries no volume to divide by and the capability report already flags it.

• sbml_to_bngl keeps a synthesized rate-law parameter symbolic (#513).
  The .bngl writer had the defect #497 fixed in write_net: every parameter
  was emitted as its evaluated literal, so a synthesized _rateLaw_* lost its
  reference to the model parameters it was built from, and a forward
  sensitivity with respect to one of them came back identically zero from
  Model.from_bngl — finite, no warning — while Model.from_sbml and, since
  #497, Model.from_net gave the right answer. A parameter that is not
  constant is now written as its expression, through the same ExprTk→BNGL
  normalization the functions block uses; everything else stays a literal.
  BNG2.pl carries the expression through to its .net, so the reload keeps
  the reference.

• sbml_to_net keeps a synthesized rate-law parameter symbolic (#496).
  write_net emitted every parameter as its evaluated literal, on the stated
  assumption that a derived parameter is a constant expression. A synthesized
  _rateLaw_* is not: it refers to the model parameters it was built from, and
  folding it to a number dropped that reference. The parameter stayed declared,
  so nothing downstream objected, but the ODE right-hand side no longer
  mentioned the model parameter and a forward sensitivity with respect to it
  came back identically zero — finite, the right shape, no warning — while the
  same model through Model.from_sbml gave the right answer. A parameter that
  is not constant is now written as its expression and everything else stays a
  literal; codegen_data lists parameters in dependency order, so anything an
  expression refers to is already declared above it.

• A .net parameter expression keeps the whitespace it was written with
  (#498). The .net reader split a parameter line on whitespace and then joined
  the value tokens back with nothing between them, so a space separating two
  word characters was deleted and the tokens welded together:
  if(k > 0.1 and thr > 0.5, k, 0.0) was read as if(k>0.1andthr>0.5,k,0.0).
  ExprTk parses 0.1andthr without complaint, so the parameter took a wrong
  value that was finite, with no warning and no error — the same rate law
  written with parentheses instead of spaces gave a different number. The
  functions block twenty lines away in the same file has always joined with a
  space; the parameters block now matches it.

• Every benchmark script now answers --help instead of running (#488,
  #489). A script that parses no argv answers "what does this do?" by doing
  it — no usage text, no refusal for a mistyped flag, and a results file
  overwritten by the command typed to avoid exactly that. #488 fixed the two
  cases the orchestrator can reach; #489 swept the rest of
  benchmarks/suites/, where the hazard turned out to be demonstrated rather
  than argued:

  script	what --help used to do
  psa/run_bngsim_timing.py	started a ~20-minute PSA sweep onto a fixed path
  nf/run.py	started the NFsim correctness sweep
  jacobian/probe_attach.py	probed a model named --help, overwriting the committed results/attach_probe.json with the one-row failure
  ode_engines_s4_sbml/check_sbml_engine_agreement.py	ran three engines, then overwrote its committed report
  ode_fullnet/recover_s4_points.py	regenerated two networks through BNG2.pl and copied the result into another repository's latex/generated/
  ssa_table5/convert_all.py, emit_ssa_table.py, merge_jobout.py	reconverted the corpus / rewrote the Table 5 result set
  jacobian/diagnose_divergence.py	read --help as a model id

  All of them now build an ArgumentParser, so --help exits 0 with usage and
  an unrecognized flag exits 2. Where a probe needs somewhere harmless to write,
  it now has a flag: --out on psa/run_bngsim_timing.py (alongside
  --warmup / --runs / --effort, the names psa/run.py already uses, which
  turn the sweep into a one-second smoke test and are recorded in its payload),
  --out on probe_attach.py and check_sbml_engine_agreement.py, and
  --paper-dir / --no-copy on recover_s4_points.py — which also skips
  rather than crashes when the paper checkout is absent.

  Two scripts are exempt by decision: ssa_table5/_ssa_cell.py and the SBML
  test suite's bngsim_wrapper.py are launched by a caller on a fixed argv, so
  a parser there would change a measurement path's contract to buy a usage line
  for an invocation nobody makes. Each now names its caller instead.

  python/tests/test_benchmark_runner_help.py pins the contract by structure —
  every suites/ script with a __main__ guard, minus that exemption list — so
  a script added to any suite inherits it without an edit. No measured value
  changes: every runner behaves identically when invoked as the orchestrator
  invokes it.

• Eight benchmark scripts imported only on a machine with the checkout at
  ~/Code/bngsim (#489). ode_engines_s3/run_s3_timing.py,
  ode_engines_s4_sbml/{run_s4_timing,check_sbml_engine_agreement}.py and the
  five ode_fullnet/ scripts resolved BNGSIM_ROOT to Path.home() / "Code" / "bngsim" and inserted <root>/parity_checks on sys.path, so the
  module-scope import _bng_common raised ModuleNotFoundError in any clone
  living somewhere else — a fresh checkout, a worktree, CI. The default is now
  the checkout the file itself lives in; the BNGSIM_ROOT override is unchanged
  and still covers the case it was written for, a different venv running against
  a canonical checkout. On a machine where the two already coincided, nothing
  moves. Every --help probe in the test above now runs with $HOME pointed at
  an empty directory, which is what surfaced this.

• gen_networks.py did not create the .net cache it exists to fill (#492).
  Phase 1 of the full-network ODE benchmark writes every generated network into
  benchmarks/suites/ode_fullnet/nets/ and never created that directory, so on
  a fresh checkout — or one where the cache had been cleaned to reclaim disk —
  every model failed. Expensively: generation happens first and only the write
  fails, so six models spent 28 s to produce nothing and the full corpus would
  spend hours. The script whose job is to build the cache was the one script
  that could not create it, which left a missing cache with no recovery path.
  The pass now creates the directory before it reports it, so the cache: line
  names somewhere that exists. A write failure that survives that — the
  directory removed while a multi-hour sweep is in flight — is recorded as
  cache_write_failed naming the directory, rather than as a bare
  FileNotFoundError on a .net path, which was indistinguishable from a
  genuine BNG2.pl netgen failure.

• A scoped benchmark run could not avoid the committed report it replaced
  (#493). Seven scripts under benchmarks/suites/ accept a flag that narrows
  what they measure — --engines, --models, --limit, --effort, --quick,
  --stride, --max-models — and wrote a tracked artifact at a path fixed in
  the source, with no --out and no environment override. All seven now take
  --out, defaulting to today's path so an unflagged run is unchanged; where a
  script resumes from its own report, --out moves the resume read too, so a
  probe does not inherit the committed run's rows. run_forced.py's default
  derives from --mode, so it resolves through one helper rather than an
  argparse default. Each report also records the filters it was scoped by
  (models_filter, limit, max_models, quick), which several did not.

  The sharp case is --engines: re-timing a model records a row holding only
  the engines that run measured, so --redo --engines bngsim replaced six
  models' cross-engine rows with one-engine ones — 513 lines of committed
  measurement, written atomically, so it reads as a finished run. Without
  --redo the rows survive but _meta.engines was still rewritten to the
  request, leaving a report that claimed one engine while holding six. The Table
  S3/S4 reports now also carry engines_present, derived from the rows
  themselves; engines keeps its existing meaning for anything already reading
  it.