A gradient column assembled at exactly 2x its central difference, once, and has not reproduced (Raia, IC axis)

[wshlavacek]

Summary

While sweeping every gntr slug's assembled gradient against finite differences for #535, one
column on Raia_CancerResearch2011 came back at exactly twice its central difference. It has
not recurred in five subsequent runs at that same point, two of which had byte-identical routing.

I am filing it because a factor of exactly two on a sensitivity column is not the shape numerical
noise takes, and because a non-reproducible wrong gradient is strictly worse than a reproducible
one: the #535 sweep only caught it by chance, and nothing in the product would have.

The observation

Slug Raia_CancerResearch2011 (Grein-2026 subset-I corpus), evaluated at its PEtab nominal point
displaced 2% of each box width along sampling space. ‖∇‖∞ = 2.274e+04, so this is a
well-conditioned point, and the central difference is stable to five digits across three step sizes:

param            assembled   fd(h=1e-3)  fd(h=3e-4)  fd(h=1e-4)
init_Rec_i         -45475       -22738      -22738      -22739

-45475 / -22737.6 = 2.00000. Every other one of Raia's 39 columns agreed to 1e-4 or better on that
same run.

Probing that column directly at the same point, the FD side is sound and the perturbation is
applied exactly (no bound clamping, no silent reflection):

theta(u0-h) = 162.50845020924407     f(-h) = 2683.0361725734565
theta(u0)   = 162.62074585080055     f(0)  = 2676.2079742098376
theta(u0+h) = 162.73311909023639     f(+h) = 2669.393618322428

forward  diff: -22714.5
backward diff: -22760.7
central  diff: -22737.6

Forward, backward and central all agree, so the true derivative is ≈ -22738 and the assembled
value was the wrong one.

Why it is interesting

init_Rec_i is Raia's only initial-condition-axis free parameter. Its route is a single
contribution and there is nothing in it to double:

seed terms: (SeedTerm(target='ic', key='Rec_i', node=('num', 1.0)),)
route:      [('ic', 'Rec_i', 1.0)]

Raia also puts prediction_formula σ on all eight observables (ADR-0079/0080), so the assembled
column runs through the sigma_sensitivity chain rule — the one place where a term could plausibly
be added twice. assembled - truth = -22737, i.e. exactly the true value again, which is consistent
with one contribution being counted twice rather than with a scale factor being wrong. That is a
hypothesis, not a finding; I could not make it happen again to test it.

What I ruled out

• Not the Verify the assembled gradient against finite differences on every real gntr problem, not just fixtures #535 fix. The anomalous run predates ADR-0097, but Raia's routing is byte-identical
  either side of it — verified by monkeypatching the old rule and diffing sensitivity_params,
  sensitivity_ic, and every per-parameter route. Nothing changed, so that fix neither caused nor
  cured this.
• Not a box-bound artefact. set_value clamps out-of-box values, which makes a central
  difference a half-step and was behind four other red rows in the Verify the assembled gradient against finite differences on every real gntr problem, not just fixtures #535 sweep. init_Rec_i sat at
  u = 2.211 in a [-5, 3] box, nowhere near either bound, and the probe above confirms both
  perturbations landed exactly where intended.
• Not FD noise. Stable to five digits across h ∈ {1e-3, 3e-4, 1e-4}, and forward/backward/
  central agree.
• Not the evaluation point. The same point at --disp 0.01 and 0.04 agrees to 8e-06 and 1e-06.

Reproduction attempts

Five runs at the identical point since: one in the post-fix sweep, three explicit back-to-back
repeats, and one in the final bounds-clear sweep. All five read -22737 against -22738. The corpus
sweep for #535 now reports Raia clean.

Suspects

No confirmed cause. Recording the one lead: lanl/bngsim#117 established that a module-level printer
in that stack was not thread-safe, and it is what was ultimately behind #525's per-dose sensitivity
mis-stacking. It is closed, and its symptom was a derivative reported as "not representable in
C" rather than a silently scaled sensitivity column, so I want to be clear this is a precedent for
"non-determinism has reached the sensitivity path before", not evidence that it happened here.

Suggested next step

