Add frozen, reproducible KEGG reaction benchmark (Phase 1)

.gitattributes

@@ -0,0 +1,7 @@
+# Frozen benchmark artifacts carry SHA-256 digests in benchmark/data/VERSION.json.
+# Those digests are computed over LF-terminated files, so line endings must survive
+# checkout unchanged on every platform or the digests will not verify.
+benchmark/data/*.csv text eol=lf
+benchmark/data/*.json text eol=lf
+benchmark/manifest/*.txt text eol=lf
+benchmark/manifest/*.json text eol=lf

.gitignore

@@ -29,4 +29,11 @@ analysis/**
 # testing files
 **/scratch.ipynb
 **/participants_likelihood*.csv
-**/reaction_likelihood*.csv
+**/reaction_likelihood*.csv
+
+# benchmark artifacts
+# Downloaded SBML is reproducible from the manifest + checksums, so it is not vendored.
+benchmark/models/
+# Transient build scratch and logs (leading underscore = not a frozen artifact).
+benchmark/data/_*
+**/BioModels_251106

benchmark/README.md

@@ -0,0 +1,148 @@
+# Reaction Annotation Benchmark (Phase 1)
+
+Frozen, checksum-verified benchmark for evaluating KEGG reaction annotation under
+imperfect metabolite identifiers.
+
+Current version: **`phase1-v1`** — 74 models, 5,838 ground-truth reactions,
+5,816 evaluable. See [`data/RECONCILIATION.md`](data/RECONCILIATION.md) for the
+observed counts and why they differ from the historical 68 / 4,379 report.
+
+## Layout
+
+```
+benchmark/
+  manifest/
+    models.txt            # 75 BioModels accessions (canonical manifest)
+    model_paths.txt       # Repo-relative paths for evaluation scripts
+    model_registry.json   # Generated: URLs, upstream + local SHA-256, revisions
+  models/                 # Downloaded SBML (gitignored, reproducible)
+  data/                   # Frozen benchmark tables (generated)
+  scripts/
+    download_biomodels.py
+    build_benchmark.py
+```
+
+## Quick start
+
+```bash
+python benchmark/scripts/download_biomodels.py   # resolve, download, verify
+python benchmark/scripts/build_benchmark.py      # build frozen tables
+python -m pytest tests/test_benchmark_build.py   # regression tests
+```
+
+Add `--with-candidates` to attach rule-based candidates and heuristic scores.
+This is slower and requires the KEGG reference data under `data/kegg/`.
+
+## Provenance
+
+Main SBML filenames are **not** predictable from the accession
+(`BIOMD0000000013_url.xml` vs `iJN1463.xml`), so the downloader resolves each
+model's file list from the BioModels API rather than guessing. Every file is
+verified against the upstream SHA-256 published by the API; `model_registry.json`
+records the download URL, upstream and local digests, SBML format version,
+curation status, and the latest upstream revision.
+
+## Pipeline failures are not scientific exclusions
+
+This distinction is enforced in code and in the outputs.
+
+**Pipeline failures** (`pipeline_failures.csv`) indicate a defect in our tooling
+or in the snapshot. They must be resolved or explicitly justified, never folded
+into exclusion statistics.
+
+| Failure type | Meaning |
+|---|---|
+| `file_missing` | Download absent — a tooling failure, not a property of the data |
+| `parse_error` | libSBML cannot parse the file; examined individually |
+| `checksum_unverified` | Local bytes disagree with the upstream digest |
+| `provenance_missing` | Registry has no entry for an accession |
+| `registry_missing` | No registry at all, so nothing is provenance-verified |
+
+**Scientific exclusions** (`exclusions.csv`) are properties of the data.
+
+| Level | Reason | Effect |
+|---|---|---|
+| model | `no_kegg_reaction_annotations` | Model carries no reaction-level ground truth |
+| reaction | `exchange_ssx` | Empty reactant or product side; record kept, marked not evaluable |
+| reaction | `invalid_ground_truth_id` | Identifier fails `R#####` validation |
+
+SSX exclusions reduce **reaction** counts, never model counts. A model whose
+reactions are all SSX stays in the benchmark and contributes zero evaluable
+reactions; that case is reported as `models_with_zero_eval_reactions`.
+
+## Parser problems are surfaced, not absorbed
+
+`parser_diagnostics.csv` compares, per model, the number of reactions whose raw
+annotation mentions `kegg.reaction` against the number from which an identifier
+was actually extracted. Any shortfall sets `parser_discrepancy` — the signature
+of a URI-form or parser defect rather than a data property.
+
+It also counts species carrying `kegg.reaction` annotations. The manifest was
+assembled by grepping raw file bytes, so models are selected as "KEGG-annotated"
+even when the only mentions are misplaced species annotations
+(`misplaced_annotations_only`). `BIOMD0000000579` is exactly this case.
+
+## Ground truth policy
+
+**All** valid KEGG reaction IDs are preserved in `ground_truth_kegg_all`;
+`ground_truth_kegg_primary` is the first for backward-compatible scoring. 91
+reactions carry more than one valid ID, one of them 10. Scoring against a single
+arbitrary ID systematically under-credits correct predictions.
+
+## Train/test partitioning
+
+Split on `cluster_id` from `model_clusters.csv`, never on individual reactions.
+Clusters group near-duplicate model variants by ground-truth overlap (Jaccard, or
+containment gated on comparable size and a minimum shared-identifier count) so
+closely related models always land in the same partition.
+
+Cluster IDs are `CLU_<smallest member accession>` and stay stable across rebuilds
+while membership is unchanged. Every included model has a cluster ID, singletons
+included: 74 models across 54 clusters, 9 clusters holding more than one model.
+
+## Determinism and versioning
+
+Tables are sorted on stable keys and written with a fixed line terminator; no
+timestamps enter the data tables. Two consecutive builds produce byte-identical
+output across all 11 artifacts. `VERSION.json` freezes the version with SHA-256
+digests of every artifact plus the manifest and registry, and records whether all
+invariants passed.
+
+## Dataset invariants
+
+`invariants.json` records ten checks, all currently passing:
+
+1. Manifest holds 75 unique accessions
+2. Every registry entry carries a download URL and local SHA-256
+3. Every downloaded file matches the upstream SHA-256
+4. Every included model parses
+5. Every evaluable reaction has at least one valid ground-truth ID
+6. All included ground-truth IDs are well formed
+7. Reaction records reconcile: total = evaluable + reaction-level exclusions
+8. Model records reconcile: manifest = included + excluded + pipeline failures
+9. Every included model has a cluster ID
+10. `(model_id, reaction_id)` keys are unique
+
+## Outputs
+
+| File | Contents |
+|---|---|
+| `reactions.csv` | One row per (model, reaction) with ground truth and flags |
+| `model_context.csv` | Model title and notes; join on `model_id` |
+| `model_summary.csv` | Per-model status, counts, cluster assignment |
+| `model_clusters.csv` | `model_id` → `cluster_id` for partitioning |
+| `exclusions.csv` | Scientific exclusions only |
+| `pipeline_failures.csv` | Tooling and provenance failures |
+| `duplicate_groups.csv` | Multi-member clusters with linkage rules |
+| `species_annotations.csv` | Species-level annotation inventory |
+| `parser_diagnostics.csv` | Raw-vs-extracted KEGG mention comparison |
+| `benchmark_summary.json` | Observed counts |
+| `invariants.json` | Invariant check results |
+| `VERSION.json` | Frozen version with artifact digests |
+| `RECONCILIATION.md` | Discrepancy analysis vs the historical report |
+
+## Branch integration
+
+Phase 1 integrates evaluation code from `upstream/test/no-rule-evals`, including
+the KEGG URI slash-form fix that `main` lacks, KEGG-compound species support, SSX
+detection, and the ontology/cofactor ablation toggles.

