2026001

.gitignore

@@ -36,7 +36,7 @@ secrets*
 # Data (large feature files — regenerate via scripts)
 features/*.npz
 features/*.npy
-
+features/*
 # Outputs
 
 

classes/__init__.py

@@ -1,3 +1,5 @@
 from .features import FeatureExtractor
 from .patient_identifier import PatientIdentifier
 from .classifier import PatientLeakageClassifier
+from .nifti_dataset import NiftiSliceDataset
+from .siamese import SiamesePatientMatcher, SiameseCNN

classes/classifier.py

@@ -61,7 +61,7 @@ class PatientLeakageClassifier:
     # ------------------------------------------------------------------
 
     def run(self, model_name, seed, split_type="patient",
-            nfeatures_list=None, gamma_logspace=None):
+            nfeatures_list=None, gamma_logspace=None, return_predictions=False):
         """Run full pipeline: train/test split → RF importance →
         GridSearchCV over nfeatures×gamma → retrain → test.
 
@@ -72,6 +72,9 @@ class PatientLeakageClassifier:
         split_type : str       "image" or "patient".
         nfeatures_list : list  Feature counts to search over.
         gamma_logspace : tuple Args for np.logspace (start, stop, n).
+        return_predictions : bool
+            If True, also include the test labels/predictions ("y_true",
+            "y_pred") so callers can build a confusion matrix without refitting.
 
         Returns
         -------
@@ -84,55 +87,12 @@ class PatientLeakageClassifier:
             gamma_logspace = DEFAULT_GAMMA_RANGE
 
         X, Y, groups = self._load_data(model_name)
+        (X_tr, X_te, y_tr, y_te, g_tr, cv,
+         n_tr_patients, n_te_patients) = self._split(X, Y, groups, seed, split_type)
 
-        # --- Train/test split ---
-        if split_type == "image":
-            X_tr, X_te, y_tr, y_te = train_test_split(
-                X, Y, test_size=0.2, random_state=seed, stratify=Y)
-            cv = 5
-            g_tr = None
-            n_tr_patients = None
-            n_te_patients = None
-        else:
-            unique_g = np.unique(groups)
-            g_labels = np.array([Y[groups == g][0] for g in unique_g])
-            tr_g, te_g = train_test_split(
-                unique_g, test_size=0.2, random_state=seed, stratify=g_labels)
-            tr_mask = np.isin(groups, tr_g)
-            te_mask = np.isin(groups, te_g)
-            X_tr, X_te = X[tr_mask], X[te_mask]
-            y_tr, y_te = Y[tr_mask], Y[te_mask]
-            g_tr = groups[tr_mask]
-            n_tr_patients = len(np.unique(g_tr))
-            n_te_patients = len(np.unique(groups[te_mask]))
-            cv = StratifiedGroupKFold(n_splits=5, shuffle=True,
-                                      random_state=seed)
-
-        # --- RF feature importance ---
-        fs = RandomForestClassifier(n_estimators=500, random_state=15,
-                                    class_weight="balanced", n_jobs=self.n_jobs)
-        fs.fit(X_tr, y_tr)
-        sorted_idx = np.argsort(fs.feature_importances_)[::-1]
-
-        # --- Grid search over nfeatures × gamma ---
-        pipe = Pipeline([("scaler", StandardScaler()),
-                         ("svm", SVC(kernel="rbf", C=10))])
-        valid_nfeats = [n for n in nfeatures_list if n <= X_tr.shape[1]]
-        best_cv, best_n, best_g = 0, None, None
-
-        for nfeat in valid_nfeats:
-            sel = sorted_idx[:nfeat]
-            pg = {"svm__gamma": (1.0 / nfeat) * np.logspace(*gamma_logspace)}
-            grid = GridSearchCV(pipe, pg, cv=cv, scoring="accuracy",
-                                n_jobs=self.n_jobs, verbose=0)
-            if g_tr is not None:
-                grid.fit(X_tr[:, sel], y_tr, groups=g_tr)
-            else:
-                grid.fit(X_tr[:, sel], y_tr)
-            if grid.best_score_ > best_cv:
-                best_cv = grid.best_score_
-                best_n = nfeat
-                best_g = grid.best_params_["svm__gamma"]
+        sorted_idx = self._rf_ranking(X_tr, y_tr)
+        best_cv, best_n, best_g = self._grid_search(
+            X_tr, y_tr, g_tr, cv, sorted_idx, nfeatures_list, gamma_logspace)
 
         # --- Retrain + test ---
         final = Pipeline([("scaler", StandardScaler()),
@@ -140,7 +100,7 @@ class PatientLeakageClassifier:
         final.fit(X_tr[:, sorted_idx[:best_n]], y_tr)
         y_pred = final.predict(X_te[:, sorted_idx[:best_n]])
 
-        return {
+        result = {
             "model": model_name,
             "seed": seed,
             "split_type": split_type,
@@ -153,11 +113,107 @@ class PatientLeakageClassifier:
             "n_train_patients": n_tr_patients,
             "n_test_patients": n_te_patients,
         }
+        if return_predictions:
+            result["y_true"] = [str(v) for v in y_te]
+            result["y_pred"] = [str(v) for v in y_pred]
+        return result
+
+    def cv_curve(self, model_name, seed, split_type="image",
+                 nfeatures_list=None, gamma_logspace=None):
+        """CV accuracy as a function of the number of RF-selected features.
+
+        Same split → RF ranking → per-nfeatures γ-grid as ``run``, but returns
+        the whole curve (used by the Figure 3 / S1 / S4 CV-vs-features plots)
+        instead of only the best point.
+
+        Returns
+        -------
+        dict: {"curves": {nfeat: cv_acc}, "best_nfeat", "best_cv", "best_gamma"}
+        """
+        if nfeatures_list is None:
+            nfeatures_list = DEFAULT_NFEATURES
+        if gamma_logspace is None:
+            gamma_logspace = DEFAULT_GAMMA_RANGE
+
+        X, Y, groups = self._load_data(model_name)
+        X_tr, _, y_tr, _, g_tr, cv, _, _ = self._split(
+            X, Y, groups, seed, split_type)
+        sorted_idx = self._rf_ranking(X_tr, y_tr)
+
+        pipe = Pipeline([("scaler", StandardScaler()),
+                         ("svm", SVC(kernel="rbf", C=10))])
+        valid_nfeats = [n for n in nfeatures_list if n <= X_tr.shape[1]]
+        curves, best_cv, best_n, best_g = {}, 0, None, None
+        for nfeat in valid_nfeats:
+            pg = {"svm__gamma": (1.0 / nfeat) * np.logspace(*gamma_logspace)}
+            grid = GridSearchCV(pipe, pg, cv=cv, scoring="accuracy",
+                                n_jobs=self.n_jobs, verbose=0)
+            if g_tr is not None:
+                grid.fit(X_tr[:, sorted_idx[:nfeat]], y_tr, groups=g_tr)
+            else:
+                grid.fit(X_tr[:, sorted_idx[:nfeat]], y_tr)
+            curves[int(nfeat)] = float(grid.best_score_)
+            if grid.best_score_ > best_cv:
+                best_cv = float(grid.best_score_)
+                best_n = int(nfeat)
+                best_g = float(grid.best_params_["svm__gamma"])
+        return {"curves": curves, "best_nfeat": best_n,
+                "best_cv": best_cv, "best_gamma": best_g}
 
     # ------------------------------------------------------------------
     # Internals
     # ------------------------------------------------------------------
 
+    def _split(self, X, Y, groups, seed, split_type):
+        """Image-level or patient-level 80:20 train/test split.
+
+        Returns (X_tr, X_te, y_tr, y_te, g_tr, cv, n_tr_patients, n_te_patients);
+        g_tr is None and cv is 5 for image-level splits.
+        """
+        if split_type == "image":
+            X_tr, X_te, y_tr, y_te = train_test_split(
+                X, Y, test_size=0.2, random_state=seed, stratify=Y)
+            return X_tr, X_te, y_tr, y_te, None, 5, None, None
+
+        unique_g = np.unique(groups)
+        g_labels = np.array([Y[groups == g][0] for g in unique_g])
+        tr_g, te_g = train_test_split(
+            unique_g, test_size=0.2, random_state=seed, stratify=g_labels)
+        tr_mask = np.isin(groups, tr_g)
+        te_mask = np.isin(groups, te_g)
+        cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=seed)
+        return (X[tr_mask], X[te_mask], Y[tr_mask], Y[te_mask],
+                groups[tr_mask], cv, len(np.unique(groups[tr_mask])),
+                len(np.unique(groups[te_mask])))
+
+    def _rf_ranking(self, X_tr, y_tr):
+        """Feature indices ranked by RF Gini importance (descending)."""
+        fs = RandomForestClassifier(n_estimators=500, random_state=15,
+                                    class_weight="balanced", n_jobs=self.n_jobs)
+        fs.fit(X_tr, y_tr)
+        return np.argsort(fs.feature_importances_)[::-1]
+
+    def _grid_search(self, X_tr, y_tr, g_tr, cv, sorted_idx,
+                     nfeatures_list, gamma_logspace):
+        """Search nfeatures × gamma; return (best_cv, best_nfeat, best_gamma)."""
+        pipe = Pipeline([("scaler", StandardScaler()),
+                         ("svm", SVC(kernel="rbf", C=10))])
+        valid_nfeats = [n for n in nfeatures_list if n <= X_tr.shape[1]]
+        best_cv, best_n, best_g = 0, None, None
+        for nfeat in valid_nfeats:
+            pg = {"svm__gamma": (1.0 / nfeat) * np.logspace(*gamma_logspace)}
+            grid = GridSearchCV(pipe, pg, cv=cv, scoring="accuracy",
+                                n_jobs=self.n_jobs, verbose=0)
+            if g_tr is not None:
+                grid.fit(X_tr[:, sorted_idx[:nfeat]], y_tr, groups=g_tr)
+            else:
+                grid.fit(X_tr[:, sorted_idx[:nfeat]], y_tr)
+            if grid.best_score_ > best_cv:
+                best_cv = grid.best_score_
+                best_n = nfeat
+                best_g = grid.best_params_["svm__gamma"]
+        return best_cv, best_n, best_g
+
     def _load_data(self, model_name):
         data = np.load(os.path.join(self.features_dir,
                                      f"{model_name}_features.npz"),

classes/features.py

@@ -96,8 +96,11 @@ class FeatureExtractor:
             self.model = full_model.features
 
         elif name == "EfficientNetB1":
-            # Keep conv stack, drop avgpool + classifier
-            # → output (B, 1280, 7, 7)  (at 224 px; 8×8 at 240 px)
+            # Keep conv stack, drop avgpool + classifier.
+            # Input 240×240 → 8×8 spatial (240/32 = 7.5 → 8 with padding).
+            # Output: (B, 1280, 8, 8) → 81,920-d.
+            # NOTE: manuscript incorrectly reports 62,720 (7×7×1280);
+            # they used the 224px spatial size for the calculation.
             self.model = full_model.features
 
         elif name == "MobileNetV2":
@@ -106,10 +109,16 @@ class FeatureExtractor:
             self.model = full_model.features
 
         elif name == "ResNet50":
-            # Drop fc, keep everything *including* the adaptive avg pool
-            # → output (B, 2048, 1, 1)  — matches TF include_top=False
-            full_model.fc = nn.Identity()
-            self.model = full_model
+            # Drop BOTH avgpool and FC → (B, 2048, 7, 7) → 100,352-d.
+            # This matches TF ResNet50(include_top=False) with no pooling and
+            # the manuscript's reported dimension, and keeps ResNet50 consistent
+            # with the other four models (all use flattened spatial conv maps).
+            # Empirically the spatial features give ~0.99 image-level accuracy,
+            # vs ~0.97 for the GAP-pooled 2,048-d vector, which discards the
+            # 7×7 grid the RF+SVM pipeline relies on.
+            self.model = nn.Sequential(
+                *list(full_model.children())[:-2]
+            )
 
         self.model.to(self.device)
         self.model.eval()
@@ -162,25 +171,17 @@ class FeatureExtractor:
         features = []
         for start in range(0, n_images, batch_size):
             end = min(start + batch_size, n_images)
-            batch_paths = image_paths[start:end]
             batch_tensors = []
 
-            for path in batch_paths:
+            for path in image_paths[start:end]:
                 try:
                     img = Image.open(path).convert("RGB")
-                    tensor = self.transform(img)
-                    batch_tensors.append(tensor)
+                    batch_tensors.append(self.transform(img))
                 except Exception as e:
                     print(f"  Error loading {path}: {e}")
                     batch_tensors.append(torch.zeros(3, self.input_size, self.input_size))
 
-            batch = torch.stack(batch_tensors).to(self.device)
-
-            with torch.no_grad():
-                batch_features = self.model(batch)
-                batch_features = batch_features.view(batch_features.size(0), -1)
-
-            features.append(batch_features.cpu().numpy())
+            features.append(self._forward(batch_tensors))
 
             if (start // batch_size) % 20 == 0:
                 print(f"  Processed {end}/{n_images} images...")
@@ -192,6 +193,74 @@ class FeatureExtractor:
         print(f"  Feature matrix shape: {X.shape}")
         return X, Y, filenames
 
+    def extract_from_images(self, images, labels=None, filenames=None,
+                            batch_size=32):
+        """Extract features from in-memory images (dataset-layout agnostic).
+
+        Unlike ``extract``, this does not assume a directory tree of JPEGs
+        grouped into class folders, so it works for any source that can
+        produce image arrays (e.g. axial slices decoded from NIfTI volumes).
+
+        Parameters
+        ----------
+        images : sequence
+            Each item may be a PIL.Image, an (H, W) grayscale array, or an
+            (H, W, 3) RGB array (uint8 or float). Grayscale is expanded to
+            RGB and everything is resized/normalised via ``self.transform``.
+        labels, filenames : sequence or None
+            Optional per-image class labels / names, returned unchanged so the
+            output matches ``extract``'s (X, Y, filenames) contract.
+
+        Returns
+        -------
+        X : np.ndarray  (n_images, n_features)
+        Y : np.ndarray  (n_images,)
+        filenames : np.ndarray  (n_images,)
+        """
+        n_images = len(images)
+        print(f"  Extracting features from {n_images} in-memory images "
+              f"({self.model_name})")
+
+        features = []
+        for start in range(0, n_images, batch_size):
+            end = min(start + batch_size, n_images)
+            batch_tensors = [self.transform(self._to_pil(im))
+                             for im in images[start:end]]
+            features.append(self._forward(batch_tensors))
+
+            if (start // batch_size) % 20 == 0:
+                print(f"  Processed {end}/{n_images} images...")
+
+        X = np.concatenate(features, axis=0).astype(np.float32)
+        Y = np.array(labels) if labels is not None else np.array([None] * n_images)
+        fnames = (np.array(filenames) if filenames is not None
+                  else np.arange(n_images))
+        print(f"  Feature matrix shape: {X.shape}")
+        return X, Y, fnames
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    def _forward(self, batch_tensors):
+        """Run one batch of preprocessed tensors through the model."""
+        batch = torch.stack(batch_tensors).to(self.device)
+        with torch.no_grad():
+            f = self.model(batch)
+            f = f.view(f.size(0), -1)
+        return f.cpu().numpy()
+
+    @staticmethod
+    def _to_pil(im):
+        """Coerce a PIL image or numpy array (gray/RGB) to an RGB PIL image."""
+        if isinstance(im, Image.Image):
+            return im.convert("RGB")
+        arr = np.asarray(im)
+        if arr.dtype != np.uint8:
+            lo, hi = float(arr.min()), float(arr.max())
+            arr = (255.0 * (arr - lo) / (hi - lo + 1e-8)).astype(np.uint8)
+        return Image.fromarray(arr).convert("RGB")
+
     def save_features(self, X, Y, filenames, output_dir):
         """Save extracted features to a compressed .npz file."""
         os.makedirs(output_dir, exist_ok=True)

classes/nifti_dataset.py

@@ -0,0 +1,247 @@
+"""
+Load axial slices from a directory of NIfTI CT volumes.
+
+Each ``.nii.gz`` file is treated as one patient (the filename stem is the
+ground-truth patient ID). A handful of evenly-spaced axial slices are sampled
+per volume, HU-windowed, and returned both as grayscale thumbnails (for the
+manuscript's clustering method) and as full-resolution slices (for CNN feature
+extraction).
+
+This is used to validate the patient-clustering method on a dataset where the
+true patient identity IS known — e.g. the Medical Segmentation Decathlon
+Task06_Lung set — unlike IQ-OTH/NCCD, where patient IDs are unavailable.
+"""
+
+import os
+import glob
+import numpy as np
+
+
+class NiftiSliceDataset:
+    """Sample HU-windowed axial slices from a folder of NIfTI CT volumes.
+
+    Parameters
+    ----------
+    volume_dir : str
+        Directory containing ``*.nii.gz`` volumes, one per patient.
+    n_slices_per_patient : int
+        Number of evenly-spaced axial slices to sample per volume.
+    thumbnail_size : tuple of int
+        Grayscale thumbnail size for clustering (matches the manuscript's 64x64).
+    hu_window : tuple of float
+        (low, high) Hounsfield-unit clip range applied before scaling to 8-bit.
+        The default (-1000, 400) keeps the body outline, soft tissue, and lung
+        parenchyma — the cues that identify a patient.
+    central_fraction : float
+        Fraction of the volume (centred on the middle slice) to sample from,
+        avoiding empty/partial slices at the superior/inferior ends.
+    slice_axis : int
+        Axis along which axial slices are indexed (-1 = last, correct for the
+        (H, W, Z) layout used by the Decathlon volumes).
+    """
+
+    def __init__(self, volume_dir, n_slices_per_patient=10,
+                 thumbnail_size=(64, 64), hu_window=(-1000.0, 400.0),
+                 central_fraction=0.6, slice_axis=-1, random_slices=True,
+                 seed=42, flip_vertical=False, rotate_deg=0):
+        self.volume_dir = volume_dir
+        self.n_slices_per_patient = n_slices_per_patient
+        self.thumbnail_size = thumbnail_size
+        self.hu_window = hu_window
+        self.central_fraction = central_fraction
+        self.slice_axis = slice_axis
+        self.random_slices = random_slices
+        self.seed = seed
+        self.flip_vertical = flip_vertical
+        self.rotate_deg = rotate_deg  # 0, 90, 180, or 270
+
+        # Populated by load()
+        self.thumbnails = None      # (n_slices, H*W) float32 in [0, 1]
+        self.slices_rgb = None      # list of (H, W) uint8 arrays
+        self.patient_labels = None  # (n_slices,) ground-truth patient IDs
+        self.filenames = None       # (n_slices,) per-slice names
+
+    def load(self, verbose=True):
+        """Load and slice every volume in ``volume_dir``."""
+        import nibabel as nib
+        from PIL import Image
+
+        # Prefer decompressed .nii (fast memmap random-slice access) and fall
+        # back to .nii.gz, skipping any .gz that has an uncompressed twin.
+        nii = glob.glob(os.path.join(self.volume_dir, "*.nii"))
+        stems = {os.path.basename(p)[:-4] for p in nii}
+        gz = [p for p in glob.glob(os.path.join(self.volume_dir, "*.nii.gz"))
+              if os.path.basename(p)[:-7] not in stems]
+        paths = sorted(nii + gz)
+        if not paths:
+            raise FileNotFoundError(
+                f"No .nii/.nii.gz volumes found in {self.volume_dir}")
+
+        thumbs, slices_rgb, patient_labels, filenames = [], [], [], []
+        lo, hi = self.hu_window
+
+        for path in paths:
+            patient_id = os.path.basename(path).replace(".nii.gz", "").replace(".nii", "")
+            img = nib.load(path)
+            n_z = img.shape[self.slice_axis]
+            axis = self.slice_axis % img.ndim
+
+            for k in self._slice_indices(n_z):
+                # Lazy slice: pull only this plane off disk instead of
+                # materialising the whole float64 volume with get_fdata().
+                slicer = [slice(None)] * img.ndim
+                slicer[axis] = k
+                sl = np.asarray(img.dataobj[tuple(slicer)], dtype=np.float32)
+                # HU window -> 8-bit grayscale
+                sl = np.clip(sl, lo, hi)
+                sl = (255.0 * (sl - lo) / (hi - lo)).astype(np.uint8)
+                if self.flip_vertical:
+                    sl = np.flipud(sl)
+                if self.rotate_deg:
+                    sl = np.rot90(sl, k=self.rotate_deg // 90)
+
+                thumb = Image.fromarray(sl).convert("L").resize(
+                    self.thumbnail_size)
+                thumbs.append(np.asarray(thumb, dtype=np.float32).flatten() / 255.0)
+                slices_rgb.append(sl)
+                patient_labels.append(patient_id)
+                filenames.append(f"{patient_id}_slice{k:03d}")
+
+            if verbose:
+                print(f"  {patient_id}: {n_z} slices -> "
+                      f"sampled {self.n_slices_per_patient}")
+
+        self.thumbnails = np.array(thumbs, dtype=np.float32)
+        self.slices_rgb = slices_rgb
+        self.patient_labels = np.array(patient_labels)
+        self.filenames = np.array(filenames)
+
+        if verbose:
+            print(f"\n  Loaded {len(self.filenames)} slices from "
+                  f"{len(paths)} patients")
+        return self
+
+    def load_all_slices(self, stride=3, verbose=True):
+        """Load ALL slices from the central fraction (not just a sample).
+
+        Every ``stride``-th slice is taken to keep feature extraction manageable
+        while still giving K-means enough data to cluster per patient.
+
+        Returns self (populates thumbnails, slices_rgb, patient_labels, filenames).
+        """
+        import nibabel as nib
+        from PIL import Image
+
+        nii = glob.glob(os.path.join(self.volume_dir, "*.nii"))
+        stems = {os.path.basename(p)[:-4] for p in nii}
+        gz = [p for p in glob.glob(os.path.join(self.volume_dir, "*.nii.gz"))
+              if os.path.basename(p)[:-7] not in stems]
+        paths = sorted(nii + gz)
+        if not paths:
+            raise FileNotFoundError(
+                f"No .nii/.nii.gz volumes found in {self.volume_dir}")
+
+        thumbs, slices_rgb, patient_labels, filenames = [], [], [], []
+        lo, hi = self.hu_window
+        total_slices = 0
+
+        for path in paths:
+            patient_id = os.path.basename(path).replace(".nii.gz", "").replace(".nii", "")
+            img = nib.load(path)
+            n_z = img.shape[self.slice_axis]
+            axis = self.slice_axis % img.ndim
+
+            # All slices in central fraction, stepped
+            half = self.central_fraction / 2.0
+            start = int(n_z * (0.5 - half))
+            end = int(n_z * (0.5 + half))
+            indices = list(range(start, end, stride))
+
+            for k in indices:
+                slicer = [slice(None)] * img.ndim
+                slicer[axis] = k
+                sl = np.asarray(img.dataobj[tuple(slicer)], dtype=np.float32)
+                sl = np.clip(sl, lo, hi)
+                sl = (255.0 * (sl - lo) / (hi - lo)).astype(np.uint8)
+                if self.flip_vertical:
+                    sl = np.flipud(sl)
+                if self.rotate_deg:
+                    sl = np.rot90(sl, k=self.rotate_deg // 90)
+
+                thumb = Image.fromarray(sl).convert("L").resize(
+                    self.thumbnail_size)
+                thumbs.append(np.asarray(thumb, dtype=np.float32).flatten() / 255.0)
+                slices_rgb.append(sl)
+                patient_labels.append(patient_id)
+                filenames.append(f"{patient_id}_slice{k:03d}")
+                total_slices += 1
+
+            if verbose:
+                print(f"  {patient_id}: {n_z} slices -> {len(indices)} sampled")
+
+        self.thumbnails = np.array(thumbs, dtype=np.float32)
+        self.slices_rgb = slices_rgb
+        self.patient_labels = np.array(patient_labels)
+        self.filenames = np.array(filenames)
+
+        if verbose:
+            print(f"\n  Loaded {total_slices} slices from "
+                  f"{len(paths)} patients (stride={stride})")
+        return self
+
+    def export_pngs(self, cache_dir):
+        """Export loaded slices as PNG files for fast reloading.
+
+        Creates ``cache_dir/`` with one PNG per slice and a manifest.json
+        mapping paths → patient IDs and z-indices.
+
+        Returns (paths, patient_ids, z_indices).
+        """
+        import json
+        from PIL import Image as PILImage
+
+        os.makedirs(cache_dir, exist_ok=True)
+        manifest_path = os.path.join(cache_dir, "manifest.json")
+
+        if os.path.exists(manifest_path):
+            print(f"  Loading cached PNGs from {cache_dir}")
+            with open(manifest_path) as f:
+                manifest = json.load(f)
+            return (manifest["paths"], np.array(manifest["patient_ids"]),
+                    np.array(manifest["z_indices"]))
+
+        if self.slices_rgb is None or self.patient_labels is None:
+            raise RuntimeError("Call load() or load_all_slices() first.")
+
+        print(f"  Exporting {len(self.slices_rgb)} PNGs to {cache_dir} ...")
+        paths, pids, zs = [], [], []
+        for i, (sl, pid, fname) in enumerate(
+            zip(self.slices_rgb, self.patient_labels, self.filenames)):
+            out_path = os.path.join(cache_dir, f"{fname}.png")
+            PILImage.fromarray(sl).save(out_path)
+            paths.append(out_path)
+            pids.append(pid)
+            zs.append(int(fname.split("_slice")[-1]))
+            if (i + 1) % 500 == 0:
+                print(f"    {i + 1}/{len(self.slices_rgb)} ...", flush=True)
+
+        manifest = {"paths": paths, "patient_ids": pids, "z_indices": zs}
+        with open(manifest_path, "w") as f:
+            json.dump(manifest, f)
+        print(f"  Exported {len(paths)} PNGs")
+        return paths, np.array(pids), np.array(zs)
+
+    @property
+    def n_patients(self):
+        return len(np.unique(self.patient_labels))
+
+    def _slice_indices(self, n_z):
+        """Slice indices within the central fraction (even or random)."""
+        half = self.central_fraction / 2.0
+        start = int(n_z * (0.5 - half))
+        end = int(n_z * (0.5 + half))
+        if self.random_slices:
+            rng = np.random.default_rng(self.seed)
+            return np.sort(rng.integers(start, end, size=self.n_slices_per_patient))
+        else:
+            return np.linspace(start, end - 1, self.n_slices_per_patient).astype(int)

classes/patient_identifier.py

@@ -53,8 +53,8 @@ class PatientIdentifier:
         self.labels_ = labels
         groups = self._run_clustering(
             labels, filenames,
-            data_loader_fn=lambda class_name, files: self._load_thumbnails(
-                base_path, class_name, files),
+            data_loader_fn=lambda class_name, idx: self._load_thumbnails(
+                base_path, class_name, filenames[idx]),
             valid_extensions=valid_extensions,
         )
         self.groups_ = groups
@@ -82,21 +82,28 @@ class PatientIdentifier:
         self.groups_ = groups
         return groups
 
-    def build_assignment_dict(self, groups, filenames, labels):
+    def build_assignment_dict(self, groups, filenames, labels,
+                              class_prefixes=None):
         """Convert flat group array into a human-readable dict.
 
+        Parameters
+        ----------
+        class_prefixes : dict or None
+            Optional mapping from class_name -> short prefix used in the
+            patient IDs (e.g. {"Bengin cases": "benign"}). If None, a prefix
+            is derived automatically from each class name, so the method works
+            for any dataset rather than only IQ-OTH's three classes.
+
         Returns
         -------
         dict mapping patient_id (e.g. "benign_0") -> list of image filenames.
         """
         self.labels_ = labels
-        class_map = {
-            "Bengin cases":    "benign",
-            "Malignant cases": "malig",
-            "Normal cases":    "normal",
-        }
+        if class_prefixes is None:
+            class_prefixes = {c: self._slugify(c)
+                              for c in self.patient_estimates}
         assignment = {}
-        for class_name, prefix in class_map.items():
+        for class_name, prefix in class_prefixes.items():
             idx_class = np.where(labels == class_name)[0]
             # Get unique cluster IDs within this class (preserve order)
             seen = []
@@ -115,6 +122,12 @@ class PatientIdentifier:
         self.assignments_ = assignment
         return assignment
 
+    @staticmethod
+    def _slugify(class_name):
+        """Turn a class name into a short lowercase prefix for patient IDs."""
+        token = str(class_name).strip().lower().split()[0]
+        return "".join(ch for ch in token if ch.isalnum()) or "group"
+
     def get_unique_groups(self):
         if self.groups_ is None:
             raise RuntimeError("Call identify_*() first.")

classes/siamese.py

@@ -0,0 +1,599 @@
+"""
+SiamesePatientMatcher — learned patient identification via connected components.
+
+Trains a siamese CNN on Task06_Lung (known patient IDs), then applies the
+trained model to IQ-OTH/NCCD to build a patient manifest without K-means.
+
+Architecture:
+  Two CT slices → shared backbone → [f_A, f_B, |f_A-f_B|] → MLP head → same/different
+
+For inference, we build a graph: edge between slices i,j if the siamese
+confidence exceeds a threshold, then find connected components.
+Each component = one estimated patient.
+"""
+
+import os
+import numpy as np
+from collections import defaultdict
+
+import torch
+import torch.nn as nn
+from torchvision import transforms, models
+
+
+# ---------------------------------------------------------------------------
+# Backbone registry
+# ---------------------------------------------------------------------------
+
+BACKBONES = {
+    "resnet18": (models.resnet18, models.ResNet18_Weights.IMAGENET1K_V1, 512),
+    "resnet34": (models.resnet34, models.ResNet34_Weights.IMAGENET1K_V1, 512),
+    "efficientnet_b0": (models.efficientnet_b0, models.EfficientNet_B0_Weights.IMAGENET1K_V1, 1280),
+}
+
+
+# ---------------------------------------------------------------------------
+# Model definition
+# ---------------------------------------------------------------------------
+
+
+class SiameseCNN(nn.Module):
+    """Shared CNN backbone → concatenate → MLP head → binary classification."""
+
+    def __init__(self, backbone_name="resnet18", hidden_dims=None):
+        super().__init__()
+        if hidden_dims is None:
+            hidden_dims = [512, 128]
+
+        fn, weights, feat_dim = BACKBONES[backbone_name]
+        cnn = fn(weights=weights)
+
+        if backbone_name.startswith("resnet"):
+            self.backbone = nn.Sequential(
+                cnn.conv1, cnn.bn1, cnn.relu, cnn.maxpool,
+                cnn.layer1, cnn.layer2, cnn.layer3, cnn.layer4,
+                nn.AdaptiveAvgPool2d((1, 1)), nn.Flatten())
+        elif backbone_name.startswith("efficientnet"):
+            self.backbone = nn.Sequential(
+                cnn.features, nn.AdaptiveAvgPool2d((1, 1)), nn.Flatten())
+        else:
+            raise ValueError(f"Unknown backbone: {backbone_name}")
+
+        # Head: input = backbone_dim * 3
+        head_input = feat_dim * 3
+        layers = []
+        prev = head_input
+        for h in hidden_dims:
+            layers.append(nn.Linear(prev, h))
+            layers.append(nn.BatchNorm1d(h))
+            layers.append(nn.ReLU())
+            layers.append(nn.Dropout(0.3))
+            prev = h
+        layers.append(nn.Linear(prev, 1))
+        self.head = nn.Sequential(*layers)
+
+    def forward(self, img_a, img_b):
+        fa = self.backbone(img_a)
+        fb = self.backbone(img_b)
+        combined = torch.cat([fa, fb, torch.abs(fa - fb)], dim=-1)
+        return self.head(combined).squeeze(-1)
+
+    def embed(self, images, device):
+        """Extract backbone feature vectors for a batch of images."""
+        return self.backbone(images)
+
+
+# ---------------------------------------------------------------------------
+# Patient matcher
+# ---------------------------------------------------------------------------
+
+
+class SiamesePatientMatcher:
+    """Apply a trained siamese model to identify patient groups in new data.
+
+    Parameters
+    ----------
+    model_path : str
+        Path to saved model weights (.pt file).
+    backbone : str
+        Backbone name matching the saved model.
+    device : str or None
+        Torch device. Auto-detected if None.
+    input_size : int
+        Input resolution for the backbone (224 for ResNet, 240 for EfficientNet).
+    """
+
+    def __init__(self, model_path, backbone="resnet18", device=None,
+                 input_size=224):
+        self.device = torch.device(
+            device or ("cuda" if torch.cuda.is_available() else "cpu"))
+        self.input_size = input_size
+        self.model = SiameseCNN(backbone).to(self.device)
+        self.model.load_state_dict(
+            torch.load(model_path, map_location=self.device, weights_only=True))
+        self.model.eval()
+
+        self.transform = transforms.Compose([
+            transforms.Resize((input_size, input_size)),
+            transforms.ToTensor(),
+            transforms.Normalize(mean=[0.485, 0.456, 0.406],
+                               std=[0.229, 0.224, 0.225]),
+        ])
+
+    def identify_patients(self, images, filenames=None, threshold=0.9,
+                          top_k=20, batch_size=64, k=None,
+                          cluster_method="edge_rank",
+                          min_size=None, max_size=None, keep_k=False):
+        """Build a patient manifest.
+
+        If ``k`` is None (default): connected-components clustering with a
+        hard similarity threshold. No prior knowledge of patient count needed.
+
+        If ``k`` is provided: spectral clustering into exactly ``k`` groups
+        using the siamese similarity graph. Weak/spurious edges get cut to
+        respect the known patient count. Use this when you know K a priori
+        (e.g. IQ-OTH has 15 Benign, 40 Malignant, 55 Normal patients).
+
+        Parameters
+        ----------
+        images : list
+            PIL Images, numpy arrays, or paths to image files.
+        filenames : list of str or None
+            Image filenames for the output manifest.
+        threshold : float
+            Min siamese probability for an edge (connected-components mode only).
+        top_k : int
+            Top-K similar candidates to verify per slice.
+        batch_size : int
+            Batch size for embedding extraction.
+        k : int or None
+            If provided, partition into exactly k groups via spectral clustering.
+
+        Returns
+        -------
+        dict mapping patient_id → list of filenames.
+        """
+        if filenames is None:
+            filenames = [str(i) for i in range(len(images))]
+
+        # ---- Load and embed all images ----
+        print(f"Embedding {len(images)} images ...", flush=True)
+        all_embeddings = self.embed_images(images, batch_size)
+        print(f"  Embeddings: {all_embeddings.shape}", flush=True)
+
+        # ---- Build similarity graph ----
+        sim_matrix, verified_edges = self._build_similarity_graph(
+            all_embeddings, threshold, top_k)
+
+        # ---- Cluster ----
+        if k is not None:
+            if cluster_method in ("complete", "average"):
+                # Agglomerative clustering on the dense siamese P(same) matrix.
+                # Complete/average linkage resist the single-linkage "chaining"
+                # that edge-ranking/connected-components suffer when the model
+                # is over-confident (e.g. out-of-distribution on IQ-OTH), where
+                # a few cross-patient edges merge many patients into one blob.
+                manifest = self._agglomerative_cluster(
+                    all_embeddings, filenames, k, linkage=cluster_method)
+            else:
+                # Default: edge-ranking (preserves natural clusters), falling
+                # back to spectral if not enough edges to reach k.
+                manifest = self._edge_rank_cluster(
+                    verified_edges, len(filenames), filenames, k)
+                if manifest is None:
+                    print("  Edge-ranking couldn't reach k, falling back to "
+                          "spectral", flush=True)
+                    manifest = self._spectral_cluster(
+                        sim_matrix, all_embeddings, filenames, k)
+        else:
+            manifest = self._connected_components(
+                verified_edges, len(filenames), filenames)
+
+        if min_size or max_size:
+            if keep_k and k is not None:
+                manifest = self._rebalance_keep_k(
+                    manifest, filenames, all_embeddings, k, min_size, max_size)
+            else:
+                manifest = self._rebalance(
+                    manifest, filenames, all_embeddings, min_size, max_size)
+        return manifest
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    def _build_similarity_graph(self, embeddings, threshold, top_k):
+        """Build a sparse similarity graph verified by the siamese head.
+
+        Returns
+        -------
+        sim_matrix : (n, n) ndarray  cosine similarity (all pairs)
+        edges : list of (i, j, prob)  verified edges above threshold
+        """
+        n = len(embeddings)
+
+        # Cosine similarity
+        emb_norm = embeddings / (np.linalg.norm(embeddings, axis=1,
+                                                 keepdims=True) + 1e-8)
+        sim_matrix = emb_norm @ emb_norm.T
+        top_k_idx = np.argsort(-sim_matrix, axis=1)[:, 1:top_k + 1]
+
+        # Verify top-K candidates with siamese head
+        print(f"Verifying top-{top_k} candidates per slice "
+              f"(threshold={threshold}) ...", flush=True)
+        edges = []
+        emb_t = torch.from_numpy(embeddings.astype(np.float32)).to(self.device)
+
+        for i in range(n):
+            for j in top_k_idx[i]:
+                if j <= i:
+                    continue
+                with torch.no_grad():
+                    fa = emb_t[i:i + 1]
+                    fb = emb_t[j:j + 1]
+                    combined = torch.cat([fa, fb, torch.abs(fa - fb)], dim=-1)
+                    prob = torch.sigmoid(self.model.head(combined)).item()
+                    if prob > threshold:
+                        edges.append((i, j, prob))
+            if (i + 1) % 500 == 0:
+                print(f"  {i + 1}/{n} slices, {len(edges)} edges",
+                      flush=True)
+
+        print(f"  Total verified edges: {len(edges)}", flush=True)
+        return sim_matrix, edges
+
+    def _connected_components(self, edges, n, filenames):
+        """Cluster via connected components on verified edges."""
+        from scipy.sparse.csgraph import connected_components
+        from scipy.sparse import csr_matrix
+
+        if edges:
+            row, col = zip(*[(e[0], e[1]) for e in edges])
+            # Symmetric: each edge goes both ways, double the data
+            all_row = list(row) + list(col)
+            all_col = list(col) + list(row)
+            data = np.ones(len(all_row))
+            adj = csr_matrix((data, (all_row, all_col)), shape=(n, n))
+            n_components, labels = connected_components(adj, directed=False)
+        else:
+            n_components, labels = n, np.arange(n)
+
+        return self._labels_to_manifest(labels, n_components, filenames,
+                                        "siamese")
+
+    def embed_images(self, images, batch_size=64):
+        """Backbone feature vectors for paths / PIL images / arrays."""
+        from PIL import Image as PILImage
+        embs = []
+        for start in range(0, len(images), batch_size):
+            batch = []
+            for item in images[start:start + batch_size]:
+                if isinstance(item, str):
+                    img = PILImage.open(item).convert("RGB")
+                elif isinstance(item, PILImage.Image):
+                    img = item.convert("RGB")
+                else:
+                    arr = np.asarray(item)
+                    if arr.dtype != np.uint8:
+                        arr = (255.0 * (arr - arr.min()) /
+                               (arr.max() - arr.min() + 1e-8)).astype(np.uint8)
+                    img = PILImage.fromarray(arr).convert("RGB")
+                batch.append(self.transform(img))
+            batch_t = torch.stack(batch).to(self.device)
+            with torch.no_grad():
+                embs.append(self.model.backbone(batch_t).cpu().numpy())
+            if (start // batch_size) % 20 == 0:
+                print(f"  {start + len(batch)}/{len(images)}", flush=True)
+        return np.concatenate(embs, axis=0)
+
+    def _dense_prob_matrix(self, embeddings):
+        """Full symmetric P(same-patient) matrix from the siamese head.
+
+        Scores every pair with the head (cat[fa, fb, |fa-fb|] -> sigmoid),
+        one query row at a time to avoid materialising all n^2 vectors, then
+        symmetrises since the head is not exactly order-invariant.
+        """
+        n = len(embeddings)
+        emb = torch.from_numpy(embeddings.astype(np.float32)).to(self.device)
+        P = np.zeros((n, n), dtype=np.float32)
+        with torch.no_grad():
+            for i in range(n):
+                fa = emb[i:i + 1].expand(n, -1)
+                combined = torch.cat([fa, emb, torch.abs(fa - emb)], dim=-1)
+                P[i] = torch.sigmoid(
+                    self.model.head(combined).squeeze(-1)).cpu().numpy()
+        return 0.5 * (P + P.T)
+
+    def _agglomerative_cluster(self, embeddings, filenames, k, linkage="complete"):
+        """Cluster into k groups via agglomerative clustering on 1 - P(same).
+
+        Unlike single-linkage (edge-ranking / connected components), complete
+        and average linkage will not merge two groups on the strength of a
+        single confident cross-patient edge, so they resist chaining.
+        """
+        from sklearn.cluster import AgglomerativeClustering
+
+        print(f"  Agglomerative clustering ({linkage} linkage) into k={k} ...",
+              flush=True)
+        P = self._dense_prob_matrix(embeddings)
+        dist = 1.0 - P
+        np.fill_diagonal(dist, 0.0)
+        dist[dist < 0] = 0.0
+        try:
+            model = AgglomerativeClustering(
+                n_clusters=k, metric="precomputed", linkage=linkage)
+        except TypeError:  # scikit-learn < 1.2
+            model = AgglomerativeClustering(
+                n_clusters=k, affinity="precomputed", linkage=linkage)
+        labels = model.fit_predict(dist)
+        return self._labels_to_manifest(labels, k, filenames, "siamese")
+
+    def _rebalance(self, manifest, filenames, embeddings,
+                   min_size=None, max_size=None):
+        """Enforce group-size bounds using the siamese distances.
+
+        Oversized groups (> max_size) are split at their natural gaps via
+        complete-linkage on the members' 1 - P(same) submatrix — appropriate
+        when a group is really several distinct patients chained together.
+        Undersized groups (< min_size, e.g. orphaned singletons that are a
+        leakage hazard) are absorbed into their nearest group. Splitting first,
+        then absorbing, keeps sizes within bounds where possible.
+        """
+        from sklearn.cluster import AgglomerativeClustering
+
+        f2i = {f: i for i, f in enumerate(filenames)}
+        dist = 1.0 - self._dense_prob_matrix(embeddings)
+        np.fill_diagonal(dist, 0.0)
+        dist[dist < 0] = 0.0
+        groups = {pid: [f2i[f] for f in fs if f in f2i]
+                  for pid, fs in manifest.items()}
+
+        # 1) Split oversized groups. Complete-linkage into ceil(size/max_size)
+        #    sub-groups can still leave one child over the cap (uneven splits),
+        #    so repeat until every group is <= max_size.
+        if max_size:
+            changed = True
+            while changed:
+                changed = False
+                for pid in list(groups):
+                    idx = groups[pid]
+                    if len(idx) <= max_size:
+                        continue
+                    n_sub = int(np.ceil(len(idx) / max_size))
+                    sub = dist[np.ix_(idx, idx)]
+                    try:
+                        model = AgglomerativeClustering(
+                            n_clusters=n_sub, metric="precomputed", linkage="complete")
+                    except TypeError:
+                        model = AgglomerativeClustering(
+                            n_clusters=n_sub, affinity="precomputed", linkage="complete")
+                    lab = model.fit_predict(sub)
+                    del groups[pid]
+                    for s in range(n_sub):
+                        members = [idx[j] for j in range(len(idx)) if lab[j] == s]
+                        if members:
+                            groups[f"{pid}_s{s}"] = members
+                    changed = True
+
+        # 2) Absorb undersized groups into their nearest remaining group.
+        if min_size:
+            changed = True
+            while changed:
+                changed = False
+                for pid in sorted((p for p in groups if len(groups[p]) < min_size),
+                                  key=lambda p: len(groups[p])):
+                    if pid not in groups or len(groups) == 1:
+                        continue
+                    idx = groups[pid]
+                    best, best_d = None, np.inf
+                    for opid, oidx in groups.items():
+                        if opid == pid:
+                            continue
+                        d = dist[np.ix_(idx, oidx)].min()
+                        if d < best_d:
+                            best_d, best = d, opid
+                    if best is not None:
+                        groups[best] = groups[best] + idx
+                        del groups[pid]
+                        changed = True
+
+        return {pid: [filenames[i] for i in idx] for pid, idx in groups.items()}
+
+    def _rebalance_keep_k(self, manifest, filenames, embeddings, k,
+                          min_size=None, max_size=None, max_iter=1000):
+        """Rebalance group sizes while keeping exactly k groups.
+
+        Each pass splits the largest group into two balanced halves — bisecting
+        k-means on the siamese embeddings, which cuts at the natural density gap
+        and favours an even split rather than shaving off one or two points —
+        and merges the smallest group into its nearest neighbour, so the group
+        count is preserved. Repeats until every group is within
+        [min_size, max_size]. It is a no-op when no group violates the bounds
+        (e.g. Task06, where edge-ranking already gives clean per-patient sizes).
+        """
+        from sklearn.cluster import KMeans
+
+        f2i = {f: i for i, f in enumerate(filenames)}
+        dist = 1.0 - self._dense_prob_matrix(embeddings)
+        np.fill_diagonal(dist, 0.0)
+        dist[dist < 0] = 0.0
+        clusters = [[f2i[f] for f in fs if f in f2i] for fs in manifest.values()]
+
+        def bisect(idx):
+            lab = KMeans(n_clusters=2, n_init=10, random_state=42).fit_predict(
+                embeddings[idx])
+            a = [idx[i] for i in range(len(idx)) if lab[i] == 0]
+            b = [idx[i] for i in range(len(idx)) if lab[i] == 1]
+            if not a or not b:                 # degenerate: even halves
+                half = len(idx) // 2
+                a, b = idx[:half], idx[half:]
+            return a, b
+
+        def merge_smallest():
+            si = min(range(len(clusters)), key=lambda i: len(clusters[i]))
+            small = clusters.pop(si)
+            ji = min(range(len(clusters)),
+                     key=lambda j: dist[np.ix_(small, clusters[j])].min())
+            clusters[ji].extend(small)
+
+        # Start from exactly k groups (edge-rank already yields k; be safe).
+        while len(clusters) > k:
+            merge_smallest()
+        while len(clusters) < k:
+            bi = max(range(len(clusters)), key=lambda i: len(clusters[i]))
+            clusters += list(bisect(clusters.pop(bi)))
+
+        hi = max_size if max_size else float("inf")
+        lo = min_size if min_size else 0
+        for _ in range(max_iter):
+            sizes = [len(c) for c in clusters]
+            if max(sizes) <= hi and min(sizes) >= lo:
+                break
+            bi = max(range(len(clusters)), key=lambda i: len(clusters[i]))
+            clusters += list(bisect(clusters.pop(bi)))   # +1 group
+            merge_smallest()                              # -1 group -> keeps k
+
+        return {f"siamese_{i}": [filenames[j] for j in c]
+                for i, c in enumerate(clusters)}
+
+    def _edge_rank_cluster(self, edges, n, filenames, k):
+        """Cluster into k groups by adding edges in descending confidence order.
+
+        Starts with n isolated nodes, adds edges from most to least confident
+        until exactly k connected components form. Natural clusters stay intact;
+        only the weakest-split cluster gets divided.
+        Returns None if we can't reach k components.
+        """
+        from scipy.sparse.csgraph import connected_components
+        from scipy.sparse import csr_matrix
+
+        if len(edges) == 0:
+            return None
+
+        # Sort edges by probability descending
+        sorted_edges = sorted(edges, key=lambda e: e[2], reverse=True)
+
+        # Binary search for the threshold that gives exactly k components
+        # Start with all edges → see how many components
+        all_row = [e[0] for e in sorted_edges] + [e[1] for e in sorted_edges]
+        all_col = [e[1] for e in sorted_edges] + [e[0] for e in sorted_edges]
+        all_data = np.ones(len(all_row))
+        adj_full = csr_matrix((all_data, (all_row, all_col)), shape=(n, n))
+        n_min, _ = connected_components(adj_full, directed=False)
+
+        if n_min > k:
+            print(f"  Edge ranking: even with all {len(edges)} edges, "
+                  f"only {n_min} components (need {k})", flush=True)
+            return None
+        if n_min == k:
+            # Perfect — all edges give exactly k components
+            _, labels = connected_components(adj_full, directed=False)
+            return self._labels_to_manifest(labels, k, filenames, "siamese")
+
+        # Binary search: find edge index where components == k
+        lo, hi = 0, len(sorted_edges)
+        best_labels = None
+        while lo < hi:
+            mid = (lo + hi) // 2
+            # Build graph with first `mid` edges
+            sub_edges = sorted_edges[:mid]
+            row = [e[0] for e in sub_edges] + [e[1] for e in sub_edges]
+            col = [e[1] for e in sub_edges] + [e[0] for e in sub_edges]
+            data = np.ones(len(row))
+            adj = csr_matrix((data, (row, col)), shape=(n, n))
+            n_comp, labels = connected_components(adj, directed=False)
+
+            if n_comp > k:
+                lo = mid + 1  # need more edges
+            elif n_comp < k:
+                hi = mid  # too many edges
+            else:
+                best_labels = labels
+                hi = mid  # try to find the earliest edge that achieves k
+
+        if best_labels is None:
+            return None
+
+        # Report the confidence at the split point
+        split_conf = sorted_edges[lo - 1][2] if lo > 0 else 1.0
+        print(f"  Edge ranking: k={k} reached at confidence={split_conf:.4f} "
+              f"(edge {lo}/{len(sorted_edges)})", flush=True)
+        return self._labels_to_manifest(best_labels, k, filenames, "siamese")
+
+    def _spectral_cluster(self, sim_matrix, embeddings, filenames, k):
+        """Cluster into exactly k groups via spectral clustering on the
+        siamese similarity graph.
+
+        Builds a weighted adjacency matrix from cosine similarity, then uses
+        spectral clustering (normalized cut) to partition into k groups.
+        Weak/spurious edges get cut to respect the known group count.
+        """
+        from sklearn.cluster import SpectralClustering
+        from scipy.sparse import csr_matrix
+
+        n = len(embeddings)
+
+        # Build sparse weighted adjacency from top similarities
+        # Use top_k=50 for denser graph (spectral clustering needs connectivity)
+        top_k_dense = min(50, n - 1)
+        top_idx = np.argsort(-sim_matrix, axis=1)[:, 1:top_k_dense + 1]
+
+        row, col, data = [], [], []
+        for i in range(n):
+            for j in top_idx[i]:
+                if j <= i:
+                    continue
+                # Weight = cosine similarity (in [0,1] after ReLU)
+                w = max(0.0, float(sim_matrix[i, j]))
+                if w > 0:
+                    row.append(i); col.append(j); data.append(w)
+                    row.append(j); col.append(i); data.append(w)
+
+        adj = csr_matrix((data, (row, col)), shape=(n, n))
+        print(f"  Adjacency: {len(data) // 2} edges (top-{top_k_dense})",
+              flush=True)
+
+        # Spectral clustering
+        print(f"  Spectral clustering into k={k} groups ...", flush=True)
+        sc = SpectralClustering(
+            n_clusters=k, affinity="precomputed",
+            random_state=42, n_init=20,
+            assign_labels="kmeans")  # kmeans discretization is more stable
+        labels = sc.fit_predict(adj.toarray())
+
+        # If spectral clustering fails (disconnected graph), fall back to
+        # adding a small epsilon to connect components
+        n_unique = len(np.unique(labels))
+        if n_unique < k:
+            print(f"  Warning: only {n_unique}/{k} clusters found; "
+                  f"graph may be disconnected. Adding background connectivity.",
+                  flush=True)
+            # Add weak background edges
+            adj_dense = adj.toarray()
+            adj_dense += 0.001 * (1.0 - np.eye(n))
+            labels = SpectralClustering(
+                n_clusters=k, affinity="precomputed",
+                random_state=42, n_init=20,
+                assign_labels="kmeans").fit_predict(adj_dense)
+
+        return self._labels_to_manifest(labels, k, filenames, "spectral")
+
+    def _labels_to_manifest(self, labels, n_groups, filenames, prefix):
+        """Convert flat cluster labels to a manifest dict."""
+        components = defaultdict(list)
+        for i, c in enumerate(labels):
+            components[int(c)].append(filenames[i])
+
+        manifest = {}
+        for comp_id, comp_files in sorted(components.items()):
+            manifest[f"{prefix}_{comp_id:03d}"] = sorted(comp_files)
+
+        print(f"  {len(manifest)} groups from {len(filenames)} slices",
+              flush=True)
+
+        sizes = [len(v) for v in manifest.values()]
+        if sizes:
+            print(f"  Group sizes: min={min(sizes)}, max={max(sizes)}, "
+                  f"mean={np.mean(sizes):.1f}", flush=True)
+
+        return manifest

plan.md

@@ -0,0 +1,116 @@
+# Implementation Plan — Potter et al. Data Leakage Reproduction
+
+## Manuscript pipeline (from docx)
+
+1. **Feature extraction** — 5 CNNs (VGG16, DenseNet121, EfficientNetB1, MobileNetV2, ResNet50), frozen, include_top=False
+2. **PCA + t-SNE** — visualize feature space (Figure 5 / S2)
+3. **Patient clustering** — K-means on 64×64 grayscale thumbnails, per-class (15/40/55)
+4. **Random Forest** — feature importance ranking (2000 trees → we use 500, equivalent)
+5. **SVM with GridSearchCV** — grid over γ per nfeatures, C=10 fixed, 5-fold CV
+6. **Image-level vs patient-level** — compare test accuracy with both split types
+7. **20-seed repetition** — boxplots showing distribution (Figure 6 / S4)
+
+## Current status
+
+### ✅ Implemented & Working
+
+| Component | Location | Notes |
+|---|---|---|
+| Feature extraction (5 models, PyTorch) | `classes/features.py` | ResNet50 GAP stripped (100K-d). EfficientNetB1 dim error documented |
+| Patient clustering (K-means) | `classes/patient_identifier.py` | Thumbnail + feature-based modes |
+| Classification pipeline | `classes/classifier.py` | `PatientLeakageClassifier.run()` — RF→GridSearchCV→SVM |
+| t-SNE visualization | `scripts/visualizations/simple_patient_tsne.py` | Per-class patient coloring |
+| Figure 6 (20-seed boxplots) | `scripts/visualizations/figure6.py` | Image vs patient distributions |
+| Figure S4 (CV accuracy curves) | `scripts/visualizations/figure_s4.py` | Per-model, log₁₀ x-scale, patient-level |
+| C×γ grid search | archived | Confirmed C=10 is adequate |
+| SVM class-weight analysis | archived | Confirmed no benefit |
+| Task06 NIfTI dataset loader | `classes/nifti_dataset.py` | HU windowing, orientation control |
+| Clustering validation (Task06) | `scripts/validate_clustering_task06.py` | Thumbnail + feature K-means vs ground truth |
+| Clustering comparison (IQ-OTH) | `scripts/analysis/compare_clustering_methods.py` | Thumbnail vs feature vs v8 reference |
+| Siamese patient matching | `classes/siamese.py` | `SiamesePatientMatcher` class |
+| Siamese training | `scripts/train_siamese.py` | Task06 train/test split, hard negatives |
+
+### 🚧 In Progress
+
+| Component | Status |
+|---|---|
+| Siamese-based patient manifest for IQ-OTH | `scripts/siamese_identify.py` written, needs trained model |
+
+### ❌ Still Needed
+
+| Manuscript Figure/Table | What we need | Priority |
+|---|---|---|
+| **Figure 5** (PCA + t-SNE) | Per-model t-SNE plots, matching manuscript style | Medium |
+| **Figure S1** (image-level CV curves) | 5 panels, image-level split, log₁₀ x-scale. Nearly identical to S4 but with image split | Low |
+| **Figure S2** (PCA/t-SNE per model) | Like Figure 5 but for all 5 models | Medium |
+| **Figure S3** (example cluster images) | 5 example slices from a single K-means cluster per class | Low |
+| **Figure S5** (confusion matrix) | VGG16 patient-level confusion matrix | Low |
+| **Table 2 comparison** | Run classification with thumbnail K-means manifest to match manuscript numbers | High |
+| **Siamese results integration** | Once trained: run siamese_identify, compare manifest against K-means manifests | High |
+| **Final patient manifest** | Choose best method, produce canonical patient assignments | High |
+
+## Proposed final structure
+
+```
+patient_leakage_detector/
+├── classes/
+│   ├── features.py              # FeatureExtractor (5 CNNs)
+│   ├── patient_identifier.py    # PatientIdentifier (K-means clustering)
+│   ├── classifier.py            # PatientLeakageClassifier (RF→SVM pipeline)
+│   ├── nifti_dataset.py         # NiftiSliceDataset (Task06 loader)
+│   └── siamese.py               # SiamesePatientMatcher (learned matching)
+│
+├── scripts/
+│   ├── classification.py        # Core: single-seed classification
+│   ├── train_siamese.py         # Core: siamese training on Task06
+│   ├── siamese_identify.py      # Core: IQ-OTH patient manifest via siamese
+│   ├── validate_clustering_task06.py  # Core: Task06 validation
+│   │
+│   ├── visualizations/
+│   │   ├── figure5.py           # [TODO] PCA + t-SNE plots
+│   │   ├── figure6.py           # ✅ 20-seed boxplots
+│   │   ├── figure_s1.py         # [TODO] Image-level CV curves
+│   │   ├── figure_s2.py         # [TODO] PCA/t-SNE per model
+│   │   ├── figure_s3.py         # [TODO] Example cluster images
+│   │   ├── figure_s4.py         # ✅ Patient-level CV curves
+│   │   ├── figure_s5.py         # [TODO] Confusion matrix
+│   │   └── simple_patient_tsne.py  # ✅ t-SNE with patient coloring
+│   │
+│   └── analysis/
+│       ├── compare_clustering_methods.py  # ✅ Thumbnail vs feature vs v8
+│       └── verify_feature_dims.py        # ✅ TF dimension verification
+│
+├── features/
+│   ├── VGG16_features.npz       # Extracted features (all 5 models)
+│   └── task06_lung/             # Cached Task06 features
+│
+├── models/
+│   └── siamese_resnet18.pt      # Trained siamese model
+│
+├── results/
+│   ├── simple_patient_manifest.csv     # Feature-based K-means manifest
+│   ├── siamese_manifest.csv            # [TODO] Siamese-based manifest
+│   ├── classification_results.json     # Single-seed results (all 5 models)
+│   ├── figure6_data.json              # 20-seed data
+│   ├── task06_clustering_validation.json  # Task06 validation results
+│   └── ...
+│
+├── plots/
+│   ├── figure6.png
+│   ├── figure_s4.png
+│   └── ...
+│
+├── plan.md                     # This file
+├── notes.md                    # Review notes for supervisor
+├── .gitignore
+└── .archive/                   # Superseded scripts and results
+```
+
+## Immediate next steps
+
+1. **Finish siamese training** → run `python scripts/siamese_identify.py` → produce `siamese_manifest.csv`
+2. **Run classification with thumbnail K-means manifest** to reproduce manuscript Table 2 numbers
+3. **Compare manifests**: K-means thumbnail vs K-means feature vs siamese — which gives the best patient-level split?
+4. **Choose canonical manifest** and produce final classification results
+5. **Generate remaining figures** (Figure 5, S1, S2, S3, S5)
+6. **Draft findings for supervisor discussion**

results/benchmark_results.json

@@ -0,0 +1,90 @@
+[
+  {
+    "split": "image",
+    "variant": "A_baseline",
+    "cv": 0.9909,
+    "test": 1.0,
+    "nfeat": 5000,
+    "C": 10,
+    "gamma": 2.667e-05,
+    "svm_weighted": false,
+    "c_gridded": false
+  },
+  {
+    "split": "image",
+    "variant": "B_svm_balanced",
+    "cv": 0.9875,
+    "test": 1.0,
+    "nfeat": 2000,
+    "C": 10,
+    "gamma": 0.00028117,
+    "svm_weighted": true,
+    "c_gridded": false
+  },
+  {
+    "split": "image",
+    "variant": "C_c_gamma_grid",
+    "cv": 0.992,
+    "test": 1.0,
+    "nfeat": 1000,
+    "C": 100,
+    "gamma": 3.162e-05,
+    "svm_weighted": false,
+    "c_gridded": true
+  },
+  {
+    "split": "image",
+    "variant": "D_combined",
+    "cv": 0.9886,
+    "test": 1.0,
+    "nfeat": 2000,
+    "C": 100,
+    "gamma": 6.668e-05,
+    "svm_weighted": true,
+    "c_gridded": true
+  },
+  {
+    "split": "patient",
+    "variant": "A_baseline",
+    "cv": 0.91,
+    "test": 0.8978,
+    "nfeat": 300,
+    "C": 10,
+    "gamma": 0.00010541,
+    "svm_weighted": false,
+    "c_gridded": false
+  },
+  {
+    "split": "patient",
+    "variant": "B_svm_balanced",
+    "cv": 0.9042,
+    "test": 0.8578,
+    "nfeat": 300,
+    "C": 10,
+    "gamma": 0.00010541,
+    "svm_weighted": true,
+    "c_gridded": false
+  },
+  {
+    "split": "patient",
+    "variant": "C_c_gamma_grid",
+    "cv": 0.91,
+    "test": 0.8978,
+    "nfeat": 300,
+    "C": 10,
+    "gamma": 0.00010541,
+    "svm_weighted": false,
+    "c_gridded": true
+  },
+  {
+    "split": "patient",
+    "variant": "D_combined",
+    "cv": 0.9052,
+    "test": 0.9022,
+    "nfeat": 300,
+    "C": 1000,
+    "gamma": 0.00010541,
+    "svm_weighted": true,
+    "c_gridded": true
+  }
+]

results/c_gamma_grid_results.json

@@ -0,0 +1,202 @@
+[
+  {
+    "cv": 0.9920324675324675,
+    "test": 1.0,
+    "nfeat": 5000,
+    "C": 10.0,
+    "gamma": 4.308869380063768e-05,
+    "model": "VGG16",
+    "split_type": "image",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9920324675324675,
+    "test": 1.0,
+    "nfeat": 1000,
+    "C": 100.0,
+    "gamma": 3.1622776601683795e-05,
+    "model": "VGG16",
+    "split_type": "image",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9121525126693927,
+    "test": 0.8933333333333333,
+    "nfeat": 200,
+    "C": 10.0,
+    "gamma": 0.00015811388300841897,
+    "model": "VGG16",
+    "split_type": "patient",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9121525126693927,
+    "test": 0.8933333333333333,
+    "nfeat": 200,
+    "C": 10.0,
+    "gamma": 0.00015811388300841897,
+    "model": "VGG16",
+    "split_type": "patient",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9874675324675325,
+    "test": 0.9954545454545455,
+    "nfeat": 5000,
+    "C": 10.0,
+    "gamma": 4.308869380063768e-05,
+    "model": "DenseNet121",
+    "split_type": "image",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9874675324675325,
+    "test": 0.9954545454545455,
+    "nfeat": 5000,
+    "C": 10.0,
+    "gamma": 4.308869380063768e-05,
+    "model": "DenseNet121",
+    "split_type": "image",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9367063276432266,
+    "test": 0.8711111111111111,
+    "nfeat": 1000,
+    "C": 10.0,
+    "gamma": 3.1622776601683795e-05,
+    "model": "DenseNet121",
+    "split_type": "patient",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9367063276432266,
+    "test": 0.8711111111111111,
+    "nfeat": 1000,
+    "C": 10.0,
+    "gamma": 3.1622776601683795e-05,
+    "model": "DenseNet121",
+    "split_type": "patient",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9829155844155844,
+    "test": 0.9954545454545455,
+    "nfeat": 5000,
+    "C": 10.0,
+    "gamma": 4.308869380063768e-05,
+    "model": "EfficientNetB1",
+    "split_type": "image",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9829155844155844,
+    "test": 0.9954545454545455,
+    "nfeat": 5000,
+    "C": 10.0,
+    "gamma": 4.308869380063768e-05,
+    "model": "EfficientNetB1",
+    "split_type": "image",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9158173382359915,
+    "test": 0.8844444444444445,
+    "nfeat": 1000,
+    "C": 10.0,
+    "gamma": 8.254041852680186e-05,
+    "model": "EfficientNetB1",
+    "split_type": "patient",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9158173382359915,
+    "test": 0.8844444444444445,
+    "nfeat": 1000,
+    "C": 10.0,
+    "gamma": 8.254041852680186e-05,
+    "model": "EfficientNetB1",
+    "split_type": "patient",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9874610389610389,
+    "test": 0.9954545454545455,
+    "nfeat": 2000,
+    "C": 10.0,
+    "gamma": 0.0001077217345015942,
+    "model": "MobileNetV2",
+    "split_type": "image",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9874610389610389,
+    "test": 0.9954545454545455,
+    "nfeat": 2000,
+    "C": 10.0,
+    "gamma": 0.0001077217345015942,
+    "model": "MobileNetV2",
+    "split_type": "image",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9129576751848534,
+    "test": 0.8888888888888888,
+    "nfeat": 750,
+    "C": 10.0,
+    "gamma": 4.216370213557839e-05,
+    "model": "MobileNetV2",
+    "split_type": "patient",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9129576751848534,
+    "test": 0.8888888888888888,
+    "nfeat": 750,
+    "C": 10.0,
+    "gamma": 4.216370213557839e-05,
+    "model": "MobileNetV2",
+    "split_type": "patient",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9897597402597402,
+    "test": 0.9954545454545455,
+    "nfeat": 1500,
+    "C": 10.0,
+    "gamma": 0.00014362897933545892,
+    "model": "ResNet50",
+    "split_type": "image",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9908961038961038,
+    "test": 0.9954545454545455,
+    "nfeat": 1500,
+    "C": 100.0,
+    "gamma": 2.1081851067789193e-05,
+    "model": "ResNet50",
+    "split_type": "image",
+    "variant": "C_grid"
+  },
+  {
+    "cv": 0.9156075191986195,
+    "test": 0.88,
+    "nfeat": 750,
+    "C": 10.0,
+    "gamma": 4.216370213557839e-05,
+    "model": "ResNet50",
+    "split_type": "patient",
+    "variant": "baseline_C10"
+  },
+  {
+    "cv": 0.9156075191986195,
+    "test": 0.88,
+    "nfeat": 750,
+    "C": 10.0,
+    "gamma": 4.216370213557839e-05,
+    "model": "ResNet50",
+    "split_type": "patient",
+    "variant": "C_grid"
+  }
+]

results/classification_pca50.json

@@ -0,0 +1,42 @@
+[
+  {
+    "model": "VGG16",
+    "image_cv": 0.9920324675324675,
+    "image_test": 1.0,
+    "patient_cv": 0.9121525126693927,
+    "patient_test": 0.8933333333333333,
+    "drop": 0.10666666666666669
+  },
+  {
+    "model": "DenseNet121",
+    "image_cv": 0.9874675324675325,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9367063276432266,
+    "patient_test": 0.8711111111111111,
+    "drop": 0.12434343434343442
+  },
+  {
+    "model": "EfficientNetB1",
+    "image_cv": 0.9829155844155844,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9158173382359915,
+    "patient_test": 0.8844444444444445,
+    "drop": 0.111010101010101
+  },
+  {
+    "model": "MobileNetV2",
+    "image_cv": 0.9874610389610389,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9129576751848534,
+    "patient_test": 0.8888888888888888,
+    "drop": 0.10656565656565664
+  },
+  {
+    "model": "ResNet50",
+    "image_cv": 0.9669675324675324,
+    "image_test": 0.9772727272727273,
+    "patient_cv": 0.8742211426908154,
+    "patient_test": 0.8622222222222222,
+    "drop": 0.11505050505050507
+  }
+]

results/classification_results.json

@@ -1,52 +0,0 @@
-[
-  {
-    "model": "VGG16",
-    "image_cv": 0.992,
-    "image_test": 1.0,
-    "image_nfeat": 5000,
-    "patient_cv": 0.9122,
-    "patient_test": 0.8933,
-    "patient_nfeat": 200,
-    "drop": 0.1067
-  },
-  {
-    "model": "DenseNet121",
-    "image_cv": 0.9875,
-    "image_test": 0.9955,
-    "image_nfeat": 5000,
-    "patient_cv": 0.9367,
-    "patient_test": 0.8711,
-    "patient_nfeat": 1000,
-    "drop": 0.1243
-  },
-  {
-    "model": "EfficientNetB1",
-    "image_cv": 0.9829,
-    "image_test": 0.9955,
-    "image_nfeat": 5000,
-    "patient_cv": 0.9158,
-    "patient_test": 0.8844,
-    "patient_nfeat": 1000,
-    "drop": 0.111
-  },
-  {
-    "model": "MobileNetV2",
-    "image_cv": 0.9875,
-    "image_test": 0.9955,
-    "image_nfeat": 2000,
-    "patient_cv": 0.913,
-    "patient_test": 0.8889,
-    "patient_nfeat": 750,
-    "drop": 0.1066
-  },
-  {
-    "model": "ResNet50",
-    "image_cv": 0.967,
-    "image_test": 0.9773,
-    "image_nfeat": 750,
-    "patient_cv": 0.8742,
-    "patient_test": 0.8622,
-    "patient_nfeat": 400,
-    "drop": 0.1151
-  }
-]

results/classification_siamese.json

@@ -0,0 +1,47 @@
+[
+  {
+    "model": "VGG16",
+    "manifest": "siamese",
+    "image_cv": 0.9920324675324675,
+    "image_test": 1.0,
+    "patient_cv": 0.9098522338879679,
+    "patient_test": 0.8473684210526315,
+    "drop": 0.15263157894736845
+  },
+  {
+    "model": "DenseNet121",
+    "manifest": "siamese",
+    "image_cv": 0.9874675324675325,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9153328830764822,
+    "patient_test": 0.8842105263157894,
+    "drop": 0.11124401913875603
+  },
+  {
+    "model": "EfficientNetB1",
+    "manifest": "siamese",
+    "image_cv": 0.9829155844155844,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9173969014142157,
+    "patient_test": 0.8526315789473684,
+    "drop": 0.14282296650717707
+  },
+  {
+    "model": "MobileNetV2",
+    "manifest": "siamese",
+    "image_cv": 0.9874610389610389,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9252064018337733,
+    "patient_test": 0.8736842105263158,
+    "drop": 0.12177033492822964
+  },
+  {
+    "model": "ResNet50",
+    "manifest": "siamese",
+    "image_cv": 0.9897597402597402,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9317341028632242,
+    "patient_test": 0.8631578947368421,
+    "drop": 0.13229665071770336
+  }
+]

results/classification_thumbnail.json

@@ -0,0 +1,47 @@
+[
+  {
+    "model": "VGG16",
+    "manifest": "thumbnail",
+    "image_cv": 0.9920324675324675,
+    "image_test": 1.0,
+    "patient_cv": 0.8856369070229704,
+    "patient_test": 0.9125,
+    "drop": 0.08750000000000002
+  },
+  {
+    "model": "DenseNet121",
+    "manifest": "thumbnail",
+    "image_cv": 0.9874675324675325,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9183953769749152,
+    "patient_test": 0.9,
+    "drop": 0.09545454545454546
+  },
+  {
+    "model": "EfficientNetB1",
+    "manifest": "thumbnail",
+    "image_cv": 0.9829155844155844,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9032050186948257,
+    "patient_test": 0.925,
+    "drop": 0.07045454545454544
+  },
+  {
+    "model": "MobileNetV2",
+    "manifest": "thumbnail",
+    "image_cv": 0.9874610389610389,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.9055959867591996,
+    "patient_test": 0.8625,
+    "drop": 0.13295454545454544
+  },
+  {
+    "model": "ResNet50",
+    "manifest": "thumbnail",
+    "image_cv": 0.9897597402597402,
+    "image_test": 0.9954545454545455,
+    "patient_cv": 0.910289858567831,
+    "patient_test": 0.9166666666666666,
+    "drop": 0.07878787878787885
+  }
+]

results/clustering_validation/validate_pca50_VGG16.json

@@ -0,0 +1,20 @@
+{
+  "method": "pca50_VGG16",
+  "model": "VGG16",
+  "n_patients": 8,
+  "n_slices": 1251,
+  "ARI": 0.3838311521430662,
+  "NMI": 0.588535400629693,
+  "cluster_purity": {
+    "overall": 0.6418864908073542,
+    "median": 0.6389547237659365,
+    "mean": 0.7006446210824684,
+    "frac_gt_70": 0.375,
+    "frac_gt_90": 0.375
+  },
+  "patient_capture": {
+    "median": 0.7111548252719277,
+    "mean": 0.6921309007056184,
+    "frac_gt_50": 0.875
+  }
+}

results/clustering_validation/validate_pca50_VGG16_full.json

@@ -0,0 +1,20 @@
+{
+  "method": "pca50_VGG16",
+  "model": "VGG16",
+  "n_patients": 63,
+  "n_slices": 10587,
+  "ARI": 0.2982860430503409,
+  "NMI": 0.6543620808930918,
+  "cluster_purity": {
+    "overall": 0.4898460375932748,
+    "median": 0.8548387096774194,
+    "mean": 0.6433234955790892,
+    "frac_gt_70": 0.5396825396825397,
+    "frac_gt_90": 0.49206349206349204
+  },
+  "patient_capture": {
+    "median": 0.4207650273224044,
+    "mean": 0.4439023854165772,
+    "frac_gt_50": 0.31746031746031744
+  }
+}

results/clustering_validation/validate_siamese_cc.json

@@ -0,0 +1,20 @@
+{
+  "ARI": 0.9049336824308681,
+  "NMI": 0.9575945510291562,
+  "cluster_purity": {
+    "overall": 0.8848920863309353,
+    "median": 1.0,
+    "mean": 0.9288185862580326,
+    "frac_gt_70": 0.8571428571428571,
+    "frac_gt_90": 0.8571428571428571
+  },
+  "patient_capture": {
+    "median": 1.0,
+    "mean": 1.0,
+    "frac_gt_50": 1.0
+  },
+  "method": "cc",
+  "threshold": 0.9,
+  "n_patients": 8,
+  "n_slices": 1251
+}

results/clustering_validation/validate_siamese_cc_full.json

@@ -0,0 +1,20 @@
+{
+  "ARI": 0.9430129669411811,
+  "NMI": 0.9865730997479965,
+  "cluster_purity": {
+    "overall": 0.9331255313119864,
+    "median": 1.0,
+    "mean": 0.9592822768032196,
+    "frac_gt_70": 0.9122807017543859,
+    "frac_gt_90": 0.9122807017543859
+  },
+  "patient_capture": {
+    "median": 1.0,
+    "mean": 1.0,
+    "frac_gt_50": 1.0
+  },
+  "method": "cc",
+  "threshold": 0.9,
+  "n_patients": 63,
+  "n_slices": 10587
+}

results/clustering_validation/validate_siamese_sc.json

@@ -0,0 +1,20 @@
+{
+  "ARI": 1.0,
+  "NMI": 1.0,
+  "cluster_purity": {
+    "overall": 1.0,
+    "median": 1.0,
+    "mean": 1.0,
+    "frac_gt_70": 1.0,
+    "frac_gt_90": 1.0
+  },
+  "patient_capture": {
+    "median": 1.0,
+    "mean": 1.0,
+    "frac_gt_50": 1.0
+  },
+  "method": "sc",
+  "threshold": 0.9,
+  "n_patients": 8,
+  "n_slices": 1251
+}

results/clustering_validation/validate_siamese_sc_full.json

@@ -0,0 +1,20 @@
+{
+  "ARI": 0.9872517295956881,
+  "NMI": 0.9973073869067276,
+  "cluster_purity": {
+    "overall": 0.9852649475772174,
+    "median": 1.0,
+    "mean": 0.9922860109775998,
+    "frac_gt_70": 0.9841269841269841,
+    "frac_gt_90": 0.9841269841269841
+  },
+  "patient_capture": {
+    "median": 1.0,
+    "mean": 0.9995767195767196,
+    "frac_gt_50": 1.0
+  },
+  "method": "sc",
+  "threshold": 0.9,
+  "n_patients": 63,
+  "n_slices": 10587
+}

results/clustering_validation/validate_thumbnail.json

@@ -0,0 +1,19 @@
+{
+  "method": "thumbnail_64x64",
+  "n_patients": 8,
+  "n_slices": 1251,
+  "ARI": 0.40028837820172775,
+  "NMI": 0.6177421217683747,
+  "cluster_purity": {
+    "overall": 0.697841726618705,
+    "median": 0.8719512195121951,
+    "mean": 0.7807546888218818,
+    "frac_gt_70": 0.625,
+    "frac_gt_90": 0.5
+  },
+  "patient_capture": {
+    "median": 0.6048850574712643,
+    "mean": 0.650025535852033,
+    "frac_gt_50": 0.75
+  }
+}

results/clustering_validation/validate_thumbnail_full.json

@@ -0,0 +1,19 @@
+{
+  "method": "thumbnail_64x64",
+  "n_patients": 63,
+  "n_slices": 10587,
+  "ARI": 0.279398847493094,
+  "NMI": 0.6671641782310866,
+  "cluster_purity": {
+    "overall": 0.4804004911684141,
+    "median": 0.5984251968503937,
+    "mean": 0.6204889734832222,
+    "frac_gt_70": 0.4603174603174603,
+    "frac_gt_90": 0.38095238095238093
+  },
+  "patient_capture": {
+    "median": 0.38620689655172413,
+    "mean": 0.4227929731937232,
+    "frac_gt_50": 0.30158730158730157
+  }
+}

results/compare_patient_groupings.json

@@ -0,0 +1,69 @@
+{
+  "n_images": 1097,
+  "structure": {
+    "siamese": {
+      "n_groups": 110,
+      "min": 1,
+      "median": 4.0,
+      "mean": 9.972727272727273,
+      "max": 112,
+      "singletons": 25
+    },
+    "pca50": {
+      "n_groups": 110,
+      "min": 2,
+      "median": 9.0,
+      "mean": 9.972727272727273,
+      "max": 32,
+      "singletons": 0
+    },
+    "thumbnail": {
+      "n_groups": 110,
+      "min": 1,
+      "median": 8.5,
+      "mean": 9.972727272727273,
+      "max": 30,
+      "singletons": 1
+    }
+  },
+  "pairwise": {
+    "siamese_vs_pca50": {
+      "ARI": 0.4181618426201391,
+      "NMI": 0.8689785640560458,
+      "V": 0.8689785640560458
+    },
+    "siamese_vs_thumbnail": {
+      "ARI": 0.4147240727370461,
+      "NMI": 0.8703177137592031,
+      "V": 0.8703177137592032
+    },
+    "pca50_vs_thumbnail": {
+      "ARI": 0.7432698225746973,
+      "NMI": 0.9226193956118858,
+      "V": 0.9226193956118858
+    }
+  },
+  "per_class_ARI": {
+    "siamese_vs_pca50": {
+      "Benign": 0.1917309601351469,
+      "Malignant": 0.387952185974747,
+      "Normal": 0.5325854014848859
+    },
+    "siamese_vs_thumbnail": {
+      "Benign": 0.23555079431438128,
+      "Malignant": 0.38636842596539783,
+      "Normal": 0.5002427880113821
+    },
+    "pca50_vs_thumbnail": {
+      "Benign": 0.5478101559111155,
+      "Malignant": 0.77928612571208,
+      "Normal": 0.6875861344042246
+    }
+  },
+  "mean_ari_vs_others": {
+    "siamese": 0.41644295767859263,
+    "pca50": 0.5807158325974182,
+    "thumbnail": 0.5789969476558717
+  },
+  "most_distinct": "siamese"
+}

results/siamese_manifest.csv

@@ -0,0 +1,111 @@
+patient_id,class,n_images,images
+Benign_0,Unknown,7,B_094;B_095;B_096;B_097;B_098;B_099;B_100
+Benign_1,Unknown,10,B_101;B_102;B_103;B_104;B_105;B_106;B_107;B_108;B_109;B_110
+Benign_10,Unknown,9,B_053;B_054;B_055;B_066;B_067;B_074;B_075;B_092;B_093
+Benign_11,Unknown,7,B_014;B_015;B_016;B_017;B_018;B_019;B_020
+Benign_12,Unknown,21,B_009;B_010;B_011;B_012;B_013;B_021;B_022;B_023;B_024;B_025;B_035;B_042;B_043;B_044;B_045;B_086;B_087;B_088;B_089;B_090;B_091
+Benign_13,Unknown,18,B_001;B_002;B_003;B_026;B_027;B_028;B_029;B_056;B_057;B_058;B_059;B_068;B_069;B_076;B_077;B_078;B_079;B_080
+Benign_14,Unknown,9,B_070;B_071;B_072;B_073;B_081;B_082;B_083;B_084;B_085
+Benign_2,Unknown,7,B_111;B_112;B_113;B_114;B_115;B_116;B_120
+Benign_3,Unknown,3,B_117;B_118;B_119
+Benign_4,Unknown,5,B_030;B_031;B_032;B_033;B_034
+Benign_5,Unknown,6,B_036;B_037;B_038;B_039;B_040;B_041
+Benign_6,Unknown,5,B_047;B_048;B_049;B_050;B_051
+Benign_7,Unknown,4,B_052;B_046;B_004;B_005
+Benign_8,Unknown,3,B_006;B_007;B_008
+Benign_9,Unknown,6,B_060;B_061;B_062;B_063;B_064;B_065
+Malignant_0,Unknown,26,M_001;M_017;M_018;M_539;M_540;M_541;M_543;M_544;M_545;M_546;M_547;M_548;M_549;M_550;M_551;M_552;M_553;M_554;M_555;M_556;M_557;M_558;M_559;M_560;M_561;M_527
+Malignant_1,Unknown,15,M_016;M_523;M_524;M_525;M_526;M_528;M_529;M_530;M_531;M_532;M_533;M_534;M_535;M_536;M_537
+Malignant_10,Unknown,5,M_415;M_416;M_417;M_418;M_419
+Malignant_11,Unknown,9,M_044;M_045;M_046;M_047;M_048;M_049;M_050;M_051;M_052
+Malignant_12,Unknown,9,M_441;M_442;M_443;M_444;M_445;M_446;M_447;M_448;M_440
+Malignant_13,Unknown,10,M_449;M_450;M_451;M_452;M_453;M_454;M_455;M_456;M_457;M_458
+Malignant_14,Unknown,3,M_459;M_460;M_461
+Malignant_15,Unknown,11,M_483;M_484;M_485;M_486;M_487;M_488;M_489;M_490;M_491;M_492;M_493
+Malignant_16,Unknown,10,M_494;M_499;M_500;M_501;M_502;M_503;M_504;M_505;M_506;M_507
+Malignant_17,Unknown,4,M_495;M_496;M_497;M_498
+Malignant_18,Unknown,15,M_508;M_509;M_510;M_511;M_512;M_513;M_514;M_515;M_516;M_517;M_518;M_519;M_520;M_521;M_522
+Malignant_19,Unknown,4,M_079;M_080;M_081;M_082
+Malignant_2,Unknown,12,M_184;M_185;M_186;M_187;M_188;M_189;M_190;M_191;M_180;M_183;M_181;M_182
+Malignant_20,Unknown,4,M_083;M_084;M_085;M_086
+Malignant_21,Unknown,22,M_089;M_090;M_091;M_092;M_093;M_094;M_095;M_096;M_097;M_098;M_099;M_100;M_101;M_102;M_103;M_104;M_105;M_106;M_107;M_108;M_087;M_088
+Malignant_22,Unknown,25,M_234;M_235;M_236;M_237;M_238;M_239;M_240;M_241;M_242;M_243;M_244;M_245;M_246;M_247;M_248;M_249;M_250;M_251;M_252;M_253;M_254;M_255;M_256;M_257;M_258
+Malignant_23,Unknown,28,M_200;M_201;M_202;M_203;M_204;M_205;M_206;M_207;M_208;M_209;M_210;M_211;M_212;M_213;M_214;M_215;M_216;M_217;M_218;M_219;M_220;M_221;M_222;M_223;M_224;M_225;M_226;M_227
+Malignant_24,Unknown,27,M_280;M_281;M_282;M_283;M_300;M_301;M_302;M_303;M_304;M_305;M_306;M_307;M_308;M_309;M_310;M_311;M_312;M_313;M_314;M_315;M_316;M_317;M_318;M_319;M_320;M_321;M_322
+Malignant_25,Unknown,12,M_288;M_289;M_290;M_291;M_292;M_293;M_294;M_295;M_296;M_297;M_298;M_299
+Malignant_26,Unknown,23,M_008;M_009;M_014;M_015;M_053;M_054;M_055;M_056;M_057;M_144;M_145;M_146;M_147;M_148;M_149;M_150;M_151;M_152;M_153;M_154;M_155;M_156;M_058
+Malignant_27,Unknown,15,M_269;M_270;M_271;M_272;M_273;M_274;M_275;M_276;M_277;M_278;M_279;M_284;M_285;M_286;M_287
+Malignant_28,Unknown,16,M_167;M_168;M_169;M_172;M_173;M_174;M_175;M_176;M_177;M_178;M_179;M_262;M_263;M_266;M_267;M_268
+Malignant_29,Unknown,20,M_037;M_038;M_039;M_040;M_041;M_042;M_043;M_427;M_428;M_429;M_430;M_431;M_432;M_433;M_434;M_435;M_436;M_437;M_438;M_439
+Malignant_3,Unknown,3,M_192;M_193;M_194
+Malignant_30,Unknown,28,M_002;M_003;M_004;M_005;M_006;M_007;M_059;M_060;M_061;M_062;M_063;M_064;M_065;M_066;M_067;M_068;M_069;M_070;M_071;M_072;M_073;M_074;M_075;M_076;M_077;M_078;M_365;M_366
+Malignant_31,Unknown,10,M_010;M_011;M_012;M_013;M_138;M_139;M_140;M_141;M_142;M_143
+Malignant_32,Unknown,13,M_030;M_031;M_032;M_033;M_034;M_035;M_420;M_421;M_422;M_423;M_424;M_425;M_426
+Malignant_33,Unknown,22,M_157;M_158;M_159;M_160;M_161;M_162;M_163;M_164;M_165;M_166;M_170;M_171;M_259;M_260;M_261;M_264;M_265;M_352;M_353;M_354;M_355;M_359
+Malignant_34,Unknown,21,M_381;M_382;M_383;M_384;M_385;M_386;M_387;M_388;M_389;M_390;M_391;M_392;M_393;M_394;M_395;M_396;M_397;M_398;M_399;M_400;M_401
+Malignant_35,Unknown,13,M_402;M_403;M_404;M_405;M_406;M_407;M_408;M_409;M_410;M_411;M_412;M_413;M_414
+Malignant_36,Unknown,21,M_019;M_020;M_021;M_022;M_027;M_028;M_462;M_463;M_464;M_465;M_466;M_467;M_468;M_469;M_470;M_471;M_472;M_473;M_474;M_538;M_542
+Malignant_37,Unknown,12,M_023;M_024;M_025;M_026;M_475;M_476;M_477;M_478;M_479;M_480;M_481;M_482
+Malignant_38,Unknown,10,M_122;M_123;M_124;M_125;M_126;M_127;M_128;M_129;M_130;M_131
+Malignant_39,Unknown,18,M_110;M_111;M_112;M_113;M_114;M_115;M_116;M_117;M_118;M_119;M_120;M_121;M_132;M_133;M_134;M_135;M_136;M_137
+Malignant_4,Unknown,5,M_195;M_196;M_197;M_198;M_199
+Malignant_5,Unknown,6,M_231;M_232;M_233;M_230;M_229;M_228
+Malignant_6,Unknown,3,M_323;M_324;M_325
+Malignant_7,Unknown,27,M_327;M_328;M_329;M_330;M_331;M_332;M_333;M_334;M_335;M_336;M_337;M_338;M_339;M_340;M_341;M_342;M_343;M_344;M_345;M_346;M_347;M_348;M_349;M_350;M_351;M_029;M_326
+Malignant_8,Unknown,10,M_036;M_356;M_357;M_358;M_360;M_361;M_362;M_363;M_364;M_109
+Malignant_9,Unknown,14,M_367;M_368;M_369;M_370;M_371;M_372;M_373;M_374;M_375;M_376;M_377;M_378;M_379;M_380
+Normal_0,Unknown,9,N_008;N_009;N_010;N_011;N_012;N_013;N_412;N_048;N_049
+Normal_1,Unknown,3,N_112;N_113;N_114
+Normal_10,Unknown,4,N_021;N_022;N_023;N_024
+Normal_11,Unknown,8,N_210;N_313;N_314;N_315;N_316;N_319;N_320;N_321
+Normal_12,Unknown,3,N_220;N_221;N_222
+Normal_13,Unknown,8,N_223;N_224;N_225;N_226;N_227;N_228;N_229;N_344
+Normal_14,Unknown,4,N_025;N_026;N_031;N_032
+Normal_15,Unknown,4,N_027;N_028;N_029;N_030
+Normal_16,Unknown,9,N_270;N_271;N_272;N_273;N_274;N_275;N_276;N_277;N_278
+Normal_17,Unknown,6,N_294;N_295;N_296;N_292;N_291;N_293
+Normal_18,Unknown,13,N_324;N_325;N_326;N_327;N_328;N_329;N_330;N_331;N_332;N_333;N_334;N_335;N_411
+Normal_19,Unknown,5,N_033;N_034;N_035;N_036;N_038
+Normal_2,Unknown,9,N_115;N_116;N_117;N_118;N_119;N_120;N_121;N_349;N_348
+Normal_20,Unknown,3,N_336;N_337;N_338
+Normal_21,Unknown,6,N_339;N_340;N_341;N_342;N_343;N_410
+Normal_22,Unknown,3,N_345;N_346;N_347
+Normal_23,Unknown,5,N_351;N_352;N_353;N_354;N_350
+Normal_24,Unknown,3,N_356;N_357;N_355
+Normal_25,Unknown,5,N_039;N_040;N_041;N_042;N_043
+Normal_26,Unknown,5,N_044;N_045;N_046;N_047;N_037
+Normal_27,Unknown,9,N_050;N_051;N_052;N_053;N_054;N_055;N_056;N_057;N_058
+Normal_28,Unknown,5,N_059;N_060;N_061;N_062;N_063
+Normal_29,Unknown,7,N_081;N_082;N_083;N_084;N_085;N_086;N_087
+Normal_3,Unknown,9,N_122;N_123;N_124;N_125;N_126;N_127;N_128;N_129;N_130
+Normal_30,Unknown,12,N_135;N_136;N_137;N_138;N_139;N_140;N_141;N_142;N_143;N_144;N_145;N_415
+Normal_31,Unknown,14,N_231;N_232;N_233;N_234;N_235;N_236;N_237;N_238;N_239;N_240;N_241;N_242;N_398;N_385
+Normal_32,Unknown,12,N_386;N_387;N_388;N_389;N_390;N_391;N_392;N_393;N_394;N_395;N_396;N_397
+Normal_33,Unknown,13,N_211;N_212;N_213;N_214;N_215;N_216;N_217;N_218;N_219;N_317;N_318;N_322;N_323
+Normal_34,Unknown,10,N_244;N_245;N_246;N_247;N_248;N_249;N_250;N_251;N_252;N_253
+Normal_35,Unknown,11,N_190;N_191;N_280;N_281;N_282;N_283;N_284;N_285;N_286;N_287;N_288
+Normal_36,Unknown,11,N_186;N_187;N_188;N_189;N_192;N_193;N_279;N_289;N_290;N_169;N_168
+Normal_37,Unknown,8,N_073;N_074;N_075;N_076;N_077;N_078;N_079;N_080
+Normal_38,Unknown,10,N_088;N_089;N_090;N_091;N_092;N_093;N_094;N_095;N_096;N_097
+Normal_39,Unknown,8,N_254;N_255;N_256;N_257;N_258;N_259;N_260;N_261
+Normal_4,Unknown,11,N_132;N_133;N_134;N_399;N_400;N_401;N_402;N_403;N_404;N_405;N_414
+Normal_40,Unknown,8,N_262;N_263;N_264;N_265;N_266;N_267;N_268;N_269
+Normal_41,Unknown,9,N_297;N_298;N_307;N_308;N_309;N_310;N_311;N_312;N_243
+Normal_42,Unknown,9,N_230;N_299;N_300;N_301;N_302;N_303;N_304;N_305;N_306
+Normal_43,Unknown,8,N_365;N_366;N_367;N_368;N_369;N_370;N_407;N_408
+Normal_44,Unknown,8,N_358;N_359;N_360;N_361;N_362;N_363;N_364;N_406
+Normal_45,Unknown,9,N_064;N_065;N_066;N_067;N_068;N_069;N_070;N_071;N_072
+Normal_46,Unknown,6,N_194;N_195;N_196;N_197;N_198;N_199
+Normal_47,Unknown,8,N_371;N_372;N_373;N_374;N_375;N_376;N_377;N_409
+Normal_48,Unknown,7,N_378;N_379;N_380;N_381;N_382;N_383;N_384
+Normal_49,Unknown,4,N_001;N_006;N_007;N_206
+Normal_5,Unknown,8,N_014;N_015;N_016;N_017;N_018;N_019;N_020;N_413
+Normal_50,Unknown,10,N_002;N_003;N_004;N_005;N_203;N_204;N_205;N_207;N_208;N_209
+Normal_51,Unknown,9,N_098;N_099;N_100;N_101;N_102;N_103;N_104;N_105;N_106
+Normal_52,Unknown,5,N_108;N_109;N_110;N_111;N_107
+Normal_53,Unknown,6,N_131;N_153;N_154;N_155;N_156;N_157
+Normal_54,Unknown,8,N_146;N_147;N_148;N_149;N_150;N_151;N_152;N_416
+Normal_6,Unknown,10,N_158;N_159;N_160;N_161;N_162;N_163;N_164;N_165;N_166;N_167
+Normal_7,Unknown,4,N_170;N_171;N_172;N_173
+Normal_8,Unknown,12,N_174;N_175;N_176;N_177;N_178;N_179;N_180;N_181;N_182;N_183;N_184;N_185
+Normal_9,Unknown,3,N_200;N_201;N_202

results/task06_clustering_validation.json

@@ -0,0 +1,151 @@
+{
+  "dataset": "Task06_Lung",
+  "n_patients": 63,
+  "n_slices": 3540,
+  "slices_per_patient": 10,
+  "all_slices": true,
+  "random_slices": true,
+  "seed": 42,
+  "methods": {
+    "thumbnail_64x64": {
+      "ARI": 0.29796172449326236,
+      "NMI": 0.6762977925745551,
+      "homogeneity": 0.6749507693602819,
+      "completeness": 0.6776502031446745,
+      "v_measure": 0.6762977925745552,
+      "cluster_purity": {
+        "median": 0.006779661016949152,
+        "mean": 0.007775087436104384,
+        "min": 0.003672316384180791,
+        "max": 0.01694915254237288,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.005649717514124294,
+        "mean": 0.006770693211371176,
+        "min": 0.0019774011299435027
+      }
+    },
+    "features_VGG16": {
+      "ARI": 0.25752977845774433,
+      "NMI": 0.6387724510004029,
+      "homogeneity": 0.6346040991748216,
+      "completeness": 0.6429959239291733,
+      "v_measure": 0.638772451000403,
+      "cluster_purity": {
+        "median": 0.0064971751412429375,
+        "mean": 0.007286342032104744,
+        "min": 0.002824858757062147,
+        "max": 0.014124293785310734,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.005084745762711864,
+        "mean": 0.006250560487848622,
+        "min": 0.0022598870056497176
+      }
+    },
+    "features_DenseNet121": {
+      "ARI": 0.2830074213473864,
+      "NMI": 0.6399794747034858,
+      "homogeneity": 0.6373804978038539,
+      "completeness": 0.6425997335109946,
+      "v_measure": 0.6399794747034859,
+      "cluster_purity": {
+        "median": 0.0064971751412429375,
+        "mean": 0.007447762532508295,
+        "min": 0.0019774011299435027,
+        "max": 0.017796610169491526,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.005649717514124294,
+        "mean": 0.006775177114160163,
+        "min": 0.0022598870056497176
+      }
+    },
+    "features_EfficientNetB1": {
+      "ARI": 0.24518407277228452,
+      "NMI": 0.6023771039464246,
+      "homogeneity": 0.6013989746958669,
+      "completeness": 0.6033584200844399,
+      "v_measure": 0.6023771039464246,
+      "cluster_purity": {
+        "median": 0.0064971751412429375,
+        "mean": 0.006918662003407767,
+        "min": 0.002542372881355932,
+        "max": 0.016101694915254237,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.004519774011299435,
+        "mean": 0.005851493139628732,
+        "min": 0.0019774011299435027
+      }
+    },
+    "features_MobileNetV2": {
+      "ARI": 0.22930407550375523,
+      "NMI": 0.5962965082158278,
+      "homogeneity": 0.5911426275134923,
+      "completeness": 0.6015410476089263,
+      "v_measure": 0.5962965082158279,
+      "cluster_purity": {
+        "median": 0.0062146892655367235,
+        "mean": 0.00685588736436194,
+        "min": 0.0019774011299435027,
+        "max": 0.01638418079096045,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.004519774011299435,
+        "mean": 0.006071204376289122,
+        "min": 0.001694915254237288
+      }
+    },
+    "features_ResNet50": {
+      "ARI": 0.31425357943429666,
+      "NMI": 0.6679639213211583,
+      "homogeneity": 0.666426843431343,
+      "completeness": 0.6695081059783261,
+      "v_measure": 0.6679639213211583,
+      "cluster_purity": {
+        "median": 0.0064971751412429375,
+        "mean": 0.007887185005829073,
+        "min": 0.002542372881355932,
+        "max": 0.01977401129943503,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.005932203389830509,
+        "mean": 0.007174244462380055,
+        "min": 0.0019774011299435027
+      }
+    },
+    "random_baseline": {
+      "ARI": 0.00022828855339970057,
+      "NMI": 0.14819729551651747,
+      "homogeneity": 0.14928435888387095,
+      "completeness": 0.1471259493213095,
+      "v_measure": 0.1481972955165175,
+      "cluster_purity": {
+        "median": 0.0011299435028248588,
+        "mean": 0.0011433952111918214,
+        "min": 0.000847457627118644,
+        "max": 0.001694915254237288,
+        "frac_gt_70": 0.0,
+        "frac_gt_90": 0.0
+      },
+      "dominant_capture": {
+        "median": 0.0011299435028248588,
+        "mean": 0.001040265447045108,
+        "min": 0.0005649717514124294
+      }
+    }
+  }
+}

results/thumbnail_patient_manifest.csv

@@ -0,0 +1,111 @@
+patient_id,class,n_images,images
+bengin_0,bengin,7,B_001;B_002;B_003;B_066;B_067;B_074;B_075
+bengin_1,bengin,12,B_010;B_011;B_012;B_013;B_042;B_043;B_044;B_045;B_053;B_054;B_055;B_009
+bengin_10,bengin,9,B_004;B_005;B_056;B_057;B_058;B_059;B_006;B_007;B_008
+bengin_11,bengin,7,B_046;B_047;B_048;B_049;B_050;B_051;B_052
+bengin_12,bengin,8,B_060;B_061;B_062;B_063;B_064;B_065;B_092;B_093
+bengin_13,bengin,16,B_068;B_069;B_070;B_071;B_072;B_073;B_076;B_077;B_078;B_079;B_080;B_081;B_082;B_083;B_084;B_085
+bengin_14,bengin,6,B_086;B_087;B_088;B_089;B_090;B_091
+bengin_2,bengin,7,B_100;B_094;B_095;B_096;B_097;B_098;B_099
+bengin_3,bengin,15,B_101;B_102;B_103;B_104;B_105;B_106;B_107;B_108;B_109;B_110;B_030;B_031;B_032;B_033;B_034
+bengin_4,bengin,6,B_111;B_112;B_113;B_114;B_115;B_116
+bengin_5,bengin,4,B_117;B_118;B_119;B_120
+bengin_6,bengin,8,B_014;B_015;B_016;B_017;B_018;B_019;B_020;B_035
+bengin_7,bengin,5,B_021;B_022;B_023;B_024;B_025
+bengin_8,bengin,4,B_026;B_027;B_028;B_029
+bengin_9,bengin,6,B_036;B_037;B_038;B_039;B_040;B_041
+malignant_0,malignant,30,M_001;M_017;M_018;M_352;M_353;M_354;M_355;M_359;M_539;M_540;M_541;M_543;M_544;M_545;M_546;M_547;M_548;M_549;M_550;M_551;M_552;M_553;M_554;M_555;M_556;M_557;M_558;M_559;M_560;M_561
+malignant_1,malignant,19,M_010;M_011;M_012;M_013;M_138;M_139;M_014;M_140;M_141;M_142;M_143;M_144;M_145;M_146;M_147;M_148;M_281;M_282;M_283
+malignant_10,malignant,16,M_180;M_508;M_509;M_510;M_511;M_512;M_513;M_514;M_515;M_516;M_517;M_518;M_519;M_520;M_521;M_522
+malignant_11,malignant,9,M_183;M_184;M_185;M_186;M_187;M_188;M_189;M_190;M_191
+malignant_12,malignant,12,M_019;M_020;M_021;M_022;M_023;M_024;M_025;M_026;M_027;M_028;M_538;M_542
+malignant_13,malignant,8,M_192;M_193;M_194;M_195;M_196;M_197;M_198;M_199
+malignant_14,malignant,26,M_002;M_003;M_004;M_005;M_059;M_006;M_060;M_061;M_062;M_063;M_064;M_065;M_066;M_067;M_068;M_069;M_007;M_070;M_071;M_072;M_073;M_074;M_075;M_076;M_077;M_078
+malignant_15,malignant,10,M_200;M_201;M_202;M_203;M_204;M_205;M_206;M_207;M_208;M_209
+malignant_16,malignant,9,M_210;M_211;M_212;M_213;M_214;M_215;M_216;M_217;M_218
+malignant_17,malignant,9,M_219;M_220;M_221;M_222;M_223;M_224;M_225;M_226;M_227
+malignant_18,malignant,6,M_228;M_229;M_230;M_231;M_232;M_233
+malignant_19,malignant,14,M_234;M_235;M_236;M_237;M_238;M_239;M_240;M_241;M_242;M_243;M_244;M_245;M_246;M_247
+malignant_2,malignant,21,M_100;M_101;M_102;M_103;M_104;M_105;M_106;M_107;M_108;M_109;M_089;M_090;M_091;M_092;M_093;M_094;M_095;M_096;M_097;M_098;M_099
+malignant_20,malignant,11,M_248;M_249;M_250;M_251;M_252;M_253;M_254;M_255;M_256;M_257;M_258
+malignant_21,malignant,15,M_262;M_263;M_266;M_267;M_268;M_356;M_357;M_358;M_036;M_360;M_361;M_362;M_363;M_364;M_527
+malignant_22,malignant,23,M_280;M_301;M_302;M_303;M_304;M_305;M_306;M_307;M_308;M_309;M_310;M_311;M_312;M_313;M_314;M_315;M_316;M_317;M_318;M_319;M_320;M_321;M_322
+malignant_23,malignant,13,M_288;M_289;M_290;M_291;M_292;M_293;M_294;M_295;M_296;M_297;M_298;M_299;M_300
+malignant_24,malignant,15,M_029;M_323;M_324;M_325;M_326;M_079;M_080;M_081;M_082;M_083;M_084;M_085;M_086;M_087;M_088
+malignant_25,malignant,13,M_030;M_031;M_032;M_033;M_034;M_035;M_037;M_038;M_039;M_040;M_041;M_042;M_043
+malignant_26,malignant,25,M_327;M_328;M_329;M_330;M_331;M_332;M_333;M_334;M_335;M_336;M_337;M_338;M_339;M_340;M_341;M_342;M_343;M_344;M_345;M_346;M_347;M_348;M_349;M_350;M_351
+malignant_27,malignant,8,M_365;M_366;M_053;M_054;M_055;M_056;M_057;M_058
+malignant_28,malignant,14,M_367;M_368;M_369;M_370;M_371;M_372;M_373;M_374;M_375;M_376;M_377;M_378;M_379;M_380
+malignant_29,malignant,13,M_382;M_383;M_384;M_385;M_386;M_387;M_388;M_389;M_390;M_391;M_392;M_393;M_394
+malignant_3,malignant,12,M_110;M_111;M_112;M_113;M_114;M_115;M_116;M_117;M_118;M_119;M_120;M_121
+malignant_30,malignant,20,M_395;M_396;M_397;M_398;M_399;M_400;M_401;M_402;M_403;M_404;M_405;M_406;M_407;M_408;M_409;M_410;M_411;M_412;M_413;M_414
+malignant_31,malignant,12,M_415;M_416;M_417;M_418;M_419;M_420;M_421;M_422;M_423;M_424;M_425;M_426
+malignant_32,malignant,13,M_427;M_428;M_429;M_430;M_431;M_432;M_433;M_434;M_435;M_436;M_437;M_438;M_439
+malignant_33,malignant,9,M_044;M_045;M_046;M_047;M_048;M_049;M_050;M_051;M_052
+malignant_34,malignant,5,M_440;M_441;M_442;M_443;M_444
+malignant_35,malignant,4,M_445;M_446;M_447;M_448
+malignant_36,malignant,10,M_449;M_450;M_451;M_452;M_453;M_454;M_455;M_456;M_457;M_458
+malignant_37,malignant,24,M_459;M_460;M_461;M_462;M_463;M_464;M_465;M_466;M_467;M_468;M_469;M_470;M_471;M_472;M_473;M_474;M_475;M_476;M_477;M_478;M_479;M_480;M_481;M_482
+malignant_38,malignant,11,M_483;M_484;M_485;M_486;M_487;M_488;M_489;M_490;M_491;M_492;M_493
+malignant_39,malignant,14,M_494;M_495;M_496;M_497;M_498;M_499;M_500;M_501;M_502;M_503;M_504;M_505;M_506;M_507
+malignant_4,malignant,9,M_122;M_123;M_124;M_125;M_126;M_127;M_128;M_129;M_381
+malignant_5,malignant,8,M_130;M_131;M_132;M_133;M_134;M_135;M_136;M_137
+malignant_6,malignant,26,M_149;M_015;M_150;M_151;M_152;M_153;M_154;M_155;M_156;M_269;M_270;M_271;M_272;M_273;M_274;M_275;M_276;M_277;M_278;M_279;M_284;M_285;M_286;M_287;M_008;M_009
+malignant_7,malignant,17,M_157;M_158;M_159;M_160;M_161;M_162;M_163;M_164;M_165;M_166;M_170;M_171;M_259;M_260;M_261;M_264;M_265
+malignant_8,malignant,17,M_016;M_181;M_182;M_523;M_524;M_525;M_526;M_528;M_529;M_530;M_531;M_532;M_533;M_534;M_535;M_536;M_537
+malignant_9,malignant,11,M_167;M_168;M_169;M_172;M_173;M_174;M_175;M_176;M_177;M_178;M_179
+normal_0,Normal,7,N_001;N_002;N_003;N_004;N_005;N_006;N_007
+normal_1,Normal,7,N_010;N_011;N_012;N_013;N_412;N_008;N_009
+normal_10,Normal,8,N_014;N_015;N_016;N_017;N_018;N_019;N_020;N_413
+normal_11,Normal,13,N_146;N_147;N_148;N_149;N_150;N_151;N_152;N_153;N_154;N_155;N_156;N_157;N_416
+normal_12,Normal,10,N_158;N_159;N_160;N_161;N_162;N_163;N_164;N_165;N_166;N_167
+normal_13,Normal,6,N_168;N_169;N_170;N_171;N_172;N_173
+normal_14,Normal,5,N_174;N_182;N_183;N_184;N_185
+normal_15,Normal,7,N_175;N_176;N_177;N_178;N_179;N_180;N_181
+normal_16,Normal,5,N_186;N_187;N_188;N_189;N_194
+normal_17,Normal,9,N_190;N_191;N_192;N_193;N_195;N_196;N_197;N_198;N_199
+normal_18,Normal,4,N_200;N_201;N_202;N_210
+normal_19,Normal,16,N_203;N_204;N_205;N_206;N_207;N_208;N_209;N_211;N_212;N_213;N_214;N_215;N_216;N_217;N_218;N_219
+normal_2,Normal,6,N_100;N_105;N_106;N_107;N_098;N_099
+normal_20,Normal,4,N_021;N_022;N_023;N_024
+normal_21,Normal,3,N_220;N_221;N_222
+normal_22,Normal,7,N_223;N_224;N_225;N_226;N_227;N_228;N_229
+normal_23,Normal,3,N_230;N_243;N_385
+normal_24,Normal,7,N_231;N_232;N_233;N_234;N_235;N_236;N_237
+normal_25,Normal,15,N_238;N_239;N_240;N_241;N_242;N_244;N_245;N_246;N_247;N_248;N_249;N_250;N_251;N_252;N_253
+normal_26,Normal,4,N_025;N_026;N_031;N_032
+normal_27,Normal,14,N_254;N_255;N_256;N_257;N_258;N_259;N_260;N_261;N_262;N_263;N_264;N_270;N_271;N_272
+normal_28,Normal,11,N_265;N_266;N_267;N_268;N_269;N_273;N_274;N_275;N_276;N_277;N_278
+normal_29,Normal,4,N_027;N_028;N_029;N_030
+normal_3,Normal,8,N_101;N_102;N_103;N_104;N_108;N_109;N_110;N_111
+normal_30,Normal,7,N_279;N_291;N_292;N_293;N_064;N_065;N_066
+normal_31,Normal,14,N_280;N_281;N_282;N_283;N_284;N_285;N_286;N_287;N_288;N_289;N_290;N_294;N_295;N_296
+normal_32,Normal,8,N_297;N_298;N_307;N_308;N_309;N_310;N_311;N_312
+normal_33,Normal,9,N_299;N_300;N_301;N_302;N_303;N_304;N_305;N_306;N_355
+normal_34,Normal,7,N_313;N_314;N_315;N_316;N_319;N_320;N_321
+normal_35,Normal,4,N_317;N_318;N_322;N_323
+normal_36,Normal,13,N_324;N_325;N_326;N_327;N_328;N_329;N_330;N_331;N_332;N_333;N_334;N_335;N_411
+normal_37,Normal,7,N_033;N_034;N_035;N_356;N_357;N_036;N_038
+normal_38,Normal,9,N_336;N_337;N_338;N_339;N_340;N_341;N_342;N_343;N_410
+normal_39,Normal,1,N_344
+normal_4,Normal,3,N_112;N_113;N_114
+normal_40,Normal,4,N_345;N_346;N_347;N_350
+normal_41,Normal,6,N_348;N_349;N_351;N_352;N_353;N_354
+normal_42,Normal,16,N_358;N_359;N_360;N_361;N_362;N_363;N_364;N_365;N_366;N_367;N_368;N_369;N_370;N_406;N_407;N_408
+normal_43,Normal,5,N_037;N_044;N_045;N_046;N_047
+normal_44,Normal,15,N_371;N_372;N_373;N_374;N_375;N_376;N_377;N_378;N_379;N_380;N_381;N_382;N_383;N_384;N_409
+normal_45,Normal,13,N_386;N_387;N_388;N_389;N_390;N_391;N_392;N_393;N_394;N_395;N_396;N_397;N_398
+normal_46,Normal,5,N_039;N_040;N_041;N_042;N_043
+normal_47,Normal,7,N_399;N_400;N_401;N_402;N_403;N_404;N_405
+normal_48,Normal,2,N_048;N_049
+normal_49,Normal,9,N_050;N_051;N_052;N_053;N_054;N_055;N_056;N_057;N_058
+normal_5,Normal,7,N_115;N_116;N_117;N_118;N_119;N_120;N_121
+normal_50,Normal,11,N_059;N_062;N_063;N_073;N_074;N_075;N_076;N_077;N_078;N_079;N_080
+normal_51,Normal,6,N_060;N_061;N_088;N_089;N_090;N_091
+normal_52,Normal,6,N_067;N_068;N_069;N_070;N_071;N_072
+normal_53,Normal,7,N_081;N_082;N_083;N_084;N_085;N_086;N_087
+normal_54,Normal,6,N_092;N_093;N_094;N_095;N_096;N_097
+normal_6,Normal,4,N_122;N_123;N_124;N_125
+normal_7,Normal,5,N_126;N_127;N_128;N_129;N_130
+normal_8,Normal,5,N_131;N_132;N_133;N_134;N_414
+normal_9,Normal,12,N_135;N_136;N_137;N_138;N_139;N_140;N_141;N_142;N_143;N_144;N_145;N_415

scripts/analysis/compare_clustering_methods.py

@@ -0,0 +1,245 @@
+#!/usr/bin/env python3
+"""
+compare_clustering_methods.py
+
+Compare two patient-clustering approaches against the v8 "reference":
+  1. Thumbnail-based (64×64 grayscale K-means — the manuscript's method)
+  2. Feature-based (PCA-50d of VGG16 features → K-means — our method)
+
+Metrics (all label-invariant):
+  - Adjusted Rand Index (ARI)
+  - Normalized Mutual Information (NMI)
+  - V-measure (homogeneity + completeness)
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/compare_clustering_methods.py
+"""
+
+import os, sys, csv, re
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+
+from sklearn.metrics import (
+    adjusted_rand_score,
+    normalized_mutual_info_score,
+    homogeneity_completeness_v_measure,
+)
+from sklearn.decomposition import PCA
+from sklearn.cluster import KMeans
+from PIL import Image
+
+# ---------------------------------------------------------------------------
+# Config
+# ---------------------------------------------------------------------------
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+DATASET_PATH = os.path.expanduser("~/Documents/data_leakage/The IQ-OTHNCCD lung cancer dataset")
+MANIFEST_V8 = os.path.join(ROOT, ".archive", "results", "patient_manifest_v8.csv")
+MANIFEST_SIMPLE = os.path.join(ROOT, "results", "simple_patient_manifest.csv")
+FEATURES_PATH = os.path.join(ROOT, "features", "VGG16_features.npz")
+
+PATIENT_COUNTS = {
+    "Bengin cases": 15,
+    "Malignant cases": 40,
+    "Normal cases": 55,
+}
+
+CLASS_ORDER = ["Bengin cases", "Malignant cases", "Normal cases"]
+
+SEED = 42
+THUMBNAIL_SIZE = (64, 64)
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+
+def f2n(fname):
+    m = re.search(r"\((\d+)\)", fname)
+    num = int(m.group(1)) if m else None
+    for cls_key, prefix in [
+        ("Bengin cases", "B"),
+        ("Malignant cases", "M"),
+        ("Normal cases", "N"),
+    ]:
+        if fname.startswith(cls_key.rstrip("s")):
+            return f"{prefix}_{num:03d}" if num else fname
+    return fname
+
+
+def load_v8_assignments(path):
+    """Load v8 manifest, return {image_name: patient_id}."""
+    mapping = {}
+    with open(path, newline="") as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            for img in row["confirmed_images"].split(";"):
+                if img:
+                    mapping[img] = row["patient_id"]
+    return mapping
+
+
+def load_simple_assignments(path):
+    """Load simple feature-based manifest, return {image_name: patient_id}."""
+    mapping = {}
+    with open(path, newline="") as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            for img in row["images"].split(";"):
+                if img:
+                    mapping[img] = row["patient_id"]
+    return mapping
+
+
+def run_thumbnail_clustering(dataset_path, filenames, labels):
+    """Run 64×64 grayscale thumbnail K-means (manuscript method)."""
+    groups = np.array([None] * len(labels), dtype=object)
+
+    for class_name, n_clusters in PATIENT_COUNTS.items():
+        idx_class = np.where(labels == class_name)[0]
+
+        # Load thumbnails
+        X = []
+        for i in idx_class:
+            fname = filenames[i]
+            img_path = os.path.join(dataset_path, class_name, fname)
+            try:
+                img = Image.open(img_path).convert("L")
+                img = img.resize(THUMBNAIL_SIZE)
+                img_arr = np.array(img, dtype=np.float32) / 255.0
+                X.append(img_arr.flatten())
+            except Exception:
+                X.append(np.zeros(64 * 64, dtype=np.float32))
+        X = np.array(X, dtype=np.float32)
+
+        # PCA → 50d
+        n_pca = min(50, X.shape[0] - 1, X.shape[1])
+        pca = PCA(n_components=n_pca, random_state=SEED)
+        X_pca = pca.fit_transform(X)
+
+        # K-means
+        kmeans = KMeans(n_clusters=n_clusters, random_state=SEED, n_init=20)
+        clusters = kmeans.fit_predict(X_pca)
+
+        for i, cluster_id in zip(idx_class, clusters):
+            groups[i] = f"{class_name}_cluster_{cluster_id}"
+
+        sizes = np.bincount(clusters)
+        print(f"  {class_name}: k={n_clusters}, sizes min={sizes.min()} max={sizes.max()} mean={sizes.mean():.1f}")
+
+    return groups
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+
+print("Loading VGG16 features for filename/label reference...")
+data = np.load(FEATURES_PATH, allow_pickle=True)
+all_filenames = data["filenames"]
+all_labels = data["Y"]
+X_vgg16 = data["X"]
+
+print(f"  {len(all_filenames)} images across {len(np.unique(all_labels))} classes")
+
+# Load reference (v8) and our feature-based assignments
+v8_map = load_v8_assignments(MANIFEST_V8)
+simple_map = load_simple_assignments(MANIFEST_SIMPLE)
+
+# Build per-image label arrays for all three methods, aligned by filename
+# We need all images that exist in ALL three
+v8_labels_list = []
+simple_labels_list = []
+thumb_labels_list = []  # filled after clustering
+common_filenames = []
+common_labels = []
+
+# First, run thumbnail clustering
+print("\nRunning thumbnail-based K-means (manuscript method)...")
+thumb_groups = run_thumbnail_clustering(DATASET_PATH, all_filenames, all_labels)
+
+# Build thumbnail mapping (using short names like B_001)
+thumb_map = {}
+for fname, group in zip(all_filenames, thumb_groups):
+    short = f2n(fname)
+    for cls in CLASS_ORDER:
+        if group.startswith(cls):
+            cluster_id = int(group.split("_cluster_")[-1])
+            prefix = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"}[cls]
+            thumb_map[short] = f"thumb_{prefix}_{cluster_id:02d}"
+            break
+
+# Filter to images present in all three
+# Feature filenames are like "Bengin case (1).jpg", manifests use "B_001"
+for fname, label in zip(all_filenames, all_labels):
+    short_name = f2n(fname)
+    v8_id = v8_map.get(short_name)
+    simple_id = simple_map.get(short_name)
+    thumb_id = thumb_map.get(short_name)
+    if v8_id and simple_id and thumb_id:
+        common_filenames.append(short_name)
+        common_labels.append(label)
+        v8_labels_list.append(v8_id)
+        simple_labels_list.append(simple_id)
+        thumb_labels_list.append(thumb_id)
+
+print(f"\nImages common to all three methods: {len(common_filenames)}")
+
+# Convert to numpy arrays
+v8_labels_arr = np.array(v8_labels_list)
+simple_labels_arr = np.array(simple_labels_list)
+thumb_labels_arr = np.array(thumb_labels_list)
+common_labels_arr = np.array(common_labels)
+
+# ---------------------------------------------------------------------------
+# Compute agreement metrics — per-class and overall
+# ---------------------------------------------------------------------------
+
+
+def compute_metrics(ref, pred, name):
+    """Compute clustering agreement metrics against reference."""
+    ari = adjusted_rand_score(ref, pred)
+    nmi = normalized_mutual_info_score(ref, pred)
+    h, c, v = homogeneity_completeness_v_measure(ref, pred)
+    return {"name": name, "ARI": ari, "NMI": nmi, "Homogeneity": h, "Completeness": c, "V_measure": v}
+
+
+print("\n" + "=" * 80)
+print("OVERALL AGREEMENT WITH v8 REFERENCE")
+print("=" * 80)
+
+results = []
+for name, pred in [("Thumbnail (manuscript)", thumb_labels_arr), ("Feature-based (ours)", simple_labels_arr)]:
+    r = compute_metrics(v8_labels_arr, pred, name)
+    results.append(r)
+
+print(f"\n{'Method':<30s} {'ARI':>8s} {'NMI':>8s} {'Homog':>8s} {'Compl':>8s} {'V_meas':>8s}")
+print("-" * 72)
+for r in results:
+    print(f"{r['name']:<30s} {r['ARI']:>8.4f} {r['NMI']:>8.4f} {r['Homogeneity']:>8.4f} {r['Completeness']:>8.4f} {r['V_measure']:>8.4f}")
+
+# Per-class breakdown
+print("\n" + "=" * 80)
+print("PER-CLASS ARI WITH v8 REFERENCE")
+print("=" * 80)
+print(f"\n{'Class':<20s} {'Thumbnail':>10s} {'Feature-based':>15s}")
+print("-" * 47)
+for cls in CLASS_ORDER:
+    mask = common_labels_arr == cls
+    if mask.sum() < 2:
+        continue
+    thumb_ari = adjusted_rand_score(v8_labels_arr[mask], thumb_labels_arr[mask])
+    feat_ari = adjusted_rand_score(v8_labels_arr[mask], simple_labels_arr[mask])
+    better = " ←" if feat_ari > thumb_ari else ""
+    print(f"{cls:<20s} {thumb_ari:>10.4f} {feat_ari:>15.4f}{better}")
+
+# Also: direct agreement between thumbnail and feature-based
+print("\n" + "=" * 80)
+print("THUMBNAIL vs FEATURE-BASED (direct agreement)")
+print("=" * 80)
+direct = compute_metrics(thumb_labels_arr, simple_labels_arr, "Thumb vs Feature")
+print(f"  ARI={direct['ARI']:.4f}  NMI={direct['NMI']:.4f}  V_measure={direct['V_measure']:.4f}")
+
+print("\nDONE")

scripts/analysis/compare_patient_groupings.py

@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+"""compare_patient_groupings.py
+
+Compare the three estimated patient groupings of the IQ-OTH/NCCD
+("lung_effnet") dataset against each other:
+
+    siamese    -> results/siamese_manifest.csv
+    pca50      -> results/simple_patient_manifest.csv   (PCA-50 CNN-feature K-means)
+    thumbnail  -> results/thumbnail_patient_manifest.csv (64x64 grayscale K-means)
+
+Downstream classification accuracy is affected similarly by all three, so the
+question this answers is: do the three methods actually partition the images
+differently? If they agree closely, the method choice is cosmetic; if they
+diverge, the groupings are genuinely different partitions — which, combined with
+the Task06 ground-truth validation (where siamese scored higher purity/capture),
+is what makes the siamese work worthwhile.
+
+Metrics are all label-invariant (cluster-id names don't matter):
+    ARI  - Adjusted Rand Index (chance-corrected pair agreement)
+    NMI  - Normalized Mutual Information
+    V    - V-measure (harmonic mean of homogeneity & completeness)
+
+Outputs:
+    results/compare_patient_groupings.json
+    plots/analysis/grouping_agreement_heatmap.png
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/analysis/compare_patient_groupings.py
+"""
+
+import os
+import csv
+import json
+from itertools import combinations
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.metrics import (
+    adjusted_rand_score,
+    normalized_mutual_info_score,
+    homogeneity_completeness_v_measure,
+)
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+RESULTS_DIR = os.path.join(ROOT, "results")
+PLOTS_DIR = os.path.join(ROOT, "plots", "analysis")
+
+METHODS = {
+    "siamese": "siamese_manifest.csv",
+    "pca50": "simple_patient_manifest.csv",
+    "thumbnail": "thumbnail_patient_manifest.csv",
+}
+METHOD_ORDER = ["siamese", "pca50", "thumbnail"]
+CLASS_OF = {"B": "Benign", "M": "Malignant", "N": "Normal"}
+
+
+def load_manifest(path):
+    """Return {image_short_name: group_id} from a manifest CSV."""
+    mapping = {}
+    with open(path, newline="") as f:
+        for row in csv.DictReader(f):
+            for img in row["images"].split(";"):
+                if img:
+                    mapping[img] = row["patient_id"]
+    return mapping
+
+
+def size_stats(labels):
+    """Group-size distribution for one method's label array."""
+    _, counts = np.unique(labels, return_counts=True)
+    return {
+        "n_groups": int(len(counts)),
+        "min": int(counts.min()),
+        "median": float(np.median(counts)),
+        "mean": float(counts.mean()),
+        "max": int(counts.max()),
+        "singletons": int((counts == 1).sum()),
+    }
+
+
+def agreement(a, b):
+    h, c, v = homogeneity_completeness_v_measure(a, b)
+    return {"ARI": float(adjusted_rand_score(a, b)),
+            "NMI": float(normalized_mutual_info_score(a, b)),
+            "V": float(v)}
+
+
+def main():
+    maps = {m: load_manifest(os.path.join(RESULTS_DIR, f))
+            for m, f in METHODS.items()}
+
+    # Align on images present in all three (should be the full 1,097).
+    common = sorted(set.intersection(*[set(mp) for mp in maps.values()]))
+    classes = np.array([CLASS_OF.get(img.split("_")[0], "?") for img in common])
+    labels = {m: np.array([maps[m][img] for img in common]) for m in METHOD_ORDER}
+    print(f"Aligned on {len(common)} images common to all three methods.\n")
+
+    # --- Structural summary --------------------------------------------------
+    print("=" * 72)
+    print("GROUP STRUCTURE PER METHOD")
+    print("=" * 72)
+    print(f"{'method':<12s} {'groups':>7s} {'min':>5s} {'median':>7s} "
+          f"{'mean':>6s} {'max':>5s} {'singletons':>11s}")
+    structure = {}
+    for m in METHOD_ORDER:
+        s = size_stats(labels[m])
+        structure[m] = s
+        print(f"{m:<12s} {s['n_groups']:>7d} {s['min']:>5d} {s['median']:>7.1f} "
+              f"{s['mean']:>6.1f} {s['max']:>5d} {s['singletons']:>11d}")
+
+    # --- Pairwise agreement --------------------------------------------------
+    print("\n" + "=" * 72)
+    print("PAIRWISE AGREEMENT (how similarly the methods partition the images)")
+    print("=" * 72)
+    print(f"{'pair':<24s} {'ARI':>8s} {'NMI':>8s} {'V':>8s}")
+    print("-" * 52)
+    pairwise = {}
+    ari_matrix = np.eye(len(METHOD_ORDER))
+    for i, j in combinations(range(len(METHOD_ORDER)), 2):
+        a, b = METHOD_ORDER[i], METHOD_ORDER[j]
+        g = agreement(labels[a], labels[b])
+        pairwise[f"{a}_vs_{b}"] = g
+        ari_matrix[i, j] = ari_matrix[j, i] = g["ARI"]
+        print(f"{a+' vs '+b:<24s} {g['ARI']:>8.3f} {g['NMI']:>8.3f} {g['V']:>8.3f}")
+
+    # --- Per-class ARI -------------------------------------------------------
+    print("\n" + "=" * 72)
+    print("PER-CLASS ARI (agreement within each diagnostic class)")
+    print("=" * 72)
+    print(f"{'pair':<24s} " + " ".join(f"{c:>10s}" for c in ["Benign", "Malignant", "Normal"]))
+    print("-" * 60)
+    per_class = {}
+    for i, j in combinations(range(len(METHOD_ORDER)), 2):
+        a, b = METHOD_ORDER[i], METHOD_ORDER[j]
+        row = {}
+        cells = []
+        for c in ["Benign", "Malignant", "Normal"]:
+            mask = classes == c
+            ari = float(adjusted_rand_score(labels[a][mask], labels[b][mask]))
+            row[c] = ari
+            cells.append(f"{ari:>10.3f}")
+        per_class[f"{a}_vs_{b}"] = row
+        print(f"{a+' vs '+b:<24s} " + " ".join(cells))
+
+    # --- Odd-one-out: mean ARI of each method vs the other two ---------------
+    print("\n" + "=" * 72)
+    print("MEAN ARI OF EACH METHOD VS THE OTHER TWO (lower = most distinct)")
+    print("=" * 72)
+    mean_ari = {}
+    for i, m in enumerate(METHOD_ORDER):
+        others = [ari_matrix[i, j] for j in range(len(METHOD_ORDER)) if j != i]
+        mean_ari[m] = float(np.mean(others))
+        print(f"  {m:<12s} {mean_ari[m]:.3f}")
+    odd = min(mean_ari, key=mean_ari.get)
+    print(f"\n  Most distinct grouping: {odd}")
+
+    # --- Save + heatmap ------------------------------------------------------
+    out = {
+        "n_images": len(common),
+        "structure": structure,
+        "pairwise": pairwise,
+        "per_class_ARI": per_class,
+        "mean_ari_vs_others": mean_ari,
+        "most_distinct": odd,
+    }
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    out_path = os.path.join(RESULTS_DIR, "compare_patient_groupings.json")
+    with open(out_path, "w") as f:
+        json.dump(out, f, indent=2)
+    print(f"\nSaved metrics → {out_path}")
+
+    fig, ax = plt.subplots(figsize=(5.5, 4.5))
+    im = ax.imshow(ari_matrix, vmin=0, vmax=1, cmap="viridis")
+    ax.set_xticks(range(len(METHOD_ORDER)))
+    ax.set_yticks(range(len(METHOD_ORDER)))
+    ax.set_xticklabels(METHOD_ORDER)
+    ax.set_yticklabels(METHOD_ORDER)
+    for i in range(len(METHOD_ORDER)):
+        for j in range(len(METHOD_ORDER)):
+            ax.text(j, i, f"{ari_matrix[i, j]:.2f}", ha="center", va="center",
+                    color="white" if ari_matrix[i, j] < 0.6 else "black", fontsize=11)
+    ax.set_title("Patient-grouping agreement (ARI)\nIQ-OTH/NCCD, 1,097 images")
+    fig.colorbar(im, ax=ax, label="Adjusted Rand Index")
+    plt.tight_layout()
+    os.makedirs(PLOTS_DIR, exist_ok=True)
+    plot_path = os.path.join(PLOTS_DIR, "grouping_agreement_heatmap.png")
+    plt.savefig(plot_path, dpi=150)
+    plt.close()
+    print(f"Saved heatmap → {plot_path}")
+
+
+if __name__ == "__main__":
+    main()

scripts/analysis/verify_feature_dims.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+"""
+verify_feature_dims.py
+
+Quick check: what does TF/Keras actually output for each model
+with include_top=False?  Compare against manuscript claims.
+"""
+
+import numpy as np
+from tensorflow.keras.applications import (
+    VGG16, DenseNet121, EfficientNetB1, MobileNetV2, ResNet50
+)
+
+# Match the manuscript's input sizes
+models = {
+    "VGG16":         (VGG16,         224),
+    "DenseNet121":   (DenseNet121,   224),
+    "EfficientNetB1":(EfficientNetB1, 240),
+    "MobileNetV2":   (MobileNetV2,   224),
+    "ResNet50":      (ResNet50,      224),
+}
+
+manuscript_claims = {
+    "VGG16":         25088,
+    "DenseNet121":   50176,
+    "EfficientNetB1":62720,
+    "MobileNetV2":   62720,
+    "ResNet50":      100352,
+}
+
+print(f"{'Model':<18s} {'Input':>5s} {'Spatial':>10s} {'Flattened':>10s} {'Manuscript':>12s} {'Match?':>7s}")
+print("-" * 70)
+
+for name, (model_fn, input_size) in models.items():
+    model = model_fn(weights="imagenet", include_top=False,
+                     input_shape=(input_size, input_size, 3))
+
+    # Pass a dummy batch through
+    dummy = np.random.randn(1, input_size, input_size, 3)
+    # Need to preprocess correctly — but shape doesn't depend on values
+    output = model.predict(dummy, verbose=0)
+
+    spatial = output.shape[1:4]  # (H, W, C) for channels_last
+    flattened = int(np.prod(spatial))
+    claimed = manuscript_claims[name]
+    match = "✅" if flattened == claimed else "❌"
+
+    print(f"{name:<18s} {input_size:>4}d  {str(spatial):>10s} {flattened:>10d} {claimed:>12d} {match:>7s}")
+
+print("\nNote: TF uses channels_last (NHWC) format.")

scripts/cache/cv_curves.json

@@ -0,0 +1,218 @@
+{
+  "VGG16::image": {
+    "curves": {
+      "32": 0.9624025974025974,
+      "54": 0.9749545454545455,
+      "93": 0.980642857142857,
+      "158": 0.9840649350649351,
+      "271": 0.9886168831168831,
+      "464": 0.9897532467532468,
+      "794": 0.9908961038961038,
+      "1359": 0.9886168831168831,
+      "2326": 0.9897532467532468,
+      "3981": 0.9920324675324675,
+      "6813": 0.987487012987013,
+      "11659": 0.9760714285714286,
+      "19953": 0.9635389610389609
+    },
+    "best_nfeat": 3981,
+    "best_cv": 0.9920324675324675,
+    "best_gamma": 5.411792740597548e-05
+  },
+  "MobileNetV2::image": {
+    "curves": {
+      "32": 0.93612987012987,
+      "54": 0.943,
+      "93": 0.9589675324675324,
+      "158": 0.9692077922077921,
+      "271": 0.9771818181818182,
+      "464": 0.9794740259740259,
+      "794": 0.9829090909090908,
+      "1359": 0.9851818181818182,
+      "2326": 0.9874610389610389,
+      "3981": 0.9886038961038961,
+      "6813": 0.9863246753246753,
+      "11659": 0.979474025974026,
+      "19953": 0.9783441558441558,
+      "34145": 0.9635194805194806,
+      "58434": 0.9566948051948053
+    },
+    "best_nfeat": 3981,
+    "best_cv": 0.9886038961038961,
+    "best_gamma": 5.411792740597548e-05
+  },
+  "DenseNet121::image": {
+    "curves": {
+      "32": 0.9475324675324674,
+      "54": 0.9577987012987013,
+      "93": 0.9760519480519481,
+      "158": 0.9772012987012987,
+      "271": 0.982896103896104,
+      "464": 0.9840454545454547,
+      "794": 0.9840389610389609,
+      "1359": 0.9840454545454544,
+      "2326": 0.9851818181818182,
+      "3981": 0.9863246753246753,
+      "6813": 0.9874675324675325,
+      "11659": 0.9840519480519481,
+      "19953": 0.9829155844155844,
+      "34145": 0.9829090909090908
+    },
+    "best_nfeat": 6813,
+    "best_cv": 0.9874675324675325,
+    "best_gamma": 3.162240848424899e-05
+  },
+  "ResNet50::image": {
+    "curves": {
+      "32": 0.9521168831168831,
+      "54": 0.9646948051948051,
+      "93": 0.9658376623376623,
+      "158": 0.974935064935065,
+      "271": 0.97837012987013,
+      "464": 0.9795,
+      "794": 0.9829155844155844,
+      "1359": 0.9863506493506493,
+      "2326": 0.9908961038961038,
+      "3981": 0.9886168831168831,
+      "6813": 0.987474025974026,
+      "11659": 0.987474025974026,
+      "19953": 0.9863311688311688,
+      "34145": 0.9840584415584415,
+      "58434": 0.9840649350649351,
+      "100000": 0.9829285714285714
+    },
+    "best_nfeat": 2326,
+    "best_cv": 0.9908961038961038,
+    "best_gamma": 9.262401934788839e-05
+  },
+  "EfficientNetB1::image": {
+    "curves": {
+      "32": 0.9373441558441558,
+      "54": 0.9247987012987012,
+      "93": 0.954422077922078,
+      "158": 0.9544415584415585,
+      "271": 0.9669805194805194,
+      "464": 0.9715324675324675,
+      "794": 0.9760649350649351,
+      "1359": 0.9772272727272726,
+      "2326": 0.982896103896104,
+      "3981": 0.9840584415584415,
+      "6813": 0.9863311688311688,
+      "11659": 0.9851818181818182,
+      "19953": 0.9851883116883118,
+      "34145": 0.9806298701298701,
+      "58434": 0.9760649350649351
+    },
+    "best_nfeat": 6813,
+    "best_cv": 0.9863311688311688,
+    "best_gamma": 3.162240848424899e-05
+  },
+  "VGG16::patient": {
+    "curves": {
+      "32": 0.8775928827840215,
+      "54": 0.8869860655394717,
+      "93": 0.9017508373724363,
+      "158": 0.899545406113065,
+      "271": 0.9168027205072213,
+      "464": 0.9111734782922494,
+      "794": 0.8892728094160605,
+      "1359": 0.8776515115789267,
+      "2326": 0.8788736781890476,
+      "3981": 0.8718631943594317,
+      "6813": 0.8729861925854772,
+      "11659": 0.8579623693039702,
+      "19953": 0.8415711387799029
+    },
+    "best_nfeat": 271,
+    "best_cv": 0.9168027205072213,
+    "best_gamma": 0.0001166892125523387
+  },
+  "MobileNetV2::patient": {
+    "curves": {
+      "32": 0.8788294912848766,
+      "54": 0.8777222470719559,
+      "93": 0.8671963468340518,
+      "158": 0.886821936334651,
+      "271": 0.9026897499150754,
+      "464": 0.9027096807901112,
+      "794": 0.9118403567490991,
+      "1359": 0.9084816014962362,
+      "2326": 0.9015949854380466,
+      "3981": 0.897043807427408,
+      "6813": 0.8878722188414819,
+      "11659": 0.8754424416978528,
+      "19953": 0.867454694839946,
+      "34145": 0.8502226180733488,
+      "58434": 0.8421961732955723
+    },
+    "best_nfeat": 794,
+    "best_cv": 0.9118403567490991,
+    "best_gamma": 3.98271745613146e-05
+  },
+  "DenseNet121::patient": {
+    "curves": {
+      "32": 0.8883307016099342,
+      "54": 0.8776358621830977,
+      "93": 0.8695492943753337,
+      "158": 0.9020030888753784,
+      "271": 0.9217021570809798,
+      "464": 0.9203788855354906,
+      "794": 0.9297944483539279,
+      "1359": 0.9321867182439819,
+      "2326": 0.9274615287727703,
+      "3981": 0.9065480925304101,
+      "6813": 0.8846197537489628,
+      "11659": 0.8697672582591949,
+      "19953": 0.865148271598542,
+      "34145": 0.8606789978555252
+    },
+    "best_nfeat": 1359,
+    "best_cv": 0.9321867182439819,
+    "best_gamma": 2.3269151288950545e-05
+  },
+  "ResNet50::patient": {
+    "curves": {
+      "32": 0.8429676740201181,
+      "54": 0.855707644929198,
+      "93": 0.876331816908905,
+      "158": 0.8754247864722885,
+      "271": 0.8913213110995345,
+      "464": 0.8961471621478985,
+      "794": 0.9121510194763156,
+      "1359": 0.9099166863641326,
+      "2326": 0.9086942852731289,
+      "3981": 0.8994412881459878,
+      "6813": 0.8971347519512485,
+      "11659": 0.8834703540651564,
+      "19953": 0.8708504896171615,
+      "34145": 0.8685060975020858,
+      "58434": 0.8604140407100207,
+      "100000": 0.861563465997377
+    },
+    "best_nfeat": 794,
+    "best_cv": 0.9121510194763156,
+    "best_gamma": 3.98271745613146e-05
+  },
+  "EfficientNetB1::patient": {
+    "curves": {
+      "32": 0.839392438341168,
+      "54": 0.8532932121895296,
+      "93": 0.8648301746424334,
+      "158": 0.87897726765122,
+      "271": 0.8822707327204625,
+      "464": 0.9066292229158499,
+      "794": 0.9124550742250127,
+      "1359": 0.9042832672872496,
+      "2326": 0.8927915426153833,
+      "3981": 0.8928568764738968,
+      "6813": 0.8767214401504766,
+      "11659": 0.8685078785862419,
+      "19953": 0.8597186612447733,
+      "34145": 0.8562044440055402,
+      "58434": 0.8596541437822435
+    },
+    "best_nfeat": 794,
+    "best_cv": 0.9124550742250127,
+    "best_gamma": 3.98271745613146e-05
+  }
+}

scripts/classification.py

@@ -1,50 +0,0 @@
-#!/usr/bin/env python3
-"""
-classification.py — single-seed image-level vs patient-level comparison.
-
-Usage:
-    conda activate fundus_imaging
-    python scripts/classification.py
-"""
-
-import os, sys, json
-sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-
-from classes import PatientLeakageClassifier
-
-ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
-MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
-SEED = 20
-
-clf = PatientLeakageClassifier(
-    os.path.join(ROOT, "results", "simple_patient_manifest.csv"),
-    os.path.join(ROOT, "features"),
-    n_jobs=8)
-
-print(f"{'='*60}")
-print(f"IMAGE-LEVEL vs PATIENT-LEVEL (seed={SEED})")
-print(f"{'='*60}")
-
-results = []
-for name in MODELS:
-    print(f"\n  {name} ...")
-    img = clf.run(name, SEED, "image")
-    pat = clf.run(name, SEED, "patient")
-    results.append({"model": name,
-                    "image_cv": img["cv"], "image_test": img["test"],
-                    "patient_cv": pat["cv"], "patient_test": pat["test"],
-                    "drop": img["test"] - pat["test"]})
-    print(f"    Image:   CV={img['cv']:.4f}  Test={img['test']:.4f}")
-    print(f"    Patient: CV={pat['cv']:.4f}  Test={pat['test']:.4f}")
-    print(f"    Drop:    {img['test'] - pat['test']:.4f}")
-
-print(f"\n  {'Model':<18s} {'Img-CV':>8s} {'Img-Test':>9s} "
-      f"{'Pat-CV':>8s} {'Pat-Test':>9s} {'Drop':>7s}")
-print(f"  {'-'*54}")
-for r in results:
-    print(f"  {r['model']:<18s} {r['image_cv']:>8.4f} {r['image_test']:>9.4f} "
-          f"{r['patient_cv']:>8.4f} {r['patient_test']:>9.4f} {r['drop']:>7.4f}")
-
-with open(os.path.join(ROOT, "results", "classification_results.json"), "w") as f:
-    json.dump(results, f, indent=2)
-print(f"\nDONE")

scripts/classification/classification_pca50.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+"""
+classification.py — single-seed image-level vs patient-level comparison.
+
+This script OWNS the shared classification-runs cache
+(scripts/cache/classification_runs.json): a flat list of per-(model, seed,
+split) run dicts. It reads any runs it needs from the cache and writes back any
+it has to compute, so the cache is maintained here. figure4_6.py references the
+same file to draw the 20-seed boxplots.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/classification.py
+"""
+
+import os, sys, json
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+
+from classes import PatientLeakageClassifier
+
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+ROOT = os.path.dirname(SCRIPT_DIR)
+MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
+SEED = 20
+CACHE_PATH = os.path.join(SCRIPT_DIR, "cache", "classification_runs.json")
+
+# Load the run cache: keep the raw list (to append to) and an index (to look up).
+runs = []
+if os.path.exists(CACHE_PATH):
+    with open(CACHE_PATH) as f:
+        runs = json.load(f)
+index = {(r["model"], r["seed"], r["split_type"]): r for r in runs}
+
+clf = None  # lazily created only if the cache is missing something
+
+
+def get_run(name, split):
+    """Return the cached run, else compute it and write it back to the cache."""
+    global clf
+    hit = index.get((name, SEED, split))
+    if hit is not None:
+        return hit
+    if clf is None:
+        clf = PatientLeakageClassifier(
+            os.path.join(ROOT, "results", "simple_patient_manifest.csv"),
+            os.path.join(ROOT, "features"), n_jobs=8)
+    print(f"    (cache miss for {name}/{split} — computing)")
+    r = clf.run(name, SEED, split)
+    runs.append(r)
+    index[(name, SEED, split)] = r
+    os.makedirs(os.path.dirname(CACHE_PATH), exist_ok=True)
+    with open(CACHE_PATH, "w") as f:
+        json.dump(runs, f, indent=2)
+    return r
+
+
+print(f"{'='*60}")
+print(f"IMAGE-LEVEL vs PATIENT-LEVEL (seed={SEED})")
+print(f"{'='*60}")
+
+results = []
+for name in MODELS:
+    print(f"\n  {name} ...")
+    img = get_run(name, "image")
+    pat = get_run(name, "patient")
+    results.append({"model": name,
+                    "image_cv": img["cv"], "image_test": img["test"],
+                    "patient_cv": pat["cv"], "patient_test": pat["test"],
+                    "drop": img["test"] - pat["test"]})
+    print(f"    Image:   CV={img['cv']:.4f}  Test={img['test']:.4f}")
+    print(f"    Patient: CV={pat['cv']:.4f}  Test={pat['test']:.4f}")
+    print(f"    Drop:    {img['test'] - pat['test']:.4f}")
+
+print(f"\n  {'Model':<18s} {'Img-CV':>8s} {'Img-Test':>9s} "
+      f"{'Pat-CV':>8s} {'Pat-Test':>9s} {'Drop':>7s}")
+print(f"  {'-'*54}")
+for r in results:
+    print(f"  {r['model']:<18s} {r['image_cv']:>8.4f} {r['image_test']:>9.4f} "
+          f"{r['patient_cv']:>8.4f} {r['patient_test']:>9.4f} {r['drop']:>7.4f}")
+
+with open(os.path.join(ROOT, "results", "classification_pca50.json"), "w") as f:
+    json.dump(results, f, indent=2)
+print(f"\nDONE")

scripts/classification/classification_siamese.py

@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+"""classification_siamese.py — Run classification with siamese patient manifest."""
+
+import os, sys, json
+import numpy as np
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+from classes import PatientLeakageClassifier
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+MANIFEST = os.path.join(ROOT, "results", "siamese_manifest.csv")
+FEATURES_DIR = os.path.join(ROOT, "features")
+RESULTS_DIR = os.path.join(ROOT, "results")
+SEED = 20
+MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
+
+clf = PatientLeakageClassifier(MANIFEST, FEATURES_DIR, n_jobs=6)
+
+print("=" * 60)
+print("Siamese-based patient classification")
+print("=" * 60)
+
+results = []
+for model_name in MODELS:
+    print(f"\n  {model_name} ...", flush=True)
+    img = clf.run(model_name, SEED, "image")
+    pat = clf.run(model_name, SEED, "patient")
+    drop = img["test"] - pat["test"]
+    results.append({
+        "model": model_name, "manifest": "siamese",
+        "image_cv": img["cv"], "image_test": img["test"],
+        "patient_cv": pat["cv"], "patient_test": pat["test"],
+        "drop": drop,
+    })
+    print(f"    Image:   CV={img['cv']:.4f}  Test={img['test']:.4f}")
+    print(f"    Patient: CV={pat['cv']:.4f}  Test={pat['test']:.4f}  "
+          f"Drop={drop:.4f}")
+
+# Save
+out = os.path.join(RESULTS_DIR, "classification_siamese.json")
+with open(out, "w") as f:
+    json.dump(results, f, indent=2)
+
+# Comparison table
+print(f"\n{'='*80}")
+print("COMPARISON — All three patient-clustering methods (seed=20)")
+print(f"{'='*80}")
+
+def safe_load(path):
+    if os.path.exists(path):
+        with open(path) as f: return {r["model"]: r for r in json.load(f)}
+    return None
+
+pca50 = safe_load(os.path.join(RESULTS_DIR, "classification_pca50.json"))
+thumbnail = safe_load(os.path.join(RESULTS_DIR, "classification_thumbnail.json"))
+siamese = {r["model"]: r for r in results}
+
+print(f"\n{'Model':<18s} {'PCA50 Pat':>10s} {'Thumb Pat':>11s} {'Siam Pat':>10s} "
+      f"{'PCA50 Drop':>11s} {'Thumb Drop':>11s} {'Siam Drop':>10s}")
+print("-" * 82)
+for m in MODELS:
+    f_pat = f"{pca50[m]['patient_test']:>10.4f}" if pca50 else "       N/A"
+    t_pat = f"{thumbnail[m]['patient_test']:>11.4f}" if thumbnail else "        N/A"
+    f_drop = f"{pca50[m]['drop']:>11.4f}" if pca50 else "        N/A"
+    t_drop = f"{thumbnail[m]['drop']:>11.4f}" if thumbnail else "        N/A"
+    print(f"{m:<18s} {f_pat} {t_pat} "
+          f"{siamese[m]['patient_test']:>10.4f} {f_drop} {t_drop} "
+          f"{siamese[m]['drop']:>10.4f}")
+
+print(f"\nSaved → {out}")
+print("DONE")

scripts/classification/classification_thumbnail.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+"""
+classification_thumbnail.py — Run the classification pipeline using the
+manuscript's thumbnail-based K-means patient manifest.
+
+Saves results alongside the feature-based results for comparison.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/classification_thumbnail.py
+"""
+
+import os, sys, csv, json, re
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+
+from classes import PatientIdentifier, PatientLeakageClassifier
+
+# ---------------------------------------------------------------------------
+# Config
+# ---------------------------------------------------------------------------
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+DATASET_PATH = os.path.join(os.path.dirname(ROOT),
+                            "The IQ-OTHNCCD lung cancer dataset")
+FEATURES_DIR = os.path.join(ROOT, "features")
+RESULTS_DIR = os.path.join(ROOT, "results")
+
+SEED = 20
+MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
+
+# Known patient counts per class
+PATIENT_COUNTS = {
+    "Bengin cases": 15,
+    "Malignant cases": 40,
+    "Normal cases": 55,
+}
+
+# ---------------------------------------------------------------------------
+# Step 1: Build thumbnail-based patient manifest
+# ---------------------------------------------------------------------------
+
+print("=" * 60)
+print("Building thumbnail-based patient manifest")
+print("=" * 60)
+
+# Collect all image filenames and class labels (same order as feature extraction)
+image_paths, all_labels, all_fnames = [], [], []
+for class_name in sorted(os.listdir(DATASET_PATH)):
+    class_path = os.path.join(DATASET_PATH, class_name)
+    if not os.path.isdir(class_path):
+        continue
+    for file in sorted(os.listdir(class_path)):
+        if file.lower().endswith((".png", ".jpg", ".jpeg")):
+            image_paths.append(os.path.join(class_path, file))
+            all_labels.append(class_name)
+            all_fnames.append(file)
+
+print(f"Found {len(image_paths)} images across "
+      f"{len(set(all_labels))} classes")
+
+# Run thumbnail K-means
+labels_arr = np.array(all_labels)
+fnames_arr = np.array(all_fnames)
+
+identifier = PatientIdentifier(
+    patient_estimates=PATIENT_COUNTS, random_state=42)
+groups = identifier.identify_from_thumbnails(
+    DATASET_PATH, fnames_arr, labels_arr)
+
+# Build manifest
+manifest = identifier.build_assignment_dict(groups, fnames_arr, labels_arr)
+
+# Save manifest
+manifest_path = os.path.join(RESULTS_DIR, "thumbnail_patient_manifest.csv")
+
+
+def f2n(fname):
+    m = re.search(r"\((\d+)\)", fname)
+    num = int(m.group(1)) if m else None
+    for cls_key, prefix in [
+        ("Bengin cases", "B"), ("Malignant cases", "M"), ("Normal cases", "N"),
+    ]:
+        if fname.startswith(cls_key.rstrip("s")):
+            return f"{prefix}_{num:03d}" if num else fname
+    return fname
+
+
+with open(manifest_path, "w", newline="") as f:
+    writer = csv.writer(f)
+    writer.writerow(["patient_id", "class", "n_images", "images"])
+    for pid in sorted(manifest.keys()):
+        imgs = manifest[pid]
+        short_names = [f2n(img) for img in imgs]
+        cls = pid.split("_")[0]  # "benign_0" → "benign"
+        cls = {"benign": "Benign", "malig": "Malignant",
+               "normal": "Normal"}.get(cls, cls)
+        writer.writerow([pid, cls, len(imgs), ";".join(short_names)])
+
+print(f"Manifest saved → {manifest_path}")
+print(f"  {len(manifest)} estimated patients")
+
+# Per-class stats
+for cls in ["Bengin cases", "Malignant cases", "Normal cases"]:
+    short_cls = {"Bengin cases": "benign", "Malignant cases": "malig",
+                 "Normal cases": "normal"}[cls]
+    cls_patients = {k: v for k, v in manifest.items()
+                    if k.startswith(short_cls)}
+    n_pat = len(cls_patients)
+    n_img = sum(len(v) for v in cls_patients.values())
+    print(f"  {cls}: {n_pat} patients, {n_img} images")
+
+# ---------------------------------------------------------------------------
+# Step 2: Run classification with thumbnail manifest
+# ---------------------------------------------------------------------------
+
+print(f"\n{'=' * 60}")
+print("Running classification (thumbnail manifest, seed=20)")
+print("=" * 60)
+
+clf = PatientLeakageClassifier(manifest_path, FEATURES_DIR, n_jobs=6)
+
+results = []
+for model_name in MODELS:
+    print(f"\n  {model_name} ...", flush=True)
+    img = clf.run(model_name, SEED, "image")
+    pat = clf.run(model_name, SEED, "patient")
+    results.append({
+        "model": model_name,
+        "manifest": "thumbnail",
+        "image_cv": img["cv"],
+        "image_test": img["test"],
+        "patient_cv": pat["cv"],
+        "patient_test": pat["test"],
+        "drop": img["test"] - pat["test"],
+    })
+    print(f"    Image:  CV={img['cv']:.4f}  Test={img['test']:.4f}")
+    print(f"    Patient: CV={pat['cv']:.4f}  Test={pat['test']:.4f}  "
+          f"Drop={img['test'] - pat['test']:.4f}")
+
+# ---------------------------------------------------------------------------
+# Save results
+# ---------------------------------------------------------------------------
+
+out_path = os.path.join(RESULTS_DIR, "classification_thumbnail.json")
+with open(out_path, "w") as f:
+    json.dump(results, f, indent=2)
+
+# Also load feature-based results for side-by-side comparison
+feat_path = os.path.join(RESULTS_DIR, "classification_results.json")
+if os.path.exists(feat_path):
+    with open(feat_path) as f:
+        feat_results = json.load(f)
+
+    print(f"\n{'=' * 80}")
+    print("COMPARISON: Feature-based (PCA-50) vs Thumbnail K-means")
+    print(f"{'=' * 80}")
+    print(f"{'Model':<18s} {'Feat Image':>11s} {'Thumb Image':>12s} "
+          f"{'Feat Pat':>9s} {'Thumb Pat':>10s} {'Feat Drop':>10s} "
+          f"{'Thumb Drop':>11s}")
+    print("-" * 78)
+    for tr, fr in zip(results, feat_results):
+        assert tr["model"] == fr["model"]
+        print(f"{tr['model']:<18s} {fr['image_test']:>11.4f} "
+              f"{tr['image_test']:>12.4f} {fr['patient_test']:>9.4f} "
+              f"{tr['patient_test']:>10.4f} {fr['drop']:>10.4f} "
+              f"{tr['drop']:>11.4f}")
+
+print(f"\nSaved → {out_path}")
+print("DONE")

scripts/clustering_validation/validate_pca50.py

@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+"""Validate PCA-50 feature K-means on the same Task06 test patients as siamese."""
+
+import os, sys, json, argparse
+import numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from PIL import Image as PILImage
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+from classes import PatientIdentifier, FeatureExtractor
+from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+PNG_DIR = os.path.join(ROOT, "features", "task06_pngs", "test")
+RESULTS_DIR = os.path.join(ROOT, "results", "clustering_validation")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+SEED = 42
+MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
+
+def cluster_purity_stats(y_true, y_pred):
+    purities = []; overall_correct = 0
+    for c in np.unique(y_pred):
+        mask = y_pred == c
+        _, counts = np.unique(y_true[mask], return_counts=True)
+        purities.append(counts.max() / mask.sum())
+        overall_correct += counts.max()
+    purities = np.array(purities)
+    return {"overall": float(overall_correct / len(y_true)),
+            "median": float(np.median(purities)), "mean": float(np.mean(purities)),
+            "frac_gt_70": float((purities > 0.7).mean()),
+            "frac_gt_90": float((purities > 0.9).mean())}
+
+def patient_capture_stats(y_true, y_pred):
+    captures = []
+    for p in np.unique(y_true):
+        mask = y_true == p
+        p_clusters = y_pred[mask]
+        _, counts = np.unique(p_clusters, return_counts=True)
+        captures.append(counts.max() / mask.sum())
+    captures = np.array(captures)
+    return {"median": float(np.median(captures)), "mean": float(np.mean(captures)),
+            "frac_gt_50": float((captures > 0.5).mean())}
+
+def plot_assignment_matrix(y_true, y_pred, out_path, title):
+    true_patients = sorted(np.unique(y_true))
+    pred_clusters = sorted(np.unique(y_pred))
+    matrix = np.zeros((len(true_patients), len(pred_clusters)))
+    for i, p in enumerate(true_patients):
+        for j, c in enumerate(pred_clusters):
+            matrix[i, j] = ((y_true == p) & (y_pred == c)).sum()
+    matrix_norm = matrix / (matrix.sum(axis=1, keepdims=True) + 1e-8)
+    fig, ax = plt.subplots(figsize=(max(14, len(pred_clusters)*0.22), max(10, len(true_patients)*0.18)))
+    cmap = plt.cm.YlOrRd.copy()
+    cmap.set_under('white')
+    ax.imshow(matrix_norm, aspect="auto", cmap=cmap, vmin=1e-6, vmax=1)
+    ax.set_xticks(range(len(pred_clusters)))
+    ax.set_xticklabels([f"c{c}" for c in pred_clusters], fontsize=6, rotation=90)
+    ax.set_yticks(range(len(true_patients)))
+    ax.set_yticklabels(true_patients, fontsize=7)
+    ax.set_xlabel("Predicted cluster"); ax.set_ylabel("True patient")
+    ax.set_title(title, fontsize=11)
+    plt.colorbar(ax.images[0], ax=ax, label="Fraction of patient's slices")
+    plt.tight_layout()
+    os.makedirs(os.path.dirname(out_path), exist_ok=True)
+    plt.savefig(out_path, dpi=150); plt.close()
+    print(f"  Heatmap → {out_path}")
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--model", default="VGG16", choices=MODELS)
+    ap.add_argument("--all-models", action="store_true")
+    ap.add_argument("--tag", default="")
+    ap.add_argument("--full", action="store_true",
+                    help="Run on full NIfTI dataset (63 patients, not just test).")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+    os.makedirs(RESULTS_DIR, exist_ok=True); os.makedirs(PLOTS_DIR, exist_ok=True)
+
+    if args.full:
+        tag = (tag or "") + "_full"
+        from classes import NiftiSliceDataset
+        VOL = os.path.join(os.path.dirname(ROOT), "Task06_Lung", "imagesTr")
+        ds = NiftiSliceDataset(VOL, random_slices=False, seed=SEED, rotate_deg=90)
+        ds.load_all_slices(stride=1)
+        test_slices = ds.slices_rgb
+        test_pids = ds.patient_labels
+        test_fnames = ds.filenames
+        y_true_map = {f: pid for f, pid in zip(test_fnames, test_pids)}
+        k = ds.n_patients
+        y_true = np.array([y_true_map[f] for f in test_fnames])
+        print(f"  FULL dataset: {len(test_fnames)} slices, {k} patients")
+    else:
+        manifest_path = os.path.join(PNG_DIR, "manifest.json")
+        if not os.path.exists(manifest_path):
+            print(f"Test PNGs not found at {PNG_DIR}. Run train_siamese.py first.")
+            sys.exit(1)
+        with open(manifest_path) as f: png = json.load(f)
+        test_paths = png["paths"]
+        test_pids = np.array(png["patient_ids"])
+        test_fnames = [p.split("/")[-1].replace(".png", "") for p in test_paths]
+        y_true_map = {f: pid for f, pid in zip(test_fnames, test_pids)}
+        k = len(np.unique(test_pids))
+        y_true = np.array([y_true_map[f] for f in test_fnames])
+        test_slices = [PILImage.open(p).convert("L") for p in test_paths]
+        print(f"  {len(test_fnames)} slices, {k} patients (same test set as siamese)")
+
+    models_to_run = MODELS if args.all_models else [args.model]
+
+    for model_name in models_to_run:
+        print(f"\n{'─'*50}")
+        print(f"Model: {model_name}")
+        print(f"{'─'*50}")
+
+        print(f"  Extracting features ...", flush=True)
+        ext = FeatureExtractor(model_name=model_name)
+        X, _, _ = ext.extract_from_images(test_slices)
+
+        print(f"  Clustering PCA-50 into k={k} ...", flush=True)
+        ident = PatientIdentifier(patient_estimates={"lung": k}, random_state=SEED)
+        groups = ident.identify_from_features(X, np.full(len(X), "lung"))
+        y_pred = np.array([str(g) for g in groups])
+
+        purity = cluster_purity_stats(y_true, y_pred)
+        capture = patient_capture_stats(y_true, y_pred)
+        ari = adjusted_rand_score(y_true, y_pred)
+        nmi = normalized_mutual_info_score(y_true, y_pred)
+
+        print(f"\n  RESULTS — PCA-50 K-means ({model_name})")
+        print(f"  {'─'*45}")
+        print(f"  ARI:              {ari:.4f}")
+        print(f"  NMI:              {nmi:.4f}")
+        print(f"  Overall purity:   {purity['overall']:.3f} ({purity['overall']*100:.1f}%)")
+        print(f"  Mean capture:     {capture['mean']:.3f} ({capture['mean']*100:.1f}%)")
+
+        results = {"method": f"pca50_{model_name}", "model": model_name,
+                   "n_patients": int(k), "n_slices": len(test_fnames),
+                   "ARI": ari, "NMI": nmi,
+                   "cluster_purity": purity, "patient_capture": capture}
+        out_json = os.path.join(RESULTS_DIR, f"validate_pca50_{model_name}{tag}.json")
+        with open(out_json, "w") as f: json.dump(results, f, indent=2)
+        print(f"  Metrics → {out_json}")
+
+        plot_assignment_matrix(y_true, y_pred,
+            os.path.join(PLOTS_DIR, "clustering_validation", "pca50",
+                         f"assignment_matrix_{model_name}{tag}.png"),
+            f"PCA-50 K-means {model_name} (ARI={ari:.3f}, purity={purity['overall']:.1%}, capture mean={capture['mean']:.1%})")
+
+    # Summary table if all models
+    if args.all_models:
+        print(f"\n{'='*60}")
+        print("SUMMARY — All models")
+        print(f"{'='*60}")
+        print(f"{'Model':<18s} {'ARI':>7s} {'Purity':>8s} {'Capture':>8s}")
+        print("-" * 43)
+        for model_name in MODELS:
+            j = os.path.join(RESULTS_DIR, f"validate_pca50_{model_name}{tag}.json")
+            if os.path.exists(j):
+                with open(j) as f: d = json.load(f)
+                print(f"{model_name:<18s} {d['ARI']:>7.3f} {d['cluster_purity']['overall']:>7.1%} {d['patient_capture']['mean']:>7.1%}")
+
+    print("DONE")
+
+if __name__ == "__main__": main()

scripts/clustering_validation/validate_siamese.py

@@ -0,0 +1,251 @@
+#!/usr/bin/env python3
+"""
+validate_siamese.py — Validate the trained siamese network against Task06
+held-out test patients.
+
+Runs both connected-components and edge-ranking clustering (known K) on the
+test set and reports metrics + assignment heatmaps.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/clustering_validation/validate_siamese.py
+    python scripts/clustering_validation/validate_siamese.py --threshold 0.95 --tag v2
+"""
+
+import os, sys, json, argparse
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+from classes import SiamesePatientMatcher
+from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+MODELS_DIR = os.path.join(ROOT, "models")
+RESULTS_DIR = os.path.join(ROOT, "results", "clustering_validation")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+PNG_DIR = os.path.join(ROOT, "features", "task06_pngs", "test")
+VOLUME_DIR = os.path.join(os.path.dirname(ROOT), "Task06_Lung", "imagesTr")
+SEED = 42
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+
+def cluster_purity_stats(y_true, y_pred):
+    purities = []
+    overall_correct = 0
+    for c in np.unique(y_pred):
+        mask = y_pred == c
+        _, counts = np.unique(y_true[mask], return_counts=True)
+        purities.append(counts.max() / mask.sum())
+        overall_correct += counts.max()
+    purities = np.array(purities)
+    return {
+        "overall": float(overall_correct / len(y_true)),
+        "median": float(np.median(purities)),
+        "mean": float(np.mean(purities)),
+        "frac_gt_70": float((purities > 0.7).mean()),
+        "frac_gt_90": float((purities > 0.9).mean()),
+    }
+
+
+def patient_capture_stats(y_true, y_pred):
+    captures = []
+    for p in np.unique(y_true):
+        mask = y_true == p
+        p_clusters = y_pred[mask]
+        _, counts = np.unique(p_clusters, return_counts=True)
+        captures.append(counts.max() / mask.sum())
+    captures = np.array(captures)
+    return {
+        "median": float(np.median(captures)),
+        "mean": float(np.mean(captures)),
+        "frac_gt_50": float((captures > 0.5).mean()),
+    }
+
+
+
+
+def score_manifest(manifest, y_true_map):
+    """Score a manifest against ground-truth patient IDs."""
+    preds, truths = [], []
+    for pred_pid, fnames in manifest.items():
+        for f in fnames:
+            preds.append(pred_pid)
+            truths.append(y_true_map.get(f, f"unknown_{f}"))
+    y_true = np.array(truths)
+    y_pred = np.array(preds)
+    purity = cluster_purity_stats(y_true, y_pred)
+    capture = patient_capture_stats(y_true, y_pred)
+    ari = adjusted_rand_score(y_true, y_pred)
+    nmi = normalized_mutual_info_score(y_true, y_pred)
+    return {"ARI": ari, "NMI": nmi,
+            "cluster_purity": purity, "patient_capture": capture,
+            "y_true": y_true, "y_pred": y_pred}
+
+
+def plot_assignment_matrix(y_true, y_pred, out_path, title):
+    true_patients = sorted(np.unique(y_true))
+    pred_clusters = sorted(np.unique(y_pred))
+    matrix = np.zeros((len(true_patients), len(pred_clusters)))
+    for i, p in enumerate(true_patients):
+        for j, c in enumerate(pred_clusters):
+            matrix[i, j] = ((y_true == p) & (y_pred == c)).sum()
+    matrix_norm = matrix / (matrix.sum(axis=1, keepdims=True) + 1e-8)
+
+    fig, ax = plt.subplots(figsize=(max(14, len(pred_clusters) * 0.22),
+                                    max(10, len(true_patients) * 0.18)))
+    cmap = plt.cm.YlOrRd.copy()
+    cmap.set_under('white')
+    im = ax.imshow(matrix_norm, aspect="auto", cmap=cmap, vmin=1e-6, vmax=1)
+    ax.set_xticks(range(len(pred_clusters)))
+    ax.set_xticklabels([f"c{c}" for c in pred_clusters], fontsize=6, rotation=90)
+    ax.set_yticks(range(len(true_patients)))
+    ax.set_yticklabels(true_patients, fontsize=7)
+    ax.set_xlabel("Predicted cluster"); ax.set_ylabel("True patient")
+    ax.set_title(title, fontsize=11)
+    plt.colorbar(im, ax=ax, label="Fraction of patient's slices")
+    plt.tight_layout()
+    os.makedirs(os.path.dirname(out_path), exist_ok=True)
+    plt.savefig(out_path, dpi=150)
+    plt.close()
+    print(f"  Heatmap → {out_path}")
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--model", default=os.path.join(MODELS_DIR, "siamese_resnet18.pt"))
+    ap.add_argument("--backbone", default="resnet18")
+    ap.add_argument("--threshold", type=float, default=0.9)
+    ap.add_argument("--tag", default="", help="Append tag to output filenames")
+    ap.add_argument("--full", action="store_true",
+                    help="Run on full NIfTI dataset (leaked training data).")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    os.makedirs(PLOTS_DIR, exist_ok=True)
+
+    # ---- Load test PNGs ----
+    if args.full:
+        tag = (tag or "") + "_full"
+        from classes import NiftiSliceDataset
+        VOL = os.path.join(os.path.dirname(ROOT), "Task06_Lung", "imagesTr")
+        ds = NiftiSliceDataset(VOL, random_slices=False, seed=SEED, rotate_deg=90)
+        ds.load_all_slices(stride=1)
+        # Export full PNGs to temp dir
+        import tempfile
+        tmpdir = tempfile.mkdtemp(prefix="task06_full_")
+        full_paths, full_pids, _ = ds.export_pngs(tmpdir)
+        test_paths = full_paths
+        test_pids = np.array(full_pids)
+        print(f"  FULL dataset: {len(test_paths)} slices, {len(np.unique(test_pids))} patients")
+        print(f"  (includes training data — leakage expected for siamese)")
+    else:
+        manifest_path = os.path.join(PNG_DIR, "manifest.json")
+        if not os.path.exists(manifest_path):
+            print(f"Test PNGs not found at {PNG_DIR}")
+            print("Run train_siamese.py first to generate the test set.")
+            sys.exit(1)
+        with open(manifest_path) as f:
+            png_manifest = json.load(f)
+        test_paths = png_manifest["paths"]
+        test_pids = np.array(png_manifest["patient_ids"])
+
+    # Short filenames for matching
+    test_fnames = [p.split("/")[-1].replace(".png", "") for p in test_paths]
+    y_true_map = {f: pid for f, pid in zip(test_fnames, test_pids)}
+    k = len(np.unique(test_pids))
+
+    print("=" * 60)
+    print(f"Siamese validation — held-out test set")
+    print("=" * 60)
+    print(f"  {len(test_paths)} slices, {k} patients")
+    print(f"  Model: {args.model}")
+
+    # ---- Load siamese model ----
+    print(f"\nLoading siamese model ...")
+    matcher = SiamesePatientMatcher(
+        args.model, backbone=args.backbone, input_size=224)
+
+    # ---- Connected components (no known K) ----
+    print(f"\n{'─'*50}")
+    print("Method 1: Connected components (threshold-based)")
+    print(f"{'─'*50}")
+    manifest_cc = matcher.identify_patients(
+        test_paths, filenames=test_fnames,
+        threshold=args.threshold, top_k=20, k=None)
+
+    # ---- Edge-ranking clustering (known K) ----
+    print(f"\n{'─'*50}")
+    print(f"Method 2: Edge-ranking clustering (k={k})")
+    print(f"{'─'*50}")
+    manifest_sc = matcher.identify_patients(
+        test_paths, filenames=test_fnames,
+        threshold=args.threshold, top_k=20, k=k)
+
+    # ---- Score both methods ----
+    results_cc = score_manifest(manifest_cc, y_true_map)
+    results_sc = score_manifest(manifest_sc, y_true_map)
+
+    print(f"\n{'='*60}")
+    print("RESULTS — Siamese patient identification")
+    print(f"{'='*60}")
+    print(f"\n{'Method':<30s} {'ARI':>7s} {'NMI':>7s} {'  Purity   Capture':>8s}")
+    print(f"{'':30s} {'':>7s} {'':>7s} {' (overall)  (mean)  ':>8s}")
+    print("-" * 60)
+
+    for name, r in [("Connected components", results_cc),
+                    ("Edge ranking (k=" + str(k) + ")", results_sc)]:
+        p = r["cluster_purity"]["overall"]
+        c = r["patient_capture"]["mean"]
+        print(f"{name:<30s} {r['ARI']:>7.3f} {r['NMI']:>7.3f} "
+              f"{p:>7.1%} {c:>7.1%}")
+
+    # ---- Save metrics ----
+    for suffix, r in [("cc", results_cc), ("sc", results_sc)]:
+        out = {k: v for k, v in r.items() if k not in ("y_true", "y_pred")}
+        out["method"] = suffix
+        out["threshold"] = args.threshold
+        out["n_patients"] = int(k)
+        out["n_slices"] = len(test_paths)
+        out_json = os.path.join(RESULTS_DIR, f"validate_siamese_{suffix}{tag}.json")
+        with open(out_json, "w") as f:
+            json.dump(out, f, indent=2)
+        print(f"  Metrics → {out_json}")
+
+        # ---- Plot siamese methods ----
+    for name, r in [("connected_components", results_cc),
+                    ("edge_rank", results_sc)]:
+        p = r["cluster_purity"]["overall"]
+        c = r["patient_capture"]["mean"]
+        title = (f"Siamese Patient-Cluster Assignment ({name})\n"
+                 f"(ARI={r['ARI']:.3f}, purity={p:.1%}, "
+                 f"capture mean={c:.1%})")
+        plot_assignment_matrix(
+            r["y_true"], r["y_pred"],
+            os.path.join(PLOTS_DIR, "clustering_validation", "siamese",
+                         f"assignment_matrix_{name}{tag}.png"),
+            title)
+
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/clustering_validation/validate_thumbnail.py

@@ -0,0 +1,156 @@
+#!/usr/bin/env python3
+"""Validate thumbnail K-means on the same Task06 test patients as siamese."""
+
+import os, sys, json, argparse
+import numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from PIL import Image as PILImage
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+from classes import PatientIdentifier
+from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
+
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+PNG_DIR = os.path.join(ROOT, "features", "task06_pngs", "test")
+RESULTS_DIR = os.path.join(ROOT, "results", "clustering_validation")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+SEED = 42
+
+def cluster_purity_stats(y_true, y_pred):
+    purities = []; overall_correct = 0
+    for c in np.unique(y_pred):
+        mask = y_pred == c
+        _, counts = np.unique(y_true[mask], return_counts=True)
+        purities.append(counts.max() / mask.sum())
+        overall_correct += counts.max()
+    purities = np.array(purities)
+    return {"overall": float(overall_correct / len(y_true)),
+            "median": float(np.median(purities)), "mean": float(np.mean(purities)),
+            "frac_gt_70": float((purities > 0.7).mean()),
+            "frac_gt_90": float((purities > 0.9).mean())}
+
+def patient_capture_stats(y_true, y_pred):
+    captures = []
+    for p in np.unique(y_true):
+        mask = y_true == p
+        p_clusters = y_pred[mask]
+        _, counts = np.unique(p_clusters, return_counts=True)
+        captures.append(counts.max() / mask.sum())
+    captures = np.array(captures)
+    return {"median": float(np.median(captures)), "mean": float(np.mean(captures)),
+            "frac_gt_50": float((captures > 0.5).mean())}
+
+def plot_assignment_matrix(y_true, y_pred, out_path, title):
+    true_patients = sorted(np.unique(y_true))
+    pred_clusters = sorted(np.unique(y_pred))
+    matrix = np.zeros((len(true_patients), len(pred_clusters)))
+    for i, p in enumerate(true_patients):
+        for j, c in enumerate(pred_clusters):
+            matrix[i, j] = ((y_true == p) & (y_pred == c)).sum()
+    matrix_norm = matrix / (matrix.sum(axis=1, keepdims=True) + 1e-8)
+    fig, ax = plt.subplots(figsize=(max(14, len(pred_clusters)*0.22), max(10, len(true_patients)*0.18)))
+    cmap = plt.cm.YlOrRd.copy()
+    cmap.set_under('white')
+    ax.imshow(matrix_norm, aspect="auto", cmap=cmap, vmin=1e-6, vmax=1)
+    ax.set_xticks(range(len(pred_clusters)))
+    ax.set_xticklabels([f"c{c}" for c in pred_clusters], fontsize=6, rotation=90)
+    ax.set_yticks(range(len(true_patients)))
+    ax.set_yticklabels(true_patients, fontsize=7)
+    ax.set_xlabel("Predicted cluster"); ax.set_ylabel("True patient")
+    ax.set_title(title, fontsize=11)
+    plt.colorbar(ax.images[0], ax=ax, label="Fraction of patient's slices")
+    plt.tight_layout()
+    os.makedirs(os.path.dirname(out_path), exist_ok=True)
+    plt.savefig(out_path, dpi=150); plt.close()
+    print(f"  Heatmap → {out_path}")
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--tag", default="")
+    ap.add_argument("--full", action="store_true",
+                    help="Run on full NIfTI dataset (63 patients, not just test).")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+    os.makedirs(RESULTS_DIR, exist_ok=True); os.makedirs(PLOTS_DIR, exist_ok=True)
+
+    if args.full:
+        tag = (tag or "") + "_full"
+        from classes import NiftiSliceDataset
+        VOL = os.path.join(os.path.dirname(ROOT), "Task06_Lung", "imagesTr")
+        ds = NiftiSliceDataset(VOL, random_slices=False, seed=SEED, rotate_deg=90)
+        ds.load_all_slices(stride=1)
+        test_paths = None
+        test_pids = ds.patient_labels
+        test_fnames = ds.filenames
+        y_true_map = {f: pid for f, pid in zip(test_fnames, test_pids)}
+        k = ds.n_patients
+        print("=" * 60)
+        print("Thumbnail K-means validation (manuscript method) — FULL")
+        print("=" * 60)
+        print(f"  FULL dataset: {len(test_fnames)} slices, {k} patients")
+        # Build thumbnails from loaded slices
+        thumbs = []
+        for sl in ds.slices_rgb:
+            img = PILImage.fromarray(sl).convert("L").resize((64, 64))
+            thumbs.append(np.asarray(img, dtype=np.float32).flatten() / 255.0)
+        thumbs = np.array(thumbs, dtype=np.float32)
+    else:
+        manifest_path = os.path.join(PNG_DIR, "manifest.json")
+        if not os.path.exists(manifest_path):
+            print(f"Test PNGs not found at {PNG_DIR}. Run train_siamese.py first.")
+            sys.exit(1)
+        with open(manifest_path) as f: png = json.load(f)
+        test_paths = png["paths"]
+        test_pids = np.array(png["patient_ids"])
+        test_fnames = [p.split("/")[-1].replace(".png", "") for p in test_paths]
+        y_true_map = {f: pid for f, pid in zip(test_fnames, test_pids)}
+        k = len(np.unique(test_pids))
+        print("=" * 60)
+        print("Thumbnail K-means validation (manuscript method)")
+        print("=" * 60)
+        print(f"  {len(test_fnames)} slices, {k} patients (same test set as siamese)")
+        # Build thumbnails from PNGs
+        print("  Building 64x64 thumbnails ...", flush=True)
+        thumbs = []
+        for p in test_paths:
+            img = PILImage.open(p).convert("L").resize((64, 64))
+            thumbs.append(np.asarray(img, dtype=np.float32).flatten() / 255.0)
+        thumbs = np.array(thumbs, dtype=np.float32)
+
+    print(f"  Clustering into k={k} ...", flush=True)
+    ident = PatientIdentifier(patient_estimates={"lung": k}, random_state=SEED)
+    groups = ident.identify_from_features(thumbs, np.full(len(thumbs), "lung"))
+    y_pred = np.array([str(g) for g in groups])
+    y_true = np.array([y_true_map[f] for f in test_fnames])
+
+    purity = cluster_purity_stats(y_true, y_pred)
+    capture = patient_capture_stats(y_true, y_pred)
+    ari = adjusted_rand_score(y_true, y_pred)
+    nmi = normalized_mutual_info_score(y_true, y_pred)
+
+    print(f"\n{'='*50}")
+    print("RESULTS — Thumbnail K-means (manuscript)")
+    print(f"{'='*50}")
+    print(f"  ARI:              {ari:.4f}")
+    print(f"  NMI:              {nmi:.4f}")
+    print(f"  Overall purity:   {purity['overall']:.3f} ({purity['overall']*100:.1f}%)")
+    print(f"  Mean capture:     {capture['mean']:.3f} ({capture['mean']*100:.1f}%)")
+
+    results = {"method": "thumbnail_64x64", "n_patients": int(k),
+               "n_slices": len(test_fnames), "ARI": ari, "NMI": nmi,
+               "cluster_purity": purity, "patient_capture": capture}
+    out_json = os.path.join(RESULTS_DIR, f"validate_thumbnail{tag}.json")
+    with open(out_json, "w") as f: json.dump(results, f, indent=2)
+    print(f"  Metrics → {out_json}")
+
+    plot_assignment_matrix(y_true, y_pred,
+        os.path.join(PLOTS_DIR, "clustering_validation", "thumbnail", f"assignment_matrix{tag}.png"),
+        f"Thumbnail K-means (ARI={ari:.3f}, purity={purity['overall']:.1%}, capture mean={capture['mean']:.1%})")
+    print("DONE")
+
+if __name__ == "__main__": main()

scripts/create_patient_groups.py

@@ -1,105 +0,0 @@
-#!/usr/bin/env python3
-"""
-create_patient_groups.py
-
-End-to-end pipeline for all five CNNs from the manuscript:
-  VGG16, DenseNet121, EfficientNetB1, MobileNetV2, ResNet50
-
-For each model:
-  1. Extract deep features and save to features/{Model}_features.npz
-  2. Generate PCA / t-SNE plot → plots/{Model}_pca_tsne.png
-  3. Estimate patient groups via K-means → features/{Model}_patient_groups.npy
-
-Usage:
-    conda activate fundus_imaging
-    python scripts/create_patient_groups.py
-"""
-
-import os
-import sys
-import numpy as np
-
-sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-
-from classes import FeatureExtractor, PatientIdentifier
-from classes.visualizations import plot_pca_tsne
-
-# ---------------------------------------------------------------------------
-# Configuration
-# ---------------------------------------------------------------------------
-
-BASE_PATH = os.path.join(
-    os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))),
-    "The IQ-OTHNCCD lung cancer dataset"
-)
-
-PROJECT_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
-FEATURES_DIR = os.path.join(PROJECT_ROOT, "features")
-PLOTS_DIR = os.path.join(PROJECT_ROOT, "plots")
-
-os.makedirs(FEATURES_DIR, exist_ok=True)
-os.makedirs(PLOTS_DIR, exist_ok=True)
-
-MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
-
-PATIENT_ESTIMATES = {
-    "Bengin cases": 15,
-    "Malignant cases": 40,
-    "Normal cases": 55,
-}
-
-LEGEND_NAMES = {
-    "Bengin cases": "Benign",
-    "Malignant cases": "Malignant",
-    "Normal cases": "Normal",
-}
-
-# ---------------------------------------------------------------------------
-# Run pipeline for each model
-# ---------------------------------------------------------------------------
-
-for model_name in MODELS:
-    print("\n" + "=" * 60)
-    print(f"MODEL: {model_name}")
-    print("=" * 60)
-
-    # --- Step 1: Extract features ---
-    print("\n  [1/3] Feature extraction ...")
-    extractor = FeatureExtractor(model_name=model_name)
-    features_path = os.path.join(FEATURES_DIR, f"{model_name}_features.npz")
-
-    if os.path.exists(features_path):
-        print(f"    Loading cached features from {features_path}")
-        X, Y, filenames = FeatureExtractor.load_features(features_path)
-    else:
-        X, Y, filenames = extractor.extract(BASE_PATH)
-        extractor.save_features(X, Y, filenames, FEATURES_DIR)
-
-    print(f"    {model_name}: X shape = {X.shape}")
-
-    # --- Step 2: PCA / t-SNE ---
-    print(f"\n  [2/3] PCA / t-SNE visualization ...")
-    plot_path = os.path.join(PLOTS_DIR, f"{model_name}_pca_tsne.png")
-    plot_pca_tsne(X, Y, legend_names=LEGEND_NAMES, output_path=plot_path)
-
-    # --- Step 3: Patient identification ---
-    print(f"\n  [3/3] Patient identification (K-means) ...")
-    identifier = PatientIdentifier(patient_estimates=PATIENT_ESTIMATES)
-    groups = identifier.identify(BASE_PATH, filenames, Y)
-
-    groups_path = os.path.join(FEATURES_DIR, f"{model_name}_patient_groups.npy")
-    np.save(groups_path, groups)
-    print(f"    Saved patient groups → {groups_path}")
-
-# ---------------------------------------------------------------------------
-# Summary
-# ---------------------------------------------------------------------------
-
-print("\n" + "=" * 60)
-print("ALL MODELS COMPLETE")
-print("=" * 60)
-for model_name in MODELS:
-    fp = os.path.join(FEATURES_DIR, f"{model_name}_features.npz")
-    gp = os.path.join(FEATURES_DIR, f"{model_name}_patient_groups.npy")
-    pp = os.path.join(PLOTS_DIR, f"{model_name}_pca_tsne.png")
-    print(f"  {model_name:18s}  features: {os.path.basename(fp):30s}  groups: {os.path.basename(gp)}")

scripts/figure6.py

@@ -1,97 +0,0 @@
-#!/usr/bin/env python3
-"""
-figure6.py — test accuracy boxplots across 20 seeds (Figure 6).
-
-Usage:
-    conda activate fundus_imaging
-    python scripts/figure6.py
-"""
-
-import os, sys, json
-import numpy as np
-import matplotlib
-matplotlib.use("Agg")
-import matplotlib.pyplot as plt
-from matplotlib.patches import Patch
-from tqdm import tqdm
-
-sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-
-from classes import PatientLeakageClassifier
-
-ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
-MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
-SEEDS = list(range(1, 21))
-
-clf = PatientLeakageClassifier(
-    os.path.join(ROOT, "results", "simple_patient_manifest.csv"),
-    os.path.join(ROOT, "features"),
-    n_jobs=8)
-
-# ---------------------------------------------------------------------------
-# Run all seeds
-# ---------------------------------------------------------------------------
-
-all_data = []  # flat list of per-run dicts for raw data file
-
-for name in tqdm(MODELS, desc="Models"):
-    for seed in tqdm(SEEDS, desc=f"  {name} seeds", leave=False):
-        for stype in ["image", "patient"]:
-            r = clf.run(name, seed, stype)
-            all_data.append(r)
-
-# Save raw data
-with open(os.path.join(ROOT, "results", "figure6_data.json"), "w") as f:
-    json.dump(all_data, f, indent=2)
-
-# ---------------------------------------------------------------------------
-# Build per-model accuracy lists
-# ---------------------------------------------------------------------------
-
-accs = {}  # model -> {image: [accs], patient: [accs]}
-for name in MODELS:
-    accs[name] = {"image": [], "patient": []}
-for r in all_data:
-    accs[r["model"]][r["split_type"]].append(r["test"])
-
-# ---------------------------------------------------------------------------
-# Plot
-# ---------------------------------------------------------------------------
-
-fig, ax = plt.subplots(1, 1, figsize=(10, 6))
-
-for i, name in enumerate(MODELS):
-    pos_img = i * 2 + 0.7
-    pos_pat = i * 2 + 1.3
-
-    for pos, stype, color in [(pos_img, "image", '#4C9BD4'),
-                               (pos_pat, "patient", '#6DBF6D')]:
-        data = accs[name][stype]
-        bp = ax.boxplot(data, positions=[pos], widths=0.5,
-                        patch_artist=True, showfliers=True,
-                        flierprops=dict(marker='o', markersize=3))
-        bp['boxes'][0].set_facecolor(color)
-        med = np.median(data)
-        ax.annotate(f"{med:.3f}", (pos, med), fontsize=6,
-                   ha='center', va='bottom')
-
-ax.legend(handles=[
-    Patch(facecolor='#4C9BD4', label='Image-level split'),
-    Patch(facecolor='#6DBF6D', label='Patient-level split'),
-], loc='lower right')
-
-ax.set_xticks([p + 1 for p in range(0, len(MODELS) * 2, 2)])
-ax.set_xticklabels(MODELS)
-ax.set_ylabel("Test accuracy")
-ax.set_title("Figure 6 — Image-level vs Patient-level test accuracy (20 seeds)")
-ax.set_ylim(0.7, 1.02)
-ax.grid(axis='y', alpha=0.3)
-
-plt.tight_layout()
-out = os.path.join(ROOT, "plots", "figure6.png")
-plt.savefig(out, dpi=150)
-plt.close()
-
-print(f"\nSaved → {out}")
-print(f"Saved → {ROOT}/results/figure6_data.json")
-print("DONE")

scripts/siamese_identify.py

@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+"""
+siamese_identify.py — Apply a trained siamese model to IQ-OTH/NCCD to build
+a patient manifest via connected-components clustering.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/siamese_identify.py
+    python scripts/siamese_identify.py --threshold 0.95
+    python scripts/siamese_identify.py --backbone resnet34
+"""
+
+import os
+import sys
+import csv
+import re
+import argparse
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from classes import SiamesePatientMatcher
+
+# ---------------------------------------------------------------------------
+# Config
+# ---------------------------------------------------------------------------
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+MODELS_DIR = os.path.join(ROOT, "models")
+RESULTS_DIR = os.path.join(ROOT, "results")
+DEFAULT_DATASET = os.path.join(os.path.dirname(ROOT), "The IQ-OTHNCCD lung cancer dataset")
+
+BACKBONE_INPUT_SIZES = {
+    "resnet18": 224,
+    "resnet34": 224,
+    "efficientnet_b0": 240,
+}
+
+
+def f2n(fname):
+    """Convert IQ-OTH filename to short form: 'Bengin case (1).jpg' → 'B_001'."""
+    m = re.search(r"\((\d+)\)", fname)
+    num = int(m.group(1)) if m else None
+    for cls_key, prefix in [
+        ("Bengin cases", "B"), ("Malignant cases", "M"), ("Normal cases", "N"),
+    ]:
+        if fname.startswith(cls_key.rstrip("s")):
+            return f"{prefix}_{num:03d}" if num else fname
+    return fname
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--model", default=os.path.join(MODELS_DIR, "siamese_resnet18.pt"),
+                    help="Path to trained siamese model.")
+    ap.add_argument("--backbone", default="resnet18",
+                    choices=list(BACKBONE_INPUT_SIZES.keys()))
+    ap.add_argument("--dataset", default=DEFAULT_DATASET,
+                    help="Path to IQ-OTH/NCCD dataset directory.")
+    ap.add_argument("--threshold", type=float, default=0.9,
+                    help="Minimum siamese probability to create an edge.")
+    ap.add_argument("--top-k", type=int, default=20,
+                    help="Top-K candidates to verify per slice.")
+    ap.add_argument("--cluster-method", default="edge_rank",
+                    choices=["edge_rank", "complete", "average"],
+                    help="Clustering on the siamese graph. 'edge_rank' is "
+                         "single-linkage (chains on OOD data); 'complete'/"
+                         "'average' are chaining-resistant agglomerative.")
+    ap.add_argument("--min-size", type=int, default=None,
+                    help="Absorb groups smaller than this into their nearest "
+                         "group (mitigates leakage-prone singletons).")
+    ap.add_argument("--max-size", type=int, default=None,
+                    help="Split groups larger than this at their natural gaps "
+                         "in siamese-distance space.")
+    ap.add_argument("--keep-k", action="store_true",
+                    help="Preserve the known patient count K while enforcing "
+                         "size bounds (balanced bisection + nearest-merge).")
+    ap.add_argument("--output", default=None,
+                    help="Output manifest path (default: results/siamese_manifest.csv).")
+    ap.add_argument("--device", default=None)
+    args = ap.parse_args()
+
+    input_size = BACKBONE_INPUT_SIZES[args.backbone]
+
+    # ---- Load model ----
+    print(f"Loading model: {args.model}")
+    print(f"  backbone={args.backbone}, input_size={input_size}")
+    matcher = SiamesePatientMatcher(
+        args.model, backbone=args.backbone,
+        device=args.device, input_size=input_size)
+
+    # ---- Collect IQ-OTH images ----
+    print(f"\nScanning dataset: {args.dataset}")
+    image_paths = []
+    class_labels = []
+    for class_name in sorted(os.listdir(args.dataset)):
+        class_path = os.path.join(args.dataset, class_name)
+        if not os.path.isdir(class_path):
+            continue
+        for fname in sorted(os.listdir(class_path)):
+            if fname.lower().endswith((".png", ".jpg", ".jpeg")):
+                image_paths.append(os.path.join(class_path, fname))
+                class_labels.append(class_name)
+
+    print(f"  Found {len(image_paths)} images across "
+          f"{len(set(class_labels))} classes")
+
+    # ---- Identify patients within each class (spectral clustering with known K) ----
+    KNOWN_K = {
+        "Bengin cases": 15,    # typo in original dataset
+        "Malignant cases": 40,
+        "Normal cases": 55,
+    }
+    all_assignments = {}
+
+    for class_name in sorted(set(class_labels)):
+        class_mask = [i for i, c in enumerate(class_labels) if c == class_name]
+        class_paths = [image_paths[i] for i in class_mask]
+        class_fnames = [os.path.basename(p) for p in class_paths]
+        k = KNOWN_K.get(class_name)
+
+        print(f"\n{'='*50}")
+        print(f"Class: {class_name} ({len(class_paths)} images, k={k})")
+        print(f"{'='*50}")
+
+        manifest = matcher.identify_patients(
+            class_paths,
+            filenames=[f2n(f) for f in class_fnames],
+            threshold=args.threshold,
+            top_k=args.top_k,
+            k=k,
+            cluster_method=args.cluster_method,
+            min_size=args.min_size,
+            max_size=args.max_size,
+            keep_k=args.keep_k,
+        )
+
+        # Prefix with class and a per-class running index. Enumerate rather than
+        # reuse the raw cluster id: rebalancing yields pids like "siamese_0_s0"
+        # whose last token ("s0") is not unique and would collide.
+        short_cls = {"Bengin cases": "Benign", "Malignant cases": "Malignant",
+                      "Normal cases": "Normal"}[class_name]
+        prefixed = {f"{short_cls}_{i}": imgs
+                    for i, (pid, imgs) in enumerate(manifest.items())}
+        all_assignments.update(prefixed)
+
+    # ---- Save manifest ----
+    output_path = args.output or os.path.join(
+        RESULTS_DIR, "siamese_manifest.csv")
+    os.makedirs(os.path.dirname(output_path), exist_ok=True)
+
+    # Determine class for each patient from the filenames
+    def get_class(fname):
+        if fname.startswith("B_"):
+            return "Benign"
+        elif fname.startswith("M_"):
+            return "Malignant"
+        elif fname.startswith("N_"):
+            return "Normal"
+        return "Unknown"
+
+    with open(output_path, "w", newline="") as f:
+        writer = csv.writer(f)
+        writer.writerow(["patient_id", "class", "n_images", "images"])
+        for pid in sorted(all_assignments.keys()):
+            imgs = all_assignments[pid]
+            cls = get_class(pid)
+            writer.writerow([pid, cls, len(imgs), ";".join(imgs)])
+
+    print(f"\nManifest saved → {output_path}")
+    print(f"  {len(all_assignments)} estimated patients, "
+          f"{sum(len(v) for v in all_assignments.values())} images")
+
+    # Summary per class
+    print(f"\n{'Class':<20s} {'Patients':>10s} {'Images':>8s} {'Mean imgs/pat':>14s}")
+    print("-" * 54)
+    for cls in ["Benign", "Malignant", "Normal"]:
+        cls_patients = {k: v for k, v in all_assignments.items()
+                        if get_class(k) == cls}
+        n_pat = len(cls_patients)
+        n_img = sum(len(v) for v in cls_patients.values())
+        mean = n_img / n_pat if n_pat > 0 else 0
+        print(f"{cls:<20s} {n_pat:>10d} {n_img:>8d} {mean:>14.1f}")
+
+    print("\nDONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/train_siamese.py

@@ -0,0 +1,499 @@
+#!/usr/bin/env python3
+"""
+train_siamese.py — Train a siamese CNN on Task06_Lung to determine if two
+CT slices come from the same patient.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/train_siamese.py
+    python scripts/train_siamese.py --backbone resnet34 --epochs 15
+"""
+
+import os
+import sys
+import json
+import time
+import argparse
+import numpy as np
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+import torch
+import torch.nn as nn
+import torch.optim as optim
+from torch.utils.data import Dataset, DataLoader
+from torchvision import transforms
+from PIL import Image
+from sklearn.metrics import accuracy_score, roc_auc_score
+from scipy.sparse.csgraph import connected_components
+from scipy.sparse import csr_matrix
+from collections import defaultdict
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from classes import NiftiSliceDataset, SiameseCNN
+
+# ---------------------------------------------------------------------------
+# Config
+# ---------------------------------------------------------------------------
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+MODELS_DIR = os.path.join(ROOT, "models")
+RESULTS_DIR = os.path.join(ROOT, "results")
+DEFAULT_CACHE_DIR = os.path.join(ROOT, "features", "task06_pngs")
+
+BACKBONES = {
+    "resnet18": (512, 224),
+    "resnet34": (512, 224),
+    "efficientnet_b0": (1280, 240),
+}
+
+SEED = 42
+BATCH_SIZE = 64
+EPOCHS = 15
+LR = 1e-4
+MIN_Z_GAP = 20
+HARD_NEG_FRAC = 0.5
+TEST_PATIENTS = 8
+N_PAIRS = 20000
+
+# ---------------------------------------------------------------------------
+# Pair dataset
+# ---------------------------------------------------------------------------
+
+
+class PairDataset(Dataset):
+    """Yield (img_A, img_B, label) pairs with hard negative mining."""
+
+    def __init__(self, paths, patient_ids, z_indices, n_pairs,
+                 transform=None, min_z_gap=MIN_Z_GAP,
+                 hard_neg_frac=HARD_NEG_FRAC, seed=SEED):
+        rng = np.random.default_rng(seed)
+        unique_patients = np.unique(patient_ids)
+
+        z_to_samples = defaultdict(list)
+        for idx in range(len(paths)):
+            z_to_samples[z_indices[idx]].append((idx, patient_ids[idx]))
+
+        n_pos = n_pairs // 2
+        n_neg = n_pairs - n_pos
+        n_hard = int(n_neg * hard_neg_frac)
+
+        pairs, labels = [], []
+        pos_z_pairs = []
+
+        # Positive pairs: same patient, z-distance >= min_z_gap
+        for _ in range(n_pos):
+            for __ in range(100):
+                pid = rng.choice(unique_patients)
+                idx = np.where(patient_ids == pid)[0]
+                if len(idx) < 2:
+                    continue
+                i, j = rng.choice(idx, size=2, replace=False)
+                if abs(int(z_indices[i]) - int(z_indices[j])) >= min_z_gap:
+                    pairs.append((i, j))
+                    labels.append(1)
+                    pos_z_pairs.append((z_indices[i], z_indices[j]))
+                    break
+        while len(labels) < n_pos:
+            pid = rng.choice(unique_patients)
+            idx = np.where(patient_ids == pid)[0]
+            if len(idx) < 2:
+                continue
+            i, j = rng.choice(idx, size=2, replace=False)
+            pairs.append((i, j))
+            labels.append(1)
+            pos_z_pairs.append((z_indices[i], z_indices[j]))
+
+        # Hard negatives: different patients, matched z-positions
+        for _ in range(n_hard):
+            z_i, z_j = pos_z_pairs[rng.integers(0, len(pos_z_pairs))]
+            for __ in range(100):
+                s_i = z_to_samples.get(z_i, [])
+                s_j = z_to_samples.get(z_j, [])
+                if len(s_i) < 1 or len(s_j) < 1:
+                    break
+                si = s_i[rng.integers(0, len(s_i))]
+                sj = s_j[rng.integers(0, len(s_j))]
+                if si[1] != sj[1]:
+                    pairs.append((si[0], sj[0]))
+                    labels.append(0)
+                    break
+
+        # Easy negatives
+        while len(labels) < n_pairs:
+            p1, p2 = rng.choice(unique_patients, size=2, replace=False)
+            i = rng.choice(np.where(patient_ids == p1)[0])
+            j = rng.choice(np.where(patient_ids == p2)[0])
+            pairs.append((i, j))
+            labels.append(0)
+
+        order = rng.permutation(len(labels))
+        self.pairs = [(pairs[o][0], pairs[o][1]) for o in order]
+        self.labels = [labels[o] for o in order]
+        self.paths = paths
+        self.transform = transform
+
+    def __len__(self):
+        return len(self.pairs)
+
+    def __getitem__(self, idx):
+        i, j = self.pairs[idx]
+        img_a = Image.open(self.paths[i]).convert("RGB")
+        img_b = Image.open(self.paths[j]).convert("RGB")
+        if self.transform:
+            img_a = self.transform(img_a)
+            img_b = self.transform(img_b)
+        return img_a, img_b, torch.tensor(self.labels[idx], dtype=torch.float32)
+
+
+# ---------------------------------------------------------------------------
+# Training
+# ---------------------------------------------------------------------------
+
+
+def train_epoch(model, loader, optimizer, criterion, device):
+    model.train()
+    total_loss, correct, n = 0.0, 0, 0
+    for a, b, y in loader:
+        a, b, y = a.to(device), b.to(device), y.to(device)
+        optimizer.zero_grad()
+        loss = criterion(model(a, b), y)
+        loss.backward()
+        optimizer.step()
+        preds = (torch.sigmoid(model(a, b)) > 0.5).float()
+        total_loss += loss.item() * len(y)
+        correct += (preds == y).sum().item()
+        n += len(y)
+    return total_loss / n, correct / n
+
+
+@torch.no_grad()
+def eval_epoch(model, loader, criterion, device):
+    model.eval()
+    total_loss, n = 0.0, 0
+    all_preds, all_labels = [], []
+    for a, b, y in loader:
+        a, b, y = a.to(device), b.to(device), y.to(device)
+        logits = model(a, b)
+        total_loss += criterion(logits, y).item() * len(y)
+        all_preds.extend(torch.sigmoid(logits).cpu().numpy())
+        all_labels.extend(y.cpu().numpy())
+        n += len(y)
+    all_labels = np.array(all_labels, dtype=int)
+    all_preds = np.array(all_preds, dtype=float)
+    acc = accuracy_score(all_labels, all_preds > 0.5)
+    try:
+        auc = roc_auc_score(all_labels, all_preds)
+    except ValueError:
+        auc = 0.5
+    return total_loss / n, acc, auc
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+
+def _export_subset(slices, patient_ids, filenames, cache_dir):
+    """Export a subset of slices as PNGs. Cached via manifest.json."""
+    manifest_path = os.path.join(cache_dir, "manifest.json")
+    # Include slice count in key to detect stale caches (e.g. different split)
+    cache_key = f"{len(slices)}_{len(patient_ids)}"
+    if os.path.exists(manifest_path):
+        with open(manifest_path) as f:
+            m = json.load(f)
+        if m.get("_key") == cache_key and len(m["paths"]) == len(slices):
+            return m["paths"], np.array(m["patient_ids"]), np.array(m["z_indices"])
+        else:
+            print(f"  Cache stale (key mismatch), re-exporting ...", flush=True)
+
+    print(f"  Exporting {len(slices)} PNGs to {cache_dir} ...", flush=True)
+    paths, pids, zs = [], [], []
+    for i, (sl, pid, fname) in enumerate(zip(slices, patient_ids, filenames)):
+        out_path = os.path.join(cache_dir, f"{fname}.png")
+        Image.fromarray(sl).save(out_path)
+        paths.append(out_path)
+        pids.append(pid)
+        zs.append(int(fname.split("_slice")[-1]))
+        if (i + 1) % 500 == 0:
+            print(f"    {i + 1}/{len(slices)} ...", flush=True)
+
+    with open(manifest_path, "w") as f:
+        json.dump({"_key": f"{len(slices)}_{len(patient_ids)}",
+                    "paths": paths, "patient_ids": [str(p) for p in pids],
+                    "z_indices": [int(z) for z in zs]}, f)
+    return paths, np.array(pids), np.array(zs)
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--backbone", default="resnet18",
+                    choices=list(BACKBONES.keys()))
+    ap.add_argument("--n-pairs", type=int, default=N_PAIRS)
+    ap.add_argument("--epochs", type=int, default=EPOCHS)
+    ap.add_argument("--lr", type=float, default=LR)
+    ap.add_argument("--batch-size", type=int, default=BATCH_SIZE)
+    ap.add_argument("--cache-dir", default=DEFAULT_CACHE_DIR)
+    ap.add_argument("--flip-vertical", action="store_true",
+                    help="Flip slices vertically to match IQ-OTH orientation.")
+    ap.add_argument("--rotate", type=int, default=0,
+                    choices=[0, 90, 180, 270],
+                    help="Rotate slices by N degrees (e.g. 90 if spine is on left).")
+    ap.add_argument("--test-patients", type=int, default=TEST_PATIENTS,
+                    help="Number of patients held out for final testing.")
+    ap.add_argument("--device", default=None)
+    args = ap.parse_args()
+
+    feat_dim, input_size = BACKBONES[args.backbone]
+    device = torch.device(args.device or (
+        "cuda" if torch.cuda.is_available() else "cpu"))
+    print(f"Device: {device}, backbone: {args.backbone}, "
+          f"input_size: {input_size}")
+
+    volume_dir = os.path.join(os.path.dirname(ROOT), "Task06_Lung", "imagesTr")
+
+    # ---- Split patients into train / held-out test ----
+    # First, get the list of patients without loading all data
+    from classes.nifti_dataset import NiftiSliceDataset as _DS
+    import glob as _glob
+    _nii = _glob.glob(os.path.join(volume_dir, "*.nii"))
+    _gz = [p for p in _glob.glob(os.path.join(volume_dir, "*.nii.gz"))
+           if os.path.basename(p)[:-7] not in
+           {os.path.basename(q)[:-4] for q in _nii}]
+    all_volumes = sorted(_nii + _gz)
+    all_patient_ids = [os.path.basename(p).replace(".nii.gz", "").replace(".nii", "")
+                       for p in all_volumes]
+    n_total = len(all_patient_ids)
+    rng = np.random.default_rng(SEED)
+    test_pids = set(rng.choice(all_patient_ids,
+                               size=min(args.test_patients, n_total - 1),
+                               replace=False))
+    train_pids_set = set(all_patient_ids) - test_pids
+    print(f"\nPatients: {n_total} total → {len(train_pids_set)} train, "
+          f"{len(test_pids)} held-out test\n")
+
+    # ---- Extract TRAIN slices (sparse, for training speed) ----
+    print("=" * 50)
+    print("Loading TRAIN slices (stride=3)")
+    print("=" * 50)
+    ds_train = _DS(volume_dir, random_slices=False, seed=SEED,
+                   flip_vertical=args.flip_vertical, rotate_deg=args.rotate)
+    ds_train.load_all_slices(stride=3)
+    # Filter to train patients only
+    train_slice_mask = np.isin(ds_train.patient_labels, list(train_pids_set))
+    train_slices = [sl for i, sl in enumerate(ds_train.slices_rgb)
+                    if train_slice_mask[i]]
+    train_pids_arr = ds_train.patient_labels[train_slice_mask]
+    train_fnames = ds_train.filenames[train_slice_mask]
+    train_zs = np.array([int(f.split("_slice")[-1]) for f in train_fnames])
+
+    # Export train PNGs
+    train_cache = os.path.join(args.cache_dir, "train")
+    os.makedirs(train_cache, exist_ok=True)
+    train_paths, _, _ = _export_subset(
+        train_slices, train_pids_arr, train_fnames, train_cache)
+
+    print(f"  Train: {len(train_paths)} slices, "
+          f"{len(np.unique(train_pids_arr))} patients\n")
+
+    # ---- Extract TEST slices (dense, for thorough eval) ----
+    print("=" * 50)
+    print("Loading TEST slices (stride=1, full central 60%)")
+    print("=" * 50)
+    ds_test = _DS(volume_dir, random_slices=False, seed=SEED,
+                  flip_vertical=args.flip_vertical)
+    ds_test.load_all_slices(stride=1)
+    test_slice_mask = np.isin(ds_test.patient_labels, list(test_pids))
+    test_slices = [sl for i, sl in enumerate(ds_test.slices_rgb)
+                   if test_slice_mask[i]]
+    test_pids_arr = ds_test.patient_labels[test_slice_mask]
+    test_fnames = ds_test.filenames[test_slice_mask]
+    test_zs = np.array([int(f.split("_slice")[-1]) for f in test_fnames])
+
+    # Export test PNGs
+    test_cache = os.path.join(args.cache_dir, "test")
+    os.makedirs(test_cache, exist_ok=True)
+    test_paths, _, _ = _export_subset(
+        test_slices, test_pids_arr, test_fnames, test_cache)
+
+    print(f"  Test:  {len(test_paths)} slices, "
+          f"{len(np.unique(test_pids_arr))} patients\n")
+
+    # ---- Transforms ----
+    train_tf = transforms.Compose([
+        transforms.Resize((input_size, input_size)),
+        transforms.RandomHorizontalFlip(),
+        transforms.ToTensor(),
+        transforms.Normalize(mean=[0.485, 0.456, 0.406],
+                           std=[0.229, 0.224, 0.225]),
+    ])
+    val_tf = transforms.Compose([
+        transforms.Resize((input_size, input_size)),
+        transforms.ToTensor(),
+        transforms.Normalize(mean=[0.485, 0.456, 0.406],
+                           std=[0.229, 0.224, 0.225]),
+    ])
+
+    # ---- Training pairs (80% of train patients for training, 20% for val monitoring) ----
+    train_unique = np.unique(train_pids_arr)
+    n_val_patients = max(2, int(len(train_unique) * 0.2))
+    val_pids_set = set(rng.choice(train_unique, size=n_val_patients, replace=False))
+    tr_pids_set = set(train_unique) - val_pids_set
+
+    tr_mask = np.isin(train_pids_arr, list(tr_pids_set))
+    val_mask = np.isin(train_pids_arr, list(val_pids_set))
+
+    tr_paths_list = [train_paths[i] for i in np.where(tr_mask)[0]]
+    val_paths_list = [train_paths[i] for i in np.where(val_mask)[0]]
+
+    print(f"  Training pairs from: {len(tr_paths_list)} slices, "
+          f"{len(tr_pids_set)} patients")
+    print(f"  Val pairs from:     {len(val_paths_list)} slices, "
+          f"{len(val_pids_set)} patients")
+
+    print(f"\nGenerating {args.n_pairs:,} training pairs ...", flush=True)
+    t0 = time.time()
+    tr_ds = PairDataset(tr_paths_list, train_pids_arr[tr_mask],
+                        train_zs[tr_mask],
+                        n_pairs=args.n_pairs, transform=train_tf, seed=SEED)
+    print(f"  ... done ({time.time() - t0:.0f}s)")
+
+    val_n = args.n_pairs // 4
+    print(f"Generating {val_n:,} validation pairs ...", flush=True)
+    val_ds = PairDataset(val_paths_list, train_pids_arr[val_mask],
+                         train_zs[val_mask],
+                         n_pairs=val_n, transform=val_tf, seed=SEED + 1)
+
+    tr_loader = DataLoader(tr_ds, batch_size=args.batch_size,
+                           shuffle=True, num_workers=4)
+    val_loader = DataLoader(val_ds, batch_size=args.batch_size, num_workers=4)
+
+    # ---- Model ----
+    model = SiameseCNN(args.backbone).to(device)
+    print(f"\nModel: {sum(p.numel() for p in model.parameters()):,} parameters")
+
+    for p in model.backbone.parameters():
+        p.requires_grad = False
+
+    criterion = nn.BCEWithLogitsLoss()
+    optimizer = optim.Adam(model.parameters(), lr=args.lr)
+
+    # ---- Train ----
+    best_val_acc = 0
+    unfreeze_epoch = max(1, args.epochs // 2)
+
+    print(f"\n{'Epoch':>6s} {'tr_loss':>8s} {'tr_acc':>8s} "
+          f"{'val_loss':>8s} {'val_acc':>8s} {'val_auc':>8s}")
+    print("-" * 52)
+
+    for epoch in range(1, args.epochs + 1):
+        if epoch == unfreeze_epoch + 1:
+            for p in model.backbone.parameters():
+                p.requires_grad = True
+            for g in optimizer.param_groups:
+                g["lr"] = args.lr * 0.1
+            print("  (unfreezing backbone, LR 0.1x)", flush=True)
+
+        tr_loss, tr_acc = train_epoch(model, tr_loader, optimizer, criterion, device)
+        val_loss, val_acc, val_auc = eval_epoch(model, val_loader, criterion, device)
+
+        star = "*" if val_acc > best_val_acc else ""
+        if val_acc > best_val_acc:
+            best_val_acc = val_acc
+            os.makedirs(MODELS_DIR, exist_ok=True)
+            torch.save(model.state_dict(),
+                       os.path.join(MODELS_DIR, f"siamese_{args.backbone}.pt"))
+
+        print(f"{epoch:>6d} {tr_loss:>8.4f} {tr_acc:>8.4f} "
+              f"{val_loss:>8.4f} {val_acc:>8.4f} {val_auc:>8.4f} {star}",
+              flush=True)
+
+    print(f"\nBest val accuracy: {best_val_acc:.4f}")
+
+    # ---- Final evaluation on held-out TEST set ----
+    print(f"\n{'='*50}")
+    print(f"Held-out TEST evaluation ({len(test_pids)} unseen patients, "
+          f"{len(test_paths)} slices)")
+    print("=" * 50)
+
+    model_path = os.path.join(MODELS_DIR, f"siamese_{args.backbone}.pt")
+    model.load_state_dict(torch.load(model_path, map_location=device, weights_only=True))
+    model.eval()
+
+    from classes.siamese import SiamesePatientMatcher as _SPM
+    matcher = _SPM(model_path, backbone=args.backbone,
+                   device=str(device), input_size=input_size)
+
+    test_pids_list = list(test_pids)
+    test_fnames_clean = np.array([f.split("/")[-1].replace(".png", "")
+                                   for f in test_paths])
+
+    # Connected components (no known K)
+    manifest_cc = matcher.identify_patients(
+        test_paths, filenames=test_fnames_clean,
+        threshold=0.9, top_k=20, k=None)
+
+    # Spectral clustering (known K)
+    manifest_sc = matcher.identify_patients(
+        test_paths, filenames=test_fnames_clean,
+        threshold=0.9, top_k=20, k=len(test_pids))
+
+    def score_manifest(manifest, y_true_map):
+        """Score a manifest against ground-truth patient IDs."""
+        # y_true_map: filename → true_patient_id
+        all_labels = []
+        all_preds = []
+        for pred_pid, fnames in manifest.items():
+            for f in fnames:
+                all_preds.append(pred_pid)
+                all_labels.append(y_true_map.get(f, f"unknown_{f}"))
+        all_labels = np.array(all_labels)
+        all_preds = np.array(all_preds)
+        purities = []
+        for c in np.unique(all_preds):
+            mask = all_preds == c
+            _, cts = np.unique(all_labels[mask], return_counts=True)
+            purities.append(cts.max() / mask.sum())
+        captures = []
+        for p in np.unique(all_labels):
+            mask = all_labels == p
+            _, cts = np.unique(all_preds[mask], return_counts=True)
+            captures.append(cts.max() / mask.sum())
+        return np.median(purities), np.median(captures), len(np.unique(all_preds))
+
+    # Build ground-truth map: filename → patient_id
+    true_map = {f: p for f, p in zip(test_fnames_clean, test_pids_arr)}
+
+    p_cc, c_cc, n_cc = score_manifest(manifest_cc, true_map)
+    p_sc, c_sc, n_sc = score_manifest(manifest_sc, true_map)
+
+    print(f"\n  {'Method':<30s} {'Groups':>7s} {'Purity':>8s} {'Capture':>8s}")
+    print(f"  {'-'*53}")
+    print(f"  {'Connected components':<30s} {n_cc:>7d} "
+          f"{p_cc:>7.1%} {c_cc:>7.1%}")
+    print(f"  {'Spectral (known K=' + str(len(test_pids)) + ')':<30s} {n_sc:>7d} "
+          f"{p_sc:>7.1%} {c_sc:>7.1%}")
+
+    # Compare to thumbnail baseline
+    results_path = os.path.join(RESULTS_DIR, "task06_clustering_validation.json")
+    if os.path.exists(results_path):
+        with open(results_path) as f:
+            prev = json.load(f)
+        tp = prev["methods"]["thumbnail_64x64"]["cluster_purity"]["median"]
+        tc = prev["methods"]["thumbnail_64x64"]["dominant_capture"]["median"]
+        print(f"  {'Thumbnail 64x64 (baseline)':<30s} {'—':>7s} "
+              f"{tp:>7.1%} {tc:>7.1%}")
+
+    print(f"\nModel saved → {model_path}")
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/figure1.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python3
+"""figure1.py — Sample CT images from the IQ-OTH/NCCD dataset, one per patient.
+Usage: python scripts/visualizations/figure1.py [--tag TAG]"""
+
+import os, sys, csv, argparse
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from PIL import Image
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+DATASET = os.path.join(os.path.dirname(ROOT), "The IQ-OTHNCCD lung cancer dataset")
+MANIFEST = os.path.join(ROOT, "results", "simple_patient_manifest.csv")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+
+ap = argparse.ArgumentParser()
+ap.add_argument("--tag", default="", help="Append tag to filename")
+ap.add_argument("--manifest", default=MANIFEST, help="Patient manifest CSV")
+args = ap.parse_args()
+tag = f"_{args.tag}" if args.tag else ""
+
+# Load manifest to get per-patient images
+patients = {"Benign": [], "Malignant": [], "Normal": []}
+with open(args.manifest, newline="") as f:
+    reader = csv.DictReader(f)
+    img_col = "confirmed_images" if "confirmed_images" in reader.fieldnames else "images"
+    for row in reader:
+        cls = row.get("class", "")
+        if cls in patients:
+            imgs = row[img_col].split(";")
+            if imgs:
+                patients[cls].append((row["patient_id"], imgs[0]))  # first image per patient
+
+CLASS_DIR = {"Benign": "Bengin cases", "Malignant": "Malignant cases",
+             "Normal": "Normal cases"}
+N_EXAMPLES = 3
+
+fig, axes = plt.subplots(3, N_EXAMPLES, figsize=(8, 9))
+for row, (cls_label, cls_dir) in enumerate(CLASS_DIR.items()):
+    # Pick first N_EXAMPLES patients for this class
+    selected = patients[cls_label][:N_EXAMPLES]
+    for col, (pid, short_name) in enumerate(selected):
+        ax = axes[row, col]
+        # Convert short name back to original filename
+        prefix = short_name[0]
+        num = int(short_name.split("_")[1])
+        cls_map = {"B": ("Bengin cases", "Bengin"), "M": ("Malignant cases", "Malignant"),
+                   "N": ("Normal cases", "Normal")}
+        dir_name, file_prefix = cls_map[prefix]
+        fname = f"{file_prefix} case ({num}).jpg"
+        img_path = os.path.join(DATASET, dir_name, fname)
+        try:
+            img = Image.open(img_path).convert("L")
+            ax.imshow(img, cmap="gray")
+        except Exception as e:
+            ax.text(0.5, 0.5, f"error: {e}", ha="center", va="center", fontsize=7)
+        ax.set_xticks([]); ax.set_yticks([])
+        if col == 0:
+            ax.set_ylabel(cls_label, fontsize=10, rotation=0,
+                          labelpad=20, va="center")
+
+fig.suptitle("Figure 1 — Sample images from the IQ-OTH/NCCD dataset",
+             fontsize=12, y=1.02)
+plt.tight_layout(rect=[0, 0, 1, 0.97])
+out = os.path.join(PLOTS_DIR, "figure1", f"figure1{tag}.png")
+os.makedirs(os.path.dirname(out), exist_ok=True)
+plt.savefig(out, dpi=150)
+plt.close()
+print(f"Saved → {out}")

scripts/visualizations/figure3_s1_s4.py

@@ -0,0 +1,179 @@
+#!/usr/bin/env python3
+"""Combined CV-accuracy-vs-features curves for Figures 3, S1, and S4.
+
+These three figures all compute the same per-model curve — RF importance
+ranking + a per-nfeatures gamma grid at seed 20 — differing only in which
+models and which split they show:
+
+    Figure 3  : VGG16, image-level      (1 panel)   == panel (a) of S1
+    Figure S1 : all 5 models, image     (5 panels)
+    Figure S4 : all 5 models, patient   (5 panels)
+
+Computing them together runs each (model, split) curve once (5 image + 5
+patient = 10 curves) and caches them, instead of recomputing VGG16's image
+curve for both Figure 3 and S1. The cache is incremental: dropping a model's
+entries (e.g. after a feature change) recomputes only that model.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/visualizations/figure3_s1_s4.py                # all three
+    python scripts/visualizations/figure3_s1_s4.py --only s4      # just Fig S4
+    python scripts/visualizations/figure3_s1_s4.py --from-cache   # re-render only
+    python scripts/visualizations/figure3_s1_s4.py --force        # recompute all
+"""
+
+import os
+import sys
+import json
+import argparse
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from tqdm import tqdm
+
+os.environ.setdefault("OMP_NUM_THREADS", "1")
+os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
+os.environ.setdefault("MKL_NUM_THREADS", "1")
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+# Shared results-cache tier (alongside classification_runs.json). These curves
+# are model-selection output (best nfeat/gamma/cv), not plotting scaffolding, so
+# they live in scripts/cache/ rather than under visualizations/.
+CACHE_DIR = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "cache")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+
+from classes import PatientLeakageClassifier
+
+# Panel order matches the published supplementary figures.
+MODELS = ["VGG16", "MobileNetV2", "DenseNet121", "ResNet50", "EfficientNetB1"]
+SEED = 20
+DATA_PATH = os.path.join(CACHE_DIR, "cv_curves.json")
+MANIFEST = os.path.join(ROOT, "results", "simple_patient_manifest.csv")
+
+# 16-point log grid of feature counts (each model caps at its own max dim).
+NFEATS = np.unique(np.round(np.logspace(1.5, 5, 16)).astype(int)).tolist()
+GAMMA_LOGSPACE = (-1.5, 1, 7)
+PANELS = ["(a)", "(b)", "(c)", "(d)", "(e)"]
+
+
+def key(model, split):
+    return f"{model}::{split}"
+
+
+def load_cache():
+    if os.path.exists(DATA_PATH):
+        with open(DATA_PATH) as f:
+            return json.load(f)
+    return {}
+
+
+def compute(cache, splits, n_jobs, manifest):
+    """Fill any missing (model, split) curves in the cache; return updated cache."""
+    todo = [(m, s) for s in splits for m in MODELS if key(m, s) not in cache]
+    if not todo:
+        print(f"Cache complete for requested figures ({len(cache)} curves).")
+        return cache
+    clf = PatientLeakageClassifier(manifest, os.path.join(ROOT, "features"),
+                                   n_jobs=n_jobs)
+    for m, s in tqdm(todo, desc="Computing missing curves"):
+        cache[key(m, s)] = clf.cv_curve(m, SEED, s, nfeatures_list=NFEATS,
+                                        gamma_logspace=GAMMA_LOGSPACE)
+    os.makedirs(CACHE_DIR, exist_ok=True)
+    with open(DATA_PATH, "w") as f:
+        json.dump(cache, f, indent=2)
+    print(f"Computed {len(todo)} curves; cache now {len(cache)}.")
+    return cache
+
+
+def _draw_curve(ax, curve, color, ylim, title=None):
+    # Keys are ints on a fresh compute but strings once round-tripped through
+    # JSON; iterate items so either works.
+    items = sorted((int(n), v) for n, v in curve["curves"].items())
+    nfeats = [n for n, _ in items]
+    accs = [v * 100 for _, v in items]
+    ax.plot(nfeats, accs, "o-", color=color, markersize=4, linewidth=1.2)
+    ax.set_xscale("log")
+    ax.set_xlabel("Number of selected features")
+    ax.set_ylabel("CV accuracy (%)")
+    ax.set_ylim(*ylim)
+    ax.grid(alpha=0.3)
+    bn, ba = curve["best_nfeat"], curve["best_cv"] * 100
+    ax.axvline(bn, color="red", linestyle="--", linewidth=0.8, alpha=0.6)
+    ax.plot(bn, ba, "r*", markersize=11)
+    ax.annotate(f"n={bn}\n{ba:.2f}%", (bn, ba), fontsize=8,
+                xytext=(10, -10), textcoords="offset points", color="red")
+    if title:
+        ax.set_title(title, fontsize=10)
+
+
+def render_fig3(cache, tag):
+    fig, ax = plt.subplots(figsize=(8, 5))
+    _draw_curve(ax, cache[key("VGG16", "image")], "#D62728", (92, 100))
+    ax.set_title("Figure 3 — VGG16 Image-level CV accuracy vs number of "
+                 "selected features", fontsize=12)
+    _save(fig, "figure3", f"figure3{tag}.png")
+
+
+def render_grid(cache, split, color, ylim, suptitle, subdir, fname):
+    fig, axes = plt.subplots(2, 3, figsize=(14, 9))
+    axes = axes.flatten()
+    for i, model in enumerate(MODELS):
+        _draw_curve(axes[i], cache[key(model, split)], color, ylim,
+                    title=f"{PANELS[i]} {model}")
+    axes[5].set_visible(False)
+    fig.suptitle(suptitle, fontsize=13, y=1.01)
+    _save(fig, subdir, fname, tight=True)
+
+
+def _save(fig, subdir, fname, tight=False):
+    plt.tight_layout()
+    out = os.path.join(PLOTS_DIR, subdir, fname)
+    os.makedirs(os.path.dirname(out), exist_ok=True)
+    fig.savefig(out, dpi=150, bbox_inches="tight" if tight else None)
+    plt.close(fig)
+    print(f"Saved → {out}")
+
+
+def main():
+    ap = argparse.ArgumentParser(description="Combined Figures 3, S1 & S4.")
+    ap.add_argument("--only", choices=["3", "s1", "s4"],
+                    help="Render only one figure (still computes its curves).")
+    ap.add_argument("--from-cache", action="store_true",
+                    help="Re-render from cache without computing.")
+    ap.add_argument("--force", action="store_true",
+                    help="Recompute all curves, ignoring the cache.")
+    ap.add_argument("--n-jobs", type=int, default=6)
+    ap.add_argument("--manifest", default=MANIFEST)
+    ap.add_argument("--tag", default="")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    # Which splits are needed for the requested figure(s)?
+    need_image = args.only in (None, "3", "s1")
+    need_patient = args.only in (None, "s4")
+    splits = (["image"] if need_image else []) + (["patient"] if need_patient else [])
+
+    cache = {} if args.force else load_cache()
+    if args.from_cache:
+        print(f"Loaded {len(cache)} cached curves ← {DATA_PATH}")
+    else:
+        cache = compute(cache, splits, args.n_jobs, args.manifest)
+
+    if args.only in (None, "3"):
+        render_fig3(cache, tag)
+    if args.only in (None, "s1"):
+        render_grid(cache, "image", "#D62728", (90, 100),
+                    "Figure S1 — Image-level CV accuracy vs number of selected "
+                    "features", "figure_s1", f"figure_s1{tag}.png")
+    if args.only in (None, "s4"):
+        render_grid(cache, "patient", "#2C7BB6", (82, 95),
+                    "Figure S4 — Patient-level CV accuracy vs number of selected "
+                    "features (C=10)", "figure_s4", f"figure_s4{tag}.png")
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/figure4_6.py

@@ -0,0 +1,122 @@
+#!/usr/bin/env python3
+"""Combined 20-seed classification for Figures 4 and 6 — one compute pass.
+
+Usage:
+    python scripts/visualizations/figure4_6.py              # both figures
+    python scripts/visualizations/figure4_6.py --only 4     # Fig 4 only
+    python scripts/visualizations/figure4_6.py --only 6     # Fig 6 only
+    python scripts/visualizations/figure4_6.py --from-cache # use saved data
+"""
+
+import os, sys, json, argparse
+import numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib.patches import Patch
+from tqdm import tqdm
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from classes import PatientLeakageClassifier
+
+PLOTS_DIR = os.path.join(ROOT, "plots")
+RESULTS_DIR = os.path.join(ROOT, "results")
+# Shared classification-runs cache, owned by scripts/classification.py. Both
+# scripts fill it incrementally (same model/seed/split keys), merge-safe.
+CACHE_DIR = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "cache")
+MODELS = ["VGG16", "DenseNet121", "EfficientNetB1", "MobileNetV2", "ResNet50"]
+SEEDS = list(range(1, 21))
+DATA_PATH = os.path.join(CACHE_DIR, "classification_runs.json")
+MANIFEST = os.path.join(RESULTS_DIR, "simple_patient_manifest.csv")
+
+ap = argparse.ArgumentParser()
+ap.add_argument("--only", default=None, choices=["4", "6"])
+ap.add_argument("--from-cache", action="store_true")
+ap.add_argument("--force", action="store_true", help="Recompute all runs, ignoring cache.")
+ap.add_argument("--manifest", default=MANIFEST)
+ap.add_argument("--tag", default="")
+args = ap.parse_args()
+tag = f"_{args.tag}" if args.tag else ""
+
+# Run or load. The cache is a flat list of per-run dicts; compute only the
+# (model, seed, split) combos it's missing, so dropping a stale model's rows
+# (e.g. ResNet50 after a feature change) recomputes just that model.
+all_data = []
+if os.path.exists(DATA_PATH) and not args.force:
+    with open(DATA_PATH) as f: all_data = json.load(f)
+
+if args.from_cache:
+    print(f"Loaded {len(all_data)} cached runs ← {DATA_PATH}")
+else:
+    have = {(r["model"], r["seed"], r["split_type"]) for r in all_data}
+    todo = [(m, s, st) for m in MODELS for s in SEEDS
+            for st in ("image", "patient") if (m, s, st) not in have]
+    if todo:
+        clf = PatientLeakageClassifier(args.manifest, os.path.join(ROOT, "features"), n_jobs=6)
+        for m, s, st in tqdm(todo, desc="Computing missing runs"):
+            all_data.append(clf.run(m, s, st))
+        os.makedirs(CACHE_DIR, exist_ok=True)
+        with open(DATA_PATH, "w") as f: json.dump(all_data, f, indent=2)
+        print(f"Computed {len(todo)} missing runs; cache now {len(all_data)}")
+    else:
+        print(f"Cache complete ({len(all_data)} runs); nothing to compute.")
+
+# Build accs. Only plot models that actually have data, so a partially filled
+# cache (e.g. ResNet50 dropped pending recompute) still renders without error.
+accs = {m: {"image": [], "patient": []} for m in MODELS}
+for r in all_data:
+    if r["model"] in accs: accs[r["model"]][r["split_type"]].append(r["test"])
+PLOT_MODELS = [m for m in MODELS if accs[m]["image"] or accs[m]["patient"]]
+missing = [m for m in MODELS if m not in PLOT_MODELS]
+if missing:
+    print(f"WARNING: no cached runs for {missing}; run without --from-cache "
+          "to compute them. Plotting remaining models only.")
+
+# Figure 4: image-level only
+if args.only is None or args.only == "4":
+    fig, ax = plt.subplots(figsize=(8, 5))
+    pos = list(range(1, len(PLOT_MODELS) + 1))
+    bp = ax.boxplot([accs[m]["image"] for m in PLOT_MODELS], positions=pos,
+                    widths=0.5, patch_artist=True, showfliers=True,
+                    flierprops=dict(marker='o', markersize=3))
+    for i, b in enumerate(bp['boxes']):
+        b.set_facecolor('#4C9BD4')
+        ax.annotate(f"{np.median(accs[PLOT_MODELS[i]]['image']):.3f}",
+                    (pos[i], np.median(accs[PLOT_MODELS[i]]['image'])),
+                    fontsize=6, ha='center', va='bottom')
+    ax.set_xticks(pos); ax.set_xticklabels(PLOT_MODELS)
+    ax.set_ylabel("Test accuracy"); ax.set_ylim(0.96, 1.00); ax.grid(axis='y', alpha=0.3)
+    ax.set_title("Figure 4 — Image-level test accuracy across 20 seeds", fontsize=12)
+    plt.tight_layout()
+    out = os.path.join(PLOTS_DIR, "figure4", f"figure4{tag}.png")
+    os.makedirs(os.path.dirname(out), exist_ok=True)
+    plt.savefig(out, dpi=150); plt.close()
+    print(f"Saved → {out}")
+
+# Figure 6: image vs patient
+if args.only is None or args.only == "6":
+    fig, ax = plt.subplots(figsize=(10, 6))
+    for i, name in enumerate(PLOT_MODELS):
+        for pos, stype, color in [(i*2+0.7, "image", '#4C9BD4'),
+                                   (i*2+1.3, "patient", '#6DBF6D')]:
+            data = accs[name][stype]
+            bp = ax.boxplot(data, positions=[pos], widths=0.5,
+                            patch_artist=True, showfliers=True,
+                            flierprops=dict(marker='o', markersize=3))
+            bp['boxes'][0].set_facecolor(color)
+            ax.annotate(f"{np.median(data):.3f}", (pos, np.median(data)),
+                        fontsize=6, ha='center', va='bottom')
+    ax.legend(handles=[Patch(facecolor='#4C9BD4', label='Image-level split'),
+                       Patch(facecolor='#6DBF6D', label='Patient-level split')],
+              loc='lower right')
+    ax.set_xticks([p+1 for p in range(0, len(PLOT_MODELS)*2, 2)])
+    ax.set_xticklabels(PLOT_MODELS)
+    ax.set_ylabel("Test accuracy"); ax.set_ylim(0.70, 1.00); ax.grid(axis='y', alpha=0.3)
+    ax.set_title("Figure 6 — Image vs Patient-level test accuracy (20 seeds)", fontsize=13)
+    plt.tight_layout()
+    out = os.path.join(PLOTS_DIR, "figure6", f"figure6{tag}.png")
+    os.makedirs(os.path.dirname(out), exist_ok=True)
+    plt.savefig(out, dpi=150); plt.close()
+    print(f"Saved → {out}")
+
+print("DONE")

scripts/visualizations/figure5.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+"""figure5.py — PCA and t-SNE of VGG16 features, colored by class.
+
+Usage:
+    python scripts/visualizations/figure5.py
+    python scripts/visualizations/figure5.py --tag v2
+"""
+
+import os, sys, argparse
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.decomposition import PCA
+from sklearn.manifold import TSNE
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+FEATURES_DIR = os.path.join(ROOT, "features")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+SEED = 42
+CLASS_COLORS = {"Bengin cases": "#2CA02C", "Malignant cases": "#D62728",
+                "Normal cases": "#1F77B4"}
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--tag", default="", help="Append tag to filename")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    data = np.load(os.path.join(FEATURES_DIR, "VGG16_features.npz"),
+                   allow_pickle=True)
+    X, Y = data["X"], data["Y"]
+    pca50 = PCA(n_components=50, random_state=SEED).fit_transform(X)
+    pca2 = PCA(n_components=2, random_state=SEED).fit_transform(pca50)
+    tsne = TSNE(n_components=2, random_state=SEED, perplexity=30,
+                max_iter=1000).fit_transform(pca50)
+
+    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
+    for ax, coords, title in [(axes[0], pca2, "(a) PCA"),
+                               (axes[1], tsne, "(b) t-SNE")]:
+        for cls in sorted(np.unique(Y)):
+            mask = Y == cls
+            label = cls.replace("Bengin cases", "Benign")
+            ax.scatter(coords[mask, 0], coords[mask, 1],
+                       c=CLASS_COLORS[cls], label=label,
+                       alpha=0.6, s=15, edgecolors="none")
+        ax.set_title(title, fontsize=12)
+        ax.set_xlabel("Component 1")
+        ax.set_ylabel("Component 2")
+        ax.legend(markerscale=2, fontsize=9)
+        ax.grid(alpha=0.2)
+
+    fig.suptitle("Figure 5 — PCA and t-SNE of VGG16 Features (by class)",
+                 fontsize=13)
+    plt.tight_layout()
+    out = os.path.join(PLOTS_DIR, "figure5", f"figure5{tag}.png")
+    os.makedirs(os.path.dirname(out), exist_ok=True)
+    plt.savefig(out, dpi=150)
+    plt.close()
+    print(f"Saved → {out}")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/figure_s2.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+"""figure_s2.py — PCA and t-SNE for all 5 CNN models, colored by class.
+Usage: python scripts/visualizations/figure_s2.py [--tag TAG]"""
+
+import os, sys, argparse
+import numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.decomposition import PCA
+from sklearn.manifold import TSNE
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+FEATURES_DIR = os.path.join(ROOT, "features")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+SEED = 42
+MODELS = ["VGG16", "MobileNetV2", "DenseNet121", "ResNet50", "EfficientNetB1"]
+CLASS_COLORS = {"Bengin cases": "#2CA02C", "Malignant cases": "#D62728", "Normal cases": "#1F77B4"}
+
+ap = argparse.ArgumentParser()
+ap.add_argument("--tag", default="", help="Append tag to filename")
+args = ap.parse_args()
+tag = f"_{args.tag}" if args.tag else ""
+
+fig, axes = plt.subplots(2, 5, figsize=(22, 9))
+for mi, model_name in enumerate(MODELS):
+    data = np.load(os.path.join(FEATURES_DIR, f"{model_name}_features.npz"), allow_pickle=True)
+    X, Y = data["X"], data["Y"]
+    print(f"{model_name}: {X.shape}")
+    pca50 = PCA(n_components=50, random_state=SEED).fit_transform(X)
+    pca2 = PCA(n_components=2, random_state=SEED).fit_transform(pca50)
+    tsne = TSNE(n_components=2, random_state=SEED, perplexity=30, max_iter=800).fit_transform(pca50)
+    for row, coords, title in [(0, pca2, f"{model_name} PCA"), (1, tsne, f"{model_name} t-SNE")]:
+        ax = axes[row, mi]
+        for cls in sorted(np.unique(Y)):
+            mask = Y == cls; label = cls.replace("Bengin cases", "Benign")
+            ax.scatter(coords[mask,0], coords[mask,1], c=CLASS_COLORS[cls], label=label, alpha=0.5, s=8, edgecolors="none")
+        ax.set_title(title, fontsize=9)
+        ax.set_xlabel("Component 1" if row==0 else "t-SNE 1")
+        ax.set_ylabel("Component 2" if row==0 else "t-SNE 2")
+        ax.grid(alpha=0.2)
+
+handles = [plt.Line2D([0],[0], marker='o', color='w', markerfacecolor=c, markersize=8, label=l)
+           for l,c in zip(["Benign","Malignant","Normal"], ["#2CA02C","#D62728","#1F77B4"])]
+fig.legend(handles=handles, loc='lower center', ncol=3, fontsize=10, bbox_to_anchor=(0.5,-0.02))
+fig.suptitle("Figure S2 — PCA (top) and t-SNE (bottom) per model", fontsize=13, y=1.01)
+plt.tight_layout()
+out = os.path.join(PLOTS_DIR, "figure_s2", f"figure_s2{tag}.png")
+os.makedirs(os.path.dirname(out), exist_ok=True)
+plt.savefig(out, dpi=150, bbox_inches="tight")
+plt.close()
+print(f"Saved → {out}")

scripts/visualizations/figure_s3.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+"""figure_s3.py — Example CT images from estimated patient clusters.
+Usage: python scripts/visualizations/figure_s3.py [--tag TAG]"""
+
+import os, sys, csv, argparse
+import numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from PIL import Image
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+MANIFEST = os.path.join(ROOT, "results", "simple_patient_manifest.csv")
+DATASET = os.path.join(os.path.dirname(ROOT), "The IQ-OTHNCCD lung cancer dataset")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+N_EXAMPLES = 5
+CLASS_MAP = {"Benign": ("Bengin cases", "B"), "Malignant": ("Malignant cases", "M"), "Normal": ("Normal cases", "N")}
+
+def f2n_back(fname):
+    prefix = fname[0]; num = int(fname.split("_")[1])
+    cls = {"B": ("Bengin cases", "Bengin"), "M": ("Malignant cases", "Malignant"), "N": ("Normal cases", "Normal")}[prefix]
+    return cls[0], f"{cls[1]} case ({num}).jpg"
+
+ap = argparse.ArgumentParser()
+ap.add_argument("--tag", default="", help="Append tag to filename")
+args = ap.parse_args()
+tag = f"_{args.tag}" if args.tag else ""
+
+patients = {}
+with open(MANIFEST, newline="") as f:
+    reader = csv.DictReader(f)
+    img_col = "confirmed_images" if "confirmed_images" in reader.fieldnames else "images"
+    for row in reader:
+        imgs = row[img_col].split(";")
+        if imgs: patients[row["patient_id"]] = imgs
+
+fig, axes = plt.subplots(3, N_EXAMPLES, figsize=(12, 8))
+for row, (cls_label, (cls_dir, _)) in enumerate(CLASS_MAP.items()):
+    cls_patients = [(p, imgs) for p, imgs in patients.items() if p.lower().startswith(cls_label.lower()) and len(imgs) >= N_EXAMPLES]
+    if not cls_patients: continue
+    pid, imgs = cls_patients[0]
+    for col in range(N_EXAMPLES):
+        ax = axes[row, col]
+        try:
+            cls_dir_name, orig_fname = f2n_back(imgs[col])
+            img = Image.open(os.path.join(DATASET, cls_dir_name, orig_fname)).convert("L")
+            ax.imshow(img, cmap="gray")
+        except Exception as e:
+            ax.text(0.5, 0.5, f"error: {e}", ha="center", va="center", fontsize=7)
+        ax.set_xticks([]); ax.set_yticks([])
+        if col == 0: ax.set_ylabel(f"{cls_label}\nPatient {pid}", fontsize=9, rotation=0, labelpad=40, va="center")
+
+fig.suptitle("Figure S3 — Example CT images from estimated patient clusters", fontsize=12, y=1.01)
+plt.tight_layout()
+out = os.path.join(PLOTS_DIR, "figure_s3", f"figure_s3{tag}.png")
+os.makedirs(os.path.dirname(out), exist_ok=True)
+plt.savefig(out, dpi=150)
+plt.close()
+print(f"Saved → {out}")

scripts/visualizations/figure_s5.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+"""figure_s5.py — Confusion matrix for VGG16 patient-level classification.
+
+Reuses PatientLeakageClassifier.run(return_predictions=True) so the RF ranking,
+gamma grid, and final fit are not duplicated here — the same code path that
+produces the Figure 4/6 numbers also produces these predictions.
+
+Usage: python scripts/visualizations/figure_s5.py [--tag TAG]
+"""
+
+import os
+import sys
+import argparse
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
+
+os.environ.setdefault("OMP_NUM_THREADS", "1")
+os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
+os.environ.setdefault("MKL_NUM_THREADS", "1")
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+PLOTS_DIR = os.path.join(ROOT, "plots")
+
+from classes import PatientLeakageClassifier
+
+SEED = 20
+N_JOBS = 6
+CLASS_NAMES = ["Benign", "Malignant", "Normal"]
+
+
+def display_name(raw):
+    """Map a raw dataset label ('Bengin cases', ...) to a display class name."""
+    if raw.startswith("Bengin"):
+        return "Benign"
+    if raw.startswith("Malignant"):
+        return "Malignant"
+    return "Normal"
+
+
+ap = argparse.ArgumentParser()
+ap.add_argument("--tag", default="")
+args = ap.parse_args()
+tag = f"_{args.tag}" if args.tag else ""
+
+clf = PatientLeakageClassifier(
+    os.path.join(ROOT, "results", "simple_patient_manifest.csv"),
+    os.path.join(ROOT, "features"), n_jobs=N_JOBS)
+r = clf.run("VGG16", SEED, "patient", return_predictions=True)
+print(f"Best: n={r['nfeat']}, gamma={r['gamma']:.6e}, "
+      f"CV={r['cv']:.4f}, Test={r['test']:.4f}")
+
+y_true = [display_name(c) for c in r["y_true"]]
+y_pred = [display_name(c) for c in r["y_pred"]]
+cm = confusion_matrix(y_true, y_pred, labels=CLASS_NAMES, normalize="true")
+
+fig, ax = plt.subplots(figsize=(6, 5))
+ConfusionMatrixDisplay(cm, display_labels=CLASS_NAMES).plot(
+    cmap="Blues", ax=ax, colorbar=True, values_format=".2f")
+ax.set_title("Figure S5 — Patient-level Confusion Matrix (VGG16)", fontsize=12)
+plt.tight_layout()
+out = os.path.join(PLOTS_DIR, "figure_s5", f"figure_s5{tag}.png")
+os.makedirs(os.path.dirname(out), exist_ok=True)
+plt.savefig(out, dpi=150)
+plt.close()
+print(f"Saved → {out}")

scripts/visualizations/manifest_html.py

@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+"""
+manifest_html.py — Generate HTML pages for visually inspecting patient manifests.
+
+One HTML file per class, each patient's assigned images shown in a row.
+
+Usage:
+    python scripts/visualizations/manifest_html.py
+    python scripts/visualizations/manifest_html.py --method siamese
+    python scripts/visualizations/manifest_html.py --manifest results/simple_patient_manifest.csv --method pca50
+"""
+
+import os, sys, csv, argparse, base64
+from io import BytesIO
+from PIL import Image
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__))))
+
+DATASET = os.path.join(os.path.dirname(ROOT), "The IQ-OTHNCCD lung cancer dataset")
+RESULTS_DIR = os.path.join(ROOT, "results")
+HTML_DIR = os.path.join(ROOT, "plots", "html")
+
+MANIFESTS = {
+    "pca50":    os.path.join(RESULTS_DIR, "simple_patient_manifest.csv"),
+    "thumbnail": os.path.join(RESULTS_DIR, "thumbnail_patient_manifest.csv"),
+    "siamese":  os.path.join(RESULTS_DIR, "siamese_manifest.csv"),
+}
+
+THUMB_SIZE = 150  # px, display width
+
+
+def f2n_back(short_name):
+    """B_009 → ('Bengin cases', 'Bengin case (9).jpg')"""
+    prefix = short_name[0]
+    num = int(short_name.split("_")[1])
+    cls_map = {"B": ("Bengin cases", "Bengin"),
+               "M": ("Malignant cases", "Malignant"),
+               "N": ("Normal cases", "Normal")}
+    dir_name, file_prefix = cls_map[prefix]
+    return dir_name, f"{file_prefix} case ({num}).jpg"
+
+
+def img_to_b64(path, size=THUMB_SIZE):
+    """Load an image and return a base64 data URI."""
+    try:
+        img = Image.open(path).convert("L")
+        img.thumbnail((size, size), Image.LANCZOS)
+        buf = BytesIO()
+        img.save(buf, format="PNG")
+        return f"data:image/png;base64,{base64.b64encode(buf.getvalue()).decode()}"
+    except Exception:
+        return ""
+
+
+def build_html(manifest_path, method, class_name, class_dir, patients):
+    """Generate an HTML string for one class."""
+    rows = []
+    for pid, short_names in patients.items():
+        # Build image cards
+        cards = []
+        for sn in short_names:
+            dir_name, fname = f2n_back(sn)
+            img_path = os.path.join(DATASET, dir_name, fname)
+            b64 = img_to_b64(img_path)
+            if b64:
+                cards.append(
+                    f'<div class="card">'
+                    f'<img src="{b64}" alt="{sn}">'
+                    f'<div class="label">{sn}</div>'
+                    f'</div>')
+        if cards:
+            rows.append(
+                f'<div class="patient">'
+                f'<h3>{pid} <span class="count">({len(cards)} images)</span></h3>'
+                f'<div class="images">{"".join(cards)}</div>'
+                f'</div>')
+
+    return f"""<!DOCTYPE html>
+<html lang="en">
+<head>
+<meta charset="UTF-8">
+<title>{method} — {class_name}</title>
+<style>
+  body {{ font-family: -apple-system, sans-serif; background: #1a1a2e; color: #eee; margin: 20px; }}
+  h1 {{ color: #e94560; }}
+  .patient {{ margin-bottom: 30px; border-bottom: 1px solid #333; padding-bottom: 15px; }}
+  .patient h3 {{ margin: 0 0 8px 0; color: #0f3460; background: #16213e; padding: 6px 12px; border-radius: 4px; display: inline-block; }}
+  .count {{ font-weight: normal; color: #888; font-size: 0.85em; }}
+  .images {{ display: flex; flex-wrap: wrap; gap: 8px; }}
+  .card {{ background: #16213e; border-radius: 4px; overflow: hidden; width: {THUMB_SIZE + 20}px; }}
+  .card img {{ display: block; width: {THUMB_SIZE}px; height: {THUMB_SIZE}px; object-fit: contain; margin: 0 auto; background: #000; }}
+  .label {{ font-size: 9px; color: #aaa; text-align: center; padding: 4px; word-break: break-all; }}
+  a.nav {{ color: #e94560; margin-right: 15px; }}
+</style>
+</head>
+<body>
+<h1>{method} — {class_name} <small>({len(patients)} patients, {sum(len(v) for v in patients.values())} images)</small></h1>
+<p>
+  <a class="nav" href="{method}_Benign.html">Benign</a>
+  <a class="nav" href="{method}_Malignant.html">Malignant</a>
+  <a class="nav" href="{method}_Normal.html">Normal</a>
+</p>
+{"".join(rows)}
+</body>
+</html>"""
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--method", default=None,
+                    help="Which manifest method (pca50, thumbnail, siamese). "
+                         "Default: all.")
+    ap.add_argument("--manifest", default=None,
+                    help="Path to manifest CSV (overrides --method).")
+    args = ap.parse_args()
+
+    os.makedirs(HTML_DIR, exist_ok=True)
+
+    methods_to_run = [args.method] if args.method else list(MANIFESTS.keys())
+
+    for method in methods_to_run:
+        manifest_path = args.manifest or MANIFESTS.get(method)
+        if not manifest_path or not os.path.exists(manifest_path):
+            print(f"  {method}: manifest not found ({manifest_path})")
+            continue
+
+        # Load manifest, group by class
+        patients_by_class = {"Benign": {}, "Malignant": {}, "Normal": {}}
+        with open(manifest_path, newline="") as f:
+            reader = csv.DictReader(f)
+            img_col = ("images" if "images" in reader.fieldnames
+                       else "confirmed_images")
+            for row in reader:
+                cls = row.get("class", "").strip()
+                # Normalize class names
+                cls_lower = cls.lower()
+                if cls_lower in ("benign", "bengin"):
+                    cls = "Benign"
+                elif cls_lower in ("malignant", "malig"):
+                    cls = "Malignant"
+                elif cls_lower == "normal":
+                    cls = "Normal"
+                elif cls_lower == "unknown":
+                    # Infer from patient_id prefix
+                    pid = row.get("patient_id", "")
+                    if pid.lower().startswith("benign") or pid.lower().startswith("bengin"):
+                        cls = "Benign"
+                    elif pid.lower().startswith("malignant") or pid.lower().startswith("malig"):
+                        cls = "Malignant"
+                    elif pid.lower().startswith("normal"):
+                        cls = "Normal"
+                if cls not in patients_by_class:
+                    continue
+                pid = row["patient_id"]
+                imgs = row[img_col].split(";") if row[img_col] else []
+                if imgs:
+                    patients_by_class[cls][pid] = imgs
+
+        for cls_name in ["Benign", "Malignant", "Normal"]:
+            patients = patients_by_class[cls_name]
+            if not patients:
+                print(f"  {method}/{cls_name}: no patients, skipping")
+                continue
+
+            html = build_html(manifest_path, method, cls_name,
+                              "", patients)
+            out_path = os.path.join(HTML_DIR, f"{method}_{cls_name}.html")
+            with open(out_path, "w") as f:
+                f.write(html)
+            print(f"  Saved → plots/html/{method}_{cls_name}.html "
+                  f"({len(patients)} patients)")
+
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/manifest_tsne.py

@@ -0,0 +1,140 @@
+#!/usr/bin/env python3
+"""manifest_tsne.py — visualize any patient-grouping manifest on the VGG16 t-SNE.
+
+Plots the same VGG16 feature t-SNE used elsewhere, but colored by the groups in
+a given manifest (siamese / pca50 / thumbnail). Because every manifest is drawn
+on the *identical* layout (same features, PCA, seed, perplexity), the resulting
+per-class figures are directly comparable across methods.
+
+For the siamese grouping this is diagnostic: if a siamese "patient" is a coherent
+patient it forms a tight island; if it's an over-merged chain, its color is
+smeared across feature space (VGG16 sees images the siamese wrongly linked).
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/visualizations/manifest_tsne.py \
+        --manifest results/siamese_manifest.csv --name siamese
+    python scripts/visualizations/manifest_tsne.py \
+        --manifest results/simple_patient_manifest.csv --name pca50
+"""
+
+import os
+import sys
+import re
+import csv
+import argparse
+from collections import defaultdict
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+from sklearn.decomposition import PCA
+from sklearn.manifold import TSNE
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+FEATURES_DIR = os.path.join(ROOT, "features")
+PLOTS_DIR = os.path.join(ROOT, "plots")
+
+RANDOM_STATE = 42
+CLASS_OF = {"B": "Benign", "M": "Malignant", "N": "Normal"}
+
+
+def f2n(fname):
+    """VGG16 feature filename -> manifest short name (e.g. 'B_001')."""
+    m = re.search(r"\((\d+)\)", fname)
+    num = int(m.group(1)) if m else None
+    for cls_key, prefix in [("Bengin cases", "B"), ("Malignant cases", "M"),
+                            ("Normal cases", "N")]:
+        if fname.startswith(cls_key.rstrip("s")):
+            return f"{prefix}_{num:03d}" if num else fname
+    return fname
+
+
+def load_manifest(path):
+    """Return {image_short_name: group_id}."""
+    mapping = {}
+    with open(path, newline="") as f:
+        for row in csv.DictReader(f):
+            for img in row["images"].split(";"):
+                if img:
+                    mapping[img] = row["patient_id"]
+    return mapping
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--manifest", required=True, help="Path to a grouping manifest CSV.")
+    ap.add_argument("--name", required=True, help="Short method name for titles/filenames.")
+    ap.add_argument("--features", default="VGG16", help="CNN feature set for the layout.")
+    ap.add_argument("--no-centroids", dest="centroids", action="store_false",
+                    help="Do not draw per-group centroid diamonds.")
+    ap.add_argument("--tag", default="")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    print(f"Loading {args.features} features ...")
+    data = np.load(os.path.join(FEATURES_DIR, f"{args.features}_features.npz"),
+                   allow_pickle=True)
+    X, filenames = data["X"], data["filenames"]
+    img_ids = np.array([f2n(f) for f in filenames])
+    img_class = np.array([CLASS_OF.get(i.split("_")[0], "?") for i in img_ids])
+
+    print("Computing PCA-50 + t-SNE (shared layout) ...")
+    n_pca = min(50, X.shape[0] - 1, X.shape[1])
+    X_pca = PCA(n_components=n_pca, random_state=RANDOM_STATE).fit_transform(X)
+    X_tsne = TSNE(n_components=2, perplexity=35, learning_rate="auto",
+                  init="pca", random_state=RANDOM_STATE).fit_transform(X_pca)
+
+    groups = load_manifest(args.manifest)
+    img_group = np.array([groups.get(i, "unassigned") for i in img_ids])
+
+    for cls in ["Benign", "Malignant", "Normal"]:
+        fig, ax = plt.subplots(figsize=(14, 10))
+        ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c="lightgray", s=3, alpha=0.15)
+
+        mask = img_class == cls
+        class_groups = sorted(set(img_group[mask]))
+        n = len(class_groups)
+        # Order groups by size so the biggest (most likely over-merged) is obvious.
+        sizes = {g: int((img_group[mask] == g).sum()) for g in class_groups}
+        class_groups = sorted(class_groups, key=lambda g: -sizes[g])
+        cmap = plt.cm.tab20 if n <= 20 else plt.cm.gist_ncar
+
+        for gi, g in enumerate(class_groups):
+            color = cmap(gi % 20) if n <= 20 else cmap(gi / max(n - 1, 1))
+            gm = mask & (img_group == g)
+            ax.scatter(X_tsne[gm, 0], X_tsne[gm, 1], c=[color], s=18, alpha=0.8)
+            if args.centroids:
+                cx, cy = X_tsne[gm, 0].mean(), X_tsne[gm, 1].mean()
+                ax.scatter(cx, cy, c=[color], s=60, marker="D", edgecolors="black",
+                           linewidths=0.6, zorder=5)
+
+        biggest = class_groups[0]
+        legend_handles = [
+            Line2D([0], [0], marker="o", color="w", markerfacecolor="gray",
+                   markersize=8, label="Group images (dots)"),
+        ]
+        if args.centroids:
+            legend_handles.append(
+                Line2D([0], [0], marker="D", color="w", markerfacecolor="gray",
+                       markersize=8, label="Group centroids (diamonds)"))
+        ax.legend(handles=legend_handles, loc="lower right")
+        ax.set_title(f"{cls} — {args.name} groups on {args.features} t-SNE "
+                     f"({n} groups; largest={sizes[biggest]} imgs)")
+        ax.set_xlabel("t-SNE dim 1")
+        ax.set_ylabel("t-SNE dim 2")
+        plt.tight_layout()
+        out = os.path.join(PLOTS_DIR, "tsne", f"tsne_{args.name}_{cls}{tag}.png")
+        os.makedirs(os.path.dirname(out), exist_ok=True)
+        plt.savefig(out, dpi=150)
+        plt.close()
+        print(f"  Saved {out}")
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/siamese_similarity_tsne.py

@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+"""siamese_similarity_tsne.py — t-SNE of the siamese's OWN similarity space.
+
+Unlike manifest_tsne.py (which recolors the VGG16 feature layout), this builds
+the layout directly from the siamese pairwise distance 1 - P(same-patient), so
+proximity reflects how the siamese model itself relates images. Points are
+colored by the siamese edge-rank groups.
+
+This is the diagnostic view for the over-merge: if the siamese collapses several
+patients together (over-confidence on IQ-OTH), the largest group forms one dense
+mass in its own space; coherent patients form tight, separated islands.
+
+Usage:
+    conda activate fundus_imaging
+    python scripts/visualizations/siamese_similarity_tsne.py
+"""
+
+import os
+import sys
+import re
+import csv
+import argparse
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.manifold import TSNE
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.dirname(
+    os.path.abspath(__file__)))))
+ROOT = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+PLOTS_DIR = os.path.join(ROOT, "plots")
+DEFAULT_DATASET = os.path.join(os.path.dirname(ROOT),
+                               "The IQ-OTHNCCD lung cancer dataset")
+
+from classes import SiamesePatientMatcher
+
+CLASS_DIRS = {"Bengin cases": "Benign", "Malignant cases": "Malignant",
+              "Normal cases": "Normal"}
+VALID_EXT = (".png", ".jpg", ".jpeg", ".tif", ".tiff", ".bmp")
+RANDOM_STATE = 42
+
+
+def f2n(fname):
+    m = re.search(r"\((\d+)\)", fname)
+    num = int(m.group(1)) if m else None
+    for cls_key, prefix in [("Bengin cases", "B"), ("Malignant cases", "M"),
+                            ("Normal cases", "N")]:
+        if fname.startswith(cls_key.rstrip("s")):
+            return f"{prefix}_{num:03d}" if num else fname
+    return fname
+
+
+def load_manifest(path):
+    mapping = {}
+    with open(path, newline="") as f:
+        for row in csv.DictReader(f):
+            for img in row["images"].split(";"):
+                if img:
+                    mapping[img] = row["patient_id"]
+    return mapping
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--model", default=os.path.join(ROOT, "models", "siamese_resnet18.pt"))
+    ap.add_argument("--backbone", default="resnet18")
+    ap.add_argument("--dataset", default=DEFAULT_DATASET)
+    ap.add_argument("--manifest", default=os.path.join(ROOT, "results", "siamese_manifest.csv"))
+    ap.add_argument("--name", default="siamese_sim")
+    ap.add_argument("--highlight-largest", type=int, default=0,
+                    help="Grey all points and bold only the N largest groups.")
+    ap.add_argument("--tag", default="")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    matcher = SiamesePatientMatcher(args.model, backbone=args.backbone, input_size=224)
+    groups = load_manifest(args.manifest)
+
+    for class_dir, cls in CLASS_DIRS.items():
+        cpath = os.path.join(args.dataset, class_dir)
+        files = sorted(f for f in os.listdir(cpath) if f.lower().endswith(VALID_EXT))
+        paths = [os.path.join(cpath, f) for f in files]
+        ids = [f2n(f) for f in files]
+        print(f"\n{cls}: {len(paths)} images")
+
+        # Siamese pairwise distance -> t-SNE on precomputed distances.
+        emb = matcher.embed_images(paths)
+        P = matcher._dense_prob_matrix(emb)
+        dist = np.clip(1.0 - P, 0.0, None)
+        np.fill_diagonal(dist, 0.0)
+        perp = max(5, min(30, (len(paths) - 1) // 3))
+        X = TSNE(n_components=2, metric="precomputed", init="random",
+                 perplexity=perp, random_state=RANDOM_STATE).fit_transform(dist)
+
+        img_group = np.array([groups.get(i, "unassigned") for i in ids])
+        g_order = sorted(set(img_group), key=lambda g: -(img_group == g).sum())
+        n = len(g_order)
+
+        fig, ax = plt.subplots(figsize=(14, 10))
+        if args.highlight_largest > 0:
+            # Grey everything, then bold only the N largest groups.
+            ax.scatter(X[:, 0], X[:, 1], c="lightgray", s=12, alpha=0.5)
+            hl = g_order[:args.highlight_largest]
+            hl_cmap = plt.cm.tab10
+            for gi, g in enumerate(hl):
+                gm = img_group == g
+                ax.scatter(X[gm, 0], X[gm, 1], c=[hl_cmap(gi)], s=28, alpha=0.9,
+                           edgecolors="black", linewidths=0.3,
+                           label=f"{g} ({gm.sum()} imgs)")
+            ax.legend(loc="lower right", title="Largest siamese groups")
+            ax.set_title(f"{cls} — siamese similarity t-SNE "
+                         f"(largest {len(hl)} of {n} groups highlighted)")
+        else:
+            cmap = plt.cm.tab20 if n <= 20 else plt.cm.gist_ncar
+            for gi, g in enumerate(g_order):
+                color = cmap(gi % 20) if n <= 20 else cmap(gi / max(n - 1, 1))
+                gm = img_group == g
+                ax.scatter(X[gm, 0], X[gm, 1], c=[color], s=18, alpha=0.8)
+            biggest = (img_group == g_order[0]).sum()
+            ax.set_title(f"{cls} — siamese similarity t-SNE "
+                         f"({n} groups; largest={biggest} imgs)")
+        ax.set_xlabel("t-SNE dim 1 (siamese distance)")
+        ax.set_ylabel("t-SNE dim 2 (siamese distance)")
+        plt.tight_layout()
+        suffix = "_highlight" if args.highlight_largest > 0 else ""
+        out = os.path.join(PLOTS_DIR, "tsne", f"tsne_{args.name}_{cls}{suffix}{tag}.png")
+        os.makedirs(os.path.dirname(out), exist_ok=True)
+        plt.savefig(out, dpi=150)
+        plt.close()
+        print(f"  Saved {out}")
+    print("DONE")
+
+
+if __name__ == "__main__":
+    main()

scripts/visualizations/simple_patient_tsne.py

@@ -1,49 +1,74 @@
-#!/usr/bin/env python3
-"""
-simple_patient_tsne.py
 
-The short path:
-  1. Load VGG16 features
-  2. K-means per class (15/40/55 patients) in PCA-50d space
-  3. Optional: iterative centroid refinement
-  4. Plot t-SNE colored by cluster, with centroid labels
 
-Usage:
-    conda activate fundus_imaging
-    python scripts/simple_patient_tsne.py
-"""
 
 import os, sys, re, csv, json
 from collections import defaultdict
 import numpy as np
 import matplotlib
-matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 from sklearn.decomposition import PCA
 from sklearn.manifold import TSNE
 from sklearn.cluster import KMeans
 
-sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+import argparse
+#!/usr/bin/env python3
+"""
+simple_patient_tsne.py
+The short path:
+  1. Load VGG16 features
+  2. K-means per class (15/40/55 patients) in PCA-50d space
+  3. Optional: iterative centroid refinement
+  4. Plot t-SNE colored by cluster, with centroid labels
+Usage:
+    conda activate fundus_imaging
+    python scripts/simple_patient_tsne.py
+"""
+matplotlib.use("Agg")
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
 PROJECT_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 FEATURES_DIR = os.path.join(PROJECT_ROOT, "features")
 PLOTS_DIR = os.path.join(PROJECT_ROOT, "plots")
 RESULTS_DIR = os.path.join(PROJECT_ROOT, "results")
 os.makedirs(PLOTS_DIR, exist_ok=True)
 os.makedirs(RESULTS_DIR, exist_ok=True)
-
 PATIENT_COUNTS = {"Bengin cases": 15, "Malignant cases": 40, "Normal cases": 55}
 CLASS_NAMES = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"}
 RANDOM_STATE = 42
-
 # ---------------------------------------------------------------------------
 # 1. Load VGG16 features
 # ---------------------------------------------------------------------------
-
 print("Loading VGG16 features ...")
 data = np.load(os.path.join(FEATURES_DIR, "VGG16_features.npz"), allow_pickle=True)
 X, Y, filenames = data["X"], data["Y"], data["filenames"]
-
 def f2n(fname):
     m = re.search(r'\((\d+)\)', fname)
     num = int(m.group(1)) if m else None
@@ -51,41 +76,31 @@ def f2n(fname):
         if fname.startswith(cls_key.rstrip("s")):
             return f"{prefix}_{num:03d}" if num else fname
     return fname
-
 img_nums = np.array([f2n(f) for f in filenames])
-
 # ---------------------------------------------------------------------------
 # 2. K-means per class in PCA-50d space
 # ---------------------------------------------------------------------------
-
 print("Clustering in PCA-50d space ...")
 n_pca = min(50, X.shape[0] - 1, X.shape[1])
 X_pca = PCA(n_components=n_pca, random_state=RANDOM_STATE).fit_transform(X)
-
 image_to_patient = {}
 patient_to_images = defaultdict(list)
-
 for class_name, k in PATIENT_COUNTS.items():
     mask = Y == class_name
     X_class = X_pca[mask]
     idx_class = np.where(mask)[0]
-
     kmeans = KMeans(n_clusters=k, random_state=RANDOM_STATE, n_init=20)
     labels = kmeans.fit_predict(X_class)
-
     prefix = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"}[class_name]
     for i, cluster_id in enumerate(labels):
         pid = f"{prefix}_{cluster_id:02d}"
         img = img_nums[idx_class[i]]
         image_to_patient[img] = pid
         patient_to_images[pid].append(img)
-
 print(f"  {len(patient_to_images)} patients, {len(image_to_patient)} images")
-
 # ---------------------------------------------------------------------------
 # 3. Iterative centroid refinement (optional, 3 passes)
 # ---------------------------------------------------------------------------
-
 print("Refining assignments (nearest-centroid, 5 passes) ...")
 for iteration in range(5):
     # Compute centroids
@@ -93,7 +108,6 @@ for iteration in range(5):
     for pid, imgs in patient_to_images.items():
         idxs = [np.where(img_nums == img)[0][0] for img in imgs]
         centroids[pid] = X_pca[idxs].mean(axis=0)
-
     # Reassign
     moves = 0
     for class_name in PATIENT_COUNTS:
@@ -118,28 +132,22 @@ for iteration in range(5):
     print(f"  Pass {iteration+1}: {moves} moves")
     if moves == 0:
         break
-
 # ---------------------------------------------------------------------------
 # 4. t-SNE
 # ---------------------------------------------------------------------------
-
 print("Computing t-SNE ...")
 X_tsne = TSNE(n_components=2, perplexity=35, learning_rate="auto",
               init="pca", random_state=RANDOM_STATE).fit_transform(X_pca)
-
 # ---------------------------------------------------------------------------
 # 5. Plot — one figure per class
 # ---------------------------------------------------------------------------
-
 for class_name, display_name in CLASS_NAMES.items():
     fig, ax = plt.subplots(1, 1, figsize=(14, 10))
     ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c="lightgray", s=3, alpha=0.15)
-
     mask = Y == class_name
     class_pids = sorted([p for p in patient_to_images if p.startswith(display_name)])
     n_patients = len(class_pids)
     cmap = plt.cm.tab20 if n_patients <= 20 else plt.cm.gist_ncar
-
     for pi, pid in enumerate(class_pids):
         color = cmap(pi % 20) if n_patients <= 20 else cmap(pi / max(n_patients-1, 1))
         pts_x, pts_y = [], []
@@ -152,7 +160,6 @@ for class_name, display_name in CLASS_NAMES.items():
         cx, cy = np.mean(pts_x), np.mean(pts_y)
         ax.scatter(cx, cy, c=[color], s=60, marker='D', edgecolors='black',
                   linewidths=0.6, zorder=5, label='_nolegend_')
-
     # Legend elements
     from matplotlib.lines import Line2D
     legend_elements = [
@@ -162,31 +169,32 @@ for class_name, display_name in CLASS_NAMES.items():
                markersize=8, label='Patient centroids (diamonds)'),
     ]
     ax.legend(handles=legend_elements, loc='lower right')
-
     ax.set_title(f"{display_name} — VGG16 t-SNE ({n_patients} patients)")
     ax.set_xlabel("t-SNE dim 1")
     ax.set_ylabel("t-SNE dim 2")
     plt.tight_layout()
-    out = os.path.join(PLOTS_DIR, f"tsne_{display_name}.png")
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--tag", default="", help="Append tag to filename")
+    args = ap.parse_args()
+    tag = f"_{args.tag}" if args.tag else ""
+
+    out = os.path.join(PLOTS_DIR, "tsne", f"tsne_{display_name}{tag}.png")
+    os.makedirs(os.path.dirname(out), exist_ok=True)
     plt.savefig(out, dpi=150)
     plt.close()
     print(f"  Saved {out}")
-
 # ---------------------------------------------------------------------------
 # 6. Save assignments
 # ---------------------------------------------------------------------------
-
 manifest = []
 for pid in sorted(patient_to_images.keys()):
     imgs = sorted(patient_to_images[pid])
     manifest.append({"patient_id": pid, "class": pid.split("_")[0],
                      "n_images": len(imgs), "images": ";".join(imgs)})
-
 with open(os.path.join(RESULTS_DIR, "simple_patient_manifest.csv"), "w", newline="") as f:
     w = csv.DictWriter(f, fieldnames=["patient_id", "class", "n_images", "images"])
     w.writeheader()
     w.writerows(manifest)
-
 print(f"\nSaved simple_patient_manifest.csv ({len(manifest)} patients, "
       f"{sum(m['n_images'] for m in manifest)} images)")
 print("DONE")