Nothing to fix until it reproduces. The cheap standing guard would be an assertion, or a debug-level
check, that a single-contribution ic-axis route contributes exactly once per scored point — the
narrow invariant this violated. Reproducing it may need the run repeated under load, since every
attempt above was on an otherwise idle machine.

Found during #535; the corpus-side tooling is tools/fd_check.py in wshlavacek/BNGL-Models.

[wshlavacek]

Landed the standing guard this issue asked for, in b598e2d. Not closing — the guard covers the
routing/assembly half of the hypothesis, and the half that fits the evidence best is still open. Details
below, including one hypothesis narrowed and one negative result worth recording.

What went in

Two halves, because the invariant could not be checked as stated.

The obstacle: a repeated (axis, key) column is sometimes legitimate. Two condition targets whose
initialAssignments seed one species' initial value both land on (ic, that species), and their sum
is the derivative. So "this route reads one column twice" could not be treated as a defect — which is
precisely why the failure was undiagnosable.

1. route_experiment now folds same-column terms into a single contribution carrying their summed
   derivative tree (_fold_same_column). Arithmetically neutral, one tensor read instead of two, and
   still symbolic — so a point-dependent path (Extend ADR-0076 condition routing to a non-unit seed derivative (d(IC)/d(target) != 1) -- Bertozzi_PNAS2020 #530) keeps being refreshed by at_point across the
   whole sum. One-contribution-per-column becomes structural rather than incidental.

2. ExperimentRouting.check_column_multiplicity() then asserts it where the columns are actually
   summed: once per experiment from both _raw_sensitivity_accessor and
   _constraint_sensitivity_accessor, so the gradient, EFIM and constraint assemblers are all covered.
   It raises a PybnfError naming the free parameter, the axis and key of the repeated column, and each
   duplicate's factor.

Cost is O(n_params) per experiment per evaluation, so it is always on rather than debug-gated. Nine
tests; the real-bngsim FD acceptance gates are unchanged.

Why there is no fix

Nothing localized to fix. Every deterministic path that could double an ic column is already correct:

• classify_bound_id already guards the one known self-double — the SBML backend reports each species as
  its own initializer ([(s, s)]), so without the name not in species_names guard a species-named free
  parameter would collect both the seed term and its own IC term
• _compute_ic_seed_map / _net_species_ic_seed_map emit one SeedTerm per assignment symbol over
  sorted(inputs)
• tensor_sens assigns vec[index[name]] = total rather than accumulating, so two route names cannot
  collide into one column
• enable_output_sensitivities replaces the request rather than appending; _request_keys and
  apply_routings de-duplicate

That matches the byte-identical-routing finding above. A "fix" here would mean editing correct code
against an unlocalized 2x.

Narrowing: the σ chain rule is not it

The sigma_sensitivity hypothesis above ("the one place where a term could plausibly be added twice")
does not fit the observation on reflection. assembled - truth = truth means the whole column
doubled — both the residual-Jacobian term and the noise term. PredictionFormulaSigma.sigma_sensitivity
contributes only the latter, so double-counting there could not produce an exact 2.000000 on the total.
It also could not be confined to one column: the prediction-coupling branch chains through raw_sens,
which perturbs every model-parameter column, and the other 38 agreed to 1e-4.

What does fit exactly is raw_sens[..., init_Rec_i] being doubled at every point — i.e. the d_ic
tensor itself. The sensitivity system is linear in S, so a seed of 2 rather than 1 doubles the
whole column at every time and every observable, giving exactly 2.000000 on the assembled total, on the
IC axis alone, with every other column untouched. That is the shape.

Negative result: the seed is solid on a small model

At t=0, dx_j/dx_k(0) must be exactly delta_jk, so a doubled seed is directly readable off row 0.
Probed against bngsim 0.12.1 on e2e_ode_decay.net, serially and under thread load (#117 being the
precedent for concurrency reaching this stack):

serial      : n=200  seed min/max=1/1  distinct seeds=1  max|col-truth|=5.983e-11  non-unit seeds=0
8 threads   : n=400  seed min/max=1/1  distinct seeds=1  max|col-truth|=5.983e-11  non-unit seeds=0

600 runs, every seed exactly 1.0, every column matching the analytic e^{-kt} to 6e-11. This does not
clear the hypothesis — a one-species decay model is a weak test against 39 parameters, 8 observables and
a prediction_formula sigma on every column — but it does mean there is currently nothing to hand
bngsim.

What is still open, and what would close it

A doubled column arriving from the backend's own d_ic tensor still passes this guard. The only exact
PyBNF-side check for it is the t0 unit seed above, and it does not generalize: it needs
species-labelled selectors (SBML only — the net backend's tensor is over observables) and a row that
really is t0, which pre-equilibration breaks. So it was not added.

No bngsim issue reports a doubled IC seed. The nearest are the same seeding machinery with different
symptoms — lanl/bngsim#43 (missing seed), #111 (silently zeroed seed), #113 (stale seed: a column that
should be 0 reported as e^{-kf t}) — plus lanl/bngsim#51, where a warm codegen cache made #41/#43
silently inert. That last one is the interesting shape for a once-in-six-runs anomaly: cache-state-
dependent sensitivity seeding, rather than a race. Worth filing upstream if this recurs with a
bngsim-side fingerprint, or if a Raia-scale version of the probe above reproduces it; filing now would be
a hypothesis with no reproducer.

Keeping this open for that half.

[wshlavacek]

Recording the search, since it came up empty and the next person should not repeat it. No bug found
— not in PyBNF, not demonstrably in bngsim. b598e2d adds a guard on the one invariant that could be
checked; it is not a fix, and nothing here explains the observation. Deliberately not filing at
lanl/bngsim: a hypothesis with no reproducer is not something those devs can act on.

The localization argument

Worth stating explicitly, because it constrains the search hard. Only one column moved — the other
38 agreed to 1e-4 on the same run. So whatever happened is confined to the per-axis read for
init_Rec_i, which is:

total += c.factor * _sensitivity(sens, selector, c, row, name)   # one contribution, factor 1.0

Two degrees of freedom: c.factor, and the tensor value. c.factor comes from a route verified
byte-identical across runs, whose sole contribution is a constant 1.0 with no symbolic node — so
at_point returns the routing unchanged and there is no per-point arithmetic to get wrong. That leaves
the tensor value. Every remaining hypothesis has to explain a doubled d_ic[:, :, Rec_i], and any
hypothesis that would perturb more than one column is already excluded by the other 38.

Ruled out (PyBNF side)

Each of these was traced to the code, not reasoned about abstractly:

• Route construction. classify_bound_id already guards the one known self-double — the SBML
  backend reports each species as its own initializer ([(s, s)]), so without its
  name not in species_names guard a species-named free parameter would collect the seed term and
  its own IC term. Both seed-map builders (_compute_ic_seed_map, _net_species_ic_seed_map) emit one
  SeedTerm per assignment symbol over sorted(inputs).
• Per-parameter accumulation. tensor_sens assigns (vec[index[name]] = total), so two route
  names cannot collide into one column. _iter_scored_points yields each (row, col) exactly once.
• Request plumbing. enable_output_sensitivities replaces the request rather than appending;
  _request_keys and apply_routings both de-duplicate; _setup_gradient_path is idempotent.
• Tensor addressing. _sensitivity addresses by label (labels.index(key)), so neither a permuted
  nor a duplicated sensitivity_ic_species echo can double a read on the PyBNF side — .index() takes
  the first match.
• The whole experiment scored twice. Would double every column, not one. Excluded by the 38.
• The σ chain rule (this issue's own lead). Excluded by shape — see the previous comment.
• The t0 unit seed on a small model. 600 runs, serial and 8-thread, seed exactly 1.0. See the
  previous comment.

Where I would look next, in order

1. Run the t0 seed check at Raia scale. dx_j/dx_k(0) must be exactly delta_jk, and on the SBML
   backend the tensor's selectors are species, so the seed is directly readable off row 0 —
   d_ic[0, j, k] == unit_factor * delta_jk. The decay-fixture probe was clean but is a one-species
   model; Raia is 39 parameters, 8 observables and a prediction_formula sigma on every column. This is
   the cheapest test that could actually fire, and it needs the corpus rather than this repo. Note it
   does not generalize into PyBNF as a standing check: it needs species-labelled selectors (SBML
   only — the net backend's tensor is over observables) and a row that really is t0, which
   pre-equilibration breaks. That is why b598e2d does not include it.

2. Vary cache state, not just load. All six attempts were back-to-back at the same point on an idle
   machine — i.e. always the same codegen-cache state. The anomalous run was inside a corpus sweep,
   where many other models had passed through that cache first. Codegen cache does not invalidate on sensitivity-RHS changes (#41 and #43 were silently inert on warm caches) bngsim#51 is precedent that warm-
   cache state silently changes sensitivity seeding behaviour (it made Differential evolution typo #41/Mutate parameters #43 inert). "Run Raia after a
   sweep of other slugs in the same process" is a different variable from "run it under load", and
   nobody has varied it. I would try this before trying harder on concurrency.

3. Establish whether the anomalous evaluation was a fresh worker. PyBNF scatters models to workers
   and the sensitivity request rides the pickle. Worth confirming that a reused worker cannot end up
   applying a request twice, or holding a model whose request was enabled under a different routing.
   enable_output_sensitivities replaces rather than appends, so the single-call path is clean — the
   open question is whether the scatter can deliver two live model objects for one (model, suffix).
   This one is answerable in this repo.

4. Make a recurrence diagnosable in one shot. The strongest practical change is out-of-tree: have
   tools/fd_check.py dump the raw d_ic slice, the resolved routing, and sensitivity_ic_species
   whenever a column disagrees with its finite difference beyond tolerance. This anomaly cost a
   re-run-that-never-reproduced precisely because the artifact was gone by the time anyone looked. With
   the dump, one recurrence settles hypothesis 1 immediately.

What is now automatic

b598e2d means the routing half can no longer produce this silently: any route reaching one native
column twice raises by name, on every gradient, EFIM and constraint assembly. If it recurs without
that error firing, hypothesis 1 is the answer and the tensor is where to look.

[wshlavacek]

Two branches the previous comments left out. Neither explains the observation, but one of them was an
entire unexamined category and the other bears on whether hunting further is worth it.

1. The measuring instrument

Everything above verified the thing being measured. Nothing had verified the thing doing the
measuring — a wrong tools/fd_check.py would produce exactly this signature: one bad column, no
reproduction, and nothing wrong in the product. So I read all 172 lines of it.

No mechanism there doubles a single column. Specifically ruled out:

• FreeParameter.set_value returns a copy, so the FD loop cannot alias or mutate the free list
  the assembled gradient is later built from — the two sides genuinely evaluate at the same u0.
• The report's zip(g.param_names, g.gradient, grad_fd) is order-aligned: param_names and grad_fd
  both derive from the same free ordering.
• LikelihoodObjective.evaluate_multiple resets _pset_values = {} on entry, so the ~78 FD
  evaluations that run before the assembly cannot leave stale bindings behind for
  PredictionFormulaSigma to misclassify as coefficients.
• The routing construction mirrors _setup_gradient_path faithfully, including _condition_for_suffix
  and the per-experiment at_point.

But it did surface two real divergences from gradient_at. Neither can be Raia's — checking the
conf, Raia declares no normalization and no constraints — but both are live ways for this instrument
to lie on other slugs, and belong in tools/README.md alongside the five gotchas already there:

• The harness omits data_key. It appends 3-tuples, (sim, exp, routing), where gradient_at
  appends (sim, exp, routing, suffix). On a slug with an analytic per-series scale (ADR-0066, Track the deferred gradient for analytic per-series scaling (ADR-0066/ADR-0475): no gradient job_type can fit a scaled column #533)
  that has two consequences, and the second is the dangerous one. The scaled column is refused —
  loud, fine. But _iter_scored_points only repoints objective._scale_factors at the current
  experiment when data_key is supplied; without it, _scale_factors keeps whatever the last
  evaluate() inside evaluate_multiple left there, so every experiment's residual is scored with
  the last experiment's profiled c*. That is silent, and the ratio of two experiments' c* is
  exactly the kind of clean number that can present as a small integer factor.
• Constraints are in the FD side only. loss_at calls
  evaluate_multiple(..., config.constraints), which adds cset.total_penalty(...); the assembled
  side calls assemble_gaussian_gradient alone and never assemble_constraint_gradient. On a
  constrained slug every constraint-touched column reads red as an artefact.

So the instrument moves from unexamined to examined, no single-column doubling mechanism, two
unrelated latent defects found. It does not explain Raia.

2. Is it even a software defect? The exponent check

Worth doing before spending more on the hunt, because "doubled" has a hardware reading: a soft error
that increments an IEEE-754 exponent field. If that were the story, there is nothing to find.

It mostly does not hold here. The true value is ≈ -22737.6, which lies in [2^14, 2^15), so its
biased exponent is 1037 = 0b100_0000_1101. Doubling means 1038 = 0b100_0000_1110 — the two differ
in two bits, not one. A single-bit flip in this value's exponent LSB would have halved it
(1036, exponent 13), not doubled it. So the tidy one-bit story does not apply, which pushes weight
back onto a real defect.

Two related notes:

• A bit flip in a single tensor entry could not produce this either: the assembled column is a sum
  over scored points, so corrupting one entry perturbs the total by a fraction, not by exactly 2. Only
  a doubled seed (which scales every point identically, the sensitivity system being linear in S)
  or a corrupted final scalar gives a clean 2.
• Those two remain distinguishable in principle, and the distinguishing test was not run: the report
  gives -45475 / -22737.6 = 2.00000 to six digits, which cannot separate "exactly 2.0 bitwise" from
  "2.0000x". A doubled seed gives 2.0 to all 16 digits; anything accumulating through the sum
  generally does not. Recommendation 4 in the previous comment — dumping the raw tensor and routing on
  any FD mismatch — would settle that on the first recurrence, and printing the ratio at full
  precision costs nothing.

Net: after examining the instrument and discounting the hardware reading, the doubled-d_ic hypothesis
is still the one standing, and hypothesis 1 (the t0 unit-seed check at Raia scale) is still the test
that would confirm or kill it.

[wshlavacek]

Instrument audit finished and fixed. Three defects in tools/fd_check.py, none of them this anomaly,
plus two more non-reproductions. Corpus side is wshlavacek/BNGL-Models 2e3bfd9 and 2ac53da; the PyBNF
side is a docstring, 3e42e67.

The two divergences, fixed

Neither can be Raia's: it declares no normalization and no constraints. In fact no slug in this
corpus declares either, so both fixes are no-ops across the corpus today and are forward protection
only. Both are now in tools/README.md beside the five gotchas already there, because both produce a
plausible wrong answer rather than an error.

• data_key was omitted. The script appended (sim, exp, routing) where gradient_at appends
  (sim, exp, routing, suffix). The loud consequence is mild — a column under an analytic per-series
  scale is refused. The silent one is not: _iter_scored_points repoints
  objective._scale_factors at the current experiment only when data_key is supplied, so without
  it every experiment's residual is scored with whatever profiled c* the last evaluate() inside
  evaluate_multiple left behind. A ratio of two experiments' c* presents as a clean multiplicative
  factor — the same shape as this issue, on a slug that normalizes.
• The constraint penalty was on one side only. loss_at passes config.constraints to
  evaluate_multiple, which adds cset.total_penalty(...), so the finite difference has always
  differentiated the penalty; the assembled side never called assemble_constraint_gradient. Every
  constraint-touched column would read red by construction.

Verified as an exact no-op where it must be: same point, same step, pre- and post-fix on Raia agree
on every assembled value and every central difference to all printed digits.

A third defect: the default-point path never ran

With no param-values.json the script died before simulating anything, with a traceback that points at
the wrong place:

AttributeError: 'float' object has no attribute 'name'
  pset.py:2401 in set_value  ->  if new_value < self.lower_bound
  pset.py:2660 in __lt__     ->  return self.name < other.name

FreeParameter.initial_value_from_quantile returns a FreeParameter, not a value. The name reads
otherwise, but that is the contract every caller in PyBNF depends on — algorithms/base.py builds its
latin-hypercube PSet straight from the return value. The script wrapped it in a second set_value, so
a FreeParameter arrived where a float belongs and new_value < self.lower_bound resolved to
FreeParameter.__lt__(float), which reaches for .name on the float. PyBNF's behaviour is correct;
only its docstrings said "value", so 3e42e67 documents the return type and names this symptom so the
next person searching that traceback lands on it.

Two more non-reproductions

Both at --disp 0.02, h = 3e-4, all 39 columns:

point	init_Rec_i assembled	central diff	rel err	worst column
SBML parameter values	-93181.4	-93182.4	1.1e-05	1.05e-03
median-quantile default	-3.09057e+06	-3.09057e+06	1.1e-07	1.73e-04

That is now eight evaluations of this column at four distinct points with no anomaly.

Incidental, and relevant to how the non-reproduction was established

Two independent invocations of the identical command at the identical point do not agree bitwise.
Assembled and central-difference values print the same to six figures, but the relative-error column
moves in its last digits on five rows — a ~1e-9 relative run-to-run difference, consistent with an
adaptive solver's step sequence varying between runs.

Not remotely a factor of two, and it does not explain anything here. But it does mean byte-identical
routing does not imply byte-identical numerics on this stack, so "five subsequent runs read -22737"
constrains the anomaly less tightly than it appears to: those runs agreed to printed precision, not to
the last bit. Worth knowing before treating any future non-reproduction as decisive.

Still open, still on hypothesis 1 (the t0 unit-seed check at Raia scale).

[wshlavacek]

Found the mechanism, and reproduced this issue's exact signature on demand. d660ba9.

init_Rec_i's parameter axis is not zero — it is the whole derivative

ADR-0097 drops the sensitivity_params axis of a parameter that only seeds species initial values,
and documents the reason as economy: that axis "would be identically zero", so requesting it merely
wastes a forward-sensitivity vector. That premise is false. Requesting both axes at once on Raia:

max|d_param[init_Rec_i] - d_ic[Rec_i]| / max|d_ic[Rec_i]|  ==  0.000e+00      (all 20 outputs)

Bit-for-bit identical. bngsim seeds ∂x(0)/∂θ into the parameter axis as well (lanl/bngsim#43), so
for a parameter the right-hand side never otherwise reads, that axis carries the entire seeding
derivative — the same numbers the initial-condition axis carries. Injecting the dropped axis back into
the route:

param                              assembled    central diff    rel err
init_Rec_i                      -6.18114e+06    -3.09057e+06   1.00e+00  <--
worst relative error: 1.00e+00

2.00000x, on that column, every other column untouched. This issue's signature, deterministically,
from a routing that reaches one derivative on two axes.

What that does and does not settle

The mechanism is settled: the 2x is one seeding derivative counted on both axes. The trigger on the
original run is not. Both the pre-#535 rule (not _seeds_only_initial_conditions(seeds)) and the
post-#535 rule route exactly one axis for Raia, and the degenerate states do not help — an empty
ic_seed_map yields a lone (param, init_Rec_i) term, which by the measurement above is also the
correct value. So I cannot yet produce the state that would have routed both.

Three fixes

1. The claim is corrected wherever it appeared — the routing module docstring,
   classify_free_param, classify_bound_id, and both backends' ode_rhs_symbols. The drop is
   load-bearing for correctness, not an optimization.
2. The unanswerable-model fallback is now a refusal. ADR-0097 kept the axis when
   ode_rhs_symbols() returned None, on the strength of "worst case, one redundant vector". With
   that false, the two errors are deleting half a derivative (Verify the assembled gradient against finite differences on every real gntr problem, not just fixtures #535) and doubling it (this issue), and
   there is no safe default left — so the router names the parameter and points at a gradient-free
   job_type. Both shipped backends always answer, so no fit changes; it closes the hole.
3. The assembly checks numerically, once per experiment: a route holding both its own parameter
   axis and an IC axis whose tensor slices are bit-identical is refused by name. Exact equality is the
   discriminator on purpose — two independent derivatives of a live model do not agree to the last bit.

Verified: the injected-both-axes Raia run now refuses instead of reporting 2x; plain Raia and
Fiedler_BMCSystBiol2016 are unchanged (Fiedler legitimately routes both axes for seven parameters,
whose columns genuinely differ, and still agrees at ≤ 3.7e-06); gradient suite 383 passed / 4 skipped.

The guard from b598e2d would not have caught this

Worth stating plainly. Those two contributions sit on different columns — (ic, Rec_i) and
(param, init_Rec_i) — so check_column_multiplicity passes them. It guards a shape that cannot
currently occur while the shape that reproduces this issue went straight through. Both guards now stand,
but the first one was aimed at the wrong thing.

A forward regression, and a contract question — lanl/bngsim#155

Neither of bngsim's two quantities is wrong; summing them was ours. But which of them the parameter axis
carries is version- and expression-shape-dependent, and lanl/bngsim#147 (merged today, on main,
not in the installed v0.12.1) just changed it for a whole class:

shape	example	d_param[p] on v0.12.1	PyBNF routes
bare <ci>	Rec_i = init_Rec_i	RHS + seeding	ic axis only — correct
compound	Fiedler's steady-state ICs over K_1, k2…k11	RHS only	ic terms + param axis — correct

#147 lowers a compound parameter-only initialAssignment so it seeds too. Fiedler's six initial
assignments are compound expressions over ten symbols, every one a constant="true" parameter —
precisely that predicate. So on the first build carrying #147, Fiedler's ic terms will duplicate the
seeding and the seeding half of seven columns will be counted twice.

Fix 3 above will not catch that: once the parameter also drives the RHS, d_param[p] is
RHS + seeding, which is not bit-identical to any single ic column. Detecting it needs the contract
stated, which is what lanl/bngsim#155 asks for.

Keeping this open: the trigger for the original observation is still unexplained, and the #147 exposure
needs closing once bngsim answers.

[wshlavacek]

Trigger hunt. I have a verified trigger that produces this exact signature, but I cannot yet show
it was reachable on the original run. Also one correction to what I shipped.

The routing layer, enumerated exhaustively

Rather than inspect, I enumerated every reachable input state of classify_bound_id for
init_Rec_i — seed map (absent / IC seed / IC+PARAM seed / PARAM seed / non-routable) x
rhs_symbols (None / absent / present) x param_ids (member / not). 30 states; exactly 4 route
both axes, in two families:

• Family A — ic seed + init_Rec_i present in rhs_symbols.
• Family B — the seed map hands init_Rec_i a PARAM seed term alongside the IC one.

Both are unreachable for Raia, checked against the SBML rather than assumed:

• Family B needs a second initialAssignment targeting a parameter. Raia has exactly one, on the
  species Rec_i, so _compute_ic_seed_map can only ever emit (IC, Rec_i) and
  _seeds_only_initial_conditions is unconditionally true.
• Family A needs init_Rec_i in a kinetic law, rule, event, or function body — the four
  over-inclusive sources _compute_rhs_symbols walks. Raia has 16 reactions whose 60 <ci> nodes
  never name it, and zero rules, zero events, zero function definitions.

Two more paths cleared the same way:

• Condition param-refs. Raia's four conditions are numeric = perturbations of il13_level
  (il13_level = 0 / 20 / 4 / 80). None is is_param_ref, so ref_contribs is empty and the
  condition path cannot add a contribution to any route.
• Request-set interference. d_ic[Rec_i] is invariant across four parameter-axis request sets
  (nothing / init_Rec_i alone / the real 17-parameter fit request / that plus init_Rec_i):
  ratio range [1.000000, 1.000000] on all 20 outputs. Asking for more does not perturb the IC seed,
  so "two seeding paths collide in one seed matrix" is not supported on this build.

The trigger that does work: rhs_symbols is None

_absent_from_rhs returned False for None, so "the model cannot say" kept the axis. Verified
directly against aaa9d7aa, the shipped commit before the fix:

aaa9d7aa (pre-fix)  rhs known, absent  -> [('ic', 'Rec_i')]
aaa9d7aa (pre-fix)  rhs = None         -> [('ic', 'Rec_i'), ('param', 'init_Rec_i')]

and both axes on Raia is exactly 2.00000x. So the chain is complete and every link is measured.

What makes this more than hypothetical is how None is reachable. _rhs_symbols (SBML) and
_net_rhs_symbols (net) are class attributes defaulting to None, and the comment on each says
why:

A class attribute so an instance built via object.__new__ (the test fakes, unpickling) answers
"cannot say" rather than raising — and "cannot say" keeps every parameter axis, which is the safe
direction.

That is the false premise, written at the point where the default is chosen. "Cannot say" is reachable
by omission — any instance that reaches the router without __init__ having run — and getattr
succeeds silently because the class attribute is there.

Against Raia's actual path, though: _extract_sbml_structure() is called unconditionally in
__init__, so a normally-constructed model always has a real frozenset, and fd_check.py builds its
models in-process. So I can show the door exists and exactly what comes through it, but not that
anything went through it on that run.

A timeline tension worth recording

The report says the run "predates ADR-0097". It also says the point was the PEtab nominal displaced 2%
of each box width — that is --disp, which landed in corpus commit ef52ef5 at
2026-08-02 18:13:35, which is 102 seconds after ADR-0097 (2d111a38, 18:11:53). Both cannot
hold unless the displacement was applied by hand through a param-values.json. If the attribution is
off by those two minutes, Family A is live for that run and the class-attribute default is the
mechanism.

I could not settle it by re-running: 21 of Raia's 39 nominal values are observable/noise parameters
that live in upstream parameters.tsv, which is deliberately not vendored, so the original point is
not reconstructible from this checkout.

Correction: the fold undercuts the multiplicity guard

_fold_same_column (b598e2d) merges two contributions on one column into a single contribution
carrying their summed factor. That means a route that reaches (ic, Rec_i) twice now arrives at the
assembly as one contribution with factor 2.0 — and check_column_multiplicity cannot see it. I
made the guard blind to exactly the router-produced duplicates it was meant to catch; it now fires only
on hand-built routings.

The fold is still right (it is what makes the invariant structural, and two condition targets seeding
one species genuinely do meet). The fix is to have it record what it folded, and have the check read
that, rather than inferring multiplicity from the contribution count after the fact. I will do that
unless there is a reason to prefer a different shape.

The cross-axis guard added in d660ba9 is unaffected — it compares tensor slices, not contribution
counts, and it is the one that catches the demonstrated shape.

[wshlavacek]

Fixed the fold weakness reported above. 153c974.

Counting could not fix it. After a fold, a doubled factor is indistinguishable from a legitimate
single term — and same-column terms genuinely do arise, so the fold cannot simply be dropped: two
condition targets whose initialAssignments seed one species both land on that species' ic column,
and their sum is the derivative. What separates the cases is not how many terms met but which
paths they came from. One path arriving twice is the defect; two paths meeting is arithmetic.

So each contribution now carries its provenance — 'bind' for the bind-by-id classification,
'ref:<target>' per condition parameter-reference. The fold concatenates the labels,
check_column_multiplicity refuses a contribution whose labels repeat (naming the path and the folded
factor), and at_point carries them across, which matters because a point-dependent routing is rebuilt
every optimizer step and would otherwise go blind exactly where it is re-resolved.

Why a repeated label is a defect by construction, not a heuristic: within one classify_bound_id call
the columns are already distinct — seeds come one per initialAssignment symbol, and the id's own axis
cannot collide with them — so a repeated 'bind' means the bind path was applied twice. Two
param-refs through the same target are likewise a redundant re-application rather than a second path.

origins is metadata, not identity: compare=False, repr=False, so equality, hashing and the
dataclass repr are untouched and every existing assertion about what a routing computes still holds.

Gradient suite 432 passed / 4 skipped; Fiedler_BMCSystBiol2016 unchanged end to end at 3.68e-06.

Where this leaves #537

Three guards now stand, and they cover different things:

guard	catches	would it have caught the original?
check_column_multiplicity + provenance	one chain-rule path reaching one native column twice	no — the reported route was single-contribution
_check_axes_are_not_the_same_derivative	ic axis and own param axis holding the same numbers	yes — this is the demonstrated 2.00000x shape
classify_bound_id refusal on rhs_symbols is None	the one reachable state that routes both axes	yes — closes the door rather than detecting it

Still open for the reason it has been: the mechanism is settled and reproducible, the trigger for that
one run is not, and the lanl/bngsim#147 exposure needs closing once lanl/bngsim#155 answers.