fix(core+engine): three correctness defects — cross-process network determinism, observable count filters, batch-SSA seed aliasing; one report-only

[akutuva21]

Three independently-verified correctness fixes and one report-only item, as three commits that can be read and merged one at a time. Base 6889fba.

Merge hazard, stated first because it is not obvious. Commit 3 changes batch-SSA output values by design. A byte-identity or md5 comparison of a batch-SSA artifact across this PR will show a diff that is not a regression. cpp/engine/OdeIntegrator.cpp now carries eight declared lanes; review against the function list, not git status. The invariant that still holds is that a base seed determines a batch exactly.

Defect 0 (P0, systemic) — network generation is not reproducible across processes

Mechanism. ullmann_M_t is std::map<Node*, node_container_t*> (cpp/core/BNGcore.hpp:92), so iterating it yields ascending address order. UllmannSGIso::find_maps seeded its recursion at M.begin() and next_node advanced with ++row_iter, so the order in which pattern nodes were assigned — and therefore the order subgraph isomorphisms are emitted — was a function of heap addresses, i.e. ASLR. Downstream: ReactionRule.cpp:1546,1620 push embeddings in emission order, NetworkGenerator builds reactions in embedding order, RxnList::add preserves first-insertion order (RxnList.cpp:59-60). Two reactions tied on everything the .net writer prints came out in an address-dependent order.

Fix. Row order is now Ga order. Ga is a std::vector<Node*>, so its iteration order is insertion order and is stable across processes. Each level resolves its own row with M.find(rowOrder[d]); the map and its find/insert call sites are untouched, so copy_M's two-map zip is unaffected. next_node lost its row_iter_t& parameter, leaving the 3-arg copy_M overload dead — removed rather than left in place.

Proving test. 20 fresh-directory runs each, md5 of the emitted .net:

model	before	after
Motivating_example_cBNGL	3 variants	1
SHP2_base_model	3 variants	1
tlbr	2 variants	1

Re-measured after rebasing onto origin/main at 1c78a72 (12 runs each): 1, 1, 1.

One note on instrument strength, from swarmMemory's independent characterisation: their 40-run tlbr check compared species/reaction counts via generateNative, not file bytes, so identical counts there were never in tension with differing bytes here — a reorder of tied reactions moves one and not the other. The byte-level md5 above is the stronger instrument and is the one this fix is verified with.

All three share this single root cause — tlbr and SHP2 are not two bugs. Measured on main's shared binary at 6889fba (HEAD 6889fba; git status --porcelain --untracked-files=no | wc -l → 0; git log --since='2026-09-29 15:08:00' --format=%h -- cpp/ → 75b22a7, NFsim-only, never on this path) and on my own build after.

Defect 1 — observable count filters compared the wrong quantity

Mechanism. OdeIntegrator::compileGroups applied a stoichiometric comparison to matchCount, the number of pattern embeddings in a species graph, not a molecule count. On one species A() holding 10 molecules, A()>5 compared 1 > 5, failed, and contributed 0.

Semantics decision, stated rather than chosen silently. The grammar comment (BNGParser.g4:227) intends "stoichiometry comparison", but BNG2 2.9.3 is the oracle and it drops the threshold for Molecules entirely: Perl2/Observable.pm tests $patt->Quantifier only inside the Type eq "Species" branch (lines 278-292); the Molecules branch (lines 241-256) accumulates raw match counts and never inspects it. I matched BNG2. The other reading — comparing a true molecule count — would have returned 10 for A()>5 and A() alike too, but 0 for A()>10 and A()>99 where BNG2 returns 10; it is a different fix with different failure modes, and it would have been a deliberate divergence from the oracle on models that currently "work" in BNG2.

Where the 0 came from. The .net begin groups block carries no threshold in either engine — that is BNG2's own output shape, not a writer bug — so the number is produced at observable-evaluation time in compileGroups, not at lowering. A writer-only fix would not have changed it.

What a user gets, before and after. Model: A() 10, B() 0, reversible A() <-> B(), SSA seed 7, t=1. Columns are the four observables in source order. Measured, not reasoned:

	Plain A()	Gt5 A()>5	Gt10 A()>10	Gt99 A()>99
BNG2 2.9.3	10	10	10	10
BNG3 before	10	0	0	0
BNG3 after	10	10	10	10

So the concrete user-visible change is: a Molecules observable with a count predicate used to report 0 whenever the predicate was not satisfied by the per-species embedding count, and now reports the same value as the unfiltered observable. Before, A()>5 on 10 molecules returned 0 — a model that ran and reported a wrong number with no diagnostic. After, it returns 10, which is what BNG2 returns and what the pattern plainly contains.

Two consequences worth stating rather than leaving implicit:

1. This makes BNG2's own arguably-wrong behaviour the compatible one. The other reading — comparing a true molecule count — would also return 10 for A()>5, but 0 for A()>10 and A()>99 where BNG2 returns 10. That is a different fix with different failure modes, and it would have been a deliberate divergence from the oracle on models that currently work in BNG2. A user porting a model from BNG2 gets identical numbers either way for thresholds the model satisfies, and identical numbers for thresholds it does not.
2. Species observables are unaffected. The comparison is honoured there, in both engines, before and after. Measured term-for-term on the same models.

Proving test. test_correctness_regressions.cpp pins all six Molecules columns and the three Species columns, so a future change cannot quietly reintroduce this.

Defect 2 — batch-SSA base seeds aliased

