began work on v3

2026-03-19 16:58:29 +01:00
parent 786457b30d
commit eb9eafe715
42 changed files with 11214 additions and 1 deletions
@@ -0,0 +1,240 @@
+from __future__ import annotations
+
+from typing import Callable, Dict, List, Optional
+
+import numpy as np
+import pandas as pd
+
+from v3.classes.data_bundle import DataBundle
+
+# ---- Pachymetry → IOP correction (per PAPILA Table 3) ----
+_PACHY_TABLE: Dict[int, int] = {
+    475: +5,
+    485: +4,
+    495: +4,
+    505: +3,
+    515: +2,
+    525: +1,
+    535: +1,
+    545: 0,
+    555: -1,
+    565: -1,
+    575: -2,
+    585: -3,
+    595: -4,
+    605: -4,
+    615: -5,
+}
+_PACHY_KEYS = np.array(sorted(_PACHY_TABLE.keys()))
+
+
+def _nearest_pachy_key(x: float) -> int:
+    idx = int(np.argmin(np.abs(_PACHY_KEYS - float(x))))
+    return int(_PACHY_KEYS[idx])
+
+
+def _fit_perkins_converter(
+    frames: List[pd.DataFrame], method: str
+) -> Callable[[float, Optional[float]], float]:
+    """
+    Fit a Perkins→Pneumatic converter from pooled paired observations across all frames.
+    Returns a callable: converter(perkins_value, pachymetry_value) -> float.
+    Supported methods: "ratio", "ols", "lad", "multi".
+    """
+    combined = pd.concat(frames, ignore_index=True)
+    paired = combined.dropna(subset=["Pneumatic", "Perkins"])
+    pneumatic = paired["Pneumatic"].values.astype(float)
+    perkins   = paired["Perkins"].values.astype(float)
+
+    if len(paired) == 0:
+        raise ValueError("No paired Pneumatic+Perkins observations found; cannot fit converter.")
+
+    if method == "ratio":
+        ratio = float((pneumatic / perkins).mean())
+        def converter_ratio(p: float, pachy: Optional[float] = None) -> float:
+            return p * ratio
+        return converter_ratio
+
+    elif method == "ols":
+        from scipy import stats as _stats
+        slope, intercept, *_ = _stats.linregress(perkins, pneumatic)
+        slope, intercept = float(slope), float(intercept)
+        def converter_ols(p: float, pachy: Optional[float] = None) -> float:
+            return p * slope + intercept
+        return converter_ols
+
+    elif method == "lad":
+        from scipy import stats as _stats
+        from scipy.optimize import minimize as _minimize
+        slope0, intercept0, *_ = _stats.linregress(perkins, pneumatic)
+        def _lad_loss(params):
+            a, b = params
+            return np.abs(pneumatic - (a * perkins + b)).mean()
+        res = _minimize(_lad_loss, x0=[slope0, intercept0], method="Nelder-Mead")
+        slope, intercept = float(res.x[0]), float(res.x[1])
+        def converter_lad(p: float, pachy: Optional[float] = None) -> float:
+            return p * slope + intercept
+        return converter_lad
+
+    elif method == "multi":
+        from numpy.linalg import lstsq as _lstsq
+        paired_multi = combined.dropna(subset=["Pneumatic", "Perkins", "Pachymetry"])
+        if len(paired_multi) == 0:
+            raise ValueError("No paired Pneumatic+Perkins+Pachymetry rows; cannot fit multi method.")
+        pneu  = paired_multi["Pneumatic"].values.astype(float)
+        perk  = paired_multi["Perkins"].values.astype(float)
+        pachy_vals = paired_multi["Pachymetry"].values.astype(float)
+        X = np.column_stack([perk, pachy_vals, np.ones(len(perk))])
+        coeffs, *_ = _lstsq(X, pneu, rcond=None)
+        slope, pachy_coef, intercept = float(coeffs[0]), float(coeffs[1]), float(coeffs[2])
+        pachy_fallback = float(pachy_vals.mean())
+        def converter_multi(p: float, pachy: Optional[float] = None) -> float:
+            pv = pachy if (pachy is not None and not np.isnan(pachy)) else pachy_fallback
+            return p * slope + pachy_coef * pv + intercept
+        return converter_multi
+
+    else:
+        raise ValueError(f"Unknown iop_corr_method: {method!r}. Choose ratio/ols/lad/multi.")
+
+
+def _pick_iop(row: pd.Series, converter: Callable) -> float:
+    """Prefer Pneumatic; convert Perkins to Pneumatic scale if Pneumatic is absent."""
+    pneumatic = row.get("Pneumatic", np.nan)
+    if not pd.isna(pneumatic):
+        return float(pneumatic)
+    perkins = row.get("Perkins", np.nan)
+    if pd.isna(perkins):
+        return np.nan
+    pachy = row.get("Pachymetry", np.nan)
+    return converter(float(perkins), None if pd.isna(pachy) else float(pachy))
+
+
+def _correct_iop(raw_iop: float, pachy: float) -> float:
+    """Return corrected IOP using nearest pachymetry bin; if pachy missing, return raw."""
+    if pd.isna(raw_iop):
+        return np.nan
+    if pd.isna(pachy):
+        return float(raw_iop)
+    key = _nearest_pachy_key(float(pachy))
+    return float(raw_iop) + float(_PACHY_TABLE[key])
+
+
+def _apply_iop_and_drop_md(
+    df: pd.DataFrame,
+    converter: Callable,
+    drop_raw: bool = False,
+) -> pd.DataFrame:
+    """Add IOP_raw/IOP_corr and drop source IOP columns + VF_MD if present (in-place safe)."""
+    df["IOP_raw"] = df.apply(lambda row: _pick_iop(row, converter), axis=1)
+    pachy = df.get("Pachymetry", pd.Series(np.nan, index=df.index))
+    df["IOP_corr"] = [
+        _correct_iop(r, p) for r, p in zip(df["IOP_raw"].values, pachy.values)
+    ]
+    drop_cols = [c for c in ("Pneumatic", "Perkins", "VF_MD") if c in df.columns]
+    if drop_raw:
+        drop_cols.append("IOP_raw")
+    if drop_cols:
+        df.drop(columns=drop_cols, inplace=True)
+    return df
+
+
+def _canonicalize_eye_column(df: pd.DataFrame) -> None:
+    if "eyeID" in df.columns:
+        src = "eyeID"
+    else:
+        src = None
+        for c in df.columns:
+            if "eye" in c.lower():
+                src = c
+                break
+        if src is None:
+            df["eyeID"] = "OS"
+            return
+
+    s = df[src]
+
+    def norm(v):
+        if pd.isna(v):
+            return None
+        x = str(v).strip().upper()
+        if x in {"OS", "L", "LEFT", "0"}:
+            return "OS"
+        if x in {"OD", "R", "RIGHT", "1"}:
+            return "OD"
+        try:
+            num = int(float(x))
+            return "OD" if num % 2 == 1 else "OS"
+        except Exception:
+            return None
+
+    mapped = s.map(norm)
+    uniq = {u for u in mapped.dropna().unique().tolist()}
+    if not uniq.issubset({"OS", "OD"}):
+        raise ValueError(f"eyeID must be binary; found values {sorted(uniq)}")
+    df["eyeID"] = mapped.fillna("OS")
+
+
+def build_papila_data(
+    *,
+    image_dir: str,
+    clinical_dir: str,
+    label_col: str,
+    cat_cols: List[str],
+    n_splits: int = 5,
+    random_seed: int = 42,
+    iop_corr_method: str = "ratio",
+    iop_drop_raw: bool = False,
+    exclude_cols: Optional[List[str]] = None,
+) -> DataBundle:
+    """
+    Build a DataBundle for PAPILA with dataset-specific preprocessing:
+      - load OD/OS Excel sheets
+      - normalize Patient ID
+      - canonicalize eyeID
+      - compute IOP_raw / IOP_corr, drop VF_MD
+      - build feature typing & folds
+    """
+    _exclude = list(exclude_cols) if exclude_cols else []
+
+    # Remove excluded cols from cat_cols too so the bundle doesn't try to encode them
+    effective_cat_cols = [c for c in cat_cols if c not in _exclude]
+
+    bundle = DataBundle(
+        image_dir=image_dir,
+        clinical_dir=clinical_dir,
+        label_col=label_col,
+        patient_col="Patient ID",
+        cat_cols=effective_cat_cols,
+        n_splits=n_splits,
+        random_seed=random_seed,
+        filename_template="RET{pid:03d}{eye}.jpg",
+    )
+
+    od = pd.read_excel(f"{clinical_dir}/patient_data_od.xlsx", header=1)
+    od["eyeID"] = "OD"
+    os = pd.read_excel(f"{clinical_dir}/patient_data_os.xlsx", header=1)
+    os["eyeID"] = "OS"
+
+    for frame in (od, os):
+        if "Patient ID" not in frame.columns and "ID" in frame.columns:
+            frame.rename(columns={"ID": "Patient ID"}, inplace=True)
+        frame["Patient ID"] = frame["Patient ID"].astype(str).str.extract(r"(\d+)")[0].astype(int)
+        _canonicalize_eye_column(frame)
+
+    bundle.add_df(od, id_column="ID", exclude_cols=_exclude or None)
+    bundle.add_df(os, id_column="ID", exclude_cols=_exclude or None)
+
+    converter = _fit_perkins_converter(bundle.frames, method=iop_corr_method)
+    for i in range(len(bundle.frames)):
+        bundle.frames[i] = _apply_iop_and_drop_md(
+            bundle.frames[i], converter=converter, drop_raw=iop_drop_raw
+        )
+
+    bundle._refresh_master_df(exclude_cols=_exclude or None)
+    bundle._infer_or_validate_feature_types(exclude_cols=_exclude or None)
+    bundle._compute_numeric_stats()
+    bundle._build_cat_maps()
+    bundle._compute_feature_dim()
+    bundle._build_kfold_indices()
+
+    return bundle