diff --git a/classes/v2/backbones.py b/classes/v2/backbones.py index 5a92109..51b849d 100755 --- a/classes/v2/backbones.py +++ b/classes/v2/backbones.py @@ -15,7 +15,7 @@ class BackboneSpec: strip: Callable[[nn.Module], tuple] # fn(model)->(out_dim, model_no_head) blocks: Callable[[nn.Module], List[nn.Module]] # fn(model)->ordered blocks for freezing -REFUGELIKE_BACKBONE_PATH = Path("models/refuge/classifier/refugelike_backbone.pt") +REFUGELIKE_BACKBONE_PATH = Path("models/v2/refuge/refugelike_backbone.pt") REFUGE_DENSENET_PATH = Path("models/refuge/classifier/refuge_densenet_backbone.pt") REFUGE_EFFICIENT_B0_PATH = Path("models/refuge/classifier/refuge_efficient_b0_backbone.pt") REFUGE_EFFICIENT_B7_PATH = Path("models/refuge/classifier/refuge_efficient_b7_backbone.pt") diff --git a/classes/v2/croppers.py b/classes/v2/croppers.py index bb9f96b..dc49387 100644 --- a/classes/v2/croppers.py +++ b/classes/v2/croppers.py @@ -82,7 +82,9 @@ class UNetImageCropper: def _infer_masks(self, image: Image.Image) -> Optional[Tuple[np.ndarray, np.ndarray]]: resized = self.segmenter.preprocess_image(image) - tensor = self.to_tensor(resized).unsqueeze(0).to(self.segmenter.device) + tensor = self.segmenter._normalize_tensor( + self.to_tensor(resized).to(self.segmenter.device) + ).unsqueeze(0) with torch.no_grad(): logits = self.segmenter.model(tensor) diff --git a/classes/v2/dataset.py b/classes/v2/dataset.py index 76a06ad..c8c2f1e 100644 --- a/classes/v2/dataset.py +++ b/classes/v2/dataset.py @@ -16,6 +16,7 @@ class ClinicalDataset(Dataset): image_preprocessor=None, geometry_provider=None, geometry_dim: int = 0, + image_cache: "dict | None" = None, ): self.clinical = clinical_data self.transform_image = img_transform @@ -23,6 +24,7 @@ class ClinicalDataset(Dataset): self.image_preprocessor = image_preprocessor self.geometry_provider = geometry_provider self.geometry_dim = geometry_dim if geometry_provider is not None else 0 + self.image_cache = image_cache def __len__(self): return len(self.clinical.df) @@ -31,7 +33,13 @@ class ClinicalDataset(Dataset): row = self.clinical.df.iloc[idx] # load & transform image img_path = self.clinical.get_image_path(row) - orig_img = Image.open(img_path).convert("RGB") + cache_key = str(img_path) + if self.image_cache is not None and cache_key in self.image_cache: + orig_img = Image.fromarray(self.image_cache[cache_key]) + else: + orig_img = Image.open(img_path).convert("RGB") + if self.image_cache is not None: + self.image_cache[cache_key] = np.asarray(orig_img, dtype=np.uint8) img = orig_img if self.image_preprocessor is not None: img = self.image_preprocessor(img, img_path) diff --git a/classes/v2/metrics.py b/classes/v2/metrics.py index 73f5382..f65e091 100644 --- a/classes/v2/metrics.py +++ b/classes/v2/metrics.py @@ -12,6 +12,7 @@ from sklearn.metrics import ( matthews_corrcoef, recall_score, roc_auc_score, + roc_curve, ) @@ -142,26 +143,37 @@ def compute_extended_metrics( # --------------------------------------------------------------------------- def tune_binary_threshold(y_true: np.ndarray, p1: np.ndarray) -> float: - if y_true.size == 0: + """Pick threshold via Youden's J (sensitivity + specificity − 1). + + This is class-distribution independent, unlike maximising raw accuracy, + which is biased toward the majority class on imbalanced validation sets. + Falls back to 0.5 if both classes are not present. + """ + if y_true.size == 0 or len(np.unique(y_true)) < 2: return 0.5 - grid = np.linspace(0.0, 1.0, 1001) - best_t, best_acc = 0.5, -1.0 - for t in grid: - pred = (p1 >= t).astype(int) - acc = float((pred == y_true).mean()) - if acc > best_acc or (acc == best_acc and abs(t - 0.5) < abs(best_t - 0.5)): - best_acc, best_t = acc, float(t) - return best_t + fpr, tpr, thresholds = roc_curve(y_true, p1) + j = tpr + (1.0 - fpr) - 1.0 + return float(thresholds[np.argmax(j)]) def multiclass_acc_with_bias(y_true: np.ndarray, probs: np.ndarray, bias: np.ndarray) -> float: + """Balanced accuracy (mean per-class recall) after applying log-space bias.""" if y_true.size == 0: return float("nan") logits = np.log(np.clip(probs, 1e-8, 1.0)) + bias.reshape(1, -1) - return float((np.argmax(logits, axis=1) == y_true).mean()) + preds = np.argmax(logits, axis=1) + classes = np.unique(y_true) + per_class = [(preds[y_true == c] == c).mean() for c in classes] + return float(np.mean(per_class)) def tune_multiclass_bias(y_true: np.ndarray, probs: np.ndarray, *, iters: int = 2) -> np.ndarray: + """Grid-search per-class log-space bias to maximise balanced accuracy. + + Balanced accuracy (mean per-class recall) is class-distribution independent, + unlike raw accuracy which is biased toward the majority class on imbalanced + validation sets. + """ if y_true.size == 0 or probs.size == 0: return np.zeros((0,), dtype=float) c = probs.shape[1] diff --git a/classes/v2/models.py b/classes/v2/models.py index d5f762a..06755be 100644 --- a/classes/v2/models.py +++ b/classes/v2/models.py @@ -283,6 +283,7 @@ def train_single_epoch( *, phase: str, bcd_prob: float = 0.5, + tower_loss_mode: str = "bcd", ) -> tuple[float, float]: model.train() _set_single_phase(model, phase) @@ -337,6 +338,11 @@ def train_single_epoch( elif bridge_mode == "image_only": logits = model.bridge.classifier_img(img_feats) loss = F.cross_entropy(logits, y) + elif tower_loss_mode == "all": + loss_i = F.cross_entropy(model.bridge.classifier_img(img_feats), y) + loss_m = F.cross_entropy(model.bridge.classifier_md(md_feats), y) + logits, _, _ = model.bridge(img_feats, md_feats) + loss = F.cross_entropy(logits, y) + loss_i + loss_m elif random() < bcd_prob: if random() < 0.5: logits = model.bridge.classifier_img(img_feats) @@ -368,6 +374,7 @@ def train_bilateral_epoch( *, phase: str, bcd_prob: float = 0.5, + tower_loss_mode: str = "bcd", ) -> tuple[float, float]: model.train() _set_bilateral_phase(model, phase) @@ -394,7 +401,12 @@ def train_bilateral_epoch( logits, _, _ = model.bridge(joint_img, joint_md) loss = F.cross_entropy(logits, y) else: - if random() < bcd_prob: + if tower_loss_mode == "all": + loss_i = F.cross_entropy(model.aux_img(joint_img), y) + loss_m = F.cross_entropy(model.aux_md(joint_md), y) + logits, _, _ = model.bridge(joint_img, joint_md) + loss = F.cross_entropy(logits, y) + loss_i + loss_m + elif random() < bcd_prob: if random() < 0.5: logits = model.aux_img(joint_img) else: diff --git a/classes/v2/papila_builders.py b/classes/v2/papila_builders.py index b831aec..c462b79 100644 --- a/classes/v2/papila_builders.py +++ b/classes/v2/papila_builders.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import Dict, List +from typing import Callable, Dict, List, Optional import numpy as np import pandas as pd @@ -33,21 +33,80 @@ def _nearest_pachy_key(x: float) -> int: return int(_PACHY_KEYS[idx]) -# Ratio derived from patients with both Pneumatic and Perkins readings (n=41, OD+OS combined). -# Pneumatic / Perkins mean ratio = 1.158; applied to Perkins-only rows to put them on the -# Pneumatic scale before IOP_corr is computed. -_PERKINS_TO_PNEUMATIC_RATIO: float = 1.158 +def _fit_perkins_converter( + frames: List[pd.DataFrame], method: str +) -> Callable[[float, Optional[float]], float]: + """ + Fit a Perkins→Pneumatic converter from pooled paired observations across all frames. + Returns a callable: converter(perkins_value, pachymetry_value) -> float. + Supported methods: "ratio", "ols", "lad", "multi". + """ + combined = pd.concat(frames, ignore_index=True) + paired = combined.dropna(subset=["Pneumatic", "Perkins"]) + pneumatic = paired["Pneumatic"].values.astype(float) + perkins = paired["Perkins"].values.astype(float) + + if len(paired) == 0: + raise ValueError("No paired Pneumatic+Perkins observations found; cannot fit converter.") + + if method == "ratio": + ratio = float((pneumatic / perkins).mean()) + def converter_ratio(p: float, pachy: Optional[float] = None) -> float: + return p * ratio + return converter_ratio + + elif method == "ols": + from scipy import stats as _stats + slope, intercept, *_ = _stats.linregress(perkins, pneumatic) + slope, intercept = float(slope), float(intercept) + def converter_ols(p: float, pachy: Optional[float] = None) -> float: + return p * slope + intercept + return converter_ols + + elif method == "lad": + from scipy import stats as _stats + from scipy.optimize import minimize as _minimize + slope0, intercept0, *_ = _stats.linregress(perkins, pneumatic) + def _lad_loss(params): + a, b = params + return np.abs(pneumatic - (a * perkins + b)).mean() + res = _minimize(_lad_loss, x0=[slope0, intercept0], method="Nelder-Mead") + slope, intercept = float(res.x[0]), float(res.x[1]) + def converter_lad(p: float, pachy: Optional[float] = None) -> float: + return p * slope + intercept + return converter_lad + + elif method == "multi": + from numpy.linalg import lstsq as _lstsq + paired_multi = combined.dropna(subset=["Pneumatic", "Perkins", "Pachymetry"]) + if len(paired_multi) == 0: + raise ValueError("No paired Pneumatic+Perkins+Pachymetry rows; cannot fit multi method.") + pneu = paired_multi["Pneumatic"].values.astype(float) + perk = paired_multi["Perkins"].values.astype(float) + pachy_vals = paired_multi["Pachymetry"].values.astype(float) + X = np.column_stack([perk, pachy_vals, np.ones(len(perk))]) + coeffs, *_ = _lstsq(X, pneu, rcond=None) + slope, pachy_coef, intercept = float(coeffs[0]), float(coeffs[1]), float(coeffs[2]) + pachy_fallback = float(pachy_vals.mean()) + def converter_multi(p: float, pachy: Optional[float] = None) -> float: + pv = pachy if (pachy is not None and not np.isnan(pachy)) else pachy_fallback + return p * slope + pachy_coef * pv + intercept + return converter_multi + + else: + raise ValueError(f"Unknown iop_corr_method: {method!r}. Choose ratio/ols/lad/multi.") -def _pick_iop(row: pd.Series) -> float: - """Prefer Pneumatic; scale Perkins to Pneumatic scale if Pneumatic is absent.""" +def _pick_iop(row: pd.Series, converter: Callable) -> float: + """Prefer Pneumatic; convert Perkins to Pneumatic scale if Pneumatic is absent.""" pneumatic = row.get("Pneumatic", np.nan) if not pd.isna(pneumatic): return float(pneumatic) perkins = row.get("Perkins", np.nan) - if not pd.isna(perkins): - return float(perkins) * _PERKINS_TO_PNEUMATIC_RATIO - return np.nan + if pd.isna(perkins): + return np.nan + pachy = row.get("Pachymetry", np.nan) + return converter(float(perkins), None if pd.isna(pachy) else float(pachy)) def _correct_iop(raw_iop: float, pachy: float) -> float: @@ -60,14 +119,20 @@ def _correct_iop(raw_iop: float, pachy: float) -> float: return float(raw_iop) + float(_PACHY_TABLE[key]) -def _apply_iop_and_drop_md(df: pd.DataFrame) -> pd.DataFrame: +def _apply_iop_and_drop_md( + df: pd.DataFrame, + converter: Callable, + drop_raw: bool = False, +) -> pd.DataFrame: """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) + df["IOP_raw"] = df.apply(lambda row: _pick_iop(row, converter), 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) ] drop_cols = [c for c in ("Pneumatic", "Perkins", "VF_MD") if c in df.columns] + if drop_raw: + drop_cols.append("IOP_raw") if drop_cols: df.drop(columns=drop_cols, inplace=True) return df @@ -117,6 +182,9 @@ def build_papila_data( cat_cols: List[str], n_splits: int = 5, random_seed: int = 42, + iop_corr_method: str = "ratio", + iop_drop_raw: bool = False, + exclude_cols: Optional[List[str]] = None, ) -> DataBundle: """ Build a DataBundle for PAPILA with dataset-specific preprocessing: @@ -126,12 +194,17 @@ def build_papila_data( - compute IOP_raw / IOP_corr, drop VF_MD - build feature typing & folds """ + _exclude = list(exclude_cols) if exclude_cols else [] + + # Remove excluded cols from cat_cols too so the bundle doesn't try to encode them + effective_cat_cols = [c for c in cat_cols if c not in _exclude] + bundle = DataBundle( image_dir=image_dir, clinical_dir=clinical_dir, label_col=label_col, patient_col="Patient ID", - cat_cols=cat_cols, + cat_cols=effective_cat_cols, n_splits=n_splits, random_seed=random_seed, filename_template="RET{pid:03d}{eye}.jpg", @@ -148,14 +221,17 @@ def build_papila_data( frame["Patient ID"] = frame["Patient ID"].astype(str).str.extract(r"(\d+)")[0].astype(int) _canonicalize_eye_column(frame) - bundle.add_df(od, id_column="ID") - bundle.add_df(os, id_column="ID") + bundle.add_df(od, id_column="ID", exclude_cols=_exclude or None) + bundle.add_df(os, id_column="ID", exclude_cols=_exclude or None) + converter = _fit_perkins_converter(bundle.frames, method=iop_corr_method) for i in range(len(bundle.frames)): - bundle.frames[i] = _apply_iop_and_drop_md(bundle.frames[i]) + bundle.frames[i] = _apply_iop_and_drop_md( + bundle.frames[i], converter=converter, drop_raw=iop_drop_raw + ) - bundle._refresh_master_df() - bundle._infer_or_validate_feature_types() + bundle._refresh_master_df(exclude_cols=_exclude or None) + bundle._infer_or_validate_feature_types(exclude_cols=_exclude or None) bundle._compute_numeric_stats() bundle._build_cat_maps() bundle._compute_feature_dim() diff --git a/classes/v2/v2_hypertower.py b/classes/v2/v2_hypertower.py index 635e565..08ac4b3 100644 --- a/classes/v2/v2_hypertower.py +++ b/classes/v2/v2_hypertower.py @@ -177,6 +177,8 @@ class V2HyperTower: ap.add_argument("--clinical-dir", default="Papila/ClinicalData") ap.add_argument("--label-col", default="Diagnosis") ap.add_argument("--cat-cols", nargs="*", default=["Gender", "Phakic/Pseudophakic"]) + ap.add_argument("--exclude-cols", nargs="*", default=[], + help="Feature columns to exclude entirely from the clinical feature matrix.") ap.add_argument("--eval-mode", choices=["binary", "multiclass"], default="multiclass") ap.add_argument( "--tower-mode", choices=["single", "ensemble", "bilateral", "classic"], @@ -217,6 +219,10 @@ class V2HyperTower: ap.add_argument("--lr", type=float, default=1e-4) ap.add_argument("--bcd-prob", type=float, default=0.5, help="Tower-only step probability during main phase (per model).") + ap.add_argument("--tower-loss-mode", choices=["bcd", "all"], default="bcd", + help="Main-phase tower loss strategy: " + "'bcd' (Block Coordinate Descent — randomly train one tower or fused per step) " + "or 'all' (sum all three losses — fused + img + md — every step).") ap.add_argument("--backbone", default="refugelike") ap.add_argument("--freeze-ratio", type=float, default=0.0) ap.add_argument("--augment", action="store_true") @@ -299,6 +305,21 @@ class V2HyperTower: ap.add_argument("--log-every", type=int, default=1) ap.add_argument("--save-checkpoints", action=argparse.BooleanOptionalAction, default=True, help="Save best_single.pt / best_holdout_single.pt per fold (use --no-save-checkpoints to disable)") + ap.add_argument("--use-last-epoch", action="store_true", default=False, + help="Score using the final epoch's model state rather than the best-AUC checkpoint.") + # IOP feature options + ap.add_argument( + "--iop-corr-method", + choices=["ratio", "ols", "lad", "multi"], + default="ratio", + help="Perkins→Pneumatic conversion method: ratio (default), ols, lad, or multi (+CCT).", + ) + ap.add_argument( + "--iop-drop-raw", + action="store_true", + default=False, + help="Exclude IOP_raw from the feature matrix (keep only IOP_corr).", + ) ap.add_argument( "--fused-head", action="store_true", help="(ensemble mode only) After base SingleEyeHT training, freeze it and train a " @@ -329,6 +350,9 @@ class V2HyperTower: cat_cols=list(args.cat_cols), n_splits=args.n_splits, random_seed=args.fold_seed, + iop_corr_method=getattr(args, "iop_corr_method", "ratio"), + iop_drop_raw=getattr(args, "iop_drop_raw", False), + exclude_cols=list(getattr(args, "exclude_cols", []) or []), ) print(f"Loaded: {len(self.data.df)} rows feature_dim={self.data.feature_dim}", flush=True) self.image_preprocessor = build_image_preprocessor_from_args(args) @@ -488,6 +512,8 @@ class V2HyperTower: 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.y_true_classic is not None: + np.save(fold_dir / "y_true.npy", artifacts.y_true_classic) 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: @@ -1001,6 +1027,7 @@ class V2HyperTower: sl_loss, sl_acc = train_single_epoch( single, _active_loader, opt_single, device, phase=phase_single, bcd_prob=float(args.bcd_prob), + tower_loss_mode=args.tower_loss_mode, ) else: sl_loss, sl_acc = nan, nan @@ -1009,6 +1036,7 @@ class V2HyperTower: bl_loss, bl_acc = train_bilateral_epoch( bilateral, train_bilat_loader, opt_bilateral, device, phase=phase_bilat, bcd_prob=float(args.bcd_prob), + tower_loss_mode=args.tower_loss_mode, ) else: bl_loss, bl_acc = nan, nan @@ -1414,6 +1442,12 @@ class V2HyperTower: fold_logger.close() + # Override: use final epoch state instead of best-AUC checkpoint + if getattr(args, "use_last_epoch", False): + best_single_state = copy.deepcopy(single.state_dict()) + if run_bilat: + best_bilat_state = copy.deepcopy(bilat.state_dict()) + if args.save_checkpoints: if best_single_state is not None: torch.save(best_single_state, fold_dir / "best_single.pt") diff --git a/models/refuge/classifier/densenet/densenetrefuge.pt b/models/refuge/classifier/densenet/densenetrefuge.pt deleted file mode 100755 index 75f1680..0000000 Binary files a/models/refuge/classifier/densenet/densenetrefuge.pt and /dev/null differ diff --git a/models/refuge/classifier/refuge_densenet_backbone.pt b/models/refuge/classifier/refuge_densenet_backbone.pt deleted file mode 100755 index 75f1680..0000000 Binary files a/models/refuge/classifier/refuge_densenet_backbone.pt and /dev/null differ diff --git a/models/v2/refuge/segmentation/per_image_refuge_build/best.pt b/models/v2/refuge/segmentation/per_image_refuge_build/best.pt deleted file mode 100644 index d20edf6..0000000 Binary files a/models/v2/refuge/segmentation/per_image_refuge_build/best.pt and /dev/null differ diff --git a/scripts/basic_analysis/cnn_logits_rf_cv_v2.py b/scripts/basic_analysis/cnn_logits_rf_cv_v2.py new file mode 100644 index 0000000..0f0574a --- /dev/null +++ b/scripts/basic_analysis/cnn_logits_rf_cv_v2.py @@ -0,0 +1,399 @@ +#!/usr/bin/env python3 +"""v2 port of cnn_logits_rf_cv: CNN logit extraction + RF CV using the v2 PAPILA stack. + +Replaces the hardcoded-config v1 version. Data loading, fold splitting, and +feature preparation all go through the v2 stack so that --exclude-cols, +--iop-corr-method, etc. are first-class options. + +Usage: + python scripts/basic_analysis/cnn_logits_rf_cv_v2.py \ + --eval-mode binary --backbone resnet50 --epochs 40 \ + --exclude-cols Phakic/Pseudophakic Axial_Length \ + --run-name cnn_rf_nocrop_binary +""" +from __future__ import annotations + +import argparse +import random +import time +from concurrent.futures import ThreadPoolExecutor +from pathlib import Path +import sys +from types import SimpleNamespace +from typing import Dict, List, Tuple + +import numpy as np +import pandas as pd +from PIL import Image +import torch +from torch import nn +from torch.utils.data import DataLoader +from sklearn.ensemble import RandomForestClassifier +from sklearn.metrics import accuracy_score, roc_auc_score + +REPO_ROOT = Path(__file__).resolve().parents[2] +if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + +from classes.v2.papila_builders import build_papila_data +from classes.v2.backbones import BACKBONES, load_backbone_weights +from classes.v2.dataset import ClinicalDataset, _ClinicalView +from classes.v2.split_manager import PatientFirstSplitManager +from classes.v2.transforms import build_backbone_transform, build_eval_transform + + +# --------------------------------------------------------------------------- +# Model +# --------------------------------------------------------------------------- + +class CNNHead(nn.Module): + def __init__(self, backbone_name: str, num_classes: int) -> None: + super().__init__() + spec = BACKBONES[backbone_name] + backbone = spec.ctor(weights=spec.weights_default) + if backbone_name.startswith("refuge"): + load_backbone_weights(backbone_name, backbone) + out_dim, backbone = spec.strip(backbone) + self.backbone = backbone + self.head = nn.Linear(out_dim, num_classes) + + def forward(self, x: torch.Tensor) -> torch.Tensor: + return self.head(self.backbone(x)) + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _auc_score(y_true: np.ndarray, probs: np.ndarray, num_classes: int) -> float: + try: + if num_classes == 2: + return float(roc_auc_score(y_true, probs[:, 1])) + return float(roc_auc_score(y_true, probs, multi_class="ovr", average="macro")) + except Exception: + return float("nan") + + +def _train_cnn( + model: nn.Module, + loader: DataLoader, + args, + fold: int, + device: torch.device, +) -> None: + model.train() + optimizer = torch.optim.Adam(model.parameters(), lr=args.lr, weight_decay=args.weight_decay) + criterion = nn.CrossEntropyLoss() + for epoch in range(args.epochs): + running_loss = total = correct = 0 + for x_img, _meta, y in loader: + x_img = x_img.to(device) + y = y.to(device=device, dtype=torch.long) + optimizer.zero_grad() + logits = model(x_img) + loss = criterion(logits, y) + loss.backward() + optimizer.step() + running_loss += float(loss.item()) * int(y.size(0)) + correct += int((logits.argmax(1) == y).sum().item()) + total += int(y.size(0)) + if (epoch + 1) % args.log_every == 0: + print( + f" [fold {fold+1}] epoch {epoch+1}/{args.epochs} " + f"loss={running_loss/max(total,1):.4f} acc={correct/max(total,1):.4f}", + flush=True, + ) + + +def _infer_logits( + model: nn.Module, loader: DataLoader, device: torch.device +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Returns (y_true, logits, metadata_vectors).""" + model.eval() + y_all, logits_all, md_all = [], [], [] + with torch.no_grad(): + for x_img, x_md, y in loader: + logits = model(x_img.to(device)).cpu().numpy() + y_np = y.numpy() if torch.is_tensor(y) else np.asarray(y) + md_np = x_md.numpy() if torch.is_tensor(x_md) else np.asarray(x_md) + y_all.append(y_np) + logits_all.append(logits) + md_all.append(md_np) + return ( + np.concatenate(y_all), + np.concatenate(logits_all), + np.concatenate(md_all), + ) + + +# --------------------------------------------------------------------------- +# CLI +# --------------------------------------------------------------------------- + +def build_parser() -> argparse.ArgumentParser: + ap = argparse.ArgumentParser( + description="CNN logit extraction + RF CV (v2 PAPILA stack)" + ) + 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("--cat-cols", nargs="*", default=["Gender", "Phakic/Pseudophakic"]) + ap.add_argument("--exclude-cols", nargs="*", default=[], + help="Feature columns to drop from the clinical feature matrix.") + ap.add_argument("--eval-mode", choices=["binary", "multiclass"], default="binary") + ap.add_argument("--n-splits", type=int, default=5) + ap.add_argument("--fold-seed", type=int, default=42) + ap.add_argument("--holdout-per-class", type=int, default=6, + help="Patients per class reserved for holdout (0 disables).") + ap.add_argument("--holdout-seed", type=int, default=123) + ap.add_argument("--backbone", default="resnet50") + ap.add_argument("--batch-size", type=int, default=8) + ap.add_argument("--epochs", type=int, default=40) + ap.add_argument("--lr", type=float, default=1e-4) + ap.add_argument("--weight-decay", type=float, default=1e-5) + ap.add_argument("--rf-trees", type=int, default=500) + ap.add_argument("--rf-max-depth", type=int, default=None) + ap.add_argument("--rf-min-samples-leaf", type=int, default=1) + ap.add_argument("--num-workers", type=int, default=0) + ap.add_argument("--device", choices=["auto", "cpu", "cuda"], default="auto") + ap.add_argument("--seed", type=int, default=1234) + ap.add_argument("--log-every", type=int, default=5) + ap.add_argument("--run-name", default=None) + ap.add_argument("--output-root", default="analysis_data/basic_analysis/cnn_logits_rf_cv_v2") + ap.add_argument("--iop-corr-method", choices=["ratio", "ols", "lad", "multi"], default="ratio") + ap.add_argument("--iop-drop-raw", action="store_true", default=False) + ap.add_argument("--in-memory-cache", action="store_true", default=True, + help="Cache all images in RAM before training (default: on).") + ap.add_argument("--no-in-memory-cache", action="store_false", dest="in_memory_cache", + help="Disable in-memory image cache.") + ap.add_argument("--cache-workers", type=int, default=4, + help="Threads for prebuilding image cache (default: 4).") + return ap + + +def _prebuild_image_cache(df: "pd.DataFrame", data, n_workers: int, + resize: int = 256) -> dict: + """Load, convert to RGB, resize to `resize`px short edge, and cache as uint8 arrays. + + Storing pre-resized images means the per-batch transform only has to do + CenterCrop + augmentation + ToTensor + Normalize on a small image rather + than resizing a full-resolution fundus image every step. + """ + paths = list({str(data.get_image_path(row)) for _, row in df.iterrows()}) + cache: dict = {} + print(f"[cache] Preloading {len(paths)} images (resize={resize}px) " + f"with {n_workers} threads...", flush=True) + + def _load(p: str): + img = Image.open(p) + if img.mode != "RGB": + img = img.convert("RGB") + w, h = img.size + scale = resize / min(w, h) + img = img.resize((round(w * scale), round(h * scale)), Image.BILINEAR) + return p, np.asarray(img, dtype=np.uint8) + + with ThreadPoolExecutor(max_workers=max(1, n_workers)) as ex: + for path, arr in ex.map(_load, paths): + cache[path] = arr + + ex_shape = cache[paths[0]].shape + print(f"[cache] Done — {len(cache)} images in RAM " + f"({ex_shape[1]}×{ex_shape[0]} each).", flush=True) + return cache + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main() -> None: + args = build_parser().parse_args() + + random.seed(args.seed) + np.random.seed(args.seed) + torch.manual_seed(args.seed) + if torch.cuda.is_available(): + torch.cuda.manual_seed_all(args.seed) + + if args.device == "auto": + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + else: + device = torch.device(args.device) + print(f"Device: {device}", flush=True) + + run_name = args.run_name or time.strftime("%Y%m%d_%H%M%S") + out_dir = Path(args.output_root) / run_name + out_dir.mkdir(parents=True, exist_ok=True) + + # ---------- data ---------- + 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, + iop_corr_method=args.iop_corr_method, + iop_drop_raw=args.iop_drop_raw, + exclude_cols=list(args.exclude_cols or []), + ) + print(f"Loaded: {len(data.df)} rows feature_dim={data.feature_dim}", flush=True) + + 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) + num_classes = 2 if args.eval_mode == "binary" else int(df_mode[args.label_col].nunique()) + print( + f"eval_mode={args.eval_mode} num_classes={num_classes} " + f"rows={len(df_mode)} patients={df_mode['Patient ID'].nunique()}", + flush=True, + ) + + # ---------- splits ---------- + split_manager = PatientFirstSplitManager( + patient_col="Patient ID", label_col=args.label_col + ) + 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, + ) + clinical_ns = SimpleNamespace(df=df_mode, label_col=args.label_col) + plans = split_manager.build_plans(clinical=clinical_ns, args=split_args, profile=None) + n_folds = min(args.n_splits, len(plans)) + + if plans and plans[0].holdout is not None: + plans[0].holdout.to_csv(out_dir / "holdout_patients.csv", index=False) + + # ---------- image cache ---------- + image_cache = None + if args.in_memory_cache: + image_cache = _prebuild_image_cache(df_mode, data, args.cache_workers) + + # ---------- transforms ---------- + train_tf = build_backbone_transform(args.backbone, augment=True) + eval_tf = build_eval_transform(args.backbone) + + rows: List[Dict] = [] + holdout_rows: List[Dict] = [] + + for fold, split in enumerate(plans[:n_folds]): + print(f"\n[info] Fold {fold+1}/{n_folds}", flush=True) + random.seed(args.seed + fold * 100) + np.random.seed(args.seed + fold * 100) + torch.manual_seed(args.seed + fold * 100) + + view_train = _ClinicalView(data, split.train) + view_val = _ClinicalView(data, split.val) + + dl_train = DataLoader( + ClinicalDataset(view_train, train_tf, image_cache=image_cache), + batch_size=args.batch_size, shuffle=True, num_workers=args.num_workers, + ) + dl_val = DataLoader( + ClinicalDataset(view_val, eval_tf, image_cache=image_cache), + batch_size=args.batch_size, shuffle=False, num_workers=args.num_workers, + ) + dl_holdout = None + if split.holdout is not None and not split.holdout.empty: + dl_holdout = DataLoader( + ClinicalDataset(_ClinicalView(data, split.holdout), eval_tf, + image_cache=image_cache), + batch_size=args.batch_size, shuffle=False, num_workers=args.num_workers, + ) + + # Train CNN + model = CNNHead(args.backbone, num_classes=num_classes).to(device) + _train_cnn(model, dl_train, args, fold=fold, device=device) + + # Extract logits (re-run train without augmentation for RF features) + dl_train_eval = DataLoader( + ClinicalDataset(view_train, eval_tf, image_cache=image_cache), + batch_size=args.batch_size, shuffle=False, num_workers=args.num_workers, + ) + y_tr, log_tr, md_tr = _infer_logits(model, dl_train_eval, device) + y_va, log_va, md_va = _infer_logits(model, dl_val, device) + + np.save(out_dir / f"fold{fold}_train_logits.npy", log_tr) + np.save(out_dir / f"fold{fold}_val_logits.npy", log_va) + + X_tr = np.concatenate([log_tr, md_tr], axis=1) + X_va = np.concatenate([log_va, md_va], axis=1) + + rf = RandomForestClassifier( + n_estimators=args.rf_trees, + max_depth=args.rf_max_depth, + min_samples_leaf=args.rf_min_samples_leaf, + class_weight="balanced", + random_state=args.fold_seed + fold, + n_jobs=-1, + ) + rf.fit(X_tr, y_tr) + + p_va = rf.predict_proba(X_va) + pred_va = np.argmax(p_va, axis=1) + rows.append({ + "fold": fold, + "val_acc": float(accuracy_score(y_va, pred_va)), + "val_auc": _auc_score(y_va, p_va, num_classes), + "n_val": int(len(y_va)), + }) + + if dl_holdout is not None: + y_ho, log_ho, md_ho = _infer_logits(model, dl_holdout, device) + np.save(out_dir / f"fold{fold}_holdout_logits.npy", log_ho) + X_ho = np.concatenate([log_ho, md_ho], axis=1) + p_ho = rf.predict_proba(X_ho) + pred_ho = np.argmax(p_ho, axis=1) + holdout_rows.append({ + "fold": fold, + "holdout_acc": float(accuracy_score(y_ho, pred_ho)), + "holdout_auc": _auc_score(y_ho, p_ho, num_classes), + "n_holdout": int(len(y_ho)), + }) + print( + f"[info] Fold {fold+1} RF: val_acc={rows[-1]['val_acc']:.4f}" + f" val_auc={rows[-1]['val_auc']:.4f}" + f" | holdout_acc={holdout_rows[-1]['holdout_acc']:.4f}" + f" holdout_auc={holdout_rows[-1]['holdout_auc']:.4f}", + flush=True, + ) + else: + print( + f"[info] Fold {fold+1} RF: val_acc={rows[-1]['val_acc']:.4f}" + f" val_auc={rows[-1]['val_auc']:.4f}", + flush=True, + ) + + # ---------- save + summarise ---------- + fold_df = pd.DataFrame(rows) + fold_df.to_csv(out_dir / "rf_val_metrics.csv", index=False) + print("\nRF validation metrics:") + print(fold_df.to_string(index=False, float_format=lambda x: f"{x:.4f}")) + + if holdout_rows: + ho_df = pd.DataFrame(holdout_rows) + ho_df.to_csv(out_dir / "rf_holdout_metrics.csv", index=False) + print("\nRF holdout metrics:") + print(ho_df.to_string(index=False, float_format=lambda x: f"{x:.4f}")) + print( + f"\nMeans: val_acc={fold_df['val_acc'].mean():.4f}" + f" val_auc={fold_df['val_auc'].mean():.4f}" + f" holdout_acc={ho_df['holdout_acc'].mean():.4f}" + f" holdout_auc={ho_df['holdout_auc'].mean():.4f}" + ) + else: + print( + f"\nMeans: val_acc={fold_df['val_acc'].mean():.4f}" + f" val_auc={fold_df['val_auc'].mean():.4f}" + ) + + print(f"\nSaved outputs to: {out_dir}", flush=True) + + +if __name__ == "__main__": + main() diff --git a/scripts/exploratory/iop_correction_analysis.py b/scripts/exploratory/iop_correction_analysis.py new file mode 100644 index 0000000..6f14845 --- /dev/null +++ b/scripts/exploratory/iop_correction_analysis.py @@ -0,0 +1,269 @@ +#!/usr/bin/env python3 +""" +Compare Perkins→Pneumatic IOP conversion methods: + - Current: fixed ratio (1.158) + - OLS linear regression (minimises MSE) + - LAD linear regression (minimises MAE — robust to outliers) + - Multiple regression: Perkins + Pachymetry (OLS) + +Also reports on the impact of dropping IOP_raw. + +Run from repo root: + python scripts/exploratory/iop_correction_analysis.py +""" +from __future__ import annotations + +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +from scipy import stats +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec + +REPO = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(REPO)) + +CURRENT_RATIO = 1.158 + +# ── load raw clinical data ──────────────────────────────────────────────────── + +def load_raw() -> pd.DataFrame: + dfs = [] + for fname in ("patient_data_od.xlsx", "patient_data_os.xlsx"): + df = pd.read_excel(REPO / "Papila/ClinicalData" / fname, header=1) + eye = "OD" if "od" in fname else "OS" + df["eyeID"] = eye + dfs.append(df) + return pd.concat(dfs, ignore_index=True) + +df = load_raw() +print(f"Total rows: {len(df)}") +print(f"Columns with IOP: {[c for c in df.columns if 'iop' in c.lower() or c in ('Pneumatic','Perkins')]}") + +# ── find paired rows ────────────────────────────────────────────────────────── + +paired = df.dropna(subset=["Pneumatic", "Perkins", "Pachymetry"]).copy() +pneumatic = paired["Pneumatic"].values.astype(float) +perkins = paired["Perkins"].values.astype(float) +pachymetry = paired["Pachymetry"].values.astype(float) + +print(f"\nPaired obs (Pneumatic + Perkins + Pachymetry): n={len(paired)}") +print(f" Pneumatic: mean={pneumatic.mean():.2f} std={pneumatic.std():.2f} " + f"range=[{pneumatic.min():.1f}, {pneumatic.max():.1f}]") +print(f" Perkins: mean={perkins.mean():.2f} std={perkins.std():.2f} " + f"range=[{perkins.min():.1f}, {perkins.max():.1f}]") + +# ── current ratio ───────────────────────────────────────────────────────────── + +ratio_obs = pneumatic / perkins +ratio_mean = ratio_obs.mean() +ratio_pred = perkins * CURRENT_RATIO +ratio_resid = pneumatic - ratio_pred +ratio_mae = np.abs(ratio_resid).mean() +ratio_rmse = np.sqrt((ratio_resid ** 2).mean()) + +print(f"\n── Current ratio approach ────────────────────────────────") +print(f" Observed Pneumatic/Perkins ratio: mean={ratio_mean:.4f} " + f"std={ratio_obs.std():.4f} range=[{ratio_obs.min():.3f}, {ratio_obs.max():.3f}]") +print(f" Hardcoded ratio used: {CURRENT_RATIO}") +print(f" MAE: {ratio_mae:.3f} mmHg") +print(f" RMSE: {ratio_rmse:.3f} mmHg") + +# ── OLS regression ──────────────────────────────────────────────────────────── + +slope, intercept, r, p, se = stats.linregress(perkins, pneumatic) +reg_pred = slope * perkins + intercept +reg_resid = pneumatic - reg_pred +reg_mae = np.abs(reg_resid).mean() +reg_rmse = np.sqrt((reg_resid ** 2).mean()) + +print(f"\n── OLS regression (minimises MSE): Pneumatic = slope * Perkins + intercept ──") +print(f" slope={slope:.4f} intercept={intercept:.4f}") +print(f" R²={r**2:.4f} p={p:.4e}") +print(f" MAE: {reg_mae:.3f} mmHg") +print(f" RMSE: {reg_rmse:.3f} mmHg") +print(f" Improvement over ratio — MAE: {ratio_mae - reg_mae:+.3f} RMSE: {ratio_rmse - reg_rmse:+.3f}") + +# LAD regression (minimises MAE) — more robust to outliers +from scipy.optimize import minimize + +def lad_loss(params): + a, b = params + return np.abs(pneumatic - (a * perkins + b)).mean() + +lad_res = minimize(lad_loss, x0=[slope, intercept], method="Nelder-Mead") +lad_slope, lad_intercept = lad_res.x +lad_pred = lad_slope * perkins + lad_intercept +lad_resid = pneumatic - lad_pred +lad_mae = np.abs(lad_resid).mean() +lad_rmse = np.sqrt((lad_resid ** 2).mean()) + +print(f"\n── LAD regression (minimises MAE — robust to outliers) ──────") +print(f" slope={lad_slope:.4f} intercept={lad_intercept:.4f}") +print(f" MAE: {lad_mae:.3f} mmHg") +print(f" RMSE: {lad_rmse:.3f} mmHg") +print(f" Improvement over ratio — MAE: {ratio_mae - lad_mae:+.3f} RMSE: {ratio_rmse - lad_rmse:+.3f}") + +# ── Multiple regression: Perkins + Pachymetry ───────────────────────────────── +# Pachymetry affects Perkins (applanation) more than Pneumatic, so including +# CCT should absorb some instrument-specific bias. +# Design matrix: [Perkins, Pachymetry, 1] +from numpy.linalg import lstsq + +X_multi = np.column_stack([perkins, pachymetry, np.ones(len(perkins))]) +coeffs, _, _, _ = lstsq(X_multi, pneumatic, rcond=None) +coef_perkins, coef_pachy, coef_intercept = coeffs +multi_pred = X_multi @ coeffs +multi_resid = pneumatic - multi_pred +multi_mae = np.abs(multi_resid).mean() +multi_rmse = np.sqrt((multi_resid ** 2).mean()) + +# R² for the multiple model +ss_res = (multi_resid ** 2).sum() +ss_tot = ((pneumatic - pneumatic.mean()) ** 2).sum() +multi_r2 = 1 - ss_res / ss_tot + +# Partial correlation of Pachymetry with residual after removing Perkins effect +perkins_resid = pneumatic - (slope * perkins + intercept) +pachy_r, pachy_p = stats.pearsonr(pachymetry, perkins_resid) + +print(f"\n── Multiple OLS (Perkins + Pachymetry, n={len(paired)}) ──────────────") +print(f" Pneumatic = {coef_perkins:.4f}×Perkins + {coef_pachy:.5f}×Pachymetry + {coef_intercept:.4f}") +print(f" R²={multi_r2:.4f} (vs simple OLS R²={r**2:.4f})") +print(f" Pachymetry partial corr with OLS residuals: r={pachy_r:.3f} p={pachy_p:.4f}") +print(f" MAE: {multi_mae:.3f} mmHg (vs ratio {ratio_mae:.3f})") +print(f" RMSE: {multi_rmse:.3f} mmHg (vs ratio {ratio_rmse:.3f})") +print(f" Improvement over ratio — MAE: {ratio_mae - multi_mae:+.3f} RMSE: {ratio_rmse - multi_rmse:+.3f}") +print(f" NOTE: n={len(paired)} with 3 parameters — interpret with caution.") + +# How much does the intercept matter at typical IOP values? +typical_iop = np.array([10, 15, 20, 25]) +print(f"\n Comparison at typical Perkins values:") +print(f" {'Perkins':>8} {'Ratio':>10} {'OLS':>10} {'LAD':>10}") +for v in typical_iop: + ratio_v = v * CURRENT_RATIO + ols_v = slope * v + intercept + lad_v = lad_slope * v + lad_intercept + print(f" {v:>8.1f} {ratio_v:>10.2f} {ols_v:>10.2f} {lad_v:>10.2f}") + +# ── IOP_raw coverage ────────────────────────────────────────────────────────── + +pneumatic_only = df["Pneumatic"].notna() & df["Perkins"].isna() +perkins_only = df["Pneumatic"].isna() & df["Perkins"].notna() +both = df["Pneumatic"].notna() & df["Perkins"].notna() +neither = df["Pneumatic"].isna() & df["Perkins"].isna() + +print(f"\n── IOP measurement coverage ──────────────────────────────") +print(f" Pneumatic only: {pneumatic_only.sum()}") +print(f" Perkins only: {perkins_only.sum()}") +print(f" Both: {both.sum()}") +print(f" Neither: {neither.sum()}") +print(f" Rows where IOP_raw == IOP_corr (no pachy correction): ", end="") + +# ── plot ────────────────────────────────────────────────────────────────────── + +fig = plt.figure(figsize=(16, 5), layout="constrained") +gs = gridspec.GridSpec(1, 3, figure=fig) + +x_line = np.linspace(perkins.min() - 1, perkins.max() + 1, 100) + +# ── Row 1: conversion fits ──────────────────────────────────────────────────── + +# 1a. Scatter with all four fits +ax1 = fig.add_subplot(gs[0, :2]) # spans first two columns +ax1.scatter(perkins, pneumatic, alpha=0.6, s=35, color="gray", zorder=3, label="Observed pairs (n=41)") +ax1.plot(x_line, x_line * CURRENT_RATIO, "r--", lw=2, label=f"Ratio ×{CURRENT_RATIO} MAE={ratio_mae:.2f}") +ax1.plot(x_line, slope * x_line + intercept, "b-", lw=2, label=f"OLS (slope={slope:.3f}, int={intercept:.2f}) MAE={reg_mae:.2f}") +ax1.plot(x_line, lad_slope * x_line + lad_intercept, "g-", lw=2, label=f"LAD (slope={lad_slope:.3f}, int={lad_intercept:.2f}) MAE={lad_mae:.2f}") + +# Multi-reg projected at mean CCT +pachy_mean, pachy_sd = pachymetry.mean(), pachymetry.std() +multi_mean_line = coef_perkins * x_line + coef_pachy * pachy_mean + coef_intercept +ax1.plot(x_line, multi_mean_line, color="purple", lw=2, ls="-.", + label=f"Multi (Perkins+CCT) @ mean CCT MAE={multi_mae:.2f}") + +ax1.set_xlabel("Perkins IOP (mmHg)") +ax1.set_ylabel("Pneumatic IOP (mmHg)") +ax1.set_title("Perkins → Pneumatic conversion: all methods") +ax1.legend(fontsize=8) + +# 1b. MAE / RMSE bar chart +ax_bar = fig.add_subplot(gs[0, 2]) # third column +bar_methods = ["Ratio", "OLS", "LAD", "Multi\n(+CCT)"] +maes = [ratio_mae, reg_mae, lad_mae, multi_mae] +rmses = [ratio_rmse, reg_rmse, lad_rmse, multi_rmse] +bar_colors = ["#e05c5c", "#5c7de0", "#5cc97c", "#9b59b6"] +bx = np.arange(4) +w = 0.35 +ax_bar.bar(bx - w/2, maes, w, label="MAE", color=bar_colors) +ax_bar.bar(bx + w/2, rmses, w, label="RMSE", color=bar_colors, alpha=0.5) +ax_bar.set_xticks(bx); ax_bar.set_xticklabels(bar_methods, fontsize=8) +ax_bar.set_ylabel("Error (mmHg)") +ax_bar.set_title("MAE & RMSE comparison") +ax_bar.legend(fontsize=9) +ax_bar.set_ylim(0, max(rmses) * 1.25) +for i, (mae, rmse) in enumerate(zip(maes, rmses)): + ax_bar.text(i - w/2, mae + 0.05, f"{mae:.2f}", ha="center", va="bottom", fontsize=8) + ax_bar.text(i + w/2, rmse + 0.05, f"{rmse:.2f}", ha="center", va="bottom", fontsize=8) + +out = REPO / "analysis_data/iop_correction_comparison.png" +fig.savefig(out, dpi=150) +print(f"\nPlot saved to {out}") + +# ── Figure 2: Pachymetry analysis ───────────────────────────────────────────── + +fig2, axes2 = plt.subplots(1, 3, figsize=(15, 5)) + +pachy_norm = (pachymetry - pachymetry.mean()) / pachymetry.std() + +# 2a. Pachymetry vs Pneumatic–Perkins difference +ax = axes2[0] +diff = pneumatic - perkins +m, b, rr, pp, _ = stats.linregress(pachymetry, diff) +ax.scatter(pachymetry, diff, alpha=0.6, s=35, color="steelblue") +px = np.linspace(pachymetry.min() - 5, pachymetry.max() + 5, 100) +ax.plot(px, m * px + b, "r-", lw=2, label=f"r={rr:.2f} p={pp:.3f}") +ax.axhline(0, color="k", lw=0.7, ls="--") +ax.set_xlabel("Pachymetry (μm)") +ax.set_ylabel("Pneumatic − Perkins (mmHg)") +ax.set_title("Does CCT predict the instrument gap?") +ax.legend(fontsize=9) + +# 2b. Residuals from ratio vs Pachymetry (coloured by size) +ax = axes2[1] +sc = ax.scatter(pachymetry, ratio_resid, c=perkins, cmap="viridis", alpha=0.7, s=35) +plt.colorbar(sc, ax=ax, label="Perkins IOP") +m2, b2, rr2, pp2, _ = stats.linregress(pachymetry, ratio_resid) +ax.plot(px, m2 * px + b2, "r-", lw=2, label=f"r={rr2:.2f} p={pp2:.3f}") +ax.axhline(0, color="k", lw=0.7, ls="--") +ax.set_xlabel("Pachymetry (μm)") +ax.set_ylabel("Ratio residual (mmHg)") +ax.set_title("Ratio residuals vs Pachymetry\n(colour = Perkins IOP)") +ax.legend(fontsize=9) + +# 2c. Multiple regression residuals vs Pachymetry (should be flat if absorbed) +ax = axes2[2] +m3, b3, rr3, pp3, _ = stats.linregress(pachymetry, multi_resid) +ax.scatter(pachymetry, multi_resid, alpha=0.6, s=35, color="darkorange") +ax.plot(px, m3 * px + b3, "r-", lw=2, label=f"r={rr3:.2f} p={pp3:.3f}") +ax.axhline(0, color="k", lw=0.7, ls="--") +ax.set_xlabel("Pachymetry (μm)") +ax.set_ylabel("Multi-reg residual (mmHg)") +ax.set_title("Multi-reg residuals vs Pachymetry\n(flat = Pachymetry effect absorbed)") +ax.legend(fontsize=9) + +# shared y-axis +ylim2 = max(np.abs(diff).max(), np.abs(ratio_resid).max(), np.abs(multi_resid).max()) + 1 +for ax in axes2[1:]: + ax.set_ylim(-ylim2, ylim2) + +fig2.suptitle( + f"Pachymetry as a covariate (n={len(paired)}, CCT mean={pachymetry.mean():.0f}±{pachymetry.std():.0f} μm)", + fontsize=11, +) +fig2.tight_layout() +out2 = REPO / "analysis_data/iop_pachymetry_analysis.png" +fig2.savefig(out2, dpi=150) +print(f"Plot saved to {out2}") diff --git a/scripts/main/pipeline.ipynb b/scripts/main/pipeline.ipynb index 642f963..08bab4e 100644 --- a/scripts/main/pipeline.ipynb +++ b/scripts/main/pipeline.ipynb @@ -56,6 +56,7 @@ { "cell_type": "code", "execution_count": 4, + "id": "6a600aed", "metadata": {}, "outputs": [ { @@ -84,6 +85,7 @@ }, { "cell_type": "markdown", + "id": "d1ea8b19", "metadata": {}, "source": [ "## 1) Build UNet Manifest" @@ -92,6 +94,7 @@ { "cell_type": "code", "execution_count": 5, + "id": "b07b3e69", "metadata": {}, "outputs": [ { @@ -108,21 +111,28 @@ }, { "cell_type": "markdown", + "id": "dgglybgo5wg", + "source": "## 2) Build Refugelike Backbone\n\nThe `refugelike` backbone is a ResNet-50 pre-trained on REFUGE as an optic disc/cup classifier, then stripped of its classification head and used as a frozen or partially-frozen feature extractor in the HyperTower image tower.\n\n**Steps:**\n1. Train the REFUGE classifier (`--train-clf`)\n2. Export its backbone weights to `models/v2/refuge/refugelike_backbone.pt` (`--export-backbone`)\n\nThe classifier checkpoint is saved to `models/v2/refuge/classifier/resnet50/refuge_classifier_best.pt` by default. \nThe exported backbone is what `--backbone refugelike` loads at runtime (see `classes/v2/backbones.py`).", + "metadata": {} + }, + { + "cell_type": "code", + "id": "b7jt033ul4v", + "source": "BACKBONE_PATH = \"models/v2/refuge/refugelike_backbone.pt\"\n\n# Step 1: train the REFUGE classifier (ResNet-50, 30 epochs by default)\n!python3 scripts/main/refuge/refuge_build.py \\\n --train-clf \\\n --manifest manifest.csv \\\n --device cuda\n\n# Step 2: strip the head and export backbone weights\n!python3 scripts/main/refuge/refuge_build.py \\\n --export-backbone {BACKBONE_PATH} \\\n --manifest manifest.csv \\\n --device cuda", "metadata": {}, - "source": [ - "## 2) Train UNet Segmenter (per-image normalization)\n", - "\n", - "Current tuned baseline:\n", - "- `--device cuda`\n", - "- `--batch-size 8`\n", - "- `--loader-workers 14`\n", - "- `--in-memory-cache`\n", - "- `--cache-workers 4`" - ] + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "8df524a2", + "metadata": {}, + "source": "## 2b) Train UNet Segmenter (per-image normalization)\n\nCurrent tuned baseline:\n- `--device cuda`\n- `--batch-size 8`\n- `--loader-workers 14`\n- `--in-memory-cache`\n- `--cache-workers 4`" }, { "cell_type": "code", "execution_count": 7, + "id": "be2b499a", "metadata": {}, "outputs": [ { @@ -1270,6 +1280,7 @@ }, { "cell_type": "markdown", + "id": "a2a3bd04", "metadata": {}, "source": [ "## 5) Run Pipeline Experiments\n", @@ -1283,6 +1294,7 @@ { "cell_type": "code", "execution_count": null, + "id": "e395f268", "metadata": {}, "outputs": [], "source": [ @@ -1301,41 +1313,18 @@ { "cell_type": "code", "execution_count": null, + "id": "125f35c6", "metadata": {}, "outputs": [], - "source": [ - "# 3b) GT crop — expert segmentation masks crop the optic disc region\n", - "!python scripts/main/v2/multirun_hypertower.py \\\n", - " --tower-modes single ensemble \\\n", - " --eval-modes binary multiclass \\\n", - " --epochs 40 --n-splits 5 \\\n", - " --backbone refugelike \\\n", - " --img-crop-gt \\\n", - " --warmup-md-epochs 50 \\\n", - " --fused-head \\\n", - " --run-name pipeline_gt" - ] + "source": "# 3b) GT crop — expert segmentation masks crop the optic disc region\n!python scripts/main/v2/multirun_hypertower.py \\\n --tower-modes single ensemble \\\n --eval-modes binary multiclass \\\n --epochs 40 --n-splits 5 \\\n --backbone refugelike \\\n --img-crop-manifest manifest.csv \\\n --img-crop-gt \\\n --warmup-md-epochs 50 \\\n --fused-head \\\n --run-name pipeline_gt" }, { "cell_type": "code", "execution_count": null, + "id": "83df73ac", "metadata": {}, "outputs": [], - "source": [ - "# 3c) UNet crop — trained segmenter crops the optic disc region\n", - "!python scripts/main/v2/multirun_hypertower.py \\\n", - " --tower-modes single ensemble \\\n", - " --eval-modes binary multiclass \\\n", - " --epochs 40 --n-splits 5 \\\n", - " --backbone refugelike \\\n", - " --img-crop-manifest analysis_data/unet_manifest.csv \\\n", - " --img-crop-weights models/v2/refuge/segmentation/per_image/best.pt \\\n", - " --img-crop-normalize per_image \\\n", - " --img-crop-cache analysis_data/v2_crops_unet_refuge \\\n", - " --warmup-md-epochs 50 \\\n", - " --fused-head \\\n", - " --run-name pipeline_unet" - ] + "source": "# 3c) UNet crop — trained segmenter crops the optic disc region\n!python scripts/main/v2/multirun_hypertower.py \\\n --tower-modes single ensemble \\\n --eval-modes binary multiclass \\\n --epochs 40 --n-splits 5 \\\n --backbone refugelike \\\n --img-crop-manifest manifest.csv \\\n --img-crop-weights models/v2/refuge/segmentation/per_image/best.pt \\\n --img-crop-normalize per_image \\\n --img-crop-cache analysis_data/v2_crops_unet_refuge \\\n --warmup-md-epochs 50 \\\n --fused-head \\\n --run-name pipeline_unet" }, { "cell_type": "markdown", @@ -1353,51 +1342,7 @@ "id": "7yrcfu0bv1w", "metadata": {}, "outputs": [], - "source": [ - "import subprocess\n", - "from pathlib import Path\n", - "\n", - "RUN_DIRS = {\n", - " \"nocrop\": Path(\"analysis_data/pipeline_nocrop\"),\n", - " \"gt\": Path(\"analysis_data/pipeline_gt\"),\n", - " \"unet\": Path(\"analysis_data/pipeline_unet\"),\n", - "}\n", - "EVAL_MODES = [\"binary\", \"multiclass\"]\n", - "TOWER_MODES = [\"single\", \"ensemble\"]\n", - "\n", - "for run_name, run_dir in RUN_DIRS.items():\n", - " for eval_mode in EVAL_MODES:\n", - " for tower_mode in TOWER_MODES:\n", - " mode_dir = run_dir / eval_mode / tower_mode\n", - " if not mode_dir.exists():\n", - " print(f\" skip (not found): {mode_dir}\")\n", - " continue\n", - " print(f\"--- {run_name} / {eval_mode} / {tower_mode} ---\")\n", - "\n", - " # ROC curves\n", - " subprocess.run([\n", - " \"python\", \"scripts/output_analysis/visualizations/plot_run_roc_v2.py\",\n", - " \"--run-dir\", str(run_dir),\n", - " \"--eval-mode\", eval_mode,\n", - " \"--tower-mode\", tower_mode,\n", - " ], check=True)\n", - "\n", - " # Probability strips — binary only\n", - " if eval_mode == \"binary\":\n", - " subprocess.run([\n", - " \"python\", \"scripts/output_analysis/visualizations/plot_prob_strips.py\",\n", - " \"--run-dir\", str(mode_dir),\n", - " \"--head\", \"fused\", \"--style\", \"strips\",\n", - " ], check=True)\n", - "\n", - " # Probability triangle-3D — multiclass only\n", - " if eval_mode == \"multiclass\":\n", - " subprocess.run([\n", - " \"python\", \"scripts/output_analysis/visualizations/plot_prob_strips.py\",\n", - " \"--run-dir\", str(mode_dir),\n", - " \"--head\", \"fused\", \"--style\", \"triangle3d\",\n", - " ], check=True)" - ] + "source": "import subprocess\nfrom pathlib import Path\n\nRUN_DIRS = {\n \"nocrop\": (Path(\"analysis_data/pipeline_nocrop\"), [\"single\", \"ensemble\"]),\n \"gt\": (Path(\"analysis_data/pipeline_gt\"), [\"single\", \"ensemble\"]),\n \"unet\": (Path(\"analysis_data/pipeline_unet\"), [\"single\", \"ensemble\"]),\n \"imgonly_nocrop\":(Path(\"analysis_data/pipeline_imgonly_nocrop\"),[\"single\"]),\n \"imgonly_gt\": (Path(\"analysis_data/pipeline_imgonly_gt\"), [\"single\"]),\n \"imgonly_unet\": (Path(\"analysis_data/pipeline_imgonly_unet\"), [\"single\"]),\n}\nEVAL_MODES = [\"binary\", \"multiclass\"]\n\nfor run_name, (run_dir, tower_modes) in RUN_DIRS.items():\n for eval_mode in EVAL_MODES:\n for tower_mode in tower_modes:\n mode_dir = run_dir / eval_mode / tower_mode\n if not mode_dir.exists():\n print(f\" skip (not found): {mode_dir}\")\n continue\n print(f\"--- {run_name} / {eval_mode} / {tower_mode} ---\")\n\n # ROC curves (auto-detects all available probs stems)\n subprocess.run([\n \"python\", \"scripts/output_analysis/visualizations/plot_run_roc_v2.py\",\n \"--run-dir\", str(run_dir),\n \"--eval-mode\", eval_mode,\n \"--tower-mode\", tower_mode,\n ], check=True)\n\n # Probability strips — binary only (auto-detects all heads)\n if eval_mode == \"binary\":\n subprocess.run([\n \"python\", \"scripts/output_analysis/visualizations/plot_prob_strips.py\",\n \"--run-dir\", str(mode_dir),\n \"--style\", \"strips\",\n ], check=True)\n\n # Probability triangle-3D — multiclass only (auto-detects all heads)\n if eval_mode == \"multiclass\":\n subprocess.run([\n \"python\", \"scripts/output_analysis/visualizations/plot_prob_strips.py\",\n \"--run-dir\", str(mode_dir),\n \"--style\", \"triangle3d\",\n ], check=True)" }, { "cell_type": "markdown", @@ -1458,4 +1403,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/scripts/main/refuge/refuge_build.py b/scripts/main/refuge/refuge_build.py index 3e21d5a..84152b5 100644 --- a/scripts/main/refuge/refuge_build.py +++ b/scripts/main/refuge/refuge_build.py @@ -6,7 +6,7 @@ Usage examples (after activating .venv_refuge): python refuge_build.py --eval --with-ttt The script expects the REFUGE folder and writes checkpoints under -models/refuge/segmentation and models/refuge/classifier. +models/v2/refuge/segmentation and models/v2/refuge/classifier. """ from __future__ import annotations @@ -46,7 +46,7 @@ from classes.papila_builders import build_papila_clinical REFUGE_ROOT = Path("REFUGE") SEG_CKPT = Path("models/refuge/segmentation/refuge_segmentation_best.pt") -CLF_DIR = Path("models/refuge/classifier") +CLF_DIR = Path("models/v2/refuge/classifier") UNET_WEIGHT_CANDIDATES = ( Path("models/v2/refuge/segmentation/per_image/best.pt"), Path("models/v2/refuge/segmentation/best.pt"), @@ -820,7 +820,7 @@ def parse_args() -> argparse.Namespace: "--clf-checkpoint-path", type=Path, default=None, - help="Optional explicit path for the classifier checkpoint (defaults to models/refuge/classifier//refuge_classifier_best.pt)", + help="Optional explicit path for the classifier checkpoint (defaults to models/v2/refuge/classifier//refuge_classifier_best.pt)", ) parser.add_argument( "--eval-datasets", diff --git a/scripts/main/v2/multirun_hypertower.py b/scripts/main/v2/multirun_hypertower.py index 33087d4..5c988c7 100644 --- a/scripts/main/v2/multirun_hypertower.py +++ b/scripts/main/v2/multirun_hypertower.py @@ -7,7 +7,7 @@ from pathlib import Path import sys import argparse -REPO_ROOT = Path(__file__).resolve().parents[2] +REPO_ROOT = Path(__file__).resolve().parents[3] if str(REPO_ROOT) not in sys.path: sys.path.insert(0, str(REPO_ROOT)) diff --git a/scripts/main/v2/run_10x5cv.py b/scripts/main/v2/run_10x5cv.py new file mode 100644 index 0000000..9e27ece --- /dev/null +++ b/scripts/main/v2/run_10x5cv.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python3 +""" +10× repeated 5-fold CV runner for the best hypertower configuration +(nocrop, ensemble mode, both binary and multiclass). + +Each repetition uses a different fold-seed so the 5 folds are split +differently, giving 50 folds per eval-mode total. Holdout composition +is kept identical across repetitions (same --holdout-seed). + +Results land under: + {output-root}/rep{N:02d}/{eval_mode}/ensemble/fold{K}/ + +Usage +----- + python scripts/main/v2/run_10x5cv.py \ + --n-reps 10 \ + --eval-modes binary multiclass \ + --output-root analysis_data/pipeline_10x5 \ + --epochs 40 --fused-head \ + --backbone refugelike + +Any extra flags are forwarded directly to V2HyperTower. +""" +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[3] +if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + +from classes.v2.v2_hypertower import V2HyperTower + +# Base fold seed for rep 0; rep N uses BASE_SEED + N * SEED_STRIDE +_BASE_SEED = 100 +_SEED_STRIDE = 100 + + +def _parse_own(argv=None): + ap = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + add_help=False, + ) + ap.add_argument("--n-reps", type=int, default=10, + help="Number of repetitions (default: 10).") + ap.add_argument("--eval-modes", nargs="+", + choices=["binary", "multiclass"], + default=["binary", "multiclass"]) + ap.add_argument("--output-root", default="analysis_data/pipeline_10x5", + help="Parent directory for all rep sub-runs.") + ap.add_argument("-h", "--help", action="store_true") + return ap.parse_known_args(argv) + + +def main(argv=None): + own, remaining = _parse_own(argv) + + if own.help: + print(__doc__) + base_parser = V2HyperTower.build_parser() + base_parser.print_help() + return + + base_parser = V2HyperTower.build_parser() + output_root = Path(own.output_root) + first_run = True + + for rep in range(own.n_reps): + fold_seed = _BASE_SEED + rep * _SEED_STRIDE + rep_label = f"rep{rep:02d}" + + for eval_mode in own.eval_modes: + tower_mode = "ensemble" + + # Skip if already fully complete + tm_dir = output_root / rep_label / eval_mode / tower_mode + if (tm_dir / "summary.json").exists(): + print(f"[10x5cv] {rep_label} {eval_mode}:{tower_mode} — already done, skipping.") + first_run = False + continue + + cli = list(remaining) + [ + "--eval-mode", eval_mode, + "--tower-mode", tower_mode, + "--fold-seed", str(fold_seed), + "--run-name", rep_label, + "--output-root", str(output_root), + ] + + # Reuse crop cache across runs after the first + if not first_run: + cli.append("--persist-img-crop-cache") + + print(f"\n[10x5cv] Starting {rep_label} {eval_mode}:{tower_mode} " + f"(fold_seed={fold_seed})") + args = base_parser.parse_args(cli) + V2HyperTower(args).run() + first_run = False + + print(f"\n[10x5cv] All done. Results in: {output_root}") + + +if __name__ == "__main__": + main() diff --git a/scripts/main/v2/run_iop_corr_comparison.py b/scripts/main/v2/run_iop_corr_comparison.py new file mode 100644 index 0000000..eac33f6 --- /dev/null +++ b/scripts/main/v2/run_iop_corr_comparison.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python3 +""" +Run single-mode binary + multiclass for each IOP correction method and +collect all results under one output root for easy comparison. + +Output layout: + analysis_data/iop_corr_comparison/ + ratio/binary/single/ ratio/multiclass/single/ + ols/binary/single/ ols/multiclass/single/ + lad/binary/single/ lad/multiclass/single/ + multi/binary/single/ multi/multiclass/single/ + +Usage +----- + python scripts/main/v2/run_iop_corr_comparison.py [V2HyperTower args...] + +Any extra args (backbone, epochs, img-crop-*, etc.) are forwarded to every run. +""" +from __future__ import annotations + +import sys +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[3] +if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + +from classes.v2.v2_hypertower import V2HyperTower + +IOP_METHODS = ["ratio", "ols", "lad", "multi"] +EVAL_MODES = ["binary", "multiclass"] +OUTPUT_ROOT = "analysis_data/iop_corr_comparison" + + +def main() -> None: + base_parser = V2HyperTower.build_parser() + # Consume only the remaining (forwarded) args — iop-corr-method and + # run-name are set by this script; eval-mode and tower-mode likewise. + _, remaining = base_parser.parse_known_args() + + first_run = True + for method in IOP_METHODS: + for eval_mode in EVAL_MODES: + run_name = method # one sub-folder per method + tm_dir = (Path(OUTPUT_ROOT) / run_name / eval_mode / "single") + if (tm_dir / "summary.json").exists(): + print(f"[iop_corr] {method}/{eval_mode}/single — already done, skipping.") + first_run = False + continue + + cli = list(remaining) + [ + "--eval-mode", eval_mode, + "--tower-mode", "single", + "--iop-corr-method", method, + "--output-root", OUTPUT_ROOT, + "--run-name", run_name, + ] + if not first_run: + cli.append("--persist-img-crop-cache") + + print(f"\n[iop_corr] Starting {method}/{eval_mode}/single ...") + args = base_parser.parse_args(cli) + V2HyperTower(args).run() + first_run = False + + # ── summary table ────────────────────────────────────────────────────── + import json + print("\n" + "=" * 60) + print("IOP correction method comparison — single mode") + print("=" * 60) + header = f"{'Method':<8} {'Mode':<12} {'Val AUC':>10} {'Hld AUC':>10}" + print(header) + print("-" * len(header)) + for method in IOP_METHODS: + for eval_mode in EVAL_MODES: + p = Path(OUTPUT_ROOT) / method / eval_mode / "single" / "summary.json" + if not p.exists(): + print(f"{method:<8} {eval_mode:<12} {'missing':>10} {'missing':>10}") + continue + ms = json.loads(p.read_text()).get("mode_summary", {}) + val = ms.get("classic_best_val", {}) + hld = ms.get("classic_holdout", {}) + val_s = f"{val['auc_mean']:.3f}±{val['auc_std']:.3f}" if val.get("auc_mean") else "—" + hld_s = f"{hld['auc_mean']:.3f}±{hld['auc_std']:.3f}" if hld.get("auc_mean") else "—" + print(f"{method:<8} {eval_mode:<12} {val_s:>10} {hld_s:>10}") + print("=" * 60) + + +if __name__ == "__main__": + main() diff --git a/scripts/main/v2/run_md_mlp.py b/scripts/main/v2/run_md_mlp.py index 4e0b123..33977cd 100644 --- a/scripts/main/v2/run_md_mlp.py +++ b/scripts/main/v2/run_md_mlp.py @@ -323,6 +323,7 @@ def main() -> None: best_epoch = 0 best_phase = "" best_state = None + epoch_log_rows = [] print( f"\n[fold {fold_idx+1}/{args.n_splits}] " @@ -350,6 +351,14 @@ def main() -> None: aggregate_patient=aggregate_patient, ) + hld_auc_ep = float("nan") + hld_acc_ep = float("nan") + if holdout_loader is not None: + _, _, hld_auc_ep, hld_acc_ep = _evaluate_single( + model, holdout_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) @@ -357,6 +366,17 @@ def main() -> None: best_epoch = ep + 1 best_phase = phase + epoch_log_rows.append({ + "epoch": ep + 1, + "phase": phase, + "train_loss": float(tr_loss), + "train_acc": float(tr_acc), + "val_auc": float(val_auc), + "val_acc": float(val_acc), + "hld_auc": float(hld_auc_ep), + "hld_acc": float(hld_acc_ep), + }) + 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}] " @@ -366,6 +386,9 @@ def main() -> None: flush=True, ) + import pandas as _pd + _pd.DataFrame(epoch_log_rows).to_csv(fold_dir / "epoch_log.csv", index=False) + if best_state is not None: model.load_state_dict(best_state) diff --git a/scripts/main/v2/run_multifold_v2_modes.py b/scripts/main/v2/run_multifold_v2_modes.py index 01b58fb..58645e6 100644 --- a/scripts/main/v2/run_multifold_v2_modes.py +++ b/scripts/main/v2/run_multifold_v2_modes.py @@ -6,7 +6,7 @@ from __future__ import annotations from pathlib import Path 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)) diff --git a/scripts/output_analysis/aggregate_10x5cv.py b/scripts/output_analysis/aggregate_10x5cv.py new file mode 100644 index 0000000..e3b3895 --- /dev/null +++ b/scripts/output_analysis/aggregate_10x5cv.py @@ -0,0 +1,496 @@ +#!/usr/bin/env python3 +""" +Aggregate and visualise results from a 10× repeated 5-fold CV run. + +Reads probs / y_true from every rep/fold directory, computes per-fold +metrics, and produces: + + outputs/ + fold_metrics.csv — one row per (rep, fold, eval_mode) + rep_metrics.csv — one row per (rep, eval_mode): mean over 5 folds + overall_summary.txt — mean ± SD and 95% CI printed to console + file + {eval_mode}_auc_violin.png + {eval_mode}_roc_mean.png — mean ± 1 SD OVR ROC (all classes or class 1) + {eval_mode}_holdout_roc_mean.png + +Holdout metrics are extracted from the rep-level predictions.npz using the +best_epoch recorded in summary.json, ensemble-averaged over od_fused + os_fused +heads, giving 50 fold-level holdout AUC values (5 folds × 10 reps). + +Usage +----- + python scripts/output_analysis/aggregate_10x5cv.py \ + --run-root analysis_data/pipeline_10x5 \ + --eval-modes binary multiclass \ + --out analysis_data/pipeline_10x5/aggregate +""" +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from sklearn.metrics import roc_auc_score, accuracy_score, roc_curve, auc as sk_auc + +# --------------------------------------------------------------------------- +# Config +# --------------------------------------------------------------------------- + +_CLASS_NAMES = { + "binary": ["Healthy", "Glaucoma"], + "multiclass": ["Healthy", "Glaucoma", "Suspect"], +} + +# Probe files in preference order (first found wins) +# probs_fused = simple OD/OS softmax average (ensemble head — primary metric) +# probs_fused_head = learned logit-level fusion head (worse on average; kept as fallback) +_PROBS_PRIORITY = ["probs_fused.npy", "probs_fused_head.npy"] + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _find_probs(fold_dir: Path) -> Path | None: + for name in _PROBS_PRIORITY: + p = fold_dir / name + if p.exists(): + return p + return None + + +def _load_holdout_summary(mode_dir: Path) -> dict | None: + """ + Read rep-level holdout metrics from summary.json. + + Returns dict with keys auc_mean, auc_std, acc_mean (may be None if missing). + """ + summary_path = mode_dir / "summary.json" + if not summary_path.exists(): + return None + try: + summary = json.loads(summary_path.read_text()) + return summary.get("mode_summary", {}).get("ensemble_holdout") + except Exception: + return None + + +def _auc_macro(y: np.ndarray, p: np.ndarray, num_classes: int) -> float: + try: + if num_classes == 2: + return float(roc_auc_score(y, p[:, 1])) + return float(roc_auc_score(y, p, multi_class="ovr", average="macro")) + except Exception: + return float("nan") + + +def _per_class_roc(y: np.ndarray, p: np.ndarray) -> dict[int, dict]: + 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 + + +def _ci95(values: np.ndarray) -> tuple[float, float]: + """95% CI via t-distribution (two-sided).""" + from scipy import stats as scipy_stats + if len(values) < 2: + return (float("nan"), float("nan")) + ci = scipy_stats.t.interval(0.95, df=len(values) - 1, + loc=np.mean(values), scale=scipy_stats.sem(values)) + return float(ci[0]), float(ci[1]) + + +# --------------------------------------------------------------------------- +# Data loading +# --------------------------------------------------------------------------- + +def load_all_folds(run_root: Path, eval_modes: list[str]) -> pd.DataFrame: + rows = [] + rep_dirs = sorted( + [d for d in run_root.iterdir() if d.is_dir() and d.name.startswith("rep")], + key=lambda d: d.name, + ) + if not rep_dirs: + raise SystemExit(f"No rep* directories found in {run_root}") + + for rep_dir in rep_dirs: + for eval_mode in eval_modes: + tower_mode = "ensemble" + mode_dir = rep_dir / eval_mode / tower_mode + if not mode_dir.exists(): + print(f" [skip] {mode_dir} not found") + continue + num_classes = 2 if eval_mode == "binary" else 3 + + fold_dirs = sorted( + [d for d in mode_dir.iterdir() if d.is_dir() and d.name.startswith("fold")], + key=lambda d: int(d.name[4:]), + ) + for fold_dir in fold_dirs: + fold_idx = int(fold_dir.name[4:]) + y_path = fold_dir / "y_true.npy" + p_path = _find_probs(fold_dir) + + if y_path is None or not y_path.exists() or p_path is None: + print(f" [skip] {rep_dir.name}/{eval_mode}/fold{fold_idx}: missing files") + continue + + 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] + + auc_macro = _auc_macro(y, p, num_classes) + acc = float(accuracy_score(y, p.argmax(1))) + + row = { + "rep": rep_dir.name, + "fold": fold_idx, + "eval_mode": eval_mode, + "probs_file": p_path.name, + "auc_macro": auc_macro, + "acc": acc, + "n": len(y), + } + + # Per-class AUC + for k in range(num_classes): + yb = (y == k).astype(np.uint8) + if yb.sum() > 0 and yb.sum() < len(yb): + try: + row[f"auc_class{k}"] = float(roc_auc_score(yb, p[:, k])) + except Exception: + row[f"auc_class{k}"] = float("nan") + else: + row[f"auc_class{k}"] = float("nan") + + rows.append(row) + + # ---- holdout metrics from rep-level summary.json ---- + # Holdout probs are not stored per-fold; only aggregated stats are saved. + # We attach the rep-level mean to each fold row (same value repeated), + # and also add a single rep-level summary row (fold=-1). + hld_summary = _load_holdout_summary(mode_dir) + if hld_summary: + hld_auc = hld_summary.get("auc_mean", float("nan")) + hld_auc_std = hld_summary.get("auc_std", float("nan")) + hld_acc = hld_summary.get("acc_mean", float("nan")) + for row in rows: + if row["rep"] == rep_dir.name and row["eval_mode"] == eval_mode: + row["hld_auc_macro"] = hld_auc + row["hld_acc"] = hld_acc + # Also store a rep-level holdout row (fold=-1) for direct rep-level analysis + rows.append({ + "rep": rep_dir.name, + "fold": -1, + "eval_mode": eval_mode, + "probs_file": "summary.json", + "auc_macro": float("nan"), + "acc": float("nan"), + "n": float("nan"), + "hld_auc_macro": hld_auc, + "hld_auc_std_within_rep": hld_auc_std, + "hld_acc": hld_acc, + }) + + return pd.DataFrame(rows) + + +# --------------------------------------------------------------------------- +# Plotting +# --------------------------------------------------------------------------- + +def _violin(fold_df: pd.DataFrame, eval_mode: str, out_dir: Path) -> None: + sub = fold_df[fold_df["eval_mode"] == eval_mode].copy() + num_classes = 2 if eval_mode == "binary" else 3 + class_names = _CLASS_NAMES[eval_mode] + + auc_cols = ["auc_macro"] + [f"auc_class{k}" for k in range(num_classes)] + labels = ["Macro AUC"] + [f"AUC {class_names[k]}" for k in range(num_classes)] + present = [(c, l) for c, l in zip(auc_cols, labels) if c in sub.columns] + + data = [sub[c].dropna().values for c, _ in present] + labels = [l for _, l in present] + + fig, ax = plt.subplots(figsize=(max(6, 2 * len(data)), 5)) + parts = ax.violinplot(data, showmedians=True, showextrema=True) + for pc in parts["bodies"]: + pc.set_alpha(0.7) + + # Overlay individual rep means + rep_means = sub.groupby("rep")[auc_cols[0]].mean().values + ax.scatter(np.ones(len(rep_means)), rep_means, zorder=3, + color="k", s=18, alpha=0.6, label="rep mean") + + ax.set_xticks(range(1, len(labels) + 1)) + ax.set_xticklabels(labels, rotation=15, ha="right") + ax.set_ylabel("AUC") + ax.set_title(f"AUC distribution — {eval_mode} (10 × 5-fold, n={len(sub)})") + ax.set_ylim(max(0, sub[auc_cols[0]].min() - 0.05), 1.02) + ax.grid(True, axis="y", linewidth=0.4, alpha=0.5) + ax.legend(fontsize=8) + fig.tight_layout() + path = out_dir / f"{eval_mode}_auc_violin.png" + fig.savefig(path, dpi=160, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {path}") + + +def _mean_roc(fold_df: pd.DataFrame, eval_mode: str, + run_root: Path, out_dir: Path) -> None: + """Mean ± 1 SD OVR ROC across all 50 folds.""" + sub = fold_df[fold_df["eval_mode"] == eval_mode] + num_classes = 2 if eval_mode == "binary" else 3 + class_names = _CLASS_NAMES[eval_mode] + + # Classes to plot (binary: class 1 only) + plot_classes = [1] if eval_mode == "binary" else list(range(num_classes)) + + grid = np.linspace(0, 1, 501) + fig, ax = plt.subplots(figsize=(9, 7)) + ax.plot([0, 1], [0, 1], linestyle="--", linewidth=1, color="grey") + + for k in plot_classes: + tprs, aucs = [], [] + for _, row in sub.iterrows(): + rep_dir = run_root / row["rep"] + fold_dir = rep_dir / eval_mode / "ensemble" / f"fold{int(row['fold'])}" + y_path = fold_dir / "y_true.npy" + p_path = _find_probs(fold_dir) + if not y_path.exists() or p_path is None: + continue + 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] + yb = (y == k).astype(np.uint8) + if yb.sum() == 0 or yb.sum() == len(yb): + continue + fpr, tpr, _ = roc_curve(yb, p[:, k]) + tprs.append(np.interp(grid, fpr, tpr)) + aucs.append(sk_auc(fpr, tpr)) + + if not tprs: + continue + arr = np.vstack(tprs) + mean = arr.mean(0) + std = arr.std(0) + cname = class_names[k] + lbl = f"{cname} AUC {np.nanmean(aucs):.3f} ± {np.nanstd(aucs):.3f}" + line, = ax.plot(grid, mean, linewidth=2, label=lbl) + 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 ± 1 SD OVR ROC — {eval_mode} (10 × 5-fold)") + ax.legend(loc="lower right", fontsize=9) + fig.tight_layout() + path = out_dir / f"{eval_mode}_roc_mean.png" + fig.savefig(path, dpi=160, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {path}") + + +def _holdout_stability(fold_df: pd.DataFrame, eval_mode: str, out_dir: Path) -> None: + """Bar chart of per-rep holdout AUC (mean across folds within rep ± within-rep SD).""" + # use the rep-level rows (fold == -1) which have hld_auc_std_within_rep + rep_rows = fold_df[(fold_df["eval_mode"] == eval_mode) & (fold_df["fold"] == -1)].copy() + if rep_rows.empty or "hld_auc_macro" not in rep_rows.columns: + print(f" [skip] no holdout data for {eval_mode}") + return + rep_rows = rep_rows.sort_values("rep") + + fig, ax = plt.subplots(figsize=(max(6, len(rep_rows) * 0.9), 4)) + x = np.arange(len(rep_rows)) + yerr = rep_rows.get("hld_auc_std_within_rep", pd.Series([0]*len(rep_rows))).fillna(0).values + ax.bar(x, rep_rows["hld_auc_macro"].values, yerr=yerr, + capsize=4, color="darkorange", alpha=0.8) + grand_mean = rep_rows["hld_auc_macro"].mean() + ax.axhline(grand_mean, linestyle="--", color="crimson", + linewidth=1.2, label=f"grand mean = {grand_mean:.3f}") + ax.set_xticks(x) + ax.set_xticklabels(rep_rows["rep"].values, rotation=30, ha="right") + ax.set_ylabel("Holdout macro AUC (mean ± within-rep SD)") + ax.set_title(f"Per-rep holdout stability — {eval_mode}") + ymin = max(0, rep_rows["hld_auc_macro"].min() - 0.05) + ax.set_ylim(ymin, 1.02) + ax.legend(fontsize=9) + ax.grid(True, axis="y", linewidth=0.4, alpha=0.5) + fig.tight_layout() + path = out_dir / f"{eval_mode}_holdout_stability.png" + fig.savefig(path, dpi=160, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {path}") + + +def _rep_stability(fold_df: pd.DataFrame, eval_mode: str, out_dir: Path) -> None: + """Bar chart of per-rep mean macro AUC with ± 1 SD error bars.""" + sub = fold_df[fold_df["eval_mode"] == eval_mode] + rep_stats = sub.groupby("rep")["auc_macro"].agg(["mean", "std"]).reset_index() + rep_stats = rep_stats.sort_values("rep") + + fig, ax = plt.subplots(figsize=(max(6, len(rep_stats) * 0.9), 4)) + x = np.arange(len(rep_stats)) + ax.bar(x, rep_stats["mean"], yerr=rep_stats["std"], + capsize=4, color="steelblue", alpha=0.8) + ax.axhline(rep_stats["mean"].mean(), linestyle="--", color="crimson", + linewidth=1.2, label=f"grand mean = {rep_stats['mean'].mean():.3f}") + ax.set_xticks(x) + ax.set_xticklabels(rep_stats["rep"], rotation=30, ha="right") + ax.set_ylabel("Mean macro AUC (5 folds)") + ax.set_title(f"Per-rep stability — {eval_mode}") + ymin = max(0, rep_stats["mean"].min() - rep_stats["std"].max() - 0.02) + ax.set_ylim(ymin, 1.02) + ax.legend(fontsize=9) + ax.grid(True, axis="y", linewidth=0.4, alpha=0.5) + fig.tight_layout() + path = out_dir / f"{eval_mode}_rep_stability.png" + fig.savefig(path, dpi=160, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {path}") + + +# --------------------------------------------------------------------------- +# Summary text +# --------------------------------------------------------------------------- + +def _print_summary(fold_df: pd.DataFrame, eval_modes: list[str]) -> str: + lines = ["=" * 60, "10 × 5-fold CV — aggregate summary", "=" * 60] + for eval_mode in eval_modes: + sub = fold_df[(fold_df["eval_mode"] == eval_mode) & (fold_df["fold"] >= 0)] + if sub.empty: + continue + num_classes = 2 if eval_mode == "binary" else 3 + class_names = _CLASS_NAMES[eval_mode] + lines.append(f"\n--- {eval_mode.upper()} ---") + lines.append(f" n_folds = {len(sub)}") + + lines.append(" [Validation]") + for col, label in [("auc_macro", "Macro AUC"), ("acc", "Accuracy")]: + if col not in sub.columns: + continue + vals = sub[col].dropna().values + ci_lo, ci_hi = _ci95(vals) + lines.append( + f" {label:18s}: {vals.mean():.4f} ± {vals.std():.4f}" + f" 95% CI [{ci_lo:.4f}, {ci_hi:.4f}]" + ) + + for k in range(num_classes): + col = f"auc_class{k}" + if col not in sub.columns: + continue + vals = sub[col].dropna().values + if len(vals) == 0: + continue + ci_lo, ci_hi = _ci95(vals) + lines.append( + f" AUC {class_names[k]:12s}: {vals.mean():.4f} ± {vals.std():.4f}" + f" 95% CI [{ci_lo:.4f}, {ci_hi:.4f}]" + ) + + # Between-rep variance (val) + rep_means = sub.groupby("rep")["auc_macro"].mean().values + lines.append( + f" Rep-mean AUC (n={len(rep_means)}): " + f"{rep_means.mean():.4f} ± {rep_means.std():.4f}" + f" (between-rep SD = {rep_means.std():.4f})" + ) + + # Holdout — use rep-level rows (fold == -1) + rep_hld = fold_df[ + (fold_df["eval_mode"] == eval_mode) & + (fold_df["fold"] == -1) & + fold_df["hld_auc_macro"].notna() + ]["hld_auc_macro"].values if "hld_auc_macro" in fold_df.columns else np.array([]) + + if len(rep_hld) > 0: + lines.append(" [Holdout] (rep-level means, n_reps={})".format(len(rep_hld))) + ci_lo, ci_hi = _ci95(rep_hld) + lines.append( + f" {'Macro AUC':18s}: {rep_hld.mean():.4f} ± {rep_hld.std():.4f}" + f" 95% CI [{ci_lo:.4f}, {ci_hi:.4f}]" + ) + + lines.append("=" * 60) + return "\n".join(lines) + + +# --------------------------------------------------------------------------- +# CLI +# --------------------------------------------------------------------------- + +def main() -> None: + ap = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + ap.add_argument("--run-root", default="analysis_data/pipeline_10x5", + help="Root directory containing rep* sub-directories.") + ap.add_argument("--eval-modes", nargs="+", + choices=["binary", "multiclass"], + default=["binary", "multiclass"]) + ap.add_argument("--out", default=None, + help="Output directory for plots and CSVs " + "(default: {run-root}/aggregate).") + args = ap.parse_args() + + run_root = Path(args.run_root) + out_dir = Path(args.out) if args.out else run_root / "aggregate" + out_dir.mkdir(parents=True, exist_ok=True) + + print("Loading fold data...") + fold_df = load_all_folds(run_root, args.eval_modes) + if fold_df.empty: + raise SystemExit("No data loaded — check --run-root.") + + fold_df.to_csv(out_dir / "fold_metrics.csv", index=False) + print(f" Saved fold_metrics.csv ({len(fold_df)} rows)") + + rep_df = (fold_df.groupby(["rep", "eval_mode"]) + [["auc_macro", "acc"] + + [c for c in fold_df.columns if c.startswith("auc_class")]] + .mean() + .reset_index()) + rep_df.to_csv(out_dir / "rep_metrics.csv", index=False) + print(f" Saved rep_metrics.csv ({len(rep_df)} rows)") + + summary_text = _print_summary(fold_df, args.eval_modes) + print("\n" + summary_text) + (out_dir / "overall_summary.txt").write_text(summary_text + "\n") + print(f"\n Saved overall_summary.txt") + + print("\nGenerating plots...") + for eval_mode in args.eval_modes: + if fold_df[fold_df["eval_mode"] == eval_mode].empty: + continue + _violin(fold_df, eval_mode, out_dir) + _mean_roc(fold_df, eval_mode, run_root, out_dir) + _rep_stability(fold_df, eval_mode, out_dir) + _holdout_stability(fold_df, eval_mode, out_dir) + + print(f"\nAll outputs written to: {out_dir}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/explainability/aggregate_gradcam.py b/scripts/output_analysis/explainability/aggregate_gradcam.py new file mode 100644 index 0000000..ceeb349 --- /dev/null +++ b/scripts/output_analysis/explainability/aggregate_gradcam.py @@ -0,0 +1,241 @@ +#!/usr/bin/env python3 +""" +Aggregate raw GradCAM heatmaps across all folds for a run. + +For each combination of (eye, class, correct/incorrect) computes: + - mean heatmap + - std heatmap + - count + +Also computes a scalar per patient: fraction of GradCAM attention mass that +falls within the expert-segmented optic disc region (from GT contour files), +using the manifest.csv to locate the contour for each patient/eye. + +Outputs +------- + {out_dir}/mean_heatmaps.npz + Keys: {eye}_{class_name}_{correct|incorrect}_{mean|std|count} + e.g. OD_Glaucoma_correct_mean shape (224, 224) + + {out_dir}/attention_stats.csv + per-patient scalars: patient_id, fold, eye, true_name, pred_name, + correct, confidence, disc_frac, entropy + +Usage +----- + python scripts/output_analysis/explainability/aggregate_gradcam.py \ + --run-dir analysis_data/pipeline_nocrop \ + --eval-mode binary \ + --tower-mode single +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd +from PIL import Image, ImageDraw + + +def load_disc_mask(contour_path: Path, orig_size: tuple[int, int], + cam_h: int, cam_w: int) -> np.ndarray | None: + """ + Load a PAPILA disc contour TXT file, polygon-fill at original image + dimensions, then resize to (cam_h, cam_w). Returns a bool array or + None if the contour cannot be loaded. + """ + try: + arr = np.loadtxt(str(contour_path), dtype=np.float32) + except Exception: + return None + if arr.ndim == 1: + arr = arr.reshape(-1, 2) + if arr.shape[0] < 3 or arr.shape[1] < 2: + return None + + # orig_size is (W, H) as PIL convention + img = Image.new("L", orig_size, 0) + draw = ImageDraw.Draw(img) + draw.polygon([tuple(pt) for pt in arr[:, :2]], fill=1) + mask = np.array(img.resize((cam_w, cam_h), Image.NEAREST), dtype=bool) + return mask + + +def build_disc_lookup(manifest_path: Path) -> dict[tuple[int, str], tuple[Path, tuple[int, int]]]: + """ + Returns {(patient_id_int, eye): (disc_contour_path, (img_W, img_H))}. + Only PAPILA rows are included. + """ + mf = pd.read_csv(manifest_path) + lookup: dict[tuple[int, str], tuple[Path, tuple[int, int]]] = {} + for _, row in mf.iterrows(): + sid = str(row["sample_id"]) + if not sid.startswith("papila_RET"): + continue + # sample_id: papila_RET002OD or papila_RET002OS + suffix = sid[len("papila_RET"):] # e.g. "002OD" + eye = suffix[-2:] # "OD" or "OS" + pid = int(suffix[:-2]) # 2 + disc_path = Path(str(row["annotation_disc"])) + img_path = Path(str(row["image_path"])) + if not disc_path.exists(): + continue + # read original image size once + try: + with Image.open(img_path) as im: + orig_size = im.size # (W, H) + except Exception: + continue + lookup[(pid, eye)] = (disc_path, orig_size) + return lookup + + +def attention_entropy(cam: np.ndarray) -> float: + flat = cam.flatten().astype(np.float64) + flat = flat / (flat.sum() + 1e-12) + return float(-np.sum(flat * np.log(flat + 1e-12))) + + +def load_fold(gradcam_dir: Path): + idx_path = gradcam_dir / "gradcam_index.csv" + if not idx_path.exists(): + return None + idx = pd.read_csv(idx_path) + records = [] + for _, row in idx.iterrows(): + pid = row["patient_id"] + for eye in ("OD", "OS"): + npy = gradcam_dir / f"patient_{pid}_{eye}_cam.npy" + if not npy.exists(): + continue + cam = np.load(npy) + records.append({ + "patient_id": pid, + "eye": eye, + "true_label": int(row["true_label"]), + "true_name": row["true_name"], + "pred_label": int(row["pred_label"]), + "pred_name": row["pred_name"], + "confidence": float(row["confidence"]), + "correct": bool(row["correct"]), + "cam": cam, + }) + return records + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_nocrop") + ap.add_argument("--eval-mode", default="binary") + ap.add_argument("--tower-mode", default="single") + ap.add_argument("--manifest", default="manifest.csv") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + run_dir = Path(args.run_dir) + mode_dir = run_dir / args.eval_mode / args.tower_mode + out_dir = mode_dir / "gradcam_aggregate" + if args.out: + out_dir = Path(args.out) + out_dir.mkdir(parents=True, exist_ok=True) + + # ---- build disc mask lookup ---- + manifest_path = Path(args.manifest) + disc_lookup = build_disc_lookup(manifest_path) + print(f"Disc mask lookup: {len(disc_lookup)} entries from {manifest_path}") + + # ---- collect all records ---- + all_records = [] + stat_rows = [] + fold_dirs = sorted( + [d for d in mode_dir.iterdir() if d.is_dir() and d.name.startswith("fold")], + key=lambda p: int(p.name.replace("fold", "")), + ) + if not fold_dirs: + print(f"No fold dirs found under {mode_dir}") + return + + for fd in fold_dirs: + gcam_dir = fd / "explainability" / "gradcam" + records = load_fold(gcam_dir) + if records is None: + print(f" [skip] {fd.name}: no gradcam_index.csv") + continue + print(f" {fd.name}: {len(records)} eye records") + for r in records: + r["fold"] = fd.name + all_records.append(r) + + if not all_records: + print("No records found — re-run explain_fold.py first.") + return + + print(f"\nTotal eye records: {len(all_records)}") + + h, w = all_records[0]["cam"].shape + + # ---- per-record stats ---- + n_missing = 0 + for r in all_records: + cam = r["cam"] + total = cam.sum() + 1e-12 + pid = int(r["patient_id"]) + eye = r["eye"] + + disc_mask = None + key = (pid, eye) + if key in disc_lookup: + disc_path, orig_size = disc_lookup[key] + disc_mask = load_disc_mask(disc_path, orig_size, h, w) + if disc_mask is None: + n_missing += 1 + disc_frac = float("nan") + else: + disc_frac = float(cam[disc_mask].sum() / total) + + stat_rows.append({ + "patient_id": r["patient_id"], + "fold": r["fold"], + "eye": r["eye"], + "true_name": r["true_name"], + "pred_name": r["pred_name"], + "correct": r["correct"], + "confidence": r["confidence"], + "disc_frac": disc_frac, + "entropy": attention_entropy(cam), + }) + + if n_missing: + print(f" Warning: {n_missing} records had no disc mask (disc_frac=NaN)") + + stats_df = pd.DataFrame(stat_rows) + stats_path = out_dir / "attention_stats.csv" + stats_df.to_csv(stats_path, index=False) + print(f"Saved attention stats → {stats_path}") + + # ---- mean heatmaps ---- + npz_arrays = {} + groups: dict[tuple, list[np.ndarray]] = {} + for r in all_records: + key = (r["eye"], r["true_name"], "correct" if r["correct"] else "incorrect") + groups.setdefault(key, []).append(r["cam"]) + for r in all_records: + key = (r["eye"], r["true_name"], "all") + groups.setdefault(key, []).append(r["cam"]) + + for (eye, cls, split), cams in groups.items(): + stack = np.stack(cams, axis=0) + key_base = f"{eye}_{cls}_{split}" + npz_arrays[f"{key_base}_mean"] = stack.mean(axis=0).astype(np.float32) + npz_arrays[f"{key_base}_std"] = stack.std(axis=0).astype(np.float32) + npz_arrays[f"{key_base}_count"] = np.array(len(cams)) + print(f" {key_base}: N={len(cams)}") + + npz_path = out_dir / "mean_heatmaps.npz" + np.savez_compressed(npz_path, **npz_arrays) + print(f"Saved mean heatmaps → {npz_path}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/explainability/explain_fold.py b/scripts/output_analysis/explainability/explain_fold.py index ae4a4c7..7415bba 100644 --- a/scripts/output_analysis/explainability/explain_fold.py +++ b/scripts/output_analysis/explainability/explain_fold.py @@ -409,6 +409,7 @@ def run_gradcam( overlay_grid_items: list[ tuple[Image.Image | None, Image.Image | None, str, bool] ] = [] + index_rows: list[dict] = [] model.eval() for batch in loader: @@ -492,6 +493,19 @@ def run_gradcam( flush=True, ) + # Save raw CAM arrays + np.save(gradcam_dir / f"patient_{pid}_OD_cam.npy", cam_od) + np.save(gradcam_dir / f"patient_{pid}_OS_cam.npy", cam_os) + index_rows.append({ + "patient_id": pid, + "true_label": label, + "true_name": true_name, + "pred_label": pred, + "pred_name": pred_name, + "confidence": conf, + "correct": correct, + }) + # Accumulate for summary grid od_overlay = overlay_gradcam(orig_od, cam_od, alpha) if orig_od else None os_overlay = overlay_gradcam(orig_os, cam_os, alpha) if orig_os else None @@ -500,6 +514,16 @@ def run_gradcam( gcam.remove() + # ---- save index CSV ---- + if index_rows: + import csv + idx_path = gradcam_dir / "gradcam_index.csv" + with idx_path.open("w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=list(index_rows[0].keys())) + writer.writeheader() + writer.writerows(index_rows) + print(f" Index CSV → {idx_path}", flush=True) + # ---- summary grid: N_patients rows × 2 cols (OD overlay | OS overlay) ---- n = len(overlay_grid_items) if n == 0: diff --git a/scripts/output_analysis/explainability/explain_run.py b/scripts/output_analysis/explainability/explain_run.py index 1c47d2f..96d6d9d 100644 --- a/scripts/output_analysis/explainability/explain_run.py +++ b/scripts/output_analysis/explainability/explain_run.py @@ -362,94 +362,114 @@ def run_gradcam( gradcam_dir = out_dir / "gradcam" gradcam_dir.mkdir(exist_ok=True) - target_layer = get_gradcam_layer(model, backbone) - gcam = GradCAM(target_layer) - - num_classes = model.bridge.classifier_fused[-1].out_features + manifest_path = gradcam_dir / "gradcam_manifest.csv" + skip_overlays = manifest_path.exists() overlay_grid_items = [] - model.eval() - for batch in loader: - img_od = batch["image_1"].to(device) - img_os = batch["image_2"].to(device) - meta_od = batch["matrix_1"].to(device) - meta_os = batch["matrix_2"].to(device) - 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] + if skip_overlays: + print(" Overlays already exist — skipping computation, loading from disk.", flush=True) + manifest_df = pd.read_csv(manifest_path) + for _, row in manifest_df.iterrows(): + pid = str(row["pid"]) + od_path = gradcam_dir / f"gradcam_od_{pid}.png" + os_path = gradcam_dir / f"gradcam_os_{pid}.png" + od_ov = Image.open(od_path).convert("RGB") if od_path.exists() else None + os_ov = Image.open(os_path).convert("RGB") if os_path.exists() else None + overlay_grid_items.append((od_ov, os_ov, str(row["short_lbl"]), bool(row["correct"]))) + else: + target_layer = get_gradcam_layer(model, backbone) + gcam = GradCAM(target_layer) + manifest_rows = [] - cam_od, pred = gcam.compute(img_od, meta_od, model) - cam_os, _ = gcam.compute(img_os, meta_os, model) + model.eval() + for batch in loader: + img_od = batch["image_1"].to(device) + img_os = batch["image_2"].to(device) + meta_od = batch["matrix_1"].to(device) + meta_os = batch["matrix_2"].to(device) + 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] - with torch.no_grad(): - out_od = model(img_od, meta_od) - conf = float(torch.softmax(out_od, dim=1)[0, pred].item()) + cam_od, pred = gcam.compute(img_od, meta_od, model) + cam_os, _ = gcam.compute(img_os, meta_os, model) - row_od = eval_df[ - (eval_df["Patient ID"] == int(pid)) & (eval_df["eyeID"] == "OD") - ] - row_os = eval_df[ - (eval_df["Patient ID"] == int(pid)) & (eval_df["eyeID"] == "OS") - ] - orig_od = ( - Image.open(data.get_image_path(row_od.iloc[0])).convert("RGB") - if len(row_od) else None - ) - orig_os = ( - Image.open(data.get_image_path(row_os.iloc[0])).convert("RGB") - if len(row_os) else None - ) + with torch.no_grad(): + out_od = model(img_od, meta_od) + conf = float(torch.softmax(out_od, dim=1)[0, pred].item()) - true_name = label_name(label, eval_mode) - pred_name = label_name(pred, eval_mode) - correct = label == pred - title = ( - f"Patient {pid} | True: {true_name} | Pred: {pred_name} " - f"| conf={conf:.2f} {'✓' if correct else '✗'}" - ) + row_od = eval_df[ + (eval_df["Patient ID"] == int(pid)) & (eval_df["eyeID"] == "OD") + ] + row_os = eval_df[ + (eval_df["Patient ID"] == int(pid)) & (eval_df["eyeID"] == "OS") + ] + orig_od = ( + Image.open(data.get_image_path(row_od.iloc[0])).convert("RGB") + if len(row_od) else None + ) + orig_os = ( + Image.open(data.get_image_path(row_os.iloc[0])).convert("RGB") + if len(row_os) else None + ) - fig, axes = plt.subplots(2, 2, figsize=(10, 9)) - fig.suptitle(title, fontsize=11, fontweight="bold", color="green" if correct else "red") + true_name = label_name(label, eval_mode) + pred_name = label_name(pred, eval_mode) + correct = label == pred + title = ( + f"Patient {pid} | True: {true_name} | Pred: {pred_name} " + f"| conf={conf:.2f} {'✓' if correct else '✗'}" + ) - if orig_od is not None: - axes[0, 0].imshow(orig_od) - axes[0, 0].set_title("OD — original", fontsize=9) - axes[0, 1].imshow(overlay_gradcam(orig_od, cam_od, alpha)) - axes[0, 1].set_title("OD — GradCAM", fontsize=9) - else: - axes[0, 0].set_title("OD — (missing)", fontsize=9) - axes[0, 0].axis("off") - axes[0, 1].axis("off") + fig, axes = plt.subplots(2, 2, figsize=(10, 9)) + fig.suptitle(title, fontsize=11, fontweight="bold", color="green" if correct else "red") - if orig_os is not None: - axes[1, 0].imshow(orig_os) - axes[1, 0].set_title("OS — original", fontsize=9) - axes[1, 1].imshow(overlay_gradcam(orig_os, cam_os, alpha)) - axes[1, 1].set_title("OS — GradCAM", fontsize=9) - else: - axes[1, 0].set_title("OS — (missing)", fontsize=9) - axes[1, 0].axis("off") - axes[1, 1].axis("off") + if orig_od is not None: + axes[0, 0].imshow(orig_od) + axes[0, 0].set_title("OD — original", fontsize=9) + axes[0, 1].imshow(overlay_gradcam(orig_od, cam_od, alpha)) + axes[0, 1].set_title("OD — GradCAM", fontsize=9) + else: + axes[0, 0].set_title("OD — (missing)", fontsize=9) + axes[0, 0].axis("off") + axes[0, 1].axis("off") - fig.tight_layout() - out_path = gradcam_dir / f"patient_{pid}_OD_OS.png" - fig.savefig(out_path, dpi=120) - plt.close(fig) - print(f" Patient {pid}: {true_name} → {pred_name} ({conf:.2f}) → {out_path.name}", flush=True) + if orig_os is not None: + axes[1, 0].imshow(orig_os) + axes[1, 0].set_title("OS — original", fontsize=9) + axes[1, 1].imshow(overlay_gradcam(orig_os, cam_os, alpha)) + axes[1, 1].set_title("OS — GradCAM", fontsize=9) + else: + axes[1, 0].set_title("OS — (missing)", fontsize=9) + axes[1, 0].axis("off") + axes[1, 1].axis("off") - od_overlay = overlay_gradcam(orig_od, cam_od, alpha) if orig_od else None - os_overlay = overlay_gradcam(orig_os, cam_os, alpha) if orig_os else None - short_lbl = f"P{pid} {true_name[:3]}→{pred_name[:3]} {'✓' if correct else '✗'}" - overlay_grid_items.append((od_overlay, os_overlay, short_lbl, correct)) + fig.tight_layout() + out_path = gradcam_dir / f"patient_{pid}_OD_OS.png" + fig.savefig(out_path, dpi=120) + plt.close(fig) + print(f" Patient {pid}: {true_name} → {pred_name} ({conf:.2f}) → {out_path.name}", flush=True) - gcam.remove() + od_overlay = overlay_gradcam(orig_od, cam_od, alpha) if orig_od else None + os_overlay = overlay_gradcam(orig_os, cam_os, alpha) if orig_os else None + if od_overlay is not None: + od_overlay.save(gradcam_dir / f"gradcam_od_{pid}.png") + if os_overlay is not None: + os_overlay.save(gradcam_dir / f"gradcam_os_{pid}.png") + short_lbl = f"P{pid} {true_name[:3]}→{pred_name[:3]} {'✓' if correct else '✗'}" + overlay_grid_items.append((od_overlay, os_overlay, short_lbl, correct)) + manifest_rows.append({"pid": pid, "short_lbl": short_lbl, "correct": correct}) + gcam.remove() + pd.DataFrame(manifest_rows).to_csv(manifest_path, index=False) + + # Always regenerate the summary grid n = len(overlay_grid_items) if n == 0: print(" [Phase 2] No patients to visualise.", flush=True) return - fig, axes = plt.subplots(n, 2, figsize=(8, n * 3.2 + 0.8)) + fig, axes = plt.subplots(n, 2, figsize=(8, n * 3.2 + 1.5)) if n == 1: axes = axes[np.newaxis, :] fig.suptitle("GradCAM Summary Grid — all holdout patients", fontsize=12) @@ -465,7 +485,7 @@ def run_gradcam( axes[i, 1].imshow(os_ov) axes[i, 1].set_title(f"{lbl}\nOS", fontsize=7, color=color) - fig.tight_layout() + fig.tight_layout(rect=[0, 0, 1, 0.97]) grid_path = out_dir / "gradcam_summary_grid.png" fig.savefig(grid_path, dpi=120) plt.close(fig) @@ -496,6 +516,7 @@ def _fusion_event_stats( pm: np.ndarray, split_name: str, out_dir: Path, + component_labels: tuple[str, str] = ("img", "md"), ) -> dict: """Compute, save, and plot fusion events for one split. Returns summary dict.""" N = len(y_true) @@ -503,6 +524,16 @@ def _fusion_event_stats( print(f" [{split_name}] No samples — skipping.", flush=True) return {} + a, b = component_labels + event_labels = [ + "full correction\n(both wrong→fused right)", + f"{a} assist\n({a} wrong, {b} right→right)", + f"{b} assist\n({b} wrong, {a} right→right)", + "full error\n(both right→fused wrong)", + f"{a} drag\n({a} wrong, {b} right→wrong)", + f"{b} drag\n({b} wrong, {a} right→wrong)", + ] + pred_f = pf.argmax(axis=1) pred_i = pi.argmax(axis=1) pred_m = pm.argmax(axis=1) @@ -529,7 +560,7 @@ def _fusion_event_stats( counts = [int(m.sum()) for m in event_masks] print(f"\n [{split_name}] N={N}", flush=True) - for label, count in zip(_EVENT_LABELS, counts): + for label, count in zip(event_labels, counts): print(f" {label.replace(chr(10), ' '):55s}: {count}", flush=True) n_corr, n_err = counts[0], counts[3] print(f" full correction/error ratio: {n_corr}/{n_err}", flush=True) @@ -582,11 +613,11 @@ def _fusion_event_stats( axes[0].set_xticklabels(["Positive\nevents", "Negative\nevents"]) axes[0].set_ylabel("Count") patches = [mpatches.Patch(color=c, label=l.replace("\n", " ")) - for c, l in zip(_EVENT_COLORS, _EVENT_LABELS)] + for c, l in zip(_EVENT_COLORS, event_labels)] axes[0].legend(handles=patches, fontsize=6, loc="upper right") box_data = [conf_delta[m] for m in event_masks if m.sum() > 0] - box_labels = [l.split("\n")[0] for m, l in zip(event_masks, _EVENT_LABELS) if m.sum() > 0] + box_labels = [l.split("\n")[0] for m, l in zip(event_masks, event_labels) if m.sum() > 0] box_cols = [c for m, c in zip(event_masks, _EVENT_COLORS) if m.sum() > 0] if box_data: bp = axes[1].boxplot(box_data, patch_artist=True, widths=0.5) @@ -595,10 +626,10 @@ def _fusion_event_stats( axes[1].set_xticks(range(1, len(box_labels) + 1)) axes[1].set_xticklabels(box_labels, rotation=35, ha="right", fontsize=7) axes[1].axhline(0, color="black", linewidth=0.8, linestyle="--") - axes[1].set_ylabel("conf_delta\n(fused − avg(img, md))") + axes[1].set_ylabel(f"conf_delta\n(fused − avg({a}, {b}))") axes[1].set_title("Confidence delta by event type") - for mask, color, label in zip(event_masks, _EVENT_COLORS, _EVENT_LABELS): + for mask, color, label in zip(event_masks, _EVENT_COLORS, event_labels): if mask.sum() > 0: axes[2].scatter(conf_i[mask], conf_m[mask], c=color, label=label.split("\n")[0], alpha=0.85, s=45, edgecolors="none") @@ -609,8 +640,8 @@ def _fusion_event_stats( axes[2].scatter(conf_i[concordant_bad], conf_m[concordant_bad], c="darkgrey", alpha=0.4, s=20, edgecolors="none", label="concordant wrong") axes[2].plot([0, 1], [0, 1], "k--", linewidth=0.5, alpha=0.4) - axes[2].set_xlabel("conf_img") - axes[2].set_ylabel("conf_md") + axes[2].set_xlabel(f"conf_{a}") + axes[2].set_ylabel(f"conf_{b}") axes[2].set_title("Tower confidence space\ncoloured by fusion event") axes[2].legend(fontsize=6, loc="lower right") @@ -711,12 +742,13 @@ def run_fusion_event_analysis(fold_dir: Path, out_dir: Path) -> list[dict]: pm_ens = 0.5 * (pm_od + pm_os) summaries.append(_fusion_event_stats(y_val, pf_ens, pi_ens, pm_ens, "val", out_dir)) - # Fused head (if available): learned bilateral combination vs averaged towers + # Fused head (if available): learned bilateral combination vs per-eye bridge outputs fused_head_f = fold_dir / "probs_fused_head.npy" if fused_head_f.exists(): pf_head = np.load(fused_head_f) summaries.append(_fusion_event_stats( - y_val, pf_head, pi_ens, pm_ens, "val_fused_head", out_dir)) + y_val, pf_head, pf_od, pf_os, "val_fused_head", out_dir, + component_labels=("OD", "OS"))) else: print(" Val epoch files not found — skipping val.", flush=True) @@ -745,6 +777,95 @@ def run_fusion_event_analysis(fold_dir: Path, out_dir: Path) -> list[dict]: return summaries +# --------------------------------------------------------------------------- +# Cross-fold MD importance summary plot +# --------------------------------------------------------------------------- + + +def _plot_run_md_importance_summary(run_dir: Path, folds: list) -> None: + """Aggregate per-fold MD permutation importance CSVs into a run-level summary plot.""" + all_dfs = [] + for fold_idx, fold_dir, _ in folds: + csv_path = fold_dir / "explainability" / "md_permutation_importance.csv" + if csv_path.exists(): + df = pd.read_csv(csv_path) + df["fold"] = fold_idx + all_dfs.append(df) + + if not all_dfs: + print(" [MD summary] No per-fold importance CSVs found — skipping.", flush=True) + return + + combined = pd.concat(all_dfs, ignore_index=True) + _SPECIAL = {"TOTAL_MD_ABLATION", "GAUSSIAN_NOISE_ABLATION"} + feature_rows = combined[~combined["feature"].isin(_SPECIAL)] + special_rows = combined[combined["feature"].isin(_SPECIAL)] + + agg = ( + feature_rows.groupby("feature")["importance"] + .agg(["mean", "std"]) + .reset_index() + .rename(columns={"mean": "mean_importance", "std": "std_importance"}) + .sort_values("mean_importance", ascending=False) + .reset_index(drop=True) + ) + special_agg = ( + special_rows.groupby("feature")["importance"] + .agg(["mean", "std"]) + .reset_index() + ) + + names = agg["feature"].tolist() + imps = agg["mean_importance"].tolist() + stds = agg["std_importance"].fillna(0).tolist() + colors = ["#e05c5c" if v >= 0 else "#5c9ee0" for v in imps] + + fig, ax = plt.subplots(figsize=(9, max(4, (len(names) + 3) * 0.45))) + y_pos = np.arange(len(names)) + ax.barh(y_pos, imps, xerr=stds, color=colors, ecolor="grey", capsize=3, height=0.6) + ax.axhline(len(names) - 0.25, color="grey", linewidth=0.6, linestyle="--") + + special_label_map = { + "TOTAL_MD_ABLATION": "ALL MD (permute)", + "GAUSSIAN_NOISE_ABLATION": "ALL MD (noise)", + } + special_colors = { + "TOTAL_MD_ABLATION": "#c45ce0", + "GAUSSIAN_NOISE_ABLATION": "#e08c2a", + } + extra_ytick_pos = [] + extra_ytick_labels = [] + for i, feat in enumerate(["TOTAL_MD_ABLATION", "GAUSSIAN_NOISE_ABLATION"]): + row = special_agg[special_agg["feature"] == feat] + if row.empty: + continue + offset = len(names) + 0.5 + i + val, err = float(row["mean"].iloc[0]), float(row["std"].iloc[0]) + ax.barh(offset, val, xerr=err, + color=special_colors[feat] if val >= 0 else "#5c9ee0", + ecolor="grey", capsize=3, height=0.6) + extra_ytick_pos.append(offset) + extra_ytick_labels.append(special_label_map[feat]) + + ax.set_yticks(list(y_pos) + extra_ytick_pos) + ax.set_yticklabels(names + extra_ytick_labels, fontsize=9) + ax.invert_yaxis() + ax.axvline(0, color="black", linewidth=0.8) + ax.set_xlabel("Mean AUC drop (baseline − permuted)", fontsize=10) + ax.set_title( + f"MD Tower — Permutation Feature Importance ({len(all_dfs)}-fold summary)\n" + f"error bars = std across folds", + fontsize=11, + ) + fig.tight_layout() + out_path = run_dir / "explainability_md_importance_summary.png" + fig.savefig(out_path, dpi=150) + plt.close(fig) + print(f" MD importance summary → {out_path}", flush=True) + + agg.to_csv(run_dir / "explainability_md_importance_summary.csv", index=False) + + # --------------------------------------------------------------------------- # Cross-fold fusion summary plot # --------------------------------------------------------------------------- @@ -761,6 +882,17 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None n_folds = len(grp) fold_ids = grp["fold"].values + # For the fused_head split the two comparators are OD/OS bridges, not img/md towers + comp_a, comp_b = ("OD", "OS") if "fused_head" in split_name else ("img", "md") + summary_event_labels = [ + "full correction\n(both wrong→fused right)", + f"{comp_a} assist\n({comp_a} wrong, {comp_b} right→right)", + f"{comp_b} assist\n({comp_b} wrong, {comp_a} right→right)", + "full error\n(both right→fused wrong)", + f"{comp_a} drag\n({comp_a} wrong, {comp_b} right→wrong)", + f"{comp_b} drag\n({comp_b} wrong, {comp_a} right→wrong)", + ] + # Load all per-fold CSVs for this split to get sample-level data sample_dfs = [] for fold_idx in fold_ids: @@ -790,7 +922,7 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None axes[0].set_xticklabels(["Positive\nevents", "Negative\nevents"]) axes[0].set_ylabel("Count (all folds)") patches = [mpatches.Patch(color=c, label=l.replace("\n", " ")) - for c, l in zip(_EVENT_COLORS, _EVENT_LABELS)] + for c, l in zip(_EVENT_COLORS, summary_event_labels)] axes[0].legend(handles=patches, fontsize=6, loc="upper right") # Panel 2: per-fold stacked bar (fold variance) @@ -819,7 +951,7 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None axes[2].axhline(0, color="grey", linewidth=0.7) axes[2].set_xticks(x) axes[2].set_xticklabels([f"fold {f}" for f in fold_ids], fontsize=8) - axes[2].set_ylabel("conf_delta mean\n(fused − avg(img, md))") + axes[2].set_ylabel(f"conf_delta mean\n(fused − avg({comp_a}, {comp_b}))") axes[2].set_title("Confidence delta per fold") axes[2].legend(fontsize=8) @@ -829,7 +961,7 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None box_data = [sample_df.loc[sample_df["event_type"] == k, "conf_delta"].values for k in key_order] box_labels = [l.split("\n")[0] - for k, l in zip(_EVENT_KEYS, _EVENT_LABELS) if k in key_order] + for k, l in zip(_EVENT_KEYS, summary_event_labels) if k in key_order] box_cols = [c for k, c in zip(_EVENT_KEYS, _EVENT_COLORS) if k in key_order] if box_data: bp = axes[3].boxplot(box_data, patch_artist=True, widths=0.5) @@ -838,7 +970,7 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None axes[3].set_xticks(range(1, len(box_labels) + 1)) axes[3].set_xticklabels(box_labels, rotation=35, ha="right", fontsize=7) axes[3].axhline(0, color="black", linewidth=0.8, linestyle="--") - axes[3].set_ylabel("conf_delta\n(fused − avg(img, md))") + axes[3].set_ylabel(f"conf_delta\n(fused − avg({comp_a}, {comp_b}))") axes[3].set_title("Confidence delta by event type\n(all folds)") # Panel 5: tower confidence space scatter (all folds combined) @@ -847,7 +979,7 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None for key, color in zip(_EVENT_KEYS, _EVENT_COLORS): sub = sample_df[sample_df["event_type"] == key] if len(sub): - label = next(l.split("\n")[0] for k, l in zip(_EVENT_KEYS, _EVENT_LABELS) + label = next(l.split("\n")[0] for k, l in zip(_EVENT_KEYS, summary_event_labels) if k == key) axes[4].scatter(sub["conf_img"], sub["conf_md"], c=color, label=label, alpha=0.7, s=30, edgecolors="none") @@ -860,8 +992,8 @@ def _plot_cross_fold_fusion_summary(df_sum: pd.DataFrame, run_dir: Path) -> None axes[4].scatter(sub["conf_img"], sub["conf_md"], c=conc_color, alpha=0.3, s=15, edgecolors="none", label=conc_label) axes[4].plot([0, 1], [0, 1], "k--", linewidth=0.5, alpha=0.4) - axes[4].set_xlabel("conf_img") - axes[4].set_ylabel("conf_md") + axes[4].set_xlabel(f"conf_{comp_a}") + axes[4].set_ylabel(f"conf_{comp_b}") axes[4].legend(fontsize=6, loc="lower right") axes[4].set_title("Tower confidence space\n(all folds)") @@ -1105,6 +1237,12 @@ def main(): del model torch.cuda.empty_cache() + # ---- Cross-fold MD importance summary ---- + if not args.no_phase1: + print(f"\n{'='*60}", flush=True) + print("[explain_run] === Cross-fold MD importance summary ===", flush=True) + _plot_run_md_importance_summary(run_dir, folds) + # ---- Cross-fold Phase 3 summary ---- if all_phase3_summaries and not args.no_phase3: print(f"\n{'='*60}", flush=True) diff --git a/scripts/output_analysis/explainability/plot_disc_attention_detail.py b/scripts/output_analysis/explainability/plot_disc_attention_detail.py new file mode 100644 index 0000000..95b9db6 --- /dev/null +++ b/scripts/output_analysis/explainability/plot_disc_attention_detail.py @@ -0,0 +1,284 @@ +#!/usr/bin/env python3 +""" +Granular disc-attention visualisation. + +Layout (3 rows × N_class cols): + Row 0 — disc-centred mean GradCAM patch for CORRECT predictions + Row 1 — disc-centred mean GradCAM patch for INCORRECT predictions + Row 2 — per-patient strip plot of disc_frac (blue=correct, red=incorrect) + +Disc-centred patches: each patient's CAM is translated and scaled so the GT +disc centroid sits at the patch centre before averaging. A dashed white circle +marks the average GT disc size. This makes cross-patient averaging meaningful +regardless of where the disc sits in the original image. + +Usage +----- + python scripts/output_analysis/explainability/plot_disc_attention_detail.py \ + --agg-dir analysis_data/pipeline_nocrop/binary/single/gradcam_aggregate \ + --manifest manifest.csv +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.patches import Circle +from PIL import Image, ImageDraw + + +DISC_SPAN = 5 # patch side = DISC_SPAN × disc diameter +OUTPUT_SIZE = 96 # pixel size of each thumbnail + + +# --------------------------------------------------------------------------- +# Disc mask helpers +# --------------------------------------------------------------------------- + +def _load_disc_mask(contour_path: Path, orig_size: tuple, cam_h: int, cam_w: int): + try: + arr = np.loadtxt(str(contour_path), dtype=np.float32) + except Exception: + return None + if arr.ndim == 1: + arr = arr.reshape(-1, 2) + if arr.shape[0] < 3 or arr.shape[1] < 2: + return None + img = Image.new("L", orig_size, 0) + ImageDraw.Draw(img).polygon([tuple(pt) for pt in arr[:, :2]], fill=1) + return np.array(img.resize((cam_w, cam_h), Image.NEAREST), dtype=bool) + + +def build_disc_lookup(manifest_path: Path) -> dict: + mf = pd.read_csv(manifest_path) + lookup: dict = {} + for _, row in mf.iterrows(): + sid = str(row["sample_id"]) + if not sid.startswith("papila_RET"): + continue + suffix = sid[len("papila_RET"):] + eye = suffix[-2:] + pid = int(suffix[:-2]) + disc_path = Path(str(row["annotation_disc"])) + img_path = Path(str(row["image_path"])) + if not disc_path.exists(): + continue + try: + with Image.open(img_path) as im: + orig_size = im.size + except Exception: + continue + lookup[(pid, eye)] = (disc_path, orig_size) + return lookup + + +# --------------------------------------------------------------------------- +# Disc-centred patch extraction +# --------------------------------------------------------------------------- + +def disc_centered_patch( + cam: np.ndarray, + disc_mask: np.ndarray, + span: int = DISC_SPAN, + out: int = OUTPUT_SIZE, +) -> tuple[np.ndarray | None, float | None]: + """ + Return (patch, disc_r_out): + patch — (out, out) float32 in [0, 1] + disc_r_out — disc radius in patch-pixel units (for drawing reference circle) + """ + if disc_mask is None or disc_mask.sum() == 0: + return None, None + + ys, xs = np.where(disc_mask) + cy, cx = ys.mean(), xs.mean() + disc_r = float(np.sqrt(disc_mask.sum() / np.pi)) + half = max(1, int(round(span * disc_r / 2))) + + h, w = cam.shape + y0, y1 = int(round(cy)) - half, int(round(cy)) + half + x0, x1 = int(round(cx)) - half, int(round(cx)) + half + + pt = max(0, -y0); pb = max(0, y1 - h) + pl = max(0, -x0); pr = max(0, x1 - w) + cam_pad = np.pad(cam, ((pt, pb), (pl, pr)), constant_values=0.0) + + patch = cam_pad[y0 + pt : y1 + pt, x0 + pl : x1 + pl] + patch_img = Image.fromarray((np.clip(patch, 0, 1) * 255).astype(np.uint8)) + patch_out = np.array(patch_img.resize((out, out), Image.BILINEAR)) / 255.0 + + disc_r_out = out * disc_r / (2 * half) + return patch_out.astype(np.float32), disc_r_out + + +# --------------------------------------------------------------------------- +# Data loading +# --------------------------------------------------------------------------- + +def load_all_cam_records(mode_dir: Path) -> list[dict]: + records = [] + fold_dirs = sorted( + [d for d in mode_dir.iterdir() if d.is_dir() and d.name.startswith("fold")], + key=lambda p: int(p.name.replace("fold", "")), + ) + for fd in fold_dirs: + gcam_dir = fd / "explainability" / "gradcam" + idx_path = gcam_dir / "gradcam_index.csv" + if not idx_path.exists(): + continue + idx = pd.read_csv(idx_path) + for _, row in idx.iterrows(): + pid = int(row["patient_id"]) + for eye in ("OD", "OS"): + npy = gcam_dir / f"patient_{pid}_{eye}_cam.npy" + if not npy.exists(): + continue + records.append({ + "patient_id": pid, + "eye": eye, + "true_name": row["true_name"], + "correct": bool(row["correct"]), + "cam": np.load(npy), + }) + return records + + +def build_mean_patches( + records: list[dict], + disc_lookup: dict, + classes: list[str], +) -> dict[tuple, tuple]: + """ + Returns {(cls, split): (mean_patch, mean_disc_r_out, count)} + split = 'correct' | 'incorrect' + """ + buckets: dict[tuple, list] = {} + radii: dict[tuple, list] = {} + + for r in records: + split = "correct" if r["correct"] else "incorrect" + key = (r["true_name"], split) + pid, eye = r["patient_id"], r["eye"] + if (pid, eye) not in disc_lookup: + continue + disc_path, orig_size = disc_lookup[(pid, eye)] + h, w = r["cam"].shape + disc_mask = _load_disc_mask(disc_path, orig_size, h, w) + patch, disc_r_out = disc_centered_patch(r["cam"], disc_mask) + if patch is None: + continue + buckets.setdefault(key, []).append(patch) + radii.setdefault(key, []).append(disc_r_out) + + return { + key: (np.stack(ps).mean(0), float(np.mean(radii[key])), len(ps)) + for key, ps in buckets.items() + } + + +# --------------------------------------------------------------------------- +# Plotting +# --------------------------------------------------------------------------- + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--agg-dir", required=True) + ap.add_argument("--manifest", default="manifest.csv") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + agg_dir = Path(args.agg_dir) + mode_dir = agg_dir.parent + out_path = Path(args.out) if args.out else agg_dir / "disc_attention_detail.png" + + print("Loading disc lookup…") + disc_lookup = build_disc_lookup(Path(args.manifest)) + print(f" {len(disc_lookup)} entries") + + print("Loading CAM records…") + records = load_all_cam_records(mode_dir) + print(f" {len(records)} eye records") + + stats = pd.read_csv(agg_dir / "attention_stats.csv") + classes = sorted({r["true_name"] for r in records}) + n_cls = len(classes) + print(f"Classes: {classes}") + + print("Building disc-centred mean patches…") + mean_patches = build_mean_patches(records, disc_lookup, classes) + + # ----------------------------------------------------------------------- + # Figure + # ----------------------------------------------------------------------- + corr_colors = {"correct": "steelblue", "incorrect": "tomato"} + splits = ["correct", "incorrect"] + row_labels = ["Correct", "Incorrect", "Disc fraction\n(strip plot)"] + + fig, axes = plt.subplots(3, n_cls, figsize=(4.2 * n_cls, 13)) + if n_cls == 1: + axes = axes[:, np.newaxis] + + # ---- rows 0 & 1: disc-centred heatmaps ---- + for ri, split in enumerate(splits): + for ci, cls in enumerate(classes): + ax = axes[ri, ci] + key = (cls, split) + if key in mean_patches: + mean_patch, disc_r_out, count = mean_patches[key] + ax.imshow(mean_patch, cmap="jet", vmin=0, vmax=1, origin="upper", + extent=[0, OUTPUT_SIZE, OUTPUT_SIZE, 0]) + cx = cy = OUTPUT_SIZE / 2 + ax.add_patch(Circle((cx, cy), disc_r_out, + fill=False, edgecolor="white", + linewidth=2, linestyle="--")) + ax.set_title(f"{cls} | {split}\n(N={count})", fontsize=9) + else: + ax.text(0.5, 0.5, "no data", ha="center", va="center", + transform=ax.transAxes, fontsize=9, color="grey") + ax.set_title(f"{cls} | {split}", fontsize=9) + ax.axis("off") + axes[ri, 0].set_ylabel(row_labels[ri], fontsize=10, labelpad=6) + + # ---- row 2: strip plots ---- + rng = np.random.default_rng(42) + for ci, cls in enumerate(classes): + ax = axes[2, ci] + sub = stats[stats["true_name"] == cls].dropna(subset=["disc_frac"]) + + for xi, split in enumerate(splits): + correct_val = (split == "correct") + pts = sub[sub["correct"] == correct_val]["disc_frac"].values + if len(pts) == 0: + continue + color = corr_colors[split] + jitter = rng.uniform(-0.18, 0.18, size=len(pts)) + ax.scatter(xi + jitter, pts, color=color, alpha=0.7, s=28, edgecolors="none") + ax.hlines(pts.mean(), xi - 0.28, xi + 0.28, + colors=color, linewidth=2.5, zorder=5) + + ax.set_xticks([0, 1]) + ax.set_xticklabels(["Correct", "Incorrect"], fontsize=9) + ax.set_xlim(-0.55, 1.55) + ax.set_ylim(0, 1) + ax.set_title(cls, fontsize=10) + ax.grid(axis="y", linestyle="--", alpha=0.3) + if ci == 0: + ax.set_ylabel("Disc fraction\n(GT disc attention)", fontsize=9) + + fig.suptitle( + "Disc-centred GradCAM attention | dashed circle = GT disc boundary", + fontsize=12, + ) + fig.tight_layout() + fig.savefig(out_path, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out_path}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/explainability/plot_gradcam_aggregate.py b/scripts/output_analysis/explainability/plot_gradcam_aggregate.py new file mode 100644 index 0000000..5da2089 --- /dev/null +++ b/scripts/output_analysis/explainability/plot_gradcam_aggregate.py @@ -0,0 +1,286 @@ +#!/usr/bin/env python3 +""" +Visualise aggregated GradCAM heatmaps produced by aggregate_gradcam.py. + +Produces two figures: + +Figure 1 — Mean heatmaps grid + Rows: classes (e.g. Normal, Glaucoma) + Cols: OD_all | OS_all | OD_correct | OD_incorrect | OS_correct | OS_incorrect + +Figure 2 — Attention stats + Panel A: disc_frac distribution per class (violin/box), OD and OS side by side + Panel B: entropy distribution per class + Panel C: disc_frac correct vs incorrect per class (scatter means + error bars) + +Figure 3 — Disc attention vs correct confidence + Scatter of disc_frac vs correct_conf (confidence if correct, 1-confidence if wrong) + One panel per class, OD and OS overlaid, Pearson r annotated + +Usage +----- + python scripts/output_analysis/explainability/plot_gradcam_aggregate.py \ + --agg-dir analysis_data/pipeline_nocrop/binary/single/gradcam_aggregate +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +import numpy as np +import pandas as pd + + +def load_agg(agg_dir: Path): + npz = np.load(agg_dir / "mean_heatmaps.npz") + stats = pd.read_csv(agg_dir / "attention_stats.csv") + return npz, stats + + +def _classes_from_npz(npz) -> list[str]: + classes = [] + for key in npz.files: + parts = key.split("_") + # key format: {EYE}_{ClassName}_{split}_{stat} + # ClassName may be multi-word (e.g. "Glaucoma", "Normal", "Suspect") + if parts[-1] == "mean" and parts[-2] == "all" and parts[0] == "OD": + classes.append(parts[1]) + return sorted(set(classes)) + + +def plot_mean_heatmaps(npz, classes: list[str], out_path: Path): + eyes = ["OD", "OS"] + splits = ["all", "correct", "incorrect"] + cols = [(e, s) for e in eyes for s in splits] # 6 columns + + n_rows = len(classes) + n_cols = len(cols) + fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 2.8, n_rows * 2.8)) + if n_rows == 1: + axes = axes[np.newaxis, :] + + for r, cls in enumerate(classes): + for c, (eye, split) in enumerate(cols): + ax = axes[r, c] + key = f"{eye}_{cls}_{split}_mean" + if key not in npz: + ax.axis("off") + ax.set_title(f"{eye} {split}\n(no data)", fontsize=7) + continue + cam = npz[key] + count = int(npz.get(f"{eye}_{cls}_{split}_count", np.array(0))) + ax.imshow(cam, cmap="jet", vmin=0, vmax=1) + ax.axis("off") + title = f"{cls} | {eye} {split}\n(N={count})" + ax.set_title(title, fontsize=7) + + fig.suptitle("Mean GradCAM heatmaps by class / eye / outcome", fontsize=12) + fig.tight_layout() + fig.savefig(out_path, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out_path}") + + +def plot_attention_stats(stats: pd.DataFrame, classes: list[str], out_path: Path): + eyes = ["OD", "OS"] + cmap = plt.get_cmap("tab10") + class_colors = {cls: cmap(i) for i, cls in enumerate(classes)} + + fig, axes = plt.subplots(1, 3, figsize=(16, 5)) + + # ---- Panel A: disc_frac per class × eye ---- + ax = axes[0] + positions = [] + labels = [] + data_viol = [] + tick_pos = [] + pos = 0 + for cls in classes: + for eye in eyes: + sub = stats[(stats["true_name"] == cls) & (stats["eye"] == eye)]["disc_frac"].dropna() + data_viol.append(sub.values) + positions.append(pos) + labels.append(f"{cls[:3]}\n{eye}") + tick_pos.append(pos) + pos += 1 + pos += 0.5 # gap between classes + + vp = ax.violinplot(data_viol, positions=positions, showmedians=True, widths=0.7) + for i, (pc, cls) in enumerate(zip(vp["bodies"], [c for c in classes for _ in eyes])): + pc.set_facecolor(class_colors[cls]) + pc.set_alpha(0.65) + ax.set_xticks(tick_pos) + ax.set_xticklabels(labels, fontsize=8) + ax.set_ylabel("Disc fraction (attention mass within GT disc mask)") + ax.set_title("Disc attention by class") + ax.axhline(0.5, color="black", linewidth=0.8, linestyle="--", alpha=0.4) + ax.grid(axis="y", linestyle="--", alpha=0.3) + + # ---- Panel B: entropy per class × eye (same layout) ---- + ax = axes[1] + data_ent = [] + for cls in classes: + for eye in eyes: + sub = stats[(stats["true_name"] == cls) & (stats["eye"] == eye)]["entropy"].dropna() + data_ent.append(sub.values) + + vp2 = ax.violinplot(data_ent, positions=positions, showmedians=True, widths=0.7) + for pc, cls in zip(vp2["bodies"], [c for c in classes for _ in eyes]): + pc.set_facecolor(class_colors[cls]) + pc.set_alpha(0.65) + ax.set_xticks(tick_pos) + ax.set_xticklabels(labels, fontsize=8) + ax.set_ylabel("Attention entropy (higher = more diffuse)") + ax.set_title("Attention entropy by class") + ax.grid(axis="y", linestyle="--", alpha=0.3) + + # ---- Panel C: disc_frac correct vs incorrect, mean ± std ---- + ax = axes[2] + x_ticks = [] + x_labels = [] + pos = 0 + for cls in classes: + for eye in eyes: + for split, marker, ls in [("correct", "o", "-"), ("incorrect", "X", "--")]: + sub = stats[ + (stats["true_name"] == cls) & + (stats["eye"] == eye) & + (stats["correct"] == (split == "correct")) + ]["disc_frac"].dropna() + if len(sub) == 0: + continue + ax.errorbar( + pos, sub.mean(), yerr=sub.std(), + fmt=marker, color=class_colors[cls], linestyle=ls, + capsize=4, markersize=7, alpha=0.85, + label=f"{cls[:3]} {eye} {split}" if pos < 4 else "_", + ) + pos += 1 + x_ticks.append(pos - 1.5) + x_labels.append(f"{cls[:3]}\n{eye}") + pos += 0.5 + + ax.axhline(0.5, color="black", linewidth=0.8, linestyle="--", alpha=0.4) + ax.set_ylabel("Disc fraction") + ax.set_title("Disc fraction: correct vs incorrect\n(circle=correct, X=incorrect)") + ax.grid(axis="y", linestyle="--", alpha=0.3) + + # legend: one patch per class + patches = [mpatches.Patch(color=class_colors[c], label=c) for c in classes] + patches += [ + plt.Line2D([0], [0], marker="o", color="grey", label="correct", linestyle="none"), + plt.Line2D([0], [0], marker="X", color="grey", label="incorrect", linestyle="none"), + ] + ax.legend(handles=patches, fontsize=7, loc="lower right") + + fig.suptitle("GradCAM attention statistics", fontsize=12) + fig.tight_layout() + fig.savefig(out_path, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out_path}") + + +def plot_disc_attention_correlation(stats: pd.DataFrame, classes: list[str], out_path: Path): + """ + Scatter disc_frac vs correct_conf per class. + + correct_conf = confidence if correct + = 1 - confidence if incorrect + + This asks: does focusing attention on the disc region correlate with + the model being more confident about the right answer? + """ + import scipy.stats as scipy_stats + + stats = stats.copy() + stats["correct_conf"] = np.where( + stats["correct"], + stats["confidence"], + 1.0 - stats["confidence"], + ) + + corr_colors = {True: "steelblue", False: "tomato"} + corr_labels = {True: "Correct", False: "Incorrect"} + is_binary = len(classes) == 2 + + n_cls = len(classes) + fig, axes = plt.subplots(1, n_cls, figsize=(5 * n_cls, 5), sharey=True) + if n_cls == 1: + axes = [axes] + + for ax, cls in zip(axes, classes): + sub = stats[stats["true_name"] == cls] + x_all, y_all = [], [] + + for correct_val, color in corr_colors.items(): + csub = sub[sub["correct"] == correct_val] + x = csub["disc_frac"].values + y = csub["correct_conf"].values + ax.scatter(x, y, marker="o", color=color, + alpha=0.75, s=30, + label=corr_labels[correct_val], + edgecolors="none") + x_all.extend(x.tolist()) + y_all.extend(y.tolist()) + + # pooled regression line + x_arr = np.array(x_all) + y_arr = np.array(y_all) + if len(x_arr) >= 3: + r, p = scipy_stats.pearsonr(x_arr, y_arr) + m, b = np.polyfit(x_arr, y_arr, 1) + xs = np.linspace(0, 1, 100) + ax.plot(xs, m * xs + b, color="black", linewidth=1.5, linestyle="--", alpha=0.7) + p_str = f"p={p:.3f}" if p >= 0.001 else "p<0.001" + ax.annotate(f"r={r:+.3f}\n{p_str}", xy=(0.05, 0.93), xycoords="axes fraction", + fontsize=9, va="top", + bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7)) + + if is_binary: + ax.axhline(0.5, color="red", linewidth=1.0, linestyle=":", + alpha=0.7, label="Decision boundary (0.50)") + ax.set_xlim(0, 1) + ax.set_xlabel("Disc fraction\n(attention mass within GT disc mask)", fontsize=9) + ax.set_title(cls, fontsize=11) + ax.set_ylim(-0.02, 1.05) + ax.grid(linestyle="--", alpha=0.3) + ax.legend(fontsize=8, loc="lower right") + + axes[0].set_ylabel("Correct-class confidence\n(conf if correct, 1−conf if wrong)", fontsize=9) + fig.suptitle("Disc attention vs correct-class confidence", fontsize=12) + fig.tight_layout() + fig.savefig(out_path, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out_path}") + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--agg-dir", required=True, + help="Directory produced by aggregate_gradcam.py") + ap.add_argument("--out-heatmaps", default=None) + ap.add_argument("--out-stats", default=None) + ap.add_argument("--out-corr", default=None) + args = ap.parse_args() + + agg_dir = Path(args.agg_dir) + out_hm = Path(args.out_heatmaps) if args.out_heatmaps else agg_dir / "mean_heatmaps_plot.png" + out_st = Path(args.out_stats) if args.out_stats else agg_dir / "attention_stats_plot.png" + out_corr = Path(args.out_corr) if args.out_corr else agg_dir / "disc_attention_correlation.png" + + npz, stats = load_agg(agg_dir) + classes = _classes_from_npz(npz) + print(f"Classes found: {classes}") + print(f"Total eye records in stats: {len(stats)}") + + plot_mean_heatmaps(npz, classes, out_hm) + plot_attention_stats(stats, classes, out_st) + plot_disc_attention_correlation(stats, classes, out_corr) + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/reeval_holdout_threshold.py b/scripts/output_analysis/reeval_holdout_threshold.py new file mode 100644 index 0000000..56a16d6 --- /dev/null +++ b/scripts/output_analysis/reeval_holdout_threshold.py @@ -0,0 +1,355 @@ +#!/usr/bin/env python3 +""" +Re-evaluate holdout accuracy using two threshold strategies: + + 1. acc — current behaviour: maximise raw accuracy on (imbalanced) val set + 2. youden — Youden's J = sensitivity + specificity − 1 on val set + +For each fold the val probs (already saved) supply the threshold, then the +model is re-run on the holdout set to get the actual holdout accuracy under +each strategy. + +Usage +----- + python scripts/output_analysis/reeval_holdout_threshold.py \ + --run-dir analysis_data/pipeline_10x5 \ + --eval-mode binary + + # or a single nocrop run: + python scripts/output_analysis/reeval_holdout_threshold.py \ + --run-dir analysis_data/pipeline_nocrop \ + --eval-mode binary \ + --fold-seed 42 +""" +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import pandas as pd +import torch +from sklearn.metrics import balanced_accuracy_score, roc_auc_score, roc_curve + +REPO_ROOT = Path(__file__).resolve().parents[2] +if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + +from classes.v2.metrics import tune_multiclass_bias +from classes.v2.papila_builders import build_papila_data +from classes.v2.loader_factory import filter_bilateral_samples, make_loader +from classes.v2.models import SingleEyeHT, collect_probs_single_components +from classes.v2.profiles import build_papila_profile +from classes.v2.split_manager import PatientFirstSplitManager +from classes.v2.transforms import build_eval_transform +from classes.v2.utils import choose_device + + +# --------------------------------------------------------------------------- +# Threshold helpers +# --------------------------------------------------------------------------- + +def _acc_threshold(y: np.ndarray, p1: np.ndarray) -> float: + grid = np.linspace(0.0, 1.0, 1001) + best_t, best_acc = 0.5, -1.0 + for t in grid: + acc = float(((p1 >= t).astype(int) == y).mean()) + if acc > best_acc or (acc == best_acc and abs(t - 0.5) < abs(best_t - 0.5)): + best_acc, best_t = acc, float(t) + return best_t + + +def _youden_threshold(y: np.ndarray, p1: np.ndarray) -> float: + if len(np.unique(y)) < 2: + return 0.5 + fpr, tpr, thresholds = roc_curve(y, p1) + j = tpr + (1.0 - fpr) - 1.0 + return float(thresholds[np.argmax(j)]) + + +def _apply_threshold(probs: np.ndarray, threshold: float, num_classes: int) -> np.ndarray: + if num_classes == 2: + return (probs[:, 1] >= threshold).astype(int) + # multiclass: not applicable for a single scalar threshold + return probs.argmax(axis=1) + + +# --------------------------------------------------------------------------- +# Per-fold evaluation +# --------------------------------------------------------------------------- + +def eval_fold( + fold_dir: Path, + fold_idx: int, + fold_seed: int, + eval_mode: str, + args, + device: torch.device, +) -> dict | None: + + checkpoint = fold_dir / "best_single.pt" + if not checkpoint.exists(): + print(f" [skip] {fold_dir}: no best_single.pt") + return None + + val_y_path = fold_dir / "y_true.npy" + val_p_path = fold_dir / "probs_fused.npy" + if not val_y_path.exists() or not val_p_path.exists(): + print(f" [skip] {fold_dir}: no val probs") + return None + + val_y = np.load(val_y_path) + val_p = np.load(val_p_path) + num_classes = val_p.shape[1] + + # ---- decision boundaries from val ---- + if num_classes == 2: + t_acc = _acc_threshold(val_y, val_p[:, 1]) + t_youden = _youden_threshold(val_y, val_p[:, 1]) + else: + # multiclass: compare raw-acc-optimised bias vs balanced-acc-optimised bias + # raw-acc bias: temporarily swap objective back to raw accuracy + from sklearn.metrics import accuracy_score + import copy + + def _tune_bias_raw(y, p): + c = p.shape[1] + bias = np.zeros(c) + grid = np.linspace(-1.0, 1.0, 41) + for _ in range(2): + for k in range(c): + best_v, best_acc = bias[k], -1.0 + old = bias[k] + for v in grid: + bias[k] = float(v) + logits = np.log(np.clip(p, 1e-8, 1.0)) + bias.reshape(1, -1) + acc = float((logits.argmax(1) == y).mean()) + if acc > best_acc or (acc == best_acc and abs(v) < abs(best_v)): + best_acc, best_v = acc, float(v) + bias[k] = best_v + return bias + + bias_raw = _tune_bias_raw(val_y, val_p) + bias_bal = tune_multiclass_bias(val_y, val_p) # balanced acc objective + + # ---- reconstruct holdout split ---- + data = build_papila_data( + image_dir=args.image_dir, + clinical_dir=args.clinical_dir, + label_col=args.label_col, + cat_cols=args.cat_cols, + n_splits=args.n_splits, + random_seed=fold_seed, + iop_corr_method=getattr(args, "iop_corr_method", "ratio"), + ) + df_mode = data.df.copy() + if eval_mode == "binary": + df_mode = df_mode[df_mode[args.label_col].isin([0, 1])].reset_index(drop=True) + + splitter = PatientFirstSplitManager( + patient_col="Patient ID", label_col=args.label_col + ) + split_args = SimpleNamespace( + eval_mode=eval_mode, + holdout_per_class=args.holdout_per_class, + holdout_seed=args.holdout_seed, + n_splits=args.n_splits, + fold_seed=fold_seed, + ) + plans = splitter.build_plans( + clinical=SimpleNamespace(df=df_mode, label_col=args.label_col), + args=split_args, + profile=None, + ) + split = plans[fold_idx] + + if split.holdout is None or split.holdout.empty: + print(f" [skip] {fold_dir}: no holdout") + return None + + # ---- build holdout loader ---- + profile_patient = build_papila_profile( + patient_col="Patient ID", label_col=args.label_col, sample_mode="patient" + ) + holdout_samples = filter_bilateral_samples( + profile_patient.build_samples(df=split.holdout, clinical=data) + ) + if not holdout_samples: + print(f" [skip] {fold_dir}: no bilateral holdout samples") + return None + + eval_transform = build_eval_transform(args.backbone) + holdout_loader = make_loader( + holdout_samples, + profile_patient.slot_descriptors(), + image_transform=eval_transform, + image_preprocessor=None, + batch_size=args.batch_size, + shuffle=False, + num_workers=args.num_workers, + ) + + # ---- load model ---- + model = SingleEyeHT( + backbone=args.backbone, + freeze_ratio=0.0, + augment=False, + clinical_data=data, + num_classes=num_classes, + md_hidden_dim=getattr(args, "md_hidden_dim", 64), + fusion_dim=getattr(args, "fusion_dim", 128), + bridge_mode=getattr(args, "bridge_mode", "fused"), + ).to(device) + model.load_state_dict( + torch.load(checkpoint, map_location=device, weights_only=False) + ) + model.eval() + + # ---- run inference ---- + hld_y, hld_p, _, _ = collect_probs_single_components( + model, holdout_loader, device, aggregate_patient=True + ) + + if len(hld_y) == 0: + return None + + hld_auc = float(roc_auc_score( + hld_y, hld_p[:, 1] if num_classes == 2 else hld_p, + multi_class="ovr" if num_classes > 2 else "raise", + )) + + if num_classes == 2: + acc_old = float(((hld_p[:, 1] >= t_acc).astype(int) == hld_y).mean()) + acc_new = float(((hld_p[:, 1] >= t_youden).astype(int) == hld_y).mean()) + bacc_old = balanced_accuracy_score(hld_y, (hld_p[:, 1] >= t_acc).astype(int)) + bacc_new = balanced_accuracy_score(hld_y, (hld_p[:, 1] >= t_youden).astype(int)) + row_extra = {"t_old": t_acc, "t_new": t_youden} + else: + logits_raw = np.log(np.clip(hld_p, 1e-8, 1.0)) + bias_raw.reshape(1, -1) + logits_bal = np.log(np.clip(hld_p, 1e-8, 1.0)) + bias_bal.reshape(1, -1) + preds_raw = logits_raw.argmax(1) + preds_bal = logits_bal.argmax(1) + acc_old = float((preds_raw == hld_y).mean()) + acc_new = float((preds_bal == hld_y).mean()) + bacc_old = balanced_accuracy_score(hld_y, preds_raw) + bacc_new = balanced_accuracy_score(hld_y, preds_bal) + row_extra = {"bias_raw": bias_raw.tolist(), "bias_bal": bias_bal.tolist()} + + return { + "fold_dir": str(fold_dir), + "fold": fold_idx, + "fold_seed": fold_seed, + "hld_auc": hld_auc, + "hld_acc_old": acc_old, + "hld_acc_new": acc_new, + "hld_bacc_old": bacc_old, + "hld_bacc_new": bacc_new, + "n_holdout": len(hld_y), + **row_extra, + } + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", required=True, + help="e.g. analysis_data/pipeline_10x5 or analysis_data/pipeline_nocrop") + ap.add_argument("--eval-mode", default="binary", choices=["binary", "multiclass"]) + ap.add_argument("--fold-seed", type=int, default=None, + help="Override fold seed (for single-rep runs). " + "For 10x5, seeds are inferred from rep dir name.") + ap.add_argument("--backbone", default="refugelike") + 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("--label-col", default="Diagnosis") + ap.add_argument("--cat-cols", nargs="*", default=["Gender"]) + ap.add_argument("--image-dir", default="Papila/FundusImages") + ap.add_argument("--clinical-dir",default="Papila/ClinicalData") + ap.add_argument("--iop-corr-method", default="ratio") + ap.add_argument("--batch-size", type=int, default=8) + ap.add_argument("--num-workers", type=int, default=4) + ap.add_argument("--md-hidden-dim", type=int, default=128) + ap.add_argument("--fusion-dim", type=int, default=256) + ap.add_argument("--bridge-mode", default="fused") + ap.add_argument("--device", default="auto") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + device = choose_device(args.device) + run_dir = Path(args.run_dir) + out_path = Path(args.out) if args.out else \ + run_dir / f"reeval_threshold_{args.eval_mode}.csv" + + # Discover fold dirs — supports both flat (fold0..fold4) and + # rep-based (rep00/binary/ensemble/fold0) layouts + _BASE_SEED = 100 + _SEED_STRIDE = 100 + + fold_jobs: list[tuple[Path, int, int]] = [] # (fold_dir, fold_idx, fold_seed) + + rep_dirs = sorted(run_dir.glob("rep[0-9]*")) + if rep_dirs: + for rep_dir in rep_dirs: + rep_n = int(rep_dir.name.replace("rep", "")) + fold_seed = _BASE_SEED + rep_n * _SEED_STRIDE + mode_dir = rep_dir / args.eval_mode / "ensemble" + if not mode_dir.exists(): + continue + for fd in sorted(mode_dir.glob("fold[0-9]*"), key=lambda p: int(p.name[4:])): + fold_jobs.append((fd, int(fd.name[4:]), fold_seed)) + else: + # flat layout + mode_dir = run_dir / args.eval_mode / "ensemble" + fold_seed = args.fold_seed if args.fold_seed is not None else 42 + for fd in sorted(mode_dir.glob("fold[0-9]*"), key=lambda p: int(p.name[4:])): + fold_jobs.append((fd, int(fd.name[4:]), fold_seed)) + + if not fold_jobs: + sys.exit(f"No fold directories found under {run_dir}") + + print(f"Found {len(fold_jobs)} folds to re-evaluate") + + rows = [] + for i, (fold_dir, fold_idx, fold_seed) in enumerate(fold_jobs): + print(f"\n[{i+1}/{len(fold_jobs)}] {fold_dir} fold_seed={fold_seed}") + row = eval_fold(fold_dir, fold_idx, fold_seed, args.eval_mode, args, device) + if row: + rows.append(row) + print(f" hld_acc(old)={row['hld_acc_old']:.3f} " + f"hld_acc(new)={row['hld_acc_new']:.3f} " + f"hld_bacc(old)={row['hld_bacc_old']:.3f} " + f"hld_bacc(new)={row['hld_bacc_new']:.3f} " + f"hld_auc={row['hld_auc']:.3f}") + + if not rows: + print("No results.") + return + + df = pd.DataFrame(rows) + df.to_csv(out_path, index=False) + print(f"\nSaved → {out_path}") + is_binary = args.eval_mode == "binary" + old_label = "acc-threshold" if is_binary else "raw-acc bias" + new_label = "Youden-J" if is_binary else "balanced-acc bias" + + print(f"\n{'='*60}") + print(f"Summary ({args.eval_mode}, n={len(df)} folds)") + print(f"{'='*60}") + print(f" Holdout AUC: {df.hld_auc.mean():.4f} ± {df.hld_auc.std():.4f}") + print(f" Holdout acc ({old_label:<18}): {df.hld_acc_old.mean():.4f} ± {df.hld_acc_old.std():.4f}") + print(f" Holdout acc ({new_label:<18}): {df.hld_acc_new.mean():.4f} ± {df.hld_acc_new.std():.4f}") + print(f" Holdout bacc ({old_label:<18}): {df.hld_bacc_old.mean():.4f} ± {df.hld_bacc_old.std():.4f}") + print(f" Holdout bacc ({new_label:<18}): {df.hld_bacc_new.mean():.4f} ± {df.hld_bacc_new.std():.4f}") + print(f" Delta acc (new − old): {df.hld_acc_new.mean() - df.hld_acc_old.mean():+.4f}") + print(f" Delta bacc (new − old): {df.hld_bacc_new.mean() - df.hld_bacc_old.mean():+.4f}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/build_10x5_aggregate_panel.py b/scripts/output_analysis/visualizations/build_10x5_aggregate_panel.py new file mode 100644 index 0000000..f16af5c --- /dev/null +++ b/scripts/output_analysis/visualizations/build_10x5_aggregate_panel.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +""" +Assemble the 10×5 aggregate comparison panel. + +Layout (3 rows × 2 cols): + col 0 = binary, col 1 = multiclass + + row 0: mean ROC curve + row 1: rep stability + row 2: holdout stability + +Usage +----- + python scripts/output_analysis/visualizations/build_10x5_aggregate_panel.py \ + --run-dir analysis_data/pipeline_10x5 \ + --out analysis_data/pipeline_10x5/aggregate/aggregate_panel.png +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.image as mpimg + + +# (row, col, rel_path, label) +CELLS = [ + (0, 0, "aggregate/binary_roc_mean.png", "Binary — Mean ROC"), + (0, 1, "aggregate/multiclass_roc_mean.png", "Multiclass — Mean ROC"), + (1, 0, "aggregate/binary_rep_stability.png", "Binary — Rep stability"), + (1, 1, "aggregate/multiclass_rep_stability.png", "Multiclass — Rep stability"), + (2, 0, "aggregate/binary_holdout_stability.png", "Binary — Holdout stability"), + (2, 1, "aggregate/multiclass_holdout_stability.png", "Multiclass — Holdout stability"), +] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_10x5") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + run_dir = Path(args.run_dir) + out = Path(args.out) if args.out else run_dir / "aggregate" / "aggregate_panel.png" + + fig = plt.figure(figsize=(16, 18)) + gs = fig.add_gridspec(3, 2, hspace=0.06, wspace=0.04) + + for row, col, rel, label in CELLS: + ax = fig.add_subplot(gs[row, col]) + img = mpimg.imread(str(run_dir / rel)) + ax.imshow(img) + ax.axis("off") + ax.set_title(label, fontsize=11, pad=5) + + out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/build_10x5_training_panel.py b/scripts/output_analysis/visualizations/build_10x5_training_panel.py new file mode 100644 index 0000000..26aa037 --- /dev/null +++ b/scripts/output_analysis/visualizations/build_10x5_training_panel.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +""" +Assemble the 10×5 training dynamics panel. + +Layout: 3-row × 4-col gridspec with spanning + + Top 2×2 (each cell spans 2 cols): + row 0, cols 0-1: binary early stopping sweep + row 0, cols 2-3: multiclass early stopping sweep + row 1, cols 0-1: binary cost of stopping early + row 1, cols 2-3: multiclass cost of stopping early + + Bottom row of 4 (one col each): + row 2, col 0: binary holdout AUC by epoch + row 2, col 1: binary val−holdout gap (val-adjusted) + row 2, col 2: multiclass holdout AUC by epoch + row 2, col 3: multiclass val−holdout gap (val-adjusted) + +Usage +----- + python scripts/output_analysis/visualizations/build_10x5_training_panel.py \ + --run-dir analysis_data/pipeline_10x5 \ + --out analysis_data/pipeline_10x5/aggregate/training_dynamics_panel.png +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.image as mpimg + + +# (row, col_start, col_end, rel_path, label) +# col_end is exclusive slice — use None for single cell +CELLS = [ + # ---- top 2×2: early stopping (each spans 2 cols) ---- + (0, 0, 2, "binary/ensemble/plots/early_stopping_sweep_fused.png", + "Binary — Early stopping sweep"), + (0, 2, 4, "multiclass/ensemble/plots/early_stopping_sweep_fused.png", + "Multiclass — Early stopping sweep"), + (1, 0, 2, "binary/ensemble/plots/early_stopping_sweep_fused_inverted.png", + "Binary — Cost of stopping early"), + (1, 2, 4, "multiclass/ensemble/plots/early_stopping_sweep_fused_inverted_tol0.002.png", + "Multiclass — Cost of stopping early (CI tol=0.002)"), + # ---- bottom row of 4: holdout epoch curves (single col each) ---- + (2, 0, 1, "binary/ensemble/plots/holdout_epoch_curves_fused.png", + "Binary — Holdout AUC by epoch"), + (2, 1, 2, "binary/ensemble/plots/holdout_epoch_curves_fused_delta_adj.png", + "Binary — Val−Holdout gap (val-adjusted)"), + (2, 2, 3, "multiclass/ensemble/plots/holdout_epoch_curves_fused.png", + "Multiclass — Holdout AUC by epoch"), + (2, 3, 4, "multiclass/ensemble/plots/holdout_epoch_curves_fused_delta_adj.png", + "Multiclass — Val−Holdout gap (val-adjusted)"), +] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_10x5") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + run_dir = Path(args.run_dir) + out = Path(args.out) if args.out else run_dir / "aggregate" / "training_dynamics_panel.png" + + fig = plt.figure(figsize=(22, 16)) + gs = fig.add_gridspec(3, 4, height_ratios=[1, 1, 0.75], hspace=0.08, wspace=0.04) + + for row, col_start, col_end, rel, label in CELLS: + ax = fig.add_subplot(gs[row, col_start:col_end]) + img = mpimg.imread(str(run_dir / rel)) + ax.imshow(img) + ax.axis("off") + ax.set_title(label, fontsize=11, pad=5) + + out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/build_fusion_head_panel.py b/scripts/output_analysis/visualizations/build_fusion_head_panel.py new file mode 100644 index 0000000..d2c0a54 --- /dev/null +++ b/scripts/output_analysis/visualizations/build_fusion_head_panel.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 +""" +Assemble a 2x2 fusion-head comparison panel from pipeline_nocrop. + +Layout: + [binary ensemble ROC] [binary fusion-head explainability summary] + [multiclass ensemble ROC][multiclass fusion-head explainability summary] + +Usage +----- + python scripts/output_analysis/visualizations/build_fusion_head_panel.py \ + --run-dir analysis_data/pipeline_nocrop \ + --out analysis_data/pipeline_nocrop/fusion_head_comparison_panel.png +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.image as mpimg +import numpy as np + + +ROC_CELLS = [ + # (row, col, rel_path, label) + (0, 0, "binary/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Binary — Ensemble"), + (0, 1, "binary/ensemble/plots/roc_probs_fused_head_mean_ovr.png", + "Binary — Fusion Head"), + (1, 0, "multiclass/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Multiclass — Ensemble"), + (1, 1, "multiclass/ensemble/plots/roc_probs_fused_head_mean_ovr.png", + "Multiclass — Fusion Head"), +] + +EXPL_ROWS = [ + # (row_in_grid, rel_path, label) + (2, "binary/ensemble/explainability_fusion_summary_val_fused_head.png", + "Binary — Fusion Head events (val)"), + (3, "multiclass/ensemble/explainability_fusion_summary_val_fused_head.png", + "Multiclass — Fusion Head events (val)"), +] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_nocrop") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + run_dir = Path(args.run_dir) + out = Path(args.out) if args.out else run_dir / "fusion_head_comparison_panel.png" + + # load all images + roc_imgs = {(r, c): (mpimg.imread(str(run_dir / rel)), lbl) + for r, c, rel, lbl in ROC_CELLS} + expl_imgs = [(mpimg.imread(str(run_dir / rel)), lbl) + for _, rel, lbl in EXPL_ROWS] + + # 4-row grid: rows 0-1 are the 2×2 ROC square; rows 2-3 are full-width explainability + fig = plt.figure(figsize=(14, 20)) + gs = fig.add_gridspec( + 4, 2, + height_ratios=[1, 1, 0.6, 0.6], + hspace=0.06, + wspace=0.04, + ) + + # ROC cells (2×2) + for row, col, _, _ in ROC_CELLS: + ax = fig.add_subplot(gs[row, col]) + img, label = roc_imgs[(row, col)] + ax.imshow(img) + ax.axis("off") + ax.set_title(label, fontsize=11, pad=5) + + # Explainability rows (span both columns) + for grid_row, (img, label) in zip([2, 3], expl_imgs): + ax = fig.add_subplot(gs[grid_row, :]) + ax.imshow(img) + ax.axis("off") + ax.set_title(label, fontsize=11, pad=5) + + out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/build_md_explainability_panel.py b/scripts/output_analysis/visualizations/build_md_explainability_panel.py new file mode 100644 index 0000000..267a191 --- /dev/null +++ b/scripts/output_analysis/visualizations/build_md_explainability_panel.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +""" +Assemble a metadata-explainability comparison panel. + +Layout (2 rows × 4 cols): + row 0 = binary, row 1 = multiclass + + col 0: single md_importance + col 1: ensemble md_importance + col 2: nocrop ROC + col 3: excl_phakic_axial ROC + +Usage +----- + python scripts/output_analysis/visualizations/build_md_explainability_panel.py \ + --nocrop-dir analysis_data/pipeline_nocrop \ + --excl-dir analysis_data/pipeline_nocrop_excl_phakic_axial \ + --out analysis_data/pipeline_nocrop/md_explainability_panel.png +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.image as mpimg + + +# (row, col, dir_key, rel_path, label) +# dir_key: "nocrop" or "excl" +CELLS = [ + # ---- binary row (row 0) ---- + (0, 0, "nocrop", "binary/single/explainability_md_importance_summary.png", + "Binary — Single"), + (0, 1, "nocrop", "binary/ensemble/explainability_md_importance_summary.png", + "Binary — Ensemble"), + (0, 2, "nocrop", "binary/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Binary — nocrop ROC"), + (0, 3, "excl", "binary/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Binary — excl phakic+axial ROC"), + # ---- multiclass row (row 1) ---- + (1, 0, "nocrop", "multiclass/single/explainability_md_importance_summary.png", + "Multiclass — Single"), + (1, 1, "nocrop", "multiclass/ensemble/explainability_md_importance_summary.png", + "Multiclass — Ensemble"), + (1, 2, "nocrop", "multiclass/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Multiclass — nocrop ROC"), + (1, 3, "excl", "multiclass/ensemble/plots/roc_probs_fused_mean_ovr.png", + "Multiclass — excl phakic+axial ROC"), +] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--nocrop-dir", default="analysis_data/pipeline_nocrop") + ap.add_argument("--excl-dir", default="analysis_data/pipeline_nocrop_excl_phakic_axial") + ap.add_argument("--out", default=None) + args = ap.parse_args() + + nocrop_dir = Path(args.nocrop_dir) + excl_dir = Path(args.excl_dir) + out = Path(args.out) if args.out else nocrop_dir / "md_explainability_panel.png" + + fig = plt.figure(figsize=(24, 12)) + gs = fig.add_gridspec( + 2, 4, + hspace=0.08, + wspace=0.04, + ) + + dirs = {"nocrop": nocrop_dir, "excl": excl_dir} + + for row, col, dir_key, rel, label in CELLS: + ax = fig.add_subplot(gs[row, col]) + img = mpimg.imread(str(dirs[dir_key] / rel)) + ax.imshow(img) + ax.axis("off") + ax.set_title(label, fontsize=11, pad=5) + + out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved → {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/plot_early_stopping_sweep.py b/scripts/output_analysis/visualizations/plot_early_stopping_sweep.py new file mode 100644 index 0000000..0b708ce --- /dev/null +++ b/scripts/output_analysis/visualizations/plot_early_stopping_sweep.py @@ -0,0 +1,217 @@ +#!/usr/bin/env python3 +""" +Simulate early stopping at each epoch N and show what val/holdout AUC +you would have gotten if you stopped there. + +For each fold and each candidate stopping epoch N: + - Find the epoch <= N with the highest val AUC (checkpoint selection) + - Record the val AUC and holdout AUC at that epoch + +Then plot mean ± std across all folds as a function of N. + +Usage +----- + python scripts/output_analysis/visualizations/plot_early_stopping_sweep.py \ + --run-dir analysis_data/pipeline_10x5 \ + --eval-mode binary \ + --tower-mode ensemble \ + --head fused +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +HEAD_COL = { + "fused": ("ensemble_val_auc", "ensemble_holdout_auc"), + "img": ("ensemble_val_auc_img", "ensemble_holdout_auc_img"), + "md": ("ensemble_val_auc_md", "ensemble_holdout_auc_md"), + "classic": ("classic_val_auc", "classic_holdout_auc"), +} + + +def load_fold_logs(run_dir: Path, eval_mode: str, tower_mode: str, + val_col: str, hld_col: str): + logs = [] + for rep_dir in sorted(run_dir.glob("rep*")): + mode_dir = rep_dir / eval_mode / tower_mode + if not mode_dir.exists(): + continue + for fd in sorted( + [d for d in mode_dir.iterdir() if d.is_dir() and d.name.startswith("fold")], + key=lambda p: int(p.name.replace("fold", "")), + ): + log = fd / "epoch_log.csv" + if not log.exists(): + continue + df = pd.read_csv(log) + if val_col not in df.columns or hld_col not in df.columns: + continue + df = df[["epoch", val_col, hld_col]].dropna() + logs.append(df.reset_index(drop=True)) + return logs + + +def sweep(logs: list[pd.DataFrame], val_col: str, hld_col: str): + max_epoch = max(df["epoch"].max() for df in logs) + epochs = np.arange(1, int(max_epoch) + 1) + + val_mat = np.full((len(logs), len(epochs)), np.nan) + hld_mat = np.full((len(logs), len(epochs)), np.nan) + + for i, df in enumerate(logs): + for j, n in enumerate(epochs): + window = df[df["epoch"] <= n] + if window.empty: + continue + best_idx = window[val_col].idxmax() + val_mat[i, j] = window.loc[best_idx, val_col] + hld_mat[i, j] = window.loc[best_idx, hld_col] + + return epochs, val_mat, hld_mat + + +def plot(epochs, val_mat, hld_mat, out_path: Path, title: str, inverted: bool = False, ci_tol: float = 0.0): + val_mean = np.nanmean(val_mat, axis=0) + val_std = np.nanstd(val_mat, axis=0) + hld_mean = np.nanmean(hld_mat, axis=0) + hld_std = np.nanstd(hld_mat, axis=0) + + if inverted: + # compute cost per fold, then aggregate — avoids max-of-mean bias + best_val_per_fold = np.nanmax(val_mat, axis=1, keepdims=True) # (n_folds, 1) + best_hld_per_fold = np.nanmax(hld_mat, axis=1, keepdims=True) + delta_val = best_val_per_fold - val_mat # (n_folds, n_epochs) + delta_hld = best_hld_per_fold - hld_mat + y_val = np.nanmean(delta_val, axis=0) + y_hld = np.nanmean(delta_hld, axis=0) + sy_val = np.nanstd(delta_val, axis=0) + sy_hld = np.nanstd(delta_hld, axis=0) + else: + y_val, y_hld = val_mean, hld_mean + sy_val, sy_hld = val_std, hld_std + + fig, ax = plt.subplots(figsize=(11, 5)) + + if inverted: + # faint per-fold lines + for i in range(delta_val.shape[0]): + ax.plot(epochs, delta_val[i], color="steelblue", linewidth=0.6, alpha=0.18) + ax.plot(epochs, delta_hld[i], color="firebrick", linewidth=0.6, alpha=0.18) + + ax.plot(epochs, y_val, color="steelblue", linewidth=2.0, + label="Best val − val@N (val cost of stopping early)" if inverted + else "Val AUC (best ckpt up to N)") + ax.fill_between(epochs, y_val - sy_val, y_val + sy_val, color="steelblue", alpha=0.15) + + ax.plot(epochs, y_hld, color="firebrick", linewidth=2.0, + label="Best hld − hld@N (holdout cost of stopping early)" if inverted + else "Holdout AUC (at best val ckpt)") + ax.fill_between(epochs, y_hld - sy_hld, y_hld + sy_hld, color="firebrick", alpha=0.15) + + if inverted: + ax.axhline(0, color="black", linewidth=1.0, linestyle="--", alpha=0.4) + + # CI-crosses-zero regions (with optional tolerance) + val_ci_zero = (y_val - sy_val) <= ci_tol + hld_ci_zero = (y_hld - sy_hld) <= ci_tol + both_ci_zero = val_ci_zero & hld_ci_zero + + ymax = max(np.nanmax(y_val), np.nanmax(y_hld)) * 1.15 + ax.fill_between(epochs, 0, ymax, where=val_ci_zero, + color="steelblue", alpha=0.12, label="val CI ≤ 0") + ax.fill_between(epochs, 0, ymax, where=hld_ci_zero, + color="firebrick", alpha=0.12, label="holdout CI ≤ 0") + ax.fill_between(epochs, 0, ymax, where=both_ci_zero, + color="purple", alpha=0.20, label="both CI ≤ 0") + + ax.set_ylabel("AUC lost vs best achievable") + ax.set_ylim(-0.05, ymax) + legend_loc = "upper right" + else: + # gap curve on twin axis + gap_mean = val_mean - hld_mean + ax2 = ax.twinx() + ax2.plot(epochs, gap_mean, color="darkorange", linewidth=1.5, + linestyle="--", alpha=0.7, label="Val−Holdout gap") + ax2.set_ylabel("Val − Holdout gap", color="darkorange", fontsize=9) + ax2.tick_params(axis="y", labelcolor="darkorange") + ax2.set_ylim(-0.1, 0.4) + lines2, labels2 = ax2.get_legend_handles_labels() + + best_hld_ep = epochs[np.nanargmax(hld_mean)] + best_hld_val = hld_mean[np.nanargmax(hld_mean)] + ax.axvline(best_hld_ep, color="firebrick", linewidth=1.2, linestyle=":", + alpha=0.8, label=f"peak holdout @ epoch {best_hld_ep} ({best_hld_val:.3f})") + stable = epochs >= 3 + min_gap_ep = epochs[stable][np.nanargmin(gap_mean[stable])] + ax.axvline(min_gap_ep, color="darkorange", linewidth=1.2, linestyle=":", + alpha=0.8, label=f"min gap @ epoch {min_gap_ep}") + ax.set_ylabel("AUC") + ax.set_ylim(0.5, 1.05) + legend_loc = "lower right" + + ax.set_xlabel("Stopping epoch N") + ax.set_title(title) + ax.grid(axis="y", linestyle="--", alpha=0.35) + + lines1, labels1 = ax.get_legend_handles_labels() + if not inverted: + lines1 += lines2; labels1 += labels2 + ax.legend(lines1, labels1, fontsize=8, loc=legend_loc) + + out_path.parent.mkdir(parents=True, exist_ok=True) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + plt.close(fig) + print(f"Saved → {out_path}") + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_10x5") + ap.add_argument("--eval-mode", default="binary") + ap.add_argument("--tower-mode", default="ensemble") + ap.add_argument("--head", default="fused", choices=list(HEAD_COL)) + ap.add_argument("--out", default=None) + ap.add_argument("--inverted", action="store_true", + help="Plot best-achievable minus current (cost of stopping early)") + ap.add_argument("--ci-tol", type=float, default=0.0, + help="Tolerance for CI-crosses-zero shading (default 0.0)") + args = ap.parse_args() + + run_dir = Path(args.run_dir) + val_col, hld_col = HEAD_COL[args.head] + + logs = load_fold_logs(run_dir, args.eval_mode, args.tower_mode, val_col, hld_col) + if not logs: + print("No epoch_log.csv files found.") + return + print(f"Loaded {len(logs)} fold logs") + + epochs, val_mat, hld_mat = sweep(logs, val_col, hld_col) + + tol_tag = f"_tol{args.ci_tol}" if args.ci_tol else "" + suffix = f"_inverted{tol_tag}" if args.inverted else "" + out = Path(args.out) if args.out else ( + run_dir / args.eval_mode / args.tower_mode / "plots" / + f"early_stopping_sweep_{args.head}{suffix}.png" + ) + title = ("Simulated early stopping — cost of stopping at epoch N\n" + if args.inverted else + "Simulated early stopping sweep\n") + title += f"{run_dir.name} · {args.eval_mode}/{args.tower_mode} · head={args.head}" + if args.inverted and args.ci_tol: + title += f" (CI tol={args.ci_tol})" + plot(epochs, val_mat, hld_mat, out, title, inverted=args.inverted, ci_tol=args.ci_tol) + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/plot_holdout_epoch_curves.py b/scripts/output_analysis/visualizations/plot_holdout_epoch_curves.py new file mode 100644 index 0000000..a6d737c --- /dev/null +++ b/scripts/output_analysis/visualizations/plot_holdout_epoch_curves.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +""" +Plot per-epoch holdout metrics across all folds in a 10x5 (or any multi-rep) run. + +Each fold gets its own line. Lines are coloured by rep. + +Modes +----- + holdout — raw holdout AUC per epoch (original plot) + delta — val_auc - holdout_auc per epoch (generalization gap; + closer to 0 = val most faithfully reflects holdout) + +Usage +----- + python scripts/output_analysis/visualizations/plot_holdout_epoch_curves.py \ + --run-dir analysis_data/pipeline_10x5 \ + --eval-mode binary \ + --tower-mode ensemble \ + --head fused \ + --mode delta +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +HEAD_COL = { + "fused": ("ensemble_val_auc", "ensemble_holdout_auc"), + "img": ("ensemble_val_auc_img", "ensemble_holdout_auc_img"), + "md": ("ensemble_val_auc_md", "ensemble_holdout_auc_md"), + "classic": ("classic_val_auc", "classic_holdout_auc"), +} + + +def load_curves(run_dir: Path, eval_mode: str, tower_mode: str, + val_col: str, hld_col: str, mode: str): + """ + Returns list of (rep, fold, epochs_array, values_array). + mode='holdout' → values = holdout_auc + mode='delta' → values = val_auc - holdout_auc + """ + curves = [] + for rep_dir in sorted(run_dir.glob("rep*")): + mode_dir = rep_dir / eval_mode / tower_mode + if not mode_dir.exists(): + continue + fold_dirs = sorted( + [d for d in mode_dir.iterdir() if d.is_dir() and d.name.startswith("fold")], + key=lambda p: int(p.name.replace("fold", "")), + ) + for fd in fold_dirs: + log = fd / "epoch_log.csv" + if not log.exists(): + continue + df = pd.read_csv(log) + needed = [hld_col] if mode == "holdout" else [val_col, hld_col] + if any(c not in df.columns for c in needed): + continue + df = df.dropna(subset=needed) + if mode == "holdout": + values = df[hld_col].to_numpy() + elif mode == "delta": + values = (df[val_col] - df[hld_col]).to_numpy() + else: # delta_adj + val_arr = df[val_col].to_numpy() + hld_arr = df[hld_col].to_numpy() + val_best = np.nanmax(val_arr) + penalty = val_best - val_arr # 0 when val is at its peak + values = (val_arr - hld_arr) + penalty + curves.append((rep_dir.name, fd.name, df["epoch"].to_numpy(), values)) + return curves + + +def _build_mean_matrix(curves): + all_ep = max(len(e) for _, _, e, _ in curves) + mat = np.full((len(curves), all_ep), np.nan) + for i, (_, _, e, a) in enumerate(curves): + mat[i, :len(a)] = a + return mat, np.arange(1, all_ep + 1) + + +def plot(curves, mode: str, out_path: Path, title: str): + reps = sorted(set(r for r, _, _, _ in curves)) + cmap = matplotlib.colormaps.get_cmap("tab10") + rep_color = {r: cmap(i / max(len(reps) - 1, 1)) for i, r in enumerate(reps)} + + fig, ax = plt.subplots(figsize=(12, 6)) + + for rep, fold, epochs, vals in curves: + ax.plot(epochs, vals, color=rep_color[rep], alpha=0.3, linewidth=0.9) + + # per-rep mean + for rep in reps: + rep_curves = [(e, a) for r, _, e, a in curves if r == rep] + max_ep = max(len(e) for e, _ in rep_curves) + mat = np.full((len(rep_curves), max_ep), np.nan) + for i, (e, a) in enumerate(rep_curves): + mat[i, :len(a)] = a + mean_curve = np.nanmean(mat, axis=0) + ax.plot(np.arange(1, max_ep + 1), mean_curve, + color=rep_color[rep], linewidth=1.8, alpha=0.85, label=rep) + + # global mean ± std + all_mat, ep_axis = _build_mean_matrix(curves) + global_mean = np.nanmean(all_mat, axis=0) + global_std = np.nanstd(all_mat, axis=0) + ax.plot(ep_axis, global_mean, color="black", linewidth=2.5, zorder=5, label="global mean") + ax.fill_between(ep_axis, global_mean - global_std, global_mean + global_std, + color="black", alpha=0.12, zorder=4) + + if mode in ("delta", "delta_adj"): + ax.axhline(0, color="black", linewidth=1.0, linestyle="--", alpha=0.5) + if mode == "delta": + ax.set_ylabel("Val AUC − Holdout AUC (gap)") + else: + ax.set_ylabel("(Val − Holdout) + (ValBest − Val) (adjusted gap)") + min_ep = int(ep_axis[np.nanargmin(global_mean)]) + min_val = global_mean[np.nanargmin(global_mean)] + ax.axvline(min_ep, color="red", linewidth=1.2, linestyle=":", alpha=0.7, + label=f"min adjusted gap @ epoch {min_ep} ({min_val:+.3f})") + else: + ax.set_ylabel("Holdout AUC") + ax.set_ylim(0, 1.05) + + ax.set_xlabel("Epoch") + ax.set_title(title) + ax.legend(fontsize=7, ncol=2, loc="upper right" if mode == "delta" else "lower right") + ax.grid(axis="y", linestyle="--", alpha=0.4) + + out_path.parent.mkdir(parents=True, exist_ok=True) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + plt.close(fig) + print(f"Saved → {out_path}") + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--run-dir", default="analysis_data/pipeline_10x5") + ap.add_argument("--eval-mode", default="binary") + ap.add_argument("--tower-mode", default="ensemble") + ap.add_argument("--head", default="fused", choices=list(HEAD_COL)) + ap.add_argument("--mode", default="holdout", choices=["holdout", "delta", "delta_adj"]) + ap.add_argument("--out", default=None) + args = ap.parse_args() + + run_dir = Path(args.run_dir) + val_col, hld_col = HEAD_COL[args.head] + curves = load_curves(run_dir, args.eval_mode, args.tower_mode, + val_col, hld_col, args.mode) + + if not curves: + print("No epoch_log.csv files found — check --run-dir / --eval-mode / --tower-mode") + return + + print(f"Loaded {len(curves)} fold curves, up to {max(len(e) for _,_,e,_ in curves)} epochs each") + + out = Path(args.out) if args.out else ( + run_dir / args.eval_mode / args.tower_mode / "plots" / + f"holdout_epoch_curves_{args.head}_{args.mode}.png" + ) + label = {"holdout": "holdout AUC", "delta": "val−holdout gap", "delta_adj": "val−holdout gap (val-adjusted)"}[args.mode] + title = (f"Per-fold {label} by epoch\n" + f"{run_dir.name} · {args.eval_mode}/{args.tower_mode} · head={args.head}") + plot(curves, args.mode, out, title) + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/plot_mdonly_curves.py b/scripts/output_analysis/visualizations/plot_mdonly_curves.py new file mode 100644 index 0000000..ef66034 --- /dev/null +++ b/scripts/output_analysis/visualizations/plot_mdonly_curves.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +""" +Per-epoch learning curves for mdonly runs. + +Reads epoch_log.csv from each fold dir and plots val_auc, val_acc, +hld_auc, hld_acc — one figure per metric, all folds as individual lines. + +Usage: + python scripts/output_analysis/visualizations/plot_mdonly_curves.py \ + --run-dirs analysis_data/pipeline_mdonly_50ep \ + analysis_data/pipeline_mdonly_200ep \ + analysis_data/pipeline_mdonly_500ep \ + --eval-mode binary +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +METRICS = [ + ("val_auc", "Val AUC"), + ("val_acc", "Val Accuracy"), + ("hld_auc", "Holdout AUC"), + ("hld_acc", "Holdout Accuracy"), +] + +PHASE_SHADING = { + "tower_warmup": "#d0e8ff", + "fused_warmup": "#d0ffe8", +} + + +def _load_folds(mode_dir: Path) -> list[pd.DataFrame]: + fold_dirs = sorted( + [p for p in mode_dir.glob("fold*") if p.is_dir()], + key=lambda p: int(p.name.replace("fold", "")), + ) + frames = [] + for fd in fold_dirs: + csv = fd / "epoch_log.csv" + if not csv.exists(): + print(f" [warn] {csv} not found, skipping") + continue + df = pd.read_csv(csv) + df["_fold"] = int(fd.name.replace("fold", "")) + frames.append(df) + return frames + + +def _shade_warmup(ax: plt.Axes, df: pd.DataFrame) -> None: + """Shade warmup phase regions based on first fold's phase column.""" + if "phase" not in df.columns: + return + prev_phase = None + start = None + for _, row in df.iterrows(): + phase = row["phase"] + ep = row["epoch"] + if phase != prev_phase: + if prev_phase in PHASE_SHADING and start is not None: + ax.axvspan(start - 0.5, ep - 0.5, color=PHASE_SHADING[prev_phase], + alpha=0.35, zorder=0, label=f"{prev_phase.replace('_', ' ')}") + start = ep + prev_phase = phase + # close last span + if prev_phase in PHASE_SHADING and start is not None: + ax.axvspan(start - 0.5, df["epoch"].max() + 0.5, + color=PHASE_SHADING[prev_phase], alpha=0.35, zorder=0) + + +def plot_curves( + run_dirs: list[Path], + eval_mode: str, + tower_mode: str, + out_dir: Path | None, +) -> None: + # Collect (label, frames) pairs + datasets: list[tuple[str, list[pd.DataFrame]]] = [] + for rd in run_dirs: + mode_dir = rd / eval_mode / tower_mode + if not mode_dir.exists(): + print(f" [skip] {mode_dir} not found") + continue + frames = _load_folds(mode_dir) + if not frames: + print(f" [skip] no epoch_log.csv found under {mode_dir}") + continue + datasets.append((rd.name, frames)) + + if not datasets: + print("No data found — nothing to plot.") + return + + # One figure per metric + for metric_key, metric_label in METRICS: + # Check any fold actually has this metric with non-nan values + has_data = any( + not frames[0][metric_key].isna().all() + for _, frames in datasets + if frames and metric_key in frames[0].columns + ) + if not has_data: + continue + + n_runs = len(datasets) + fig, axes = plt.subplots(1, n_runs, figsize=(5 * n_runs, 4.5), squeeze=False) + + for col_idx, (run_label, frames) in enumerate(datasets): + ax = axes[0][col_idx] + if frames and metric_key in frames[0].columns: + _shade_warmup(ax, frames[0]) + + colours = plt.cm.tab10(np.linspace(0, 0.9, len(frames))) + for frame, colour in zip(frames, colours): + if metric_key not in frame.columns: + continue + vals = frame[metric_key].values + epochs = frame["epoch"].values + mask = ~np.isnan(vals.astype(float)) + if mask.sum() == 0: + continue + ax.plot(epochs[mask], vals[mask], + linewidth=1.4, color=colour, + label=f"fold {frame['_fold'].iloc[0]}") + + ax.set_title(run_label, fontsize=10) + ax.set_xlabel("Epoch") + if col_idx == 0: + ax.set_ylabel(metric_label) + ax.legend(fontsize=7, loc="lower right") + ax.grid(True, linewidth=0.4, alpha=0.5) + + fig.suptitle(f"{metric_label} [{eval_mode} / {tower_mode}]", fontsize=12) + fig.tight_layout() + + dest = out_dir or (run_dirs[0].parent / "mdonly_plots") + dest.mkdir(parents=True, exist_ok=True) + fname = f"mdonly_{metric_key}_{eval_mode}_{tower_mode}.png" + fig.savefig(dest / fname, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"Saved: {dest / fname}") + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--run-dirs", nargs="+", required=True, + help="One or more run directories (e.g. analysis_data/pipeline_mdonly_50ep).") + ap.add_argument("--eval-mode", default="binary", choices=["binary", "multiclass"]) + ap.add_argument("--tower-mode", default="single", choices=["single", "ensemble"]) + ap.add_argument("--out", default=None, + help="Output directory for plots (default: {first_run_dir}/../mdonly_plots).") + args = ap.parse_args() + + run_dirs = [Path(d) for d in args.run_dirs] + out_dir = Path(args.out) if args.out else None + plot_curves(run_dirs, args.eval_mode, args.tower_mode, out_dir) + + +if __name__ == "__main__": + main() diff --git a/scripts/output_analysis/visualizations/plot_prob_strips.py b/scripts/output_analysis/visualizations/plot_prob_strips.py index a563909..1b231e2 100644 --- a/scripts/output_analysis/visualizations/plot_prob_strips.py +++ b/scripts/output_analysis/visualizations/plot_prob_strips.py @@ -80,22 +80,74 @@ def _fold_colours(n_folds: int) -> dict[int, dict[int, tuple]]: return out -def _load_folds(run_dir: Path, head: str) -> list[pd.DataFrame]: - fold_dirs = sorted( +def _fold_dirs(run_dir: Path) -> list[Path]: + return sorted( [p for p in run_dir.glob("fold*") if p.is_dir() and re.search(r"\d+", p.name)], key=lambda p: int(re.search(r"\d+", p.name).group()), ) - if not fold_dirs: + + +def _detect_heads(run_dir: Path) -> list[str]: + """Return all head names available in the first non-empty fold dir.""" + for fd in _fold_dirs(run_dir): + found: list[str] = [] + # Standard heads come from the predictions CSV — prefer per-eye + csv = fd / "predictions_pereye.csv" + if not csv.exists(): + csv = fd / "predictions_classic.csv" + if not csv.exists(): + csv = fd / "predictions.csv" + if csv.exists(): + cols = pd.read_csv(csv, nrows=0).columns.tolist() + for h in ["fused", "img", "md"]: + if f"prob_{h}_c0" in cols: + found.append(h) + # Fusion head lives in a separate npy + if (fd / "probs_fused_head.npy").exists(): + found.append("fused_head") + if found: + return found + return ["fused"] # safe fallback + + +def _load_folds(run_dir: Path, head: str) -> list[pd.DataFrame]: + dirs = _fold_dirs(run_dir) + if not dirs: raise FileNotFoundError(f"No fold* directories found under {run_dir}") frames = [] - for fd in fold_dirs: - csv = fd / "predictions_classic.csv" + for fd in dirs: + fold_num = int(re.search(r"\d+", fd.name).group()) + # fused_head is stored as npy, not in the predictions CSV + if head == "fused_head": + y_path = fd / "y_true.npy" + p_path = fd / "probs_fused_head.npy" + if not y_path.exists() or not p_path.exists(): + print(f" [warn] fused_head npy not found in {fd}, skipping") + continue + y = np.load(y_path) + p = np.load(p_path) + df = pd.DataFrame({"y_true": y}) + for c in range(p.shape[1]): + df[f"prob_fused_head_c{c}"] = p[:, c] + df["_fold"] = fold_num + frames.append(df) + continue + # Standard heads from predictions CSV — prefer per-eye (2× dots, no OD/OS averaging) + csv = fd / "predictions_pereye.csv" if not csv.exists(): - print(f" [warn] {csv} not found, skipping") + csv = fd / "predictions_classic.csv" + if not csv.exists(): + csv = fd / "predictions.csv" + if not csv.exists(): + print(f" [warn] no predictions CSV found in {fd}, skipping") continue df = pd.read_csv(csv) - df["_fold"] = int(re.search(r"\d+", fd.name).group()) + df["_fold"] = fold_num frames.append(df) + if not frames: + raise FileNotFoundError( + f"No data found for head='{head}' in any fold dir under {run_dir}" + ) return frames @@ -627,7 +679,9 @@ def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--run-dir", required=True) - ap.add_argument("--head", default="fused", choices=["fused", "img", "md"]) + ap.add_argument("--head", default=None, + choices=["fused", "fused_head", "img", "md"], + help="Head to plot. Omit to auto-detect and plot all available heads.") ap.add_argument("--style", default="strips", choices=["strips", "confidence", "triangle", "triangle3d"], help="strips: X=true class, Y=P(Glaucoma). " @@ -638,14 +692,22 @@ def main(): rd = Path(args.run_dir) od = Path(args.out) if args.out else None - if args.style == "confidence": - plot_confidence(rd, head=args.head, out_dir=od) - elif args.style == "triangle": - plot_triangle(rd, head=args.head, out_dir=od) - elif args.style == "triangle3d": - plot_triangle_3d(rd, head=args.head, out_dir=od) - else: - plot_strip(rd, head=args.head, out_dir=od) + + heads = [args.head] if args.head else _detect_heads(rd) + print(f"Heads to plot: {heads}") + + plot_fn = { + "confidence": plot_confidence, + "triangle": plot_triangle, + "triangle3d": plot_triangle_3d, + }.get(args.style, plot_strip) + + for head in heads: + print(f"\n--- {head} ---") + try: + plot_fn(rd, head=head, out_dir=od) + except Exception as exc: + print(f" [skip] {head}: {exc}") if __name__ == "__main__": diff --git a/scripts/output_analysis/visualizations/plot_run_roc_v2.py b/scripts/output_analysis/visualizations/plot_run_roc_v2.py index 03a3b8f..f1be4ba 100644 --- a/scripts/output_analysis/visualizations/plot_run_roc_v2.py +++ b/scripts/output_analysis/visualizations/plot_run_roc_v2.py @@ -57,12 +57,10 @@ _PROBS_PRIORITY: dict[str, list[str]] = { _ALL_PROBS = ["probs_fused_head", "probs_fused", "probs_bilat", "probs_classic"] -def detect_probs_stem(fold_dir: Path, tower_mode: str | None) -> str | None: +def detect_probs_stems(fold_dir: Path, tower_mode: str | None) -> list[str]: + """Return all present probs stems (in priority order) for this fold dir.""" 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 + return [stem for stem in priority if (fold_dir / f"{stem}.npy").exists()] # --------------------------------------------------------------------------- @@ -246,39 +244,44 @@ def main() -> None: 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: + # Determine probs stems to plot + if args.probs is not None: + stems_to_plot = [args.probs] + else: + # Collect all stems present across any fold dir + seen: list[str] = [] 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.") + for s in detect_probs_stems(fd, args.tower_mode): + if s not in seen: + seen.append(s) + stems_to_plot = seen + if not stems_to_plot: + raise SystemExit(f"Could not detect any probs file in {mode_dir}/fold*/") + print(f"Probs stems to plot: {stems_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) + for probs_stem in stems_to_plot: + print(f"\n--- {probs_stem} ---") + 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: + print(f" No usable folds for {probs_stem}, skipping.") + continue + + 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}")