diff --git a/v4/classes/v4_hypertower.py b/v4/classes/v4_hypertower.py index 8eba28a..6a773d0 100644 --- a/v4/classes/v4_hypertower.py +++ b/v4/classes/v4_hypertower.py @@ -266,6 +266,7 @@ def main(): save_features = cfg.get("save_features", False) fold_results = [] eval_stage_preds = [] # list[dict] — one per fold, only for eval_stage + all_phase_preds: dict[str, list[dict]] = {} # phase → list[dict] across folds t0 = time.time() for fold in range(cfg.get("folds", 5)): @@ -277,6 +278,11 @@ def main(): fold_results.append(result) if save_predictions and eval_stage in fold_preds: eval_stage_preds.append(fold_preds[eval_stage]) + if save_features: + for ph, pdata in fold_preds.items(): + if pdata.get("val_z") is None: + continue + all_phase_preds.setdefault(ph, []).append(pdata) print( f" fold{fold+1} DONE" f" val_auc={result.get(f'{eval_stage}_val_auc', float('nan')):.4f}" @@ -346,42 +352,43 @@ def main(): store.save(pred_path) print(f"Predictions saved: {pred_path}", flush=True) - if save_features and eval_stage_preds: - emb_dim = eval_stage_preds[0]["val_z"].shape[-1] - fstore = FeatureStore(n_folds=len(eval_stage_preds)) - - # Build entity_id / y_true universe (same as predictions). - seen, all_ids, id_to_y = set(), [], {} - for fp in eval_stage_preds: - for eid, y in zip(fp["val_ids"], fp["val_y"]): - k = str(eid) - if k not in seen: - seen.add(k); all_ids.append(eid) - id_to_y[k] = int(y) - if fp.get("test_ids"): - for eid, y in zip(fp["test_ids"], fp["test_y"]): + if save_features and all_phase_preds: + # One FeatureStore covers all phases; each phase gets its own group. + n_folds_any = max(len(v) for v in all_phase_preds.values()) + fstore = FeatureStore(n_folds=n_folds_any) + for phase, phase_preds in all_phase_preds.items(): + emb_dim = phase_preds[0]["val_z"].shape[-1] + seen, all_ids, id_to_y = set(), [], {} + for fp in phase_preds: + for eid, y in zip(fp["val_ids"], fp["val_y"]): k = str(eid) if k not in seen: seen.add(k); all_ids.append(eid) id_to_y[k] = int(y) + if fp.get("test_ids"): + for eid, y in zip(fp["test_ids"], fp["test_y"]): + k = str(eid) + if k not in seen: + seen.add(k); all_ids.append(eid) + id_to_y[k] = int(y) + y_true = np.array([id_to_y.get(str(e), -1) for e in all_ids], dtype=np.int64) + fstore.register_phase(phase=phase, entity_ids=all_ids, y_true=y_true) + fstore.register_head(phase=phase, head=f"{phase}_embedding", + n_epochs=1, embedding_dim=emb_dim) - y_true = np.array([id_to_y.get(str(e), -1) for e in all_ids], dtype=np.int64) - fstore.register_phase(phase=eval_stage, entity_ids=all_ids, y_true=y_true) - fstore.register_head(phase=eval_stage, head=f"{eval_stage}_embedding", - n_epochs=1, embedding_dim=emb_dim) - - for fold_idx, fp in enumerate(eval_stage_preds): - fstore.record(eval_stage, fold_idx, 0, fp["val_ids"], - f"{eval_stage}_embedding", fp["val_z"]) - fstore.set_split(eval_stage, fold_idx, fp["val_ids"], "val") - if fp.get("test_ids") and fp.get("test_z") is not None: - fstore.record(eval_stage, fold_idx, 0, fp["test_ids"], - f"{eval_stage}_embedding", fp["test_z"]) - fstore.set_split(eval_stage, fold_idx, fp["test_ids"], "test") + for fold_idx, fp in enumerate(phase_preds): + fstore.record(phase, fold_idx, 0, fp["val_ids"], + f"{phase}_embedding", fp["val_z"]) + fstore.set_split(phase, fold_idx, fp["val_ids"], "val") + if fp.get("test_ids") and fp.get("test_z") is not None: + fstore.record(phase, fold_idx, 0, fp["test_ids"], + f"{phase}_embedding", fp["test_z"]) + fstore.set_split(phase, fold_idx, fp["test_ids"], "test") feat_path = out_dir / "features.h5" fstore.save(feat_path) - print(f"Features saved: {feat_path}", flush=True) + print(f"Features saved (phases: {sorted(all_phase_preds.keys())}): {feat_path}", + flush=True) if __name__ == "__main__": diff --git a/v4/configs/cd_geom_duo.json b/v4/configs/cd_geom_duo.json new file mode 100644 index 0000000..41f3ef9 --- /dev/null +++ b/v4/configs/cd_geom_duo.json @@ -0,0 +1,143 @@ +{ + "_notes": [ + "Two-tower ensemble: cd + geom seg-CNN (UNet seg source).", + "Mirrors ensemble_fused.json but swaps the img tower for the geom tower.", + "Tells us how close the seg-CNN-over-masks gets to the image tower's contribution", + "when paired with clinical features." + ], + "run_name": "v4/cd_geom_duo", + "num_classes": 2, + "label_filter": [0, 1], + "split_identity_level": 1, + "eval_stage": "hb", + "save_predictions": true, + "save_features": true, + "seed": 1234, + "folds": 5, + "fold_seed": 100, + "output_root": "v4/results", + "out_dir_tags": ["binary"], + + "data": { + "module": "v4.classes.profiles.v4papila", + "args": { + "image_dir": "Papila/FundusImages", + "clinical_dir": "Papila/ClinicalData", + "label_col": "Diagnosis", + "iop_corr_method": "ratio", + "iop_drop_raw": true, + "exclude_cols": ["Axial_Length"], + "in_memory_cache": false + } + }, + + "towers": [ + { + "name": "cd", + "module": "v4.classes.towers.clinical_tower", + "class": "ClinicalEncoder", + "data_source": "matrix", + "args": { + "hidden_dim": 128 + } + }, + { + "name": "geom", + "module": "v4.classes.towers.geometry_tower", + "class": "GeometrySegEncoder", + "data_source": "image", + "args": { + "backbone": "resnet18", + "channels": 3, + "target_size": 224, + "augment": true, + "freeze_ratio": 0.0, + "seg_source": "unet", + "weights_path": "models/v2/refuge/segmentation/per_image/best.pt", + "contour_dir": "Papila/ExpertsSegmentations/Contours", + "unet_size": 512, + "normalize": "per_image", + "threshold": 0.5, + "crop_to_disc": true, + "finetune_epochs": 10, + "finetune_lr": 1e-5, + "finetune_batch_size": 4 + } + } + ], + + "stages": [ + { + "name": "cd_warm", + "type": "warm", + "tower": "cd", + "head_name": "cd_aux", + "level": "eye", + "epochs": 40 + }, + { + "name": "cd_aux", + "type": "head", + "input": "cd", + "train_with": "nt", + "bcd": true + }, + { + "name": "geom_aux", + "type": "head", + "input": "geom", + "train_with": "nt", + "bcd": true + }, + { + "name": "nt", + "type": "fusion", + "module": "v4.classes.bridges.fusion_bridge", + "class": "FusionBridge", + "inputs": ["cd", "geom"], + "level": "eye", + "epochs": 36, + "train_towers": true, + "warmup": { + "tower_epochs": 3, + "fused_epochs": 3 + }, + "args": { + "fusion_dim": 256 + } + }, + { + "name": "nt_head", + "type": "head", + "input": "nt", + "train_with": "nt" + }, + { + "name": "hb", + "type": "fusion", + "module": "v4.classes.bridges.hyperbridge", + "class": "HyperBridge", + "inputs": { "a": "nt", "b": "nt" }, + "level": "patient", + "epochs": 10, + "args": { + "hidden_dim": 256, + "mode": "embedding_mlp" + } + }, + { + "name": "hb_head", + "type": "head", + "input": "hb", + "train_with": "hb", + "args": { "dropout": 0.3 } + } + ], + + "training": { + "lr": 1e-4, + "batch_size": 8, + "bcd_prob": 0.5, + "tune_binary_threshold": true + } +} diff --git a/v4/distributed/jobs.db b/v4/distributed/jobs.db index 5f5bc10..f53c305 100644 Binary files a/v4/distributed/jobs.db and b/v4/distributed/jobs.db differ diff --git a/v4/scripts/analysis/inspect_embeddings.py b/v4/scripts/analysis/inspect_embeddings.py index 49dc494..99c469b 100644 --- a/v4/scripts/analysis/inspect_embeddings.py +++ b/v4/scripts/analysis/inspect_embeddings.py @@ -35,8 +35,20 @@ def load_phase(features_path: Path, phase: str | None): f"Phase {phase!r} not in {features_path} (available: {phases})" ) g = f[phase] - z = g["z"][:] # (n_folds, n_samples, n_dim) - split = g["split"][:].astype(str) # (n_folds, n_samples) + emb_keys = [k for k in g.keys() if k.endswith("_embedding")] + if not emb_keys: + raise SystemExit( + f"No '*_embedding' dataset in phase {phase!r} of {features_path}" + ) + # Prefer the embedding named after the phase if present, else first match. + emb_key = f"{phase}_embedding" if f"{phase}_embedding" in emb_keys else emb_keys[0] + z_raw = g[emb_key][:] + # Stored shape is (n_folds, n_epochs, n_samples, n_dim). Use last epoch. + if z_raw.ndim == 4: + z = z_raw[:, -1, :, :] + else: + z = z_raw + split = g["split"][:].astype(str) y_true = g["y_true"][:] return phase, z, split, y_true, phases @@ -58,9 +70,11 @@ def main(): phase, z, split, _, all_phases = load_phase(feat, args.phase) val_mask = (split == "val") - z_val = z[val_mask] # (n_val_total, n_dim) + z_val = z[val_mask] z_train = z[(split == "train")] - n_dim = z_val.shape[-1] + n_dim = z_val.shape[-1] if z_val.size else 0 + if n_dim == 0: + raise SystemExit(f"No val embeddings found for phase {phase!r}") print(f"Run: {args.run_dir}") print(f"File: {feat}") diff --git a/v4/scripts/analysis/logreg_cdr_compare.py b/v4/scripts/analysis/logreg_cdr_compare.py new file mode 100644 index 0000000..695f0a0 --- /dev/null +++ b/v4/scripts/analysis/logreg_cdr_compare.py @@ -0,0 +1,129 @@ +"""logreg_cdr_compare — compare LogReg AUCs across GT / base-UNet / fine-tuned-UNet CDR sources. + +Lets us answer: how much of the apparent +0.11 AUC from the LogReg(clinical+CDR-GT) +result is from the GT source quality vs the CDR features themselves? Tritower +uses per-fold-fine-tuned UNet — this script reproduces that pipeline as a 5-fold +LogReg baseline so the numbers are directly comparable to v4 cd_solo_geom and +tritower runs. + +5-fold StratifiedGroupKFold (patient-grouped), LogReg with StandardScaler. + +Run: + python -m v4.scripts.analysis.logreg_cdr_compare +""" +from __future__ import annotations + +import sys, time +from pathlib import Path + +import numpy as np +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import make_pipeline +from sklearn.model_selection import StratifiedGroupKFold +from sklearn.metrics import roc_auc_score + +REPO_ROOT = Path(__file__).resolve().parents[3] +sys.path.insert(0, str(REPO_ROOT)) + +from v4.classes.profiles.v4papila import build_data +from v4.classes.profiles.fundus_images import ( + build_geometry_loader, + _PapilaUNetMaskPipeline, + compute_geometry_features, +) + +CONTOUR_DIR = REPO_ROOT / "Papila/ExpertsSegmentations/Contours" +UNET_WEIGHTS = REPO_ROOT / "models/v2/refuge/segmentation/per_image/best.pt" + + +def main(): + data = build_data({ + "image_dir": "Papila/FundusImages", + "clinical_dir": "Papila/ClinicalData", + "iop_corr_method": "ratio", + "iop_drop_raw": True, + "exclude_cols": ["Axial_Length"], + }) + + df = data.df[data.df[data.label_col].isin([0, 1])].reset_index(drop=True) + samples = data.collect_samples(df) + y = df[data.label_col].astype(int).values + groups = df["Patient ID"].astype(int).values + view = data.matrix + X_clin = np.array([view.vectorize_entity(int(r["Patient ID"]), str(r["eyeID"])) + for _, r in df.iterrows()]) + print(f"n eyes={len(df)} n patients={len(set(groups))} " + f"label balance={np.bincount(y).tolist()}") + + # ─── GT CDR ─────────────────────────────────────────────────────────────── + print("\n[1/3] GT CDR (rasterise expert contours; split-independent)") + gt = build_geometry_loader("gt", contour_dir=str(CONTOUR_DIR)) + gt.precompute(samples) + gt_vecs = gt.all_vectors() + X_cdr_gt = np.array([gt_vecs[(int(r["Patient ID"]), str(r["eyeID"]))] + for _, r in df.iterrows()]) + + # ─── Base UNet CDR (no fine-tune; same masks every fold) ────────────────── + print("\n[2/3] base UNet CDR (REFUGE weights, no fine-tune)") + t0 = time.time() + unet_base = build_geometry_loader("unet", + weights_path=str(UNET_WEIGHTS), + contour_dir=str(CONTOUR_DIR), + finetune_epochs=0, + ) + unet_base.precompute(samples) + base_vecs = unet_base.all_vectors() + X_cdr_unet_base = np.array([base_vecs[(int(r["Patient ID"]), str(r["eyeID"]))] + for _, r in df.iterrows()]) + print(f" base inference done in {time.time()-t0:.1f}s") + + # ─── Build the 5 folds (used for both per-fold UNet ft AND for LogReg CV) ─ + sgkf = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=42) + splits = list(sgkf.split(np.arange(len(df)), y, groups)) + + # ─── Per-fold fine-tuned UNet CDR ───────────────────────────────────────── + print("\n[3/3] per-fold-fine-tuned UNet CDR (matches tritower / cd_solo_geom)") + X_cdr_unet_ft = np.zeros((len(df), 5), dtype=np.float32) + pipe = _PapilaUNetMaskPipeline( + str(UNET_WEIGHTS), + contour_dir=str(CONTOUR_DIR), + finetune_epochs=10, + finetune_lr=1e-5, + finetune_batch_size=4, + ) + for fold, (tr, te) in enumerate(splits): + t1 = time.time() + train_samples = [samples[i] for i in tr] + test_samples = [samples[i] for i in te] + pipe.reset_weights() + pipe.finetune(train_samples) + masks = pipe.predict(test_samples) + for j, (pid, eye, _) in zip(te, test_samples): + disc, cup = masks[(pid, eye)] + X_cdr_unet_ft[j] = compute_geometry_features(disc, cup) + print(f" fold {fold+1}/5 done ({time.time()-t1:.0f}s)") + + # ─── LogReg CV ──────────────────────────────────────────────────────────── + def cv(X, label): + aucs = [] + for tr, te in splits: + clf = make_pipeline(StandardScaler(), + LogisticRegression(max_iter=2000, C=1.0)) + clf.fit(X[tr], y[tr]) + aucs.append(roc_auc_score(y[te], clf.predict_proba(X[te])[:, 1])) + print(f" {label:32s}: {np.mean(aucs):.4f} ± {np.std(aucs):.4f}") + + print("\n=== StratifiedGroupKFold LogReg AUCs ===") + cv(X_clin, "clinical only") + cv(X_cdr_gt, "CDR-GT only") + cv(X_cdr_unet_base, "CDR-UNet (base) only") + cv(X_cdr_unet_ft, "CDR-UNet (per-fold ft) only") + print() + cv(np.concatenate([X_clin, X_cdr_gt], axis=1), "clinical + CDR-GT") + cv(np.concatenate([X_clin, X_cdr_unet_base], axis=1), "clinical + CDR-UNet (base)") + cv(np.concatenate([X_clin, X_cdr_unet_ft], axis=1), "clinical + CDR-UNet (ft)") + + +if __name__ == "__main__": + main() diff --git a/v4/scripts/experiments/tri_v1/baseline_tri_features_all.json b/v4/scripts/experiments/tri_v1/baseline_tri_features_all.json new file mode 100644 index 0000000..048449d --- /dev/null +++ b/v4/scripts/experiments/tri_v1/baseline_tri_features_all.json @@ -0,0 +1,10 @@ +[ + { + "_note": "Single rep of tritower default with save_features=true and ALL fusion phases captured (nt + hb). New run_name so it doesn't skip the prior single-phase result.", + "run_name": "experiments/tri_v1/baseline_tri_features_all", + "reps": 1, + "overrides": { + "save_features": true + } + } +] diff --git a/v4/scripts/experiments/tri_v1/cd_geom_duo.json b/v4/scripts/experiments/tri_v1/cd_geom_duo.json new file mode 100644 index 0000000..4a548a7 --- /dev/null +++ b/v4/scripts/experiments/tri_v1/cd_geom_duo.json @@ -0,0 +1,7 @@ +[ + { + "_note": "cd + geom (seg-CNN) two-tower ensemble — img-tower removed. 3 reps as a pilot; promote to 10 if it lands near img+cd's 0.896.", + "run_name": "experiments/tri_v1/cd_geom_duo", + "reps": 3 + } +]