diff --git a/classes/v2/loader_factory.py b/classes/v2/loader_factory.py
index 734b2ec..e515ed7 100644
--- a/classes/v2/loader_factory.py
+++ b/classes/v2/loader_factory.py
@@ -20,6 +20,13 @@ def _default_slot_descriptors(patient_col: str, label_col: str) -> dict[str, Slo
             required=True,
             shape_hint="scalar",
         ),
+        "eye_id_1": SlotDescriptor(
+            key="eye_id_1",
+            kind="id",
+            description="Eye side identifier (OD/OS)",
+            required=False,
+            shape_hint="scalar",
+        ),
         "label_1": SlotDescriptor(
             key="label_1",
             kind="label",
@@ -53,6 +60,7 @@ def _row_to_sample(
 ) -> dict[str, Any]:
     return {
         "id_1": row[patient_col],
+        "eye_id_1": str(row.get("eyeID", "")),
         "label_1": row[label_col],
         "image_1": clinical.get_image_path(row) if hasattr(clinical, "get_image_path") else None,
         "matrix_1": clinical.vectorize_row(row) if hasattr(clinical, "vectorize_row") else None,
diff --git a/classes/v2/models.py b/classes/v2/models.py
index 5744951..26d443a 100644
--- a/classes/v2/models.py
+++ b/classes/v2/models.py
@@ -34,6 +34,7 @@ class SingleEyeHT(nn.Module):
         num_classes: int,
         md_hidden_dim: int = 128,
         fusion_dim: int = 256,
+        bridge_mode: str = "fused",
     ):
         super().__init__()
         self.img_tower = ImageTower(
@@ -52,7 +53,7 @@ class SingleEyeHT(nn.Module):
             meta_dim=self.md_tower.out_dim,
             num_classes=num_classes,
             fusion_dim=fusion_dim,
-            mode="fused",
+            mode=bridge_mode,
             use_se=False,
         )
 
@@ -61,8 +62,8 @@ class SingleEyeHT(nn.Module):
         return self.img_tower.transform
 
     def forward(self, x: torch.Tensor, meta: torch.Tensor) -> torch.Tensor:
-        img_feats = self.img_tower(x)
-        md_feats  = self.md_tower(meta)
+        img_feats = None if self.bridge.mode == "metadata_only" else self.img_tower(x)
+        md_feats  = None if self.bridge.mode == "image_only"    else self.md_tower(meta)
         out_f, _, _ = self.bridge(img_feats, md_feats)
         return out_f
 
@@ -214,11 +215,15 @@ def _set_requires_grad(module: nn.Module, enabled: bool) -> None:
 
 
 def _set_single_phase(model: SingleEyeHT, phase: str) -> None:
+    bridge_mode = model.bridge.mode
+    # Ablation modes have no fusion bridge; fused_warmup is meaningless — treat as tower_warmup
+    if bridge_mode in ("image_only", "metadata_only") and phase == "fused_warmup":
+        phase = "tower_warmup"
     if phase == "tower_warmup":
-        _set_requires_grad(model.img_tower, True)
-        _set_requires_grad(model.md_tower, True)
-        _set_requires_grad(model.bridge.classifier_img, True)
-        _set_requires_grad(model.bridge.classifier_md, True)
+        _set_requires_grad(model.img_tower, bridge_mode != "metadata_only")
+        _set_requires_grad(model.md_tower, bridge_mode != "image_only")
+        _set_requires_grad(model.bridge.classifier_img, bridge_mode != "metadata_only")
+        _set_requires_grad(model.bridge.classifier_md, bridge_mode != "image_only")
         _set_requires_grad(model.bridge.W_img, False)
         _set_requires_grad(model.bridge.W_md, False)
         _set_requires_grad(model.bridge.classifier_fused, False)
@@ -282,19 +287,33 @@ def train_single_epoch(
         x = x.to(device)
         m = m.to(device)
         y = _to_label_tensor(y, device)
-        img_feats = model.img_tower(x)
-        md_feats = model.md_tower(m)
+        bridge_mode = model.bridge.mode
+        img_feats = None if bridge_mode == "metadata_only" else model.img_tower(x)
+        md_feats  = None if bridge_mode == "image_only"    else model.md_tower(m)
 
         if phase == "tower_warmup":
-            logits_i = model.bridge.classifier_img(img_feats)
-            logits_m = model.bridge.classifier_md(md_feats)
-            loss = 0.5 * (F.cross_entropy(logits_i, y) + F.cross_entropy(logits_m, y))
-            logits = 0.5 * (F.softmax(logits_i, dim=1) + F.softmax(logits_m, dim=1))
+            if bridge_mode == "metadata_only":
+                logits = model.bridge.classifier_md(md_feats)
+                loss = F.cross_entropy(logits, y)
+            elif bridge_mode == "image_only":
+                logits = model.bridge.classifier_img(img_feats)
+                loss = F.cross_entropy(logits, y)
+            else:
+                logits_i = model.bridge.classifier_img(img_feats)
+                logits_m = model.bridge.classifier_md(md_feats)
+                loss = 0.5 * (F.cross_entropy(logits_i, y) + F.cross_entropy(logits_m, y))
+                logits = 0.5 * (F.softmax(logits_i, dim=1) + F.softmax(logits_m, dim=1))
         elif phase == "fused_warmup":
             logits, _, _ = model.bridge(img_feats, md_feats)
             loss = F.cross_entropy(logits, y)
         else:
-            if random() < bcd_prob:
+            if bridge_mode == "metadata_only":
+                logits = model.bridge.classifier_md(md_feats)
+                loss = F.cross_entropy(logits, y)
+            elif bridge_mode == "image_only":
+                logits = model.bridge.classifier_img(img_feats)
+                loss = F.cross_entropy(logits, y)
+            elif random() < bcd_prob:
                 if random() < 0.5:
                     logits = model.bridge.classifier_img(img_feats)
                 else:
@@ -447,6 +466,79 @@ def collect_probs_classic(
     return np.concatenate(y_chunks), np.concatenate(p_chunks, axis=0)
 
 
+def collect_probs_ensemble_pereye(
+    model: "SingleEyeHT",
+    loader: DataLoader,
+    device: torch.device,
+    *,
+    return_ids: bool = False,
+):
+    """
+    Per-patient, per-eye probs for all 3 heads from a bilateral loader (ensemble mode).
+
+    OD corresponds to image_1/matrix_1; OS to image_2/matrix_2.
+    Arrays are in patient order (not interleaved at sample level).
+
+    Returns:
+        (y, pf_od, pi_od, pm_od, pf_os, pi_os, pm_os)
+        or, when return_ids=True:
+        (y, pf_od, pi_od, pm_od, pf_os, pi_os, pm_os, patient_ids)
+
+    Patient-level averaged ensemble probs can be recovered as:
+        p_en = 0.5 * (pf_od + pf_os)
+    """
+    model.eval()
+    y_chunks: list = []
+    pf_od_c, pi_od_c, pm_od_c = [], [], []
+    pf_os_c, pi_os_c, pm_os_c = [], [], []
+    id_chunks: list[str] = []
+
+    with torch.no_grad():
+        for batch in loader:
+            x1 = batch.get("image_1"); m1 = batch.get("matrix_1")
+            x2 = batch.get("image_2"); m2 = batch.get("matrix_2")
+            y  = batch.get("label_1")
+            if not (torch.is_tensor(x1) and torch.is_tensor(m1) and
+                    torch.is_tensor(x2) and torch.is_tensor(m2)):
+                continue
+            y_t = _to_label_tensor(y, device)
+
+            def _fwd(x, m):
+                img_feats = None if model.bridge.mode == "metadata_only" else model.img_tower(x.to(device))
+                md_feats  = None if model.bridge.mode == "image_only"    else model.md_tower(m.to(device))
+                out_f, out_i, out_m = model.bridge(img_feats, md_feats)
+                pf = F.softmax(out_f, dim=1)
+                pi = F.softmax(out_i, dim=1) if out_i is not None else pf
+                pm = F.softmax(out_m, dim=1) if out_m is not None else pf
+                return pf, pi, pm
+
+            pf_od, pi_od, pm_od = _fwd(x1, m1)
+            pf_os, pi_os, pm_os = _fwd(x2, m2)
+
+            y_chunks.append(y_t.cpu().numpy())
+            pf_od_c.append(pf_od.cpu().numpy()); pi_od_c.append(pi_od.cpu().numpy()); pm_od_c.append(pm_od.cpu().numpy())
+            pf_os_c.append(pf_os.cpu().numpy()); pi_os_c.append(pi_os.cpu().numpy()); pm_os_c.append(pm_os.cpu().numpy())
+
+            if return_ids:
+                ids = batch.get("id_1", [""] * len(y_t))
+                if torch.is_tensor(ids):
+                    ids = ids.tolist()
+                id_chunks.extend([str(i) for i in ids])
+
+    if not y_chunks:
+        z = np.zeros((0, 0), dtype=np.float32)
+        empty_i = np.array([], dtype=np.int64)
+        base = (empty_i, z, z, z, z, z, z)
+        return base + (np.array([], dtype=object),) if return_ids else base
+
+    y     = np.concatenate(y_chunks)
+    pf_od = np.concatenate(pf_od_c, axis=0); pi_od = np.concatenate(pi_od_c, axis=0); pm_od = np.concatenate(pm_od_c, axis=0)
+    pf_os = np.concatenate(pf_os_c, axis=0); pi_os = np.concatenate(pi_os_c, axis=0); pm_os = np.concatenate(pm_os_c, axis=0)
+    if return_ids:
+        return y, pf_od, pi_od, pm_od, pf_os, pi_os, pm_os, np.array(id_chunks, dtype=object)
+    return y, pf_od, pi_od, pm_od, pf_os, pi_os, pm_os
+
+
 def collect_probs_ensemble(
     model: SingleEyeHT,
     loader: DataLoader,
@@ -532,15 +624,22 @@ def collect_probs_single_components(
     device: torch.device,
     *,
     aggregate_patient: bool,
-) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    return_logits: bool = False,
+):
     """
-    Collect fused/img/md probabilities for SingleEyeHT.
+    Collect fused/img/md probabilities (and optionally raw logits) for SingleEyeHT.
     - aggregate_patient=False: eye-level (OD/OS as independent samples)
     - aggregate_patient=True : patient-level (average OD/OS per head)
+    - return_logits=False: returns (y, probs_f, probs_i, probs_m)
+    - return_logits=True:  returns (y, probs_f, probs_i, probs_m,
+                                       logits_f, logits_i, logits_m)
+      Note: logits are averaged across eyes when aggregate_patient=True,
+      which is equivalent to averaging in logit space (before softmax).
     """
     model.eval()
     y_chunks = []
     pf_chunks, pi_chunks, pm_chunks = [], [], []
+    lf_chunks, li_chunks, lm_chunks = [], [], []
     with torch.no_grad():
         for batch in loader:
             x1 = batch.get("image_1"); m1 = batch.get("matrix_1")
@@ -550,40 +649,118 @@ def collect_probs_single_components(
                 continue
             y_t = _to_label_tensor(y, device)
 
-            def _per_eye_probs(x, m):
-                img_feats = model.img_tower(x.to(device))
-                md_feats = model.md_tower(m.to(device))
+            def _per_eye(x, m):
+                img_feats = None if model.bridge.mode == "metadata_only" else model.img_tower(x.to(device))
+                md_feats  = None if model.bridge.mode == "image_only"    else model.md_tower(m.to(device))
                 out_f, out_i, out_m = model.bridge(img_feats, md_feats)
-                return (
-                    F.softmax(out_f, dim=1),
-                    F.softmax(out_i, dim=1),
-                    F.softmax(out_m, dim=1),
-                )
+                pf = F.softmax(out_f, dim=1)
+                pi = F.softmax(out_i, dim=1) if out_i is not None else pf
+                pm = F.softmax(out_m, dim=1) if out_m is not None else pf
+                lf = out_f
+                li = out_i if out_i is not None else out_f
+                lm = out_m if out_m is not None else out_f
+                return pf, pi, pm, lf, li, lm
 
-            pf_od, pi_od, pm_od = _per_eye_probs(x1, m1)
-            pf_os, pi_os, pm_os = _per_eye_probs(x2, m2)
+            pf_od, pi_od, pm_od, lf_od, li_od, lm_od = _per_eye(x1, m1)
+            pf_os, pi_os, pm_os, lf_os, li_os, lm_os = _per_eye(x2, m2)
 
             if aggregate_patient:
                 y_chunks.append(y_t.cpu().numpy())
                 pf_chunks.append((0.5 * (pf_od + pf_os)).cpu().numpy())
                 pi_chunks.append((0.5 * (pi_od + pi_os)).cpu().numpy())
                 pm_chunks.append((0.5 * (pm_od + pm_os)).cpu().numpy())
+                lf_chunks.append((0.5 * (lf_od + lf_os)).cpu().numpy())
+                li_chunks.append((0.5 * (li_od + li_os)).cpu().numpy())
+                lm_chunks.append((0.5 * (lm_od + lm_os)).cpu().numpy())
             else:
                 y_np = y_t.cpu().numpy()
                 y_chunks += [y_np, y_np]
                 pf_chunks += [pf_od.cpu().numpy(), pf_os.cpu().numpy()]
                 pi_chunks += [pi_od.cpu().numpy(), pi_os.cpu().numpy()]
                 pm_chunks += [pm_od.cpu().numpy(), pm_os.cpu().numpy()]
+                lf_chunks += [lf_od.cpu().numpy(), lf_os.cpu().numpy()]
+                li_chunks += [li_od.cpu().numpy(), li_os.cpu().numpy()]
+                lm_chunks += [lm_od.cpu().numpy(), lm_os.cpu().numpy()]
 
     if not y_chunks:
         z = np.zeros((0, 0), dtype=np.float32)
+        if return_logits:
+            return np.array([], dtype=np.int64), z, z, z, z, z, z
         return np.array([], dtype=np.int64), z, z, z
-    return (
-        np.concatenate(y_chunks),
-        np.concatenate(pf_chunks, axis=0),
-        np.concatenate(pi_chunks, axis=0),
-        np.concatenate(pm_chunks, axis=0),
-    )
+
+    y   = np.concatenate(y_chunks)
+    pf  = np.concatenate(pf_chunks, axis=0)
+    pi  = np.concatenate(pi_chunks, axis=0)
+    pm  = np.concatenate(pm_chunks, axis=0)
+    if return_logits:
+        lf = np.concatenate(lf_chunks, axis=0)
+        li = np.concatenate(li_chunks, axis=0)
+        lm = np.concatenate(lm_chunks, axis=0)
+        return y, pf, pi, pm, lf, li, lm
+    return y, pf, pi, pm
+
+
+def collect_probs_eye_level(
+    model: "SingleEyeHT",
+    loader: DataLoader,
+    device: torch.device,
+    *,
+    return_ids: bool = False,
+):
+    """
+    Collect fused/img/md probabilities from a single-eye loader (image_1/matrix_1 only).
+    Used for eval-mode passes over the training set.
+
+    Returns (y, probs_f, probs_i, probs_m) or, when return_ids=True,
+    (y, probs_f, probs_i, probs_m, sample_ids) where sample_ids is an
+    array of strings like "2OD", "4OS".
+    """
+    model.eval()
+    y_chunks, pf_chunks, pi_chunks, pm_chunks, id_chunks = [], [], [], [], []
+    with torch.no_grad():
+        for batch in loader:
+            x = batch.get("image_1")
+            m = batch.get("matrix_1")
+            y = batch.get("label_1")
+            if not (torch.is_tensor(x) and torch.is_tensor(m)):
+                continue
+            y_t = _to_label_tensor(y, device)
+            img_feats = None if model.bridge.mode == "metadata_only" else model.img_tower(x.to(device))
+            md_feats  = None if model.bridge.mode == "image_only"    else model.md_tower(m.to(device))
+            out_f, out_i, out_m = model.bridge(img_feats, md_feats)
+            pf = F.softmax(out_f, dim=1)
+            pi = F.softmax(out_i, dim=1) if out_i is not None else pf
+            pm = F.softmax(out_m, dim=1) if out_m is not None else pf
+            y_chunks.append(y_t.cpu().numpy())
+            pf_chunks.append(pf.cpu().numpy())
+            pi_chunks.append(pi.cpu().numpy())
+            pm_chunks.append(pm.cpu().numpy())
+            if return_ids:
+                ids   = batch.get("id_1", [""] * len(y_t))
+                eyes  = batch.get("eye_id_1", [""] * len(y_t))
+                # ids/eyes may be tensors (int) or lists of strings
+                if torch.is_tensor(ids):
+                    ids = ids.tolist()
+                if torch.is_tensor(eyes):
+                    eyes = eyes.tolist()
+                id_chunks.extend(
+                    [f"{pid}{eye}" for pid, eye in zip(ids, eyes)]
+                )
+
+    if not y_chunks:
+        z = np.zeros((0, 0), dtype=np.float32)
+        empty_ids = np.array([], dtype=object)
+        if return_ids:
+            return np.array([], dtype=np.int64), z, z, z, empty_ids
+        return np.array([], dtype=np.int64), z, z, z
+
+    y  = np.concatenate(y_chunks)
+    pf = np.concatenate(pf_chunks, axis=0)
+    pi = np.concatenate(pi_chunks, axis=0)
+    pm = np.concatenate(pm_chunks, axis=0)
+    if return_ids:
+        return y, pf, pi, pm, np.array(id_chunks, dtype=object)
+    return y, pf, pi, pm
 
 
 def collect_probs_bilateral_components(
diff --git a/classes/v2/papila_builders.py b/classes/v2/papila_builders.py
index bcc2ce6..f787d93 100644
--- a/classes/v2/papila_builders.py
+++ b/classes/v2/papila_builders.py
@@ -50,14 +50,15 @@ def _correct_iop(raw_iop: float, pachy: float) -> float:
 
 
 def _apply_iop_and_drop_md(df: pd.DataFrame) -> pd.DataFrame:
-    """Add IOP_raw/IOP_corr and drop VF_MD if present (in-place safe)."""
+    """Add IOP_raw/IOP_corr and drop source IOP columns + VF_MD if present (in-place safe)."""
     df["IOP_raw"] = df.apply(_pick_iop, axis=1)
     pachy = df.get("Pachymetry", pd.Series(np.nan, index=df.index))
     df["IOP_corr"] = [
         _correct_iop(r, p) for r, p in zip(df["IOP_raw"].values, pachy.values)
     ]
-    if "VF_MD" in df.columns:
-        df.drop(columns=["VF_MD"], inplace=True)
+    drop_cols = [c for c in ("Pneumatic", "Perkins", "VF_MD") if c in df.columns]
+    if drop_cols:
+        df.drop(columns=drop_cols, inplace=True)
     return df
 
 
diff --git a/classes/v2/results.py b/classes/v2/results.py
index 4378f6f..28c4cdd 100644
--- a/classes/v2/results.py
+++ b/classes/v2/results.py
@@ -109,3 +109,23 @@ class FoldArtifacts:
     probs_bilat:     Optional[np.ndarray]
     y_true_fused:    Optional[np.ndarray] = None
     probs_fused:     Optional[np.ndarray] = None
+    probs_ensemble_img: Optional[np.ndarray] = None
+    probs_ensemble_md:  Optional[np.ndarray] = None
+    probs_classic_img:  Optional[np.ndarray] = None
+    probs_classic_md:   Optional[np.ndarray] = None
+    # per-eye (pre-averaged) versions for ensemble mode
+    y_true_ensemble_pereye:    Optional[np.ndarray] = None
+    probs_ensemble_pereye:     Optional[np.ndarray] = None
+    probs_ensemble_img_pereye: Optional[np.ndarray] = None
+    probs_ensemble_md_pereye:  Optional[np.ndarray] = None
+    # raw logits (before softmax) — patient-level
+    logits_ensemble:     Optional[np.ndarray] = None
+    logits_ensemble_img: Optional[np.ndarray] = None
+    logits_ensemble_md:  Optional[np.ndarray] = None
+    logits_classic:      Optional[np.ndarray] = None
+    logits_classic_img:  Optional[np.ndarray] = None
+    logits_classic_md:   Optional[np.ndarray] = None
+    # raw logits — per-eye
+    logits_ensemble_pereye:     Optional[np.ndarray] = None
+    logits_ensemble_img_pereye: Optional[np.ndarray] = None
+    logits_ensemble_md_pereye:  Optional[np.ndarray] = None
diff --git a/classes/v2/v2_hypertower.py b/classes/v2/v2_hypertower.py
index a6b8257..81baf80 100644
--- a/classes/v2/v2_hypertower.py
+++ b/classes/v2/v2_hypertower.py
@@ -39,6 +39,8 @@ from classes.v2.models import (
     collect_probs_bilateral_components,
     collect_probs_classic,
     collect_probs_ensemble,
+    collect_probs_ensemble_pereye,
+    collect_probs_eye_level,
     collect_probs_fused,
     collect_probs_single_components,
     train_bilateral_epoch,
@@ -59,6 +61,98 @@ from classes.v2.utils import (
 from classes.v2.hypertower_logger import HypertowerLogger
 
 
+# ---------------------------------------------------------------------------
+# Fusion-event helper
+# ---------------------------------------------------------------------------
+
+def _fusion_events(
+    y: np.ndarray,
+    pf: np.ndarray,
+    pi: np.ndarray,
+    pm: np.ndarray,
+) -> tuple[int, int]:
+    """
+    Count fusion corrections and errors.
+    - correction: fused correct, both img and md wrong
+    - error:      fused wrong, both img and md correct
+    Returns (n_corrections, n_errors).
+    """
+    pred_f = pf.argmax(1); pred_i = pi.argmax(1); pred_m = pm.argmax(1)
+    corr = int(((pred_f == y) & (pred_i != y) & (pred_m != y)).sum())
+    err  = int(((pred_f != y) & (pred_i == y) & (pred_m == y)).sum())
+    return corr, err
+
+
+def _cm_cells(y: np.ndarray, p: np.ndarray, num_classes: int) -> dict[str, int]:
+    """
+    Return confusion matrix cells as a flat dict.
+    Binary:     keys tn/fp/fn/tp
+    Multiclass: keys cm_{i}_{j} for true class i, predicted class j
+    Returns empty dict if arrays are empty or wrong shape.
+    """
+    if not y.size or p.ndim < 2 or p.shape[1] != num_classes:
+        return {}
+    pred = p.argmax(1)
+    if num_classes == 2:
+        tn = int(((pred == 0) & (y == 0)).sum())
+        fp = int(((pred == 1) & (y == 0)).sum())
+        fn = int(((pred == 0) & (y == 1)).sum())
+        tp = int(((pred == 1) & (y == 1)).sum())
+        return {"tn": tn, "fp": fp, "fn": fn, "tp": tp}
+    # multiclass: full NxN matrix
+    out: dict[str, int] = {}
+    for i in range(num_classes):
+        for j in range(num_classes):
+            out[f"cm_{i}_{j}"] = int(((y == i) & (pred == j)).sum())
+    return out
+
+
+# ---------------------------------------------------------------------------
+# Per-sample prediction logging
+# ---------------------------------------------------------------------------
+
+def _save_predictions_csv(
+    fold_dir: Path,
+    eval_mode: str,
+    y_true: np.ndarray,
+    heads: dict,          # {"fused": probs_array, "img": probs_array, "md": probs_array, ...}
+    suffix: str = "",     # e.g. "_pereye"
+) -> None:
+    """
+    Save a per-sample CSV with predicted class, per-class probabilities,
+    and TP/FP/TN/FN (binary) or correct flag (multiclass) for every head.
+    """
+    N = len(y_true)
+    num_classes = next(p.shape[1] for p in heads.values() if p is not None)
+    rows = []
+    for i in range(N):
+        true = int(y_true[i])
+        row: dict = {"idx": i, "y_true": true}
+        for head_name, probs in heads.items():
+            if probs is None:
+                continue
+            pred = int(probs[i].argmax())
+            row[f"pred_{head_name}"] = pred
+            for c in range(num_classes):
+                row[f"prob_{head_name}_c{c}"] = float(probs[i, c])
+            if eval_mode == "binary":
+                row[f"tp_{head_name}"] = int(pred == 1 and true == 1)
+                row[f"fp_{head_name}"] = int(pred == 1 and true == 0)
+                row[f"tn_{head_name}"] = int(pred == 0 and true == 0)
+                row[f"fn_{head_name}"] = int(pred == 0 and true == 1)
+            else:
+                row[f"correct_{head_name}"] = int(pred == true)
+        rows.append(row)
+
+    if not rows:
+        return
+    csv_path = fold_dir / f"predictions{suffix}.csv"
+    with csv_path.open("w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
+        writer.writeheader()
+        writer.writerows(rows)
+
+
 # ---------------------------------------------------------------------------
 # V2HyperTower
 # ---------------------------------------------------------------------------
@@ -156,6 +250,9 @@ class V2HyperTower:
                         help="MDTower hidden dimension.")
         ap.add_argument("--fusion-dim",    type=int, default=256,
                         help="Bridge/BilateralBridge fusion dimension.")
+        ap.add_argument("--bridge-mode", default="fused",
+                        choices=["fused", "image_only", "metadata_only"],
+                        help="Bridge fusion mode: fused (default), image_only, or metadata_only.")
         # Mixed patients
         ap.add_argument(
             "--exclude-mixed-patients",
@@ -332,10 +429,73 @@ class V2HyperTower:
                 np.save(fold_dir / "y_true.npy", artifacts.y_true_ensemble)
             if artifacts.probs_ensemble is not None:
                 np.save(fold_dir / "probs_fused.npy", artifacts.probs_ensemble)
+            if artifacts.probs_ensemble_img is not None:
+                np.save(fold_dir / "probs_img.npy", artifacts.probs_ensemble_img)
+            if artifacts.probs_ensemble_md is not None:
+                np.save(fold_dir / "probs_md.npy", artifacts.probs_ensemble_md)
             if artifacts.probs_classic is not None:
                 np.save(fold_dir / "probs_classic.npy", artifacts.probs_classic)
+            if artifacts.probs_classic_img is not None:
+                np.save(fold_dir / "probs_classic_img.npy", artifacts.probs_classic_img)
+            if artifacts.probs_classic_md is not None:
+                np.save(fold_dir / "probs_classic_md.npy", artifacts.probs_classic_md)
+            if artifacts.y_true_ensemble_pereye is not None:
+                np.save(fold_dir / "y_true_pereye.npy", artifacts.y_true_ensemble_pereye)
+            if artifacts.probs_ensemble_pereye is not None:
+                np.save(fold_dir / "probs_fused_pereye.npy", artifacts.probs_ensemble_pereye)
+            if artifacts.probs_ensemble_img_pereye is not None:
+                np.save(fold_dir / "probs_img_pereye.npy", artifacts.probs_ensemble_img_pereye)
+            if artifacts.probs_ensemble_md_pereye is not None:
+                np.save(fold_dir / "probs_md_pereye.npy", artifacts.probs_ensemble_md_pereye)
+            if artifacts.logits_ensemble is not None:
+                np.save(fold_dir / "logits_fused.npy", artifacts.logits_ensemble)
+            if artifacts.logits_ensemble_img is not None:
+                np.save(fold_dir / "logits_img.npy", artifacts.logits_ensemble_img)
+            if artifacts.logits_ensemble_md is not None:
+                np.save(fold_dir / "logits_md.npy", artifacts.logits_ensemble_md)
+            if artifacts.logits_classic is not None:
+                np.save(fold_dir / "logits_classic.npy", artifacts.logits_classic)
+            if artifacts.logits_classic_img is not None:
+                np.save(fold_dir / "logits_classic_img.npy", artifacts.logits_classic_img)
+            if artifacts.logits_classic_md is not None:
+                np.save(fold_dir / "logits_classic_md.npy", artifacts.logits_classic_md)
+            if artifacts.logits_ensemble_pereye is not None:
+                np.save(fold_dir / "logits_fused_pereye.npy", artifacts.logits_ensemble_pereye)
+            if artifacts.logits_ensemble_img_pereye is not None:
+                np.save(fold_dir / "logits_img_pereye.npy", artifacts.logits_ensemble_img_pereye)
+            if artifacts.logits_ensemble_md_pereye is not None:
+                np.save(fold_dir / "logits_md_pereye.npy", artifacts.logits_ensemble_md_pereye)
             if artifacts.probs_bilat is not None:
                 np.save(fold_dir / "probs_bilat.npy", artifacts.probs_bilat)
+            if artifacts.probs_fused is not None:
+                np.save(fold_dir / "probs_fused_head.npy", artifacts.probs_fused)
+            # y_true is shared across all heads for the same fold
+            if artifacts.y_true_bilat is not None and artifacts.y_true_ensemble is None:
+                np.save(fold_dir / "y_true.npy", artifacts.y_true_bilat)
+            # per-sample prediction CSVs
+            if artifacts.y_true_ensemble is not None:
+                _save_predictions_csv(
+                    fold_dir, mode, artifacts.y_true_ensemble,
+                    {"fused": artifacts.probs_ensemble,
+                     "img":   artifacts.probs_ensemble_img,
+                     "md":    artifacts.probs_ensemble_md},
+                )
+            if artifacts.y_true_ensemble_pereye is not None:
+                _save_predictions_csv(
+                    fold_dir, mode, artifacts.y_true_ensemble_pereye,
+                    {"fused": artifacts.probs_ensemble_pereye,
+                     "img":   artifacts.probs_ensemble_img_pereye,
+                     "md":    artifacts.probs_ensemble_md_pereye},
+                    suffix="_pereye",
+                )
+            if artifacts.y_true_classic is not None:
+                _save_predictions_csv(
+                    fold_dir, mode, artifacts.y_true_classic,
+                    {"fused": artifacts.probs_classic,
+                     "img":   artifacts.probs_classic_img,
+                     "md":    artifacts.probs_classic_md},
+                    suffix="_classic",
+                )
 
         fold_csv = tm_dir / "fold_results.csv"
         csv_fields = list(FoldResult.__dataclass_fields__.keys())
@@ -509,6 +669,7 @@ class V2HyperTower:
                 augment=args.augment, clinical_data=data,
                 num_classes=num_classes,
                 md_hidden_dim=args.md_hidden_dim, fusion_dim=args.fusion_dim,
+                bridge_mode=getattr(args, "bridge_mode", "fused"),
             ).to(device)
         if run_bilat:
             bilateral = BilateralHT(
@@ -525,6 +686,7 @@ class V2HyperTower:
         # ---- loaders ---------------------------------------------------
         use_balanced = bool(getattr(args, "balanced_sampling", False))
         train_single_loader = None
+        train_eval_loader   = None   # non-shuffled, no sampler — for per-epoch train logging
         train_bilat_loader  = None
         if run_single:
             single_sampler = build_balanced_sampler(eye_train) if use_balanced else None
@@ -536,6 +698,13 @@ class V2HyperTower:
                 sampler=single_sampler,
                 **loader_kw,
             )
+            train_eval_loader = make_loader(
+                eye_train, slots_eye,
+                image_transform=build_eval_transform(args.backbone),
+                image_preprocessor=image_preprocessor,
+                shuffle=False,
+                **loader_kw,
+            )
         if run_bilat:
             bilat_sampler = build_balanced_sampler(bilat_train) if use_balanced else None
             train_bilat_loader = make_loader(
@@ -593,18 +762,68 @@ class V2HyperTower:
             "main_epoch_single", "main_epoch_bilat",
             "single_active", "bilat_active",
             "single_train_loss", "single_train_acc",
+            # val — fused head (existing)
             "classic_val_auc",      "classic_val_acc",      "classic_val_n",
             "ensemble_val_auc",     "ensemble_val_acc",     "ensemble_val_n",
             "bilat_train_loss",     "bilat_train_acc",
             "bilat_val_auc",        "bilat_val_acc",        "bilat_val_n",
+            # val — img/md heads + fusion events
+            "classic_val_auc_img",  "classic_val_acc_img",
+            "classic_val_auc_md",   "classic_val_acc_md",
+            "classic_val_fe_corr",  "classic_val_fe_err",
+            "ensemble_val_auc_img", "ensemble_val_acc_img",
+            "ensemble_val_auc_md",  "ensemble_val_acc_md",
+            "ensemble_val_fe_corr", "ensemble_val_fe_err",
+            "bilat_val_auc_img",    "bilat_val_acc_img",
+            "bilat_val_auc_md",     "bilat_val_acc_md",
+            "bilat_val_fe_corr",    "bilat_val_fe_err",
+            # holdout — fused head (existing)
             "classic_holdout_auc",  "classic_holdout_acc",
             "ensemble_holdout_auc", "ensemble_holdout_acc",
             "bilat_holdout_auc",    "bilat_holdout_acc",
+            # holdout — img/md heads + fusion events
+            "classic_holdout_auc_img",  "classic_holdout_acc_img",
+            "classic_holdout_auc_md",   "classic_holdout_acc_md",
+            "classic_holdout_fe_corr",  "classic_holdout_fe_err",
+            "ensemble_holdout_auc_img", "ensemble_holdout_acc_img",
+            "ensemble_holdout_auc_md",  "ensemble_holdout_acc_md",
+            "ensemble_holdout_fe_corr", "ensemble_holdout_fe_err",
+            # train-set eval pass (eval mode, all 3 heads)
+            "train_auc_fused",  "train_acc_fused",
+            "train_auc_img",    "train_acc_img",
+            "train_auc_md",     "train_acc_md",
+            "train_fe_corr",    "train_fe_err",
+            "train_n",
             "is_best_single", "is_best_bilat",
             "is_best_holdout_single", "is_best_holdout_bilat",
         ]
+        # CM columns — named by num_classes so binary and multiclass both work
+        if num_classes == 2:
+            _cm_keys = ["tn", "fp", "fn", "tp"]
+        else:
+            _cm_keys = [f"cm_{i}_{j}" for i in range(num_classes) for j in range(num_classes)]
+        for _split in ("classic_val", "ensemble_val", "classic_holdout", "ensemble_holdout", "train"):
+            for _head in ("fused", "img", "md"):
+                for _k in _cm_keys:
+                    epoch_fields.append(f"{_split}_{_head}_{_k}")
         fold_logger = HypertowerLogger(run_dir=fold_dir)
 
+        # per-epoch accumulation for npy tensors
+        _epoch_train_pf:  list[np.ndarray] = []
+        _epoch_train_pi:  list[np.ndarray] = []
+        _epoch_train_pm:  list[np.ndarray] = []
+        _epoch_train_ids: list[np.ndarray] = []
+        _epoch_train_y:   list[np.ndarray] = []
+        # per-eye val accumulators (ensemble mode: OD and OS separate)
+        _epoch_val_pf_od: list[np.ndarray] = []
+        _epoch_val_pi_od: list[np.ndarray] = []
+        _epoch_val_pm_od: list[np.ndarray] = []
+        _epoch_val_pf_os: list[np.ndarray] = []
+        _epoch_val_pi_os: list[np.ndarray] = []
+        _epoch_val_pm_os: list[np.ndarray] = []
+        _epoch_val_y:     list[np.ndarray] = []
+        _epoch_val_ids:   list[np.ndarray] = []
+
         # ---- best-epoch trackers ---------------------------------------
         best_single_auc = -1.0
         best_bilat_auc  = -1.0
@@ -704,9 +923,16 @@ class V2HyperTower:
                 en_n = 0
                 en_acc_img = en_acc_md = en_auc_img = en_auc_md = nan
             elif run_single and tower_mode == "ensemble":
-                y_en, p_en, p_en_img, p_en_md = collect_probs_single_components(
-                    single, val_loader, device, aggregate_patient=True
+                (y_en,
+                 _p_en_f_od, _p_en_i_od, _p_en_m_od,
+                 _p_en_f_os, _p_en_i_os, _p_en_m_os,
+                 _en_pat_ids) = collect_probs_ensemble_pereye(
+                    single, val_loader, device, return_ids=True
                 )
+                # patient-level averages (used for metrics, same as before)
+                p_en     = 0.5 * (_p_en_f_od + _p_en_f_os)
+                p_en_img = 0.5 * (_p_en_i_od + _p_en_i_os)
+                p_en_md  = 0.5 * (_p_en_m_od + _p_en_m_os)
                 en_acc, en_auc, en_n = _score_arrays(y_en, p_en, num_classes)
                 en_acc_img = float((p_en_img.argmax(1) == y_en).mean()) if y_en.size else nan
                 en_acc_md  = float((p_en_md.argmax(1) == y_en).mean()) if y_en.size else nan
@@ -742,23 +968,47 @@ class V2HyperTower:
                 bi_acc_img = bi_acc_md = bi_auc_img = bi_auc_md = nan
 
             # --- holdout evaluation ------------------------------------
+            # defaults (overwritten below when holdout_loader is not None)
+            _z2 = np.zeros((0, num_classes), dtype=np.float32)
+            _e2 = np.array([], dtype=np.int64)
+            y_cl_h = y_en_h = _e2
+            p_cl_h = p_cl_h_img = p_cl_h_md = _z2
+            p_en_h = p_en_h_img = p_en_h_md = _z2
+
             if holdout_loader is not None:
                 if run_single and tower_mode == "single":
-                    y_cl_h, p_cl_h, _, _ = collect_probs_single_components(
+                    y_cl_h, p_cl_h, p_cl_h_img, p_cl_h_md = collect_probs_single_components(
                         single, holdout_loader, device, aggregate_patient=False
                     )
                     _, cl_auc_h, _ = _score_arrays(y_cl_h, p_cl_h, num_classes)
                     cl_acc_h = float((p_cl_h.argmax(1) == y_cl_h).mean()) if y_cl_h.size else nan
+                    _, cl_auc_h_img, _ = _score_arrays(y_cl_h, p_cl_h_img, num_classes)
+                    cl_acc_h_img = float((p_cl_h_img.argmax(1) == y_cl_h).mean()) if y_cl_h.size else nan
+                    _, cl_auc_h_md, _  = _score_arrays(y_cl_h, p_cl_h_md,  num_classes)
+                    cl_acc_h_md  = float((p_cl_h_md.argmax(1) == y_cl_h).mean()) if y_cl_h.size else nan
+                    cl_fe_h_corr, cl_fe_h_err = _fusion_events(y_cl_h, p_cl_h, p_cl_h_img, p_cl_h_md)
                     en_auc_h = en_acc_h = nan
+                    en_auc_h_img = en_acc_h_img = en_auc_h_md = en_acc_h_md = nan
+                    en_fe_h_corr = en_fe_h_err = 0
                 elif run_single and tower_mode == "ensemble":
-                    y_en_h, p_en_h, _, _ = collect_probs_single_components(
+                    y_en_h, p_en_h, p_en_h_img, p_en_h_md = collect_probs_single_components(
                         single, holdout_loader, device, aggregate_patient=True
                     )
                     _, en_auc_h, _ = _score_arrays(y_en_h, p_en_h, num_classes)
                     en_acc_h = float((p_en_h.argmax(1) == y_en_h).mean()) if y_en_h.size else nan
+                    _, en_auc_h_img, _ = _score_arrays(y_en_h, p_en_h_img, num_classes)
+                    en_acc_h_img = float((p_en_h_img.argmax(1) == y_en_h).mean()) if y_en_h.size else nan
+                    _, en_auc_h_md, _  = _score_arrays(y_en_h, p_en_h_md,  num_classes)
+                    en_acc_h_md  = float((p_en_h_md.argmax(1) == y_en_h).mean()) if y_en_h.size else nan
+                    en_fe_h_corr, en_fe_h_err = _fusion_events(y_en_h, p_en_h, p_en_h_img, p_en_h_md)
                     cl_auc_h = cl_acc_h = nan
+                    cl_auc_h_img = cl_acc_h_img = cl_auc_h_md = cl_acc_h_md = nan
+                    cl_fe_h_corr = cl_fe_h_err = 0
                 else:
                     cl_auc_h = cl_acc_h = en_auc_h = en_acc_h = nan
+                    cl_auc_h_img = cl_acc_h_img = cl_auc_h_md = cl_acc_h_md = nan
+                    en_auc_h_img = en_acc_h_img = en_auc_h_md = en_acc_h_md = nan
+                    cl_fe_h_corr = cl_fe_h_err = en_fe_h_corr = en_fe_h_err = 0
                 if run_bilat:
                     y_bi_h, p_bi_h, _, _ = collect_probs_bilateral_components(
                         bilateral, holdout_loader, device
@@ -769,6 +1019,59 @@ class V2HyperTower:
                     bi_auc_h = bi_acc_h = nan
             else:
                 cl_auc_h = cl_acc_h = en_auc_h = en_acc_h = bi_auc_h = bi_acc_h = nan
+                cl_auc_h_img = cl_acc_h_img = cl_auc_h_md = cl_acc_h_md = nan
+                en_auc_h_img = en_acc_h_img = en_auc_h_md = en_acc_h_md = nan
+                cl_fe_h_corr = cl_fe_h_err = en_fe_h_corr = en_fe_h_err = 0
+
+            # --- fusion-event helpers for val sets ----------------------
+            cl_fe_corr, cl_fe_err = _fusion_events(y_cl, p_cl, p_cl_img, p_cl_md) if y_cl.size else (0, 0)
+            en_fe_corr, en_fe_err = _fusion_events(y_en, p_en, p_en_img, p_en_md) if y_en.size else (0, 0)
+            bi_fe_corr, bi_fe_err = (0, 0)  # bilateral components not separated the same way
+
+            # --- train eval pass (eval mode, all 3 heads) ----------------
+            tr_auc_f = tr_acc_f = tr_auc_i = tr_acc_i = tr_auc_m = tr_acc_m = nan
+            tr_fe_corr = tr_fe_err = tr_n = 0
+            y_tr = np.array([], dtype=np.int64)
+            p_tr_f = p_tr_i = p_tr_m = np.zeros((0, num_classes), dtype=np.float32)
+            if run_single and train_eval_loader is not None:
+                y_tr, p_tr_f, p_tr_i, p_tr_m, tr_ids = collect_probs_eye_level(
+                    single, train_eval_loader, device, return_ids=True
+                )
+                if y_tr.size:
+                    _, tr_auc_f, _ = _score_arrays(y_tr, p_tr_f, num_classes)
+                    tr_acc_f = float((p_tr_f.argmax(1) == y_tr).mean())
+                    _, tr_auc_i, _ = _score_arrays(y_tr, p_tr_i, num_classes)
+                    tr_acc_i = float((p_tr_i.argmax(1) == y_tr).mean())
+                    _, tr_auc_m, _ = _score_arrays(y_tr, p_tr_m, num_classes)
+                    tr_acc_m = float((p_tr_m.argmax(1) == y_tr).mean())
+                    tr_fe_corr, tr_fe_err = _fusion_events(y_tr, p_tr_f, p_tr_i, p_tr_m)
+                    tr_n = int(y_tr.size)
+                    # accumulate for npy tensors
+                    _epoch_train_pf.append(p_tr_f)
+                    _epoch_train_pi.append(p_tr_i)
+                    _epoch_train_pm.append(p_tr_m)
+                    _epoch_train_ids.append(tr_ids)
+                    _epoch_train_y.append(y_tr)
+
+            # accumulate val for npy tensors
+            if run_single and tower_mode == "ensemble" and y_en.size:
+                _epoch_val_pf_od.append(_p_en_f_od)
+                _epoch_val_pi_od.append(_p_en_i_od)
+                _epoch_val_pm_od.append(_p_en_m_od)
+                _epoch_val_pf_os.append(_p_en_f_os)
+                _epoch_val_pi_os.append(_p_en_i_os)
+                _epoch_val_pm_os.append(_p_en_m_os)
+                _epoch_val_y.append(y_en)
+                _epoch_val_ids.append(_en_pat_ids)
+            elif run_single and tower_mode == "single" and y_cl.size:
+                # single mode: no per-eye split, reuse same array for both slots
+                _epoch_val_pf_od.append(p_cl)
+                _epoch_val_pi_od.append(p_cl_img)
+                _epoch_val_pm_od.append(p_cl_md)
+                _epoch_val_pf_os.append(p_cl)
+                _epoch_val_pi_os.append(p_cl_img)
+                _epoch_val_pm_os.append(p_cl_md)
+                _epoch_val_y.append(y_cl)
 
             # Best-epoch checks (restricted to main phase).
             target_single_auc         = cl_auc if tower_mode == "single" else en_auc
@@ -830,23 +1133,65 @@ class V2HyperTower:
                 best_epoch_holdout_bilat = epoch + 1
                 best_holdout_bilat_state = copy.deepcopy(bilateral.state_dict())
 
+            # --- confusion matrix cells per split × head -------------------
+            def _prefixed_cm(prefix: str, y: np.ndarray, pf: np.ndarray,
+                              pi: np.ndarray, pm: np.ndarray) -> dict:
+                out: dict = {}
+                for head, p in (("fused", pf), ("img", pi), ("md", pm)):
+                    for k, v in _cm_cells(y, p, num_classes).items():
+                        out[f"{prefix}_{head}_{k}"] = v
+                return out
+
+            cm_row: dict = {}
+            cm_row.update(_prefixed_cm("classic_val",        y_cl,   p_cl,     p_cl_img,   p_cl_md))
+            cm_row.update(_prefixed_cm("ensemble_val",       y_en,   p_en,     p_en_img,   p_en_md))
+            cm_row.update(_prefixed_cm("classic_holdout",    y_cl_h, p_cl_h,   p_cl_h_img, p_cl_h_md))
+            cm_row.update(_prefixed_cm("ensemble_holdout",   y_en_h, p_en_h,   p_en_h_img, p_en_h_md))
+            cm_row.update(_prefixed_cm("train",              y_tr,   p_tr_f,   p_tr_i,     p_tr_m))
+
             fold_logger.write_epoch_row({
                 "fold": fold, "epoch": epoch + 1,
                 "phase_single": phase_single, "phase_bilat": phase_bilat,
                 "main_epoch_single": main_epoch_single, "main_epoch_bilat": main_epoch_bilat,
                 "single_active": int(single_active), "bilat_active": int(bilat_active),
                 "single_train_loss": _f(sl_loss), "single_train_acc": _f(sl_acc),
+                # val — fused
                 "classic_val_auc":  _f(cl_auc), "classic_val_acc":  _f(cl_acc), "classic_val_n":  cl_n,
                 "ensemble_val_auc": _f(en_auc), "ensemble_val_acc": _f(en_acc), "ensemble_val_n": en_n,
                 "bilat_train_loss": _f(bl_loss), "bilat_train_acc": _f(bl_acc),
                 "bilat_val_auc":    _f(bi_auc), "bilat_val_acc":    _f(bi_acc), "bilat_val_n":    bi_n,
+                # val — img/md + fusion events
+                "classic_val_auc_img":  _f(cl_auc_img),  "classic_val_acc_img":  _f(cl_acc_img),
+                "classic_val_auc_md":   _f(cl_auc_md),   "classic_val_acc_md":   _f(cl_acc_md),
+                "classic_val_fe_corr":  cl_fe_corr,       "classic_val_fe_err":   cl_fe_err,
+                "ensemble_val_auc_img": _f(en_auc_img),  "ensemble_val_acc_img": _f(en_acc_img),
+                "ensemble_val_auc_md":  _f(en_auc_md),   "ensemble_val_acc_md":  _f(en_acc_md),
+                "ensemble_val_fe_corr": en_fe_corr,       "ensemble_val_fe_err":  en_fe_err,
+                "bilat_val_auc_img":    _f(bi_auc_img),  "bilat_val_acc_img":    _f(bi_acc_img),
+                "bilat_val_auc_md":     _f(bi_auc_md),   "bilat_val_acc_md":     _f(bi_acc_md),
+                "bilat_val_fe_corr":    bi_fe_corr,       "bilat_val_fe_err":     bi_fe_err,
+                # holdout — fused
                 "classic_holdout_auc":  _f(cl_auc_h), "classic_holdout_acc":  _f(cl_acc_h),
                 "ensemble_holdout_auc": _f(en_auc_h), "ensemble_holdout_acc": _f(en_acc_h),
                 "bilat_holdout_auc":    _f(bi_auc_h), "bilat_holdout_acc":    _f(bi_acc_h),
+                # holdout — img/md + fusion events
+                "classic_holdout_auc_img":  _f(cl_auc_h_img),  "classic_holdout_acc_img":  _f(cl_acc_h_img),
+                "classic_holdout_auc_md":   _f(cl_auc_h_md),   "classic_holdout_acc_md":   _f(cl_acc_h_md),
+                "classic_holdout_fe_corr":  cl_fe_h_corr,       "classic_holdout_fe_err":   cl_fe_h_err,
+                "ensemble_holdout_auc_img": _f(en_auc_h_img),  "ensemble_holdout_acc_img": _f(en_acc_h_img),
+                "ensemble_holdout_auc_md":  _f(en_auc_h_md),   "ensemble_holdout_acc_md":  _f(en_acc_h_md),
+                "ensemble_holdout_fe_corr": en_fe_h_corr,       "ensemble_holdout_fe_err":  en_fe_h_err,
+                # train eval pass
+                "train_auc_fused": _f(tr_auc_f), "train_acc_fused": _f(tr_acc_f),
+                "train_auc_img":   _f(tr_auc_i), "train_acc_img":   _f(tr_acc_i),
+                "train_auc_md":    _f(tr_auc_m), "train_acc_md":    _f(tr_acc_m),
+                "train_fe_corr":   tr_fe_corr,    "train_fe_err":    tr_fe_err,
+                "train_n":         tr_n,
                 "is_best_single":         int(is_best_single),
                 "is_best_bilat":          int(is_best_bilat),
                 "is_best_holdout_single": int(is_best_holdout_single),
                 "is_best_holdout_bilat":  int(is_best_holdout_bilat),
+                **cm_row,
             }, optional_cols=epoch_fields)
 
             if args.log_every > 0 and (epoch + 1) % args.log_every == 0:
@@ -898,6 +1243,27 @@ class V2HyperTower:
             if best_holdout_bilat_state is not None:
                 torch.save(best_holdout_bilat_state, fold_dir / "best_holdout_bilateral.pt")
 
+        # ---- Save per-epoch per-patient npy tensors -------------------------
+        if _epoch_train_pf:
+            # Use the order from the first epoch (consistent since loader is non-shuffled)
+            ids_ref = _epoch_train_ids[0]
+            y_ref   = _epoch_train_y[0]
+            np.save(fold_dir / "train_patient_ids.npy", ids_ref)
+            np.save(fold_dir / "train_y_true.npy",      y_ref)
+            np.save(fold_dir / "train_probs_fused.npy", np.stack(_epoch_train_pf))  # (n_ep, n_pts, n_cls)
+            np.save(fold_dir / "train_probs_img.npy",   np.stack(_epoch_train_pi))
+            np.save(fold_dir / "train_probs_md.npy",    np.stack(_epoch_train_pm))
+        if _epoch_val_pf_od:
+            np.save(fold_dir / "val_y_true_epochs.npy",         np.stack(_epoch_val_y))     # (n_ep, N)
+            np.save(fold_dir / "val_probs_fused_od_epochs.npy", np.stack(_epoch_val_pf_od)) # (n_ep, N, C)
+            np.save(fold_dir / "val_probs_img_od_epochs.npy",   np.stack(_epoch_val_pi_od))
+            np.save(fold_dir / "val_probs_md_od_epochs.npy",    np.stack(_epoch_val_pm_od))
+            np.save(fold_dir / "val_probs_fused_os_epochs.npy", np.stack(_epoch_val_pf_os))
+            np.save(fold_dir / "val_probs_img_os_epochs.npy",   np.stack(_epoch_val_pi_os))
+            np.save(fold_dir / "val_probs_md_os_epochs.npy",    np.stack(_epoch_val_pm_os))
+            if _epoch_val_ids:  # only ensemble mode populates this
+                np.save(fold_dir / "val_patient_ids.npy", _epoch_val_ids[0])
+
         # ---- Phase 2: fused head training (ensemble + --fused-head only) ----
         best_fused_auc        = -1.0
         best_fused_state: Optional[dict] = None
@@ -999,15 +1365,30 @@ class V2HyperTower:
             single.load_state_dict(best_single_state)
         if run_bilat and best_bilat_state is not None:
             bilateral.load_state_dict(best_bilat_state)
+        y_en_pe_best = p_en_pe_best = p_en_pe_best_img = p_en_pe_best_md = None
+        l_en_best = l_en_best_img = l_en_best_md = None
+        l_cl_best = l_cl_best_img = l_cl_best_md = None
+        l_en_pe_best = l_en_pe_best_img = l_en_pe_best_md = None
         if run_single and tower_mode == "single":
-            y_cl_best, p_cl_best = collect_probs_classic(single, val_loader, device)
-            y_en_best = p_en_best = None
+            y_cl_best, p_cl_best, p_cl_best_img, p_cl_best_md, \
+                l_cl_best, l_cl_best_img, l_cl_best_md = collect_probs_single_components(
+                single, val_loader, device, aggregate_patient=False, return_logits=True
+            )
+            y_en_best = p_en_best = p_en_best_img = p_en_best_md = None
         elif run_single and tower_mode == "ensemble":
-            y_en_best, p_en_best = collect_probs_ensemble(single, val_loader, device)
-            y_cl_best = p_cl_best = None
+            y_en_best, p_en_best, p_en_best_img, p_en_best_md, \
+                l_en_best, l_en_best_img, l_en_best_md = collect_probs_single_components(
+                single, val_loader, device, aggregate_patient=True, return_logits=True
+            )
+            y_en_pe_best, p_en_pe_best, p_en_pe_best_img, p_en_pe_best_md, \
+                l_en_pe_best, l_en_pe_best_img, l_en_pe_best_md = collect_probs_single_components(
+                single, val_loader, device, aggregate_patient=False, return_logits=True
+            )
+            y_cl_best = p_cl_best = p_cl_best_img = p_cl_best_md = None
         else:
             y_cl_best = y_en_best = None
-            p_cl_best = p_en_best = None
+            p_cl_best = p_en_best = p_en_best_img = p_en_best_md = None
+            p_cl_best_img = p_cl_best_md = None
         if run_bilat:
             y_bi_best, p_bi_best = collect_probs_bilateral(bilateral, val_loader, device)
         else:
@@ -1077,6 +1458,23 @@ class V2HyperTower:
             y_true_ensemble=y_en_best, probs_ensemble=p_en_best,
             y_true_bilat=y_bi_best,    probs_bilat=p_bi_best,
             y_true_fused=y_fu_best,    probs_fused=p_fu_best,
+            probs_ensemble_img=p_en_best_img,
+            probs_ensemble_md=p_en_best_md,
+            probs_classic_img=p_cl_best_img,
+            probs_classic_md=p_cl_best_md,
+            y_true_ensemble_pereye=y_en_pe_best,
+            probs_ensemble_pereye=p_en_pe_best,
+            probs_ensemble_img_pereye=p_en_pe_best_img,
+            probs_ensemble_md_pereye=p_en_pe_best_md,
+            logits_ensemble=l_en_best,
+            logits_ensemble_img=l_en_best_img,
+            logits_ensemble_md=l_en_best_md,
+            logits_classic=l_cl_best,
+            logits_classic_img=l_cl_best_img,
+            logits_classic_md=l_cl_best_md,
+            logits_ensemble_pereye=l_en_pe_best,
+            logits_ensemble_img_pereye=l_en_pe_best_img,
+            logits_ensemble_md_pereye=l_en_pe_best_md,
         )
 
     # ------------------------------------------------------------------
diff --git a/results/basic_analysis/ensemble_hypertower_binary.png b/results/basic_analysis/ensemble_hypertower_binary.png
new file mode 100644
index 0000000..284a569
Binary files /dev/null and b/results/basic_analysis/ensemble_hypertower_binary.png differ
diff --git a/results/basic_analysis/ensemble_hypertower_multiclass.png b/results/basic_analysis/ensemble_hypertower_multiclass.png
new file mode 100644
index 0000000..02ba1ef
Binary files /dev/null and b/results/basic_analysis/ensemble_hypertower_multiclass.png differ
diff --git a/results/basic_analysis/fused_head_binary.png b/results/basic_analysis/fused_head_binary.png
new file mode 100644
index 0000000..678ceb2
Binary files /dev/null and b/results/basic_analysis/fused_head_binary.png differ
diff --git a/results/basic_analysis/fused_head_multiclass.png b/results/basic_analysis/fused_head_multiclass.png
new file mode 100644
index 0000000..be073a9
Binary files /dev/null and b/results/basic_analysis/fused_head_multiclass.png differ
diff --git a/results/basic_analysis/gradcam_summary_grid.png b/results/basic_analysis/gradcam_summary_grid.png
new file mode 100644
index 0000000..7ea6dfc
Binary files /dev/null and b/results/basic_analysis/gradcam_summary_grid.png differ
diff --git a/results/basic_analysis/image_only_binary.png b/results/basic_analysis/image_only_binary.png
new file mode 100755
index 0000000..cb274d8
Binary files /dev/null and b/results/basic_analysis/image_only_binary.png differ
diff --git a/results/basic_analysis/image_only_binary_correct.png b/results/basic_analysis/image_only_binary_correct.png
new file mode 100755
index 0000000..72de4e9
Binary files /dev/null and b/results/basic_analysis/image_only_binary_correct.png differ
diff --git a/results/basic_analysis/image_only_fused.png b/results/basic_analysis/image_only_fused.png
new file mode 100755
index 0000000..6e7c6da
Binary files /dev/null and b/results/basic_analysis/image_only_fused.png differ
diff --git a/results/basic_analysis/md_only_binary.png b/results/basic_analysis/md_only_binary.png
new file mode 100644
index 0000000..813ebee
Binary files /dev/null and b/results/basic_analysis/md_only_binary.png differ
diff --git a/results/basic_analysis/md_only_multiclass.png b/results/basic_analysis/md_only_multiclass.png
new file mode 100644
index 0000000..7cf3a0e
Binary files /dev/null and b/results/basic_analysis/md_only_multiclass.png differ
diff --git a/results/basic_analysis/md_permutation_importance copy.png b/results/basic_analysis/md_permutation_importance copy.png
new file mode 100644
index 0000000..1271631
Binary files /dev/null and b/results/basic_analysis/md_permutation_importance copy.png differ
diff --git a/results/basic_analysis/md_permutation_importance.png b/results/basic_analysis/md_permutation_importance.png
new file mode 100644
index 0000000..8a7f8f5
Binary files /dev/null and b/results/basic_analysis/md_permutation_importance.png differ
diff --git a/results/basic_analysis/papila_knn_roc_mean.png b/results/basic_analysis/papila_knn_roc_mean.png
new file mode 100755
index 0000000..d945d9d
Binary files /dev/null and b/results/basic_analysis/papila_knn_roc_mean.png differ
diff --git a/results/basic_analysis/papila_logistic_regression_roc_mean.png b/results/basic_analysis/papila_logistic_regression_roc_mean.png
new file mode 100755
index 0000000..38104b9
Binary files /dev/null and b/results/basic_analysis/papila_logistic_regression_roc_mean.png differ
diff --git a/results/basic_analysis/papila_random_forest_feature_importance_top.png b/results/basic_analysis/papila_random_forest_feature_importance_top.png
new file mode 100755
index 0000000..23a8ac9
Binary files /dev/null and b/results/basic_analysis/papila_random_forest_feature_importance_top.png differ
diff --git a/results/basic_analysis/papila_random_forest_roc_mean.png b/results/basic_analysis/papila_random_forest_roc_mean.png
new file mode 100755
index 0000000..205e224
Binary files /dev/null and b/results/basic_analysis/papila_random_forest_roc_mean.png differ
diff --git a/results/basic_analysis/papila_svm_roc_mean.png b/results/basic_analysis/papila_svm_roc_mean.png
new file mode 100755
index 0000000..ff305dd
Binary files /dev/null and b/results/basic_analysis/papila_svm_roc_mean.png differ
diff --git a/results/basic_analysis/patient_213_OD_OS.png b/results/basic_analysis/patient_213_OD_OS.png
new file mode 100644
index 0000000..e620f87
Binary files /dev/null and b/results/basic_analysis/patient_213_OD_OS.png differ
diff --git a/results/basic_analysis/patient_275_OD_OS.png b/results/basic_analysis/patient_275_OD_OS.png
new file mode 100644
index 0000000..4605494
Binary files /dev/null and b/results/basic_analysis/patient_275_OD_OS.png differ
diff --git a/results/basic_analysis/single_hypertower_binary.png b/results/basic_analysis/single_hypertower_binary.png
new file mode 100644
index 0000000..7dd73d2
Binary files /dev/null and b/results/basic_analysis/single_hypertower_binary.png differ
diff --git a/results/basic_analysis/single_hypertower_binary_correct.png b/results/basic_analysis/single_hypertower_binary_correct.png
new file mode 100644
index 0000000..1781871
Binary files /dev/null and b/results/basic_analysis/single_hypertower_binary_correct.png differ
diff --git a/results/basic_analysis/single_hypertower_multiclass.png b/results/basic_analysis/single_hypertower_multiclass.png
new file mode 100755
index 0000000..60d7f18
Binary files /dev/null and b/results/basic_analysis/single_hypertower_multiclass.png differ
diff --git a/scripts/main/v2/run_binary_v2.2.sh b/scripts/main/v2/run_binary_v2.2.sh
new file mode 100755
index 0000000..2b18228
--- /dev/null
+++ b/scripts/main/v2/run_binary_v2.2.sh
@@ -0,0 +1,74 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+# Binary runs v2.2 (6 total):
+#   UNet crop:  single | ensemble | fused head
+#   GT crop:    single | ensemble | fused head
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
+cd "$ROOT_DIR"
+
+MANIFEST="manifest.csv"
+UNET_WEIGHTS="models/v2/refuge/segmentation/per_image/best.pt"
+
+COMMON=(
+  --epochs 40
+  --n-splits 5
+  --batch-size 8
+  --backbone refugelike
+  --eval-mode binary
+  --img-crop-manifest "$MANIFEST"
+)
+
+UNET_CROP=(
+  --img-crop-weights "$UNET_WEIGHTS"
+)
+
+GT_CROP=(
+  --img-crop-gt
+)
+
+# ── UNet crop ────────────────────────────────────────────────────────────────
+
+echo "[1/6] UNet crop — binary, single..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode single \
+  --run-name v2.2_single_binary_unet_40ep_5fold
+
+echo "[2/6] UNet crop — binary, ensemble..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode ensemble \
+  --run-name v2.2_ensemble_binary_unet_40ep_5fold
+
+echo "[3/6] UNet crop — binary, ensemble + fused head..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode ensemble \
+  --fused-head --fusion-epochs 10 \
+  --run-name v2.2_fused_binary_unet_40ep_5fold
+
+# ── GT crop ──────────────────────────────────────────────────────────────────
+
+echo "[4/6] GT crop — binary, single..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode single \
+  --run-name v2.2_single_binary_gt_40ep_5fold
+
+echo "[5/6] GT crop — binary, ensemble..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode ensemble \
+  --run-name v2.2_ensemble_binary_gt_40ep_5fold
+
+echo "[6/6] GT crop — binary, ensemble + fused head..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode ensemble \
+  --img-crop-gt \
+  --fused-head --fusion-epochs 10 \
+  --run-name v2.2_fused_binary_gt_40ep_5fold
+
+echo "Binary v2.2 runs complete."
diff --git a/scripts/main/v2/run_fused_head_overnight.sh b/scripts/main/v2/run_fused_head_overnight.sh
deleted file mode 100755
index dda1dc0..0000000
--- a/scripts/main/v2/run_fused_head_overnight.sh
+++ /dev/null
@@ -1,41 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Fused-head ensemble runs (2 total) — UNet ROI crop:
-#   binary     × ensemble + fused head
-#   multiclass × ensemble + fused head
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-weights models/v2/refuge/segmentation/per_image_refuge_build/best.pt
-  --img-crop-normalize per_image
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-  --tower-mode ensemble
-  --fused-head
-  --fusion-epochs 10
-)
-
-echo "[1/2] UNet ROI — binary, ensemble + fused head..."
-python3 scripts/main/v2/run_multifold_v2.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-mode binary \
-  --run-name v2_ensemble_fused_binary_unet_40ep_5fold_v1
-
-echo "[2/2] UNet ROI — multiclass, ensemble + fused head..."
-python3 scripts/main/v2/run_multifold_v2.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-mode multiclass \
-  --run-name v2_ensemble_fused_multiclass_unet_40ep_5fold_v1
-
-echo "Fused-head runs complete."
diff --git a/scripts/main/v2/run_md_mlp.py b/scripts/main/v2/run_md_mlp.py
new file mode 100644
index 0000000..4e0b123
--- /dev/null
+++ b/scripts/main/v2/run_md_mlp.py
@@ -0,0 +1,461 @@
+#!/usr/bin/env python3
+"""
+V2-parity metadata-only runner.
+
+Goal:
+- Match V2HyperTower single-model metadata-only behavior as closely as possible.
+- Avoid image tower/image IO overhead in forward/training.
+
+How parity is achieved:
+- Uses PatientFirstSplitManager (same split policy).
+- Uses PAPILA profile builders + V2 filters:
+  - eye_train      = filter_eye_samples(...)
+  - bilat_val/test = filter_bilateral_samples(...)
+- Uses V2 training/eval helpers directly:
+  - train_single_epoch(...)
+  - collect_probs_single_components(...)
+- Uses bridge_mode="metadata_only".
+
+Implementation detail:
+- Batch dictionaries still include image slots to satisfy shared V2 helpers,
+  but these are tiny dummy tensors and are never consumed in metadata-only mode.
+
+Outputs:
+  analysis_data/{run_name}/{eval_mode}/{tower_mode}/fold{N}/
+    y_true.npy
+    probs_classic.npy or probs_ensemble.npy
+    y_true_holdout.npy
+    probs_classic_holdout.npy or probs_ensemble_holdout.npy
+"""
+
+from __future__ import annotations
+
+import argparse
+import copy
+import json
+import sys
+import time
+from pathlib import Path
+from types import SimpleNamespace
+
+# ensure repo root is on sys.path when run directly
+_REPO_ROOT = Path(__file__).resolve().parents[3]
+if str(_REPO_ROOT) not in sys.path:
+    sys.path.insert(0, str(_REPO_ROOT))
+
+import numpy as np
+import torch
+from torch import nn
+from torch.utils.data import DataLoader, Dataset
+
+from classes.v2.bridges import Bridge
+from classes.v2.loader_factory import (
+    build_balanced_sampler,
+    filter_bilateral_samples,
+    filter_eye_samples,
+)
+from classes.v2.metrics import _score_arrays
+from classes.v2.models import collect_probs_single_components, train_single_epoch
+from classes.v2.papila_builders import build_papila_data
+from classes.v2.profiles import build_papila_profile
+from classes.v2.split_manager import PatientFirstSplitManager
+from classes.v2.towers import MDTower
+from classes.v2.utils import choose_device, seed_everything
+
+
+class MetadataOnlySingleHT(nn.Module):
+    """SingleEyeHT-compatible shell without real image tower usage."""
+
+    def __init__(
+        self,
+        *,
+        clinical_data,
+        num_classes: int,
+        md_hidden_dim: int,
+        fusion_dim: int,
+        dropout: float,
+        use_se: bool,
+        se_reduction: int,
+        se_pre_norm: bool,
+    ):
+        super().__init__()
+        # Placeholder module to satisfy phase toggling logic.
+        self.img_tower = nn.Identity()
+        self.md_tower = MDTower(
+            clinical_data=clinical_data,
+            hidden_dim=md_hidden_dim,
+            dropout=dropout,
+            use_se=use_se,
+            se_reduction=se_reduction,
+            se_pre_norm=se_pre_norm,
+        )
+        # img_dim is irrelevant in metadata_only mode, but Bridge defines img head params.
+        self.bridge = Bridge(
+            img_dim=1,
+            meta_dim=self.md_tower.out_dim,
+            num_classes=num_classes,
+            fusion_dim=fusion_dim,
+            mode="metadata_only",
+            use_se=False,
+            se_reduction=16,
+            se_pre_norm=True,
+        )
+
+
+class EyeMetaDataset(Dataset):
+    """Eye-level dataset for V2 train_single_epoch input contract."""
+
+    def __init__(self, samples: list[dict]):
+        self.samples = samples
+
+    def __len__(self) -> int:
+        return len(self.samples)
+
+    def __getitem__(self, idx: int) -> dict[str, torch.Tensor]:
+        s = self.samples[idx]
+        return {
+            "image_1": torch.zeros(1, dtype=torch.float32),
+            "matrix_1": torch.as_tensor(s["matrix_1"], dtype=torch.float32),
+            "label_1": torch.tensor(int(s["label_1"]), dtype=torch.long),
+        }
+
+
+class BilatMetaDataset(Dataset):
+    """Patient-level bilateral dataset for collect_probs_single_components."""
+
+    def __init__(self, samples: list[dict]):
+        self.samples = samples
+
+    def __len__(self) -> int:
+        return len(self.samples)
+
+    def __getitem__(self, idx: int) -> dict[str, torch.Tensor]:
+        s = self.samples[idx]
+        return {
+            "image_1": torch.zeros(1, dtype=torch.float32),
+            "image_2": torch.zeros(1, dtype=torch.float32),
+            "matrix_1": torch.as_tensor(s["matrix_1"], dtype=torch.float32),
+            "matrix_2": torch.as_tensor(s["matrix_2"], dtype=torch.float32),
+            "label_1": torch.tensor(int(s["label_1"]), dtype=torch.long),
+        }
+
+
+def _phase_for_epoch(epoch_idx: int, warm_tower: int, warm_fused: int, main_epochs: int) -> tuple[str, int]:
+    if epoch_idx < warm_tower:
+        return "tower_warmup", 0
+    if epoch_idx < (warm_tower + warm_fused):
+        return "fused_warmup", 0
+    if epoch_idx < (warm_tower + warm_fused + main_epochs):
+        main_ep = epoch_idx - warm_tower - warm_fused + 1
+        return "main", main_ep
+    return "done", main_epochs
+
+
+def _evaluate_single(
+    model: nn.Module,
+    loader: DataLoader,
+    device: torch.device,
+    num_classes: int,
+    aggregate_patient: bool,
+) -> tuple[np.ndarray, np.ndarray, float, float]:
+    y, p_fused, _, p_md = collect_probs_single_components(
+        model, loader, device, aggregate_patient=aggregate_patient
+    )
+    # In metadata_only mode p_fused == p_md; keep md explicitly for clarity.
+    probs = p_md if p_md.size else p_fused
+    acc, auc, _ = _score_arrays(y, probs, num_classes)
+    return y, probs, float(auc), float(acc)
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description="V2-parity metadata-only runner.")
+    ap.add_argument("--eval-mode", required=True, choices=["binary", "multiclass"])
+    ap.add_argument("--tower-mode", default="single", choices=["single", "ensemble"])
+    ap.add_argument("--run-name", required=True)
+
+    ap.add_argument("--epochs", type=int, default=40, help="Main-phase epochs.")
+    ap.add_argument("--warmup-tower-epochs", type=int, default=None)
+    ap.add_argument("--warmup-fused-epochs", type=int, default=None)
+    ap.add_argument("--batch-size", type=int, default=8)
+    ap.add_argument("--lr", type=float, default=1e-4)
+    ap.add_argument("--weight-decay", type=float, default=0.0)
+    ap.add_argument("--bcd-prob", type=float, default=0.5)
+
+    ap.add_argument("--md-hidden-dim", type=int, default=128)
+    ap.add_argument("--fusion-dim", type=int, default=256)
+    ap.add_argument("--dropout", type=float, default=0.1)
+    ap.add_argument("--use-se", action="store_true")
+    ap.add_argument("--se-reduction", type=int, default=16)
+    ap.add_argument("--se-pre-norm", action="store_true")
+
+    ap.add_argument("--n-splits", type=int, default=5)
+    ap.add_argument("--holdout-per-class", type=int, default=5)
+    ap.add_argument("--holdout-seed", type=int, default=123)
+    ap.add_argument("--fold-seed", type=int, default=42)
+    ap.add_argument("--seed", type=int, default=1234)
+    ap.add_argument("--balanced-sampling", action=argparse.BooleanOptionalAction, default=False)
+
+    ap.add_argument("--analysis-dir", default="analysis_data")
+    ap.add_argument("--image-dir", default="Papila/FundusImages")
+    ap.add_argument("--clinical-dir", default="Papila/ClinicalData")
+    ap.add_argument("--label-col", default="Diagnosis")
+    ap.add_argument("--patient-col", default="Patient ID")
+    ap.add_argument("--cat-cols", nargs="*", default=["Gender", "Phakic/Pseudophakic"])
+
+    args = ap.parse_args()
+
+    seed_everything(args.seed)
+    device = choose_device(None)
+    print(f"Device: {device}", flush=True)
+
+    print("Loading PAPILA data...", flush=True)
+    data = build_papila_data(
+        image_dir=args.image_dir,
+        clinical_dir=args.clinical_dir,
+        label_col=args.label_col,
+        cat_cols=list(args.cat_cols),
+        n_splits=args.n_splits,
+        random_seed=args.fold_seed,
+    )
+    print(f"Loaded: {len(data.df)} rows  feature_dim={data.feature_dim}", flush=True)
+
+    num_classes = 2 if args.eval_mode == "binary" else 3
+    df_mode = data.df.copy()
+    if args.eval_mode == "binary":
+        df_mode = df_mode[df_mode[args.label_col].isin([0, 1])].reset_index(drop=True)
+    print(f"[{args.eval_mode}] rows={len(df_mode)}", flush=True)
+
+    class _ClinicalShim:
+        label_col = args.label_col
+
+        def __init__(self, df):
+            self.df = df
+
+    split_args = SimpleNamespace(
+        eval_mode=args.eval_mode,
+        holdout_per_class=args.holdout_per_class,
+        holdout_seed=args.holdout_seed,
+        n_splits=args.n_splits,
+        fold_seed=args.fold_seed,
+    )
+    splitter = PatientFirstSplitManager(patient_col=args.patient_col, label_col=args.label_col)
+    plans = splitter.build_plans(clinical=_ClinicalShim(df_mode), args=split_args)
+
+    profile_eye = build_papila_profile(
+        patient_col=args.patient_col,
+        label_col=args.label_col,
+        sample_mode="eye",
+    )
+    profile_patient = build_papila_profile(
+        patient_col=args.patient_col,
+        label_col=args.label_col,
+        sample_mode="patient",
+    )
+
+    warm_tower = int(args.warmup_tower_epochs) if args.warmup_tower_epochs is not None else 2
+    warm_fused = int(args.warmup_fused_epochs) if args.warmup_fused_epochs is not None else 2
+    total_epochs = warm_tower + warm_fused + int(args.epochs)
+
+    out_root = Path(args.analysis_dir) / args.run_name / args.eval_mode / args.tower_mode
+    out_root.mkdir(parents=True, exist_ok=True)
+
+    fold_metrics = []
+    aggregate_patient = args.tower_mode == "ensemble"
+
+    for fold_idx, split in enumerate(plans[: args.n_splits]):
+        fold_dir = out_root / f"fold{fold_idx}"
+        fold_dir.mkdir(parents=True, exist_ok=True)
+
+        eye_train = filter_eye_samples(profile_eye.build_samples(df=split.train, clinical=data))
+        bilat_val = filter_bilateral_samples(profile_patient.build_samples(df=split.val, clinical=data))
+
+        holdout_bilat = []
+        if split.holdout is not None and not split.holdout.empty:
+            holdout_bilat = filter_bilateral_samples(
+                profile_patient.build_samples(df=split.holdout, clinical=data)
+            )
+
+        if not eye_train or not bilat_val:
+            print(f"[fold {fold_idx+1}] skipped (eye_train={len(eye_train)} bilat_val={len(bilat_val)})", flush=True)
+            fold_metrics.append(
+                {
+                    "fold": fold_idx,
+                    "best_epoch": None,
+                    "best_phase": None,
+                    "val_auc": float("nan"),
+                    "val_acc": float("nan"),
+                    "hld_auc": float("nan"),
+                    "hld_acc": float("nan"),
+                    "eye_train_n": len(eye_train),
+                    "bilat_val_n": len(bilat_val),
+                    "bilat_holdout_n": len(holdout_bilat),
+                }
+            )
+            continue
+
+        sampler = build_balanced_sampler(eye_train) if args.balanced_sampling else None
+        train_loader = DataLoader(
+            EyeMetaDataset(eye_train),
+            batch_size=args.batch_size,
+            shuffle=(sampler is None),
+            sampler=sampler,
+        )
+        val_loader = DataLoader(BilatMetaDataset(bilat_val), batch_size=args.batch_size, shuffle=False)
+        holdout_loader = (
+            DataLoader(BilatMetaDataset(holdout_bilat), batch_size=args.batch_size, shuffle=False)
+            if holdout_bilat
+            else None
+        )
+
+        model = MetadataOnlySingleHT(
+            clinical_data=data,
+            num_classes=num_classes,
+            md_hidden_dim=args.md_hidden_dim,
+            fusion_dim=args.fusion_dim,
+            dropout=args.dropout,
+            use_se=bool(args.use_se),
+            se_reduction=int(args.se_reduction),
+            se_pre_norm=bool(args.se_pre_norm),
+        ).to(device)
+        opt = torch.optim.Adam(model.parameters(), lr=args.lr, weight_decay=args.weight_decay)
+
+        best_auc = -1.0
+        best_epoch = 0
+        best_phase = ""
+        best_state = None
+
+        print(
+            f"\n[fold {fold_idx+1}/{args.n_splits}] "
+            f"eye_train_n={len(eye_train)} bilat_val_n={len(bilat_val)} "
+            f"holdout_n={len(holdout_bilat)} warmup={warm_tower}+{warm_fused} total={total_epochs}",
+            flush=True,
+        )
+
+        for ep in range(total_epochs):
+            phase, main_ep = _phase_for_epoch(ep, warm_tower, warm_fused, int(args.epochs))
+            tr_loss, tr_acc = train_single_epoch(
+                model,
+                train_loader,
+                opt,
+                device,
+                phase=phase,
+                bcd_prob=float(args.bcd_prob),
+            )
+
+            _, p_val, val_auc, val_acc = _evaluate_single(
+                model,
+                val_loader,
+                device,
+                num_classes,
+                aggregate_patient=aggregate_patient,
+            )
+
+            is_main = phase == "main"
+            if is_main and (not np.isnan(val_auc)) and val_auc > best_auc:
+                best_auc = float(val_auc)
+                best_state = copy.deepcopy(model.state_dict())
+                best_epoch = ep + 1
+                best_phase = phase
+
+            if ep == 0 or (ep + 1) % 10 == 0 or (ep + 1) == total_epochs:
+                print(
+                    f"  ep {ep+1:>3}/{total_epochs} [{phase}:{main_ep}/{args.epochs}] "
+                    f"loss={tr_loss:.4f} acc={tr_acc:.4f} "
+                    f"val_auc={val_auc:.4f} val_acc={val_acc:.4f} "
+                    f"best_auc={best_auc:.4f}",
+                    flush=True,
+                )
+
+        if best_state is not None:
+            model.load_state_dict(best_state)
+
+        y_val, p_val, val_auc, val_acc = _evaluate_single(
+            model,
+            val_loader,
+            device,
+            num_classes,
+            aggregate_patient=aggregate_patient,
+        )
+
+        if args.tower_mode == "single":
+            probs_name = "probs_classic.npy"
+            probs_h_name = "probs_classic_holdout.npy"
+        else:
+            probs_name = "probs_ensemble.npy"
+            probs_h_name = "probs_ensemble_holdout.npy"
+
+        np.save(fold_dir / "y_true.npy", y_val)
+        np.save(fold_dir / probs_name, p_val)
+
+        hld_auc = float("nan")
+        hld_acc = float("nan")
+        if holdout_loader is not None:
+            y_h, p_h, hld_auc, hld_acc = _evaluate_single(
+                model,
+                holdout_loader,
+                device,
+                num_classes,
+                aggregate_patient=aggregate_patient,
+            )
+            np.save(fold_dir / "y_true_holdout.npy", y_h)
+            np.save(fold_dir / probs_h_name, p_h)
+
+        print(
+            f"  [fold {fold_idx+1}] best_epoch={best_epoch} best_auc={best_auc:.4f} "
+            f"val_auc={val_auc:.4f} val_acc={val_acc:.4f} "
+            f"hld_auc={hld_auc:.4f} hld_acc={hld_acc:.4f}",
+            flush=True,
+        )
+
+        fold_metrics.append(
+            {
+                "fold": fold_idx,
+                "best_epoch": best_epoch,
+                "best_phase": best_phase,
+                "val_auc": float(val_auc),
+                "val_acc": float(val_acc),
+                "hld_auc": float(hld_auc),
+                "hld_acc": float(hld_acc),
+                "eye_train_n": len(eye_train),
+                "bilat_val_n": len(bilat_val),
+                "bilat_holdout_n": len(holdout_bilat),
+            }
+        )
+
+    val_aucs = [m["val_auc"] for m in fold_metrics if not np.isnan(m["val_auc"])]
+    hld_aucs = [m["hld_auc"] for m in fold_metrics if not np.isnan(m["hld_auc"])]
+    if val_aucs:
+        print(f"\nMean val AUC: {np.mean(val_aucs):.4f} ± {np.std(val_aucs):.4f}", flush=True)
+    if hld_aucs:
+        print(f"Mean hld AUC: {np.mean(hld_aucs):.4f} ± {np.std(hld_aucs):.4f}", flush=True)
+
+    summary = {
+        "run_name": args.run_name,
+        "eval_mode": args.eval_mode,
+        "tower_mode": args.tower_mode,
+        "bridge_mode": "metadata_only",
+        "model": "MetadataOnlySingleHT",
+        "epochs": int(args.epochs),
+        "warmup_tower_epochs": warm_tower,
+        "warmup_fused_epochs": warm_fused,
+        "md_hidden_dim": int(args.md_hidden_dim),
+        "fusion_dim": int(args.fusion_dim),
+        "dropout": float(args.dropout),
+        "lr": float(args.lr),
+        "weight_decay": float(args.weight_decay),
+        "bcd_prob": float(args.bcd_prob),
+        "balanced_sampling": bool(args.balanced_sampling),
+        "feature_dim": int(data.feature_dim),
+        "timestamp": time.strftime("%Y%m%d_%H%M%S"),
+        "fold_metrics": fold_metrics,
+        "val_auc_mean": float(np.mean(val_aucs)) if val_aucs else None,
+        "val_auc_std": float(np.std(val_aucs)) if val_aucs else None,
+        "hld_auc_mean": float(np.mean(hld_aucs)) if hld_aucs else None,
+        "hld_auc_std": float(np.std(hld_aucs)) if hld_aucs else None,
+    }
+    (out_root / "summary.json").write_text(json.dumps(summary, indent=2))
+    print(f"\nOutputs written to: {out_root}", flush=True)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/main/v2/run_modes_roi_overnight.sh b/scripts/main/v2/run_modes_roi_overnight.sh
deleted file mode 100755
index f19e461..0000000
--- a/scripts/main/v2/run_modes_roi_overnight.sh
+++ /dev/null
@@ -1,45 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Six V2HyperTower runs using GT ROI crop:
-#   binary  × {ensemble, bilateral}               (runs 1-2)
-#   multiclass × {ensemble, bilateral}             (runs 3-4)
-#   multiclass × {ensemble, bilateral} + balanced  (runs 5-6)
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-gt
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-)
-
-# Runs 1-4: binary + multiclass, ensemble + bilateral, no balanced sampling
-echo "[1/2] GT ROI — binary + multiclass, ensemble + bilateral..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes binary multiclass \
-  --tower-modes ensemble bilateral \
-  --run-name v2_modes_gt_40ep_5fold_no_single_v2
-
-# Runs 5-6: multiclass only, ensemble + bilateral, balanced sampling
-# (reuse the crop cache built during runs 1-4)
-echo "[2/2] GT ROI — multiclass, ensemble + bilateral, balanced sampling..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes multiclass \
-  --tower-modes ensemble bilateral \
-  --balanced-sampling \
-  --persist-img-crop-cache \
-  --run-name v2_modes_gt_40ep_5fold_multiclass_balanced_v2
-
-echo "All runs complete."
diff --git a/scripts/main/v2/run_modes_roi_overnight_bilateral.sh b/scripts/main/v2/run_modes_roi_overnight_bilateral.sh
deleted file mode 100644
index 475c18c..0000000
--- a/scripts/main/v2/run_modes_roi_overnight_bilateral.sh
+++ /dev/null
@@ -1,44 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Bilateral tower runs (3 total):
-#   binary     × bilateral
-#   multiclass × bilateral
-#   multiclass × bilateral + balanced sampling
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-gt
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-  --tower-modes bilateral
-)
-
-# Runs 1-2: binary + multiclass bilateral (no balanced sampling)
-echo "[1/2] GT ROI — binary + multiclass, bilateral..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes binary multiclass \
-  --run-name v2_modes_gt_40ep_5fold_bilateral_v2
-
-# Run 3: multiclass bilateral + balanced sampling
-# (reuse the crop cache built above)
-echo "[2/2] GT ROI — multiclass, bilateral, balanced sampling..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes multiclass \
-  --balanced-sampling \
-  --persist-img-crop-cache \
-  --run-name v2_modes_gt_40ep_5fold_bilateral_balanced_v2
-
-echo "Bilateral runs complete."
diff --git a/scripts/main/v2/run_modes_roi_overnight_bilateral_unet.sh b/scripts/main/v2/run_modes_roi_overnight_bilateral_unet.sh
deleted file mode 100755
index 46d889c..0000000
--- a/scripts/main/v2/run_modes_roi_overnight_bilateral_unet.sh
+++ /dev/null
@@ -1,45 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Bilateral tower runs (3 total) — UNet ROI crop:
-#   binary     × bilateral
-#   multiclass × bilateral
-#   multiclass × bilateral + balanced sampling
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-weights models/v2/refuge/segmentation/per_image_refuge_build/best.pt
-  --img-crop-normalize per_image
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-  --tower-modes bilateral
-)
-
-# Runs 1-2: binary + multiclass bilateral (no balanced sampling)
-echo "[1/2] UNet ROI — binary + multiclass, bilateral..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes binary multiclass \
-  --run-name v2_modes_unet_40ep_5fold_bilateral_v2
-
-# Run 3: multiclass bilateral + balanced sampling
-# (reuse the crop cache built above)
-echo "[2/2] UNet ROI — multiclass, bilateral, balanced sampling..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes multiclass \
-  --balanced-sampling \
-  --persist-img-crop-cache \
-  --run-name v2_modes_unet_40ep_5fold_bilateral_balanced_v2
-
-echo "Bilateral UNet runs complete."
diff --git a/scripts/main/v2/run_modes_roi_overnight_ensemble.sh b/scripts/main/v2/run_modes_roi_overnight_ensemble.sh
deleted file mode 100644
index 0def8a0..0000000
--- a/scripts/main/v2/run_modes_roi_overnight_ensemble.sh
+++ /dev/null
@@ -1,44 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Ensemble tower runs (3 total):
-#   binary     × ensemble
-#   multiclass × ensemble
-#   multiclass × ensemble + balanced sampling
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-gt
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-  --tower-modes ensemble
-)
-
-# Runs 1-2: binary + multiclass ensemble (no balanced sampling)
-echo "[1/2] GT ROI — binary + multiclass, ensemble..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes binary multiclass \
-  --run-name v2_modes_gt_40ep_5fold_ensemble_v2
-
-# Run 3: multiclass ensemble + balanced sampling
-# (reuse the crop cache built above)
-echo "[2/2] GT ROI — multiclass, ensemble, balanced sampling..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes multiclass \
-  --balanced-sampling \
-  --persist-img-crop-cache \
-  --run-name v2_modes_gt_40ep_5fold_ensemble_balanced_v2
-
-echo "Ensemble runs complete."
diff --git a/scripts/main/v2/run_modes_roi_overnight_ensemble_unet.sh b/scripts/main/v2/run_modes_roi_overnight_ensemble_unet.sh
deleted file mode 100755
index 7e725b8..0000000
--- a/scripts/main/v2/run_modes_roi_overnight_ensemble_unet.sh
+++ /dev/null
@@ -1,45 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Ensemble tower runs (3 total) — UNet ROI crop:
-#   binary     × ensemble
-#   multiclass × ensemble
-#   multiclass × ensemble + balanced sampling
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-CROP_ARGS=(
-  --img-crop-manifest manifest.csv
-  --img-crop-weights models/v2/refuge/segmentation/per_image_refuge_build/best.pt
-  --img-crop-normalize per_image
-)
-
-COMMON_ARGS=(
-  --epochs 40
-  --n-splits 5
-  --batch-size 8
-  --backbone refugelike
-  --tower-modes ensemble
-)
-
-# Runs 1-2: binary + multiclass ensemble (no balanced sampling)
-echo "[1/2] UNet ROI — binary + multiclass, ensemble..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes binary multiclass \
-  --run-name v2_modes_unet_40ep_5fold_ensemble_v2
-
-# Run 3: multiclass ensemble + balanced sampling
-# (reuse the crop cache built above)
-echo "[2/2] UNet ROI — multiclass, ensemble, balanced sampling..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  "${CROP_ARGS[@]}" \
-  --eval-modes multiclass \
-  --balanced-sampling \
-  --persist-img-crop-cache \
-  --run-name v2_modes_unet_40ep_5fold_ensemble_balanced_v2
-
-echo "Ensemble UNet runs complete."
diff --git a/scripts/main/v2/run_modes_roi_smoke.sh b/scripts/main/v2/run_modes_roi_smoke.sh
deleted file mode 100755
index ca64b4c..0000000
--- a/scripts/main/v2/run_modes_roi_smoke.sh
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/usr/bin/env bash
-set -euo pipefail
-
-# Quick smoke test for ROI mode runs:
-#   1) GT masks
-#   2) UNet masks
-# Uses 1 epoch and 1 fold for fast validation.
-
-ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
-cd "$ROOT_DIR"
-
-COMMON_ARGS=(
-  --eval-modes binary multiclass
-  --tower-modes single ensemble bilateral
-  --epochs 1
-  --n-splits 2
-  --folds 1
-  --batch-size 8
-  --backbone refugelike
-  --img-crop-manifest manifest.csv
-  --warmup-tower-epochs 0
-  --warmup-fused-epochs 0
-  --single-warmup-tower-epochs 0
-  --single-warmup-fused-epochs 0
-  --bilat-warmup-tower-epochs 0
-  --bilat-warmup-fused-epochs 0
-)
-
-echo "[smoke 1/2] Starting GT ROI run..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  --img-crop-gt \
-  --run-name smoke_v2_modes_roi_gt
-
-echo "[smoke 2/2] Starting UNet ROI run..."
-python3 scripts/basic_analysis/compare_hypertower_modes.py \
-  "${COMMON_ARGS[@]}" \
-  --img-crop-weights models/v2/refuge/segmentation/per_image_refuge_build/best.pt \
-  --img-crop-normalize per_image \
-  --run-name smoke_v2_modes_roi_unet_perimage
-
-echo "Smoke runs complete."
diff --git a/scripts/main/v2/run_multiclass_v2.2.sh b/scripts/main/v2/run_multiclass_v2.2.sh
new file mode 100755
index 0000000..64a8be1
--- /dev/null
+++ b/scripts/main/v2/run_multiclass_v2.2.sh
@@ -0,0 +1,73 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+# Multiclass runs v2.2 (6 total):
+#   UNet crop:  single | ensemble | fused head
+#   GT crop:    single | ensemble | fused head
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)"
+cd "$ROOT_DIR"
+
+MANIFEST="manifest.csv"
+UNET_WEIGHTS="models/v2/refuge/segmentation/per_image/best.pt"
+
+COMMON=(
+  --epochs 40
+  --n-splits 5
+  --batch-size 8
+  --backbone refugelike
+  --eval-mode multiclass
+  --img-crop-manifest "$MANIFEST"
+)
+
+UNET_CROP=(
+  --img-crop-weights "$UNET_WEIGHTS"
+)
+
+GT_CROP=(
+  --img-crop-gt
+)
+
+# ── UNet crop ────────────────────────────────────────────────────────────────
+
+echo "[1/6] UNet crop — multiclass, single..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode single \
+  --run-name v2.2_single_multiclass_unet_40ep_5fold
+
+echo "[2/6] UNet crop — multiclass, ensemble..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode ensemble \
+  --run-name v2.2_ensemble_multiclass_unet_40ep_5fold
+
+echo "[3/6] UNet crop — multiclass, ensemble + fused head..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${UNET_CROP[@]}" \
+  --tower-mode ensemble \
+  --fused-head --fusion-epochs 10 \
+  --run-name v2.2_fused_multiclass_unet_40ep_5fold
+
+# ── GT crop ──────────────────────────────────────────────────────────────────
+
+echo "[4/6] GT crop — multiclass, single..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode single \
+  --run-name v2.2_single_multiclass_gt_40ep_5fold
+
+echo "[5/6] GT crop — multiclass, ensemble..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode ensemble \
+  --run-name v2.2_ensemble_multiclass_gt_40ep_5fold
+
+echo "[6/6] GT crop — multiclass, ensemble + fused head..."
+python3 scripts/main/v2/run_multifold_v2.py \
+  "${COMMON[@]}" "${GT_CROP[@]}" \
+  --tower-mode ensemble \
+  --fused-head --fusion-epochs 10 \
+  --run-name v2.2_fused_multiclass_gt_40ep_5fold
+
+echo "Multiclass v2.2 runs complete."
diff --git a/scripts/output_analysis/explainability/explain_fold.py b/scripts/output_analysis/explainability/explain_fold.py
index 95ee358..6122961 100644
--- a/scripts/output_analysis/explainability/explain_fold.py
+++ b/scripts/output_analysis/explainability/explain_fold.py
@@ -197,20 +197,30 @@ def run_permutation_importance(
 ) -> None:
     print("\n[Phase 1] MD permutation importance ...", flush=True)
 
-    # ---- cache image embeddings + collect meta tensors + labels ----
-    img_feats_list, md_list, label_list = [], [], []
+    # ---- cache bilateral image embeddings + metadata tensors + labels ----
+    img1_feats_list, img2_feats_list = [], []
+    md1_list, md2_list, label_list = [], [], []
     model.eval()
     with torch.no_grad():
         for batch in loader:
-            imgs = batch["image_1"].to(device)
-            meta = batch["matrix_1"].to(device)
+            img1 = batch["image_1"].to(device)
+            img2 = batch["image_2"].to(device)
+            md1 = batch["matrix_1"].to(device)
+            md2 = batch["matrix_2"].to(device)
             labels = batch["label_1"]
-            img_feats_list.append(model.img_tower(imgs))
-            md_list.append(meta)
-            label_list.append(labels)
+            img1_feats_list.append(model.img_tower(img1))
+            img2_feats_list.append(model.img_tower(img2))
+            md1_list.append(md1)
+            md2_list.append(md2)
+            if isinstance(labels, torch.Tensor):
+                label_list.append(labels)
+            else:
+                label_list.append(torch.tensor(labels, dtype=torch.long))
 
-    img_feats = torch.cat(img_feats_list)  # [N, img_dim]
-    md_tensor = torch.cat(md_list)  # [N, feature_dim]
+    img1_feats = torch.cat(img1_feats_list)  # [N, img_dim]
+    img2_feats = torch.cat(img2_feats_list)  # [N, img_dim]
+    md1_tensor = torch.cat(md1_list)  # [N, feature_dim]
+    md2_tensor = torch.cat(md2_list)  # [N, feature_dim]
     y_true = torch.cat(label_list).numpy()
     N = len(y_true)
 
@@ -218,11 +228,15 @@ def run_permutation_importance(
         print("  [Phase 1] No samples — skipping.", flush=True)
         return
 
-    # ---- baseline AUC ----
+    # ---- baseline AUC (patient-level: average OD/OS fused probabilities) ----
     with torch.no_grad():
-        md_feats = model.md_tower(md_tensor)
-        fused, _, _ = model.bridge(img_feats, md_feats)
-        probs_baseline = torch.softmax(fused, dim=1).cpu().numpy()
+        md1_feats = model.md_tower(md1_tensor)
+        md2_feats = model.md_tower(md2_tensor)
+        fused1, _, _ = model.bridge(img1_feats, md1_feats)
+        fused2, _, _ = model.bridge(img2_feats, md2_feats)
+        probs_baseline = (
+            0.5 * (torch.softmax(fused1, dim=1) + torch.softmax(fused2, dim=1))
+        ).cpu().numpy()
     _, baseline_auc, _ = _score_arrays(y_true, probs_baseline, num_classes)
     print(f"  Baseline AUC: {baseline_auc:.4f}  (N={N})", flush=True)
 
@@ -235,13 +249,25 @@ def run_permutation_importance(
         all_dims = dims["value_dims"] + dims["missing_dims"]
         drops = []
         for _ in range(n_permutations):
-            perm = md_tensor.clone()
+            perm1 = md1_tensor.clone()
+            perm2 = md2_tensor.clone()
             perm_idx = torch.from_numpy(rng.permutation(N)).to(device)
-            perm[:, all_dims] = perm[perm_idx][:, all_dims]
+            # Apply the same donor patient permutation to both eyes to preserve
+            # within-patient coherence while breaking feature-label association.
+            perm1[:, all_dims] = perm1[perm_idx][:, all_dims]
+            perm2[:, all_dims] = perm2[perm_idx][:, all_dims]
             with torch.no_grad():
-                md_p = model.md_tower(perm)
-                fused_p, _, _ = model.bridge(img_feats, md_p)
-                probs_p = torch.softmax(fused_p, dim=1).cpu().numpy()
+                md1_p = model.md_tower(perm1)
+                md2_p = model.md_tower(perm2)
+                fused1_p, _, _ = model.bridge(img1_feats, md1_p)
+                fused2_p, _, _ = model.bridge(img2_feats, md2_p)
+                probs_p = (
+                    0.5
+                    * (
+                        torch.softmax(fused1_p, dim=1)
+                        + torch.softmax(fused2_p, dim=1)
+                    )
+                ).cpu().numpy()
             _, auc_p, _ = _score_arrays(y_true, probs_p, num_classes)
             drops.append(baseline_auc - auc_p)
 
@@ -254,6 +280,56 @@ def run_permutation_importance(
 
     results.sort(key=lambda r: r["importance"], reverse=True)
 
+    # ---- total MD ablation (all features permuted simultaneously) ----
+    print("  Running total MD ablation ...", flush=True)
+    total_drops = []
+    for _ in range(n_permutations):
+        perm_idx = torch.from_numpy(rng.permutation(N)).to(device)
+        perm1_all = md1_tensor[perm_idx]
+        perm2_all = md2_tensor[perm_idx]
+        with torch.no_grad():
+            md1_all = model.md_tower(perm1_all)
+            md2_all = model.md_tower(perm2_all)
+            f1, _, _ = model.bridge(img1_feats, md1_all)
+            f2, _, _ = model.bridge(img2_feats, md2_all)
+            probs_all = (
+                0.5 * (torch.softmax(f1, dim=1) + torch.softmax(f2, dim=1))
+            ).cpu().numpy()
+        _, auc_all, _ = _score_arrays(y_true, probs_all, num_classes)
+        total_drops.append(baseline_auc - auc_all)
+    total_mean = float(np.mean(total_drops))
+    total_std = float(np.std(total_drops))
+    print(
+        f"  Total MD ablation  Δ AUC = {total_mean:+.4f} ± {total_std:.4f}", flush=True
+    )
+
+    # ---- Gaussian noise ablation (tests architectural vs informational benefit) ----
+    print("  Running Gaussian noise ablation ...", flush=True)
+    noise_drops = []
+    for _ in range(n_permutations):
+        noise1 = torch.randn_like(md1_tensor)
+        noise2 = torch.randn_like(md2_tensor)
+        with torch.no_grad():
+            md1_noise = model.md_tower(noise1)
+            md2_noise = model.md_tower(noise2)
+            f1, _, _ = model.bridge(img1_feats, md1_noise)
+            f2, _, _ = model.bridge(img2_feats, md2_noise)
+            probs_noise = (
+                0.5 * (torch.softmax(f1, dim=1) + torch.softmax(f2, dim=1))
+            ).cpu().numpy()
+        _, auc_noise, _ = _score_arrays(y_true, probs_noise, num_classes)
+        noise_drops.append(baseline_auc - auc_noise)
+    noise_mean = float(np.mean(noise_drops))
+    noise_std = float(np.std(noise_drops))
+    print(
+        f"  Gaussian noise ablation  Δ AUC = {noise_mean:+.4f} ± {noise_std:.4f}", flush=True
+    )
+    print(
+        f"  [interpretation] permutation Δ={total_mean:+.4f}  noise Δ={noise_mean:+.4f}  "
+        f"informational gain = {total_mean - noise_mean:+.4f}",
+        flush=True,
+    )
+
     # ---- save CSV ----
     import csv
 
@@ -262,6 +338,8 @@ def run_permutation_importance(
         writer = csv.DictWriter(f, fieldnames=["feature", "importance", "std"])
         writer.writeheader()
         writer.writerows(results)
+        writer.writerow({"feature": "TOTAL_MD_ABLATION", "importance": total_mean, "std": total_std})
+        writer.writerow({"feature": "GAUSSIAN_NOISE_ABLATION", "importance": noise_mean, "std": noise_std})
 
     # ---- bar chart ----
     names = [r["feature"] for r in results]
@@ -269,13 +347,26 @@ def run_permutation_importance(
     stds = [r["std"] for r in results]
     colors = ["#e05c5c" if v >= 0 else "#5c9ee0" for v in imps]
 
-    fig, ax = plt.subplots(figsize=(9, max(4, len(names) * 0.45)))
+    fig, ax = plt.subplots(figsize=(9, max(4, (len(names) + 3) * 0.45)))
     y_pos = np.arange(len(names))
-    bars = ax.barh(
+    ax.barh(
         y_pos, imps, xerr=stds, color=colors, ecolor="grey", capsize=3, height=0.6
     )
-    ax.set_yticks(y_pos)
-    ax.set_yticklabels(names, fontsize=9)
+    ax.axhline(len(names) - 0.25, color="grey", linewidth=0.6, linestyle="--")
+    # total ablation
+    ax.barh(
+        len(names) + 0.5, total_mean, xerr=total_std,
+        color="#c45ce0" if total_mean >= 0 else "#5c9ee0",
+        ecolor="grey", capsize=3, height=0.6,
+    )
+    # gaussian noise ablation
+    ax.barh(
+        len(names) + 1.5, noise_mean, xerr=noise_std,
+        color="#e08c2a" if noise_mean >= 0 else "#5c9ee0",
+        ecolor="grey", capsize=3, height=0.6,
+    )
+    ax.set_yticks(list(y_pos) + [len(names) + 0.5, len(names) + 1.5])
+    ax.set_yticklabels(names + ["ALL MD (permute)", "ALL MD (noise)"], fontsize=9)
     ax.invert_yaxis()
     ax.axvline(0, color="black", linewidth=0.8)
     ax.set_xlabel("Mean AUC drop (baseline − permuted)", fontsize=10)
@@ -325,7 +416,8 @@ def run_gradcam(
         img_os = batch["image_2"].to(device)  # [1, 3, H, W]
         meta_od = batch["matrix_1"].to(device)  # [1, feature_dim]
         meta_os = batch["matrix_2"].to(device)
-        label = int(batch["label_1"][0].item())
+        lbl_raw = batch["label_1"][0]
+        label = int(lbl_raw.item() if isinstance(lbl_raw, torch.Tensor) else lbl_raw)
         pid = batch["id_1"][0]
 
         # GradCAM for each eye (OD drives the prediction label)
@@ -488,9 +580,149 @@ def parse_args():
     ap.add_argument("--batch-size", type=int, default=1)
     ap.add_argument("--no-phase1", action="store_true", help="Skip MD importance")
     ap.add_argument("--no-phase2", action="store_true", help="Skip GradCAM")
+    ap.add_argument("--no-phase3", action="store_true", help="Skip fusion event analysis")
     return ap.parse_args()
 
 
+# ---------------------------------------------------------------------------
+# Phase 3 — Fusion event analysis
+# ---------------------------------------------------------------------------
+
+
+def run_fusion_event_analysis(
+    model: SingleEyeHT,
+    loader,
+    device: torch.device,
+    out_dir: Path,
+) -> None:
+    print("\n[Phase 3] Fusion event analysis ...", flush=True)
+
+    from classes.v2.models import collect_probs_single_components
+
+    y_true, pf, pi, pm = collect_probs_single_components(
+        model, loader, device, aggregate_patient=True
+    )
+    N = len(y_true)
+    if N == 0:
+        print("  [Phase 3] No samples — skipping.", flush=True)
+        return
+
+    pred_f = pf.argmax(axis=1)
+    pred_i = pi.argmax(axis=1)
+    pred_m = pm.argmax(axis=1)
+
+    corrections = (pred_f == y_true) & (pred_i != y_true) & (pred_m != y_true)
+    errors      = (pred_f != y_true) & (pred_i == y_true) & (pred_m == y_true)
+    n_corr = corrections.sum()
+    n_err  = errors.sum()
+    both_wrong   = ((pred_i != y_true) & (pred_m != y_true)).sum()
+    both_correct = ((pred_i == y_true) & (pred_m == y_true)).sum()
+
+    print(f"  N={N}  corrections={n_corr}  errors={n_err}  ratio={n_corr}/{n_err}", flush=True)
+    print(f"  correction rate: {n_corr}/{both_wrong} = {n_corr/max(both_wrong,1):.2%} of both-wrong cases", flush=True)
+    print(f"  error rate:      {n_err}/{both_correct} = {n_err/max(both_correct,1):.2%} of both-correct cases", flush=True)
+
+    # ---- cache intermediate hm/hi vectors for all patients ----
+    model.eval()
+    hm1_list, hm2_list, hi1_list, hi2_list = [], [], [], []
+    with torch.no_grad():
+        for batch in loader:
+            x1 = batch.get("image_1"); m1 = batch.get("matrix_1")
+            x2 = batch.get("image_2"); m2 = batch.get("matrix_2")
+            if not (torch.is_tensor(x1) and torch.is_tensor(m1)):
+                continue
+            hi1 = model.bridge.ln_img(model.bridge.W_img(model.img_tower(x1.to(device))))
+            hi2 = model.bridge.ln_img(model.bridge.W_img(model.img_tower(x2.to(device))))
+            hm1 = model.bridge.ln_md(model.bridge.W_md(model.md_tower(m1.to(device))))
+            hm2 = model.bridge.ln_md(model.bridge.W_md(model.md_tower(m2.to(device))))
+            hi1_list.append(hi1.cpu()); hi2_list.append(hi2.cpu())
+            hm1_list.append(hm1.cpu()); hm2_list.append(hm2.cpu())
+
+    hi1 = torch.cat(hi1_list)  # [N, fusion_dim]
+    hi2 = torch.cat(hi2_list)
+    hm1 = torch.cat(hm1_list)  # [N, fusion_dim]
+    hm2 = torch.cat(hm2_list)
+
+    hm1_mean = hm1.mean(dim=0, keepdim=True)
+    hm2_mean = hm2.mean(dim=0, keepdim=True)
+
+    # ---- for each patient: compare logit[true_class] with real hm vs mean hm ----
+    gains = []
+    with torch.no_grad():
+        for idx in range(N):
+            true_cls = int(y_true[idx])
+            # patient-level average of OD/OS fused vectors (SE skipped: hard to replicate outside forward)
+            fused_real = (hi1[idx:idx+1] * hm1[idx:idx+1] + hi2[idx:idx+1] * hm2[idx:idx+1]) * 0.5
+            fused_mean = (hi1[idx:idx+1] * hm1_mean + hi2[idx:idx+1] * hm2_mean) * 0.5
+            logit_real = model.bridge.classifier_fused(fused_real.to(device))
+            logit_mean = model.bridge.classifier_fused(fused_mean.to(device))
+            gain = (logit_real[0, true_cls] - logit_mean[0, true_cls]).item()
+            gains.append(gain)
+
+    gains = np.array(gains)
+
+    if n_corr > 0:
+        corr_gains = gains[corrections]
+        helped = (corr_gains > 0).sum()
+        print(f"\n  Fusion corrections — MD gate gain vs mean gate:", flush=True)
+        print(f"    mean gain = {corr_gains.mean():+.4f}  median = {np.median(corr_gains):+.4f}", flush=True)
+        print(f"    real MD helped {helped}/{n_corr} correction patients ({helped/n_corr:.0%})", flush=True)
+
+    if n_err > 0:
+        err_gains = gains[errors]
+        print(f"\n  Fusion errors — MD gate gain vs mean gate:", flush=True)
+        print(f"    mean gain = {err_gains.mean():+.4f}  median = {np.median(err_gains):+.4f}", flush=True)
+
+    # ---- save CSV ----
+    import csv
+    rows = []
+    for idx in range(N):
+        rows.append({
+            "patient_idx": idx,
+            "y_true": int(y_true[idx]),
+            "pred_fused": int(pred_f[idx]),
+            "pred_img": int(pred_i[idx]),
+            "pred_md": int(pred_m[idx]),
+            "conf_fused": float(pf[idx].max()),
+            "conf_img": float(pi[idx].max()),
+            "conf_md": float(pm[idx].max()),
+            "is_correction": bool(corrections[idx]),
+            "is_error": bool(errors[idx]),
+            "md_gate_gain": float(gains[idx]),
+        })
+    csv_path = out_dir / "fusion_events.csv"
+    with csv_path.open("w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
+        writer.writeheader()
+        writer.writerows(rows)
+    print(f"  Saved → {csv_path}", flush=True)
+
+    # ---- chart ----
+    fig, axes = plt.subplots(1, 2, figsize=(11, 4))
+
+    categories = ["corrections\n(both wrong→fused right)", "errors\n(both right→fused wrong)"]
+    counts = [int(n_corr), int(n_err)]
+    axes[0].bar(categories, counts, color=["#e05c5c", "#5c9ee0"], width=0.5)
+    axes[0].set_ylabel("Count")
+    axes[0].set_title(f"Fusion Events  (N={N})")
+    for i, v in enumerate(counts):
+        axes[0].text(i, v + 0.1, str(v), ha="center", fontsize=11)
+
+    if n_corr > 0:
+        axes[1].hist(gains[corrections], bins=10, alpha=0.7, color="#e05c5c", label=f"corrections (n={n_corr})")
+    if n_err > 0:
+        axes[1].hist(gains[errors], bins=10, alpha=0.7, color="#5c9ee0", label=f"errors (n={n_err})")
+    axes[1].axvline(0, color="black", linewidth=0.8)
+    axes[1].set_xlabel("MD gate gain vs mean gate\n(logit[true class]: real − mean)")
+    axes[1].set_title("Does real MD help the fused prediction?")
+    axes[1].legend(fontsize=9)
+
+    fig.tight_layout()
+    fig.savefig(out_dir / "fusion_events.png", dpi=150)
+    plt.close(fig)
+    print(f"  Saved → {out_dir / 'fusion_events.png'}", flush=True)
+
+
 def main():
     args = parse_args()
     fold_dir = args.fold_dir.resolve()
@@ -639,6 +871,15 @@ def main():
             out_dir=out_dir,
         )
 
+    # ---- Phase 3 ----
+    if not args.no_phase3:
+        run_fusion_event_analysis(
+            model=model,
+            loader=loader,
+            device=device,
+            out_dir=out_dir,
+        )
+
     print("\n[explain_fold] Done.", flush=True)
 
 
diff --git a/scripts/output_analysis/visualizations/aggregate_roc_perclass_all_models_v2.py b/scripts/output_analysis/visualizations/aggregate_roc_perclass_all_models_v2.py
new file mode 100755
index 0000000..eedcf82
--- /dev/null
+++ b/scripts/output_analysis/visualizations/aggregate_roc_perclass_all_models_v2.py
@@ -0,0 +1,281 @@
+#!/usr/bin/env python3
+"""
+Per-class ROC curves for v2 HyperTower runs.
+
+Plots multiple runs as separate lines on the same axes — one figure per class
+(multiclass) or one figure total (binary).
+
+v2 directory layout
+-------------------
+  analysis_data/{run_name}/{eval_mode}/{tower_mode}/
+    fold0/  y_true.npy  probs_fused.npy | probs_bilat.npy | probs_classic.npy | probs_fused_head.npy
+    fold1/  ...
+
+Usage examples
+--------------
+  # Compare UNet ensemble vs bilateral vs fused head (binary)
+  python scripts/output_analysis/visualizations/aggregate_roc_perclass_all_models_v2.py \\
+    --mode binary --tag unet_binary_comparison \\
+    --runs \\
+      analysis_data/v2_modes_full_40ep_5fold_roi_unet_perimage_refugebuild_holdout/binary/ensemble:"UNet Ensemble" \\
+      analysis_data/v2_modes_full_40ep_5fold_roi_unet_perimage_refugebuild_holdout/binary/bilateral:"UNet Bilateral" \\
+      analysis_data/v2_ensemble_fused_binary_unet_40ep_5fold_v1/binary/ensemble:"UNet Fused Head"
+
+  Each --runs entry is  <path>:<label>  where <path> points directly to the
+  {eval_mode}/{tower_mode} subdirectory and <label> is shown in the legend.
+
+Outputs (written to --output-dir, default: analysis_data/roc_plots/)
+  {tag}_binary_roc.png          (binary mode)
+  {tag}_class{k}_roc.png        (multiclass mode, one file per class)
+  {tag}_roc_summary.json
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.metrics import roc_curve, auc as sk_auc
+
+
+# ---------------------------------------------------------------------------
+# Probs filename auto-detection priority per tower mode
+# ---------------------------------------------------------------------------
+_PROBS_PRIORITY: dict[str, list[str]] = {
+    "ensemble":  ["probs_fused_head", "probs_fused"],
+    "bilateral": ["probs_bilat"],
+    "single":    ["probs_classic"],
+    "classic":   ["probs_classic"],
+}
+_ALL_PROBS = ["probs_fused_head", "probs_fused", "probs_bilat", "probs_classic"]
+
+
+def _detect_probs_stem(fold_dir: Path, tower_mode: str | None) -> str | None:
+    priority = _PROBS_PRIORITY.get(tower_mode, _ALL_PROBS) if tower_mode else _ALL_PROBS
+    for stem in priority:
+        if (fold_dir / f"{stem}.npy").exists():
+            return stem
+    return None
+
+
+# ---------------------------------------------------------------------------
+# Fold discovery and loading
+# ---------------------------------------------------------------------------
+
+def find_fold_dirs(mode_dir: Path) -> list[Path]:
+    return sorted(
+        [p for p in mode_dir.iterdir() if p.is_dir() and p.name.startswith("fold")],
+        key=lambda p: int(p.name.replace("fold", "")),
+    )
+
+
+def load_fold(fold_dir: Path, probs_stem: str) -> tuple[np.ndarray, np.ndarray] | None:
+    y_path = fold_dir / "y_true.npy"
+    p_path = fold_dir / f"{probs_stem}.npy"
+    if not y_path.exists() or not p_path.exists():
+        return None
+    return np.load(y_path), np.load(p_path)
+
+
+# ---------------------------------------------------------------------------
+# Per-class ROC helpers
+# ---------------------------------------------------------------------------
+
+def per_class_roc(y: np.ndarray, p: np.ndarray) -> dict[int, tuple]:
+    K = p.shape[1]
+    out: dict[int, tuple] = {}
+    for k in range(K):
+        yb = (y == k).astype(np.uint8)
+        if yb.sum() == 0 or yb.sum() == len(yb):
+            continue
+        fpr, tpr, _ = roc_curve(yb, p[:, k])
+        out[k] = (fpr, tpr, sk_auc(fpr, tpr))
+    return out
+
+
+def build_mean_curve(mode_dir: Path, tower_mode: str | None, mode: str) -> dict | None:
+    fold_dirs = find_fold_dirs(mode_dir)
+    if not fold_dirs:
+        return None
+
+    probs_stem: str | None = None
+    for fd in fold_dirs:
+        probs_stem = _detect_probs_stem(fd, tower_mode)
+        if probs_stem:
+            break
+    if probs_stem is None:
+        return None
+
+    grid = np.linspace(0, 1, 501)
+    per_fold: list[dict] = []
+    for fd in fold_dirs:
+        result = load_fold(fd, probs_stem)
+        if result is None:
+            continue
+        y, p = result
+        if mode == "binary":
+            mask = np.isin(y, [0, 1])
+            y, p = y[mask], p[mask]
+            if p.shape[1] > 2:
+                p = p[:, :2]
+        per_fold.append(per_class_roc(y, p))
+
+    if not per_fold:
+        return None
+
+    K = max(max(d.keys()) for d in per_fold) + 1
+    class_curves: dict[int, dict] = {}
+    for k in range(K):
+        tprs, aucs = [], []
+        for d in per_fold:
+            if k not in d:
+                continue
+            fpr, tpr, a = d[k]
+            tprs.append(np.interp(grid, fpr, tpr))
+            aucs.append(a)
+        if not tprs:
+            continue
+        tprs_arr = np.vstack(tprs)
+        class_curves[k] = {
+            "fpr":      grid,
+            "tpr_mean": tprs_arr.mean(axis=0),
+            "tpr_std":  tprs_arr.std(axis=0),
+            "auc_mean": float(np.nanmean(aucs)),
+            "auc_std":  float(np.nanstd(aucs)),
+        }
+
+    return {"probs_stem": probs_stem, "class_curves": class_curves}
+
+
+# ---------------------------------------------------------------------------
+# CLI
+# ---------------------------------------------------------------------------
+
+def parse_run_entry(entry: str) -> tuple[Path, str]:
+    """Parse  path:label  or  path  (label defaults to last two dir components)."""
+    if ":" in entry:
+        raw_path, label = entry.rsplit(":", 1)
+    else:
+        raw_path = entry
+        p = Path(entry)
+        label = f"{p.parent.name}/{p.name}"
+    return Path(raw_path), label
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(
+        description="Per-class ROC curves comparing multiple v2 HyperTower runs.",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=__doc__,
+    )
+    ap.add_argument(
+        "--runs", nargs="+", required=True, metavar="PATH[:LABEL]",
+        help=(
+            "Mode-level directories to compare, each optionally followed by :label. "
+            "Path should point to the {eval_mode}/{tower_mode} subdirectory."
+        ),
+    )
+    ap.add_argument("--mode", required=True, choices=["binary", "multiclass"])
+    ap.add_argument("--tag", required=True, help="Output filename prefix.")
+    ap.add_argument(
+        "--output-dir", default="analysis_data/roc_plots",
+        help="Directory for PNG and JSON output (default: analysis_data/roc_plots).",
+    )
+    ap.add_argument(
+        "--class-names", nargs="*", default=["Healthy", "Glaucoma", "Suspect"],
+    )
+    ap.add_argument("--shade", action="store_true", help="Shade ±1 SD bands.")
+    args = ap.parse_args()
+
+    out_dir = Path(args.output_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    per_model: list[tuple[str, dict]] = []
+    for entry in args.runs:
+        mode_dir, label = parse_run_entry(entry)
+        if not mode_dir.exists():
+            print(f"  WARNING: {mode_dir} not found — skipping.")
+            continue
+        tower_mode = mode_dir.name
+        result = build_mean_curve(mode_dir, tower_mode, args.mode)
+        if result is None:
+            print(f"  WARNING: no usable folds in {mode_dir} — skipping.")
+            continue
+        auc_str = "  ".join(
+            f"class{k} AUC={v['auc_mean']:.3f}±{v['auc_std']:.3f}"
+            for k, v in result["class_curves"].items()
+        )
+        print(f"  {label}  [{result['probs_stem']}]  {auc_str}")
+        per_model.append((label, result["class_curves"]))
+
+    if not per_model:
+        raise SystemExit("No usable runs — nothing to plot.")
+
+    if args.mode == "binary":
+        classes_to_plot = [1]
+        out_names = [f"{args.tag}_binary_roc.png"]
+        titles = ["Binary — Glaucoma (positive class)"]
+    else:
+        max_k = max(max(curves.keys()) for _, curves in per_model)
+        classes_to_plot = list(range(min(3, max_k + 1)))
+        out_names = [f"{args.tag}_class{k}_roc.png" for k in classes_to_plot]
+        titles = [
+            f"Multiclass OVR — "
+            f"{args.class_names[k] if k < len(args.class_names) else f'class {k}'}"
+            for k in classes_to_plot
+        ]
+
+    out_json: dict = {"tag": args.tag, "mode": args.mode, "figures": []}
+
+    for k, out_name, title in zip(classes_to_plot, out_names, titles):
+        fig, ax = plt.subplots(figsize=(9, 7))
+        ax.plot([0, 1], [0, 1], linestyle="--", linewidth=1, color="grey")
+        ax.set_xlabel("False Positive Rate")
+        ax.set_ylabel("True Positive Rate")
+        ax.set_title(f"{title}\n{args.tag}")
+
+        entries = []
+        for label, curves in per_model:
+            if k not in curves:
+                continue
+            c = curves[k]
+            ax.plot(
+                c["fpr"], c["tpr_mean"], linewidth=2,
+                label=f"{label}  (AUC {c['auc_mean']:.3f} ± {c['auc_std']:.3f})",
+            )
+            if args.shade:
+                ax.fill_between(
+                    c["fpr"],
+                    np.clip(c["tpr_mean"] - c["tpr_std"], 0, 1),
+                    np.clip(c["tpr_mean"] + c["tpr_std"], 0, 1),
+                    alpha=0.10,
+                )
+            entries.append({
+                "label": label,
+                "auc_mean": c["auc_mean"],
+                "auc_std": c["auc_std"],
+            })
+
+        ax.legend(loc="lower right")
+        fig.tight_layout()
+        out_path = out_dir / out_name
+        fig.savefig(out_path, dpi=160)
+        plt.close(fig)
+        print(f"  Saved: {out_path}")
+
+        out_json["figures"].append({
+            "class_index": k,
+            "output_png": str(out_path),
+            "models": entries,
+        })
+
+    summary_path = out_dir / f"{args.tag}_roc_summary.json"
+    summary_path.write_text(json.dumps(out_json, indent=2))
+    print(f"  Summary: {summary_path}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/output_analysis/visualizations/plot_run_roc_v2.py b/scripts/output_analysis/visualizations/plot_run_roc_v2.py
new file mode 100644
index 0000000..03a3b8f
--- /dev/null
+++ b/scripts/output_analysis/visualizations/plot_run_roc_v2.py
@@ -0,0 +1,286 @@
+#!/usr/bin/env python3
+"""
+Per-fold and mean OVR ROC plots for a single v2 HyperTower run.
+
+Reads the saved .npy artifacts from a completed run and produces:
+  - One per-fold ROC figure per class (all folds as individual lines)
+  - One mean ± SD OVR ROC figure (all classes on the same axes)
+
+Outputs are written to {mode_dir}/plots/.
+
+v2 directory layout expected
+-----------------------------
+  {run_dir}/{eval_mode}/{tower_mode}/
+    fold0/  y_true.npy  probs_fused.npy | probs_bilat.npy | probs_classic.npy | probs_fused_head.npy
+    fold1/  ...
+
+Usage examples
+--------------
+  # Ensemble binary run
+  python scripts/output_analysis/visualizations/plot_run_roc_v2.py \\
+      --run-dir analysis_data/v2_modes_full_40ep_5fold_roi_unet_perimage_refugebuild_holdout \\
+      --eval-mode binary --tower-mode ensemble
+
+  # Bilateral multiclass run
+  python scripts/output_analysis/visualizations/plot_run_roc_v2.py \\
+      --run-dir analysis_data/v2_modes_full_40ep_5fold_roi_unet_perimage_refugebuild_holdout \\
+      --eval-mode multiclass --tower-mode bilateral
+
+  # Fused head — explicitly select probs file
+  python scripts/output_analysis/visualizations/plot_run_roc_v2.py \\
+      --run-dir analysis_data/v2_ensemble_fused_binary_unet_40ep_5fold_v1 \\
+      --eval-mode binary --tower-mode ensemble --probs probs_fused_head
+"""
+
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.metrics import roc_curve, auc as sk_auc
+
+
+# ---------------------------------------------------------------------------
+# Probs auto-detection
+# ---------------------------------------------------------------------------
+
+_PROBS_PRIORITY: dict[str, list[str]] = {
+    "ensemble":  ["probs_fused_head", "probs_fused"],
+    "bilateral": ["probs_bilat"],
+    "single":    ["probs_classic"],
+    "classic":   ["probs_classic"],
+}
+_ALL_PROBS = ["probs_fused_head", "probs_fused", "probs_bilat", "probs_classic"]
+
+
+def detect_probs_stem(fold_dir: Path, tower_mode: str | None) -> str | None:
+    priority = _PROBS_PRIORITY.get(tower_mode, _ALL_PROBS) if tower_mode else _ALL_PROBS
+    for stem in priority:
+        if (fold_dir / f"{stem}.npy").exists():
+            return stem
+    return None
+
+
+# ---------------------------------------------------------------------------
+# Data loading
+# ---------------------------------------------------------------------------
+
+def find_fold_dirs(mode_dir: Path) -> list[Path]:
+    return sorted(
+        [p for p in mode_dir.iterdir() if p.is_dir() and p.name.startswith("fold")],
+        key=lambda p: int(p.name.replace("fold", "")),
+    )
+
+
+def _find_y_true(fold_dir: Path) -> Path | None:
+    """
+    Return path to y_true.npy for this fold.  If missing from fold_dir
+    (can happen for bilateral-only old runs), fall back to the same fold
+    index under sibling tower-mode directories (ensemble → single → classic).
+    Labels are shared across tower modes within the same fold.
+    """
+    local = fold_dir / "y_true.npy"
+    if local.exists():
+        return local
+    fold_name = fold_dir.name           # e.g. "fold0"
+    tower_dir = fold_dir.parent         # e.g. .../binary/bilateral
+    eval_dir  = tower_dir.parent        # e.g. .../binary
+    for fallback in ("ensemble", "single", "classic"):
+        candidate = eval_dir / fallback / fold_name / "y_true.npy"
+        if candidate.exists():
+            return candidate
+    return None
+
+
+def load_fold(
+    fold_dir: Path,
+    probs_stem: str,
+    eval_mode: str,
+) -> tuple[np.ndarray, np.ndarray] | None:
+    y_path = _find_y_true(fold_dir)
+    p_path = fold_dir / f"{probs_stem}.npy"
+    if y_path is None or not p_path.exists():
+        return None
+    y = np.load(y_path)
+    p = np.load(p_path)
+    if eval_mode == "binary":
+        mask = np.isin(y, [0, 1])
+        y, p = y[mask], p[mask]
+        if p.shape[1] > 2:
+            p = p[:, :2]
+    return y, p
+
+
+# ---------------------------------------------------------------------------
+# ROC helpers
+# ---------------------------------------------------------------------------
+
+def per_class_roc(y: np.ndarray, p: np.ndarray) -> dict[int, dict]:
+    """OVR ROC for each class. Returns {k: {fpr, tpr, auc}}."""
+    out: dict[int, dict] = {}
+    for k in range(p.shape[1]):
+        yb = (y == k).astype(np.uint8)
+        if yb.sum() == 0 or yb.sum() == len(yb):
+            continue
+        fpr, tpr, _ = roc_curve(yb, p[:, k])
+        out[k] = {"fpr": fpr, "tpr": tpr, "auc": sk_auc(fpr, tpr)}
+    return out
+
+
+# ---------------------------------------------------------------------------
+# Plotting
+# ---------------------------------------------------------------------------
+
+def plot_perfold(
+    per_fold: list[tuple[int, dict]],
+    out_dir: Path,
+    class_names: list[str],
+    probs_stem: str,
+    eval_mode: str,
+) -> None:
+    """One figure per class: each fold as a separate line."""
+    all_classes = sorted({k for _, curves in per_fold for k in curves})
+    if eval_mode == "binary":
+        all_classes = [k for k in all_classes if k == 1]
+    out_dir.mkdir(parents=True, exist_ok=True)
+    for k in all_classes:
+        cname = class_names[k] if k < len(class_names) else f"class_{k}"
+        fig, ax = plt.subplots(figsize=(9, 7))
+        ax.plot([0, 1], [0, 1], linestyle="--", linewidth=1, color="grey")
+        for fold_idx, curves in per_fold:
+            if k not in curves:
+                continue
+            c = curves[k]
+            auc_val = c["auc"]
+            ax.plot(c["fpr"], c["tpr"], linewidth=1.5,
+                    label=f"Fold {fold_idx}  (AUC {auc_val:.3f})")
+        ax.set_xlabel("False Positive Rate")
+        ax.set_ylabel("True Positive Rate")
+        ax.set_title(f"Per-fold ROC — {cname}  [{probs_stem}]")
+        ax.legend(loc="lower right")
+        fig.tight_layout()
+        safe = cname.replace(" ", "_")
+        fig.savefig(out_dir / f"roc_{probs_stem}_{safe}_perfold.png", dpi=160)
+        plt.close(fig)
+        print(f"  Saved per-fold ROC ({cname})")
+
+
+def plot_mean_ovr(
+    per_fold: list[tuple[int, dict]],
+    out_dir: Path,
+    class_names: list[str],
+    probs_stem: str,
+    eval_mode: str,
+) -> None:
+    """Mean ± SD OVR ROC — all classes on one figure."""
+    all_classes = sorted({k for _, curves in per_fold for k in curves})
+    if eval_mode == "binary":
+        all_classes = [k for k in all_classes if k == 1]
+    grid = np.linspace(0, 1, 501)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    fig, ax = plt.subplots(figsize=(9, 7))
+    ax.plot([0, 1], [0, 1], linestyle="--", linewidth=1, color="grey")
+    for k in all_classes:
+        cname = class_names[k] if k < len(class_names) else f"class_{k}"
+        tprs, aucs = [], []
+        for _, curves in per_fold:
+            if k not in curves:
+                continue
+            c = curves[k]
+            tprs.append(np.interp(grid, c["fpr"], c["tpr"]))
+            aucs.append(c["auc"])
+        if not tprs:
+            continue
+        tprs_arr = np.vstack(tprs)
+        mean = tprs_arr.mean(axis=0)
+        std  = tprs_arr.std(axis=0)
+        label = f"{cname}  (AUC {np.nanmean(aucs):.3f} ± {np.nanstd(aucs):.3f})"
+        line, = ax.plot(grid, mean, linewidth=2, label=label)
+        ax.fill_between(grid,
+                        np.clip(mean - std, 0, 1),
+                        np.clip(mean + std, 0, 1),
+                        alpha=0.15, color=line.get_color())
+    ax.set_xlabel("False Positive Rate")
+    ax.set_ylabel("True Positive Rate")
+    ax.set_title(f"Mean OVR ROC (± 1 SD)  [{probs_stem}]")
+    ax.legend(loc="lower right")
+    fig.tight_layout()
+    out_path = out_dir / f"roc_{probs_stem}_mean_ovr.png"
+    fig.savefig(out_path, dpi=160)
+    plt.close(fig)
+    print(f"  Saved mean OVR ROC: {out_path}")
+
+
+# ---------------------------------------------------------------------------
+# CLI
+# ---------------------------------------------------------------------------
+
+def main() -> None:
+    ap = argparse.ArgumentParser(
+        description="Per-fold and mean OVR ROC plots for a single v2 run.",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=__doc__,
+    )
+    ap.add_argument("--run-dir", required=True, type=Path,
+                    help="Top-level run directory (e.g. analysis_data/v2_my_run).")
+    ap.add_argument("--eval-mode", required=True, choices=["binary", "multiclass"])
+    ap.add_argument("--tower-mode", required=True,
+                    choices=["single", "classic", "ensemble", "bilateral"],
+                    help="Tower mode subdirectory to read from.")
+    ap.add_argument("--probs", default=None,
+                    help="Probs file stem to use (e.g. probs_fused, probs_bilat, "
+                         "probs_fused_head). Auto-detected if omitted.")
+    ap.add_argument("--class-names", nargs="*",
+                    default=["Healthy", "Glaucoma", "Suspect"])
+    args = ap.parse_args()
+
+    mode_dir = args.run_dir / args.eval_mode / args.tower_mode
+    if not mode_dir.exists():
+        raise SystemExit(f"Directory not found: {mode_dir}")
+
+    fold_dirs = find_fold_dirs(mode_dir)
+    if not fold_dirs:
+        raise SystemExit(f"No fold subdirectories found in {mode_dir}")
+
+    # Determine probs stem
+    probs_stem = args.probs
+    if probs_stem is None:
+        for fd in fold_dirs:
+            probs_stem = detect_probs_stem(fd, args.tower_mode)
+            if probs_stem:
+                break
+    if probs_stem is None:
+        raise SystemExit(f"Could not detect a probs file in {mode_dir}/fold*/")
+    print(f"Using probs: {probs_stem}.npy")
+
+    # Load all folds
+    per_fold: list[tuple[int, dict]] = []
+    for fd in fold_dirs:
+        fold_idx = int(fd.name.replace("fold", ""))
+        result = load_fold(fd, probs_stem, args.eval_mode)
+        if result is None:
+            print(f"  [skip] fold {fold_idx}: missing y_true or {probs_stem}.npy")
+            continue
+        y, p = result
+        curves = per_class_roc(y, p)
+        per_fold.append((fold_idx, curves))
+        auc_str = "  ".join(
+            f"class{k}={v['auc']:.3f}" for k, v in curves.items()
+        )
+        print(f"  fold {fold_idx}: {auc_str}")
+
+    if not per_fold:
+        raise SystemExit("No usable folds — nothing to plot.")
+
+    out_dir = mode_dir / "plots"
+    plot_perfold(per_fold, out_dir, args.class_names, probs_stem, args.eval_mode)
+    plot_mean_ovr(per_fold, out_dir, args.class_names, probs_stem, args.eval_mode)
+    print(f"\nPlots written to {out_dir}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/output_analysis/visualizations/rebuild_run_best_plots.py b/scripts/output_analysis/visualizations/rebuild_run_best_plots.py
index c81bbe9..7f7d693 100755
--- a/scripts/output_analysis/visualizations/rebuild_run_best_plots.py
+++ b/scripts/output_analysis/visualizations/rebuild_run_best_plots.py
@@ -19,7 +19,7 @@ from types import SimpleNamespace
 import matplotlib
 import sys
 
-REPO_ROOT = Path(__file__).resolve().parents[1]
+REPO_ROOT = Path(__file__).resolve().parents[3]
 if str(REPO_ROOT) not in sys.path:
     sys.path.insert(0, str(REPO_ROOT))
 matplotlib.use("Agg")
@@ -70,10 +70,34 @@ def load_summary(run_dir: Path) -> dict:
         return json.load(fh)
 
 
+def resolve_data_dir(raw_dir: str) -> str:
+    """
+    Resolve dataset paths saved in legacy cli_args.json.
+
+    Older runs often store "ClinicalData"/"FundusImages" relative to a
+    dataset root, while current repo layout uses "Papila/<dir>".
+    """
+    p = Path(raw_dir)
+    if p.exists():
+        return str(p)
+
+    candidates = [
+        REPO_ROOT / p,
+        REPO_ROOT / "Papila" / p,
+    ]
+    for cand in candidates:
+        if cand.exists():
+            return str(cand)
+
+    return str(p)
+
+
 def prepare_clinical(cli_args: dict, run_dir: Path) -> tuple:
+    image_dir = resolve_data_dir(cli_args["image_dir"])
+    clinical_dir = resolve_data_dir(cli_args["clinical_dir"])
     clinical = build_papila_clinical(
-        cli_args["image_dir"],
-        cli_args["clinical_dir"],
+        image_dir,
+        clinical_dir,
         cli_args["label_col"],
         cli_args["cat_cols"],
         n_splits=cli_args["n_splits"],
@@ -103,10 +127,12 @@ def prepare_clinical(cli_args: dict, run_dir: Path) -> tuple:
 
 
 def build_ht_args(cli_args: dict, fold: int, run_dir: Path, models_dir: Path, holdout_df):
+    image_dir = resolve_data_dir(cli_args["image_dir"])
+    clinical_dir = resolve_data_dir(cli_args["clinical_dir"])
     # Copy of the training-time namespace so HyperTower can be re-instantiated.
     return SimpleNamespace(
-        image_dir=cli_args["image_dir"],
-        clinical_dir=cli_args["clinical_dir"],
+        image_dir=image_dir,
+        clinical_dir=clinical_dir,
         label_col=cli_args["label_col"],
         cat_cols=cli_args["cat_cols"],
         batch_size=cli_args["batch_size"],
@@ -253,6 +279,14 @@ def choose_head_probs(head: str, probs_f, probs_i, probs_m):
     return probs_i
 
 
+def legacy_probs_suffix(head: str) -> str:
+    if head == "image":
+        return "img"
+    if head == "metadata":
+        return "md"
+    return "fused"
+
+
 def ensure_binary_slice(y_true, *arrays):
     mask = np.isin(y_true, [0, 1])
     filtered = [y_true[mask]]
@@ -264,8 +298,17 @@ def ensure_binary_slice(y_true, *arrays):
     return filtered
 
 
-def plot_overlays(per_fold_curves, out_dir: Path, class_names: list[str], head: str, suffix: str = ""):
+def plot_overlays(
+    per_fold_curves,
+    out_dir: Path,
+    class_names: list[str],
+    head: str,
+    eval_mode: str,
+    suffix: str = "",
+):
     keys = sorted({k for _, curves in per_fold_curves for k in curves.keys()})
+    if eval_mode == "binary":
+        keys = [k for k in keys if k == 1]
     if not keys:
         return
     name_map = {k: (class_names[k] if k < len(class_names) else f"class_{k}") for k in keys}
@@ -293,8 +336,17 @@ def plot_overlays(per_fold_curves, out_dir: Path, class_names: list[str], head:
         plt.close(fig)
 
 
-def plot_mean_sd(per_fold_curves, out_dir: Path, class_names: list[str], head: str, suffix: str = ""):
+def plot_mean_sd(
+    per_fold_curves,
+    out_dir: Path,
+    class_names: list[str],
+    head: str,
+    eval_mode: str,
+    suffix: str = "",
+):
     keys = sorted({k for _, curves in per_fold_curves for k in curves.keys()})
+    if eval_mode == "binary":
+        keys = [k for k in keys if k == 1]
     if not keys:
         return
     grid = np.linspace(0, 1, 501)
@@ -363,10 +415,43 @@ def main():
             print(f"[skip] Fold {fold_idx}: no best_epoch recorded.")
             continue
 
+        # Prefer saved fold arrays when available. This avoids reconstructing
+        # HyperTower for legacy runs whose external weight paths no longer exist.
+        base = run_dir / f"fold{fold_idx}{file_suffix}"
+        y_path = Path(f"{base}_y_true.npy")
+        p_path = Path(f"{base}_probs_{legacy_probs_suffix(head)}.npy")
+        if y_path.exists() and p_path.exists():
+            y_true = np.load(y_path)
+            head_probs = np.load(p_path)
+            if cli_args["eval_mode"] == "binary" and head_probs.shape[1] >= 2:
+                head_probs = head_probs[:, :2]
+            curves = compute_per_class_curves(y_true, head_probs)
+            per_fold_curves.append((fold_idx, curves))
+            try:
+                if head_probs.shape[1] > 2:
+                    fold_auc = roc_auc_score(y_true, head_probs, multi_class="ovr", average="macro")
+                else:
+                    target_scores = head_probs[:, 1] if head_probs.shape[1] > 1 else head_probs[:, 0]
+                    fold_auc = roc_auc_score(y_true, target_scores)
+                fold_aucs.append(fold_auc)
+                print(
+                    f"[info] Fold {fold_idx}: using saved arrays "
+                    f"({y_path.name}, {p_path.name}), AUC={fold_auc:.4f}"
+                )
+            except Exception:
+                print(
+                    f"[warning] Fold {fold_idx}: using saved arrays "
+                    f"({y_path.name}, {p_path.name}) but AUC failed."
+                )
+            continue
+
         fold_models_dir = base_models_dir / f"fold{fold_idx}"
         best_checkpoint = fold_models_dir / "model_best.pt"
         if not best_checkpoint.exists():
-            print(f"[warning] Fold {fold_idx}: missing model_best.pt at {best_checkpoint}")
+            print(
+                f"[warning] Fold {fold_idx}: missing model_best.pt at {best_checkpoint} "
+                f"and missing fallback arrays {y_path.name}/{p_path.name}"
+            )
             continue
 
         ht_args = build_ht_args(cli_args, fold_idx, run_dir, fold_models_dir, holdout_df)
@@ -435,8 +520,8 @@ def main():
         raise SystemExit("No folds processed; nothing to plot.")
 
     plots_dir = run_dir / "plots"
-    plot_overlays(per_fold_curves, plots_dir, class_names, head, file_suffix)
-    plot_mean_sd(per_fold_curves, plots_dir, class_names, head, file_suffix)
+    plot_overlays(per_fold_curves, plots_dir, class_names, head, cli_args["eval_mode"], file_suffix)
+    plot_mean_sd(per_fold_curves, plots_dir, class_names, head, cli_args["eval_mode"], file_suffix)
 
     if fold_aucs:
         print(f"[info] {head} head mean AUC across folds: {np.mean(fold_aucs):.4f} ± {np.std(fold_aucs):.4f}")