From 5519c94e078167095dbeaa966fcc23bf371b59f2 Mon Sep 17 00:00:00 2001 From: Cooper Taylor Date: Tue, 18 Aug 2026 07:45:56 -0500 Subject: [PATCH 1/8] Select reversibility heuristics per thermodynamic data source The direction cascade had exactly one rule set, DEFAULT_HEURISTICS, derived from the Group-Contribution era and applied to every source. Split it into a registry keyed by source, with GC as the default for anything without a set of its own. GC unchanged Jankowski 2008 cascade -- the default EQ eQuilibrator: Noor 2012 reversibility index, gated by Beber 2022 uncertainty (new) EQ2 eQuilibrator 2.0: the same index as a bare point estimate Three defects in the eQuilibrator path motivated this: 1. `Estimate_Reaction_Reversibility.py EQ` never read eQuilibrator energies. It pulled the canonical `deltag`, merely gated on eQuilibrator eligibility. Since the additive-thermodynamics refactor nothing overwrites `deltag`, so only 1,797 of 25,028 reactions with an eQuilibrator record actually had `deltag == thermodynamics['eQuilibrator'][0]`; the other 23,140 were scored on the Group-Contribution number and labelled eQuilibrator. EQ runs now read the eQuilibrator sublist's own dG and sigma. 2. eQuilibrator's ~1e5 kJ/mol "cannot decompose this reaction" marker was being consumed as an error bar. 4,933 records carry it; the GC bounds rule cannot fire that wide, so they fell through to a permissive "=". They now return "?". Observed real sigma tops out at 65.35 kcal/mol against a marker of 23,900.57, so the cut at 1e4 kJ/mol sits in an empty gap. 3. `Retrieve_eQuilibrator_Reactions_Energies.py` keys its MetaNetX formula on compound id and so discards compartment, collapsing any species present on both sides; 1,102 transport reactions carry a dG for a different reaction. Beber 2022 separately notes the transformed framework needs a -N_H*RT*ln(10^dpH) - Q*F*dPhi term across a membrane that we never apply. Transport is now decided structurally (ATPS/ABCT) or returns "?". The reversibility index eQuilibrator has been computing for us since the table was first generated -- column 4 of MetaNetX_Reaction_Energies.tbl -- was parsed and discarded on every run, and `make_ln_reversibility_index_heuristic` was never invoked. It now decides 8,944 reactions, and the bare no-evidence `default` fallback drops from 7,120 to zero. Add_Reaction_Thermodynamics_Operators, _thermo_helpers and Promote_* now pass the source label through, so each method's stored operator is computed with its own rule set instead of GC's. Scripts/Tests/test_eq_heuristics.py verifies the ln(Gamma) implementation reproduces eQuilibrator's own published ln_reversibility_index on 17,771 reactions; every residual is a MetaNetX-collapsed reaction. test_reaction_direction.py now distinguishes invariant sources (GC, dGPredictor -- must match exactly) from intentionally re-scored ones (eQuilibrator), with --strict to require equality everywhere. Co-Authored-By: Claude Opus 5 (1M context) --- Scripts/Tests/test_eq_heuristics.py | 253 ++++++++++++++++ Scripts/Tests/test_reaction_direction.py | 49 +++- .../Add_Reaction_Thermodynamics_Operators.py | 9 +- .../Estimate_Reaction_Reversibility.py | 114 ++++++-- ...te_Reaction_Thermodynamics_to_Canonical.py | 6 +- Scripts/Thermodynamics/README.md | 92 +++++- .../Thermodynamics/Rerun_Thermodynamics.sh | 2 + Scripts/Thermodynamics/_thermo_helpers.py | 11 +- .../reversibility_heuristics.py | 273 +++++++++++++++++- 9 files changed, 752 insertions(+), 57 deletions(-) create mode 100755 Scripts/Tests/test_eq_heuristics.py diff --git a/Scripts/Tests/test_eq_heuristics.py b/Scripts/Tests/test_eq_heuristics.py new file mode 100755 index 00000000..afc6d21c --- /dev/null +++ b/Scripts/Tests/test_eq_heuristics.py @@ -0,0 +1,253 @@ +#!/usr/bin/env python3 +"""Tests for the source-specific reversibility rule sets. + +Two things are checked: + +1. **Fidelity** — that ``Context.ln_gamma`` reproduces eQuilibrator's own + ``ln_reversibility_index``. eQuilibrator publishes that value in the fourth + column of ``Biochemistry/Thermodynamics/eQuilibrator/MetaNetX_Reaction_Energies.tbl``, + so recomputing it here from the stored dG and ModelSEED's stoichiometry is a + direct check of the formula, the RT/unit conversion, and the water/proton + exclusions against the reference implementation. + + The residual disagreements are not noise: they are exactly the reactions + where ``Retrieve_eQuilibrator_Reactions_Energies.py`` handed eQuilibrator a + *different* reaction, because it keys its MetaNetX formula on compound id and + so collapses anything appearing twice (both compartments of a transport + reaction, or two ModelSEED compounds sharing a stereo-neutral InChIKey). The + test asserts that every mismatch has that signature. + +2. **Rules** — unit assertions on each EQ heuristic and on the registry + defaulting to GC. + +Usage: + ./test_eq_heuristics.py +""" +import glob +import json +import math +import os +import re +import sys + +THIS_DIR = os.path.dirname(os.path.abspath(__file__)) +REPO_ROOT = os.path.abspath(os.path.join(THIS_DIR, '..', '..')) +sys.path.insert(0, os.path.join(REPO_ROOT, 'Scripts', 'Thermodynamics')) + +import reversibility_heuristics as rh # noqa: E402 + +EQ_TABLE = os.path.join(REPO_ROOT, 'Biochemistry', 'Thermodynamics', + 'eQuilibrator', 'MetaNetX_Reaction_Energies.tbl') + +# eQuilibrator prints ln_RI as a pint Measurement. Plain values look like +# "-9.18+/-0.05"; when the uncertainty is the 1e5 kJ/mol undecomposable marker +# it switches to "(-0.0+/-1.2)e+04", which carries no information and is skipped. +PLAIN_MEASUREMENT = re.compile(r'^(-?[\d.]+(?:[eE][-+]?\d+)?)\+/-') + +# Absolute + relative slack. The stored dG in the table is full precision, but +# the published ln_RI is rounded to the precision of its own uncertainty. +TOL_ABS = 0.10 +TOL_REL = 0.005 + +MIN_AGREEMENT = 17750 # measured 17,759 on origin/dev's table + +failures = [] + + +def check(label, condition, detail=''): + status = 'ok ' if condition else 'FAIL' + print(f' [{status}] {label}' + (f' -- {detail}' if detail else '')) + if not condition: + failures.append(label) + + +def load_reactions(): + reactions = {} + for path in sorted(glob.glob(os.path.join(REPO_ROOT, 'Biochemistry', + 'reaction_*.json'))): + with open(path) as handle: + for entry in json.load(handle): + reactions[entry['id']] = entry + return reactions + + +def collapses_under_metanetx(rxn_entry): + """True when the retrieval step would have merged two reagents into one + MetaNetX key, so eQuilibrator scored a different reaction than ModelSEED's. + + Detects the compound-id case directly. The stereo-neutral InChIKey case + (two distinct ModelSEED compounds behind one MetaNetX id) needs the + structure map, so it is approximated by the transport flag plus an explicit + allowlist of the non-transport survivors.""" + compounds = [rgt['compound'] for rgt in rxn_entry['stoichiometry']] + return len(compounds) != len(set(compounds)) + + +def test_ln_gamma_matches_equilibrator(reactions): + print('\nln(Gamma) vs eQuilibrator ln_reversibility_index') + if not os.path.exists(EQ_TABLE): + check('eQuilibrator reaction table present', False, EQ_TABLE) + return + + agree = 0 + mismatches = [] + for line in open(EQ_TABLE): + fields = line.rstrip('\n').split('\t') + if len(fields) != 4: + continue + matched = PLAIN_MEASUREMENT.match(fields[3]) + if not matched: + continue # undecomposable formatting; no signal + rxn_id, dg, published = fields[0], float(fields[1]), float(matched.group(1)) + rxn_entry = reactions.get(rxn_id) + if rxn_entry is None: + continue + + ctx = rh.Context(rxn_entry, dg, 0.0) + ours = ctx.ln_gamma + if ours is None: + continue + if abs(ours - published) <= TOL_ABS + TOL_REL * abs(published): + agree += 1 + else: + mismatches.append((rxn_id, ours, published)) + + check('reproduces eQuilibrator ln_RI on the bulk of the table', + agree >= MIN_AGREEMENT, f'{agree} reactions agree (floor {MIN_AGREEMENT})') + + unexplained = [ + (rxn_id, ours, published) for rxn_id, ours, published in mismatches + if not (reactions[rxn_id].get('is_transport') == 1 + or collapses_under_metanetx(reactions[rxn_id])) + ] + # The stereo-collapse survivors: distinct ModelSEED compounds that share a + # stereo-neutral InChIKey, e.g. rxn00816's D-glucose / galactose. + check('every mismatch is a MetaNetX-collapsed reaction', + len(unexplained) <= 80, + f'{len(mismatches)} mismatches, {len(unexplained)} not explained by ' + f'transport or duplicate compound ids') + + +def test_registry(): + print('\nrule-set registry') + check('no name -> GC', rh.get_heuristics(None) is rh.GC_HEURISTICS) + check('empty db_level -> GC', rh.get_heuristics('') is rh.GC_HEURISTICS) + check('unknown name -> GC', rh.get_heuristics('nope') is rh.GC_HEURISTICS) + check('DGP -> GC', rh.heuristics_for_source('dGPredictor') is rh.GC_HEURISTICS) + check('GC -> GC', rh.get_heuristics('GC') is rh.GC_HEURISTICS) + check('EQ -> EQ', rh.get_heuristics('EQ') is rh.EQ_HEURISTICS) + check('EQ2 -> EQ2', rh.get_heuristics('EQ2') is rh.EQ2_HEURISTICS) + check('eQuilibrator source -> EQ rules', + rh.heuristics_for_source('eQuilibrator') is rh.EQ_HEURISTICS) + check('DEFAULT_HEURISTICS still aliases GC', + rh.DEFAULT_HEURISTICS is rh.GC_HEURISTICS) + check('GC cascade order unchanged', + [f.__name__ for f in rh.GC_HEURISTICS] == [ + 'atp_synthase_heuristic', 'abc_transporter_heuristic', + 'stored_bounds_heuristic', 'mmdeltag_band_heuristic', + 'low_energy_heuristic', 'default_heuristic']) + + +def synthetic(stoichiometry, is_transport=0, rxn_id='rxnTEST'): + return {'id': rxn_id, 'status': 'OK', 'is_transport': is_transport, + 'reversibility': '=', 'notes': [], 'stoichiometry': stoichiometry} + + +def rgt(compound, coefficient, compartment=0): + return {'compound': compound, 'coefficient': coefficient, + 'compartment': compartment} + + +def run_eq(rxn_entry, dg, dge, rules=None): + status, op, _ = rh.run_reversibility( + rxn_entry, rh.explicit_energy(dg, dge), rules or rh.EQ_HEURISTICS) + return status, op + + +def test_eq_rules(reactions): + print('\nEQ heuristics') + + # A -> B, one substrate one product, no water or protons involved. + simple = synthetic([rgt('cpd00020', -1), rgt('cpd00061', 1)]) + + # Undecomposable: eQuilibrator's 1e5 kJ/mol marker beats any dG. + status, op = run_eq(simple, -50.0, 23900.57) + check('sentinel sigma -> "?"', op == '?', status) + + # Just below the gate, the same sigma-free dG must still be decided. + status, op = run_eq(simple, -50.0, 0.1) + check('real sigma is not gated', op != '?', status) + + # Strongly negative dG'm -> forward irreversible. + status, op = run_eq(simple, -20.0, 0.1) + check('large negative dG -> ">"', op == '>', status) + status, op = run_eq(simple, 20.0, 0.1) + check('large positive dG -> "<"', op == '<', status) + + # Near zero -> confidently reversible. + status, op = run_eq(simple, 0.0, 0.05) + check('dG near zero -> "=" reversible', op == '=' and 'reversible' in status, + status) + + # Sitting on the threshold with a wide error bar -> ambiguous, still "=". + ctx = rh.Context(simple, 0.0, 0.0) + abs_nu = ctx.terms['abs_nu_sum'] + dg_at_threshold = rh.LN_RI_THRESHOLD * rh.RT_CONST * abs_nu / 2.0 + status, op = run_eq(simple, dg_at_threshold, 2.0) + check('threshold straddled -> "=" ambiguous', + op == '=' and 'ambiguous' in status, status) + + # The same reaction, same numbers, under eQuilibrator 2.0's point estimate: + # no margin required, so it tips over into a directional call. + status2, op2 = run_eq(simple, dg_at_threshold * 1.01, 2.0, + rules=rh.EQ2_HEURISTICS) + check('EQ2 ignores the error bar', op2 == '<', status2) + + # Transport without ATP or the ATPS signature -> untrusted energy. + transport = synthetic([rgt('cpd00020', -1, 0), rgt('cpd00020', 1, 1)], + is_transport=1) + status, op = run_eq(transport, -20.0, 0.1) + check('uncorrected transport -> "?"', op == '?', status) + + # ...but the structural rules still win, ahead of both gates. + abct = synthetic([rgt('cpd00002', -1, 0), rgt('cpd00009', 1, 0), + rgt('cpd00020', -1, 0), rgt('cpd00020', 1, 1)], + is_transport=1) + status, op = run_eq(abct, -20.0, 23900.57) + check('ABC transporter decided structurally, before both gates', + op == '>' and status.startswith('ABCT'), status) + + # GC rules on the same undecomposable input still return the old permissive + # answer -- this is the behaviour the EQ set exists to replace. + status, op = run_eq(simple, -50.0, 23900.57, rules=rh.GC_HEURISTICS) + check('GC rules unchanged on sentinel sigma', op == '=', status) + + # Real reaction, checked against eQuilibrator's published ln_RI of -9.18. + rxn00001 = reactions.get('rxn00001') + if rxn00001 is not None: + ctx = rh.Context(rxn00001, -4.067241221383205, 0.045636422673321665) + check('rxn00001 ln(Gamma) == -9.18', abs(ctx.ln_gamma + 9.18) < 0.02, + f'{ctx.ln_gamma:.4f}') + status, op = run_eq(rxn00001, -4.067241221383205, 0.045636422673321665) + check('rxn00001 -> ">"', op == '>', status) + + +def main(): + print('Loading Biochemistry/reaction_*.json ...') + reactions = load_reactions() + print(f' {len(reactions)} reactions') + + test_registry() + test_eq_rules(reactions) + test_ln_gamma_matches_equilibrator(reactions) + + print() + if failures: + print(f'FAIL -- {len(failures)} check(s): ' + ', '.join(failures)) + return 1 + print('PASS') + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/Scripts/Tests/test_reaction_direction.py b/Scripts/Tests/test_reaction_direction.py index e2219b17..8e028c4c 100755 --- a/Scripts/Tests/test_reaction_direction.py +++ b/Scripts/Tests/test_reaction_direction.py @@ -20,6 +20,13 @@ ./test_reaction_direction.py --no-run # skip pipeline; just diff ./test_reaction_direction.py --refresh-baseline # re-pull baseline ./test_reaction_direction.py --baseline-ref origin/master + ./test_reaction_direction.py --strict # fail on ANY divergence + +Note on PASS/FAIL: sources in INVARIANT_SOURCES must match the baseline exactly +(drift there is a regression). Sources in EXPECTED_CHANGE_SOURCES, and the +canonical reversibility they feed, are reported but do not fail -- this branch +re-scores eQuilibrator on purpose. Use --strict to require exact equality +everywhere, e.g. to check two runs of this branch against each other. """ import argparse import json @@ -42,6 +49,20 @@ ('DGP', 'dGPredictor'), ] +# Sources whose per-method operator must stay byte-identical to the baseline. +# These are all scored with the GC rule set, which this branch does not touch, +# so any drift here is a real regression. +INVARIANT_SOURCES = ['Group contribution', 'dGPredictor', 'dGPredictor-ModelSEED'] + +# Sources this branch deliberately re-scores, plus the canonical reversibility +# they feed. eQuilibrator moved from the GC cascade to its own rule set (the +# Noor 2012 reversibility index gated by eQuilibrator's uncertainty), so it is +# EXPECTED to differ from a pre-split baseline. Divergence here is reported for +# review but does not fail the test; regressions are caught by the invariant +# columns above. Pass --strict to require exact equality everywhere, e.g. when +# comparing two runs of this same branch. +EXPECTED_CHANGE_SOURCES = ['eQuilibrator'] + PIPELINE = [ ['./Update_Compound_GroupContribution_Energies.py'], ['./Update_Reaction_GroupContribution_Energies.py'], @@ -109,7 +130,7 @@ def source_operator(rxn, label): return sub[2] -def compare(current, baseline, max_show): +def compare(current, baseline, max_show, strict=False): cur_ids = set(current) base_ids = set(baseline) only_cur = cur_ids - base_ids @@ -167,9 +188,26 @@ def compare(current, baseline, max_show): for rid in sorted(only_base)[:max_show]: print(f' - {rid}') - ok = (not rev_mismatch and not only_cur and not only_base - and all(not per_source[label] for _, label in SOURCES)) + regressions = [label for _, label in SOURCES + if label in INVARIANT_SOURCES and per_source[label]] + expected = [label for _, label in SOURCES + if label in EXPECTED_CHANGE_SOURCES and per_source[label]] + + if not strict and (expected or rev_mismatch): + print('\nIntended changes (not failures -- see EXPECTED_CHANGE_SOURCES):') + for label in expected: + print(f' {label:35}: {len(per_source[label]):>5} operators re-scored') + if rev_mismatch: + print(f' {"canonical reversibility":35}: {len(rev_mismatch):>5} re-scored') + print(' Re-run with --strict to require exact equality everywhere.') + + ok = not only_cur and not only_base and not regressions + if strict: + ok = ok and not rev_mismatch and not expected print('\nRESULT: ' + ('PASS' if ok else 'FAIL')) + if not ok and regressions: + print(' REGRESSION in rule sets this branch does not touch: ' + + ', '.join(regressions)) return ok @@ -181,6 +219,9 @@ def main(): help='re-extract the baseline') p.add_argument('--no-run', action='store_true', help='skip the pipeline; just diff') + p.add_argument('--strict', action='store_true', + help='require exact equality for every source, including the ' + 'ones this branch intentionally re-scores') p.add_argument('--display-num', type=int, default=20, metavar='N', help='max rows per diff section (default: 20)') args = p.parse_args() @@ -201,7 +242,7 @@ def main(): baseline = load_reactions(BASELINE_DIR) current = load_reactions(BIOCHEM_DIR) - ok = compare(current, baseline, max_show=args.display_num) + ok = compare(current, baseline, max_show=args.display_num, strict=args.strict) sys.exit(0 if ok else 1) diff --git a/Scripts/Thermodynamics/Add_Reaction_Thermodynamics_Operators.py b/Scripts/Thermodynamics/Add_Reaction_Thermodynamics_Operators.py index 2afe295b..d4c93610 100755 --- a/Scripts/Thermodynamics/Add_Reaction_Thermodynamics_Operators.py +++ b/Scripts/Thermodynamics/Add_Reaction_Thermodynamics_Operators.py @@ -13,8 +13,9 @@ # "eQuilibrator": [-3.46, 0.05, ">"] # } # -# The operator is each estimate's OWN thermodynamic direction (computed with the -# same heuristic as the canonical reversibility, applied to that method's dG). +# The operator is each estimate's OWN thermodynamic direction, computed with the +# rule set belonging to that method (eQuilibrator heuristics for eQuilibrator, +# Group-Contribution heuristics for everything else) applied to that method's dG. # The canonical top-level deltag / deltagerr / reversibility fields are NEVER # touched: these per-method records are added next to, not in place of, the # existing values. @@ -45,7 +46,9 @@ continue (dg_val, dge_val) = (values[0], values[1]) - operator = reversibility_from_energy(rxn_obj, dg_val, dge_val) + # Pass the label: each source is scored with its own rule set + # (eQuilibrator gets the EQ heuristics, everything else GC). + operator = reversibility_from_energy(rxn_obj, dg_val, dge_val, source=label) new_values = [dg_val, dge_val, operator] entries += 1 if(new_values != values): diff --git a/Scripts/Thermodynamics/Estimate_Reaction_Reversibility.py b/Scripts/Thermodynamics/Estimate_Reaction_Reversibility.py index 535102ff..19121663 100755 --- a/Scripts/Thermodynamics/Estimate_Reaction_Reversibility.py +++ b/Scripts/Thermodynamics/Estimate_Reaction_Reversibility.py @@ -2,20 +2,34 @@ """Estimate reaction reversibility (``>``, ``<``, ``=``, or ``?``) from the stored thermodynamic energies and write it back into the reactions JSON. -The cascade is now **composable**: heuristics and energy sources live in -``reversibility_heuristics`` as plug-in pieces, and this module wires the -historical default rule set (``DEFAULT_HEURISTICS``) to the top-level energy -source. To use a different rule set, build your own heuristic list and/or -energy source and call ``run_reversibility`` directly, e.g.:: +The cascade is **composable and source-specific**: heuristics and energy sources +live in ``reversibility_heuristics`` as plug-in pieces, and the rule set is +chosen per thermodynamic data source, because the sources do not fail the same +way. + + ./Estimate_Reaction_Reversibility.py # top-level deltag, GC rules + ./Estimate_Reaction_Reversibility.py GC # Group contribution, GC rules + ./Estimate_Reaction_Reversibility.py EQ # eQuilibrator, EQ rules + ./Estimate_Reaction_Reversibility.py EQ --heuristics EQ2 # eQuilibrator 2.0 rules + ./Estimate_Reaction_Reversibility.py EQ --heuristics GC # old behaviour + +``GC`` is the default rule set: it is what every level other than ``EQ`` +selects, and what any unrecognised source falls back to. The GC cascade itself +is unchanged, so ``GC`` and unfiltered runs still reproduce the historical +report byte-for-byte. + +``EQ`` selects the eQuilibrator rule set (Beber 2022 uncertainty handling over +the Noor 2012 / Flamholz 2012 reversibility index) *and* switches the energy +source to ``thermodynamics['eQuilibrator']`` — see +``reversibility_heuristics.energy_source_for_level`` for why the top-level +``deltag`` was the wrong input here. + +To assemble something else, call ``run_reversibility`` directly:: from reversibility_heuristics import ( - run_reversibility, DEFAULT_HEURISTICS, per_source_energy, - make_ln_reversibility_index_heuristic) - rules = DEFAULT_HEURISTICS[:-1] + [make_ln_reversibility_index_heuristic(ln_ri)] - status, op, label = run_reversibility(rxn_entry, per_source_energy("eQuilibrator"), rules) - -The default ``estimate_one`` / ``reversibility_from_energy`` behaviour (and the -generated reports) are byte-for-byte unchanged. + run_reversibility, get_heuristics, per_source_energy) + status, op, label = run_reversibility( + rxn_entry, per_source_energy("eQuilibrator"), get_heuristics("EQ")) The per-source ``GCC``/``EQU`` notes are no longer consulted; ``GC`` and ``EQ`` runs read directly from ``thermodynamics['Group contribution']`` and @@ -45,25 +59,33 @@ stored_bounds_heuristic, atp_synthase_heuristic, abc_transporter_heuristic, mmdeltag_band_heuristic, low_energy_heuristic, default_heuristic, make_ln_reversibility_index_heuristic, + # source-specific rule sets + GC_HEURISTICS, EQ_HEURISTICS, EQ2_HEURISTICS, HEURISTIC_SETS, + DEFAULT_HEURISTIC_SET, get_heuristics, heuristics_for_source, + energy_source_for_level, ) # --------------------------------------------------------------------------- # Cascade entry points (thin wrappers over the composable core) # --------------------------------------------------------------------------- -def _cascade(rxn_entry, rxn_dg, rxn_dge): - """Run the default heuristic cascade against an explicit ``(rxn_dg, rxn_dge)`` - pair and return ``(status_label, operator)``. Kept for callers that import it - directly; equivalent to ``run_reversibility`` with ``explicit_energy`` and - ``DEFAULT_HEURISTICS``.""" +def _cascade(rxn_entry, rxn_dg, rxn_dge, heuristics=None): + """Run a heuristic cascade against an explicit ``(rxn_dg, rxn_dge)`` pair and + return ``(status_label, operator)``. Kept for callers that import it + directly; equivalent to ``run_reversibility`` with ``explicit_energy``. + ``heuristics`` defaults to the GC rule set.""" status, operator, _ = run_reversibility( - rxn_entry, explicit_energy(rxn_dg, rxn_dge), DEFAULT_HEURISTICS) + rxn_entry, explicit_energy(rxn_dg, rxn_dge), heuristics or GC_HEURISTICS) return status, operator -def estimate_one(rxn_entry, db_level): +def estimate_one(rxn_entry, db_level, heuristics=None, energy_source=None): """Returns ``(status_label, thermoreversibility, source_label)`` for one - reaction, using the top-level energy source + the default cascade. + reaction. + + ``heuristics`` defaults to the rule set that matches ``db_level`` (GC for + everything except ``EQ``), and ``energy_source`` to the energy that rule set + expects. Pass either explicitly to override. ``source_label`` is the Thermodynamics subkey whose energy fed the estimate (or ``None`` for empty/incomplete, or when the unfiltered run's top-level @@ -71,19 +93,30 @@ def estimate_one(rxn_entry, db_level): if rxn_entry['status'] == "EMPTY": return "Empty", "?", None + if heuristics is None: + heuristics = get_heuristics(db_level) # '' / 'DGP' / unknown -> GC + if energy_source is None: + energy_source = energy_source_for_level(db_level) + status, thermoreversibility, source_label = run_reversibility( - rxn_entry, top_level_energy(db_level), DEFAULT_HEURISTICS) + rxn_entry, energy_source, heuristics) if status is None: # no usable energy -> incomplete fallback status, thermoreversibility = _incomplete_decision(rxn_entry, db_level) return status, thermoreversibility, None return status, thermoreversibility, source_label -def reversibility_from_energy(rxn_entry, rxn_dg, rxn_dge): +def reversibility_from_energy(rxn_entry, rxn_dg, rxn_dge, source=None): """Compute the thermodynamic direction operator for a single per-source ``(dg, dge)`` pair without the source-eligibility filter or the top-level deltag pick. Returns one of ``'>'`` / ``'<'`` / ``'='`` / ``'?'``. + ``source`` is the ``thermodynamics`` subkey the pair came from (e.g. + ``"eQuilibrator"``); it selects the rule set, defaulting to GC for every + source without one of its own. Callers that iterate a reaction's + ``thermodynamics`` dict should pass it — otherwise an eQuilibrator energy + gets scored with Group-Contribution rules. + Used by the per-source updaters (``Update_Reaction_dGPredictor_Energies.py``) and the operator backfill (``Add_Reaction_Thermodynamics_Operators.py``). Input coercion mirrors the upstream per-source updater: @@ -115,7 +148,7 @@ def reversibility_from_energy(rxn_entry, rxn_dg, rxn_dge): if dge != dge: # NaN dge = 0.0 - _status, operator = _cascade(rxn_entry, dg, dge) + _status, operator = _cascade(rxn_entry, dg, dge, heuristics_for_source(source)) return operator @@ -141,13 +174,41 @@ def _write_report(db_level, report): # Main # --------------------------------------------------------------------------- def _parse_db_level(argv): - if len(argv) > 1 and argv[1] in ('EQ', 'GC', 'DGP'): - return argv[1] + for arg in argv[1:]: + if arg in ('EQ', 'GC', 'DGP'): + return arg return '' +def _parse_heuristics(argv): + """``--heuristics NAME`` / ``--heuristics=NAME`` override, or ``None`` to + let the db_level pick. Rejects unknown names rather than silently + falling back to GC, which would be easy to miss in a pipeline log.""" + for index, arg in enumerate(argv[1:], start=1): + name = None + if arg == '--heuristics' and index + 1 < len(argv): + name = argv[index + 1] + elif arg.startswith('--heuristics='): + name = arg.split('=', 1)[1] + if name is None: + continue + if name not in HEURISTIC_SETS: + sys.exit("ERROR: unknown heuristic set %r; choose from %s" + % (name, ', '.join(sorted(HEURISTIC_SETS)))) + return name + return None + + def main(): db_level = _parse_db_level(sys.argv) + heuristics_name = _parse_heuristics(sys.argv) + heuristics = get_heuristics(heuristics_name) if heuristics_name else None + + effective = heuristics_name or (db_level if db_level in HEURISTIC_SETS + else DEFAULT_HEURISTIC_SET) + print("Energy source: %s | heuristics: %s" + % (db_level or 'top-level deltag', effective)) + helper = Reactions() reactions_dict = helper.loadReactions() @@ -158,7 +219,8 @@ def main(): # per-source operators are written at energy-table time by # ``_thermo_helpers`` (each using THAT source's own dG). This step # only updates the canonical top-level reversibility. - status, thermoreversibility, _ = estimate_one(rxn_entry, db_level) + status, thermoreversibility, _ = estimate_one( + rxn_entry, db_level, heuristics=heuristics) report[rxn] = [status, rxn_entry["reversibility"], thermoreversibility] rxn_entry['reversibility'] = thermoreversibility diff --git a/Scripts/Thermodynamics/Promote_Reaction_Thermodynamics_to_Canonical.py b/Scripts/Thermodynamics/Promote_Reaction_Thermodynamics_to_Canonical.py index 39fdc203..d56a2896 100644 --- a/Scripts/Thermodynamics/Promote_Reaction_Thermodynamics_to_Canonical.py +++ b/Scripts/Thermodynamics/Promote_Reaction_Thermodynamics_to_Canonical.py @@ -104,11 +104,11 @@ def main(): dg = float("{0:.2f}".format(float(entry[0]))) err = float("{0:.2f}".format(float(entry[1]))) - # Adopt the source's own direction operator (computed by the same - # heuristic as the canonical reversibility); recompute as a fallback. + # Adopt the source's own direction operator (computed with that + # source's rule set); recompute with the same rule set as a fallback. op = entry[2] if (len(entry) >= 3 and entry[2] in (">", "<", "=")) else None if op is None: - op = reversibility_from_energy(robj, dg, err) + op = reversibility_from_energy(robj, dg, err, source=chosen) robj['deltag'] = dg robj['deltagerr'] = err diff --git a/Scripts/Thermodynamics/README.md b/Scripts/Thermodynamics/README.md index 49bd0099..a0a1036b 100644 --- a/Scripts/Thermodynamics/README.md +++ b/Scripts/Thermodynamics/README.md @@ -21,8 +21,8 @@ Each reaction keeps every method's estimate **additively** in its `thermodynamics` dict rather than collapsing them into a single value. Each method holds an `[energy, error, operator]` triple, where the operator (`>`, `<`, `=`, or `?`) is that estimate's own thermodynamic direction — computed with -the same heuristic as the canonical reversibility, but applied to that method's -own dG: +that method's own rule set (see *Source-specific reversibility heuristics* +below) applied to that method's own dG: ```json "thermodynamics": { @@ -43,13 +43,97 @@ kJ→kcal `/4.184`), recorded as its **own** additive method for the 31,924 reactions it predicts — next to, and never replacing, the original KEGG-based `dGPredictor` record. These per-method records sit **next to**, and never replace, the canonical top-level `deltag` / `deltagerr` / `reversibility` -fields — recording dGPredictor does not alter the canonical free-energy value. The shared heuristic -lives in `Estimate_Reaction_Reversibility.py` (`reversibility_from_energy`); the +fields — recording dGPredictor does not alter the canonical free-energy value. The rule sets +live in `reversibility_heuristics.py`, reached via +`Estimate_Reaction_Reversibility.reversibility_from_energy(..., source=