{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.12","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","sourceId":15485469,"datasetId":9907068,"databundleVersionId":16409804},{"sourceType":"datasetVersion","sourceId":1488548,"datasetId":873734,"databundleVersionId":1522473}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# RSNA PE: Preprocessing, Augmentation, and Detection Pipeline\n\nThis notebook implements an end-to-end workflow in one place:\n1. Download + validate dataset\n2. No-leak study-level split\n3. Image preprocessing + enhancement + augmentation\n4. ViT image model training + evaluation\n5. GradCAM explainability\n6. XGBoost tabular baseline\n7. Quality-aware ensemble fusion","metadata":{}},{"cell_type":"code","source":"import os\nimport json\nimport random\nfrom pathlib import Path\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\n\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nfrom tqdm.auto import tqdm\n\nimport albumentations as A\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import (\n    roc_auc_score, accuracy_score, precision_score, recall_score,\n    f1_score, confusion_matrix\n)\n\nimport timm\nfrom xgboost import XGBClassifier\n\n# SEED = 42\n# random.seed(SEED)\n# np.random.seed(SEED)\n# torch.manual_seed(SEED)\n# torch.cuda.manual_seed_all(SEED)\n\n# device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n# print('Device:', device)\n\n#  safer CUDA seed\nif torch.cuda.is_available():\n    torch.cuda.manual_seed(SEED)\n    torch.cuda.manual_seed_all(SEED)\n\n#  ensure reproducibility + performance balance\ntorch.backends.cudnn.deterministic = False\ntorch.backends.cudnn.benchmark = True\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint('Device:', device)\n\n#  confirm GPU\nif device.type == 'cuda':\n    print('GPU:', torch.cuda.get_device_name(0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:00:57.497764Z","iopub.execute_input":"2026-04-02T09:00:57.498091Z","iopub.status.idle":"2026-04-02T09:00:57.505925Z","shell.execute_reply.started":"2026-04-02T09:00:57.498067Z","shell.execute_reply":"2026-04-02T09:00:57.505166Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1) Resolve Kaggle Input Paths (Separate CSV and Images)","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\n\n# Default Kaggle input paths (as requested)\nimage_root = \"/kaggle/input/datasets/vaillant/rsna-str-pe-detection-jpeg-256\"\ntrain_csv_path = \"/kaggle/input/datasets/hansajpatidar2/hp-train/train.csv\"\n\n# Convert image root to Path for downstream rglob usage\nimage_root = Path(image_root)\n\nprint('Resolved train.csv path:', train_csv_path)\nprint('Resolved image root:', image_root)\nprint('Image root exists:', image_root.exists())\nprint('Train CSV exists:', Path(train_csv_path).exists())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:01:20.275523Z","iopub.execute_input":"2026-04-02T09:01:20.275786Z","iopub.status.idle":"2026-04-02T09:01:20.28588Z","shell.execute_reply.started":"2026-04-02T09:01:20.275764Z","shell.execute_reply":"2026-04-02T09:01:20.284873Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2) Load CSV and (Optional) Fast Debug Sampling","metadata":{}},{"cell_type":"code","source":"if not Path(train_csv_path).exists():\n    raise FileNotFoundError(f'train.csv path not found: {train_csv_path}')\n\ndf = pd.read_csv(train_csv_path)\nprint('Loaded CSV from:', train_csv_path)\nprint('Rows:', len(df))\ndisplay(df.head())\n\nUSE_SAMPLE = True\nSAMPLE_ROWS = 120000\n\nif USE_SAMPLE and len(df) > SAMPLE_ROWS:\n    df = df.sample(SAMPLE_ROWS, random_state=SEED).reset_index(drop=True)\n    print(f'Using sampled dataframe: {len(df)} rows')\nelse:\n    print('Using full dataframe')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:01:23.024825Z","iopub.execute_input":"2026-04-02T09:01:23.025641Z","iopub.status.idle":"2026-04-02T09:01:26.216017Z","shell.execute_reply.started":"2026-04-02T09:01:23.0256Z","shell.execute_reply":"2026-04-02T09:01:26.215251Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 3) Build JPEG Index and Map CSV to File Paths","metadata":{}},{"cell_type":"code","source":"def index_jpegs(root: Path):\n    sop_to_path = {}\n    jpg_count = 0\n    duplicate_suffix_count = 0\n\n    for jpg in tqdm(root.rglob('*.jpg'), desc='Indexing JPGs'):\n        jpg_count += 1\n        stem = jpg.stem\n        parts = stem.split('_')\n        if len(parts) < 2:\n            continue\n\n        sop_suffix = parts[-1]\n        if sop_suffix in sop_to_path:\n            duplicate_suffix_count += 1\n        sop_to_path[sop_suffix] = str(jpg)\n\n    print('Image dataset root:', root)\n    print('Indexed JPG files:', jpg_count)\n    print('Indexed SOP suffix keys:', len(sop_to_path))\n    print('Duplicate suffix collisions resolved by latest path:', duplicate_suffix_count)\n    return sop_to_path\n\nsop_to_path = index_jpegs(image_root)\n\ndf['sop_suffix'] = df['SOPInstanceUID'].astype(str).apply(lambda x: x.split('.')[-1])\ndf['file_path'] = df['sop_suffix'].map(sop_to_path)\ndf['file_found'] = df['file_path'].notna()\n\nmissing = (~df['file_found']).sum()\nprint(f'Found files for {df[\"file_found\"].sum()} / {len(df)} rows')\nprint('Missing:', missing)\n\ndf_valid = df[df['file_found']].copy().reset_index(drop=True)\nprint('Valid rows:', len(df_valid))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T06:56:01.677782Z","iopub.execute_input":"2026-04-02T06:56:01.678093Z","iopub.status.idle":"2026-04-02T08:00:55.532565Z","shell.execute_reply.started":"2026-04-02T06:56:01.678065Z","shell.execute_reply":"2026-04-02T08:00:55.531819Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4) No-Leak Study-Level Stratified Split (70/15/15)","metadata":{}},{"cell_type":"code","source":"study_col = 'StudyInstanceUID'\nlabel_col = 'pe_present_on_image'\n\nstudy_df = df_valid.groupby(study_col)[label_col].max().reset_index(name='has_pe')\n\ntrain_studies, temp_studies = train_test_split(\n    study_df,\n    test_size=0.30,\n    stratify=study_df['has_pe'],\n    random_state=SEED\n)\n\nval_studies, test_studies = train_test_split(\n    temp_studies,\n    test_size=0.50,\n    stratify=temp_studies['has_pe'],\n    random_state=SEED\n)\n\ntrain_set = set(train_studies[study_col])\nval_set = set(val_studies[study_col])\ntest_set = set(test_studies[study_col])\n\nassert len(train_set & val_set) == 0\nassert len(train_set & test_set) == 0\nassert len(val_set & test_set) == 0\n\ndf_train = df_valid[df_valid[study_col].isin(train_set)].copy().reset_index(drop=True)\ndf_val = df_valid[df_valid[study_col].isin(val_set)].copy().reset_index(drop=True)\ndf_test = df_valid[df_valid[study_col].isin(test_set)].copy().reset_index(drop=True)\n\nprint('Train rows:', len(df_train), 'PE ratio:', df_train[label_col].mean())\nprint('Val rows:', len(df_val), 'PE ratio:', df_val[label_col].mean())\nprint('Test rows:', len(df_test), 'PE ratio:', df_test[label_col].mean())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:01:32.386949Z","iopub.execute_input":"2026-04-02T09:01:32.387514Z","iopub.status.idle":"2026-04-02T09:01:32.493403Z","shell.execute_reply.started":"2026-04-02T09:01:32.387485Z","shell.execute_reply":"2026-04-02T09:01:32.492733Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"processed_dir = Path('./processed')\nprocessed_dir.mkdir(parents=True, exist_ok=True)\n\ndf_train.to_csv(processed_dir / 'train_metadata.csv', index=False)\ndf_val.to_csv(processed_dir / 'val_metadata.csv', index=False)\ndf_test.to_csv(processed_dir / 'test_metadata.csv', index=False)\n\nsplit_info = {\n    'seed': SEED,\n    'train_studies': list(train_set),\n    'val_studies': list(val_set),\n    'test_studies': list(test_set)\n}\nwith open(processed_dir / 'split_info.json', 'w', encoding='utf-8') as f:\n    json.dump(split_info, f, indent=2)\n\nprint('Saved split metadata to', processed_dir.resolve())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T08:28:37.392842Z","iopub.execute_input":"2026-04-02T08:28:37.39356Z","iopub.status.idle":"2026-04-02T08:28:38.415365Z","shell.execute_reply.started":"2026-04-02T08:28:37.3935Z","shell.execute_reply":"2026-04-02T08:28:38.414589Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5) Compute Medical Normalization Statistics from Training Images","metadata":{}},{"cell_type":"code","source":"def compute_mean_std(image_paths, max_samples=2000):\n    paths = list(image_paths)\n    if len(paths) > max_samples:\n        paths = random.sample(paths, max_samples)\n\n    channel_sum = np.zeros(3, dtype=np.float64)\n    channel_sum_sq = np.zeros(3, dtype=np.float64)\n    total_pixels = 0\n\n    for p in tqdm(paths, desc='Computing mean/std'):\n        img = np.array(Image.open(p).convert('RGB'), dtype=np.float32) / 255.0\n        h, w, _ = img.shape\n        total_pixels += h * w\n        channel_sum += img.reshape(-1, 3).sum(axis=0)\n        channel_sum_sq += (img.reshape(-1, 3) ** 2).sum(axis=0)\n\n    mean = channel_sum / total_pixels\n    std = np.sqrt((channel_sum_sq / total_pixels) - (mean ** 2))\n    return mean.tolist(), std.tolist()\n\nMEDICAL_MEAN, MEDICAL_STD = compute_mean_std(df_train['file_path'].tolist(), max_samples=1500)\nprint('MEDICAL_MEAN =', MEDICAL_MEAN)\nprint('MEDICAL_STD  =', MEDICAL_STD)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:01:52.618003Z","iopub.execute_input":"2026-04-02T09:01:52.618284Z","iopub.status.idle":"2026-04-02T09:02:15.389139Z","shell.execute_reply.started":"2026-04-02T09:01:52.61826Z","shell.execute_reply":"2026-04-02T09:02:15.388328Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6) Image Enhancement and Augmentation (as requested)","metadata":{}},{"cell_type":"code","source":"def apply_clahe(image_np, clip_limit=2.0, tile_size=8):\n    clahe = cv2.createCLAHE(clipLimit=clip_limit, tileGridSize=(tile_size, tile_size))\n    enhanced = np.stack([clahe.apply(image_np[:, :, c]) for c in range(3)], axis=2)\n    return enhanced\n\ndef apply_unsharp_mask(image_np, sigma=1.0, strength=0.5):\n    blurred = cv2.GaussianBlur(image_np, ksize=(0, 0), sigmaX=sigma)\n    sharpened = cv2.addWeighted(image_np, 1.0 + strength, blurred, -strength, 0)\n    return np.clip(sharpened, 0, 255).astype(np.uint8)\n\ndef enhance_image(image_np, clahe_clip=2.0, unsharp_sigma=1.0, unsharp_strength=0.5):\n    x = apply_clahe(image_np, clip_limit=clahe_clip, tile_size=8)\n    x = apply_unsharp_mask(x, sigma=unsharp_sigma, strength=unsharp_strength)\n    return x\n\ndef get_train_transforms(img_size=224):\n    return A.Compose([\n        A.Resize(img_size, img_size),\n        A.HorizontalFlip(p=0.5),\n        A.Rotate(limit=10, p=0.5, border_mode=cv2.BORDER_REFLECT_101),\n        A.RandomBrightnessContrast(brightness_limit=0.2, contrast_limit=0.2, p=0.5),\n        A.GaussianBlur(blur_limit=3, p=0.3),\n        A.GaussNoise(std_range=(0.03, 0.10), p=0.3),\n    ])\n\ndef get_eval_transforms(img_size=224):\n    return A.Compose([A.Resize(img_size, img_size)])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:02:15.390406Z","iopub.execute_input":"2026-04-02T09:02:15.390717Z","iopub.status.idle":"2026-04-02T09:02:15.397985Z","shell.execute_reply.started":"2026-04-02T09:02:15.390693Z","shell.execute_reply":"2026-04-02T09:02:15.397233Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7) Dataset, Dataloaders, and Training Utilities","metadata":{}},{"cell_type":"code","source":"class PEImageDataset(Dataset):\n    def __init__(self, df, transform, mean, std, enhance=True):\n        self.df = df.reset_index(drop=True)\n        self.transform = transform\n        self.mean = np.array(mean, dtype=np.float32).reshape(1, 1, 3)\n        self.std = np.array(std, dtype=np.float32).reshape(1, 1, 3)\n        self.enhance = enhance\n\n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, idx):\n        row = self.df.iloc[idx]\n        image = np.array(Image.open(row['file_path']).convert('RGB'), dtype=np.uint8)\n\n        if self.enhance:\n            image = enhance_image(image, clahe_clip=2.0, unsharp_sigma=1.0, unsharp_strength=0.5)\n\n        image = self.transform(image=image)['image']\n        image = image.astype(np.float32) / 255.0\n        image = (image - self.mean) / (self.std + 1e-7)\n        image = np.transpose(image, (2, 0, 1))\n\n        label = np.float32(row['pe_present_on_image'])\n        sop_suffix = row['sop_suffix']\n\n        return torch.tensor(image, dtype=torch.float32), torch.tensor(label, dtype=torch.float32), sop_suffix\n\ntrain_ds = PEImageDataset(df_train, get_train_transforms(224), MEDICAL_MEAN, MEDICAL_STD, enhance=True)\nval_ds = PEImageDataset(df_val, get_eval_transforms(224), MEDICAL_MEAN, MEDICAL_STD, enhance=True)\ntest_ds = PEImageDataset(df_test, get_eval_transforms(224), MEDICAL_MEAN, MEDICAL_STD, enhance=True)\n\nBATCH_SIZE = 32 if device.type == 'cuda' else 8\ntrain_loader = DataLoader(train_ds, batch_size=BATCH_SIZE, shuffle=True, num_workers=2, pin_memory=True)\nval_loader = DataLoader(val_ds, batch_size=BATCH_SIZE, shuffle=False, num_workers=2, pin_memory=True)\ntest_loader = DataLoader(test_ds, batch_size=BATCH_SIZE, shuffle=False, num_workers=2, pin_memory=True)\n\nneg = (df_train['pe_present_on_image'] == 0).sum()\npos = (df_train['pe_present_on_image'] == 1).sum()\npos_weight = torch.tensor([max(1.0, neg / max(1, pos))], device=device, dtype=torch.float32)\nprint('pos_weight:', float(pos_weight.item()))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T09:02:20.067295Z","iopub.execute_input":"2026-04-02T09:02:20.068029Z","iopub.status.idle":"2026-04-02T09:02:20.111827Z","shell.execute_reply.started":"2026-04-02T09:02:20.067997Z","shell.execute_reply":"2026-04-02T09:02:20.11119Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from torch.amp import autocast, GradScaler\n\nscaler = GradScaler(device='cuda')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T12:30:54.603621Z","iopub.execute_input":"2026-04-02T12:30:54.604313Z","iopub.status.idle":"2026-04-02T12:30:54.607569Z","shell.execute_reply.started":"2026-04-02T12:30:54.604281Z","shell.execute_reply":"2026-04-02T12:30:54.606752Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model = timm.create_model('vit_base_patch16_224', pretrained=True, num_classes=1).to(device)\ncriterion = nn.BCEWithLogitsLoss(pos_weight=pos_weight)\noptimizer = torch.optim.AdamW(model.parameters(), lr=1e-5, weight_decay=1e-4)\nscheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=10, eta_min=1e-6)\n\ndef run_epoch(model, loader, optimizer=None):\n    train_mode = optimizer is not None\n    model.train(train_mode)\n\n    losses = []\n    y_true = []\n    y_prob = []\n\n    for images, labels, _ in tqdm(loader, disable=False):\n        images = images.to(device, non_blocking=True)\n        labels = labels.to(device, non_blocking=True).unsqueeze(1)\n\n        with torch.set_grad_enabled(train_mode):\n            # logits = model(images)\n            # loss = criterion(logits, labels)\n            with autocast(device_type='cuda',enabled=train_mode):\n                logits = model(images)\n                loss = criterion(logits, labels)\n\n            if train_mode:\n                optimizer.zero_grad()\n                # loss.backward()\n                scaler.scale(loss).backward()\n                torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n                # optimizer.step()\n                scaler.step(optimizer)\n                scaler.update()\n\n        probs = torch.sigmoid(logits).detach().cpu().numpy().ravel()\n        y_prob.extend(probs.tolist())\n        y_true.extend(labels.detach().cpu().numpy().ravel().tolist())\n        losses.append(float(loss.item()))\n\n    y_true = np.array(y_true)\n    y_prob = np.array(y_prob)\n    y_pred = (y_prob >= 0.5).astype(int)\n\n    metrics = {\n        'loss': float(np.mean(losses)) if losses else np.nan,\n        'auc': float(roc_auc_score(y_true, y_prob)) if len(np.unique(y_true)) > 1 else np.nan,\n        'acc': float(accuracy_score(y_true, y_pred)),\n        'precision': float(precision_score(y_true, y_pred, zero_division=0)),\n        'recall': float(recall_score(y_true, y_pred, zero_division=0)),\n        'f1': float(f1_score(y_true, y_pred, zero_division=0)),\n    }\n    return metrics\n\ndef evaluate_with_probs(model, loader):\n    model.eval()\n    y_true, y_prob, sop_keys = [], [], []\n\n    with torch.no_grad():\n        for images, labels, sops in tqdm(loader, disable=False):\n            images = images.to(device, non_blocking=True)\n            logits = model(images)\n            probs = torch.sigmoid(logits).cpu().numpy().ravel()\n\n            y_prob.extend(probs.tolist())\n            y_true.extend(labels.numpy().ravel().tolist())\n            sop_keys.extend(list(sops))\n\n    y_true = np.array(y_true)\n    y_prob = np.array(y_prob)\n    y_pred = (y_prob >= 0.5).astype(int)\n\n    metrics = {\n        'auc': float(roc_auc_score(y_true, y_prob)) if len(np.unique(y_true)) > 1 else np.nan,\n        'acc': float(accuracy_score(y_true, y_pred)),\n        'precision': float(precision_score(y_true, y_pred, zero_division=0)),\n        'recall': float(recall_score(y_true, y_pred, zero_division=0)),\n        'f1': float(f1_score(y_true, y_pred, zero_division=0)),\n    }\n\n    pred_df = pd.DataFrame({\n        'sop_suffix': sop_keys,\n        'y_true': y_true,\n        'p_image': y_prob\n    })\n\n    return metrics, pred_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T12:31:46.869352Z","iopub.execute_input":"2026-04-02T12:31:46.869783Z","iopub.status.idle":"2026-04-02T12:31:58.639042Z","shell.execute_reply.started":"2026-04-02T12:31:46.869753Z","shell.execute_reply":"2026-04-02T12:31:58.638459Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"EPOCHS = 10\nPATIENCE = 2\n\nbest_auc = -np.inf\nbest_state = None\npatience_counter = 0\n\nhistory = []\nfor epoch in range(1, EPOCHS + 1):\n    print(f'\\nEpoch {epoch}/{EPOCHS}')\n\n    train_metrics = run_epoch(model, train_loader, optimizer=optimizer)\n    val_metrics = run_epoch(model, val_loader, optimizer=None)\n    scheduler.step()\n\n    row = {'epoch': epoch, **{f'train_{k}': v for k, v in train_metrics.items()}, **{f'val_{k}': v for k, v in val_metrics.items()}}\n    history.append(row)\n    print(row)\n\n    if np.isfinite(val_metrics['auc']) and val_metrics['auc'] > best_auc:\n        best_auc = val_metrics['auc']\n        torch.save(model.state_dict(), '/kaggle/working/best_model.pth')\n        patience_counter = 0\n    else:\n        patience_counter += 1\n\n    if patience_counter >= PATIENCE:\n        print('Early stopping triggered')\n        break\n\nhistory_df = pd.DataFrame(history)\ndisplay(history_df)\n\nif best_state is not None:\n    model.load_state_dict(torch.load('/kaggle/working/best_model.pth'))\n    print('Loaded best model state with val AUC:', best_auc)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T12:34:12.263446Z","iopub.execute_input":"2026-04-02T12:34:12.264171Z","iopub.status.idle":"2026-04-02T13:39:55.14445Z","shell.execute_reply.started":"2026-04-02T12:34:12.264137Z","shell.execute_reply":"2026-04-02T13:39:55.143574Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8) Evaluate Image Model on Validation and Test","metadata":{}},{"cell_type":"code","source":"val_metrics_img, val_pred_img = evaluate_with_probs(model, val_loader)\ntest_metrics_img, test_pred_img = evaluate_with_probs(model, test_loader)\n\nprint('Image model validation metrics:', val_metrics_img)\nprint('Image model test metrics:', test_metrics_img)\n\ndisplay(test_pred_img.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T13:41:23.232677Z","iopub.execute_input":"2026-04-02T13:41:23.233591Z","iopub.status.idle":"2026-04-02T13:48:27.694764Z","shell.execute_reply.started":"2026-04-02T13:41:23.233545Z","shell.execute_reply":"2026-04-02T13:48:27.693943Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 9) GradCAM Explainability (ExAI section)","metadata":{}},{"cell_type":"markdown","source":"NOTE: below cell is one time run only","metadata":{}},{"cell_type":"code","source":"# # One-time setup for GradCAM cell (safe to run independently).\n# import os\n# import subprocess\n# import sys\n# from pathlib import Path\n\n# def gradcam_preflight_setup():\n#     kaggle_working = Path('/kaggle/working')\n#     if kaggle_working.exists():\n#         try:\n#             os.chdir(kaggle_working)\n#             print('Working directory set to writable path:', Path.cwd())\n#         except Exception as e:\n#             print(f'Could not switch to /kaggle/working: {e}')\n\n#     ip = get_ipython()\n#     if ip is not None and hasattr(ip, 'history_manager'):\n#         try:\n#             # Avoid readonly sqlite history errors in Kaggle mounted datasets.\n#             ip.history_manager.enabled = False\n#             print('IPython history DB writes disabled for this session')\n#         except Exception as e:\n#             print(f'Could not disable history manager: {e}')\n\n# def ensure_package(module_name, pip_name=None):\n#     try:\n#         __import__(module_name)\n#         print(f'{module_name} is already available')\n#         return True\n#     except ModuleNotFoundError:\n#         pkg = pip_name or module_name\n#         print(f'Installing missing package: {pkg}')\n#         try:\n#             subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', pkg])\n#             __import__(module_name)\n#             print(f'Installed and imported: {module_name}')\n#             return True\n#         except Exception as e:\n#             print(f'Package install failed ({pkg}): {e}')\n#             return False\n\n# gradcam_preflight_setup()\n# GRADCAM_READY = ensure_package('pytorch_grad_cam', 'grad-cam')\n# print('GradCAM ready:', GRADCAM_READY)\n# if not GRADCAM_READY:\n#     print('If this is Kaggle, enable Internet in notebook settings and rerun this cell.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T13:49:24.565988Z","iopub.execute_input":"2026-04-02T13:49:24.566996Z","iopub.status.idle":"2026-04-02T13:49:24.575538Z","shell.execute_reply.started":"2026-04-02T13:49:24.566957Z","shell.execute_reply":"2026-04-02T13:49:24.574966Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ntry:\n    from pytorch_grad_cam import GradCAM\n    from pytorch_grad_cam.utils.image import show_cam_on_image\n\n    def vit_reshape_transform(tensor, height=14, width=14):\n        # ViT tokens -> feature map for CAM\n        result = tensor[:, 1:, :].reshape(tensor.size(0), height, width, tensor.size(2))\n        result = result.permute(0, 3, 1, 2)\n        return result\n\n    sample_image, sample_label, _ = test_ds[0]\n    input_tensor = sample_image.unsqueeze(0).to(device)\n\n    target_layers = [model.blocks[-1].norm1]\n    cam = GradCAM(model=model, target_layers=target_layers, reshape_transform=vit_reshape_transform)\n    grayscale_cam = cam(input_tensor=input_tensor)[0]\n\n    img_vis = sample_image.permute(1, 2, 0).numpy()\n    img_vis = (img_vis - img_vis.min()) / (img_vis.max() - img_vis.min() + 1e-8)\n    overlay = show_cam_on_image(img_vis.astype(np.float32), grayscale_cam, use_rgb=True)\n\n    plt.figure(figsize=(10, 4))\n    plt.subplot(1, 2, 1)\n    plt.imshow(img_vis)\n    plt.title(f'Input (label={int(sample_label.item())})')\n    plt.axis('off')\n\n    plt.subplot(1, 2, 2)\n    plt.imshow(overlay)\n    plt.title('GradCAM')\n    plt.axis('off')\n    plt.show()\n\nexcept Exception as e:\n    print('GradCAM skipped:', e)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T13:49:59.33461Z","iopub.execute_input":"2026-04-02T13:49:59.335054Z","iopub.status.idle":"2026-04-02T13:51:19.764989Z","shell.execute_reply.started":"2026-04-02T13:49:59.335012Z","shell.execute_reply":"2026-04-02T13:51:19.764232Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 10) Tabular Model (XGBoost)","metadata":{}},{"cell_type":"code","source":"candidate_features = [\n    'negative_exam_for_pe', 'qa_motion', 'qa_contrast', 'flow_artifact',\n    'rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1', 'leftsided_pe', 'rightsided_pe',\n    'central_pe', 'chronic_pe', 'acute_and_chronic_pe',\n    'true_filling_defect_not_pe', 'indeterminate'\n]\n\ntab_features = [c for c in candidate_features if c in df_valid.columns]\nprint('Tabular features used:', tab_features)\n\nif len(tab_features) == 0:\n    print('No tabular features available in CSV. Skipping XGBoost + ensemble.')\n    xgb_model = None\nelse:\n    for frame in [df_train, df_val, df_test]:\n        frame['rv_enlarged'] = frame.get('rv_lv_ratio_gte_1', 0).fillna(0)\n        frame['rv_normal'] = frame.get('rv_lv_ratio_lt_1', 0).fillna(0)\n        frame['rv_uncertainty'] = 1 - (frame['rv_enlarged'] + frame['rv_normal'])\n        frame['bilateral_pe'] = frame.get('leftsided_pe', 0).fillna(0) * frame.get('rightsided_pe', 0).fillna(0)\n        frame['any_pe_detected'] = ((frame.get('leftsided_pe', 0).fillna(0) + frame.get('rightsided_pe', 0).fillna(0)) > 0).astype(int)\n\n    derived = ['rv_enlarged', 'rv_normal', 'rv_uncertainty', 'bilateral_pe', 'any_pe_detected']\n    tab_features = tab_features + derived\n\n    X_train = df_train[tab_features].fillna(0)\n    y_train = df_train['pe_present_on_image'].astype(int)\n    X_val = df_val[tab_features].fillna(0)\n    y_val = df_val['pe_present_on_image'].astype(int)\n    X_test = df_test[tab_features].fillna(0)\n    y_test = df_test['pe_present_on_image'].astype(int)\n\n    neg_tab = (y_train == 0).sum()\n    pos_tab = (y_train == 1).sum()\n    scale_pos_weight = float(neg_tab / max(1, pos_tab))\n\n    xgb_model = XGBClassifier(\n        max_depth=6,\n        learning_rate=0.05,\n        n_estimators=150,\n        subsample=0.8,\n        colsample_bytree=0.8,\n        scale_pos_weight=scale_pos_weight,\n        objective='binary:logistic',\n        random_state=SEED,\n        n_jobs=-1\n    )\n\n    xgb_model.fit(X_train, y_train)\n    p_val_tab = xgb_model.predict_proba(X_val)[:, 1]\n    p_test_tab = xgb_model.predict_proba(X_test)[:, 1]\n\n    print('XGBoost val AUC:', roc_auc_score(y_val, p_val_tab) if len(np.unique(y_val)) > 1 else np.nan)\n    print('XGBoost test AUC:', roc_auc_score(y_test, p_test_tab) if len(np.unique(y_test)) > 1 else np.nan)\n\n    val_pred_tab = pd.DataFrame({'sop_suffix': df_val['sop_suffix'].values, 'p_tabular': p_val_tab})\n    test_pred_tab = pd.DataFrame({'sop_suffix': df_test['sop_suffix'].values, 'p_tabular': p_test_tab})","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T13:55:17.426651Z","iopub.execute_input":"2026-04-02T13:55:17.427336Z","iopub.status.idle":"2026-04-02T13:55:18.152349Z","shell.execute_reply.started":"2026-04-02T13:55:17.427305Z","shell.execute_reply":"2026-04-02T13:55:18.151653Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 11) Quality-Aware Ensemble (Image + Tabular)","metadata":{}},{"cell_type":"code","source":"def compute_metrics(y_true, p):\n    y_true = np.array(y_true)\n    p = np.array(p)\n    y_pred = (p >= 0.5).astype(int)\n    tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()\n    return {\n        'auc': float(roc_auc_score(y_true, p)) if len(np.unique(y_true)) > 1 else np.nan,\n        'accuracy': float(accuracy_score(y_true, y_pred)),\n        'precision': float(precision_score(y_true, y_pred, zero_division=0)),\n        'recall': float(recall_score(y_true, y_pred, zero_division=0)),\n        'f1': float(f1_score(y_true, y_pred, zero_division=0)),\n        'specificity': float(tn / max(1, (tn + fp)))\n    }\n\nif 'xgb_model' in globals() and xgb_model is not None:\n    ens_df = test_pred_img.merge(test_pred_tab, on='sop_suffix', how='inner')\n    meta_cols = ['sop_suffix', 'qa_motion', 'qa_contrast', 'pe_present_on_image']\n    meta_subset = df_test[[c for c in meta_cols if c in df_test.columns]].copy()\n    ens_df = ens_df.merge(meta_subset, on='sop_suffix', how='left')\n\n    qa_motion = ens_df['qa_motion'] if 'qa_motion' in ens_df.columns else 0\n    qa_contrast = ens_df['qa_contrast'] if 'qa_contrast' in ens_df.columns else 0\n\n    alpha = np.where((qa_motion == 1) | (qa_contrast == 1), 0.5, 0.2)\n    ens_df['p_final'] = alpha * ens_df['p_tabular'] + (1 - alpha) * ens_df['p_image']\n\n    y_true = ens_df['pe_present_on_image'].values\n    ensemble_metrics = compute_metrics(y_true, ens_df['p_final'].values)\n    image_only_metrics = compute_metrics(y_true, ens_df['p_image'].values)\n    tab_only_metrics = compute_metrics(y_true, ens_df['p_tabular'].values)\n\n    print('Image-only metrics:', image_only_metrics)\n    print('Tabular-only metrics:', tab_only_metrics)\n    print('Ensemble metrics:', ensemble_metrics)\n\n    ens_df['predicted_label'] = (ens_df['p_final'] >= 0.5).astype(int)\n    ens_df['actual_label'] = ens_df['pe_present_on_image'].astype(int)\n    ens_df['predicted_output'] = ens_df.apply(\n        lambda r: (\n            f\"Predicted: {'PE' if r['predicted_label'] == 1 else 'Not PE'} | \"\n            f\"Actual: {'PE' if r['actual_label'] == 1 else 'Not PE'}\"\n        ),\n        axis=1\n    )\n\n    sample_n = min(5, len(ens_df))\n    display(\n        ens_df[\n            ['sop_suffix', 'p_image', 'p_tabular', 'p_final', 'pe_present_on_image', 'predicted_output']\n        ].sample(n=sample_n).reset_index(drop=True)\n    )\nelse:\n    print('Ensemble skipped because tabular model was not trained (features unavailable).')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T14:10:14.371636Z","iopub.execute_input":"2026-04-02T14:10:14.371924Z","iopub.status.idle":"2026-04-02T14:10:14.550249Z","shell.execute_reply.started":"2026-04-02T14:10:14.371901Z","shell.execute_reply":"2026-04-02T14:10:14.549456Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 12) Save Artifacts","metadata":{}},{"cell_type":"code","source":"models_dir = Path('./models')\nresults_dir = Path('./results')\nmodels_dir.mkdir(parents=True, exist_ok=True)\nresults_dir.mkdir(parents=True, exist_ok=True)\n\ntorch.save(model.state_dict(), models_dir / 'vit_image_model.pth')\nhistory_df.to_csv(results_dir / 'image_training_history.csv', index=False)\ntest_pred_img.to_csv(results_dir / 'image_test_predictions.csv', index=False)\n\nartifacts = {\n    'medical_mean': MEDICAL_MEAN,\n    'medical_std': MEDICAL_STD,\n    'image_val_metrics': val_metrics_img,\n    'image_test_metrics': test_metrics_img\n}\n\nif 'ensemble_metrics' in globals():\n    artifacts['ensemble_test_metrics'] = ensemble_metrics\n\nwith open(results_dir / 'metrics.json', 'w', encoding='utf-8') as f:\n    json.dump(artifacts, f, indent=2)\n\nprint('Saved model + metrics + predictions')\nprint('models:', models_dir.resolve())\nprint('results:', results_dir.resolve())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T13:57:46.003573Z","iopub.execute_input":"2026-04-02T13:57:46.004189Z","iopub.status.idle":"2026-04-02T13:57:46.425171Z","shell.execute_reply.started":"2026-04-02T13:57:46.004158Z","shell.execute_reply":"2026-04-02T13:57:46.424438Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Comparison plots: Image model (ViT) vs XGBoost\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\nfrom sklearn.metrics import (\n    roc_auc_score, roc_curve,\n    precision_recall_curve, average_precision_score,\n    confusion_matrix, classification_report\n\n)\n\nif 'test_pred_img' not in globals():\n    raise RuntimeError(\"test_pred_img not found. Run the image-model evaluation cell first.\")\n\ncmp_df = test_pred_img[['sop_suffix', 'y_true', 'p_image']].copy()\ncmp_df['y_true'] = cmp_df['y_true'].astype(int)\n\nhas_tabular = 'test_pred_tab' in globals() and test_pred_tab is not None\nif has_tabular:\n    cmp_df = cmp_df.merge(\n        test_pred_tab[['sop_suffix', 'p_tabular']],\n        on='sop_suffix',\n        how='inner'\n    )\n    if len(cmp_df) == 0:\n        raise RuntimeError('No overlapping sop_suffix keys between image and tabular predictions.')\nelse:\n    print('test_pred_tab is unavailable. Showing image-model plots only.')\n\nmodels = {'Image (ViT)': 'p_image'}\nif has_tabular:\n    models['XGBoost'] = 'p_tabular'\n\nif cmp_df['y_true'].nunique() < 2:\n    print('Warning: only one class present in y_true; ROC/PR AUC may be undefined.')\n\ndef plot_confusion(ax, y_true, y_pred, title):\n    cm = confusion_matrix(y_true, y_pred, labels=[0, 1])\n    im = ax.imshow(cm, cmap='Blues')\n    ax.set_title(title)\n    ax.set_xticks([0, 1])\n    ax.set_yticks([0, 1])\n    ax.set_xticklabels(['Not PE', 'PE'])\n    ax.set_yticklabels(['Not PE', 'PE'])\n    ax.set_xlabel('Predicted')\n    ax.set_ylabel('True')\n\n    for i in range(cm.shape[0]):\n        for j in range(cm.shape[1]):\n            ax.text(j, i, str(cm[i, j]), ha='center', va='center', color='black')\n\nsummary_rows = []\nclass_rows = []\ny_true = cmp_df['y_true'].values\n\nfor model_name, prob_col in models.items():\n    probs = cmp_df[prob_col].values.astype(float)\n    y_pred = (probs >= 0.5).astype(int)\n\n    auc_val = roc_auc_score(y_true, probs) if np.unique(y_true).size > 1 else np.nan\n    ap_val = average_precision_score(y_true, probs) if np.unique(y_true).size > 1 else np.nan\n\n    report = classification_report(\n        y_true, y_pred, labels=[0, 1],\n        target_names=['Not PE', 'PE'],\n        output_dict=True, zero_division=0\n    )\n\n    summary_rows.append({\n        'model': model_name,\n        'roc_auc': float(auc_val) if np.isfinite(auc_val) else np.nan,\n        'pr_auc': float(ap_val) if np.isfinite(ap_val) else np.nan,\n        'accuracy': float(report['accuracy'])\n    })\n\n    for cls_name in ['Not PE', 'PE']:\n        class_rows.append({\n            'model': model_name,\n            'class': cls_name,\n            'precision': float(report[cls_name]['precision']),\n            'recall': float(report[cls_name]['recall']),\n            'f1': float(report[cls_name]['f1-score']),\n            'support': int(report[cls_name]['support'])\n        })\n\nsummary_df = pd.DataFrame(summary_rows)\nclass_df = pd.DataFrame(class_rows)\n\n# Figure 1: ROC + PR + confusion matrices\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\n\nax_roc = axes[0, 0]\nfor model_name, prob_col in models.items():\n    probs = cmp_df[prob_col].values.astype(float)\n    if np.unique(y_true).size > 1:\n        fpr, tpr, _ = roc_curve(y_true, probs)\n        auc_val = roc_auc_score(y_true, probs)\n        ax_roc.plot(fpr, tpr, label=f'{model_name} (AUC={auc_val:.4f})')\nax_roc.plot([0, 1], [0, 1], '--', color='gray')\nax_roc.set_title('ROC Curve')\nax_roc.set_xlabel('False Positive Rate')\nax_roc.set_ylabel('True Positive Rate')\nax_roc.legend(loc='lower right')\nax_roc.grid(alpha=0.3)\n\nax_pr = axes[0, 1]\nfor model_name, prob_col in models.items():\n    probs = cmp_df[prob_col].values.astype(float)\n    if np.unique(y_true).size > 1:\n        precision, recall, _ = precision_recall_curve(y_true, probs)\n        ap_val = average_precision_score(y_true, probs)\n        ax_pr.plot(recall, precision, label=f'{model_name} (AP={ap_val:.4f})')\nax_pr.set_title('Precision-Recall Curve')\nax_pr.set_xlabel('Recall')\nax_pr.set_ylabel('Precision')\nax_pr.legend(loc='lower left')\nax_pr.grid(alpha=0.3)\n\nplot_confusion(\n    axes[1, 0],\n    y_true,\n    (cmp_df['p_image'].values >= 0.5).astype(int),\n    'Confusion Matrix: Image (ViT)'\n)\n\nif 'p_tabular' in cmp_df.columns:\n    plot_confusion(\n        axes[1, 1],\n        y_true,\n        (cmp_df['p_tabular'].values >= 0.5).astype(int),\n        'Confusion Matrix: XGBoost'\n    )\nelse:\n    axes[1, 1].axis('off')\n    axes[1, 1].text(0.5, 0.5, 'XGBoost predictions unavailable', ha='center', va='center', fontsize=12)\n\nplt.tight_layout()\nplt.show()\n\n# Figure 2: Class-wise precision and recall (PE / Not PE)\nfig2, axes2 = plt.subplots(1, 2, figsize=(14, 4))\nfor idx, metric in enumerate(['precision', 'recall']):\n    pivot_df = class_df.pivot(index='model', columns='class', values=metric)\n    pivot_df.plot(kind='bar', ax=axes2[idx], rot=0)\n    axes2[idx].set_title(f'Class-wise {metric.capitalize()}')\n    axes2[idx].set_ylabel(metric.capitalize())\n    axes2[idx].set_ylim(0, 1.05)\n    axes2[idx].grid(axis='y', alpha=0.3)\n    axes2[idx].legend(title='Class')\n\nplt.tight_layout()\nplt.show()\n\nprint('Overall model summary:')\ndisplay(summary_df.sort_values('roc_auc', ascending=False))\n\nprint('Class-wise metrics (includes PE and Not PE precision/recall):')\ndisplay(class_df.sort_values(['model', 'class']).reset_index(drop=True))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-02T14:00:07.870286Z","iopub.execute_input":"2026-04-02T14:00:07.870559Z","iopub.status.idle":"2026-04-02T14:00:08.775016Z","shell.execute_reply.started":"2026-04-02T14:00:07.870537Z","shell.execute_reply":"2026-04-02T14:00:08.774356Z"}},"outputs":[],"execution_count":null}]}