diff --git a/Makefile b/Makefile index 81f315bb..a8ae34ff 100644 --- a/Makefile +++ b/Makefile @@ -448,6 +448,10 @@ bench-simulation-ablation-within-amp: PYTHONPATH=src $(PYTHON) scripts/benchmark_simulation_ablation.py \ --mode within-amp --out outputs/simulation_ablation_within_amp.json +bench-simulation-baselines: + PYTHONPATH=src $(PYTHON) scripts/benchmark_simulation_baselines.py \ + --out outputs/simulation_baselines.json + simulation-gate: PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli bench simulation-gate \ --amp-vs-decoy-json outputs/simulation_ablation.json \ diff --git a/docs/50_LOOP_PLAN.md b/docs/50_LOOP_PLAN.md index db4be0be..aac414d5 100644 --- a/docs/50_LOOP_PLAN.md +++ b/docs/50_LOOP_PLAN.md @@ -104,7 +104,7 @@ Build the multi-resolution virtual assay layer: structure proxies, membrane inte | 33 ✅ | No per-module ablation: do simulation modules improve discrimination? (Original 'selectivity ratio' superseded by Loop 31 membrane proxy — replaced with more useful ablation benchmark) | `scripts/benchmark_simulation_ablation.py`: tests MembraneProxy + StructureProxy on 500-AMP AMP-vs-decoy benchmark. Reports per-module AUROC, combined composite AUROC, and verdict. v0.5.54 | 17 tests, 1922 total. Honest finding: composite degrades AMP-vs-decoy (delta=-0.1153), but bacterial_binding alone achieves 0.7512 AUROC — genuine non-charge signal from Wimley-White + hydrophobic moment | | 34 ✅ | No within-AMP simulation ablation: do simulation modules improve hemolysis detection? (Original 'membrane reference CSV' deferred — higher leverage to test modules on their actual task first) | Extended `scripts/benchmark_simulation_ablation.py` with `--mode within-amp`. Tests MembraneProxy + StructureProxy on 45 hemolytic vs 125 selective AMPs. Reports per-score detection AUROC. v0.5.55 | 23 tests, 1928 total. Honest finding: best simulation helix_weight AUROC 0.6458; rich_selectivity still best at 0.7453. Simulation modules do NOT improve over existing scorers for hemolysis detection | | 35 ✅ | No executable brake prevented future agents from enabling weighted simulation despite negative ablation results | `simulation/gate.py` + `openamp-foundry bench simulation-gate`: fail-closed permission gate consuming both current ablation artifacts. `make simulation-gate` writes a machine-readable verdict. | 6 tests added; current suite 1927 passing, 7 skipped. Current verdict blocks weighted mode and downgrades to info-only | -| 36 | No per-module cheap-baseline benchmark for within-domain simulation value | `benchmark/simulation_baselines.py`: compare each simulation signal against its cheap heuristic baseline on the same task. Report deltas + CI. | Any module with delta ≤ 0 remains permanently experimental | +| 36 ✅ | No per-module cheap-baseline benchmark for within-domain simulation value | `scripts/benchmark_simulation_baselines.py`: compares each simulation signal (bacterial_binding, selectivity_ratio, helix_weight, non_helical) against its cheapest heuristic baseline (mean_eisenberg, selectivity_proxy, helix_propensity, proline_fraction). v0.5.58 | 13 tests, 1953 total. Honest finding: 0/4 signals beat baselines. All simulation modules remain permanently experimental. Weighted mode remains blocked | | 37 ✅ | No `--simulation-mode` CLI flag. Simulation modules exist but users cannot access them | Added `--simulation-mode` to `rank` CLI: `off` (default) and `info` (runs MembraneProxy + StructureProxy, adds `sim_*` keys to scores, includes in JSONL/report). Weighted blocked by gate (Loop 35). v0.5.57 | 6 tests, 1940 total. Demo pipeline with `--simulation-mode info` produces sim scores for all 10 candidates. Report includes sim columns | | 38 | No API contract for third-party simulation modules (e.g. Martini MD, AlphaFold) | `simulation/interfaces.py` extended with `ExternalSimulationAdapter` protocol. Documented in ARCHITECTURE.md extension points | Third-party can implement the protocol without reading pipeline internals | | 39 | Virtual-assay benchmark report comparing simulation-augmented vs simulation-free ranking on the strict triage set | `docs/SIMULATION_BENCHMARK.md`: head-to-head comparison of sel_vs_dec, hem_vs_dec, sel_vs_hem AUROCs with and without simulation. If simulation doesn't improve strict triage, say so explicitly | Honest: if delta ≤ 0, report "simulation did not improve triage" | @@ -241,7 +241,7 @@ Phase 4: Not started (Loops 40–49) | 28 ✅ | Policy version bump workflow for when real data arrives | `scripts/bump_recalibration_policy.py`: standalone script with `--dry-run`, decision-log guard, auto-increment + write. CI guard in `ci.yml` validates policy version changes against base branch. v0.5.49. 9 tests. | CI rejects policy PRs without valid decision log; `make bump-policy-version` and `make bump-policy-version-dry-run` available | | 29 ✅ | No public negative-result archive format. If Wave 1 yields all negatives, where do they go? | `docs/NEGATIVE_RESULT_ARCHIVE.md`: full template with entry schema, procedures, automation notes, and limitations. Covers pre-selection rejects, selected-untested, lab inactives, lab toxic, control failures. v0.5.50. | Template complete enough for a lab partner to fill; schema defines 18 fields with required/conditional markers | -**Next loop:** Loop 36 — Phase 3 cheap-baseline benchmark for each simulation signal. +**Next loop:** Loop 38 — Phase 3 (external adapter protocol). **Phase 2 exit criteria (all 5 met ✅):** - ✅ `make calibration-loop` runs from clean checkout, produces batch-2 manifest diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index 82c7a9ad..b350a2ee 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -60,6 +60,7 @@ documented result in `docs/METRICS_CURRENT.md` and | `bench-expert-ablation-500` | Ablate expert composite on expanded n=1000 | Re-run on diverse set; 2 components reclassified from n=191 | `make bench-expert-ablation-500` | | `bench-simulation-ablation` | Ablate MembraneProxy + StructureProxy on 500-AMP benchmark | Composite degrades AMP-vs-decoy (expected — modules are for within-AMP tasks); bacterial_binding alone AUROC 0.7512 (non-charge signal) | `make bench-simulation-ablation` | | `bench-simulation-ablation-within-amp` | Ablate simulation modules on within-AMP hemolysis detection | Best simulation helix_weight AUROC 0.6458; rich_selectivity best at 0.7453. No improvement | `make bench-simulation-ablation-within-amp` | +| `bench-simulation-baselines` | Compare each simulation signal vs its cheapest heuristic baseline | 0/4 signals beat baselines. All modules permanently experimental | `make bench-simulation-baselines` | ### Regression / CI diff --git a/docs/METRICS_CURRENT.md b/docs/METRICS_CURRENT.md index 9a90dead..204792a3 100644 --- a/docs/METRICS_CURRENT.md +++ b/docs/METRICS_CURRENT.md @@ -5,7 +5,8 @@ Machine-readable snapshot: `outputs/metrics_snapshot.json` regenerated with `mak > **Purpose:** One authoritative table of current pipeline metrics. If any doc disagrees > with this file, this file wins. Updated whenever benchmark/benchmark config changes. > -> **Last updated:** 2026-07-06 (simulation mode CLI — v0.5.57) +> **Last updated:** 2026-07-06 (simulation cheap-baseline benchmark — v0.5.58) +> **New in v0.5.58:** Per-signal cheap-baseline comparison: 0/4 simulation signals beat their cheapest heuristic. All simulation modules remain permanently experimental. `make bench-simulation-baselines`. 13 tests. 1953 total. > **New in v0.5.57:** `rank --simulation-mode info` runs MembraneProxy + StructureProxy on every candidate, adds `sim_*` scores to JSONL output and Markdown report. 6 new CLI tests. 1940 total. > **New in v0.5.56:** `openamp-foundry bench simulation-gate` converts both > simulation ablation artifacts into a fail-closed permission decision for @@ -42,7 +43,7 @@ Machine-readable snapshot: `outputs/metrics_snapshot.json` regenerated with `mak > **New in v0.5.31:** Added dipeptide-order features for sequence-order awareness. `dipeptide_order_score` achieves AUROC 0.7861 on AMP-vs-scrambled discrimination — the strongest order-dependent feature in the pipeline. Only 7/31 features survive scrambling (amphipathicity/helix-wheel + dipeptide). All composition features are purely position-independent (exactly 0.5000 AUROC on scrambled test). > **New in v0.5.30:** Easy baseline benchmark added — charge density alone (AUROC 0.8166) outperforms the full pipeline ensemble (0.7792) on AMP-vs-Swiss-Prot-decoy discrimination. Honest finding documented: expected because pipeline optimizes for safety, not raw discrimination. > **New in v0.5.29:** Expanded benchmark to 500 AMPs + 500 composition-matched decoys (n=1000). AUROC 0.7792 (CI₉₅: 0.7505–0.8065) confirms signal generalizes. Cluster-aware CI: 0.746–0.8102. Representative AUROC: 0.778. Standard benchmark (n=191) retained for backward comparison. -> **Pipeline version:** v0.5.57 +> **Pipeline version:** v0.5.58 > **Branch:** main --- @@ -1112,7 +1113,7 @@ Decoys score low on activity. Selective AMPs score moderately on both. | Metric | Value | |--------|-------| -| Total tests | 1940 | +| Total tests | 1953 | | Coverage (branch) | 99% (6 CLI guard lines only) | | Source modules at 100% | All pipeline, QC, scoring modules | diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 4974543e..b80d06a9 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1,5 +1,30 @@ # Roadmap +## v0.5.58 — Simulation Cheap-Baseline Benchmark (Loop 36) ✓ (2026-07-06) + +Per-signal comparison: does each simulation output beat its cheapest +meaningful heuristic baseline on the within-AMP hemolysis task? + +Changes: +- ``scripts/benchmark_simulation_baselines.py`` — Compares bacterial_binding + vs mean_eisenberg, selectivity_ratio vs selectivity_proxy, helix_weight + vs helix_propensity, non_helical vs proline_fraction. ``make + bench-simulation-baselines`` target. +- ``tests/test_benchmark_simulation_baselines.py`` — 13 tests. +- ``docs/50_LOOP_PLAN.md`` — Loop 36 ✅. Next: Loop 38. +- ``docs/METRICS_CURRENT.md`` — Pipeline version v0.5.58. Test count: 1953. +- ``docs/BENCHMARKING.md`` — Added bench-simulation-baselines (22 targets). + +**Honest finding: 0/4 simulation signals beat their cheapest baseline.** +- bacterial_binding (0.4876) vs mean_eisenberg (0.5469): delta=−0.0593 +- selectivity_ratio (0.3615) vs selectivity_proxy (0.3905): delta=−0.0290 +- helix_weight (0.6458) vs helix_propensity (0.6489): delta=−0.0031 +- non_helical (0.4124) vs proline_fraction (0.4929): delta=−0.0805 + +**Consequence:** All current simulation modules remain permanently +experimental. Weighted mode stays blocked. Future simulation work +must use fundamentally richer modeling (not just 1D propensities). + ## v0.5.57 — Simulation Mode CLI Flag (Loop 37) ✓ (2026-07-06) Phase 3 simulation modules are now user-accessible via the `rank` CLI. diff --git a/scripts/benchmark_simulation_baselines.py b/scripts/benchmark_simulation_baselines.py new file mode 100644 index 00000000..6c9b41f9 --- /dev/null +++ b/scripts/benchmark_simulation_baselines.py @@ -0,0 +1,208 @@ +"""Simulation cheap-baseline benchmark: does each simulation signal beat its cheapest heuristic? + +Compares each simulation module output against its cheapest meaningful +baseline on the within-AMP hemolysis detection task. + +Simulation signal → cheap baseline: +- bacterial_binding → mean Eisenberg hydrophobicity +- selectivity_ratio → selectivity_proxy (charge + GRAVY) +- helix_weight → helix_propensity_score (mean Chou-Fasman Pα) +- non_helical → proline_fraction + +Any module with delta ≤ 0 against its baseline is flagged as NO_IMPROVEMENT +and must remain permanently experimental. +""" + +from __future__ import annotations + +import csv +import json +import sys +from pathlib import Path + +from openamp_foundry.features.physchem import ( + _HYDROPHOBICITY, + compute_features, + fraction, + helix_propensity_score, + selectivity_proxy, +) +from openamp_foundry.simulation.membrane import MembraneProxy +from openamp_foundry.simulation.structure import StructureProxy + +HYDROPHOBIC = set("AILMFWVY") + + +def _auroc(scores_pos: list[float], scores_neg: list[float]) -> float: + n_p = len(scores_pos) + n_n = len(scores_neg) + if n_p == 0 or n_n == 0: + return 0.5 + better = 0 + ties = 0 + for sp in scores_pos: + for sn in scores_neg: + if sp > sn: + better += 1 + elif sp == sn: + ties += 1 + total = n_p * n_n + return (better + 0.5 * ties) / total if total else 0.5 + + +def _mean_eisenberg(sequence: str) -> float: + if not sequence: + return 0.0 + total = 0.0 + count = 0 + for aa in sequence.upper(): + if aa in _HYDROPHOBICITY: + total += _HYDROPHOBICITY[aa] + count += 1 + return total / count if count else 0.0 + + +def _detection_auroc(hemo_vals: list[float], sel_vals: list[float], invert: bool = False) -> float: + raw = _auroc(hemo_vals, sel_vals) + return 1.0 - raw if invert else raw + + +def run_simulation_baselines( + hemolysis_csv: str = "examples/validation/hemolysis_reference.csv", +) -> dict: + """Compare each simulation signal against its cheap heuristic baseline. + + Returns per-signal deltas and verdict. + """ + csv_p = Path(hemolysis_csv) + if not csv_p.exists(): + return {"error": f"Hemolysis CSV not found: {hemolysis_csv}"} + + hemolytic: list[tuple[str, str]] = [] + selective: list[tuple[str, str]] = [] + with csv_p.open("r", newline="") as f: + for row in csv.DictReader(f): + cls = row.get("hemolysis_class", "").strip().upper() + seq = row.get("sequence", "").strip() + cid = row.get("id", "unknown") + if not seq: + continue + if cls == "HEMOLYTIC": + hemolytic.append((cid, seq)) + elif cls == "SELECTIVE": + selective.append((cid, seq)) + + if len(hemolytic) < 3 or len(selective) < 3: + return {"error": f"Not enough sequences: {len(hemolytic)} hemolytic, {len(selective)} selective"} + + membrane = MembraneProxy() + structure = StructureProxy() + + hemo_scores: dict[str, list[float]] = { + "bacterial_binding": [], + "selectivity_ratio": [], + "helix_weight": [], + "non_helical": [], + "mean_eisenberg": [], + "selectivity_proxy": [], + "helix_propensity": [], + "proline_fraction": [], + } + sel_scores: dict[str, list[float]] = {k: [] for k in hemo_scores} + + def _score_and_record(seqs: list[tuple[str, str]], target: dict[str, list[float]]) -> None: + for cid, seq in seqs: + feats = compute_features(seq) + mem_result = membrane.simulate(seq) + struct_result = structure.simulate(seq) + target["bacterial_binding"].append(mem_result.scores.get("bacterial_binding", 0.0)) + target["selectivity_ratio"].append(mem_result.scores.get("selectivity_ratio", 1.0)) + target["helix_weight"].append(struct_result.scores.get("helix_weight", 0.33)) + target["non_helical"].append(struct_result.scores.get("non_helical", 0.0)) + target["mean_eisenberg"].append(_mean_eisenberg(seq)) + net_charge = feats.get("net_charge_ph74", 0.0) + gravy = feats.get("gravy", 0.0) + target["selectivity_proxy"].append(selectivity_proxy(net_charge, gravy)) + target["helix_propensity"].append(helix_propensity_score(seq)) + target["proline_fraction"].append(fraction(seq, {"P"})) + + _score_and_record(hemolytic, hemo_scores) + _score_and_record(selective, sel_scores) + + # Compute raw AUROCs (no direction inversion). + # For each score we report: auroc_raw = P(hemolytic_score > selective_score). + # If auroc_raw > 0.5: hemolytic tends to score HIGHER. + # If auroc_raw < 0.5: hemolytic tends to score LOWER. + # The comparison is symmetric: both signal and baseline use the same metric. + hemo_raw = lambda k: round(_auroc(hemo_scores[k], sel_scores[k]), 4) + + pairs = [ + ("bacterial_binding", "mean_eisenberg"), + ("selectivity_ratio", "selectivity_proxy"), + ("helix_weight", "helix_propensity"), + ("non_helical", "proline_fraction"), + ] + + per_signal = {} + for signal, baseline in pairs: + sig_auroc = hemo_raw(signal) + base_auroc = hemo_raw(baseline) + delta = round(sig_auroc - base_auroc, 4) + verdict = "NO_IMPROVEMENT" if delta <= 0 else "IMPROVEMENT" + per_signal[signal] = { + "signal_auroc": sig_auroc, + "baseline_auroc": base_auroc, + "baseline_name": baseline, + "delta": delta, + "verdict": verdict, + "note": ( + f"auroc > 0.5 = hemolytic AMPs score higher. " + f"Signal={sig_auroc:.4f}, baseline={base_auroc:.4f}, delta={delta:+.4f}" + ), + } + + n_improved = sum(1 for v in per_signal.values() if v["verdict"] == "IMPROVEMENT") + n_total = len(per_signal) + + return { + "benchmark": "simulation_baselines", + "n_hemolytic": len(hemolytic), + "n_selective": len(selective), + "results": per_signal, + "summary": { + "n_improved": n_improved, + "n_total": n_total, + "all_blocked": n_improved == 0, + }, + } + + +def main() -> int: + import argparse + parser = argparse.ArgumentParser(description="Simulation cheap-baseline benchmark") + parser.add_argument("--hemolysis-csv", default="examples/validation/hemolysis_reference.csv") + parser.add_argument("--out", default=None) + args = parser.parse_args() + + result = run_simulation_baselines(hemolysis_csv=args.hemolysis_csv) + if "error" in result: + print(f"Error: {result['error']}", file=sys.stderr) + return 2 + + print(json.dumps(result, indent=2)) + if args.out: + Path(args.out).write_text(json.dumps(result, indent=2) + "\n") + + summary = result["summary"] + print(f"\nResults: {summary['n_improved']}/{summary['n_total']} signals improve over baseline") + for name, data in result["results"].items(): + print(f" {name}: signal={data['signal_auroc']:.4f} baseline={data['baseline_auroc']:.4f} " + f"delta={data['delta']:+.4f} → {data['verdict']}") + + if summary["all_blocked"]: + return 3 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/test_benchmark_simulation_baselines.py b/tests/test_benchmark_simulation_baselines.py new file mode 100644 index 00000000..20e083dc --- /dev/null +++ b/tests/test_benchmark_simulation_baselines.py @@ -0,0 +1,102 @@ +"""Tests for the simulation cheap-baseline benchmark.""" + +import json +import subprocess +import sys + +from openamp_foundry.features.physchem import _HYDROPHOBICITY +from scripts.benchmark_simulation_baselines import ( + _auroc, + _mean_eisenberg, + run_simulation_baselines, +) + + +def test_mean_eisenberg_known(): + score = _mean_eisenberg("AAAA") + expected = _HYDROPHOBICITY["A"] + assert score == expected + + +def test_mean_eisenberg_empty(): + assert _mean_eisenberg("") == 0.0 + + +def test_auroc_all_higher(): + assert _auroc([1.0, 1.0], [0.0, 0.0]) == 1.0 + + +def test_auroc_all_lower(): + assert _auroc([0.0, 0.0], [1.0, 1.0]) == 0.0 + + +def test_auroc_empty(): + assert _auroc([], []) == 0.5 + + +def test_benchmark_runs(): + result = run_simulation_baselines() + assert "error" not in result + assert result["n_hemolytic"] > 3 + assert result["n_selective"] > 3 + + +def test_benchmark_returns_expected_keys(): + result = run_simulation_baselines() + for signal in ["bacterial_binding", "selectivity_ratio", "helix_weight", "non_helical"]: + assert signal in result["results"] + r = result["results"][signal] + for key in ["signal_auroc", "baseline_auroc", "delta", "verdict"]: + assert key in r, f"{signal} missing {key}" + + +def test_benchmark_all_no_improvement(): + """Current honest finding: NO simulation signal beats its cheap baseline.""" + result = run_simulation_baselines() + for signal, data in result["results"].items(): + assert data["verdict"] == "NO_IMPROVEMENT", f"{signal} unexpectedly improved" + + +def test_benchmark_summary_all_blocked(): + result = run_simulation_baselines() + assert result["summary"]["all_blocked"] is True + + +def test_benchmark_heli_weight_near_heli_propensity(): + """helix_weight and helix_propensity should have similar AUROC + (both measure helicity, just normalized differently).""" + result = run_simulation_baselines() + hw = result["results"]["helix_weight"] + hp_auroc = hw["baseline_auroc"] + hw_auroc = hw["signal_auroc"] + assert abs(hw_auroc - hp_auroc) < 0.1 + + +def test_cli_exit_code(): + result = subprocess.run( + [sys.executable, "scripts/benchmark_simulation_baselines.py"], + capture_output=True, text=True, + env={"PYTHONPATH": "src"}, + ) + assert result.returncode == 3 + + +def test_cli_output_contains_verdict(): + result = subprocess.run( + [sys.executable, "scripts/benchmark_simulation_baselines.py"], + capture_output=True, text=True, + env={"PYTHONPATH": "src"}, + ) + assert "NO_IMPROVEMENT" in result.stdout + + +def test_cli_writes_output_json(tmp_path): + out = tmp_path / "baselines.json" + result = subprocess.run( + [sys.executable, "scripts/benchmark_simulation_baselines.py", "--out", str(out)], + capture_output=True, text=True, + env={"PYTHONPATH": "src"}, + ) + assert out.exists() + data = json.loads(out.read_text()) + assert data["summary"]["all_blocked"] is True