benchmark/data/RECONCILIATION.md

@@ -0,0 +1,179 @@
+# Reconciliation: observed counts vs the historical 68 models / 4,379 reactions
+
+Benchmark version `phase1-v1`. Built from a checksum-verified snapshot of all 75
+manifest accessions. Every number below is read from the generated artifacts,
+not assumed.
+
+## Observed counts
+
+| Quantity | Observed |
+|---|---|
+| Manifest accessions (unique) | 75 |
+| Models downloaded and SHA-256 verified | 75 / 75 |
+| Models parsed by libSBML | 75 / 75 |
+| Pipeline failures | 0 |
+| Models included in benchmark | 74 |
+| Models scientifically excluded | 1 |
+| Ground-truth reactions (all records) | 5,838 |
+| Evaluable reactions | 5,816 |
+| Reactions excluded, exchange/SSX | 19 |
+| Reactions excluded, malformed ground-truth ID | 3 |
+| Reactions carrying more than one KEGG ID | 91 |
+| Distinct KEGG reaction identifiers | 1,906 |
+| Duplicate clusters with more than one member | 9 (29 models) |
+| Distinct clusters (partition units) | 54 |
+
+`reactions.csv` SHA-256: `f11ebc6f4d0734deb14c054caf277b81695b22668017d70af9c43cd6146004f1`
+
+All ten dataset invariants pass; see `invariants.json`.
+
+## The single scientific exclusion
+
+**BIOMD0000000579** carries no reaction-level KEGG annotation. It was
+investigated individually rather than accepted at face value, because a
+"KEGG-annotated" model yielding zero ground truth is exactly the signature of a
+parser defect.
+
+Finding: the model has 175 reactions, **all** of which carry annotations, but
+none reference `kegg.reaction`. The 16 raw `kegg.reaction` mentions in the file
+belong to **13 species**, for example species `s308` annotated with
+`kegg.reaction/R00289`. Annotating a metabolite with a reaction identifier is a
+curation error in the source model.
+
+This is a genuine data-level exclusion, not a tooling failure. It is recorded in
+`parser_diagnostics.csv` under `misplaced_annotations_only`, and it exposes a
+manifest-construction weakness: the original `find_kegg_in_biomodels.py` selected
+models by grepping raw file bytes for `kegg.reaction`, which cannot distinguish a
+reaction annotation from a misplaced species annotation. The manifest therefore
+over-selects.
+
+## Why the counts differ from 68 / 4,379
+
+The historical figures are treated as a clue. Exclusions were **not** adjusted to
+reproduce them. Four candidate explanations were tested.
+
+### 1. Upstream file drift — ruled out
+
+Every one of the 75 models has its latest upstream revision dated **before** the
+2025-11-06 snapshot date, so no model has changed at BioModels since the
+historical run. The corpus we built is the same set of model versions the
+historical run would have seen.
+
+One caveat worth recording: the copy of `BIOMD0000000013` previously vendored in
+`tests/test_models/` is 99,166 bytes, whereas the verified upstream main file is
+96,608 bytes. The historical `tests/BioModels_251106/` directory is absent from
+the repository and from this machine, so its byte-level contents cannot be
+compared. Local snapshots obtained by bulk download were evidently not identical
+to the API's `main` SBML file, which is precisely why this benchmark now pins
+upstream checksums.
+
+### 2. A colon-only KEGG URI pattern — a real defect, and the likely reason this
+was hard to reproduce
+
+`main` matched only `kegg.reaction:R00024`. Real BioModels files overwhelmingly
+use the slash form `kegg.reaction/R00024`, so `main`'s pattern extracts **zero**
+ground-truth reactions from this corpus. The unmerged evaluation branch carried
+the fix; `main` did not. The fix is now on this branch and locked down by
+parametrised regression tests over slash, colon, and URN forms.
+
+This means the historical numbers could not have been produced by `main`'s
+extraction code, and any attempt to reproduce them from `main` would have yielded
+nothing at all.
+
+### 3. Ground-truth policy — affects scoring, not counts
+
+The historical helper kept only the first KEGG ID per reaction. That changes
+measured accuracy but not the number of reaction records, so it cannot account
+for the 1,459-reaction gap. It does matter scientifically: **91 reactions** carry
+multiple valid IDs, one of them (`BIOMD0000000015`, reaction `den`) carrying 10.
+Scoring against a single arbitrary ID under-credits correct predictions on those
+reactions.
+
+### 4. Which models were actually evaluated — the most plausible cause
+
+The corpus is extremely top-heavy. Seven genome-scale models contribute **4,556
+of 5,838** ground-truth reactions (78%); the other 67 models contribute 1,282.
+Median reactions per model is 10; the largest single model contributes 1,047.
+Any difference in how the largest models were handled swings the total enormously:
+
+| Scenario | Models | Reactions |
+|---|---|---|
+| This build (all included) | 74 | 5,838 |
+| Excluding the 2 genome-scale models with no species annotations | 72 | 4,418 |
+| Excluding all 7 genome-scale models | 67 | 1,282 |
+| Historical report | 68 | 4,379 |
+
+No subset rule tested reproduces 68 / 4,379 exactly. The nearest simple rule —
+dropping the two genome-scale models that have no usable species annotations
+(`BIOMD0000001090`, `BIOMD0000001091`, 1,420 reactions between them) — yields
+4,418 reactions across 72 models, still not a match.
+
+Two concrete pieces of evidence point at the model list rather than the rules:
+
+- The committed evaluation script defaulted to
+  `tests/kegg_annotated_files-test.txt`, a file that **does not exist** in the
+  repository, rather than the 75-entry `kegg_annotated_files.txt`. The published
+  numbers therefore likely came from a different, probably truncated, list.
+- The historical outputs (`per_reaction_results.csv`, `results_summary.csv`) were
+  never committed, so the excluded set cannot be recovered directly.
+
+**Conclusion.** The difference is attributable to which models were fed to the
+historical run, not to the corpus having changed and not to the exclusion rules
+applied here. 68 / 4,379 is not reproducible from the 75-model manifest under any
+defensible rule, and it is not adopted as an acceptance criterion. The current
+counts are fully reconciled and reproducible: 75 = 74 included + 1 excluded + 0
+pipeline failures, and 5,838 = 5,816 evaluable + 22 reaction-level exclusions.
+
+## Reaction-level exclusions in detail
+
+**Exchange/SSX (19 reactions across 9 models).** Reactions with an empty
+reactant or product side, such as `PRPP =>` in `BIOMD0000000015`. The rule-based
+matcher generates no candidates for these when `include_exchange_reactions=False`,
+so they are retained as records and marked `included_in_eval=False`. They reduce
+reaction counts only; no model is excluded on this basis, and
+`models_with_zero_eval_reactions` is 0.
+
+**Malformed ground-truth ID (3 reactions).** Reaction `gluconeogenesis_ser` in
+`BIOMD0000000268`, `BIOMD0000000450`, and `BIOMD0000000674` is annotated
+`R0006565`, which has seven digits where KEGG uses five. A curation typo in three
+related models, all of which fall in cluster `CLU_BIOMD0000000268`.
+
+## Duplicate clusters and partitioning
+
+Clustering links models by ground-truth overlap using two rules, both recorded
+per cluster in `duplicate_groups.csv`:
+
+- **Jaccard** ≥ 0.9 catches same-scope variants, such as
+  `Yamada2003_JAK_STAT_pathway` and its SOCS1 knockout.
+- **Containment** ≥ 0.9 catches models that extend another, gated on comparable
+  set size (ratio ≥ 0.5) and at least 5 shared identifiers. Without containment
+  the three Smallbone2013 yeast variants split apart: `BIOMD0000000473` shares
+  189 reactions with `BIOMD0000000471`/`472` yet scores only 0.86 by Jaccard
+  because it adds 23 more. Without the size and overlap gates, clustering
+  collapses — a two-reaction model is trivially contained in a genome-scale one,
+  which chained unrelated models into one cluster and cut distinct clusters from
+  65 to 29. Both failure modes are covered by regression tests.
+
+Cluster IDs are `CLU_<lexicographically smallest member>`, so they are stable
+across rebuilds for unchanged membership. Every included model has a cluster ID,
+singletons included, so train/test splits partition on `cluster_id` and never on
+individual reactions.
+
+One judgment call is worth flagging. `CLU_BIOMD0000000042` groups seven yeast
+glycolysis models from seven different papers (Nielsen1998, Hynne2001,
+Galazzo1990, Teusink2000, Bakker2001, Albert2005, Ralser2007). They are not
+variants of one another in provenance terms, but they annotate substantially the
+same KEGG reactions, so separating them across a split would leak labels. The
+grouping is deliberately conservative; thresholds are exposed as
+`--duplicate-threshold`, `--containment-threshold` if a stricter provenance-based
+definition is preferred later.
+
+## Reproducing this build
+
+```bash
+python benchmark/scripts/download_biomodels.py   # 75/75 checksum-verified
+python benchmark/scripts/build_benchmark.py      # deterministic artifacts
+python -m pytest tests/test_benchmark_build.py   # 21 regression tests
+```
+
+Two consecutive builds produced byte-identical output across all 11 artifacts.

