{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.12.13"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":154281},{"sourceType":"datasetVersion","sourceId":18956429},{"sourceType":"datasetVersion","sourceId":18673646},{"sourceType":"datasetVersion","sourceId":18229736},{"sourceType":"datasetVersion","sourceId":18842180},{"sourceType":"datasetVersion","sourceId":18839182},{"sourceType":"datasetVersion","sourceId":18757740},{"sourceType":"datasetVersion","sourceId":18673450},{"sourceType":"datasetVersion","sourceId":18875869},{"sourceType":"datasetVersion","sourceId":18879001},{"sourceType":"datasetVersion","sourceId":18716507},{"sourceType":"kernelVersion","sourceId":342671664},{"sourceType":"kernelVersion","sourceId":342849430},{"sourceType":"modelInstanceVersion","sourceId":4533},{"sourceType":"datasetVersion","sourceId":18706996},{"sourceType":"datasetVersion","sourceId":18715672},{"sourceType":"datasetVersion","sourceId":19003959},{"sourceType":"modelInstanceVersion","sourceId":4534}],"dockerImageVersionId":31430,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"rsna_optimization":{"official_source_score":0.891,"revision":"v66-v65-parent-legacy-dino-002","source":"pilkwang/rsna-knee-baseline-v1"},"dinosaurs":{"name":"RSNA Knee | DINOsaur V17 TRAIN","output":"rsna-knee-axial-ae-v17","strategy":"denoising axial AE bottleneck + meniscus-only OOF heads"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from pathlib import Path\nimport os, gc, json, math, time, random, hashlib, glob, warnings, zipfile\nfrom concurrent.futures import ThreadPoolExecutor, as_completed\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom sklearn.decomposition import PCA\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.ensemble import HistGradientBoostingClassifier\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.preprocessing import StandardScaler\n\nwarnings.filterwarnings(\"ignore\")\n\nSEED = 20260821\nrandom.seed(SEED)\nnp.random.seed(SEED)\ntorch.manual_seed(SEED)\n\nROOT = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")\nOUT = Path(\"/kaggle/working/rsna-knee-axial-ae-v17-improved\")\nOUT.mkdir(parents=True, exist_ok=True)\n\nLABELS = [\n    \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n    \"Medial OA\", \"Lateral OA\", \"PF OA\", \"Effusion\",\n    \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\",\n]\n\nSIZE = 128\nSLICES_PER_STUDY = 12\nAE_EPOCHS = 15\nAE_BATCH = 96\n\ntrain = pd.read_csv(\n    ROOT / \"train.csv\",\n    dtype={\"StudyInstanceUID\": str},\n)\nseries = pd.read_csv(\n    ROOT / \"train_series.csv\",\n    dtype={\n        \"StudyInstanceUID\": str,\n        \"SeriesInstanceUID\": str,\n    },\n)\nseries = series.loc[:, ~series.columns.duplicated()]\ntrain_ids = train[\"StudyInstanceUID\"].astype(str).tolist()\nprint(f\"train studies: {len(train_ids)}\")\nprint(f\"label columns: {LABELS}\")\nprint(f\"label coverage: {train[LABELS].notna().all(axis=1).sum()}/{len(train)} gold\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def _find_artifact(filename):\n    matches = []\n    root = Path(\"/kaggle/input\")\n    if root.exists():\n        for current, dirs, files in os.walk(root):\n            dirs[:] = [\n                d for d in dirs\n                if d not in (\"train_series\", \"test_series\", \".git\", \"__pycache__\")\n            ]\n            if filename in files:\n                candidate = Path(current) / filename\n                if candidate.is_file():\n                    matches.append(candidate)\n    matches.sort(key=lambda p: (len(p.parts), str(p)))\n    return matches[0] if matches else None\n\n\ndef _series_score(row):\n    score = 0.0\n    try:\n        score += 3.0 * float(pd.to_numeric(row.get(\"Fat_Suppression\", 0), errors=\"coerce\") or 0)\n    except Exception:\n        pass\n    try:\n        score += 2.0 * float(pd.to_numeric(row.get(\"Fluid_Sensitive\", 0), errors=\"coerce\") or 0)\n    except Exception:\n        pass\n    description = (\n        str(row.get(\"SeriesDescription\", \"\"))\n        + \" \"\n        + str(row.get(\"SequenceName\", \"\"))\n    ).lower()\n    if \"fs\" in description or \"fat\" in description or \"stir\" in description:\n        score += 1.5\n    if \"pd\" in description or \"proton\" in description or \"t2\" in description:\n        score += 0.8\n    return score\n\n\ndef _choose_evenly(paths, n):\n    paths = sorted(paths)\n    if len(paths) <= n:\n        return paths\n    idx = np.linspace(0, len(paths) - 1, n).round().astype(int)\n    return [paths[i] for i in idx]\n\n\ndef _study_paths():\n    axial = series[\n        series[\"Anatomical_Plane\"].astype(str).str.lower().eq(\"axial\")\n    ].copy()\n    axial[\"_score\"] = axial.apply(_series_score, axis=1)\n    grouped = {\n        str(uid): frame.sort_values(\"_score\", ascending=False)\n        for uid, frame in axial.groupby(\"StudyInstanceUID\", sort=False)\n    }\n    mapping = {}\n    missing = 0\n    for uid in train_ids:\n        frame = grouped.get(str(uid))\n        if frame is None or frame.empty:\n            mapping[str(uid)] = []\n            missing += 1\n            continue\n        selected = []\n        for _, row in frame.head(2).iterrows():\n            series_uid = str(row[\"SeriesInstanceUID\"])\n            folder = ROOT / \"train_series\" / str(uid) / series_uid\n            selected.extend(folder.glob(\"*.dcm\"))\n        mapping[str(uid)] = _choose_evenly(\n            [str(p) for p in selected],\n            SLICES_PER_STUDY,\n        )\n    print(f\"axial studies={len(train_ids)-missing}/{len(train_ids)} missing={missing}\")\n    return mapping\n\n\nPATHS = _study_paths()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def _read_slice(path):\n    ds = pydicom.dcmread(path)\n    image = ds.pixel_array.astype(np.float32)\n    slope = float(getattr(ds, \"RescaleSlope\", 1.0) or 1.0)\n    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0) or 0.0)\n    image = image * slope + intercept\n    finite = np.isfinite(image)\n    if not finite.any():\n        raise ValueError(\"non-finite image\")\n    lo, hi = np.percentile(image[finite], [1.0, 99.0])\n    if hi <= lo:\n        lo = float(np.min(image[finite]))\n        hi = float(np.max(image[finite]))\n    if hi > lo:\n        image = np.clip((image - lo) / (hi - lo), 0.0, 1.0)\n    else:\n        image = np.zeros_like(image, dtype=np.float32)\n    if str(getattr(ds, \"PhotometricInterpretation\", \"\")).upper() == \"MONOCHROME1\":\n        image = 1.0 - image\n    image = cv2.resize(image, (SIZE, SIZE), interpolation=cv2.INTER_AREA)\n    return np.clip(image * 255.0, 0, 255).astype(np.uint8)\n\n\ncache_path = Path(\"/kaggle/working/v17_improved_axial_uint8.npy\")\nowner_path = Path(\"/kaggle/working/v17_improved_axial_owner.npy\")\n\nif cache_path.exists() and owner_path.exists():\n    ALL_SLICES = np.load(cache_path, mmap_mode=\"r\")\n    OWNERS = np.load(owner_path)\n    print(\"loaded axial cache\", ALL_SLICES.shape)\nelse:\n    jobs = []\n    for study_index, uid in enumerate(train_ids):\n        for path in PATHS[str(uid)]:\n            jobs.append((study_index, path))\n    arrays = []\n    owners = []\n    t0 = time.time()\n\n    def worker(item):\n        study_index, path = item\n        try:\n            return study_index, _read_slice(path)\n        except Exception:\n            return None\n\n    with ThreadPoolExecutor(max_workers=12) as pool:\n        futures = [pool.submit(worker, item) for item in jobs]\n        for completed, future in enumerate(as_completed(futures), 1):\n            result = future.result()\n            if result is not None:\n                study_index, image = result\n                arrays.append(image)\n                owners.append(study_index)\n            if completed % 4000 == 0:\n                elapsed = time.time() - t0\n                print(f\"read {completed}/{len(futures)} {elapsed/60:.1f}m\", flush=True)\n\n    ALL_SLICES = np.stack(arrays).astype(np.uint8)\n    OWNERS = np.asarray(owners, dtype=np.int32)\n    np.save(cache_path, ALL_SLICES)\n    np.save(owner_path, OWNERS)\n    print(\"saved axial cache\", ALL_SLICES.shape)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class SEBlock(nn.Module):\n    \"\"\"Squeeze-and-Excitation: channel attention.\"\"\"\n    def __init__(self, channels, reduction=4):\n        super().__init__()\n        self.pool = nn.AdaptiveAvgPool2d(1)\n        self.fc = nn.Sequential(\n            nn.Linear(channels, channels // reduction, bias=False),\n            nn.ReLU(inplace=True),\n            nn.Linear(channels // reduction, channels, bias=False),\n            nn.Sigmoid(),\n        )\n\n    def forward(self, x):\n        b, c, _, _ = x.size()\n        w = self.pool(x).view(b, c)\n        w = self.fc(w).view(b, c, 1, 1)\n        return x * w\n\n\nclass ResidualBlock(nn.Module):\n    def __init__(self, channels):\n        super().__init__()\n        self.block = nn.Sequential(\n            nn.Conv2d(channels, channels, 3, padding=1),\n            nn.BatchNorm2d(channels),\n            nn.LeakyReLU(0.2, inplace=True),\n            nn.Conv2d(channels, channels, 3, padding=1),\n            nn.BatchNorm2d(channels),\n        )\n        self.se = SEBlock(channels)\n\n    def forward(self, x):\n        return x + self.se(self.block(x))\n\n\nclass AxialDenoisingAE(nn.Module):\n    \"\"\"\n    Improved denoising autoencoder with:\n    - Squeeze-and-Excitation attention in residual blocks\n    - Deeper bottleneck (4 residual blocks instead of 3)\n    - 128 x 8x8 = 8192 bottleneck activations\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.init_conv = nn.Sequential(\n            nn.Conv2d(1, 32, 3, padding=1),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.enc1 = nn.Sequential(\n            nn.Conv2d(32, 64, 3, stride=2, padding=1),\n            nn.BatchNorm2d(64),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.res1 = ResidualBlock(64)\n        self.enc2 = nn.Sequential(\n            nn.Conv2d(64, 128, 3, stride=2, padding=1),\n            nn.BatchNorm2d(128),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.res2 = ResidualBlock(128)\n        self.enc3 = nn.Sequential(\n            nn.Conv2d(128, 256, 3, stride=2, padding=1),\n            nn.BatchNorm2d(256),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.res3 = ResidualBlock(256)\n        self.enc4 = nn.Sequential(\n            nn.Conv2d(256, 128, 3, stride=2, padding=1),\n            nn.BatchNorm2d(128),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.res4 = ResidualBlock(128)\n\n        self.dec4 = nn.Sequential(\n            nn.ConvTranspose2d(128, 256, 4, stride=2, padding=1),\n            nn.BatchNorm2d(256),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.dec3 = nn.Sequential(\n            nn.ConvTranspose2d(256, 128, 4, stride=2, padding=1),\n            nn.BatchNorm2d(128),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.dec2 = nn.Sequential(\n            nn.ConvTranspose2d(128, 64, 4, stride=2, padding=1),\n            nn.BatchNorm2d(64),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.dec1 = nn.Sequential(\n            nn.ConvTranspose2d(64, 32, 4, stride=2, padding=1),\n            nn.BatchNorm2d(32),\n            nn.LeakyReLU(0.2, inplace=True),\n        )\n        self.final = nn.Sequential(\n            nn.Conv2d(32, 1, 3, padding=1),\n            nn.Sigmoid(),\n        )\n\n    def encode(self, x):\n        x = self.init_conv(x)\n        x = self.res1(self.enc1(x))\n        x = self.res2(self.enc2(x))\n        x = self.res3(self.enc3(x))\n        x = self.res4(self.enc4(x))\n        return x\n\n    def decode(self, z):\n        x = self.dec4(z)\n        x = self.dec3(x)\n        x = self.dec2(x)\n        x = self.dec1(x)\n        return self.final(x)\n\n    def forward(self, x):\n        return self.decode(self.encode(x))\n\n\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nmodel = AxialDenoisingAE().to(device)\nn_params = sum(p.numel() for p in model.parameters())\nprint(f\"model params: {n_params:,}\")\nprint(f\"bottleneck: 128 x 8x8 = {128*8*8} activations\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"rng = np.random.default_rng(SEED)\nstudy_order = np.arange(len(train_ids))\nrng.shuffle(study_order)\nn_val = max(1, int(0.10 * len(study_order)))\nval_studies = set(study_order[:n_val].tolist())\n\nis_val = np.asarray([owner in val_studies for owner in OWNERS])\ntrain_idx = np.flatnonzero(~is_val)\nval_idx = np.flatnonzero(is_val)\n\noptimizer = torch.optim.AdamW(model.parameters(), lr=3e-4, weight_decay=1e-4)\nscheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=AE_EPOCHS, eta_min=1e-5)\n\nbest_val = float(\"inf\")\nbest_state = None\n\n\ndef _batch(indices):\n    array = np.asarray(ALL_SLICES[indices], dtype=np.float32) / 255.0\n    return torch.from_numpy(array[:, None]).to(device)\n\n\nfor epoch in range(AE_EPOCHS):\n    rng = np.random.default_rng(SEED + epoch)\n    order = train_idx.copy()\n    rng.shuffle(order)\n\n    model.train()\n    total = 0.0\n    n = 0\n\n    for start in range(0, len(order), AE_BATCH):\n        idx = order[start:start + AE_BATCH]\n        if len(idx) < 8:\n            continue\n        clean = _batch(idx)\n        noise_scale = 0.02 + 0.02 * (epoch / max(AE_EPOCHS - 1, 1))\n        noisy = clean + noise_scale * torch.randn_like(clean)\n        pixel_keep = (torch.rand_like(noisy) > 0.08).float()\n        noisy = torch.clamp(noisy * pixel_keep, 0.0, 1.0)\n\n        optimizer.zero_grad(set_to_none=True)\n        pred = model(noisy)\n\n        mse = F.mse_loss(pred, clean)\n        l1 = F.l1_loss(pred, clean)\n        loss = mse + 0.08 * l1\n\n        loss.backward()\n        torch.nn.utils.clip_grad_norm_(model.parameters(), 2.0)\n        optimizer.step()\n\n        total += float(loss.detach().cpu()) * len(idx)\n        n += len(idx)\n\n    scheduler.step()\n\n    model.eval()\n    val_loss = 0.0\n    val_n = 0\n    with torch.no_grad():\n        for start in range(0, len(val_idx), AE_BATCH):\n            idx = val_idx[start:start + AE_BATCH]\n            if len(idx) == 0:\n                continue\n            clean = _batch(idx)\n            pred = model(clean)\n            loss = F.mse_loss(pred, clean)\n            val_loss += float(loss.cpu()) * len(idx)\n            val_n += len(idx)\n\n    train_loss = total / max(n, 1)\n    val_loss = val_loss / max(val_n, 1)\n    lr_now = scheduler.get_last_lr()[0]\n    print(f\"AE epoch {epoch+1:02d}/{AE_EPOCHS} train={train_loss:.6f} val_mse={val_loss:.6f} lr={lr_now:.2e}\")\n\n    if val_loss < best_val:\n        best_val = val_loss\n        best_state = {\n            key: value.detach().cpu().clone()\n            for key, value in model.state_dict().items()\n        }\n\nif best_state is None:\n    raise RuntimeError(\"AE best checkpoint missing\")\n\nmodel.load_state_dict(best_state, strict=True)\nmodel.eval()\n\nencoder_path = OUT / \"v17_improved_axial_ae_encoder.pt\"\ntorch.save(\n    {\n        \"version\": \"17-improved\",\n        \"size\": SIZE,\n        \"state_dict\": best_state,\n        \"best_val_mse\": float(best_val),\n        \"n_params\": n_params,\n    },\n    encoder_path,\n)\nprint(f\"\\nsaved encoder: {encoder_path}\")\nprint(f\"best val MSE: {best_val:.6f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"@torch.no_grad()\ndef _slice_feature(images_uint8):\n    x = torch.from_numpy(\n        images_uint8.astype(np.float32)[:, None] / 255.0\n    ).to(device)\n\n    z = model.encode(x)\n\n    pooled_4x4 = F.adaptive_avg_pool2d(z, (4, 4)).flatten(1)\n    pooled_2x2 = F.adaptive_avg_pool2d(z, (2, 2)).flatten(1)\n    gmax = z.amax(dim=(2, 3))\n    recon = model.decode(z)\n    err = ((recon - x) ** 2).mean(dim=(1, 2, 3), keepdim=False).unsqueeze(1)\n    h, w = recon.shape[2], recon.shape[3]\n    qh, qw = h // 2, w // 2\n    q_errors = []\n    for y0, y1 in [(0, qh), (qh, h)]:\n        for x0, x1 in [(0, qw), (qw, w)]:\n            qe = ((recon[:, :, y0:y1, x0:x1] - x[:, :, y0:y1, x0:x1]) ** 2).mean(dim=(1, 2, 3))\n            q_errors.append(qe.unsqueeze(1))\n    quadrant_err = torch.cat(q_errors, dim=1)\n    z_mean = z.mean(dim=(2, 3))\n    z_std = z.std(dim=(2, 3))\n\n    return torch.cat(\n        [pooled_4x4, pooled_2x2, gmax, err, quadrant_err, z_mean, z_std],\n        dim=1,\n    ).cpu().numpy().astype(np.float32)\n\n\nFEATURE_DIM = None\nstudy_features = []\n\nfor study_index in range(len(train_ids)):\n    idx = np.flatnonzero(OWNERS == study_index)\n\n    if len(idx) == 0:\n        study_features.append(None)\n        continue\n\n    per_slice = []\n    for start in range(0, len(idx), 64):\n        feat = _slice_feature(np.asarray(ALL_SLICES[idx[start:start + 64]]))\n        per_slice.append(feat)\n\n    per_slice = np.concatenate(per_slice, axis=0)\n\n    aggregate = np.concatenate([\n        per_slice.mean(axis=0),\n        per_slice.std(axis=0),\n        per_slice.max(axis=0),\n    ]).astype(np.float32)\n\n    if FEATURE_DIM is None:\n        FEATURE_DIM = aggregate.shape[0]\n        print(f\"feature dim: {FEATURE_DIM}\")\n\n    study_features.append(aggregate)\n\n    if (study_index + 1) % 500 == 0:\n        print(f\"study features {study_index+1}/{len(train_ids)}\", flush=True)\n\nstudy_features = [\n    f if f is not None else np.zeros(FEATURE_DIM, np.float32)\n    for f in study_features\n]\n\nFEATURES = np.stack(study_features)\nprint(\"study feature matrix\", FEATURES.shape)\n\nfeature_path = OUT / \"v17_improved_train_study_features.npy\"\nnp.save(feature_path, FEATURES)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"oof_path = _find_artifact(\"oof.npz\")\nfold_path = _find_artifact(\"v52_e11_oof.csv\")\npublic_path = _find_artifact(\"v52_oof.csv\")\n\nif oof_path is None or fold_path is None or public_path is None:\n    raise FileNotFoundError(\"V17 needs oof.npz + v52_oof.csv + v52_e11_oof.csv\")\n\nwith np.load(oof_path, allow_pickle=False) as data:\n    ids = data[\"ids\"].astype(str)\n    targets = data[\"targets\"].astype(str).tolist()\n    y = data[\"y_derived\"].astype(np.float64)\n    base_oof = data[\"pred\"].astype(np.float64)\n    gold = data[\"gold_mask\"].astype(bool)\n\nif targets != LABELS:\n    raise RuntimeError(\"target contract mismatch\")\nif not np.array_equal(ids, np.asarray(train_ids, dtype=str)):\n    raise RuntimeError(\"train ID contract mismatch\")\n\nexact = train[LABELS].apply(pd.to_numeric, errors=\"coerce\")\nofficial_gold = exact.notna().all(axis=1).to_numpy()\n\nif not np.array_equal(gold, official_gold):\n    raise RuntimeError(\"gold contract mismatch\")\n\ny = np.where(np.isfinite(y), y, 0.5)\ny[gold] = exact.loc[gold, LABELS].to_numpy(np.float64)\n\nfold_frame = pd.read_csv(fold_path, dtype={\"StudyInstanceUID\": str})\nfold_frame = train[[\"StudyInstanceUID\"]].merge(\n    fold_frame[[\"StudyInstanceUID\", \"fold\"] + LABELS],\n    on=\"StudyInstanceUID\", how=\"left\", validate=\"one_to_one\",\n)\nfold = pd.to_numeric(fold_frame[\"fold\"], errors=\"coerce\").to_numpy()\nif not np.isfinite(fold).all():\n    raise RuntimeError(\"fold assignment incomplete\")\nfold = fold.astype(np.int64)\n\npublic_frame = pd.read_csv(public_path, dtype={\"StudyInstanceUID\": str})\npublic_frame = train[[\"StudyInstanceUID\"]].merge(\n    public_frame[[\"StudyInstanceUID\"] + LABELS],\n    on=\"StudyInstanceUID\", how=\"left\", validate=\"one_to_one\",\n)\n\n\ndef rank_cols(values):\n    return pd.DataFrame(values).rank(method=\"average\", pct=True).to_numpy(np.float64)\n\n\nbase_rank = rank_cols(base_oof)\npublic_rank = rank_cols(public_frame[LABELS].to_numpy(np.float64))\npass2_rank = rank_cols(fold_frame[LABELS].to_numpy(np.float64))\noof_anchor = rank_cols(\n    0.85 * rank_cols(0.50 * base_rank + 0.50 * public_rank)\n    + 0.15 * pass2_rank\n)\n\nprint(f\"oof anchor shape: {oof_anchor.shape}\")\nprint(f\"gold labels: {gold.sum()}/{len(gold)}\")\nprint(f\"folds: {sorted(np.unique(fold).tolist())}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"scaler = StandardScaler()\nx_scaled = scaler.fit_transform(FEATURES)\n\npca = PCA(n_components=192, random_state=SEED)\nX = pca.fit_transform(x_scaled).astype(np.float64)\n\nprint(f\"PCA explained variance: {float(pca.explained_variance_ratio_.sum()):.6f}\")\nprint(f\"PCA output shape: {X.shape}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def target_auc(truth, prediction):\n    if len(np.unique(truth)) < 2:\n        return np.nan\n    return float(roc_auc_score(truth, prediction))\n\n\ndef bootstrap_support(truth, anchor, candidate, seed, n_boot=400):\n    pos = np.flatnonzero(truth > 0.5)\n    neg = np.flatnonzero(truth <= 0.5)\n    if len(pos) < 2 or len(neg) < 2:\n        return 0.0\n    rng = np.random.default_rng(seed)\n    gains = []\n    for _ in range(n_boot):\n        idx = np.concatenate([\n            rng.choice(pos, len(pos), replace=True),\n            rng.choice(neg, len(neg), replace=True),\n        ])\n        a0 = target_auc(truth[idx], anchor[idx])\n        a1 = target_auc(truth[idx], candidate[idx])\n        if np.isfinite(a0) and np.isfinite(a1):\n            gains.append(a1 - a0)\n    return float((np.asarray(gains) > 0).mean()) if gains else 0.0\n\n\ndef fit_target_lr(target_index, train_mask):\n    target_y = y[:, target_index]\n    weak_high = (~gold) & ((target_y <= 0.15) | (target_y >= 0.85))\n    use = train_mask & (gold | weak_high)\n    labels = (target_y[use] >= 0.5).astype(int)\n    confidence = np.clip(np.abs(target_y[use] - 0.5) * 2.0, 0.25, 1.0)\n    sample_weight = 0.75 + 1.25 * confidence\n    sample_weight[gold[use]] = 12.0\n\n    models = []\n    for c_value in (0.025, 0.08, 0.20):\n        clf = LogisticRegression(\n            C=float(c_value), solver=\"liblinear\", penalty=\"l2\",\n            class_weight=\"balanced\", max_iter=3000,\n            random_state=SEED + target_index,\n        )\n        clf.fit(X[use], labels, sample_weight=sample_weight)\n        models.append(clf)\n    return models\n\n\ndef fit_target_hgb(target_index, train_mask):\n    target_y = y[:, target_index]\n    weak_high = (~gold) & ((target_y <= 0.15) | (target_y >= 0.85))\n    use = train_mask & (gold | weak_high)\n    labels = (target_y[use] >= 0.5).astype(int)\n    confidence = np.clip(np.abs(target_y[use] - 0.5) * 2.0, 0.25, 1.0)\n    sample_weight = 0.75 + 1.25 * confidence\n    sample_weight[gold[use]] = 12.0\n\n    models = []\n    for n_iter, lr, depth in [(128, 0.05, 5), (200, 0.03, 6), (64, 0.08, 4)]:\n        clf = HistGradientBoostingClassifier(\n            max_iter=n_iter, learning_rate=lr, max_depth=depth,\n            min_samples_leaf=20, l2_regularization=1.0,\n            random_state=SEED + target_index,\n        )\n        clf.fit(X[use], labels, sample_weight=sample_weight)\n        models.append(clf)\n    return models\n\n\ndef predict_ensemble(lr_models, hgb_models, indices):\n    lr_pred = np.stack([m.predict_proba(X[indices])[:, 1] for m in lr_models], axis=0).mean(axis=0)\n    hgb_pred = np.stack([m.predict_proba(X[indices])[:, 1] for m in hgb_models], axis=0).mean(axis=0)\n    return 0.5 * lr_pred + 0.5 * hgb_pred\n\n\ngold_index = np.flatnonzero(gold)\nunique_folds = sorted(int(v) for v in np.unique(fold[gold]) if int(v) >= 0)\n\noof_special = np.full((len(train), len(LABELS)), np.nan, np.float64)\n\nfor target_pos, target in enumerate(LABELS):\n    target_index = LABELS.index(target)\n    t0 = time.time()\n\n    for outer in unique_folds:\n        valid = gold & (fold == outer)\n        train_mask = np.ones(len(train), dtype=bool)\n        train_mask[valid] = False\n\n        lr_models = fit_target_lr(target_index, train_mask)\n        hgb_models = fit_target_hgb(target_index, train_mask)\n        valid_idx = np.flatnonzero(valid)\n\n        oof_special[valid_idx, target_pos] = predict_ensemble(\n            lr_models, hgb_models, valid_idx\n        )\n\n    elapsed = time.time() - t0\n    valid_auc = target_auc(y[gold_index, target_index], oof_special[gold_index, target_pos])\n    anchor_auc = target_auc(y[gold_index, target_index], oof_anchor[gold_index, target_index])\n    print(f\"{target:20s} OOF AUC={valid_auc:.4f} anchor={anchor_auc:.4f} delta={valid_auc-anchor_auc:+.4f} [{elapsed:.0f}s]\")\n\nif not np.isfinite(oof_special[gold_index]).all():\n    raise RuntimeError(\"V17 OOF incomplete\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def nested_deploy(truth, anchor, specialist, folds, seed):\n    specialist_rank = pd.Series(specialist).rank(method=\"average\", pct=True).to_numpy(np.float64)\n    grid = (0.00, 0.03, 0.05, 0.08, 0.10, 0.12)\n    nested = np.full(len(truth), np.nan, np.float64)\n\n    for outer in sorted(np.unique(folds)):\n        choose = folds != outer\n        valid = folds == outer\n        best_weight = 0.0\n        best_obj = target_auc(truth[choose], anchor[choose])\n        for weight in grid[1:]:\n            candidate = (1.0 - weight) * anchor[choose] + weight * specialist_rank[choose]\n            auc = target_auc(truth[choose], candidate)\n            objective = auc - 0.0025 * weight\n            if np.isfinite(objective) and objective > best_obj:\n                best_obj = objective\n                best_weight = weight\n        nested[valid] = (1.0 - best_weight) * anchor[valid] + best_weight * specialist_rank[valid]\n\n    if not np.isfinite(nested).all():\n        return 0.0, 0.0, 0.0\n\n    base_auc = target_auc(truth, anchor)\n    nested_auc = target_auc(truth, nested)\n    gain = nested_auc - base_auc\n    support = bootstrap_support(truth, anchor, nested, seed)\n\n    best_weight = 0.0\n    best_obj = base_auc\n    for weight in grid[1:]:\n        candidate = (1.0 - weight) * anchor + weight * specialist_rank\n        auc = target_auc(truth, candidate)\n        objective = auc - 0.0025 * weight\n        if np.isfinite(objective) and objective > best_obj:\n            best_obj = objective\n            best_weight = weight\n\n    if gain >= 0.006 and support >= 0.66:\n        deploy = min(best_weight, 0.12)\n    elif gain >= 0.003 and support >= 0.60:\n        deploy = min(best_weight, 0.06)\n    else:\n        deploy = 0.0\n\n    return float(deploy), float(gain), float(support)\n\n\ndeploy_weights = {}\nvalidation = {}\nfinal_lr_coefficients = {}\nfinal_hgb_models = {}\n\nfor target_pos, target in enumerate(LABELS):\n    target_index = LABELS.index(target)\n    truth = y[gold_index, target_index]\n    anchor = oof_anchor[gold_index, target_index]\n    specialist = oof_special[gold_index, target_pos]\n\n    deploy, gain, support = nested_deploy(\n        truth, anchor, specialist, fold[gold_index], SEED + target_index,\n    )\n\n    deploy_weights[target] = deploy\n    validation[target] = {\n        \"nested_gain\": gain,\n        \"bootstrap_support\": support,\n        \"specialist_auc\": target_auc(truth, specialist),\n        \"anchor_auc\": target_auc(truth, anchor),\n        \"deploy_weight\": deploy,\n    }\n\n    all_mask = np.ones(len(train), dtype=bool)\n    lr_final = fit_target_lr(target_index, all_mask)\n    hgb_final = fit_target_hgb(target_index, all_mask)\n\n    final_lr_coefficients[target] = {\n        \"coef\": np.stack([m.coef_[0] for m in lr_final]).astype(np.float32),\n        \"intercept\": np.asarray([m.intercept_[0] for m in lr_final], np.float32),\n    }\n    final_hgb_models[target] = hgb_final\n\nprint(\"\\n\" + \"=\" * 80)\nprint(\"DEPLOY DECISIONS\")\nprint(\"=\" * 80)\nn_deployed = 0\nfor target in LABELS:\n    info = validation[target]\n    status = \"DEPLOY\" if deploy_weights[target] > 0.01 else \"skip\"\n    if deploy_weights[target] > 0.01:\n        n_deployed += 1\n    print(\n        f\"{target:20s} w={deploy_weights[target]:.3f} \"\n        f\"gain={info['nested_gain']:+.4f} sup={info['bootstrap_support']:.3f} \"\n        f\"sp_auc={info['specialist_auc']:.4f} anc={info['anchor_auc']:.4f} [{status}]\"\n    )\nprint(f\"\\ntotal deployed: {n_deployed}/{len(LABELS)}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"npz_data = {\n    \"scaler_mean\": scaler.mean_.astype(np.float32),\n    \"scaler_scale\": scaler.scale_.astype(np.float32),\n    \"pca_mean\": pca.mean_.astype(np.float32),\n    \"pca_components\": pca.components_.astype(np.float32),\n}\nfor t in LABELS:\n    safe = t.replace(\" \", \"_\").replace(\"'\", \"\").lower()\n    npz_data[f\"{safe}_coef\"] = final_lr_coefficients[t][\"coef\"]\n    npz_data[f\"{safe}_intercept\"] = final_lr_coefficients[t][\"intercept\"]\nnp.savez_compressed(OUT / \"v17_improved_projection_and_heads.npz\", **npz_data)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"manifest = {\n    \"artifact_type\": \"rsna-knee-axial-ae-v17-improved\",\n    \"version\": \"17-improved\",\n    \"targets\": LABELS,\n    \"image_size\": SIZE,\n    \"slices_per_study\": SLICES_PER_STUDY,\n    \"feature_layout\": (\n        \"per-slice: 4x4 avg pool + 2x2 avg pool + channel max + \"\n        \"recon MSE + 4-quadrant MSE + latent mean + latent std; \"\n        \"study: mean/std/max; StandardScaler + PCA192\"\n    ),\n    \"classifier\": \"LogisticRegression(3x) + HistGradientBoostingClassifier(3x) ensemble\",\n    \"deploy_weights\": deploy_weights,\n    \"validation\": validation,\n    \"inference_shrink\": 0.90,\n    \"best_reconstruction_val_mse\": float(best_val),\n    \"n_ae_params\": n_params,\n}\n\n(OUT / \"v17_improved_manifest.json\").write_text(\n    json.dumps(manifest, indent=2), encoding=\"utf-8\",\n)\n\nimport pickle\nfor target in LABELS:\n    safe = target.replace(\" \", \"_\").replace(\"'\", \"\").lower()\n    with open(OUT / f\"hgb_{safe}.pkl\", \"wb\") as f:\n        pickle.dump(final_hgb_models[target], f)\n\noof_frame = train[[\"StudyInstanceUID\"]].copy()\nfor target_pos, target in enumerate(LABELS):\n    oof_frame[target] = oof_special[:, target_pos]\noof_frame[\"gold\"] = gold.astype(np.uint8)\noof_frame[\"fold\"] = fold\noof_frame.to_csv(OUT / \"v17_improved_oof.csv\", index=False)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"try:\n    feature_path.unlink()\nexcept OSError:\n    pass\n\nzip_path = Path(\"/kaggle/working/rsna-knee-axial-ae-v17-improved.zip\")\nwith zipfile.ZipFile(zip_path, \"w\", compression=zipfile.ZIP_DEFLATED) as archive:\n    for path in sorted(OUT.iterdir()):\n        if path.is_file():\n            archive.write(path, arcname=path.name)\n\nprint(\"\\nV17 IMPROVED deployment:\")\nfor target in LABELS:\n    info = validation[target]\n    print(\n        f\"  {target:18s} \"\n        f\"weight={deploy_weights[target]:.3f} \"\n        f\"gain={info['nested_gain']:+.4f} \"\n        f\"support={info['bootstrap_support']:.3f} \"\n        f\"specialist_auc={info['specialist_auc']:.4f}\"\n    )\nprint(f\"\\nweights folder: {OUT}\")\nprint(f\"weights zip: {zip_path}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}