Benchmark assign_compartments transport pruning and gap-filling

CHANGELOG.md

@@ -4,6 +4,19 @@ Milestones in the raven-toolbox port. For function-level status see
 [docs/raven_migration.md](https://github.com/SysBioChalmers/raven-toolbox/blob/develop/docs/reference/migration.md); for open work see
 [docs/todo.md](https://github.com/SysBioChalmers/raven-toolbox/blob/develop/docs/reference/todo.md).
 
+## Unreleased
+
+* **Benchmarked two previously-unmeasured `assign_compartments` features.** A new ablation study
+  ([`assignment_ablations.md`](docs/studies/assignment_ablations.md),
+  `scripts/benchmark_assignment_ablations.py`) measures transport pruning and gap-filling on yeast-GEM,
+  which were shipped with correctness tests but no performance numbers. Transport pruning removes 21.6 %
+  of provisioned transports (a strict subset) at no reaction-accuracy cost. Gap-filling never fires on a
+  well-connected draft (no gratuitous additions at genome scale), and under ground-truthed
+  knockout-recovery it re-adds the *exact* removed reaction in 100 % of recoveries with zero wrong
+  additions — its recall (~45 %) is bounded by `cobra.flux_analysis.gapfill`'s numerical tolerance (a
+  sharp optimum at cobra's default, confirmed by an `integer_threshold` sweep), and it fails safe by
+  declining rather than mis-filling.
+
 ## 0.3.0 — 2026-07-16
 
 Compartment localisation and per-reaction confidence tracking, new gap-filling and flux-sampling

docs/studies/assignment_ablations.md

@@ -0,0 +1,101 @@
+# Compartment-assignment ablations: transport pruning & gap-filling
+
+Two `assign_compartments` features shipped with correctness tests but no performance measurement. This
+study measures each on curated *S. cerevisiae* yeast-GEM, flattened to one compartment and reassigned —
+the same draft the other yeast studies use (2569 draft reactions, 2296 to place, `min_growth` = 50 % of
+curated). Regenerate with `scripts/benchmark_assignment_ablations.py`.
+
+## 1. Transport-reaction pruning (`prune_transports`)
+
+After placement, `assign_compartments` adds the inter-compartment transports the network needs to stay
+functional. Pruning then removes the ones that turned out redundant — a transport whose removal leaves the
+model still certified. The question: how much does it remove, and does removing it cost accuracy?
+
+| | transports added | reaction agreement vs curated | runtime |
+|---|--:|--:|--:|
+| `prune_transports=True` (default) | **994** | 0.721 | 209 s |
+| `prune_transports=False` | 1268 | 0.721 | 74 s |
+
+- **Pruning removes 274 transports — 21.6 % — at no accuracy cost.** Reaction-level agreement is identical
+  to four decimals, and the pruned set is a **strict subset** of the unpruned one (0 transports appear only
+  when pruning is on). So pruning only ever removes genuinely redundant transports; it never trades one
+  placement for another.
+- **The cost is runtime**, not quality: pruning re-certifies the model after each candidate removal, so the
+  default is ~3× slower. That is the price of a leaner network — 274 fewer artificial inter-organelle
+  shuttles that a curator would otherwise have to inspect.
+
+This is the intra-method complement to the [CarveFungi head-to-head](carvefungi_milp_benchmark.md), which
+showed the *transport-minimisation objective term* yields ~41 % fewer transports than CarveFungi's carve;
+here the *post-hoc pruning pass* removes a further fifth of what placement provisionally added.
+
+## 2. Gap-filling on the natural draft (does it over-add?)
+
+`assign_compartments` can pull reactions from a `universal` model to restore biomass when a
+compartmentalised placement cannot grow. The first question is whether it fires *gratuitously* — adds
+reactions a well-connected model did not actually need.
+
+Reassigning the flattened draft with a universal available (the draft itself) vs without:
+
+| | added reactions | certified | growth |
+|---|--:|:--:|--:|
+| no `universal` | 0 | yes | 0.14 |
+| `universal=` draft | **0** | yes | 0.14 |
+
+**Gap-fill adds nothing on the natural draft.** Transport addition alone restores growth well above the
+floor (0.14 vs the 0.04 floor), so the growth-failure feedback that triggers gap-fill never fires. This is
+the genome-scale confirmation of the `test_no_gratuitous_gapfill` unit test: the feature is a safety net
+that stays inert when the model is already functional, rather than a source of spurious additions.
+
+## 3. Gap-filling under real gaps (does it add the *right* reaction?)
+
+To measure gap-fill *working*, it needs real gaps. Ground-truthed knockout-recovery on the certified
+compartmentalised model: take each of the 352 internal reactions whose single removal drops growth below
+5 % of optimum, remove it, and gap-fill from a universal that contains it (a copy of the model — so the
+removed reaction is a candidate and the ground truth is exact). This is the same
+`cobra.flux_analysis.gapfill` call `assign_compartments._gapfill` wraps. A 60-reaction seeded sample:
+
+| outcome | count | of sample |
+|---|--:|--:|
+| recovered growth, **re-added the exact removed reaction** | 27 | 45 % |
+| `cobra.gapfill` numerical failure (declined) | 33 | 55 % |
+| recovered growth with a *wrong* reaction | **0** | 0 % |
+
+Two things stand out, and they are the honest headline:
+
+- **Precision is perfect. Of every recovery, 100 % re-added the exact reaction that was removed** — and
+  zero knockouts were "fixed" with a different or superfluous reaction. When gap-fill answers, it answers
+  correctly.
+- **The only failure mode is declining to answer.** The 55 % shortfall is entirely
+  `cobra.flux_analysis.gapfill`'s own numerical-validation limit ("Failed to validate gap filled model, try
+  lowering the integer threshold"), not a wrong addition. `_gapfill` catches it and reports no gap-fill, so
+  a placement that cobra cannot solve stays *uncertified and visible* rather than silently mis-filled. The
+  recovery rate is therefore a property of cobra's MILP tolerance, not of the assignment logic — and it is
+  bounded below the true recoverable set, never above it.
+
+cobra's failure message suggests "try lowering the integer threshold", so we swept it. It does not help
+— the default is a **sharp optimum**:
+
+| `integer_threshold` | recovered (of 30) | exact |
+|---|--:|--:|
+| 1e-9 | 0 % | — |
+| **1e-6 (cobra default)** | **~47 %** | **100 %** |
+| 1e-5 / 1e-4 / 1e-3 | 0 % | — |
+
+Every value other than the default collapses recovery to zero, in both directions. So there is no tuning
+win to be had, the documented remedy is misleading, and `_gapfill` correctly leaves the threshold at
+cobra's default. The ~45–47 % ceiling is cobra's, and it is not cheaply liftable — which makes the
+fail-safe design (decline, never mis-fill) the right call rather than a workaround.
+
+## 4. Verdict
+
+- **Transport pruning:** a strict win — 21.6 % fewer transports at zero accuracy cost, paid for in runtime.
+  Correctly a default.
+- **Gap-filling:** precise and conservative. It never fires gratuitously (§2) and never adds a wrong
+  reaction (§3); when it recovers a gap it recovers it *exactly*. Its recall is capped by cobra's gapfill
+  tolerance, and it fails safe — an unsolved gap is reported, not papered over.
+
+## 5. Reproducing
+
+```
+python scripts/benchmark_assignment_ablations.py   # all three parts -> .research_tmp/assignment_ablations.json
+```

docs/studies/index.md

@@ -21,6 +21,10 @@ its equivalence claims against MATLAB RAVEN.
 - **[Yeast-GEM validation](yeast_validation.md)** — `assign_compartments` on curated yeast-GEM:
   recovery of the compartmentalisation, a CarveFungi head-to-head (McNemar), and a biological check
   (transporter connectivity, pathway localisation, dual-localised enzymes).
+- **[Compartment-assignment ablations](assignment_ablations.md)** — measuring two `assign_compartments`
+  features that were only correctness-tested: transport pruning (21.6 % fewer transports at no accuracy
+  cost) and gap-filling (never gratuitous; 100 % exact when it recovers, capped by cobra's gapfill
+  tolerance).
 - **[Multi-organism validation](multiorganism_validation.md)** — the same method across four kingdoms
   (yeast, Human-GEM, AraCore, iCre1355), including the chloroplast, with no per-organism changes.
 - **[Confidence tracking](confidence_tracking.md)** — per-reaction, multi-facet confidence persisted in
@@ -72,6 +76,7 @@ kegg_hmm_cutoff_calibration
 localization_redesign
 curation_priority_signals
 yeast_validation
+assignment_ablations
 multiorganism_validation
 confidence_tracking
 ```

scripts/benchmark_assignment_ablations.py

@@ -0,0 +1,221 @@
+#!/usr/bin/env python3
+"""Ablation benchmarks for two `assign_compartments` features that were previously only
+correctness-tested, never measured on a real model: **transport-reaction pruning** and **gap-filling**.
+
+Run on curated *S. cerevisiae* yeast-GEM, flattened to one compartment and reassigned (the same draft the
+other yeast studies use). Three parts, one JSON output:
+
+1. **Transport pruning** (`prune_transports=True` vs `False`). How many transport reactions pruning
+   removes, whether the pruned set is a strict subset of the unpruned one, whether reaction-level
+   placement accuracy moves, and the runtime cost.
+
+2. **Gap-fill, natural draft** (`universal=` passed vs not, no induced gaps). Whether gap-fill fires at
+   all when the flattened draft is simply reassigned — i.e. does the feature add reactions gratuitously
+   at genome scale, or does transport addition alone restore growth?
+
+3. **Gap-fill, knockout-recovery** (ground-truthed). Remove known growth-essential reactions from the
+   certified model one at a time and gap-fill from a universal that contains them (the same
+   `cobra.flux_analysis.gapfill` call `assign_compartments` uses internally). Of the reactions gap-fill
+   adds, how many restore growth, and how many are the *exact* removed reaction — the precision of the
+   additions against ground truth. cobra's gapfill has a known numerical-validation failure mode; its
+   rate is reported honestly, since `assign_compartments` degrades to "no gap-fill" when it trips.
+
+ASCII-only output. Deterministic (seeded knockout sample).
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+import time
+import warnings
+from pathlib import Path
+
+warnings.filterwarnings("ignore")
+import cobra  # noqa: E402
+from cobra.flux_analysis import gapfill as cobra_gapfill  # noqa: E402
+from cobra.flux_analysis import single_reaction_deletion  # noqa: E402
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+from benchmark_replicate_yeast_gem import (  # noqa: E402
+    _name,
+    _norm,
+    build_draft,
+    curated_reaction_compartments,
+    load_yeast_scores,
+)
+
+from raven_toolbox.localization import apply_assignment, assign_compartments  # noqa: E402
+
+
+def _run(draft, biomass_id, scores, relocate, *, min_growth, time_limit, **kw):
+    """One assign_compartments -> apply_assignment pass; returns (proposal, applied, growth, wall)."""
+    t0 = time.monotonic()
+    prop = assign_compartments(draft, scores, relocate, default_compartment="c", base_metabolite=_name,
+                               biomass_reaction=biomass_id, min_growth=min_growth,
+                               time_limit=time_limit, **kw)
+    wall = round(time.monotonic() - t0, 1)
+    applied = apply_assignment(draft, prop, default_compartment="c", base_metabolite=_name,
+                               universal=kw.get("universal"))
+    applied.objective = biomass_id
+    growth = round(applied.slim_optimize(error_value=0.0) or 0.0, 4)
+    return prop, applied, growth, wall
+
+
+def _agreement(prop, curated_rxn):
+    rc = {rid: _norm(cs[0]) for rid, cs in prop.placements.items() if cs}
+    common = set(rc) & set(curated_rxn)
+    return (round(sum(rc[r] == curated_rxn[r] for r in common) / len(common), 4) if common else None,
+            len(common))
+
+
+def part_prune(draft, biomass_id, scores, relocate, curated_rxn, *, min_growth, time_limit):
+    """prune_transports True vs False: transport count, subset relation, accuracy, runtime."""
+    on, _, g_on, w_on = _run(draft, biomass_id, scores, relocate, min_growth=min_growth,
+                             time_limit=time_limit, prune_transports=True)
+    off, _, g_off, w_off = _run(draft, biomass_id, scores, relocate, min_growth=min_growth,
+                                time_limit=time_limit, prune_transports=False)
+    t_on, t_off = set(on.added_transports), set(off.added_transports)
+    acc_on, n_on = _agreement(on, curated_rxn)
+    acc_off, n_off = _agreement(off, curated_rxn)
+    removed = len(t_off) - len(t_on)
+    return {
+        "transports_pruned": len(t_on), "transports_unpruned": len(t_off),
+        "removed_by_pruning": removed,
+        "removed_fraction": round(removed / len(t_off), 4) if t_off else None,
+        "pruned_is_strict_subset": t_on <= t_off,          # pruning only ever removes
+        "only_in_pruned": len(t_on - t_off),               # expect 0
+        "agreement_pruned": acc_on, "agreement_unpruned": acc_off, "agreement_n": n_on,
+        "accuracy_unchanged": acc_on == acc_off,
+        "growth_pruned": g_on, "growth_unpruned": g_off,
+        "wall_pruned_s": w_on, "wall_unpruned_s": w_off,
+    }
+
+
+def part_gapfill_natural(draft, biomass_id, scores, relocate, universal, *, min_growth, time_limit):
+    """Does gap-fill fire when the flattened draft is simply reassigned (no induced gaps)?"""
+    no_u, _, g_no, _ = _run(draft, biomass_id, scores, relocate, min_growth=min_growth,
+                            time_limit=time_limit)
+    with_u, _, g_u, _ = _run(draft, biomass_id, scores, relocate, min_growth=min_growth,
+                             time_limit=time_limit, universal=universal)
+    return {
+        "added_reactions_without_universal": len(no_u.added_reactions),
+        "added_reactions_with_universal": len(with_u.added_reactions),
+        "certified_without_universal": no_u.certified,
+        "certified_with_universal": with_u.certified,
+        "growth_without_universal": g_no, "growth_with_universal": g_u,
+        "gratuitous_gapfill": len(with_u.added_reactions) > 0,   # expect False
+    }
+
+
+def part_gapfill_knockout(applied, *, floor_fraction, sample_size, seed):
+    """Ground-truthed recovery: remove each essential reaction, gap-fill from a universal that has it.
+
+    The universal is a copy of the certified model, so every removed reaction is a candidate and the
+    ground truth is exact. Uses the same cobra.flux_analysis.gapfill call assign_compartments._gapfill
+    wraps. Reports recovery rate, exact-match rate, and cobra's numerical-failure rate."""
+    import random
+
+    g0 = applied.slim_optimize()
+    floor = floor_fraction * g0
+    dl = single_reaction_deletion(applied, processes=1).reset_index(drop=True)
+    essential = []
+    for _, row in dl.iterrows():
+        ids = row["ids"]
+        if len(ids) != 1:
+            continue
+        rid = next(iter(ids))
+        gr = row["growth"]
+        if (gr is None or gr != gr or gr < floor) and not applied.reactions.get_by_id(rid).boundary:
+            essential.append(rid)
+
+    universal = applied.copy()
+    sample = random.Random(seed).sample(essential, min(sample_size, len(essential)))
+    recovered = exact = cobra_failure = added_nothing = 0
+    for rid in sample:
+        with applied:
+            applied.remove_reactions([applied.reactions.get_by_id(rid)])
+            try:
+                sols = cobra_gapfill(applied, universal, lower_bound=max(floor, 1e-4),
+                                     demand_reactions=False, iterations=1)
+                added = [r.id for r in sols[0]] if sols else []
+            except Exception:  # noqa: BLE001 — cobra tolerance/backend quirk; _gapfill swallows it too
+                cobra_failure += 1
+                continue
+        if not added:
+            added_nothing += 1
+            continue
+        recovered += 1
+        exact += rid in added
+    return {
+        "internal_essential_reactions": len(essential),
+        "sampled": len(sample),
+        "recovered_growth": recovered,
+        "recovery_rate": round(recovered / len(sample), 4) if sample else None,
+        "exact_reaction_readded": exact,
+        "exact_rate_of_recovered": round(exact / recovered, 4) if recovered else None,
+        "cobra_gapfill_numerical_failures": cobra_failure,
+        "added_nothing": added_nothing,
+        "floor_fraction": floor_fraction,
+    }
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--yeast-gem", type=Path, default=Path("C:/Work/GitHub/yeast-GEM/model/yeast-GEM.xml"))
+    ap.add_argument("--data-dir", type=Path, default=Path("data/deeploc"))
+    ap.add_argument("--min-growth-fraction", type=float, default=0.5)
+    ap.add_argument("--time-limit", type=float, default=120.0)
+    ap.add_argument("--knockout-sample", type=int, default=50)
+    ap.add_argument("--seed", type=int, default=1234)
+    ap.add_argument("--out", type=Path, default=Path(".research_tmp/assignment_ablations.json"))
+    args = ap.parse_args(argv)
+
+    yeast = cobra.io.read_sbml_model(str(args.yeast_gem))
+    curated_growth = yeast.slim_optimize()
+    curated_rxn = curated_reaction_compartments(yeast)
+    draft, biomass_id = build_draft(yeast)
+    scores = load_yeast_scores(args.data_dir)
+    relocate = [r.id for r in draft.reactions if not r.boundary and r.id != biomass_id]
+    min_growth = args.min_growth_fraction * curated_growth
+    print(f"curated growth {curated_growth:.4f}; draft {len(draft.reactions)} reactions; "
+          f"relocate {len(relocate)}; min_growth {min_growth:.4f}", flush=True)
+
+    out: dict = {"min_growth": round(min_growth, 4), "curated_growth": round(curated_growth, 4)}
+
+    print("\n[1/3] transport-pruning ablation ...", flush=True)
+    out["pruning"] = part_prune(draft, biomass_id, scores, relocate, curated_rxn,
+                                min_growth=min_growth, time_limit=args.time_limit)
+    p = out["pruning"]
+    print(f"  pruned {p['transports_pruned']} vs unpruned {p['transports_unpruned']} transports "
+          f"({p['removed_by_pruning']} removed, {p['removed_fraction']:.1%}); "
+          f"strict subset={p['pruned_is_strict_subset']}")
+    print(f"  reaction agreement pruned {p['agreement_pruned']} vs unpruned {p['agreement_unpruned']} "
+          f"(unchanged={p['accuracy_unchanged']}); runtime {p['wall_pruned_s']}s vs {p['wall_unpruned_s']}s")
+
+    print("\n[2/3] gap-fill on the natural draft (no induced gaps) ...", flush=True)
+    out["gapfill_natural"] = part_gapfill_natural(draft, biomass_id, scores, relocate, universal=draft,
+                                                  min_growth=min_growth, time_limit=args.time_limit)
+    g = out["gapfill_natural"]
+    print(f"  added reactions: without universal {g['added_reactions_without_universal']}, "
+          f"with universal {g['added_reactions_with_universal']} "
+          f"(gratuitous gap-fill={g['gratuitous_gapfill']})")
+
+    print("\n[3/3] gap-fill knockout-recovery (ground-truthed) ...", flush=True)
+    _base, applied, _g, _w = _run(draft, biomass_id, scores, relocate, min_growth=min_growth,
+                                  time_limit=args.time_limit)
+    out["gapfill_knockout"] = part_gapfill_knockout(applied, floor_fraction=0.05,
+                                                    sample_size=args.knockout_sample, seed=args.seed)
+    k = out["gapfill_knockout"]
+    print(f"  {k['sampled']} essential knockouts: recovered {k['recovered_growth']} "
+          f"({k['recovery_rate']:.1%}); of those, {k['exact_reaction_readded']} re-added the exact "
+          f"reaction ({k['exact_rate_of_recovered']:.1%})")
+    print(f"  cobra.gapfill numerical failures: {k['cobra_gapfill_numerical_failures']}/{k['sampled']}")
+
+    args.out.parent.mkdir(parents=True, exist_ok=True)
+    args.out.write_text(json.dumps(out, indent=2, default=str))
+    print(f"\nwritten -> {args.out}")
+
+
+if __name__ == "__main__":
+    main()