Add scripts for patient classification and visualization

- Created `classification.py` for comparing image-level and patient-level classification results using various CNN models.
- Implemented `create_patient_groups.py` to extract features, generate PCA/t-SNE plots, and identify patient groups via K-means clustering.
- Added `figure6.py` to generate boxplots for test accuracy across multiple seeds.
- Developed `simple_patient_tsne.py` to perform t-SNE visualization of patient groups and save results in a manifest file.
- Introduced `simple_patient_manifest.csv` to store patient IDs, classes, image counts, and associated images.

.gitignore

@@ -0,0 +1,47 @@
+# Archive
+.archive/
+
+# Python
+__pycache__/
+*.py[cod]
+*.egg-info/
+dist/
+build/
+
+# Environments
+.env
+.venv
+env/
+venv/
+
+# IDE
+.idea/
+.vscode/
+*.swp
+*.swo
+
+# OS
+.DS_Store
+Thumbs.db
+
+# Sensitive
+*.pem
+*.key
+*.p12
+*.pfx
+credentials*
+secrets*
+*.env.*
+
+# Data (large feature files — regenerate via scripts)
+features/*.npz
+features/*.npy
+
+# Outputs
+plots/*.png
+
+# Jupyter
+.ipynb_checkpoints/
+
+# Conda
+.python-version

classes/__init__.py

@@ -0,0 +1,3 @@
+from .features import FeatureExtractor
+from .patient_identifier import PatientIdentifier
+from .classifier import PatientLeakageClassifier

classes/alignment.py

@@ -0,0 +1,199 @@
+"""
+Align patient clusters across models using the Hungarian algorithm.
+
+Each model clusters independently, so its cluster-3 for benign may not
+correspond to another model's cluster-3.  The Hungarian algorithm finds
+the optimal 1-to-1 mapping between cluster labels by maximizing
+image overlap.
+"""
+
+import numpy as np
+from scipy.optimize import linear_sum_assignment
+
+
+class PatientAligner:
+    """Align patient clusters across multiple models.
+
+    Parameters
+    ----------
+    assignments : dict
+        Nested dict:  model_name -> {patient_id: [filenames, ...]}
+        e.g. {"VGG16": {"benign_0": [...], ...}, "DenseNet121": {...}, ...}
+    """
+
+    def __init__(self, assignments):
+        self.assignments = assignments
+        self.model_names = sorted(assignments.keys())
+        self.n_models = len(self.model_names)
+
+    def align_to_reference(self, reference_model=None):
+        """Align every model's clusters to a reference model.
+
+        Parameters
+        ----------
+        reference_model : str or None
+            Model to use as reference.  Defaults to the first model.
+
+        Returns
+        -------
+        mapping : dict
+            model_name -> {model_patient_id: reference_patient_id}
+        """
+        if reference_model is None:
+            reference_model = self.model_names[0]
+
+        print(f"\nAligning all models to reference: {reference_model}")
+        mapping = {reference_model: {pid: pid for pid in self.assignments[reference_model]}}
+
+        for model in self.model_names:
+            if model == reference_model:
+                continue
+            print(f"  {model} -> {reference_model} ...")
+            mapping[model] = {}
+            # Align class by class
+            for prefix in self._class_prefixes():
+                mapping[model].update(
+                    self._align_class(model, reference_model, prefix)
+                )
+
+        return mapping
+
+    def build_consensus(self, mapping):
+        """For each image, measure model agreement via co-clustering overlap.
+
+        Instead of comparing arbitrary cluster IDs (even aligned ones),
+        we ask: *which other images* does each model place in this image's
+        cluster?  High Jaccard overlap = models agree on patient membership
+        regardless of what they call the cluster.
+
+        Returns
+        -------
+        consensus : dict
+            image_filename -> {
+                "consensus_pid": <reference-model patient id>,
+                "co_cluster_agreement": <mean pairwise Jaccard>,
+                "cluster_size": {model: n_images_in_cluster},
+                "raw_pids": {model: original_patient_id},
+            }
+        """
+        # Build cluster lookup: model -> pid -> set of filenames
+        clusters = {}
+        for model, model_assignments in self.assignments.items():
+            clusters[model] = {pid: set(files) for pid, files in model_assignments.items()}
+
+        # Per-image: find which cluster it's in per model
+        image_to_cluster = {}
+        for model, model_clusters in clusters.items():
+            for pid, files in model_clusters.items():
+                for f in files:
+                    image_to_cluster.setdefault(f, {})[model] = pid
+
+        consensus = {}
+        for img, model_pids in image_to_cluster.items():
+            # Co-cluster set for each model (images in same cluster, excluding self)
+            co_sets = {}
+            for model in self.model_names:
+                pid = model_pids.get(model)
+                if pid is not None:
+                    co_sets[model] = clusters[model][pid] - {img}
+                else:
+                    co_sets[model] = set()
+
+            # Pairwise Jaccard between all model pairs
+            jaccards = []
+            for i, ma in enumerate(self.model_names):
+                for j, mb in enumerate(self.model_names):
+                    if i < j:
+                        sa, sb = co_sets[ma], co_sets[mb]
+                        union = len(sa | sb)
+                        jac = len(sa & sb) / union if union > 0 else 1.0
+                        jaccards.append(jac)
+
+            mean_jaccard = np.mean(jaccards) if jaccards else 1.0
+
+            # Consensus PID via aligned majority vote
+            aligned_votes = {}
+            for model, pid in model_pids.items():
+                aligned_votes[model] = mapping.get(model, {}).get(pid, pid)
+            pid_counts = {}
+            for pid in aligned_votes.values():
+                pid_counts[pid] = pid_counts.get(pid, 0) + 1
+            consensus_pid = max(pid_counts, key=pid_counts.get)
+
+            consensus[img] = {
+                "consensus_pid": consensus_pid,
+                "co_cluster_agreement": round(mean_jaccard, 3),
+                "cluster_size": {m: len(co_sets[m]) + 1 for m in self.model_names},
+                "raw_pids": model_pids,
+            }
+
+        return consensus
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    def _class_prefixes(self):
+        """Return the set of class prefixes (benign, malig, normal)."""
+        prefixes = set()
+        for assignments in self.assignments.values():
+            for pid in assignments:
+                prefixes.add(pid.rsplit("_", 1)[0])
+        return sorted(prefixes)
+
+    def _align_class(self, model_a, model_b, prefix):
+        """Hungarian-align patient IDs for one class between two models.
+
+        Returns dict: model_a_patient_id -> model_b_patient_id
+        """
+        # Gather patient IDs for this class prefix
+        pids_a = sorted([p for p in self.assignments[model_a] if p.startswith(prefix)])
+        pids_b = sorted([p for p in self.assignments[model_b] if p.startswith(prefix)])
+
+        n_a, n_b = len(pids_a), len(pids_b)
+        if n_a == 0 or n_b == 0:
+            return {}
+
+        # Build overlap cost matrix (negate overlap → minimize = maximize overlap)
+        cost = np.zeros((n_a, n_b))
+        for i, pa in enumerate(pids_a):
+            files_a = set(self.assignments[model_a][pa])
+            for j, pb in enumerate(pids_b):
+                files_b = set(self.assignments[model_b][pb])
+                overlap = len(files_a & files_b)
+                cost[i, j] = -overlap  # Hungarian *minimizes* cost
+
+        row_ind, col_ind = linear_sum_assignment(cost)
+
+        # Build mapping
+        mapping = {}
+        zero_matches = 0
+        missed_opportunities = 0
+
+        for i, j in zip(row_ind, col_ind):
+            overlap = -cost[i, j]
+            total_a = len(self.assignments[model_a][pids_a[i]])
+
+            # Best possible match for this cluster (even if Hungarian didn't pick it)
+            best_overlap = int(-cost[i].min())  # most negative → best match
+            best_j = int(np.argmin(cost[i]))
+
+            mapping[pids_a[i]] = pids_b[j]
+
+            if overlap == 0:
+                zero_matches += 1
+                if best_overlap > 0:
+                    missed_opportunities += 1
+                    best_pid = pids_b[best_j]
+                    print(f"    LOW-OVERLAP: {pids_a[i]} ({total_a} imgs) matched to "
+                          f"{pids_b[j]} (0 overlap) — best would be {best_pid} "
+                          f"({best_overlap}/{total_a} overlap)")
+                else:
+                    print(f"    UNMATCHED: {pids_a[i]} ({total_a} imgs) — "
+                          f"no cluster in {model_b} shares ANY images with this group")
+
+        if zero_matches:
+            print(f"    {zero_matches} zero-overlap matches ({missed_opportunities} "
+                  f"forced by global optimum, {zero_matches - missed_opportunities} "
+                  f"truly unmatched)")
+        return mapping

classes/classifier.py

@@ -0,0 +1,170 @@
+"""
+PatientLeakageClassifier — runs the RF→GridSearchCV→SVM pipeline
+with patient-aware or image-level train/test splitting.
+"""
+
+import os, csv, re
+import numpy as np
+from sklearn.model_selection import (
+    train_test_split, GridSearchCV, StratifiedGroupKFold
+)
+from sklearn.ensemble import RandomForestClassifier
+from sklearn.svm import SVC
+from sklearn.pipeline import Pipeline
+from sklearn.preprocessing import StandardScaler
+from sklearn.metrics import accuracy_score
+
+
+# ---------------------------------------------------------------------------
+# 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
+
+DEFAULT_NFEATURES = [50, 100, 200, 300, 400, 500, 750, 1000,
+                     1500, 2000, 3000, 5000]
+DEFAULT_GAMMA_RANGE = (-1.5, 1, 7)  # logspace args
+
+
+class PatientLeakageClassifier:
+    """Run the classification pipeline for one model, seed, and split type.
+
+    Parameters
+    ----------
+    manifest_path : str   Path to patient manifest CSV.
+    features_dir : str    Directory with {Model}_features.npz files.
+    n_jobs : int          Parallelism for RF and GridSearchCV.
+    """
+
+    def __init__(self, manifest_path, features_dir, n_jobs=8):
+        self.features_dir = features_dir
+        self.n_jobs = n_jobs
+        self._image_to_patient = {}
+        with open(manifest_path, newline="") as f:
+            reader = csv.DictReader(f)
+            img_col = ("confirmed_images" if "confirmed_images" in reader.fieldnames
+                       else "images")
+            for row in reader:
+                pid = row["patient_id"]
+                for img in row[img_col].split(";"):
+                    if img:
+                        self._image_to_patient[img] = pid
+
+    # ------------------------------------------------------------------
+    # Public API
+    # ------------------------------------------------------------------
+
+    def run(self, model_name, seed, split_type="patient",
+            nfeatures_list=None, gamma_logspace=None):
+        """Run full pipeline: train/test split → RF importance →
+        GridSearchCV over nfeatures×gamma → retrain → test.
+
+        Parameters
+        ----------
+        model_name : str       e.g. "VGG16"
+        seed : int             Random seed for the train/test split.
+        split_type : str       "image" or "patient".
+        nfeatures_list : list  Feature counts to search over.
+        gamma_logspace : tuple Args for np.logspace (start, stop, n).
+
+        Returns
+        -------
+        dict with keys: model, seed, split_type, cv, test, nfeat, gamma,
+                        n_train, n_test, n_train_patients, n_test_patients
+        """
+        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)
+
+        # --- 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"]
+
+        # --- Retrain + test ---
+        final = Pipeline([("scaler", StandardScaler()),
+                          ("svm", SVC(kernel="rbf", C=10, gamma=best_g))])
+        final.fit(X_tr[:, sorted_idx[:best_n]], y_tr)
+        y_pred = final.predict(X_te[:, sorted_idx[:best_n]])
+
+        return {
+            "model": model_name,
+            "seed": seed,
+            "split_type": split_type,
+            "cv": float(best_cv),
+            "test": float(accuracy_score(y_te, y_pred)),
+            "nfeat": best_n,
+            "gamma": float(best_g),
+            "n_train": X_tr.shape[0],
+            "n_test": X_te.shape[0],
+            "n_train_patients": n_tr_patients,
+            "n_test_patients": n_te_patients,
+        }
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    def _load_data(self, model_name):
+        data = np.load(os.path.join(self.features_dir,
+                                     f"{model_name}_features.npz"),
+                       allow_pickle=True)
+        X = data["X"]
+        Y = data["Y"]
+        filenames = data["filenames"]
+        groups = np.array([self._image_to_patient.get(_f2n(f), f"unk_{i}")
+                           for i, f in enumerate(filenames)])
+        return X, Y, groups

classes/cluster_confidence.py

@@ -0,0 +1,95 @@
+"""
+Cluster confidence scoring for patient identification.
+
+For each K-means cluster (in PCA-50d space), compute how well-separated
+it is from its nearest neighbor.  Well-separated clusters are more likely
+to represent real individual patients.
+"""
+
+import numpy as np
+from sklearn.decomposition import PCA
+from sklearn.cluster import KMeans
+from scipy.spatial.distance import cdist
+
+
+def compute_confidence(X, labels, patient_estimates, random_state=42):
+    """Score each patient cluster by its separation from neighbors.
+
+    Works in PCA-50d space (same space K-means operates in).
+
+    Parameters
+    ----------
+    X : np.ndarray  (n_images, n_features)
+        Full feature matrix (all classes together).
+    labels : np.ndarray  (n_images,)
+        Class label for each image.
+    patient_estimates : dict
+        class_name -> n_patients.
+    random_state : int
+
+    Returns
+    -------
+    cluster_scores : dict
+        patient_id -> {
+            "confidence": float,        # separation / internal_spread
+            "intra_spread": float,      # mean distance to centroid
+            "nearest_neighbor_dist": float,
+            "nearest_neighbor_id": str,
+            "images": list of filenames,  # set by caller
+            "n_images": int,
+        }
+    """
+    n_pca = min(50, X.shape[0] - 1, X.shape[1])
+    pca = PCA(n_components=n_pca, random_state=random_state)
+    X_pca = pca.fit_transform(X)
+
+    cluster_scores = {}
+
+    for class_name, n_clusters in patient_estimates.items():
+        idx_class = np.where(labels == class_name)[0]
+        X_class = X_pca[idx_class]
+
+        kmeans = KMeans(n_clusters=n_clusters,
+                        random_state=random_state, n_init=20)
+        cluster_ids = kmeans.fit_predict(X_class)
+        centroids = kmeans.cluster_centers_   # (n_clusters, n_pca)
+
+        # Intra-cluster spread: mean distance from centroid
+        intra_spreads = np.zeros(n_clusters)
+        for k in range(n_clusters):
+            mask = cluster_ids == k
+            if mask.sum() > 0:
+                dists = np.linalg.norm(X_class[mask] - centroids[k], axis=1)
+                intra_spreads[k] = dists.mean()
+
+        # Inter-cluster: distance to nearest OTHER centroid
+        centroid_dists = cdist(centroids, centroids, metric='euclidean')
+        np.fill_diagonal(centroid_dists, np.inf)
+        nearest_dists = centroid_dists.min(axis=1)
+        nearest_idx = centroid_dists.argmin(axis=1)
+
+        # Confidence = separation / internal spread
+        # (add epsilon to avoid division by zero)
+        with np.errstate(divide='ignore', invalid='ignore'):
+            confidences = nearest_dists / (intra_spreads + 1e-8)
+        confidences = np.nan_to_num(confidences, nan=0.0)
+
+        # Build result dict
+        class_prefix = {
+            "Bengin cases": "benign",
+            "Malignant cases": "malig",
+            "Normal cases": "normal",
+        }[class_name]
+
+        for k in range(n_clusters):
+            pid = f"{class_prefix}_{k}"
+            cluster_scores[pid] = {
+                "confidence": float(confidences[k]),
+                "intra_spread": float(intra_spreads[k]),
+                "nearest_neighbor_dist": float(nearest_dists[k]),
+                "nearest_neighbor_id": f"{class_prefix}_{int(nearest_idx[k])}",
+                "class_name": class_name,
+                "cluster_id": k,
+            }
+
+    return cluster_scores, X_pca, pca

classes/features.py

@@ -0,0 +1,212 @@
+"""
+Feature extraction using pretrained CNNs (PyTorch).
+
+Saves feature maps as .npz files with the model name in the filename,
+e.g. "VGG16_features.npz" containing X, Y, and filenames arrays.
+
+Supports all five models from the manuscript:
+    VGG16, DenseNet121, EfficientNetB1, MobileNetV2, ResNet50
+"""
+
+import os
+import numpy as np
+from PIL import Image
+import torch
+import torch.nn as nn
+from torchvision import transforms
+from torchvision.models import (
+    vgg16, VGG16_Weights,
+    densenet121, DenseNet121_Weights,
+    efficientnet_b1, EfficientNet_B1_Weights,
+    mobilenet_v2, MobileNet_V2_Weights,
+    resnet50, ResNet50_Weights,
+)
+
+
+# ---------------------------------------------------------------------------
+# Model registry
+# ---------------------------------------------------------------------------
+
+MODEL_CONFIGS = {
+    "VGG16": {
+        "fn": vgg16,
+        "weights": VGG16_Weights.IMAGENET1K_V1,
+        "input_size": 224,
+    },
+    "DenseNet121": {
+        "fn": densenet121,
+        "weights": DenseNet121_Weights.IMAGENET1K_V1,
+        "input_size": 224,
+    },
+    "EfficientNetB1": {
+        "fn": efficientnet_b1,
+        "weights": EfficientNet_B1_Weights.IMAGENET1K_V1,
+        "input_size": 240,
+    },
+    "MobileNetV2": {
+        "fn": mobilenet_v2,
+        "weights": MobileNet_V2_Weights.IMAGENET1K_V1,
+        "input_size": 224,
+    },
+    "ResNet50": {
+        "fn": resnet50,
+        "weights": ResNet50_Weights.IMAGENET1K_V1,
+        "input_size": 224,
+    },
+}
+
+
+class FeatureExtractor:
+    """Extract deep features from images using a pretrained CNN.
+
+    Parameters
+    ----------
+    model_name : str
+        One of: VGG16, DenseNet121, EfficientNetB1, MobileNetV2, ResNet50.
+    device : str or None
+        Torch device string. Auto-detected if None.
+    """
+
+    def __init__(self, model_name="VGG16", device=None):
+        if model_name not in MODEL_CONFIGS:
+            raise ValueError(
+                f"Unsupported model '{model_name}'. "
+                f"Choose from: {list(MODEL_CONFIGS.keys())}"
+            )
+        self.model_name = model_name
+        self.cfg = MODEL_CONFIGS[model_name]
+        self.input_size = self.cfg["input_size"]
+        self.device = device or ("cuda" if torch.cuda.is_available() else "cpu")
+        self.model = None
+        self._load_model()
+
+    def _load_model(self):
+        """Load pretrained model and strip the classifier head."""
+        full_model = self.cfg["fn"](weights=self.cfg["weights"])
+
+        name = self.model_name
+
+        if name == "VGG16":
+            # Drop the classifier Sequential → output (B, 512, 7, 7)
+            self.model = nn.Sequential(*list(full_model.children())[:-1])
+
+        elif name == "DenseNet121":
+            # Keep conv stack (dense blocks), drop classifier Linear
+            # → output (B, 1024, 7, 7)
+            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)
+            self.model = full_model.features
+
+        elif name == "MobileNetV2":
+            # Keep conv stack, drop adaptive pool + classifier
+            # → output (B, 1280, 7, 7)
+            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
+
+        self.model.to(self.device)
+        self.model.eval()
+
+        # Preprocessing
+        self.transform = transforms.Compose([
+            transforms.Resize((self.input_size, self.input_size)),
+            transforms.ToTensor(),
+            transforms.Normalize(
+                mean=[0.485, 0.456, 0.406],
+                std=[0.229, 0.224, 0.225]
+            ),
+        ])
+
+    # ------------------------------------------------------------------
+    # extract / save / load — unchanged public API
+    # ------------------------------------------------------------------
+
+    def extract(self, image_dir, valid_extensions=None, batch_size=32):
+        """Extract features from all images in a directory tree.
+
+        Expects subdirectories named by class (e.g. "Bengin cases/").
+
+        Returns
+        -------
+        X : np.ndarray  (n_images, n_features)
+        Y : np.ndarray  (n_images,)  class labels
+        filenames : np.ndarray  (n_images,)  image filenames
+        """
+        if valid_extensions is None:
+            valid_extensions = (".png", ".jpg", ".jpeg",
+                                ".tif", ".tiff", ".bmp")
+
+        image_paths, labels, fnames = [], [], []
+        for class_name in sorted(os.listdir(image_dir)):
+            class_path = os.path.join(image_dir, class_name)
+            if not os.path.isdir(class_path):
+                continue
+            print(f"  Scanning class: {class_name}")
+            for file in sorted(os.listdir(class_path)):
+                if not file.lower().endswith(valid_extensions):
+                    continue
+                image_paths.append(os.path.join(class_path, file))
+                labels.append(class_name)
+                fnames.append(file)
+
+        n_images = len(image_paths)
+        print(f"  Found {n_images} images across all classes")
+
+        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:
+                try:
+                    img = Image.open(path).convert("RGB")
+                    tensor = self.transform(img)
+                    batch_tensors.append(tensor)
+                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())
+
+            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)
+        filenames = np.array(fnames)
+
+        print(f"  Feature matrix shape: {X.shape}")
+        return X, Y, filenames
+
+    def save_features(self, X, Y, filenames, output_dir):
+        """Save extracted features to a compressed .npz file."""
+        os.makedirs(output_dir, exist_ok=True)
+        filepath = os.path.join(output_dir, f"{self.model_name}_features.npz")
+        np.savez_compressed(filepath, X=X, Y=Y, filenames=filenames)
+        print(f"  Saved features to {filepath}")
+        return filepath
+
+    @staticmethod
+    def load_features(filepath):
+        """Load saved features from a .npz file.
+
+        Returns
+        -------
+        X, Y, filenames : np.ndarray
+        """
+        data = np.load(filepath, allow_pickle=True)
+        return data["X"], data["Y"], data["filenames"]

classes/patient_identifier.py

@@ -0,0 +1,190 @@
+"""
+Patient identification via K-means clustering.
+
+Two modes:
+  - Thumbnail-based: cluster on 64×64 grayscale images (Andreas's approach).
+  - Feature-based:   cluster on PCA-reduced CNN features.
+
+Since the IQ-OTH/NCCD dataset does not provide patient IDs,
+we estimate which images belong to the same patient by clustering
+within each class.
+"""
+
+import os
+import numpy as np
+from PIL import Image
+from sklearn.decomposition import PCA
+from sklearn.cluster import KMeans
+
+
+class PatientIdentifier:
+    """Estimate patient groups via K-means clustering within each class.
+
+    Parameters
+    ----------
+    patient_estimates : dict
+        Mapping from class_name -> number_of_patients, e.g.
+        {"Bengin cases": 15, "Malignant cases": 40, "Normal cases": 55}.
+    cluster_image_size : tuple of int
+        Grayscale thumbnail size for clustering (default (64, 64)).
+    random_state : int
+        Seed for PCA and K-means reproducibility.
+    """
+
+    def __init__(self, patient_estimates, cluster_image_size=(64, 64),
+                 random_state=42):
+        self.patient_estimates = patient_estimates
+        self.cluster_image_size = cluster_image_size
+        self.random_state = random_state
+        self.groups_ = None
+        self.labels_ = None
+        self.assignments_ = None  # dict: class_name -> list of cluster_ids per sample
+
+    # ------------------------------------------------------------------
+    # Public API
+    # ------------------------------------------------------------------
+
+    def identify_from_thumbnails(self, base_path, filenames, labels,
+                                 valid_extensions=None):
+        """Cluster images within each class using grayscale thumbnails.
+
+        (Original approach from Andreas — model-agnostic.)
+        """
+        self.labels_ = labels
+        groups = self._run_clustering(
+            labels, filenames,
+            data_loader_fn=lambda class_name, files: self._load_thumbnails(
+                base_path, class_name, files),
+            valid_extensions=valid_extensions,
+        )
+        self.groups_ = groups
+        return groups
+
+    def identify_from_features(self, X, labels):
+        """Cluster images within each class using pre-extracted CNN features.
+
+        Parameters
+        ----------
+        X : np.ndarray  (n_images, n_features)
+            Feature matrix from a CNN (e.g. VGG16 output).
+        labels : np.ndarray  (n_images,)
+            Class label strings in the same order as X.
+
+        Returns
+        -------
+        groups : np.ndarray  (n_images,)  group label strings
+        """
+        self.labels_ = labels
+        groups = self._run_clustering(
+            labels, None,
+            data_loader_fn=lambda class_name, idx_class: X[idx_class],
+        )
+        self.groups_ = groups
+        return groups
+
+    def build_assignment_dict(self, groups, filenames, labels):
+        """Convert flat group array into a human-readable dict.
+
+        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",
+        }
+        assignment = {}
+        for class_name, prefix in class_map.items():
+            idx_class = np.where(labels == class_name)[0]
+            # Get unique cluster IDs within this class (preserve order)
+            seen = []
+            cluster_ids = []
+            for g in groups[idx_class]:
+                cid = int(g.split("_cluster_")[-1])
+                if cid not in seen:
+                    seen.append(cid)
+                cluster_ids.append(seen.index(cid))
+
+            for local_id in seen:
+                patient_id = f"{prefix}_{local_id}"
+                mask = np.array(cluster_ids) == local_id
+                assignment[patient_id] = filenames[idx_class][mask].tolist()
+
+        self.assignments_ = assignment
+        return assignment
+
+    def get_unique_groups(self):
+        if self.groups_ is None:
+            raise RuntimeError("Call identify_*() first.")
+        return np.unique(self.groups_)
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    def _run_clustering(self, labels, filenames, data_loader_fn,
+                        valid_extensions=None):
+        """Shared K-means pipeline used by both thumbnail and feature modes."""
+        if valid_extensions is None:
+            valid_extensions = (".png", ".jpg", ".jpeg",
+                                ".tif", ".tiff", ".bmp")
+
+        n_samples = len(labels)
+        groups = np.array([None] * n_samples, dtype=object)
+
+        for class_name, n_clusters in self.patient_estimates.items():
+            print(f"\n  Clustering class: {class_name} "
+                  f"(k={n_clusters})")
+
+            idx_class = np.where(labels == class_name)[0]
+
+            X_cluster = data_loader_fn(class_name, idx_class)
+            if isinstance(X_cluster, list):
+                X_cluster = np.array(X_cluster, dtype=np.float32)
+
+            # Dimensionality reduction
+            n_pca = min(50, X_cluster.shape[0] - 1, X_cluster.shape[1])
+            pca = PCA(n_components=n_pca, random_state=self.random_state)
+            X_pca = pca.fit_transform(X_cluster)
+
+            # K-means
+            kmeans = KMeans(
+                n_clusters=n_clusters,
+                random_state=self.random_state,
+                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}"
+
+            cluster_sizes = np.bincount(clusters)
+            print(f"    Cluster sizes: min={cluster_sizes.min()}, "
+                  f"max={cluster_sizes.max()}, "
+                  f"mean={cluster_sizes.mean():.1f}")
+
+        n_unique = len(np.unique(groups))
+        print(f"\n  Estimated {n_unique} patient groups total")
+        return groups
+
+    def _load_thumbnails(self, base_path, class_name, files):
+        """Load grayscale 64×64 thumbnails for a class."""
+        class_path = os.path.join(base_path, class_name)
+        X = []
+        for file in files:
+            img_path = os.path.join(class_path, file)
+            try:
+                img = Image.open(img_path).convert("L")
+                img = img.resize(self.cluster_image_size)
+                img = np.array(img, dtype=np.float32) / 255.0
+                X.append(img.flatten())
+            except Exception as e:
+                print(f"    Error loading {file}: {e}")
+                X.append(np.zeros(
+                    self.cluster_image_size[0] *
+                    self.cluster_image_size[1],
+                    dtype=np.float32
+                ))
+        return np.array(X, dtype=np.float32)

classes/visualizations.py

@@ -0,0 +1,148 @@
+"""
+Visualization utilities for feature space exploration.
+
+PCA and t-SNE plots to inspect class separability in the
+extracted feature representations.
+"""
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from sklearn.decomposition import PCA
+from sklearn.manifold import TSNE
+
+
+def plot_pca_tsne(X, Y, legend_names=None, perplexity=35,
+                  random_state=42, output_path=None, figsize=(12, 5),
+                  dpi=150):
+    """Generate side-by-side PCA and t-SNE plots of feature vectors.
+
+    Parameters
+    ----------
+    X : np.ndarray of shape (n_samples, n_features)
+        Feature matrix (e.g. VGG16 outputs).
+    Y : np.ndarray of shape (n_samples,)
+        Class labels.
+    legend_names : dict, optional
+        Mapping from raw label -> display label, e.g.
+        {"Bengin cases": "Benign", ...}.
+    perplexity : int
+        t-SNE perplexity (default 35).
+    random_state : int
+        Seed for reproducibility.
+    output_path : str, optional
+        Path to save the figure. If None, saved as "pca_tsne.png".
+    figsize : tuple
+        Figure dimensions.
+    dpi : int
+        Output resolution.
+
+    Returns
+    -------
+    fig : matplotlib Figure
+    """
+    if legend_names is None:
+        legend_names = {}
+
+    n_samples, n_features = X.shape
+
+    # --- PCA: 50D → 2D ---
+    n_pca = min(50, n_samples - 1, n_features)
+    pca_pre = PCA(n_components=n_pca, random_state=random_state)
+    X_pca50 = pca_pre.fit_transform(X)
+
+    pca_2d = PCA(n_components=2, random_state=random_state)
+    X_pca2 = pca_2d.fit_transform(X_pca50)
+    explained_variance = pca_2d.explained_variance_ratio_ * 100
+
+    # --- t-SNE: 50D → 2D ---
+    tsne = TSNE(
+        n_components=2,
+        perplexity=perplexity,
+        learning_rate="auto",
+        init="pca",
+        random_state=random_state
+    )
+    X_tsne = tsne.fit_transform(X_pca50)
+
+    # --- Plot ---
+    fig, axes = plt.subplots(1, 2, figsize=figsize)
+    unique_classes = np.unique(Y)
+
+    for class_name in unique_classes:
+        idx = Y == class_name
+        label = legend_names.get(class_name, class_name)
+
+        axes[0].scatter(
+            X_pca2[idx, 0], X_pca2[idx, 1],
+            label=label, alpha=0.7, s=8
+        )
+    axes[0].set_xlabel(f"PC1 ({explained_variance[0]:.1f}%)")
+    axes[0].set_ylabel(f"PC2 ({explained_variance[1]:.1f}%)")
+    axes[0].set_title("PCA")
+    axes[0].legend(markerscale=2)
+
+    for class_name in unique_classes:
+        idx = Y == class_name
+        label = legend_names.get(class_name, class_name)
+
+        axes[1].scatter(
+            X_tsne[idx, 0], X_tsne[idx, 1],
+            label=label, alpha=0.7, s=8
+        )
+    axes[1].set_xlabel("t-SNE dimension 1")
+    axes[1].set_ylabel("t-SNE dimension 2")
+    axes[1].set_title("t-SNE")
+    axes[1].legend(markerscale=2)
+
+    plt.tight_layout()
+
+    if output_path is None:
+        output_path = "pca_tsne.png"
+    plt.savefig(output_path, dpi=dpi)
+    plt.close()
+    print(f"  Saved PCA/t-SNE plot to {output_path}")
+
+    return fig
+
+
+def plot_cv_accuracy(cv_results, split_type="", output_path=None,
+                     figsize=(7, 5), dpi=150):
+    """Plot cross-validation accuracy vs number of selected features.
+
+    Parameters
+    ----------
+    cv_results : list of dict
+        Each dict has keys "nfeatures" and "best_cv_accuracy".
+    split_type : str
+        Label for the plot title, e.g. "Image-level split".
+    output_path : str, optional
+        Output file path.
+    figsize : tuple
+    dpi : int
+
+    Returns
+    -------
+    fig : matplotlib Figure
+    """
+    fig, ax = plt.subplots(figsize=figsize)
+
+    nfeats = [r["nfeatures"] for r in cv_results]
+    accs = [r["best_cv_accuracy"] for r in cv_results]
+
+    ax.plot(nfeats, accs, marker="o")
+    ax.set_xscale("log")
+    ax.set_xlabel("Number of selected features")
+    ax.set_ylabel("Mean 5-fold validation accuracy")
+    if split_type:
+        ax.set_title(split_type)
+    plt.tight_layout()
+
+    if output_path is None:
+        output_path = f"{split_type.replace(' ', '_').replace('-', '_')}_cv.png"
+    plt.savefig(output_path, dpi=dpi)
+    plt.close()
+    print(f"  Saved CV accuracy plot to {output_path}")
+
+    return fig

results/classification_results.json

@@ -0,0 +1,52 @@
+[
+  {
+    "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/simple_patient_manifest.csv

@@ -0,0 +1,111 @@
+patient_id,class,n_images,images
+Benign_00,Benign,13,B_021;B_022;B_023;B_024;B_025;B_035;B_066;B_067;B_074;B_075;B_086;B_087;B_088
+Benign_01,Benign,16,B_014;B_015;B_016;B_017;B_018;B_019;B_020;B_036;B_037;B_038;B_039;B_040;B_041;B_089;B_090;B_091
+Benign_02,Benign,9,B_070;B_071;B_072;B_073;B_081;B_082;B_083;B_084;B_085
+Benign_03,Benign,10,B_101;B_102;B_103;B_104;B_105;B_106;B_107;B_108;B_109;B_110
+Benign_04,Benign,6,B_111;B_112;B_113;B_114;B_115;B_116
+Benign_05,Benign,4,B_026;B_027;B_028;B_029
+Benign_06,Benign,12,B_030;B_031;B_042;B_043;B_044;B_045;B_060;B_061;B_062;B_063;B_064;B_065
+Benign_07,Benign,3,B_001;B_002;B_003
+Benign_08,Benign,7,B_032;B_033;B_034;B_117;B_118;B_119;B_120
+Benign_09,Benign,7,B_068;B_069;B_076;B_077;B_078;B_079;B_080
+Benign_10,Benign,7,B_006;B_007;B_008;B_056;B_057;B_058;B_059
+Benign_11,Benign,10,B_009;B_010;B_011;B_012;B_013;B_053;B_054;B_055;B_092;B_093
+Benign_12,Benign,5,B_047;B_048;B_049;B_050;B_051
+Benign_13,Benign,4,B_004;B_005;B_046;B_052
+Benign_14,Benign,7,B_094;B_095;B_096;B_097;B_098;B_099;B_100
+Malignant_00,Malignant,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_01,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_02,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_03,Malignant,26,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
+Malignant_04,Malignant,8,M_192;M_193;M_194;M_195;M_196;M_197;M_198;M_199
+Malignant_05,Malignant,12,M_029;M_326;M_327;M_328;M_329;M_330;M_331;M_332;M_341;M_342;M_343;M_344
+Malignant_06,Malignant,9,M_036;M_356;M_357;M_358;M_360;M_361;M_362;M_363;M_364
+Malignant_07,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_08,Malignant,6,M_228;M_229;M_230;M_231;M_232;M_233
+Malignant_09,Malignant,15,M_019;M_020;M_021;M_022;M_023;M_024;M_025;M_026;M_027;M_028;M_086;M_087;M_088;M_538;M_542
+Malignant_10,Malignant,18,M_030;M_037;M_038;M_039;M_040;M_041;M_042;M_043;M_387;M_388;M_389;M_390;M_391;M_392;M_393;M_512;M_513;M_514
+Malignant_11,Malignant,25,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
+Malignant_12,Malignant,8,M_079;M_080;M_081;M_082;M_323;M_324;M_325;M_507
+Malignant_13,Malignant,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_14,Malignant,19,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_15,Malignant,21,M_394;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_16,Malignant,32,M_008;M_009;M_015;M_053;M_054;M_055;M_056;M_057;M_058;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_365;M_366
+Malignant_17,Malignant,5,M_016;M_523;M_524;M_525;M_526
+Malignant_18,Malignant,13,M_446;M_447;M_448;M_449;M_450;M_451;M_452;M_453;M_454;M_455;M_456;M_457;M_458
+Malignant_19,Malignant,6,M_440;M_441;M_442;M_443;M_444;M_445
+Malignant_20,Malignant,26,M_280;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_21,Malignant,13,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_137
+Malignant_22,Malignant,8,M_240;M_241;M_242;M_243;M_244;M_245;M_246;M_247
+Malignant_23,Malignant,12,M_180;M_181;M_182;M_183;M_184;M_185;M_186;M_187;M_188;M_189;M_190;M_191
+Malignant_24,Malignant,21,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_109
+Malignant_25,Malignant,15,M_122;M_123;M_124;M_125;M_126;M_127;M_128;M_129;M_130;M_131;M_132;M_133;M_134;M_135;M_136
+Malignant_26,Malignant,12,M_508;M_509;M_510;M_511;M_515;M_516;M_517;M_518;M_519;M_520;M_521;M_522
+Malignant_27,Malignant,9,M_044;M_045;M_046;M_047;M_048;M_049;M_050;M_051;M_052
+Malignant_28,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_29,Malignant,8,M_459;M_460;M_461;M_462;M_463;M_464;M_465;M_466
+Malignant_30,Malignant,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_31,Malignant,9,M_200;M_201;M_202;M_203;M_204;M_205;M_206;M_207;M_208
+Malignant_32,Malignant,17,M_031;M_032;M_033;M_034;M_035;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_33,Malignant,7,M_381;M_382;M_383;M_384;M_385;M_386;M_527
+Malignant_34,Malignant,10,M_528;M_529;M_530;M_531;M_532;M_533;M_534;M_535;M_536;M_537
+Malignant_35,Malignant,16,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_36,Malignant,16,M_083;M_084;M_085;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
+Malignant_37,Malignant,6,M_234;M_235;M_236;M_237;M_238;M_239
+Malignant_38,Malignant,15,M_333;M_334;M_335;M_336;M_337;M_338;M_339;M_340;M_345;M_346;M_347;M_348;M_349;M_350;M_351
+Malignant_39,Malignant,19,M_010;M_011;M_012;M_013;M_014;M_138;M_139;M_140;M_141;M_142;M_143;M_144;M_145;M_146;M_147;M_148;M_149;M_150;M_281
+Normal_00,Normal,12,N_014;N_015;N_016;N_017;N_018;N_019;N_020;N_046;N_047;N_396;N_397;N_413
+Normal_01,Normal,6,N_132;N_133;N_134;N_400;N_401;N_414
+Normal_02,Normal,4,N_317;N_318;N_322;N_323
+Normal_03,Normal,19,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_371;N_372;N_373;N_406;N_407;N_408
+Normal_04,Normal,8,N_324;N_328;N_329;N_330;N_331;N_333;N_335;N_411
+Normal_05,Normal,5,N_169;N_170;N_171;N_172;N_173
+Normal_06,Normal,8,N_065;N_066;N_098;N_099;N_100;N_106;N_107;N_292
+Normal_07,Normal,9,N_345;N_346;N_347;N_348;N_350;N_351;N_352;N_353;N_354
+Normal_08,Normal,5,N_174;N_182;N_183;N_184;N_185
+Normal_09,Normal,8,N_297;N_298;N_307;N_308;N_309;N_310;N_311;N_312
+Normal_10,Normal,17,N_131;N_144;N_145;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_415;N_416
+Normal_11,Normal,12,N_280;N_281;N_282;N_283;N_284;N_285;N_286;N_287;N_288;N_294;N_295;N_296
+Normal_12,Normal,9,N_336;N_337;N_338;N_339;N_340;N_341;N_342;N_343;N_410
+Normal_13,Normal,4,N_033;N_355;N_356;N_357
+Normal_14,Normal,8,N_115;N_116;N_117;N_118;N_119;N_120;N_121;N_349
+Normal_15,Normal,6,N_203;N_204;N_205;N_207;N_208;N_209
+Normal_16,Normal,7,N_223;N_224;N_225;N_226;N_227;N_228;N_229
+Normal_17,Normal,5,N_158;N_164;N_165;N_166;N_167
+Normal_18,Normal,5,N_240;N_241;N_242;N_243;N_398
+Normal_19,Normal,15,N_255;N_256;N_257;N_258;N_259;N_260;N_261;N_262;N_263;N_264;N_265;N_266;N_267;N_268;N_269
+Normal_20,Normal,7,N_050;N_051;N_052;N_053;N_054;N_055;N_056
+Normal_21,Normal,9,N_039;N_040;N_041;N_042;N_043;N_044;N_045;N_206;N_254
+Normal_22,Normal,9,N_135;N_136;N_137;N_138;N_139;N_140;N_141;N_142;N_143
+Normal_23,Normal,12,N_073;N_074;N_075;N_076;N_077;N_078;N_079;N_080;N_190;N_191;N_192;N_193
+Normal_24,Normal,7,N_008;N_009;N_010;N_011;N_012;N_013;N_412
+Normal_25,Normal,10,N_386;N_387;N_388;N_389;N_390;N_391;N_392;N_393;N_394;N_395
+Normal_26,Normal,8,N_210;N_313;N_314;N_315;N_316;N_319;N_320;N_321
+Normal_27,Normal,12,N_064;N_067;N_068;N_069;N_070;N_071;N_072;N_195;N_196;N_197;N_198;N_199
+Normal_28,Normal,5,N_034;N_035;N_036;N_038;N_230
+Normal_29,Normal,6,N_004;N_005;N_127;N_128;N_129;N_130
+Normal_30,Normal,6,N_021;N_022;N_023;N_024;N_037;N_385
+Normal_31,Normal,9,N_231;N_232;N_233;N_234;N_235;N_236;N_237;N_238;N_239
+Normal_32,Normal,5,N_059;N_060;N_061;N_062;N_063
+Normal_33,Normal,4,N_027;N_028;N_029;N_030
+Normal_34,Normal,6,N_188;N_189;N_194;N_279;N_291;N_293
+Normal_35,Normal,9,N_270;N_271;N_272;N_273;N_274;N_275;N_276;N_277;N_278
+Normal_36,Normal,9,N_211;N_212;N_213;N_214;N_215;N_216;N_217;N_218;N_219
+Normal_37,Normal,3,N_220;N_221;N_222
+Normal_38,Normal,9,N_101;N_102;N_103;N_104;N_105;N_108;N_109;N_110;N_111
+Normal_39,Normal,10,N_244;N_245;N_246;N_247;N_248;N_249;N_250;N_251;N_252;N_253
+Normal_40,Normal,7,N_175;N_176;N_177;N_178;N_179;N_180;N_181
+Normal_41,Normal,5,N_168;N_402;N_403;N_404;N_405
+Normal_42,Normal,4,N_025;N_026;N_031;N_032
+Normal_43,Normal,6,N_299;N_300;N_301;N_302;N_303;N_304
+Normal_44,Normal,5,N_325;N_326;N_327;N_332;N_334
+Normal_45,Normal,12,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_46,Normal,7,N_001;N_002;N_003;N_006;N_007;N_057;N_058
+Normal_47,Normal,5,N_048;N_049;N_200;N_201;N_202
+Normal_48,Normal,7,N_081;N_082;N_083;N_084;N_085;N_086;N_087
+Normal_49,Normal,9,N_088;N_089;N_090;N_091;N_186;N_187;N_289;N_290;N_399
+Normal_50,Normal,2,N_305;N_306
+Normal_51,Normal,4,N_112;N_113;N_114;N_344
+Normal_52,Normal,5,N_159;N_160;N_161;N_162;N_163
+Normal_53,Normal,6,N_092;N_093;N_094;N_095;N_096;N_097
+Normal_54,Normal,5,N_122;N_123;N_124;N_125;N_126

scripts/classification.py

@@ -0,0 +1,50 @@
+#!/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/create_patient_groups.py

@@ -0,0 +1,105 @@
+#!/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

@@ -0,0 +1,97 @@
+#!/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/simple_patient_tsne.py

@@ -0,0 +1,192 @@
+#!/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__))))
+
+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
+    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
+
+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
+    centroids = {}
+    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:
+        mask = Y == class_name
+        for i in np.where(mask)[0]:
+            img = img_nums[i]
+            old_pid = image_to_patient[img]
+            # Find nearest centroid in same class
+            best_pid = old_pid
+            best_dist = float('inf')
+            for pid, c in centroids.items():
+                if pid.startswith(CLASS_NAMES[class_name]):
+                    d = float(np.linalg.norm(X_pca[i] - c))
+                    if d < best_dist:
+                        best_dist = d
+                        best_pid = pid
+            if best_pid != old_pid:
+                patient_to_images[old_pid].remove(img)
+                patient_to_images[best_pid].append(img)
+                image_to_patient[img] = best_pid
+                moves += 1
+    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 = [], []
+        for img in patient_to_images[pid]:
+            i = np.where(img_nums == img)[0][0]
+            pts_x.append(X_tsne[i, 0])
+            pts_y.append(X_tsne[i, 1])
+        ax.scatter(pts_x, pts_y, c=[color], s=18, alpha=0.8, label='_nolegend_')
+        # Centroid diamond (no label)
+        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 = [
+        Line2D([0], [0], marker='o', color='w', markerfacecolor='gray',
+               markersize=8, label='Patient images (dots)'),
+        Line2D([0], [0], marker='D', color='w', markerfacecolor='gray',
+               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")
+    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")