[wshlavacek]

Closes #109. Stacked on #116, which refuses every Arrhenius rule the load-time expansion cannot price exactly. This makes those rules run at BioNetGen's rate instead. The base is fix/arrhenius-context-refusal, so the diff shows this change alone; GitHub retargets it to main once #116 merges.

What changes

The expansion stays wherever it is exact, and those models run byte-identical. Every other Arrhenius rule now runs as itself, with its own reactant and product patterns, and its dG is counted at each firing the way BioNetGen's network counts it per reaction.

What BNG2 counts. BNG2 weighs an energy pattern in a species by matches / automorphisms (EnergyPattern.pm update; the automorphisms are the pattern matched on itself, SpeciesGraph.pm). Both are counted with BNG2's own component test, under which a !+ lands on a !? and a ~? needs some state, so L(r!+,r!?) has two automorphisms. A reaction's dG is Σ G × (weight in products − weight in reactants).

Why the count is local. For a connected pattern, the matches a bond changes are exactly those that put some pattern component on one of the bond's two components; any other match reads only components the bond leaves alone. So a firing's dG is the change in those anchored matches, counted on the complexes the bond joins or parts, with the bond made or broken virtually. No species, no network.

The matcher. energy_match.{hpp,cpp} is a matcher templated over a small graph view. The engine runs it on its pool, the automorphism count runs it on the pattern itself, and the unit tests run it on a hand-built graph.

Thinning. The rule rides the local-rate machinery, as a Function (unbinding) or FunctionProduct (binding) rate law with per-molecule weights. The engine draws a firing against those weights, counts the drawn instance's rate, and keeps the firing with probability rate / weight.

• Unbinding: a molecule's weight is the largest counted rate among its bonds of the rule's kind.
• Binding:
  • dG is cut into a term each end's own complex decides, plus, for each pattern bond of the rule's kind, a product of the matches on the two sides of that bond.
  • Where one side is the cut site alone, the product is exactly the other side's matches, and that end's weight carries it. So the weights are exact unless an energy pattern has context on both sides of the bond.
  • For such a pattern, end 0's weight takes the far side at its most matches where that raises the rate, and none where it lowers it.
  • Where that bound could make a pair be drawn more than e times as often as it fires, the rule is split into one rule per class of far-side matches, up to 8. Each class rule keeps only the pairs whose far-side site is in its class, so each pair is priced exactly. The class rules are named with _arr_cK appended, as expanded rules take _arr_vK.
  • A bound left larger is reported with a warning. Splitting costs a weight computation per class, so small bounds stay bounded.
• Same complex: a binding pair already in one complex never fires, as BNG2 applies a rule with two reactant patterns, including under -no-bscb.
• set_param: re-resolves phi, Ea0 and the energies.

Still refused, with the rule and reason named: a state change, two bonds at once, ring closure, a homodimer, include/exclude constraints, and, off the expansion, an energy pattern joined only by . that the rule's bond can reach. BNG2 counts such a pattern per complex, so joining or parting two complexes can change it. It is let through for a pattern naming a molecule type no bond can bring near the rule's, and, beside an unbinding written with . (which never parts a complex), for a pattern that asks nothing of the bond's two sites.

Verification