Mechanism. integrateBatchSSA's CPU pool derived each trajectory's seed as base + traj. Addition aliases: two batch runs whose base seeds differ by delta share B - delta trajectories.

CPU/GPU coupling. The derivation moved to engine/BatchSsa.hpp::batchTrajectorySeed because every backend must agree — tests/test_batch_ssa_statistical_parity.py compares CPU against GPU at the same base seed, so a per-backend derivation would silently invalidate it. It is a pure function of (baseSeed, trajectory): no thread id, device id, or loop order. perfBatch holds BatchSsa.cpp:235,365 and is mirroring it there.

Acceptance, all four met. Across-seed spread of the batch mean against the analytic Monte-Carlo sd, 20 base seeds per arm:

criterion	before	after
base seeds 5000/5001/5002	0.79835 / 0.79835 / 0.79835	0.79855 / 0.80195 / 0.80585
near-adjacent spread ratio	0.051	1.48
far-apart spread ratio (reference)	1.276	1.22
per-batch observable sd (exact 0.4)	0.3988	0.3977-0.3988

GPU parity is NOT verified and I do not claim it. Every model available on this host is below the Metal launch-overhead threshold — this binary logs metal backend skipped, batch of 2000 over 2 reactions is below the launch-overhead threshold — so the GPU path was never exercised. The risk I am accepting is stated rather than hidden: the GPU kernels keep their existing counter-based PCG32 derivation (MetalSsaBackend.mm:123, CudaSsaBackend.cu:215), so CPU and GPU now derive seeds by different formulas. Both are pure functions of (base, traj) and both decorrelate adjacent seeds, so the parity test compares two decorrelated samples as it intends — but the streams differ, and only a maintainer with a working accelerator can confirm.

Report-only — untagged complex product

Reproduced, then declined to fix, because no landing is unambiguous.

S() + E() -> S().E() generates an empty .net and simulates with nothing happening. But BNG2 2.9.3 also generates nothing here — with the complex as a seed species it aborts (ABORT: Species S().E() 1 is not connected); with no seed species it warns (WARNING: No transformations detected in reaction rule (S() + E() -> S().E() _rateLaw1). Please verify that this was your intent.) and generates zero reactions. So this is not an entire class of mass-action models silently generating nothing; it is one spelling that neither engine supports, where BNG3 is silent and BNG2 is not.

I attempted the fail-closed refusal and reverted it: no syntactic discriminator exists. S().E() reaches buildPatternGraph as two disconnected molecule nodes, but so do R.R(tf~Y) and TF.R.R.TF(d~pY) (Motivating_example.bngl:141,152) and @EM:R.R (Motivating_example_cBNGL.bngl:109) — and all three generate normally in both engines (356 reaction lines for Motivating_example under BNG2 and here). Every scoping I tried over-fired on shipped validation models; the only variant that fired nowhere over-refused nothing useful. Shipping it would have been a regression against the oracle, so it is reverted and reported.

What would settle it: a BNG2 model that forms a hetero-oligomer through a component-free . and generates ≥1 reaction. I looked for one and did not find it; sciMetabolic reported the same gap. Without such an oracle, "instantiate the identity mapping" would invent semantics BNG2 does not have, and "refuse" would reject models BNG2 accepts. The narrow, safe landing that remains open is emitting BNG2's warning verbatim when a reaction rule yields no transformation — a diagnostic parity fix, not a semantics decision, and it needs no new semantics.

Verification tokens, per commit

commit	gate	result
53d2e60 determinism	20 + 12 fresh-dir runs x 3 models	3→1, 3→1, 2→1 (both before and after rebase)
53d2e60 determinism	ctest -j4	479/479
f083d7e observables	BNG2 vs before vs after, 4 columns	exact agreement with BNG2 after; 0s before
f083d7e observables	ctest -j4	479/479
d165ff5 seeds	4 acceptance criteria from the reporter	all met (table above)
d165ff5 seeds	independent: orchScope ran PR #31's closed-form harness against this binary	10 passed, both paths to this worktree
d165ff5 seeds	ctest -j4	479/479

orchScope's run is the one instrument that can judge the seed commit independently, because none of its ten assertions compares against a recorded value — exact integer mass balance, pool conservation, Poisson mean and variance, the multinomial cross-moment, binomial thinning, the two-sided k*C(n,2) convention, and bit-identical reproducibility across threads = 0/1/3.

Verification

Failing-first demonstrated: reverting all five source files to 6889fba and rebuilding gives 3 failed / 5 passed (8 cases, 83 assertions); with the fixes, all 9 pass (94 assertions). (The defect-1 tests were removed with the fix.)

• ctest --test-dir build -j4 → 100% tests passed out of 479 (post-rebase, on origin/main at 1c78a72)
• python -m pytest tests/python → 958 passed, 28 skipped (post-rebase on 1c78a72; 862 before the rebase, the difference being other lanes' merged tests), with both paths asserted inside a pytest session using the sys.modules probe:
  • UNDER TEST pkg: /Users/akutuva/Documents/BioNetGen/BNG3-correctness/python/bionetgen/__init__.py
  • UNDER TEST ext: /Users/akutuva/Documents/BioNetGen/BNG3-correctness/build/cpp/_bionetgen_cpp.cpython-314-darwin.so
• Build: 643 steps, 183.55 s wall / 448.36 s user, 477,265,920 B peak RSS. No timing claim is made — this is a correctness PR.

Co-Authored-By: Claude Opus 4.8 (1M context) noreply@anthropic.com