{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceType":"datasetVersion","datasetSlug":"chest-xray-pneumonia"},{"sourceType":"competition","competitionSlug":"rsna-pneumonia-detection-challenge"}]},"accelerator":"GPU"},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"5720721e","cell_type":"markdown","source":"# Final Manuscript Experiment: ROOF-DARS with Full RSNA External Validation\n\nThis notebook generates the final manuscript-ready results for the proposed\n**Robust Out-of-Fold Difficulty-Aware Resampling (ROOF-DARS)** framework.\n\n## Experimental design\n\n- **Development:** complete Kermany pediatric chest X-ray training and validation folders.\n- **Internal test:** untouched original Kermany test folder.\n- **External validation:** complete eligible RSNA Pneumonia Detection Challenge cohort.\n- **RSNA positive endpoint:** `Lung Opacity`.\n- **RSNA negative endpoint:** `Normal`.\n- **Excluded RSNA endpoint:** `No Lung Opacity / Not Normal`.\n- **External prevalence:** preserved naturally; no class balancing is applied to the RSNA test cohort.\n- **Repeated training:** five prespecified seeds using identical patient/group-disjoint splits for every method.\n- **Calibration and thresholding:** learned only from Kermany development data.\n- **Primary threshold:** maximum MCC subject to sensitivity ≥ 0.90 on the Kermany calibration subset.\n- **Checkpointing:** every completed seed–method experiment is saved and can be resumed.\n\n## Compared methods\n\n1. Unbalanced BCE  \n2. Class-weighted BCE  \n3. Focal loss  \n4. Online hard-example mining  \n5. Random undersampling  \n6. SMOTE  \n7. Original in-sample DARS  \n8. OOF loss-only ablation  \n9. Proposed ROOF-DARS  \n\nThe notebook produces per-seed results, mean ± standard deviation, 95% confidence\nintervals, ensemble predictions, paired statistical tests, calibration analysis,\ngeneralization gaps, runtime measurements, and manuscript-ready tables and figures.\n","metadata":{}},{"id":"52dd5e69","cell_type":"markdown","source":"## 1. Kaggle runtime and required datasets\n\nUse a Kaggle GPU runtime and attach:\n\n- `paultimothymooney/chest-xray-pneumonia`\n- `rsna-pneumonia-detection-challenge`\n\nAccept the RSNA competition rules before attaching the competition data.\n","metadata":{}},{"id":"9f8c5245","cell_type":"code","source":"import importlib.util\nimport subprocess\nimport sys\n\nrequired = []\nif importlib.util.find_spec(\"torchxrayvision\") is None:\n    required.append(\"torchxrayvision==1.4.0\")\nif importlib.util.find_spec(\"pydicom\") is None:\n    required.append(\"pydicom==3.0.1\")\n\nif required:\n    subprocess.check_call(\n        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"--no-deps\", *required]\n    )\n\nif importlib.util.find_spec(\"imblearn\") is None:\n    subprocess.check_call(\n        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"imbalanced-learn\"]\n    )\n\nprint(\"Dependency check completed.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:02:48.314532Z","iopub.execute_input":"2026-07-03T10:02:48.315225Z","iopub.status.idle":"2026-07-03T10:02:51.583943Z","shell.execute_reply.started":"2026-07-03T10:02:48.315192Z","shell.execute_reply":"2026-07-03T10:02:51.583279Z"}},"outputs":[],"execution_count":null},{"id":"d6d74ba5","cell_type":"markdown","source":"## 2. Imports and environment\n","metadata":{}},{"id":"83d44e49","cell_type":"code","source":"import gc\nimport hashlib\nimport json\nimport math\nimport os\nimport platform\nimport random\nimport re\nimport time\nimport warnings\nfrom dataclasses import asdict, dataclass\nfrom pathlib import Path\nfrom typing import Any, Dict, Iterable, Optional, Sequence, Tuple\n\nos.environ.setdefault(\"TF_CPP_MIN_LOG_LEVEL\", \"3\")\nos.environ.setdefault(\"ABSL_MIN_LOG_LEVEL\", \"3\")\n\nimport cv2\nimport joblib\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport scipy\nimport sklearn\nimport statsmodels\nimport torch\nimport torch.nn.functional as torch_F\nimport torchvision\nimport torchxrayvision as xrv\nimport tensorflow as tf\n\nfrom imblearn.over_sampling import SMOTE\nfrom imblearn.under_sampling import RandomUnderSampler\nfrom scipy.optimize import minimize\nfrom scipy.special import expit\nfrom scipy.stats import binomtest, norm, t\nfrom sklearn.calibration import calibration_curve\nfrom sklearn.metrics import (\n    accuracy_score,\n    average_precision_score,\n    balanced_accuracy_score,\n    brier_score_loss,\n    confusion_matrix,\n    f1_score,\n    matthews_corrcoef,\n    precision_recall_curve,\n    precision_score,\n    recall_score,\n    roc_auc_score,\n    roc_curve,\n)\nfrom sklearn.model_selection import StratifiedGroupKFold\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.utils.class_weight import compute_class_weight\nfrom statsmodels.stats.multitest import multipletests\nfrom tensorflow import keras\nfrom tensorflow.keras import layers, models\nfrom tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau\nfrom tensorflow.keras.optimizers import Adam\nfrom torch.utils.data import DataLoader, Dataset\nfrom IPython.display import display\n\nfor _gpu in tf.config.list_physical_devices(\"GPU\"):\n    try:\n        tf.config.experimental.set_memory_growth(_gpu, True)\n    except Exception:\n        pass\n\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\n\nprint(\"Python:\", platform.python_version())\nprint(\"TensorFlow:\", tf.__version__)\nprint(\"PyTorch:\", torch.__version__)\nprint(\"TorchXRayVision:\", getattr(xrv, \"__version__\", \"unknown\"))\nprint(\"scikit-learn:\", sklearn.__version__)\nprint(\"SciPy:\", scipy.__version__)\nprint(\"statsmodels:\", statsmodels.__version__)\nprint(\"PyTorch CUDA available:\", torch.cuda.is_available())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:02:51.585218Z","iopub.execute_input":"2026-07-03T10:02:51.585528Z","iopub.status.idle":"2026-07-03T10:03:17.248961Z","shell.execute_reply.started":"2026-07-03T10:02:51.585504Z","shell.execute_reply":"2026-07-03T10:03:17.248116Z"}},"outputs":[],"execution_count":null},{"id":"38a780f3","cell_type":"markdown","source":"## 3. Final manuscript configuration\n","metadata":{}},{"id":"07b2ab66","cell_type":"code","source":"@dataclass\nclass Config:\n    seed: int = 42\n    seeds: Tuple[int, ...] = (42, 123, 2025, 31415, 27182)\n    img_size: int = 224\n    feature_batch_size: int = 32\n    mlp_batch_size: int = 64\n    mlp_epochs: int = 60\n    tuning_epochs: int = 16\n    scout_epochs: int = 18\n    scout_folds: int = 3\n    scout_replicates: int = 2\n    dropout: float = 0.35\n    learning_rate: float = 3e-4\n    hidden_units: int = 128\n    weight_cap: float = 4.0\n    noise_quantile: float = 0.97\n    noise_penalty: float = 1.25\n    target_sensitivity: float = 0.90\n    ohem_top_fraction: float = 0.50\n    focal_gamma: float = 2.0\n    xrv_weights: str = \"densenet121-res224-chex\"\n    xrv_feature_dim: int = 1024\n    num_workers: int = 2\n    bootstrap_iterations: int = 2000\n    output_dir: str = \"/kaggle/working/roof_dars_rsna_final_manuscript\"\n    kermany_root: Optional[str] = None\n    rsna_root: Optional[str] = None\n    max_dev: Optional[int] = None\n    max_internal_test: Optional[int] = None\n    max_rsna_per_class: Optional[int] = None\n    audit_exact_duplicates: bool = True\n    verbose: int = 0\n\n\ncfg = Config()\n\nMETHODS = [\n    \"Unbalanced-BCE\",\n    \"ClassWeighted-BCE\",\n    \"Focal-Loss\",\n    \"OHEM\",\n    \"Random-Under-Sampling\",\n    \"SMOTE\",\n    \"DARS-Original\",\n    \"OOF-Loss\",\n    \"ROOF-DARS\",\n]\n\nREFERENCE_METHOD = \"ROOF-DARS\"\nOUTPUT_DIR = Path(cfg.output_dir)\nCACHE_DIR = OUTPUT_DIR / \"feature_cache\"\nMODEL_DIR = OUTPUT_DIR / \"models\"\nPREDICTION_DIR = OUTPUT_DIR / \"predictions\"\nMETADATA_DIR = OUTPUT_DIR / \"metadata\"\nTABLE_DIR = OUTPUT_DIR / \"tables\"\nFIGURE_DIR = OUTPUT_DIR / \"figures\"\n\nfor directory in [\n    OUTPUT_DIR,\n    CACHE_DIR,\n    MODEL_DIR,\n    PREDICTION_DIR,\n    METADATA_DIR,\n    TABLE_DIR,\n    FIGURE_DIR,\n]:\n    directory.mkdir(parents=True, exist_ok=True)\n\nprint(\"Seeds:\", cfg.seeds)\nprint(\"Methods:\", METHODS)\nprint(\"Output directory:\", OUTPUT_DIR)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.250229Z","iopub.execute_input":"2026-07-03T10:03:17.250776Z","iopub.status.idle":"2026-07-03T10:03:17.261844Z","shell.execute_reply.started":"2026-07-03T10:03:17.250754Z","shell.execute_reply":"2026-07-03T10:03:17.26102Z"}},"outputs":[],"execution_count":null},{"id":"7b0cb01d","cell_type":"code","source":"def set_all_seeds(seed: int) -> None:\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    random.seed(seed)\n    np.random.seed(seed)\n    keras.utils.set_random_seed(seed)\n    torch.manual_seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed_all(seed)\n    try:\n        torch.use_deterministic_algorithms(True, warn_only=True)\n    except Exception:\n        pass\n    try:\n        tf.config.experimental.enable_op_determinism()\n    except Exception:\n        pass\n\nset_all_seeds(cfg.seeds[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.263692Z","iopub.execute_input":"2026-07-03T10:03:17.264097Z","iopub.status.idle":"2026-07-03T10:03:17.302036Z","shell.execute_reply.started":"2026-07-03T10:03:17.264035Z","shell.execute_reply":"2026-07-03T10:03:17.30133Z"}},"outputs":[],"execution_count":null},{"id":"0c84f313","cell_type":"markdown","source":"## 4. Kermany and full RSNA cohort construction\n","metadata":{}},{"id":"4c80dd2f","cell_type":"code","source":"IMAGE_EXTENSIONS = {\".jpg\", \".jpeg\", \".png\", \".bmp\", \".tif\", \".tiff\"}\n\n\ndef find_kermany_root(explicit: Optional[str]) -> Path:\n    if explicit:\n        root = Path(explicit)\n        if (root / \"train\").exists() and (root / \"test\").exists():\n            return root\n        if (root / \"chest_xray\" / \"train\").exists():\n            return root / \"chest_xray\"\n        raise FileNotFoundError(f\"Invalid Kermany root: {root}\")\n\n    candidates = [\n        Path(\"/kaggle/input/chest-xray-pneumonia/chest_xray\"),\n        Path(\"/kaggle/input/chest-xray-pneumonia/chest_xray/chest_xray\"),\n        Path(\"/kaggle/input/chest-xray-pneumonia\"),\n    ]\n    for root in candidates:\n        if (root / \"train\").exists() and (root / \"test\").exists():\n            return root\n\n    for root, dirs, _ in os.walk(\"/kaggle/input\"):\n        root_path = Path(root)\n        if (\n            \"train\" in dirs\n            and \"test\" in dirs\n            and (root_path / \"train\" / \"NORMAL\").exists()\n            and (root_path / \"train\" / \"PNEUMONIA\").exists()\n        ):\n            return root_path\n\n    raise FileNotFoundError(\n        \"Add the Kaggle dataset paultimothymooney/chest-xray-pneumonia.\"\n    )\n\n\ndef kermany_local_group(path: str) -> str:\n    stem = Path(path).stem.lower()\n    match = re.search(r\"person(\\d+)\", stem)\n    if match:\n        return f\"PNEU_person{match.group(1)}\"\n    match = re.search(r\"normal2-im-(\\d+)\", stem)\n    if match:\n        return f\"NORMAL2_IM_{match.group(1)}\"\n    match = re.search(r\"(?:^|_)im-(\\d+)\", stem)\n    if match:\n        return f\"NORMAL_IM_{match.group(1)}\"\n    return f\"FALLBACK_{Path(path).stem}\"\n\n\ndef collect_kermany_split(\n    split_dir: Path,\n    split_name: str,\n    namespace: str,\n) -> pd.DataFrame:\n    records = []\n    for class_name, label in [(\"NORMAL\", 0), (\"PNEUMONIA\", 1)]:\n        class_dir = split_dir / class_name\n        if not class_dir.exists():\n            continue\n        for path in sorted(class_dir.iterdir()):\n            if path.is_file() and path.suffix.lower() in IMAGE_EXTENSIONS:\n                local_group = kermany_local_group(str(path))\n                records.append(\n                    {\n                        \"path\": str(path),\n                        \"label\": label,\n                        \"class_name\": class_name,\n                        \"patient_id\": f\"KERMANY::{namespace}::{local_group}\",\n                        \"split\": split_name,\n                        \"source\": \"Kermany\",\n                    }\n                )\n    return pd.DataFrame(records)\n\n\ndef stratified_limit(\n    df: pd.DataFrame,\n    maximum: Optional[int],\n    seed: int,\n) -> pd.DataFrame:\n    if maximum is None or maximum >= len(df):\n        return df.reset_index(drop=True)\n    limited, _ = train_test_split(\n        df,\n        train_size=maximum,\n        stratify=df[\"label\"],\n        random_state=seed,\n    )\n    return limited.reset_index(drop=True)\n\n\ndef load_kermany(cfg: Config) -> Tuple[pd.DataFrame, pd.DataFrame]:\n    root = find_kermany_root(cfg.kermany_root)\n    train_df = collect_kermany_split(root / \"train\", \"train\", \"development\")\n    val_df = (\n        collect_kermany_split(root / \"val\", \"val\", \"development\")\n        if (root / \"val\").exists()\n        else pd.DataFrame(columns=train_df.columns)\n    )\n    test_df = collect_kermany_split(\n        root / \"test\", \"internal_test\", \"internal_test\"\n    )\n\n    dev_df = pd.concat([train_df, val_df], ignore_index=True)\n    dev_df = stratified_limit(dev_df, cfg.max_dev, cfg.seed)\n    test_df = stratified_limit(test_df, cfg.max_internal_test, cfg.seed)\n\n    print(f\"Kermany root: {root}\")\n    print(\n        f\"Development: {len(dev_df)} images | \"\n        f\"{dev_df['patient_id'].nunique()} groups | \"\n        f\"Normal={(dev_df.label == 0).sum()} | \"\n        f\"Pneumonia={(dev_df.label == 1).sum()}\"\n    )\n    print(\n        f\"Internal test: {len(test_df)} images | \"\n        f\"{test_df['patient_id'].nunique()} groups | \"\n        f\"Normal={(test_df.label == 0).sum()} | \"\n        f\"Pneumonia={(test_df.label == 1).sum()}\"\n    )\n    return dev_df, test_df\n\n\ndef find_rsna_root(explicit: Optional[str]) -> Path:\n    candidates = []\n    if explicit:\n        candidates.append(Path(explicit))\n    candidates.extend(\n        [\n            Path(\"/kaggle/input/rsna-pneumonia-detection-challenge\"),\n            Path(\"/kaggle/input/rsna-pneumonia-detection-challenge/rsna-pneumonia-detection-challenge\"),\n        ]\n    )\n\n    for root in candidates:\n        if (\n            root.exists()\n            and (root / \"stage_2_detailed_class_info.csv\").exists()\n            and (root / \"stage_2_train_labels.csv\").exists()\n            and (root / \"stage_2_train_images\").exists()\n        ):\n            return root\n\n    for csv_path in Path(\"/kaggle/input\").rglob(\"stage_2_detailed_class_info.csv\"):\n        root = csv_path.parent\n        if (\n            (root / \"stage_2_train_labels.csv\").exists()\n            and (root / \"stage_2_train_images\").exists()\n        ):\n            return root\n\n    raise FileNotFoundError(\n        \"RSNA data not found. Add the Kaggle competition \"\n        \"'rsna-pneumonia-detection-challenge' and accept its rules.\"\n    )\n\n\ndef load_rsna_external(cfg: Config) -> pd.DataFrame:\n    root = find_rsna_root(cfg.rsna_root)\n    class_info_path = root / \"stage_2_detailed_class_info.csv\"\n    labels_path = root / \"stage_2_train_labels.csv\"\n    images_dir = root / \"stage_2_train_images\"\n\n    class_info = pd.read_csv(class_info_path)\n    labels = pd.read_csv(labels_path)\n\n    required_class_cols = {\"patientId\", \"class\"}\n    required_label_cols = {\"patientId\", \"Target\"}\n    if not required_class_cols.issubset(class_info.columns):\n        raise ValueError(\n            f\"Unexpected RSNA class-info columns: {list(class_info.columns)}\"\n        )\n    if not required_label_cols.issubset(labels.columns):\n        raise ValueError(\n            f\"Unexpected RSNA label columns: {list(labels.columns)}\"\n        )\n\n    class_info = class_info[[\"patientId\", \"class\"]].drop_duplicates(\"patientId\")\n    target_by_patient = (\n        labels.groupby(\"patientId\", as_index=False)[\"Target\"].max()\n    )\n    bbox_count = (\n        labels.assign(has_box=labels[\"Target\"].eq(1).astype(int))\n        .groupby(\"patientId\", as_index=False)[\"has_box\"]\n        .sum()\n        .rename(columns={\"has_box\": \"bbox_count\"})\n    )\n\n    bbox_map: Dict[str, list[dict[str, float]]] = {}\n    positive_box_rows = labels[labels[\"Target\"].eq(1)].copy()\n    for patient_id, group in positive_box_rows.groupby(\"patientId\"):\n        boxes = []\n        for row in group.itertuples(index=False):\n            coordinates = [\n                getattr(row, \"x\", np.nan),\n                getattr(row, \"y\", np.nan),\n                getattr(row, \"width\", np.nan),\n                getattr(row, \"height\", np.nan),\n            ]\n            if all(pd.notna(value) for value in coordinates):\n                boxes.append(\n                    {\n                        \"x\": float(coordinates[0]),\n                        \"y\": float(coordinates[1]),\n                        \"width\": float(coordinates[2]),\n                        \"height\": float(coordinates[3]),\n                    }\n                )\n        bbox_map[str(patient_id)] = boxes\n\n    cohort = (\n        class_info.merge(target_by_patient, on=\"patientId\", how=\"left\")\n        .merge(bbox_count, on=\"patientId\", how=\"left\")\n    )\n    cohort[\"bbox_json\"] = cohort[\"patientId\"].astype(str).map(\n        lambda patient_id: json.dumps(bbox_map.get(patient_id, []))\n    )\n    cohort[\"Target\"] = cohort[\"Target\"].fillna(0).astype(int)\n    cohort[\"bbox_count\"] = cohort[\"bbox_count\"].fillna(0).astype(int)\n\n    cohort = cohort[\n        cohort[\"class\"].isin([\"Lung Opacity\", \"Normal\"])\n    ].copy()\n    cohort[\"label\"] = cohort[\"class\"].map(\n        {\"Normal\": 0, \"Lung Opacity\": 1}\n    ).astype(int)\n\n    mismatch = cohort[\n        ((cohort[\"class\"] == \"Lung Opacity\") & (cohort[\"Target\"] != 1))\n        | ((cohort[\"class\"] == \"Normal\") & (cohort[\"Target\"] != 0))\n    ]\n    if len(mismatch):\n        print(\n            f\"WARNING: excluding {len(mismatch)} RSNA rows with inconsistent \"\n            \"class and Target labels.\"\n        )\n        cohort = cohort.drop(index=mismatch.index)\n\n    cohort[\"path\"] = cohort[\"patientId\"].map(\n        lambda pid: str(images_dir / f\"{pid}.dcm\")\n    )\n    missing_mask = ~cohort[\"path\"].map(lambda value: Path(value).exists())\n    missing_count = int(missing_mask.sum())\n    cohort = cohort[~missing_mask].copy()\n\n    if cfg.max_rsna_per_class is not None:\n        parts = []\n        for label in [0, 1]:\n            class_df = cohort[cohort[\"label\"] == label]\n            n = min(cfg.max_rsna_per_class, len(class_df))\n            parts.append(class_df.sample(n=n, random_state=cfg.seed + label))\n        cohort = pd.concat(parts, ignore_index=True)\n\n    cohort[\"class_name\"] = cohort[\"class\"].map(\n        {\"Normal\": \"NORMAL\", \"Lung Opacity\": \"LUNG_OPACITY\"}\n    )\n    cohort[\"patient_id\"] = \"RSNA::\" + cohort[\"patientId\"].astype(str)\n    cohort[\"split\"] = \"external\"\n    cohort[\"source\"] = \"RSNA_Pneumonia_Detection_Challenge\"\n    cohort[\"protocol\"] = \"Lung Opacity versus Normal; abnormal-no-opacity excluded\"\n\n    output = cohort[\n        [\n            \"path\",\n            \"label\",\n            \"class_name\",\n            \"patient_id\",\n            \"split\",\n            \"source\",\n            \"protocol\",\n            \"patientId\",\n            \"class\",\n            \"Target\",\n            \"bbox_count\",\n            \"bbox_json\",\n        ]\n    ].reset_index(drop=True)\n\n    if output[\"patient_id\"].duplicated().any():\n        raise RuntimeError(\"RSNA external cohort is not patient unique.\")\n\n    print(f\"RSNA root: {root}\")\n    print(\n        f\"RSNA external: {len(output)} patients | \"\n        f\"Normal={(output.label == 0).sum()} | \"\n        f\"Lung Opacity={(output.label == 1).sum()} | \"\n        f\"Missing DICOM={missing_count}\"\n    )\n    print(\"Excluded class: No Lung Opacity / Not Normal\")\n    return output\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.303019Z","iopub.execute_input":"2026-07-03T10:03:17.303858Z","iopub.status.idle":"2026-07-03T10:03:17.333474Z","shell.execute_reply.started":"2026-07-03T10:03:17.303836Z","shell.execute_reply":"2026-07-03T10:03:17.332657Z"}},"outputs":[],"execution_count":null},{"id":"6e3ee473","cell_type":"markdown","source":"## 5. DICOM-aware preprocessing and CheXpert DenseNet121 features\n","metadata":{}},{"id":"0425628c","cell_type":"code","source":"def read_xray_gray(path: str) -> np.ndarray:\n    path_obj = Path(path)\n\n    if path_obj.suffix.lower() == \".dcm\":\n        dataset = pydicom.dcmread(str(path_obj))\n        image = dataset.pixel_array.astype(np.float32)\n\n        if image.ndim == 3:\n            image = image[0]\n\n        slope = float(getattr(dataset, \"RescaleSlope\", 1.0))\n        intercept = float(getattr(dataset, \"RescaleIntercept\", 0.0))\n        image = image * slope + intercept\n\n        finite = image[np.isfinite(image)]\n        if finite.size == 0:\n            raise ValueError(f\"DICOM contains no finite pixels: {path}\")\n\n        lower, upper = np.percentile(finite, [0.5, 99.5])\n        if upper <= lower:\n            lower, upper = float(finite.min()), float(finite.max())\n\n        image = np.clip(image, lower, upper)\n        image = (image - lower) / (upper - lower + 1e-8)\n\n        if str(getattr(dataset, \"PhotometricInterpretation\", \"\")).upper() == \"MONOCHROME1\":\n            image = 1.0 - image\n\n        return np.round(image * 255.0).astype(np.uint8)\n\n    image = cv2.imread(str(path_obj), cv2.IMREAD_GRAYSCALE)\n    if image is None:\n        raise ValueError(f\"Could not read image: {path}\")\n    return image\n\ndef crop_black_borders(\n    gray: np.ndarray,\n    threshold: int = 8,\n    pad: int = 8,\n) -> np.ndarray:\n    coords = np.argwhere(gray > threshold)\n    if coords.size == 0:\n        return gray\n    y0, x0 = coords.min(axis=0)\n    y1, x1 = coords.max(axis=0) + 1\n    return gray[\n        max(0, y0 - pad) : min(gray.shape[0], y1 + pad),\n        max(0, x0 - pad) : min(gray.shape[1], x1 + pad),\n    ]\n\n\ndef remove_corner_markers(gray: np.ndarray) -> np.ndarray:\n    h, w = gray.shape\n    region = np.zeros_like(gray, dtype=bool)\n    region[: int(0.22 * h), : int(0.28 * w)] = True\n    region[: int(0.22 * h), int(0.72 * w) :] = True\n    mask = np.zeros_like(gray, dtype=np.uint8)\n    mask[(gray > 220) & region] = 255\n    mask = cv2.morphologyEx(\n        mask, cv2.MORPH_OPEN, np.ones((3, 3), np.uint8)\n    )\n    if not mask.any():\n        return gray\n    return cv2.inpaint(gray, mask, 3, cv2.INPAINT_TELEA)\n\n\ndef deterministic_gray_preprocess(path: str) -> np.ndarray:\n    gray = read_xray_gray(path)\n    gray = crop_black_borders(gray)\n    gray = remove_corner_markers(gray)\n    h, w = gray.shape\n    gray = gray[\n        int(0.02 * h) : int(0.99 * h),\n        int(0.03 * w) : int(0.97 * w),\n    ]\n    gray = cv2.createCLAHE(\n        clipLimit=2.0, tileGridSize=(8, 8)\n    ).apply(gray)\n    return gray.astype(np.uint8)\n\n\nclass XRayPathDataset(Dataset):\n    def __init__(\n        self,\n        paths: Sequence[str],\n        labels: Sequence[int],\n        transform: Any,\n    ):\n        self.paths = list(paths)\n        self.labels = np.asarray(labels, dtype=np.float32)\n        self.transform = transform\n\n    def __len__(self) -> int:\n        return len(self.paths)\n\n    def __getitem__(self, index: int):\n        gray = deterministic_gray_preprocess(self.paths[index])\n        image = xrv.datasets.normalize(gray, 255)\n        image = image[None, ...]\n        image = self.transform(image)\n        image = np.ascontiguousarray(image, dtype=np.float32)\n        return (\n            torch.from_numpy(image),\n            torch.tensor(self.labels[index], dtype=torch.float32),\n        )\n\n\nclass XRVFeatureExtractor:\n    def __init__(self, cfg: Config):\n        self.cfg = cfg\n        self.device = torch.device(\n            \"cuda\" if torch.cuda.is_available() else \"cpu\"\n        )\n        self.transform = torchvision.transforms.Compose(\n            [\n                xrv.datasets.XRayCenterCrop(),\n                xrv.datasets.XRayResizer(cfg.img_size),\n            ]\n        )\n        if not hasattr(xrv.models, \"DenseNet\"):\n            available = [\n                name for name in dir(xrv.models)\n                if not name.startswith(\"_\")\n            ]\n            raise AttributeError(\n                \"TorchXRayVision does not expose xrv.models.DenseNet. \"\n                f\"Installed version: {getattr(xrv, '__version__', 'unknown')}. \"\n                f\"Available model symbols: {available}\"\n            )\n\n        self.model = xrv.models.DenseNet(weights=cfg.xrv_weights)\n        self.model = self.model.to(self.device)\n        self.model.eval()\n        for parameter in self.model.parameters():\n            parameter.requires_grad = False\n        print(\n            f\"Feature extractor: {cfg.xrv_weights} on {self.device} | \"\n            f\"TorchXRayVision={getattr(xrv, '__version__', 'unknown')}\"\n        )\n\n    def _cache_key(self, paths: Sequence[str], name: str) -> str:\n        payload = json.dumps(\n            {\n                \"name\": name,\n                \"paths\": list(map(str, paths)),\n                \"weights\": self.cfg.xrv_weights,\n                \"img_size\": self.cfg.img_size,\n            },\n            sort_keys=True,\n        )\n        return hashlib.sha256(payload.encode()).hexdigest()[:20]\n\n    def extract(\n        self,\n        paths: Sequence[str],\n        labels: Sequence[int],\n        cache_dir: Path,\n        cache_name: str,\n    ) -> np.ndarray:\n        cache_dir.mkdir(parents=True, exist_ok=True)\n        key = self._cache_key(paths, cache_name)\n        cache_path = cache_dir / f\"{cache_name}_{key}.npy\"\n        if cache_path.exists():\n            print(f\"Loading cached features: {cache_path}\")\n            return np.load(cache_path)\n\n        dataset = XRayPathDataset(paths, labels, self.transform)\n        loader = DataLoader(\n            dataset,\n            batch_size=self.cfg.feature_batch_size,\n            shuffle=False,\n            num_workers=self.cfg.num_workers,\n            pin_memory=self.device.type == \"cuda\",\n            persistent_workers=self.cfg.num_workers > 0,\n        )\n        matrices = []\n        with torch.inference_mode():\n            for batch_index, (images, _) in enumerate(loader, start=1):\n                images = images.to(self.device, non_blocking=True)\n                maps = self.model.features(images)\n                maps = torch_F.relu(maps, inplace=False)\n                pooled = torch_F.adaptive_avg_pool2d(\n                    maps, (1, 1)\n                ).flatten(1)\n                matrices.append(\n                    pooled.detach().cpu().numpy().astype(np.float32)\n                )\n                if batch_index % 25 == 0:\n                    print(\n                        f\"Feature batches completed: \"\n                        f\"{batch_index}/{len(loader)}\"\n                    )\n\n        features = np.concatenate(matrices, axis=0)\n        if features.shape[1] != self.cfg.xrv_feature_dim:\n            raise RuntimeError(\n                f\"Expected {self.cfg.xrv_feature_dim} features, \"\n                f\"received {features.shape[1]}.\"\n            )\n        np.save(cache_path, features)\n        print(f\"Saved features: {cache_path}\")\n        return features\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.334435Z","iopub.execute_input":"2026-07-03T10:03:17.334705Z","iopub.status.idle":"2026-07-03T10:03:17.358417Z","shell.execute_reply.started":"2026-07-03T10:03:17.334675Z","shell.execute_reply":"2026-07-03T10:03:17.357543Z"}},"outputs":[],"execution_count":null},{"id":"ae91c925","cell_type":"markdown","source":"## 6. Duplicate-audit utilities\n","metadata":{}},{"id":"43879e2d","cell_type":"code","source":"def standardized_pixel_hash(path: str) -> str:\n    gray = deterministic_gray_preprocess(path)\n    gray = cv2.resize(gray, (256, 256), interpolation=cv2.INTER_AREA)\n    return hashlib.sha256(gray.tobytes()).hexdigest()\n\n\ndef exact_duplicate_audit(\n    reference_df: pd.DataFrame,\n    candidate_df: pd.DataFrame,\n    output_csv: Path,\n) -> pd.DataFrame:\n    reference_map: Dict[str, str] = {}\n    for path in reference_df[\"path\"]:\n        reference_map.setdefault(standardized_pixel_hash(path), str(path))\n\n    records = []\n    for path in candidate_df[\"path\"]:\n        digest = standardized_pixel_hash(path)\n        if digest in reference_map:\n            records.append(\n                {\n                    \"reference_path\": reference_map[digest],\n                    \"candidate_path\": str(path),\n                    \"standardized_pixel_sha256\": digest,\n                }\n            )\n\n    report = pd.DataFrame(records)\n    report.to_csv(output_csv, index=False)\n    return report\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.359359Z","iopub.execute_input":"2026-07-03T10:03:17.359844Z","iopub.status.idle":"2026-07-03T10:03:17.373845Z","shell.execute_reply.started":"2026-07-03T10:03:17.359823Z","shell.execute_reply":"2026-07-03T10:03:17.373313Z"}},"outputs":[],"execution_count":null},{"id":"de4f8cb7","cell_type":"markdown","source":"## 7. Classification head, calibration, thresholds, and metrics\n","metadata":{}},{"id":"3686600d","cell_type":"code","source":"def binary_focal_loss(gamma: float = 2.0):\n    def loss(y_true, y_pred):\n        y_true = tf.cast(y_true, tf.float32)\n        y_pred = tf.clip_by_value(tf.cast(y_pred, tf.float32), 1e-7, 1.0 - 1e-7)\n        positive = -y_true * tf.pow(1.0 - y_pred, gamma) * tf.math.log(y_pred)\n        negative = -(1.0 - y_true) * tf.pow(y_pred, gamma) * tf.math.log(1.0 - y_pred)\n        return tf.reduce_mean(positive + negative)\n    return loss\n\n\ndef ohem_binary_crossentropy(top_fraction: float = 0.50):\n    def loss(y_true, y_pred):\n        y_true = tf.cast(y_true, tf.float32)\n        per_example = tf.keras.backend.binary_crossentropy(y_true, y_pred)\n        per_example = tf.reshape(per_example, [-1])\n        count = tf.size(per_example)\n        k = tf.maximum(\n            1,\n            tf.cast(\n                tf.math.ceil(tf.cast(count, tf.float32) * top_fraction),\n                tf.int32,\n            ),\n        )\n        hardest, _ = tf.math.top_k(per_example, k=k, sorted=False)\n        return tf.reduce_mean(hardest)\n    return loss\n\n\ndef build_mlp(\n    cfg: Config,\n    input_dim: int,\n    loss_function: Any = \"binary_crossentropy\",\n) -> keras.Model:\n    inputs = layers.Input(shape=(input_dim,), name=\"features\")\n    x = layers.BatchNormalization(name=\"head_batch_norm\")(inputs)\n    x = layers.Dense(cfg.hidden_units, activation=\"relu\", name=\"dense_1\")(x)\n    x = layers.Dropout(cfg.dropout, name=\"dropout_1\")(x)\n    x = layers.Dense(\n        max(32, cfg.hidden_units // 2),\n        activation=\"relu\",\n        name=\"dense_2\",\n    )(x)\n    x = layers.Dropout(cfg.dropout / 2.0, name=\"dropout_2\")(x)\n    logits = layers.Dense(1, activation=None, name=\"logit\")(x)\n    probabilities = layers.Activation(\"sigmoid\", name=\"probability\")(logits)\n\n    model = models.Model(inputs, probabilities, name=\"roof_dars_mlp\")\n    model.compile(\n        optimizer=Adam(learning_rate=cfg.learning_rate),\n        loss=loss_function,\n        metrics=[\n            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n            keras.metrics.AUC(name=\"roc_auc\", curve=\"ROC\"),\n            keras.metrics.AUC(name=\"pr_auc\", curve=\"PR\"),\n        ],\n    )\n    return model\n\n\ndef class_weight_dict(y: np.ndarray) -> Dict[int, float]:\n    classes = np.array([0, 1])\n    weights = compute_class_weight(\n        class_weight=\"balanced\",\n        classes=classes,\n        y=np.asarray(y).astype(int),\n    )\n    return {int(cls): float(weight) for cls, weight in zip(classes, weights)}\n\n\ndef fit_mlp(\n    cfg: Config,\n    X_train: np.ndarray,\n    y_train: np.ndarray,\n    X_valid: np.ndarray,\n    y_valid: np.ndarray,\n    seed: int,\n    epochs: Optional[int] = None,\n    class_weight: Optional[Dict[int, float]] = None,\n    loss_function: Any = \"binary_crossentropy\",\n    fixed_epochs: bool = False,\n) -> Tuple[keras.Model, pd.DataFrame]:\n    set_all_seeds(seed)\n    model = build_mlp(cfg, X_train.shape[1], loss_function=loss_function)\n\n    callbacks = []\n    if not fixed_epochs:\n        callbacks = [\n            EarlyStopping(\n                monitor=\"val_pr_auc\",\n                mode=\"max\",\n                patience=7,\n                restore_best_weights=True,\n                verbose=0,\n            ),\n            ReduceLROnPlateau(\n                monitor=\"val_pr_auc\",\n                mode=\"max\",\n                factor=0.4,\n                patience=3,\n                min_lr=1e-7,\n                verbose=0,\n            ),\n        ]\n\n    history = model.fit(\n        X_train,\n        y_train,\n        validation_data=(X_valid, y_valid),\n        epochs=epochs or cfg.mlp_epochs,\n        batch_size=cfg.mlp_batch_size,\n        callbacks=callbacks,\n        class_weight=class_weight,\n        verbose=cfg.verbose,\n    )\n    return model, pd.DataFrame(history.history)\n\n\ndef predict_logits(model: keras.Model, X: np.ndarray) -> np.ndarray:\n    logit_model = models.Model(model.input, model.get_layer(\"logit\").output)\n    return logit_model.predict(X, batch_size=256, verbose=0).reshape(-1)\n\n\ndef fit_affine_calibrator(\n    y_true: np.ndarray,\n    logits: np.ndarray,\n) -> Tuple[float, float]:\n    y = np.asarray(y_true, dtype=np.float64)\n    z = np.asarray(logits, dtype=np.float64)\n\n    def objective(parameters: np.ndarray) -> float:\n        log_slope, intercept = parameters\n        slope = np.exp(log_slope)\n        probabilities = expit(slope * z + intercept)\n        probabilities = np.clip(probabilities, 1e-7, 1.0 - 1e-7)\n        return float(\n            -np.mean(\n                y * np.log(probabilities)\n                + (1.0 - y) * np.log(1.0 - probabilities)\n            )\n        )\n\n    result = minimize(\n        objective,\n        x0=np.array([0.0, 0.0]),\n        method=\"L-BFGS-B\",\n        bounds=[(-4.0, 4.0), (-10.0, 10.0)],\n    )\n    if not result.success:\n        return 1.0, 0.0\n    return float(np.exp(result.x[0])), float(result.x[1])\n\n\ndef calibrated_probabilities(\n    logits: np.ndarray,\n    slope: float,\n    intercept: float,\n) -> np.ndarray:\n    return expit(slope * np.asarray(logits) + intercept)\n\n\ndef expected_calibration_error(\n    y_true: np.ndarray,\n    probabilities: np.ndarray,\n    bins: int = 15,\n) -> float:\n    y = np.asarray(y_true)\n    p = np.asarray(probabilities)\n    edges = np.linspace(0.0, 1.0, bins + 1)\n    value = 0.0\n    for lower, upper in zip(edges[:-1], edges[1:]):\n        mask = (p > lower) & (p <= upper)\n        if mask.any():\n            value += mask.mean() * abs(p[mask].mean() - y[mask].mean())\n    return float(value)\n\n\ndef tune_threshold_mcc(\n    y_true: np.ndarray,\n    probabilities: np.ndarray,\n    target_sensitivity: float,\n) -> Tuple[float, Dict[str, float]]:\n    y = np.asarray(y_true).astype(int)\n    p = np.asarray(probabilities, dtype=float)\n    candidates = []\n\n    for threshold in np.linspace(0.01, 0.99, 393):\n        predictions = (p >= threshold).astype(int)\n        sensitivity = recall_score(y, predictions, zero_division=0)\n        specificity = recall_score(y, predictions, pos_label=0, zero_division=0)\n        mcc = matthews_corrcoef(y, predictions)\n        balanced = balanced_accuracy_score(y, predictions)\n        candidates.append(\n            {\n                \"threshold\": float(threshold),\n                \"sensitivity\": float(sensitivity),\n                \"specificity\": float(specificity),\n                \"mcc\": float(mcc),\n                \"balanced_accuracy\": float(balanced),\n            }\n        )\n\n    eligible = [\n        item for item in candidates\n        if item[\"sensitivity\"] >= target_sensitivity\n    ]\n    pool = eligible if eligible else candidates\n    best = max(\n        pool,\n        key=lambda item: (\n            item[\"mcc\"],\n            item[\"balanced_accuracy\"],\n            item[\"specificity\"],\n        ),\n    )\n    return float(best[\"threshold\"]), best\n\n\ndef compute_metrics(\n    y_true: np.ndarray,\n    probabilities: np.ndarray,\n    threshold: float,\n) -> Dict[str, float]:\n    y = np.asarray(y_true).astype(int)\n    p = np.asarray(probabilities, dtype=float)\n    predictions = (p >= threshold).astype(int)\n    tn, fp, fn, tp = confusion_matrix(y, predictions, labels=[0, 1]).ravel()\n    sensitivity = tp / (tp + fn) if tp + fn else 0.0\n    specificity = tn / (tn + fp) if tn + fp else 0.0\n\n    return {\n        \"accuracy\": float(accuracy_score(y, predictions)),\n        \"balanced_accuracy\": float(\n            balanced_accuracy_score(y, predictions)\n        ),\n        \"sensitivity\": float(sensitivity),\n        \"specificity\": float(specificity),\n        \"precision\": float(\n            precision_score(y, predictions, zero_division=0)\n        ),\n        \"f1\": float(f1_score(y, predictions, zero_division=0)),\n        \"mcc\": float(matthews_corrcoef(y, predictions)),\n        \"roc_auc\": float(roc_auc_score(y, p)),\n        \"pr_auc\": float(average_precision_score(y, p)),\n        \"brier\": float(brier_score_loss(y, p)),\n        \"ece\": float(expected_calibration_error(y, p)),\n        \"threshold\": float(threshold),\n        \"tn\": int(tn),\n        \"fp\": int(fp),\n        \"fn\": int(fn),\n        \"tp\": int(tp),\n        \"predicted_positive_rate\": float(predictions.mean()),\n        \"true_positive_prevalence\": float(y.mean()),\n    }\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.374742Z","iopub.execute_input":"2026-07-03T10:03:17.37515Z","iopub.status.idle":"2026-07-03T10:03:17.402391Z","shell.execute_reply.started":"2026-07-03T10:03:17.375129Z","shell.execute_reply":"2026-07-03T10:03:17.401432Z"}},"outputs":[],"execution_count":null},{"id":"f5caaefb","cell_type":"markdown","source":"## 8. Original DARS, OOF ablation, and proposed ROOF-DARS\n","metadata":{}},{"id":"8930d286","cell_type":"code","source":"def bce_numpy(y: np.ndarray, p: np.ndarray) -> np.ndarray:\n    p = np.clip(np.asarray(p), 1e-7, 1.0 - 1e-7)\n    y = np.asarray(y)\n    return -(y * np.log(p) + (1.0 - y) * np.log(1.0 - p))\n\n\ndef class_quantile_normalize(\n    values: np.ndarray,\n    y: np.ndarray,\n    lower_q: float = 0.10,\n    upper_q: float = 0.90,\n) -> np.ndarray:\n    output = np.zeros_like(values, dtype=float)\n    for cls in [0, 1]:\n        mask = y == cls\n        low, high = np.quantile(values[mask], [lower_q, upper_q])\n        output[mask] = np.clip(\n            (values[mask] - low) / (high - low + 1e-8),\n            0.0,\n            1.0,\n        )\n    return output\n\n\ndef generate_oof_scout_predictions(\n    cfg: Config,\n    X_train: np.ndarray,\n    y_train: np.ndarray,\n    groups_train: np.ndarray,\n    seed: int,\n) -> Tuple[np.ndarray, np.ndarray]:\n    cache_path = CACHE_DIR / f\"oof_scout_seed_{seed}.npz\"\n    if cache_path.exists():\n        cached = np.load(cache_path)\n        return cached[\"mean_probability\"], cached[\"disagreement\"]\n\n    replicate_predictions = []\n    for replicate in range(cfg.scout_replicates):\n        scout_seed = seed + 10000 * (replicate + 1)\n        print(f\"OOF scout replicate {replicate + 1}/{cfg.scout_replicates}\")\n        splitter = StratifiedGroupKFold(\n            n_splits=cfg.scout_folds,\n            shuffle=True,\n            random_state=scout_seed,\n        )\n        oof = np.full(len(y_train), np.nan, dtype=float)\n\n        for fold, (fit_idx, eval_idx) in enumerate(\n            splitter.split(X_train, y_train, groups_train),\n            start=1,\n        ):\n            print(\n                f\"  scout seed={scout_seed}, \"\n                f\"fold={fold}/{cfg.scout_folds}\"\n            )\n            model, _ = fit_mlp(\n                cfg,\n                X_train[fit_idx],\n                y_train[fit_idx],\n                X_train[eval_idx],\n                y_train[eval_idx],\n                seed=scout_seed + fold,\n                epochs=cfg.scout_epochs,\n                class_weight=class_weight_dict(y_train[fit_idx]),\n                fixed_epochs=True,\n            )\n            oof[eval_idx] = expit(predict_logits(model, X_train[eval_idx]))\n            del model\n            keras.backend.clear_session()\n            gc.collect()\n\n        if np.isnan(oof).any():\n            raise RuntimeError(\"Incomplete OOF scout predictions.\")\n        replicate_predictions.append(oof)\n\n    matrix = np.vstack(replicate_predictions)\n    mean_probability = matrix.mean(axis=0)\n    disagreement = matrix.var(axis=0)\n    np.savez_compressed(\n        cache_path,\n        mean_probability=mean_probability,\n        disagreement=disagreement,\n    )\n    return mean_probability, disagreement\n\n\ndef generate_in_sample_scout_predictions(\n    cfg: Config,\n    X_train: np.ndarray,\n    y_train: np.ndarray,\n    X_tune: np.ndarray,\n    y_tune: np.ndarray,\n    seed: int,\n) -> np.ndarray:\n    model, _ = fit_mlp(\n        cfg,\n        X_train,\n        y_train,\n        X_tune,\n        y_tune,\n        seed=seed + 7000,\n        epochs=cfg.scout_epochs,\n        class_weight=class_weight_dict(y_train),\n    )\n    probabilities = expit(predict_logits(model, X_train))\n    del model\n    keras.backend.clear_session()\n    gc.collect()\n    return probabilities\n\n\n@dataclass(frozen=True)\nclass RoofParameters:\n    alpha_loss: float\n    beta_misclassified: float\n    gamma_uncertainty: float\n    delta_disagreement: float\n    target_ratio: float\n    noise_penalty: float\n    weight_cap: float\n\n\ndef roof_parameter_grid(cfg: Config) -> list[RoofParameters]:\n    coefficients = [\n        (1.50, 0.25, 1.00, 0.50),\n        (1.50, 0.50, 1.00, 0.50),\n        (2.00, 0.25, 0.75, 0.50),\n        (2.00, 0.50, 0.75, 0.25),\n    ]\n    ratios = [0.35, 0.50, 0.65, 0.80]\n    return [\n        RoofParameters(\n            alpha_loss=a,\n            beta_misclassified=b,\n            gamma_uncertainty=g,\n            delta_disagreement=d,\n            target_ratio=ratio,\n            noise_penalty=cfg.noise_penalty,\n            weight_cap=cfg.weight_cap,\n        )\n        for a, b, g, d in coefficients\n        for ratio in ratios\n    ]\n\n\ndef roof_difficulty(\n    cfg: Config,\n    y: np.ndarray,\n    probabilities: np.ndarray,\n    disagreement: np.ndarray,\n    params: RoofParameters,\n) -> Tuple[np.ndarray, Dict[str, Any]]:\n    y = np.asarray(y).astype(int)\n    p = np.asarray(probabilities, dtype=float)\n    disagreement = np.asarray(disagreement, dtype=float)\n\n    loss = bce_numpy(y, p)\n    loss_norm = class_quantile_normalize(loss, y)\n    misclassified = ((p >= 0.5).astype(int) != y).astype(float)\n    uncertainty = np.clip(1.0 - 2.0 * np.abs(p - 0.5), 0.0, 1.0)\n    disagreement_norm = class_quantile_normalize(\n        disagreement, y, 0.05, 0.95\n    )\n\n    noise_flag = np.zeros(len(y), dtype=float)\n    for cls in [0, 1]:\n        mask = y == cls\n        extreme_cutoff = np.quantile(loss[mask], cfg.noise_quantile)\n        extreme_loss = loss > extreme_cutoff\n        confidently_wrong = (\n            ((y == 1) & (p < 0.15))\n            | ((y == 0) & (p > 0.85))\n        )\n        noise_flag[mask] = (\n            extreme_loss[mask] & confidently_wrong[mask]\n        ).astype(float)\n\n    raw = (\n        1.0\n        + params.alpha_loss * loss_norm\n        + params.beta_misclassified * misclassified\n        + params.gamma_uncertainty * uncertainty\n        + params.delta_disagreement * disagreement_norm\n    )\n    weights = np.clip(\n        raw * np.exp(-params.noise_penalty * noise_flag),\n        0.5,\n        params.weight_cap,\n    )\n\n    details = {\n        \"mean_oof_loss\": float(loss.mean()),\n        \"oof_misclassified\": int(misclassified.sum()),\n        \"mean_uncertainty\": float(uncertainty.mean()),\n        \"mean_disagreement\": float(disagreement.mean()),\n        \"probable_noise_count\": int(noise_flag.sum()),\n        \"mean_sampling_weight\": float(weights.mean()),\n        \"max_sampling_weight\": float(weights.max()),\n    }\n    return weights, details\n\n\ndef roof_resample(\n    X: np.ndarray,\n    y: np.ndarray,\n    weights: np.ndarray,\n    target_ratio: float,\n    seed: int,\n) -> Tuple[np.ndarray, np.ndarray, Dict[str, int]]:\n    rng = np.random.default_rng(seed)\n    y = np.asarray(y).astype(int)\n    counts = {cls: int((y == cls).sum()) for cls in [0, 1]}\n    majority_class = max(counts, key=counts.get)\n    minority_class = min(counts, key=counts.get)\n\n    majority_idx = np.where(y == majority_class)[0]\n    minority_idx = np.where(y == minority_class)[0]\n    target_minority = min(\n        counts[majority_class],\n        max(\n            counts[minority_class],\n            int(round(target_ratio * counts[majority_class])),\n        ),\n    )\n\n    selected = list(majority_idx) + list(minority_idx)\n    extras_needed = target_minority - len(minority_idx)\n    if extras_needed > 0:\n        probabilities = weights[minority_idx].astype(float)\n        probabilities /= probabilities.sum()\n        extras = rng.choice(\n            minority_idx,\n            size=extras_needed,\n            replace=True,\n            p=probabilities,\n        )\n        selected.extend(extras.tolist())\n\n    selected = np.asarray(selected)\n    rng.shuffle(selected)\n    details = {\n        \"majority_class\": int(majority_class),\n        \"minority_class\": int(minority_class),\n        \"majority_original\": counts[majority_class],\n        \"minority_original\": counts[minority_class],\n        \"minority_after_resampling\": int(target_minority),\n        \"total_after_resampling\": int(len(selected)),\n    }\n    return X[selected], y[selected], details\n\n\ndef composite_tuning_score(metrics: Dict[str, float]) -> float:\n    if metrics[\"sensitivity\"] < cfg.target_sensitivity:\n        return -1e6\n    return float(\n        0.35 * metrics[\"mcc\"]\n        + 0.25 * metrics[\"balanced_accuracy\"]\n        + 0.15 * metrics[\"specificity\"]\n        + 0.10 * metrics[\"pr_auc\"]\n        + 0.10 * metrics[\"roc_auc\"]\n        - 0.025 * metrics[\"ece\"]\n        - 0.025 * metrics[\"brier\"]\n    )\n\n\ndef tune_resampling_candidate(\n    cfg: Config,\n    candidate_records: list[tuple[str, np.ndarray, Dict[str, Any]]],\n    X_train: np.ndarray,\n    y_train: np.ndarray,\n    X_tune: np.ndarray,\n    y_tune: np.ndarray,\n    seed: int,\n) -> Tuple[str, pd.DataFrame]:\n    records = []\n    best_key = None\n    best_score = -np.inf\n\n    for index, (key, selected_indices, metadata) in enumerate(candidate_records):\n        model, _ = fit_mlp(\n            cfg,\n            X_train[selected_indices],\n            y_train[selected_indices],\n            X_tune,\n            y_tune,\n            seed=seed + 3000 + index,\n            epochs=cfg.tuning_epochs,\n        )\n        probabilities = expit(predict_logits(model, X_tune))\n        threshold, _ = tune_threshold_mcc(\n            y_tune,\n            probabilities,\n            cfg.target_sensitivity,\n        )\n        metrics = compute_metrics(y_tune, probabilities, threshold)\n        score = composite_tuning_score(metrics)\n        record = {\n            \"candidate\": key,\n            **metadata,\n            **metrics,\n            \"composite_score\": score,\n        }\n        records.append(record)\n\n        if score > best_score:\n            best_score = score\n            best_key = key\n\n        del model\n        keras.backend.clear_session()\n        gc.collect()\n\n    if best_key is None:\n        raise RuntimeError(\"No valid resampling candidate was selected.\")\n    return best_key, pd.DataFrame(records)\n\n\ndef indices_from_weighted_resampling(\n    y: np.ndarray,\n    weights: np.ndarray,\n    target_ratio: float,\n    seed: int,\n) -> Tuple[np.ndarray, Dict[str, int]]:\n    index_feature = np.arange(len(y), dtype=np.int64).reshape(-1, 1)\n    resampled_indices, _, details = roof_resample(\n        index_feature,\n        y,\n        weights,\n        target_ratio,\n        seed,\n    )\n    return resampled_indices.reshape(-1).astype(int), details\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.403348Z","iopub.execute_input":"2026-07-03T10:03:17.403606Z","iopub.status.idle":"2026-07-03T10:03:17.43465Z","shell.execute_reply.started":"2026-07-03T10:03:17.403579Z","shell.execute_reply":"2026-07-03T10:03:17.433995Z"}},"outputs":[],"execution_count":null},{"id":"1ed91184","cell_type":"markdown","source":"## 9. Repeated-seed training, checkpointing, and evaluation\n","metadata":{}},{"id":"savefix1","cell_type":"markdown","source":"### Model checkpoint serialization note\n\nFocal-loss and OHEM models are trained with local custom loss closures. For reliable Keras checkpoint saving, the notebook creates an identical BCE-compiled inference model, transfers the trained weights, and saves that copy. Predictions and learned weights are unchanged.\n","metadata":{}},{"id":"7ae2ea0e","cell_type":"code","source":"def group_three_way_split(\n    y: np.ndarray,\n    groups: np.ndarray,\n    seed: int,\n) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:\n    first = StratifiedGroupKFold(\n        n_splits=5,\n        shuffle=True,\n        random_state=seed,\n    )\n    train_idx, temporary_idx = next(\n        first.split(np.zeros(len(y)), y, groups)\n    )\n    temp_y = y[temporary_idx]\n    temp_groups = groups[temporary_idx]\n\n    second = StratifiedGroupKFold(\n        n_splits=2,\n        shuffle=True,\n        random_state=seed + 7919,\n    )\n    tune_rel, calibration_rel = next(\n        second.split(\n            np.zeros(len(temporary_idx)),\n            temp_y,\n            temp_groups,\n        )\n    )\n    return (\n        train_idx,\n        temporary_idx[tune_rel],\n        temporary_idx[calibration_rel],\n    )\n\n\ndef safe_method_name(method: str) -> str:\n    return re.sub(r\"[^A-Za-z0-9_-]+\", \"_\", method)\n\n\ndef experiment_paths(seed: int, method: str) -> Dict[str, Path]:\n    name = safe_method_name(method)\n    return {\n        \"model\": MODEL_DIR / f\"{name}_seed_{seed}.keras\",\n        \"scaler\": MODEL_DIR / f\"{name}_seed_{seed}_scaler.joblib\",\n        \"metadata\": METADATA_DIR / f\"{name}_seed_{seed}.json\",\n        \"internal_predictions\": (\n            PREDICTION_DIR / f\"internal_test_{name}_seed_{seed}.csv\"\n        ),\n        \"external_predictions\": (\n            PREDICTION_DIR / f\"rsna_external_{name}_seed_{seed}.csv\"\n        ),\n        \"history\": METADATA_DIR / f\"{name}_seed_{seed}_history.csv\",\n        \"tuning\": METADATA_DIR / f\"{name}_seed_{seed}_tuning.csv\",\n    }\n\n\ndef experiment_complete(seed: int, method: str) -> bool:\n    paths = experiment_paths(seed, method)\n    required = [\n        paths[\"model\"],\n        paths[\"scaler\"],\n        paths[\"metadata\"],\n        paths[\"internal_predictions\"],\n        paths[\"external_predictions\"],\n    ]\n    return all(path.exists() for path in required)\n\n\ndef save_predictions(\n    manifest: pd.DataFrame,\n    probabilities: np.ndarray,\n    threshold: float,\n    path: Path,\n) -> None:\n    output = manifest.reset_index(drop=True).copy()\n    output[\"probability\"] = probabilities\n    output[\"prediction\"] = (\n        np.asarray(probabilities) >= threshold\n    ).astype(int)\n    output.to_csv(path, index=False)\n\n\ndef build_method_training_data(\n    method: str,\n    cfg: Config,\n    X_train: np.ndarray,\n    y_train: np.ndarray,\n    groups_train: np.ndarray,\n    X_tune: np.ndarray,\n    y_tune: np.ndarray,\n    seed: int,\n    shared_oof: Optional[Tuple[np.ndarray, np.ndarray]],\n) -> Tuple[\n    np.ndarray,\n    np.ndarray,\n    Optional[Dict[int, float]],\n    Any,\n    Dict[str, Any],\n    pd.DataFrame,\n]:\n    class_weight = None\n    loss_function: Any = \"binary_crossentropy\"\n    metadata: Dict[str, Any] = {}\n    tuning_df = pd.DataFrame()\n\n    if method == \"Unbalanced-BCE\":\n        return (\n            X_train,\n            y_train,\n            class_weight,\n            loss_function,\n            metadata,\n            tuning_df,\n        )\n\n    if method == \"ClassWeighted-BCE\":\n        class_weight = class_weight_dict(y_train)\n        return (\n            X_train,\n            y_train,\n            class_weight,\n            loss_function,\n            {\"class_weight\": class_weight},\n            tuning_df,\n        )\n\n    if method == \"Focal-Loss\":\n        loss_function = binary_focal_loss(cfg.focal_gamma)\n        return (\n            X_train,\n            y_train,\n            class_weight,\n            loss_function,\n            {\"gamma\": cfg.focal_gamma},\n            tuning_df,\n        )\n\n    if method == \"OHEM\":\n        loss_function = ohem_binary_crossentropy(cfg.ohem_top_fraction)\n        return (\n            X_train,\n            y_train,\n            class_weight,\n            loss_function,\n            {\"top_fraction\": cfg.ohem_top_fraction},\n            tuning_df,\n        )\n\n    if method == \"Random-Under-Sampling\":\n        sampler = RandomUnderSampler(random_state=seed)\n        X_resampled, y_resampled = sampler.fit_resample(X_train, y_train)\n        return (\n            X_resampled,\n            y_resampled,\n            class_weight,\n            loss_function,\n            {\"sampled_count\": int(len(y_resampled))},\n            tuning_df,\n        )\n\n    if method == \"SMOTE\":\n        minority = int(min((y_train == 0).sum(), (y_train == 1).sum()))\n        sampler = SMOTE(\n            random_state=seed,\n            k_neighbors=max(1, min(5, minority - 1)),\n        )\n        X_resampled, y_resampled = sampler.fit_resample(X_train, y_train)\n        return (\n            X_resampled,\n            y_resampled,\n            class_weight,\n            loss_function,\n            {\n                \"sampled_count\": int(len(y_resampled)),\n                \"synthetic_feature_space_sampling\": True,\n            },\n            tuning_df,\n        )\n\n    if method == \"DARS-Original\":\n        in_sample_probability = generate_in_sample_scout_predictions(\n            cfg,\n            X_train,\n            y_train,\n            X_tune,\n            y_tune,\n            seed,\n        )\n        loss = class_quantile_normalize(\n            bce_numpy(y_train, in_sample_probability),\n            y_train,\n        )\n        misclassified = (\n            (in_sample_probability >= 0.5).astype(int) != y_train\n        ).astype(float)\n        uncertainty = np.clip(\n            1.0 - 2.0 * np.abs(in_sample_probability - 0.5),\n            0.0,\n            1.0,\n        )\n        weights = np.clip(\n            1.0 + 2.0 * loss + 3.0 * misclassified + uncertainty,\n            0.5,\n            cfg.weight_cap,\n        )\n        indices, details = indices_from_weighted_resampling(\n            y_train,\n            weights,\n            target_ratio=1.0,\n            seed=seed + 900,\n        )\n        metadata = {\n            \"alpha_loss\": 2.0,\n            \"beta_misclassified\": 3.0,\n            \"gamma_uncertainty\": 1.0,\n            \"target_ratio\": 1.0,\n            **details,\n        }\n        return (\n            X_train[indices],\n            y_train[indices],\n            class_weight,\n            loss_function,\n            metadata,\n            tuning_df,\n        )\n\n    if shared_oof is None:\n        raise RuntimeError(f\"OOF predictions are required for {method}.\")\n\n    oof_probability, disagreement = shared_oof\n\n    if method == \"OOF-Loss\":\n        loss_norm = class_quantile_normalize(\n            bce_numpy(y_train, oof_probability),\n            y_train,\n        )\n        weights = np.clip(\n            1.0 + 1.5 * loss_norm,\n            0.5,\n            cfg.weight_cap,\n        )\n        candidates = []\n        candidate_indices = {}\n        for ratio in [0.35, 0.50, 0.65, 0.80]:\n            key = f\"ratio_{ratio:.2f}\"\n            indices, details = indices_from_weighted_resampling(\n                y_train,\n                weights,\n                target_ratio=ratio,\n                seed=seed + int(1000 * ratio),\n            )\n            candidate_indices[key] = indices\n            candidates.append(\n                (\n                    key,\n                    indices,\n                    {\n                        \"alpha_loss\": 1.5,\n                        \"target_ratio\": ratio,\n                        **details,\n                    },\n                )\n            )\n        best_key, tuning_df = tune_resampling_candidate(\n            cfg,\n            candidates,\n            X_train,\n            y_train,\n            X_tune,\n            y_tune,\n            seed,\n        )\n        selected = candidate_indices[best_key]\n        best_record = (\n            tuning_df.loc[tuning_df[\"candidate\"] == best_key]\n            .iloc[0]\n            .to_dict()\n        )\n        return (\n            X_train[selected],\n            y_train[selected],\n            class_weight,\n            loss_function,\n            best_record,\n            tuning_df,\n        )\n\n    if method == \"ROOF-DARS\":\n        candidates = []\n        candidate_indices = {}\n        candidate_details = {}\n\n        for candidate_index, params in enumerate(roof_parameter_grid(cfg)):\n            weights, difficulty_details = roof_difficulty(\n                cfg,\n                y_train,\n                oof_probability,\n                disagreement,\n                params,\n            )\n            indices, sampling_details = indices_from_weighted_resampling(\n                y_train,\n                weights,\n                target_ratio=params.target_ratio,\n                seed=seed + 2000 + candidate_index,\n            )\n            key = f\"candidate_{candidate_index:02d}\"\n            candidate_indices[key] = indices\n            candidate_details[key] = {\n                **asdict(params),\n                **difficulty_details,\n                **sampling_details,\n            }\n            candidates.append(\n                (\n                    key,\n                    indices,\n                    candidate_details[key],\n                )\n            )\n\n        best_key, tuning_df = tune_resampling_candidate(\n            cfg,\n            candidates,\n            X_train,\n            y_train,\n            X_tune,\n            y_tune,\n            seed,\n        )\n        selected = candidate_indices[best_key]\n        best_record = (\n            tuning_df.loc[tuning_df[\"candidate\"] == best_key]\n            .iloc[0]\n            .to_dict()\n        )\n        return (\n            X_train[selected],\n            y_train[selected],\n            class_weight,\n            loss_function,\n            best_record,\n            tuning_df,\n        )\n\n    raise ValueError(f\"Unknown method: {method}\")\n\n\ndef train_and_evaluate_method(\n    method: str,\n    seed: int,\n    dev_df: pd.DataFrame,\n    internal_df: pd.DataFrame,\n    rsna_df: pd.DataFrame,\n    X_dev_raw: np.ndarray,\n    X_internal_raw: np.ndarray,\n    X_rsna_raw: np.ndarray,\n    shared_oof: Optional[Tuple[np.ndarray, np.ndarray]],\n) -> list[Dict[str, Any]]:\n    paths = experiment_paths(seed, method)\n    if experiment_complete(seed, method):\n        print(f\"Skipping completed experiment: seed={seed}, method={method}\")\n        metadata = json.loads(paths[\"metadata\"].read_text(encoding=\"utf-8\"))\n        return metadata[\"metric_records\"]\n\n    y_dev = dev_df[\"label\"].to_numpy(dtype=int)\n    groups = dev_df[\"patient_id\"].to_numpy()\n    train_idx, tune_idx, calibration_idx = group_three_way_split(\n        y_dev,\n        groups,\n        seed,\n    )\n\n    scaler = StandardScaler()\n    X_train = scaler.fit_transform(X_dev_raw[train_idx])\n    X_tune = scaler.transform(X_dev_raw[tune_idx])\n    X_calibration = scaler.transform(X_dev_raw[calibration_idx])\n    X_internal = scaler.transform(X_internal_raw)\n    X_rsna = scaler.transform(X_rsna_raw)\n\n    y_train = y_dev[train_idx]\n    y_tune = y_dev[tune_idx]\n    y_calibration = y_dev[calibration_idx]\n    groups_train = groups[train_idx]\n\n    method_oof = shared_oof\n    if method_oof is not None:\n        # shared_oof was computed in the seed-specific standardized feature space.\n        pass\n\n    (\n        X_method,\n        y_method,\n        class_weight,\n        loss_function,\n        method_metadata,\n        tuning_df,\n    ) = build_method_training_data(\n        method,\n        cfg,\n        X_train,\n        y_train,\n        groups_train,\n        X_tune,\n        y_tune,\n        seed,\n        method_oof,\n    )\n\n    print(\n        f\"Training seed={seed}, method={method}, \"\n        f\"samples={len(y_method)}\"\n    )\n    start = time.perf_counter()\n    model, history = fit_mlp(\n        cfg,\n        X_method,\n        y_method,\n        X_tune,\n        y_tune,\n        seed=seed + 50000 + METHODS.index(method),\n        epochs=cfg.mlp_epochs,\n        class_weight=class_weight,\n        loss_function=loss_function,\n    )\n    training_seconds = time.perf_counter() - start\n\n    calibration_logits = predict_logits(model, X_calibration)\n    slope, intercept = fit_affine_calibrator(\n        y_calibration,\n        calibration_logits,\n    )\n    calibration_probabilities = calibrated_probabilities(\n        calibration_logits,\n        slope,\n        intercept,\n    )\n    threshold, threshold_details = tune_threshold_mcc(\n        y_calibration,\n        calibration_probabilities,\n        cfg.target_sensitivity,\n    )\n\n    metric_records = []\n    cohort_inputs = [\n        (\n            \"internal_test\",\n            internal_df,\n            X_internal,\n            paths[\"internal_predictions\"],\n        ),\n        (\n            \"rsna_external\",\n            rsna_df,\n            X_rsna,\n            paths[\"external_predictions\"],\n        ),\n    ]\n\n    for cohort_name, cohort_df, X_cohort, prediction_path in cohort_inputs:\n        inference_start = time.perf_counter()\n        logits = predict_logits(model, X_cohort)\n        inference_seconds = time.perf_counter() - inference_start\n        probabilities = calibrated_probabilities(logits, slope, intercept)\n        metrics = compute_metrics(\n            cohort_df[\"label\"].to_numpy(dtype=int),\n            probabilities,\n            threshold,\n        )\n        save_predictions(\n            cohort_df,\n            probabilities,\n            threshold,\n            prediction_path,\n        )\n        metric_records.append(\n            {\n                \"seed\": seed,\n                \"method\": method,\n                \"cohort\": cohort_name,\n                \"training_seconds\": training_seconds,\n                \"inference_seconds\": inference_seconds,\n                \"inference_ms_per_image\": (\n                    1000.0 * inference_seconds / len(cohort_df)\n                ),\n                \"trainable_parameters\": int(model.count_params()),\n                \"training_sample_count\": int(len(y_method)),\n                \"calibration_slope\": slope,\n                \"calibration_intercept\": intercept,\n                **metrics,\n            }\n        )\n\n    # Save an inference-compatible copy using a standard serializable loss.\n    # The architecture and trained weights remain unchanged; only compile\n    # metadata is replaced so focal/OHEM closure losses do not break saving.\n    serializable_model = build_mlp(\n        cfg,\n        input_dim=int(model.input_shape[-1]),\n        loss_function=\"binary_crossentropy\",\n    )\n    serializable_model.set_weights(model.get_weights())\n    serializable_model.save(paths[\"model\"])\n    del serializable_model\n\n    joblib.dump(scaler, paths[\"scaler\"])\n    history.to_csv(paths[\"history\"], index=False)\n    if len(tuning_df):\n        tuning_df.to_csv(paths[\"tuning\"], index=False)\n\n    metadata = {\n        \"seed\": seed,\n        \"method\": method,\n        \"split\": {\n            \"train_indices\": train_idx.tolist(),\n            \"tune_indices\": tune_idx.tolist(),\n            \"calibration_indices\": calibration_idx.tolist(),\n        },\n        \"method_metadata\": method_metadata,\n        \"calibration\": {\n            \"slope\": slope,\n            \"intercept\": intercept,\n        },\n        \"threshold\": threshold,\n        \"threshold_details\": threshold_details,\n        \"metric_records\": metric_records,\n        \"config\": asdict(cfg),\n    }\n    paths[\"metadata\"].write_text(\n        json.dumps(metadata, indent=2, default=float),\n        encoding=\"utf-8\",\n    )\n\n    del model\n    keras.backend.clear_session()\n    gc.collect()\n    return metric_records\n\n\ndef run_repeated_experiments(\n    dev_df: pd.DataFrame,\n    internal_df: pd.DataFrame,\n    rsna_df: pd.DataFrame,\n    X_dev_raw: np.ndarray,\n    X_internal_raw: np.ndarray,\n    X_rsna_raw: np.ndarray,\n) -> pd.DataFrame:\n    all_records = []\n\n    for seed in cfg.seeds:\n        print(\"\\n\" + \"=\" * 80)\n        print(\"SEED\", seed)\n        print(\"=\" * 80)\n        set_all_seeds(seed)\n\n        y_dev = dev_df[\"label\"].to_numpy(dtype=int)\n        groups = dev_df[\"patient_id\"].to_numpy()\n        train_idx, tune_idx, _ = group_three_way_split(y_dev, groups, seed)\n\n        scaler_for_scout = StandardScaler()\n        X_train_for_scout = scaler_for_scout.fit_transform(\n            X_dev_raw[train_idx]\n        )\n        X_tune_for_scout = scaler_for_scout.transform(\n            X_dev_raw[tune_idx]\n        )\n        y_train_for_scout = y_dev[train_idx]\n        groups_train = groups[train_idx]\n\n        needs_oof = any(\n            not experiment_complete(seed, method)\n            for method in [\"OOF-Loss\", \"ROOF-DARS\"]\n        )\n        shared_oof = None\n        if needs_oof:\n            shared_oof = generate_oof_scout_predictions(\n                cfg,\n                X_train_for_scout,\n                y_train_for_scout,\n                groups_train,\n                seed,\n            )\n\n        for method in METHODS:\n            try:\n                records = train_and_evaluate_method(\n                    method,\n                    seed,\n                    dev_df,\n                    internal_df,\n                    rsna_df,\n                    X_dev_raw,\n                    X_internal_raw,\n                    X_rsna_raw,\n                    shared_oof,\n                )\n                all_records.extend(records)\n            except Exception as error:\n                failure = {\n                    \"seed\": seed,\n                    \"method\": method,\n                    \"error\": repr(error),\n                }\n                failure_path = (\n                    METADATA_DIR\n                    / f\"FAILED_{safe_method_name(method)}_seed_{seed}.json\"\n                )\n                failure_path.write_text(\n                    json.dumps(failure, indent=2),\n                    encoding=\"utf-8\",\n                )\n                print(\"FAILED:\", failure)\n                raise\n\n        pd.DataFrame(all_records).to_csv(\n            TABLE_DIR / \"per_seed_results_partial.csv\",\n            index=False,\n        )\n\n    results = pd.DataFrame(all_records)\n    results.to_csv(TABLE_DIR / \"per_seed_results.csv\", index=False)\n    return results\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.436653Z","iopub.execute_input":"2026-07-03T10:03:17.43702Z","iopub.status.idle":"2026-07-03T10:03:17.47687Z","shell.execute_reply.started":"2026-07-03T10:03:17.436999Z","shell.execute_reply":"2026-07-03T10:03:17.476029Z"}},"outputs":[],"execution_count":null},{"id":"3732435f","cell_type":"markdown","source":"## 10. Statistical analysis and manuscript tables\n","metadata":{}},{"id":"262cb202","cell_type":"code","source":"def summarize_across_seeds(results: pd.DataFrame) -> pd.DataFrame:\n    metrics = [\n        \"accuracy\",\n        \"balanced_accuracy\",\n        \"sensitivity\",\n        \"specificity\",\n        \"precision\",\n        \"f1\",\n        \"mcc\",\n        \"roc_auc\",\n        \"pr_auc\",\n        \"brier\",\n        \"ece\",\n        \"training_seconds\",\n        \"inference_ms_per_image\",\n    ]\n    rows = []\n\n    for (cohort, method), group in results.groupby([\"cohort\", \"method\"]):\n        for metric in metrics:\n            values = group[metric].dropna().to_numpy(dtype=float)\n            n = len(values)\n            mean = float(values.mean())\n            sd = float(values.std(ddof=1)) if n > 1 else 0.0\n            half_width = (\n                float(t.ppf(0.975, df=n - 1) * sd / math.sqrt(n))\n                if n > 1\n                else 0.0\n            )\n            rows.append(\n                {\n                    \"cohort\": cohort,\n                    \"method\": method,\n                    \"metric\": metric,\n                    \"n_seeds\": n,\n                    \"mean\": mean,\n                    \"sd\": sd,\n                    \"ci_lower\": mean - half_width,\n                    \"ci_upper\": mean + half_width,\n                }\n            )\n    return pd.DataFrame(rows)\n\n\ndef build_manuscript_wide_table(summary: pd.DataFrame) -> pd.DataFrame:\n    selected_metrics = [\n        \"accuracy\",\n        \"balanced_accuracy\",\n        \"sensitivity\",\n        \"specificity\",\n        \"f1\",\n        \"mcc\",\n        \"roc_auc\",\n        \"pr_auc\",\n        \"brier\",\n        \"ece\",\n    ]\n    selected = summary[summary[\"metric\"].isin(selected_metrics)].copy()\n    selected[\"formatted\"] = selected.apply(\n        lambda row: f\"{row['mean']:.4f} ± {row['sd']:.4f}\",\n        axis=1,\n    )\n    return selected.pivot_table(\n        index=[\"cohort\", \"method\"],\n        columns=\"metric\",\n        values=\"formatted\",\n        aggfunc=\"first\",\n    ).reset_index()\n\n\ndef load_seed_predictions(\n    cohort: str,\n    method: str,\n) -> Tuple[pd.DataFrame, np.ndarray, np.ndarray]:\n    frames = []\n    thresholds = []\n    for seed in cfg.seeds:\n        paths = experiment_paths(seed, method)\n        prediction_path = (\n            paths[\"internal_predictions\"]\n            if cohort == \"internal_test\"\n            else paths[\"external_predictions\"]\n        )\n        frame = pd.read_csv(prediction_path)\n        metadata = json.loads(\n            paths[\"metadata\"].read_text(encoding=\"utf-8\")\n        )\n        frames.append(frame)\n        thresholds.append(float(metadata[\"threshold\"]))\n\n    reference = frames[0].copy()\n    reference_paths = reference[\"path\"].astype(str).to_numpy()\n    matrix = []\n\n    for frame in frames:\n        if not np.array_equal(\n            frame[\"path\"].astype(str).to_numpy(),\n            reference_paths,\n        ):\n            raise RuntimeError(\n                f\"Prediction order mismatch for {cohort}, {method}.\"\n            )\n        matrix.append(frame[\"probability\"].to_numpy(dtype=float))\n\n    mean_probability = np.vstack(matrix).mean(axis=0)\n    mean_threshold = float(np.mean(thresholds))\n    return reference, mean_probability, np.asarray(thresholds)\n\n\ndef build_ensemble_results() -> Tuple[pd.DataFrame, Dict[str, Dict[str, Any]]]:\n    records = []\n    store: Dict[str, Dict[str, Any]] = {}\n\n    for cohort in [\"internal_test\", \"rsna_external\"]:\n        store[cohort] = {}\n        for method in METHODS:\n            frame, probabilities, thresholds = load_seed_predictions(\n                cohort,\n                method,\n            )\n            threshold = float(thresholds.mean())\n            y = frame[\"label\"].to_numpy(dtype=int)\n            metrics = compute_metrics(y, probabilities, threshold)\n            predictions = (probabilities >= threshold).astype(int)\n\n            output = frame.copy()\n            output[\"ensemble_probability\"] = probabilities\n            output[\"ensemble_prediction\"] = predictions\n            output.to_csv(\n                PREDICTION_DIR\n                / f\"{cohort}_{safe_method_name(method)}_ensemble.csv\",\n                index=False,\n            )\n\n            records.append(\n                {\n                    \"cohort\": cohort,\n                    \"method\": method,\n                    \"ensemble_threshold\": threshold,\n                    **metrics,\n                }\n            )\n            store[cohort][method] = {\n                \"frame\": frame,\n                \"y\": y,\n                \"probabilities\": probabilities,\n                \"predictions\": predictions,\n                \"threshold\": threshold,\n            }\n\n    return pd.DataFrame(records), store\n\n\ndef compute_midrank(x: np.ndarray) -> np.ndarray:\n    order = np.argsort(x)\n    sorted_x = x[order]\n    result = np.zeros(len(x), dtype=float)\n    index = 0\n    while index < len(x):\n        end = index\n        while end < len(x) and sorted_x[end] == sorted_x[index]:\n            end += 1\n        result[index:end] = 0.5 * (index + end - 1) + 1\n        index = end\n    output = np.empty(len(x), dtype=float)\n    output[order] = result\n    return output\n\n\ndef fast_delong(\n    predictions_sorted_transposed: np.ndarray,\n    label_1_count: int,\n) -> Tuple[np.ndarray, np.ndarray]:\n    m = label_1_count\n    n = predictions_sorted_transposed.shape[1] - m\n    positive = predictions_sorted_transposed[:, :m]\n    negative = predictions_sorted_transposed[:, m:]\n    k = predictions_sorted_transposed.shape[0]\n\n    tx = np.empty((k, m))\n    ty = np.empty((k, n))\n    tz = np.empty((k, m + n))\n    for classifier in range(k):\n        tx[classifier] = compute_midrank(positive[classifier])\n        ty[classifier] = compute_midrank(negative[classifier])\n        tz[classifier] = compute_midrank(\n            predictions_sorted_transposed[classifier]\n        )\n\n    aucs = (\n        tz[:, :m].sum(axis=1) / m / n\n        - (m + 1.0) / 2.0 / n\n    )\n    v01 = (tz[:, :m] - tx) / n\n    v10 = 1.0 - (tz[:, m:] - ty) / m\n    sx = np.cov(v01)\n    sy = np.cov(v10)\n    covariance = sx / m + sy / n\n    return aucs, np.atleast_2d(covariance)\n\n\ndef delong_roc_test(\n    y_true: np.ndarray,\n    probabilities_a: np.ndarray,\n    probabilities_b: np.ndarray,\n) -> Tuple[float, float, float, float]:\n    y = np.asarray(y_true).astype(int).ravel()\n    a = np.asarray(probabilities_a, dtype=float).ravel()\n    b = np.asarray(probabilities_b, dtype=float).ravel()\n    order = np.argsort(-y)\n    predictions = np.vstack([a, b])[:, order]\n    aucs, covariance = fast_delong(\n        predictions,\n        int(y.sum()),\n    )\n    contrast = np.array([[1.0, -1.0]])\n    variance_matrix = contrast @ covariance @ contrast.T\n    variance = float(np.asarray(variance_matrix).squeeze().item())\n    z = abs(float(aucs[0] - aucs[1])) / math.sqrt(\n        max(variance, 1e-12)\n    )\n    p_value = 2.0 * (1.0 - norm.cdf(z))\n    return float(aucs[0]), float(aucs[1]), float(z), float(p_value)\n\n\ndef paired_bootstrap_difference(\n    y_true: np.ndarray,\n    probabilities_reference: np.ndarray,\n    probabilities_competitor: np.ndarray,\n    threshold_reference: float,\n    threshold_competitor: float,\n    metric: str,\n    iterations: int,\n    seed: int,\n) -> Tuple[float, float, float, float]:\n    y = np.asarray(y_true).astype(int)\n    ref = np.asarray(probabilities_reference)\n    comp = np.asarray(probabilities_competitor)\n    rng = np.random.default_rng(seed)\n    class_indices = {\n        label: np.where(y == label)[0]\n        for label in [0, 1]\n    }\n\n    def metric_value(\n        labels: np.ndarray,\n        probabilities: np.ndarray,\n        threshold: float,\n    ) -> float:\n        if metric == \"pr_auc\":\n            return float(average_precision_score(labels, probabilities))\n        if metric == \"mcc\":\n            predictions = (probabilities >= threshold).astype(int)\n            return float(matthews_corrcoef(labels, predictions))\n        if metric == \"balanced_accuracy\":\n            predictions = (probabilities >= threshold).astype(int)\n            return float(\n                balanced_accuracy_score(labels, predictions)\n            )\n        raise ValueError(metric)\n\n    observed = (\n        metric_value(y, ref, threshold_reference)\n        - metric_value(y, comp, threshold_competitor)\n    )\n    differences = []\n\n    for _ in range(iterations):\n        sampled = np.concatenate(\n            [\n                rng.choice(indices, size=len(indices), replace=True)\n                for indices in class_indices.values()\n            ]\n        )\n        rng.shuffle(sampled)\n        differences.append(\n            metric_value(\n                y[sampled],\n                ref[sampled],\n                threshold_reference,\n            )\n            - metric_value(\n                y[sampled],\n                comp[sampled],\n                threshold_competitor,\n            )\n        )\n\n    low, high = np.percentile(differences, [2.5, 97.5])\n    p_value = 2.0 * min(\n        np.mean(np.asarray(differences) <= 0.0),\n        np.mean(np.asarray(differences) >= 0.0),\n    )\n    return float(observed), float(low), float(high), float(min(1.0, p_value))\n\n\ndef pairwise_statistics(\n    ensemble_store: Dict[str, Dict[str, Any]],\n    cohort: str = \"rsna_external\",\n) -> pd.DataFrame:\n    reference = ensemble_store[cohort][REFERENCE_METHOD]\n    y = reference[\"y\"]\n    records = []\n\n    for competitor in METHODS:\n        if competitor == REFERENCE_METHOD:\n            continue\n        other = ensemble_store[cohort][competitor]\n\n        b = int(\n            np.sum(\n                (reference[\"predictions\"] == y)\n                & (other[\"predictions\"] != y)\n            )\n        )\n        c = int(\n            np.sum(\n                (reference[\"predictions\"] != y)\n                & (other[\"predictions\"] == y)\n            )\n        )\n        discordant = b + c\n        mcnemar_p = (\n            float(\n                binomtest(\n                    min(b, c),\n                    discordant,\n                    p=0.5,\n                    alternative=\"two-sided\",\n                ).pvalue\n            )\n            if discordant\n            else 1.0\n        )\n\n        _, _, delong_z, delong_p = delong_roc_test(\n            y,\n            reference[\"probabilities\"],\n            other[\"probabilities\"],\n        )\n        pr_diff, pr_low, pr_high, pr_p = paired_bootstrap_difference(\n            y,\n            reference[\"probabilities\"],\n            other[\"probabilities\"],\n            reference[\"threshold\"],\n            other[\"threshold\"],\n            metric=\"pr_auc\",\n            iterations=cfg.bootstrap_iterations,\n            seed=7000 + METHODS.index(competitor),\n        )\n        mcc_diff, mcc_low, mcc_high, mcc_p = paired_bootstrap_difference(\n            y,\n            reference[\"probabilities\"],\n            other[\"probabilities\"],\n            reference[\"threshold\"],\n            other[\"threshold\"],\n            metric=\"mcc\",\n            iterations=cfg.bootstrap_iterations,\n            seed=9000 + METHODS.index(competitor),\n        )\n\n        records.append(\n            {\n                \"cohort\": cohort,\n                \"reference\": REFERENCE_METHOD,\n                \"competitor\": competitor,\n                \"reference_correct_competitor_wrong\": b,\n                \"reference_wrong_competitor_correct\": c,\n                \"mcnemar_p\": mcnemar_p,\n                \"delong_z\": delong_z,\n                \"delong_p\": delong_p,\n                \"pr_auc_difference\": pr_diff,\n                \"pr_auc_ci_lower\": pr_low,\n                \"pr_auc_ci_upper\": pr_high,\n                \"pr_auc_p\": pr_p,\n                \"mcc_difference\": mcc_diff,\n                \"mcc_ci_lower\": mcc_low,\n                \"mcc_ci_upper\": mcc_high,\n                \"mcc_p\": mcc_p,\n            }\n        )\n\n    result = pd.DataFrame(records)\n    for column in [\n        \"mcnemar_p\",\n        \"delong_p\",\n        \"pr_auc_p\",\n        \"mcc_p\",\n    ]:\n        result[f\"{column}_holm\"] = multipletests(\n            result[column].fillna(1.0),\n            method=\"holm\",\n        )[1]\n    return result\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.477889Z","iopub.execute_input":"2026-07-03T10:03:17.478155Z","iopub.status.idle":"2026-07-03T10:03:17.510422Z","shell.execute_reply.started":"2026-07-03T10:03:17.478136Z","shell.execute_reply":"2026-07-03T10:03:17.509704Z"}},"outputs":[],"execution_count":null},{"id":"2e24f238","cell_type":"markdown","source":"## 11. Manuscript figures\n","metadata":{}},{"id":"bce6a523","cell_type":"code","source":"def plot_metric_comparison(\n    summary: pd.DataFrame,\n    cohort: str,\n    metric: str,\n) -> None:\n    selected = summary[\n        (summary[\"cohort\"] == cohort)\n        & (summary[\"metric\"] == metric)\n    ].copy()\n    selected = selected.sort_values(\"mean\", ascending=False)\n\n    plt.figure(figsize=(10, 5.5))\n    plt.bar(\n        selected[\"method\"],\n        selected[\"mean\"],\n        yerr=selected[\"sd\"],\n        capsize=4,\n    )\n    plt.ylabel(metric.replace(\"_\", \" \").title())\n    plt.xlabel(\"Method\")\n    plt.title(f\"{cohort}: {metric.replace('_', ' ').title()} across seeds\")\n    plt.xticks(rotation=40, ha=\"right\")\n    plt.tight_layout()\n    plt.savefig(\n        FIGURE_DIR / f\"{cohort}_{metric}_comparison.png\",\n        dpi=300,\n    )\n    plt.show()\n\n\ndef plot_ensemble_curves(\n    ensemble_store: Dict[str, Dict[str, Any]],\n    cohort: str,\n) -> None:\n    plt.figure(figsize=(9, 7))\n    for method in METHODS:\n        data = ensemble_store[cohort][method]\n        fpr, tpr, _ = roc_curve(\n            data[\"y\"],\n            data[\"probabilities\"],\n        )\n        auc = roc_auc_score(data[\"y\"], data[\"probabilities\"])\n        plt.plot(fpr, tpr, label=f\"{method}: {auc:.3f}\")\n    plt.plot([0, 1], [0, 1], \"--\")\n    plt.xlabel(\"False-positive rate\")\n    plt.ylabel(\"True-positive rate\")\n    plt.title(f\"{cohort}: ensemble ROC curves\")\n    plt.legend(fontsize=7)\n    plt.tight_layout()\n    plt.savefig(FIGURE_DIR / f\"{cohort}_ensemble_roc.png\", dpi=300)\n    plt.show()\n\n    plt.figure(figsize=(9, 7))\n    prevalence = ensemble_store[cohort][REFERENCE_METHOD][\"y\"].mean()\n    for method in METHODS:\n        data = ensemble_store[cohort][method]\n        precision, recall, _ = precision_recall_curve(\n            data[\"y\"],\n            data[\"probabilities\"],\n        )\n        auc = average_precision_score(\n            data[\"y\"],\n            data[\"probabilities\"],\n        )\n        plt.plot(recall, precision, label=f\"{method}: {auc:.3f}\")\n    plt.axhline(prevalence, linestyle=\"--\", label=f\"Prevalence: {prevalence:.3f}\")\n    plt.xlabel(\"Recall\")\n    plt.ylabel(\"Precision\")\n    plt.title(f\"{cohort}: ensemble precision-recall curves\")\n    plt.legend(fontsize=7)\n    plt.tight_layout()\n    plt.savefig(FIGURE_DIR / f\"{cohort}_ensemble_pr.png\", dpi=300)\n    plt.show()\n\n\ndef plot_reference_calibration(\n    ensemble_store: Dict[str, Dict[str, Any]],\n    cohort: str,\n) -> None:\n    data = ensemble_store[cohort][REFERENCE_METHOD]\n    fraction_positive, mean_predicted = calibration_curve(\n        data[\"y\"],\n        data[\"probabilities\"],\n        n_bins=12,\n        strategy=\"quantile\",\n    )\n    plt.figure(figsize=(6, 5.5))\n    plt.plot(\n        mean_predicted,\n        fraction_positive,\n        marker=\"o\",\n        label=REFERENCE_METHOD,\n    )\n    plt.plot([0, 1], [0, 1], \"--\", label=\"Perfect calibration\")\n    plt.xlabel(\"Mean predicted probability\")\n    plt.ylabel(\"Observed positive fraction\")\n    plt.title(f\"{cohort}: {REFERENCE_METHOD} calibration\")\n    plt.legend()\n    plt.tight_layout()\n    plt.savefig(\n        FIGURE_DIR / f\"{cohort}_{safe_method_name(REFERENCE_METHOD)}_calibration.png\",\n        dpi=300,\n    )\n    plt.show()\n\n\ndef generalization_gap_table(\n    ensemble_results: pd.DataFrame,\n) -> pd.DataFrame:\n    metrics = [\n        \"balanced_accuracy\",\n        \"sensitivity\",\n        \"specificity\",\n        \"f1\",\n        \"mcc\",\n        \"roc_auc\",\n        \"pr_auc\",\n        \"brier\",\n        \"ece\",\n    ]\n    rows = []\n    for method in METHODS:\n        internal = ensemble_results[\n            (ensemble_results[\"cohort\"] == \"internal_test\")\n            & (ensemble_results[\"method\"] == method)\n        ].iloc[0]\n        external = ensemble_results[\n            (ensemble_results[\"cohort\"] == \"rsna_external\")\n            & (ensemble_results[\"method\"] == method)\n        ].iloc[0]\n        for metric in metrics:\n            rows.append(\n                {\n                    \"method\": method,\n                    \"metric\": metric,\n                    \"internal_test\": internal[metric],\n                    \"rsna_external\": external[metric],\n                    \"external_minus_internal\": (\n                        external[metric] - internal[metric]\n                    ),\n                }\n            )\n    return pd.DataFrame(rows)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.511276Z","iopub.execute_input":"2026-07-03T10:03:17.511633Z","iopub.status.idle":"2026-07-03T10:03:17.527469Z","shell.execute_reply.started":"2026-07-03T10:03:17.511603Z","shell.execute_reply":"2026-07-03T10:03:17.526819Z"}},"outputs":[],"execution_count":null},{"id":"fdfeb7a1","cell_type":"markdown","source":"## 12. Load the complete cohorts\n","metadata":{}},{"id":"9c620c9d","cell_type":"code","source":"dev_df, internal_test_df = load_kermany(cfg)\nrsna_df = load_rsna_external(cfg)\n\ndev_df.to_csv(OUTPUT_DIR / \"kermany_development_manifest.csv\", index=False)\ninternal_test_df.to_csv(\n    OUTPUT_DIR / \"kermany_internal_test_manifest.csv\",\n    index=False,\n)\nrsna_df.to_csv(OUTPUT_DIR / \"rsna_external_manifest.csv\", index=False)\n\ncohort_summary = pd.DataFrame(\n    [\n        {\n            \"cohort\": \"Kermany development\",\n            \"images\": len(dev_df),\n            \"patients_or_groups\": dev_df[\"patient_id\"].nunique(),\n            \"negative\": int((dev_df.label == 0).sum()),\n            \"positive\": int((dev_df.label == 1).sum()),\n            \"positive_prevalence\": float(dev_df.label.mean()),\n        },\n        {\n            \"cohort\": \"Kermany internal test\",\n            \"images\": len(internal_test_df),\n            \"patients_or_groups\": internal_test_df[\"patient_id\"].nunique(),\n            \"negative\": int((internal_test_df.label == 0).sum()),\n            \"positive\": int((internal_test_df.label == 1).sum()),\n            \"positive_prevalence\": float(internal_test_df.label.mean()),\n        },\n        {\n            \"cohort\": \"RSNA external\",\n            \"images\": len(rsna_df),\n            \"patients_or_groups\": rsna_df[\"patient_id\"].nunique(),\n            \"negative\": int((rsna_df.label == 0).sum()),\n            \"positive\": int((rsna_df.label == 1).sum()),\n            \"positive_prevalence\": float(rsna_df.label.mean()),\n        },\n    ]\n)\ncohort_summary.to_csv(TABLE_DIR / \"cohort_summary.csv\", index=False)\ndisplay(cohort_summary)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:03:17.528273Z","iopub.execute_input":"2026-07-03T10:03:17.528519Z","iopub.status.idle":"2026-07-03T10:04:22.056259Z","shell.execute_reply.started":"2026-07-03T10:03:17.528494Z","shell.execute_reply":"2026-07-03T10:04:22.055508Z"}},"outputs":[],"execution_count":null},{"id":"d42c54a6","cell_type":"markdown","source":"## 13. Run the mandatory duplicate audit\n","metadata":{}},{"id":"e26dfdb4","cell_type":"code","source":"if cfg.audit_exact_duplicates:\n    print(\"Running standardized-pixel exact duplicate audits...\")\n\n    internal_duplicates = exact_duplicate_audit(\n        dev_df,\n        internal_test_df,\n        OUTPUT_DIR / \"dev_vs_internal_exact_duplicates.csv\",\n    )\n    rsna_duplicates = exact_duplicate_audit(\n        pd.concat([dev_df, internal_test_df], ignore_index=True),\n        rsna_df,\n        OUTPUT_DIR / \"kermany_vs_rsna_exact_duplicates.csv\",\n    )\n\n    print(\"Development/internal duplicates:\", len(internal_duplicates))\n    print(\"Kermany/RSNA duplicates:\", len(rsna_duplicates))\n\n    if len(internal_duplicates):\n        raise RuntimeError(\n            \"Exact duplicates were detected between development and internal testing.\"\n        )\n\n    if len(rsna_duplicates):\n        duplicate_paths = set(rsna_duplicates[\"candidate_path\"])\n        rsna_df = rsna_df[\n            ~rsna_df[\"path\"].isin(duplicate_paths)\n        ].reset_index(drop=True)\n        rsna_df.to_csv(\n            OUTPUT_DIR / \"rsna_external_manifest_after_duplicate_removal.csv\",\n            index=False,\n        )\n        print(\"RSNA after duplicate removal:\", len(rsna_df))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:04:22.057194Z","iopub.execute_input":"2026-07-03T10:04:22.057559Z","iopub.status.idle":"2026-07-03T10:25:17.344762Z","shell.execute_reply.started":"2026-07-03T10:04:22.057537Z","shell.execute_reply":"2026-07-03T10:25:17.344022Z"}},"outputs":[],"execution_count":null},{"id":"de52af77","cell_type":"markdown","source":"## 14. Extract or load cached DenseNet121 features\n","metadata":{}},{"id":"e15d700e","cell_type":"code","source":"extractor = XRVFeatureExtractor(cfg)\n\nX_dev_raw = extractor.extract(\n    dev_df[\"path\"].tolist(),\n    dev_df[\"label\"].tolist(),\n    CACHE_DIR,\n    \"kermany_development_full\",\n)\nX_internal_raw = extractor.extract(\n    internal_test_df[\"path\"].tolist(),\n    internal_test_df[\"label\"].tolist(),\n    CACHE_DIR,\n    \"kermany_internal_test_full\",\n)\nX_rsna_raw = extractor.extract(\n    rsna_df[\"path\"].tolist(),\n    rsna_df[\"label\"].tolist(),\n    CACHE_DIR,\n    \"rsna_external_full\",\n)\n\nprint(\"Feature shapes:\")\nprint(\" Development:\", X_dev_raw.shape)\nprint(\" Internal test:\", X_internal_raw.shape)\nprint(\" RSNA external:\", X_rsna_raw.shape)\n\n# DenseNet features are cached. Move the PyTorch backbone to CPU so TensorFlow\n# can use the GPU for repeated MLP training without unnecessary memory pressure.\nextractor.model.to(\"cpu\")\ntorch.cuda.empty_cache()\ngc.collect()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:25:17.34591Z","iopub.execute_input":"2026-07-03T10:25:17.346233Z","iopub.status.idle":"2026-07-03T10:25:18.992516Z","shell.execute_reply.started":"2026-07-03T10:25:17.346202Z","shell.execute_reply":"2026-07-03T10:25:18.991782Z"}},"outputs":[],"execution_count":null},{"id":"6c25fbb3","cell_type":"markdown","source":"## 15. Run or resume the full repeated-seed experiment\n","metadata":{}},{"id":"1f28c205","cell_type":"code","source":"per_seed_results = run_repeated_experiments(\n    dev_df,\n    internal_test_df,\n    rsna_df,\n    X_dev_raw,\n    X_internal_raw,\n    X_rsna_raw,\n)\n\ndisplay(\n    per_seed_results.sort_values(\n        [\"cohort\", \"method\", \"seed\"]\n    )\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:25:18.993505Z","iopub.execute_input":"2026-07-03T10:25:18.994136Z","iopub.status.idle":"2026-07-03T10:25:23.041885Z","shell.execute_reply.started":"2026-07-03T10:25:18.994102Z","shell.execute_reply":"2026-07-03T10:25:23.041117Z"}},"outputs":[],"execution_count":null},{"id":"f99feda9","cell_type":"markdown","source":"## 16. Generate final manuscript tables and statistical tests\n","metadata":{}},{"id":"eae32296","cell_type":"code","source":"seed_summary = summarize_across_seeds(per_seed_results)\nseed_summary.to_csv(\n    TABLE_DIR / \"seed_summary_long.csv\",\n    index=False,\n)\n\nmanuscript_table = build_manuscript_wide_table(seed_summary)\nmanuscript_table.to_csv(\n    TABLE_DIR / \"manuscript_results_mean_sd.csv\",\n    index=False,\n)\n\nensemble_results, ensemble_store = build_ensemble_results()\nensemble_results.to_csv(\n    TABLE_DIR / \"ensemble_results.csv\",\n    index=False,\n)\n\npairwise_results = pairwise_statistics(\n    ensemble_store,\n    cohort=\"rsna_external\",\n)\npairwise_results.to_csv(\n    TABLE_DIR / \"rsna_pairwise_statistics_vs_roof_dars.csv\",\n    index=False,\n)\n\ngeneralization_gap = generalization_gap_table(ensemble_results)\ngeneralization_gap.to_csv(\n    TABLE_DIR / \"generalization_gap.csv\",\n    index=False,\n)\n\nprint(\"Manuscript mean ± SD table\")\ndisplay(manuscript_table)\n\nprint(\"Ensemble results\")\ndisplay(\n    ensemble_results.sort_values(\n        [\"cohort\", \"mcc\"],\n        ascending=[True, False],\n    )\n)\n\nprint(\"RSNA paired statistics versus ROOF-DARS\")\ndisplay(pairwise_results)\n\nprint(\"Generalization gap\")\ndisplay(generalization_gap)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:25:23.042858Z","iopub.execute_input":"2026-07-03T10:25:23.043648Z","iopub.status.idle":"2026-07-03T10:30:19.904503Z","shell.execute_reply.started":"2026-07-03T10:25:23.043614Z","shell.execute_reply":"2026-07-03T10:30:19.903721Z"}},"outputs":[],"execution_count":null},{"id":"706a964d","cell_type":"markdown","source":"## 17. Generate manuscript figures\n","metadata":{}},{"id":"04b4e541","cell_type":"code","source":"for cohort in [\"internal_test\", \"rsna_external\"]:\n    for metric in [\n        \"balanced_accuracy\",\n        \"specificity\",\n        \"mcc\",\n        \"roc_auc\",\n        \"pr_auc\",\n        \"brier\",\n        \"ece\",\n    ]:\n        plot_metric_comparison(seed_summary, cohort, metric)\n\n    plot_ensemble_curves(ensemble_store, cohort)\n    plot_reference_calibration(ensemble_store, cohort)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:30:19.905522Z","iopub.execute_input":"2026-07-03T10:30:19.905806Z","iopub.status.idle":"2026-07-03T10:30:31.533853Z","shell.execute_reply.started":"2026-07-03T10:30:19.90578Z","shell.execute_reply":"2026-07-03T10:30:31.533036Z"}},"outputs":[],"execution_count":null},{"id":"5d68ba8f","cell_type":"markdown","source":"## Manuscript use\n\nThe final manuscript should prioritize:\n\n- the **mean ± standard deviation across five seeds**;\n- the untouched full RSNA external results using natural prevalence;\n- 95% confidence intervals across seeds;\n- ensemble paired statistics only as a supplementary robustness analysis;\n- McNemar, DeLong, and paired-bootstrap tests with Holm correction;\n- calibration and generalization-gap results;\n- explicit reporting that SMOTE synthesizes **1024-dimensional feature vectors**, not images;\n- explicit reporting that RSNA `Lung Opacity` is a radiographic endpoint and not microbiologically confirmed pneumonia.\n\nDo not copy any results from a partially completed run. Use the final tables only after all\n45 seed–method experiments have completed successfully.\n","metadata":{}},{"id":"3bda9ad8","cell_type":"code","source":"manifest = {\n    \"study\": \"ROOF-DARS full RSNA final-manuscript experiment\",\n    \"seeds\": list(cfg.seeds),\n    \"methods\": METHODS,\n    \"reference_method\": REFERENCE_METHOD,\n    \"development_dataset\": \"Kermany pediatric chest X-ray\",\n    \"internal_test\": \"Original Kermany test partition\",\n    \"external_dataset\": \"RSNA Pneumonia Detection Challenge\",\n    \"external_positive\": \"Lung Opacity\",\n    \"external_negative\": \"Normal\",\n    \"external_excluded\": \"No Lung Opacity / Not Normal\",\n    \"external_prevalence_preserved\": True,\n    \"external_threshold_tuning\": False,\n    \"calibration\": \"Positive-slope affine logit calibration using Kermany calibration subset\",\n    \"threshold_strategy\": (\n        f\"Maximum MCC subject to sensitivity >= {cfg.target_sensitivity:.2f} \"\n        \"on Kermany calibration subset\"\n    ),\n    \"backbone\": \"TorchXRayVision DenseNet121\",\n    \"backbone_weights\": cfg.xrv_weights,\n    \"backbone_pretraining\": \"CheXpert\",\n    \"backbone_frozen\": True,\n    \"feature_dimension\": cfg.xrv_feature_dim,\n    \"image_size\": cfg.img_size,\n    \"duplicate_audit\": cfg.audit_exact_duplicates,\n    \"bootstrap_iterations\": cfg.bootstrap_iterations,\n    \"python\": platform.python_version(),\n    \"tensorflow\": tf.__version__,\n    \"pytorch\": torch.__version__,\n    \"torchxrayvision\": getattr(xrv, \"__version__\", \"unknown\"),\n    \"scikit_learn\": sklearn.__version__,\n    \"scipy\": scipy.__version__,\n    \"statsmodels\": statsmodels.__version__,\n}\n\n(OUTPUT_DIR / \"reproducibility_manifest.json\").write_text(\n    json.dumps(manifest, indent=2),\n    encoding=\"utf-8\",\n)\n\nprint(\"Final manuscript experiment completed.\")\nprint(\"Output directory:\", OUTPUT_DIR)\nprint(\"\\nPrimary manuscript files:\")\nfor path in [\n    TABLE_DIR / \"cohort_summary.csv\",\n    TABLE_DIR / \"per_seed_results.csv\",\n    TABLE_DIR / \"manuscript_results_mean_sd.csv\",\n    TABLE_DIR / \"ensemble_results.csv\",\n    TABLE_DIR / \"rsna_pairwise_statistics_vs_roof_dars.csv\",\n    TABLE_DIR / \"generalization_gap.csv\",\n    OUTPUT_DIR / \"reproducibility_manifest.json\",\n]:\n    print(\" -\", path)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:30:31.534896Z","iopub.execute_input":"2026-07-03T10:30:31.535281Z","iopub.status.idle":"2026-07-03T10:30:31.543925Z","shell.execute_reply.started":"2026-07-03T10:30:31.535253Z","shell.execute_reply":"2026-07-03T10:30:31.543076Z"}},"outputs":[],"execution_count":null},{"id":"e2e3c4d1-7fd5-4485-9197-8d90ce3f9db2","cell_type":"code","source":"# ---------------------------------------------------------\n# 2. Direct paired comparison on the RSNA external cohort\n# ---------------------------------------------------------\nMETHOD_A = \"DARS-Original\"\nMETHOD_B = \"Random-Under-Sampling\"\nCOHORT = \"rsna_external\"\n\nmetrics = [\n    \"accuracy\",\n    \"balanced_accuracy\",\n    \"sensitivity\",\n    \"specificity\",\n    \"precision\",\n    \"f1\",\n    \"mcc\",\n    \"roc_auc\",\n    \"pr_auc\",\n    \"brier\",\n    \"ece\",\n]\n\nexternal = per_seed_results[\n    per_seed_results[\"cohort\"] == COHORT\n].copy()\n\navailable_methods = set(external[\"method\"].unique())\n\nif METHOD_A not in available_methods:\n    raise ValueError(f\"{METHOD_A} was not found in the results.\")\n\nif METHOD_B not in available_methods:\n    raise ValueError(f\"{METHOD_B} was not found in the results.\")\n\nrows = []\n\nfor metric in metrics:\n    if metric not in external.columns:\n        continue\n\n    pivot = external.pivot_table(\n        index=\"seed\",\n        columns=\"method\",\n        values=metric,\n        aggfunc=\"first\"\n    )\n\n    pivot = pivot.dropna(subset=[METHOD_A, METHOD_B])\n\n    dars = pivot[METHOD_A].to_numpy(dtype=float)\n    rus = pivot[METHOD_B].to_numpy(dtype=float)\n\n    if len(dars) < 2:\n        continue\n\n    differences = dars - rus\n    n = len(differences)\n\n    # Paired t-test\n    t_result = ttest_rel(dars, rus)\n\n    # Wilcoxon signed-rank test\n    try:\n        w_result = wilcoxon(\n            dars,\n            rus,\n            alternative=\"two-sided\",\n            zero_method=\"wilcox\"\n        )\n        wilcoxon_stat = float(w_result.statistic)\n        wilcoxon_p = float(w_result.pvalue)\n    except ValueError:\n        wilcoxon_stat = np.nan\n        wilcoxon_p = np.nan\n\n    # 95% confidence interval for paired mean difference\n    mean_difference = float(np.mean(differences))\n    difference_sd = float(np.std(differences, ddof=1))\n    standard_error = difference_sd / np.sqrt(n)\n    critical_t = t.ppf(0.975, df=n - 1)\n\n    ci_lower = mean_difference - critical_t * standard_error\n    ci_upper = mean_difference + critical_t * standard_error\n\n    # Paired effect size: Cohen's dz\n    cohens_dz = (\n        mean_difference / difference_sd\n        if difference_sd > 0\n        else np.nan\n    )\n\n    rows.append({\n        \"metric\": metric,\n        \"n_pairs\": n,\n        \"DARS_mean\": float(np.mean(dars)),\n        \"DARS_sd\": float(np.std(dars, ddof=1)),\n        \"RUS_mean\": float(np.mean(rus)),\n        \"RUS_sd\": float(np.std(rus, ddof=1)),\n        \"difference_DARS_minus_RUS\": mean_difference,\n        \"difference_ci_lower\": ci_lower,\n        \"difference_ci_upper\": ci_upper,\n        \"paired_t_statistic\": float(t_result.statistic),\n        \"paired_t_p\": float(t_result.pvalue),\n        \"wilcoxon_statistic\": wilcoxon_stat,\n        \"wilcoxon_p\": wilcoxon_p,\n        \"cohens_dz\": cohens_dz,\n    })\n\ncomparison_df = pd.DataFrame(rows)\n\n# ---------------------------------------------------------\n# 3. Correct for testing several metrics\n# ---------------------------------------------------------\ncomparison_df[\"paired_t_p_holm\"] = multipletests(\n    comparison_df[\"paired_t_p\"].fillna(1.0),\n    method=\"holm\"\n)[1]\n\ncomparison_df[\"wilcoxon_p_holm\"] = multipletests(\n    comparison_df[\"wilcoxon_p\"].fillna(1.0),\n    method=\"holm\"\n)[1]\n\npd.set_option(\"display.max_columns\", None)\ndisplay(comparison_df.round(6))\n\n# ---------------------------------------------------------\n# 4. Save the table\n# ---------------------------------------------------------\nOUTPUT_PATH = RESULTS_PATH.parent / (\n    \"dars_vs_random_undersampling_seed_statistics.csv\"\n)\n\ncomparison_df.to_csv(OUTPUT_PATH, index=False)\nprint(\"Saved comparison to:\", OUTPUT_PATH)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-03T10:30:31.545017Z","iopub.execute_input":"2026-07-03T10:30:31.545345Z","iopub.status.idle":"2026-07-03T10:30:32.277241Z","shell.execute_reply.started":"2026-07-03T10:30:31.545325Z","shell.execute_reply":"2026-07-03T10:30:32.276327Z"}},"outputs":[],"execution_count":null},{"id":"7c58deee-1744-4a7f-9a6d-dcc702a54411","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}