diff --git a/.gitignore b/.gitignore index 1ba52c9..ad54bdb 100644 --- a/.gitignore +++ b/.gitignore @@ -10,4 +10,5 @@ cache_data/ models/refuge/ models/v2/refuge/ **/.archive/ -.archive/ \ No newline at end of file +.archive/ +scripts/deprecated/ \ No newline at end of file diff --git a/classes/__init__.py b/classes/__init__.py index 81c93d7..e69de29 100755 --- a/classes/__init__.py +++ b/classes/__init__.py @@ -1,19 +0,0 @@ -from .clinical_data import ClinicalData -from .dataset import ClinicalDataset -from .image_tower import ImageTower -from .md_tower import MDTower -from .bridge import Bridge, VoteBridge -# from .hypertower import HyperTower -from .backbones import list_names, BackboneSpec, BACKBONES -from .papila_builders import build_papila_clinical -from .SE_attention import SEBlock, SEGateLogger -from .early_stop import EarlyStopper -__all__ = [ - "ClinicalData", - "ClinicalDataset", - "ImageTower", - "MDTower", - "Bridge", - "VoteBridge", -# "HyperTower", -] diff --git a/classes/refuge_classification.py b/classes/refuge_classification.py new file mode 100755 index 0000000..b36d85c --- /dev/null +++ b/classes/refuge_classification.py @@ -0,0 +1,849 @@ +"""REFUGE glaucoma classification with rotation-based TTT.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path +from typing import Callable, Dict, Iterable, List, Optional, Sequence, Tuple +import random + +import numpy as np +from PIL import Image +import torch +from torch import nn +from torch.utils.data import DataLoader, Dataset +from torchvision import models, transforms +from torchvision.transforms import functional as TF +import torch.nn.functional as F +from sklearn.metrics import roc_auc_score +from skimage.transform import warp_polar +from tqdm import tqdm + +from classes.geometry_features import ( + FEATURE_DIM, + EPS, + compute_geometry_features, + disc_cup_from_mask_image, +) +from classes.refuge_preprocessing import RefugePreprocessing, RefugeSample +from classes.refuge_segmentation import RefugeSegmentation +from classes.unet_segmenter import UNetSegmenter + + +# --------------------------------------------------------------------------- +# Dataset utilities +# --------------------------------------------------------------------------- + + +def _default_image_transform(size: int = 256) -> transforms.Compose: + return transforms.Compose( + [ + transforms.Resize((size, size)), + transforms.ToTensor(), + transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]), + ] + ) + + +def _augment_image_transform(size: int = 256) -> transforms.Compose: + return transforms.Compose( + [ + transforms.Resize((size, size)), + transforms.RandomHorizontalFlip(), + transforms.RandomRotation(10), + transforms.ColorJitter(0.1, 0.1, 0.1, 0.05), + transforms.ToTensor(), + transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]), + ] + ) + + +def _crop_from_geometry(image: Image.Image, geometry: Dict[str, float], size: int = 256) -> Image.Image: + cx, cy = geometry["centre_x"], geometry["centre_y"] + r = geometry["crop_radius"] + left = max(0.0, cx - r) + upper = max(0.0, cy - r) + right = min(image.width, cx + r) + lower = min(image.height, cy + r) + crop = image.crop((left, upper, right, lower)) + return crop.resize((size, size), Image.BILINEAR) + + +def _geometry_from_mask(mask: np.ndarray, scale: float) -> Dict[str, float]: + mask = np.asarray(mask) > 0 + coords = np.argwhere(mask) + if coords.size == 0: + raise RuntimeError("Empty mask; cannot derive geometry") + ys, xs = coords[:, 0], coords[:, 1] + centre_x = float(xs.mean()) + centre_y = float(ys.mean()) + width = float(xs.max() - xs.min()) + height = float(ys.max() - ys.min()) + diameter = max(width, height) + radius = diameter / 2.0 + crop_radius = radius * scale + return { + "centre_x": centre_x, + "centre_y": centre_y, + "radius": radius, + "crop_radius": crop_radius, + "crop_size": crop_radius * 2.0, + } + + +def _compute_feature_vector(disc_mask: np.ndarray, cup_mask: np.ndarray) -> np.ndarray: + return compute_geometry_features(disc_mask, cup_mask) + + +def _compute_polar_image(crop: Image.Image, size: int) -> Image.Image: + arr = np.asarray(crop).astype(np.float32) / 255.0 + radius = min(arr.shape[0], arr.shape[1]) / 2.0 + polar = warp_polar( + arr, + radius=radius, + scaling="linear", + channel_axis=-1, + ) + polar = np.clip(polar, 0.0, 1.0) + polar_img = Image.fromarray((polar * 255).astype(np.uint8)) + return polar_img.resize((size, size), Image.BILINEAR) + + +def _crop_mask_from_geometry(mask: np.ndarray, geometry: Dict[str, float], size: int) -> np.ndarray: + mask_img = Image.fromarray((mask > 0).astype(np.uint8) * 255) + cx, cy = geometry["centre_x"], geometry["centre_y"] + r = geometry["crop_radius"] + left = max(0.0, cx - r) + upper = max(0.0, cy - r) + right = min(mask_img.width, cx + r) + lower = min(mask_img.height, cy + r) + crop = mask_img.crop((left, upper, right, lower)).resize((size, size), Image.NEAREST) + return (np.asarray(crop) > 0).astype(np.uint8) + + +@dataclass +class RefugeClassificationRecord: + sample: RefugeSample + geometry: Dict[str, float] + disc_mask: Optional[np.ndarray] = None + cup_mask: Optional[np.ndarray] = None + + +class RefugeClassificationDataset(Dataset): + def __init__( + self, + records: Sequence[RefugeClassificationRecord], + transform: transforms.Compose, + polar_transform: transforms.Compose, + size: int = 256, + ) -> None: + self.records = list(records) + self.transform = transform + self.polar_transform = polar_transform + self.size = size + + def __len__(self) -> int: + return len(self.records) + + def __getitem__(self, idx: int) -> Dict[str, torch.Tensor]: + rec = self.records[idx] + image = Image.open(rec.sample.image_path).convert("RGB") + crop = _crop_from_geometry(image, rec.geometry, size=self.size) + polar_image = _compute_polar_image(crop, size=self.size) + tensor = self.transform(crop) + polar_tensor = self.polar_transform(polar_image) + + features = np.zeros((FEATURE_DIM,), dtype=np.float32) + if rec.disc_mask is not None and rec.cup_mask is not None: + disc_crop = _crop_mask_from_geometry(rec.disc_mask, rec.geometry, self.size) + cup_crop = _crop_mask_from_geometry(rec.cup_mask, rec.geometry, self.size) + features = _compute_feature_vector(disc_crop, cup_crop) + + feature_tensor = torch.from_numpy(features).float() + label = rec.sample.label + if label is None: + raise ValueError(f"Sample {rec.sample.sample_id} is missing glaucoma label") + return { + "image": tensor, + "polar": polar_tensor, + "features": feature_tensor, + "label": torch.tensor(label, dtype=torch.long), + "sample_id": rec.sample.sample_id, + } + + +class RefugeTTTDataset(Dataset): + """Dataset providing unlabeled crops for test-time training.""" + + def __init__(self, records: Sequence[RefugeClassificationRecord], transform: transforms.Compose, size: int = 256) -> None: + self.records = list(records) + self.transform = transform + self.size = size + + def __len__(self) -> int: + return len(self.records) + + def __getitem__(self, idx: int) -> torch.Tensor: + rec = self.records[idx] + image = Image.open(rec.sample.image_path).convert("RGB") + crop = _crop_from_geometry(image, rec.geometry, size=self.size) + return self.transform(crop) + + +class UNetGeometryProvider: + """Callable wrapper that derives disc geometry using a trained UNetSegmenter.""" + + def __init__( + self, + segmenter: UNetSegmenter, + threshold: float = 0.5, + tta: bool = False, + ) -> None: + self.segmenter = segmenter + self.threshold = threshold + self.tta = tta + self.segmenter.model.eval() + + def __call__(self, sample: RefugeSample, scale: float) -> Tuple[Dict[str, float], np.ndarray, np.ndarray]: + image = Image.open(sample.image_path).convert("RGB") + resized = self.segmenter.preprocess_image(image) + tensor = transforms.ToTensor()(resized) + tensor = self.segmenter._normalize_tensor(tensor) + tensor = tensor.unsqueeze(0).to(self.segmenter.device) + with torch.no_grad(): + logits = self.segmenter.model(tensor) + if self.tta: + t_h = torch.flip(tensor, dims=[3]) + log_h = self.segmenter.model(t_h) + log_h = torch.flip(log_h, dims=[3]) + t_v = torch.flip(tensor, dims=[2]) + log_v = self.segmenter.model(t_v) + log_v = torch.flip(log_v, dims=[2]) + logits = (logits + log_h + log_v) / 3.0 + probs = torch.sigmoid(logits)[0].cpu().numpy() + + disc_pred = (probs[0] > self.threshold).astype(np.uint8) * 255 + cup_pred = (probs[1] > self.threshold).astype(np.uint8) * 255 + disc_img = Image.fromarray(disc_pred, mode="L").resize(image.size, Image.NEAREST) + cup_img = Image.fromarray(cup_pred, mode="L").resize(image.size, Image.NEAREST) + disc_mask = (np.array(disc_img, dtype=np.uint8) > 0).astype(np.uint8) + cup_mask = (np.array(cup_img, dtype=np.uint8) > 0).astype(np.uint8) + cup_mask = (cup_mask > 0) & (disc_mask > 0) + cup_mask = cup_mask.astype(np.uint8) + geom = _geometry_from_mask(disc_mask, scale) + return geom, disc_mask, cup_mask + + +# --------------------------------------------------------------------------- +# Classification module +# --------------------------------------------------------------------------- + + +class ArcMarginProduct(nn.Module): + """Additive angular margin (ArcFace) head.""" + + def __init__( + self, + in_features: int, + out_features: int, + s: float = 30.0, + m: float = 0.5, + easy_margin: bool = False, + ) -> None: + super().__init__() + self.in_features = in_features + self.out_features = out_features + self.s = float(s) + self.m = float(m) + self.easy_margin = easy_margin + self.weight = nn.Parameter(torch.empty(out_features, in_features)) + nn.init.xavier_uniform_(self.weight) + + self.cos_m = math.cos(m) + self.sin_m = math.sin(m) + self.th = math.cos(math.pi - m) + self.mm = math.sin(math.pi - m) * m + + def forward(self, input: torch.Tensor, label: Optional[torch.Tensor] = None) -> torch.Tensor: + cosine = F.linear(F.normalize(input), F.normalize(self.weight)) + if label is None: + return cosine * self.s + + sine = torch.sqrt(torch.clamp(1.0 - cosine.pow(2), min=0.0)) + phi = cosine * self.cos_m - sine * self.sin_m + if self.easy_margin: + phi = torch.where(cosine > 0, phi, cosine) + else: + phi = torch.where(cosine > self.th, phi, cosine - self.mm) + + one_hot = torch.zeros_like(cosine) + one_hot.scatter_(1, label.view(-1, 1), 1.0) + logits = (one_hot * phi) + ((1.0 - one_hot) * cosine) + logits *= self.s + return logits + + +class RefugeClassification: + """Train and evaluate REFUGE glaucoma classifiers with TTT support.""" + + def __init__( + self, + preprocessing: RefugePreprocessing, + segmentation: RefugeSegmentation, + backbone: Optional[nn.Module] = None, + geometry_fn: Optional[ + Callable[ + [RefugeSample, float], + Tuple[Dict[str, float], Optional[np.ndarray], Optional[np.ndarray]], + ] + ] = None, + cache_dir: Optional[Path] = None, + use_all_labeled: bool = False, + auto_val_ratio: float = 0.1, + use_margin: bool = False, + margin_s: float = 30.0, + margin_m: float = 0.5, + ) -> None: + self.preprocessing = preprocessing + self.segmentation = segmentation + if backbone is not None: + self.backbone = backbone + in_features = getattr(self.backbone, "_feature_dim", None) + if in_features is None: + if hasattr(self.backbone, "fc") and hasattr(self.backbone.fc, "in_features"): + in_features = self.backbone.fc.in_features # type: ignore[attr-defined] + self.backbone.fc = nn.Identity() # type: ignore[attr-defined] + else: + raise ValueError( + "Provided backbone must have '_feature_dim' or expose fc.in_features" + ) + else: + self.backbone = self._default_backbone() + in_features = getattr(self.backbone, "_feature_dim", None) + if in_features is None: + in_features = self.backbone.fc.in_features # type: ignore[attr-defined] + self.backbone.fc = nn.Identity() # type: ignore[attr-defined] + self.feature_dim = in_features + self.use_polar = True + self.extra_feature_dim = FEATURE_DIM + combined_dim = self.feature_dim * (1 + int(self.use_polar)) + self.extra_feature_dim + self.margin_s = float(margin_s) + self.margin_m = float(margin_m) + self.use_margin = bool(use_margin) + if self.use_margin: + self.classifier_head = ArcMarginProduct( + combined_dim, 2, s=self.margin_s, m=self.margin_m + ) + else: + self.classifier_head = nn.Linear(combined_dim, 2) + self.rotation_head = nn.Linear(self.feature_dim, 4) + + self.train_dataset: Optional[RefugeClassificationDataset] = None + self.val_dataset: Optional[RefugeClassificationDataset] = None + self.train_loader: Optional[DataLoader] = None + self.val_loader: Optional[DataLoader] = None + self.ttt_transform = _default_image_transform() + self.train_transform = _augment_image_transform() + self.eval_transform = _default_image_transform() + self.polar_transform = _default_image_transform() + self.crop_scale = 2.5 + self.crop_size = 256 + self.geometry_cache: Dict[ + str, Tuple[Dict[str, float], Optional[np.ndarray], Optional[np.ndarray]] + ] = {} + self.train_records: List[RefugeClassificationRecord] = [] + self.val_records: List[RefugeClassificationRecord] = [] + self._geometry_fn = geometry_fn + self.cache_dir = cache_dir + if self.cache_dir is not None: + self.cache_dir.mkdir(parents=True, exist_ok=True) + self.use_all_labeled = use_all_labeled + self.auto_val_ratio = auto_val_ratio + + # ------------------------------------------------------------------ + @staticmethod + def _default_backbone() -> nn.Module: + weights = models.ResNet50_Weights.IMAGENET1K_V2 + model = models.resnet50(weights=weights) + in_features = model.fc.in_features + model.fc = nn.Identity() + setattr(model, "_feature_dim", in_features) + return model + + # ------------------------------------------------------------------ + def build_datasets( + self, + crop_scale: float = 2.5, + crop_size: int = 256, + batch_size: int = 16, + num_workers: int = 4, + ) -> None: + self.crop_scale = crop_scale + self.crop_size = crop_size + self.train_transform = _augment_image_transform(crop_size) + self.eval_transform = _default_image_transform(crop_size) + self.ttt_transform = _default_image_transform(crop_size) + self.polar_transform = _default_image_transform(crop_size) + + manifest = list(self.preprocessing.build_manifest()) + train_records: List[RefugeClassificationRecord] = [] + val_records: List[RefugeClassificationRecord] = [] + + allowed_splits = {"train", "val"} + candidates = [ + sample + for sample in manifest + if sample.label is not None and sample.split in allowed_splits + ] + + print( + f"[classifier] Building datasets from {len(candidates)} labelled samples (train/val)" + ) + + skipped: List[str] = [] + for sample in tqdm( + candidates, + desc="Preparing records", + unit="sample", + leave=False, + ): + try: + geom, disc_mask, cup_mask = self._resolve_geometry(sample, crop_scale) + except RuntimeError: + skipped.append(sample.sample_id) + continue + record = RefugeClassificationRecord( + sample=sample, + geometry=geom, + disc_mask=disc_mask, + cup_mask=cup_mask, + ) + if sample.split == "train" or ( + self.use_all_labeled and sample.split == "val" + ): + train_records.append(record) + else: + val_records.append(record) + + if skipped: + print( + f"[classifier] WARNING: {len(skipped)}/{len(candidates)} samples skipped " + f"due to empty segmentation mask: {skipped}" + ) + + if (not val_records or self.use_all_labeled) and train_records and self.auto_val_ratio > 0.0: + rng = random.Random(42) + label_groups: Dict[int, List[RefugeClassificationRecord]] = {} + for rec in train_records: + label = int(rec.sample.label or 0) + label_groups.setdefault(label, []).append(rec) + + new_train: List[RefugeClassificationRecord] = [] + new_val: List[RefugeClassificationRecord] = [] + for recs in label_groups.values(): + rng.shuffle(recs) + if len(recs) <= 1: + new_train.extend(recs) + continue + val_count = max(1, int(round(len(recs) * self.auto_val_ratio))) + if val_count >= len(recs): + val_count = len(recs) - 1 + new_val.extend(recs[:val_count]) + new_train.extend(recs[val_count:]) + + if not new_val: + # Fallback: ensure at least one validation sample if possible + if len(new_train) > 1: + new_val.append(new_train.pop()) + + if new_val: + val_records = new_val + train_records = new_train + + self.train_records = train_records + self.val_records = val_records + + print( + f"[classifier] Records ready → train: {len(train_records)}, val: {len(val_records)}" + ) + + self.train_dataset = RefugeClassificationDataset( + train_records, + transform=self.train_transform, + polar_transform=self.polar_transform, + size=crop_size, + ) + self.val_dataset = RefugeClassificationDataset( + val_records, + transform=self.eval_transform, + polar_transform=self.polar_transform, + size=crop_size, + ) + + self.train_loader = DataLoader( + self.train_dataset, + batch_size=batch_size, + shuffle=True, + num_workers=num_workers, + pin_memory=True, + ) + self.val_loader = DataLoader( + self.val_dataset, + batch_size=batch_size, + shuffle=False, + num_workers=num_workers, + pin_memory=True, + ) + + print( + "[classifier] DataLoaders prepared — training batches will start shortly" + ) + + # ------------------------------------------------------------------ + def _resolve_geometry( + self, sample: RefugeSample, scale: float + ) -> Tuple[Dict[str, float], Optional[np.ndarray], Optional[np.ndarray]]: + key = self._cache_key(sample.sample_id, scale) + cached = self.geometry_cache.get(key) + if cached is not None: + return cached + + cache_path = self._cache_path(sample.sample_id, scale) + if cache_path is not None and cache_path.exists(): + data = np.load(cache_path, allow_pickle=False) + geom = { + "centre_x": float(data["centre_x"]), + "centre_y": float(data["centre_y"]), + "radius": float(data["radius"]), + "crop_radius": float(data["crop_radius"]), + "crop_size": float(data["crop_size"]), + } + disc_mask = None + cup_mask = None + if int(data["has_disc"]): + disc_mask = data["disc_mask"].astype(np.uint8) + if int(data["has_cup"]): + cup_mask = data["cup_mask"].astype(np.uint8) + self.geometry_cache[key] = (geom, disc_mask, cup_mask) + return geom, disc_mask, cup_mask + + disc_mask: Optional[np.ndarray] = None + cup_mask: Optional[np.ndarray] = None + + if sample.mask_path and sample.mask_path.exists(): + mask_img = Image.open(sample.mask_path).convert("RGB") + disc_mask, cup_mask = disc_cup_from_mask_image(mask_img) + geom = _geometry_from_mask(disc_mask, scale) + elif self._geometry_fn is not None: + geom, disc_mask, cup_mask = self._geometry_fn(sample, scale) + else: + geom = self.segmentation.infer_disc_geometry(sample, scale=scale) + try: + pred_mask = self.segmentation.predict_mask(sample).numpy() + disc_mask = pred_mask.astype(np.uint8) + except Exception: + disc_mask = None + cup_mask = None + + if cache_path is not None: + try: + np.savez_compressed( + cache_path, + centre_x=geom["centre_x"], + centre_y=geom["centre_y"], + radius=geom["radius"], + crop_radius=geom["crop_radius"], + crop_size=geom.get("crop_size", geom["crop_radius"] * 2.0), + disc_mask=disc_mask if disc_mask is not None else np.array([], dtype=np.uint8), + cup_mask=cup_mask if cup_mask is not None else np.array([], dtype=np.uint8), + has_disc=int(disc_mask is not None), + has_cup=int(cup_mask is not None), + ) + except Exception: + pass + + self.geometry_cache[key] = (geom, disc_mask, cup_mask) + return geom, disc_mask, cup_mask + + def set_geometry_fn( + self, + geometry_fn: Optional[ + Callable[ + [RefugeSample, float], + Tuple[Dict[str, float], Optional[np.ndarray], Optional[np.ndarray]], + ] + ], + ) -> None: + self._geometry_fn = geometry_fn + self.geometry_cache.clear() + + def build_records_for_samples( + self, + samples: Sequence[RefugeSample], + crop_scale: Optional[float] = None, + progress_prefix: Optional[str] = None, + ) -> List[RefugeClassificationRecord]: + scale = crop_scale if crop_scale is not None else self.crop_scale + records: List[RefugeClassificationRecord] = [] + skipped: List[str] = [] + iterator: Iterable[RefugeSample] + if progress_prefix is not None: + iterator = tqdm(samples, desc=progress_prefix, unit="sample", leave=False) + else: + iterator = samples + labeled = [s for s in samples if s.label is not None] + for sample in iterator: + if sample.label is None: + continue + try: + geom, disc_mask, cup_mask = self._resolve_geometry(sample, scale) + except RuntimeError: + skipped.append(sample.sample_id) + continue + records.append( + RefugeClassificationRecord( + sample=sample, + geometry=geom, + disc_mask=disc_mask, + cup_mask=cup_mask, + ) + ) + prefix = f"[{progress_prefix}]" if progress_prefix else "[classifier]" + if skipped: + print( + f"{prefix} WARNING: {len(skipped)}/{len(labeled)} samples skipped " + f"due to empty segmentation mask: {skipped}" + ) + else: + print(f"{prefix} All {len(labeled)} samples processed successfully.") + return records + + def clear_disk_cache(self) -> None: + """Delete all cached geometry/mask .npz files in cache_dir.""" + if self.cache_dir is None or not self.cache_dir.exists(): + return + removed = 0 + for f in self.cache_dir.glob("*.npz"): + f.unlink() + removed += 1 + self.geometry_cache.clear() + print(f"[classifier] Cleared {removed} cached geometry files from {self.cache_dir}") + + def _cache_key(self, sample_id: str, scale: float) -> str: + scale_tag = int(round(scale * 100)) + return f"{sample_id}_s{scale_tag}" + + def _cache_path(self, sample_id: str, scale: float) -> Optional[Path]: + if self.cache_dir is None: + return None + return self.cache_dir / f"{self._cache_key(sample_id, scale)}.npz" + + # ------------------------------------------------------------------ + def train( + self, + epochs: int = 30, + lr: float = 1e-4, + weight_decay: float = 1e-4, + device: Optional[str] = None, + rotation_weight: float = 0.5, + checkpoint_dir: Optional[Path] = None, + ) -> Dict[str, float]: + if self.train_loader is None or self.val_loader is None: + self.build_datasets() + + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.backbone.to(device) + self.classifier_head.to(device) + self.rotation_head.to(device) + + params = list(self.backbone.parameters()) + list(self.classifier_head.parameters()) + list(self.rotation_head.parameters()) + optimizer = torch.optim.Adam(params, lr=lr, weight_decay=weight_decay) + clf_loss = nn.CrossEntropyLoss() + rot_loss = nn.CrossEntropyLoss() + + best_auc = 0.0 + history: Dict[str, float] = {} + + epoch_iter = tqdm(range(1, epochs + 1), desc="Epochs", unit="epoch") + + print( + f"[classifier] Starting training for {epochs} epochs with batch size {self.train_loader.batch_size}" + ) + + for epoch in epoch_iter: + self.backbone.train() + self.classifier_head.train() + self.rotation_head.train() + running_loss = 0.0 + + batch_iter = tqdm( + self.train_loader, # type: ignore[arg-type] + desc=f"Train {epoch}/{epochs}", + leave=False, + unit="batch", + ) + + for batch in batch_iter: + images = batch["image"].to(device) + polars = batch["polar"].to(device) + extra_feats = batch["features"].to(device) + labels = batch["label"].to(device) + optimizer.zero_grad() + + feats_img = self.backbone(images) + feats = feats_img + if self.use_polar: + feats_polar = self.backbone(polars) + feats = torch.cat([feats, feats_polar], dim=1) + if self.extra_feature_dim > 0: + feats = torch.cat([feats, extra_feats], dim=1) + if self.use_margin: + logits = self.classifier_head(feats, labels) + else: + logits = self.classifier_head(feats) + loss_cls = clf_loss(logits, labels) + + rot_imgs, rot_labels = self._build_rotation_batch(images) + feats_rot = self.backbone(rot_imgs) + logits_rot = self.rotation_head(feats_rot) + loss_rot = rot_loss(logits_rot, rot_labels) + + loss = loss_cls + rotation_weight * loss_rot + loss.backward() + optimizer.step() + running_loss += loss.item() * images.size(0) + + train_loss = running_loss / len(self.train_loader.dataset) # type: ignore[arg-type] + metrics = self.evaluate(device=device) + history[f"epoch_{epoch}_loss"] = train_loss + history[f"epoch_{epoch}_auc"] = metrics.get("auc", float("nan")) + + auc_val = metrics.get("auc", 0.0) + epoch_iter.set_postfix(loss=f"{train_loss:.4f}", auc=f"{auc_val:.4f}") + + if auc_val > best_auc: + best_auc = metrics["auc"] + if checkpoint_dir is not None: + checkpoint_dir.mkdir(parents=True, exist_ok=True) + torch.save({ + "backbone": self.backbone.state_dict(), + "classifier": self.classifier_head.state_dict(), + "rotation": self.rotation_head.state_dict(), + }, checkpoint_dir / "refuge_classifier_best.pt") + + return {"best_auc": best_auc, **history} + + # ------------------------------------------------------------------ + def evaluate( + self, + split: str = "val", + apply_ttt: bool = False, + device: Optional[str] = None, + ) -> Dict[str, float]: + if split != "val": + raise ValueError("Only validation split supported currently") + if self.val_loader is None: + self.build_datasets() + + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.backbone.to(device) + self.classifier_head.to(device) + self.rotation_head.to(device) + + if apply_ttt: + ttt_ds = RefugeTTTDataset(self.val_records, transform=self.ttt_transform, size=self.crop_size) + ttt_loader = DataLoader(ttt_ds, batch_size=32, shuffle=False) + self.apply_ttt(ttt_loader, device=device) + + self.backbone.eval() + self.classifier_head.eval() + preds: List[float] = [] + targets: List[int] = [] + + with torch.no_grad(): + val_iter = tqdm(self.val_loader, desc="Validate", leave=False, unit="batch") + for batch in val_iter: # type: ignore[arg-type] + images = batch["image"].to(device) + labels = batch["label"].to(device) + polars = batch["polar"].to(device) + extra_feats = batch["features"].to(device) + feats_img = self.backbone(images) + feats = feats_img + if self.use_polar: + feats_polar = self.backbone(polars) + feats = torch.cat([feats, feats_polar], dim=1) + if self.extra_feature_dim > 0: + feats = torch.cat([feats, extra_feats], dim=1) + if self.use_margin: + logits = self.classifier_head(feats) + else: + logits = self.classifier_head(feats) + probs = torch.softmax(logits, dim=1)[:, 1] + preds.extend(probs.cpu().numpy().tolist()) + targets.extend(labels.cpu().numpy().tolist()) + + auc = 0.0 + try: + if len(set(targets)) > 1: + auc = float(roc_auc_score(targets, preds)) + except ValueError: + auc = 0.0 + + return {"auc": auc} + + # ------------------------------------------------------------------ + def apply_ttt(self, loader: DataLoader, device: Optional[str] = None, steps: int = 1, lr: float = 1e-5) -> None: + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.backbone.to(device) + self.rotation_head.to(device) + self.backbone.train() + self.rotation_head.train() + + optimizer = torch.optim.Adam(list(self.backbone.parameters()) + list(self.rotation_head.parameters()), lr=lr) + criterion = nn.CrossEntropyLoss() + + for _ in range(steps): + for batch in tqdm(loader, desc="TTT adapt", leave=False, unit="batch"): + if isinstance(batch, dict): + images = batch["image"].to(device) + else: + images = batch.to(device) + optimizer.zero_grad() + rot_imgs, rot_labels = self._build_rotation_batch(images) + feats = self.backbone(rot_imgs) + logits = self.rotation_head(feats) + loss = criterion(logits, rot_labels) + loss.backward() + optimizer.step() + + # ------------------------------------------------------------------ + def extract_backbone(self) -> nn.Module: + return self.backbone + + def save_checkpoint(self, output_dir: Path) -> None: + output_dir.mkdir(parents=True, exist_ok=True) + torch.save({ + "backbone": self.backbone.state_dict(), + "classifier": self.classifier_head.state_dict(), + "rotation": self.rotation_head.state_dict(), + }, output_dir / "refuge_classifier.pt") + + def load_checkpoint(self, checkpoint_path: Path) -> None: + payload = torch.load(checkpoint_path, map_location="cpu") + self.backbone.load_state_dict(payload["backbone"]) + self.classifier_head.load_state_dict(payload["classifier"]) + self.rotation_head.load_state_dict(payload["rotation"]) + + # ------------------------------------------------------------------ + def _build_rotation_batch(self, images: torch.Tensor) -> Tuple[torch.Tensor, torch.Tensor]: + rotations = [0, 90, 180, 270] + rotated = [] + labels = [] + for idx, angle in enumerate(rotations): + rot = TF.rotate(images, angle) + rotated.append(rot) + labels.append(torch.full((images.size(0),), idx, dtype=torch.long, device=images.device)) + batch = torch.cat(rotated, dim=0) + batch_labels = torch.cat(labels, dim=0) + return batch, batch_labels diff --git a/classes/refuge_preprocessing.py b/classes/refuge_preprocessing.py new file mode 100755 index 0000000..acab6d0 --- /dev/null +++ b/classes/refuge_preprocessing.py @@ -0,0 +1,306 @@ +"""Utilities for preparing REFUGE (REFUGE1/REFUGE2) datasets. + +Builds a unified manifest across all provided splits (REFUGE1 train/val/test +and REFUGE2 validation/test), exposing image paths, glaucoma labels, disc/cup +masks, and fovea coordinates so downstream segmentation/classification modules +can operate without additional bookkeeping. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Tuple + +import pandas as pd + + +@dataclass +class RefugeSample: + """Lightweight container describing a REFUGE sample.""" + + sample_id: str + dataset: str + split: str + image_path: Path + label: Optional[int] + device: Optional[str] + mask_path: Optional[Path] + fovea_coord: Optional[Tuple[float, float]] + + +class RefugePreprocessing: + """Builds manifests and provides shared helpers for REFUGE workflows. + + Responsibilities: + * scan the REFUGE directory structure and build a consistent manifest + (train/val/test, device vendor, ground-truth labels) + * expose convenience loaders for raw RGB frames, OD/OC masks, and + optional fovea landmarks + * compute geometric metadata (disc centres, diameters) so downstream + stages can crop ROIs lazily instead of storing pre-rendered tiles + """ + + def __init__(self, root_dir: Path | str) -> None: + self.root_dir = Path(root_dir) + self._manifest = None # populated by build_manifest() + + # ------------------------------------------------------------------ + # Manifest handling + # ------------------------------------------------------------------ + def build_manifest(self, refresh: bool = False) -> Iterable[RefugeSample]: + """Return an iterable of :class:`RefugeSample` records. + + Parameters + ---------- + refresh: + when True, force a rescan of the filesystem instead of reusing the + cached manifest. + + Returns + ------- + Iterable[RefugeSample] + A sequence containing one entry per sample in the REFUGE datasets. + + Notes + ----- + The actual manifest-building logic will live here: parsing the + directory structure, reading any provided CSV/Excel metadata, and + aligning masks/labels. For now, this method raises ``NotImplementedError`` + so callers are reminded to hook it up before use. + """ + + if self._manifest is not None and not refresh: + return self._manifest + + manifest: List[RefugeSample] = [] + + manifest.extend(self._collect_refuge1_train()) + manifest.extend(self._collect_refuge1_val()) + manifest.extend(self._collect_refuge1_test()) + manifest.extend(self._collect_refuge2_val()) + manifest.extend(self._collect_refuge2_test()) + + self._manifest = manifest + return self._manifest + + # ------------------------------------------------------------------ + # Accessors for downstream modules + # ------------------------------------------------------------------ + def load_image(self, sample: RefugeSample): + """Return the RGB fundus image for ``sample``. + + Implementors should handle color-space consistency (e.g., ensure RGB vs + BGR) and any global normalisation desired across devices. + """ + + raise NotImplementedError("Image loading to be implemented") + + def load_mask(self, sample: RefugeSample): + """Return the optic disc/cup mask for ``sample`` if available.""" + + raise NotImplementedError("Mask loading to be implemented") + + def disc_geometry(self, sample: RefugeSample) -> Dict[str, float]: + """Compute disc centre and diameter from the mask. + + The segmentation module will rely on this to crop 2.5–3× disc-diameter + ROIs at training time. + """ + + raise NotImplementedError("Disc geometry helper to be implemented") + + # ------------------------------------------------------------------ + # Internal helpers + # ------------------------------------------------------------------ + def _collect_refuge1_train(self) -> List[RefugeSample]: + base = self.root_dir / "Train" / "REFUGE1-train" + if not base.exists(): + return [] + + fovea_path = base / "Fovea_location.xlsx" + fovea_map = self._read_fovea_table(fovea_path, img_col="ImgName") + + samples: List[RefugeSample] = [] + image_root = base / "Training400" + mask_root = base / "Disc_Cup_Masks" + + for label_name, label_val in ("Glaucoma", 1), ("Non-Glaucoma", 0): + img_dir = image_root / label_name + mask_dir = mask_root / label_name + if not img_dir.exists(): + continue + for image_path in sorted(img_dir.glob("*.jpg")): + img_name = image_path.name + mask_path = (mask_dir / image_path.with_suffix(".bmp").name) + fovea = fovea_map.get(img_name) + sample_id = f"refuge1_train_{image_path.stem}" + samples.append( + RefugeSample( + sample_id=sample_id, + dataset="refuge1", + split="train", + image_path=image_path, + label=label_val, + device=None, + mask_path=mask_path if mask_path.exists() else None, + fovea_coord=fovea, + ) + ) + return samples + + def _collect_refuge1_val(self) -> List[RefugeSample]: + base = self.root_dir / "Train" / "REFUGE1-val" + if not base.exists(): + return [] + + fovea_path = base / "Fovea_locations.xlsx" + df = pd.read_excel(fovea_path) + samples: List[RefugeSample] = [] + image_root = base / "REFUGE-Validation400" + mask_root = base / "Disc_Cup_Masks" + + for _, row in df.iterrows(): + img_name = row["ImgName"] + image_path = image_root / img_name + mask_path = mask_root / Path(img_name).with_suffix(".bmp").name + fovea = self._extract_fovea(row, x_key="Fovea_X", y_key="Fovea_Y") + label = int(row.get("Glaucoma Label", 0)) if not pd.isna(row.get("Glaucoma Label", 0)) else None + sample_id = f"refuge1_val_{Path(img_name).stem}" + samples.append( + RefugeSample( + sample_id=sample_id, + dataset="refuge1", + split="val", + image_path=image_path, + label=label, + device=None, + mask_path=mask_path if mask_path.exists() else None, + fovea_coord=fovea, + ) + ) + return samples + + def _collect_refuge1_test(self) -> List[RefugeSample]: + base = self.root_dir / "Train" / "REFUGE1-test" + if not base.exists(): + return [] + + df = pd.read_excel(base / "Glaucoma_label_and_Fovea_location.xlsx") + image_root = base / "Test400" + mask_root = base / "Disc_Cup_Masks" + samples: List[RefugeSample] = [] + + for _, row in df.iterrows(): + img_name = row["ImgName"] + image_path = image_root / img_name + mask_path = mask_root / Path(img_name).with_suffix(".bmp").name + fovea = self._extract_fovea(row, x_key="Fovea_X", y_key="Fovea_Y") + label = int(row.get("Label(Glaucoma=1)", 0)) if not pd.isna(row.get("Label(Glaucoma=1)", 0)) else None + sample_id = f"refuge1_test_{Path(img_name).stem}" + samples.append( + RefugeSample( + sample_id=sample_id, + dataset="refuge1", + split="test", + image_path=image_path, + label=label, + device=None, + mask_path=mask_path if mask_path.exists() else None, + fovea_coord=fovea, + ) + ) + return samples + + def _collect_refuge2_val(self) -> List[RefugeSample]: + base = self.root_dir / "Validation" + if not base.exists(): + return [] + + label_df = pd.read_csv(base / "glaucoma.csv") + fovea_df = pd.read_csv(base / "fovea.csv") + fovea_map = { + row["ImageName"]: (float(row["Fovea_X"]), float(row["Fovea_Y"])) + for _, row in fovea_df.iterrows() + } + samples: List[RefugeSample] = [] + image_root = base / "Images" + mask_root = base / "Disc_Masks" + + for _, row in label_df.iterrows(): + img_name = row["FileName"] + image_path = image_root / img_name + mask_path = mask_root / Path(img_name).with_suffix(".png").name + label = row.get("Glaucoma Risk") + label = int(label) if label == label else None + sample_id = f"refuge2_val_{Path(img_name).stem}" + samples.append( + RefugeSample( + sample_id=sample_id, + dataset="refuge2", + split="val", + image_path=image_path, + label=label, + device=None, + mask_path=mask_path if mask_path.exists() else None, + fovea_coord=fovea_map.get(img_name), + ) + ) + return samples + + def _collect_refuge2_test(self) -> List[RefugeSample]: + base = self.root_dir / "Test" + if not base.exists(): + return [] + + label_df = pd.read_excel(base / "task1.xls", header=None, names=["ImgName", "Glaucoma"]) + fovea_df = pd.read_excel(base / "fovea.xlsx") + fovea_map = { + row["ImageName"]: (float(row["Fovea_X"]), float(row["Fovea_Y"])) + for _, row in fovea_df.iterrows() + } + samples: List[RefugeSample] = [] + image_root = base / "refuge2-test" + mask_root = base / "Disc_Mask" + + for _, row in label_df.iterrows(): + img_name = row["ImgName"] + image_path = image_root / img_name + mask_path = mask_root / Path(img_name).with_suffix(".png").name + label = row.get("Glaucoma") + label = int(label) if label == label else None + sample_id = f"refuge2_test_{Path(img_name).stem}" + samples.append( + RefugeSample( + sample_id=sample_id, + dataset="refuge2", + split="test", + image_path=image_path, + label=label, + device=None, + mask_path=mask_path if mask_path.exists() else None, + fovea_coord=fovea_map.get(img_name), + ) + ) + return samples + + @staticmethod + def _read_fovea_table(path: Path, img_col: str) -> Dict[str, Tuple[float, float]]: + if not path.exists(): + return {} + df = pd.read_excel(path) + mapping: Dict[str, Tuple[float, float]] = {} + for _, row in df.iterrows(): + mapping[row[img_col]] = ( + float(row.get("Fovea_X", float("nan"))), + float(row.get("Fovea_Y", float("nan"))), + ) + return mapping + + @staticmethod + def _extract_fovea(row: pd.Series, x_key: str, y_key: str) -> Optional[Tuple[float, float]]: + x_val = row.get(x_key) + y_val = row.get(y_key) + if pd.isna(x_val) or pd.isna(y_val): + return None + return float(x_val), float(y_val) diff --git a/classes/refuge_segmentation.py b/classes/refuge_segmentation.py new file mode 100755 index 0000000..7d300c9 --- /dev/null +++ b/classes/refuge_segmentation.py @@ -0,0 +1,383 @@ +"""REFUGE optic disc / cup segmentation utilities.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Sequence, Tuple + +import numpy as np +from PIL import Image +import torch +from torch import nn +from torch.utils.data import DataLoader, Dataset +from torchvision import transforms + +from classes.refuge_preprocessing import RefugePreprocessing, RefugeSample + + +# --------------------------------------------------------------------------- +# Dataset helpers +# --------------------------------------------------------------------------- + + +def _load_rgb(path: Path) -> Image.Image: + img = Image.open(path) + if img.mode != "RGB": + img = img.convert("RGB") + return img + + +def _load_mask_array(path: Path) -> np.ndarray: + mask_img = Image.open(path).convert("L") + mask = np.array(mask_img, dtype=np.float32) + # REFUGE masks encode disc/cup with different intensities; treat any + # positive value as disc for coarse localisation. + mask = np.where(mask > 0, 1.0, 0.0) + return mask + + +@dataclass +class RefugeSegmentationSample: + sample: RefugeSample + image_path: Path + mask_path: Path + + +class RefugeSegmentationDataset(Dataset): + """Simple segmentation dataset returning tensors.""" + + def __init__( + self, + samples: Sequence[RefugeSegmentationSample], + image_size: int = 512, + ) -> None: + self.samples = list(samples) + self.image_size = image_size + self.to_tensor = transforms.ToTensor() + + def __len__(self) -> int: + return len(self.samples) + + def __getitem__(self, idx: int) -> Dict[str, torch.Tensor]: + rec = self.samples[idx] + image = _load_rgb(rec.image_path) + mask_arr = _load_mask_array(rec.mask_path) + + if self.image_size is not None: + image = image.resize((self.image_size, self.image_size), Image.BILINEAR) + mask_img = Image.fromarray(mask_arr).resize( + (self.image_size, self.image_size), Image.NEAREST + ) + mask_arr = np.array(mask_img, dtype=np.float32) + + image_tensor = self.to_tensor(image) + mask_tensor = torch.from_numpy(mask_arr).unsqueeze(0) # [1,H,W] + return { + "image": image_tensor, + "mask": mask_tensor, + "sample_id": rec.sample.sample_id, + } + + +# --------------------------------------------------------------------------- +# Model definition (lightweight U-Net) +# --------------------------------------------------------------------------- + + +class DoubleConv(nn.Module): + def __init__(self, in_channels: int, out_channels: int): + super().__init__() + self.net = nn.Sequential( + nn.Conv2d(in_channels, out_channels, 3, padding=1, bias=False), + nn.BatchNorm2d(out_channels), + nn.ReLU(inplace=True), + nn.Conv2d(out_channels, out_channels, 3, padding=1, bias=False), + nn.BatchNorm2d(out_channels), + nn.ReLU(inplace=True), + ) + + def forward(self, x: torch.Tensor) -> torch.Tensor: + return self.net(x) + + +class UNet(nn.Module): + def __init__(self, in_channels: int = 3, base_channels: int = 64): + super().__init__() + self.enc1 = DoubleConv(in_channels, base_channels) + self.enc2 = DoubleConv(base_channels, base_channels * 2) + self.enc3 = DoubleConv(base_channels * 2, base_channels * 4) + self.enc4 = DoubleConv(base_channels * 4, base_channels * 8) + + self.pool = nn.MaxPool2d(2) + self.bottleneck = DoubleConv(base_channels * 8, base_channels * 16) + + self.up4 = nn.ConvTranspose2d(base_channels * 16, base_channels * 8, 2, stride=2) + self.dec4 = DoubleConv(base_channels * 16, base_channels * 8) + self.up3 = nn.ConvTranspose2d(base_channels * 8, base_channels * 4, 2, stride=2) + self.dec3 = DoubleConv(base_channels * 8, base_channels * 4) + self.up2 = nn.ConvTranspose2d(base_channels * 4, base_channels * 2, 2, stride=2) + self.dec2 = DoubleConv(base_channels * 4, base_channels * 2) + self.up1 = nn.ConvTranspose2d(base_channels * 2, base_channels, 2, stride=2) + self.dec1 = DoubleConv(base_channels * 2, base_channels) + + self.out = nn.Conv2d(base_channels, 1, 1) + + def forward(self, x: torch.Tensor) -> torch.Tensor: + e1 = self.enc1(x) + e2 = self.enc2(self.pool(e1)) + e3 = self.enc3(self.pool(e2)) + e4 = self.enc4(self.pool(e3)) + b = self.bottleneck(self.pool(e4)) + + d4 = self.up4(b) + d4 = torch.cat([d4, e4], dim=1) + d4 = self.dec4(d4) + d3 = self.up3(d4) + d3 = torch.cat([d3, e3], dim=1) + d3 = self.dec3(d3) + d2 = self.up2(d3) + d2 = torch.cat([d2, e2], dim=1) + d2 = self.dec2(d2) + d1 = self.up1(d2) + d1 = torch.cat([d1, e1], dim=1) + d1 = self.dec1(d1) + return self.out(d1) + + +# --------------------------------------------------------------------------- +# Segmentation manager +# --------------------------------------------------------------------------- + + +class RefugeSegmentation: + """Train and run coarse-to-fine OD/OC segmentation for REFUGE.""" + + def __init__( + self, + preprocessing: RefugePreprocessing, + model: Optional[nn.Module] = None, + ) -> None: + self.preprocessing = preprocessing + self.model = model or UNet() + self.train_dataset: Optional[RefugeSegmentationDataset] = None + self.val_dataset: Optional[RefugeSegmentationDataset] = None + self.train_loader: Optional[DataLoader] = None + self.val_loader: Optional[DataLoader] = None + + # ------------------------------------------------------------------ + def build_datasets( + self, + image_size: int = 512, + batch_size: int = 8, + num_workers: int = 4, + ) -> None: + manifest = self.preprocessing.build_manifest() + + train_samples: List[RefugeSegmentationSample] = [] + val_samples: List[RefugeSegmentationSample] = [] + + for sample in manifest: + if not sample.mask_path or not sample.mask_path.exists(): + continue + rec = RefugeSegmentationSample(sample=sample, image_path=sample.image_path, mask_path=sample.mask_path) + if sample.split == "train": + train_samples.append(rec) + elif sample.split in {"val", "validation"}: + val_samples.append(rec) + + if not val_samples: + # Fall back to using a subset of training data for validation + split = max(1, int(0.1 * len(train_samples))) + val_samples = train_samples[:split] + train_samples = train_samples[split:] + + self.train_dataset = RefugeSegmentationDataset(train_samples, image_size=image_size) + self.val_dataset = RefugeSegmentationDataset(val_samples, image_size=image_size) + self.train_loader = DataLoader( + self.train_dataset, + batch_size=batch_size, + shuffle=True, + num_workers=num_workers, + pin_memory=True, + ) + self.val_loader = DataLoader( + self.val_dataset, + batch_size=batch_size, + shuffle=False, + num_workers=num_workers, + pin_memory=True, + ) + + # ------------------------------------------------------------------ + def train( + self, + epochs: int = 40, + lr: float = 1e-3, + weight_decay: float = 1e-5, + device: Optional[str] = None, + checkpoint_dir: Optional[Path] = None, + ) -> Dict[str, float]: + if self.train_loader is None or self.val_loader is None: + self.build_datasets() + + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.model.to(device) + criterion = nn.BCEWithLogitsLoss() + optimizer = torch.optim.Adam(self.model.parameters(), lr=lr, weight_decay=weight_decay) + + best_dice = 0.0 + history: Dict[str, float] = {} + + for epoch in range(1, epochs + 1): + print(f"[Seg] Processing epoch {epoch}/{epochs}") + self.model.train() + running_loss = 0.0 + for batch in self.train_loader: # type: ignore[arg-type] + images = batch["image"].to(device) + masks = batch["mask"].to(device) + optimizer.zero_grad() + logits = self.model(images) + loss = criterion(logits, masks) + loss.backward() + optimizer.step() + running_loss += loss.item() * images.size(0) + + train_loss = running_loss / len(self.train_loader.dataset) # type: ignore[arg-type] + val_metrics = self.evaluate(device=device) + history[f"epoch_{epoch}_loss"] = train_loss + history[f"epoch_{epoch}_dice"] = val_metrics.get("dice", float("nan")) + + if val_metrics.get("dice", 0.0) > best_dice: + best_dice = val_metrics["dice"] + if checkpoint_dir is not None: + checkpoint_dir.mkdir(parents=True, exist_ok=True) + torch.save(self.model.state_dict(), checkpoint_dir / "refuge_segmentation_best.pt") + + return {"best_dice": best_dice, **history} + + # ------------------------------------------------------------------ + def evaluate(self, split: str = "val", device: Optional[str] = None) -> Dict[str, float]: + if split != "val": + raise ValueError("Only validation split supported currently") + if self.val_loader is None: + self.build_datasets() + + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.model.to(device) + self.model.eval() + + dices: List[float] = [] + criterion = nn.BCEWithLogitsLoss() + losses: List[float] = [] + + with torch.no_grad(): + for batch in self.val_loader: # type: ignore[arg-type] + images = batch["image"].to(device) + masks = batch["mask"].to(device) + logits = self.model(images) + loss = criterion(logits, masks) + losses.append(loss.item() * images.size(0)) + probs = torch.sigmoid(logits) + preds = (probs > 0.5).float() + dice = self._dice_coefficient(preds, masks) + dices.extend(dice) + + mean_dice = float(np.mean(dices)) if dices else 0.0 + mean_loss = float(np.sum(losses) / len(self.val_loader.dataset)) # type: ignore[arg-type] + return {"dice": mean_dice, "loss": mean_loss} + + # ------------------------------------------------------------------ + def predict_mask(self, sample: RefugeSample, device: Optional[str] = None) -> torch.Tensor: + if self.train_dataset is None: + self.build_datasets() + device = device or ("cuda" if torch.cuda.is_available() else "cpu") + self.model.to(device) + self.model.eval() + + image = _load_rgb(sample.image_path) + original_size = image.size # (width, height) + image_resized = image.resize((self.train_dataset.image_size, self.train_dataset.image_size), Image.BILINEAR) # type: ignore[union-attr] + tensor = transforms.ToTensor()(image_resized).unsqueeze(0).to(device) + + with torch.no_grad(): + logits = self.model(tensor) + mask_resized = torch.sigmoid(logits)[0, 0] + + mask_np = mask_resized.cpu().numpy() + mask_np = (mask_np > 0.5).astype(np.float32) + mask_img = Image.fromarray(mask_np) + mask_img = mask_img.resize(original_size, Image.NEAREST) + return torch.from_numpy(np.array(mask_img, dtype=np.float32)) + + def infer_disc_geometry( + self, + sample: RefugeSample, + scale: float = 2.5, + ) -> Dict[str, float]: + if sample.mask_path and sample.mask_path.exists(): + mask = _load_mask_array(sample.mask_path) + else: + mask = self.predict_mask(sample).numpy() + + coords = np.argwhere(mask > 0.5) + if coords.size == 0: + raise RuntimeError(f"Unable to locate disc for sample {sample.sample_id}") + + ys, xs = coords[:, 0], coords[:, 1] + centre_x = float(xs.mean()) + centre_y = float(ys.mean()) + width = float(xs.max() - xs.min()) + height = float(ys.max() - ys.min()) + diameter = max(width, height) + radius = diameter / 2.0 + crop_radius = radius * scale + return { + "centre_x": centre_x, + "centre_y": centre_y, + "radius": radius, + "crop_radius": crop_radius, + "crop_size": crop_radius * 2.0, + } + + def batch_crops( + self, + samples: Iterable[RefugeSample], + scale: float = 2.5, + output_dir: Optional[Path] = None, + size: int = 256, + ) -> Dict[str, Path]: + output_paths: Dict[str, Path] = {} + if output_dir is not None: + output_dir.mkdir(parents=True, exist_ok=True) + + for sample in samples: + geom = self.infer_disc_geometry(sample, scale=scale) + image = _load_rgb(sample.image_path) + cx, cy = geom["centre_x"], geom["centre_y"] + r = geom["crop_radius"] + left = max(0.0, cx - r) + upper = max(0.0, cy - r) + right = min(image.width, cx + r) + lower = min(image.height, cy + r) + crop = image.crop((left, upper, right, lower)).resize((size, size), Image.BILINEAR) + if output_dir is not None: + out_path = output_dir / f"{sample.sample_id}_crop.png" + crop.save(out_path) + output_paths[sample.sample_id] = out_path + return output_paths + + # ------------------------------------------------------------------ + @staticmethod + def _dice_coefficient(preds: torch.Tensor, targets: torch.Tensor) -> List[float]: + eps = 1e-6 + dices = [] + preds = preds.view(preds.size(0), -1) + targets = targets.view(targets.size(0), -1) + for p, t in zip(preds, targets): + intersection = float((p * t).sum().item()) + union = float(p.sum().item() + t.sum().item()) + dice = (2.0 * intersection + eps) / (union + eps) + dices.append(dice) + return dices diff --git a/classes/v2/loader_factory.py b/classes/v2/loader_factory.py index e515ed7..f27d940 100644 --- a/classes/v2/loader_factory.py +++ b/classes/v2/loader_factory.py @@ -206,6 +206,7 @@ def make_loader( *, image_transform, image_preprocessor=None, + image_cache=None, batch_size: int, shuffle: bool, num_workers: int, @@ -216,6 +217,7 @@ def make_loader( slots, image_transform=image_transform, image_preprocessor=image_preprocessor, + image_cache=image_cache, ) return DataLoader( ds, diff --git a/classes/v2/slot_dataset.py b/classes/v2/slot_dataset.py index 9335113..d083880 100644 --- a/classes/v2/slot_dataset.py +++ b/classes/v2/slot_dataset.py @@ -1,5 +1,6 @@ from __future__ import annotations +from concurrent.futures import ThreadPoolExecutor, as_completed from typing import Any, Callable, Optional from pathlib import Path @@ -45,12 +46,14 @@ class SlotDataset(Dataset): image_transform: Optional[Callable[[Image.Image], torch.Tensor]] = None, matrix_transform: Optional[Callable[[Any], torch.Tensor]] = None, image_preprocessor: Optional[Callable[..., Image.Image]] = None, + image_cache: Optional[dict[str, np.ndarray]] = None, ) -> None: self.samples = samples self.slot_descriptors = slot_descriptors self.image_transform = image_transform or transforms.ToTensor() self.matrix_transform = matrix_transform or self._default_matrix_transform self.image_preprocessor = image_preprocessor + self.image_cache = image_cache def __len__(self) -> int: return len(self.samples) @@ -74,14 +77,72 @@ class SlotDataset(Dataset): raise ValueError("Missing required image slot") return None path = Path(value) + cache_key = str(value) + + if self.image_cache is not None: + cached = self.image_cache.get(cache_key) + if cached is not None: + return self.image_transform(Image.fromarray(cached, mode="RGB")) + img = Image.open(path).convert("RGB") if self.image_preprocessor is not None: try: img = self.image_preprocessor(img, path) except TypeError: img = self.image_preprocessor(img) + + if self.image_cache is not None: + self.image_cache[cache_key] = np.asarray(img, dtype=np.uint8) + return self.image_transform(img) + def prebuild_image_cache(self, cache_workers: int = 0) -> None: + """Pre-populate image_cache for all samples in this dataset.""" + if self.image_cache is None: + return + paths = list({ + str(record[key]) + for record in self.samples + for key, desc in self.slot_descriptors.items() + if desc.kind == "image" and record.get(key) is not None + }) + to_warm = [p for p in paths if p not in self.image_cache] + if not to_warm: + return + print( + f"[image_cache] warming {len(to_warm)} images " + f"({len(paths) - len(to_warm)} already cached)", + flush=True, + ) + + def _warm_one(path_str: str) -> None: + if path_str in self.image_cache: + return + p = Path(path_str) + img = Image.open(p).convert("RGB") + if self.image_preprocessor is not None: + try: + img = self.image_preprocessor(img, p) + except TypeError: + img = self.image_preprocessor(img) + self.image_cache[path_str] = np.asarray(img, dtype=np.uint8) + + try: + from tqdm import tqdm + except ImportError: + tqdm = None + + if cache_workers <= 1: + it = tqdm(to_warm, desc="Warm image cache", unit="img") if tqdm else to_warm + for path_str in it: + _warm_one(path_str) + else: + with ThreadPoolExecutor(max_workers=cache_workers) as ex: + futures = {ex.submit(_warm_one, p): p for p in to_warm} + it = tqdm(as_completed(futures), total=len(futures), desc="Warm image cache", unit="img") if tqdm else as_completed(futures) + for fut in it: + fut.result() + def _load_matrix(self, value: Any, *, required: bool) -> Optional[torch.Tensor]: if value is None: if required: diff --git a/classes/v2/v2_hypertower.py b/classes/v2/v2_hypertower.py index a662bf0..635e565 100644 --- a/classes/v2/v2_hypertower.py +++ b/classes/v2/v2_hypertower.py @@ -222,7 +222,13 @@ class V2HyperTower: ap.add_argument("--augment", action="store_true") ap.add_argument("--balanced-sampling", action="store_true", help="Use WeightedRandomSampler during training to equalise class frequency (default: off).") - ap.add_argument("--num-workers", type=int, default=0) + ap.add_argument("--num-workers", type=int, default=4) + ap.add_argument("--in-memory-cache", action="store_true", default=True, + help="Cache preprocessed images in RAM (default: on).") + ap.add_argument("--no-in-memory-cache", action="store_false", dest="in_memory_cache", + help="Disable in-memory image cache.") + ap.add_argument("--cache-workers", type=int, default=4, + help="Threads for prebuilding in-memory image cache (default: 4).") ap.add_argument("--device", choices=["auto", "cpu", "cuda"], default="auto") ap.add_argument("--seed", type=int, default=1234) ap.add_argument("--run-name", default=None) @@ -452,6 +458,8 @@ class V2HyperTower: n_classes=num_classes, ) + image_cache: dict | None = {} if getattr(args, "in_memory_cache", False) else None + for fold in range(n_folds): seed_everything(args.seed + fold * 100) fold_dir = tm_dir / f"fold{fold}" @@ -469,6 +477,7 @@ class V2HyperTower: fold_dir=fold_dir, tower_mode=tower_mode, pred_store=pred_store, + image_cache=image_cache, ) fold_results.append(result) if artifacts.y_true_ensemble is not None: @@ -632,6 +641,7 @@ class V2HyperTower: fold_dir: Path, tower_mode: str, pred_store: "PredictionStore | None" = None, + image_cache: "dict | None" = None, ) -> tuple[FoldResult, FoldArtifacts]: args = self.args device = self.device @@ -743,7 +753,8 @@ class V2HyperTower: slots_eye = profile_eye.slot_descriptors() slots_patient = profile_patient.slot_descriptors() - loader_kw = dict(batch_size=args.batch_size, num_workers=args.num_workers) + loader_kw = dict(batch_size=args.batch_size, num_workers=args.num_workers, + image_cache=image_cache) # ---- loaders --------------------------------------------------- use_balanced = bool(getattr(args, "balanced_sampling", False)) @@ -831,6 +842,16 @@ class V2HyperTower: if pred_store is not None: pred_store.set_split(fold, [str(s["id_1"]) for s in holdout_bilat], "holdout") + # ---- prebuild in-memory image cache (fold 0 only; shared dict fills for later folds) ---- + if image_cache is not None: + cache_workers = int(getattr(args, "cache_workers", 4)) + _loaders_to_warm = [ + train_single_loader, train_bilat_loader, val_loader, holdout_loader, + ] + for _ldr in _loaders_to_warm: + if _ldr is not None: + _ldr.dataset.prebuild_image_cache(cache_workers=cache_workers) + opt_single = torch.optim.Adam(single.parameters(), lr=args.lr) if run_single else None opt_bilateral = torch.optim.Adam(bilateral.parameters(), lr=args.lr) if run_bilat else None @@ -942,6 +963,7 @@ class V2HyperTower: # ---- epoch loop ------------------------------------------------ _prev_phase_single = "inactive" # used to detect md_warmup → next phase transition for epoch in range(total_epochs): + _epoch_t0 = time.time() if not run_single: phase_single, main_epoch_single, single_active = "inactive", 0, False elif epoch < single_warmup_md: @@ -1329,6 +1351,7 @@ class V2HyperTower: print() # seal the progress bar line if args.log_every > 0 and (epoch + 1) % args.log_every == 0: + _epoch_secs = time.time() - _epoch_t0 hld_auc = target_holdout_single_auc if run_single else bi_auc_h hld_suffix = f" hld_auc={hld_auc:.4f}" if holdout_loader is not None else "" @@ -1356,7 +1379,7 @@ class V2HyperTower: if run_single: if tower_mode == "single": msg = ( - f" ep {epoch+1:>3}/{total_epochs} " + f" ep {epoch+1:>3}/{total_epochs} ({_epoch_secs:.1f}s) " f"[single:{phase_single} {single_phase_epoch}/{single_phase_total}] " f"fused(acc={cl_acc:.4f},auc={cl_auc:.4f}) " f"img(acc={cl_acc_img:.4f},auc={cl_auc_img:.4f}) " @@ -1366,7 +1389,7 @@ class V2HyperTower: ) else: msg = ( - f" ep {epoch+1:>3}/{total_epochs} " + f" ep {epoch+1:>3}/{total_epochs} ({_epoch_secs:.1f}s) " f"[single:{phase_single} {single_phase_epoch}/{single_phase_total}] " f"fused(acc={en_acc:.4f},auc={en_auc:.4f}) " f"img(acc={en_acc_img:.4f},auc={en_auc_img:.4f}) " @@ -1376,7 +1399,7 @@ class V2HyperTower: ) else: msg = ( - f" ep {epoch+1:>3}/{total_epochs} " + f" ep {epoch+1:>3}/{total_epochs} ({_epoch_secs:.1f}s) " f"[bilat:{phase_bilat} {bilat_phase_epoch}/{bilat_phase_total}] " f"fused(acc={bi_acc:.4f},auc={bi_auc:.4f}) " f"img(acc={bi_acc_img:.4f},auc={bi_auc_img:.4f}) " diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..7242fe3 --- /dev/null +++ b/environment.yml @@ -0,0 +1,216 @@ +name: fundus_imaging +channels: + - conda-forge + - defaults +dependencies: + - _libgcc_mutex=0.1=conda_forge + - _openmp_mutex=4.5=2_gnu + - alsa-lib=1.2.11=hd590300_1 + - asttokens=2.4.1=pyhd8ed1ab_0 + - attr=2.5.1=h166bdaf_1 + - blas=1.0=openblas + - bottleneck=1.3.7=py312ha883a20_0 + - brotli=1.0.9=h5eee18b_8 + - brotli-bin=1.0.9=h5eee18b_8 + - bzip2=1.0.8=hd590300_5 + - ca-certificates=2025.9.9=h06a4308_0 + - cairo=1.18.0=h3faef2a_0 + - comm=0.2.2=pyhd8ed1ab_0 + - contourpy=1.2.0=py312hdb19cb5_0 + - cycler=0.11.0=pyhd3eb1b0_0 + - dbus=1.13.18=hb2f20db_0 + - debugpy=1.8.1=py312h30efb56_0 + - decorator=5.1.1=pyhd8ed1ab_0 + - exceptiongroup=1.2.0=pyhd8ed1ab_2 + - executing=2.0.1=pyhd8ed1ab_0 + - expat=2.6.2=h6a678d5_0 + - font-ttf-dejavu-sans-mono=2.37=hd3eb1b0_0 + - font-ttf-inconsolata=2.001=hcb22688_0 + - font-ttf-source-code-pro=2.030=hd3eb1b0_0 + - font-ttf-ubuntu=0.83=h8b1ccd4_0 + - fontconfig=2.14.2=h14ed4e7_0 + - fonts-anaconda=1=h8fa9717_0 + - fonts-conda-ecosystem=1=hd3eb1b0_0 + - fonttools=4.51.0=py312h5eee18b_0 + - freetype=2.12.1=h4a9f257_0 + - gettext=0.22.5=h59595ed_2 + - gettext-tools=0.22.5=h59595ed_2 + - glib=2.80.2=hf974151_0 + - glib-tools=2.80.2=hb6ce0ca_0 + - graphite2=1.3.14=h295c915_1 + - gst-plugins-base=1.14.1=h6a678d5_1 + - gstreamer=1.14.1=h5eee18b_1 + - harfbuzz=8.5.0=hfac3d4d_0 + - icu=73.2=h59595ed_0 + - imageio=2.37.0=py312h06a4308_0 + - importlib-metadata=7.1.0=pyha770c72_0 + - importlib_metadata=7.1.0=hd8ed1ab_0 + - ipykernel=6.29.3=pyhd33586a_0 + - ipython=8.24.0=pyh707e725_0 + - ipywidgets=8.1.2=pyhd8ed1ab_1 + - jedi=0.19.1=pyhd8ed1ab_0 + - joblib=1.4.0=py312h06a4308_0 + - jpeg=9e=h5eee18b_1 + - jupyter_client=8.6.1=pyhd8ed1ab_0 + - jupyter_core=5.7.2=py312h7900ff3_0 + - jupyterlab_widgets=3.0.10=py312h06a4308_0 + - keyutils=1.6.1=h166bdaf_0 + - kiwisolver=1.4.4=py312h6a678d5_0 + - krb5=1.20.1=h81ceb04_0 + - lame=3.100=h7b6447c_0 + - lazy_loader=0.4=py312h06a4308_0 + - lcms2=2.12=h3be6417_0 + - ld_impl_linux-64=2.40=h55db66e_0 + - lerc=3.0=h295c915_0 + - libasprintf=0.22.5=h661eb56_2 + - libasprintf-devel=0.22.5=h661eb56_2 + - libbrotlicommon=1.0.9=h5eee18b_8 + - libbrotlidec=1.0.9=h5eee18b_8 + - libbrotlienc=1.0.9=h5eee18b_8 + - libcap=2.69=h0f662aa_0 + - libclang=14.0.6=default_hc6dbbc7_1 + - libclang-cpp15=15.0.7=default_h127d8a8_5 + - libclang13=14.0.6=default_he11475f_1 + - libcups=2.4.2=h2d74bed_1 + - libdeflate=1.17=h5eee18b_1 + - libedit=3.1.20191231=he28a2e2_2 + - libevent=2.1.12=hdbd6064_1 + - libexpat=2.6.2=h59595ed_0 + - libffi=3.4.2=h7f98852_5 + - libflac=1.4.3=h59595ed_0 + - libgcc-ng=13.2.0=h77fa898_7 + - libgcrypt=1.10.3=hd590300_0 + - libgettextpo=0.22.5=h59595ed_2 + - libgettextpo-devel=0.22.5=h59595ed_2 + - libgfortran=3.0.0=1 + - libgfortran-ng=11.2.0=h00389a5_1 + - libgfortran5=11.2.0=h1234567_1 + - libglib=2.80.2=hf974151_0 + - libgomp=13.2.0=h77fa898_7 + - libgpg-error=1.49=h4f305b6_0 + - libiconv=1.17=hd590300_2 + - libjpeg-turbo=2.1.4=h166bdaf_0 + - libllvm14=14.0.6=hdb19cb5_3 + - libllvm15=15.0.7=hb3ce162_4 + - libllvm18=18.1.5=hb77312f_0 + - libnsl=2.0.1=hd590300_0 + - libogg=1.3.5=h27cfd23_1 + - libopenblas=0.3.21=h043d6bf_0 + - libopus=1.3.1=h7b6447c_0 + - libpng=1.6.43=h2797004_0 + - libpq=12.17=hdbd6064_0 + - libsndfile=1.2.2=hc60ed4a_1 + - libsodium=1.0.18=h36c2ea0_1 + - libsqlite=3.45.3=h2797004_0 + - libstdcxx-ng=13.2.0=hc0a3c3a_7 + - libsystemd0=255=h3516f8a_1 + - libtiff=4.5.1=h6a678d5_0 + - libuuid=2.38.1=h0b41bf4_0 + - libvorbis=1.3.7=h7b6447c_0 + - libwebp-base=1.3.2=h5eee18b_0 + - libxcb=1.15=h7f8727e_0 + - libxcrypt=4.4.36=hd590300_1 + - libxkbcommon=1.7.0=h662e7e4_0 + - libxml2=2.12.7=hc051c1a_0 + - libzlib=1.2.13=hd590300_5 + - lz4-c=1.9.4=h6a678d5_1 + - matplotlib=3.8.4=py312h06a4308_0 + - matplotlib-base=3.8.4=py312h526ad5a_0 + - matplotlib-inline=0.1.7=pyhd8ed1ab_0 + - mpg123=1.32.6=h59595ed_0 + - mysql=5.7.20=hf484d3e_1001 + - mysql-common=8.3.0=hf1915f5_4 + - mysql-libs=8.3.0=hca2cd23_4 + - ncurses=6.5=h59595ed_0 + - nest-asyncio=1.6.0=pyhd8ed1ab_0 + - networkx=3.4.2=py312h06a4308_0 + - nspr=4.35=h6a678d5_0 + - nss=3.100=hca3bf56_0 + - numexpr=2.8.7=py312he7dcb8a_0 + - numpy=1.26.4=py312h2809609_0 + - numpy-base=1.26.4=py312he1a6c75_0 + - openblas=0.3.4=ha44fe06_0 + - openjpeg=2.4.0=h3ad879b_0 + - openssl=3.3.0=hd590300_0 + - packaging=24.0=pyhd8ed1ab_0 + - pandas=2.2.1=py312h526ad5a_0 + - parso=0.8.4=pyhd8ed1ab_0 + - pcre2=10.43=hcad00b1_0 + - pexpect=4.9.0=pyhd8ed1ab_0 + - pickleshare=0.7.5=py_1003 + - pillow=10.3.0=py312h5eee18b_0 + - pixman=0.43.2=h59595ed_0 + - platformdirs=4.2.1=pyhd8ed1ab_0 + - ply=3.11=py312h06a4308_1 + - prompt-toolkit=3.0.42=pyha770c72_0 + - psutil=5.9.8=py312h98912ed_0 + - ptyprocess=0.7.0=pyhd3deb0d_0 + - pulseaudio-client=17.0=hb77b528_0 + - pure_eval=0.2.2=pyhd8ed1ab_0 + - pybind11-abi=5=hd3eb1b0_0 + - pygments=2.18.0=pyhd8ed1ab_0 + - pyparsing=3.0.9=py312h06a4308_0 + - pyqt=5.15.10=py312h6a678d5_0 + - pyqt5-sip=12.13.0=py312h5eee18b_0 + - python=3.12.3=hab00c5b_0_cpython + - python-dateutil=2.9.0=pyhd8ed1ab_0 + - python-tzdata=2023.3=pyhd3eb1b0_0 + - python_abi=3.12=4_cp312 + - pytz=2024.1=py312h06a4308_0 + - pyzmq=26.0.3=py312h8fd38d8_0 + - qt-main=5.15.2=h53bd1ea_10 + - readline=8.2=h8228510_1 + - scikit-image=0.25.2=py312hc74f9fe_0 + - scikit-learn=1.4.2=py312h526ad5a_1 + - scipy=1.13.0=py312h2809609_0 + - setuptools=69.5.1=pyhd8ed1ab_0 + - sip=6.7.12=py312h6a678d5_0 + - six=1.16.0=pyh6c4a22f_0 + - sqlite=3.45.3=h5eee18b_0 + - stack_data=0.6.2=pyhd8ed1ab_0 + - threadpoolctl=2.2.0=pyh0d69192_0 + - tifffile=2024.12.12=py312h06a4308_0 + - tk=8.6.13=noxft_h4845f30_101 + - tornado=6.4=py312h98912ed_0 + - traitlets=5.14.3=pyhd8ed1ab_0 + - typing_extensions=4.11.0=pyha770c72_0 + - tzdata=2024a=h0c530f3_0 + - unicodedata2=15.1.0=py312h5eee18b_0 + - wcwidth=0.2.13=pyhd8ed1ab_0 + - widgetsnbextension=4.0.10=py312h06a4308_0 + - xcb-util=0.4.0=hd590300_1 + - xcb-util-image=0.4.0=h8ee46fc_1 + - xcb-util-keysyms=0.4.0=h8ee46fc_1 + - xcb-util-renderutil=0.3.9=hd590300_1 + - xcb-util-wm=0.4.1=h8ee46fc_1 + - xkeyboard-config=2.41=hd590300_0 + - xlrd=2.0.1=pyhd3eb1b0_1 + - xorg-kbproto=1.0.7=h7f98852_1002 + - xorg-libice=1.1.1=hd590300_0 + - xorg-libsm=1.2.4=h7391055_0 + - xorg-libx11=1.8.9=h8ee46fc_0 + - xorg-libxau=1.0.11=hd590300_0 + - xorg-libxext=1.3.4=h0b41bf4_2 + - xorg-libxrender=0.9.11=hd590300_0 + - xorg-renderproto=0.11.1=h7f98852_1002 + - xorg-xextproto=7.3.0=h0b41bf4_1003 + - xorg-xf86vidmodeproto=2.3.1=h7f98852_1002 + - xorg-xproto=7.0.31=h27cfd23_1007 + - xz=5.4.6=h5eee18b_1 + - zeromq=4.3.5=h6a678d5_0 + - zipp=3.17.0=pyhd8ed1ab_0 + - zlib=1.2.13=hd590300_5 + - zstd=1.5.6=ha6fb4c9_0 + - pip: + - et-xmlfile==2.0.0 + - opencv-python==4.9.0.80 + - openpyxl==3.1.5 + - pip==25.2 + - pytorch-triton-rocm==3.2.0+rocm6.4.1.git6da9e660 + - sympy==1.13.1 + - torch==2.6.0+rocm6.4.1.git1ded221d + - torchaudio==2.6.0+rocm6.4.1.gitd8831425 + - torchvision==0.21.0+rocm6.4.1.git4040d51f + - tqdm==4.66.4 + - wheel==0.45.1 +prefix: /home/rpotter/miniconda3/envs/fundus_imaging diff --git a/scripts/basic_analysis/compare_dual_eye_towers.py b/scripts/exploratory/compare_dual_eye_towers.py similarity index 100% rename from scripts/basic_analysis/compare_dual_eye_towers.py rename to scripts/exploratory/compare_dual_eye_towers.py diff --git a/scripts/basic_analysis/compare_siamese_tower.py b/scripts/exploratory/compare_siamese_tower.py similarity index 100% rename from scripts/basic_analysis/compare_siamese_tower.py rename to scripts/exploratory/compare_siamese_tower.py diff --git a/scripts/basic_analysis/compare_siamese_v2.py b/scripts/exploratory/compare_siamese_v2.py similarity index 100% rename from scripts/basic_analysis/compare_siamese_v2.py rename to scripts/exploratory/compare_siamese_v2.py diff --git a/scripts/grid_search_analytics/best_holdout_multiclass.py b/scripts/exploratory/grid_search_analytics/best_holdout_multiclass.py similarity index 100% rename from scripts/grid_search_analytics/best_holdout_multiclass.py rename to scripts/exploratory/grid_search_analytics/best_holdout_multiclass.py diff --git a/scripts/grid_search_analytics/check_crop_cache_vs_gt.py b/scripts/exploratory/grid_search_analytics/check_crop_cache_vs_gt.py similarity index 100% rename from scripts/grid_search_analytics/check_crop_cache_vs_gt.py rename to scripts/exploratory/grid_search_analytics/check_crop_cache_vs_gt.py diff --git a/scripts/grid_search_analytics/compare_suspect_geometry_vs_image.py b/scripts/exploratory/grid_search_analytics/compare_suspect_geometry_vs_image.py similarity index 100% rename from scripts/grid_search_analytics/compare_suspect_geometry_vs_image.py rename to scripts/exploratory/grid_search_analytics/compare_suspect_geometry_vs_image.py diff --git a/scripts/grid_search_analytics/derived_analysis.py b/scripts/exploratory/grid_search_analytics/derived_analysis.py similarity index 100% rename from scripts/grid_search_analytics/derived_analysis.py rename to scripts/exploratory/grid_search_analytics/derived_analysis.py diff --git a/scripts/grid_search_analytics/derived_statistics.py b/scripts/exploratory/grid_search_analytics/derived_statistics.py similarity index 100% rename from scripts/grid_search_analytics/derived_statistics.py rename to scripts/exploratory/grid_search_analytics/derived_statistics.py diff --git a/scripts/grid_search_analytics/grid_search_heatmap.py b/scripts/exploratory/grid_search_analytics/grid_search_heatmap.py similarity index 100% rename from scripts/grid_search_analytics/grid_search_heatmap.py rename to scripts/exploratory/grid_search_analytics/grid_search_heatmap.py diff --git a/scripts/grid_search_analytics/holdout_roc_for_run.py b/scripts/exploratory/grid_search_analytics/holdout_roc_for_run.py similarity index 100% rename from scripts/grid_search_analytics/holdout_roc_for_run.py rename to scripts/exploratory/grid_search_analytics/holdout_roc_for_run.py diff --git a/scripts/grid_search_analytics/param_perf_correlations.py b/scripts/exploratory/grid_search_analytics/param_perf_correlations.py similarity index 100% rename from scripts/grid_search_analytics/param_perf_correlations.py rename to scripts/exploratory/grid_search_analytics/param_perf_correlations.py diff --git a/scripts/grid_search_analytics/rerun_grid_item.py b/scripts/exploratory/grid_search_analytics/rerun_grid_item.py similarity index 100% rename from scripts/grid_search_analytics/rerun_grid_item.py rename to scripts/exploratory/grid_search_analytics/rerun_grid_item.py diff --git a/scripts/grid_search_analytics/rerun_grid_item_v2.py b/scripts/exploratory/grid_search_analytics/rerun_grid_item_v2.py similarity index 100% rename from scripts/grid_search_analytics/rerun_grid_item_v2.py rename to scripts/exploratory/grid_search_analytics/rerun_grid_item_v2.py diff --git a/scripts/grid_search_analytics/run_derived_analysis.py b/scripts/exploratory/grid_search_analytics/run_derived_analysis.py similarity index 100% rename from scripts/grid_search_analytics/run_derived_analysis.py rename to scripts/exploratory/grid_search_analytics/run_derived_analysis.py diff --git a/scripts/legacy/classic_ml_roc.py b/scripts/legacy/classic_ml_roc.py deleted file mode 100755 index b8695dd..0000000 --- a/scripts/legacy/classic_ml_roc.py +++ /dev/null @@ -1,326 +0,0 @@ -#!/usr/bin/env python3 -"""Run classical ML models and plot *combined* ROC curves (multimodel overlays). - -Keeps your original workflow for folds/tests exactly the same. -Only changes: collects predictions per test and makes: - • One ROC plot per class (OvR), overlaying all models - • One binary ROC plot (Healthy vs Glaucoma), overlaying all models -""" -import os -import re -from pathlib import Path -from typing import Iterable, List, Tuple, Dict, Optional - -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt - -from sklearn.preprocessing import StandardScaler, label_binarize -from sklearn.pipeline import Pipeline -from sklearn.linear_model import LogisticRegression -from sklearn.neighbors import KNeighborsClassifier -from sklearn.ensemble import RandomForestClassifier -from sklearn.svm import SVC -from sklearn.metrics import roc_curve, auc - -from classes import build_papila_clinical - -# --------------------------------------------------------------------------- -# Config -# --------------------------------------------------------------------------- -SPLIT_ROOT = Path("HelpCode/kfold") -TRUST_INDEX_COL = False - -# --------------------------------------------------------------------------- -# Feature matrix -# --------------------------------------------------------------------------- -def build_feature_matrix(clinical): - df = clinical.df.copy() - scalars = ["Age", "dioptre_1", "dioptre_2", "astigmatism", "Pachymetry", "Axial_Length", "IOP_corr"] - cats = ["Gender", "Phakic/Pseudophakic"] - X = pd.concat([df[scalars], pd.get_dummies(df[cats].astype("category"), drop_first=False, prefix=cats)], axis=1) - y = df[clinical.label_col].astype(int).values - return X, y, scalars, df # X keeps NaNs; we impute per-fold - -# --------------------------------------------------------------------------- -# Models with tuned hyper-parameters (unchanged) -# --------------------------------------------------------------------------- -def make_models() -> Dict[str, Pipeline]: - return { - "LogReg": Pipeline([ - ("scaler", StandardScaler()), - ("clf", LogisticRegression( - C=1, - class_weight="balanced", - max_iter=200, - solver="lbfgs", - multi_class="auto")), - ]), - "kNN": Pipeline([ - ("scaler", StandardScaler()), - ("clf", KNeighborsClassifier( - n_neighbors=11, weights="distance")), - ]), - "RF": Pipeline([ - ("clf", RandomForestClassifier(n_estimators=200, max_depth=8, - min_samples_split=4, random_state=42)), - ]), - "SVM": Pipeline([ - ("scaler", StandardScaler()), - ("clf", SVC(C=10, kernel="rbf", gamma=0.1, probability=True)), - ]), - } - -# --------------------------------------------------------------------------- -# Split helpers copied from paper_clinical_baselines_official.py (unchanged) -# --------------------------------------------------------------------------- -_FNAME_RE = re.compile(r"RET\s*(\d+)\s*([Oo][DSs])\.jpg$", re.IGNORECASE) - -def _read_sheet_any(p: Path) -> pd.DataFrame: - if p.suffix.lower() == ".xlsx": - return pd.read_excel(p) - if p.suffix.lower() == ".csv": - return pd.read_csv(p) - if p.suffix.lower() == ".txt": - lines = [ln.strip() for ln in p.read_text(encoding="utf-8", errors="ignore").splitlines() if ln.strip()] - return pd.DataFrame({"filename": lines}) - raise ValueError(f"Unsupported split file type: {p.suffix}") - -def _normcols(cols: List[str]) -> Dict[str, str]: - def norm(s: str) -> str: - return re.sub(r"[^a-z0-9]", "", s.lower()) - return {norm(c): c for c in cols} - -def _parse_fname_to_pid_eye(fname: str) -> Optional[Tuple[int, str]]: - base = os.path.basename(str(fname)) - m = _FNAME_RE.search(base.replace(" ", "")) - if not m: - return None - return int(m.group(1)), m.group(2).upper() - -def _rows_from_sheet(sheet: pd.DataFrame, df_master: pd.DataFrame) -> List[int]: - cols = _normcols(list(sheet.columns)) - if "filename" in cols: - fn_col = cols["filename"] - lookup: Dict[str, List[int]] = {} - for i, (pid, eye) in enumerate(zip(df_master["Patient ID"].astype(int), df_master["eyeID"].astype(str))): - lookup.setdefault(f"{pid}|{eye.upper()}", []).append(i) - rows: List[int] = [] - for fn in sheet[fn_col].astype(str).tolist(): - pe = _parse_fname_to_pid_eye(fn) - if pe is None: - continue - pid, eye = pe - rows.extend(lookup.get(f"{pid}|{eye}", [])) - return rows - if "patientid" in cols and "eyeid" in cols: - pid_col, eye_col = cols["patientid"], cols["eyeid"] - lookup = {} - for i, (pid, eye) in enumerate(zip(df_master["Patient ID"].astype(int), df_master["eyeID"].astype(str))): - lookup.setdefault(f"{pid}|{eye.upper()}", []).append(i) - rows = [] - for pid, eye in zip(sheet[pid_col], sheet[eye_col]): - rows.extend(lookup.get(f"{int(pid)}|{str(eye).upper()}", [])) - return rows - if TRUST_INDEX_COL and "index" in cols: - idx = sheet[cols["index"]].astype(int).tolist() - n = len(df_master) - return [i for i in idx if 0 <= i < n] - raise RuntimeError("Split sheet missing usable columns") - -def _pair_train_test_files(dir_train: Path, dir_test: Path) -> List[Tuple[Path, Path]]: - def fold_key(p: Path) -> str: - m = re.search(r"(\d+)", p.stem) - return m.group(1) if m else p.stem.lower() - trains = sorted([p for p in dir_train.iterdir() if p.is_file() and p.suffix.lower() in (".xlsx", ".csv", ".txt")], key=fold_key) - tests = sorted([p for p in dir_test.iterdir() if p.is_file() and p.suffix.lower() in (".xlsx", ".csv", ".txt")], key=fold_key) - return [(trains[i], tests[i]) for i in range(min(len(trains), len(tests)))] - -def iter_official_folds_xlsx(clinical, split_root: Path, test_name: str) -> Iterable[Tuple[pd.DataFrame, pd.DataFrame]]: - df_master = clinical.df.copy() - test_dir = split_root / test_name - dir_train = test_dir / "Train" - dir_test = test_dir / "Test" - if not dir_train.exists() or not dir_test.exists(): - raise FileNotFoundError(f"Expected: {dir_train} and {dir_test}") - for train_file, test_file in _pair_train_test_files(dir_train, dir_test): - sh_tr, sh_te = _read_sheet_any(train_file), _read_sheet_any(test_file) - tr_rows, te_rows = _rows_from_sheet(sh_tr, df_master), _rows_from_sheet(sh_te, df_master) - tr_df, te_df = df_master.iloc[tr_rows].copy(), df_master.iloc[te_rows].copy() - yield tr_df, te_df - -# --------------------------------------------------------------------------- -# Utilities (unchanged) -# --------------------------------------------------------------------------- -def _prepare_fold_X(X: pd.DataFrame, scalars: List[str], tr_idx: np.ndarray, te_idx: np.ndarray): - Xtr, Xte = X.iloc[tr_idx].copy(), X.iloc[te_idx].copy() - med = Xtr[scalars].median(numeric_only=True) - Xtr[scalars] = Xtr[scalars].fillna(med) - Xte[scalars] = Xte[scalars].fillna(med) - return Xtr.values.astype(np.float32), Xte.values.astype(np.float32) - -# --------------------------------------------------------------------------- -# NEW: combined plotting helpers (multimodel overlays) -# --------------------------------------------------------------------------- -def _plot_multiclass_overlay(y_true: np.ndarray, prob_dict: Dict[str, np.ndarray], out_dir: Path, test_tag: str): - """One figure per class (OvR), overlaying all models.""" - n_classes = next(iter(prob_dict.values())).shape[1] - class_names = [f"Class{k}" for k in range(n_classes)] - y_bin = label_binarize(y_true, classes=list(range(n_classes))) - - for k in range(n_classes): - fig, ax = plt.subplots(figsize=(6, 5)) - for model_name, proba in prob_dict.items(): - fpr, tpr, _ = roc_curve(y_bin[:, k], proba[:, k]) - auc_val = auc(fpr, tpr) - ax.plot(fpr, tpr, lw=1.8, label=f"{model_name} (AUC={auc_val:.3f})") - ax.plot([0, 1], [0, 1], "k--", lw=1) - ax.set_xlabel("False Positive Rate") - ax.set_ylabel("True Positive Rate") - ax.set_title(f"{class_names[k]} vs Rest — {test_tag}") - ax.legend(loc="lower right") - ax.grid(True, alpha=0.3, linestyle="--") - fig.tight_layout() - fig.savefig(out_dir / f"{test_tag}_{class_names[k]}.png", dpi=170) - plt.close(fig) - -def _plot_binary_overlay(y_true: np.ndarray, prob1d_dict: Dict[str, np.ndarray], out_dir: Path, test_tag: str): - """One figure (Healthy vs Glaucoma), overlaying all models. Assumes y_true ∈ {0,1}.""" - fig, ax = plt.subplots(figsize=(6, 5)) - any_curve = False - for model_name, scores in prob1d_dict.items(): - if scores.size == 0: - continue - fpr, tpr, _ = roc_curve(y_true, scores, pos_label=1) - auc_val = auc(fpr, tpr) - ax.plot(fpr, tpr, lw=1.8, label=f"{model_name} (AUC={auc_val:.3f})") - any_curve = True - ax.plot([0, 1], [0, 1], "k--", lw=1) - ax.set_xlabel("False Positive Rate") - ax.set_ylabel("True Positive Rate") - ax.set_title(f"Binary Healthy vs Glaucoma — {test_tag}") - if any_curve: - ax.legend(loc="lower right") - ax.grid(True, alpha=0.3, linestyle="--") - fig.tight_layout() - fig.savefig(out_dir / f"{test_tag}_binary.png", dpi=170) - plt.close(fig) - -# --------------------------------------------------------------------------- -# Main (same folds/tests flow; only result collation & plotting changed) -# --------------------------------------------------------------------------- -def main(): - clinical = build_papila_clinical( - image_dir="Papila/FundusImages", - clinical_dir="Papila/ClinicalData", - label_col="Diagnosis", - cat_cols=["Gender", "Phakic/Pseudophakic"], - ) - X, y, scalars, _ = build_feature_matrix(clinical) - models = make_models() - out_dir = Path("analysis_data/roc_baselines") - out_dir.mkdir(parents=True, exist_ok=True) - - tests = [ - ("Test 3", False), ("Test 4", True) - ] if (SPLIT_ROOT / "Test 3").exists() else [ - ("Test 1", False), ("Test 2", True) - ] - - for test_name, is_binary in tests: - # Collect per-model probabilities following your original per-model loop. - # For multiclass: dict[model] -> (N, C) - # For binary: dict[model] -> (N,) (probability of class 1) - prob_dict_multi: Dict[str, np.ndarray] = {} - prob_dict_bin: Dict[str, np.ndarray] = {} - y_ref_multi: Optional[np.ndarray] = None - y_ref_bin: Optional[np.ndarray] = None - - for model_name, model in models.items(): - y_all: List[np.ndarray] = [] - p_all: List[np.ndarray] = [] - - for fold_idx, (train_df, test_df) in enumerate(iter_official_folds_xlsx(clinical, SPLIT_ROOT, test_name), 1): - # Keep your exact masking/handling - dup_rows = set(train_df.index).intersection(set(test_df.index)) - shared_pids = set(train_df["Patient ID"]).intersection(set(test_df["Patient ID"])) - if test_name in ("Test 1", "Test 2") and shared_pids: - train_df = train_df[~train_df["Patient ID"].isin(shared_pids)].copy() - dup_rows = set(train_df.index).intersection(set(test_df.index)) - shared_pids = set(train_df["Patient ID"]).intersection(set(test_df["Patient ID"])) - - tr_idx, te_idx = train_df.index.values, test_df.index.values - - if is_binary: - # original binary handling: drop Suspects on both sets - mask_tr = np.isin(y[tr_idx], [0, 1]) - mask_te = np.isin(y[te_idx], [0, 1]) - if not mask_tr.any() or not mask_te.any(): - # skip empty fold (keeps behavior safe without changing fold logic) - continue - Xtr, Xte = _prepare_fold_X(X, scalars, tr_idx[mask_tr], te_idx[mask_te]) - ytr, yte = y[tr_idx][mask_tr], y[te_idx][mask_te] - else: - Xtr, Xte = _prepare_fold_X(X, scalars, tr_idx, te_idx) - ytr, yte = y[tr_idx], y[te_idx] - - # Fit and score (unchanged approach) - model.fit(Xtr, ytr) - if is_binary: - if hasattr(model[-1], "predict_proba"): - prob = model.predict_proba(Xte)[:, 1] - else: - dec = model.decision_function(Xte) - prob = 1.0 / (1.0 + np.exp(-dec)) if np.ptp(dec) > 0 else np.full_like(dec, 0.5) - y_all.append(yte) - p_all.append(prob) - else: - if hasattr(model[-1], "predict_proba"): - prob = model.predict_proba(Xte) - else: - dec = model.decision_function(Xte) - if dec.ndim == 1: - dec = np.stack([-dec, dec], axis=1) - e = np.exp(dec - dec.max(axis=1, keepdims=True)) - prob = e / e.sum(axis=1, keepdims=True) - y_all.append(yte) - p_all.append(prob) - - if not y_all: - # No valid folds for this model under this test (e.g., all-bad after mask); skip - continue - - y_cat = np.concatenate(y_all) - p_cat = np.concatenate(p_all) - - if is_binary: - # Store 1D scores per model - prob_dict_bin[model_name] = p_cat - if y_ref_bin is None: - y_ref_bin = y_cat - else: - # Align lengths defensively (should match in normal use) - n = min(len(y_ref_bin), len(y_cat)) - y_ref_bin = y_ref_bin[:n] - prob_dict_bin[model_name] = prob_dict_bin[model_name][:n] - else: - # Store (N, C) per model - prob_dict_multi[model_name] = p_cat - if y_ref_multi is None: - y_ref_multi = y_cat - else: - # Align lengths defensively (should match in normal use) - n = min(len(y_ref_multi), len(y_cat)) - y_ref_multi = y_ref_multi[:n] - prob_dict_multi[model_name] = prob_dict_multi[model_name][:n, :] - - tag = test_name.replace(" ", "") - - # Produce overlays - if prob_dict_multi and y_ref_multi is not None: - _plot_multiclass_overlay(y_ref_multi, prob_dict_multi, out_dir, tag) - if prob_dict_bin and y_ref_bin is not None: - _plot_binary_overlay(y_ref_bin, prob_dict_bin, out_dir, tag) - -if __name__ == "__main__": - main() diff --git a/scripts/legacy/eval_papila_classifier.py b/scripts/legacy/eval_papila_classifier.py deleted file mode 100755 index ba93e40..0000000 --- a/scripts/legacy/eval_papila_classifier.py +++ /dev/null @@ -1,384 +0,0 @@ -"""Evaluate REFUGE-trained classifier on Papila images using UNet crops.""" - -from __future__ import annotations - -import argparse -import csv -from pathlib import Path -from typing import Dict, List, Optional, Sequence, Set - -import numpy as np -import torch -from torch.utils.data import DataLoader -from tqdm import tqdm -from PIL import Image, ImageDraw - -import sys - -REPO_ROOT = Path(__file__).resolve().parents[1] -if str(REPO_ROOT) not in sys.path: - sys.path.insert(0, str(REPO_ROOT)) - -from classes.refuge_preprocessing import RefugePreprocessing, RefugeSample -from classes.refuge_segmentation import RefugeSegmentation -from classes.refuge_classification import ( - RefugeClassification, - RefugeClassificationDataset, - RefugeClassificationRecord, - UNetGeometryProvider, - _default_image_transform, - _geometry_from_mask, -) -from classes.backbones import BACKBONES, load_backbone_weights -from classes.unet_segmenter import UNetSegmenter -from classes.papila_builders import build_papila_clinical - - -def parse_args() -> argparse.Namespace: - parser = argparse.ArgumentParser(description="Evaluate classifier on Papila with UNet crops") - parser.add_argument("--filtered-metrics", type=Path, required=True, help="CSV of Papila samples with acceptable Dice") - parser.add_argument("--segmenter-manifest", type=Path, required=True, help="Manifest used to train the UNet segmenter") - parser.add_argument("--segmenter-weights", type=Path, required=True, help="Path to trained UNet weights (best.pt)") - parser.add_argument("--classifier-weights", type=Path, required=False, help="Path to classifier checkpoint (refuge_classifier_best.pt)") - parser.add_argument("--refuge-root", type=Path, default=Path("REFUGE")) - parser.add_argument("--image-dir", type=Path, default=Path("Papila/FundusImages")) - parser.add_argument("--clinical-dir", type=Path, default=Path("Papila/ClinicalData")) - parser.add_argument("--label-col", type=str, default="Diagnosis", help="Column name holding Papila labels") - parser.add_argument( - "--positive-labels", - nargs="*", - default=["glaucoma", "glaucoma suspect", "suspect"], - help="Values treated as glaucoma-positive when labels are non-numeric", - ) - parser.add_argument("--dice-threshold", type=float, default=0.01, help="Minimum Dice (disc or cup) to keep a sample") - parser.add_argument("--segmenter-threshold", type=float, default=0.5, help="Probability threshold for UNet geometry") - parser.add_argument("--segmenter-normalize", choices=["none", "imagenet", "per_image"], default="per_image") - parser.add_argument("--segmenter-tta", action="store_true", help="Enable TTA (H/V flips) when deriving geometry") - parser.add_argument("--crop-scale", type=float, default=2.5) - parser.add_argument("--crop-size", type=int, default=224) - parser.add_argument("--batch-size", type=int, default=32) - parser.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu") - parser.add_argument("--output", type=Path, default=None, help="Optional CSV to store per-sample probabilities") - parser.add_argument( - "--cache-dir", - type=Path, - default=Path("analysis_data/classifier_cache"), - help="Directory to reuse classifier preprocessing cache", - ) - parser.add_argument( - "--use-gt-masks", - action="store_true", - help="Use ground truth Papila contours instead of UNet predictions", - ) - parser.add_argument( - "--gt-contours-dir", - type=Path, - default=Path("Papila/ExpertsSegmentations/Contours"), - help="Directory containing Papila contour text files", - ) - parser.add_argument( - "--backbone", - type=str, - default=None, - help="Optional backbone name (e.g. inception_v3, densenet121). Requires matching classifier weights.", - ) - return parser.parse_args() - - -def load_allowed_ids(path: Path, dice_threshold: float) -> Set[str]: - allowed: Set[str] = set() - with path.open(newline="") as fp: - reader = csv.DictReader(fp) - for row in reader: - sample_id = row.get("sample_id") - if not sample_id or sample_id == "__mean__": - continue - try: - disc = float(row.get("dice_disc", "nan")) - cup = float(row.get("dice_cup", "nan")) - except ValueError: - continue - if disc < dice_threshold and cup < dice_threshold: - continue - allowed.add(sample_id) - return allowed - - -def build_papila_samples( - image_dir: Path, - clinical_dir: Path, - label_col: str, - positive_labels: Sequence[str], - allowed_ids: Set[str], -) -> List[RefugeSample]: - clinical = build_papila_clinical( - image_dir=str(image_dir), - clinical_dir=str(clinical_dir), - label_col=label_col, - cat_cols=[], - ) - positives = {lbl.lower() for lbl in positive_labels} - samples: Dict[str, RefugeSample] = {} - for _, row in clinical.df.iterrows(): - image_path = clinical.get_image_path(row) - sample_id = f"papila_{image_path.stem}" - if sample_id not in allowed_ids or sample_id in samples: - continue - value = row.get(label_col) - if value is None or (isinstance(value, float) and np.isnan(value)): - continue - label: Optional[int] - try: - label_int = int(value) - if label_int == 2: - continue - label = 1 if label_int > 0 else 0 - except (TypeError, ValueError): - label = 1 if str(value).strip().lower() in positives else 0 - samples[sample_id] = RefugeSample( - sample_id=sample_id, - dataset="papila", - split="eval", - image_path=Path(image_path), - label=label, - device=None, - mask_path=None, - fovea_coord=None, - ) - return list(samples.values()) - - -def load_contour(path: Path) -> np.ndarray: - coords = np.loadtxt(path) - if coords.ndim == 1: - coords = coords.reshape(-1, 2) - return coords - - -def contour_to_mask(coords: np.ndarray, size: Sequence[int]) -> np.ndarray: - if coords is None or coords.size == 0: - return np.zeros((size[1], size[0]), dtype=np.uint8) - img = Image.new("L", size, 0) - draw = ImageDraw.Draw(img) - points = [tuple(map(float, pt)) for pt in coords] - draw.polygon(points, outline=1, fill=1) - return np.array(img, dtype=np.uint8) - - -class PapilaGTGeometryProvider: - def __init__(self, contours_dir: Path) -> None: - self.contours_dir = contours_dir - - def _pick(self, base: str, kind: str) -> Optional[Path]: - for exp in ("exp2", "exp1"): - cand = self.contours_dir / f"{base}_{kind}_{exp}.txt" - if cand.exists(): - return cand - return None - - def __call__(self, sample: RefugeSample, scale: float): - base = Path(sample.image_path).stem - disc_path = self._pick(base, "disc") - cup_path = self._pick(base, "cup") - if disc_path is None or cup_path is None: - raise RuntimeError(f"Missing ground-truth contours for {sample.sample_id}") - - image = Image.open(sample.image_path).convert("RGB") - disc_coords = load_contour(disc_path) - cup_coords = load_contour(cup_path) - disc_mask = contour_to_mask(disc_coords, image.size) - cup_mask = contour_to_mask(cup_coords, image.size) - cup_mask = ((cup_mask > 0) & (disc_mask > 0)).astype(np.uint8) - geom = _geometry_from_mask(disc_mask, scale) - return geom, disc_mask.astype(np.uint8), cup_mask.astype(np.uint8) - - -def build_backbone(name: Optional[str]) -> Optional[torch.nn.Module]: - if not name: - return None - key = name.lower() - if key not in BACKBONES: - raise ValueError(f"Unknown backbone '{name}'. Available: {', '.join(sorted(BACKBONES.keys()))}") - spec = BACKBONES[key] - model = spec.ctor(weights=spec.weights_default) - out_dim, model = spec.strip(model) - setattr(model, "_feature_dim", out_dim) - if key == "refugelike": - load_backbone_weights(key, model) - return model - - -def evaluate_records( - clf: RefugeClassification, - records: Sequence[RefugeClassificationRecord], - device: str, - batch_size: int, -) -> Dict[str, float]: - dataset = RefugeClassificationDataset( - records, - transform=clf.eval_transform, - polar_transform=clf.polar_transform, - size=clf.crop_size, - ) - loader = DataLoader(dataset, batch_size=batch_size, shuffle=False, num_workers=0) - clf.backbone.to(device).eval() - clf.classifier_head.to(device).eval() - preds: List[float] = [] - targets: List[int] = [] - with torch.no_grad(): - for batch in tqdm(loader, desc="Papila Eval", leave=False, unit="batch"): - images = batch["image"].to(device) - polars = batch["polar"].to(device) - extra_feats = batch["features"].to(device) - labels = batch["label"].cpu().numpy().tolist() - feats_img = clf.backbone(images) - feats = feats_img - if clf.use_polar: - feats_polar = clf.backbone(polars) - feats = torch.cat([feats, feats_polar], dim=1) - if clf.extra_feature_dim > 0: - feats = torch.cat([feats, extra_feats], dim=1) - logits = clf.classifier_head(feats) - probs = torch.softmax(logits, dim=1)[:, 1].cpu().numpy().tolist() - preds.extend(probs) - targets.extend(labels) - metrics: Dict[str, float] = {"count": float(len(targets))} - unique_labels = set(targets) - if len(unique_labels) >= 2: - metrics["auc"] = float(torchmetrics_auc(targets, preds)) - else: - metrics["auc"] = float("nan") - preds_bin = [1 if p >= 0.5 else 0 for p in preds] - accuracy = sum(int(p == t) for p, t in zip(preds_bin, targets)) / max(1, len(targets)) - metrics["accuracy"] = float(accuracy) - metrics["mean_prob"] = float(np.mean(preds)) if preds else float("nan") - metrics["labels_pos"] = float(sum(targets)) - if preds: - metrics["probs_std"] = float(np.std(preds)) - return metrics - - -def torchmetrics_auc(targets: Sequence[int], preds: Sequence[float]) -> float: - try: - from sklearn.metrics import roc_auc_score - except ImportError as exc: - raise RuntimeError("scikit-learn is required to compute AUC") from exc - - return float(roc_auc_score(targets, preds)) - - -def main() -> None: - args = parse_args() - device = args.device - - allowed_ids = load_allowed_ids(args.filtered_metrics, args.dice_threshold) - if not allowed_ids: - raise SystemExit("No Papila samples passed the Dice threshold.") - - papila_samples = build_papila_samples( - args.image_dir, - args.clinical_dir, - args.label_col, - args.positive_labels, - allowed_ids, - ) - if not papila_samples: - raise SystemExit("No Papila samples with labels matched the filtered metrics.") - - cache_dir = args.cache_dir - if args.use_gt_masks and cache_dir is not None: - cache_dir = cache_dir / "gt" - - if args.use_gt_masks: - geometry_provider = PapilaGTGeometryProvider(args.gt_contours_dir) - segmenter = None - else: - segmenter = UNetSegmenter( - manifest_path=args.segmenter_manifest, - device=device, - normalize=args.segmenter_normalize, - ) - seg_state = torch.load(args.segmenter_weights, map_location=device) - seg_state_dict = seg_state.get("model", seg_state) - segmenter.model.load_state_dict(seg_state_dict) - segmenter.model.to(device) - geometry_provider = UNetGeometryProvider( - segmenter=segmenter, - threshold=args.segmenter_threshold, - tta=args.segmenter_tta, - ) - - pre = RefugePreprocessing(args.refuge_root) - dummy_seg = RefugeSegmentation(pre) - backbone = build_backbone(args.backbone) - clf = RefugeClassification( - pre, - dummy_seg, - geometry_fn=geometry_provider, - cache_dir=cache_dir, - backbone=backbone, - ) - clf.crop_scale = args.crop_scale - clf.crop_size = args.crop_size - clf.eval_transform = _default_image_transform(args.crop_size) - clf.ttt_transform = clf.eval_transform - - if args.classifier_weights is not None: - clf_state = torch.load(args.classifier_weights, map_location=device) - clf.backbone.load_state_dict(clf_state["backbone"]) - clf.classifier_head.load_state_dict(clf_state["classifier"]) - clf.rotation_head.load_state_dict(clf_state["rotation"]) - if "feature_reg" in clf_state and getattr(clf, "feature_reg_head", None) is not None: - clf.feature_reg_head.load_state_dict(clf_state["feature_reg"]) - - records = clf.build_records_for_samples( - papila_samples, - crop_scale=args.crop_scale, - progress_prefix="papila_eval", - ) - if not records: - raise SystemExit("Unable to build any records; check geometry predictions or labels.") - - metrics = evaluate_records(clf, records, device=device, batch_size=args.batch_size) - print(f"Samples evaluated: {int(metrics['count'])}") - print(f"AUC: {metrics['auc']:.4f}" if not np.isnan(metrics['auc']) else "AUC: NaN") - print(f"Accuracy @0.5: {metrics['accuracy']:.4f}") - print(f"Mean glaucoma prob: {metrics['mean_prob']:.4f}") - - if args.output: - args.output.parent.mkdir(parents=True, exist_ok=True) - with args.output.open("w", newline="") as fp: - writer = csv.writer(fp) - writer.writerow(["sample_id", "prob_glaucoma", "label"]) - clf.backbone.eval() - clf.classifier_head.eval() - dataset = RefugeClassificationDataset( - records, - transform=clf.eval_transform, - polar_transform=clf.polar_transform, - size=clf.crop_size, - ) - loader = DataLoader(dataset, batch_size=args.batch_size, shuffle=False, num_workers=0) - with torch.no_grad(): - for batch in tqdm(loader, desc="Papila Output", leave=False, unit="batch"): - images = batch["image"].to(device) - polars = batch["polar"].to(device) - extra_feats = batch["features"].to(device) - ids = batch["sample_id"] - labels = batch["label"].tolist() - feats_img = clf.backbone(images) - feats = feats_img - if clf.use_polar: - feats_polar = clf.backbone(polars) - feats = torch.cat([feats, feats_polar], dim=1) - if clf.extra_feature_dim > 0: - feats = torch.cat([feats, extra_feats], dim=1) - logits = clf.classifier_head(feats) - probs = torch.softmax(logits, dim=1)[:, 1].cpu().numpy().tolist() - for sid, prob, label in zip(ids, probs, labels): - writer.writerow([sid, prob, label]) - print(f"Per-sample probabilities written to {args.output}") - - -if __name__ == "__main__": - main() diff --git a/scripts/legacy/extract_best_auc.py b/scripts/legacy/extract_best_auc.py deleted file mode 100755 index db1e68d..0000000 --- a/scripts/legacy/extract_best_auc.py +++ /dev/null @@ -1,115 +0,0 @@ -#!/usr/bin/env python3 -""" -Quick utility to recover the best epoch metrics from HyperTower run folders. - -Example: - python scripts/extract_best_auc.py analysis_data/img_only_densenet_gt_bin/img_only_densenet_gt_bin_20251028_112733 - -By default it looks for columns named like `auc_fused` (set via --metric) inside each -`fold{n}_epoch_log.csv`, returning the epoch with the highest value plus the holdout -metrics, if present. -""" - -from __future__ import annotations - -import argparse -import csv -import json -import math -from pathlib import Path -from typing import Dict, Optional, Tuple - - -def to_float(value: Optional[str]) -> Optional[float]: - if value is None: - return None - value = value.strip() - if not value: - return None - try: - out = float(value) - except ValueError: - return None - if math.isnan(out): - return None - return out - - -def best_row(path: Path, metric: str) -> Optional[Dict[str, str]]: - if not path.exists(): - return None - best: Optional[Tuple[float, int, Dict[str, str]]] = None - with path.open("r", newline="") as fp: - reader = csv.DictReader(fp) - for row in reader: - val = to_float(row.get(metric)) - if val is None: - continue - epoch = int(row.get("epoch", reader.line_num)) - if best is None or val > best[0]: - best = (val, epoch, row) - return best[2] if best else None - - -def summarize_fold(row: Dict[str, str], metric: str) -> Dict[str, float]: - data: Dict[str, float] = {} - for key in (metric, f"holdout_{metric.split('_', 1)[-1]}", "holdout_auc_img", "holdout_auc_fused"): - val = to_float(row.get(key)) - if val is not None: - data[key] = val - epoch_val = to_float(row.get("epoch")) - if epoch_val is not None: - data["epoch"] = int(epoch_val) - return data - - -def main() -> None: - ap = argparse.ArgumentParser(description="Extract best-per-fold metric from HyperTower runs.") - ap.add_argument("run_dir", type=Path, help="Run directory (contains fold*_epoch_log.csv)") - ap.add_argument("--metric", default="auc_fused", help="Metric column to maximise (default: auc_fused)") - ap.add_argument("--json", type=Path, default=None, help="Optional path to dump JSON summary") - args = ap.parse_args() - - run_dir: Path = args.run_dir - metric: str = args.metric - - if not run_dir.exists(): - raise SystemExit(f"Run directory not found: {run_dir}") - - fold_summaries: Dict[str, Dict[str, float]] = {} - metric_values = [] - - for csv_path in sorted(run_dir.glob("fold*_epoch_log.csv")): - best = best_row(csv_path, metric) - fold_name = csv_path.stem.replace("_epoch_log", "") - if best is None: - print(f"{fold_name}: no valid '{metric}' values found") - continue - summary = summarize_fold(best, metric) - fold_summaries[fold_name] = summary - val = summary.get(metric) - if val is not None: - metric_values.append(val) - holdout_val = summary.get(f"holdout_{metric.split('_', 1)[-1]}") - print(f"{fold_name}: epoch={summary.get('epoch')} {metric}={val:.4f}" if val is not None else f"{fold_name}: epoch={summary.get('epoch')}") - if holdout_val is not None: - print(f" holdout_{metric.split('_', 1)[-1]}={holdout_val:.4f}") - - if metric_values: - mean_val = sum(metric_values) / len(metric_values) - print(f"\nMean best {metric}: {mean_val:.4f}") - - if args.json: - payload = { - "run_dir": str(run_dir), - "metric": metric, - "folds": fold_summaries, - "mean_metric": (sum(metric_values) / len(metric_values)) if metric_values else None, - } - args.json.parent.mkdir(parents=True, exist_ok=True) - args.json.write_text(json.dumps(payload, indent=2)) - print(f"Summary written to {args.json}") - - -if __name__ == "__main__": - main() diff --git a/scripts/legacy/paper_metadata.py b/scripts/legacy/paper_metadata.py deleted file mode 100755 index 5e5e5cd..0000000 --- a/scripts/legacy/paper_metadata.py +++ /dev/null @@ -1,493 +0,0 @@ - -import pandas as pd -from classes import HyperTower, ClinicalData, list_names, build_papila_clinical -from pathlib import Path -import shutil, json, textwrap -from datetime import datetime -import numpy as np -from typing import Dict, List, Tuple -from sklearn.model_selection import GroupKFold -from sklearn.preprocessing import StandardScaler -from sklearn.pipeline import Pipeline -from sklearn.metrics import roc_auc_score -from sklearn.linear_model import LogisticRegression -from sklearn.neighbors import KNeighborsClassifier -from sklearn.ensemble import RandomForestClassifier -from sklearn.svm import SVC - -from sklearn.preprocessing import label_binarize - -def _proba_from_model(model, X): - if hasattr(model[-1], "predict_proba"): - return model.predict_proba(X) - dec = model.decision_function(X) - if dec.ndim == 1: # binary margins -> make 2-col - dec = np.stack([-dec, dec], axis=1) - e = np.exp(dec - dec.max(axis=1, keepdims=True)) - return e / e.sum(axis=1, keepdims=True) - -def _cv_auc_multiclass_per_class(X, y, groups, model, n_splits=5) -> np.ndarray: - """ - Returns a length-3 array of mean OvR AUCs for Class0/1/2 across GroupKFold. - Uses nan-safe means if a class is absent in a fold's test split. - """ - gkf = GroupKFold(n_splits=n_splits) - per_class_lists = [[], [], []] - for tr, te in gkf.split(X, y, groups): - model.fit(X[tr], y[tr]) - proba = _proba_from_model(model, X[te]) - y_te = y[te] - y_bin = label_binarize(y_te, classes=[0, 1, 2]) # (n,3) - for k in range(3): - yk = y_bin[:, k] - if yk.min() != yk.max(): # both classes present - per_class_lists[k].append(roc_auc_score(yk, proba[:, k])) - else: - per_class_lists[k].append(np.nan) - return np.array([np.nanmean(per_class_lists[k]) for k in range(3)], dtype=float) - -def _cv_auc_binary(X, y, groups, model, n_splits=5) -> float: - mask = np.isin(y, [0, 1]) - Xb, yb, gb = X[mask], y[mask], groups[mask] - gkf = GroupKFold(n_splits=n_splits) - aucs = [] - for tr, te in gkf.split(Xb, yb, gb): - model.fit(Xb[tr], yb[tr]) - if hasattr(model[-1], "predict_proba"): - p = model.predict_proba(Xb[te])[:, 1] - else: - p = model.decision_function(Xb[te]) - # logistic squash for safety - if np.ptp(p) > 0: - p = 1.0 / (1.0 + np.exp(-p)) - else: - p = np.full_like(p, 0.5, dtype=float) - # only compute if both classes present - if len(np.unique(yb[te])) == 2: - aucs.append(roc_auc_score(yb[te], p)) - else: - aucs.append(np.nan) - return float(np.nanmean(aucs)) - - - -# ----------------------------------- -# 1) Build Clinical Data (paper-faithful) -# ----------------------------------- -IMAGE_DIR = "Papila/FundusImages" -CLINICAL_DIR = "Papila/ClinicalData" -LABEL_COL = "Diagnosis" -CAT_COLS = ["Gender", "Phakic/Pseudophakic"] - -paper_auc = { - "TEST3_multiclass": { # Class0=Healthy, Class1=Glaucoma, Class2=Suspect - "LogReg": {"Class0": 0.67, "Class1": 0.66, "Class2": 0.67}, # from Fig. 7 (rounded) - "kNN": {"Class0": 0.72, "Class1": 0.70, "Class2": 0.76}, # your read of Fig. 7 - "RF": {"Class0": 0.66, "Class1": 0.66, "Class2": 0.67}, # from Fig. 7 (rounded) - "SVM": {"Class0": 0.66, "Class1": 0.65, "Class2": 0.66}, # from Fig. 7 (rounded) - }, - "TEST4_binary": { # Healthy vs Glaucoma (Suspects removed) - "LogReg": 0.71, # from text/Fig. 7 range midpoint - "kNN": 0.75, # your read of Fig. 7 - "RF": 0.70, # from Fig. 7 (rounded) - "SVM": 0.69, # from Fig. 7 (rounded) - } -} - - -clinical = build_papila_clinical( - image_dir=IMAGE_DIR, - clinical_dir=CLINICAL_DIR, - label_col=LABEL_COL, - cat_cols=CAT_COLS, -) - -# ----------------------------------- -# 2) Feature matrix (no MD; IOP_corr already present) -# ----------------------------------- -def build_feature_matrix(clinical) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[str]]: - """ - Returns: - X: features (N x D) - y: labels (Diagnosis: 0 healthy, 1 glaucoma, 2 suspect) - groups: patient IDs for GroupKFold - feat_names: list of feature names in X order - """ - df = clinical.df.copy() - - # Scalars used in paper-style baselines (no VF_MD) - scalars = ["Age", "dioptre_1", "dioptre_2", "astigmatism", - "Pachymetry", "Axial_Length", "IOP_corr"] - - # Categorical one-hot - cats = ["Gender", "Phakic/Pseudophakic"] - df_cats = pd.get_dummies(df[cats].astype("category"), drop_first=False, prefix=cats) - - # Combine - X = pd.concat([df[scalars], df_cats], axis=1) - - # Median impute numerics (simple, consistent) - for c in scalars: - med = pd.to_numeric(X[c], errors="coerce").median() - X[c] = pd.to_numeric(X[c], errors="coerce").fillna(med) - - y = df[LABEL_COL].astype(int).values - groups = df["Patient ID"].astype(int).values - feat_names = list(X.columns) - return X.values.astype(np.float32), y, groups, feat_names - -# ---------------------------- -# 3) Model zoo (the four methods used in the paper) -# ---------------------------- -def make_models(best_params: dict | None = None, random_state: int = 42) -> dict: - """ - Build paper-like baseline models. If best_params is provided (a dict mapping - model-name -> param dict with pipeline-style keys like 'clf__C'), those - params are applied to the corresponding pipelines. - """ - models = { - "LogReg": Pipeline([ - ("scaler", StandardScaler()), - ("clf", LogisticRegression( - max_iter=100, - solver="lbfgs", - multi_class="auto" - )) - ]), - "kNN": Pipeline([ - ("scaler", StandardScaler()), - ("clf", KNeighborsClassifier( - n_neighbors=5, - weights="uniform", - metric="minkowski", - p=2 - )) - ]), - "RF": Pipeline([ - ("clf", RandomForestClassifier( - n_estimators=100, - criterion="gini", - max_depth=None, - min_samples_split=2, - min_samples_leaf=1, - max_features="sqrt", - bootstrap=True, - # random_state left as default; set via best_params if desired - )) - ]), - "SVM": Pipeline([ - ("scaler", StandardScaler()), - ("clf", SVC( - C=1.0, - kernel="rbf", - gamma="scale", - probability=False - )) - ]), - } - - # Apply overrides if provided - if best_params: - for name, params in best_params.items(): - if name in models and params: - models[name].set_params(**params) - - return models - - -# ----------------------------------- -# 4) CV AUCs (mean over 5 folds; GroupKFold by patient) -# ----------------------------------- -def _cv_auc_multiclass(X, y, groups, model, n_splits=5) -> float: - gkf = GroupKFold(n_splits=n_splits) - aucs = [] - for tr, te in gkf.split(X, y, groups): - model.fit(X[tr], y[tr]) - if hasattr(model[-1], "predict_proba"): - proba = model.predict_proba(X[te]) - else: - dec = model.decision_function(X[te]) - if dec.ndim == 1: - dec = np.stack([-dec, dec], axis=1) - e = np.exp(dec - dec.max(axis=1, keepdims=True)) - proba = e / e.sum(axis=1, keepdims=True) - aucs.append(roc_auc_score(y[te], proba, multi_class="ovr", average="macro")) - return float(np.mean(aucs)) - - -def _cv_auc_binary(X, y, groups, model, n_splits=5) -> float: - # Keep classes 0 (healthy) and 1 (glaucoma); drop suspects (2) - mask = np.isin(y, [0, 1]) - Xb, yb, gb = X[mask], y[mask], groups[mask] - - gkf = GroupKFold(n_splits=n_splits) - aucs = [] - for tr, te in gkf.split(Xb, yb, gb): - model.fit(Xb[tr], yb[tr]) - if hasattr(model[-1], "predict_proba"): - p = model.predict_proba(Xb[te])[:, 1] - else: - p = model.decision_function(Xb[te]) - # simple logistic squashing if needed - if np.ptp(p) > 0: - p = 1.0 / (1.0 + np.exp(-p)) - else: - p = np.full_like(p, 0.5, dtype=float) - aucs.append(roc_auc_score(yb[te], p)) - return float(np.mean(aucs)) - -# ----------------------------------- -# 5) Run both tests (multiclass + binary) and print table -# ----------------------------------- -def run_papila_clinical_baselines(clinical, n_splits: int = 5, - random_state: int = 42, - best_params: dict | None = None) -> pd.DataFrame: - X, y, groups, feat_names = build_feature_matrix(clinical) - models = make_models(best_params=best_params, random_state=random_state) - - rows = [] - for name, model in models.items(): - c0, c1, c2 = _cv_auc_multiclass_per_class(X, y, groups, model, n_splits=n_splits) - auc_bin = _cv_auc_binary(X, y, groups, model, n_splits=n_splits) - rows.append({"model": name, "Class0": c0, "Class1": c1, "Class2": c2, "Binary": auc_bin}) - - df = pd.DataFrame(rows).set_index("model").sort_index() - return df - - -results = run_papila_clinical_baselines(clinical, n_splits=5) -# print(results.to_string(float_format=lambda x: f"{x:.3f}")) - - - - - - - - - - - - - - - -############################## -from sklearn.model_selection import ParameterGrid -from sklearn.base import clone -from sklearn.preprocessing import label_binarize -from sklearn.utils import check_random_state - -# ============================== -# Helper: per-class & binary AUC with GroupKFold -# ============================== -def _proba_from_model(model, X): - if hasattr(model[-1], "predict_proba"): - return model.predict_proba(X) - # decision_function fallback - dec = model.decision_function(X) - if dec.ndim == 1: # binary margin -> 2-col probs - dec = np.stack([-dec, dec], axis=1) - e = np.exp(dec - dec.max(axis=1, keepdims=True)) - return e / e.sum(axis=1, keepdims=True) - -def _cv_auc_perclass_and_binary(X, y, groups, model, n_splits=5): - """ - Returns: - per_class_auc: length-3 array (Class0, Class1, Class2) averaged over folds - binary_auc: scalar (0 vs 1) averaged over folds - """ - gkf = GroupKFold(n_splits=n_splits) - - # Hold fold-wise per-class AUCs (list of arrays of length 3) - perclass_fold_scores = [] - binary_fold_scores = [] - - for tr, te in gkf.split(X, y, groups): - y_te = y[te] - # Multiclass per-class (OvR) - model.fit(X[tr], y[tr]) - proba = _proba_from_model(model, X[te]) - - # One-vs-rest per-class AUCs (skip a class if absent in test fold) - y_bin = label_binarize(y_te, classes=[0, 1, 2]) # shape (n, 3) - perclass_scores = [] - for k in range(3): - yk = y_bin[:, k] - # Only compute if both 0 and 1 are present - if yk.min() != yk.max(): - perclass_scores.append(roc_auc_score(yk, proba[:, k])) - else: - perclass_scores.append(np.nan) - perclass_fold_scores.append(perclass_scores) - - # Binary AUC (0 vs 1; drop class 2) - mask = np.isin(y_te, [0, 1]) - if mask.sum() > 0 and len(np.unique(y_te[mask])) == 2: - # we need probabilities/margins for class 1 among (0,1) - # Map proba[:, 1] if the model was trained 3-way; we restrict te samples to 0/1 - binary_p = proba[mask, 1] - binary_y = y_te[mask] - binary_fold_scores.append(roc_auc_score(binary_y, binary_p)) - else: - binary_fold_scores.append(np.nan) - - # Average over folds (ignore NaNs if a class was missing in a fold) - perclass_arr = np.array(perclass_fold_scores, dtype=float) # (n_folds, 3) - per_class_auc = np.nanmean(perclass_arr, axis=0) - binary_auc = float(np.nanmean(np.array(binary_fold_scores, dtype=float))) - return per_class_auc, binary_auc - -# ============================== -# Distance-to-paper objective -# ============================== -def _distance_to_paper(model_name: str, - per_class_auc: np.ndarray, - binary_auc: float, - paper_auc: Dict, - w_mc: float = 1.0, - w_bin: float = 1.0) -> float: - mc_targets = paper_auc["TEST3_multiclass"][model_name] - tvec = np.array([mc_targets["Class0"], mc_targets["Class1"], mc_targets["Class2"]], dtype=float) - mc_diff = np.nanmean(np.abs(per_class_auc - tvec)) # mean absolute difference over 3 classes - - bin_target = paper_auc["TEST4_binary"][model_name] - bin_diff = abs(binary_auc - bin_target) - - return float(w_mc * mc_diff + w_bin * bin_diff) - -# ============================== -# Parameter grids (paper-ish, not crazy-large) -# ============================== -def get_param_grids() -> Dict[str, List[dict]]: - return { - "LogReg": [ - { - "clf__C": [0.01, 0.1, 1.0, 3.0, 10.0], - "clf__class_weight": [None, "balanced"], - "clf__max_iter": [200, 500], - # lbfgs + l2 is implied - } - ], - "kNN": [ - { - "clf__n_neighbors": [3, 5, 7, 9, 11], - "clf__weights": ["uniform", "distance"], - "clf__p": [1, 2], # Manhattan vs Euclidean - } - ], - "RF": [ - { - "clf__n_estimators": [200, 500, 1000], - "clf__max_depth": [None, 5, 10, 20], - "clf__max_features": ["sqrt", "log2", 0.5], - "clf__min_samples_leaf": [1, 2, 5], - "clf__class_weight": [None, "balanced"], - # If you want determinism add: "clf__random_state": [42] - } - ], - "SVM": [ - { - "clf__C": [0.1, 1.0, 3.0, 10.0], - "clf__gamma": ["scale", "auto", 0.1, 0.01, 0.001], - "clf__kernel": ["rbf"], # fixed to rbf as in paper-like default - } - ], - } - -# ============================== -# Grid search loop minimizing distance-to-paper -# ============================== -def search_params_to_match_paper( - clinical, - models: Dict[str, Pipeline], - paper_auc: Dict, - n_splits: int = 5, - w_mc: float = 1.0, - w_bin: float = 1.0, - verbose: bool = True, -) -> Tuple[pd.DataFrame, Dict[str, dict]]: - X, y, groups, feat_names = build_feature_matrix(clinical) - grids = get_param_grids() - - summary_rows = [] - best_params_by_model = {} - - for name, base_model in models.items(): - if name not in grids: - if verbose: - print(f"[warn] No grid for {name}, skipping.") - continue - - best_loss = np.inf - best_params = None - best_mc = None - best_bin = None - - for param_set in ParameterGrid(grids[name]): - model = clone(base_model).set_params(**param_set) - per_class_auc, binary_auc = _cv_auc_perclass_and_binary( - X, y, groups, model, n_splits=n_splits - ) - loss = _distance_to_paper( - name, per_class_auc, binary_auc, paper_auc, w_mc=w_mc, w_bin=w_bin - ) - - if verbose: - mc_str = " / ".join(f"{a:.3f}" if np.isfinite(a) else "nan" for a in per_class_auc) - print(f"[{name}] params={param_set} | mc per-class={mc_str} | bin={binary_auc:.3f} | loss={loss:.4f}") - - if loss < best_loss: - best_loss = loss - best_params = param_set - best_mc = per_class_auc - best_bin = binary_auc - - # store - best_params_by_model[name] = best_params - summary_rows.append({ - "model": name, - "best_loss": best_loss, - "best_params": json.dumps(best_params), - "mc_Class0": float(best_mc[0]), - "mc_Class1": float(best_mc[1]), - "mc_Class2": float(best_mc[2]), - "binary_auc": float(best_bin), - "paper_mc_Class0": paper_auc["TEST3_multiclass"][name]["Class0"], - "paper_mc_Class1": paper_auc["TEST3_multiclass"][name]["Class1"], - "paper_mc_Class2": paper_auc["TEST3_multiclass"][name]["Class2"], - "paper_binary": paper_auc["TEST4_binary"][name], - }) - - df = pd.DataFrame(summary_rows).set_index("model").sort_values("best_loss") - return df, best_params_by_model - -# ============================== -# Run the search -# ============================== -models = make_models(random_state=42) -df_match, best_params = search_params_to_match_paper( - clinical=clinical, - models=models, - paper_auc=paper_auc, - n_splits=5, - w_mc=1.0, # weight multiclass distance - w_bin=1.0, # weight binary distance - verbose=True -) - -# print("\n=== Best params found (by minimal distance-to-paper) ===") -# print(df_match[["best_loss","best_params","mc_Class0","mc_Class1","mc_Class2","binary_auc", -# "paper_mc_Class0","paper_mc_Class1","paper_mc_Class2","paper_binary"]]) - -# print("\nBest param dicts:") -for k, v in best_params.items(): - print(k, "->", v) - -results2 = run_papila_clinical_baselines(clinical, n_splits=5, random_state=42, best_params=best_params) -print(f"Default Settings: {results.round(2)}") -print(f"Best Params Settings: {results2.round(2)}") -print(f" Paper Results: {pd.DataFrame({ - model: {**vals, "Binary": paper_auc["TEST4_binary"][model]} - for model, vals in paper_auc["TEST3_multiclass"].items() -}).T[["Class0","Class1","Class2","Binary"]]}") \ No newline at end of file diff --git a/scripts/legacy/run_all_sweep.sh b/scripts/legacy/run_all_sweep.sh deleted file mode 100755 index 4bfc813..0000000 --- a/scripts/legacy/run_all_sweep.sh +++ /dev/null @@ -1,119 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Usage: -# bash scripts/run_all_sweep.sh --epochs 25 --n-splits 5 --batch-size 8 [extra args] -# -# Merged sweep: runs the SE attention grid (bridge/tower/both × R=8/16/32), -# skipping tower-only non-normalized variants (tower normalization is a no-op), -# and includes binary eval counterparts for each baseline run. It also submits -# the full gradual-thaw grid (multiclass + binary variants). - -ARGS=("$@") - -run() { - local SHORT="$1"; shift - echo "=== Running: $SHORT ===" - # Skip if a summary for this shortname already exists - if ls "analysis_data/${SHORT}_"*.md >/dev/null 2>&1; then - echo "… skipping ${SHORT} (summary already present)" - return 0 - fi - python3 scripts/run_multifold.py \ - --shortname "$SHORT" \ - "$@" \ - "${ARGS[@]}" || true -} - -echo "--- SE Grid (bridge/tower/both × R=8/16/32; tower nonorm skipped) ---" -for R in 8 16 32; do - # Bridge-only - run "se_bridge_R${R}_norm" --se-where bridge --se-reduction ${R} --se-pre-norm --checkpoint-best - run "se_bridge_R${R}_norm_bin" --se-where bridge --se-reduction ${R} --se-pre-norm --checkpoint-best --eval_mode binary - run "se_bridge_R${R}_nonorm" --se-where bridge --se-reduction ${R} --no-se-pre-norm --checkpoint-best - run "se_bridge_R${R}_nonorm_bin" --se-where bridge --se-reduction ${R} --no-se-pre-norm --checkpoint-best --eval_mode binary - - # Tower-only - run "se_tower_R${R}_norm" --se-where tower --se-reduction-tower ${R} --se-pre-norm-tower --checkpoint-best - run "se_tower_R${R}_norm_bin" --se-where tower --se-reduction-tower ${R} --se-pre-norm-tower --checkpoint-best --eval_mode binary - - # Tower+Bridge - run "se_tower_bridge_R${R}_norm" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --se-pre-norm --se-pre-norm-tower --checkpoint-best - run "se_tower_bridge_R${R}_norm_bin" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --se-pre-norm --se-pre-norm-tower --checkpoint-best --eval_mode binary - run "se_tower_bridge_R${R}_nonorm" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --no-se-pre-norm --no-se-pre-norm-tower --checkpoint-best - run "se_tower_bridge_R${R}_nonorm_bin" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --no-se-pre-norm --no-se-pre-norm-tower --checkpoint-best --eval_mode binary -done - -THAW_COMMON_ARGS=( - --gradual-thaw - --thaw-phase-duration 5 - --thaw-ratio 0.33 - --thaw-start-epoch 5 - --early-stop - --early-patience 5 -) - -echo "--- Gradual Thaw Grid (multiclass + binary) ---" - -# Bridge-only thaw runs (norm and nonorm) -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=(--se-where bridge --se-reduction "$R" --se-pre-norm --checkpoint-best) - else - FLAGS=(--se-where bridge --se-reduction "$R" --no-se-pre-norm --checkpoint-best) - fi - run "thaw_se_bridge_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" - run "thawbin_se_bridge_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" --eval_mode binary - done -done - -# Tower-only thaw runs (norm and nonorm) -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=(--se-where tower --se-reduction-tower "$R" --se-pre-norm-tower --checkpoint-best) - else - FLAGS=(--se-where tower --se-reduction-tower "$R" --no-se-pre-norm-tower --checkpoint-best) - fi - run "thaw_se_tower_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" - run "thawbin_se_tower_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" --eval_mode binary - done -done - -# Tower+bridge thaw runs (norm and nonorm) -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=( - --se-where both - --se-reduction "$R" - --se-reduction-tower "$R" - --se-pre-norm - --se-pre-norm-tower - --checkpoint-best - ) - else - FLAGS=( - --se-where both - --se-reduction "$R" - --se-reduction-tower "$R" - --no-se-pre-norm - --no-se-pre-norm-tower - --checkpoint-best - ) - fi - run "thaw_se_tower_bridge_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" - run "thawbin_se_tower_bridge_R${R}_${MODE}" "${FLAGS[@]}" "${THAW_COMMON_ARGS[@]}" --eval_mode binary - done -done - -echo "Merged sweep submitted. Check analysis_data/* and models/* for outputs." diff --git a/scripts/legacy/run_gradual_thaw_top5.sh b/scripts/legacy/run_gradual_thaw_top5.sh deleted file mode 100755 index 5109c24..0000000 --- a/scripts/legacy/run_gradual_thaw_top5.sh +++ /dev/null @@ -1,79 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Usage: -# bash scripts/run_gradual_thaw_top5.sh --epochs 20 --n-splits 5 --batch-size 8 [extra args] -# -# Runs the full gradual-thaw grid aligned with the SE sweep (bridge/tower/both × R=8/16/32 × norm vs nonorm). - -ARGS=("$@") - -run() { - local SHORT="$1"; shift - echo "=== Running: $SHORT ===" - if ls "analysis_data/${SHORT}_"*.md >/dev/null 2>&1; then - echo "… skipping ${SHORT} (summary already present)" - return 0 - fi - python3 scripts/run_multifold.py \ - --shortname "$SHORT" \ - --gradual-thaw --thaw-phase-duration 5 --thaw-ratio 0.33 --thaw-start-epoch 5 \ - --early-stop --early-patience 5 \ - "$@" \ - "${ARGS[@]}" || true -} - -# Bridge-only thaw runs -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=(--se-where bridge --se-reduction "$R" --se-pre-norm --checkpoint-best) - else - FLAGS=(--se-where bridge --se-reduction "$R" --no-se-pre-norm --checkpoint-best) - fi - run "thaw_se_bridge_R${R}_${MODE}" "${FLAGS[@]}" - run "thawbin_se_bridge_R${R}_${MODE}" "${FLAGS[@]}" --eval_mode binary - done -done - -# Tower-only thaw runs -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=(--se-where tower --se-reduction-tower "$R" --se-pre-norm-tower --checkpoint-best) - else - FLAGS=(--se-where tower --se-reduction-tower "$R" --no-se-pre-norm-tower --checkpoint-best) - fi - run "thaw_se_tower_R${R}_${MODE}" "${FLAGS[@]}" - run "thawbin_se_tower_R${R}_${MODE}" "${FLAGS[@]}" --eval_mode binary - done -done - -# Tower+bridge thaw runs -for R in 8 16 32; do - for MODE in norm nonorm; do - if [[ "$MODE" == "norm" ]]; then - FLAGS=( - --se-where both - --se-reduction "$R" - --se-reduction-tower "$R" - --se-pre-norm - --se-pre-norm-tower - --checkpoint-best - ) - else - FLAGS=( - --se-where both - --se-reduction "$R" - --se-reduction-tower "$R" - --no-se-pre-norm - --no-se-pre-norm-tower - --checkpoint-best - ) - fi - run "thaw_se_tower_bridge_R${R}_${MODE}" "${FLAGS[@]}" - run "thawbin_se_tower_bridge_R${R}_${MODE}" "${FLAGS[@]}" --eval_mode binary - done -done - -echo "Gradual thaw grid submitted. Check analysis_data/* and models/* for outputs." diff --git a/scripts/legacy/run_multifold_gui.py b/scripts/legacy/run_multifold_gui.py deleted file mode 100755 index faf9e30..0000000 --- a/scripts/legacy/run_multifold_gui.py +++ /dev/null @@ -1,15 +0,0 @@ -#!/usr/bin/env python3 -"""Launch the Tkinter front-end for run_multifold.""" - -import sys -from pathlib import Path - -REPO_ROOT = Path(__file__).resolve().parents[1] -if str(REPO_ROOT) not in sys.path: - sys.path.insert(0, str(REPO_ROOT)) - -from classes.frontend import launch_frontend - - -if __name__ == "__main__": - launch_frontend() diff --git a/scripts/legacy/run_multimodel.py b/scripts/legacy/run_multimodel.py deleted file mode 100755 index 726faa0..0000000 --- a/scripts/legacy/run_multimodel.py +++ /dev/null @@ -1,145 +0,0 @@ -#!/usr/bin/env python3 -import argparse, subprocess, sys, time, json -from pathlib import Path - -# Backbones in the paper that torchvision supports -BACKBONES = [ - "efficientnet_b0", - "resnet50", - "densenet121", - "vgg16", - "mobilenet_v2", - "inception_v3", - # (Xception omitted; not in torchvision — add via timm later if needed) -] - -MODES = [ - ("multiclass", ["Healthy", "Glaucoma", "Suspect"]), - ("binary", ["Healthy", "Glaucoma"]), -] - -def run(cmd): - print("\n$ " + " ".join(map(str, cmd))) - res = subprocess.run(cmd, check=True) - return res.returncode - -def main(): - ap = argparse.ArgumentParser(description="Run all paper CNNs across folds in multiclass + binary, then compile plots.") - ap.add_argument("--epochs", type=int, default=5, help="Epochs per fold (fast sanity first).") - ap.add_argument("--shortname", type=str, default="papergrid", help="Prefix for run IDs.") - ap.add_argument("--n-splits", type=int, default=5, help="Number of folds.") - ap.add_argument("--fusion-mode", type=str, default="fused", choices=["image_only","fused","metadata_only","vote"], - help="Paper CNNs are image-only; leave as image_only unless you’re testing others.") - ap.add_argument("--freeze-ratio", type=float, default=0.0, help="0.0 = full fine-tune (as in the paper).") - # You can override data roots if needed - ap.add_argument("--image-dir", default="Papila/FundusImages") - ap.add_argument("--clinical-dir", default="Papila/ClinicalData") - ap.add_argument("--label-col", default="Diagnosis") - ap.add_argument("--cat-cols", nargs="*", default=["Gender", "Phakic/Pseudophakic"]) - args = ap.parse_args() - - ts = time.strftime("%Y%m%d_%H%M%S") - master_tag = f"{args.shortname}_{ts}" - master_dir = Path("analysis_data") / master_tag - master_dir.mkdir(parents=True, exist_ok=True) - - # Keep a log of all subruns for the master report - index = [] - - for backbone in BACKBONES: - for eval_mode, class_names in MODES: - # build a child shortname per (backbone, mode) - sub_prefix = f"{args.shortname}_{backbone}_{eval_mode}" - cmd = [ - sys.executable, "scripts/run_multifold.py", - "--backbone", backbone, - "--freeze-ratio", str(args.freeze_ratio), - "--fusion-mode", args.fusion_mode, - "--epochs", str(args.epochs), - "--n-splits", str(args.n_splits), - "--shortname", sub_prefix, - "--eval_mode", eval_mode, - "--image-dir", args.image_dir, - "--clinical-dir", args.clinical_dir, - "--label-col", args.label_col, - ] - - # class names by mode (ensures plot legends are correct) - cmd += ["--class-names", *class_names] - - plot_head_map = { - "image_only": "image", - "fused" : "fused", - "metadata_only": "metadata", - "vote": "fused", - } - # We always aggregate/plot the image head for paper CNNs - cmd += ["--plot-head", plot_head_map.get(args.fusion_mode)] - - # Delegate the whole run to run_multifold.py - run(cmd) - - # Discover the child run folder (the newest folder matching the shortname prefix) - # We do this because run_multifold appends its own timestamp. - adir = Path("analysis_data") - children = sorted([p for p in adir.glob(f"{sub_prefix}_*") if p.is_dir()]) - if not children: - print(f"[WARN] No analysis_data folder found for {sub_prefix}; skipping index entry.") - continue - run_dir = children[-1] - summary_json = run_dir / "summary.json" - plots_dir = run_dir / "plots" - - # Record entry - entry = { - "backbone": backbone, - "eval_mode": eval_mode, - "run_dir": str(run_dir), - "summary_json": str(summary_json) if summary_json.exists() else None, - "plots": { - "mean": str(plots_dir / "roc_image_mean_ovr.png"), - "overlay": str(plots_dir / "roc_image_perfold_overlay.png"), - } - } - # Try to read AUCs - try: - if summary_json.exists(): - entry.update(json.loads(summary_json.read_text())) - except Exception: - pass - index.append(entry) - - # Write a master JSON + markdown report - (master_dir / "index.json").write_text(json.dumps(index, indent=2), encoding="utf-8") - - # Simple markdown table of results with links - lines = [ - f"# Multimodel grid — {master_tag}", - "", - f"- Epochs per fold: **{args.epochs}**", - f"- Folds: **{args.n_splits}**", - f"- Fusion mode: **{args.fusion_mode}** (paper CNNs = image-only)", - f"- Freeze ratio: **{args.freeze_ratio}**", - "", - "| Backbone | Mode | Mean AUC (macro/mc or ROC-AUC/bin) | Plots | Run folder |", - "|---|---|---:|---|---|", - ] - for e in index: - auc_mean = e.get("macro_ovr_auc_mean", None) - if auc_mean is not None: - auc_str = f"{auc_mean:.3f}" - else: - auc_str = "—" - mean_png = e["plots"]["mean"] - overlay_png = e["plots"]["overlay"] - plots_md = f"[mean]({mean_png}) / [overlay]({overlay_png})" - lines.append( - f"| `{e['backbone']}` | `{e['eval_mode']}` | {auc_str} | {plots_md} | `{e['run_dir']}` |" - ) - (master_dir / "README.md").write_text("\n".join(lines) + "\n", encoding="utf-8") - - print(f"\nAll done.\n- Master index: {master_dir/'index.json'}\n- Report: {master_dir/'README.md'}") - print(f"- Individual runs live under analysis_data/ with plots and summaries.") - -if __name__ == "__main__": - main() diff --git a/scripts/legacy/run_se_sweep.sh b/scripts/legacy/run_se_sweep.sh deleted file mode 100755 index 1bf307c..0000000 --- a/scripts/legacy/run_se_sweep.sh +++ /dev/null @@ -1,56 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Usage: -# bash scripts/run_se_sweep.sh --epochs 50 --n-splits 5 --batch-size 8 --eval_mode multiclass [extra args] -# -# This will launch a series of runs covering the grid from the slide: -# - Bridge-only R8/R16/R32 (normalized and non-normalized) -# - Tower-only R8/R16/R32 (normalized and non-normalized) -# - Tower+Bridge R8/R16/R32 (normalized and non-normalized) - -ARGS=("$@") - -run() { - local SHORT="$1"; shift - echo "=== Running: $SHORT ===" - # Skip if a summary for this shortname already exists - if ls "analysis_data/${SHORT}_"*.md >/dev/null 2>&1; then - echo "… skipping ${SHORT} (summary already present)" - return 0 - fi - python3 scripts/run_multifold.py \ - --shortname "$SHORT" \ - "$@" \ - "${ARGS[@]}" || true -} - -# Bridge-only (normalized + non-normalized) -for R in 8 16 32; do - run "se_bridge_R${R}_norm" --se-where bridge --se-reduction ${R} --se-pre-norm --checkpoint-best - run "se_bridge_R${R}_nonorm" --se-where bridge --se-reduction ${R} --no-se-pre-norm --checkpoint-best -done - -# Tower-only (normalized + non-normalized) -for R in 8 16 32; do - run "se_tower_R${R}_norm" \ - --se-where tower --se-reduction-tower ${R} --se-pre-norm-tower \ - --checkpoint-best - run "se_tower_R${R}_nonorm" \ - --se-where tower --se-reduction-tower ${R} --no-se-pre-norm-tower \ - --checkpoint-best -done - -# Tower+Bridge (normalized + non-normalized) -for R in 8 16 32; do - # normalized (both pre-norm on) - run "se_tower_bridge_R${R}_norm" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --se-pre-norm --se-pre-norm-tower --checkpoint-best - # non-normalized (both pre-norm off) - run "se_tower_bridge_R${R}_nonorm" \ - --se-where both --se-reduction ${R} --se-reduction-tower ${R} \ - --no-se-pre-norm --no-se-pre-norm-tower --checkpoint-best -done - -echo "Sweep submitted. Check analysis_data/* and models/* for outputs." diff --git a/scripts/main/pipeline.ipynb b/scripts/main/pipeline.ipynb index dc910c5..642f963 100644 --- a/scripts/main/pipeline.ipynb +++ b/scripts/main/pipeline.ipynb @@ -1287,7 +1287,7 @@ "outputs": [], "source": [ "# 3a) No crop — original full-size images\n", - "!python scripts/basic_analysis/compare_hypertower_modes.py \\\n", + "!python scripts/main/v2/multirun_hypertower.py \\\n", " --tower-modes single ensemble \\\n", " --eval-modes binary multiclass \\\n", " --epochs 40 --n-splits 5 \\\n", @@ -1305,7 +1305,7 @@ "outputs": [], "source": [ "# 3b) GT crop — expert segmentation masks crop the optic disc region\n", - "!python scripts/basic_analysis/compare_hypertower_modes.py \\\n", + "!python scripts/main/v2/multirun_hypertower.py \\\n", " --tower-modes single ensemble \\\n", " --eval-modes binary multiclass \\\n", " --epochs 40 --n-splits 5 \\\n", @@ -1323,7 +1323,7 @@ "outputs": [], "source": [ "# 3c) UNet crop — trained segmenter crops the optic disc region\n", - "!python scripts/basic_analysis/compare_hypertower_modes.py \\\n", + "!python scripts/main/v2/multirun_hypertower.py \\\n", " --tower-modes single ensemble \\\n", " --eval-modes binary multiclass \\\n", " --epochs 40 --n-splits 5 \\\n", diff --git a/scripts/basic_analysis/compare_hypertower_modes.py b/scripts/main/v2/multirun_hypertower.py similarity index 100% rename from scripts/basic_analysis/compare_hypertower_modes.py rename to scripts/main/v2/multirun_hypertower.py diff --git a/scripts/main/v2/run_binary_v2.2.sh b/scripts/main/v2/run_binary_v2.2.sh deleted file mode 100755 index 1a1cc2e..0000000 --- a/scripts/main/v2/run_binary_v2.2.sh +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Binary runs v2.2 (4 total): -# UNet crop: single | fused head -# GT crop: single | fused head - -ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)" -cd "$ROOT_DIR" - -MANIFEST="manifest.csv" -UNET_WEIGHTS="models/v2/refuge/segmentation/per_image/best.pt" - -COMMON=( - --epochs 40 - --n-splits 5 - --batch-size 8 - --backbone refugelike - --eval-mode binary - --single-warmup-tower-epochs 4 - --single-warmup-fused-epochs 4 - --img-crop-manifest "$MANIFEST" -) - -UNET_CROP=( - --img-crop-weights "$UNET_WEIGHTS" -) - -GT_CROP=( - --img-crop-gt -) - -# ── UNet crop ──────────────────────────────────────────────────────────────── - -echo "[1/4] UNet crop — binary, single..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${UNET_CROP[@]}" \ - --tower-mode single \ - --run-name v2.2_single_binary_unet_40ep_5fold - -echo "[2/4] UNet crop — binary, ensemble + fused head..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${UNET_CROP[@]}" \ - --tower-mode ensemble \ - --fused-head --fusion-epochs 20 \ - --run-name v2.2_fused_binary_unet_40ep_5fold - -# ── GT crop ────────────────────────────────────────────────────────────────── - -echo "[3/4] GT crop — binary, single..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${GT_CROP[@]}" \ - --tower-mode single \ - --run-name v2.2_single_binary_gt_40ep_5fold - -echo "[4/4] GT crop — binary, ensemble + fused head..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${GT_CROP[@]}" \ - --tower-mode ensemble \ - --fused-head --fusion-epochs 20 \ - --run-name v2.2_fused_binary_gt_40ep_5fold - -echo "Binary v2.2 runs complete." diff --git a/scripts/main/v2/run_imgonly_v2.21.sh b/scripts/main/v2/run_imgonly_v2.21.sh deleted file mode 100755 index c0f76c9..0000000 --- a/scripts/main/v2/run_imgonly_v2.21.sh +++ /dev/null @@ -1,68 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Image-only runs v2.21 (6 total): -# No crop: binary | multiclass -# GT crop: binary | multiclass -# UNet crop: binary | multiclass -# -# Purpose: isolate the effect of ROI cropping at the single-CNN level, -# without any MD tower contribution (bridge-mode=image_only). - -ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)" -cd "$ROOT_DIR" - -MANIFEST="manifest.csv" -UNET_WEIGHTS="models/v2/refuge/segmentation/per_image/best.pt" - -COMMON=( - --epochs 40 - --n-splits 5 - --batch-size 8 - --backbone refugelike - --tower-mode single - --bridge-mode image_only - --single-warmup-tower-epochs 4 - --single-warmup-fused-epochs 0 - --img-crop-manifest "$MANIFEST" -) - -# ── No crop ────────────────────────────────────────────────────────────────── - -echo "[1/6] No crop — binary, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode binary \ - --run-name v2.21_imgonly_binary_nocrop_40ep_5fold - -echo "[2/6] No crop — multiclass, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode multiclass \ - --run-name v2.21_imgonly_multiclass_nocrop_40ep_5fold - -# ── GT crop ────────────────────────────────────────────────────────────────── - -echo "[3/6] GT crop — binary, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode binary --img-crop-gt \ - --run-name v2.21_imgonly_binary_gt_40ep_5fold - -echo "[4/6] GT crop — multiclass, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode multiclass --img-crop-gt \ - --run-name v2.21_imgonly_multiclass_gt_40ep_5fold - -# ── UNet crop ──────────────────────────────────────────────────────────────── - -echo "[5/6] UNet crop — binary, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode binary \ - --img-crop-weights "$UNET_WEIGHTS" \ - --run-name v2.21_imgonly_binary_unet_40ep_5fold - -echo "[6/6] UNet crop — multiclass, image-only..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" --eval-mode multiclass \ - --img-crop-weights "$UNET_WEIGHTS" \ - --run-name v2.21_imgonly_multiclass_unet_40ep_5fold - -echo "Image-only v2.21 runs complete." diff --git a/scripts/main/v2/run_multiclass_v2.2.sh b/scripts/main/v2/run_multiclass_v2.2.sh deleted file mode 100755 index 8e0302e..0000000 --- a/scripts/main/v2/run_multiclass_v2.2.sh +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env bash -set -euo pipefail - -# Multiclass runs v2.2 (4 total): -# UNet crop: single | fused head -# GT crop: single | fused head - -ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/../../.." && pwd)" -cd "$ROOT_DIR" - -MANIFEST="manifest.csv" -UNET_WEIGHTS="models/v2/refuge/segmentation/per_image/best.pt" - -COMMON=( - --epochs 40 - --n-splits 5 - --batch-size 8 - --backbone refugelike - --eval-mode multiclass - --single-warmup-tower-epochs 4 - --single-warmup-fused-epochs 4 - --img-crop-manifest "$MANIFEST" -) - -UNET_CROP=( - --img-crop-weights "$UNET_WEIGHTS" -) - -GT_CROP=( - --img-crop-gt -) - -# ── UNet crop ──────────────────────────────────────────────────────────────── - -echo "[1/4] UNet crop — multiclass, single..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${UNET_CROP[@]}" \ - --tower-mode single \ - --run-name v2.2_single_multiclass_unet_40ep_5fold - -echo "[2/4] UNet crop — multiclass, ensemble + fused head..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${UNET_CROP[@]}" \ - --tower-mode ensemble \ - --fused-head --fusion-epochs 20 \ - --run-name v2.2_fused_multiclass_unet_40ep_5fold - -# ── GT crop ────────────────────────────────────────────────────────────────── - -echo "[3/4] GT crop — multiclass, single..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${GT_CROP[@]}" \ - --tower-mode single \ - --run-name v2.2_single_multiclass_gt_40ep_5fold - -echo "[4/4] GT crop — multiclass, ensemble + fused head..." -python3 scripts/main/v2/run_multifold_v2.py \ - "${COMMON[@]}" "${GT_CROP[@]}" \ - --tower-mode ensemble \ - --fused-head --fusion-epochs 20 \ - --run-name v2.2_fused_multiclass_gt_40ep_5fold - -echo "Multiclass v2.2 runs complete." diff --git a/scripts/main/v2/run_v2.3_ensemble_nocrop.sh b/scripts/main/v2/run_v2.3_ensemble_nocrop.sh deleted file mode 100755 index 7e3f2b9..0000000 --- a/scripts/main/v2/run_v2.3_ensemble_nocrop.sh +++ /dev/null @@ -1,20 +0,0 @@ -#!/usr/bin/env bash -set -e - -echo "=== v2.3 ensemble binary nocrop ===" -python scripts/basic_analysis/compare_hypertower_modes.py \ - --tower-modes ensemble --eval-modes binary \ - --epochs 40 --n-splits 5 \ - --backbone refugelike \ - --img-crop-manifest analysis_data/unet_manifest.csv \ - --run-name v2.3_ensemble_binary_nocrop - -echo "=== v2.3 ensemble multiclass nocrop ===" -python scripts/basic_analysis/compare_hypertower_modes.py \ - --tower-modes ensemble --eval-modes multiclass \ - --epochs 40 --n-splits 5 \ - --backbone refugelike \ - --img-crop-manifest analysis_data/unet_manifest.csv \ - --run-name v2.3_ensemble_multiclass_nocrop - -echo "=== done ===" diff --git a/scripts/main/v2/run_v2.3_fused_nocrop.sh b/scripts/main/v2/run_v2.3_fused_nocrop.sh deleted file mode 100755 index f8181ec..0000000 --- a/scripts/main/v2/run_v2.3_fused_nocrop.sh +++ /dev/null @@ -1,24 +0,0 @@ -#!/usr/bin/env bash -set -e - -echo "=== v2.3 fused binary nocrop ===" -python scripts/basic_analysis/compare_hypertower_modes.py \ - --tower-modes ensemble --eval-modes binary \ - --epochs 40 --n-splits 5 \ - --backbone refugelike \ - --img-crop-manifest analysis_data/unet_manifest.csv \ - --warmup-md-epochs 50 \ - --fused-head \ - --run-name v2.3_fused_binary_nocrop - -echo "=== v2.3 fused multiclass nocrop ===" -python scripts/basic_analysis/compare_hypertower_modes.py \ - --tower-modes ensemble --eval-modes multiclass \ - --epochs 40 --n-splits 5 \ - --backbone refugelike \ - --img-crop-manifest analysis_data/unet_manifest.csv \ - --warmup-md-epochs 50 \ - --fused-head \ - --run-name v2.3_fused_multiclass_nocrop - -echo "=== done ==="