{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"3a144037","cell_type":"markdown","source":"# Phase 8 — Robust Benchmarking (real + simAIRR + LIgO implanted motifs)\n\nEvaluates each branch on (1) real AIRR challenge data, (2) simAIRR-style realistic sharing benchmark, (3) LIgO-style implanted low-witness-rate benchmark. Reveals which branch captures public signal vs. recovers rare implanted motifs vs. exploits simulation artifacts.\n\n---\n\n## Kaggle inputs to add before running this notebook\n\nAdd the following as Kaggle inputs (via the \"Add Input\" button on the right\npanel of the notebook editor):\n- AIRR-ML competition dataset\n- Phase 1 notebook output\n- Phase 6 notebook output (ranked attribution tables)\n- Phase 7 notebook output (fusion_oof.csv)\n\n## Workflow\n\n1. Run the **SETUP** cell — it creates `/kaggle/working/project/` and auto-merges\n   any previous-phase notebook outputs found under `/kaggle/input/`.\n2. Run each subsequent cell in order. The **Run** cell executes the phase's\n   training script; the **Inspect** cell prints a quick summary of the outputs.\n3. When the run completes, click **Save Version → Save & Run All (Commit)** so\n   the next phase can pick this phase's outputs up via \"Add Input\".\n","metadata":{}},{"id":"cdb12981","cell_type":"code","source":"# ============================================================\n# SETUP — Initialize project + merge previous-phase inputs\n# ============================================================\n# This cell:\n#   1. Creates /kaggle/working/project/ fresh (idempotent re-runs).\n#   2. Auto-detects previous-phase notebook outputs under /kaggle/input/.\n#   3. Merges their project/ contents (src/, artifacts/, configs/) into\n#      /kaggle/working/project/ so this phase can build on them.\n#\n# Kaggle \"Add Input\" workflow:\n#   - Phase 1: add the AIRR-ML competition dataset only.\n#   - Phase N (N>=2): add the AIRR-ML competition dataset AND the previous\n#     phase notebook output(s). For Phase 9 add ALL of Phases 1..8.\n# ============================================================\n\nimport os\nimport shutil\nfrom pathlib import Path\n\nPROJECT_ROOT = Path(\"/kaggle/working/project\")\n\n# Reset working project dir (idempotent re-runs)\nif PROJECT_ROOT.exists():\n    shutil.rmtree(PROJECT_ROOT)\nPROJECT_ROOT.mkdir(parents=True, exist_ok=True)\n(PROJECT_ROOT / \"src\").mkdir(parents=True, exist_ok=True)\n(PROJECT_ROOT / \"src\" / \"__init__.py\").write_text(\"\", encoding=\"utf-8\")\n(PROJECT_ROOT / \"configs\").mkdir(parents=True, exist_ok=True)\n(PROJECT_ROOT / \"artifacts\").mkdir(parents=True, exist_ok=True)\n\n# Auto-detect and merge previous-phase inputs.\n# Each previous-phase notebook output should contain a top-level `project/`\n# directory (created by Phase 1 and propagated through every later phase).\nprev_inputs_found = []\ninput_root = Path(\"/kaggle/input\")\nif input_root.exists():\n    for entry in sorted(input_root.iterdir()):\n        if not entry.is_dir():\n            continue\n        prev_project = entry / \"project\"\n        if not prev_project.is_dir():\n            continue\n        prev_inputs_found.append(entry.name)\n        for sub in [\"src\", \"artifacts\", \"configs\"]:\n            src_dir = prev_project / sub\n            if not src_dir.exists():\n                continue\n            for path in src_dir.rglob(\"*\"):\n                if path.is_file():\n                    rel = path.relative_to(src_dir)\n                    target = PROJECT_ROOT / sub / rel\n                    target.parent.mkdir(parents=True, exist_ok=True)\n                    shutil.copy2(path, target)\n\nprint(\"PROJECT_ROOT :\", PROJECT_ROOT)\nif prev_inputs_found:\n    print(\"Previous-phase inputs merged:\")\n    for name in prev_inputs_found:\n        print(\"  -\", name)\nelse:\n    print(\"Previous-phase inputs: (none — running fresh)\")\nprint(\"\\nExisting artifacts:\")\nart_dir = PROJECT_ROOT / \"artifacts\"\nif art_dir.exists():\n    for p in sorted(art_dir.glob(\"*\")):\n        print(\"  -\", p.name)\nelse:\n    print(\"  (none)\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:29:59.570063Z","iopub.execute_input":"2026-07-24T11:29:59.570565Z","iopub.status.idle":"2026-07-24T11:29:59.583020Z","shell.execute_reply.started":"2026-07-24T11:29:59.570536Z","shell.execute_reply":"2026-07-24T11:29:59.582040Z"}},"outputs":[],"execution_count":null},{"id":"9c797094","cell_type":"code","source":"# Cell 1: Phase 8 benchmark config file write\n\nfrom pathlib import Path\n\nPROJECT_ROOT = Path(\"/kaggle/working/project\")\n\nbenchmark_yaml = \"\"\"\nenabled: true\n\npaths:\n  output_root: /kaggle/working/project/artifacts/phase8_benchmarks\n  canonical_metadata: /kaggle/working/project/artifacts/phase1/canonical_metadata.csv\n  fusion_oof: /kaggle/working/project/artifacts/phase7_fusion/fusion_oof.csv\n  phase6_dir: /kaggle/working/project/artifacts/phase6_attribution\n\nreal:\n  topk_list: [5, 10, 20]\n\nsimairr:\n  n_pos: 120\n  n_neg: 120\n  n_background_per_rep: 220\n  n_public_background: 180\n  n_disease_public: 60\n  public_common_prob: 0.18\n  disease_prob_pos: 0.22\n  disease_prob_neg: 0.03\n  random_state: 42\n\nligo:\n  n_pos: 120\n  n_neg: 120\n  n_background_per_rep: 240\n  implanted_motifs: [\"CASSLG\", \"ASSIR\", \"GGQET\", \"QYFDT\"]\n  implants_per_positive: 3\n  implant_prob_positive: 1.0\n  implant_prob_negative: 0.02\n  recovery_k_list: [1, 5, 10]\n  random_state: 42\n\ndictionary:\n  max_repertoires_per_class: 150\n  min_freq: 0.10\n  enrichment: 3.0\n  top_exact: 3000\n  top_clusters: 5000\n\"\"\"\n\n(PROJECT_ROOT / \"configs\" / \"benchmark.yaml\").write_text(benchmark_yaml.strip() + \"\\n\", encoding=\"utf-8\")\nprint(\"Written:\", PROJECT_ROOT / \"configs\" / \"benchmark.yaml\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:29:59.584735Z","iopub.execute_input":"2026-07-24T11:29:59.585313Z","iopub.status.idle":"2026-07-24T11:29:59.598870Z","shell.execute_reply.started":"2026-07-24T11:29:59.585289Z","shell.execute_reply":"2026-07-24T11:29:59.598290Z"}},"outputs":[],"execution_count":null},{"id":"05319e94","cell_type":"code","source":"# Cell 2: Write src/simulation.py, src/evaluation.py, and train_benchmarks.py\n\nfrom pathlib import Path\nfrom textwrap import dedent\n\nPROJECT_ROOT = Path(\"/kaggle/working/project\")\n\nsimulation_py = dedent(\"\"\"\nfrom __future__ import annotations\n\nfrom collections import Counter\nfrom pathlib import Path\nfrom typing import Dict, List, Tuple\n\nimport numpy as np\nimport pandas as pd\n\nfrom .branch_a_features import read_repertoire\nfrom .branch_b_cluster import approx_cluster_keys\nfrom .utils import stable_hash\n\n\nAA_LIST = list(\"ARNDCQEGHILKMFPSTWYV\")\n\n\ndef estimate_empirical_background(\n    canonical_train_meta: pd.DataFrame,\n    train_root: Path,\n    max_repertoires: int = 16,\n    max_sequences_per_file: int = 4000,\n    random_state: int = 42,\n):\n    rng = np.random.RandomState(random_state)\n    if len(canonical_train_meta) == 0:\n        aa_probs = np.ones(len(AA_LIST), dtype=float) / len(AA_LIST)\n        lengths = np.arange(10, 21)\n        return aa_probs, lengths\n\n    sample_meta = canonical_train_meta.sample(\n        n=min(max_repertoires, len(canonical_train_meta)),\n        random_state=random_state,\n        replace=False,\n    )\n\n    aa_counter = Counter()\n    lengths = []\n\n    for _, row in sample_meta.iterrows():\n        ds_path = train_root / row[\"dataset_name\"]\n        df = read_repertoire(ds_path / row[\"filename\"], max_seqs=max_sequences_per_file, random_state=random_state)\n        if len(df) == 0:\n            continue\n        seqs = df[\"junction_aa\"].fillna(\"\").astype(str).tolist()\n        for s in seqs:\n            if not s:\n                continue\n            lengths.append(len(s))\n            for ch in s:\n                if ch in AA_LIST:\n                    aa_counter[ch] += 1\n\n    if len(lengths) == 0:\n        lengths = list(range(10, 21))\n\n    aa_probs = np.array([aa_counter.get(aa, 1) for aa in AA_LIST], dtype=float)\n    aa_probs = aa_probs / aa_probs.sum()\n    lengths = np.asarray(lengths, dtype=int)\n    return aa_probs, lengths\n\n\ndef random_seq(length: int, aa_probs: np.ndarray, rng: np.random.RandomState) -> str:\n    toks = rng.choice(AA_LIST, size=int(length), replace=True, p=aa_probs)\n    return \"\".join(toks.tolist())\n\n\ndef random_lengths(length_source: np.ndarray, n: int, rng: np.random.RandomState) -> np.ndarray:\n    if len(length_source) == 0:\n        length_source = np.arange(10, 21)\n    idx = rng.choice(np.arange(len(length_source)), size=n, replace=True)\n    return np.asarray(length_source[idx], dtype=int)\n\n\ndef make_public_pool(n_seq: int, aa_probs: np.ndarray, length_source: np.ndarray, rng: np.random.RandomState) -> List[str]:\n    lengths = random_lengths(length_source, n_seq, rng)\n    return [random_seq(int(L), aa_probs, rng) for L in lengths]\n\n\ndef build_counter_from_seq_list(seqs: List[str], weights: List[float]) -> Counter:\n    c = Counter()\n    for s, w in zip(seqs, weights):\n        if s:\n            c[s] += float(w)\n    return c\n\n\ndef simulate_simairr_dataset(\n    n_pos: int,\n    n_neg: int,\n    n_background_per_rep: int,\n    n_public_background: int,\n    n_disease_public: int,\n    public_common_prob: float,\n    disease_prob_pos: float,\n    disease_prob_neg: float,\n    aa_probs: np.ndarray,\n    length_source: np.ndarray,\n    random_state: int = 42,\n):\n    rng = np.random.RandomState(random_state)\n\n    public_bg_pool = make_public_pool(n_public_background, aa_probs, length_source, rng)\n    disease_pool = make_public_pool(n_disease_public, aa_probs, length_source, rng)\n\n    records = []\n    cache = {}\n\n    def make_one(rep_id: str, label: int):\n        seqs = []\n        weights = []\n\n        # background unique-ish sequences\n        bg_lengths = random_lengths(length_source, n_background_per_rep, rng)\n        for L in bg_lengths:\n            s = random_seq(int(L), aa_probs, rng)\n            seqs.append(s)\n            weights.append(float(rng.lognormal(mean=0.4, sigma=0.8)))\n\n        # realistic public sharing in both classes\n        for s in public_bg_pool:\n            if rng.rand() < public_common_prob:\n                seqs.append(s)\n                weights.append(float(rng.lognormal(mean=0.8, sigma=0.5)))\n\n        # disease-enriched public signatures\n        p = disease_prob_pos if label == 1 else disease_prob_neg\n        for s in disease_pool:\n            if rng.rand() < p:\n                seqs.append(s)\n                weights.append(float(rng.lognormal(mean=1.0, sigma=0.6)))\n\n        seq_counter = build_counter_from_seq_list(seqs, weights)\n        cache[rep_id] = {\n            \"seq_counter\": seq_counter,\n            \"total_weight\": float(sum(seq_counter.values())),\n        }\n        records.append({\n            \"repertoire_id\": rep_id,\n            \"dataset_name\": \"simairr_synth\",\n            \"label_positive\": int(label),\n        })\n\n    for i in range(n_pos):\n        make_one(f\"simairr_pos_{i:04d}\", 1)\n    for i in range(n_neg):\n        make_one(f\"simairr_neg_{i:04d}\", 0)\n\n    meta = pd.DataFrame(records)\n    truth = {\n        \"public_background_pool\": public_bg_pool,\n        \"disease_pool\": disease_pool,\n    }\n    return meta, cache, truth\n\n\ndef implant_motif_in_sequence(seq: str, motif: str, rng: np.random.RandomState) -> str:\n    if len(seq) < len(motif):\n        seq = seq + motif\n        return seq[: max(len(motif), len(seq))]\n\n    start_max = max(0, len(seq) - len(motif))\n    pos = int(rng.randint(0, start_max + 1))\n    out = seq[:pos] + motif + seq[pos + len(motif):]\n    return out\n\n\ndef simulate_ligo_dataset(\n    n_pos: int,\n    n_neg: int,\n    n_background_per_rep: int,\n    implanted_motifs: List[str],\n    implants_per_positive: int,\n    implant_prob_positive: float,\n    implant_prob_negative: float,\n    aa_probs: np.ndarray,\n    length_source: np.ndarray,\n    random_state: int = 42,\n):\n    rng = np.random.RandomState(random_state)\n\n    records = []\n    cache = {}\n    truth_rows = []\n\n    def make_one(rep_id: str, label: int):\n        seqs = []\n        weights = []\n        implanted = []\n\n        bg_lengths = random_lengths(length_source, n_background_per_rep, rng)\n        bg_lengths = np.maximum(bg_lengths, max(len(m) for m in implanted_motifs) + 2)\n\n        bg_sequences = [random_seq(int(L), aa_probs, rng) for L in bg_lengths]\n\n        # choose implants\n        local_sequences = list(bg_sequences)\n        n_implants = implants_per_positive if (label == 1 and rng.rand() < implant_prob_positive) else 0\n        if label == 0 and rng.rand() < implant_prob_negative:\n            n_implants = 1\n\n        chosen_positions = rng.choice(np.arange(len(local_sequences)), size=min(n_implants, len(local_sequences)), replace=False).tolist() if n_implants > 0 else []\n\n        for pos_idx in chosen_positions:\n            motif = implanted_motifs[pos_idx % len(implanted_motifs)]\n            s = implant_motif_in_sequence(local_sequences[pos_idx], motif, rng)\n            local_sequences[pos_idx] = s\n            implanted.append(s)\n\n        for s in local_sequences:\n            seqs.append(s)\n            weights.append(float(rng.lognormal(mean=0.45, sigma=0.85)))\n\n        seq_counter = build_counter_from_seq_list(seqs, weights)\n        cache[rep_id] = {\n            \"seq_counter\": seq_counter,\n            \"total_weight\": float(sum(seq_counter.values())),\n        }\n        records.append({\n            \"repertoire_id\": rep_id,\n            \"dataset_name\": \"ligo_synth\",\n            \"label_positive\": int(label),\n        })\n\n        for s in implanted:\n            truth_rows.append({\n                \"repertoire_id\": rep_id,\n                \"sequence\": s,\n                \"is_implanted\": 1,\n                \"label_positive\": int(label),\n            })\n\n    for i in range(n_pos):\n        make_one(f\"ligo_pos_{i:04d}\", 1)\n    for i in range(n_neg):\n        make_one(f\"ligo_neg_{i:04d}\", 0)\n\n    meta = pd.DataFrame(records)\n    truth_df = pd.DataFrame(truth_rows)\n    return meta, cache, truth_df\n\n\ndef sequence_scores_from_bundle(seq_counter: Counter, bundle: dict) -> pd.DataFrame:\n    exact_catalog = bundle.get(\"exact_catalog\", {})\n    cluster_catalog = bundle.get(\"cluster_catalog\", {})\n\n    rows = []\n    for seq, cw in seq_counter.items():\n        public_score = 0.0\n        cluster_score = 0.0\n\n        if seq in exact_catalog:\n            public_score = float(exact_catalog[seq][\"score\"])\n\n        for key in approx_cluster_keys(seq):\n            if key in cluster_catalog:\n                cluster_score = max(cluster_score, float(cluster_catalog[key][\"score\"]))\n\n        total_score = public_score + cluster_score\n        rows.append({\n            \"sequence\": seq,\n            \"weight\": float(cw),\n            \"public_score\": float(public_score),\n            \"cluster_score\": float(cluster_score),\n            \"score\": float(total_score),\n        })\n\n    out = pd.DataFrame(rows)\n    if len(out):\n        out = out.sort_values([\"score\", \"public_score\", \"cluster_score\", \"weight\"], ascending=[False, False, False, False]).reset_index(drop=True)\n        out[\"rank\"] = np.arange(1, len(out) + 1)\n    return out\n\"\"\")\n\nevaluation_py = dedent(\"\"\"\nfrom __future__ import annotations\n\nfrom pathlib import Path\nfrom typing import List\n\nimport numpy as np\nimport pandas as pd\n\n\ndef brier_score(y_true, y_prob):\n    y_true = np.asarray(y_true, dtype=float)\n    y_prob = np.asarray(y_prob, dtype=float)\n    valid = ~np.isnan(y_true) & ~np.isnan(y_prob)\n    y_true = y_true[valid]\n    y_prob = y_prob[valid]\n    if len(y_true) == 0:\n        return np.nan\n    return float(np.mean((y_prob - y_true) ** 2))\n\n\ndef expected_calibration_error(y_true, y_prob, n_bins: int = 10):\n    y_true = np.asarray(y_true, dtype=float)\n    y_prob = np.asarray(y_prob, dtype=float)\n    valid = ~np.isnan(y_true) & ~np.isnan(y_prob)\n    y_true = y_true[valid]\n    y_prob = y_prob[valid]\n\n    if len(y_true) == 0:\n        return np.nan\n\n    bins = np.linspace(0, 1, n_bins + 1)\n    ece = 0.0\n    n = len(y_true)\n\n    for i in range(n_bins):\n        lo, hi = bins[i], bins[i + 1]\n        if i == n_bins - 1:\n            mask = (y_prob >= lo) & (y_prob <= hi)\n        else:\n            mask = (y_prob >= lo) & (y_prob < hi)\n\n        if mask.sum() == 0:\n            continue\n\n        acc = y_true[mask].mean()\n        conf = y_prob[mask].mean()\n        ece += (mask.sum() / n) * abs(acc - conf)\n\n    return float(ece)\n\n\ndef safe_auc(y_true, y_prob):\n    from sklearn.metrics import roc_auc_score\n\n    y_true = np.asarray(y_true)\n    y_prob = np.asarray(y_prob)\n    valid = ~pd.isna(y_true) & ~pd.isna(y_prob)\n    y_true = y_true[valid]\n    y_prob = y_prob[valid]\n\n    if len(y_true) == 0 or len(np.unique(y_true)) < 2:\n        return np.nan\n    return float(roc_auc_score(y_true, y_prob))\n\n\ndef load_ranked_attribution_tables(phase6_dir: Path) -> pd.DataFrame:\n    parts = []\n\n    parquets = sorted(phase6_dir.glob(\"ranked_sequence_attribution_*.parquet\"))\n    for p in parquets:\n        try:\n            parts.append(pd.read_parquet(p))\n        except Exception:\n            pass\n\n    fallbacks = sorted(phase6_dir.glob(\"ranked_sequence_attribution_*_fallback.csv\"))\n    for p in fallbacks:\n        try:\n            parts.append(pd.read_csv(p))\n        except Exception:\n            pass\n\n    if len(parts) == 0:\n        return pd.DataFrame()\n    return pd.concat(parts, ignore_index=True)\n\n\ndef real_benchmark_table(\n    fusion_oof: pd.DataFrame,\n    canonical_meta: pd.DataFrame,\n    ranked_attr_df: pd.DataFrame,\n    topk_list: List[int],\n):\n    rows = []\n\n    meta_train = canonical_meta[canonical_meta[\"source\"] == \"train\"].copy()\n    meta_train = meta_train.rename(columns={\"dataset_name\": \"dataset\"})\n    keep_meta = [c for c in [\"dataset\", \"repertoire_id\", \"sex\", \"race\", \"age\"] if c in meta_train.columns]\n    if \"repertoire_id\" in keep_meta:\n        meta_use = meta_train[keep_meta].copy().rename(columns={\"repertoire_id\": \"ID\"})\n    else:\n        meta_use = pd.DataFrame(columns=[\"ID\", \"dataset\"])\n\n    fx = fusion_oof.merge(meta_use, on=[\"ID\", \"dataset\"], how=\"left\")\n\n    method_cols = [c for c in [\"SIMPLE_MEAN\", \"RANK_MEAN\", \"STACK_GLOBAL\", \"STACK_DATASET\"] if c in fx.columns]\n\n    for ds_name, g in fx.groupby(\"dataset\"):\n        y = g[\"label_positive\"].values\n        for m in method_cols:\n            rows.append({\n                \"benchmark\": \"real\",\n                \"scope\": \"dataset\",\n                \"dataset\": ds_name,\n                \"metric\": \"auc\",\n                \"method\": m,\n                \"value\": safe_auc(y, g[m].values),\n            })\n            rows.append({\n                \"benchmark\": \"real\",\n                \"scope\": \"dataset\",\n                \"dataset\": ds_name,\n                \"metric\": \"brier\",\n                \"method\": m,\n                \"value\": brier_score(y, g[m].values),\n            })\n            rows.append({\n                \"benchmark\": \"real\",\n                \"scope\": \"dataset\",\n                \"dataset\": ds_name,\n                \"metric\": \"ece\",\n                \"method\": m,\n                \"value\": expected_calibration_error(y, g[m].values),\n            })\n\n    y_all = fx[\"label_positive\"].values\n    for m in method_cols:\n        rows.append({\n            \"benchmark\": \"real\",\n            \"scope\": \"global\",\n            \"dataset\": \"ALL\",\n            \"metric\": \"auc\",\n            \"method\": m,\n            \"value\": safe_auc(y_all, fx[m].values),\n        })\n        rows.append({\n            \"benchmark\": \"real\",\n            \"scope\": \"global\",\n            \"dataset\": \"ALL\",\n            \"metric\": \"brier\",\n            \"method\": m,\n            \"value\": brier_score(y_all, fx[m].values),\n        })\n        rows.append({\n            \"benchmark\": \"real\",\n            \"scope\": \"global\",\n            \"dataset\": \"ALL\",\n            \"metric\": \"ece\",\n            \"method\": m,\n            \"value\": expected_calibration_error(y_all, fx[m].values),\n        })\n\n    # subgroup performance if metadata exists\n    for subgroup_col in [\"sex\", \"race\"]:\n        if subgroup_col not in fx.columns:\n            continue\n        tmp = fx[pd.notna(fx[subgroup_col])].copy()\n        if len(tmp) == 0:\n            continue\n        for subgroup_val, g in tmp.groupby(subgroup_col):\n            if len(g) < 12:\n                continue\n            for m in method_cols:\n                rows.append({\n                    \"benchmark\": \"real\",\n                    \"scope\": f\"subgroup_{subgroup_col}\",\n                    \"dataset\": str(subgroup_val),\n                    \"metric\": \"auc\",\n                    \"method\": m,\n                    \"value\": safe_auc(g[\"label_positive\"].values, g[m].values),\n                })\n\n    # ranking proxy from phase6 attribution outputs\n    if len(ranked_attr_df):\n        if \"weak_target\" in ranked_attr_df.columns:\n            z = ranked_attr_df.copy()\n            z[\"weak_positive\"] = (pd.to_numeric(z[\"weak_target\"], errors=\"coerce\").fillna(0.0) > 0).astype(int)\n            for ds_name, g in z.groupby(\"dataset\"):\n                for k in topk_list:\n                    topk = g[g[\"rank\"] <= k].copy()\n                    if len(topk) == 0:\n                        continue\n                    rep_level = (\n                        topk.groupby(\"ID\", as_index=False)[\"weak_positive\"]\n                        .mean()\n                        .rename(columns={\"weak_positive\": \"support\"})\n                    )\n                    rows.append({\n                        \"benchmark\": \"real\",\n                        \"scope\": \"dataset\",\n                        \"dataset\": ds_name,\n                        \"metric\": f\"weak_support_at_{k}\",\n                        \"method\": \"ATTRIBUTION_PROXY\",\n                        \"value\": float(rep_level[\"support\"].mean()),\n                    })\n\n            for k in topk_list:\n                topk = z[z[\"rank\"] <= k].copy()\n                if len(topk) == 0:\n                    continue\n                rep_level = (\n                    topk.groupby(\"ID\", as_index=False)[\"weak_positive\"]\n                    .mean()\n                    .rename(columns={\"weak_positive\": \"support\"})\n                )\n                rows.append({\n                    \"benchmark\": \"real\",\n                    \"scope\": \"global\",\n                    \"dataset\": \"ALL\",\n                    \"metric\": f\"weak_support_at_{k}\",\n                    \"method\": \"ATTRIBUTION_PROXY\",\n                    \"value\": float(rep_level[\"support\"].mean()),\n                })\n\n    return pd.DataFrame(rows)\n\"\"\")\n\ntrain_benchmarks_py = dedent(\"\"\"\nfrom __future__ import annotations\n\nimport sys\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport yaml\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.model_selection import train_test_split\n\nPROJECT_ROOT = Path(__file__).resolve().parent\nif str(PROJECT_ROOT) not in sys.path:\n    sys.path.insert(0, str(PROJECT_ROOT))\n\nfrom src.branch_b_public import (\n    build_feature_table_from_meta,\n    build_fold_signal_catalog,\n)\nfrom src.evaluation import real_benchmark_table, safe_auc\nfrom src.simulation import (\n    estimate_empirical_background,\n    sequence_scores_from_bundle,\n    simulate_ligo_dataset,\n    simulate_simairr_dataset,\n)\nfrom src.utils import ensure_dir, seed_everything\n\n\ndef fit_lr_scores(train_feat: pd.DataFrame, val_feat: pd.DataFrame, feature_cols):\n    X_tr = train_feat[feature_cols].fillna(0.0).values\n    y_tr = train_feat[\"label_positive\"].astype(int).values\n    X_va = val_feat[feature_cols].fillna(0.0).values\n\n    clf = LogisticRegression(max_iter=3000, solver=\"liblinear\", class_weight=\"balanced\", random_state=42)\n    clf.fit(X_tr, y_tr)\n    pred = clf.predict_proba(X_va)[:, 1]\n    return pred\n\n\ndef build_synth_branch_b_features(meta_train, meta_val, cache, bundle):\n    tr_feat = build_feature_table_from_meta(meta_train, cache, bundle)\n    va_feat = build_feature_table_from_meta(meta_val, cache, bundle)\n    return tr_feat, va_feat\n\n\ndef run_simairr_benchmark(cfg, aa_probs, length_source):\n    rows = []\n\n    sim_cfg = cfg[\"simairr\"]\n    dict_cfg = cfg[\"dictionary\"]\n\n    meta, cache, truth = simulate_simairr_dataset(\n        n_pos=sim_cfg[\"n_pos\"],\n        n_neg=sim_cfg[\"n_neg\"],\n        n_background_per_rep=sim_cfg[\"n_background_per_rep\"],\n        n_public_background=sim_cfg[\"n_public_background\"],\n        n_disease_public=sim_cfg[\"n_disease_public\"],\n        public_common_prob=sim_cfg[\"public_common_prob\"],\n        disease_prob_pos=sim_cfg[\"disease_prob_pos\"],\n        disease_prob_neg=sim_cfg[\"disease_prob_neg\"],\n        aa_probs=aa_probs,\n        length_source=length_source,\n        random_state=sim_cfg[\"random_state\"],\n    )\n\n    tr_meta, va_meta = train_test_split(\n        meta,\n        test_size=0.30,\n        stratify=meta[\"label_positive\"],\n        random_state=sim_cfg[\"random_state\"],\n    )\n\n    bundle, exact_df, cluster_df = build_fold_signal_catalog(\n        fold_train_meta=tr_meta,\n        cache=cache,\n        min_freq=dict_cfg[\"min_freq\"],\n        enrichment=dict_cfg[\"enrichment\"],\n        top_exact=dict_cfg[\"top_exact\"],\n        top_clusters=dict_cfg[\"top_clusters\"],\n        max_repertoires_per_class=dict_cfg[\"max_repertoires_per_class\"],\n        random_state=sim_cfg[\"random_state\"],\n    )\n\n    tr_feat, va_feat = build_synth_branch_b_features(tr_meta, va_meta, cache, bundle)\n    y_va = va_feat[\"label_positive\"].astype(int).values\n\n    feature_cols = [c for c in tr_feat.columns if c not in [\"ID\", \"dataset\", \"label_positive\"]]\n\n    pred_full = fit_lr_scores(tr_feat, va_feat, feature_cols)\n    auc_full = safe_auc(y_va, pred_full)\n\n    exact_signal = va_feat[\"public_seq_weighted_score\"].values if \"public_seq_weighted_score\" in va_feat.columns else np.zeros(len(va_feat))\n    cluster_signal = va_feat[\"cluster_enrichment_sum\"].values if \"cluster_enrichment_sum\" in va_feat.columns else np.zeros(len(va_feat))\n    combo_signal = exact_signal + cluster_signal\n\n    auc_exact = safe_auc(y_va, exact_signal)\n    auc_cluster = safe_auc(y_va, cluster_signal)\n    auc_combo = safe_auc(y_va, combo_signal)\n\n    neg_mask = y_va == 0\n    pos_mask = y_va == 1\n\n    neg_public_hit_rate = float((va_feat.loc[neg_mask, \"public_seq_count\"] > 0).mean()) if neg_mask.sum() > 0 else np.nan\n    pos_public_hit_rate = float((va_feat.loc[pos_mask, \"public_seq_count\"] > 0).mean()) if pos_mask.sum() > 0 else np.nan\n    neg_cluster_hit_rate = float((va_feat.loc[neg_mask, \"cluster_hit_count\"] > 0).mean()) if neg_mask.sum() > 0 else np.nan\n    pos_cluster_hit_rate = float((va_feat.loc[pos_mask, \"cluster_hit_count\"] > 0).mean()) if pos_mask.sum() > 0 else np.nan\n\n    metrics = {\n        \"branch_b_full_auc\": auc_full,\n        \"exact_only_auc\": auc_exact,\n        \"cluster_only_auc\": auc_cluster,\n        \"exact_plus_cluster_auc\": auc_combo,\n        \"branch_b_advantage_over_exact\": auc_full - auc_exact if pd.notna(auc_full) and pd.notna(auc_exact) else np.nan,\n        \"branch_b_advantage_over_combo\": auc_full - auc_combo if pd.notna(auc_full) and pd.notna(auc_combo) else np.nan,\n        \"neg_public_hit_rate\": neg_public_hit_rate,\n        \"pos_public_hit_rate\": pos_public_hit_rate,\n        \"neg_cluster_hit_rate\": neg_cluster_hit_rate,\n        \"pos_cluster_hit_rate\": pos_cluster_hit_rate,\n        \"disease_catalog_size\": float(len(exact_df)),\n        \"cluster_catalog_size\": float(len(cluster_df)),\n    }\n\n    for metric, value in metrics.items():\n        rows.append({\n            \"benchmark\": \"simairr\",\n            \"scenario\": \"realistic_public_sharing\",\n            \"metric\": metric,\n            \"value\": value,\n        })\n\n    return pd.DataFrame(rows)\n\n\ndef run_ligo_benchmark(cfg, aa_probs, length_source):\n    rows = []\n\n    ligo_cfg = cfg[\"ligo\"]\n    dict_cfg = cfg[\"dictionary\"]\n\n    meta, cache, truth_df = simulate_ligo_dataset(\n        n_pos=ligo_cfg[\"n_pos\"],\n        n_neg=ligo_cfg[\"n_neg\"],\n        n_background_per_rep=ligo_cfg[\"n_background_per_rep\"],\n        implanted_motifs=ligo_cfg[\"implanted_motifs\"],\n        implants_per_positive=ligo_cfg[\"implants_per_positive\"],\n        implant_prob_positive=ligo_cfg[\"implant_prob_positive\"],\n        implant_prob_negative=ligo_cfg[\"implant_prob_negative\"],\n        aa_probs=aa_probs,\n        length_source=length_source,\n        random_state=ligo_cfg[\"random_state\"],\n    )\n\n    tr_meta, va_meta = train_test_split(\n        meta,\n        test_size=0.30,\n        stratify=meta[\"label_positive\"],\n        random_state=ligo_cfg[\"random_state\"],\n    )\n\n    # repeated catalogs for stability\n    bundle_list = []\n    seed_list = [42, 52, 62]\n    for seed in seed_list:\n        bundle, _, _ = build_fold_signal_catalog(\n            fold_train_meta=tr_meta,\n            cache=cache,\n            min_freq=dict_cfg[\"min_freq\"],\n            enrichment=dict_cfg[\"enrichment\"],\n            top_exact=dict_cfg[\"top_exact\"],\n            top_clusters=dict_cfg[\"top_clusters\"],\n            max_repertoires_per_class=dict_cfg[\"max_repertoires_per_class\"],\n            random_state=seed,\n        )\n        bundle_list.append((seed, bundle))\n\n    # AUROC on first bundle with branch-B-style features\n    main_bundle = bundle_list[0][1]\n    tr_feat = build_feature_table_from_meta(tr_meta, cache, main_bundle)\n    va_feat = build_feature_table_from_meta(va_meta, cache, main_bundle)\n\n    feature_cols = [c for c in tr_feat.columns if c not in [\"ID\", \"dataset\", \"label_positive\"]]\n    pred = fit_lr_scores(tr_feat, va_feat, feature_cols)\n    auc = safe_auc(va_feat[\"label_positive\"].astype(int).values, pred)\n\n    rows.append({\n        \"benchmark\": \"ligo\",\n        \"scenario\": \"low_witness_implant\",\n        \"metric\": \"repertoire_auc\",\n        \"value\": auc,\n    })\n\n    # sequence recovery@k\n    pos_val_ids = set(va_meta.loc[va_meta[\"label_positive\"] == 1, \"repertoire_id\"].astype(str).tolist())\n    truth_pos = truth_df[truth_df[\"repertoire_id\"].astype(str).isin(pos_val_ids)].copy()\n\n    recovery_by_k = {k: [] for k in ligo_cfg[\"recovery_k_list\"]}\n    topk_sets_by_seed = {seed: {} for seed in seed_list}\n\n    for seed, bundle in bundle_list:\n        for rep_id in sorted(pos_val_ids):\n            seq_counter = cache[rep_id][\"seq_counter\"]\n            score_df = sequence_scores_from_bundle(seq_counter, bundle)\n\n            if len(score_df) == 0:\n                continue\n\n            true_implants = set(truth_pos.loc[truth_pos[\"repertoire_id\"] == rep_id, \"sequence\"].astype(str).tolist())\n            if len(true_implants) == 0:\n                continue\n\n            top10 = set(score_df.head(10)[\"sequence\"].astype(str).tolist())\n            topk_sets_by_seed[seed][rep_id] = top10\n\n            for k in ligo_cfg[\"recovery_k_list\"]:\n                pred_topk = set(score_df.head(k)[\"sequence\"].astype(str).tolist())\n                rec = len(pred_topk & true_implants) / max(1, len(true_implants))\n                recovery_by_k[k].append(rec)\n\n    for k, vals in recovery_by_k.items():\n        rows.append({\n            \"benchmark\": \"ligo\",\n            \"scenario\": \"low_witness_implant\",\n            \"metric\": f\"recovery_at_{k}\",\n            \"value\": float(np.mean(vals)) if len(vals) else np.nan,\n        })\n\n    # attribution stability: average pairwise jaccard of top10 across seeds\n    jac_vals = []\n    rep_ids_common = sorted(pos_val_ids)\n    for rep_id in rep_ids_common:\n        sets = []\n        for seed in seed_list:\n            if rep_id in topk_sets_by_seed[seed]:\n                sets.append(topk_sets_by_seed[seed][rep_id])\n\n        if len(sets) < 2:\n            continue\n\n        for i in range(len(sets)):\n            for j in range(i + 1, len(sets)):\n                a, b = sets[i], sets[j]\n                denom = len(a | b)\n                if denom == 0:\n                    continue\n                jac_vals.append(len(a & b) / denom)\n\n    rows.append({\n        \"benchmark\": \"ligo\",\n        \"scenario\": \"low_witness_implant\",\n        \"metric\": \"attribution_stability_top10_jaccard\",\n        \"value\": float(np.mean(jac_vals)) if len(jac_vals) else np.nan,\n    })\n\n    return pd.DataFrame(rows)\n\n\ndef main():\n    cfg = yaml.safe_load((PROJECT_ROOT / \"configs\" / \"benchmark.yaml\").read_text())\n    out_dir = ensure_dir(cfg[\"paths\"][\"output_root\"])\n\n    seed_everything(42)\n\n    canonical = pd.read_csv(cfg[\"paths\"][\"canonical_metadata\"])\n    fusion_oof = pd.read_csv(cfg[\"paths\"][\"fusion_oof\"])\n    phase6_dir = Path(cfg[\"paths\"][\"phase6_dir\"])\n\n    from src.evaluation import load_ranked_attribution_tables\n    ranked_attr = load_ranked_attribution_tables(phase6_dir)\n\n    train_meta = canonical[canonical[\"source\"] == \"train\"].copy()\n    train_root = Path(\"/kaggle/input/competitions/adaptive-immune-profiling-challenge-2025/train_datasets/train_datasets\")\n\n    aa_probs, length_source = estimate_empirical_background(\n        canonical_train_meta=train_meta,\n        train_root=train_root,\n        max_repertoires=16,\n        max_sequences_per_file=4000,\n        random_state=42,\n    )\n\n    # Benchmark 1: real data (OOF + attribution proxy)\n    real_df = real_benchmark_table(\n        fusion_oof=fusion_oof,\n        canonical_meta=canonical,\n        ranked_attr_df=ranked_attr,\n        topk_list=cfg[\"real\"][\"topk_list\"],\n    )\n    real_df.to_csv(out_dir / \"benchmark_real.csv\", index=False)\n\n    # Benchmark 2: simAIRR-style\n    sim_df = run_simairr_benchmark(cfg, aa_probs=aa_probs, length_source=length_source)\n    sim_df.to_csv(out_dir / \"benchmark_simairr.csv\", index=False)\n\n    # Benchmark 3: LIgO-style implanted signal\n    ligo_df = run_ligo_benchmark(cfg, aa_probs=aa_probs, length_source=length_source)\n    ligo_df.to_csv(out_dir / \"benchmark_ligo.csv\", index=False)\n\n    print(\"=\" * 80)\n    print(\"PHASE 8: Benchmarking section\")\n    print(\"=\" * 80)\n\n    print(\"\\\\nReal benchmark:\")\n    print(real_df.head(20).to_string(index=False))\n\n    print(\"\\\\nsimAIRR benchmark:\")\n    print(sim_df.to_string(index=False))\n\n    print(\"\\\\nLIgO benchmark:\")\n    print(ligo_df.to_string(index=False))\n\n    print(\"\\\\nSaved:\")\n    print(out_dir / \"benchmark_real.csv\")\n    print(out_dir / \"benchmark_simairr.csv\")\n    print(out_dir / \"benchmark_ligo.csv\")\n    print(\"=\" * 80)\n\n\nif __name__ == \"__main__\":\n    main()\n\"\"\")\n\n(PROJECT_ROOT / \"src\" / \"simulation.py\").write_text(simulation_py, encoding=\"utf-8\")\n(PROJECT_ROOT / \"src\" / \"evaluation.py\").write_text(evaluation_py, encoding=\"utf-8\")\n(PROJECT_ROOT / \"train_benchmarks.py\").write_text(train_benchmarks_py, encoding=\"utf-8\")\n\nprint(\"Written:\")\nprint(\"-\", PROJECT_ROOT / \"src\" / \"simulation.py\")\nprint(\"-\", PROJECT_ROOT / \"src\" / \"evaluation.py\")\nprint(\"-\", PROJECT_ROOT / \"train_benchmarks.py\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:29:59.697389Z","iopub.execute_input":"2026-07-24T11:29:59.697594Z","iopub.status.idle":"2026-07-24T11:29:59.717347Z","shell.execute_reply.started":"2026-07-24T11:29:59.697576Z","shell.execute_reply":"2026-07-24T11:29:59.716492Z"}},"outputs":[],"execution_count":null},{"id":"53a56c9f","cell_type":"code","source":"# Cell 3: Benchmark imports verify\n\nfrom pathlib import Path\np = Path(\"/kaggle/working/project/src\")\nprint(\"src/ exists:\", p.exists())\nprint(\"contents:\", sorted(f.name for f in p.glob(\"*\")) if p.exists() else \"N/A\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:29:59.718753Z","iopub.execute_input":"2026-07-24T11:29:59.719088Z","iopub.status.idle":"2026-07-24T11:29:59.733250Z","shell.execute_reply.started":"2026-07-24T11:29:59.719068Z","shell.execute_reply":"2026-07-24T11:29:59.732675Z"}},"outputs":[],"execution_count":null},{"id":"336ddaf9","cell_type":"code","source":"# Cell 4: Run Phase 8 benchmarking\n\n!python /kaggle/working/project/train_benchmarks.py","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:29:59.734197Z","iopub.execute_input":"2026-07-24T11:29:59.734516Z","iopub.status.idle":"2026-07-24T11:30:01.509542Z","shell.execute_reply.started":"2026-07-24T11:29:59.734496Z","shell.execute_reply":"2026-07-24T11:30:01.508807Z"}},"outputs":[],"execution_count":null},{"id":"f1a7a183","cell_type":"code","source":"# Cell 5: Inspect benchmark outputs\n\n#Cell 5 — outputs inspect করো\nfrom pathlib import Path\nimport pandas as pd\n\nOUT_DIR = Path(\"/kaggle/working/project/artifacts/phase3_branch_b\")\n\nprint(\"OUT_DIR exists:\", OUT_DIR.exists())\nprint(\"Files:\")\nfound = sorted(OUT_DIR.glob(\"*\")) if OUT_DIR.exists() else []\nfor p in found:\n    print(\"-\", p.name)\nif not found:\n    print(\"(none — Branch B training/inference did not write any outputs here.)\")\n\ndef show(name, path):\n    print(f\"\\n{name}\")\n    if path.exists():\n        display(pd.read_csv(path).head())\n    else:\n        print(f\"  -> missing: {path}\")\n\nshow(\"branch_b_oof.csv\", OUT_DIR / \"branch_b_oof.csv\")\n\nsummary_path = OUT_DIR / \"branch_b_summary.csv\"\nif summary_path.exists():\n    print(\"\\nbranch_b_summary.csv\")\n    display(pd.read_csv(summary_path))\nelse:\n    print(\"\\nbranch_b_summary.csv -> missing\")\n\nexact_files = sorted(OUT_DIR.glob(\"disease_enriched_sequences_*.csv\"))\nif exact_files:\n    print(f\"\\n{exact_files[0].name}\")\n    display(pd.read_csv(exact_files[0]).head())\nelse:\n    print(\"\\ndisease_enriched_sequences_*.csv -> none found\")\n\ncluster_files = sorted(OUT_DIR.glob(\"cluster_catalog_*.csv\"))\nif cluster_files:\n    print(f\"\\n{cluster_files[0].name}\")\n    display(pd.read_csv(cluster_files[0]).head())\nelse:\n    print(\"\\ncluster_catalog_*.csv -> none found\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-24T11:30:01.510762Z","iopub.execute_input":"2026-07-24T11:30:01.511206Z","iopub.status.idle":"2026-07-24T11:30:01.520887Z","shell.execute_reply.started":"2026-07-24T11:30:01.511176Z","shell.execute_reply":"2026-07-24T11:30:01.520267Z"}},"outputs":[],"execution_count":null}]}