• BNG2 oracle. Two new models pin 29 shapes against BioNetGen's network: the issue's nine, and twenty more.
  • arrhenius_context (ode) covers unidirectional unbinding, one column pair per shape. The shapes beyond the issue's are:
    • a context partner that only a rule makes;
    • a context site declared twice;
    • a !? counted twice;
    • ring opening by a rule written with .;
    • a third molecule and a free site in the rule's reactant pattern;
    • a molecule holding two bonds of different rates;
    • BNG2's automorphisms for !+ beside !?, and for ~?.
  • arrhenius_context_bind (ssa, 1000 reps) covers reversible binding, both halves on the exact path:
    • cooperativity across the bond, favourable and not;
    • context on the far end only;
    • a far side that can match twice;
    • wildcard-only and repeated-site patterns, states and free sites;
    • state and bond context in the rule's own reactant pattern;
    • a two-sided pattern split by class;
    • a split pattern on a molecule with a free site in each class.
  • Main misses the bind model by 61 standard errors, and the unbinding model by an unbounded z on the ring opening, which never fires there. Refuse Arrhenius rules the load-time expansion cannot price exactly (#109) #116 refuses both. Here: max z 2.70 and 1.99.
• Unit tests. energy_match_test checks the counting against BNG2's weights by hand, for the issue's cases and !? C(m, k). It also checks automorphism counts, the binding cut against the direct count, cycles, molecule injectivity and the structural bound.
• Other tests.
  • arrhenius_context_test now pins what is still refused, that the exact path never binds a pair inside one complex under -no-bscb, and that an XML naming its AddBond sites the other way round runs bit-identical.
  • arrhenius_set_param_test runs both new models under no-op overrides. It also runs two overrides that change the answer: the context energy switched off, scored against the bare column, and phi = 0, scored in closed form.
• Bias check at high precision. The simplest exact-path rule, which has only reactant context, run at 5000 seeds against the same rule on the expansion path, agrees within 1 SE at every time point. A wildcard-only binding pattern at 2000 seeds per side against BNG2's SSA gives max z 1.33.
• Mutation check. 47 single-point mutations of the matcher, the weights, the class split, the accept step, the loader wiring and the set_param re-resolution: 40 killed. The other seven:
  • two are equivalent: a fit check a later check repeats, and a factor that is 1 wherever it applies;
  • five change only how often draws are thinned, never the rate: folding a one-sided product on both ends, not recognising a single-site side, admitting patterns that cannot match, splitting by class, and a looser structural bound. No test can see those; the null-event count is not in the API.
• Independent review. A reviewer scored its own probe models against BNG2.
  • It found two automorphism miscounts (a !+ beside a !?, z = 39 against BNG2's SSA; a ~?, z = 23).
  • It found a crash on an exact-path rule with no energy patterns, a heap read past the end on a disconnected pattern, and a two-sided bound loose enough that a stress model drew a million null events without firing.
  • A second round found two more: a split rule taking a drawn site outside its class (z = 9.2), and a ring-opening exemption that dropped a disconnected pattern constraining a bond site (z = 97).
  • All of these are fixed and pinned above. Its other probes ran at BNG2's rate: rings, cyclic patterns, multi-molecule reactant patterns on both sides, reactant patterns written B-first, phi = 1.4 and −0.3, reversible halves, split reverse halves, and splits over two cuts.
• Suite. ctest 66/66; basicmodels.py 29/29.
• Blast radius. 251 of 259 XMLs under tests/ run byte-identical to main at 2 seeds, stdout and exit code. The 8 that differ are:
  • the two new oracle models;
  • arrhenius_formula, extended in Refuse Arrhenius rules the load-time expansion cannot price exactly (#109) #116;
  • four new fixtures main misprices: bivalent, no_energy_patterns, ep_disconnected_apart and ep_disconnected_ring;
  • ep_disconnected, which this refuses.
• Cost. The exact path recomputes the weights of every molecule in each complex an event changes, so its cost per event grows with the size of the complexes events touch.
  • On a bivalent ligand–receptor system that keeps forming and breaking chains (2,000 to 20,000 molecules), it measured 20 to 60 µs per event against 2.5 µs for the same rule expanded, with 2 null events in 22,709.
  • On the reviewer's stress models (trivalent and bivalent molecules under six multi-bond patterns), 15 to 70% of draws are thinned.

Not covered

The reverse-half ambiguity documented in #116 remains. Two unidirectional rules, X and _reverse_X, that mirror each other and name the same parameters write the same XML as one reversible rule. Telling them apart needs BNG2's writeXML to mark reverse halves.