From ce9d2074e2dc99bdbfc55e0d630232f588eb5dc3 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Thu, 22 Oct 2026 21:56:48 -0600 Subject: [PATCH] rm: run the Arrhenius rules RuleMonkey 3.11.0 misread (#1058) RuleMonkey 3.12.0 (vendored in #1220) gives every Arrhenius rule it runs BNG2's rate: it expands a rule at load where that is exact, and elsewhere counts dG at each firing as BNG2's network does per reaction (richardposner/RuleMonkey#109). Shapes it cannot run, it refuses itself, naming the rule and the reason. So bngsim's own gate, refuse_rulemonkey_energy_misreads, has nothing left to catch. Deleted from _network_free_checks.py: refuse_rulemonkey_energy_misreads, rulemonkey_energy_misreads, and the helpers only they used (_arrhenius_centers, _expanded_exactly, _bond_partners, _component_names, _plain_bond_pattern, _graph, _components, Site), with the module docstring's paragraph on them. Their calls in Simulator and RuleMonkeySession.__init__ are gone too. _local and _children stay; other checks use them. All fifteen misread/ models load on 3.12.0. None is a shape RuleMonkey still refuses (a state change, a homodimer, a ring closure, two bonds, constraints, a bond that can reach a '.'-joined energy pattern). So test_rm_refuses_an_energy_rule_it_would_misread is now three value tests per model: - method='rm', a RuleMonkeySession, and network SSA on the .net, each against the closed form the model's .bngl header gives: each observable is binomial(1000, exp(-k t)), or fixed at 1000 where BNG2 never fires the rule (rule_context's s~0, ringopen's open chain); - method='rm' against network SSA, two-sample z per time point. Each replicate gets a fresh Model or session, 40 for rm, 200 for SSA. On 3.12.0, swapping each closed-form rate for the misread one, the rate with the context dropped or misapplied, fails every observable the misread changes, at |z| from 50 to 515 on rm's runs. New model state_pattern: a reversible binding rule with the energy pattern A(b!1).B(a!1,s~1), 100 each of A, B(s~0) and B(s~1), generated with BNG2 2.9.3. Detailed balance gives the exact stationary means, AB0 = 22.33 and AB1 = 68.20; dropping the state would put both at 39.33. rm settles on them at t = 5 and matches network SSA along the way, and network SSA settles on them too. test_rulemonkey_names_an_energy_rule_shape_it_refuses_itself now matches RuleMonkey's own message naming the rule ("rule '_R1' is outside the shapes RM runs: it changes a state"). The network-free guide's "Arrhenius (energy) rules" section says what rm runs, what RuleMonkey refuses, and one case RuleMonkey reads differently from BNG2: two unidirectional mirror-image rules X and _reverse_X with the same phi and Ea0 parameters write the XML of one reversible rule, and the second gets 1-phi. The 18 files that drive RuleMonkey: 496 passed, 3 skipped (448 before). --- changelog.d/1058.fixed.md | 9 + docs/user-guide/network-free.md | 30 +- python/bngsim/_network_free_checks.py | 247 -------------- python/bngsim/_rulemonkey_session.py | 2 - python/bngsim/_simulator.py | 2 - python/tests/test_rulemonkey_arrhenius.py | 174 +++++++--- tests/data/nfsim/arrhenius/state_change.bngl | 5 +- tests/data/nfsim/arrhenius/state_pattern.bngl | 37 +++ tests/data/nfsim/arrhenius/state_pattern.net | 28 ++ tests/data/nfsim/arrhenius/state_pattern.xml | 306 ++++++++++++++++++ 10 files changed, 534 insertions(+), 306 deletions(-) create mode 100644 changelog.d/1058.fixed.md create mode 100644 tests/data/nfsim/arrhenius/state_pattern.bngl create mode 100644 tests/data/nfsim/arrhenius/state_pattern.net create mode 100644 tests/data/nfsim/arrhenius/state_pattern.xml diff --git a/changelog.d/1058.fixed.md b/changelog.d/1058.fixed.md new file mode 100644 index 000000000..5280e52da --- /dev/null +++ b/changelog.d/1058.fixed.md @@ -0,0 +1,9 @@ +- **`method="rm"` runs Arrhenius rules whose energy patterns carry context + (issue #1058).** RuleMonkey 3.12.0 gives them BNG2's rate, expanding a rule at + load where that is exact and counting `dG` at each firing elsewhere + (richardposner/RuleMonkey#109), so `rm` and `RuleMonkeySession` no longer + refuse an energy pattern with a state, a free or wildcard site, a site named + twice or context beyond one bond, context in the rule's own reactant pattern, + a ring-opening unbinding, or a unidirectional rule named `_reverse_...`. + Shapes RuleMonkey cannot run (a state change, a ring closure, a homodimer, + two bonds, constraints) are still refused, by RuleMonkey, naming the rule. diff --git a/docs/user-guide/network-free.md b/docs/user-guide/network-free.md index 525ae9040..983c4923d 100644 --- a/docs/user-guide/network-free.md +++ b/docs/user-guide/network-free.md @@ -200,15 +200,27 @@ reverse half of a reversible rule. Energies are in units of RT: there is no RT in the formula, and a parameter named `RT` is an ordinary parameter, which a model in kJ/mol divides by where it uses it (issue #1022). -RuleMonkey expands each energy rule at load, and reads only some energy -patterns as BNG2 does: the bonded pair of sites a rule joins or parts, alone or -with one more bond to a molecule of a third type whose site binds nothing else, -on a rule that names only those two sites. That covers pair and cooperativity -patterns. It drops or misapplies a state, a site that must be free, a wildcard -bond, a site named twice, context further out, and context in the rule itself -(richardposner/RuleMonkey#109), so `method="rm"` refuses such a model by name -until that is fixed upstream (issue #1058). Simulate it with `method="ode"` or -`method="ssa"` on the generated network. +RuleMonkey runs an energy rule that makes one bond between two reactant +patterns, or breaks one within one, between molecules of two types, with no +other operation. It expands the rule at load where that is exact, as for pair +and cooperativity patterns, and everywhere else counts `dG` at each firing as +BNG2's network does per reaction: an energy pattern with a state, a site that +must be free, wildcard bonds, a site named twice or context further out, and +context in the rule's own reactant pattern all run at BNG2's rate +(richardposner/RuleMonkey#109; `method="rm"` refused these until issue #1058). +Counting at each firing costs more per event, and more the larger the complexes +an event touches. + +RuleMonkey refuses at load, naming the rule and the reason, an energy rule that +changes a state or adds or deletes a molecule, makes or breaks more than one +bond, closes a ring, binds two molecules of one type, has more than two +reactants, or carries include or exclude constraints, and a rule it cannot +expand whose bond can reach an energy pattern joined only by `.`. Simulate such +a model with `method="ode"` or `method="ssa"` on the generated network. BNG2 +writes the reverse half of a reversible rule `X` as a rule named `_reverse_X`, +so two unidirectional rules, `X` and a mirror image named `_reverse_X` with the +same `phi` and `Ea0` parameters, read as one reversible rule, and RuleMonkey +gives the second `1-phi` where BNG2 gives it `phi`. ### Compartments diff --git a/python/bngsim/_network_free_checks.py b/python/bngsim/_network_free_checks.py index dd91ae9e5..45168ec08 100644 --- a/python/bngsim/_network_free_checks.py +++ b/python/bngsim/_network_free_checks.py @@ -24,15 +24,6 @@ A RuleMonkey session uses :func:`clock_readers` to refuse a segment whose labels are not the clock the model reads. -Issue #1058: RuleMonkey 3.11.0 gives an Arrhenius rule BNG2's rate, -exp(-(Ea0 + phi*dG)), where every energy pattern the rule's bond is part of is -that bonded pair of sites, alone or with one more bond to a third molecule -whose site binds nothing else, and the rule names only those two sites. It -expands each energy rule at load, and outside that it drops or misapplies the -context (richardposner/RuleMonkey#109), so the rule fires at the wrong rate -with no warning. :func:`refuse_rulemonkey_energy_misreads` refuses such a model -by name until the upstream fix is vendored. - An Arrhenius (energy-based) rate law is refused on NFsim (issue #854's review). The vendored NFsim expands energy rules at load with three silent defects (RuleWorld/nfsim#102): it divides every energy by an ``RT`` that BioNetGen does @@ -260,244 +251,6 @@ def _children(el: ET.Element, name: str) -> list[ET.Element]: return [c for c in el if _local(c.tag) == name] -def _graph(el: ET.Element) -> tuple[list[ET.Element], dict[str, tuple[int, str, str]], list]: - """The molecules of a pattern or species, its sites (component id -> - molecule index, molecule type, component name) and its bonds (pairs of - component ids).""" - molecules = [m for lst in _children(el, "ListOfMolecules") for m in lst] - sites = {} - for i, mol in enumerate(molecules): - for lst in _children(mol, "ListOfComponents"): - for comp in lst: - sites[comp.get("id") or ""] = (i, mol.get("name") or "", comp.get("name") or "") - bonds = [ - (b.get("site1") or "", b.get("site2") or "") - for lst in _children(el, "ListOfBonds") - for b in lst - ] - return molecules, sites, bonds - - -def _components(mol: ET.Element) -> list[ET.Element]: - return [c for lst in _children(mol, "ListOfComponents") for c in lst] - - -def _plain_bond_pattern(pattern: ET.Element) -> bool: - """Whether a pattern names nothing but explicit bonds between its own - molecules: every component bound once to a partner the pattern lists, with - no state, no site named twice on one molecule, and no compartment.""" - molecules, sites, bonds = _graph(pattern) - if not molecules or pattern.get("compartment"): - return False - ends = [site for bond in bonds for site in bond] - for mol in molecules: - comps = _components(mol) - names = [c.get("name") for c in comps] - if mol.get("compartment") or not comps or len(set(names)) != len(names): - return False - for comp in comps: - if comp.get("state") is not None or comp.get("numberOfBonds") != "1": - return False - if ends.count(comp.get("id") or "") != 1: - return False - return all(a in sites and b in sites for a, b in bonds) - - -Site = tuple[str, str] # (molecule type, component name) - - -def _bond_partners(root: ET.Element) -> dict[Site, set[Site]]: - """For every site, the sites a bond in the model can join it to: the bonds - of the seed species and of every rule's patterns. A bond a rule adds, to a - molecule it synthesizes too, is in its product pattern.""" - partners: dict[Site, set[Site]] = {} - for el in root.iter(): - if _local(el.tag) in ("Species", "ReactantPattern", "ProductPattern"): - _, sites, bonds = _graph(el) - for a, b in bonds: - if a in sites and b in sites: - sa, sb = sites[a][1:], sites[b][1:] - partners.setdefault(sa, set()).add(sb) - partners.setdefault(sb, set()).add(sa) - return partners - - -def _component_names(root: ET.Element) -> dict[str, list[str]]: - """Each molecule type's component names, a name repeated as declared.""" - names: dict[str, list[str]] = {} - for el in root.iter(): - if _local(el.tag) == "MoleculeType": - names[el.get("id") or ""] = [ - c.get("id") or "" for lst in _children(el, "ListOfComponentTypes") for c in lst - ] - return names - - -def _arrhenius_centers(root: ET.Element) -> tuple[list[tuple[Site, Site]], list[str], bool]: - """The bond each Arrhenius rule makes or breaks, as a pair of sites; the - rules RuleMonkey would run on more than those two sites, or with the wrong - products; and whether some Arrhenius rule has a shape RuleMonkey refuses - itself at load (try_expand_arrhenius in third_party/rulemonkey), which then - names it better than this check could.""" - centers: list[tuple[Site, Site]] = [] - with_context: list[str] = [] - rule_names = {r.get("name") for r in root.iter() if _local(r.tag) == "ReactionRule"} - for rule in root.iter(): - if _local(rule.tag) != "ReactionRule" or not any( - (rl.get("type") or "").lower() == "arrhenius" for rl in _children(rule, "RateLaw") - ): - continue - reactants = [rp for lst in _children(rule, "ListOfReactantPatterns") for rp in lst] - products = [pp for lst in _children(rule, "ListOfProductPatterns") for pp in lst] - molecules: list[ET.Element] = [] - sites: dict[str, tuple[int, str, str]] = {} - pattern_of: dict[str, int] = {} - for k, rp in enumerate(reactants): - mols, rp_sites, _ = _graph(rp) - sites.update({c: (len(molecules) + v[0], v[1], v[2]) for c, v in rp_sites.items()}) - pattern_of.update(dict.fromkeys(rp_sites, k)) - molecules += mols - ops = [op for lst in _children(rule, "ListOfOperations") for op in lst] - # RuleMonkey refuses anything but one AddBond joining two reactant - # patterns, or one DeleteBond inside one, between two molecule types. - if len(ops) != 1 or _local(ops[0].tag) not in ("AddBond", "DeleteBond"): - return [], [], True - a, b = ops[0].get("site1") or "", ops[0].get("site2") or "" - if a not in sites or b not in sites or sites[a][1] == sites[b][1]: - return [], [], True - binding = _local(ops[0].tag) == "AddBond" - if binding and (len(reactants) != 2 or pattern_of[a] == pattern_of[b]): - return [], [], True - if not binding and len(reactants) != 1: - return [], [], True - centers.append((sites[a][1:], sites[b][1:])) - # It builds the expanded rule from the two sites alone: one complex - # from two when binding, two from one when unbinding (a bond inside a - # ring, whose products stay one complex, would part it). - bare = len(products) == (1 if binding else 2) and len(molecules) == 2 - bare = bare and not any(rp.get("compartment") for rp in reactants) - for mol in molecules: - comps = _components(mol) - if ( - mol.get("compartment") - or len(comps) != 1 - or comps[0].get("state") is not None - or (comps[0].get("id") or "") not in (a, b) - ): - bare = False - # It takes a rule named _reverse_X for the reverse half of rule X and - # gives it 1-phi, which a one-way rule of that name is not. - name = rule.get("name") or "" - if name.startswith("_reverse_") and name[len("_reverse_") :] not in rule_names: - bare = False - if not bare: - with_context.append(name or rule.get("id") or "?") - return centers, with_context, False - - -def _expanded_exactly( - pattern: ET.Element, - center: tuple[Site, Site], - partners: dict[Site, set[Site]], - names: dict[str, list[str]], -) -> bool: - """Whether RuleMonkey's expansion of an energy rule with bond *center* - reads the plain pattern *pattern* as BNG2 does: it is the center's two - sites bonded, with at most one more bond, from one of them to a molecule of - a third type that names only that site, which nothing else binds.""" - molecules, sites, bonds = _graph(pattern) - kind = set(center) - on_center = [bd for bd in bonds if {sites[bd[0]][1:], sites[bd[1]][1:]} == kind] - if not on_center: - return True # it holds no bond the rule makes or breaks - if len(on_center) > 1: - return False - core = {sites[on_center[0][0]][0], sites[on_center[0][1]][0]} - others = [bd for bd in bonds if bd != on_center[0]] - leaves = [i for i in range(len(molecules)) if i not in core] - if not leaves: - return not others - if len(leaves) > 1 or len(others) != 1: - return False - a, b = others[0] - leaf_end, core_end = (a, b) if sites[a][0] == leaves[0] else (b, a) - if sites[leaf_end][0] != leaves[0] or sites[core_end][0] not in core: - return False - if sites[leaf_end][1] in {site[0] for site in center}: - return False - if len(_components(molecules[leaves[0]])) != 1: - return False - # RuleMonkey adds one "bound" or "free" copy of the context site to the - # expanded rule, so a type with two sites of its name, or of the center - # site's, matches both variants or counts one embedding for two. - core_type, context = sites[core_end][1], sites[core_end][2] - center_site = next(n for t, n in center if t == core_type) - declared = names.get(core_type, []) - if declared.count(context) != 1 or declared.count(center_site) != 1: - return False - # RuleMonkey tests the context site as bound to anything (energy_expand.cpp). - return partners.get(sites[core_end][1:], set()) <= {sites[leaf_end][1:]} - - -def rulemonkey_energy_misreads(xml_path: str | Path) -> tuple[list[str], list[str]]: - """The energy patterns (as BNGL) and the Arrhenius rules (by name) of - *xml_path* that RuleMonkey 3.11.0 expands differently from BNG2 (issue - #1058). Both empty for a model with no Arrhenius rule, for one with an - Arrhenius rule RuleMonkey refuses itself, or when the file cannot be - parsed.""" - try: - root = ET.parse(str(xml_path)).getroot() - except (ET.ParseError, OSError): - return [], [] - centers, with_context, left_to_rulemonkey = _arrhenius_centers(root) - if left_to_rulemonkey or (not centers and not with_context): - return [], [] - partners = _bond_partners(root) - names = _component_names(root) - patterns = [] - for ep in root.iter(): - if _local(ep.tag) != "EnergyPattern": - continue - graphs = _children(ep, "Pattern") or [ep] - if not all( - _plain_bond_pattern(g) - and all(_expanded_exactly(g, c, partners, names) for c in centers) - for g in graphs - ): - patterns.append(ep.get("pattern") or ep.get("id") or "?") - return patterns, with_context - - -def refuse_rulemonkey_energy_misreads(xml_path: str | Path) -> None: - """Raise :class:`ModelError` if RuleMonkey would expand an Arrhenius rule of - *xml_path* differently from BNG2 (issue #1058). - - RuleMonkey 3.11.0 drops or misapplies an energy pattern's context, and an - Arrhenius rule's own reactant context, when it expands the rule at load - (richardposner/RuleMonkey#109). Delete this check, its calls in - ``Simulator`` and ``RuleMonkeySession``, and the helpers only it uses, once - the upstream fix is vendored. - """ - patterns, rules = rulemonkey_energy_misreads(xml_path) - if not patterns and not rules: - return - parts = [] - if patterns: - parts.append("energy patterns " + ", ".join(f"'{n}'" for n in patterns)) - if rules: - parts.append("rules " + ", ".join(f"'{n}'" for n in rules)) - raise ModelError( - "RuleMonkey (method='rm' / 'nf_exact') cannot run this Arrhenius model: when it " - f"expands the energy rules at load it would misread {' and '.join(parts)} " - "(richardposner/RuleMonkey#109), so they would fire at the wrong rate. It reads an " - "energy pattern that is the bonded pair of sites a rule joins or parts, alone or " - "with one more bond to a molecule of a third type whose site binds nothing else, on " - "rules that name only those two sites and join two complexes or part one in two. " - "Simulate it with method='ode' or method='ssa' on the generated network " - "(issue #1058)." - ) - - def compartments(xml_path: str | Path) -> list[str]: """The ids of the compartments *xml_path* declares. Empty when the file cannot be parsed.""" diff --git a/python/bngsim/_rulemonkey_session.py b/python/bngsim/_rulemonkey_session.py index 797f20bef..69b6ea364 100644 --- a/python/bngsim/_rulemonkey_session.py +++ b/python/bngsim/_rulemonkey_session.py @@ -101,12 +101,10 @@ def __init__( from bngsim._network_free_checks import ( refuse_rulemonkey_compartments, - refuse_rulemonkey_energy_misreads, rulemonkey_clock, ) refuse_rulemonkey_compartments(self._xml_path) - refuse_rulemonkey_energy_misreads(self._xml_path) # until issue #1058 # The functions that read the clock, for which a segment's labels must # be the session clock (see simulate). Where the model reads it through # its functions alone, the core reads it as time() + __bngsim_t0, set diff --git a/python/bngsim/_simulator.py b/python/bngsim/_simulator.py index 54605f563..72e463dcc 100644 --- a/python/bngsim/_simulator.py +++ b/python/bngsim/_simulator.py @@ -1045,13 +1045,11 @@ def __init__( from bngsim._bngsim_core import RuleMonkeySimulator from bngsim._network_free_checks import ( refuse_rulemonkey_compartments, - refuse_rulemonkey_energy_misreads, rulemonkey_clock, ) from bngsim._rulemonkey_session import warn_rulemonkey_diagnostics refuse_rulemonkey_compartments(xml_path) - refuse_rulemonkey_energy_misreads(xml_path) # until issue #1058 # Elapsed time, which keeps event times small at a large t_start. A # model that reads the clock reads model time (issue #1022): through # time(), from a copy of its XML whose clock is time() + __bngsim_t0, diff --git a/python/tests/test_rulemonkey_arrhenius.py b/python/tests/test_rulemonkey_arrhenius.py index b4b8015a2..e028a73e2 100644 --- a/python/tests/test_rulemonkey_arrhenius.py +++ b/python/tests/test_rulemonkey_arrhenius.py @@ -13,17 +13,20 @@ stationary distribution for the equilibrium (detailed balance gives it a product form), and the closed form of a one-way unbinding rule, which is first order. -RuleMonkey 3.11.0 still expands some energy rules differently from BNG2: it -drops or misapplies an energy pattern's context, and an Arrhenius rule's own -reactant context (richardposner/RuleMonkey#109). rm refuses those by name until -issue #1058 lifts the refusal; the models are in arrhenius/misread/. +RuleMonkey 3.11.0 expanded some energy rules differently from BNG2: it dropped +or misapplied an energy pattern's context, and an Arrhenius rule's own reactant +context, and rm refused those by name. RuleMonkey 3.12.0 runs each at BNG2's +rate, expanded at load where that is exact and counted at each firing elsewhere +(richardposner/RuleMonkey#109), and issue #1058 lifted the refusal. The models +are in arrhenius/misread/, each one-way unbinding with a closed form, and +state_pattern, a reversible rule with a state in its energy pattern, which +settles on its exact stationary distribution. """ from __future__ import annotations import functools import math -import re import warnings from pathlib import Path @@ -51,14 +54,20 @@ def _has_rulemonkey() -> bool: @functools.cache def _runs(method: str, name: str, n: int, observables: tuple[str, ...]) -> np.ndarray: - """replicates x observables x times.""" + """replicates x observables x times, a fresh Model or session per replicate. + *method* is a Simulator method or ``"session"``, a RuleMonkeySession.""" out = [] for seed in range(1, n + 1): - m = bngsim.Model.from_net(str(ARR / f"{name}.net")) - kw = {"xml_path": str(ARR / f"{name}.xml")} if method == "rm" else {} with warnings.catch_warnings(): warnings.simplefilter("error", bngsim.RuleMonkeyWarning) - r = bngsim.Simulator(m, method=method, **kw).run(sample_times=T, seed=seed) + if method == "session": + with bngsim.RuleMonkeySession(str(ARR / f"{name}.xml")) as s: + s.initialize(seed=seed) + r = s.simulate(sample_times=T) + else: + m = bngsim.Model.from_net(str(ARR / f"{name}.net")) + kw = {"xml_path": str(ARR / f"{name}.xml")} if method == "rm" else {} + r = bngsim.Simulator(m, method=method, **kw).run(sample_times=T, seed=seed) out.append([np.asarray(r.observables[o]) for o in observables]) return np.asarray(out, dtype=float) @@ -156,47 +165,127 @@ def test_one_way_energy_rules_fire_one_way_at_bng2s_rate(): assert np.all(np.abs(z) < 4.5), z -# Each .bngl in misread/ says what BNG2's rate is. RuleMonkey 3.11.0 gets all -# but plain_two_a and samename_explicit wrong (|z| from 24 to 142 against BNG2's -# closed form, 40 seeds); those two lie outside the shapes the refusal can tell -# are safe from the XML. Issue #1058 turns these into value tests. +# Each .bngl in misread/ unbinds one-way and its header gives BNG2's rate for each +# kind of complex it seeds, 1000 of each, so an observable that counts one kind +# is binomial(1000, exp(-k t)) at t, and k = 0 holds it at 1000. RuleMonkey 3.11.0 +# read all but plain_two_a and samename_explicit wrong (|z| from 24 to 142 against +# these, 40 seeds), and rm refused all fifteen until issue #1058. +E09, E125 = math.exp(-0.9), math.exp(-1.25) MISREADS = { - "state_ctx": "energy patterns 'A(b!1).B(a!1,s~1)'", - "freectx": "energy patterns 'A(b!1,c).B(a!1)'", - "wild": "energy patterns 'L(r!+,r!+)'", - "samename": "energy patterns 'L(r!1,r!+).R(l!1)'", - "samename_partner": "energy patterns 'L(r!1,r!2).R(l!1).S(l!2)'", - "partner": "energy patterns 'A(b!1,c!2).B(a!1).C(a!2)'", - "chain2ctx": "energy patterns 'A(b!1).B(a!1,c!2).C(b!2,d!3).D(c!3)'", - "two_leaves": "energy patterns 'A(b!1,c!2,d!3).B(a!1).C(a!2).D(a!3)'", - # A(x,c,c): one "c bound" variant matches an A with one c of each kind. - "samectx": "energy patterns 'A(x!1,c!2).B(y!1).C(a!2)'", + # A state in the energy pattern: A(b!1).B(a!1,s~1). + "state_ctx": {"AB1": E125, "AB0": E09}, + # A site the pattern needs free: A(b!1,c).B(a!1). + "freectx": {"AB1": E125, "AB0": E09}, + # Wildcard bonds only, L(r!+,r!+): both bonds of L.R.R break at e^-1.25. + "wild": {"AB1": 2 * E125}, + # L(r!1,r!+).R(l!1) needs the other r bound, which L(r!1,r).R(l!1) lacks. + "samename": {"AB1": E09}, + # L(r!1,r!2).R(l!1).S(l!2) applies only where the other r holds an S. + "samename_partner": {"LRS": E125, "LRR": 2 * E09}, + # A's c binds C or D, and A(b!1,c!2).B(a!1).C(a!2) only C. + "partner": {"ABC": E125, "ABD": E09}, + # Context two bonds out: A(b!1).B(a!1,c!2).C(b!2,d!3).D(c!3). + "chain2ctx": {"ABCD": E125, "ABC": E09}, + # Two context bonds, both needed: A(b!1,c!2,d!3).B(a!1).C(a!2).D(a!3). + "two_leaves": {"ABCD": E125, "ABC": E09}, + # A(x,c,c): one embedding of A(x!1,c!2).B(y!1).C(a!2) per bound c. + "samectx": {"AB1": E125, "AB2": math.exp(-1.6)}, # Z, not C, binds A.c, through a rule that synthesizes Z bound. - "synth": "energy patterns 'A(x!1,c!2).B(y!1).C(a!2)'", - "rule_context": "rules '_R1'", - # Unbinding inside a ring: the product stays one complex. - "ringopen": "rules '_R1'", - # A one-way rule named like a reverse half, which RuleMonkey gives 1-phi. - "revname": "rules '_reverse_u'", - "plain_two_a": "energy patterns 'A(b!1).B(a1!1,a2!2).A(b!2)'", - "samename_explicit": "energy patterns 'L(r!1,r!2).R(l!1).R(l!2)'", + "synth": {"AB": E09}, + # A state in the rule's own reactant pattern: only A(s~1).B unbinds. + "rule_context": {"AB1": E09, "AB0": 0.0}, + # Unbinding written with '.': the ring opens, the chain never parts. + "ringopen": {"Ring": E09, "Open": 0.0}, + # A unidirectional rule named like a reverse half, which BNG2 gives phi. + "revname": {"AB": E09}, + # Two A on one B, A(b!1).B(a1!1,a2!2).A(b!2); only the a1 bond breaks. + "plain_two_a": {"Full": E125}, + # A site named twice, with explicit bonds: L(r!1,r!2).R(l!1).R(l!2). + "samename_explicit": {"Full": 2 * E125}, } +def _assert_decays(counts: np.ndarray, k: float) -> None: + """*counts* (replicates x times) of 1000 complexes that each leave at rate *k*.""" + if k == 0: + assert np.all(counts == 1000), counts.min() + return + q = np.exp(-k * np.asarray(T)) + assert np.all(counts[:, 0] == 1000) + assert np.all(np.diff(counts, axis=1) <= 0) # one-way: nothing comes back + se = np.sqrt(1000 * q[1:] * (1 - q[1:]) / len(counts)) + z = (counts.mean(axis=0)[1:] - 1000 * q[1:]) / se + assert np.all(np.abs(z) < 4.5), z + + +@pytest.mark.parametrize("method", ["rm", "session", "ssa"]) +@pytest.mark.parametrize("name", sorted(MISREADS)) +def test_energy_context_fires_at_bng2s_closed_form(name, method): + # ssa, on the generated network, checks the closed form itself. + obs = tuple(MISREADS[name]) + runs = _runs(method, f"misread/{name}", N_SSA if method == "ssa" else N_RM, obs) + for j, o in enumerate(obs): + _assert_decays(runs[:, j, :], MISREADS[name][o]) + + @pytest.mark.parametrize("name", sorted(MISREADS)) -def test_rm_refuses_an_energy_rule_it_would_misread(name): - xml = ARR / "misread" / f"{name}.xml" - match = re.escape(MISREADS[name]) + r".*RuleMonkey#109.*issue #1058" - m = bngsim.Model.from_net(str(ARR / "misread" / f"{name}.net")) - with pytest.raises(bngsim.ModelError, match=match): - bngsim.Simulator(m, method="rm", xml_path=str(xml)) - with pytest.raises(bngsim.ModelError, match=match): - bngsim.RuleMonkeySession(str(xml)) +def test_energy_context_matches_network_ssa(name): + obs = tuple(MISREADS[name]) + rm = _runs("rm", f"misread/{name}", N_RM, obs) + ssa = _runs("ssa", f"misread/{name}", N_SSA, obs) + for j in range(len(obs)): + _assert_same_mean(rm[:, j, :], ssa[:, j, :]) + + +@functools.cache +def _state_pattern_equilibrium() -> dict[str, float]: + """Exact stationary means of AB0 and AB1 for state_pattern. + + 100 A, 100 B(s~0) and 100 B(s~1), and nothing changes s. Detailed balance + makes the stationary distribution prod_s w_s^n_s / n_s! with w = 1 for a + monomer, exp(-3.5) for A.B(s~0) and exp(-(3.5 - 2)) for A.B(s~1). + """ + n = 100 + x = np.arange(n + 1)[:, None] # AB0 + y = np.arange(n + 1)[None, :] # AB1 + a = n - x - y + ok = a >= 0 + log_fact = np.array([math.lgamma(k + 1.0) for k in range(n + 1)]) + + def lgam(k): + return log_fact[np.where(ok, k, 0)] + + logw = -3.5 * x - 1.5 * y - lgam(a) - lgam(n - x) - lgam(n - y) - lgam(x) - lgam(y) + logw = np.where(ok, logw, -np.inf) + p = np.exp(logw - logw.max()) + p /= p.sum() + return {"AB0": float((p * x).sum()), "AB1": float((p * y).sum())} + + +def test_a_reversible_rule_with_a_state_energy_pattern_matches_network_ssa(): + obs = ("AB0", "AB1") + rm = _runs("rm", "state_pattern", N_RM, obs) + ssa = _runs("ssa", "state_pattern", N_SSA, obs) + for j in range(len(obs)): + _assert_same_mean(rm[:, j, :], ssa[:, j, :]) + + +@pytest.mark.parametrize("method", ["rm", "ssa"]) +def test_a_reversible_rule_with_a_state_energy_pattern_settles_on_the_exact_equilibrium( + method, +): + # ssa, on the generated network, checks the oracle itself. + obs = ("AB0", "AB1") + runs = _runs(method, "state_pattern", N_SSA if method == "ssa" else N_RM, obs)[:, :, -1] + exact = _state_pattern_equilibrium() + for j, o in enumerate(obs): + se = runs[:, j].std(ddof=1) / math.sqrt(len(runs)) + assert abs(runs[:, j].mean() - exact[o]) < 4.5 * se, (o, runs[:, j].mean(), exact[o]) def test_rulemonkey_names_an_energy_rule_shape_it_refuses_itself(): """A state-change Arrhenius rule is refused by RuleMonkey's own load check, - whose message says what it is; this check stands aside.""" + whose message names the rule and says what it is.""" xml = str(ARR / "state_change.xml") m = bngsim.Model.from_net(str(ARR / "state_change.net")) for build in ( @@ -204,10 +293,9 @@ def test_rulemonkey_names_an_energy_rule_shape_it_refuses_itself(): lambda: bngsim.RuleMonkeySession(xml), ): with pytest.raises( - RuntimeError, match="outside the shapes RM runs: it changes a state" - ) as err: + RuntimeError, match="rule '_R1' is outside the shapes RM runs: it changes a state" + ): build() - assert "#109" not in str(err.value) def test_cooperativity_on_both_reaction_centers_matches_network_ssa(): diff --git a/tests/data/nfsim/arrhenius/state_change.bngl b/tests/data/nfsim/arrhenius/state_change.bngl index 8392cdb6d..d513ac6e9 100644 --- a/tests/data/nfsim/arrhenius/state_change.bngl +++ b/tests/data/nfsim/arrhenius/state_change.bngl @@ -1,6 +1,5 @@ -# Issue #1022: an Arrhenius state-change rule, which RuleMonkey 3.11.0 refuses -# itself at load ("not a supported 2-reactant binding/unbinding shape"). bngsim -# leaves the refusal to it, rather than blaming richardposner/RuleMonkey#109. +# Issue #1022: an Arrhenius state-change rule, which RuleMonkey refuses itself at +# load ("outside the shapes RM runs: it changes a state"), naming the rule. begin model begin parameters Ea0 0.2 diff --git a/tests/data/nfsim/arrhenius/state_pattern.bngl b/tests/data/nfsim/arrhenius/state_pattern.bngl new file mode 100644 index 000000000..297d5663e --- /dev/null +++ b/tests/data/nfsim/arrhenius/state_pattern.bngl @@ -0,0 +1,37 @@ +# Issue #1058: a reversible Arrhenius binding rule with an energy pattern that +# holds a state. B(s~1) binds A with G_AB + G_S1 = 1.5, B(s~0) with G_AB = 3.5, +# and nothing changes s. Detailed balance gives the stationary distribution +# prod_s w_s^n_s / n_s!, w = 1 for a monomer and exp(-G) for a complex. +# RuleMonkey 3.11.0's load-time expansion misread a state in an energy pattern +# (richardposner/RuleMonkey#109); 3.12.0 counts dG at each firing. +begin model +begin parameters + Ea0 0.5 + phi 0.3 + G_AB 3.5 + G_S1 -2 +end parameters +begin molecule types + A(b) + B(a,s~0~1) +end molecule types +begin seed species + A(b) 100 + B(a,s~0) 100 + B(a,s~1) 100 +end seed species +begin energy patterns + A(b!1).B(a!1) G_AB + A(b!1).B(a!1,s~1) G_S1 +end energy patterns +begin observables + Molecules AB0 A(b!1).B(a!1,s~0) + Molecules AB1 A(b!1).B(a!1,s~1) +end observables +begin reaction rules + A(b) + B(a) <-> A(b!1).B(a!1) Arrhenius(phi,Ea0) +end reaction rules +end model + +generate_network({overwrite=>1}) +writeXML() diff --git a/tests/data/nfsim/arrhenius/state_pattern.net b/tests/data/nfsim/arrhenius/state_pattern.net new file mode 100644 index 000000000..42499a12f --- /dev/null +++ b/tests/data/nfsim/arrhenius/state_pattern.net @@ -0,0 +1,28 @@ +# Created by BioNetGen 2.9.3 +begin parameters + 1 Ea0 0.5 # Constant + 2 phi 0.3 # Constant + 3 G_AB 3.5 # Constant + 4 G_S1 -2 # Constant + 5 __R1_local1 exp((-(Ea0+(phi*G_AB)))) # ConstantExpression + 6 __R1_local2 exp((-(Ea0+(phi*(G_AB+G_S1))))) # ConstantExpression + 7 __reverse__R1_local1 exp((-(Ea0+((1-phi)*(-G_AB))))) # ConstantExpression + 8 __reverse__R1_local2 exp((-(Ea0+((1-phi)*((-G_AB)+(-G_S1)))))) # ConstantExpression +end parameters +begin species + 1 A(b) 100 + 2 B(a,s~0) 100 + 3 B(a,s~1) 100 + 4 A(b!1).B(a!1,s~0) 0 + 5 A(b!1).B(a!1,s~1) 0 +end species +begin reactions + 1 1,2 4 __R1_local1 #_R1 + 2 1,3 5 __R1_local2 #_R1 + 3 4 1,2 __reverse__R1_local1 #_reverse__R1 + 4 5 1,3 __reverse__R1_local2 #_reverse__R1 +end reactions +begin groups + 1 AB0 4 + 2 AB1 5 +end groups diff --git a/tests/data/nfsim/arrhenius/state_pattern.xml b/tests/data/nfsim/arrhenius/state_pattern.xml new file mode 100644 index 000000000..e0d5704f4 --- /dev/null +++ b/tests/data/nfsim/arrhenius/state_pattern.xml @@ -0,0 +1,306 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +