{"cells":[{"cell_type":"markdown","id":"89092460","metadata":{},"source":"# RSNA 2026 Knee Abnormality Detection — Advanced Production Submission Notebook (v2)\n**Multimodal Learning (MRI + Radiology Reports) with DINOv3 Backbone, Temporal 1D ConvNeXt Aggregation, Asymmetric Focal Loss & Rank-Averaged Ensembling**\n\n---\n\n### Key System & Architectural Highlights (v2):\n1. **High-Throughput Decoupled Ingestion Pipeline (3x Speedup)**:\n   - Multi-threaded CPU decoding (`max_workers=10`) for DICOM pixel reading.\n   - Sequential main-thread GPU forward passes (`batch_size=64`) for full T4 VRAM utilization without CUDA thread stalls.\n   - Unbuffered progress logs every 100 studies (`flush=True`).\n2. **Temporal 1D ConvNeXt Sequence Aggregator**:\n   - 1D depth convolutions ($k=3$) to explicitly capture multi-slice continuity (meniscal tears & ACL ruptures spanning 2–3 contiguous slices).\n3. **Asymmetric Focal Loss for Imbalanced Labels**:\n   - Multi-label loss with negative suppression ($\\gamma_{neg} = 2.0, \\gamma_{pos} = 0.0$) preventing easy negative dominance on rare classes like Fracture ($6.9\\%$).\n4. **Clinical DICOM Pixel Handling & Geometry**:\n   - `MONOCHROME1` contrast inversion & `RescaleSlope`/`RescaleIntercept` intensity scaling.\n   - Laterality normalization mapping Right knees onto Left-knee convention (horizontal mirroring for Coronal/Axial, stack reversal for Sagittal).\n   - Header-only physical z-position extraction (`stop_before_pixels=True`) yielding **20x speedups** in slice ordering.\n5. **Rank-Averaged Multi-Model Ensemble**: Closed-form Ridge Probe + PyTorch `TemporalConvNeXtHead` + Fine-Tuned DINOv3 blocks.\n"},{"cell_type":"code","execution_count":null,"id":"8343bfc0","metadata":{},"outputs":[],"source":"# Cell 1: Environment Setup & Package Installation\nimport os\nimport sys\nimport subprocess\n\n# Package Verification\ntry:\n    import transformers\n    import pydicom\n    import cv2\nexcept ImportError:\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"-U\", \"transformers\", \"huggingface_hub\", \"pydicom\", \"opencv-python-headless\", \"scikit-learn\"])\n\nimport time\nimport math\nimport glob\nimport re\nimport gc\nimport json\nimport warnings\nimport unicodedata\nfrom pathlib import Path\nfrom concurrent.futures import ThreadPoolExecutor\n\nimport numpy as np\nimport pandas as pd\nfrom scipy.stats import rankdata\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.preprocessing import StandardScaler\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# No HuggingFace client is used. The backbone is attached as a Kaggle model input\n# and loaded from disk (Cell 6). Competition reruns execute with internet\n# disabled, so a hub download cannot succeed at scoring time; the previous\n# version spent about two minutes in DNS retries before failing. Removing the\n# dependency also removes the API token that was hardcoded in this cell.\n\n# Abort in seconds rather than in hours when the accelerator is missing.\n#\n# Version 10 requested a T4, reported \"Running on device: cpu\", and still spent\n# 4 h 25 m encoding 4,407 studies before anyone could see the problem. The cause\n# was not quota: it was the environment pin in kernel-metadata.json, which named\n# gcr.io/kaggle-images/python -- the CPU-only image. enable_gpu allocates the\n# hardware, but a CPU-only container has a CPU-only torch build, so\n# torch.cuda.is_available() is False no matter what hardware is attached. The\n# GPU image is gcr.io/kaggle-gpu-images/python; leaving docker_image unset lets\n# Kaggle choose correctly from enable_gpu.\n#\n# Set REQUIRE_GPU = False to take the CPU path deliberately. It does finish and\n# it does score -- Cell 7b drops the input size to keep it inside the budget --\n# but at 224 px with a ViT-S rather than 336 px with a ViT-B.\nREQUIRE_GPU = True\n\nDEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nprint(f\"Running on device: {DEVICE}\", flush=True)\nif DEVICE == \"cuda\":\n    print(f\"GPU Model: {torch.cuda.get_device_name(0)}\", flush=True)\n    print(f\"VRAM Available: {torch.cuda.get_device_properties(0).total_memory / 1e9:.2f} GB\", flush=True)\n    print(f\"CUDA devices visible: {torch.cuda.device_count()}\", flush=True)\nelif REQUIRE_GPU:\n    raise RuntimeError(\n        \"No CUDA device is visible, so this run would take about 4.5 hours \"\n        \"instead of about 40 minutes and would produce weaker features.\\n\"\n        \"  Checked first: torch \" + torch.__version__ + \" reports \"\n        \"torch.cuda.is_available() == False.\\n\"\n        \"  Most likely cause: kernel-metadata.json pins docker_image to \"\n        \"gcr.io/kaggle-images/python, which is the CPU-only image. Remove the \"\n        \"docker_image key entirely, or point it at \"\n        \"gcr.io/kaggle-gpu-images/python, and push again.\\n\"\n        \"  Also worth checking: the accelerator dropdown in the notebook editor \"\n        \"(a Quick Save uses the editor's setting, not the pushed metadata), and \"\n        \"GPU quota at kaggle.com/settings.\\n\"\n        \"  To run on CPU on purpose, set REQUIRE_GPU = False in this cell.\")\n"},{"cell_type":"code","execution_count":null,"id":"e96f5446","metadata":{},"outputs":[],"source":"# Cell 2: Paths, Configuration & Hyperparameters\nSEED = 20260805\nnp.random.seed(SEED)\ntorch.manual_seed(SEED)\n\nLABELS = [\n    \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n    \"Medial OA\", \"Lateral OA\", \"PF OA\",\n    \"Effusion\", \"Synovitis\", \"Baker's\",\n    \"Contusion\", \"Fracture\"\n]\nN_LABELS = len(LABELS)\nFOV_MM = 150.0  # Physical field of view\n\n# Slot specification: (Name, Plane, Fluid_Flag, Slice_Budget)\nSLOT_SPEC = [\n    (\"Sagittal_Fluid\", \"Sagittal\", 1, 12),\n    (\"Sagittal_Structural\", \"Sagittal\", 0, 12),\n    (\"Coronal_Fluid\", \"Coronal\", 1, 8),\n    (\"Coronal_Structural\", \"Coronal\", 0, 6),\n    (\"Axial_Fluid\", \"Axial\", 1, 8),\n]\nSLOT_NAMES = [s[0] for s in SLOT_SPEC]\nN_SLOTS = len(SLOT_SPEC)\n\nPLANE_WINDOW = {\n    \"Sagittal\": (0.12, 0.88),\n    \"Coronal\":  (0.20, 0.80),\n    \"Axial\":    (0.12, 0.88),\n}\n\n# Image & Model Parameters\nIMG_SIZE = 336\nPATCH_PX = 16\nGPU_BATCH_SIZE = 64  # per-device encode batch; scaled by the device count\n# Version 12 reported \"CUDA devices visible: 2\" and used one of them. A\n# ViT-B/16 costs about three times a ViT-S/14 per study, so the second card is\n# worth having now. The encode probe in Cell 7b measures whatever this\n# actually delivers before the full ingest commits to it.\nUSE_DATA_PARALLEL = True\n# Kaggle drops a model_sources entry it cannot resolve without erroring, so a\n# wrong version number looks exactly like a successful push and then costs a\n# full ingest against the wrong backbone. Versions 10 and 12 both lost their\n# DINOv3 mount this way. Fail in the first minute instead.\nREQUIRE_BACKBONE_FAMILY = True\nCPU_WORKERS = 10     # Parallel CPU threads for DICOM loading\n\n# ---- Vision backbone (Kaggle model input, never downloaded) ----------------\n# Attach via Add Input -> Models. The mount lands at\n# /kaggle/input/<name>/... or /kaggle/input/models/<owner>/<name>/<framework>/\n# <variation>/<version>, and Cell 6 discovers it by walking for a config.json\n# next to a weights file, so the exact layout does not have to be hardcoded.\n#\n# BACKBONE_FAMILY filters which mount to use when several are attached.\n# BACKBONE_VARIANT_ORDER breaks the tie between variations of the same model.\n# Bigger is better per image and costs more per study; the encode probe in Cell\n# 7b measures the actual cost and downshifts resolution if the run will not fit.\nBACKBONE_FAMILY = \"dinov3\"\n\n# Which size to prefer when several mounts of that family are attached, given as\n# the transformer hidden size read from each checkpoint's config.json:\n#   384 = ViT-S    768 = ViT-B    1024 = ViT-L    1280 = ViT-H\n# Ranking on the mount's path name was tried first and is unreliable: a repo\n# called \"dinov3-vit-large\" advertises its size in the folder name while the\n# ViT-B mount does not, so the larger model won on a name match alone.\n#\n# 768 (ViT-B) is the default because this runs on a single T4 across 4,407\n# studies and ViT-L is roughly 3.5x the compute. Cell 7b measures the real cost\n# and warns if the run will not fit. Set 1024 to prefer ViT-L.\nBACKBONE_TARGET_HIDDEN = 768\n\n# Head Hyperparameters\nSVD_COMPONENTS = 512\nHEAD_HIDDEN = 256\nHEAD_DROPOUT = 0.2\nHEAD_LR = 2.5e-3\nHEAD_WD = 1e-2\nHEAD_EPOCHS = 18\nHEAD_SEEDS = [20260805, 20260806, 20260807]\nGOLD_WEIGHT = 8.0\n# BLEND_W is gone: a fixed 50/50 average shipped a submission worse than\n# its own best component on the last run. Cell 11 now picks the weights by\n# greedy forward selection on the out-of-fold predictions instead.\nENSEMBLE_ROUNDS = 12\nN_FOLDS = 5\n\n# ---- Weak-label source -----------------------------------------------------\n# The rule extractor is the floor. A stronger label table, if one is mounted as a\n# dataset input, overwrites it; the 58 gold studies overwrite whatever wins.\n#\n# Measured on those 58 gold studies (companion notebook, same corpus):\n#   rule extractor      agreement 0.7960  precision 0.6802  recall 0.7648\n#   LLM-extracted       agreement 0.8233  precision 0.6910  recall 0.8602\n#   paired delta +0.0273 agreement, 95% CI [+0.0029, +0.0532], McNemar p=0.0204\n#\n# COMPLIANCE: produce any such table with an open-weights model run locally or\n# inside the kernel. Public write-ups of this competition state the rules forbid\n# sending report text to a hosted LLM API -- confirm that against the official\n# rules page before mounting a table produced that way. This notebook only reads\n# a file you supply; it never calls an external API.\nLABEL_SOURCE = \"auto\"  # \"auto\" | \"file\" | \"rule_inline\"\nLABEL_FILE_PATTERNS = [\"labels_llm_*.parquet\", \"labels_llm_*.csv\",\n                       \"labels_rule_*.parquet\", \"labels_rule_*.csv\"]\n\n# ---- Fold grouping ---------------------------------------------------------\n# PatientID is 1:1 with StudyInstanceUID across all 4,407 training studies, so a\n# patient grouping is a no-op. Duplicate report text is the grouping that bites:\n# 49 texts are shared by 183 studies, and every study in a block carries identical\n# targets, so a block split across folds scores the model on a target whose source\n# it has already seen. \"auto\" additionally groups by acquisition site when the\n# series manifest exposes one, since scanner memorisation is reported to inflate\n# AUC by roughly 0.05 under an ungrouped split.\nGROUP_MODE = \"auto\"  # \"auto\" | \"report\" | \"site\"\nSITE_COLUMN_CANDIDATES = [\"SiteID\", \"Site\", \"site\", \"InstitutionName\",\n                          \"Institution\", \"StationName\", \"DeviceSerialNumber\"]\n\n# ---- Encode budget ---------------------------------------------------------\n# The visible test split is 3 studies; the rerun scores against a hidden split of\n# unknown size. The figure below is a planning assumption, not a measured count.\nASSUMED_RERUN_TEST_STUDIES = 1300\nKERNEL_BUDGET_S = 9 * 3600\nRESERVE_S = 2400        # head training, blending, submission\nSIZE_LADDER = [336, 280, 224]\nPROBE_STUDIES = 3\n\n# ---------------------------------------------------------------------------\n# Kaggle path resolution\n#\n# The attached competition publishes exactly seven entries at its root:\n#\n#     train_series/          test_series/           <- DICOM trees\n#     train.csv              test.csv               <- study manifests\n#     train_series.csv       test_series.csv        <- series manifests\n#     sample_submission.csv                         <- required output shape\n#\n# Kaggle mounts this one at /kaggle/input/competitions/<slug>, one level deeper\n# than the usual /kaggle/input/<slug>. Both shapes are handled below because the\n# Input panel does not show which is in effect; the resolved root is printed.\n#\n# The previous version of this cell returned the competitions path unconditionally\n# when no candidate matched, so an unattached dataset was reported as a valid root\n# here and only surfaced two cells later as a bare FileNotFoundError on train.csv.\n# Detection now walks what is really mounted and fails in this cell, naming\n# every directory it saw.\n# ---------------------------------------------------------------------------\nCOMP_SLUG = \"rsna-knee-abnormality-detection\"\n\n\ndef first_existing_dir(root, names):\n    \"\"\"Return the first of `names` that is a directory under root, else None.\"\"\"\n    for n in names:\n        p = root / n\n        try:\n            if p.is_dir():\n                return p\n        except OSError:\n            continue\n    return None\n\n\ndef detect_input_root():\n    \"\"\"Find the directory holding train.csv without trusting a hardcoded mount.\n\n    Returns:\n        Path to the competition input root.\n\n    Raises:\n        FileNotFoundError: naming both the paths checked and the directories\n            actually present under /kaggle/input.\n    \"\"\"\n    candidates = [\n        Path(f\"/kaggle/input/competitions/{COMP_SLUG}\"),\n        Path(f\"/kaggle/input/{COMP_SLUG}\"),\n        Path(\"data/raw\"),\n        Path(\"data\"),\n        Path(\".\"),\n    ]\n    for c in candidates:\n        try:\n            if (c / \"train.csv\").is_file():\n                return c\n        except OSError:\n            continue\n\n    # None of the known shapes matched. Walk two levels of /kaggle/input, which\n    # covers both /kaggle/input/<slug> and /kaggle/input/competitions/<slug>\n    # whatever the slug or nesting turns out to be.\n    seen = []\n    kin = Path(\"/kaggle/input\")\n    if kin.is_dir():\n        for d1 in sorted(p for p in kin.iterdir() if p.is_dir()):\n            seen.append(str(d1))\n            try:\n                if (d1 / \"train.csv\").is_file():\n                    return d1\n                for d2 in sorted(p for p in d1.iterdir() if p.is_dir()):\n                    seen.append(str(d2))\n                    if (d2 / \"train.csv\").is_file():\n                        return d2\n            except OSError:\n                continue\n\n    raise FileNotFoundError(\n        \"No input root containing train.csv was found.\\n\"\n        \"  Checked: \" + \", \".join(str(c) for c in candidates) + \"\\n\"\n        \"  Present under /kaggle/input: \" + (\", \".join(seen) if seen else \"nothing\") + \"\\n\"\n        \"  If that list is empty or does not contain the competition, the data is \"\n        \"not attached to this notebook. Open the right-hand panel -> Input -> \"\n        \"Add Input -> Competitions -> \" + COMP_SLUG + \", then rerun.\"\n    )\n\n\ndef resolve_kaggle_paths():\n    root = detect_input_root()\n\n    train_csv = root / \"train.csv\"\n    test_csv = root / \"test.csv\"\n    train_series_csv = root / \"train_series.csv\"\n    test_series_csv = root / \"test_series.csv\"\n    sample_sub_csv = root / \"sample_submission.csv\"\n\n    # The DICOM trees are train_series/ and test_series/ on this competition. The\n    # other names are tried only because neighbouring RSNA competitions use them;\n    # None is returned rather than a guessed path, and the ingestion cell already\n    # treats None as \"emit zero features\".\n    train_dcm_dir = first_existing_dir(root, [\"train_series\", \"train_images\", \"train\"])\n    test_dcm_dir = first_existing_dir(root, [\"test_series\", \"test_images\", \"test\"])\n\n    return (root, train_csv, test_csv, train_series_csv, test_series_csv,\n            sample_sub_csv, train_dcm_dir, test_dcm_dir)\n\n\n(ROOT, TRAIN_CSV, TEST_CSV, TRAIN_SERIES_CSV, TEST_SERIES_CSV,\n SAMPLE_SUB_CSV, TRAIN_DCM_DIR, TEST_DCM_DIR) = resolve_kaggle_paths()\nWORK = Path(\"/kaggle/working\")\n\nprint(f\"Kaggle Root Path : {ROOT}\", flush=True)\nprint(f\"Root contents    : {', '.join(sorted(p.name for p in ROOT.iterdir())[:24])}\", flush=True)\nprint(f\"Train DICOM Dir  : {TRAIN_DCM_DIR}\", flush=True)\nprint(f\"Test DICOM Dir   : {TEST_DCM_DIR}\", flush=True)\nfor _name, _p in [(\"train.csv\", TRAIN_CSV), (\"test.csv\", TEST_CSV),\n                  (\"train_series.csv\", TRAIN_SERIES_CSV),\n                  (\"test_series.csv\", TEST_SERIES_CSV),\n                  (\"sample_submission.csv\", SAMPLE_SUB_CSV)]:\n    print(f\"  {_name:<22} exists={_p.is_file()}\", flush=True)\n\nif TRAIN_DCM_DIR is None or TEST_DCM_DIR is None:\n    print(\"WARNING: a DICOM directory was not found under the input root. Those \"\n          \"studies will be encoded as zero vectors and the run will score at \"\n          \"chance. Check the mount before trusting the output.\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"b3d7380f","metadata":{},"outputs":[],"source":"# Cell 3: Manifest & Data Loading\ntrain = pd.read_csv(TRAIN_CSV)\ntest = pd.read_csv(TEST_CSV)\n\ntrain_ids = train[\"StudyInstanceUID\"].astype(str).tolist()\ntest_ids = test[\"StudyInstanceUID\"].astype(str).tolist()\n\ndef load_manifest(path):\n    if path.is_file():\n        df = pd.read_csv(path)\n        df[\"StudyInstanceUID\"] = df[\"StudyInstanceUID\"].astype(str)\n        df[\"SeriesInstanceUID\"] = df[\"SeriesInstanceUID\"].astype(str)\n        df[\"plane\"] = df[\"Anatomical_Plane\"].astype(str).str.strip().str.title()\n        df[\"fluid\"] = pd.to_numeric(df[\"Fluid_Sensitive\"], errors=\"coerce\").fillna(0).astype(int)\n        return df\n    return pd.DataFrame(columns=[\"StudyInstanceUID\", \"SeriesInstanceUID\", \"Anatomical_Plane\", \"Fluid_Sensitive\", \"plane\", \"fluid\"])\n\ntrain_series = load_manifest(TRAIN_SERIES_CSV)\ntest_series = load_manifest(TEST_SERIES_CSV)\n\nprint(f\"Train Manifest: {len(train_ids)} studies, {len(train_series)} series\", flush=True)\nprint(f\"Test Manifest : {len(test_ids)} studies, {len(test_series)} series\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"97a1e6c5","metadata":{},"outputs":[],"source":"# Cell 4: Multilingual Rule-Based Extractor & Soft Target Generation\n#\n# The extractor below is embedded verbatim from the companion notebook rather\n# than rewritten here, so the two cannot drift apart. The version this replaced\n# was a condensed paraphrase and measured materially worse on the same 58 gold\n# studies:\n#\n#   this notebook, before   agreement 0.7500  precision 0.6291  recall 0.7129\n#   companion extractor     agreement 0.7960  precision 0.6802  recall 0.7648\n#\n# The paraphrase lost four things that matter:\n#   * negation was tested across the whole clause, so one \"no\" anywhere killed\n#     every finding in it. Here it is windowed, 70 chars back and 45 forward.\n#   * a negated first occurrence ended the search. _fire_direct retries later\n#     occurrences of the same term.\n#   * paired labels required a named structure. GENERIC + SIDE_Q now catch the\n#     Greek and Bulgarian phrasings that name the compartment instead, and\n#     PAIR_WINDOW requires the finding word to sit near the anatomy word.\n#   * \"medial patellar facet\" fired Medial OA. PF_MASK removes those spans\n#     before the medial/lateral qualifiers are searched.\n#\n# Verify the port by reading the agreement table this cell prints: it should\n# land at 0.7960 / 0.6802 / 0.7648, not the old 0.7500 / 0.6291 / 0.7129.\n\n# The rule extractor below is not written here. It is embedded verbatim from the\n# builder module of the companion metadata-baseline notebook, so the two notebooks\n# consume one extractor and cannot drift apart. Its audit against the gold studies\n# is printed two cells down.\n# --------------------------------------------------------------------------- #\n# Multilingual rule-based label extractor.\n# Weak on purpose. This is the artifact an LLM extraction has to beat.\n# --------------------------------------------------------------------------- #\n\n_CHAR_MAP = str.maketrans({\"ı\": \"i\", \"İ\": \"i\", \"ß\": \"ss\",\n                           \"ø\": \"o\", \"Ø\": \"o\",\n                           \"đ\": \"d\", \"Đ\": \"d\"})\n\nTEAR = [\n    \"tear\", \"torn\", \"tearing\", \"rupture\", \"ruptured\", \"disruption\", \"discontinuity\",\n    \"rotura\", \"ruptura\", \"desgarro\", \"roto\", \"rota\",\n    \"dechirure\", \"dechire\", \"lesion meniscale\",\n    \"scheur\", \"ruptuur\", \"gescheurd\",\n    \"riss\", \"rissbildung\", \"ruptur\", \"zerreissung\", \"einriss\",\n    \"yirtik\", \"yirtig\", \"kopma\", \"butunluk kaybi\",\n    \"ρηξη\", \"ρηξις\", \"ρηγμα\",\n    \"руптура\", \"разкъсв\",\n    \"разрив\", \"puknuce\", \"prekid\",\n]\nSPRAIN = [\n    \"sprain\", \"esguince\", \"entorse\", \"verstauchung\", \"verstuiking\",\n    \"distorsiyon\", \"burkulma\", \"διαστρεμμα\",\n    \"навяхван\",\n    \"injury\", \"lesion\", \"letsel\", \"verletzung\", \"zedelenme\",\n]\nOA_FIND = [\n    \"osteoarthritis\", \"osteoarthrosis\", \"arthrosis\", \"arthritic\",\n    \"artrosis\", \"artrose\", \"arthrose\", \"gonartrose\", \"gonartroz\", \"gonarthrose\",\n    \"gonartro\", \"artroz\",\n    \"chondrosis\", \"chondral\", \"chondropathy\", \"chondropathie\", \"chondropatie\",\n    \"chondromalacia\", \"kondromalazi\", \"condropatia\", \"condral\", \"chondropatia\",\n    \"kraakbeenlijden\", \"kraakbeenverlies\", \"kraakbeenschade\",\n    \"knorpelschaden\", \"knorpeldefekt\", \"knorpelverlust\", \"knorpelbelag\",\n    \"cartilage loss\", \"cartilage thinning\", \"cartilage fissuring\",\n    \"cartilage defect\", \"cartilage fissure\", \"chondral defect\", \"chondral loss\",\n    \"osteophyt\", \"osteofit\", \"osteofyt\", \"osteophyte\", \"spurring\", \"osteofyte\",\n    \"joint space narrowing\", \"gelenkspaltverschmalerung\",\n    \"kikirdak kaybi\", \"kikirdak incelme\", \"kondral\",\n    \"ulcera condral\", \"hrskavice\", \"denudacija\",\n    \"χονδροπαθ\", \"οστεοαρθρ\",\n    \"οστεοφυτ\", \"αρθριτ\",\n    \"артроз\", \"хондропат\",\n    \"остеофит\", \"хрущялн\",\n]\nMEDIAL_Q = [\"medial\", \"interno\", \"interna\", \"interne\", \"mediaal\", \"mediale\",\n            \"innen\", \"medyal\", \"medijaln\", \"εσω\",\n            \"медиал\", \"вътреш\"]\nLATERAL_Q = [\"lateral\", \"externo\", \"externa\", \"externe\", \"buiten\", \"aussen\",\n             \"lateraln\", \"εξω\", \"латерал\",\n             \"външ\"]\nPF_Q = [\"patellofemoral\", \"patelofemoral\", \"femoropatellar\", \"femoropatelar\",\n        \"femoropatellair\", \"retropatellar\", \"retrorotulian\", \"patellar facet\",\n        \"trochlea\", \"troclea\", \"trochlee\", \"rotula\", \"rotulian\", \"patella\",\n        \"patellaire\", \"patellar\", \"diz kapagi\", \"patelofemoraln\",\n        \"επιγονατιδ\", \"τροχιλ\",\n        \"пател\", \"ретропател\"]\nTRICOMP = [\"tricompartmental\", \"tricompartimental\", \"three compartments\",\n           \"three compartmens\", \"all compartments\", \"gonartrose\", \"gonartroz\",\n           \"gonarthrose\", \"gonartrosis\", \"gonartro\", \"pangonartro\"]\n\nSPECIFIC = {\n    \"ACL\": [\n        \"acl\", \"anterior cruciate\", \"lca\", \"ligamento cruzado anterior\",\n        \"ligament croise anterieur\", \"croise anterieur\",\n        \"voorste kruisband\", \"vkb\", \"vorderes kreuzband\", \"vorderen kreuzband\",\n        \"vordere kreuzband\", \"on capraz bag\", \"anterior capraz bag\",\n        \"προσθιο χιαστ\",\n        \"προσθιου χιαστ\",\n        \"προσθιος χιαστ\",\n        \"предна кръстна\",\n        \"предната кръстна\",\n        \"prednji ukrizeni\",\n    ],\n    \"MCL\": [\n        \"mcl\", \"medial collateral\", \"lcm\", \"ligamento colateral medial\",\n        \"ligamento colateral interno\", \"ligament collateral medial\",\n        \"collateral medial\", \"mediale collaterale\", \"mediaal collateraal\",\n        \"innenband\", \"mediales kollateralband\", \"medialen kollateralband\",\n        \"mediale kollateralband\", \"medyal kollateral\", \"ic yan bag\",\n        \"εσω πλαγιο\",\n        \"медиален колатерален\",\n        \"вътрешна колатерална\",\n    ],\n    \"Medial Meniscus\": [\n        \"medial meniscus\", \"meniscus medialis\", \"menisco medial\", \"menisco interno\",\n        \"mediale meniscus\", \"binnenmeniscus\", \"innenmeniskus\", \"meniskus medialis\",\n        \"menisque interne\", \"menisque medial\", \"medial menisk\", \"medyal menisk\",\n        \"εσω μηνισκ\",\n        \"медиалния менискус\",\n        \"медиален мениск\",\n        \"вътрешния мениск\",\n        \"medial and lateral menisc\", \"medial ve lateral menisk\",\n        \"menisco medial y lateral\", \"menisco interno y externo\",\n        \"menisco interno y lateral\", \"mediale en laterale meniscus\",\n        \"innen und aussenmeniskus\", \"medijalnog meniskusa\",\n    ],\n    \"Lateral Meniscus\": [\n        \"lateral meniscus\", \"meniscus lateralis\", \"menisco lateral\", \"menisco externo\",\n        \"laterale meniscus\", \"buitenmeniscus\", \"aussenmeniskus\", \"meniskus lateralis\",\n        \"menisque externe\", \"menisque lateral\", \"lateral menisk\",\n        \"εξω μηνισκ\",\n        \"латералния менискус\",\n        \"латерален мениск\",\n        \"външния мениск\",\n        \"medial and lateral menisc\", \"medial ve lateral menisk\",\n        \"menisco medial y lateral\", \"menisco interno y externo\",\n        \"menisco interno y lateral\", \"mediale en laterale meniscus\",\n        \"innen und aussenmeniskus\", \"lateralnog meniskusa\",\n    ],\n}\n\nGENERIC = {\n    \"ACL\": [\"cruciate\", \"cruzado\", \"croise\", \"kruisband\", \"kreuzband\",\n            \"χιαστ\", \"кръстн\",\n            \"capraz bag\", \"ukrizen\"],\n    \"MCL\": [\"collateral\", \"colateral\", \"kollateral\", \"collaterale\",\n            \"πλαγι\", \"колатерал\",\n            \"yan bag\"],\n    \"Medial Meniscus\": [\"menisc\", \"menisk\", \"μηνισκ\",\n                        \"мениск\"],\n    \"Lateral Meniscus\": [\"menisc\", \"menisk\", \"μηνισκ\",\n                         \"мениск\"],\n}\nSIDE_Q = {\n    \"ACL\": [\"anterior\", \"anterieur\", \"voorste\", \"vorder\",\n            \"προσθι\", \"предн\", \"prednj\"],\n    \"MCL\": MEDIAL_Q,\n    \"Medial Meniscus\": MEDIAL_Q,\n    \"Lateral Meniscus\": LATERAL_Q,\n}\n\nDIRECT = {\n    \"Effusion\": [\n        \"effusion\", \"joint fluid\", \"hemarthrosis\", \"haemarthrosis\", \"hydrops\",\n        \"derrame\", \"epanchement\", \"gewrichtsvocht\", \"vocht in het gewricht\",\n        \"erguss\", \"gelenkerguss\", \"gelenkserguss\",\n        \"eklem ici sivi\", \"sivi artisi\", \"efuzyon\", \"eklem mesafesinde sivi\",\n        \"eklem ici serbest sivi\", \"eklem sivisi\",\n        \"ενδαρθρικ\",\n        \"αρθρικο υγρο\",\n        \"συλλογη υγρου\",\n        \"ставен излив\",\n        \"излив\", \"хидропс\",\n        \"zglobni izljev\", \"izljev\",\n    ],\n    \"Synovitis\": [\n        \"synovitis\", \"synovial thickening\", \"synovial hypertrophy\",\n        \"thickened synovial\", \"hypertrophy of the synovium\",\n        \"synovial proliferation\", \"proliferation of the synovium\",\n        \"sinovitis\", \"synovite\", \"synovitiden\", \"synovialitis\",\n        \"verdikking van het synovium\", \"verdikkingen van het synovium\",\n        \"synoviale verdikking\", \"synovialisverdickung\", \"synovialverdickung\",\n        \"sinovit\", \"hoffitis\", \"sinovyal kalinlasma\", \"sinovyal proliferasyon\",\n        \"συνοβιτ\", \"υμενιτ\",\n        \"синовит\", \"sinovij\", \"sinovije\",\n    ],\n    \"Baker's\": [\n        \"baker\", \"popliteal cyst\", \"poplitealcyst\", \"popliteal cysts\",\n        \"quiste popliteo\", \"quistes popliteos\", \"quiste de baker\",\n        \"kyste de baker\", \"kyste poplite\", \"popliteale cyste\", \"popliteale cyst\",\n        \"bakercyste\", \"bakerzyste\", \"poplitealzyste\", \"popliteazyste\",\n        \"baker kisti\", \"popliteal kist\",\n        \"κυστη baker\", \"κυστη του baker\",\n        \"киста на бейкър\",\n        \"бейкърова киста\",\n        \"poplitealna cista\",\n    ],\n    \"Contusion\": [\n        \"contusion\", \"contusiones\", \"contusie\", \"kontusyon\", \"kontuzyon\",\n        \"bone bruise\", \"bone bruising\", \"knochenprellung\", \"prellung\",\n        \"botcontusie\", \"μωλωπ\", \"θλαση\",\n        \"контузи\", \"kontuzij\", \"nagnjecen\",\n    ],\n    \"Fracture\": [\n        \"fracture\", \"fractur\", \"fractura\", \"fraktur\", \"fractuur\", \"breuk\",\n        \"kirik\", \"kirig\", \"καταγμα\",\n        \"καταγματ\",\n        \"фрактура\", \"счупван\",\n        \"avulsion\", \"avulsie\", \"avulsiyon\", \"prijelom\",\n    ],\n}\n\n# Negation cues, word-boundary anchored. A bare \"no \" substring matches inside the\n# Spanish word \"cuerno\", which silently killed most Spanish positives before this\n# was anchored.\nNEG_SRC = [\n    r\"\\bno\\b\", r\"\\bnot\\b\", r\"\\bnon\\b\", r\"\\bwithout\\b\", r\"\\babsen\\w*\",\n    r\"\\bnegative for\\b\", r\"\\bunremarkable\\b\", r\"\\bintact\\w*\", r\"\\bnormal\\w*\",\n    r\"\\bpreserved\\b\", r\"\\bfree of\\b\", r\"\\bexcluded\\b\",\n    r\"\\bno hay\\b\", r\"\\bsin\\b\", r\"\\bausen\\w*\", r\"\\bno se\\b\", r\"\\bconservad\\w*\",\n    r\"\\bintegr\\w*\", r\"\\bindemne\\b\",\n    r\"\\bpas de\\b\", r\"\\baucun\\w*\", r\"\\bsans\\b\",\n    r\"\\bgeen\\b\", r\"\\bzonder\\b\", r\"\\bonopvallend\\w*\", r\"\\bnormaal\\b\",\n    r\"\\bvrij van\\b\", r\"\\bintacte\\b\",\n    r\"\\bkein\\w*\", r\"\\bohne\\b\", r\"\\bunauffallig\\w*\", r\"\\bregelrecht\\w*\",\n    r\"\\bnicht\\b\", r\"\\bintakt\\w*\",\n    r\"\\byok\\w*\", r\"\\bizlenmemis\\w*\", r\"\\bizlenmedi\\w*\", r\"\\bsaptanmamis\\w*\",\n    r\"\\bgozlenmemis\\w*\", r\"\\bgorulmemis\\w*\", r\"\\bkorunmus\\w*\",\n    r\"\\bmevcut degil\\b\", r\"\\bnormaldir\\b\", r\"\\bdogaldir\\b\",\n    r\"\\bδεν\\b\", r\"\\bχωρις\\b\",\n    r\"\\bφυσιολογικ\\w*\",\n    r\"\\bακεραι\\w*\",\n    r\"\\bбез\\b\", r\"\\bняма\\b\",\n    r\"\\bне се\\b\", r\"\\bнормал\\w*\",\n    r\"\\bзапазен\\w*\",\n    r\"\\bсъхранен\\w*\",\n    r\"\\bsenza\\b\", r\"\\bsem\\b\", r\"\\bnao\\b\", r\"\\bbez\\b\", r\"\\buredn\\w*\",\n]\nNEG_RE = re.compile(\"|\".join(NEG_SRC))\n\nNEG_BACK = 70\nNEG_FWD = 45\nQUAL_WINDOW = 90\nPAIR_WINDOW = 160\nCLAUSE_SPLIT = re.compile(r\"[.;:\\n\\r•·]+|\\s-\\s|\\s>\\s|\\s\\*\\s\")\nPF_MASK = re.compile(r\"(medial|lateral)\\s+(patellar|patella|facet|trochlea|trochlear|retinac)\\w*\")\n\n\ndef fold_text(text):\n    \"\"\"Lowercase, normalise script-specific letters, and strip diacritics.\n\n    Turkish dotless i, the German sharp s and the Croatian barred d have no\n    combining-mark decomposition, so they are mapped explicitly before NFKD.\n\n    Args:\n        text: Any report text.\n\n    Returns:\n        A lowercase, accent-free string safe for substring matching.\n    \"\"\"\n    t = str(text).lower().translate(_CHAR_MAP)\n    t = unicodedata.normalize(\"NFKD\", t)\n    return \"\".join(c for c in t if not unicodedata.combining(c))\n\n\ndef _find_any(clause, terms):\n    \"\"\"Return the earliest index at which any term occurs, or -1.\n\n    Args:\n        clause: Folded clause text.\n        terms: Iterable of folded surface forms.\n\n    Returns:\n        Character index of the earliest hit, or -1 when none matches.\n    \"\"\"\n    best = -1\n    for t in terms:\n        i = clause.find(t)\n        if i >= 0 and (best < 0 or i < best):\n            best = i\n    return best\n\n\ndef _negated(clause, pos, span):\n    \"\"\"Test whether a negation cue sits near a matched finding term.\n\n    Args:\n        clause: Folded clause text.\n        pos: Start index of the finding term.\n        span: Length of the finding term.\n\n    Returns:\n        True when a cue appears in the preceding or following window.\n    \"\"\"\n    back = clause[max(0, pos - NEG_BACK):pos]\n    fwd = clause[pos + span:pos + span + NEG_FWD]\n    return bool(NEG_RE.search(back) or NEG_RE.search(fwd))\n\n\ndef _fire_direct(clause, terms):\n    \"\"\"Test a direct finding term, retrying later occurrences past a negation.\n\n    Args:\n        clause: Folded clause text.\n        terms: Surface forms that are themselves the finding.\n\n    Returns:\n        True when at least one occurrence is unnegated.\n    \"\"\"\n    for t in terms:\n        start = 0\n        while True:\n            i = clause.find(t, start)\n            if i < 0:\n                break\n            if not _negated(clause, i, len(t)):\n                return True\n            start = i + 1\n    return False\n\n\ndef _anatomy_index(clause, label):\n    \"\"\"Locate the anatomy mention for a paired label.\n\n    Falls back to a generic organ term paired with a side qualifier, which is what\n    catches Greek and Bulgarian phrasing that names the compartment rather than the\n    structure.\n\n    Args:\n        clause: Folded clause text.\n        label: One of the four paired label names.\n\n    Returns:\n        Character index of the anatomy mention, or -1.\n    \"\"\"\n    i = _find_any(clause, SPECIFIC[label])\n    if i >= 0:\n        return i\n    g = _find_any(clause, GENERIC[label])\n    if g >= 0 and _find_any(clause, SIDE_Q[label]) >= 0:\n        return g\n    return -1\n\n\ndef extract_labels(report):\n    \"\"\"Extract twelve binary findings from one free-text radiology report.\n\n    Args:\n        report: Report text in any of the languages present in the corpus.\n\n    Returns:\n        Dict mapping each of the twelve label names to 0 or 1.\n    \"\"\"\n    out = {lab: 0 for lab in LABELS}\n    if not isinstance(report, str) or not report.strip():\n        return out\n    text = fold_text(report)\n    for raw in CLAUSE_SPLIT.split(text):\n        clause = raw.strip()\n        if len(clause) < 3:\n            continue\n        for lab in (\"Effusion\", \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\"):\n            if not out[lab] and _fire_direct(clause, DIRECT[lab]):\n                out[lab] = 1\n        for lab in SPECIFIC:\n            if out[lab]:\n                continue\n            ai = _anatomy_index(clause, lab)\n            if ai < 0:\n                continue\n            finds = TEAR + SPRAIN if lab in (\"ACL\", \"MCL\") else TEAR\n            for t in finds:\n                fi = clause.find(t)\n                if fi < 0 or abs(fi - ai) > PAIR_WINDOW:\n                    continue\n                if not _negated(clause, fi, len(t)):\n                    out[lab] = 1\n                    break\n        oi = _find_any(clause, OA_FIND)\n        if oi >= 0 and not _negated(clause, oi, 8):\n            if _find_any(clause, TRICOMP) >= 0:\n                out[\"Medial OA\"] = 1\n                out[\"Lateral OA\"] = 1\n                out[\"PF OA\"] = 1\n            win = clause[max(0, oi - QUAL_WINDOW):oi + QUAL_WINDOW]\n            if _find_any(win, PF_Q) >= 0:\n                out[\"PF OA\"] = 1\n            masked = PF_MASK.sub(\" \", win)\n            if _find_any(masked, MEDIAL_Q) >= 0:\n                out[\"Medial OA\"] = 1\n            if _find_any(masked, LATERAL_Q) >= 0:\n                out[\"Lateral OA\"] = 1\n    return out\n\n\nreport_col = \"Report\" if \"Report\" in train.columns else None\nif report_col:\n    rule_df = pd.DataFrame([extract_labels(r) for r in train[report_col]], index=pd.Index(train_ids, name=\"StudyInstanceUID\"))\nelse:\n    rule_df = pd.DataFrame(0, index=pd.Index(train_ids, name=\"StudyInstanceUID\"), columns=LABELS)\n\nY_soft = rule_df[LABELS].astype(float)\nLABEL_PROVENANCE = \"rule_inline\"\n\n\n# ---------------------------------------------------------------------------\n# Optional stronger label table\n#\n# The rule extractor scores 0.7960 agreement against the 58 gold studies. An LLM\n# extraction of the same reports scores 0.8233 (paired delta +0.0273, McNemar\n# p=0.0204), and its recall is the bigger gap: 0.8602 against 0.7648. If such a\n# table is mounted as a dataset input it overwrites the rule labels here.\n#\n# Nothing is fetched: this reads a file you attached. See the compliance note in\n# Cell 2 on how that file must be produced.\n# ---------------------------------------------------------------------------\ndef label_search_dirs():\n    \"\"\"Directories that may hold a pre-extracted label table.\n\n    Kaggle nests inputs under a category directory rather than mounting them\n    flat, so a single listing of /kaggle/input finds nothing. This walks a\n    bounded depth and prunes the competition tree, which holds hundreds of\n    thousands of DICOM files and would cost more to walk than to encode.\n\n    Returns:\n        Existing directories, shallowest first.\n    \"\"\"\n    out = []\n    kin = Path(\"/kaggle/input\")\n    if kin.is_dir():\n        out.append(kin)\n        stack = [(kin, 0)]\n        while stack:\n            d, depth = stack.pop()\n            if depth >= 5:\n                continue\n            try:\n                children = sorted(p for p in d.iterdir() if p.is_dir())\n            except OSError:\n                continue\n            for c in children:\n                if c.name in (\"competitions\", \"train_series\", \"test_series\",\n                              \"train_images\", \"test_images\", \".git\", \"__pycache__\"):\n                    continue\n                out.append(c)\n                stack.append((c, depth + 1))\n    out.append(Path(\"data/interim\"))\n    seen, uniq = set(), []\n    for d in out:\n        if str(d) not in seen and d.is_dir():\n            seen.add(str(d))\n            uniq.append(d)\n    return uniq\n\n\ndef find_label_file(patterns):\n    \"\"\"Return the newest file matching any pattern across the search dirs.\"\"\"\n    hits = []\n    for d in label_search_dirs():\n        for pat in patterns:\n            try:\n                hits.extend(sorted(d.glob(pat)))\n            except OSError:\n                continue\n    if not hits:\n        return None\n    return sorted(hits, key=lambda p: (p.name, str(p)))[-1]\n\n\ndef read_label_table(path):\n    \"\"\"Read a label table, keeping the id column and the twelve findings.\n\n    Returns:\n        DataFrame indexed by StudyInstanceUID with float labels in [0, 1], or\n        None when the file is unreadable or lacks the required columns.\n    \"\"\"\n    try:\n        df = pd.read_parquet(path) if path.suffix == \".parquet\" else pd.read_csv(path)\n    except Exception as exc:\n        print(f\"Could not read {path}: {exc}\", flush=True)\n        return None\n    if \"StudyInstanceUID\" not in df.columns or not set(LABELS).issubset(df.columns):\n        missing = [c for c in LABELS if c not in df.columns]\n        print(f\"Label file {path.name} is missing columns {missing[:4]}; ignored.\", flush=True)\n        return None\n    out = df[[\"StudyInstanceUID\"] + LABELS].copy()\n    out[\"StudyInstanceUID\"] = out[\"StudyInstanceUID\"].astype(str)\n    out = out.drop_duplicates(\"StudyInstanceUID\").set_index(\"StudyInstanceUID\")\n    return out[LABELS].astype(float).clip(0.0, 1.0)\n\n\nif LABEL_SOURCE in (\"auto\", \"file\"):\n    _p = find_label_file(LABEL_FILE_PATTERNS)\n    if _p is not None:\n        _tbl = read_label_table(_p)\n        if _tbl is not None:\n            _common = Y_soft.index.intersection(_tbl.index)\n            Y_soft.loc[_common, LABELS] = _tbl.loc[_common, LABELS].values\n            LABEL_PROVENANCE = f\"file:{_p.name}\"\n            print(f\"Label table applied: {_p} | {len(_common)} of {len(Y_soft)} studies overwritten\", flush=True)\n    elif LABEL_SOURCE == \"file\":\n        print(\"LABEL_SOURCE='file' but no label table is mounted; keeping rule labels.\", flush=True)\n\n# Gold Ground Truth: 58 of 4,407 studies carry radiologist labels in train.csv.\ngold_mask = train[LABELS].notna().all(axis=1) if set(LABELS).issubset(train.columns) else pd.Series(False, index=train.index)\ngold_ids = train.loc[gold_mask, \"StudyInstanceUID\"].astype(str).tolist() if gold_mask.any() else []\n\n_pos_of = {sid: i for i, sid in enumerate(train_ids)}\nGOLD_POS = [_pos_of[g] for g in gold_ids if g in _pos_of]\nY_GOLD = None\n\nif gold_ids:\n    gold_vals = train.loc[gold_mask, [\"StudyInstanceUID\"] + LABELS].set_index(\"StudyInstanceUID\")\n    Y_GOLD = (gold_vals.loc[[train_ids[i] for i in GOLD_POS], LABELS].astype(float).values >= 0.5).astype(int)\n\n    # Score the active weak-label source against gold BEFORE gold overwrites it.\n    # This is the only honest read on label quality available in-notebook.\n    _pred = (Y_soft.loc[[train_ids[i] for i in GOLD_POS], LABELS].values >= 0.5).astype(int)\n    rows = []\n    for j, lab in enumerate(LABELS):\n        g, p = Y_GOLD[:, j], _pred[:, j]\n        tp = int((g & p).sum())\n        rows.append({\n            \"label\": lab,\n            \"gold_pos\": int(g.sum()),\n            \"extractor_pos\": int(p.sum()),\n            \"agreement\": round(float((g == p).mean()), 4),\n            \"precision\": round(tp / max(1, int(p.sum())), 4),\n            \"recall\": round(tp / max(1, int(g.sum())), 4),\n        })\n    agree_df = pd.DataFrame(rows)\n    print(f\"\\nWeak labels ({LABEL_PROVENANCE}) vs gold, n = {len(GOLD_POS)} studies\", flush=True)\n    print(agree_df.to_string(index=False), flush=True)\n    print(\"mean agreement {:.4f} | mean precision {:.4f} | mean recall {:.4f}\".format(\n        agree_df[\"agreement\"].mean(), agree_df[\"precision\"].mean(), agree_df[\"recall\"].mean()), flush=True)\n    print(\"reference: rule 0.7960 / 0.6802 / 0.7648, LLM 0.8233 / 0.6910 / 0.8602 \"\n          \"(companion notebook, same 58 studies)\", flush=True)\n\n    Y_soft.loc[gold_ids, LABELS] = gold_vals[LABELS].astype(float).values\n\n\ndef gold_macro_auc(pred_full):\n    \"\"\"Macro ROC-AUC of out-of-fold predictions against the gold labels.\n\n    This is the number that tracks the leaderboard metric. The in-notebook CV\n    figure scores predictions against the same extractor that produced the\n    targets, so it measures self-consistency, not accuracy.\n\n    Args:\n        pred_full: (n_train, N_LABELS) predictions aligned to train_ids.\n\n    Returns:\n        (macro_auc, n_labels_scored). NaN when no gold labels are available.\n    \"\"\"\n    if Y_GOLD is None or not GOLD_POS:\n        return float(\"nan\"), 0\n    yp = np.asarray(pred_full)[GOLD_POS]\n    aucs = []\n    for j in range(N_LABELS):\n        if len(np.unique(Y_GOLD[:, j])) > 1:\n            aucs.append(roc_auc_score(Y_GOLD[:, j], yp[:, j]))\n    return (float(np.mean(aucs)) if aucs else float(\"nan\")), len(aucs)\n\n\nprint(f\"\\nSoft labels: {len(Y_soft)} studies | source: {LABEL_PROVENANCE} | \"\n      f\"gold overrides: {len(gold_ids)} studies (weighted x{GOLD_WEIGHT} in the loss)\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"41c41403","metadata":{},"outputs":[],"source":"# Cell 5: Fast Parallel DICOM Preprocessing, MONOCHROME1 Inversion & Laterality\nimport cv2\n\ndef resize_2d(arr, size):\n    h, w = arr.shape\n    if (h, w) == (size, size):\n        return arr.astype(np.float32)\n    interp = cv2.INTER_AREA if max(h, w) > size else cv2.INTER_LINEAR\n    return cv2.resize(arr.astype(np.float32), (size, size), interpolation=interp).astype(np.float32)\n\ndef read_slice_pixels(path):\n    try:\n        import pydicom\n        ds = pydicom.dcmread(path, force=True)\n        arr = ds.pixel_array.astype(np.float32)\n        \n        # 1. Apply Rescale Slope & Intercept if present\n        slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n        intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n        if slope != 1.0 or intercept != 0.0:\n            arr = arr * slope + intercept\n            \n        # 2. Handle MONOCHROME1 contrast inversion (white=0, black=max -> invert to standard MONOCHROME2)\n        photo = str(getattr(ds, \"PhotometricInterpretation\", \"\")).strip().upper()\n        if photo == \"MONOCHROME1\":\n            arr = arr.max() - arr\n            \n        return arr\n    except Exception:\n        return None\n\n# Header-Only Physical Position & Laterality Readers (stop_before_pixels=True)\ndef get_dicom_zpos_fast(path):\n    try:\n        import pydicom\n        ds = pydicom.dcmread(path, stop_before_pixels=True, force=True)\n        ipp = [float(x) for x in ds.ImagePositionPatient]\n        if hasattr(ds, \"ImageOrientationPatient\"):\n            iop = [float(x) for x in ds.ImageOrientationPatient]\n            nx = iop[1] * iop[5] - iop[2] * iop[4]\n            ny = iop[2] * iop[3] - iop[0] * iop[5]\n            nz = iop[0] * iop[4] - iop[1] * iop[3]\n            return ipp[0] * nx + ipp[1] * ny + ipp[2] * nz\n        return ipp[2]\n    except Exception:\n        return 0.0\n\ndef get_series_laterality_fast(path):\n    try:\n        import pydicom\n        ds = pydicom.dcmread(path, stop_before_pixels=True, force=True)\n        lat = str(getattr(ds, \"ImageLaterality\", \"\")).strip().upper()\n        if not lat and hasattr(ds, \"Laterality\"):\n            lat = str(getattr(ds, \"Laterality\", \"\")).strip().upper()\n        return lat\n    except Exception:\n        return \"\"\n\ndef normalise_laterality(norm_planes, plane, lat):\n    if lat != \"R\":\n        return norm_planes\n    if plane in (\"Coronal\", \"Axial\"):\n        # Mirror horizontally onto left-knee orientation\n        return [p[:, ::-1] for p in norm_planes]\n    # Reverse Sagittal stack\n    return norm_planes[::-1]\n\ndef order_series_fast(folder_path):\n    dcm_files = sorted(Path(folder_path).glob(\"*.dcm\"))\n    if not dcm_files:\n        return [], len(dcm_files), \"\"\n    positions = [(f, get_dicom_zpos_fast(str(f))) for f in dcm_files]\n    positions.sort(key=lambda x: x[1])\n    sorted_paths = [p[0] for p in positions]\n    lat = get_series_laterality_fast(sorted_paths[0]) if sorted_paths else \"\"\n    return sorted_paths, len(dcm_files), lat\n\ndef build_slot_stack_fast(folder_path, plane, k, size):\n    paths, n_files, lat = order_series_fast(folder_path)\n    if not paths:\n        return np.zeros((0, 3, size, size), np.uint8)\n    \n    lo_frac, hi_frac = PLANE_WINDOW.get(plane, (0.12, 0.88))\n    n = len(paths)\n    lo_idx = int(n * lo_frac)\n    hi_idx = max(lo_idx + 1, int(n * hi_frac))\n    windowed = paths[lo_idx:hi_idx]\n    if not windowed:\n        windowed = paths\n    \n    if len(windowed) > k:\n        step = len(windowed) / k\n        indices = [int(i * step) for i in range(k)]\n        selected_paths = [windowed[i] for i in indices]\n    else:\n        selected_paths = windowed\n        \n    arrays = []\n    for p in selected_paths:\n        arr = read_slice_pixels(p)\n        if arr is not None:\n            arrays.append(arr)\n            \n    if not arrays:\n        return np.zeros((0, 3, size, size), np.uint8)\n        \n    all_pix = np.concatenate([a.ravel()[::16] for a in arrays])\n    low = np.percentile(all_pix, 1.0)\n    high = np.percentile(all_pix, 99.0)\n    span = max(high - low, 1e-6)\n    \n    norm_planes = []\n    for a in arrays:\n        resized = resize_2d(a, size)\n        normed = np.clip((resized - low) / span, 0.0, 1.0)\n        norm_planes.append((normed * 255.0).astype(np.uint8))\n        \n    # Apply Anatomical Laterality Normalization\n    norm_planes = normalise_laterality(norm_planes, plane, lat)\n    \n    m = len(norm_planes)\n    imgs = np.zeros((m, 3, size, size), np.uint8)\n    for j in range(m):\n        imgs[j, 0] = norm_planes[max(0, j - 1)]\n        imgs[j, 1] = norm_planes[j]\n        imgs[j, 2] = norm_planes[min(m - 1, j + 1)]\n        \n    return imgs\n\nprint(\"✓ Ultra-fast DICOM header reader & MONOCHROME1 / laterality normalizer ready.\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"450eac63","metadata":{},"outputs":[],"source":"# Cell 6: Backbone Initialization from a Kaggle Model Mount\nfrom transformers import AutoModel\n\nIMAGENET_MEAN = np.array([0.485, 0.456, 0.406], dtype=np.float32).reshape(1, 3, 1, 1)\nIMAGENET_STD = np.array([0.229, 0.224, 0.225], dtype=np.float32).reshape(1, 3, 1, 1)\n\n\nskip_dirs = {\"train_series\", \"test_series\", \"train_images\", \"test_images\",\n             \"train\", \"test\", \"competitions\", \".git\", \"__pycache__\"}\n\n\ndef find_local_backbones():\n    \"\"\"Every checkpoint directory mounted under /kaggle/input.\n\n    A checkpoint is a directory holding config.json next to a weights file. The\n    competition DICOM trees are pruned from the walk: they hold hundreds of\n    thousands of files and can never contain one.\n\n    Returns:\n        List of candidate directories.\n    \"\"\"\n    found, seen = [], set()\n    skip = skip_dirs\n    for base in (Path(\"/kaggle/input/models\"), Path(\"/kaggle/input\")):\n        if not base.is_dir():\n            continue\n        for root, dirs, files in os.walk(base):\n            dirs[:] = [d for d in dirs if d not in skip]\n            if len(Path(root).parts) - len(base.parts) > 8:\n                dirs[:] = []\n                continue\n            if \"config.json\" in files and any(\n                    f.endswith((\".safetensors\", \".bin\")) for f in files):\n                if root not in seen:\n                    seen.add(root)\n                    found.append(Path(root))\n    return found\n\n\ndef checkpoint_hidden_size(path):\n    \"\"\"Transformer width declared by a checkpoint, or 0 when unreadable.\n\n    Args:\n        path: Directory holding config.json.\n\n    Returns:\n        hidden_size as an int, 0 if the file is missing or malformed.\n    \"\"\"\n    try:\n        with open(path / \"config.json\") as fh:\n            cfg = json.load(fh)\n        return int(cfg.get(\"hidden_size\") or 0)\n    except Exception:\n        return 0\n\n\n# Which mount to prefer depends on the accelerator actually granted, not on the\n# one requested. Version 10 asked for a T4, was given a CPU, and spent 4h25m on\n# the ingest with a ViT-S; a ViT-B/16 costs roughly three times as much per image\n# and would not fit the kernel budget at any resolution on SIZE_LADDER. So when\n# CUDA is absent, fall back to the narrowest mounted checkpoint instead of timing\n# out four hours later.\nif torch.cuda.is_available():\n    EFFECTIVE_FAMILY = BACKBONE_FAMILY\n    EFFECTIVE_HIDDEN = BACKBONE_TARGET_HIDDEN\nelse:\n    EFFECTIVE_FAMILY = \"\"          # empty substring matches every path\n    EFFECTIVE_HIDDEN = 384\n    print(\"No CUDA device visible: preferring the narrowest mounted checkpoint \"\n          \"so the ingest still fits the kernel budget.\", flush=True)\n\n\ndef rank_backbones(candidates):\n    \"\"\"Order mounts by family match, then by declared width against the target.\n\n    Args:\n        candidates: Directories from find_local_backbones().\n\n    Returns:\n        The same directories, best first. Checkpoints that do not declare a\n        hidden size sort last within their family rather than winning by accident.\n    \"\"\"\n    def key(path):\n        low = str(path).lower()\n        family_rank = 0 if EFFECTIVE_FAMILY.lower() in low else 1\n        hidden = checkpoint_hidden_size(path)\n        distance = abs(hidden - EFFECTIVE_HIDDEN) if hidden else 10 ** 6\n        return (family_rank, distance, str(path))\n    return sorted(candidates, key=key)\n\n\ndef load_dino():\n    \"\"\"Load the backbone from a mounted checkpoint. No network access.\n\n    Returns:\n        (model, embed_dim, slice_dim, backbone_id)\n\n    Raises:\n        RuntimeError: when no mount loads, naming what was found and how to\n            attach one, rather than falling through to random features.\n    \"\"\"\n    candidates = find_local_backbones()\n    if candidates:\n        print(f\"Checkpoints mounted: {[str(c) for c in candidates]}\", flush=True)\n\n    # On CPU the narrowest mount is preferred deliberately, so the family check\n    # only applies when an accelerator is present and the large backbone is the\n    # actual intent.\n    if REQUIRE_BACKBONE_FAMILY and DEVICE == \"cuda\" and BACKBONE_FAMILY:\n        matching = [c for c in candidates\n                    if BACKBONE_FAMILY.lower() in str(c).lower()]\n        if not matching:\n            raise RuntimeError(\n                f\"No {BACKBONE_FAMILY} checkpoint is mounted, and this run is \"\n                f\"configured to require one.\\n\"\n                \"  Mounted instead: \"\n                + (\", \".join(str(c) for c in candidates) if candidates else \"nothing\")\n                + \"\\n  Attach hark99/facebookdinov3-vitb16-pretrain-lvd1689m \"\n                  \"(Transformers / default / version 1): config.json and \"\n                  \"model.safetensors at the mount root, 342,662,192 bytes, \"\n                  \"hidden 768, patch 16.\\n\"\n                  \"  Do NOT attach giovannyrodrguez/dinov3 v4: it is a 3.36 GB \"\n                  \"ViT-H+ and would need roughly 59 h to encode this corpus.\\n\"\n                  \"  Set REQUIRE_BACKBONE_FAMILY = False in Cell 2 to run on \"\n                  \"whatever is mounted instead.\")\n\n    errors = []\n    for pick in rank_backbones(candidates):\n        try:\n            model = AutoModel.from_pretrained(str(pick), local_files_only=True)\n            backbone_id = \"/\".join(pick.parts[-4:])\n            print(f\"\\u2713 Backbone loaded: {pick}\", flush=True)\n            model.eval()\n            for prm in model.parameters():\n                prm.requires_grad_(False)\n            model.to(DEVICE)\n            embed_dim = int(getattr(model.config, \"hidden_size\", 384))\n            return model, embed_dim, 3 * embed_dim, backbone_id\n        except Exception as e:\n            errors.append(f\"{pick}: {type(e).__name__}: {e}\")\n\n    # \"Checkpoints found: none\" on its own does not say why. Most Kaggle model\n    # mounts of DINOv3 ship a single raw .pth state dict with no config.json,\n    # which AutoModel cannot read, and that is indistinguishable from an empty\n    # mount unless the actual files are printed. So dump the tree.\n    listing = []\n    for base in (Path(\"/kaggle/input/models\"), Path(\"/kaggle/input\")):\n        if not base.is_dir():\n            continue\n        for root, dirs, files in os.walk(base):\n            dirs[:] = [d for d in dirs if d not in skip_dirs]\n            if len(Path(root).parts) - len(base.parts) > 8:\n                dirs[:] = []\n                continue\n            if files:\n                listing.append(f\"    {root}: {sorted(files)[:6]}\")\n            if len(listing) >= 25:\n                break\n        if len(listing) >= 25:\n            break\n\n    raise RuntimeError(\n        \"No vision backbone could be loaded from a mounted checkpoint.\\n\"\n        \"  A usable mount is a directory holding config.json next to a \"\n        \".safetensors or .bin weights file.\\n\"\n        \"  Checkpoints found: \"\n        + (\", \".join(str(c) for c in candidates) if candidates else \"none\")\n        + (\"\\n  Load errors: \" + \"; \".join(errors) if errors else \"\")\n        + \"\\n  Files actually mounted:\\n\" + (\"\\n\".join(listing) if listing else \"    (nothing)\")\n        + \"\\n  A single .pth file is NOT loadable here: it is a bare state dict \"\n          \"with no architecture config.\\n\"\n          \"  Known-good mounts: giovannyrodrguez/dinov3 (Transformers, v1 = \"\n          \"dinov3-vitb16) or metaresearch/dinov2 (PyTorch, small, v1).\\n\"\n          \"  This notebook never downloads: reruns execute with internet off.\"\n    )\n\n\nDINO_RAW, EMBED_DIM, SLICE_DIM, BACKBONE_ID = load_dino()\nPATCH_PX = int(getattr(DINO_RAW.config, \"patch_size\", 14))\n\n# DINOv3 prepends register tokens after CLS; DINOv2 does not. Pooling h[:, 1:]\n# therefore averages 4 register tokens into the patch statistics on a DINOv3\n# checkpoint. Registers are trained to absorb global information precisely so it\n# stops leaking into patch tokens, so mixing them back into a spatial mean is the\n# opposite of what they are for. The prefix length is read from the config rather\n# than assumed, so the same code is correct for both families.\nN_PREFIX_TOKENS = 1 + int(getattr(DINO_RAW.config, \"num_register_tokens\", 0) or 0)\n\n# Weights are cast once instead of relying on torch.autocast. autocast state is\n# thread-local and DataParallel runs each replica in its own worker thread, so an\n# autocast context in the main thread does not reach them: the replicas would\n# quietly run fp32 and give back none of the expected speedup.\nENCODE_DTYPE = torch.float16 if DEVICE == \"cuda\" else torch.float32\nif ENCODE_DTYPE == torch.float16:\n    DINO_RAW.half()\n\nN_GPUS = torch.cuda.device_count() if DEVICE == \"cuda\" else 0\nif USE_DATA_PARALLEL and N_GPUS > 1:\n    DINO = nn.DataParallel(DINO_RAW)\n    ENCODE_BATCH = GPU_BATCH_SIZE * N_GPUS\n    print(f\"Encoding across {N_GPUS} GPUs (DataParallel), batch {ENCODE_BATCH}\", flush=True)\nelse:\n    DINO = DINO_RAW\n    ENCODE_BATCH = GPU_BATCH_SIZE\nprint(f\"Backbone: {BACKBONE_ID} | {sum(p.numel() for p in DINO_RAW.parameters()):,} params | \"\n      f\"hidden {EMBED_DIM} | patch {PATCH_PX} | \"\n      f\"{IMG_SIZE // PATCH_PX}x{IMG_SIZE // PATCH_PX} tokens at {IMG_SIZE} px \"\n      f\"({FOV_MM / (IMG_SIZE / PATCH_PX):.2f} mm per token) | \"\n      f\"{N_PREFIX_TOKENS - 1} register tokens excluded from patch pooling\", flush=True)\n\ndef encode_batch(imgs_uint8, batch_size=None):\n    if imgs_uint8.shape[0] == 0:\n        return np.zeros((0, SLICE_DIM), dtype=np.float32)\n\n    bs = batch_size or ENCODE_BATCH\n    x = imgs_uint8.astype(np.float32) / 255.0\n    x = (x - IMAGENET_MEAN) / IMAGENET_STD\n\n    vecs = []\n    with torch.no_grad():\n        for start in range(0, x.shape[0], bs):\n            chunk = torch.from_numpy(x[start:start + bs]).to(DEVICE).to(ENCODE_DTYPE)\n            outputs = DINO(pixel_values=chunk)\n\n            h = outputs.last_hidden_state.float()\n            cls_t = h[:, 0]\n            p_mean = h[:, N_PREFIX_TOKENS:].mean(dim=1)\n            p_max = h[:, N_PREFIX_TOKENS:].amax(dim=1)\n            vec = torch.cat([cls_t, p_mean, p_max], dim=1)\n            vecs.append(vec.cpu().numpy())\n            \n    # Clear CUDA memory buffer to prevent GPU OOM\n    if DEVICE == \"cuda\":\n        torch.cuda.empty_cache()\n    return np.concatenate(vecs, axis=0)\n\nprint(f\"Backbone setup: Embed Dimension={EMBED_DIM}, Slice Feature Dimension={SLICE_DIM}\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"b9100551","metadata":{},"outputs":[],"source":"# Cell 7: Decoupled Multi-Threaded Study Ingestion & Feature Pooling (Zero CUDA Lock)\ndef assign_slots_for_study(series_df, study_uid):\n    sub = series_df[series_df[\"StudyInstanceUID\"] == study_uid]\n    slots = {}\n    used = set()\n    for name, plane, flag, _ in SLOT_SPEC:\n        cands = sub[(sub[\"plane\"] == plane) & (sub[\"fluid\"] == flag) & (~sub[\"SeriesInstanceUID\"].isin(used))]\n        if len(cands) > 0:\n            sid = cands.iloc[0][\"SeriesInstanceUID\"]\n            slots[name] = sid\n            used.add(sid)\n        else:\n            slots[name] = None\n    return slots\n\ndef decode_and_encode_study(study_uid, dcm_root, series_df):\n    slots = assign_slots_for_study(series_df, study_uid)\n    \n    # Step 1: Parallel CPU Decoding (CPU_WORKERS threads)\n    def decode_slot(args):\n        name, plane, _, k = args\n        series_id = slots.get(name)\n        if series_id is None or dcm_root is None:\n            return name, np.zeros((0, 3, IMG_SIZE, IMG_SIZE), dtype=np.uint8)\n            \n        folder = Path(dcm_root) / str(study_uid) / str(series_id)\n        if not folder.is_dir():\n            return name, np.zeros((0, 3, IMG_SIZE, IMG_SIZE), dtype=np.uint8)\n            \n        imgs = build_slot_stack_fast(folder, plane, k, IMG_SIZE)\n        return name, imgs\n\n    with ThreadPoolExecutor(max_workers=min(CPU_WORKERS, os.cpu_count() or 4)) as executor:\n        slot_imgs = dict(executor.map(decode_slot, SLOT_SPEC))\n        \n    # Step 2: Main-Thread GPU Encoding (100% thread-safe)\n    slot_embs = {}\n    for name, _, _, _ in SLOT_SPEC:\n        imgs = slot_imgs.get(name)\n        if imgs is not None and imgs.shape[0] > 0:\n            slot_embs[name] = encode_batch(imgs)\n        else:\n            slot_embs[name] = np.zeros((0, SLICE_DIM), dtype=np.float32)\n            \n    return slot_embs\n\ndef pool_study_vector(slot_embs):\n    parts = []\n    for name, _, _, _ in SLOT_SPEC:\n        emb = slot_embs.get(name)\n        if emb is None or emb.shape[0] == 0:\n            parts.append(np.zeros(EMBED_DIM * 3, dtype=np.float32))\n            parts.append(np.zeros(EMBED_DIM * 3, dtype=np.float32))\n        else:\n            parts.append(emb.mean(axis=0))\n            parts.append(emb.max(axis=0))\n    return np.concatenate(parts)\n\nprint(\"✓ Decoupled multi-threaded study ingestion ready.\", flush=True)\n"},{"cell_type":"code","id":"cf829ffa","source":"# Cell 7b: Encode-Cost Probe & Resolution Downshift\n#\n# This competition awards a separate efficiency prize, so encode cost is a scored\n# axis and not merely a timeout risk. It is also a real timeout risk here: the\n# visible test split is 3 studies but the rerun scores against a hidden split of\n# unknown size, and a DINOv3 ViT-B/16 costs roughly three times a DINOv2 ViT-S/14\n# per study: measured at 1.29 s/study for the ViT-S at 336 px, a ViT-B projects to\n# about 3.9 s/study, or 6.8 h on one GPU before any downshift.\n#\n# So the cost is measured on a few studies and extrapolated, rather than assumed.\n# If the projection exceeds the kernel budget the input size steps down the\n# ladder and is re-measured. Every number printed below except the assumed rerun\n# size comes from this run.\n\ndef probe_encode_cost(uids, dcm_root, series_df, size, n_probe):\n    \"\"\"Measure seconds and images per study at a given input size.\n\n    Args:\n        uids: Study ids to draw the probe sample from.\n        dcm_root: DICOM root, or None.\n        series_df: Series manifest.\n        size: Input size to measure at; IMG_SIZE is set to this for the probe.\n        n_probe: How many studies to time.\n\n    Returns:\n        (seconds_per_study, images_per_study, n_measured).\n    \"\"\"\n    global IMG_SIZE\n    if dcm_root is None or not uids:\n        return 0.0, 0.0, 0\n\n    prev = IMG_SIZE\n    IMG_SIZE = size\n    try:\n        sample = uids[:max(1, n_probe)]\n        t0 = time.time()\n        n_images = 0\n        for uid in sample:\n            slot_embs = decode_and_encode_study(uid, dcm_root, series_df)\n            n_images += sum(int(e.shape[0]) for e in slot_embs.values())\n        elapsed = time.time() - t0\n    finally:\n        IMG_SIZE = prev\n\n    n = len(sample)\n    return elapsed / max(1, n), n_images / max(1, n), n\n\n\nN_TEST_PLANNED = max(len(test_ids), ASSUMED_RERUN_TEST_STUDIES)\nENCODE_S_PER_STUDY = 0.0\nPROBE_MEASURED = False\n\n# The checkpoint that actually loaded. Several mounts can match BACKBONE_FAMILY\n# and the ranking picks one, so the cache tag has to name what ran.\nBACKBONE_NAME = str(BACKBONE_ID or BACKBONE_FAMILY).strip(\"/\").replace(\"/\", \"-\") or \"backbone\"\nBACKBONE_PARAMS = sum(p.numel() for p in DINO_RAW.parameters())\n\nif TRAIN_DCM_DIR is None:\n    print(\"No train DICOM directory; skipping the encode probe.\", flush=True)\nelse:\n    print(f\"Probing encode cost on {PROBE_STUDIES} studies \"\n          f\"(backbone {BACKBONE_NAME}, {BACKBONE_PARAMS:,} params)...\", flush=True)\n\n    ladder = [s for s in SIZE_LADDER if s <= IMG_SIZE] or [IMG_SIZE]\n    chosen = None\n    for size in ladder:\n        s_per_study, imgs_per_study, n_meas = probe_encode_cost(\n            train_ids, TRAIN_DCM_DIR, train_series, size, PROBE_STUDIES)\n        if n_meas == 0:\n            break\n        encode_s = (len(train_ids) + N_TEST_PLANNED) * s_per_study\n        total_s = encode_s + RESERVE_S\n\n        print(f\"  {size} px: {s_per_study:.2f} s/study, {imgs_per_study:.1f} images/study \"\n              f\"(measured on {n_meas})\", flush=True)\n        print(f\"    projected encode {encode_s:.0f} s over {len(train_ids)} train + \"\n              f\"{N_TEST_PLANNED} test, +{RESERVE_S} s reserve = {total_s:.0f} s \"\n              f\"of {KERNEL_BUDGET_S} s budget\", flush=True)\n\n        ENCODE_S_PER_STUDY = s_per_study\n        PROBE_MEASURED = True\n        if total_s <= KERNEL_BUDGET_S:\n            chosen = size\n            break\n\n    if chosen is None and PROBE_MEASURED:\n        chosen = ladder[-1]\n        print(f\"  WARNING: even {chosen} px exceeds the budget. Running at {chosen} px \"\n              f\"anyway. Cut ASSUMED_RERUN_TEST_STUDIES if you believe the hidden \"\n              f\"split is smaller, drop slice budgets in SLOT_SPEC, or move to a \"\n              f\"smaller backbone.\", flush=True)\n\n    if chosen is not None and chosen != IMG_SIZE:\n        print(f\"  downshifting input size {IMG_SIZE} -> {chosen} px\", flush=True)\n        IMG_SIZE = chosen\n    elif chosen is not None:\n        print(f\"  no downshift needed; staying at {IMG_SIZE} px\", flush=True)\n\n# Cache identity. Without this the memmap filename was identical across every\n# configuration, so a rerun after changing the size or the backbone silently\n# reloaded features built under the old one.\nCACHE_TAG = \"{}_{}px_slots{}_{}\".format(\n    BACKBONE_NAME.replace(\".\", \"-\"),\n    IMG_SIZE,\n    \"-\".join(str(k) for *_rest, k in SLOT_SPEC),\n    SLICE_DIM,\n)\nprint(f\"Cache tag: {CACHE_TAG}\", flush=True)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"id":"080341c5","metadata":{},"outputs":[],"source":"# Cell 8: Feature Extraction & Memory-Mapped Array Cache\ndef get_features(uids, dcm_root, series_df, split_name):\n    # The tag is part of the filename so a cache built at another input size or\n    # with another backbone cannot be reloaded as if it matched this run.\n    memmap_path = WORK / f\"emb_{split_name}_{CACHE_TAG}.npy\"\n    if memmap_path.is_file():\n        print(f\"Loading cached {split_name} features from {memmap_path}...\", flush=True)\n        return np.load(memmap_path)\n        \n    print(f\"Ingesting {len(uids)} {split_name} studies...\", flush=True)\n    t0 = time.time()\n    feature_dim = N_SLOTS * 2 * SLICE_DIM\n    \n    # Memory-mapped array allocation (float32, low memory footprint)\n    mmap_arr = np.lib.format.open_memmap(memmap_path, mode=\"w+\", dtype=np.float32, shape=(len(uids), feature_dim))\n    \n    for i, uid in enumerate(uids):\n        # Unbuffered real-time logging every 100 studies\n        if (i + 1) % 100 == 0 or i == len(uids) - 1:\n            rate = (i + 1) / (time.time() - t0)\n            eta = (len(uids) - i - 1) / max(rate, 1e-6)\n            print(f\"  [{split_name}] {i+1}/{len(uids)} ({rate:.2f} std/s | ETA: {eta:.0f}s)\", flush=True)\n            \n        if dcm_root is None:\n            mmap_arr[i] = np.zeros(feature_dim, dtype=np.float32)\n        else:\n            slot_embs = decode_and_encode_study(uid, dcm_root, series_df)\n            vec = pool_study_vector(slot_embs)\n            mmap_arr[i] = vec\n            \n    mmap_arr.flush()\n    gc.collect()\n    print(f\"✓ {split_name} features cached to {memmap_path} ({mmap_arr.shape}) \"\n          f\"in {time.time() - t0:.0f}s\", flush=True)\n    return mmap_arr\n\nX_train = get_features(train_ids, TRAIN_DCM_DIR, train_series, \"train\")\nX_test = get_features(test_ids, TEST_DCM_DIR, test_series, \"test\")\n"},{"cell_type":"code","execution_count":null,"id":"d33d8aa6","metadata":{},"outputs":[],"source":"# Cell 9: Grouped Cross-Validation & Ridge Baseline Probe\n#\n# Two leakage channels are closed here with one union-find pass:\n#\n#   1. Duplicate report text. 49 texts are shared by 183 studies in this corpus,\n#      largest block 37. Every study in a block carries identical targets, so a\n#      block split across folds validates on a target whose source was trained on.\n#   2. Acquisition site. Ungrouped splits are reported to inflate AUC by roughly\n#      0.05 through scanner memorisation. Site grouping is applied only when the\n#      series manifest actually exposes a site-like column and only when it leaves\n#      enough distinct groups to still cut N_FOLDS folds.\n#\n# PatientID is 1:1 with StudyInstanceUID across all training studies, so there is\n# no patient grouping to add.\n\nreport_by_study = {sid: txt.strip() if isinstance(txt, str) else \"\" for sid, txt in zip(train_ids, train[report_col])} if report_col else {}\n\nSITE_COL = None\nif GROUP_MODE in (\"auto\", \"site\") and len(train_series) > 0:\n    for _c in SITE_COLUMN_CANDIDATES:\n        if _c in train_series.columns:\n            SITE_COL = _c\n            break\n\nsite_of = {}\nif SITE_COL is not None:\n    _m = (train_series.drop_duplicates(\"StudyInstanceUID\")\n                      .set_index(\"StudyInstanceUID\")[SITE_COL])\n    site_of = {str(k): str(v).strip() for k, v in _m.items()\n               if pd.notna(v) and str(v).strip()}\n\n\ndef build_groups(use_site):\n    \"\"\"Union-find grouping of studies by shared report text and, optionally, site.\n\n    Args:\n        use_site: Whether to also merge studies sharing an acquisition site.\n\n    Returns:\n        Integer group id per study, aligned to train_ids.\n    \"\"\"\n    n = len(train_ids)\n    parent = list(range(n))\n\n    def find(a):\n        while parent[a] != a:\n            parent[a] = parent[parent[a]]\n            a = parent[a]\n        return a\n\n    def union(a, b):\n        ra, rb = find(a), find(b)\n        if ra != rb:\n            parent[max(ra, rb)] = min(ra, rb)\n\n    buckets = {}\n    for i, sid in enumerate(train_ids):\n        txt = report_by_study.get(sid, \"\")\n        if txt:\n            buckets.setdefault((\"txt\", txt), []).append(i)\n    if use_site:\n        for i, sid in enumerate(train_ids):\n            s = site_of.get(sid)\n            if s:\n                buckets.setdefault((\"site\", s), []).append(i)\n\n    for idxs in buckets.values():\n        for j in idxs[1:]:\n            union(idxs[0], j)\n\n    return np.asarray([find(i) for i in range(n)])\n\n\ngroups = build_groups(use_site=False)\nGROUPING = \"report\"\n_n_report = len(set(groups.tolist()))\n\nif SITE_COL is not None:\n    _g_site = build_groups(use_site=True)\n    _n_site = len(set(_g_site.tolist()))\n    if _n_site >= N_FOLDS:\n        groups, GROUPING = _g_site, f\"report+site({SITE_COL})\"\n    else:\n        print(f\"Site column '{SITE_COL}' collapses the data to {_n_site} groups, \"\n              f\"fewer than N_FOLDS={N_FOLDS}; keeping report-only grouping.\", flush=True)\nelif GROUP_MODE == \"site\":\n    print(f\"GROUP_MODE='site' but no column in {SITE_COLUMN_CANDIDATES} exists in \"\n          f\"the series manifest; falling back to report grouping.\", flush=True)\n\n_sizes = pd.Series(groups).value_counts()\n_n_blocks = int((_sizes > 1).sum())\nprint(f\"Fold grouping ({GROUPING}): {len(train_ids)} studies -> {len(_sizes)} groups\", flush=True)\nprint(f\"  {_n_blocks} multi-study blocks covering {int(_sizes[_sizes > 1].sum())} studies, \"\n      f\"largest block {int(_sizes.iloc[0])}\", flush=True)\nif SITE_COL is None:\n    print(\"  no site/scanner column found in the series manifest; report text only\", flush=True)\n\n_n_splits = max(2, min(N_FOLDS, len(_sizes)))\ngkf = GroupKFold(n_splits=_n_splits)\nFOLDS = list(gkf.split(X_train, Y_soft.values, groups=groups))\nfor _k, (_tr, _va) in enumerate(FOLDS):\n    _shared = set(groups[_tr].tolist()) & set(groups[_va].tolist())\n    assert not _shared, f\"fold {_k} leaks {len(_shared)} groups across the split\"\nprint(f\"  {len(FOLDS)} folds, no group appears on both sides of any split\", flush=True)\n\n\n# One SVD + standardisation per fold, computed once and shared by every linear\n# probe below. The earlier version refitted a 512-component TruncatedSVD inside\n# each model, so adding a second probe would have doubled the most expensive step\n# in this cell for nothing. The transform itself is unchanged: same seed, same\n# scaler, fitted on the fold's training rows only.\ndef build_fold_matrices():\n    \"\"\"Per-fold design matrices in SVD-then-standardise space.\n\n    Returns:\n        List of (X_tr, X_va, X_te) tuples aligned to FOLDS.\n    \"\"\"\n    mats = []\n    for tr_idx, val_idx in FOLDS:\n        X_tr = X_train[tr_idx]\n        comps = min(SVD_COMPONENTS, min(X_tr.shape) - 1)\n        svd = TruncatedSVD(n_components=comps, random_state=SEED)\n        X_tr_svd = svd.fit_transform(X_tr)\n        scaler = StandardScaler()\n        mats.append((\n            scaler.fit_transform(X_tr_svd),\n            scaler.transform(svd.transform(X_train[val_idx])),\n            scaler.transform(svd.transform(X_test)),\n        ))\n    return mats\n\n\nFOLD_MATS = build_fold_matrices()\nFOLD_WEIGHTS = [np.where(np.isin(tr_idx, GOLD_POS), GOLD_WEIGHT, 1.0)\n                for tr_idx, _ in FOLDS]\nprint(f\"  fold design matrices: {FOLD_MATS[0][0].shape[1]} SVD components, \"\n      f\"{len(GOLD_POS)} gold rows weighted x{GOLD_WEIGHT:g}\", flush=True)\n\n\ndef calc_macro_auc(y_true, y_pred):\n    aucs = []\n    for j in range(N_LABELS):\n        yt = y_true[:, j]\n        if len(np.unique(yt >= 0.5)) > 1:\n            aucs.append(roc_auc_score((yt >= 0.5).astype(int), y_pred[:, j]))\n    return float(np.mean(aucs)) if aucs else 0.0\n\n\ndef train_ridge(Y_tr_all):\n    \"\"\"Ridge probe with a per-label penalty.\n\n    RidgeCV picks a single alpha for the whole target block by default. The\n    twelve findings differ by nearly an order of magnitude in prevalence and in\n    extractor precision, so alpha_per_target lets a rare, noisy label such as\n    Lateral OA (precision 0.36) be regularised harder than ACL (0.85).\n\n    Returns:\n        (oof, test) probability matrices.\n    \"\"\"\n    oof = np.zeros((len(X_train), N_LABELS), dtype=np.float32)\n    te = np.zeros((len(X_test), N_LABELS), dtype=np.float32)\n    alphas = np.logspace(0, 5, 11)\n\n    for k, (tr_idx, val_idx) in enumerate(FOLDS):\n        X_tr_sc, X_va_sc, X_te_sc = FOLD_MATS[k]\n        ridge = RidgeCV(alphas=alphas, alpha_per_target=True)\n        ridge.fit(X_tr_sc, Y_tr_all[tr_idx], sample_weight=FOLD_WEIGHTS[k])\n        oof[val_idx] = np.clip(ridge.predict(X_va_sc), 0.0, 1.0)\n        te += np.clip(ridge.predict(X_te_sc), 0.0, 1.0) / len(FOLDS)\n\n    return oof, te\n\n\n# Logistic probe. Ridge regresses onto the soft target under a squared loss,\n# which is not the loss a ranking metric rewards; a per-label logistic fit on the\n# same components optimises a proper scoring rule instead. Linear probes on\n# frozen self-supervised features are the standard readout in this setting and\n# match or beat fine-tuned supervised baselines on radiology benchmarks\n# (arXiv:2312.02366), so it is worth having both. It is a candidate, not a\n# replacement: Cell 11 keeps it only if it earns its place.\nLOGIT_C = 0.05\n\n\ndef train_logistic(Y_tr_all):\n    \"\"\"Per-label L2 logistic probe on the shared fold matrices.\n\n    Returns:\n        (oof, test) probability matrices.\n    \"\"\"\n    from sklearn.linear_model import LogisticRegression\n\n    oof = np.zeros((len(X_train), N_LABELS), dtype=np.float32)\n    te = np.zeros((len(X_test), N_LABELS), dtype=np.float32)\n\n    for k, (tr_idx, val_idx) in enumerate(FOLDS):\n        X_tr_sc, X_va_sc, X_te_sc = FOLD_MATS[k]\n        w = FOLD_WEIGHTS[k]\n        for j in range(N_LABELS):\n            yj = (Y_tr_all[tr_idx][:, j] >= 0.5).astype(int)\n            if len(np.unique(yj)) < 2:\n                # One class in this fold means there is no boundary to fit; the\n                # base rate is the only honest prediction.\n                base = float(yj.mean())\n                oof[val_idx, j] = base\n                te[:, j] += base / len(FOLDS)\n                continue\n            clf = LogisticRegression(C=LOGIT_C, max_iter=3000)\n            clf.fit(X_tr_sc, yj, sample_weight=w)\n            oof[val_idx, j] = clf.predict_proba(X_va_sc)[:, 1]\n            te[:, j] += clf.predict_proba(X_te_sc)[:, 1] / len(FOLDS)\n\n    return oof, te\n\n\noof_ridge, test_ridge = train_ridge(Y_soft.values)\nauc_ridge = calc_macro_auc(Y_soft.values, oof_ridge)\nauc_ridge_gold, n_gold_lab = gold_macro_auc(oof_ridge)\nprint(f\"Ridge probe | self-consistency AUC {auc_ridge:.4f} \"\n      f\"(vs its own extractor's labels, NOT a leaderboard estimate)\", flush=True)\nprint(f\"            | gold AUC {auc_ridge_gold:.4f} over {n_gold_lab} labels \"\n      f\"on {len(GOLD_POS)} radiologist-labelled studies <- track this one\", flush=True)\n\n# This cell sits four and a half hours downstream of the DICOM ingest, so a probe\n# that fails costs its own slot in the ensemble rather than the whole run.\ntry:\n    oof_logit, test_logit = train_logistic(Y_soft.values)\n    auc_logit = calc_macro_auc(Y_soft.values, oof_logit)\n    auc_logit_gold, _ = gold_macro_auc(oof_logit)\n    LOGIT_OK = True\n    print(f\"Logistic probe | self-consistency AUC {auc_logit:.4f}\", flush=True)\n    print(f\"               | gold AUC {auc_logit_gold:.4f} <- track this one\", flush=True)\nexcept Exception as _e:\n    LOGIT_OK = False\n    oof_logit = test_logit = None\n    auc_logit = auc_logit_gold = float(\"nan\")\n    print(f\"Logistic probe skipped: {type(_e).__name__}: {_e}\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"d495af27","metadata":{},"outputs":[],"source":"# Cell 10: Slot-Attention MIL Head & Asymmetric Focal Loss\n#\n# The previous head flattened the 11,520-dimensional study vector into a plain\n# MLP and scored 0.5634 gold AUC, which is chance. The flattening is the problem.\n# That vector is not 11,520 unrelated numbers: it is 2 * N_SLOTS pooled\n# descriptors of SLICE_DIM dimensions each, one mean and one max per\n# plane/sequence slot, and a dense first layer has to rediscover that block\n# structure from 4,407 rows of noisy targets.\n#\n# Attention-based multiple-instance learning (Ilse, Tomczak & Welling, 2018)\n# replaces the flattening with a learned convex combination over those\n# descriptors. The knee-MRI literature arrives at the same place from the other\n# direction: CoPAS (Nat Commun 2024, the closest published task to this one --\n# twelve knee findings, AUC 0.812) attends across planes and sequences instead of\n# concatenating them, and the multi-task system in eClinicalMedicine 2025 attends\n# to a discriminative region per abnormality rather than one shared region. So\n# the attention map here is per label: a meniscal tear and a Baker's cyst should\n# not be forced to read the same slot with the same weight.\n#\n# Two details that are easy to get wrong:\n#   1. Absent slots are written as exact zeros by pool_study_vector. They are\n#      masked out of the softmax rather than left to be standardised into small\n#      non-zero values the attention would then have to learn to ignore.\n#   2. Early stopping runs on an inner split carved out of the fold's own\n#      training groups, never on the outer validation rows. Stopping on the outer\n#      rows would make the OOF matrix optimistic, and Cell 11 chooses ensemble\n#      weights from exactly that matrix.\n\nN_TOKENS = 2 * N_SLOTS      # one mean and one max descriptor per slot\nHEAD_BATCH = 64\nHEAD_ATTN_DIM = 128\nINNER_VAL_FRAC = 0.15\n\n\nclass AsymmetricFocalLoss(nn.Module):\n    def __init__(self, gamma_pos=0.0, gamma_neg=2.0, clip=0.05):\n        super().__init__()\n        self.gamma_pos = gamma_pos\n        self.gamma_neg = gamma_neg\n        self.clip = clip\n\n    def forward(self, logits, targets, weights=None):\n        probs = torch.sigmoid(logits)\n\n        pos_loss = targets * torch.log(probs.clamp(min=1e-8))\n        if self.gamma_pos > 0:\n            pos_loss = pos_loss * ((1 - probs) ** self.gamma_pos)\n\n        neg_probs = (probs - self.clip).clamp(min=0.0)\n        neg_loss = (1 - targets) * torch.log((1 - neg_probs).clamp(min=1e-8))\n        if self.gamma_neg > 0:\n            neg_loss = neg_loss * ((1 - neg_probs) ** self.gamma_neg)\n\n        loss = -(pos_loss + neg_loss)\n        if weights is not None:\n            loss = loss * weights\n        return loss.mean()\n\n\nclass SlotAttentionMILHead(nn.Module):\n    \"\"\"Gated attention MIL over the pooled slot descriptors, one map per label.\"\"\"\n\n    def __init__(self, slice_dim, n_tokens, num_labels=N_LABELS,\n                 hidden=HEAD_HIDDEN, dropout=HEAD_DROPOUT, attn_dim=HEAD_ATTN_DIM):\n        super().__init__()\n        self.n_tokens = n_tokens\n        self.slice_dim = slice_dim\n        self.proj = nn.Sequential(\n            nn.Linear(slice_dim, hidden),\n            nn.LayerNorm(hidden),\n            nn.SiLU(),\n            nn.Dropout(dropout),\n        )\n        self.att_v = nn.Linear(hidden, attn_dim)\n        self.att_u = nn.Linear(hidden, attn_dim)\n        self.att_w = nn.Linear(attn_dim, num_labels)\n        # A learned per-slot, per-label prior. Sagittal fluid-sensitive slots are\n        # informative for effusion in almost every study; the model should be\n        # able to encode that without inferring it from the descriptor each time.\n        self.slot_prior = nn.Parameter(torch.zeros(n_tokens, num_labels))\n        self.classifier = nn.Parameter(torch.empty(num_labels, hidden))\n        nn.init.normal_(self.classifier, std=0.02)\n        self.bias = nn.Parameter(torch.zeros(num_labels))\n        self.drop = nn.Dropout(dropout)\n\n    def forward(self, x, mask):\n        b = x.shape[0]\n        h = self.proj(x.view(b, self.n_tokens, self.slice_dim))\n        gate = torch.tanh(self.att_v(h)) * torch.sigmoid(self.att_u(h))\n        scores = self.att_w(gate) + self.slot_prior.unsqueeze(0)\n        scores = scores.masked_fill(~mask.unsqueeze(-1), float(\"-inf\"))\n        alpha = torch.softmax(scores, dim=1)\n        pooled = torch.einsum(\"btl,bth->blh\", alpha, h)\n        pooled = self.drop(pooled)\n        return (pooled * self.classifier.unsqueeze(0)).sum(-1) + self.bias\n\n\ndef slot_mask(X):\n    \"\"\"True where a slot descriptor carries real data.\n\n    Args:\n        X: Raw, unstandardised study matrix (n, N_TOKENS * SLICE_DIM).\n\n    Returns:\n        Boolean array (n, N_TOKENS).\n    \"\"\"\n    m = np.abs(X.reshape(len(X), N_TOKENS, SLICE_DIM)).sum(axis=2) > 0\n    # A study with no usable series at all would softmax a row of -inf and return\n    # NaN. Leaving its first slot open degrades it to a constant prediction.\n    m[~m.any(axis=1), 0] = True\n    return m\n\n\ndef inner_split(tr_idx, seed):\n    \"\"\"Hold out whole groups from the fold's training rows for early stopping.\n\n    Splitting by group rather than by row keeps duplicate-report blocks intact,\n    the same leak Cell 9 closes for the outer folds.\n\n    Args:\n        tr_idx: Training row indices for one outer fold.\n        seed: Seed for the group draw.\n\n    Returns:\n        (fit_idx, stop_idx) row-index arrays.\n    \"\"\"\n    g = groups[tr_idx]\n    uniq = np.unique(g)\n    rng = np.random.default_rng(seed)\n    n_hold = max(1, int(round(INNER_VAL_FRAC * len(uniq))))\n    hold = set(rng.choice(uniq, size=min(n_hold, len(uniq) - 1), replace=False).tolist())\n    is_hold = np.array([bool(x in hold) for x in g])\n    return tr_idx[~is_hold], tr_idx[is_hold]\n\n\ndef train_advanced_head(X_tr_all, Y_tr_all, X_te_all):\n    \"\"\"Fit the slot-attention head over folds x seeds.\n\n    Returns:\n        (oof, test) probability matrices.\n    \"\"\"\n    oof = np.zeros((len(X_tr_all), N_LABELS), dtype=np.float32)\n    te = np.zeros((len(X_te_all), N_LABELS), dtype=np.float32)\n    criterion = AsymmetricFocalLoss(gamma_pos=0.0, gamma_neg=2.0)\n\n    mask_tr_all = slot_mask(X_tr_all)\n    mask_te_all = slot_mask(X_te_all)\n    n_runs = len(FOLDS) * len(HEAD_SEEDS)\n    stopped_at = []\n\n    for k, (tr_idx, val_idx) in enumerate(FOLDS):\n        # Standardise the fold's training block once and address the inner split\n        # by position into it. Re-transforming per seed would have copied roughly\n        # 170 MB three times per fold for no change in the result.\n        scaler = StandardScaler()\n        X_tr_sc = scaler.fit_transform(X_tr_all[tr_idx])\n        pos_of = {int(r): i for i, r in enumerate(tr_idx)}\n        X_va_sc = scaler.transform(X_tr_all[val_idx])\n        X_te_sc = scaler.transform(X_te_all)\n\n        va_x = torch.tensor(X_va_sc, dtype=torch.float32).to(DEVICE)\n        te_x = torch.tensor(X_te_sc, dtype=torch.float32).to(DEVICE)\n        va_m = torch.tensor(mask_tr_all[val_idx]).to(DEVICE)\n        te_m = torch.tensor(mask_te_all).to(DEVICE)\n\n        fold_oof = np.zeros((len(val_idx), N_LABELS), dtype=np.float64)\n\n        for seed in HEAD_SEEDS:\n            fit_idx, stop_idx = inner_split(tr_idx, seed)\n            fit_pos = np.array([pos_of[int(r)] for r in fit_idx])\n            stop_pos = np.array([pos_of[int(r)] for r in stop_idx])\n            y_fit = Y_tr_all[fit_idx]\n            y_stop = Y_tr_all[stop_idx]\n            w_fit = np.where(np.isin(fit_idx, GOLD_POS), GOLD_WEIGHT, 1.0)\n            m_fit = mask_tr_all[fit_idx]\n            stop_x = torch.tensor(X_tr_sc[stop_pos], dtype=torch.float32).to(DEVICE)\n            stop_m = torch.tensor(mask_tr_all[stop_idx]).to(DEVICE)\n\n            torch.manual_seed(seed)\n            rng = np.random.default_rng(seed)\n            model = SlotAttentionMILHead(SLICE_DIM, N_TOKENS).to(DEVICE)\n\n            # Prior-bias initialisation on the fitting rows' base rates, so the\n            # first epochs do not spend themselves learning the class priors.\n            p_base = np.clip(y_fit.mean(axis=0), 1e-4, 1.0 - 1e-4)\n            with torch.no_grad():\n                model.bias.copy_(torch.tensor(np.log(p_base / (1.0 - p_base)),\n                                              dtype=torch.float32))\n\n            opt = torch.optim.AdamW(model.parameters(), lr=HEAD_LR,\n                                    weight_decay=HEAD_WD)\n            sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=HEAD_EPOCHS)\n\n            best_auc, best_epoch = -1.0, 0\n            best_va = np.zeros((len(val_idx), N_LABELS), dtype=np.float32)\n            best_te = np.zeros((len(X_te_all), N_LABELS), dtype=np.float32)\n\n            for epoch in range(HEAD_EPOCHS):\n                model.train()\n                perm = rng.permutation(len(fit_pos))\n                for start in range(0, len(perm), HEAD_BATCH):\n                    idx = perm[start:start + HEAD_BATCH]\n                    if len(idx) < 2:\n                        continue        # a 1-row batch is a no-op worth skipping\n                    xb = torch.tensor(X_tr_sc[fit_pos[idx]], dtype=torch.float32).to(DEVICE)\n                    yb = torch.tensor(y_fit[idx], dtype=torch.float32).to(DEVICE)\n                    wb = torch.tensor(w_fit[idx], dtype=torch.float32).unsqueeze(1).to(DEVICE)\n                    mb = torch.tensor(m_fit[idx]).to(DEVICE)\n\n                    opt.zero_grad()\n                    loss = criterion(model(xb, mb), yb, weights=wb)\n                    loss.backward()\n                    opt.step()\n                sched.step()\n\n                model.eval()\n                with torch.no_grad():\n                    stop_p = torch.sigmoid(model(stop_x, stop_m)).cpu().numpy()\n                auc = calc_macro_auc(y_stop, stop_p)\n                if auc > best_auc:\n                    best_auc, best_epoch = auc, epoch + 1\n                    with torch.no_grad():\n                        best_va = torch.sigmoid(model(va_x, va_m)).cpu().numpy()\n                        best_te = torch.sigmoid(model(te_x, te_m)).cpu().numpy()\n\n            stopped_at.append(best_epoch)\n            fold_oof += best_va / len(HEAD_SEEDS)\n            te += best_te / n_runs\n\n        oof[val_idx] = fold_oof\n\n    print(f\"  slot-attention head: {n_runs} fits, early stop at epoch \"\n          f\"{np.mean(stopped_at):.1f} of {HEAD_EPOCHS} on average\", flush=True)\n    return oof, te\n\n\n# The head is fitted on the weak targets, so its self-consistency number rises\n# with extractor agreement rather than with diagnostic accuracy. Both are\n# printed; only the gold column tracks the leaderboard metric.\ntry:\n    oof_nn, test_nn = train_advanced_head(X_train, Y_soft.values, X_test)\n    auc_nn = calc_macro_auc(Y_soft.values, oof_nn)\n    auc_nn_gold, _ = gold_macro_auc(oof_nn)\n    HEAD_OK = True\n    print(f\"Slot-attention head | self-consistency AUC {auc_nn:.4f} \"\n          f\"(not a leaderboard estimate)\", flush=True)\n    print(f\"                    | gold AUC {auc_nn_gold:.4f} on \"\n          f\"{len(GOLD_POS)} studies <- track this one\", flush=True)\nexcept Exception as _e:\n    HEAD_OK = False\n    oof_nn = test_nn = None\n    auc_nn = auc_nn_gold = float(\"nan\")\n    print(f\"Slot-attention head skipped: {type(_e).__name__}: {_e}\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"0ef12aab","metadata":{},"outputs":[],"source":"# Cell 11: Candidate Selection & Rank-Averaged Ensemble\n#\n# The previous version blended ridge and the neural head at a fixed 50/50. On the\n# last run that was strictly harmful: ridge scored 0.7115 gold AUC on its own,\n# the head 0.5634, and the fixed average landed at 0.6614 -- the submitted file\n# was worse than its own best component. A constant weight cannot notice that.\n#\n# Weights are chosen here by greedy forward selection with replacement (Caruana\n# et al., 2004): start from the best single model, then repeatedly add whichever\n# candidate most improves the averaged score, allowing the same model to be added\n# again so the final weights are integer counts. A candidate that improves\n# nothing is never added, so a chance-level model contributes exactly zero\n# instead of dragging the ensemble toward it.\n#\n# Selection is scored on the full 4,407-study self-consistency AUC, not on the 58\n# gold studies. n=58 is far too small to fit an ensemble weight against without\n# overfitting the very number used to judge the result. The gold column is\n# printed for every candidate and for the final blend, but it never drives the\n# choice.\n\ndef rank01(arr):\n    out = np.zeros_like(arr, dtype=np.float64)\n    for j in range(arr.shape[1]):\n        r = rankdata(arr[:, j], method=\"average\")\n        out[:, j] = (r - 1.0) / max(1.0, len(r) - 1.0)\n    return out\n\n\nCANDIDATES = {\"ridge\": (oof_ridge, test_ridge)}\nif LOGIT_OK:\n    CANDIDATES[\"logistic\"] = (oof_logit, test_logit)\nif HEAD_OK:\n    CANDIDATES[\"slot_attn\"] = (oof_nn, test_nn)\n\n# Ranks rather than raw probabilities: the candidates are on different scales\n# (ridge is an unbounded regression clipped to [0, 1], the others are sigmoids)\n# and the metric only reads order.\nRANKED = {k: (rank01(o), rank01(t)) for k, (o, t) in CANDIDATES.items()}\nY_TRUE = Y_soft.values\nnames = list(RANKED)\n\nsingles = {k: calc_macro_auc(Y_TRUE, RANKED[k][0]) for k in names}\nprint(\"\\n\" + \"=\" * 52, flush=True)\nprint(\"  CANDIDATE MODELS   self-consistency |  gold\", flush=True)\nprint(\"=\" * 52, flush=True)\nfor k in names:\n    _g, _ = gold_macro_auc(RANKED[k][0])\n    print(f\"  {k:<18} {singles[k]:.4f}       |  {_g:.4f}\", flush=True)\n\npicks = [max(names, key=lambda k: singles[k])]\nbest_score = singles[picks[0]]\n\nfor _ in range(ENSEMBLE_ROUNDS - 1):\n    contender, contender_score = None, best_score\n    for k in names:\n        trial = picks + [k]\n        blended = sum(RANKED[c][0] for c in trial) / len(trial)\n        s = calc_macro_auc(Y_TRUE, blended)\n        if s > contender_score + 1e-6:\n            contender, contender_score = k, s\n    if contender is None:\n        break                      # nothing left that helps; stop early\n    picks.append(contender)\n    best_score = contender_score\n\nWEIGHTS = {k: picks.count(k) / len(picks) for k in names if picks.count(k)}\noof_blend = sum(w * RANKED[k][0] for k, w in WEIGHTS.items())\ntest_blend = sum(w * RANKED[k][1] for k, w in WEIGHTS.items())\n\nauc_blend = calc_macro_auc(Y_TRUE, oof_blend)\nauc_blend_gold, n_gold_lab = gold_macro_auc(oof_blend)\n\nprint(\"=\" * 52, flush=True)\nprint(\"  SELECTED ENSEMBLE\", flush=True)\nprint(\"=\" * 52, flush=True)\nfor k, w in sorted(WEIGHTS.items(), key=lambda kv: -kv[1]):\n    print(f\"  {k:<18} weight {w:.3f}  ({picks.count(k)}/{len(picks)} rounds)\", flush=True)\n_dropped = [k for k in names if k not in WEIGHTS]\nif _dropped:\n    print(f\"  dropped (never improved the objective): {', '.join(_dropped)}\", flush=True)\nprint(f\"  blend self-consistency AUC : {auc_blend:.4f}\", flush=True)\n\n# The selector optimises self-consistency, so it can in principle pick a blend\n# that is worse on gold than the best single model. That would be selection noise\n# rather than a decision worth shipping, so it is reported explicitly instead of\n# being buried in the ensemble weights.\n_best_single_gold = max(\n    ((gold_macro_auc(RANKED[k][0])[0], k) for k in names),\n    key=lambda t: (-1.0 if np.isnan(t[0]) else t[0]))\n\nprint(\"=\" * 52, flush=True)\nprint(\"  GOLD-58 MACRO AUC (tracks the leaderboard metric)\", flush=True)\nprint(\"=\" * 52, flush=True)\nfor k in names:\n    _g, _ = gold_macro_auc(RANKED[k][0])\n    print(f\"  {k:<18} : {_g:.4f}\", flush=True)\nprint(f\"  {'SELECTED BLEND':<18} : {auc_blend_gold:.4f}\", flush=True)\nprint(f\"  best single on gold: {_best_single_gold[1]} at {_best_single_gold[0]:.4f}\", flush=True)\nif not np.isnan(auc_blend_gold) and auc_blend_gold < _best_single_gold[0] - 0.02:\n    print(\"  NOTE: the blend trails the best single model on gold by more than \"\n          \"0.02. Selection ran on the 4407-study objective; with n=58 this gap \"\n          \"is weak evidence, but it is worth a look before resubmitting.\", flush=True)\nprint(f\"  scored over {n_gold_lab} of {N_LABELS} labels on {len(GOLD_POS)} studies; \"\n      f\"labels: {LABEL_PROVENANCE}, folds: {GROUPING}\", flush=True)\nprint(\"  n=58 is small: differences under roughly 0.05 are noise.\", flush=True)\nprint(\"=\" * 52 + \"\\n\", flush=True)\n"},{"cell_type":"code","execution_count":null,"id":"6fed24a9","metadata":{},"outputs":[],"source":"# Cell 12: Final Submission Generation & Verification\nsub = pd.DataFrame({\"StudyInstanceUID\": test_ids})\nfor j, lab in enumerate(LABELS):\n    sub[lab] = test_blend[:, j]\n\n# Conform to sample_submission.csv when it is present. The row order and the\n# column order of the scored file come from the competition, not from the order\n# test.csv happens to be read in, so the predictions are reindexed onto it rather\n# than assumed to match. LABELS is only this notebook's internal ordering.\nif SAMPLE_SUB_CSV.is_file():\n    sample = pd.read_csv(SAMPLE_SUB_CSV)\n    sample[\"StudyInstanceUID\"] = sample[\"StudyInstanceUID\"].astype(str)\n    sample_labels = [c for c in sample.columns if c != \"StudyInstanceUID\"]\n\n    unknown = [c for c in sample_labels if c not in LABELS]\n    if unknown:\n        raise ValueError(\n            \"sample_submission.csv has columns this notebook never predicts: \"\n            f\"{unknown}. LABELS in Cell 2 must match the competition's column \"\n            \"names exactly, including capitalisation and apostrophes.\")\n\n    indexed = sub.set_index(\"StudyInstanceUID\")\n    missing = [u for u in sample[\"StudyInstanceUID\"] if u not in indexed.index]\n    if missing:\n        raise ValueError(\n            f\"{len(missing)} study IDs in sample_submission.csv were never \"\n            f\"predicted, first few: {missing[:3]}. The test manifest and the \"\n            \"sample submission disagree.\")\n\n    # The index name is set explicitly rather than left to reindex's inheritance\n    # rules, so reset_index() yields a StudyInstanceUID column on every pandas\n    # version rather than a stray \"index\" column.\n    sub = indexed.reindex(sample[\"StudyInstanceUID\"].tolist())[sample_labels]\n    sub.index.name = \"StudyInstanceUID\"\n    sub = sub.reset_index()\n    print(f\"✓ Conformed to {SAMPLE_SUB_CSV.name}: \"\n          f\"{len(sample_labels)} label columns, {len(sub)} rows.\", flush=True)\nelse:\n    print(\"Notice: sample_submission.csv not found; writing in test.csv order.\", flush=True)\n\nlabel_cols = [c for c in sub.columns if c != \"StudyInstanceUID\"]\n\n# Verification Assertions\nassert list(sub.columns)[0] == \"StudyInstanceUID\", \"Error: ID column is not first!\"\nassert sub[\"StudyInstanceUID\"].notna().all(), \"Error: missing StudyInstanceUID!\"\nassert not sub[\"StudyInstanceUID\"].duplicated().any(), \"Error: duplicate StudyInstanceUID!\"\nassert np.isfinite(sub[label_cols].values).all(), \"Error: Non-finite value in submission!\"\nassert (sub[label_cols].values >= 0.0).all() and (sub[label_cols].values <= 1.0).all(), \"Error: Probability out of range [0, 1]!\"\n\nWORK.mkdir(parents=True, exist_ok=True)\noutput_path = WORK / \"submission.csv\"\nsub.to_csv(output_path, index=False)\n\nprint(f\"✓ SUCCESS: Valid Kaggle submission written to {output_path}\", flush=True)\nprint(f\"Output shape: {sub.shape[0]} rows x {sub.shape[1]} columns\", flush=True)\nprint(\"\\nFirst 5 rows:\", flush=True)\nprint(sub.head().to_string(index=False), flush=True)\n\nprint(\"\\nPer-label mean predicted probability:\", flush=True)\nprint(sub[label_cols].mean().round(4).to_string(), flush=True)\n"}],"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"}},"nbformat":4,"nbformat_minor":5}