benchmark/data/VERSION.json

@@ -0,0 +1,24 @@
+{
+  "benchmark_version": "phase1-v1",
+  "manifest_sha256": "3fd01d100bf544af9217722b5fe4a7204cf82aea9ec391571eae11965e26a861",
+  "registry_sha256": "8ba6de7b4171309828a043fa83835dcb1c47d202952152fa945119fa8cb1e367",
+  "artifact_sha256": {
+    "reactions.csv": "f11ebc6f4d0734deb14c054caf277b81695b22668017d70af9c43cd6146004f1",
+    "model_summary.csv": "90210acf3d2e05c1a6b94f38aee25c4c0663b96383d7fdb5fdf32392a7f2f3d8",
+    "model_context.csv": "94ea74e1506183b5b1ebd231d5cf24b33824e71b6a5188c9f2b588ff16cd9f1d",
+    "model_clusters.csv": "9878d843888077e4b3ac8110260492bdbb81491f2d8197eca702ee35ce8e48fb",
+    "exclusions.csv": "d71b0614016e3c0693d6f6d860af41ea4f9d582aad3dbd697124d699afdaafb5",
+    "pipeline_failures.csv": "01ba4719c80b6fe911b091a7c05124b64eeece964e09c058ef8f9805daca546b",
+    "duplicate_groups.csv": "56e62170cee8fffe2a0270a7130cd83db43e188e9bba6b403582286a10e67460",
+    "species_annotations.csv": "99320d1da098215211d0364b70e719a6e59383b82383a30e3436f96d8abafd82",
+    "parser_diagnostics.csv": "889e72069702858817a07c275c68cb1640bf70e2ef362ea202a931a4afe56b06",
+    "benchmark_summary.json": "bcdb4b0ef07d1735e6ba31a221e8cf11a7bd1c04f8f40c0e6f8154c0c903ff5d",
+    "invariants.json": "135665dbbeb036816e18e910918ca87a664463c4a5b1b12708fad1b3b6e89918"
+  },
+  "counts": {
+    "models_included": 74,
+    "total_ground_truth_reactions": 5838,
+    "evaluable_reactions": 5816
+  },
+  "invariants_all_passed": true
+}