diff --git a/v4/classes/accessory/transforms.py b/v4/classes/accessory/transforms.py index 70145af..bd3e5dc 100644 --- a/v4/classes/accessory/transforms.py +++ b/v4/classes/accessory/transforms.py @@ -76,39 +76,75 @@ class ImageTransformConfig: return transforms.Compose(ops) -def backbone_transform_config(backbone_name: str, augment: bool = True) -> ImageTransformConfig: - """Build an ImageTransformConfig using the backbone's default normalisation stats.""" +def backbone_transform_config( + backbone_name: str, + augment: bool = True, + crop_size: int | None = None, + resize_size: int | None = None, +) -> ImageTransformConfig: + """Build an ImageTransformConfig using the backbone's default normalisation stats. + + crop_size / resize_size override the backbone's default input resolution. When + crop_size is overridden but resize_size is not, resize_size is scaled + proportionally (8/7 ratio, matching the standard 224 → 256 pattern). + """ key = (backbone_name or "").lower() if _is_timm_backbone(key): # ConvNeXt-V2 and other timm models we currently expose are all # pretrained with standard ImageNet stats at 224×224. - return ImageTransformConfig(crop_size=224, mean=IMAGENET_MEAN, - std=IMAGENET_STD, augment=augment) - if key not in BACKBONES: - raise ValueError(f"Unknown backbone '{backbone_name}'.") - spec = BACKBONES[key] - mean = getattr(spec.weights_default, "meta", {}).get("mean", IMAGENET_MEAN) - std = getattr(spec.weights_default, "meta", {}).get("std", IMAGENET_STD) - crop = 299 if key == "inception_v3" else 224 - return ImageTransformConfig(crop_size=crop, mean=mean, std=std, augment=augment) + mean, std = IMAGENET_MEAN, IMAGENET_STD + default_crop = 224 + else: + if key not in BACKBONES: + raise ValueError(f"Unknown backbone '{backbone_name}'.") + spec = BACKBONES[key] + mean = getattr(spec.weights_default, "meta", {}).get("mean", IMAGENET_MEAN) + std = getattr(spec.weights_default, "meta", {}).get("std", IMAGENET_STD) + default_crop = 299 if key == "inception_v3" else 224 + + crop = crop_size if crop_size is not None else default_crop + resize = resize_size if resize_size is not None else round(crop * 8 / 7) + return ImageTransformConfig(crop_size=crop, resize_size=resize, + mean=mean, std=std, augment=augment) -def build_backbone_transform(backbone_name: str, augment: bool = True) -> transforms.Compose: - return backbone_transform_config(backbone_name, augment=augment).build() +def build_backbone_transform( + backbone_name: str, + augment: bool = True, + crop_size: int | None = None, + resize_size: int | None = None, +) -> transforms.Compose: + return backbone_transform_config( + backbone_name, augment=augment, + crop_size=crop_size, resize_size=resize_size, + ).build() -def build_eval_transform(backbone_name: str) -> transforms.Compose: +def build_eval_transform( + backbone_name: str, + crop_size: int | None = None, + resize_size: int | None = None, +) -> transforms.Compose: """Deterministic eval transform — no augmentation, backbone-matched normalisation.""" - return build_backbone_transform(backbone_name, augment=False) + return build_backbone_transform( + backbone_name, augment=False, + crop_size=crop_size, resize_size=resize_size, + ) def build_split_transforms( - backbone_name: str, augment: bool = True + backbone_name: str, + augment: bool = True, + crop_size: int | None = None, + resize_size: int | None = None, ) -> tuple[transforms.Compose, transforms.Compose]: """Return (precache, postcache) transform pair for tensor-cached image towers. precache : PIL → CHW float32 in [0, 1] (deterministic, run once at fill) postcache : tensor → augmented + normalized tensor (run per batch) """ - cfg = backbone_transform_config(backbone_name, augment=augment) + cfg = backbone_transform_config( + backbone_name, augment=augment, + crop_size=crop_size, resize_size=resize_size, + ) return cfg.build_precache(), cfg.build_postcache() diff --git a/v4/classes/logging/prediction_store.py b/v4/classes/logging/prediction_store.py index 9874145..0fb2a91 100644 --- a/v4/classes/logging/prediction_store.py +++ b/v4/classes/logging/prediction_store.py @@ -6,19 +6,19 @@ FeatureStore — records embeddings (opt-in); same structure but per-head HDF5 layout — PredictionStore ------------------------------ -/{phase}/logits float32 (n_folds, n_epochs, n_samples, n_heads, n_classes) -/{phase}/head_names str (n_heads,) -/{phase}/y_true int64 (n_samples,) -/{phase}/entity_id_{k} int64|str (n_samples,) — one dataset per id component -/{phase}/split str (n_folds, n_samples) -/{phase}/loss float32 (n_folds, n_epochs) +/{phase}/logits float32 (n_folds, n_epochs, n_samples, n_heads, n_classes) +/{phase}/head_names str (n_heads,) +/{phase}/y_true int64 | float64 (n_samples,) — float64 for regression targets, int64 otherwise +/{phase}/entity_id_{k} int64|str (n_samples,) — one dataset per id component +/{phase}/split str (n_folds, n_samples) +/{phase}/loss float32 (n_folds, n_epochs) HDF5 layout — FeatureStore --------------------------- -/{phase}/{head_name} float32 (n_folds, n_epochs, n_samples, embedding_dim) -/{phase}/y_true int64 (n_samples,) -/{phase}/entity_id_{k} int64|str (n_samples,) -/{phase}/split str (n_folds, n_samples) +/{phase}/{head_name} float32 (n_folds, n_epochs, n_samples, embedding_dim) +/{phase}/y_true int64 | float64 (n_samples,) — float64 for regression targets, int64 otherwise +/{phase}/entity_id_{k} int64|str (n_samples,) +/{phase}/split str (n_folds, n_samples) """ from __future__ import annotations @@ -35,6 +35,18 @@ except ImportError as e: _STR_DT = h5py.string_dtype() +def _coerce_y_true(y_true) -> np.ndarray: + """Coerce y_true to int64 for integer-typed input, float64 otherwise. + + Forcing int64 unconditionally would silently round regression targets + (e.g. VF_MD), so we honour float input by storing as float64. + """ + arr = np.asarray(y_true) + if np.issubdtype(arr.dtype, np.floating): + return arr.astype(np.float64) + return arr.astype(np.int64) + + # --------------------------------------------------------------------------- # Internal phase buffer # --------------------------------------------------------------------------- @@ -52,7 +64,7 @@ class _PhaseBuffer: n_s = len(entity_ids) n_h = len(head_names) self.entity_ids = list(entity_ids) - self.y_true = np.asarray(y_true, dtype=np.int64) + self.y_true = _coerce_y_true(y_true) self.head_names = list(head_names) self.n_epochs = n_epochs self.logits = np.full((n_folds, n_epochs, n_s, n_h, n_classes), np.nan, dtype=np.float32) @@ -74,7 +86,7 @@ class _FeaturePhaseBuffer: n_folds: int, ): self.entity_ids = list(entity_ids) - self.y_true = np.asarray(y_true, dtype=np.int64) + self.y_true = _coerce_y_true(y_true) self.split = np.full((n_folds, len(entity_ids)), "", dtype=object) self._sid = {str(eid): i for i, eid in enumerate(entity_ids)} # head_name → (buffer array, n_epochs) @@ -148,7 +160,7 @@ class PredictionStore: """Register a training phase before recording begins.""" self._phases[phase] = _PhaseBuffer( entity_ids=list(entity_ids), - y_true=np.asarray(y_true, dtype=np.int64), + y_true=_coerce_y_true(y_true), head_names=list(head_names), n_epochs=n_epochs, n_folds=self.n_folds, @@ -306,7 +318,7 @@ class FeatureStore: ) -> None: self._phases[phase] = _FeaturePhaseBuffer( entity_ids=list(entity_ids), - y_true=np.asarray(y_true, dtype=np.int64), + y_true=_coerce_y_true(y_true), n_folds=self.n_folds, ) diff --git a/v4/classes/profiles/fundus_images.py b/v4/classes/profiles/fundus_images.py index 468b559..b34bd94 100644 --- a/v4/classes/profiles/fundus_images.py +++ b/v4/classes/profiles/fundus_images.py @@ -650,6 +650,177 @@ def build_geometry_loader(source: str, **kwargs): raise NotImplementedError(f"build_geometry_loader: source={source!r} not implemented") +class _GTContourBboxLoader: + """Disc bounding-box loader from PAPILA expert disc contours. + + Computes a square bbox centred on the disc, expanded by ``margin`` × + max(disc_w, disc_h). Returned bboxes are in original-image pixel coords + and may extend past image bounds (clip at crop time). + """ + + def __init__( + self, + contour_dir: str | Path, + *, + margin: float = 2.5, + expert: int = 1, + ) -> None: + self._contour_dir = Path(contour_dir) + self._margin = float(margin) + self._expert = int(expert) + self._cache: dict[tuple, tuple[int, int, int, int] | None] = {} + + def reset_cache(self) -> None: + self._cache.clear() + + def precompute(self, samples: Iterable[Tuple[int, str, Path]]) -> None: + for pid, eye, _ in list(samples): + key = (int(pid), str(eye)) + if key in self._cache: + continue + self._cache[key] = self._compute_bbox(*key) + + def bbox_for(self, pid, eye) -> tuple[int, int, int, int] | None: + key = (int(pid), str(eye)) + if key not in self._cache: + self._cache[key] = self._compute_bbox(*key) + return self._cache[key] + + def _compute_bbox(self, pid: int, eye: str) -> tuple[int, int, int, int] | None: + path = self._contour_dir / f"RET{pid:03d}{eye}_disc_exp{self._expert}.txt" + if not path.exists(): + return None + try: + arr = np.loadtxt(str(path), dtype=np.float32) + except Exception: + return None + if arr.ndim == 1: + arr = arr.reshape(-1, 2) + if arr.shape[0] < 3: + return None + return _expand_bbox_from_points(arr[:, 0], arr[:, 1], self._margin) + + +class _UNetBboxLoader: + """Disc bounding-box loader from U-Net predicted masks. + + Reuses _PapilaUNetMaskPipeline for per-fold fine-tune + inference. + """ + + def __init__( + self, + weights_path: str | Path, + *, + contour_dir: str | Path, + margin: float = 2.5, + unet_size: int = 512, + normalize: str = "per_image", + threshold: float = 0.5, + finetune_epochs: int = 0, + finetune_lr: float = 1e-5, + finetune_batch_size: int = 4, + device: str | None = None, + ) -> None: + self._margin = float(margin) + self._unet_size = int(unet_size) + self._pipeline = _PapilaUNetMaskPipeline( + weights_path, + contour_dir=contour_dir, + unet_size=unet_size, + normalize=normalize, + threshold=threshold, + finetune_epochs=finetune_epochs, + finetune_lr=finetune_lr, + finetune_batch_size=finetune_batch_size, + device=device, + ) + self._cache: dict[tuple, tuple[int, int, int, int] | None] = {} + + def reset_cache(self) -> None: + self._cache.clear() + + def reset_weights(self) -> None: + self._pipeline.reset_weights() + + def finetune(self, train_samples: list) -> None: + self._pipeline.finetune(train_samples) + + def precompute(self, samples: Iterable[Tuple[int, str, Path]]) -> None: + samples = list(samples) + if not samples: + return + # Need original-image dims to rescale mask coords back; capture per sample. + orig_sizes: dict[tuple, tuple[int, int]] = {} + for pid, eye, image_path in samples: + try: + with Image.open(image_path) as im: + orig_sizes[(int(pid), str(eye))] = im.size # (w, h) + except Exception: + continue + masks = self._pipeline.predict(samples) + for key, (disc, _cup) in masks.items(): + ow, oh = orig_sizes.get(key, (self._unet_size, self._unet_size)) + self._cache[key] = _bbox_from_mask(disc, ow, oh, self._margin) + + def bbox_for(self, pid, eye) -> tuple[int, int, int, int] | None: + return self._cache.get((int(pid), str(eye))) + + +def _expand_bbox_from_points( + xs: np.ndarray, ys: np.ndarray, margin: float +) -> tuple[int, int, int, int]: + """Square bbox centred on disc centroid, half-side = margin × max(w,h) / 2.""" + x0, x1 = float(xs.min()), float(xs.max()) + y0, y1 = float(ys.min()), float(ys.max()) + cx, cy = (x0 + x1) / 2.0, (y0 + y1) / 2.0 + half = max(x1 - x0, y1 - y0) * float(margin) / 2.0 + return (int(round(cx - half)), int(round(cy - half)), + int(round(cx + half)), int(round(cy + half))) + + +def _bbox_from_mask( + mask: np.ndarray, orig_w: int, orig_h: int, margin: float, +) -> tuple[int, int, int, int] | None: + """Compute original-image bbox from a binary mask at mask resolution.""" + ys, xs = np.where(mask > 0) + if len(xs) == 0: + return None + sx = float(orig_w) / float(mask.shape[1]) + sy = float(orig_h) / float(mask.shape[0]) + return _expand_bbox_from_points(xs * sx, ys * sy, margin) + + +def build_disc_bbox_loader(source: str, **kwargs): + """Return a disc-bbox loader for the given source. + + Parameters + ---------- + source : "gt" | "unet" + + GT kwargs: + contour_dir (required), margin=2.5, expert=1 + UNet kwargs: + weights_path (required), contour_dir (required), + margin=2.5, unet_size=512, normalize="per_image", threshold=0.5, + finetune_epochs=0, finetune_lr=1e-5, finetune_batch_size=4, device=None + """ + if source == "gt": + if "contour_dir" not in kwargs: + raise ValueError("build_disc_bbox_loader source='gt' requires contour_dir") + gt_keys = {"contour_dir", "margin", "expert"} + return _GTContourBboxLoader(**{k: v for k, v in kwargs.items() if k in gt_keys}) + if source == "unet": + if "weights_path" not in kwargs: + raise ValueError("build_disc_bbox_loader source='unet' requires weights_path") + if "contour_dir" not in kwargs: + raise ValueError( + "build_disc_bbox_loader source='unet' requires contour_dir " + "(needed for per-fold fine-tuning, even if finetune_epochs=0)" + ) + return _UNetBboxLoader(**kwargs) + raise NotImplementedError(f"build_disc_bbox_loader: source={source!r} not implemented") + + def build_seg_map_loader(source: str, **kwargs): """Return the appropriate seg-map loader for the given source string. diff --git a/v4/classes/profiles/v4papila.py b/v4/classes/profiles/v4papila.py index d2cdf91..985f114 100644 --- a/v4/classes/profiles/v4papila.py +++ b/v4/classes/profiles/v4papila.py @@ -480,6 +480,12 @@ class ImageDataView: from v4.classes.profiles.fundus_images import build_seg_map_loader as _build return _build(source, **self._resolve_paths(kwargs)) + def build_disc_bbox_loader(self, source: str, **kwargs): + """Return a disc bounding-box loader for crop-to-disc preprocessing.""" + from v4.classes.profiles.fundus_images import build_disc_bbox_loader as _build + kwargs.setdefault("contour_dir", self._DEFAULT_CONTOUR_DIR) + return _build(source, **self._resolve_paths(kwargs)) + # ── Optional explainability hooks ──────────────────────────────────────── # # These methods are consumed by v4.classes.accessory.explainability via diff --git a/v4/classes/stages/fusion.py b/v4/classes/stages/fusion.py index 5ebf03b..e28c6b7 100644 --- a/v4/classes/stages/fusion.py +++ b/v4/classes/stages/fusion.py @@ -1,6 +1,7 @@ """stages/fusion — fusion stage runner: trains a bridge + associated head stages.""" from __future__ import annotations +import contextlib import importlib from random import choice, random as _random @@ -8,6 +9,20 @@ import numpy as np import torch import torch.nn.functional as F + +def _amp_ctx(cfg: dict, device): + """Autocast context for forward+loss when training.amp is enabled. + + bf16 is the default dtype because its dynamic range matches fp32 and no + GradScaler is required. Falls back to a no-op context when amp is disabled + or the device is not CUDA/ROCm. + """ + train_cfg = cfg.get("training", {}) + if not train_cfg.get("amp", False) or getattr(device, "type", None) != "cuda": + return contextlib.nullcontext() + dtype = getattr(torch, train_cfg.get("amp_dtype", "bfloat16")) + return torch.autocast(device_type="cuda", dtype=dtype) + from v4.classes.dataset import LoaderShell, to_label_tensor from v4.classes.metrics import score_arrays, compute_extended_metrics, tune_binary_threshold from v4.classes.stages.helpers import ( @@ -223,82 +238,78 @@ def run( if y_t.numel() == 0: continue - if is_bilateral: - side_embs = { - side: encode_embedding(src, batch, side, towers, stage_models, cfg_stages, device) - for side, src in inputs.items() - } - local_embs = {name: bridge(side_embs)} - else: - local_embs = {n: towers[n](batch[n].to(device)) for n in inputs - if n in batch and torch.is_tensor(batch[n])} - if len(local_embs) != len(inputs): - continue - local_embs[name] = bridge(list(local_embs[n] for n in inputs)) + loss = None + logits = None - head_logits = { - hs["name"]: head_models[hs["name"]](local_embs[hs["input"]]) - for hs in head_stage_cfgs - if hs["input"] in local_embs - } - - if is_bilateral or phase == "fused_warmup": - logits = head_logits.get(primary_hs_cfg["name"]) - chosen_head = head_models.get(primary_hs_cfg["name"]) - elif phase == "tower_warmup" and bcd_head_cfgs: - losses = [ - head_compute_loss(head_models[hs["name"]], - head_logits[hs["name"]], batch, y_t, - class_weights=cw) - for hs in bcd_head_cfgs if hs["name"] in head_logits - ] - if not losses: - continue - loss = sum(losses) / len(losses) - if hasattr(bridge, "modify_loss"): - loss = bridge.modify_loss(loss) - opt.zero_grad(); loss.backward(); opt.step() - total_loss += loss.item() * len(y_t) - total_n += len(y_t) - continue - elif tower_loss_mode == "all_losses" and bcd_head_cfgs: - # All-losses (v3 phase 3 control): sum primary + every aux head - # loss every step. Effective LR is implicitly N× single-head BCD - # — matches v3 semantics so the comparison is apples-to-apples. - all_head_names = ([primary_hs_cfg["name"]] - + [hs["name"] for hs in bcd_head_cfgs]) - losses = [ - head_compute_loss(head_models[n], head_logits[n], batch, y_t, - class_weights=cw) - for n in all_head_names if n in head_logits - ] - if not losses: - continue - loss = sum(losses) - if hasattr(bridge, "modify_loss"): - loss = bridge.modify_loss(loss) - opt.zero_grad(); loss.backward(); opt.step() - total_loss += loss.item() * len(y_t) - total_n += len(y_t) - continue - else: - if bcd_head_cfgs and _random() < bcd_prob: - chosen_hs = choice(bcd_head_cfgs) + with _amp_ctx(cfg, device): + if is_bilateral: + side_embs = { + side: encode_embedding(src, batch, side, towers, stage_models, cfg_stages, device) + for side, src in inputs.items() + } + local_embs = {name: bridge(side_embs)} else: - chosen_hs = primary_hs_cfg - logits = head_logits.get(chosen_hs["name"]) - chosen_head = head_models.get(chosen_hs["name"]) + local_embs = {n: towers[n](batch[n].to(device)) for n in inputs + if n in batch and torch.is_tensor(batch[n])} + if len(local_embs) != len(inputs): + continue + local_embs[name] = bridge(list(local_embs[n] for n in inputs)) - if logits is None: + head_logits = { + hs["name"]: head_models[hs["name"]](local_embs[hs["input"]]) + for hs in head_stage_cfgs + if hs["input"] in local_embs + } + + if is_bilateral or phase == "fused_warmup": + logits = head_logits.get(primary_hs_cfg["name"]) + chosen_head = head_models.get(primary_hs_cfg["name"]) + if logits is not None: + loss = head_compute_loss(chosen_head, logits, batch, y_t, class_weights=cw) + elif phase == "tower_warmup" and bcd_head_cfgs: + losses = [ + head_compute_loss(head_models[hs["name"]], + head_logits[hs["name"]], batch, y_t, + class_weights=cw) + for hs in bcd_head_cfgs if hs["name"] in head_logits + ] + if losses: + loss = sum(losses) / len(losses) + elif tower_loss_mode == "all_losses" and bcd_head_cfgs: + # All-losses (v3 phase 3 control): sum primary + every aux head + # loss every step. Effective LR is implicitly N× single-head BCD + # — matches v3 semantics so the comparison is apples-to-apples. + all_head_names = ([primary_hs_cfg["name"]] + + [hs["name"] for hs in bcd_head_cfgs]) + losses = [ + head_compute_loss(head_models[n], head_logits[n], batch, y_t, + class_weights=cw) + for n in all_head_names if n in head_logits + ] + if losses: + loss = sum(losses) + else: + if bcd_head_cfgs and _random() < bcd_prob: + chosen_hs = choice(bcd_head_cfgs) + else: + chosen_hs = primary_hs_cfg + logits = head_logits.get(chosen_hs["name"]) + chosen_head = head_models.get(chosen_hs["name"]) + if logits is not None: + loss = head_compute_loss(chosen_head, logits, batch, y_t, class_weights=cw) + + if loss is not None and hasattr(bridge, "modify_loss"): + loss = bridge.modify_loss(loss) + + if loss is None: continue - loss = head_compute_loss(chosen_head, logits, batch, y_t, class_weights=cw) - if hasattr(bridge, "modify_loss"): - loss = bridge.modify_loss(loss) + opt.zero_grad(); loss.backward(); opt.step() - if logits.dim() >= 2: + + if logits is not None and logits.dim() >= 2: total_correct += int((logits.argmax(1) == y_t).sum()) - total_loss += loss.item() * len(y_t) - total_n += len(y_t) + total_loss += loss.item() * len(y_t) + total_n += len(y_t) tr_loss = total_loss / total_n if total_n else nan tr_acc = total_correct / total_n if total_n else nan diff --git a/v4/classes/towers/image_tower.py b/v4/classes/towers/image_tower.py index cadf79a..ac33c9d 100644 --- a/v4/classes/towers/image_tower.py +++ b/v4/classes/towers/image_tower.py @@ -71,6 +71,13 @@ class ImageEncoder(TowerBase): Normalize on tensors only (no PIL, no Resize, no decode). Memory: ~3 × crop_size² × 4B per cached image. Cache is rebuilt at the start of every fold via early_pass. + crop_source : if set, crop each input image to a square disc-region + bbox before the standard transform pipeline. Values: + "gt" (use GT contour file) | "unet" (use U-Net mask) + | None (disabled, full-image pipeline). + crop_kwargs : dict forwarded to image_data.build_disc_bbox_loader(). + Common keys: margin (default 2.5), expert (GT only), + weights_path / finetune_epochs (U-Net only). geometry_source : source key passed to image_data.build_geometry_loader() (e.g. "gt", "unet"). None = geometry disabled. **geom_kwargs : forwarded verbatim to build_geometry_loader() — e.g. @@ -89,6 +96,10 @@ class ImageEncoder(TowerBase): se_pre_norm: bool = True, augment: bool = True, cache_transformed: bool = False, + crop_size: int | None = None, + resize_size: int | None = None, + crop_source: str | None = None, + crop_kwargs: dict | None = None, geometry_source: str | None = None, **geom_kwargs: Any, ): @@ -98,17 +109,37 @@ class ImageEncoder(TowerBase): self.backbone, self._base_dim, self._blocks = build_backbone(backbone, freeze_ratio) self._cache_transformed = cache_transformed + tf_kw = dict(crop_size=crop_size, resize_size=resize_size) if cache_transformed: - self._precache_tf, self._post_train_tf = build_split_transforms(backbone, augment=augment) - _, self._post_eval_tf = build_split_transforms(backbone, augment=False) + self._precache_tf, self._post_train_tf = build_split_transforms( + backbone, augment=augment, **tf_kw) + _, self._post_eval_tf = build_split_transforms( + backbone, augment=False, **tf_kw) self._tensor_cache: dict[tuple, torch.Tensor] = {} else: - self.transform = build_backbone_transform(backbone, augment=augment) - self.eval_transform = build_eval_transform(backbone) + self.transform = build_backbone_transform(backbone, augment=augment, **tf_kw) + self.eval_transform = build_eval_transform(backbone, **tf_kw) self.tower_ln = nn.LayerNorm(self._base_dim) if se_pre_norm else nn.Identity() self.tower_se = SEBlock(self._base_dim, reduction=se_reduction, residual=True) if use_se else None + self._bbox_loader = None + if crop_source is not None: + if not hasattr(image_data, "build_disc_bbox_loader"): + raise TypeError( + f"ImageEncoder crop_source={crop_source!r} requires " + f"image_data to implement build_disc_bbox_loader(), " + f"but {type(image_data).__name__} does not." + ) + self._bbox_loader = image_data.build_disc_bbox_loader( + crop_source, **(crop_kwargs or {}), + ) + print( + f"[ImageEncoder] crop_source={crop_source!r} " + f"kwargs={crop_kwargs or {}}", + flush=True, + ) + self._geom_loader = None if geometry_source is not None: if not hasattr(image_data, "build_geometry_loader"): @@ -134,16 +165,32 @@ class ImageEncoder(TowerBase): def _side_map(self) -> dict[str, str]: return self.image_data.side_map + def _load_image(self, *ids): + """Load image, optionally cropped to the disc-region bbox.""" + pil = self.image_data.load_image(*ids) + if self._bbox_loader is None: + return pil + bbox = self._bbox_loader.bbox_for(*ids[:2]) + if bbox is None: + return pil + w, h = pil.size + x0, y0, x1, y1 = bbox + x0, y0 = max(0, x0), max(0, y0) + x1, y1 = min(w, x1), min(h, y1) + if x1 <= x0 or y1 <= y0: + return pil + return pil.crop((x0, y0, x1, y1)) + def _get(self, *ids) -> torch.Tensor: if self._cache_transformed: key = tuple(ids) cached = self._tensor_cache.get(key) if cached is None: - cached = self._precache_tf(self.image_data.load_image(*ids)) + cached = self._precache_tf(self._load_image(*ids)) self._tensor_cache[key] = cached tail = self._post_train_tf if self.training else self._post_eval_tf return tail(cached) - img = self.image_data.load_image(*ids) + img = self._load_image(*ids) t = self.transform if self.training else self.eval_transform return t(img) @@ -154,6 +201,24 @@ class ImageEncoder(TowerBase): data = context.require("data") split = context.require("split") + # Disc-region bbox precomputation (must run before any image load/cache). + if self._bbox_loader is not None: + train_samples = data.collect_samples(split.train) + all_samples = train_samples + data.collect_samples(split.val) + if split.test is not None: + all_samples += data.collect_samples(split.test) + if hasattr(self._bbox_loader, "reset_cache"): + self._bbox_loader.reset_cache() + if hasattr(self._bbox_loader, "reset_weights"): + self._bbox_loader.reset_weights() + if hasattr(self._bbox_loader, "finetune"): + self._bbox_loader.finetune(train_samples) + self._bbox_loader.precompute(all_samples) + print( + f"[ImageEncoder] precomputed disc bboxes for {len(all_samples)} samples", + flush=True, + ) + if self._cache_transformed: self._tensor_cache.clear() n = self._warm_tensor_cache(data, split) @@ -200,7 +265,7 @@ class ImageEncoder(TowerBase): key = (pid, eye) if key in self._tensor_cache or key in seen: continue - self._tensor_cache[key] = self._precache_tf(self.image_data.load_image(pid, eye)) + self._tensor_cache[key] = self._precache_tf(self._load_image(pid, eye)) seen.add(key) return len(self._tensor_cache) diff --git a/v4/classes/v4_hypertower.py b/v4/classes/v4_hypertower.py index f1a6313..b7b473a 100644 --- a/v4/classes/v4_hypertower.py +++ b/v4/classes/v4_hypertower.py @@ -359,21 +359,28 @@ def main(): if save_predictions and eval_stage_preds: # Collect all unique entity_ids across val+test sets of all folds. + # Preserve the natural dtype of y so regression targets keep their + # fractional values (casting to int silently rounds VF_MD). seen, all_ids, id_to_y = set(), [], {} + y_is_float = False for fp in eval_stage_preds: for eid, y in zip(fp["val_ids"], fp["val_y"]): k = str(eid) if k not in seen: seen.add(k); all_ids.append(eid) - id_to_y[k] = int(y) + y_is_float = y_is_float or np.issubdtype(np.asarray(y).dtype, np.floating) + id_to_y[k] = float(y) if y_is_float else int(y) if fp.get("test_ids"): for eid, y in zip(fp["test_ids"], fp["test_y"]): k = str(eid) if k not in seen: seen.add(k); all_ids.append(eid) - id_to_y[k] = int(y) + y_is_float = y_is_float or np.issubdtype(np.asarray(y).dtype, np.floating) + id_to_y[k] = float(y) if y_is_float else int(y) - y_true = np.array([id_to_y.get(str(e), -1) for e in all_ids], dtype=np.int64) + sentinel = float("nan") if y_is_float else -1 + dtype = np.float64 if y_is_float else np.int64 + y_true = np.array([id_to_y.get(str(e), sentinel) for e in all_ids], dtype=dtype) store = PredictionStore(n_folds=len(eval_stage_preds), n_classes=num_classes) store.register_phase( phase=eval_stage, @@ -402,19 +409,24 @@ def main(): for phase, phase_preds in all_phase_preds.items(): emb_dim = phase_preds[0]["val_z"].shape[-1] seen, all_ids, id_to_y = set(), [], {} + y_is_float = False for fp in phase_preds: for eid, y in zip(fp["val_ids"], fp["val_y"]): k = str(eid) if k not in seen: seen.add(k); all_ids.append(eid) - id_to_y[k] = int(y) + y_is_float = y_is_float or np.issubdtype(np.asarray(y).dtype, np.floating) + id_to_y[k] = float(y) if y_is_float else int(y) if fp.get("test_ids"): for eid, y in zip(fp["test_ids"], fp["test_y"]): k = str(eid) if k not in seen: seen.add(k); all_ids.append(eid) - id_to_y[k] = int(y) - y_true = np.array([id_to_y.get(str(e), -1) for e in all_ids], dtype=np.int64) + y_is_float = y_is_float or np.issubdtype(np.asarray(y).dtype, np.floating) + id_to_y[k] = float(y) if y_is_float else int(y) + sentinel = float("nan") if y_is_float else -1 + dtype = np.float64 if y_is_float else np.int64 + y_true = np.array([id_to_y.get(str(e), sentinel) for e in all_ids], dtype=dtype) fstore.register_phase(phase=phase, entity_ids=all_ids, y_true=y_true) fstore.register_head(phase=phase, head=f"{phase}_embedding", n_epochs=1, embedding_dim=emb_dim) diff --git a/v4/distributed/client.py b/v4/distributed/client.py index ff63e4a..7fd87f3 100644 --- a/v4/distributed/client.py +++ b/v4/distributed/client.py @@ -168,20 +168,28 @@ def _gpu_info() -> str: # rsync helpers # ────────────────────────────────────────────────────────────── -def _rsync(src: str, dst: str, delete: bool = False): +def _rsync(src: str, dst: str, delete: bool = False, + excludes: list[str] | None = None): cmd = ["rsync", "-az", "--info=progress2"] if delete: cmd.append("--delete") + for pat in (excludes or []): + cmd.append(f"--exclude={pat}") cmd += [src, dst] subprocess.run(cmd, check=True) def _sync_code(server_ssh: str, server_path: str, local_path: str): - """Pull v4/ source from server → local (overwrites local changes).""" + """Pull v4/ source from server → local (overwrites local changes). + + results/ is excluded so the client never overwrites or deletes its own + per-job output directory, and so it never pulls down the full corpus of + historical results from the server. + """ src = f"{server_ssh}:{server_path}/v4/" dst = f"{local_path}/v4/" Path(dst).mkdir(parents=True, exist_ok=True) - _rsync(src, dst, delete=True) + _rsync(src, dst, delete=True, excludes=["results/"]) def _upload_results(server_ssh: str, server_path: str, local_path: str, diff --git a/v4/distributed/jobs.db b/v4/distributed/jobs.db index 6eed7aa..9a44f29 100644 Binary files a/v4/distributed/jobs.db and b/v4/distributed/jobs.db differ diff --git a/v4/distributed/server.py b/v4/distributed/server.py index 955e316..0d0f9c6 100644 --- a/v4/distributed/server.py +++ b/v4/distributed/server.py @@ -413,8 +413,33 @@ def submit_job(job: JobSubmit): ) if cur.rowcount == 0: existing = conn.execute( - "SELECT job_id FROM jobs WHERE args=?", (args_json,) + "SELECT job_id, state FROM jobs WHERE args=?", (args_json,) ).fetchone() + # If the existing duplicate is a terminal failure, drop it and + # take the new submission — saves an explicit /jobs/clear round- + # trip when re-deploying after a fix. + if existing["state"] == "failed": + conn.execute( + "DELETE FROM jobs WHERE job_id=?", (existing["job_id"],) + ) + conn.execute( + "INSERT INTO jobs " + "(job_id, run_name, module, args, output_dir, priority, created_at) " + "VALUES (?,?,?,?,?,?,?)", + (job_id, job.run_name, job.module, args_json, + job.output_dir, job.priority, _now()), + ) + print( + f"[server] requeued failed {existing['job_id']} → {job_id} " + f"({job.run_name})", + flush=True, + ) + return { + "job_id": job_id, + "duplicate": False, + "requeued": True, + "previous_job_id": existing["job_id"], + } job_id = existing["job_id"] print(f"[server] duplicate ignored ({job.run_name}) → {job_id}", flush=True) return {"job_id": job_id, "duplicate": True} diff --git a/v4/figures/F2_papila_replication_and_single_mode.py b/v4/figures/F2_papila_replication_and_single_mode.py index 51c15dc..02891ad 100644 --- a/v4/figures/F2_papila_replication_and_single_mode.py +++ b/v4/figures/F2_papila_replication_and_single_mode.py @@ -1,50 +1,56 @@ -"""F2 — Backbone selection panel. +"""F2 - Backbone selection panel. Box plot in the style of v3/figures/phase2_analysis.png (black-bordered boxes, red median lines, baseline median reference). Three left-to-right sections: - Block 1 (blue) — Basic backbones (img-only, single-eye, ImageNet pretraining): + Block 1 (blue) -- Basic backbones (img-only, single-eye, ImageNet pretraining): VGG16, MobileNetV2, DenseNet121, InceptionV3, ResNet50 - Sourced from v3 phase 1 / phase 2 fold AUCs. Will be refined with v4 - 10x5 runs later; means should not move much. + Sourced from v4 experiments/backbone_replication/basic_* (10x5 = 50 fold-rep + AUCs each, img-only single-eye, patient-grouped CV). - Block 2 (blue) — ResNet50 preprocessing/CV variations: - leaky CV, GT crop, U-Net crop (all 2.5x scale; 1.1x dropped from labels) - Sourced from v3 phase 2 'classic_test_auc' (single-mode image-only). + Block 2 (blue) -- ResNet50 preprocessing/CV variations: + Anonymous CV, GT crop, U-Net crop (disc crops use margin 2.5x) + "Anonymous CV" = patient-identity-agnostic cross-validation: fold + assignment ignores PAPILA's patient IDs, allowing the same patient's + OD/OS pair to be split across train and test. Reflects the standard + protocol in benchmark reports that do not have patient-level labels + (or do not respect them). + Sourced from v4 experiments/backbone_replication/{anonymous_cv,gtcrop, + unetcrop}_refugelike. - Block 3 (orange) — Baseline reference: - "Baseline (fine-tuned ResNet50)" — what we previously called refugelike. - Sourced from v3 phase 2 imageonly_refugelike_proper. + Block 3 (orange) -- Baseline reference: + "Baseline (fine-tuned ResNet50)" -- REFUGE-pretrained R50 image-only, + sourced from v4 experiments/refuge_v2m_baseline/img_solo_single_refugelike. Each non-baseline box is labelled with a Wilcoxon two-sided p-value comparing -its fold AUCs to the baseline. +its fold-rep AUCs to the baseline fold-rep AUCs. Re-run anytime: python -m v4.figures.F2_papila_replication_and_single_mode """ from __future__ import annotations +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 scipy.stats import wilcoxon -from v4.figures.util.loaders import REPO_ROOT +from v4.figures.util.loaders import RESULTS_ROOT OUT = Path(__file__).parent / "output" / "F2_backbones.png" -# ── Colors / styling (mirrors v3 phase2_analysis) ──────────────────────────── -C_VAR = "#4c72b0" # blue — non-baseline boxes (basic backbones + variants) -C_BASE = "#dd8452" # orange — baseline reference box -C_MEDIAN = "#c44e52" # red — median line inside boxes +# Colors / styling (mirrors v3 phase2_analysis) +C_VAR = "#4c72b0" # blue - non-baseline boxes +C_BASE = "#dd8452" # orange - baseline reference box +C_MEDIAN = "#c44e52" # red - median line inside boxes ALPHA = 0.82 -V3_PHASE1_DIR = REPO_ROOT / "v3" / "results" / "phase1" -V3_PHASE2_DIR = REPO_ROOT / "v3" / "results" / "phase2" +# Stage key for img-only single-eye fusion +STAGE_KEY = "img_fuse_test_auc" def _wilcoxon_p(a: np.ndarray, b: np.ndarray) -> float: @@ -57,67 +63,64 @@ def _wilcoxon_p(a: np.ndarray, b: np.ndarray) -> float: return float("nan") -def _load_phase1_fold_aucs(subdir: str) -> np.ndarray: - fp = V3_PHASE1_DIR / subdir / "fold_metrics.csv" - if not fp.exists(): - return np.array([]) - df = pd.read_csv(fp) - return df["auc"].dropna().astype(float).values - - -def _load_phase2_classic_aucs(run_name: str) -> np.ndarray: - """Collect classic_test_auc across all rep×fold for a phase 2 run folder.""" - root = V3_PHASE2_DIR / run_name +def _load_fold_aucs(rel: str) -> np.ndarray: + """Collect STAGE_KEY across all rep x fold for a v4 results subdirectory.""" + root = RESULTS_ROOT / rel if not root.exists(): return np.array([]) out: list[float] = [] - for rep in sorted(root.glob("rep*")): - fp = rep / "binary" / "single" / "fold_results.csv" - if not fp.exists(): continue - df = pd.read_csv(fp) - if "classic_test_auc" not in df.columns: continue - out.extend(df["classic_test_auc"].dropna().astype(float).tolist()) + for s in sorted(root.glob("rep*/binary/summary.json")): + d = json.loads(s.read_text()) + for fr in d.get("fold_results", []): + v = fr.get(STAGE_KEY) + if v is None or not np.isfinite(v): + continue + out.append(float(v)) return np.array(out) -# ── Per-section data definitions ───────────────────────────────────────────── -# Each entry: (label, loader_fn, *args) +# Per-section data definitions (label, results-subdir under RESULTS_ROOT) BASIC_BACKBONES = [ - ("VGG16", _load_phase1_fold_aucs, "cnn_vgg16"), - ("MobileNetV2", _load_phase1_fold_aucs, "cnn_mobilenet_v2"), - ("DenseNet121", _load_phase1_fold_aucs, "cnn_densenet121"), - ("InceptionV3", _load_phase1_fold_aucs, "cnn_inception_v3"), - # Use phase 2 ResNet50 (50 fold AUCs) for tighter statistics on the - # backbone that we sweep variations of in block 2. - ("ResNet50", _load_phase2_classic_aucs, "imageonly_resnet50_proper"), + ("VGG16", "backbone_replication/basic_vgg16"), + ("MobileNetV2", "backbone_replication/basic_mobilenet_v2"), + ("DenseNet121", "backbone_replication/basic_densenet121"), + ("InceptionV3", "backbone_replication/basic_inception_v3"), + ("ResNet50", "backbone_replication/basic_resnet50"), ] RESNET_VARIATIONS = [ - ("leaky CV", _load_phase2_classic_aucs, "imageonly_resnet50_leaky"), - ("GT crop", _load_phase2_classic_aucs, "imageonly_resnet50_gtcrop_2.5"), - ("U-Net crop", _load_phase2_classic_aucs, "imageonly_resnet50_unetcrop_2.5"), + ("Anonymous CV", "backbone_replication/anonymous_cv_refugelike"), + ("GT crop", "backbone_replication/gtcrop_refugelike"), + ("U-Net crop", "backbone_replication/unetcrop_refugelike"), ] BASELINE_LABEL = "baseline\n(fine-tuned ResNet50)" -BASELINE_DATA = (_load_phase2_classic_aucs, "imageonly_refugelike_proper") +BASELINE_REL = "refuge_v2m_baseline/img_solo_single_refugelike" def render() -> None: - # Load everything - block1 = [(lbl, fn(arg)) for lbl, fn, arg in BASIC_BACKBONES] - block2 = [(lbl, fn(arg)) for lbl, fn, arg in RESNET_VARIATIONS] - base_fn, base_arg = BASELINE_DATA - base_aucs = base_fn(base_arg) + block1 = [(lbl, _load_fold_aucs(rel)) for lbl, rel in BASIC_BACKBONES] + block2 = [(lbl, _load_fold_aucs(rel)) for lbl, rel in RESNET_VARIATIONS] + base_aucs = _load_fold_aucs(BASELINE_REL) - print("Block 1 — Basic backbones:") + print("Block 1 - Basic backbones (ImageNet pretraining):") for lbl, a in block1: - print(f" {lbl:<14s} n={len(a):>3d} mean={a.mean():.3f}±{a.std():.3f}" if len(a) else f" {lbl:<14s} no data") - print("Block 2 — ResNet50 variations:") + if len(a): + print(f" {lbl:<14s} n={len(a):>3d} mean={a.mean():.4f} +/- {a.std():.4f}") + else: + print(f" {lbl:<14s} no data") + print("Block 2 - ResNet50 (REFUGE) variations:") for lbl, a in block2: - print(f" {lbl:<14s} n={len(a):>3d} mean={a.mean():.3f}±{a.std():.3f}" if len(a) else f" {lbl:<14s} no data") - print(f"Block 3 — Baseline: n={len(base_aucs)} " - f"mean={base_aucs.mean():.3f}±{base_aucs.std():.3f}" if len(base_aucs) else "Block 3 — no baseline data") + if len(a): + print(f" {lbl:<14s} n={len(a):>3d} mean={a.mean():.4f} +/- {a.std():.4f}") + else: + print(f" {lbl:<14s} no data") + if len(base_aucs): + print(f"Block 3 - Baseline: n={len(base_aucs)} " + f"mean={base_aucs.mean():.4f} +/- {base_aucs.std():.4f}") + else: + print("Block 3 - no baseline data") # Lay out positions gap = 0.7 @@ -135,7 +138,6 @@ def render() -> None: section3_left = p pos.append(p) section3_right = p - total_w = p + 0.6 fig, ax = plt.subplots(figsize=(13, 5.8)) fig.suptitle("Backbone Selection", fontsize=13, fontweight="bold") @@ -156,10 +158,10 @@ def render() -> None: all_labels.append(lbl); all_aucs.append(a); all_colors.append(C_VAR) all_labels.append(BASELINE_LABEL); all_aucs.append(base_aucs); all_colors.append(C_BASE) - # Draw boxes for x, aucs, color in zip(pos, all_aucs, all_colors): - if not len(aucs): continue - bp = ax.boxplot( + if not len(aucs): + continue + ax.boxplot( aucs, positions=[x], widths=box_w, patch_artist=True, manage_ticks=False, boxprops=dict(facecolor=color, alpha=ALPHA, **boxprops_kw), medianprops=medianprops, @@ -168,31 +170,27 @@ def render() -> None: flierprops=flierprops, ) - # Baseline median reference line spanning the variant blocks if len(base_aucs): ax.axhline(np.median(base_aucs), color=C_BASE, linewidth=1.2, linestyle="--", alpha=0.55, label="Baseline median") - # Dividers between sections (vertical light lines) div1 = (section1_right + section2_left) / 2 div2 = (section2_right + section3_left) / 2 for d in (div1, div2): ax.axvline(d, color="#aaa", linewidth=0.7, alpha=0.65, linestyle="-") - # Section labels just above each block y_band = 1.02 section_centers = [ - ((pos[0] + section1_right) / 2, "Basic backbones (img-only, single)"), - ((section2_left + section2_right) / 2, "ResNet50 variations"), - ((section3_left + section3_right) / 2, "Baseline"), + ((pos[0] + section1_right) / 2, "Basic backbones (img-only, single)"), + ((section2_left + section2_right) / 2, "ResNet50 variations"), + ((section3_left + section3_right) / 2, "Baseline"), ] for cx, txt in section_centers: ax.text(cx, y_band, txt, ha="center", va="bottom", fontsize=10, color="#333", fontweight="bold", transform=ax.get_xaxis_transform()) - # X-tick labels (with p-values vs baseline beneath each variant box) tick_labels = [] for lbl, aucs, color in zip(all_labels, all_aucs, all_colors): if color == C_BASE or not len(aucs) or not len(base_aucs): diff --git a/v4/figures/F3_hyperfeature_ablation.py b/v4/figures/F3_hyperfeature_ablation.py index 34de10f..3162c0c 100644 --- a/v4/figures/F3_hyperfeature_ablation.py +++ b/v4/figures/F3_hyperfeature_ablation.py @@ -58,9 +58,9 @@ SEV_COLORS = { "unknown": C_UNKNOWN, } -SEV_ORDER = ["normal", "unknown", "early", "moderate", "severe"] -SEV_ALPHA = {"normal": 0.40, "unknown": 0.35, "early": 0.55, "moderate": 0.70, "severe": 0.85} -SEV_SIZE = {"normal": 6, "unknown": 6, "early": 8, "moderate": 10, "severe": 12} +SEV_ORDER = ["severe", "moderate", "unknown", "early", "normal"] +SEV_ALPHA = {"normal": 0.55, "unknown": 0.55, "early": 0.55, "moderate": 0.55, "severe": 0.55} +SEV_SIZE = {"normal": 8, "unknown": 8, "early": 8, "moderate": 8, "severe": 8} # ── Per-panel definitions: (label, results dir, eval_stage) ────────────────── # Top row: single-modality reference runs diff --git a/v4/figures/F3b_anonymous_cv_roc.py b/v4/figures/F3b_anonymous_cv_roc.py new file mode 100644 index 0000000..62cc117 --- /dev/null +++ b/v4/figures/F3b_anonymous_cv_roc.py @@ -0,0 +1,155 @@ +"""F3b - Hadamard L1 fusion ROC: patient-grouped vs anonymous CV (overlay). + +Single panel overlaying the ROC of the same single-eye img+cd Hadamard L1 +fusion configuration evaluated under patient-grouped 5-fold CV (the headline +protocol) and patient-anonymous 5-fold CV (the prevailing benchmark protocol). +Both runs use matched seeds, matched backbones, matched training schedules; +only the fold-grouping rule changes. + +For each configuration the figure shows: + - per-fold-rep ROC curves as faint coloured lines (50 curves per condition) + - mean ROC across fold-reps with a shaded SD band + - the rep-mean test AUC ± SD in a corner annotation + +Re-run: + python -m v4.figures.F3b_anonymous_cv_roc +""" +from __future__ import annotations + +import warnings +warnings.filterwarnings("ignore") + +from pathlib import Path + +import h5py +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from sklearn.metrics import roc_auc_score, roc_curve + +from v4.figures.util.loaders import RESULTS_ROOT + + +OUT = Path(__file__).parent / "output" / "F3b_anonymous_cv_roc.png" + +CONDITIONS = [ + ("Patient-grouped CV", + RESULTS_ROOT / "refuge_v2m_baseline" / "ensemble_single_refugelike", + "nt", "#1f6fb0"), + ("Anonymous CV", + RESULTS_ROOT / "backbone_replication" / "anonymous_cv_ensemble_single_refugelike", + "nt", "#c44e52"), +] + +# Common FPR grid for per-fold-rep ROC interpolation +FPR_GRID = np.linspace(0.0, 1.0, 201) + + +def collect_per_foldrep(run_dir: Path, eval_stage: str): + """Return list of (y_true, y_score) tuples, one per (rep, fold).""" + per: list[tuple[np.ndarray, np.ndarray]] = [] + for rep in sorted(run_dir.glob("rep*")): + fp = next(iter(rep.rglob("predictions.h5")), None) + if fp is None: + continue + with h5py.File(fp, "r") as f: + if eval_stage not in f: + continue + grp = f[eval_stage] + logits = grp["logits"][:] + y_true = grp["y_true"][:].astype(int) + split = grp["split"][:] + n_folds, n_epochs, _, n_heads, n_outputs = logits.shape + if n_outputs != 2: + continue + ep, head = n_epochs - 1, n_heads - 1 + for fold in range(n_folds): + labels = np.array( + [s.decode() if isinstance(s, bytes) else str(s) for s in split[fold]] + ) + test_mask = labels == "test" + if not test_mask.any(): + continue + lg = logits[fold, ep, test_mask, head, :] + e = np.exp(lg - lg.max(axis=1, keepdims=True)) + p = e / e.sum(axis=1, keepdims=True) + y = y_true[test_mask] + s = p[:, 1] + per.append((y, s)) + return per + + +def interp_tpr(y: np.ndarray, s: np.ndarray) -> np.ndarray: + if len(np.unique(y)) < 2: + return np.full_like(FPR_GRID, np.nan, dtype=float) + fpr, tpr, _ = roc_curve(y, s) + return np.interp(FPR_GRID, fpr, tpr) + + +def render() -> None: + fig, ax = plt.subplots(figsize=(7.4, 6.4)) + + # Diagonal reference first so it sits behind everything + ax.plot([0, 1], [0, 1], color="#aaa", linewidth=0.8, linestyle=":", zorder=1) + + annotations = [] + + for label, run_dir, stage, color in CONDITIONS: + per = collect_per_foldrep(run_dir, stage) + if not per: + continue + + tprs = np.array([interp_tpr(y, s) for (y, s) in per]) + valid = ~np.isnan(tprs).any(axis=1) + tprs = tprs[valid] + + per_aucs = np.array( + [roc_auc_score(y, s) for (y, s) in per if len(np.unique(y)) >= 2] + ) + rep_means = (per_aucs.reshape(-1, 5).mean(axis=1) + if len(per_aucs) % 5 == 0 else per_aucs) + + # Per-fold-rep curves + for (y, s) in per: + if len(np.unique(y)) < 2: + continue + fpr, tpr, _ = roc_curve(y, s) + ax.plot(fpr, tpr, color=color, linewidth=0.4, alpha=0.13, zorder=2) + + # Mean ± SD band + mean_tpr = tprs.mean(axis=0) + sd_tpr = tprs.std(axis=0) + ax.fill_between( + FPR_GRID, np.clip(mean_tpr - sd_tpr, 0, 1), + np.clip(mean_tpr + sd_tpr, 0, 1), + color=color, alpha=0.20, zorder=3, + ) + ax.plot( + FPR_GRID, mean_tpr, color=color, linewidth=2.2, zorder=4, + label=f"{label} (AUC = {rep_means.mean():.3f} ± {rep_means.std():.3f})", + ) + + annotations.append((label, rep_means)) + + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.set_xlabel("False Positive Rate", fontsize=11) + ax.set_ylabel("True Positive Rate", fontsize=11) + ax.set_title( + "Hadamard L1 fusion ROC under matched architecture,\n" + "patient-grouped vs anonymous cross-validation", + fontsize=12, fontweight="bold", + ) + ax.grid(alpha=0.25, linestyle="--") + ax.legend(loc="lower right", fontsize=10, framealpha=0.94) + + fig.tight_layout() + OUT.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(OUT, dpi=180, bbox_inches="tight") + plt.close(fig) + print(f"saved {OUT}") + + +if __name__ == "__main__": + render() diff --git a/v4/figures/F4_bilateral.py b/v4/figures/F4_bilateral.py index 85b6590..f199ead 100644 --- a/v4/figures/F4_bilateral.py +++ b/v4/figures/F4_bilateral.py @@ -58,15 +58,15 @@ SEV_LABELS = { "severe": "Glaucoma — severe (VF_MD < −12)", "unknown": "Glaucoma — VF_MD not recorded", } -SEV_ORDER = ["normal", "unknown", "early", "moderate", "severe"] +SEV_ORDER = ["severe", "moderate", "unknown", "early", "normal"] SEV_ALPHA = { - "normal": 0.40, - "unknown": 0.35, + "normal": 0.55, + "unknown": 0.55, "early": 0.55, - "moderate": 0.70, - "severe": 0.85, + "moderate": 0.55, + "severe": 0.55, } -SEV_SIZE = {"normal": 6, "unknown": 6, "early": 8, "moderate": 10, "severe": 12} +SEV_SIZE = {"normal": 8, "unknown": 8, "early": 8, "moderate": 8, "severe": 8} # Panel grid: [row][col] = (label, run_dir, eval_stage) GRID = [ diff --git a/v4/figures/F6_regression.py b/v4/figures/F6_regression.py index a58057a..81f842a 100644 --- a/v4/figures/F6_regression.py +++ b/v4/figures/F6_regression.py @@ -27,7 +27,10 @@ from sklearn.metrics import roc_curve, roc_auc_score from v4.figures.util.loaders import RESULTS_ROOT OUT = Path(__file__).parent / "output" / "F6_regression.png" -RUN_DIR = RESULTS_ROOT / "reg_head" / "baseline_reg_nt50" +# Points at the post-fix run that stores VF_MD as float64. The earlier +# baseline_reg_nt50 run stored y_true as int64, silently rounding the +# regression targets; do not mix the two. +RUN_DIR = RESULTS_ROOT / "reg_head" / "baseline_reg_nt50_floaty" # Prediction-side bin boundaries NP_THRESH = -1.097 # mean of measured-healthy MD diff --git a/v4/figures/F8_attention_combined.py b/v4/figures/F8_attention_combined.py new file mode 100644 index 0000000..745cc68 --- /dev/null +++ b/v4/figures/F8_attention_combined.py @@ -0,0 +1,72 @@ +"""F8 combined explainability panel: disc-centred attention + quadrant breakdown. + +Composes a single A/B figure from two existing renderings: + A: ``F8_gradcam/disc_attention_detail.png`` — 2x2 grid of mean Grad-CAM + heatmaps for {correct, incorrect} x {Normal, Glaucoma} cells, with the + mean disc boundary annotated as a dashed circle. + B: ``F8_quadrant_attention.png`` — grouped-bar chart of mean full-image + Grad-CAM fraction per optic-disc quadrant, by cell. + +Both source panels are produced by ``v4.figures.F8_explainability`` and +``v4.figures.F8_quadrant_plot`` respectively; this script just stitches the +two PNGs into a single combined figure with A/B subfigure labels. + +Re-run: + python -m v4.figures.F8_attention_combined +""" +from __future__ import annotations + +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from PIL import Image + + +SRC_A = Path(__file__).parent / "output" / "F8_gradcam" / "disc_attention_detail.png" +SRC_B = Path(__file__).parent / "output" / "F8_quadrant_attention.png" +OUT = Path(__file__).parent / "output" / "F8_attention_combined.png" + + +def render() -> None: + for p in (SRC_A, SRC_B): + if not p.exists(): + raise SystemExit( + f"Source panel missing: {p}\n" + "Run F8_explainability (for A) and F8_quadrant_plot (for B) first." + ) + + img_a = Image.open(SRC_A) + img_b = Image.open(SRC_B) + + # Stack vertically: A on top (square), B below (wider). + fig = plt.figure(figsize=(13.0, 13.6)) + gs = fig.add_gridspec( + 2, 1, + height_ratios=[img_a.size[1] / img_a.size[0], + img_b.size[1] / img_b.size[0] * 13.0 / 13.0], + hspace=0.06, + ) + + ax_a = fig.add_subplot(gs[0]) + ax_a.imshow(img_a) + ax_a.axis("off") + ax_a.text(-0.01, 1.01, "A", transform=ax_a.transAxes, + ha="left", va="bottom", fontsize=22, fontweight="bold") + + ax_b = fig.add_subplot(gs[1]) + ax_b.imshow(img_b) + ax_b.axis("off") + ax_b.text(-0.01, 1.01, "B", transform=ax_b.transAxes, + ha="left", va="bottom", fontsize=22, fontweight="bold") + + fig.tight_layout() + OUT.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(OUT, dpi=180, bbox_inches="tight") + plt.close(fig) + print(f"saved {OUT}") + + +if __name__ == "__main__": + render() diff --git a/v4/figures/F8_explainability.py b/v4/figures/F8_explainability.py index 3df6e88..435fba9 100644 --- a/v4/figures/F8_explainability.py +++ b/v4/figures/F8_explainability.py @@ -1,8 +1,10 @@ -"""F8 - Explainability figures from the V2-M checkpointed v4 run. +"""F8 - Explainability figures from the R50 checkpointed v4 run. Sources predictions and Grad-CAM panels exclusively from the -``experiments/explainability/ensemble_v2m_ckpt`` run (img+cd ensemble with the -refuge_efficientnet_v2_m backbone, save_checkpoints=true). +``experiments/explainability/ensemble_refugelike_ckpt`` run (img+cd ensemble +with the refugelike R50 backbone, save_checkpoints=true). rep00 (seed=1234) +is the single rep used for the figure; matches the headline configuration in +section 3. GradCAM machinery lives in ``v4.classes.accessory.explainability``; PAPILA specific knowledge (disc contour rasterisation, OS→OD orientation flip) is @@ -33,7 +35,7 @@ from v4.figures.util.loaders import REPO_ROOT OUT_DIR = Path(__file__).parent / "output" GRADCAM_DIR = OUT_DIR / "F8_gradcam" -V4_CKPT_RUN = REPO_ROOT / "v4" / "results" / "experiments" / "explainability" / "ensemble_v2m_ckpt" / "binary" +V4_CKPT_RUN = REPO_ROOT / "v4" / "results" / "experiments" / "explainability" / "ensemble_refugelike_ckpt" / "rep00" / "binary" LABEL_NAMES = {0: "Normal", 1: "Glaucoma"} EVENT_ORDER = [ @@ -43,11 +45,11 @@ EVENT_ORDER = [ ] EVENT_LABELS = { "full_correction": "Both wrong -> fused right", - "img_assist": "Image right, MD wrong", - "md_assist": "MD right, image wrong", + "img_assist": "Image right, clinical wrong", + "md_assist": "Clinical right, image wrong", "full_error": "Both right -> fused wrong", - "img_drag": "MD right, image wrong -> fused wrong", - "md_drag": "Image right, MD wrong -> fused wrong", + "img_drag": "Clinical right, image wrong -> fused wrong", + "md_drag": "Image right, clinical wrong -> fused wrong", "concordant_correct": "All correct", "concordant_wrong": "All wrong", } @@ -339,6 +341,8 @@ def make_fusion_event_panel(split: str = "test") -> None: ax.invert_yaxis() ax.set_xlabel("Count") ax.set_title("Fusion event taxonomy", fontsize=10, fontweight="bold") + max_count = max(display_counts[k] for k in bars) if bars else 0 + ax.set_xlim(0, max_count * 1.10 + 1) for yi, k in enumerate(bars): ax.text(display_counts[k] + 0.8, yi, str(int(display_counts[k])), va="center", fontsize=8) @@ -346,7 +350,7 @@ def make_fusion_event_panel(split: str = "test") -> None: ax = fig.add_subplot(gs[0, 1]) per_fold = pd.DataFrame( { - "fold": [f"{r}/{f}" for (r, f), _ in fold_groups], + "fold": [str(f) for (_, f), _ in fold_groups], "positive": [sum((g["event_type"] == k).sum() for k in positive_keys) for _, g in fold_groups], "negative": [sum((g["event_type"] == k).sum() for k in negative_keys) for _, g in fold_groups], } @@ -437,10 +441,6 @@ def make_fusion_event_panel(split: str = "test") -> None: ax.legend(handles=point_handles + shade_handles, ncol=5, fontsize=7, loc="upper center", bbox_to_anchor=(0.5, -0.14), frameon=False) - fig.suptitle( - f"S8a - Checkpoint Fusion Events ({split}; AUC={auc:.3f}, n={len(df)})", - fontsize=12, fontweight="bold", - ) out = OUT_DIR / "S8a_comparison_panel.png" fig.savefig(out, dpi=180, bbox_inches="tight") plt.close(fig) @@ -450,7 +450,7 @@ def make_fusion_event_panel(split: str = "test") -> None: def make_clinical_importance(n_permutations: int = 30, seed: int = 0) -> None: """S8e clinical permutation importance via the cd-tower → cd_aux head. - Isolates the clinical-only prediction path at the V2-M ckpt run, then + Isolates the clinical-only prediction path at the R50 ckpt run, then column-shuffles the encoded clinical vector to measure per-feature AUC drop. Per-original-column grouping comes from ``ClinicalDataView.feature_groups`` (one-hot encoded dims for a @@ -590,7 +590,7 @@ def make_clinical_importance(n_permutations: int = 30, seed: int = 0) -> None: ax.set_xlabel("Mean AUC drop on shuffling (averaged across folds)", fontsize=10) ax.axvline(0, color="black", linewidth=0.7) ax.set_title( - f"S8e — Clinical permutation importance (cd-only head, V2-M ckpt run)\n" + f"S8e — Clinical permutation importance (cd-only head, R50 ckpt run)\n" f"baseline AUC = {baseline_mean:.3f}; n_permutations = {n_permutations}", fontsize=10, fontweight="bold", ) @@ -692,6 +692,60 @@ def _disc_centred_patch_array( return patch_out.astype(np.float32), out * disc_r / (2 * half) +QUAD_ORDER = ("ST", "SN", "IT", "IN") # superotemporal, superonasal, inferotemporal, inferonasal + + +def _quadrant_fractions( + cam: np.ndarray, + disc_mask: np.ndarray, + *, + peri_inner: float = 1.0, + peri_outer: float = 2.0, +) -> tuple[dict[str, float], dict[str, float], dict[str, float]] | tuple[None, None, None]: + """Per-quadrant Grad-CAM fractions in OD-oriented coordinates, for three + region scopes. + + Quadrant boundaries are the disc-mask centroid (cx, cy). In the OD-oriented + frame nasal is left (x < cx) and temporal is right (x > cx); superior is + top (y < cy) and inferior is bottom (y > cy): + + ST = x > cx, y < cy + SN = x < cx, y < cy + IT = x > cx, y > cy + IN = x < cx, y > cy + + Three region scopes are returned: + disc_q : fractions of CAM intensity that fall inside the GT disc mask + peri_q : fractions inside a peri-disc annulus of disc-radius units + (peri_inner to peri_outer, default 1x-2x), excluding the disc + full_q : fractions over the entire image + Each dict sums to 1 (within floating-point error). Returns (None, None, None) + if the disc mask is empty. + """ + if disc_mask is None or disc_mask.sum() == 0: + return None, None, None + ys, xs = np.where(disc_mask) + cy = float(ys.mean()); cx = float(xs.mean()) + disc_r = float(np.sqrt(disc_mask.sum() / np.pi)) + h, w = cam.shape + yy, xx = np.mgrid[0:h, 0:w] + dist = np.sqrt((xx - cx) ** 2 + (yy - cy) ** 2) + peri_mask = (dist >= peri_inner * disc_r) & (dist <= peri_outer * disc_r) & ~disc_mask + quads = { + "ST": (xx > cx) & (yy < cy), + "SN": (xx < cx) & (yy < cy), + "IT": (xx > cx) & (yy > cy), + "IN": (xx < cx) & (yy > cy), + } + disc_total = float(cam[disc_mask].sum()) + 1e-8 + peri_total = float(cam[peri_mask].sum()) + 1e-8 + full_total = float(cam.sum()) + 1e-8 + disc_q = {k: float(cam[disc_mask & q].sum()) / disc_total for k, q in quads.items()} + peri_q = {k: float(cam[peri_mask & q].sum()) / peri_total for k, q in quads.items()} + full_q = {k: float(cam[q].sum()) / full_total for k, q in quads.items()} + return disc_q, peri_q, full_q + + def _annotate_nasal_temporal(ax, *, fontsize: int = 9, color: str = "white", pad: float = 2.5) -> None: """Label the disc-side (nasal) and macula-side (temporal) edges of an @@ -868,6 +922,10 @@ def _make_oriented_gradcam(n_grid: int = 16, alpha: float = 0.45, disc_patch_count: dict[tuple[str, str], int] = {} disc_radius_sum: dict[tuple[str, str], float] = {} disc_frac_sum: dict[tuple[str, str], float] = {} + # Per-eye quadrant fractions for three region scopes; aggregated per cell. + quad_disc_list: dict[tuple[str, str], list[dict[str, float]]] = {} + quad_peri_list: dict[tuple[str, str], list[dict[str, float]]] = {} + quad_full_list: dict[tuple[str, str], list[dict[str, float]]] = {} examples: dict[tuple[str, str], tuple[np.ndarray, int, str, int]] = {} fold_range = range(cfg.get("folds", 5)) @@ -927,6 +985,11 @@ def _make_oriented_gradcam(n_grid: int = 16, alpha: float = 0.45, disc_frac_sum[key] = disc_frac_sum.get(key, 0.0) + float( cam_np[roi_mask].sum() / (cam_np.sum() + 1e-8) ) + disc_q, peri_q, full_q = _quadrant_fractions(cam_np, roi_mask) + if disc_q is not None: + quad_disc_list.setdefault(key, []).append(disc_q) + quad_peri_list.setdefault(key, []).append(peri_q) + quad_full_list.setdefault(key, []).append(full_q) if key not in examples: ov = overlay_gradcam(pil_oriented, cam_np, alpha) ov_small = np.array(ov.resize(cam_np.shape[::-1], Image.BILINEAR)) @@ -1020,6 +1083,32 @@ def _make_oriented_gradcam(n_grid: int = 16, alpha: float = 0.45, if mean_patches: _make_oriented_disc_detail(mean_patches, examples, GRADCAM_DIR / "disc_attention_detail.png") + # Per-quadrant CAM fractions (within-disc and full-image scopes). + # One CSV row per (class, outcome, scope, quadrant) cell, with mean and SD + # computed over the per-eye fractions in that cell. + rows = [] + for key in sorted(quad_disc_list.keys()): + cls_name, outcome = key + n_eyes = len(quad_disc_list[key]) + for scope_label, scope_list in (("disc", quad_disc_list[key]), + ("peri", quad_peri_list[key]), + ("full", quad_full_list[key])): + for q in QUAD_ORDER: + vals = np.array([d[q] for d in scope_list], dtype=np.float64) + rows.append({ + "class": cls_name, + "outcome": outcome, + "scope": scope_label, + "quadrant": q, + "n_eyes": n_eyes, + "mean": float(vals.mean()), + "sd": float(vals.std(ddof=1)) if n_eyes > 1 else float("nan"), + }) + if rows: + out_csv = OUT_DIR / "F8_quadrant_fractions.csv" + pd.DataFrame(rows).to_csv(out_csv, index=False) + print(f"saved quadrant fractions: {out_csv}") + def make_oriented_gradcam(n_grid: int = 16, alpha: float = 0.45, target_class: int | None = None) -> None: diff --git a/v4/figures/F8_quadrant_plot.py b/v4/figures/F8_quadrant_plot.py new file mode 100644 index 0000000..bc562f7 --- /dev/null +++ b/v4/figures/F8_quadrant_plot.py @@ -0,0 +1,123 @@ +"""F8 quadrant attention chart (full-image scope). + +Reads the per-cell quadrant fractions produced by +``v4.figures.F8_explainability`` (see CSV at +``output/F8_quadrant_fractions.csv``) and renders a single-panel grouped bar +chart of full-image Grad-CAM intensity by optic-disc quadrant. Within-disc +and peri-disc scopes are present in the CSV but are not plotted here; see +the CSV for the per-cell numbers if you need them. + +Re-run: + python -m v4.figures.F8_quadrant_plot +""" +from __future__ import annotations + +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +CSV = Path(__file__).parent / "output" / "F8_quadrant_fractions.csv" +OUT = Path(__file__).parent / "output" / "F8_quadrant_attention.png" + +QUAD_ORDER = ("ST", "SN", "IT", "IN") +QUAD_LABEL = { + "ST": "Superotemporal", + "SN": "Superonasal", + "IT": "Inferotemporal", + "IN": "Inferonasal", +} + +CELL_ORDER = [ + ("Normal", "correct"), + ("Normal", "incorrect"), + ("Glaucoma", "correct"), + ("Glaucoma", "incorrect"), +] +CELL_COLOR = { + ("Normal", "correct"): "#3B6FB5", + ("Normal", "incorrect"): "#8BB0DA", + ("Glaucoma", "correct"): "#c44e52", + ("Glaucoma", "incorrect"): "#e6a3a4", +} +CELL_LABEL = { + ("Normal", "correct"): "Normal correct", + ("Normal", "incorrect"): "Normal incorrect", + ("Glaucoma", "correct"): "Glaucoma correct", + ("Glaucoma", "incorrect"): "Glaucoma incorrect", +} + + +def render() -> None: + if not CSV.exists(): + raise SystemExit( + f"CSV {CSV} not found. Run `python -m v4.figures.F8_explainability " + "--only-gradcam --run-gradcam` first." + ) + df = pd.read_csv(CSV) + df_full = df[df["scope"] == "full"] + + fig, ax = plt.subplots(figsize=(9.6, 6.0)) + + n_quad = len(QUAD_ORDER) + n_cell = len(CELL_ORDER) + bar_w = 0.18 + x = np.arange(n_quad, dtype=float) + + for i, cell in enumerate(CELL_ORDER): + means, sds = [], [] + n_eyes = None + for q in QUAD_ORDER: + row = df_full[ + (df_full["class"] == cell[0]) + & (df_full["outcome"] == cell[1]) + & (df_full["quadrant"] == q) + ] + means.append(float(row["mean"].iloc[0]) if len(row) else float("nan")) + sds.append(float(row["sd"].iloc[0]) if len(row) else float("nan")) + if n_eyes is None and len(row): + n_eyes = int(row["n_eyes"].iloc[0]) + + offsets = (i - (n_cell - 1) / 2.0) * bar_w + bars = ax.bar( + x + offsets, means, width=bar_w, + yerr=sds, capsize=2, + color=CELL_COLOR[cell], alpha=0.92, + edgecolor="black", linewidth=0.6, + label=f"{CELL_LABEL[cell]} (n = {n_eyes})", + error_kw=dict(ecolor="#444", linewidth=0.8, capthick=0.8), + ) + for rect, m in zip(bars, means): + if np.isnan(m): + continue + ax.text(rect.get_x() + rect.get_width() / 2, m + 0.005, + f"{m:.2f}", ha="center", va="bottom", + fontsize=8.5, color="#222") + + ax.set_xticks(x) + ax.set_xticklabels([QUAD_LABEL[q] for q in QUAD_ORDER], fontsize=11) + ax.set_ylabel("Mean fraction of full-image Grad-CAM intensity", fontsize=11) + ax.set_title( + "Image-tower Grad-CAM by optic-disc quadrant (OD-oriented; centroid-split)", + fontsize=12.5, fontweight="bold", + ) + ax.grid(axis="y", alpha=0.3, linestyle="--") + ax.set_ylim(0, max(ax.get_ylim()[1], 0.7)) + ax.axhline(0.25, color="#888", linestyle=":", linewidth=0.8, alpha=0.6, zorder=0) + + ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.10), ncol=4, + framealpha=0.94, fontsize=10) + + fig.tight_layout(rect=(0, 0.04, 1, 0.98)) + OUT.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(OUT, dpi=180, bbox_inches="tight") + plt.close(fig) + print(f"saved {OUT}") + + +if __name__ == "__main__": + render() diff --git a/v4/figures/S2_variance_decomposition.py b/v4/figures/S2_variance_decomposition.py new file mode 100644 index 0000000..c2fa0bf --- /dev/null +++ b/v4/figures/S2_variance_decomposition.py @@ -0,0 +1,277 @@ +"""S2 - Variance decomposition of the four-bridge L1 fusion comparison. + +Two-panel supplementary figure summarising the variance decomposition +reported alongside section 3.2 of the manuscript. + +Panel A (left): paired-line plot. + x-axis : the 4 bridge variants (Concat, Pairwise, Gated, Hadamard) + y-axis : eye-level test AUC + each line : one fold-rep, connecting that fold-rep's 4 bridge AUCs + overlay : per-bridge boxplot showing the marginal AUC distribution + annotation : pooled across-architecture and across-fold-rep SDs, + and the SD ratio + +Panel B (right): centered-offset KDEs. + x-axis : AUC offset from grouping mean (centered at 0) + y-axis : density + four coloured curves : per-bridge fold-rep distributions (50 fold-reps + per bridge, centered by subtracting each bridge's + own mean) + dashed dark curve : architectural offset distribution (200 values, + centered by subtracting each fold-rep's mean) + annotation : variance ratio + +Read directly from v4/results/experiments/{phase3_v4,refuge_v2m_baseline}. + +Re-run: + python -m v4.figures.S2_variance_decomposition +""" +from __future__ import annotations + +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 scipy.stats import gaussian_kde + +from v4.figures.util.loaders import RESULTS_ROOT + + +OUT = Path(__file__).parent / "output" / "S2_variance_decomposition.png" + + +# Bridge label, results dir relative to experiments/, and stage key +BRIDGES = [ + ("Concat", "phase3_v4/single_bcd_concat"), + ("Pairwise", "phase3_v4/single_bcd_pairwise"), + ("Gated", "phase3_v4/single_bcd_gated"), + ("Hadamard", "refuge_v2m_baseline/ensemble_single_refugelike"), +] +STAGE_KEY = "nt_test_auc" + + +# Panel-A colour palette (paired lines + boxplots) +C_LINE = "#1f6fb0" # single blue for all paired-cell lines +C_LINE_ALPHA = 0.22 +C_MARKER = "#1f6fb0" +C_MEDIAN = "#c44e52" # red box median line +C_BOX_FILL = "#dbe6f0" # pale blue box fill + +# Panel-B colour palette (per-bridge KDEs) +C_ARCH = "#222" # dark grey for the architectural curve +BRIDGE_COLORS = { + "Concat": "#7f7f7f", + "Pairwise": "#ff7f0e", + "Gated": "#2ca02c", + "Hadamard": "#1f77b4", +} + + +# ── Data loading + decomposition ──────────────────────────────────────────── + +def collect() -> pd.DataFrame: + rows = [] + for label, rel in BRIDGES: + root = RESULTS_ROOT / rel + for s in sorted(root.glob("rep*/binary/summary.json")): + rep = int(s.parents[1].name.replace("rep", "")) + d = json.loads(s.read_text()) + for fr in d.get("fold_results", []): + v = fr.get(STAGE_KEY) + if v is None or not np.isfinite(v): + continue + rows.append({ + "bridge": label, + "rep": rep, + "fold": int(fr["fold"]), + "auc": float(v), + }) + return pd.DataFrame(rows) + + +def decompose(df: pd.DataFrame) -> tuple[float, float, float, dict[str, float]]: + """Return (arch_sd_pooled, fold_sd_pooled, ratio, per_bridge_sd).""" + cell_var = df.groupby(["rep", "fold"])["auc"].var(ddof=1) + bridge_var = df.groupby("bridge")["auc"].var(ddof=1) + arch_sd = float(np.sqrt(cell_var.mean())) + fold_sd = float(np.sqrt(bridge_var.mean())) + ratio = fold_sd / arch_sd if arch_sd > 0 else float("inf") + per_bridge_sd = {b: float(np.sqrt(v)) for b, v in bridge_var.items()} + return arch_sd, fold_sd, ratio, per_bridge_sd + + +# ── Panel A: paired-line plot ─────────────────────────────────────────────── + +def draw_panel_a(ax, df: pd.DataFrame, + arch_sd: float, fold_sd: float, ratio: float, + n_cells: int) -> None: + bridge_order = [b for b, _ in BRIDGES] + pivot = df.pivot_table( + index=["rep", "fold"], columns="bridge", values="auc" + )[bridge_order] + x_positions = np.arange(len(bridge_order), dtype=float) + + # Paired cell lines: one per fold-rep + for _, row in pivot.iterrows(): + if row.isna().any(): + continue + ax.plot( + x_positions, row.values, color=C_LINE, alpha=C_LINE_ALPHA, + linewidth=0.9, marker="o", markersize=2.0, + markerfacecolor=C_MARKER, markeredgecolor="none", zorder=2, + ) + + # Per-bridge boxplot + box_data = [pivot[b].dropna().values for b in bridge_order] + ax.boxplot( + box_data, + positions=x_positions, + widths=0.32, + patch_artist=True, + manage_ticks=False, + zorder=3, + boxprops=dict(facecolor=C_BOX_FILL, edgecolor="black", + linewidth=1.0, alpha=0.85), + whiskerprops=dict(color="black", linewidth=0.9), + capprops=dict(color="black", linewidth=0.9), + medianprops=dict(color=C_MEDIAN, linewidth=1.8), + flierprops=dict(marker="", markersize=0), + ) + + ax.set_xticks(x_positions) + ax.set_xticklabels(bridge_order, fontsize=11) + ax.set_xlabel("L1 fusion bridge", fontsize=11) + ax.set_ylabel("Eye-level test AUC", fontsize=11) + ax.set_xlim(-0.5, len(bridge_order) - 0.5) + ax.grid(axis="y", alpha=0.3, linestyle="--") + + txt = ( + f"n = {n_cells} fold-rep AUC values | variance ratio {ratio:.2f}\n" + f" across-architecture SD = {arch_sd:.3f}\n" + f" across-fold-rep SD = {fold_sd:.3f}" + ) + ax.text( + 0.985, 0.025, txt, transform=ax.transAxes, + ha="right", va="bottom", fontsize=9.5, family="monospace", + bbox=dict(boxstyle="round,pad=0.5", facecolor="white", + edgecolor="#888", alpha=0.92), + ) + + +# ── Panel B: centered-offset KDE curves ───────────────────────────────────── + +def draw_panel_b(ax, df: pd.DataFrame, + arch_sd: float, fold_sd: float, ratio: float, + per_bridge_sd: dict[str, float]) -> None: + df = df.copy() + fold_rep_mean = df.groupby(["rep", "fold"])["auc"].transform("mean") + bridge_mean = df.groupby("bridge")["auc"].transform("mean") + df["arch_offset"] = df["auc"] - fold_rep_mean + df["fold_offset"] = df["auc"] - bridge_mean + + all_offsets = np.concatenate( + [df["arch_offset"].values, df["fold_offset"].values] + ) + x_max = float(np.abs(all_offsets).max()) * 1.10 + xs = np.linspace(-x_max, x_max, 600) + + # Per-bridge fold-rep offset curves + for b in [name for name, _ in BRIDGES]: + offs = df.loc[df["bridge"] == b, "fold_offset"].values + kde = gaussian_kde(offs) + y = kde(xs) + sd = per_bridge_sd[b] + ax.plot(xs, y, color=BRIDGE_COLORS[b], linewidth=1.6, alpha=0.92, + zorder=3, + label=f"{b} SD = {sd:.3f}") + + # Architectural offset curve (pooled) + arch_offsets = df["arch_offset"].values + kde_arch = gaussian_kde(arch_offsets) + y_arch = kde_arch(xs) + ax.fill_between(xs, y_arch, color=C_ARCH, alpha=0.18, zorder=2) + ax.plot(xs, y_arch, color=C_ARCH, linewidth=2.2, linestyle="--", zorder=4, + label=f"Architectural SD = {arch_sd:.3f}") + + ax.axvline(0, color="#666", linewidth=0.8, linestyle=":", + alpha=0.6, zorder=0) + ax.set_xlabel("AUC offset from grouping mean", fontsize=11) + ax.set_ylabel("Probability density", fontsize=11) + ax.set_xlim(-x_max, x_max) + ax.set_yticklabels([]) + ax.tick_params(axis="y", which="both", left=True, labelleft=False) + ax.grid(axis="y", alpha=0.25, linestyle="--") + + # Headroom on the y-axis so the variance-ratio box does not crowd the + # architectural-curve peak. + ymin, ymax = ax.get_ylim() + ax.set_ylim(0, ymax * 1.10) + + # Legend below the top so it clears the variance-ratio annotation. + ax.legend( + loc="upper left", bbox_to_anchor=(0.0, 0.82), + fontsize=9, framealpha=0.92, + ) + + txt = f"variance ratio fold-SD / arch-SD = {ratio:.2f}" + ax.text( + 0.985, 0.975, txt, transform=ax.transAxes, + ha="right", va="top", fontsize=10.5, family="monospace", + bbox=dict(boxstyle="round,pad=0.45", facecolor="white", + edgecolor="#888", alpha=0.92), + ) + + +# ── Combined render ───────────────────────────────────────────────────────── + +def render() -> None: + df = collect() + if df.empty: + print("No data collected; check the source paths.") + return + + n_cells = df.groupby(["rep", "fold"]).ngroups + n_bridges = df["bridge"].nunique() + arch_sd, fold_sd, ratio, per_bridge_sd = decompose(df) + + print(f"Collected {len(df)} observations ({n_bridges} bridges x {n_cells} cells)") + print(f" across-architecture SD (within cell) : {arch_sd:.4f}") + print(f" across-fold-rep SD (within bridge) : {fold_sd:.4f}") + print(f" ratio fold-SD / arch-SD : {ratio:.2f}x") + print("Per-bridge fold-rep SDs:") + for b in [name for name, _ in BRIDGES]: + print(f" {b:<10s} SD = {per_bridge_sd[b]:.4f}") + + fig, (axA, axB) = plt.subplots( + nrows=1, ncols=2, figsize=(16.0, 6.0), + gridspec_kw=dict(wspace=0.22), + ) + + draw_panel_a(axA, df, arch_sd, fold_sd, ratio, n_cells) + draw_panel_b(axB, df, arch_sd, fold_sd, ratio, per_bridge_sd) + + # Subfigure labels + for ax, label in ((axA, "A"), (axB, "B")): + ax.text( + -0.07, 1.03, label, transform=ax.transAxes, + ha="left", va="bottom", fontsize=15, fontweight="bold", + ) + + fig.suptitle( + "Variance decomposition: fold-assignment noise vs L1 bridge choice", + fontsize=13.0, fontweight="bold", y=1.00, + ) + fig.tight_layout(rect=(0, 0, 1, 0.97)) + + OUT.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(OUT, dpi=180, bbox_inches="tight") + plt.close(fig) + print(f"saved {OUT}") + + +if __name__ == "__main__": + render() diff --git a/v4/figures/S3_variance_distributions.py b/v4/figures/S3_variance_distributions.py new file mode 100644 index 0000000..7333536 --- /dev/null +++ b/v4/figures/S3_variance_distributions.py @@ -0,0 +1,162 @@ +"""S3 - Centered-offset distributions of architectural vs fold-rep variance. + +Companion to S2, presenting the same variance decomposition as two +overlapping KDE curves on a common centered axis. + +For each of the 200 (fold-rep, bridge) AUC observations, compute two +mean-centered offsets: + + architectural offset = AUC - mean(AUC over the 4 bridges in that fold-rep) + fold-rep offset = AUC - mean(AUC over the 50 fold-reps for that bridge) + +Both sets have 200 values, both are centered at 0 by construction, and the +spread of each distribution corresponds directly to one of the two SDs in +the variance decomposition. Plotted as KDE curves on a shared x-axis, the +ratio of their widths is the variance ratio reported in S2. + +Re-run: + python -m v4.figures.S3_variance_distributions +""" +from __future__ import annotations + +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 scipy.stats import gaussian_kde + +from v4.figures.util.loaders import RESULTS_ROOT + + +OUT = Path(__file__).parent / "output" / "S3_variance_distributions.png" + +BRIDGES = [ + ("Concat", "phase3_v4/single_bcd_concat"), + ("Pairwise", "phase3_v4/single_bcd_pairwise"), + ("Gated", "phase3_v4/single_bcd_gated"), + ("Hadamard", "refuge_v2m_baseline/ensemble_single_refugelike"), +] +STAGE_KEY = "nt_test_auc" + +C_ARCH = "#222" # dark grey for the architectural curve +BRIDGE_COLORS = { + "Concat": "#7f7f7f", # grey (underperformer) + "Pairwise": "#ff7f0e", # orange + "Gated": "#2ca02c", # green + "Hadamard": "#1f77b4", # blue (default) +} + + +def collect() -> pd.DataFrame: + rows = [] + for label, rel in BRIDGES: + root = RESULTS_ROOT / rel + for s in sorted(root.glob("rep*/binary/summary.json")): + rep = int(s.parents[1].name.replace("rep", "")) + d = json.loads(s.read_text()) + for fr in d.get("fold_results", []): + v = fr.get(STAGE_KEY) + if v is None or not np.isfinite(v): + continue + rows.append({ + "bridge": label, + "rep": rep, + "fold": int(fr["fold"]), + "auc": float(v), + }) + return pd.DataFrame(rows) + + +def render() -> None: + df = collect() + if df.empty: + print("No data collected.") + return + + # Compute centering offsets per (fold-rep, bridge) cell + fold_rep_mean = df.groupby(["rep", "fold"])["auc"].transform("mean") + bridge_mean = df.groupby("bridge")["auc"].transform("mean") + df["arch_offset"] = df["auc"] - fold_rep_mean + df["fold_offset"] = df["auc"] - bridge_mean + + # Use the same SD formula as S2 (within-group SD averaged across groups) + # so the two figures report identical pooled numbers. + cell_var = df.groupby(["rep", "fold"])["auc"].var(ddof=1) + bridge_var = df.groupby("bridge")["auc"].var(ddof=1) + arch_sd_pooled = float(np.sqrt(cell_var.mean())) + fold_sd_pooled = float(np.sqrt(bridge_var.mean())) + ratio = fold_sd_pooled / arch_sd_pooled if arch_sd_pooled > 0 else float("inf") + + # Per-bridge fold-rep SDs (50 fold-reps per bridge) + per_bridge_sd = {b: float(np.sqrt(v)) for b, v in bridge_var.items()} + + print(f"n cells = {len(df)}") + print(f"architectural SD (pooled, within-fold-rep avg) = {arch_sd_pooled:.4f}") + print(f"fold-rep SD (pooled, within-bridge avg) = {fold_sd_pooled:.4f}") + print(f"ratio fold-SD / arch-SD = {ratio:.2f}") + print("Per-bridge fold-rep SDs:") + for b in [name for name, _ in BRIDGES]: + print(f" {b:<10s} SD = {per_bridge_sd[b]:.4f}") + + # KDE x-axis: cover the union of all offset ranges + all_offsets = np.concatenate([df["arch_offset"].values, df["fold_offset"].values]) + x_max = float(np.abs(all_offsets).max()) * 1.10 + xs = np.linspace(-x_max, x_max, 600) + + fig, ax = plt.subplots(figsize=(9.4, 5.6)) + + # Per-bridge fold-rep offset curves (4 curves, 50 values each) + bridge_order = [name for name, _ in BRIDGES] + for b in bridge_order: + offs = df.loc[df["bridge"] == b, "fold_offset"].values + kde = gaussian_kde(offs) + y = kde(xs) + sd = per_bridge_sd[b] + ax.plot(xs, y, color=BRIDGE_COLORS[b], linewidth=1.6, alpha=0.92, + zorder=3, + label=f"{b} fold-rep SD = {sd:.4f}") + + # Architectural offset curve (200 values pooled across cells) + arch_offsets = df["arch_offset"].values + kde_arch = gaussian_kde(arch_offsets) + y_arch = kde_arch(xs) + ax.fill_between(xs, y_arch, color=C_ARCH, alpha=0.18, zorder=2) + ax.plot(xs, y_arch, color=C_ARCH, linewidth=2.2, linestyle="--", zorder=4, + label=f"Architectural (pooled) SD = {arch_sd_pooled:.4f}") + + # Mean line at 0 (every distribution is centered there) + ax.axvline(0, color="#666", linewidth=0.8, linestyle=":", alpha=0.6, zorder=0) + + ax.set_xlabel("AUC offset from grouping mean", fontsize=11) + ax.set_ylabel("Density", fontsize=11) + ax.set_xlim(-x_max, x_max) + ax.grid(axis="y", alpha=0.25, linestyle="--") + ax.legend(loc="upper left", fontsize=9.5, framealpha=0.92) + + # Top-right annotation with the variance ratio + txt = f"variance ratio fold-SD / arch-SD = {ratio:.2f}" + ax.text( + 0.985, 0.975, txt, transform=ax.transAxes, + ha="right", va="top", fontsize=10.5, family="monospace", + bbox=dict(boxstyle="round,pad=0.45", facecolor="white", + edgecolor="#888", alpha=0.92), + ) + + fig.suptitle( + "Architectural vs fold-rep variance: centered-offset distributions", + fontsize=12.5, fontweight="bold", y=0.995, + ) + fig.tight_layout() + + OUT.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(OUT, dpi=180, bbox_inches="tight") + plt.close(fig) + print(f"saved {OUT}") + + +if __name__ == "__main__": + render() diff --git a/v4/figures/output/F2_backbones.png b/v4/figures/output/F2_backbones.png index 17254b2..4563477 100644 Binary files a/v4/figures/output/F2_backbones.png and b/v4/figures/output/F2_backbones.png differ diff --git a/v4/figures/output/F3_single_mode_strips.png b/v4/figures/output/F3_single_mode_strips.png index 7078ed6..581a15b 100644 Binary files a/v4/figures/output/F3_single_mode_strips.png and b/v4/figures/output/F3_single_mode_strips.png differ diff --git a/v4/figures/output/F3b_anonymous_cv_roc.png b/v4/figures/output/F3b_anonymous_cv_roc.png new file mode 100644 index 0000000..f2bcd44 Binary files /dev/null and b/v4/figures/output/F3b_anonymous_cv_roc.png differ diff --git a/v4/figures/output/F4_bilateral.png b/v4/figures/output/F4_bilateral.png index 736feb3..e2625b8 100644 Binary files a/v4/figures/output/F4_bilateral.png and b/v4/figures/output/F4_bilateral.png differ diff --git a/v4/figures/output/F6_regression.png b/v4/figures/output/F6_regression.png index cc2817a..a0ff225 100644 Binary files a/v4/figures/output/F6_regression.png and b/v4/figures/output/F6_regression.png differ diff --git a/v4/figures/output/F8_attention_combined.png b/v4/figures/output/F8_attention_combined.png new file mode 100644 index 0000000..4382bec Binary files /dev/null and b/v4/figures/output/F8_attention_combined.png differ diff --git a/v4/figures/output/F8_gradcam/disc_attention_detail.png b/v4/figures/output/F8_gradcam/disc_attention_detail.png new file mode 100644 index 0000000..e26ddf5 Binary files /dev/null and b/v4/figures/output/F8_gradcam/disc_attention_detail.png differ diff --git a/v4/figures/output/F8_gradcam/mean_cam_comparison.png b/v4/figures/output/F8_gradcam/mean_cam_comparison.png index 5c317fa..5da348e 100644 Binary files a/v4/figures/output/F8_gradcam/mean_cam_comparison.png and b/v4/figures/output/F8_gradcam/mean_cam_comparison.png differ diff --git a/v4/figures/output/F8_gradcam/mean_cam_glaucoma.png b/v4/figures/output/F8_gradcam/mean_cam_glaucoma.png index 75be8ef..a99ee7a 100644 Binary files a/v4/figures/output/F8_gradcam/mean_cam_glaucoma.png and b/v4/figures/output/F8_gradcam/mean_cam_glaucoma.png differ diff --git a/v4/figures/output/F8_gradcam/mean_cam_normal.png b/v4/figures/output/F8_gradcam/mean_cam_normal.png index 1be078e..a75e772 100644 Binary files a/v4/figures/output/F8_gradcam/mean_cam_normal.png and b/v4/figures/output/F8_gradcam/mean_cam_normal.png differ diff --git a/v4/figures/output/F8_gradcam/overlay_grid_glaucoma.png b/v4/figures/output/F8_gradcam/overlay_grid_glaucoma.png new file mode 100644 index 0000000..dadb118 Binary files /dev/null and b/v4/figures/output/F8_gradcam/overlay_grid_glaucoma.png differ diff --git a/v4/figures/output/F8_gradcam/overlay_grid_normal.png b/v4/figures/output/F8_gradcam/overlay_grid_normal.png new file mode 100644 index 0000000..f7510b4 Binary files /dev/null and b/v4/figures/output/F8_gradcam/overlay_grid_normal.png differ diff --git a/v4/figures/output/F8_quadrant_attention.png b/v4/figures/output/F8_quadrant_attention.png new file mode 100644 index 0000000..89a1a66 Binary files /dev/null and b/v4/figures/output/F8_quadrant_attention.png differ diff --git a/v4/figures/output/F8_quadrant_fractions.csv b/v4/figures/output/F8_quadrant_fractions.csv new file mode 100644 index 0000000..0ea6c14 --- /dev/null +++ b/v4/figures/output/F8_quadrant_fractions.csv @@ -0,0 +1,49 @@ +class,outcome,scope,quadrant,n_eyes,mean,sd +Glaucoma,correct,disc,ST,46,0.14268406172610618,0.08672064804334276 +Glaucoma,correct,disc,SN,46,0.1891374642326569,0.15141795670756453 +Glaucoma,correct,disc,IT,46,0.28698060971321787,0.1321960653720743 +Glaucoma,correct,disc,IN,46,0.3811978502755397,0.17407907371247724 +Glaucoma,correct,peri,ST,46,0.0801709241030986,0.08494386999472202 +Glaucoma,correct,peri,SN,46,0.15918948407117592,0.20257575474186357 +Glaucoma,correct,peri,IT,46,0.34826832191551693,0.2606073381849083 +Glaucoma,correct,peri,IN,46,0.41237125804664543,0.249789089870656 +Glaucoma,correct,full,ST,46,0.09742898919008372,0.08228642449987539 +Glaucoma,correct,full,SN,46,0.18659980165443252,0.20125297725836833 +Glaucoma,correct,full,IT,46,0.36286098988076154,0.26458782651352963 +Glaucoma,correct,full,IN,46,0.35311022712131146,0.2246580892242915 +Glaucoma,incorrect,disc,ST,34,0.252112440145696,0.15034014496616682 +Glaucoma,incorrect,disc,SN,34,0.1746002632983389,0.09578946023984732 +Glaucoma,incorrect,disc,IT,34,0.3293620129186265,0.17212955919145329 +Glaucoma,incorrect,disc,IN,34,0.24392528597147015,0.17660301429323946 +Glaucoma,incorrect,peri,ST,34,0.20385530932089996,0.17952535419921464 +Glaucoma,incorrect,peri,SN,34,0.15875330500366588,0.15449712019023593 +Glaucoma,incorrect,peri,IT,34,0.3862569142356727,0.2418765080584544 +Glaucoma,incorrect,peri,IN,34,0.25113447428742586,0.24028533916258976 +Glaucoma,incorrect,full,ST,34,0.21217464524321497,0.15540407279111246 +Glaucoma,incorrect,full,SN,34,0.16330973974559065,0.12739118311502012 +Glaucoma,incorrect,full,IT,34,0.37856896329311157,0.20473194944523931 +Glaucoma,incorrect,full,IN,34,0.24594665148198902,0.20319230952963316 +Normal,correct,disc,ST,313,0.21529517366946307,0.0928060891129116 +Normal,correct,disc,SN,313,0.1571976673739274,0.08983072836980605 +Normal,correct,disc,IT,313,0.3536251021009069,0.14392434591451472 +Normal,correct,disc,IN,313,0.27388205685844386,0.13935549411695708 +Normal,correct,peri,ST,313,0.16039412908580125,0.10491506604844639 +Normal,correct,peri,SN,313,0.11198831344232964,0.1038424179426872 +Normal,correct,peri,IT,313,0.4372451776941388,0.22439216417926178 +Normal,correct,peri,IN,313,0.29037237912573877,0.21697235549864913 +Normal,correct,full,ST,313,0.17585699926973924,0.09265567144523121 +Normal,correct,full,SN,313,0.1305376994147368,0.09449518811759064 +Normal,correct,full,IT,313,0.4112080738659084,0.19447938630719017 +Normal,correct,full,IN,313,0.28239722323230754,0.18195870133945977 +Normal,incorrect,disc,ST,27,0.12299851888544425,0.11261550366741999 +Normal,incorrect,disc,SN,27,0.2121276641099002,0.1995218721251311 +Normal,incorrect,disc,IT,27,0.26471045076221783,0.17809359884921286 +Normal,incorrect,disc,IN,27,0.40016336326290614,0.22355291760060642 +Normal,incorrect,peri,ST,27,0.07760201435962506,0.10155210979149742 +Normal,incorrect,peri,SN,27,0.1965738090304749,0.2495116802723263 +Normal,incorrect,peri,IT,27,0.2814199416369157,0.2308903646345705 +Normal,incorrect,peri,IN,27,0.4444042475568166,0.27092685369261976 +Normal,incorrect,full,ST,27,0.09632595685606479,0.13797078057713352 +Normal,incorrect,full,SN,27,0.20212570434383328,0.2336043948759492 +Normal,incorrect,full,IT,27,0.3081102746977408,0.2523555953549531 +Normal,incorrect,full,IN,27,0.3934380562908973,0.24633953116106946 diff --git a/v4/figures/output/S1_geometry.png b/v4/figures/output/S1_geometry.png index 519d076..6cf0511 100644 Binary files a/v4/figures/output/S1_geometry.png and b/v4/figures/output/S1_geometry.png differ diff --git a/v4/figures/output/S2_variance_decomposition.png b/v4/figures/output/S2_variance_decomposition.png new file mode 100644 index 0000000..8ce7923 Binary files /dev/null and b/v4/figures/output/S2_variance_decomposition.png differ diff --git a/v4/figures/output/S3_variance_distributions.png b/v4/figures/output/S3_variance_distributions.png new file mode 100644 index 0000000..d569fb2 Binary files /dev/null and b/v4/figures/output/S3_variance_distributions.png differ diff --git a/v4/figures/output/S8a_comparison_panel.png b/v4/figures/output/S8a_comparison_panel.png index 8478d83..d52ce61 100644 Binary files a/v4/figures/output/S8a_comparison_panel.png and b/v4/figures/output/S8a_comparison_panel.png differ diff --git a/v4/figures/output/S8e_clinical_importance.csv b/v4/figures/output/S8e_clinical_importance.csv new file mode 100644 index 0000000..26c069e --- /dev/null +++ b/v4/figures/output/S8e_clinical_importance.csv @@ -0,0 +1,10 @@ +feature,mean_drop,std_drop,baseline_auc_mean +Age,0.1487745098039216,0.07246596166317897,0.7117647058823529 +IOP_corr,0.0645588235294118,0.053595964337907635,0.7117647058823529 +Phakic/Pseudophakic,0.031004901960784353,0.05039528009700277,0.7117647058823529 +Pachymetry,0.011421568627451003,0.016547911692771797,0.7117647058823529 +Gender,0.006102941176470622,0.022596178307964714,0.7117647058823529 +eyeID,0.0,0.0,0.7117647058823529 +dioptre_2,-0.0010294117647058861,0.003622645200240852,0.7117647058823529 +astigmatism,-0.002156862745098008,0.008336072214271729,0.7117647058823529 +dioptre_1,-0.006593137254901939,0.011003619078899692,0.7117647058823529 diff --git a/v4/figures/output/S8e_clinical_importance.png b/v4/figures/output/S8e_clinical_importance.png new file mode 100644 index 0000000..e5a4691 Binary files /dev/null and b/v4/figures/output/S8e_clinical_importance.png differ diff --git a/v4/scripts/analysis/bridge_attention_ceiling_check.py b/v4/scripts/analysis/bridge_attention_ceiling_check.py new file mode 100644 index 0000000..bbbaa0a --- /dev/null +++ b/v4/scripts/analysis/bridge_attention_ceiling_check.py @@ -0,0 +1,117 @@ +"""Variance decomposition: is the high-backbone end of the bridge_attention +sweep hitting a dataset ceiling? + +For each (rep, fold) cell, we have 12 hb_test_auc measurements — one per +(bridge × backbone) condition. We decompose the variance two ways and compare +between the LOW-backbone and HIGH-backbone halves of the gradient: + + across-arch variance @ fixed (rep,fold) + = var across the conditions in this subset, for the same fold split + (small → architectures are interchangeable at this capacity) + across-fold variance @ fixed architecture + = var across the 50 fold-reps, for one condition + (small → the fold split doesn't matter much) + +If at the high-backbone end, across-fold dwarfs across-arch, the dataset's +fold-assignment noise dominates the architectural choice — i.e. all the +strong configurations are hitting the same ceiling. + +Reads summary.json files directly, no inference needed. +""" +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pandas as pd + +SWEEP_ROOT = Path("v4/results/experiments/bridge_attention") + +LOW_BACKBONES = ["mobilenet_v2", "resnet50", "efficientnet_b0"] +HIGH_BACKBONES = ["efficientnet_v2_m", "refugelike", "refuge_efficientnet_v2_m"] +BRIDGES = ["gated", "ortho"] + + +def collect_long() -> pd.DataFrame: + rows = [] + for bridge in BRIDGES: + for bb in LOW_BACKBONES + HIGH_BACKBONES: + run_dir = SWEEP_ROOT / f"{bridge}_{bb}" + if not run_dir.exists(): + continue + for s in sorted(run_dir.glob("rep*/binary/summary.json")): + rep = int(s.parents[1].name.replace("rep", "")) + d = json.loads(s.read_text()) + for fr in d.get("fold_results", []): + v = fr.get("hb_test_auc") + if v is None or not np.isfinite(v): + continue + rows.append({ + "bridge": bridge, + "backbone": bb, + "rep": rep, + "fold": fr["fold"], + "auc": float(v), + }) + return pd.DataFrame(rows) + + +def decompose(df: pd.DataFrame, label: str) -> None: + # condition = bridge × backbone tuple + df = df.copy() + df["condition"] = df["bridge"] + "/" + df["backbone"] + n_cond = df["condition"].nunique() + n_cells = df.groupby(["rep", "fold"]).ngroups + + # Across-architecture variance @ fixed (rep, fold) + grouped_cell = df.groupby(["rep", "fold"])["auc"] + cell_var = grouped_cell.var(ddof=1) # one var per (rep,fold) cell + cell_std_mean = float(np.sqrt(cell_var.mean())) if not cell_var.empty else float("nan") + + # Across-fold-rep variance @ fixed architecture + grouped_arch = df.groupby("condition")["auc"] + arch_var = grouped_arch.var(ddof=1) + arch_std_mean = float(np.sqrt(arch_var.mean())) if not arch_var.empty else float("nan") + + ratio = arch_std_mean / cell_std_mean if cell_std_mean > 0 else float("inf") + mean_auc = float(df["auc"].mean()) + print(f"── {label} ── (n_conditions={n_cond}, n_cells={n_cells}, n_obs={len(df)})") + print(f" mean AUC across all (cond, rep, fold) ........ {mean_auc:.4f}") + print(f" across-arch SD @ fixed (rep,fold) ............ {cell_std_mean:.4f} " + f"← architectural spread within the same fold") + print(f" across-foldrep SD @ fixed architecture ....... {arch_std_mean:.4f} " + f"← fold-assignment noise within a single architecture") + print(f" ratio arch-SD / fold-SD ..................... {ratio:>6.2f}x " + f"({'fold noise dominates' if ratio > 2 else 'arch + fold comparable'})") + print() + + +def main(): + df = collect_long() + if df.empty: + print("No data — has the full readout finished yet?") + return + + n_cond = df["bridge"].nunique() * df["backbone"].nunique() + print(f"Collected {len(df)} (cond, rep, fold) AUC observations from " + f"{n_cond} (bridge × backbone) conditions\n") + + low = df[df["backbone"].isin(LOW_BACKBONES)] + high = df[df["backbone"].isin(HIGH_BACKBONES)] + decompose(low, "LOW backbones (mobilenet_v2, resnet50, efficientnet_b0)") + decompose(high, "HIGH backbones (efficientnet_v2_m, refugelike, refuge_v2m)") + + # Headline interpretation + print("──────────────────────────────────────────────────────") + print("Interpretation:") + print(" - If both subsets show high arch-SD: architectures genuinely differ.") + print(" - If LOW shows high arch-SD but HIGH shows low arch-SD: ceiling effect") + print(" at the high-backbone end — all strong configurations hit the same wall.") + print(" - Within each subset, ratio = fold-SD / arch-SD: when fold noise") + print(" dominates by >2x, the cell-to-cell variation between architectures") + print(" is smaller than the noise floor introduced by patient assignment.") + + +if __name__ == "__main__": + main() diff --git a/v4/scripts/analysis/bridge_attention_readout.py b/v4/scripts/analysis/bridge_attention_readout.py new file mode 100644 index 0000000..1e4104a --- /dev/null +++ b/v4/scripts/analysis/bridge_attention_readout.py @@ -0,0 +1,300 @@ +"""Per-tower gate / contribution readout + per-head AUC sanity check. + +For each gated_ run (and ortho_), loads each fold's +checkpoints and runs inference on the test set, capturing in a single pass: + • per-sample per-stream bridge attention (sigmoid gates for + GatedAdditiveBridge; ||W_i·z_i|| pre-LN projection norms for OrthoBridge) + • per-fold AUCs of the eye-level aux heads (img_aux, cd_aux, nt_head) — + sanity check that the image tower really is getting stronger as the + backbone improves, and that the clinical tower stays consistent + +Outputs: + v4/results/experiments/bridge_attention/_summary/gated_gates.csv + v4/results/experiments/bridge_attention/_summary/ortho_stream_norms.csv + Both CSVs include img_aux_auc, cd_aux_auc, nt_head_auc columns. + +Run: + python -m v4.scripts.analysis.bridge_attention_readout + python -m v4.scripts.analysis.bridge_attention_readout --bridges gated + python -m v4.scripts.analysis.bridge_attention_readout --reps 3 # subsample for speed +""" +from __future__ import annotations + +import argparse +import importlib +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import torch +from tqdm import tqdm + +from v4.classes.split_manager import SplitManager +from v4.classes.v4_hypertower import _make_loader, build_towers, load_data +import v4.classes.v4_hypertower as orch + + +SWEEP_ROOT = Path("v4/results/experiments/bridge_attention") +OUT_DIR = SWEEP_ROOT / "_summary" + +BACKBONES = [ + "mobilenet_v2", + "resnet50", + "efficientnet_b0", + "efficientnet_v2_m", + "refugelike", + "refuge_efficientnet_v2_m", +] + + +# ───────────────────────────────────────────────────────────────────────────── +# Shared fold-module loader (trimmed version of F8's helper) +# ───────────────────────────────────────────────────────────────────────────── + +def _load_fold_bridge(ckpt_dir: Path, cfg: dict, data, device): + """Load towers + nt bridge + eye-level aux heads from a single fold ckpt dir.""" + from v4.classes.heads.classifier import ClassificationHead + + towers = build_towers(cfg["towers"], data) + for name, tower in towers.items(): + tower.load_state_dict( + torch.load(ckpt_dir / f"tower_{name}.pt", map_location="cpu") + ) + tower.to(device).eval() + + stage_by_name = {s["name"]: s for s in cfg["stages"]} + nt_cfg = stage_by_name["nt"] + nt_mod = importlib.import_module(nt_cfg["module"]) + nt = getattr(nt_mod, nt_cfg["class"])( + [towers[n].out_dim for n in nt_cfg["inputs"]], + **nt_cfg.get("args", {}), + ).to(device) + nt.load_state_dict(torch.load(ckpt_dir / "stage_nt.pt", map_location="cpu")) + nt.eval() + + nc = cfg["num_classes"] + img_aux = ClassificationHead(towers["img"].out_dim, nc).to(device) + cd_aux = ClassificationHead(towers["cd"].out_dim, nc).to(device) + nt_head = ClassificationHead(nt.out_dim, nc).to(device) + img_aux.load_state_dict(torch.load(ckpt_dir / "stage_img_aux.pt", map_location="cpu")) + cd_aux.load_state_dict(torch.load(ckpt_dir / "stage_cd_aux.pt", map_location="cpu")) + nt_head.load_state_dict(torch.load(ckpt_dir / "stage_nt_head.pt", map_location="cpu")) + img_aux.eval(); cd_aux.eval(); nt_head.eval() + aux = {"img_aux": img_aux, "cd_aux": cd_aux, "nt_head": nt_head} + return towers, nt, nt_cfg, aux + + +def _splits_from_cfg(cfg: dict, data): + label_filter = cfg.get("label_filter", None) + df_mode = data.df.copy() + if label_filter is not None: + df_mode = df_mode[df_mode[data.label_col].isin(label_filter)].reset_index(drop=True) + identity_cols = getattr(data, "identity_cols", []) + identity_level = cfg.get("split_identity_level", 1) + group_col = ( + identity_cols[identity_level - 1] + if identity_level and identity_cols + else None + ) + return SplitManager(group_col=group_col).build_plans( + df_mode, + label_col=data.label_col, + n_splits=cfg.get("folds", 5), + seed=cfg.get("fold_seed", 100), + ), label_filter + + +# ───────────────────────────────────────────────────────────────────────────── +# Unified inference: capture bridge attention + per-head logits in one pass +# ───────────────────────────────────────────────────────────────────────────── + +def _install_attention_hooks(nt, bridge_kind: str): + """Register the right hook(s) for the bridge kind. Returns (handles, finalize) + where finalize() reads back the captured per-sample (N, n_streams) array.""" + if bridge_kind == "gated": + captured: list[torch.Tensor] = [] + def hook(_module, _input, output): + captured.append(output.detach().cpu()) + handles = [nt.gate.register_forward_hook(hook)] + def finalize(): + return ( + torch.cat(captured, dim=0).numpy() if captured else np.empty((0, 0)) + ) + return handles, finalize + + if bridge_kind == "ortho": + # Pre-LN projection norm ||W_i z_i||_2 per stream. Post-LN norm is + # √fusion_dim by construction, so we hook the linear projection itself. + per_stream: list[list[torch.Tensor]] = [[] for _ in range(len(nt.inner.W))] + handles = [] + for i, w in enumerate(nt.inner.W): + def make_hook(idx): + def hook(_module, _input, output): + per_stream[idx].append(output.detach().norm(dim=-1).cpu()) + return hook + handles.append(w.register_forward_hook(make_hook(i))) + def finalize(): + cols = [ + torch.cat(c, dim=0).numpy() if c else np.array([]) + for c in per_stream + ] + if not all(s.size for s in cols): + return np.empty((0, 0)) + return np.stack(cols, axis=1) + return handles, finalize + + raise ValueError(f"unknown bridge_kind {bridge_kind!r}") + + +def extract_one_fold( + bridge_kind: str, + towers, nt, aux: dict, + data, split_obj, label_filter, device, +): + """Single inference pass that returns: + attention : ndarray (N, n_streams) — bridge-specific attention signal + head_aucs : dict[str, float] — AUC of softmax[:,1] for each aux head + n_samples : int + """ + from sklearn.metrics import roc_auc_score + import torch.nn.functional as F + + if split_obj.test is None: + return np.empty((0, 0)), {}, 0 + + shell = data.build_shells( + split_obj.test, level="patient", label_filter=label_filter + ) + loader = _make_loader(shell, towers, batch_size=8, shuffle=False) + + handles, finalize = _install_attention_hooks(nt, bridge_kind) + img_logits, cd_logits, nt_logits, ys = [], [], [], [] + + try: + with torch.no_grad(): + for batch in loader: + y = batch["label"].detach().cpu().numpy() + for side in ("a", "b"): + z_img = towers["img"](batch["img"][side].to(device)) + z_cd = towers["cd"](batch["cd"][side].to(device)) + img_logits.append(aux["img_aux"](z_img).cpu()) + cd_logits.append(aux["cd_aux"](z_cd).cpu()) + z_fused = nt([z_img, z_cd]) # attention hooks capture here + nt_logits.append(aux["nt_head"](z_fused).cpu()) + ys.append(y) + attention = finalize() + finally: + for h in handles: + h.remove() + + if not ys: + return attention, {}, 0 + + y_all = np.concatenate(ys) + head_aucs: dict[str, float] = {} + for name, buf in [ + ("img_aux", img_logits), + ("cd_aux", cd_logits), + ("nt_head", nt_logits), + ]: + logits = torch.cat(buf, dim=0).numpy() + probs = (np.exp(logits - logits.max(axis=1, keepdims=True)) + / np.exp(logits - logits.max(axis=1, keepdims=True)) + .sum(axis=1, keepdims=True)) + try: + head_aucs[name] = float(roc_auc_score(y_all, probs[:, 1])) + except Exception: + head_aucs[name] = float("nan") + + return attention, head_aucs, int(attention.shape[0]) + + +# ───────────────────────────────────────────────────────────────────────────── +# Driver +# ───────────────────────────────────────────────────────────────────────────── + +def run_for_bridge(bridge_kind: str, max_reps: int | None, device) -> pd.DataFrame: + rows: list[dict] = [] + for backbone in BACKBONES: + run_dir = SWEEP_ROOT / f"{bridge_kind}_{backbone}" + if not run_dir.exists(): + print(f" [{bridge_kind}/{backbone}] (missing)") + continue + rep_dirs = sorted(run_dir.glob("rep*")) + if max_reps is not None: + rep_dirs = rep_dirs[:max_reps] + for rep_dir in tqdm(rep_dirs, desc=f"{bridge_kind}/{backbone}", unit="rep"): + summary_path = rep_dir / "binary" / "summary.json" + if not summary_path.exists(): + continue + cfg = json.loads(summary_path.read_text())["config"] + orch.cfg_ref = cfg + data = load_data(cfg) + splits, label_filter = _splits_from_cfg(cfg, data) + + try: + rep_idx = int(rep_dir.name.replace("rep", "")) + except ValueError: + continue + + for fold_idx in range(cfg.get("folds", 5)): + ckpt_dir = rep_dir / "binary" / "checkpoints" / f"fold{fold_idx}" + if not ckpt_dir.exists(): + continue + towers, nt, _, aux = _load_fold_bridge(ckpt_dir, cfg, data, device) + arr, head_aucs, _ = extract_one_fold( + bridge_kind, towers, nt, aux, data, + splits[fold_idx], label_filter, device, + ) + + if arr.size: + means = arr.mean(axis=0) + rows.append({ + "bridge": bridge_kind, + "backbone": backbone, + "rep": rep_idx, + "fold": fold_idx, + "n_samples": arr.shape[0], + "img_value": float(means[0]), + "cd_value": float(means[1]), + "img_share": float(means[0] / (means.sum() + 1e-8)), + "img_aux_auc": head_aucs.get("img_aux", float("nan")), + "cd_aux_auc": head_aucs.get("cd_aux", float("nan")), + "nt_head_auc": head_aucs.get("nt_head", float("nan")), + }) + + del towers, nt, aux + if torch.cuda.is_available(): + torch.cuda.empty_cache() + return pd.DataFrame(rows) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument( + "--bridges", default="gated,ortho", + help="Comma-separated bridge kinds to process (default: gated,ortho)", + ) + ap.add_argument( + "--reps", type=int, default=None, + help="Cap reps per backbone for a quick first pass (default: all)", + ) + args = ap.parse_args() + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + print(f"device: {device}") + + OUT_DIR.mkdir(parents=True, exist_ok=True) + for bridge in args.bridges.split(","): + bridge = bridge.strip() + if not bridge: + continue + df = run_for_bridge(bridge, args.reps, device) + suffix = "gates" if bridge == "gated" else "stream_norms" + out = OUT_DIR / f"{bridge}_{suffix}.csv" + df.to_csv(out, index=False) + print(f"saved {out} ({len(df)} rows)") + + +if __name__ == "__main__": + main() diff --git a/v4/scripts/experiments/backbone_replication/anonymous_cv.json b/v4/scripts/experiments/backbone_replication/anonymous_cv.json new file mode 100644 index 0000000..e9fd862 --- /dev/null +++ b/v4/scripts/experiments/backbone_replication/anonymous_cv.json @@ -0,0 +1,10 @@ +[ + { + "_note": "F2 block 2 — Anonymous CV variant of the refugelike R50 image-only single-eye baseline. split_identity_level=0 disables patient-level fold grouping (v4_hypertower.py:254-256); fold assignment becomes patient-anonymous, allowing the two eyes of one patient to fall on opposite sides of the train/test split. Mirrors v3 phase 2 imageonly_resnet50_leaky. Same backbone, same architecture, same seed/fold_seed as the F2 block-3 baseline (refugelike img-only single-eye); only the grouping rule changes. The gap baseline→anonymous quantifies the patient-anonymous CV inflation.", + "run_name": "experiments/backbone_replication/anonymous_cv_refugelike", + "reps": 10, + "overrides": { + "split_identity_level": 0 + } + } +] diff --git a/v4/scripts/experiments/backbone_replication/anonymous_cv_ensemble_single.json b/v4/scripts/experiments/backbone_replication/anonymous_cv_ensemble_single.json new file mode 100644 index 0000000..76924fc --- /dev/null +++ b/v4/scripts/experiments/backbone_replication/anonymous_cv_ensemble_single.json @@ -0,0 +1,10 @@ +[ + { + "_note": "Anonymous CV variant of the single-eye img+cd L1 fusion configuration at the refugelike R50 backbone (the 'classic hypertower' eye-level fusion mode without L2 bilateral aggregation). split_identity_level=0 disables patient-level fold grouping; fold assignment becomes patient-anonymous, allowing the two eyes of one patient to fall on opposite sides of the train/test split. Mirrors the image-only anonymous_cv_refugelike run (F2 block 2) but with the clinical tower fused in at the L1 stage, so the patient-anonymous CV inflation can be quantified directly against the single-eye img+cd headline number (ensemble_single_refugelike rep-mean test AUC 0.874 at split_identity_level=1). Same backbone, towers, bridge, training schedule, and seeds as the single-eye headline; only the fold-grouping rule changes. Bilateral L2 is not included because patient-level aggregation is mechanically incompatible with patient-anonymous fold assignment.", + "run_name": "experiments/backbone_replication/anonymous_cv_ensemble_single_refugelike", + "reps": 10, + "overrides": { + "split_identity_level": 0 + } + } +] diff --git a/v4/scripts/experiments/backbone_replication/basic_backbones.json b/v4/scripts/experiments/backbone_replication/basic_backbones.json new file mode 100644 index 0000000..9e7b30e --- /dev/null +++ b/v4/scripts/experiments/backbone_replication/basic_backbones.json @@ -0,0 +1,46 @@ +[ + { + "_note": "F2 block 1 — basic backbones, ImageNet pretrained, image-only, single-eye, patient-grouped 5-fold CV. VGG16. Replicates v3 phase 1 imageonly_vgg16 in the v4 pipeline.", + "run_name": "experiments/backbone_replication/basic_vgg16", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "vgg16", "freeze_ratio": 0.0 } } + } + }, + + { + "_note": "F2 block 1 — MobileNetV2 (ImageNet, img-only, single-eye).", + "run_name": "experiments/backbone_replication/basic_mobilenet_v2", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "mobilenet_v2", "freeze_ratio": 0.0 } } + } + }, + + { + "_note": "F2 block 1 — DenseNet121 (ImageNet, img-only, single-eye).", + "run_name": "experiments/backbone_replication/basic_densenet121", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "densenet121", "freeze_ratio": 0.0 } } + } + }, + + { + "_note": "F2 block 1 — InceptionV3 (ImageNet, img-only, single-eye). NB: Inception expects 299x299 input — backbone_transform_config auto-overrides crop_size to 299 for inception_v3.", + "run_name": "experiments/backbone_replication/basic_inception_v3", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "inception_v3", "freeze_ratio": 0.0 } } + } + }, + + { + "_note": "F2 block 1 — ResNet50 ImageNet-pretrained (the unfrozen ImageNet R50, not refugelike). Anchor for the basic-backbone comparison.", + "run_name": "experiments/backbone_replication/basic_resnet50", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "resnet50", "freeze_ratio": 0.0 } } + } + } +] diff --git a/v4/scripts/experiments/backbone_replication/disc_crop.json b/v4/scripts/experiments/backbone_replication/disc_crop.json new file mode 100644 index 0000000..c0800f2 --- /dev/null +++ b/v4/scripts/experiments/backbone_replication/disc_crop.json @@ -0,0 +1,39 @@ +[ + { + "_note": "F2 block 2 — GT disc-contour crop. R50 refugelike, image-only, single-eye, patient-grouped 5-fold CV. Each input image is cropped to a square bbox centred on the GT disc contour with margin=2.5 (matching v3 phase 2 imageonly_resnet50_gtcrop_2.5). The crop happens in original image coords before the standard 256-resize + 224-center-crop transform pipeline runs. Eyes with no contour file fall back to the un-cropped full image.", + "run_name": "experiments/backbone_replication/gtcrop_refugelike", + "reps": 10, + "tower_overrides": { + "img": { + "args": { + "backbone": "refugelike", + "freeze_ratio": 0.0, + "crop_source": "gt", + "crop_kwargs": { "margin": 2.5, "expert": 1 } + } + } + } + }, + + { + "_note": "F2 block 2 — U-Net disc-mask crop. R50 refugelike, image-only, single-eye, patient-grouped 5-fold CV. Each input image is cropped to a square bbox centred on the U-Net-predicted disc mask with margin=2.5 (matching v3 phase 2 imageonly_resnet50_unetcrop_2.5). U-Net is loaded from the base REFUGE checkpoint and NOT fine-tuned per fold (finetune_epochs=0) to keep this an apples-to-apples preprocessing-only ablation. Eyes where U-Net predicts no disc fall back to the un-cropped full image.", + "run_name": "experiments/backbone_replication/unetcrop_refugelike", + "reps": 10, + "tower_overrides": { + "img": { + "args": { + "backbone": "refugelike", + "freeze_ratio": 0.0, + "crop_source": "unet", + "crop_kwargs": { + "margin": 2.5, + "weights_path": "models/v2/refuge/segmentation/per_image/best.pt", + "unet_size": 512, + "threshold": 0.5, + "finetune_epochs": 0 + } + } + } + } + } +] diff --git a/v4/scripts/experiments/explainability/ensemble_refugelike_ckpt.json b/v4/scripts/experiments/explainability/ensemble_refugelike_ckpt.json new file mode 100644 index 0000000..05cb61a --- /dev/null +++ b/v4/scripts/experiments/explainability/ensemble_refugelike_ckpt.json @@ -0,0 +1,10 @@ +[ + { + "_note": "10-rep checkpointed run of the production bilateral img+cd ensemble at the refugelike R50 backbone (headline configuration). Per-fold tower and stage_models state_dicts saved under each rep's checkpoints/foldN/ directory. Used for F8 explainability: reconstructing per-tower predictions at the nt (eye) and hb (patient) levels so we can compare img-tower-only, cd-tower-only, eye-fusion, and patient-fusion outputs on the headline backbone.", + "run_name": "experiments/explainability/ensemble_refugelike_ckpt", + "reps": 10, + "overrides": { + "save_checkpoints": true + } + } +] diff --git a/v4/scripts/experiments/freeze_sweep/refugelike_freeze_sweep.json b/v4/scripts/experiments/freeze_sweep/refugelike_freeze_sweep.json index 5b1587c..6cb6c6f 100644 --- a/v4/scripts/experiments/freeze_sweep/refugelike_freeze_sweep.json +++ b/v4/scripts/experiments/freeze_sweep/refugelike_freeze_sweep.json @@ -1,8 +1,8 @@ [ { - "_note": "refugelike (resnet50 + REFUGE fundus pretraining), freeze stem only (1/5 blocks). Tests whether even minimal anchoring helps stability.", + "_note": "refugelike (resnet50 + REFUGE fundus pretraining), freeze stem only (1/5 blocks). Tests whether even minimal anchoring helps stability. Bumped to 10 reps for the methods ablation table.", "run_name": "experiments/freeze_sweep/refugelike_freeze20", - "reps": 3, + "reps": 10, "tower_overrides": { "img": { "args": { "backbone": "refugelike", "freeze_ratio": 0.2 } } } @@ -18,11 +18,20 @@ }, { - "_note": "Freeze stem + layer1 + layer2 (3/5 blocks). Only the deep semantic layers adapt — most aggressive practical setting before model loses capacity.", + "_note": "Freeze stem + layer1 + layer2 (3/5 blocks). Only the deep semantic layers adapt — most aggressive practical setting before model loses capacity. Bumped to 10 reps for the methods ablation table.", "run_name": "experiments/freeze_sweep/refugelike_freeze60", - "reps": 3, + "reps": 10, "tower_overrides": { "img": { "args": { "backbone": "refugelike", "freeze_ratio": 0.6 } } } + }, + + { + "_note": "Fully frozen backbone (freeze_ratio=1.0) — image features are entirely fixed at the REFUGE-pretrained state, only the L1 bridge and downstream heads adapt. Provides the extreme end of the freeze-ratio sweep for the methods ablation table.", + "run_name": "experiments/freeze_sweep/refugelike_freeze100", + "reps": 10, + "tower_overrides": { + "img": { "args": { "backbone": "refugelike", "freeze_ratio": 1.0 } } + } } ] diff --git a/v4/scripts/experiments/reg_head/baseline_reg_nt50_floaty.json b/v4/scripts/experiments/reg_head/baseline_reg_nt50_floaty.json new file mode 100644 index 0000000..9295c8b --- /dev/null +++ b/v4/scripts/experiments/reg_head/baseline_reg_nt50_floaty.json @@ -0,0 +1,33 @@ +[ + { + "_note": "Re-run of the headline VF_MD regression configuration (baseline_reg_nt50) after the PredictionStore y_true dtype fix. The previous run stored y_true as int64, silently rounding VF_MD floats to integers (introducing ~0.3 dB of rounding noise into MAE and the regression residuals). With the fix applied, regression targets are stored as float64. Matched seeds and fold_seeds to the original baseline_reg_nt50 so fold splits are identical and rep-paired comparison is meaningful. All other settings (nt epochs = 50, all four heads as RegressionHead targeting vf_md, label_filter [0,1,2] to retain Suspect, refugelike R50 backbone, bilateral L1+L2 fusion) are unchanged.", + "run_name": "experiments/reg_head/baseline_reg_nt50_floaty", + "reps": 10, + "overrides": { + "label_filter": [0, 1, 2] + }, + "stage_overrides": { + "nt": { "epochs": 50 }, + "img_aux": { + "module": "v4.classes.heads.regression", + "class": "RegressionHead", + "args": { "dropout": 0.3, "target_key": "vf_md", "loss": "mse" } + }, + "cd_aux": { + "module": "v4.classes.heads.regression", + "class": "RegressionHead", + "args": { "dropout": 0.3, "target_key": "vf_md", "loss": "mse" } + }, + "nt_head": { + "module": "v4.classes.heads.regression", + "class": "RegressionHead", + "args": { "dropout": 0.3, "target_key": "vf_md", "loss": "mse" } + }, + "hb_head": { + "module": "v4.classes.heads.regression", + "class": "RegressionHead", + "args": { "dropout": 0.3, "target_key": "vf_md", "loss": "mse" } + } + } + } +] diff --git a/v4/scripts/experiments/sensitivity/axial_length_include.json b/v4/scripts/experiments/sensitivity/axial_length_include.json new file mode 100644 index 0000000..561b59a --- /dev/null +++ b/v4/scripts/experiments/sensitivity/axial_length_include.json @@ -0,0 +1,16 @@ +[ + { + "_note": "Sensitivity: refugelike ensemble (img + cd, bilateral) — apples-to-apples replication of tri_v1/baseline_ensemble (mean hb_test_auc ≈ 0.896) but with Axial_Length INCLUDED in the clinical feature set instead of excluded. Tests whether the pilot-era finding that Axial_Length negatively contributed to fused AUC survives the v4 architecture. Same seed/fold_seed start (1234/100) as baseline_ensemble so (rep, fold) pairs are matched for paired statistics. save_checkpoints + save_predictions enabled so the run can drive a downstream permutation-importance feature ablation if needed.", + "run_name": "experiments/sensitivity/refugelike_ensemble_with_axial_length", + "reps": 10, + "overrides": { + "save_checkpoints": true, + "save_predictions": true, + "data": { + "args": { + "exclude_cols": [] + } + } + } + } +] diff --git a/v4/scripts/experiments/sensitivity/probe_v2m_480_amp.py b/v4/scripts/experiments/sensitivity/probe_v2m_480_amp.py new file mode 100644 index 0000000..03c0dd6 --- /dev/null +++ b/v4/scripts/experiments/sensitivity/probe_v2m_480_amp.py @@ -0,0 +1,107 @@ +"""Memory probe: EfficientNetV2-M at 480x480, bilateral forward+backward, AMP bf16. + +Goal: confirm bs=8 fits in 16 GB on the available GPU before committing to the +full sensitivity experiment. + +Mimics the bilateral training step (image tower run twice on OD + OS with shared +weights, plus a small downstream head + CE loss + Adam step). The clinical tower +and L1 bridge are omitted; their memory footprint is negligible against V2-M +activations. Synthetic inputs of the correct shape — no v4 dataset needed. + +Run: python v4/scripts/experiments/sensitivity/probe_v2m_480_amp.py +Expected output: GPU name, peak memory at each phase, fit/OOM verdict. +""" +from __future__ import annotations + +import sys +import torch +from torch import nn +from torchvision.models import efficientnet_v2_m + + +BATCH = 8 +RES = 480 +DTYPE = torch.bfloat16 + + +def fmt_gb(bytes_): + return f"{bytes_ / 1024**3:.2f} GB" + + +def main(): + if not torch.cuda.is_available(): + print("No CUDA/ROCm device available; probe requires a GPU.", file=sys.stderr) + sys.exit(2) + + device = torch.device("cuda") + gpu_name = torch.cuda.get_device_name(0) + gpu_total = torch.cuda.get_device_properties(0).total_memory + print(f"GPU: {gpu_name} total VRAM: {fmt_gb(gpu_total)}") + print(f"Config: V2-M, bs={BATCH}, res={RES}, bilateral 2x forward, AMP={DTYPE}") + print("-" * 70) + + torch.cuda.empty_cache() + torch.cuda.reset_peak_memory_stats() + + backbone = efficientnet_v2_m(weights=None).to(device) + feat_dim = backbone.classifier[1].in_features + backbone.classifier = nn.Identity() + head = nn.Sequential( + nn.LayerNorm(feat_dim), + nn.Linear(feat_dim, 256), + nn.GELU(), + nn.Linear(256, 2), + ).to(device) + opt = torch.optim.Adam( + list(backbone.parameters()) + list(head.parameters()), + lr=1e-4, + ) + + print(f"After model + optimizer load: allocated={fmt_gb(torch.cuda.memory_allocated())} " + f"peak={fmt_gb(torch.cuda.max_memory_allocated())}") + + x_od = torch.randn(BATCH, 3, RES, RES, device=device) + x_os = torch.randn(BATCH, 3, RES, RES, device=device) + y = torch.randint(0, 2, (BATCH,), device=device) + + print(f"After synthetic inputs: allocated={fmt_gb(torch.cuda.memory_allocated())} " + f"peak={fmt_gb(torch.cuda.max_memory_allocated())}") + + try: + opt.zero_grad(set_to_none=True) + with torch.autocast(device_type="cuda", dtype=DTYPE): + z_od = backbone(x_od) + z_os = backbone(x_os) + z = z_od + z_os + logits = head(z) + loss = nn.functional.cross_entropy(logits, y) + + print(f"After forward: allocated={fmt_gb(torch.cuda.memory_allocated())} " + f"peak={fmt_gb(torch.cuda.max_memory_allocated())}") + + loss.backward() + print(f"After backward: allocated={fmt_gb(torch.cuda.memory_allocated())} " + f"peak={fmt_gb(torch.cuda.max_memory_allocated())}") + + opt.step() + print(f"After optimizer step: allocated={fmt_gb(torch.cuda.memory_allocated())} " + f"peak={fmt_gb(torch.cuda.max_memory_allocated())}") + + torch.cuda.synchronize() + peak = torch.cuda.max_memory_allocated() + headroom = gpu_total - peak + print("-" * 70) + print(f"VERDICT: FIT | peak={fmt_gb(peak)} of {fmt_gb(gpu_total)} " + f"headroom={fmt_gb(headroom)} ({100 * headroom / gpu_total:.1f}%)") + print(f"Loss value: {loss.item():.4f}") + + except torch.cuda.OutOfMemoryError as e: + peak = torch.cuda.max_memory_allocated() + print("-" * 70) + print(f"VERDICT: OOM | peak before OOM={fmt_gb(peak)} of {fmt_gb(gpu_total)}") + print(f"OOM details: {e}") + sys.exit(1) + + +if __name__ == "__main__": + main() diff --git a/v4/scripts/experiments/sensitivity/v2m_at_480_amp.json b/v4/scripts/experiments/sensitivity/v2m_at_480_amp.json new file mode 100644 index 0000000..9a36da8 --- /dev/null +++ b/v4/scripts/experiments/sensitivity/v2m_at_480_amp.json @@ -0,0 +1,24 @@ +[ + { + "_note": "Sensitivity: refuge V2-M ensemble (img + cd, bilateral, ortho-w0.1 inner-Hadamard not applied here — matches the plain V2-M ensemble baseline at 224 from efficientnet/refuge_efficientnetv2_m which gave 0.9132 ± 0.019). This run feeds EfficientNetV2-M at its NATIVE 480x480 input resolution instead of the pipeline default 224. Activation memory roughly 4.6x; bf16 autocast keeps bs=8 fit in 16 GB (probe peak 9.07 GB on 7800 XT). Matched seed/fold_seed (1234/100) inherited from ensemble_fused.json so reps 1-10 here pair with reps 1-10 of the 224 baseline for paired statistics. 10 reps queued — kill early if wall-clock proves prohibitive. save_checkpoints + save_predictions enabled for downstream analysis if the result is promising.", + "run_name": "experiments/sensitivity/v2m_at_480_amp", + "reps": 10, + "overrides": { + "save_checkpoints": true, + "save_predictions": true, + "training": { + "amp": true, + "amp_dtype": "bfloat16" + } + }, + "tower_overrides": { + "img": { + "args": { + "backbone": "refuge_efficientnet_v2_m", + "freeze_ratio": 0.0, + "crop_size": 480 + } + } + } + } +]