{"cells":[{"cell_type":"markdown","metadata":{},"source":"# RSNA Knee 101 — read the data before you train\n\n**A beginner's tour of this competition: what you are predicting, what is actually in the\nfiles, and what DINOv2 is doing when you use it.**\n\nIf you have never opened a DICOM, never used a vision transformer, and are not sure why\neveryone keeps saying \"the fat suppression flag is broken\" — this notebook is for you. It\nassumes you know Python and pandas, and nothing else.\n\nEverything here is **computed in front of you**. Where a number appears in the text it is\nbecause a cell above it printed that number on this run. Nothing is quoted from memory,\nincluding the findings this notebook credits to other people — those are re-measured too,\nbecause a finding you have re-measured is a finding you can build on.\n\n### How to read it\n\n| Chapter | Question it answers |\n|---|---|\n| 1 | What am I predicting, and what does the score reward? |\n| 2 | **The twelve findings in plain language** — what they are, where they live, how big |\n| 3 | Where are the labels? (Almost nowhere. This is the whole problem.) |\n| 4 | The radiology reports — your real training signal |\n| 5 | The series table, and the one column that is lying to you |\n| 6 | What is inside a DICOM file? |\n| 7 | Millimetres, not pixels: the geometry that decides your preprocessing |\n| 8 | **DINOv2 101** — what it is, what a patch token is, how to load and run it |\n| 9 | A checklist of the traps, so you do not pay for them twice |\n| 10 | **From findings to code** — the four fixes, running on the real data |\n\nRuns on **CPU or GPU** in roughly ten minutes. It reads every CSV in full and a\nconfigurable sample of the DICOM headers, so it costs you almost nothing to fork and\nre-run with the sample turned up."},{"cell_type":"markdown","metadata":{},"source":"## Credit where it is due\n\nThis notebook is a synthesis and an explainer. It did not discover the things it\nexplains — the public notebooks below did, and several of them found the traps the hard\nway, in public, before the rest of us walked into them. Each row names what that notebook\nestablished first; the chapter that uses it says so again in place.\n\n| Notebook | Author | What it contributed |\n|---|---|---|\n| [RSNA Knee baseline v1](https://www.kaggle.com/code/pilkwang/rsna-knee-baseline-v1) | **Pilkwang Kim** | The reference baseline for this competition. Recovering the acquisition axes from DICOM headers instead of trusting the CSV flags, sampling every series at a constant physical scale, normalising left and right knees onto one convention, per-diagnosis attention over series slots, and grouping folds on the report hash |\n| [RSNA Knee: EDA to 2.5D](https://www.kaggle.com/code/karnakbaevarthur/rsna-knee-eda-to-2-5d) | **Karnakbayev Artur** | The audit that named the two facts everything else has to design around: only a small subset of training studies carry per-condition labels, and `Fluid_Sensitive` is identical to `Fat_Suppression` on every row. Also the first protocol-only baseline |\n| [RSNA Knee DINOv2 at meniscus resolution](https://www.kaggle.com/code/wguesdon/rsna-knee-dinov2-at-meniscus-resolution) | **Will** (`wguesdon`) | The resolution argument: convert input size into *millimetres per DINOv2 patch token* and compare that against the size of the lesion. Also ordering slices by physical position rather than by filename |\n| [RSNA Knee EDA: duplicate reports, broken flag](https://www.kaggle.com/code/alexandremoritz/rsna-knee-eda-duplicate-reports-broken-flag) | **Alexandre Moritz** | Duplicate report text across studies, and the broken flag, called out early |\n| [RSNA Knee EDA: the reports are train-only](https://www.kaggle.com/code/wguesdon/rsna-knee-eda-the-reports-are-train-only) | **Will** (`wguesdon`) | Reports exist in `train.csv` and not in `test.csv` — which is what rules out an entire class of model |\n| [RSNA Knee \\| Data structure, EDA, baseline](https://www.kaggle.com/code/romanrozen/rsna-knee-data-structure-eda-baseline) | **Roman Rozen** | A widely-read walkthrough of the same pipeline family, with weight EMA and a pairwise ranking loss on top |\n\nChapter 8 also uses **Meta AI's DINOv2**, mounted on Kaggle as\n`metaresearch/dinov2/PyTorch/small/1`\n([paper](https://arxiv.org/abs/2304.07193), Oquab et al., 2023 — *DINOv2: Learning Robust\nVisual Features without Supervision*).\n\nIf you fork this, please keep this block. It costs you nothing and it is how a competition\nforum stays worth reading."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 1 — What am I predicting?\n\nOne **study** is one visit to the MRI scanner. In that visit the radiographer runs several\n**series** — the same knee imaged several times with different settings and from different\nangles. Each series is a stack of **slices**, and each slice is one `.dcm` file.\n\n```\nstudy (one patient visit) -> 4-8 series -> 20-40 slices each -> a .dcm file per slice\n```\n\nFor each *study* you predict **twelve probabilities**, one per finding:\n\n| | Finding | In plain English |\n|---|---|---|\n| 1 | ACL | Anterior cruciate ligament injury — the classic \"I heard a pop\" ligament |\n| 2 | MCL | Medial collateral ligament injury — the ligament on the inner side |\n| 3 | Medial Meniscus | Tear in the inner cartilage cushion |\n| 4 | Lateral Meniscus | Tear in the outer cartilage cushion |\n| 5 | Medial OA | Arthritis in the inner compartment |\n| 6 | Lateral OA | Arthritis in the outer compartment |\n| 7 | PF OA | Arthritis behind the kneecap (patellofemoral) |\n| 8 | Effusion | Fluid in the joint — a swollen knee |\n| 9 | Synovitis | Inflamed joint lining |\n| 10 | Baker's | A fluid-filled cyst behind the knee |\n| 11 | Contusion | Bone bruise |\n| 12 | Fracture | Broken bone |\n\nFour of them — the two menisci and the two side-specific arthritis labels — are **left/right\nwithin the knee** (medial = toward the other leg, lateral = away). Remember that; Chapter 7\ncomes back to it and it is the single easiest way to silently destroy four of your twelve\nscores."},{"cell_type":"markdown","metadata":{},"source":"### The score, and the two things it tells you to do\n\n$$\\text{Score} \\;=\\; \\frac{1}{12}\\sum_{i=1}^{12} \\mathrm{AUC}_i$$\n\n**AUC** (area under the ROC curve) asks one question: *pick a random positive study and a\nrandom negative study — how often does your model give the positive one the higher score?*\n1.0 is perfect, 0.5 is a coin flip.\n\nTwo consequences, and both of them delete work rather than add it:\n\n**1. Only the *order* of your predictions matters.** AUC does not look at the values, only\nat the ranking. Multiplying a whole column by 3, or squaring it, changes nothing. So:\n\n- Calibration is worth zero. Do not spend time on it.\n- Picking a threshold is worth zero. There is no threshold.\n- When you combine two models, **average their ranks, not their probabilities.** Averaging\n probabilities lets whichever model outputs more extreme numbers dominate the blend, and\n \"more extreme\" is not the same as \"more correct\".\n\n**2. Every label is worth exactly the same.** Effusion is common and fracture is rare, but\neach contributes $1/12$ of the score. A label you leave at chance costs you about\n$(0.85 - 0.5)/12 \\approx 0.029$ — no matter how good the other eleven are. **Rare findings\ndeserve more of your attention, not less**, because a rare finding is where a model most\neasily collapses to chance."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Setup. Nothing here is competition-specific except the label list.\n# ---------------------------------------------------------------------------\nimport hashlib\nimport os\nimport re\nimport time\nfrom collections import Counter\nfrom concurrent.futures import ThreadPoolExecutor\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport matplotlib.pyplot as plt\n\npd.set_option(\"display.width\", 140)\npd.set_option(\"display.max_columns\", 40)\n\nLABELS = [\n \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n \"Medial OA\", \"Lateral OA\", \"PF OA\", \"Effusion\",\n \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\",\n]\n\n# How much of the DICOM tree to open. Every CSV is read in full; only the header\n# sweep in Chapter 7 is sampled, because opening 24,000 directories is slow and\n# the answers it gives converge long before you have opened all of them.\n# Turn this up if you fork and want the numbers on the whole corpus.\nN_HEADER_STUDIES = 400\nSEED = 101\n\nrng = np.random.default_rng(SEED)"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Find the data.\n#\n# Kaggle does not always mount a competition at /kaggle/input/. This one is\n# one level deeper, under /kaggle/input/competitions/, which a shallow scan\n# misses. So try the known shapes first, then walk two levels as a fallback.\n# ---------------------------------------------------------------------------\ndef find_data_root():\n \"\"\"Return the directory that holds train.csv, or raise.\"\"\"\n known = [\n Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\"),\n Path(\"/kaggle/input/rsna-knee-abnormality-detection\"),\n Path(\"data\"),\n Path(\".\"),\n ]\n for candidate in known:\n if (candidate / \"train.csv\").is_file():\n return candidate\n\n base = Path(\"/kaggle/input\")\n if base.is_dir():\n for level1 in sorted(p for p in base.iterdir() if p.is_dir()):\n for candidate in [level1] + sorted(p for p in level1.iterdir() if p.is_dir()):\n if (candidate / \"train.csv\").is_file():\n return candidate\n raise FileNotFoundError(\"could not find train.csv - is the competition attached?\")\n\n\nROOT = find_data_root()\nprint(\"data root:\", ROOT)\nprint()\n\nfor item in sorted(ROOT.iterdir()):\n if item.is_dir():\n n = sum(1 for _ in item.iterdir())\n print(f\" {item.name + '/':22s} {n:>6,} entries\")\n else:\n print(f\" {item.name:22s} {item.stat().st_size / 1e6:>6.1f} MB\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Load the four tables and look at their shape before anything else.\n# ---------------------------------------------------------------------------\ntrain = pd.read_csv(ROOT / \"train.csv\")\ntest = pd.read_csv(ROOT / \"test.csv\")\ntrain_series = pd.read_csv(ROOT / \"train_series.csv\")\ntest_series = pd.read_csv(ROOT / \"test_series.csv\")\n\nfor name, frame in [(\"train\", train), (\"test\", test),\n (\"train_series\", train_series), (\"test_series\", test_series)]:\n print(f\"{name + '.csv':18s} {frame.shape[0]:>7,} rows x {frame.shape[1]:>2} cols\")\n print(f\"{'':18s} {list(frame.columns)}\")\n print()"},{"cell_type":"markdown","metadata":{},"source":"### Read that output again — two things should stop you\n\n**`test.csv` is tiny.** It holds a placeholder while the competition is running; the real\ntest set is swapped in when your notebook is re-run for scoring. Practical consequence:\n**you cannot compare train and test distributions**, and any \"train/test shift\" analysis\nyou do on the visible test set is measuring three rows. Your only honest signal is\ncross-validation on train.\n\n**`train.csv` has a `Report` column and `test.csv` does not** — first flagged publicly by\n[Will's EDA](https://www.kaggle.com/code/wguesdon/rsna-knee-eda-the-reports-are-train-only).\nThis is a bigger deal than it looks. Text is available when you *fit* and absent when you\n*predict*. So a model with a text branch is impossible. Reports can only be used to\nmanufacture training targets — never as an input feature."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 2 — The twelve findings, in plain language\n\nEvery other chapter treats the twelve labels as twelve columns. They are not. They are\ntwelve **different things, of different sizes, in different places in the knee**, and that\nturns out to decide preprocessing choices you would otherwise make by taste.\n\nIf you are not from a medical background, this chapter is the one to read slowly. It costs\nten minutes and it makes Chapter 7 (millimetres) and Chapter 8 (DINOv2) make sense.\n\n**A note on what kind of statement this is.** The anatomy table below is *background* — it\ncomes from radiology, not from measuring this dataset. Everything after it is measured\nfrom the competition data. The notebook keeps those two kinds of claim visibly separate,\nand so should you.\n\n## The knee in four sentences\n\nThe knee is where the **femur** (thigh bone) meets the **tibia** (shin bone), with the\n**patella** (kneecap) riding in a groove at the front.\n\nBetween femur and tibia sit two crescent-shaped shock absorbers, the **menisci** — one on\nthe inner side (**medial**) and one on the outer (**lateral**).\n\nFour ligaments hold it together; two of them are targets here — the **ACL** (anterior\ncruciate ligament, deep in the middle) and the **MCL** (medial collateral, on the inner\nedge).\n\nEverything is wrapped in a capsule lined with **synovium**, which makes lubricating fluid —\nand which, when it misbehaves, produces most of the remaining targets.\n\n## The twelve, one line each\n\n| # | Target | What it is | Where it lives | Rough size |\n|---|---|---|---|---|\n| 1 | **ACL** | tear/sprain of the central stabilising ligament | intercondylar notch, dead centre | 7–11 mm thick |\n| 2 | **MCL** | tear/sprain of the inner-side ligament | medial edge of the joint | 2–4 mm thick |\n| 3 | **Medial Meniscus** | tear or degeneration of the inner shock absorber | inner joint line | 9–12 mm wide, tear 1–3 mm |\n| 4 | **Lateral Meniscus** | same, outer side | outer joint line | 9–12 mm wide, tear 1–3 mm |\n| 5 | **Medial OA** | cartilage loss / arthritis, inner compartment | inner femur–tibia contact | cartilage 2–4 mm |\n| 6 | **Lateral OA** | same, outer compartment | outer femur–tibia contact | cartilage 2–4 mm |\n| 7 | **PF OA** | arthritis behind the kneecap | patella–femur groove, front | cartilage 2–4 mm |\n| 8 | **Effusion** | excess fluid in the joint | pools **above** the kneecap | 10–50 mm collection |\n| 9 | **Synovitis** | inflamed, thickened joint lining | throughout the capsule | diffuse, ill-defined |\n| 10 | **Baker's** | fluid-filled cyst | **behind** the knee, popliteal fossa | 20–50 mm |\n| 11 | **Contusion** | bone bruise — marrow swelling, no break | inside the bone ends | 10–30 mm patch |\n| 12 | **Fracture** | actual break in the bone | anywhere in femur/tibia/patella | line, <1 mm wide |\n\nTwo columns in that table are doing quiet work. **Size** tells you what resolution you\nneed. **Where** tells you what happens if you crop too tightly — note that Effusion lives\n*above* the joint and Baker's lives *behind* it, so both sit near the edge of the field of\nview rather than in the middle."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Real sentences, pulled live from the reports, for each of the twelve.\n#\n# Reading how radiologists actually phrase a finding is the fastest way to\n# understand it - and it is also exactly what your report labeller has to cope\n# with. Note the hedging (\"R/O\" = rule out), the grading, and the synonyms.\n#\n# Two filters matter here and are worth copying:\n# - a NEGATION filter, because \"no joint effusion\" contains \"effusion\"\n# - a non-Latin filter, only so the printed examples stay readable; the\n# corpus itself is multilingual and Chapter 4 measures that.\n# ---------------------------------------------------------------------------\nimport re\n\nDEGEN = r\"chondr|cartilage|osteophyt|joint space|arthros|arthrit|degenerat\"\nFINDING_RX = {\n \"ACL\": (r\"anterior cruciate\", r\"tear|rupture|disrupt|sprain\"),\n \"MCL\": (r\"medial collateral\", r\"tear|sprain|thicken|edema|oedema\"),\n \"Medial Meniscus\": (r\"medial meniscus\", r\"tear|degenerat|macerat|fray\"),\n \"Lateral Meniscus\": (r\"lateral meniscus\", r\"tear|degenerat|macerat|fray\"),\n \"Medial OA\": (r\"medial (?:compartment|femorotibial|femoral condyle|\"\n r\"tibial plateau)\", DEGEN),\n \"Lateral OA\": (r\"lateral (?:compartment|femorotibial|femoral condyle|\"\n r\"tibial plateau)\", DEGEN),\n \"PF OA\": (r\"patellofemoral|retropatellar|patellar facet|chondromalacia\",\n DEGEN),\n \"Effusion\": (r\"effusion|joint fluid\",\n r\"mild|moderate|large|small|trace|present|distend\"),\n \"Synovitis\": (r\"synovitis|synovial\",\n r\"thicken|hypertroph|proliferat|synovitis|enhanc\"),\n \"Baker's\": (r\"baker|popliteal cyst\", r\"cyst\"),\n \"Contusion\": (r\"bone (?:marrow )?(?:edema|oedema|bruise|bruising|contusion)\"\n r\"|contusion\", r\".\"),\n \"Fracture\": (r\"fracture\", r\".\"),\n}\n# \"preserved\" and \"maintained\" are the ones that catch people out: they are negations\n# containing no negative word at all. A radiologist writing \"the cartilage is preserved\"\n# is saying the finding is ABSENT, and a labeller that only looks for \"no\" will score\n# that sentence as a confident positive.\nNEGATION = re.compile(r\"\\b(no|not|non|without|absent|intact|normal|unremarkable|\"\n r\"negative|preserved|maintained|spared)\\b\", re.I)\nNON_LATIN = re.compile(r\"[ığşçöüİĞŞÇÖÜáéíóúñàèùâêôëïüßäöÄÖÜĆČŽŠĐ]\")\n\n\ndef sentences(text):\n return [s.strip(\" -•>\\t\")\n for s in re.split(r\"(?<=[.;])\\s+|\\n\", str(text))\n if 30 < len(s.strip()) < 175]\n\n\nreports = train[\"Report\"].fillna(\"\").tolist()\nfor target, (anat_rx, qual_rx) in FINDING_RX.items():\n anat, qual = re.compile(anat_rx, re.I), re.compile(qual_rx, re.I)\n found, seen = [], set()\n for rep in reports:\n for s in sentences(rep):\n if (anat.search(s) and qual.search(s) and not NEGATION.search(s)\n and not NON_LATIN.search(s) and s.count(\" \") > 5\n and not s.lower().startswith(\"in the\")):\n key = re.sub(r\"\\W\", \"\", s.lower())[:40]\n if key not in seen:\n seen.add(key)\n found.append(s)\n if len(found) >= 2:\n break\n print(f\"{target}\")\n for s in found[:2]:\n print(f\" {s}\")\n print()"},{"cell_type":"markdown","metadata":{},"source":"Read those and a few things jump out that matter for the labeller in Chapter 4:\n\n- **Hedging is everywhere.** `R/O tear` means *rule out* — the radiologist is raising a\n possibility, not asserting one. A naive keyword match scores that as a confident\n positive.\n- **Grading is common.** \"grade I sprain\", \"low grade chondrosis\", \"trace Baker's cyst\".\n The finding is present but minimal. Is a trace effusion a positive? The gold annotator\n and your regex may disagree, and that disagreement becomes label noise.\n- **Negation is structural, not just the word \"no\".** \"intact\", \"normal\", \"unremarkable\",\n \"within normal limits\" all mean absent.\n- **Findings co-occur in one sentence.** \"posterior horn of medial meniscus and anterior\n horn of lateral meniscus\" is two different labels in one clause, and a sloppy parser will\n assign the wrong laterality to one of them.\n\n## Why the size column decides your preprocessing\n\nChapter 8 explains that DINOv2 does not see pixels — it cuts the image into **patches** of\n14×14 pixels, and each patch becomes one token. So the question that governs everything is:\n**how many millimetres of knee land inside one token?**"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Millimetres per patch token, against the size of each finding.\n#\n# The only thing that sets mm-per-pixel is (physical extent cropped) / (pixels).\n# Two knobs, one ratio - so a tighter crop and a bigger input buy the SAME\n# resolution, but only one of them costs compute. Chapter 7 measures the field\n# of view this corpus actually has.\n# ---------------------------------------------------------------------------\nPATCH = 14 # DINOv2 patch size, in pixels\n\nFINDING_MM = {\n \"Fracture line\": 0.5, \"Meniscal tear\": 2.0, \"Cartilage (OA)\": 3.0,\n \"MCL thickness\": 3.0, \"ACL thickness\": 9.0, \"Meniscus width\": 10.0,\n \"Bone contusion\": 20.0, \"Baker's cyst\": 35.0, \"Effusion\": 30.0,\n}\n\nprint(f\"{'crop':>6} {'input':>6} {'mm/pixel':>9} {'mm/token':>9} \"\n f\"{'compute':>8}\")\nfor crop_mm, img in [(160, 224), (160, 336), (160, 448), (107, 224)]:\n n_tok = (img // PATCH) ** 2 + 1\n base = (224 // PATCH) ** 2 + 1\n # ViT cost: linear in tokens for the MLP, quadratic for attention.\n d, mlp = 384, 4\n cost = lambda n: (6 + 2 + 4 * mlp) * n * d ** 2 + 4 * n ** 2 * d\n print(f\"{crop_mm:>5}mm {img:>5}px {crop_mm / img:>9.3f} \"\n f\"{crop_mm / img * PATCH:>9.2f} {cost(n_tok) / cost(base):>7.1f}x\")\n\nprint(\"\\nHow wide is each finding, in TOKENS, at our default 160mm / 224px?\")\nmm_per_token = 160 / 224 * PATCH\nfor name, mm in sorted(FINDING_MM.items(), key=lambda kv: kv[1]):\n tokens = mm / mm_per_token\n bar = \"#\" * max(1, int(tokens * 8))\n flag = \" <- smaller than one token\" if tokens < 1 else \"\"\n print(f\" {name:<16} {mm:>5.1f} mm = {tokens:>5.2f} tokens {bar}{flag}\")"},{"cell_type":"markdown","metadata":{},"source":"That printout is the single most useful thing in this notebook for deciding preprocessing.\n\n**Seven of the nine structures are smaller than one token.** A meniscal tear is about\n**one fifth** of a token. That does not make it undetectable — the patch embedding is a\nlearned linear projection of 588 numbers into 384 dimensions, so very little is thrown\naway. But it does mean the tear cannot be expressed as a *relationship between two tokens*,\nand comparing tokens is the only thing self-attention does. At 224 px the whole meniscus is\none token, so \"bright streak through dark tissue\" has to survive inside a single vector.\n\nMeanwhile **Effusion, Baker's cyst and bone contusion are 2–5 tokens wide** and will\nsurvive almost any sensible resolution.\n\nNotice also the last row of the first table: **107 mm at 224 px gives the same mm/token as\n160 mm at 336 px, at one fortieth of the compute** — because only the *ratio* matters. The\nprice is that you have thrown away 26 mm from every edge, which is exactly where Effusion\n(above the kneecap) and Baker's cyst (behind the knee) live.\n\n**That is a real trade, not a free win, and it is cheap to test.** Under a metric that is\nthe unweighted mean of twelve AUCs, sharpening the menisci at the cost of blinding yourself\nto two other findings can easily be a net loss. Measure it; do not reason about it."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 3 — Where are the labels?\n\nThis is the defining problem of the competition, and it is worth understanding before you\nwrite a single line of model code.\n\nCredit: this was first laid out in\n[Karnakbayev Artur's EDA](https://www.kaggle.com/code/karnakbaevarthur/rsna-knee-eda-to-2-5d).\nThe cell below re-measures it."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# How many training studies actually carry the twelve labels?\n#\n# A study is \"labelled\" (\"gold\") if its twelve label columns are not null.\n# ---------------------------------------------------------------------------\nhas_labels = train[LABELS].notna().all(axis=1)\ngold = train[has_labels]\n\nprint(f\"training studies {len(train):>6,}\")\nprint(f\" with all 12 labels {len(gold):>6,} <- your only ground truth\")\nprint(f\" with no labels {(~has_labels).sum():>6,}\")\nprint(f\" labelled fraction {len(gold) / len(train):>6.1%}\")\nprint()\n\nn_report = train[\"Report\"].notna().sum() if \"Report\" in train.columns else 0\nprint(f\" with a radiology report {n_report:>6,} ({n_report / len(train):.1%})\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# What do the labelled studies contain? Two facts per label: how often it is\n# positive, and how many positive examples that actually is.\n# ---------------------------------------------------------------------------\nprevalence = pd.DataFrame({\n \"positives\": gold[LABELS].sum().astype(int),\n \"rate\": gold[LABELS].mean(),\n}).sort_values(\"rate\", ascending=False)\n\nprint(prevalence.to_string(formatters={\"rate\": \"{:.1%}\".format}))\nprint()\nprint(f\"findings per labelled study: mean {gold[LABELS].sum(axis=1).mean():.2f}, \"\n f\"min {int(gold[LABELS].sum(axis=1).min())}, \"\n f\"max {int(gold[LABELS].sum(axis=1).max())}\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Chart: positive rate per finding among the labelled studies.\n#\n# A horizontal bar chart, because the category names are long and the job of the\n# chart is magnitude comparison. Values are printed on the bars, so the reader\n# never has to measure a bar against a gridline.\n# ---------------------------------------------------------------------------\nSURFACE, INK, MUTED, GRID = \"#fcfcfb\", \"#0b0b0b\", \"#6a6862\", \"#e1e0d9\"\nBLUE, ORANGE, AQUA = \"#2a78d6\", \"#eb6834\", \"#1baf7a\"\n\nplt.rcParams.update({\n \"figure.facecolor\": SURFACE, \"axes.facecolor\": SURFACE,\n \"text.color\": INK, \"axes.labelcolor\": MUTED, \"axes.edgecolor\": GRID,\n \"xtick.color\": MUTED, \"ytick.color\": INK, \"font.size\": 11,\n \"axes.spines.top\": False, \"axes.spines.right\": False,\n})\n\norder = prevalence.sort_values(\"rate\")\nfig, ax = plt.subplots(figsize=(9, 5.6))\nax.barh(order.index, order[\"rate\"], height=0.66, color=BLUE, zorder=3)\n\nfor name, row in order.iterrows():\n ax.text(row[\"rate\"] + 0.012, name, f\"{row['rate']:.0%} (n={int(row['positives'])})\",\n va=\"center\", fontsize=10, color=MUTED)\n\nax.set_xlim(0, max(0.05, order[\"rate\"].max() * 1.42))\nax.xaxis.set_major_formatter(lambda v, _: f\"{v:.0%}\")\nax.grid(axis=\"x\", color=GRID, lw=0.8, zorder=0)\nax.set_xlabel(\"share of labelled studies that are positive\")\nax.set_title(f\"Positive rate per finding, across the {len(gold)} labelled studies\",\n loc=\"left\", fontsize=13, fontweight=\"bold\", pad=14)\nfig.text(0.005, -0.02,\n \"Every labelled study has at least one positive finding, so these rates are \"\n \"higher than the rates in the full corpus.\",\n fontsize=9.5, color=MUTED)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","metadata":{},"source":"### What this means for your model\n\nA handful of labelled studies out of thousands is not a training set. It is barely a\nvalidation set. Three things follow:\n\n**1. You cannot train on the labels.** With this many positives for a rare finding, a model\nfitted directly on them memorises them. Everyone competitive here trains on labels\n*derived from the reports* (Chapter 4) and keeps the real labels for checking.\n\n**2. The labelled studies are not a random sample.** Look at the \"findings per labelled\nstudy\" line — every labelled study has at least one positive finding. So the positive rates\nin the chart are **higher than the true rates in the corpus**, and you cannot read them as\nprevalence. They were selected for annotation *because* something was there.\n\n**3. Your validation number is noisy, and you should say so.** The standard error on an AUC\ncomputed from this many studies is roughly ±0.07. A change of +0.01 on this reference is\nnot a change. If you tune on it, you are tuning on noise — and you will find out on the\nprivate leaderboard."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# How noisy? A quick bootstrap, so the number above is not just an assertion.\n#\n# Take a label, generate random predictions, and see how far AUC wanders purely by\n# resampling which studies you happened to score. This is the floor on how precisely\n# you can measure anything with this reference.\n# ---------------------------------------------------------------------------\nfrom sklearn.metrics import roc_auc_score\n\n\ndef bootstrap_auc_spread(y_true, n_boot=2000, seed=SEED):\n \"\"\"Spread of AUC for a *random* scorer, resampling the studies with replacement.\"\"\"\n y_true = np.asarray(y_true)\n scores = np.random.default_rng(seed).random(len(y_true))\n aucs = []\n boot_rng = np.random.default_rng(seed + 1)\n for _ in range(n_boot):\n idx = boot_rng.integers(0, len(y_true), len(y_true))\n if len(np.unique(y_true[idx])) < 2: # a resample with only one class has no AUC\n continue\n aucs.append(roc_auc_score(y_true[idx], scores[idx]))\n return np.std(aucs), np.percentile(aucs, [2.5, 97.5])\n\n\nprint(f\"{'label':<18} {'n_pos':>5} {'std':>6} 95% interval for a *random* model\")\nprint(\"-\" * 66)\nfor label in LABELS:\n y = gold[label].astype(int).values\n if y.sum() < 2 or y.sum() > len(y) - 2:\n print(f\"{label:<18} {int(y.sum()):>5} {'--':>6} too few of one class to measure\")\n continue\n std, (lo, hi) = bootstrap_auc_spread(y)\n print(f\"{label:<18} {int(y.sum()):>5} {std:>6.3f} [{lo:.2f}, {hi:.2f}]\")\n\nprint()\nprint(\"Read the last column as: a model that knows nothing can land anywhere in there.\")\nprint(\"Any improvement smaller than that interval is not measurable on this reference.\")"},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 4 — The reports are the real training signal\n\nNearly every training study ships the original free-text radiology report. That text is\nwhere your labels have to come from.\n\nTwo properties of this text shape everything you build on top of it."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Look at an actual report before theorising about them.\n# ---------------------------------------------------------------------------\nreports = train[\"Report\"].fillna(\"\").astype(str)\n\nlengths = reports.str.len()\nprint(\"length in characters:\")\nprint(lengths.describe().to_string())\nprint()\n\n# Show a report of typical length rather than the first one, which may be an outlier.\ntypical = (lengths - lengths.median()).abs().idxmin()\nprint(\"=\" * 78)\nprint(f\"ONE REPORT, VERBATIM (study index {typical}, {lengths[typical]} characters)\")\nprint(\"=\" * 78)\nprint(reports[typical][:1200])"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Property 1: the reports are multilingual.\n#\n# A crude detector is enough to make the point: look at the alphabet, then at a few\n# unmistakable function words. This is not a language ID model and does not need to be;\n# it needs to show you that the corpus is not English.\n# ---------------------------------------------------------------------------\nSCRIPT_RANGES = [(\"Greek\", 0x0370, 0x03FF), (\"Cyrillic\", 0x0400, 0x04FF)]\n\n# English is scored in the SAME competition as the others, not checked afterwards as a\n# fallback. The first version of this cell tested the other languages first and fell\n# through to English, and an English report saying \"No joint effusion. No cartilage\n# defects.\" matched the Spanish list twice on the word \"no\" and was counted as Spanish.\n# It reported 445 English reports out of 4,407. Shared function words have been dropped\n# from the lists for the same reason - only distinctive markers are left.\nLANGUAGE_MARKERS = [\n (\"English\", r\"\\b(the|and|is|are|with|there|no|of|intact|tear|joint)\\b\"),\n (\"German\", r\"\\b(und|nicht|kein|keine|des|der|die|mit|nachweis|gelenk)\\b\"),\n (\"French\", r\"\\b(aucun|aucune|avec|sans|les|une|articulaire|epanchement)\\b\"),\n (\"Spanish\", r\"\\b(del|los|las|una|sin|derrame|rotura|articular|se observa)\\b\"),\n (\"Portuguese\", r\"\\b(nao|sem|dos|das|joelho|derrame|observa-se)\\b\"),\n (\"Dutch\", r\"\\b(geen|niet|het|een|van de|kruisband|gewricht)\\b\"),\n (\"Turkish\", r\"\\b(ve|yok|izlenmektedir|izlenmemektedir|mevcut|eklem)\\b\"),\n (\"Italian\", r\"\\b(non|senza|della|ginocchio|si osserva|versamento)\\b\"),\n]\n\n\ndef guess_language(text):\n \"\"\"Return a coarse language guess for one report. Indicative, not exact.\n\n Alphabet first - Greek and Cyrillic are unambiguous. Otherwise score every\n language including English and take the highest, so no language wins by default.\n \"\"\"\n low = text.lower()\n for name, lo, hi in SCRIPT_RANGES:\n if sum(lo <= ord(ch) <= hi for ch in text[:600]) > 12:\n return name\n hits = [(len(re.findall(pattern, low)), name) for name, pattern in LANGUAGE_MARKERS]\n best_count, best_name = max(hits)\n return best_name if best_count >= 2 else \"unclear\"\n\n\nlanguages = reports.map(guess_language)\ncounts = languages.value_counts()\nprint(counts.to_string())\nprint()\nprint(f\"share not detected as English: \"\n f\"{1 - counts.get('English', 0) / counts.sum():.1%}\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Chart: rough language mix. One series, so no legend - the title names it.\n# ---------------------------------------------------------------------------\ntop = counts.head(9).sort_values()\nfig, ax = plt.subplots(figsize=(9, 4.8))\ncolors = [ORANGE if name in (\"unclear\", \"English\") else BLUE for name in top.index]\nax.barh(top.index, top.values, height=0.66, color=colors, zorder=3)\nfor name, value in top.items():\n ax.text(value + top.max() * 0.012, name, f\"{value:,} ({value / counts.sum():.0%})\",\n va=\"center\", fontsize=10, color=MUTED)\nax.set_xlim(0, top.max() * 1.30)\nax.grid(axis=\"x\", color=GRID, lw=0.8, zorder=0)\nax.set_xlabel(\"reports\")\nax.set_title(\"Rough language mix of the radiology reports (heuristic detector)\",\n loc=\"left\", fontsize=13, fontweight=\"bold\", pad=14)\nfig.text(0.005, -0.03,\n \"Orange = English and the residual bucket. The exact split depends on the \"\n \"detector; the point is that a large share of the corpus is not English.\",\n fontsize=9.5, color=MUTED)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","metadata":{},"source":"### Why the language mix is a modelling problem, not a curiosity\n\nThe obvious plan is a keyword rule: search the report for `tear`, `rupture`, `effusion`,\nand emit a 1 when you find one.\n\nHere is the trap. **A rule that never fires does not raise an error — it emits a negative.**\nSo a keyword list that is thick in English and thin in Greek does not look broken. It looks\nlike *a corpus in which Greek patients have fewer injuries*.\n\nAnd that is not random noise, which would average out. Language tracks the hospital, which\ntracks the scanner, the referral pattern and the patient population. A gap in your lexicon\ntherefore becomes a bias that is **aligned with a site** — exactly the kind of bias that\nsurvives cross-validation and shows up on the private leaderboard.\n\nIf you write a rule extractor, measure two things per label, not one:\n\n- **the fire rate** — on what share of reports does *any* rule for this label match? A label\n that is silent on most reports is being trained almost entirely on your default value.\n- **agreement with the gold labels** — on the labelled studies from Chapter 3.\n\nA rule can score well on the second and still be useless because of the first."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Property 2: some reports are byte-identical across different studies.\n#\n# Credit: called out early by Alexandre Moritz's EDA, and used as the fold-grouping\n# key in Pilkwang Kim's baseline v1.\n#\n# Why it matters: if the same report appears in your training fold and your validation\n# fold, then any model that learned from the report text has effectively seen the\n# validation answer. Your CV goes up and the leaderboard does not move.\n# ---------------------------------------------------------------------------\nnormalised = reports.str.lower().str.replace(r\"\\s+\", \" \", regex=True).str.strip()\nnonempty = normalised[normalised.str.len() > 0]\n\ngroup_sizes = nonempty.value_counts()\nduplicated_text = group_sizes[group_sizes > 1]\n\nprint(f\"non-empty reports {len(nonempty):>6,}\")\nprint(f\"distinct report texts {len(group_sizes):>6,}\")\nprint(f\"texts appearing in more than 1 study {len(duplicated_text):>6,}\")\nprint(f\"studies sharing text with another {int(duplicated_text.sum()):>6,} \"\n f\"({duplicated_text.sum() / len(nonempty):.1%})\")\nprint()\nprint(\"how many studies share one report text:\")\nprint(duplicated_text.value_counts().sort_index().rename(\"report texts\")\n .rename_axis(\"studies sharing it\").to_string())\nprint()\nprint(\"Fix: hash the normalised report and pass it as `groups=` to GroupKFold, so every\")\nprint(\"study sharing a report lands in the same fold.\")"},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 5 — The series table, and the column that is lying to you\n\n`train_series.csv` describes each series with three things: which way the scanner sliced\n(`Anatomical_Plane`) and two flags about how the image was acquired.\n\nThose two flags name **physically independent** properties:\n\n- **`Fluid_Sensitive`** — is this sequence one where fluid appears bright? That is set by\n the pulse sequence timing (T2 and proton-density weighting: fluid bright; T1 weighting:\n fluid dark). It is how you see an effusion, a tear filled with fluid, a bone bruise.\n- **`Fat_Suppression`** — was a fat-cancelling preparation applied on top? Fat is bright and\n everywhere in a knee, so suppressing it makes an oedema next to it visible.\n\nYou can have either, both, or neither. Two independent flags means four combinations.\n\nCredit for finding what follows:\n[Karnakbayev Artur](https://www.kaggle.com/code/karnakbaevarthur/rsna-knee-eda-to-2-5d)\nand [Alexandre Moritz](https://www.kaggle.com/code/alexandremoritz/rsna-knee-eda-duplicate-reports-broken-flag)."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Are Fluid_Sensitive and Fat_Suppression actually two different columns?\n# ---------------------------------------------------------------------------\nboth_series = pd.concat([train_series.assign(split=\"train\"),\n test_series.assign(split=\"test\")], ignore_index=True)\n\nfs = both_series[\"Fluid_Sensitive\"]\nfat = both_series[\"Fat_Suppression\"]\n\nidentical = int((fs.fillna(-1) == fat.fillna(-1)).sum())\nprint(f\"series rows compared {len(both_series):>7,}\")\nprint(f\"rows where the two flags agree {identical:>7,} ({identical / len(both_series):.4%})\")\nprint()\nprint(\"cross-tabulation of the two columns:\")\nprint(pd.crosstab(fs, fat, dropna=False).to_string())\nprint()\nif identical == len(both_series):\n print(\">>> The two columns are IDENTICAL on every row.\")\n print(\">>> Two column names, ONE bit of information. Not two.\")"},{"cell_type":"markdown","metadata":{},"source":"### Why one duplicated bit costs you more than it sounds like\n\nIf you build your series slots as *plane x fluid_sensitive x fat_suppression* you will\nbelieve you have 3 x 2 x 2 = 12 kinds of series. You do not. You have 3 x 2 = 6, and half\nyour slots are permanently empty — and a slot that is always empty is a slot your model\nlearns to ignore, not a slot that costs nothing.\n\nWorse, you have lost a real distinction. A **fat-suppressed fluid-sensitive** series (bright\nfluid, fat cancelled — where a bone bruise screams) and a **non-suppressed fluid-sensitive**\nseries (bright fluid, fat also bright) look different and show different things. With one\nbit you cannot tell them apart.\n\n**The recovery, and this is the single highest-value preprocessing idea in this\ncompetition:** the CSV lost the distinction, but the DICOM headers never did. Every slice\ncarries `RepetitionTime`, `EchoTime`, `ScanOptions`, `SeriesDescription` — the settings\nthemselves. You can rebuild both axes from them. That is what\n[Pilkwang Kim's baseline](https://www.kaggle.com/code/pilkwang/rsna-knee-baseline-v1) does,\nand Chapter 7 shows you how."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# The rest of the series table: planes, and how many series a study gets.\n# ---------------------------------------------------------------------------\nprint(\"anatomical planes:\")\nprint(both_series[\"Anatomical_Plane\"].value_counts(dropna=False).to_string())\nprint()\n\nper_study = train_series.groupby(\"StudyInstanceUID\").size()\nprint(\"series per training study:\")\nprint(per_study.describe().to_string())\nprint()\n\nplane_by_study = pd.crosstab(train_series[\"StudyInstanceUID\"],\n train_series[\"Anatomical_Plane\"])\nprint(\"share of studies that have at least one series in each plane:\")\nfor plane in plane_by_study.columns:\n print(f\" {plane:<12} {(plane_by_study[plane] > 0).mean():.1%}\")\nprint()\n\ncombos = (plane_by_study > 0).astype(int).astype(str).agg(\"\".join, axis=1)\nprint(f\"plane availability patterns (columns = {list(plane_by_study.columns)}):\")\nprint(combos.value_counts().head(6).to_string())"},{"cell_type":"markdown","metadata":{},"source":"### What to do with that\n\nRead the plane-availability output above before assuming anything: on this corpus the\nthree planes are far more consistently present than newcomers expect, so the gaps you have\nto handle are **not** at the plane level. They appear once you subdivide a plane by\ncontrast — \"the sagittal *fat-suppressed fluid-sensitive* series\" is a much narrower ask\nthan \"a sagittal series\", and Chapter 7 measures how often each of those actually exists.\n\nWherever the gap is, you need a plan that is not \"crash\" and not \"feed zeros silently\".\nThe standard approach — from\n[Pilkwang Kim's baseline](https://www.kaggle.com/code/pilkwang/rsna-knee-baseline-v1) — is\n**slots plus a presence mask**: define a fixed list of series roles, fill the ones the study\nhas, and pass a 0/1 mask alongside so the model knows the difference between *\"this series\nwas black\"* and *\"this series does not exist\"*. An attention layer over slots then learns\nwhich role matters for which finding, which is medically real: a Baker's cyst is best seen\naxially and posteriorly; a meniscal tear sagittally."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 6 — What is inside a DICOM file?\n\nA `.dcm` file is **one slice of pixels plus a few hundred labelled fields** describing how\nthat slice was made. The fields are the interesting part: this is a scanner writing down its\nown settings, and almost nothing in the CSVs is as reliable as what is in here.\n\nDirectory layout:\n\n```\ntrain_series///.dcm\n```\n\nThe filenames are SOP Instance UIDs — **random identifiers**. Sorting files by name does\n*not* put the slices in anatomical order. Chapter 7 shows what that costs."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Open one file and read the tags that matter.\n# ---------------------------------------------------------------------------\nSERIES_ROOT = ROOT / \"train_series\"\n\nexample_study = sorted(p for p in SERIES_ROOT.iterdir() if p.is_dir())[0]\nexample_series = sorted(p for p in example_study.iterdir() if p.is_dir())[0]\nexample_files = sorted(example_series.glob(\"*.dcm\"))\n\nprint(f\"study {example_study.name}\")\nprint(f\"series {example_series.name} ({len(example_files)} slices)\")\nprint()\n\n# stop_before_pixels=True reads only the header. It is ~100x faster than decoding\n# the image, and for anything that scans the corpus it is the only sane way to do it.\nheader = pydicom.dcmread(str(example_files[len(example_files) // 2]),\n stop_before_pixels=True, force=True)\n\nFIELDS = [\n (\"SeriesDescription\", \"what the radiographer called this sequence\"),\n (\"ScanningSequence\", \"SE / GR / IR - the pulse sequence family\"),\n (\"SequenceVariant\", \"variants applied to it\"),\n (\"ScanOptions\", \"extra options; fat saturation shows up here\"),\n (\"RepetitionTime\", \"TR in ms - short TR means T1 weighting\"),\n (\"EchoTime\", \"TE in ms - long TE means T2 weighting\"),\n (\"MagneticFieldStrength\", \"1.5T or 3T scanner\"),\n (\"Manufacturer\", \"which vendor wrote these tags, and how\"),\n (\"Rows\", \"image height in pixels\"),\n (\"Columns\", \"image width in pixels\"),\n (\"PixelSpacing\", \"MILLIMETRES PER PIXEL - the key to Chapter 7\"),\n (\"SliceThickness\", \"millimetres per slice\"),\n (\"SpacingBetweenSlices\", \"centre-to-centre slice gap in mm\"),\n (\"Laterality\", \"L or R - which knee this is\"),\n (\"ImagePositionPatient\", \"xyz of this slice in patient space, in mm\"),\n (\"ImageOrientationPatient\", \"how the image plane sits in patient space\"),\n (\"RescaleSlope\", \"stored value -> real value: real = slope*stored + intercept\"),\n (\"RescaleIntercept\", \"\"),\n (\"PhotometricInterpretation\", \"MONOCHROME1 means bright and dark are INVERTED\"),\n]\n\nfor tag, note in FIELDS:\n value = getattr(header, tag, None)\n shown = str(value)[:44] if value not in (None, \"\") else \"-\"\n print(f\" {tag:<26} {shown:<46} {note}\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Now the pixels, and the three corrections you apply before looking at them.\n# ---------------------------------------------------------------------------\ndef load_slice(path):\n \"\"\"Read one DICOM into a float array with the standard corrections applied.\n\n 1. Rescale. Stored integers are not physical values until you apply\n RescaleSlope and RescaleIntercept.\n 2. MONOCHROME1. Some series store brightness inverted. Flip them, so that\n \"bright\" always means \"high signal\" everywhere in the corpus.\n 3. No windowing. Deliberately not done here - see the note below.\n \"\"\"\n ds = pydicom.dcmread(str(path), force=True)\n arr = ds.pixel_array.astype(np.float32)\n\n slope = float(getattr(ds, \"RescaleSlope\", 1.0) or 1.0)\n intercept = float(getattr(ds, \"RescaleIntercept\", 0.0) or 0.0)\n arr = arr * slope + intercept\n\n if str(getattr(ds, \"PhotometricInterpretation\", \"\")).upper() == \"MONOCHROME1\":\n arr = arr.max() - arr\n\n return arr, ds\n\n\ndef window(arr, lo=1.0, hi=99.0):\n \"\"\"Map a percentile range onto 0-1 for display.\"\"\"\n p_lo, p_hi = np.percentile(arr, [lo, hi])\n return np.clip((arr - p_lo) / max(p_hi - p_lo, 1e-6), 0, 1)\n\n\nimage, ds = load_slice(example_files[len(example_files) // 2])\nprint(f\"raw pixel array shape {image.shape}, dtype {image.dtype}\")\nprint(f\"value range {image.min():.0f} to {image.max():.0f}\")\n\nspacing = getattr(ds, \"PixelSpacing\", None)\nif spacing is not None:\n mm_per_px = float(spacing[0])\n print(f\"pixel spacing {mm_per_px:.4f} mm/pixel\")\n print(f\"physical size {image.shape[0] * mm_per_px:.1f} x \"\n f\"{image.shape[1] * mm_per_px:.1f} mm of knee\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Look at it. Six slices spread through the stack, in FILENAME order for now -\n# Chapter 7 fixes the ordering.\n# ---------------------------------------------------------------------------\npicks = np.linspace(0, len(example_files) - 1, 6).astype(int)\n\nfig, axes = plt.subplots(1, 6, figsize=(15, 2.9))\nfor ax, i in zip(axes, picks):\n img, _ = load_slice(example_files[i])\n ax.imshow(window(img), cmap=\"gray\")\n ax.set_title(f\"file #{i}\", fontsize=10, color=MUTED)\n ax.axis(\"off\")\nfig.suptitle(f\"One series, six slices - {getattr(ds, 'SeriesDescription', 'unknown')}\",\n fontsize=13, fontweight=\"bold\", x=0.008, ha=\"left\", y=1.06)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","metadata":{},"source":"### Two things about windowing, because this is where beginners lose signal\n\n**Window per *series*, not per slice.** It is tempting to normalise every slice to its own\n1st–99th percentile. Do not. A knee with a large effusion has one very bright structure in\nit; per-slice normalisation rescales that brightness away on exactly the slices where it is\nthe finding you are trying to detect. Compute the percentiles once over the whole series and\napply them to every slice in it.\n\n**Do not use a single global window either.** Intensity in MRI is not calibrated the way\nHounsfield units are in CT — the same tissue on two scanners can be numerically far apart.\nPer-series is the level that is both stable and preserves within-series contrast.\n\n---\n\n## One whole study, all at once\n\nEverything above looked at *one series*. But you do not predict on a series — you predict\ntwelve numbers for a **study**, and a study is several series of the same knee, scanned\ndifferent ways in different orientations.\n\nThis next figure is the one to sit with. It is the entire input for a single prediction."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Every series in one study, side by side, with the geometry that matters.\n#\n# This is the whole picture: one visit, one knee, several sequences. Read the\n# captions as much as the pictures - the mm figure is the physical width of the\n# image, and it is NOT constant across the rows.\n# ---------------------------------------------------------------------------\nstudy_series = sorted(p for p in example_study.iterdir() if p.is_dir())\n\n# The provided table knows the plane; the DICOM header knows everything else.\nplane_of = dict(zip(train_series.SeriesInstanceUID, train_series.Anatomical_Plane))\nfluid_of = dict(zip(train_series.SeriesInstanceUID, train_series.Fluid_Sensitive))\n\npanels = []\nfor sdir in study_series:\n files = sorted(sdir.glob(\"*.dcm\"))\n if not files:\n continue\n img, ds = load_slice(files[len(files) // 2]) # a middle slice\n px = getattr(ds, \"PixelSpacing\", None)\n mm = float(px[0]) if px is not None else float(\"nan\")\n panels.append({\n \"img\": img,\n \"plane\": plane_of.get(sdir.name, \"?\"),\n \"fluid\": fluid_of.get(sdir.name, \"?\"),\n \"desc\": str(getattr(ds, \"SeriesDescription\", \"\") or \"?\")[:22],\n \"n\": len(files),\n \"mm_px\": mm,\n \"fov\": img.shape[1] * mm,\n })\n\nn = len(panels)\nfig, axes = plt.subplots(1, n, figsize=(2.45 * n, 3.6))\naxes = np.atleast_1d(axes)\nfor ax, p in zip(axes, panels):\n ax.imshow(window(p[\"img\"]), cmap=\"gray\")\n ax.set_xticks([]); ax.set_yticks([])\n for s in ax.spines.values():\n s.set_color(GRID)\n ax.set_title(f\"{p['plane']}\\n{p['desc']}\", fontsize=10.5, color=INK, pad=6)\n ax.set_xlabel(f\"{p['n']} slices\\n{p['mm_px']:.2f} mm/px\\n{p['fov']:.0f} mm wide\",\n fontsize=9, color=MUTED, linespacing=1.5)\nfig.suptitle(f\"ONE STUDY = {n} series of the same knee - this is one prediction\",\n fontsize=13, fontweight=\"bold\", x=0.008, ha=\"left\", y=1.02)\nplt.tight_layout()\nplt.show()\n\n# The two facts the picture cannot show you.\nplanes = [p[\"plane\"] for p in panels]\nprint(f\"planes present: {dict((x, planes.count(x)) for x in sorted(set(planes)))}\")\nfovs = [p[\"fov\"] for p in panels if np.isfinite(p[\"fov\"])]\nif fovs:\n spread = max(fovs) / max(min(fovs), 1e-6)\n print(f\"field of view across these series: {min(fovs):.0f} to {max(fovs):.0f} mm \"\n f\"({spread:.2f}x spread)\")\n # State what THIS study shows, not what studies in general show. A single study\n # may well be internally consistent; the corpus is not, and Chapter 7 measures it.\n if spread > 1.05:\n print(\" -> even inside one study the physical scale is not constant\")\n else:\n print(\" -> this study is internally consistent; the variation across the\")\n print(\" corpus is what matters, and Chapter 7 measures it\")"},{"cell_type":"markdown","metadata":{},"source":"Three things to take from that figure, in order of how much they will cost you.\n\n**1. Several series share a plane.** Look at how many say *Sagittal* or *Coronal*. Across\nthe corpus, **4,406 of 4,407 studies have two or more series of the same plane** — so\n\"give me the coronal\" is an ambiguous instruction almost every time. Something has to\nchoose, and the provided `Fluid_Sensitive` column cannot do it, because it is a byte-copy\nof `Fat_Suppression` (Chapter 5). That is why the recovery in Chapter 5 exists.\n\n**2. Read the `mm wide` caption, not the picture.** Two panels can look the same size on\nscreen and cover different amounts of knee — the display stretches every array to the same\nbox, so physical scale is invisible to the eye and only the caption carries it. Whether\n*this* study is internally consistent is printed above; what is certain is that the\n**corpus** is not, and Chapter 7 measures the spread. Feed varying physical scales to a\nnetwork at a fixed pixel size and the same anatomy arrives at different sizes, so the\nencoder must learn a scale invariance it should never have needed.\n\n**3. The same anatomy looks completely different per sequence.** Fluid is bright in one\npanel and dark in another; fat is bright in one and suppressed in another. That is not\nnoise, it is the *point* — a radiologist reads effusion off the fluid-sensitive series and\nligaments off a different one. A model given only one sequence is reading with one eye shut,\nwhich is why the pipeline feeds several \"slots\" per study rather than one image."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 7 — Millimetres, not pixels\n\nThis chapter is the reason to open the DICOM headers at all. We sweep a sample of the\ncorpus, and answer four questions that each change what your preprocessing should do.\n\n1. How big is a pixel, in millimetres — and is it the same everywhere?\n2. Can we rebuild the two acquisition axes the CSV collapsed into one?\n3. Are the slices in order?\n4. Which knee is this?"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# The header sweep. One header per series, read with threads.\n#\n# We read the MIDDLE file of each series rather than the first: the first slice of an\n# MRI stack is often at the very edge of the anatomy, and some series put a localiser\n# or a derived image at one end.\n# ---------------------------------------------------------------------------\nHEADER_TAGS = [\n \"SeriesDescription\", \"SequenceName\", \"ScanOptions\", \"ScanningSequence\",\n \"RepetitionTime\", \"EchoTime\", \"Laterality\", \"ImageLaterality\",\n \"PixelSpacing\", \"Rows\", \"Columns\", \"SliceThickness\",\n \"MagneticFieldStrength\", \"Manufacturer\",\n]\n\n\ndef read_series_header(job):\n \"\"\"Read one series: how many slices it has, plus the tags above.\"\"\"\n study_id, series_id, path = job\n row = {\"StudyInstanceUID\": study_id, \"SeriesInstanceUID\": series_id}\n try:\n files = sorted(f.name for f in os.scandir(path) if f.name.endswith(\".dcm\"))\n row[\"n_slices\"] = len(files)\n if not files:\n return row\n ds = pydicom.dcmread(os.path.join(path, files[len(files) // 2]),\n stop_before_pixels=True, force=True)\n for tag in HEADER_TAGS:\n value = getattr(ds, tag, None)\n if value is None:\n row[tag] = None\n elif isinstance(value, (list, tuple)) or type(value).__name__ == \"MultiValue\":\n row[tag] = \"|\".join(str(v) for v in value)\n else:\n row[tag] = str(value)\n except Exception as exc:\n row[\"error\"] = str(exc)[:80]\n return row\n\n\nall_studies = sorted(p for p in SERIES_ROOT.iterdir() if p.is_dir())\nsampled = [all_studies[i] for i in\n rng.choice(len(all_studies), min(N_HEADER_STUDIES, len(all_studies)),\n replace=False)]\n\njobs = [(study.name, series.name, str(series))\n for study in sampled\n for series in study.iterdir() if series.is_dir()]\n\nstarted = time.time()\nwith ThreadPoolExecutor(max_workers=16) as pool:\n headers = pd.DataFrame(list(pool.map(read_series_header, jobs)))\n\nprint(f\"read {len(headers):,} series headers from {len(sampled):,} studies \"\n f\"in {time.time() - started:.0f}s\")\nprint(f\"failed: {headers.get('error', pd.Series(dtype=object)).notna().sum()}\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Question 1: how big is a pixel, and how big is the field of view?\n#\n# field of view (mm) = number of pixels across x millimetres per pixel\n#\n# This is the number that decides your resize. If two series have different mm/pixel\n# and you resize both to 224, the same anatomy ends up at two different sizes, and\n# your network has to learn to be invariant to something you could have removed.\n# ---------------------------------------------------------------------------\nheaders[\"mm_per_px\"] = pd.to_numeric(\n headers[\"PixelSpacing\"].str.split(\"|\").str[0], errors=\"coerce\")\nheaders[\"cols\"] = pd.to_numeric(headers[\"Columns\"], errors=\"coerce\")\nheaders[\"fov_mm\"] = headers[\"mm_per_px\"] * headers[\"cols\"]\n\ngeom = headers.dropna(subset=[\"mm_per_px\", \"cols\"])\n\nprint(\"millimetres per pixel:\")\nprint(geom[\"mm_per_px\"].describe().to_string())\nprint()\nprint(\"most common acquisition matrices (pixels across):\")\nprint(geom[\"cols\"].value_counts().head(8).to_string())\nprint()\nprint(\"field of view in mm (= mm/pixel x pixels across):\")\nprint(geom[\"fov_mm\"].round(1).describe().to_string())\nprint()\nprint(\"most common field-of-view values:\")\nprint(geom[\"fov_mm\"].round(1).value_counts().head(8).to_string())\nprint()\nprint(\"-\" * 70)\nshare_150 = (geom[\"fov_mm\"].round(0) == 150).mean()\nspread = geom[\"fov_mm\"].quantile([0.05, 0.95])\n\nprint(\"CHECKING A PUBLISHED CLAIM\")\nprint()\nprint(\"Will (wguesdon) reports a CONSTANT 150 mm field of view across this corpus,\")\nprint(\"and builds the millimetres-per-token argument on it. Re-measured here:\")\nprint()\nprint(f\" series at exactly 150 mm {share_150:>7.1%}\")\nprint(f\" distinct values {geom['fov_mm'].round(0).nunique():>7,}\")\nprint(f\" median {geom['fov_mm'].median():>7.1f} mm\")\nprint(f\" 5th - 95th percentile {spread.iloc[0]:>7.1f} - {spread.iloc[1]:.1f} mm\")\nprint()\nif share_150 > 0.9:\n print(\" -> The claim HOLDS on this sample. Input size alone then fixes the\")\n print(\" physical size of a pixel, and no per-series crop is needed.\")\nelse:\n print(\" -> The claim does NOT hold on this sample. Field of view VARIES.\")\n print()\n print(\" This does not weaken the resolution argument - at the median field of\")\n print(\" view a 126 px input is coarser than the 150 mm figure implies, not\")\n print(\" finer. But it does change what you have to DO about it:\")\n print()\n print(\" * you cannot infer physical scale from input size alone;\")\n print(\" * two series resized to the same pixel count are NOT at the same scale;\")\n print(\" * so the millimetre crop below is mandatory, not an optimisation.\")\n print()\n print(\" Crop a fixed physical extent per series first. THEN input size fixes\")\n print(\" millimetres per token, because you made it constant yourself.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Chart: the two geometry distributions side by side.\n# ---------------------------------------------------------------------------\ndef value_bars(ax, counts, label_fmt, color, ylabel):\n \"\"\"One horizontal bar per distinct value, ordered by the value itself.\n\n Ordering by value rather than by frequency matters here: the reader is looking\n for whether the distribution is concentrated, and a bar chart sorted by count\n hides that.\n \"\"\"\n counts = counts.sort_index(ascending=False)\n names = [label_fmt.format(v) for v in counts.index]\n ax.barh(names, counts.values, height=0.62, color=color, zorder=3)\n for i, value in enumerate(counts.values):\n ax.text(value + counts.max() * 0.015, i, f\"{value:,}\",\n va=\"center\", fontsize=10, color=MUTED)\n ax.set_xlim(0, counts.max() * 1.25)\n ax.grid(axis=\"x\", color=GRID, lw=0.8, zorder=0)\n ax.set_xlabel(\"series\")\n ax.set_ylabel(ylabel)\n\n\nspacing_counts = geom[\"mm_per_px\"].round(4).value_counts().head(8)\nfov_counts = geom[\"fov_mm\"].round(0).value_counts().head(8)\n\nn_spacings = geom[\"mm_per_px\"].round(4).nunique()\nn_fovs = geom[\"fov_mm\"].round(0).nunique()\n\nfig, (ax_left, ax_right) = plt.subplots(1, 2, figsize=(13, 4.6))\nvalue_bars(ax_left, spacing_counts, \"{:.4f}\", BLUE, \"mm per pixel\")\nvalue_bars(ax_right, fov_counts, \"{:.0f}\", ORANGE, \"field of view, mm\")\n\n# Titles report what was measured; they do not assert a conclusion the sample may\n# not support. If you re-run this on a different sample the headline follows the data.\nax_left.set_title(f\"Pixel size: {n_spacings} distinct values\", loc=\"left\",\n fontsize=12.5, fontweight=\"bold\")\nax_right.set_title(f\"Field of view: {n_fovs} distinct values\", loc=\"left\",\n fontsize=12.5, fontweight=\"bold\")\n\nif n_fovs < n_spacings:\n caption = (\"Fewer distinct fields of view than pixel sizes: the scanners differ in \"\n \"matrix size, not in how much knee they photographed. So resizing to a \"\n \"fixed pixel count is nearly a fixed physical scale already - but only \"\n \"nearly, and 'nearly' is what the millimetre crop removes.\")\nelse:\n caption = (\"Field of view is not more concentrated than pixel size on this sample, \"\n \"so a fixed-pixel resize does NOT give you a fixed physical scale. Crop \"\n \"by millimetres.\")\nfig.text(0.005, -0.06, caption, fontsize=9.5, color=MUTED, wrap=True)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","metadata":{},"source":"### What to do with pixel spacing\n\n**Crop and resize by millimetres, not by pixels.** Decide on a physical extent — a knee\nfield of view is roughly 140–180 mm — take that many millimetres out of the centre of every\nslice using its own `PixelSpacing`, and only then resize to your network's input size. Now\none pixel means the same physical distance in every study, and the network no longer has to\nspend capacity undoing the difference between a 640-pixel scanner and a 960-pixel one.\n\nThis is the *constant physical scale* idea from\n[Pilkwang Kim's baseline](https://www.kaggle.com/code/pilkwang/rsna-knee-baseline-v1). If\nyou take one preprocessing idea from this notebook, take this one."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Question 2: rebuild the two acquisition axes the CSV collapsed into one.\n#\n# AXIS A - fat suppression. Look for a fat-sat marker in the text tags. Two traps:\n#\n# * underscore is a word character in regex, so \\bwe\\b never matches inside\n# \"t2_de3d_we_tra\". Replace separators with spaces FIRST.\n# * GE writes \"SAT_GEMS\" in ScanOptions for *spatial* saturation, which is a\n# different thing. Match ScanOptions as exact tokens, not as a substring, or\n# you will mark non-suppressed series as suppressed.\n#\n# AXIS B - weighting. Read it from the description when it is stated, and fall back\n# to the physics when it is not: short TR is T1, long TE is T2, long TR with short\n# TE is proton-density.\n# ---------------------------------------------------------------------------\nFATSAT_TOKENS = {\"FS\", \"FATSAT\", \"FAT_SAT\", \"FSAT\"}\nFATSAT_WORDS = re.compile(\n r\"\\bfs\\b|fatsat|fat sat|\\bstir\\b|\\bspair\\b|\\bspir\\b|\\bwe\\b|water excit|\"\n r\"\\btirm\\b|\\bfatsup\\b\")\nT1_WORDS = re.compile(r\"\\bt1\\b|\\bt1w\\b\")\nT2_WORDS = re.compile(r\"\\bt2\\b|\\bt2w\\b\")\nPD_WORDS = re.compile(r\"\\bpd\\b|\\bpdw\\b|proton|\\bdp\\b|dens\")\n\ndescription = (headers[\"SeriesDescription\"].fillna(\"\") + \" \"\n + headers[\"SequenceName\"].fillna(\"\"))\ndescription = description.str.lower().str.replace(r\"[_\\-.]\", \" \", regex=True)\n\nscan_options = headers[\"ScanOptions\"].fillna(\"\").str.upper().str.split(\"|\")\noptions_say_fatsat = scan_options.apply(\n lambda tokens: any(t.strip() in FATSAT_TOKENS for t in tokens))\n\nheaders[\"fat_suppressed\"] = description.str.contains(FATSAT_WORDS) | options_say_fatsat\n\ntr = pd.to_numeric(headers[\"RepetitionTime\"], errors=\"coerce\")\nte = pd.to_numeric(headers[\"EchoTime\"], errors=\"coerce\")\nis_gradient_echo = headers[\"ScanningSequence\"].fillna(\"\").str.upper().str.contains(\"GR\")\n\nheaders[\"weighting\"] = np.select(\n [\n description.str.contains(T1_WORDS) & ~description.str.contains(T2_WORDS)\n & ~description.str.contains(PD_WORDS),\n description.str.contains(T2_WORDS) & ~description.str.contains(PD_WORDS),\n description.str.contains(PD_WORDS),\n is_gradient_echo,\n tr < 800,\n te > 60,\n tr >= 800,\n ],\n [\"T1\", \"T2\", \"PD\", \"GRE\", \"T1\", \"T2\", \"PD\"],\n default=\"unknown\",\n)\nheaders[\"fluid_sensitive\"] = headers[\"weighting\"].isin([\"T2\", \"PD\"])\n\nprint(\"recovered weighting:\")\nprint(headers[\"weighting\"].value_counts().to_string())\nprint()\nprint(\"recovered (fluid-sensitive, fat-suppressed) combinations:\")\nrecovered = pd.crosstab(headers[\"fluid_sensitive\"], headers[\"fat_suppressed\"])\nprint(recovered.to_string())\nprint()\nprint(f\"distinct combinations recovered from the headers: {(recovered > 0).sum().sum()}\")\nprint(\"The CSV gave you 2. Everything above 2 is contrast the CSV threw away.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Now that both axes are back, how often does each SLOT actually exist?\n#\n# A slot is (plane x fluid-sensitive x fat-suppressed). This is the number that\n# tells you how many slots to define: a slot present in 20% of studies is mostly a\n# masked-out zero, and that budget is better spent on more slices in a slot that\n# always fills.\n# ---------------------------------------------------------------------------\nplane_map = dict(zip(both_series[\"SeriesInstanceUID\"], both_series[\"Anatomical_Plane\"]))\nheaders[\"plane\"] = headers[\"SeriesInstanceUID\"].map(plane_map)\n\nslot_name = (headers[\"plane\"].fillna(\"?\") + \" / \"\n + np.where(headers[\"fluid_sensitive\"], \"fluid\", \"structural\") + \" / \"\n + np.where(headers[\"fat_suppressed\"], \"fatsat\", \"no-fatsat\"))\n\nfill = (pd.crosstab(headers[\"StudyInstanceUID\"], slot_name) > 0).mean()\nfill = fill.sort_values(ascending=False)\n\nprint(f\"share of the {headers['StudyInstanceUID'].nunique():,} sampled studies that have \"\n f\"each slot:\")\nprint()\nfor name, share in fill.items():\n bar = \"#\" * int(round(share * 40))\n print(f\" {share:>6.1%} {bar:<40} {name}\")\nprint()\nprint(f\"slots present in over 80% of studies: {(fill > 0.8).sum()} of {len(fill)}\")\nprint(\"Define your slots from this list, not from the product of every flag.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Question 3: are the slices in order?\n#\n# ImagePositionPatient is the xyz of the slice corner in millimetres, in a patient\n# coordinate system. ImageOrientationPatient gives two direction vectors lying in the\n# image plane; their cross product is the direction the stack advances along.\n#\n# Project each slice's position onto that normal, and sorting by the result puts the\n# stack in true anatomical order regardless of how the files are named.\n# ---------------------------------------------------------------------------\ndef slice_positions(series_dir):\n \"\"\"Return (filename, position-along-the-stack-normal) for each slice.\"\"\"\n files = sorted(p for p in Path(series_dir).glob(\"*.dcm\"))\n rows = []\n normal = None\n for path in files:\n ds = pydicom.dcmread(str(path), stop_before_pixels=True, force=True)\n pos = getattr(ds, \"ImagePositionPatient\", None)\n orient = getattr(ds, \"ImageOrientationPatient\", None)\n if pos is None:\n continue\n if normal is None and orient is not None:\n row_dir = np.array([float(v) for v in orient[:3]])\n col_dir = np.array([float(v) for v in orient[3:]])\n normal = np.cross(row_dir, col_dir)\n pos = np.array([float(v) for v in pos])\n offset = float(pos @ normal) if normal is not None else float(pos[2])\n rows.append((path.name, offset))\n return rows\n\n\npositions = slice_positions(example_series)\nif positions:\n names = [n for n, _ in positions]\n offsets = np.array([o for _, o in positions])\n\n filename_order = np.argsort(names)\n physical_order = np.argsort(offsets)\n agreement = float((filename_order == physical_order).mean())\n\n print(f\"slices in this series {len(positions)}\")\n print(f\"position range along the stack {offsets.min():.1f} to {offsets.max():.1f} mm \"\n f\"({offsets.max() - offsets.min():.0f} mm of knee)\")\n print(f\"median gap between slices {np.median(np.diff(np.sort(offsets))):.2f} mm\")\n print()\n print(f\"filename order == physical order for {agreement:.0%} of slices\")\n if agreement < 0.95:\n print()\n print(\">>> Sorting by filename does NOT give you anatomical order.\")\n print(\">>> Filenames are SOP Instance UIDs. They are random.\")\n print(\">>> If you take 'slices 10 to 13' by filename, you get four unrelated\")\n print(\">>> positions in the knee, not a contiguous 2.5D neighbourhood.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# See it. The same six slices, sorted correctly this time.\n# ---------------------------------------------------------------------------\nif positions:\n ordered_names = [names[i] for i in physical_order]\n picks = np.linspace(0, len(ordered_names) - 1, 6).astype(int)\n\n fig, axes = plt.subplots(2, 6, figsize=(15, 6.0))\n for ax, i in zip(axes[0], np.linspace(0, len(names) - 1, 6).astype(int)):\n img, _ = load_slice(example_series / sorted(names)[i])\n ax.imshow(window(img), cmap=\"gray\")\n ax.axis(\"off\")\n for ax, i in zip(axes[1], picks):\n img, _ = load_slice(example_series / ordered_names[i])\n ax.imshow(window(img), cmap=\"gray\")\n ax.set_xticks([])\n ax.set_yticks([])\n for spine in ax.spines.values():\n spine.set_visible(False)\n # An xlabel, not a title: a title on the lower row floats between the two\n # rows and reads as if it belonged to the upper one.\n ax.set_xlabel(f\"{offsets[physical_order][i]:.0f} mm\", fontsize=9.5, color=MUTED)\n\n fig.text(0.005, 0.97, \"Top: sorted by filename. Bottom: sorted by position in the knee.\",\n fontsize=12.5, fontweight=\"bold\")\n fig.text(0.005, 0.0,\n \"Both rows show six slices from the same series. Only the bottom row walks \"\n \"through the joint in order.\", fontsize=9.5, color=MUTED)\n plt.tight_layout(rect=[0, 0.02, 1, 0.95])\n plt.show()"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Question 4: which knee is this?\n#\n# Four of the twelve labels - both menisci, medial OA, lateral OA - are defined\n# relative to the MIDLINE OF THE BODY, not to the image. On a left knee the medial\n# compartment sits on one side of the picture; on a right knee it sits on the other.\n#\n# Feed both to a network without telling it which is which, and \"medial\" is no longer\n# a location it can learn. Four of your twelve labels degrade toward chance.\n# ---------------------------------------------------------------------------\ndef normalise_laterality(series):\n \"\"\"Fold the laterality tags onto exactly {L, R, missing}.\n\n Two traps, both of which split one answer across two buckets:\n * an EMPTY STRING is not NaN, so `.fillna()` leaves it alone. Blank it to NaN\n first or you get a \"missing\" bucket and a \"\" bucket counting the same thing.\n * some scanners write LEFT / RIGHT rather than L / R.\n \"\"\"\n out = series.astype(str).str.strip().str.upper()\n out = out.replace({\"\": np.nan, \"NAN\": np.nan, \"NONE\": np.nan})\n out = out.replace({\"LEFT\": \"L\", \"RIGHT\": \"R\"})\n return out.where(out.isin([\"L\", \"R\"]), \"missing\")\n\n\nprimary = headers[\"Laterality\"].astype(str).str.strip().replace({\"\": np.nan, \"nan\": np.nan})\nlaterality = normalise_laterality(primary.fillna(headers[\"ImageLaterality\"]))\n\ncounts = laterality.value_counts()\nprint(\"which knee, according to the DICOM header:\")\nfor value, count in counts.items():\n print(f\" {value:<9} {count:>6,} ({count / len(laterality):.1%})\")\nprint()\n\nby_study = headers.assign(lat=laterality).groupby(\"StudyInstanceUID\")[\"lat\"].nunique()\nprint(f\"studies in the sample {len(by_study):>5,}\")\nprint(f\"studies with conflicting laterality tags {(by_study > 1).sum():>5,}\")\nprint()\n\nmissing_share = (laterality == \"missing\").mean()\nif missing_share > 0.1:\n print(f\">>> {missing_share:.0%} of series do not record laterality at all.\")\n print(\">>> So 'read the tag and mirror' is not a complete plan. You need a fallback:\")\n print(\">>> take the majority tag across the study's other series, and failing that,\")\n print(\">>> infer the side from the sign of the ImagePositionPatient x-coordinate,\")\n print(\">>> which is negative on the patient's right in the DICOM patient frame.\")\n print()\nprint(\"The fix is to pick a convention - say, make every knee look like a left knee -\")\nprint(\"and mirror the ones that are not. Two cautions:\")\nprint(\" * coronal and axial images: mirror the PIXELS left-right.\")\nprint(\" * sagittal images: the medial/lateral axis runs THROUGH the stack, not across\")\nprint(\" the picture, so you reverse the SLICE ORDER instead of flipping pixels.\")\nprint(\" * once you have done this, do NOT use a horizontal flip augmentation. You would\")\nprint(\" be putting back exactly the nuisance you just removed.\")"},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 8 — DINOv2 101\n\nEveryone in this competition is using DINOv2. This chapter explains what it is, why it is a\nsensible choice here, and what the knobs actually do.\n\n### What it is, in one paragraph\n\n**DINOv2** is a vision model from Meta AI (Oquab et al., 2023,\n[arXiv:2304.07193](https://arxiv.org/abs/2304.07193)). What makes it unusual is that it was\ntrained on **142 million images with no labels at all** — no \"this is a cat\". Instead it was\nshown two different crops of the same photo and trained to produce the same representation\nfor both. To satisfy that, it has to learn what is *in* the picture: object boundaries,\nparts, textures, depth. The result is a general-purpose feature extractor, and the headline\nclaim of the paper is that its features are good enough that a **linear layer** on top of\nfrozen features competes with fully fine-tuned specialist models.\n\nThat claim is exactly what makes it right for this competition. You have very few real\nlabels. A model that works well frozen is a model you can use when you cannot afford to\ntrain one.\n\n### How a vision transformer sees an image\n\nThis is the part worth being concrete about, because every resolution decision follows\nfrom it.\n\n```\ninput image 224 x 224\n |\n | cut into a grid of non-overlapping 14 x 14 pixel PATCHES\n v\n16 x 16 = 256 patches\n |\n | each patch -> one 384-dimensional vector (a \"patch token\")\n | plus one extra learned vector (the \"CLS token\")\n v\n257 tokens -> 12 transformer blocks -> 257 output vectors\n```\n\nThe **CLS token** is the one designed to summarise the whole image; it is what you use if\nyou want a single vector per image. The **patch tokens** are one vector per location, and\nthey retain *where* things are.\n\nThe `14` is the important number. `dinov2-small` is a **ViT-S/14**: patch size 14 pixels.\nSo one token is a 14 x 14 pixel square — and how much *knee* that is depends entirely on\nwhat you resized to."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Load DINOv2 from the Kaggle model mount.\n#\n# Attach it to your notebook via Add Input -> Models -> \"metaresearch/dinov2\",\n# variation \"small\". It arrives in HuggingFace layout (config.json + weights), so\n# transformers.AutoModel reads it directly. local_files_only=True because a code\n# competition runs with internet OFF.\n# ---------------------------------------------------------------------------\nimport torch\n\n\ndef find_dinov2(variation=\"small\"):\n \"\"\"Locate a mounted DINOv2 checkpoint directory.\"\"\"\n base = Path(\"/kaggle/input\")\n if not base.is_dir():\n return None\n found = []\n for folder, subdirs, files in os.walk(base):\n # Do not descend into the 24,000 DICOM directories.\n subdirs[:] = [d for d in subdirs if d not in (\"train_series\", \"test_series\")]\n if \"config.json\" in files and \"dinov2\" in folder.lower():\n found.append(Path(folder))\n for path in found:\n if variation in str(path).lower():\n return path\n return found[0] if found else None\n\n\nDEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nDINO_DIR = find_dinov2(\"small\")\ndino = None\n\nif DINO_DIR is None:\n print(\"DINOv2 not mounted. Add Input -> Models -> metaresearch/dinov2 -> small\")\nelse:\n from transformers import AutoModel\n\n dino = AutoModel.from_pretrained(str(DINO_DIR), local_files_only=True).eval().to(DEVICE)\n for p in dino.parameters():\n p.requires_grad_(False)\n\n PATCH = int(dino.config.patch_size)\n DIM = int(dino.config.hidden_size)\n NATIVE = int(getattr(dino.config, \"image_size\", 0))\n\n print(f\"loaded from {DINO_DIR}\")\n print(f\"device {DEVICE}\")\n print(f\"parameters {sum(p.numel() for p in dino.parameters()):,}\")\n print(f\"patch size {PATCH} px\")\n print(f\"hidden size {DIM} <- the length of every token vector\")\n print(f\"blocks {dino.config.num_hidden_layers}\")\n print(f\"configured at {NATIVE} px, a {NATIVE // PATCH} x {NATIVE // PATCH} patch grid\")"},{"cell_type":"markdown","metadata":{},"source":"### The arithmetic that decides your input size\n\nChapter 7 measured how many millimetres of knee are in one slice, and crucially it measured\nthat **this varies between series**. So do the millimetre crop *first* — take a fixed\nphysical extent out of every slice — and then one number falls out:\n\n$$\\text{mm per patch token} \\;=\\; \\frac{\\text{crop extent in mm}}{\\text{input size in px}} \\times 14$$\n\nNow compare that against the thing you are trying to detect. A medial meniscus is about\n9–12 mm wide. A **tear in it is 1–3 mm.**\n\nThis framing is Will's\n([RSNA Knee DINOv2 at meniscus resolution](https://www.kaggle.com/code/wguesdon/rsna-knee-dinov2-at-meniscus-resolution)),\nand it is the most useful single idea published on this competition so far, because it turns\n\"what input size should I use\" from a hyperparameter you sweep into a question with a\nphysical answer.\n\n**One correction, offered in that spirit.** That notebook states the field of view is a\nconstant 150 mm on every series. Chapter 7 re-measured it and found a spread — check the\noutput yourself. The correction *strengthens* the argument rather than undermining it: at\nthe measured median the tokens are larger than 150 mm implies, so the resolution problem is\nslightly worse than stated. What changes is the remedy. If the field of view were constant,\ninput size alone would fix the physical scale and no crop would be needed. Because it\nvaries, the millimetre crop is load-bearing — without it, two series at the same input size\nsit at different scales and the mm-per-token figure is an average over a distribution rather\nthan a property of your pipeline."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# mm per token at each candidate input size, using the field of view we MEASURED.\n# ---------------------------------------------------------------------------\n# The physical extent you crop to. This is a CHOICE, not a measurement - which is\n# exactly the point. Chapter 7 showed the raw field of view varies, so the pipeline\n# has to impose a constant scale rather than inherit one. 160 mm is the measured\n# median here and comfortably contains a knee.\nCROP_MM = float(np.round(geom[\"fov_mm\"].median() / 10.0) * 10)\n\nFOV_MM = CROP_MM\nPATCH_PX = int(dino.config.patch_size) if dino is not None else 14\nDIM = int(dino.config.hidden_size) if dino is not None else 384\nN_BLOCKS = int(dino.config.num_hidden_layers) if dino is not None else 12\n\n\ndef vit_gflops(n_tokens, dim=DIM, blocks=N_BLOCKS, mlp_ratio=4):\n \"\"\"Forward FLOPs for one image through a plain ViT, counting multiplies and adds.\n\n Per block, with n tokens of width d:\n qkv projection 6*n*d^2\n the two attention matmuls 4*n^2*d <- the quadratic term\n output projection 2*n*d^2\n MLP 4*mlp_ratio*n*d^2\n The patch embedding and the head are under 1% and ignored.\n \"\"\"\n per_block = (6 + 2 + 4 * mlp_ratio) * n_tokens * dim ** 2 + 4 * n_tokens ** 2 * dim\n return blocks * per_block / 1e9\n\n\nprint(f\"physical crop this pipeline imposes: {CROP_MM:.0f} mm \"\n f\"(median measured field of view was {geom['fov_mm'].median():.1f} mm)\")\nprint(f\"DINOv2 patch {PATCH_PX} px, width {DIM}, {N_BLOCKS} blocks\")\nprint()\n\nbaseline_flops = vit_gflops((126 // PATCH_PX) ** 2 + 1)\nrows = []\nfor size in [126, 168, 224, 280, 336, 448]:\n n_tokens = (size // PATCH_PX) ** 2\n mm_pixel = FOV_MM / size\n rows.append({\n \"input_px\": size,\n \"tokens\": n_tokens,\n \"mm_per_pixel\": round(mm_pixel, 3),\n \"mm_per_token\": round(FOV_MM * PATCH_PX / size, 2),\n \"3mm_tear_in_px\": round(3.0 / mm_pixel, 2),\n \"gflops\": round(vit_gflops(n_tokens + 1), 1),\n \"cost_vs_126px\": round(vit_gflops(n_tokens + 1) / baseline_flops, 1),\n })\n\ntable = pd.DataFrame(rows)\nprint(table.to_string(index=False))\nprint()\nprint(\"Read the '3mm_tear_in_px' column. Below about 1 pixel, the tear is not merely\")\nprint(\"small in the image - it has been erased by the resize before the model sees it.\")\nprint(\"'cost_vs_126px' is measured FLOPs, not a guess; the n^2 attention term is why\")\nprint(\"it grows faster than the token count does.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Chart: physical size of a token against input size, with the anatomy marked.\n# ---------------------------------------------------------------------------\nsizes = np.arange(112, 476, PATCH_PX)\nmm_token = FOV_MM * PATCH_PX / sizes\n\nfig, ax = plt.subplots(figsize=(10.5, 5.4))\nax.axhspan(9, 12, color=BLUE, alpha=0.13, zorder=0)\nax.axhspan(1, 3, color=ORANGE, alpha=0.16, zorder=0)\nax.plot(sizes, mm_token, color=INK, lw=2, zorder=3)\n\nfor size in [126, 224, 336]:\n value = FOV_MM * PATCH_PX / size\n ax.plot([size], [value], \"o\", ms=8, color=INK, zorder=4)\n ax.annotate(f\"{size} px\\n{value:.1f} mm/token\", (size, value),\n textcoords=\"offset points\", xytext=(10, 12), fontsize=10, color=INK)\n\nax.text(468, 10.5, \"meniscus, 9-12 mm wide\", ha=\"right\", va=\"center\",\n fontsize=10.5, color=\"#1a5aa8\")\nax.text(468, 2.0, \"a tear in it, 1-3 mm\", ha=\"right\", va=\"center\",\n fontsize=10.5, color=\"#c04b1c\")\n\nax.set_xlim(112, 476)\nax.set_ylim(0, max(14, mm_token.max() * 1.2))\nax.grid(color=GRID, lw=0.8, zorder=0)\nax.set_xlabel(\"input size, pixels\")\nax.set_ylabel(\"millimetres of knee per patch token\")\nax.set_title(\"One DINOv2 token covers this much anatomy\", loc=\"left\",\n fontsize=13, fontweight=\"bold\", pad=14)\n\n# The caption follows the measurement rather than asserting it, so this figure can\n# never claim something its own curve contradicts.\nmm_at_126 = FOV_MM * PATCH_PX / 126\nif mm_at_126 >= 9:\n verdict = \"At 126 px one token swallows the whole meniscus.\"\nelif mm_at_126 >= 3:\n verdict = f\"At 126 px one token is {mm_at_126:.1f} mm - larger than any tear.\"\nelse:\n verdict = (f\"At 126 px one token is {mm_at_126:.1f} mm, already inside the lesion \"\n f\"range on this sample.\")\nfig.text(0.005, -0.03,\n f\"Computed from a {CROP_MM:.0f} mm physical crop and a {PATCH_PX} px patch. \"\n f\"{verdict}\",\n fontsize=9.5, color=MUTED)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","metadata":{},"source":"There is a second, quieter reason not to go too small. DINOv2 ships with a positional\nembedding grid sized for its configured resolution — printed above as the native patch grid.\nWhen you feed it a different size, it **interpolates that grid**. Going from a 37 x 37 grid\ndown to 9 x 9 (which is what 126 px does) throws away most of the positional structure the\nmodel was trained with. A larger input is not the exotic setting here; it is the *closer*\none.\n\nThe counterweight is cost — look at the `cost_vs_126px` column. Attention is quadratic in\nthe number of tokens, so going from 126 px to 336 px costs *more* than the 7x increase in\ntoken count. That is the real trade, and it is a budget decision, not a correctness one."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# Run it. One real knee slice in, tokens out.\n# ---------------------------------------------------------------------------\nIMAGENET_MEAN = np.array([0.485, 0.456, 0.406], np.float32).reshape(3, 1, 1)\nIMAGENET_STD = np.array([0.229, 0.224, 0.225], np.float32).reshape(3, 1, 1)\n\n\ndef to_model_input(arr_2d, size=224):\n \"\"\"One 2D slice -> a normalised [1, 3, size, size] tensor.\n\n DINOv2 expects 3 channels and ImageNet normalisation. Here the same greyscale\n slice is repeated across all three. In a real pipeline you would instead stack\n three NEIGHBOURING slices - that is what \"2.5D\" means, and it hands the model\n a little depth for free.\n \"\"\"\n x = torch.from_numpy(window(arr_2d).astype(np.float32))[None, None]\n x = torch.nn.functional.interpolate(x, size=(size, size),\n mode=\"bilinear\", align_corners=False)\n x = x.repeat(1, 3, 1, 1)\n x = (x - torch.from_numpy(IMAGENET_MEAN)) / torch.from_numpy(IMAGENET_STD)\n return x\n\n\nif dino is not None:\n IMG_SIZE = 224\n slice_img, _ = load_slice(example_files[len(example_files) // 2])\n batch = to_model_input(slice_img, IMG_SIZE).to(DEVICE)\n\n with torch.no_grad():\n out = dino(pixel_values=batch).last_hidden_state\n\n grid = IMG_SIZE // PATCH_PX\n print(f\"input {tuple(batch.shape)}\")\n print(f\"output {tuple(out.shape)}\")\n print(f\" = 1 image x (1 CLS + {grid}x{grid}={grid * grid} \"\n f\"patches) x {out.shape[-1]} dims\")\n print()\n\n cls_token = out[:, 0]\n patch_tokens = out[:, 1:]\n print(f\"CLS token {tuple(cls_token.shape)} <- one vector for the image\")\n print(f\"patch tokens {tuple(patch_tokens.shape)} <- one vector per location\")\n print()\n print(\"Three ways to reduce this to a fixed-length feature vector:\")\n print(f\" CLS only {cls_token.shape[-1]} dims\")\n print(f\" CLS + mean over patches {cls_token.shape[-1] * 2} dims <- a good default\")\n print(f\" CLS + mean + max over patches {cls_token.shape[-1] * 3} dims\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# What did it actually learn? The classic DINOv2 demonstration.\n#\n# Take the patch tokens, run PCA down to 3 dimensions, and paint those three as RGB\n# back onto the image grid. Nobody told this model what a femur is. If the picture\n# below comes out with coherent coloured regions that follow the anatomy, that is\n# the self-supervised representation separating tissue types on its own.\n# ---------------------------------------------------------------------------\nif dino is not None:\n from sklearn.decomposition import PCA\n\n def token_pca_image(arr_2d, size=224):\n \"\"\"Patch tokens -> 3-component PCA -> an RGB image on the patch grid.\"\"\"\n with torch.no_grad():\n tokens = dino(pixel_values=to_model_input(arr_2d, size).to(DEVICE))\n tokens = tokens.last_hidden_state[0, 1:].float().cpu().numpy()\n side = size // PATCH_PX\n reduced = PCA(n_components=3, random_state=SEED).fit_transform(tokens)\n lo, hi = np.percentile(reduced, [2, 98], axis=0)\n reduced = np.clip((reduced - lo) / np.maximum(hi - lo, 1e-6), 0, 1)\n return reduced.reshape(side, side, 3)\n\n show_idx = np.linspace(0, len(example_files) - 1, 4).astype(int)\n fig, axes = plt.subplots(2, 4, figsize=(13, 6.6))\n for column, i in enumerate(show_idx):\n img, _ = load_slice(example_files[i])\n axes[0][column].imshow(window(img), cmap=\"gray\")\n axes[0][column].axis(\"off\")\n axes[1][column].imshow(token_pca_image(img), interpolation=\"nearest\")\n axes[1][column].axis(\"off\")\n\n fig.text(0.005, 0.965, \"What DINOv2 sees without ever being trained on a knee\",\n fontsize=13.5, fontweight=\"bold\")\n fig.text(0.005, 0.93, \"Top: the MRI slice. Bottom: the first three principal \"\n \"components of its patch tokens, painted as red/green/blue.\",\n fontsize=10, color=MUTED)\n fig.text(0.005, 0.0,\n \"Each coloured square is one patch token. Regions of one colour are regions \"\n \"the model considers one kind of thing.\", fontsize=9.5, color=MUTED)\n plt.tight_layout(rect=[0, 0.02, 1, 0.91])\n plt.show()"},{"cell_type":"markdown","metadata":{},"source":"### Freeze it, or fine-tune it?\n\nBoth are defensible here and the public notebooks split on it. The honest summary:\n\n| | **Freeze** | **Fine-tune the last few blocks** |\n|---|---|---|\n| What you do | Run DINOv2 once, cache the vectors, train only a small head | Unfreeze the last ~4–6 transformer blocks at a low learning rate (~1e-5), train end to end |\n| Cost | Encode once, then every head experiment takes seconds | Every experiment is a full GPU run |\n| Risk | The features were learned on natural photographs. Nothing in that data resembles a torn meniscus on a proton-density sequence | Your training labels come from *reports*, not radiologists. Fine-tuning on noisy labels partly teaches the model the noise |\n| Do this when | You are iterating on the head, the aggregation, the label extractor — i.e. most of the time | Frozen features have stopped improving no matter what you change |\n\n**The diagnostic that tells you which situation you are in:** if raising the resolution,\nusing a bigger backbone, adding slices and improving your aggregation head *all* stop\nhelping at the same point, the bottleneck is the representation itself, and that is when\nunfreezing pays. Until then it does not.\n\nTwo rules if you do unfreeze:\n\n- **Only the last blocks.** Early transformer blocks are generic edge and texture detectors.\n You do not have enough supervision to improve them and you have plenty to damage them.\n- **Two learning rates.** The backbone wants something like 1e-5; a freshly initialised head\n wants 1e-3. One learning rate for both wrecks whichever one it was not chosen for.\n\n### The mistake that costs the most, and makes no noise\n\n```python\nmodel = timm.create_model(\"resnet18\", pretrained=False) # <-- random weights\n```\n\n`pretrained=False` gives you a **randomly initialised network**. It trains, the loss goes\ndown, and it produces a submission. It has learned nothing useful, and there is no error\nmessage anywhere. In a code competition with internet off, `pretrained=True` silently fails\nthis way too if the weights are not cached.\n\n**Assert it.** Check that the weights loaded, and raise loudly if they did not. A random\nbackbone that runs to completion is far more expensive than one that stops."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 9 — The checklist\n\nEverything above, as things to check before you spend a GPU hour. Roughly in the order the\nmistakes cost you.\n\n### Data\n\n- [ ] **Do not trust `Fluid_Sensitive` and `Fat_Suppression` as two columns.** They are one\n bit. Rebuild both axes from `RepetitionTime`, `EchoTime`, `ScanOptions` and\n `SeriesDescription`. *(Ch. 4, 6)*\n- [ ] **Crop by millimetres, not pixels.** Use `PixelSpacing` to take a fixed physical extent\n from every slice, then resize. *(Ch. 6)*\n- [ ] **Sort slices by `ImagePositionPatient`, not by filename.** Filenames are random UIDs.\n *(Ch. 6)*\n- [ ] **Normalise laterality — and have a fallback for the series that do not record it.**\n Mirror every knee onto one convention (pixels for coronal/axial, slice order for\n sagittal), then do *not* use a horizontal flip augmentation. A large share of series\n carry no laterality tag at all, so fall back to the study's other series and then to\n the sign of the `ImagePositionPatient` x-coordinate. *(Ch. 6)*\n- [ ] **Window per series, not per slice.** Per-slice normalisation rescales away the\n brightness of an effusion on exactly the slices where it is the finding. *(Ch. 5)*\n- [ ] **Apply `RescaleSlope`/`RescaleIntercept` and handle `MONOCHROME1`.** *(Ch. 5)*\n- [ ] **No vertical flip.** A knee has the femur above and the tibia below. Flipping it\n creates an anatomy that does not exist, and three of your twelve labels are defined on\n that axis. Small rotation, scale and shift are the augmentations that correspond to how\n patients are actually positioned.\n\n### Labels and validation\n\n- [ ] **Train on report-derived labels; keep the annotated studies for checking.** *(Ch. 2)*\n- [ ] **Group your folds on the report hash.** Duplicate report text across studies leaks a\n validation answer straight into training. *(Ch. 3)*\n- [ ] **Hold the annotated studies out of their own fold.** Scoring your annotated reference\n over studies the model trained on gives you a number that only goes up — and if you\n select a checkpoint on it, it selects the most overfitted one.\n- [ ] **Report the uncertainty on your validation AUC.** With this few labelled studies it is\n roughly ±0.07 per label. Do not chase +0.005. *(Ch. 2)*\n- [ ] **Measure the *fire rate* of every rule, not only its accuracy.** A rule that never\n matches emits a negative, silently, and that gap tracks language, which tracks the\n site. *(Ch. 3)*\n\n### Model\n\n- [ ] **Verify your pretrained weights actually loaded.** Raise if not. *(Ch. 7)*\n- [ ] **Choose your input size from millimetres per token, not by sweeping.** *(Ch. 7)*\n- [ ] **Combine models by rank, not by probability.** The metric reads order only. *(Ch. 1)*\n- [ ] **Give absent series a presence mask.** \"This series is black\" and \"this series does\n not exist\" must not look the same to your model. *(Ch. 4)*\n- [ ] **Spend effort on the rare labels.** Each of the twelve is worth exactly $1/12$. *(Ch. 1)*\n\n### Kaggle mechanics\n\n- [ ] **Pin the accelerator to T4.** `\"machine_shape\": \"NvidiaTeslaT4\"` in\n `kernel-metadata.json`. With `enable_gpu` alone you may be handed a P100, whose compute\n capability (6.0) is below what the pinned PyTorch build ships kernels for.\n `torch.cuda.is_available()` returns `True` anyway, so the run proceeds all the way\n through data loading and dies at the *first* kernel launch with\n `no kernel image is available for execution on the device`.\n- [ ] **Write a valid `submission.csv` before you do anything expensive**, then overwrite it.\n A run killed for memory takes a `SIGKILL` that no `try`/`except` catches, and a\n submission that never wrote scores nothing at all.\n- [ ] **Project your memory before you allocate it.** Cache size is\n studies x slots x slices x resolution². If it will not fit, drop *slices* before\n resolution: coverage is linear in memory, resolution is quadratic."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Where to go next\n\nThree things in this notebook were **stated but not measured**, and each is a genuinely open\nquestion worth a notebook of its own. If you want to contribute something useful, these are\nopen:\n\n1. **Does resolution actually help, on this data?** Chapter 8 argues from sampling theory\n that a 1–3 mm tear cannot survive a 126 px resize. That is a prediction. Measuring it —\n same head, same folds, input size as the only variable — would settle it.\n2. **How much does correct slice ordering buy?** Chapter 7 shows filename order is wrong.\n Nobody has published the score difference.\n3. **Which labels are actually reachable?** Some of the twelve may be near-invisible on the\n available sequences. Per-label OOF AUC, published, would tell everyone where the remaining\n points are.\n\nAnd one thing this notebook **did** measure that the public discussion currently has wrong:\nthe field of view is not constant across this corpus (Chapter 7). If you have been treating\ninput size as a proxy for physical scale, it is not one until you impose the millimetre crop\nyourself.\n\n### If this was useful\n\nUpvote the notebooks in the credits block at the top — several of them found these traps in\npublic, before the rest of us walked into them, and that is the only reason this notebook\ncould be written.\n\nCorrections are very welcome. If something here is wrong, say so in the comments and it gets\nfixed in the next version with credit to you."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 10 — From findings to code that runs\n\nA notebook full of findings proves nothing until the findings turn into code. This chapter\nis the bridge: the four functions from our baseline that each exist *because* of a chapter\nabove, run here on the real data so you can see them work.\n\nNone of this needs a GPU. It is the preprocessing, not the model."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# FIX 1 (Chapter 6) - order slices by geometry, not by filename.\n#\n# Filenames are SOP Instance UIDs, so sorted filename order is a random\n# permutation of the anatomy. Project each slice's ImagePositionPatient onto the\n# stack normal (the cross product of the two ImageOrientationPatient vectors) and\n# sort by that instead.\n# ---------------------------------------------------------------------------\ndef order_slices(directory, files):\n \"\"\"Return `files` sorted by physical position along the stack.\"\"\"\n keyed, missing = [], []\n normal = None\n for name in files:\n ds = pydicom.dcmread(str(Path(directory) / name),\n stop_before_pixels=True, force=True)\n pos = getattr(ds, \"ImagePositionPatient\", None)\n orient = getattr(ds, \"ImageOrientationPatient\", None)\n if pos is None:\n missing.append(name)\n continue\n if normal is None and orient is not None:\n r = np.array([float(v) for v in orient[:3]])\n c = np.array([float(v) for v in orient[3:]])\n normal = np.cross(r, c)\n p = np.array([float(v) for v in pos])\n keyed.append((float(p @ normal) if normal is not None else float(p[2]), name))\n keyed.sort()\n return [n for _, n in keyed] + sorted(missing)\n\n\nnames = [f.name for f in example_files]\nordered = order_slices(example_series, sorted(names))\nsame = sum(a == b for a, b in zip(sorted(names), ordered))\nprint(f\"slices {len(names)}\")\nprint(f\"filename order matches physical order for {same}/{len(names)} positions\")\nprint(f\"first slice by FILENAME ...{sorted(names)[0][-16:]}\")\nprint(f\"first slice by GEOMETRY ...{ordered[0][-16:]}\")\nprint(\"same file\" if sorted(names)[0] == ordered[0] else\n \">>> different files - sorting by name gives you the wrong end of the knee\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# FIX 2 (Chapter 7) - crop to a constant number of MILLIMETRES, then resize.\n#\n# The subtle bug: a guard written as `if want < min(h, w)` reduces algebraically\n# to `if FOV > crop_mm`, so every series whose field of view is at or below\n# crop_mm silently skips normalisation - and the median FOV here is 160 mm,\n# exactly the usual crop_mm. Pad instead of skipping.\n# ---------------------------------------------------------------------------\ndef crop_to_mm(vol, mm_per_px, crop_mm=160.0, out=224):\n \"\"\"One millimetre maps to a fixed number of output pixels, always.\"\"\"\n want = int(round(crop_mm / mm_per_px))\n h, w = vol.shape[-2:]\n if want < min(h, w):\n cy, cx, half = h // 2, w // 2, want // 2\n vol = vol[..., cy - half:cy + half, cx - half:cx + half]\n how = \"cropped\"\n else:\n py, px_ = max(0, want - h), max(0, want - w)\n vol = np.pad(vol, [(0, 0)] * (vol.ndim - 2) +\n [(py // 2, py - py // 2), (px_ // 2, px_ - px_ // 2)])\n how = \"PADDED (field of view smaller than crop_mm)\"\n # The resize to `out` px is what the encoder actually receives; what matters is\n # that it now always represents exactly crop_mm of knee.\n return vol, how, crop_mm / out\n\n\nimg, ds = load_slice(example_files[len(example_files) // 2])\nmm = float(ds.PixelSpacing[0])\nprint(f\"this series: {img.shape[0]}x{img.shape[1]} px at {mm:.3f} mm/px \"\n f\"= {img.shape[0] * mm:.0f} mm field of view\")\nout_vol, how, mm_out = crop_to_mm(img[None], mm)\nprint(f\" -> {how}, then resized to 224 px = {mm_out:.3f} mm/pixel\")\nprint(f\" -> one DINOv2 patch (14 px) now covers {mm_out * 14:.2f} mm, ALWAYS,\")\nprint(f\" whatever the acquisition was. That is the whole point.\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# FIX 3 (Chapter 2) - resolve which knee it is, or four targets are unlearnable.\n#\n# Medial Meniscus, Lateral Meniscus, Medial OA and Lateral OA are side-specific -\n# a third of the metric. A study with no laterality is never mirrored, so\n# \"medial\" lands on opposite sides of the image for left and right knees.\n#\n# THE TRAP: every sagittal protocol here is named `sag_...`, and `sag` is also\n# Turkish for \"right\". Match it naively and you mirror most of the corpus.\n# ---------------------------------------------------------------------------\nSIDE_RX = [(\"R\", re.compile(r\"\\b(right|rt|dexter|sağ|derech[ao]|rechts?|droite?)\\b\", re.I)),\n (\"L\", re.compile(r\"\\b(left|lt|sinister|sol|izquierd[ao]|links?|gauche)\\b\", re.I))]\n\n\ndef resolve_side(ds):\n \"\"\"Laterality tag -> ImageLaterality -> side words -> geometry (LPS: +x is left).\"\"\"\n for tag in (\"Laterality\", \"ImageLaterality\"):\n v = str(getattr(ds, tag, \"\") or \"\").strip().upper()\n if v[:1] in (\"L\", \"R\"):\n return v[0], tag\n blob = \" \".join(str(getattr(ds, t, \"\") or \"\") for t in\n (\"SeriesDescription\", \"StudyDescription\", \"BodyPartExamined\"))\n hits = {s for s, rx in SIDE_RX if rx.search(re.sub(r\"[^\\w\\s]\", \" \", blob))}\n if len(hits) == 1:\n return hits.pop(), \"description text\"\n pos = getattr(ds, \"ImagePositionPatient\", None)\n if pos is not None and abs(float(pos[0])) > 20.0:\n return (\"R\" if float(pos[0]) < 0 else \"L\"), \"geometry (ImagePositionPatient x)\"\n return None, \"nothing\"\n\n\n# Check every series in this study, and check the trap holds.\nfor sdir in sorted(p for p in example_study.iterdir() if p.is_dir()):\n f = sorted(sdir.glob(\"*.dcm\"))\n if not f:\n continue\n d = pydicom.dcmread(str(f[len(f) // 2]), stop_before_pixels=True, force=True)\n side, source = resolve_side(d)\n print(f\" {str(getattr(d, 'SeriesDescription', '?'))[:26]:<28} \"\n f\"side={side or '-'} from {source}\")\n\ntrap = re.sub(r\"[^\\w\\s]\", \" \", \"sag_pd_fs\")\nprint(f\"\\ntrap check: does 'sag_pd_fs' read as Turkish 'right'? \"\n f\"{'YES - BUG' if any(rx.search(trap) for s, rx in SIDE_RX if s == 'R') else 'no'}\")"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# ---------------------------------------------------------------------------\n# FIX 4 (Chapter 4) - fold by NORMALISED report text, not by study id.\n#\n# The derived label is a deterministic function of the report, so two studies\n# sharing a report share a target. Split them across folds and the text model\n# memorises one and predicts its twin - which inflates the very number used to\n# decide how much to trust it.\n#\n# Measured on this corpus: 131 byte-identical duplicate reports, and 25 further\n# groups (125 studies) that differ ONLY in punctuation, casing or measurement\n# digits - which a raw md5 scatters across folds.\n# ---------------------------------------------------------------------------\ndef normalise_report(s):\n s = re.sub(r\"[^\\w\\s]\", \" \", str(s).lower())\n return re.sub(r\"\\s+\", \" \", re.sub(r\"\\d+\", \"#\", s)).strip()\n\n\ndef report_fold(report, n_folds=5):\n key = normalise_report(report) or str(report)\n return int(hashlib.md5(key.encode()).hexdigest()[:8], 16) % n_folds\n\n\nreports = train[\"Report\"].fillna(\"\")\nraw = reports.map(lambda r: int(hashlib.md5(r.encode()).hexdigest()[:8], 16) % 5)\nnorm_key = reports.map(normalise_report)\nsplit = (pd.DataFrame({\"g\": norm_key, \"f\": raw}).groupby(\"g\")[\"f\"].nunique() > 1)\nn_split = int(split.sum())\nn_studies = int(norm_key.isin(split[split].index).sum())\n\nprint(f\"exact-identical report groups {reports.nunique():>6,}\")\nprint(f\"after normalising {norm_key.nunique():>6,}\")\nprint(f\"\\ngroups a RAW md5 splits across folds: {n_split} ({n_studies} studies)\")\nnew = reports.map(report_fold)\nstill = int((pd.DataFrame({\"g\": norm_key, \"f\": new}).groupby(\"g\")[\"f\"].nunique() > 1).sum())\nprint(f\"groups the NORMALISED hash splits: {still}\")\nprint(\"\\nfold sizes:\", [int((new == f).sum()) for f in range(5)])"},{"cell_type":"markdown","metadata":{},"source":"## What those four are worth\n\nHonest accounting, because this is the part people overstate.\n\n| Fix | What kind of claim is that? |\n|---|---|\n| Slice ordering | **Correctness.** Filename order matched physical order 5% of the time — the \"2.5D triplets\" were three unrelated positions in the knee. Score effect: **unmeasured.** |\n| Millimetre crop | **Correctness.** The guard skipped every series at or below `crop_mm`. Score effect: **unmeasured.** |\n| Laterality | **Correctness.** Four side-specific targets are unlearnable on unmirrored studies. Score effect: **unmeasured.** |\n| Fold by normalised report | **Leakage.** 25 groups / 125 studies were split. Effect on the reported number: it was optimistic. |\n\nEvery one is *\"this code did not do what it claimed\"*. **Not one of them is a claim that\nthe score goes up.** Those are different statements and it matters that they stay\ndifferent — the only thing that can promote a correctness fix to a performance claim is a\ncross-validated measurement.\n\nWhich needs a ruler. We built one by running the same configuration twice with only the\ntraining seed changed — same folds, same targets, same studies — so sampling noise cancels\nand what remains is training randomness alone:\n\n```\nseed 2026 : 0.7552 seed 7 : 0.7541\n noise floor = 0.0010\n```\n\nAnything below 0.0010 on the macro AUC is not a result. And per label the floor is about\n**0.03** — thirty times looser — so a single label moving a couple of points is noise, not\nprogress. Two different gates, and confusing them is the most common way a competition gets\nwasted.\n\n**The full pipeline** — the six-slot cache, the DINOv2 head, the rank fusion — is a separate\nnotebook, because it needs a GPU and a couple of hours. This one stays CPU-only and runs in\nminutes, on purpose."},{"cell_type":"markdown","metadata":{},"source":"---\n\n# Chapter 11 — Now run it\n\nEverything above is measurement. This chapter is the pipeline those measurements produced,\nend to end, and it writes a `submission.csv` you can submit.\n\n**It needs a GPU and about an hour.** Chapters 1–10 are CPU-only and take minutes — if you\nonly came for the findings, you can stop reading here and still have everything the\nnotebook measured.\n\nWhat the code below carries, and which chapter each piece came from:\n\n| From | What it does |\n|---|---|\n| Ch. 5 | Recovers fat-suppression and weighting from the DICOM headers, because `Fluid_Sensitive` is a byte-copy of `Fat_Suppression` |\n| Ch. 6 | Sorts slices by `ImagePositionPatient` projected on the stack normal — filename order matches physical order 5% of the time |\n| Ch. 7 | Crops a constant **160 mm**, padding when the field of view is smaller, so one DINOv2 patch always covers the same 10 mm of knee |\n| Ch. 2 | Mirrors every knee onto a left-knee convention, so *medial* is always on the same side — four of the twelve targets are side-specific |\n| Ch. 4 | Assigns folds by a hash of the **normalised** report, so 25 near-duplicate groups stay whole |\n| Ch. 3 | Builds targets from the reports, then reports OOF per label with `n_pos` and a standard error beside each |\n\n**One number to keep in view while you read it.** Running this twice with only the training\nseed changed — same folds, same targets, same studies — the macro AUC moved by **0.0010**.\nThat is the noise floor. A change that moves the mean by less than that has not done\nanything, and per label the floor is roughly **0.03**, thirty times looser. Two gates, not\none."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"\nimport os\n\n# ---- edit here -------------------------------------------------------------------- #\nos.environ.setdefault(\"RSNA_FOLDS\", \"1\") # 1 to validate, 5 for the real run\nos.environ.setdefault(\"RSNA_TIME_BUDGET\", \"14400\") # seconds; raise to 28800 for 5 folds\nos.environ.setdefault(\"RSNA_TRAIN_SEED\", \"2026\") # bump alone to measure training noise\n# os.environ[\"RSNA_BACKBONE\"] = \"resnet50\" # a timm name, overrides DINOv2\n# os.environ[\"SLOT_SCHEME\"] = \"public\" # ablate the header recovery\n# ------------------------------------------------------------------------------------ #\n\nfor k in (\"RSNA_FOLDS\", \"RSNA_TIME_BUDGET\", \"RSNA_TRAIN_SEED\", \"RSNA_BACKBONE\",\n \"SLOT_SCHEME\"):\n if k in os.environ:\n print(f\"{k} = {os.environ[k]}\")","id":"cell-02"},{"cell_type":"markdown","metadata":{},"source":"## 2. Where the targets come from\n\n`train.csv` has a `Report` column and `test.csv` does not. Text is available when\nfitting and absent when predicting, which rules out any fusion model with a text branch\nand leaves the reports usable only as a source of training targets and of sample\nweights.\n\nTwo facts shape the extractor:\n\n**Reports are graded, annotations are thresholded.** The reporting radiologist and the\nannotator do not share a cutoff — a report saying *small joint effusion* can sit against\na negative annotation. So a rule of the form *term present ⇒ positive* is wrong by\nconstruction; grading the mention (trace / unqualified / marked) is right, and costs\nnothing, because §1 established that only order is read.\n\n**Silence is not a negative.** A finding the report never mentions gets a low score and\na low *confidence*, and the confidence becomes the sample weight. A study whose report\nsays nothing about synovitis pulls on the synovitis head far less than one that names it.","id":"cell-03"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"\"\"\"Report -> twelve graded targets, in nine languages.\n\nProvenance: the multilingual lexicon and the assert/negate/hedge scoring are adapted\nfrom the community notebook `rsna-knee-baseline-v1`. What is added here is §2 of this\nmodule - `distill`, which fits a text model to the *rule* scores out-of-fold and blends\nit back in.\n\nThat addition exists because of a failure mode the rule extractor cannot fix from\ninside. A rule that never fires does not raise; it emits a negative. A lexicon that is\nthick in English and thin in Greek therefore looks like a corpus where Greek patients\nhave fewer findings, and because language tracks the reporting site, that bias is\naligned with a scanner and a population rather than averaging out as noise.\n\nA character n-gram model fitted to the rule scores over the *whole* corpus sees the\nGreek phrasings that co-occur with rule-positive reports and scores them, without anyone\nhaving written them into the lexicon. Fitted out-of-fold and grouped on report text, it\ncannot simply memorise the rules it was trained on.\n\"\"\"\n\nfrom __future__ import annotations\n\nimport re\nimport unicodedata\n\nimport numpy as np\n\nTARGETS = [\n \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n \"Medial OA\", \"Lateral OA\", \"PF OA\", \"Effusion\",\n \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\",\n]","id":"cell-04"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# Turkish dotted/dotless i must be folded before casefolding, otherwise \"İZLENMEZ\"\n# and \"izlenmez\" diverge. ss and the Croatian/Serbian d-with-stroke likewise.\n_PRE = str.maketrans({\n \"ı\": \"i\", \"İ\": \"i\", \"I\": \"i\", \"ß\": \"ss\", \"đ\": \"d\", \"Đ\": \"d\",\n \"ø\": \"o\", \"Ø\": \"o\", \"æ\": \"ae\", \"Æ\": \"ae\",\n})\n\n\ndef normalize(text: str) -> str:\n \"\"\"Fold case, diacritics and separators; keep Greek and Cyrillic letters.\n\n NFKD decomposition strips Latin accents and Greek tonos alike (a -> a), which is what\n we want: reports are inconsistent about accents. It also maps the MICRO SIGN U+00B5\n to a real mu, which matters because most Greek reports here use the wrong codepoint.\n \"\"\"\n if not isinstance(text, str):\n return \"\"\n text = text.translate(_PRE).lower()\n text = unicodedata.normalize(\"NFKD\", text)\n text = \"\".join(ch for ch in text if not unicodedata.combining(ch))\n text = text.replace(\"­\", \"\") # soft hyphen\n text = re.sub(r\"[_\\-/\\\\]+\", \" \", text)\n text = re.sub(r\"[ \\t]+\", \" \", text)\n return text\n\n\n_SENT_SPLIT = re.compile(r\"(?<=[.;!?])\\s+|\\n+\")\n\n\ndef clauses(text: str):\n \"\"\"Split into clauses, then attach `header:` lines to the value that follows.\n\n A report line reading `Fractures :` followed by `Aucune.` is one statement. Splitting\n on punctuation alone separates the anatomy from its negation and flips the label.\n \"\"\"\n norm = normalize(text)\n raw = [c.strip() for c in _SENT_SPLIT.split(norm) if c and c.strip()]\n\n merged = []\n for i, c in enumerate(raw):\n # A fragment ending in a colon is a heading for the next fragment. Structured\n # English reports write long ones - \"lateral compartment (meniscus, collateral\n # ligament complex, cartilage):\" is eight words - so the cap is generous.\n #\n # The heading is emitted ONLY joined to its value, never also on its own. A bare\n # `Fractures :` contains the anatomy and no polarity cue, so scoring it alone\n # reads it as an assertion - and the very next line is `Aucune.` Structured\n # reports are built out of exactly this shape, so keeping the bare heading turns\n # every negated section header into a false positive. The joined clause is a\n # superset of the heading's text, so nothing is lost by dropping it.\n if c.endswith(\":\") and len(c.split()) <= 14 and i + 1 < len(raw):\n merged.append(c + \" \" + raw[i + 1])\n continue\n merged.append(c)\n # Comma-separated enumerations inside a long clause hide separate assertions.\n out = []\n for c in merged:\n out.append(c)\n if len(c.split()) > 25:\n out.extend(p.strip() for p in c.split(\",\") if len(p.split()) > 2)\n return out\n\n\ndef _rx(*alts: str) -> re.Pattern:\n return re.compile(\"|\".join(alts))\n\n\nNEGATION = _rx(\n # `none` and `nil` matter more than their frequency suggests: a structured report\n # writes the anatomy as a heading and the finding as a one-word value beneath it,\n # so `FRACTURE:` / `None.` is the entire statement and missing that one word flips\n # the label on every section that was checked and found clear.\n r\"\\bno\\b\", r\"\\bnot\\b\", r\"\\bnone\\b\", r\"\\bnil\\b\", r\"\\babsent\\b\",\n r\"\\bwithout\\b\", r\"\\bnegative\\b\", r\"\\babsence\\b\",\n r\"\\bno evidence\\b\", r\"\\bunremarkable\\b\", r\"\\bfree of\\b\",\n r\"\\bsin\\b\", r\"\\bno hay\\b\", r\"\\bausencia\\b\", r\"\\bausentes?\\b\", r\"\\bningun[ao]?\\b\",\n r\"\\bpas de\\b\", r\"\\bsans\\b\", r\"\\baucune?\\b\", r\"\\brien\\b\",\n r\"\\bgeen\\b\", r\"\\bzonder\\b\", r\"\\bniet\\b\",\n r\"\\bkeine?\\b\", r\"\\bohne\\b\", r\"\\bnicht\\b\",\n r\"\\byok\\b\", r\"\\byoktur\\b\", r\"izlenmemekte\", r\"saptanmadi\", r\"\\bdegil\\b\",\n r\"gozlenmemekte\", r\"mevcut degil\", r\"eslik etmiyor\", r\"\\bizlenmedi\\b\",\n r\"\\bnema\\b\", r\"\\bbez\\b\", r\"\\bnisu\\b\", r\"\\bnije\\b\",\n r\"\\bδεν\\b\", r\"\\bχωρις\\b\", r\"ουδεν\",\n r\"\\bбез\\b\", r\"\\bне\\b\", r\"липсва\", r\"\\bняма\\b\",\n)\n\nNORMALITY = _rx(\n r\"\\bnormal\", r\"\\bintact\\b\", r\"\\bpreserved\\b\", r\"\\bwithin normal limits\\b\",\n r\"limites normales\", r\"\\bconservad\", r\"\\bintegr\", r\"\\bnormales\\b\",\n r\"\\bdoga(l|ll)\\b\", r\"korunmus\", r\"\\bnormaldir\\b\", r\"olagan\",\n r\"\\buredn\", r\"\\bocuvan\", r\"\\bodrzan\", r\"\\bintakt\",\n r\"φυσιολογικ\", r\"ακεραι\",\n r\"unauffallig\", r\"regelrecht\",\n r\"нормал\", r\"запазен\", r\"съхранен\", r\"\\bбез особености\\b\",\n r\"\\bgaaf\\b\", r\"\\bnormaal\\b\",\n)\n\nUNCERTAIN = _rx(\n r\"\\bpossible\\b\", r\"\\bprobable\\b\", r\"\\bsuspicious\\b\", r\"\\bsuspected\\b\",\n r\"cannot (be )?exclude\", r\"\\bmay\\b\", r\"\\bquestionable\\b\", r\"\\bequivocal\\b\",\n r\"\\bposible\\b\", r\"sin criterios categoricos\", r\"\\bdudos\",\n r\"\\bmuhtemel\\b\", r\"\\bolasi\\b\", r\"\\bsupheli\\b\", r\"\\bizlenim\",\n r\"\\bmoguce\\b\", r\"\\bvjerojatno\\b\", r\"\\bsumnja\\b\",\n r\"πιθαν\", r\"υποπτ\",\n r\"\\bmoglich\", r\"\\bverdachtig\", r\"\\bfraglich\", r\"\\bv\\.a\\.\\b\",\n r\"\\bвъзможно\\b\", r\"\\bвероятно\\b\", r\"суспект\",\n r\"\\bmogelijk\\b\", r\"\\bverdacht\\b\",\n)\n\nTEAR = _rx(\n # `disrupt`, not `\\bdisruption\\b`: \"the ACL is disrupted\" is the commonest English\n # phrasing for a complete tear and the noun form misses it entirely.\n r\"\\btear\", r\"\\btorn\\b\", r\"\\brupture\", r\"disrupt\", r\"discontinuit\",\n r\"\\bavuls\", r\"\\blacerat\",\n r\"\\brotura\\b\", r\"\\broturas\\b\", r\"\\bruptura\", r\"\\bdesgarro\", r\"\\broto\\b\",\n r\"\\bdechirure\", r\"\\bdechire\",\n r\"\\bscheur\", r\"\\bruptuur\", r\"gescheurd\",\n # German compounds the pathology onto the anatomy - Innenmeniskusriss,\n # Kreuzbandruptur, Meniskusabriss - so these stems must not carry a leading word\n # boundary or the entire German subcorpus reads as negative.\n r\"riss\\b\", r\"einriss\", r\"ruptur\", r\"zerreiss\", r\"\\blasion\",\n r\"kontinuitatsunterbrechung\",\n r\"\\byirtik\", r\"\\byirtig\", r\"\\bkopma\\b\", r\"butunluk kaybi\",\n r\"\\bpuknuce\", r\"\\bprekid\\b\", r\"\\bpukotin\",\n r\"ρηξη\", r\"ρηξις\", r\"ρηγμα\",\n r\"руптура\", r\"разкъсв\", r\"разрив\", r\"скъсв\",\n)\n\nDEGEN = _rx(\n r\"degenerat\", r\"\\bmucoid\\b\", r\"\\bmyxoid\\b\", r\"\\bfray\", r\"\\bfissur\",\n r\"dejeneratif\", r\"\\bmukoid\\b\", r\"degenerativn\", r\"εκφυλιστ\", r\"дегенерат\",\n r\"\\bmuco ?ide\\b\", r\"aufgefasert\",\n)\n\nINJURY = _rx(\n r\"\\binjur\", r\"\\bsprain\", r\"\\blesion\", r\"\\blasion\", r\"\\bedema\\b\", r\"\\boedema\\b\",\n r\"\\bodem\\b\", r\"\\bedem\\b\", r\"\\bοιδημα\", r\"\\bодем\", r\"\\bедем\", r\"\\bstrain\\b\",\n r\"\\bhigh signal\\b\", r\"\\bsignal alteration\\b\", r\"\\bhiperintens\", r\"\\bhyperintens\",\n r\"\\bthicken\", r\"\\bzadebljanje\\b\", r\"\\bverdikking\\b\", r\"\\bdistenzij\",\n r\"\\blaksite\\b\", r\"\\blaxity\\b\", r\"\\bpartial\\b\", r\"\\bparcijaln\", r\"\\bparcial\",\n r\"\\bpartiel\", r\"\\bpartiell\",\n)\n\nANAT = {\n \"ACL\": _rx(\n r\"anterior cruciate\", r\"\\bacl\\b\",\n r\"cruzado anterior\", r\"\\blca\\b\",\n r\"croise anterieur\",\n r\"voorste kruisband\", r\"\\bvkb\\b\",\n r\"vorderes kreuzband\", r\"vorderen kreuzband\", r\"vordere kreuzband\",\n r\"on capraz\", r\"\\bocb\\b\",\n r\"prednji krizni\", r\"prednjeg krizn\",\n r\"προσθι[οα][^ ]* χιαστ\", r\"προσθιου χιαστου\", r\"χιαστο[^ ]* συνδεσμ\",\n r\"предна кръстна\", r\"предната кръстна\",\n # Plural, unqualified: reports routinely clear both cruciates in one clause\n # (\"Ligamentos cruzados y colaterales dentro de limites normales\").\n r\"cruciate ligaments\", r\"ligamentos cruzados\", r\"ligaments croises\",\n r\"kruisbanden\", r\"kreuzbander\", r\"capraz baglar\", r\"krizn[a-z]* ligament[a-z]*\",\n r\"χιαστοι συνδεσμ\", r\"χιαστων συνδεσμ\", r\"кръстните връзки\", r\"кръстни връзки\",\n ),\n \"MCL\": _rx(\n r\"medial collateral\", r\"\\bmcl\\b\", r\"tibial collateral\",\n r\"colateral medial\", r\"colateral interno\", r\"\\blcm\\b\",\n r\"collateral medial\", r\"collateral interne\",\n r\"mediale collaterale\", r\"binnenband\",\n r\"innenband\", r\"mediales? kollateral\",\n r\"\\bic yan bag\", r\"medial kollateral\", r\"\\biyb\\b\",\n r\"medijalni kolateraln\", r\"medijalnog kolateraln\",\n r\"εσω πλαγι\", r\"εσωτερικο πλαγι\",\n r\"медиален колатерал\", r\"вътрешна странична\",\n # \"Ligamentos cruzados y colaterales\" separates the noun from its adjective, so\n # the adjective has to stand alone as a cue.\n r\"\\bcolaterales\\b\", r\"\\bcollateraux\\b\", r\"\\bcollateralen\\b\", r\"\\bkolateralni\\b\",\n r\"collateral ligaments\", r\"ligamentos colaterales\", r\"ligaments collateraux\",\n r\"collaterale banden\", r\"kollateralbander\", r\"seitenbander\", r\"yan baglar\",\n r\"kolateraln[a-z]* ligament[a-z]*\", r\"πλαγιοι συνδεσμ\", r\"πλαγιων συνδεσμ\",\n r\"колатерални връзки\", r\"страничните връзки\",\n ),\n \"Medial Meniscus\": _rx(\n r\"medial meniscus\", r\"\\bmm\\b(?= tear)\", r\"medial menisc\",\n r\"menisco medial\", r\"menisco interno\",\n r\"menisque medial\", r\"menisque interne\",\n r\"mediale meniscus\", r\"binnenmeniscus\",\n r\"innenmeniskus\", r\"medialen? meniskus\", r\"innenmeniskushinterhorn\",\n r\"medyal menisk\", r\"\\bic menisk\",\n r\"medijalni meniskus\", r\"medijalnog meniskusa\", r\"medijalnom meniskusu\",\n r\"εσω μηνισκ\", r\"μηνισκ[^ ]* του εσω\", r\"εσω διαμερισμα[^.]{0,40}μηνισκ\",\n r\"медиалния менискус\", r\"медиален менискус\", r\"вътрешния менискус\",\n ),\n \"Lateral Meniscus\": _rx(\n r\"lateral meniscus\", r\"lateral menisc\",\n r\"menisco lateral\", r\"menisco externo\",\n r\"menisque lateral\", r\"menisque externe\",\n r\"laterale meniscus\", r\"buitenmeniscus\",\n r\"aussenmeniskus\", r\"lateralen? meniskus\",\n r\"lateral menisk\", r\"\\bdis menisk\",\n r\"lateralni meniskus\", r\"lateralnog meniskusa\", r\"lateralnom meniskusu\",\n r\"εξω μηνισκ\", r\"μηνισκ[^ ]* του εξω\", r\"εξω διαμερισμα[^.]{0,40}μηνισκ\",\n r\"латералния менискус\", r\"латерален менискус\", r\"външния менискус\",\n ),\n}\n\n# Osteoarthritis is rarely written as \"osteoarthritis\". It is written as cartilage loss,\n# chondropathy grade, joint space narrowing, or osteophytes - scoped to a compartment.\nOA_EVIDENCE = _rx(\n r\"osteoarthrit\", r\"\\barthros\", r\"\\bgonarthros\", r\"\\bosteoarthros\",\n r\"chondropath\", r\"chondromalac\", r\"condropat\", r\"condromalac\",\n r\"cartilage loss\", r\"cartilage thinning\", r\"chondral (loss|defect|ulcer|thinning)\",\n r\"osteophyt\", r\"osteofit\", r\"osteofyt\", r\"osteofito\", r\"osteophyten\",\n r\"joint space narrowing\", r\"pinzamiento articular\",\n r\"kikirdak kayb\", r\"kikirdak incelme\", r\"kondropati\", r\"kondral\",\n r\"kraakbeen(lijden|verlies)\", r\"gonartrose\", r\"artrose\",\n r\"knorpel(verlust|schaden|defekt)\", r\"arthrose\", r\"gonarthrose\",\n r\"hrskavic\", r\"hondromalac\", r\"artroz\", r\"osteoartrit\",\n r\"χονδρ[^ ]*παθ\", r\"αρθριτ\", r\"αρθρωσ\", r\"οστεοφυτ\",\n r\"αρθρικου χονδρου\", r\"εξαλειψη του αρθρικου χονδρου\",\n r\"артроз\", r\"хондропат\", r\"остеофит\", r\"хрущял[^.]{0,30}(изтън|увред|дефект)\",\n r\"ulcera[s]? condral\", r\"cartilago[^.]{0,25}(perdida|adelgaz)\",\n r\"icrs grade\", r\"outerbridge\",\n)\n\nCOMPARTMENT = {\n \"Medial OA\": _rx(\n r\"medial (femorotibial|tibiofemoral|compartment)\",\n r\"compartimento femorotibial medial\", r\"femorotibial interno\",\n r\"mediaal femorotibiaal\", r\"mediale femorotibial\",\n r\"medial femorotibial\", r\"medialen kompartiment\", r\"innere[sn]? kompartiment\",\n r\"medyal femorotibial\", r\"ic kompartman\", r\"medyal kompartman\",\n r\"medijaln[^ ]* (femorotibi|odjelj|kompartm)\",\n r\"εσω διαμερισμα\", r\"εσω κνημιαι\", r\"εσω μηριαι\",\n r\"медиалн[^ ]* (компартм|отдел|тибиал|феморотиб)\",\n r\"medial (femoral|tibial) (condyle|plateau)\", r\"condilo femoral medial\",\n r\"medialen? (femurkondyl|tibiaplateau)\", r\"mediale femorale condyl\",\n ),\n \"Lateral OA\": _rx(\n r\"lateral (femorotibial|tibiofemoral|compartment)\",\n r\"compartimento femorotibial lateral\", r\"femorotibial externo\",\n r\"lateraal femorotibiaal\", r\"laterale femorotibial\",\n r\"lateral femorotibial\", r\"lateralen kompartiment\", r\"aussere[sn]? kompartiment\",\n r\"dis kompartman\", r\"lateral kompartman\",\n r\"lateraln[^ ]* (femorotibi|odjelj|kompartm)\",\n r\"εξω διαμερισμα\", r\"εξω κνημιαι\", r\"εξω μηριαι\",\n r\"латералн[^ ]* (компартм|отдел|тибиал|феморотиб)\",\n r\"lateral (femoral|tibial) (condyle|plateau)\", r\"condilo femoral lateral\",\n r\"lateralen? (femurkondyl|tibiaplateau)\", r\"laterale femorale condyl\",\n ),\n \"PF OA\": _rx(\n r\"patellofemoral\", r\"femoropatellar\", r\"femoropatelar\", r\"patelofemoral\",\n r\"retropatellar\", r\"retrorotulian\", r\"\\btrochlea\", r\"\\btroclea\", r\"\\btroklea\",\n r\"\\bpatella\\b\", r\"\\bpatellar\\b\", r\"\\brotulian\", r\"\\brotula\\b\", r\"\\bpatele\\b\",\n r\"\\bpatellae?\\b\", r\"patellofemoraal\", r\"femoropatellair\",\n r\"επιγονατιδ\", r\"μηροεπιγονατιδ\", r\"τροχιλ\",\n r\"пател\", r\"феморопател\", r\"тролх\",\n r\"anterior compartment\", r\"compartimento anterior\", r\"prednj[^ ]* odjeljk\",\n ),\n}\n\n# Self-declaring findings: the term itself is the finding.\nDIRECT = {\n \"Effusion\": _rx(\n r\"\\beffusion\", r\"joint fluid\", r\"intra ?articular fluid\", r\"\\bhydrops\\b\",\n r\"derrame articular\", r\"\\bderrame\\b\", r\"liquido articular\",\n r\"epanchement\",\n r\"gewrichtsvocht\", r\"\\bvocht\\b\", r\"gewrichtseffusie\",\n r\"gelenkerguss\", r\"\\berguss\\b\", r\"gelenksergu\",\n # \"diz eklemi ici sivi miktari ... artmis\" - the noun takes a possessive suffix,\n # so `eklem ` alone misses. Match the stem plus any suffix.\n r\"eklem\\w* ic\\w* sivi\", r\"efuzyon\", r\"eklem sivisi\",\n r\"sivi (miktari|artisi|birikimi)\", r\"sivi artis\", r\"\\bsivi\\b[^.]{0,25}artmis\",\n r\"\\bizljev\", r\"\\bizliv\", r\"zglobn[^ ]* tekucin\", r\"\\bhidrops\\b\",\n r\"αρθρικ[^ ]* υγρ\", r\"υγρου ενδαρθρικα\", r\"ενδαρθρικ[^ ]* υγρ\", r\"ποσοτητα υγρου\",\n r\"ενδαρθρικ\", r\"αρθρικη συλλογη\", r\"υγρο στην αρθρωση\", r\"υγρου στην αρθρωση\",\n r\"ставен излив\", r\"излив\", r\"ставна течност\", r\"синовиална течност\",\n ),\n \"Synovitis\": _rx(\n r\"synovit\", r\"sinovit\", r\"synovial (thickening|proliferation|hypertroph)\",\n r\"synoviale? (verdikking|proliferatie)\",\n r\"synovialitis\", r\"synovialis(verdickung|proliferation)\",\n r\"sinovijalitis\", r\"zadebljanje sinovij\",\n r\"υμενιτιδα\", r\"συνοβιτιδα\", r\"υμενικ[^ ]* υπερτροφ\", r\"αρθρικου υμεν\",\n r\"синовит\", r\"синовиал[^ ]* (задебел|пролифер)\",\n r\"verdikkingen van (het )?synovium\", r\"pannus\",\n ),\n \"Baker's\": _rx(\n r\"baker\", r\"popliteal cyst\", r\"quiste popliteo\", r\"quistes popliteos\",\n r\"kyste poplite\", r\"popliteale? cyst\", r\"poplitealzyste\", r\"bakerzyste\",\n r\"popliteal kist\", r\"\\bbakerova\\b\", r\"poplitealn[^ ]* cist\",\n r\"κυστη baker\", r\"πολυχωρη συνοβιακη κυστη\", r\"κυστη του baker\",\n r\"киста на бейкър\", r\"бейкърова киста\", r\"поплитеална киста\",\n r\"gastrocnemio ?semimembranos\", r\"gastrocnemius semimembranosus burs\",\n ),\n \"Contusion\": _rx(\n r\"\\bcontusion\", r\"bone bruise\", r\"bone marrow (o?edema|contusion)\",\n r\"\\bkontuz\", r\"medular bone o?edema\", r\"marrow o?edema\",\n r\"edema oseo\", r\"edema de medula osea\",\n r\"oedeme osseux\",\n r\"botcontusie\", r\"botoedeem\", r\"beenmergoedeem\", r\"botmergoedeem\",\n r\"knochenmarkodem\", r\"knochenodem\", r\"kontusion\",\n r\"kemik kontuzyonu\", r\"kemik iligi odemi\", r\"kemik odemi\",\n r\"kostani edem\", r\"edem kosti\", r\"kontuzij\",\n r\"οστεομυελικ[^ ]* οιδημα\", r\"οστικο οιδημα\", r\"μυελικο οιδημα\",\n r\"костномозъчен едем\", r\"костен едем\", r\"контузионен\",\n ),\n \"Fracture\": _rx(\n r\"\\bfractur\", r\"\\bfract\\b\",\n r\"\\bfractura\", r\"\\bfracturas\\b\",\n r\"\\bfractuur\", r\"\\bbreuk\\b\",\n r\"\\bfraktur\", r\"\\bbruch\\b\",\n r\"\\bkirik\\b\", r\"\\bkirigi\\b\",\n r\"\\bprijelom\", r\"impresijsk[^ ]* fraktur\",\n r\"καταγμα\", r\"καταγματ\",\n r\"фрактур\", r\"счупван\", r\"фисур\",\n r\"insufficiency fracture\", r\"stress fracture\", r\"avulsion fracture\",\n r\"subchondral fracture\", r\"subkondral kiri\",\n ),\n}\n\n# Terms that look like a finding but are not the finding being scored.\nDECOY = {\n \"Fracture\": _rx(r\"no fracture\", r\"microfractur\", r\"\\bfracture (risk|prophyla)\"),\n \"Baker's\": _rx(r\"meniscal cyst\", r\"quiste meniscal\", r\"ganglion\"),\n}\n\nPAIRED = {\"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\"}\nOA_TARGETS = {\"Medial OA\", \"Lateral OA\", \"PF OA\"}\n\nSTEM_MENISCUS = _rx(r\"menisc\\w*\", r\"menisk\\w*\", r\"μηνισκ\\w*\", r\"мениск\\w*\")\nSTEM_CRUCIATE = _rx(r\"cruciate\", r\"cruzado\", r\"croise\", r\"kruisband\", r\"kreuzband\",\n r\"capraz bag\\w*\", r\"krizn\\w*\", r\"χιαστ\\w*\", r\"кръстн\\w*\",\n r\"\\bacl\\b\", r\"\\bpcl\\b\", r\"\\blca\\b\", r\"\\blcp\\b\", r\"\\bvkb\\b\",\n r\"\\bhkb\\b\", r\"\\bocb\\b\", r\"\\bacb\\b\")\nSTEM_COLLATERAL = _rx(r\"collateral\\w*\", r\"colateral\\w*\", r\"kollateral\\w*\",\n r\"collaterale\\w*\", r\"kolateraln\\w*\", r\"yan bag\\w*\",\n r\"πλαγι\\w*\", r\"колатерал\\w*\", r\"странич\\w*\",\n r\"innenband\\w*\", r\"aussenband\\w*\", r\"binnenband\\w*\",\n r\"\\bmcl\\b\", r\"\\blcl\\b\", r\"\\blcm\\b\", r\"\\biyb\\b\")\n\nSIDE_MEDIAL = _rx(r\"\\bmedial\\w*\", r\"\\bmedyal\\w*\", r\"\\bmedijaln\\w*\", r\"\\bmediaal\\w*\",\n r\"\\bmediale\\w*\", r\"\\binterno\\w*\", r\"\\binterne\\w*\", r\"\\binnen\\w*\",\n r\"\\bic\\b\", r\"\\bunutarnj\\w*\", r\"\\bεσω\\w*\", r\"\\bεσωτερικ\\w*\",\n r\"\\bмедиал\\w*\", r\"\\bвътреш\\w*\", r\"\\btibial collateral\\b\")\nSIDE_LATERAL = _rx(r\"\\blateral\\w*\", r\"\\bexterno\\w*\", r\"\\bexterne\\w*\", r\"\\bdis\\b\",\n r\"\\blateraln\\w*\", r\"\\baussen\\w*\", r\"\\bbuiten\\w*\", r\"\\bεξω\\w*\",\n r\"\\bεξωτερικ\\w*\", r\"\\bлатерал\\w*\", r\"\\bвъншн\\w*\",\n r\"\\bfibular collateral\\b\", r\"\\bvanjsk\\w*\")\nSIDE_ANTERIOR = _rx(r\"\\banterior\\w*\", r\"\\bant\\b\", r\"\\bon\\b\", r\"\\bprednj\\w*\",\n r\"\\bvorder\\w*\", r\"\\bvoorste\\b\", r\"\\bπροσθι\\w*\", r\"\\bпредн\\w*\",\n r\"\\banteriyor\\w*\", r\"\\bavant\\b\", r\"\\banterieur\\w*\")\n\n# Fracture is the target whose stem varies most across the corpus.\nSTEM_FRACTURE = _rx(r\"fractur\\w*\", r\"fraktur\\w*\", r\"fractuur\\w*\", r\"\\bfract\\b\",\n r\"kiri[kg]\\w*\", r\"prijelom\\w*\", r\"lom kosti\", r\"\\bbreuk\\w*\",\n r\"\\bbruch\\w*\", r\"καταγμα\\w*\", r\"καταγματ\\w*\", r\"фрактур\\w*\",\n # NOT a bare `fissur\\w*`: \"fisuras condrales\" and \"full thickness\n # fissures in the articular cartilage\" describe cartilage, not bone.\n r\"счупван\\w*\", r\"fisur\\w* (osea|oseas|kost)\", r\"fissur\\w* kost\")\n\nSTEM_OA_COMPARTMENT = _rx(r\"compartment\\w*\", r\"compartimento\\w*\", r\"compartiment\\w*\",\n r\"kompartman\\w*\", r\"kompartiment\\w*\", r\"odjelj\\w*\",\n r\"διαμερισμα\\w*\", r\"компартм\\w*\", r\"\\bотдел\\w*\",\n r\"femorotibial\\w*\", r\"femorotibiaal\\w*\", r\"tibiofemoral\\w*\",\n r\"femoro tibial\\w*\", r\"κνημιαι\\w*\", r\"μηριαι\\w*\",\n r\"femoral condyl\\w*\", r\"tibial plateau\\w*\",\n r\"condilo femoral\", r\"platillo tibial\", r\"tibiaplateau\\w*\",\n r\"femurkondyl\\w*\", r\"femoralne? kondil\\w*\",\n r\"tibijaln\\w* plato\", r\"femoral kondil\\w*\",\n r\"tibia plato\", r\"tibyal plato\")\n\n\ndef _near(clause: str, stem_rx: re.Pattern, qual_rx: re.Pattern, window: int = 55):\n \"\"\"True if a stem match has a qualifier within `window` characters either side.\n\n Character windows rather than token windows, because word order differs: English\n puts the side before the noun, Greek and Bulgarian often after, and Turkish\n attaches it as a separate preceding adjective.\n \"\"\"\n for m in stem_rx.finditer(clause):\n lo = max(0, m.start() - window)\n hi = min(len(clause), m.end() + window)\n if qual_rx.search(clause[lo:hi]):\n return True\n return False\n\n\nSTEM_RULES = {\n \"ACL\": (STEM_CRUCIATE, SIDE_ANTERIOR),\n \"MCL\": (STEM_COLLATERAL, SIDE_MEDIAL),\n \"Medial Meniscus\": (STEM_MENISCUS, SIDE_MEDIAL),\n \"Lateral Meniscus\": (STEM_MENISCUS, SIDE_LATERAL),\n \"Medial OA\": (STEM_OA_COMPARTMENT, SIDE_MEDIAL),\n \"Lateral OA\": (STEM_OA_COMPARTMENT, SIDE_LATERAL),\n}\n\nSEV_LOW = _rx(\n r\"\\bsmall\\b\", r\"\\bminimal\\b\", r\"\\btrace\\b\", r\"\\bmild\\b\", r\"\\bslight\\b\",\n r\"\\btiny\\b\", r\"\\bscant\\b\", r\"\\bmimimal\\b\", r\"\\bdiscrete\\b\", r\"\\bfocal\\b\",\n r\"\\bleve\\b\", r\"\\bminim\", r\"\\bpeque\", r\"\\bligero\\b\", r\"\\bescaso\\b\", r\"\\bdiscreto\\b\",\n r\"\\bhafif\\b\", r\"\\baz miktarda\\b\", r\"\\bsilik\\b\",\n r\"\\bmanja\\b\", r\"\\bmanji\\b\", r\"\\bblago\\b\", r\"\\bdiskretn\", r\"\\bmalo\\b\",\n r\"\\bgering\", r\"\\bdiskret\", r\"\\bkleine?r?\\b\", r\"\\bwenig\\b\", r\"\\bzarte?\\b\",\n r\"\\bbeperkte?\\b\", r\"\\bgeringe\\b\", r\"\\bweinig\\b\", r\"\\blichte?\\b\",\n r\"\\bηπι\", r\"\\bμικρ\", r\"\\bελαχιστ\",\n r\"\\bминимал\", r\"\\bлек\", r\"\\bмалк\", r\"\\bнеголям\",\n)\n\nSEV_HIGH = _rx(\n r\"\\blarge\\b\", r\"\\bmarked\\b\", r\"\\bmassive\\b\", r\"\\bsevere\\b\", r\"\\bextensive\\b\",\n r\"\\bmoderate\\b\", r\"\\bgross\\b\", r\"\\bsignificant\\b\", r\"\\babundant\\b\", r\"\\btense\\b\",\n r\"\\bmoderad\", r\"\\bimportante\\b\", r\"\\bsevera?\\b\", r\"\\bmarcad\", r\"\\bcuantios\",\n r\"\\bbelirgin\\b\", r\"\\byaygin\\b\", r\"\\bileri\\b\", r\"\\bciddi\\b\", r\"\\bbol\\b\",\n r\"\\bopsezan\\b\", r\"\\bveliki\\b\", r\"\\bizrazit\", r\"\\bznacajn\", r\"\\bumjeren\",\n r\"\\bausgepragt\", r\"\\bdeutlich\", r\"\\bmassiv\", r\"\\bmassig\", r\"\\bgross\",\n r\"\\buitgebreid\", r\"\\bgevorderd\", r\"\\bveel\\b\", r\"\\bmatige?\\b\",\n r\"\\bμετρι\", r\"\\bμεγαλ\", r\"\\bεκτεταμεν\", r\"\\bευμεγεθ\", r\"\\bσοβαρ\",\n r\"\\bголям\", r\"\\bизразен\", r\"\\bзначим\", r\"\\bумерен\", r\"\\bобилен\",\n)\n\n# OA is often asserted for the whole joint rather than per compartment\n# (\"tricompartmental osteoarthritis\", \"gonarthrose\"). Those statements are evidence for\n# all three OA targets.\nGLOBAL_OA = _rx(\n r\"tri ?compartment\", r\"all three compartment\", r\"global(ised)? (oa|osteoarthrit)\",\n r\"\\bgonarthros\", r\"\\bgonartros\", r\"\\bgonarthrose\", r\"\\bgonartrose\",\n r\"osteoarthritis of the knee\", r\"artrosis (de |)(la )?rodilla\", r\"knee osteoarthrit\",\n r\"\\bdiz osteoartrit\", r\"\\bgonartroz\", r\"artroza koljena\",\n r\"οστεοαρθριτιδα\", r\"αρθριτιδα του γονατος\",\n r\"артроза на колянната\", r\"гонартроз\",\n r\"degenerative joint disease\", r\"\\bdjd\\b\",\n)\n\n# A bare \"bone marrow oedema\" is not a contusion when it sits under a cartilage defect:\n# subchondral oedema beneath a worn compartment is reactive degenerative signal, and\n# reading it as a bruise turns every osteoarthritic knee into a trauma case.\nDEGENERATIVE_MARROW = _rx(\n r\"subchondral\", r\"subcondral\", r\"subkondral\", r\"supkondraln\", r\"subchondraln\",\n r\"υποχονδρι\", r\"субхондрал\",\n r\"\\bcyst\", r\"\\bquist\", r\"\\bzyste\\b\", r\"\\bcistic\", r\"reactive\", r\"reactivo\",\n)\n\nTRAUMA = _rx(\n r\"\\bbruise\\b\", r\"\\bcontusion\", r\"\\bkontuz\",\n r\"\\btrauma\", r\"\\bimpaction\\b\", r\"\\bpivot shift\\b\", r\"\\bkissing\\b\",\n r\"\\bacute\\b\", r\"\\bagudo\\b\", r\"\\bakut\", r\"\\bpivot kaymasi\\b\",\n r\"\\bbone bruise\\b\", r\"\\bbotcontusie\\b\",\n r\"\\bконтузион\", r\"\\bμωλωπ\", r\"\\bkontuzij\",\n)\n\n\ndef _polarity(clause: str) -> str:\n \"\"\"Classify one clause as positive, negative or uncertain for a matched term.\n\n Scope is the whole clause. Clause segmentation already keeps statements short, and\n a window in characters mis-scopes badly across languages with different word orders -\n Turkish puts its negator at the end of the sentence, English at the front.\n \"\"\"\n if UNCERTAIN.search(clause):\n return \"uncertain\"\n if NEGATION.search(clause):\n return \"negative\"\n if NORMALITY.search(clause):\n # \"meniscus normal\" negates; \"normal ... but tear\" does not.\n if TEAR.search(clause) or re.search(r\"\\bgrade [34]\\b\", clause):\n return \"positive\"\n return \"negative\"\n return \"positive\"\n\n\nclass _Matcher:\n \"\"\"Phrase lexicon first, stem+side proximity as the fallback.\n\n Exposes `.search` so it drops into the same slot as a compiled pattern.\n \"\"\"\n\n def __init__(self, phrase_rx, stem=None, side=None, window=55):\n self.phrase_rx = phrase_rx\n self.stem = stem\n self.side = side\n self.window = window\n\n def search(self, clause):\n m = self.phrase_rx.search(clause)\n if m is not None:\n return m\n if self.stem is not None and _near(clause, self.stem, self.side, self.window):\n return self.stem.search(clause)\n return None\n\n\nANAT_MATCH = {tgt: _Matcher(ANAT[tgt], *STEM_RULES[tgt]) for tgt in PAIRED}\nCOMPARTMENT_MATCH = {\n \"Medial OA\": _Matcher(COMPARTMENT[\"Medial OA\"], *STEM_RULES[\"Medial OA\"]),\n \"Lateral OA\": _Matcher(COMPARTMENT[\"Lateral OA\"], *STEM_RULES[\"Lateral OA\"]),\n \"PF OA\": _Matcher(COMPARTMENT[\"PF OA\"]),\n}\nDIRECT_MATCH = {\n tgt: _Matcher(_rx(rx.pattern, STEM_FRACTURE.pattern) if tgt == \"Fracture\" else rx)\n for tgt, rx in DIRECT.items()\n}\n\n\ndef _severity(clause: str) -> float:\n \"\"\"Weight one positive mention by how emphatic the sentence is.\n\n Ordered, not calibrated. A \"moderate effusion\" must outrank a \"trace effusion\" and\n both must outrank silence; the absolute numbers do not matter to AUC.\n \"\"\"\n high = SEV_HIGH.search(clause) is not None\n low = SEV_LOW.search(clause) is not None\n if high and not low:\n return 1.0\n if low and not high:\n return 0.45\n return 0.75 # unqualified mention\n\n\ndef _score_clauses(cls, anat_rx, path_rx=None, decoy_rx=None, context_penalty=None,\n context_bonus=None):\n \"\"\"Accumulate graded evidence over clauses for one target.\n\n Returns (score, confidence, n_pos, n_neg). Positives are graded by severity and by\n optional context regexes; negatives only matter when nothing positive was found,\n because reports assert normality for every structure they check.\n \"\"\"\n n_pos = n_neg = n_unc = 0\n best = 0.0\n for c in cls:\n m = anat_rx.search(c)\n if not m:\n continue\n if decoy_rx is not None and decoy_rx.search(c):\n continue\n if path_rx is not None and not path_rx.search(c):\n if NORMALITY.search(c) and not NEGATION.search(c):\n n_neg += 1\n continue\n pol = _polarity(c)\n if pol == \"positive\":\n n_pos += 1\n w = _severity(c)\n if context_penalty is not None and context_penalty.search(c):\n w *= 0.45\n if context_bonus is not None and context_bonus.search(c):\n w = min(1.0, w * 1.35)\n best = max(best, w)\n elif pol == \"negative\":\n n_neg += 1\n else:\n n_unc += 1\n best = max(best, 0.30)\n\n if n_pos or n_unc:\n # 0.52 .. 0.95, ordered by the strongest single mention, nudged by repetition.\n score = min(0.95, 0.50 + 0.42 * best + 0.03 * min(n_pos, 3))\n conf = min(1.0, 0.55 + 0.15 * n_pos)\n elif n_neg:\n score = max(0.04, 0.20 - 0.04 * n_neg)\n conf = min(0.9, 0.45 + 0.12 * n_neg)\n else:\n score, conf = 0.28, 0.05 # silence sits above asserted-negative\n return score, conf, n_pos, n_neg\n\n\ndef extract(report: str) -> dict:\n \"\"\"Extract twelve (score, confidence) pairs from one report.\"\"\"\n cls = clauses(report)\n out = {}\n path_paired = _rx(TEAR.pattern, DEGEN.pattern, INJURY.pattern)\n\n for tgt in TARGETS:\n if tgt in PAIRED:\n s, c, npos, nneg = _score_clauses(cls, ANAT_MATCH[tgt], path_paired)\n elif tgt in OA_TARGETS:\n s, c, npos, nneg = _score_clauses(cls, COMPARTMENT_MATCH[tgt], OA_EVIDENCE)\n elif tgt == \"Contusion\":\n # Reactive subchondral oedema under a cartilage defect is osteoarthritis,\n # not a bruise. Explicit trauma wording pushes the other way.\n s, c, npos, nneg = _score_clauses(cls, DIRECT_MATCH[tgt], None, DECOY.get(tgt),\n context_penalty=DEGENERATIVE_MARROW,\n context_bonus=TRAUMA)\n else:\n s, c, npos, nneg = _score_clauses(cls, DIRECT_MATCH[tgt], None, DECOY.get(tgt))\n out[tgt] = s\n out[tgt + \"__conf\"] = c\n out[tgt + \"__npos\"] = npos\n out[tgt + \"__nneg\"] = nneg\n\n # --- cross-target corrections ------------------------------------------ #\n # A whole-joint osteoarthritis statement is evidence for every compartment that was\n # not separately assessed. Without this, \"incipient OA of all three compartments\"\n # scores zero on all three OA targets.\n g_hits = [c for c in cls if GLOBAL_OA.search(c) and _polarity(c) == \"positive\"]\n if g_hits:\n gscore = 0.50 + 0.42 * max(_severity(c) for c in g_hits)\n for tgt in OA_TARGETS:\n if out[tgt + \"__npos\"] == 0 and out[tgt + \"__nneg\"] == 0:\n out[tgt] = max(out[tgt], gscore * 0.92)\n out[tgt + \"__conf\"] = max(out[tgt + \"__conf\"], 0.4)\n\n # Synovitis is frequently visible on the images and absent from the text, so silence\n # is weak evidence of absence here in a way it is not for other findings. Effusion is\n # its most reliable textual proxy - the two share a mechanism - so a silent synovitis\n # inherits a fraction of the effusion evidence instead of falling to the floor.\n if out[\"Synovitis__npos\"] == 0 and out[\"Synovitis__nneg\"] == 0:\n out[\"Synovitis\"] = max(out[\"Synovitis\"], 0.28 + 0.45 * (out[\"Effusion\"] - 0.28))\n\n return out","id":"cell-05"},{"cell_type":"markdown","metadata":{},"source":"### The lexicon hole, and what to do about it\n\nThe dangerous failure of a rule extractor is silent: a rule that never fires does not\nraise, it emits a negative. A lexicon thick in English and thin in Greek does not look\nbroken — it looks like a corpus where Greek patients have fewer findings. And because\nlanguage tracks the reporting institution, which tracks the scanner and the population,\nthat gap is a bias aligned with a site rather than noise that averages out.\n\nThe fix here is to fit a character n-gram model **to the rule scores**, over all 4,407\nreports, out of fold. `char_wb` 3–5 needs no tokeniser, stemmer or stopword list, so\nGreek and Turkish are handled on the same footing as English; the model picks up the\nphrasings that co-occur with rule-positive reports and scores them without anyone having\nwritten them into the lexicon.\n\nThe precise claim matters, because the stronger one is false. The model cannot find\nfindings in a language the lexicon misses *entirely* — its only supervision is the rule\noutput, so with no variation there is nothing to learn. What it repairs is **partial**\ncoverage, which is the real situation. `tests/test_distillation.py` measures exactly\nthis: on a corpus where 64% of positive reports use phrasings outside the lexicon, the\nrules sit at chance (0.500) on those reports and the blend reaches 0.940.\n\nEach target is gated on how well the text model recovered the rules out of fold. Below\na correlation of 0.15 it contributes nothing rather than noise — the negative control in\nthat test file confirms the gate fires when coverage is zero.\n\nThe other trap is on the way back. Combining as ranks is right, but reading the combined\nrank back through the *sorted* rule scores re-quantises everything into the rules' four\natoms, and since most mass sits in one enormous \"silent\" atom, that throws away precisely\nthe ordering the text model was added to supply. Interpolating between the distinct\nlevels instead keeps the prevalence and the meaning of 0.5 while preserving the ordering.","id":"cell-06"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"#\n# This is the part the source notebooks do not have, and it exists to patch the one\n# failure the rule extractor cannot see from inside itself.\n#\n# `rsna-knee-baseline-v1` names the problem precisely: a rule that never fires emits a\n# negative rather than an error, so a lexicon thin in one language looks like a\n# population with fewer findings, and language tracks site. `rsna-knee-eda-to-2-5d`\n# fits a TF-IDF teacher, but on the 58 annotated studies only - far too few rows for a\n# 190k-feature model, which is why its reported OOF AUC sits near chance on several\n# targets.\n#\n# Fitting the text model to the *rule scores* instead uses all 4,407 reports. The model\n# then learns which character n-grams co-occur with rule-positive reports, in every\n# language present, and can score a Greek report the lexicon was silent on. Out-of-fold\n# and grouped on report text so it cannot memorise its own training signal.\n\ndef distill(reports, rule_scores, groups, n_splits=5, seed=2026, alpha=1.0,\n verbose=True):\n \"\"\"Fit a char+word TF-IDF ridge to the rule scores, out-of-fold.\n\n reports : list[str], length N\n rule_scores : (N, 12) float, the graded rule output\n groups : (N,) int, fold assignment (hash of report text - duplicates together)\n returns : (N, 12) float, OOF text predictions on the same 0-1 scale\n\n Ridge on the graded score rather than logistic on a binarisation, because §1 of the\n source notebook is right that only order is read and grading carries more signal\n than a threshold does.\n \"\"\"\n from sklearn.feature_extraction.text import TfidfVectorizer\n from sklearn.linear_model import Ridge\n from scipy import sparse\n\n texts = [normalize(r) for r in reports]\n rule_scores = np.asarray(rule_scores, dtype=np.float32)\n groups = np.asarray(groups)\n\n # char_wb 3-5 is what makes this multilingual: it needs no tokeniser, no stemmer and\n # no stopword list, so Greek and Turkish are handled on the same footing as English.\n char = TfidfVectorizer(analyzer=\"char_wb\", ngram_range=(3, 5), min_df=3,\n max_features=200_000, sublinear_tf=True)\n word = TfidfVectorizer(ngram_range=(1, 2), min_df=2, max_features=80_000,\n sublinear_tf=True)\n X = sparse.hstack([char.fit_transform(texts), word.fit_transform(texts)],\n format=\"csr\")\n\n oof = np.zeros_like(rule_scores)\n folds = sorted(set(groups.tolist()))\n for f in folds:\n va = np.flatnonzero(groups == f)\n tr = np.flatnonzero(groups != f)\n if len(tr) < 20 or len(va) == 0:\n oof[va] = rule_scores[va]\n continue\n m = Ridge(alpha=alpha, solver=\"sparse_cg\")\n m.fit(X[tr], rule_scores[tr])\n oof[va] = m.predict(X[va])\n\n oof = np.clip(oof, 0.02, 0.98)\n\n # How well the text model reproduces the rules out of fold, per target. This is the\n # gate, not a diagnostic: the text model earns its weight by disagreeing with the\n # rules where the rules are silent, but a model that cannot even recover them\n # out-of-fold has learned nothing and must not be blended in. Targets whose lexicon\n # is thin, or whose positives are too rare for the ridge to find, land here.\n corr = np.zeros(rule_scores.shape[1], np.float32)\n for j in range(rule_scores.shape[1]):\n if np.std(oof[:, j]) > 1e-8 and np.std(rule_scores[:, j]) > 1e-8:\n corr[j] = float(np.corrcoef(oof[:, j], rule_scores[:, j])[0, 1])\n if verbose:\n print(\"distilled text model, OOF corr with rules per target:\")\n for t, a in zip(TARGETS, corr):\n print(f\" {t:<18} {a:+.3f}{'' if a > 0.15 else ' (gated out)'}\")\n return oof, corr\n\n\ndef blend(rule_scores, text_scores, w_text=0.35):\n \"\"\"Combine rules and the distilled text model in rank space, on the rules' scale.\n\n Two steps, and the second is the one that is easy to get wrong.\n\n *Combine as ranks.* The two sources are on different scales - the rules emit a\n handful of discrete levels, the ridge a continuous spread - so averaging the values\n lets the rules' coarse quantisation dominate. Ranks are the common currency.\n\n *Map back onto the rules' scale, by interpolation.* The combined rank is uniform on\n [0, 1] by construction, and using it directly as a BCE target would tell the model\n that exactly half of all knees have a fracture. So it is read back through the rule\n scores - but through their *distinct levels*, interpolating between them, not\n through the sorted array.\n\n Snapping to the sorted array instead is the obvious version and it silently undoes\n the whole exercise. The rule score takes about four distinct values, so most of its\n mass sits in one enormous atom at \"silent\"; a rank that lands anywhere inside that\n atom snaps back to the same number, and every ordering the text model contributed\n within the silent reports - which is precisely the ordering it was added to supply -\n is quantised away. Interpolating between the levels keeps the prevalence and the\n meaning of 0.5 while letting the text model order the reports inside each level.\n \"\"\"\n import pandas as pd\n\n rule_scores = np.asarray(rule_scores, dtype=np.float32)\n r = pd.DataFrame(rule_scores).rank(pct=True, method=\"average\").values\n t = pd.DataFrame(text_scores).rank(pct=True, method=\"average\").values\n\n # `w_text` may be a scalar or one weight per target - see `distill`, which returns\n # the per-target OOF correlation used to gate it.\n w = np.broadcast_to(np.asarray(w_text, np.float64).ravel(),\n (rule_scores.shape[1],)) if np.ndim(w_text) else \\\n np.full(rule_scores.shape[1], float(w_text))\n combined = (1.0 - w) * r + w * t\n\n out = np.empty_like(rule_scores)\n n = len(rule_scores)\n for j in range(rule_scores.shape[1]):\n levels, counts = np.unique(rule_scores[:, j], return_counts=True)\n if len(levels) == 1:\n out[:, j] = levels[0]\n continue\n # Anchor each distinct level at the midpoint of the rank interval it occupies,\n # then interpolate. Anchoring at the midpoint rather than an edge keeps the map\n # centred on the level, so a report the rules scored at that level and the text\n # model had no opinion about comes back out where it went in.\n mid = (np.cumsum(counts) - counts / 2.0) / n\n out[:, j] = np.interp(combined[:, j], mid, levels)\n return out","id":"cell-07"},{"cell_type":"markdown","metadata":{},"source":"## 3. Reading the acquisition\n\n`train_series.csv` gives each series a plane and two flags, `Fluid_Sensitive` and\n`Fat_Suppression`. The EDA notebook found those two columns are **byte-identical across\nall 24,371 rows**.\n\nThey name two physically independent things. *Fluid sensitivity* is a property of the\ncontrast weighting, set by $T_R$ and $T_E$:\n\n$$\\text{weighting} = \\begin{cases}\nT_1 & T_R \\lesssim 800\\,\\text{ms}\\\\\nT_2 & T_R \\gtrsim 800\\,\\text{ms},\\ T_E \\gtrsim 60\\,\\text{ms}\\\\\n\\text{PD} & T_R \\gtrsim 800\\,\\text{ms},\\ T_E \\lesssim 60\\,\\text{ms}\n\\end{cases}$$\n\n*Fat suppression* is a preparation applied on top of any weighting — a chemically\nselective pulse, STIR, or water excitation — and it is what makes marrow oedema\nconspicuous. Two identical columns carry one bit between them, not two, so recovering\nboth axes from the DICOM header is the difference between six slots and three.\n\nTwo string-matching traps, both of which silently invert the answer:\n\n- underscore is a word character, so a token test for `we` (water excitation) never\n fires inside `t2_de3d_we_tra`. Separators must be normalised first.\n- `ScanOptions` must be matched as **exact tokens**: one vendor writes `SAT_GEMS` for\n *spatial* saturation, so a substring test for `SAT` marks non-suppressed series as\n suppressed.\n\n| slot | plane | weighting | fat sat | what it carries |\n|---|---|---|---|---|\n| `SAG_FLUID_FS` | sagittal | PD / T2 | yes | meniscal tears, marrow oedema, effusion |\n| `COR_FLUID_FS` | coronal | PD / T2 | yes | collateral ligaments, meniscal body |\n| `AX_FLUID_FS` | axial | PD / T2 | yes | patellofemoral joint, synovium |\n| `SAG_FLUID_NOFS` | sagittal | PD / T2 | no | meniscal morphology at high CNR |\n| `COR_T1` | coronal | T1 | no | marrow architecture, bone outline |\n| `SAG_T1` | sagittal | T1 | no | anatomy, chronic change |","id":"cell-08"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"\"\"\"RSNA knee: study-level 12-label pipeline.\n\nSynthesised from two community notebooks, with the parts that were wrong or missing in\nboth replaced. What is inherited, and from where:\n\n from `rsna-knee-baseline-v1` header-recovered slot scheme, constant-physical-scale\n sampling, laterality normalisation, uint8 cache,\n per-diagnosis slot attention, report-hash grouping\n from `rsna-knee-eda-to-2-5d` protocol-only features as a legal test-time signal\n\nWhat is new here is listed in `README.md`; the five that change results are grouped\nK-fold instead of a single holdout, gold studies held out of their own fold so the\nannotated reference is honest, anatomy-preserving augmentation, a backbone that refuses\nto train from random initialisation silently, and rank fusion of the folds and the\nprotocol model.\n\"\"\"\n\nfrom __future__ import annotations\n\nimport os\n\nfor _v in (\"OMP_NUM_THREADS\", \"OPENBLAS_NUM_THREADS\", \"MKL_NUM_THREADS\"):\n os.environ.setdefault(_v, \"4\")\n\nimport gc\nimport hashlib\nimport re\nimport time\nimport warnings\nfrom concurrent.futures import ThreadPoolExecutor\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\n\nT0 = time.time()\n\n\ndef log(msg):\n print(f\"[{time.time() - T0:7.1f}s] {msg}\", flush=True)","id":"cell-09"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"class CFG:\n seed = 2026 # folds, target construction, text distillation\n\n # Separate from `seed` on purpose. This one drives only the stochastic parts of\n # image-model training: head init, batch order, group choice, augmentation. Bump\n # it alone and the folds and the derived targets are byte-identical, so the change\n # in the OOF score is training noise and nothing else. That number is the gate an\n # experiment has to clear before its result means anything.\n train_seed = int(os.environ.get(\"RSNA_TRAIN_SEED\", \"2026\"))\n\n img = 224 # encoder input; a multiple of the DINOv2 patch size 14\n crop_mm = 160.0 # physical extent of the centre crop; a knee FOV is 140-180\n\n # The crop guard `want < min(h, w)` reduces exactly to `FOV > crop_mm`, so a series\n # whose field of view is at or below crop_mm was silently left unnormalised - and\n # the measured median FOV in this corpus is 160 mm, right at the threshold. Padding\n # to crop_mm fixes the scale for those. Off by default so it can be A/B'd against\n # the noise floor rather than smuggled in.\n pad_short_fov = os.environ.get(\"RSNA_PAD_SHORT_FOV\", \"0\") == \"1\"\n\n # Whether to accept the geometric laterality fallback (sign of ImagePositionPatient\n # x). Off by default: the tag, ImageLaterality and text sources are all explicit\n # statements by the scanner, while this one is an inference. The run logs how often\n # it agrees with the tag where both exist, so it can be switched on with evidence\n # rather than on faith.\n lat_from_geometry = os.environ.get(\"RSNA_LAT_GEOMETRY\", \"0\") == \"1\"\n group = 3 # slices per encoder input, stacked as the three channels\n n_group_max = 3 # groups per slot, before the memory budget is applied\n cache_budget_gb = 12.0\n ram_fraction = 0.45 # ceiling on the cache as a share of free RAM\n\n hdr_threads = 16\n pix_threads = 12\n\n n_folds = 5\n max_folds_to_run = int(os.environ.get(\"RSNA_FOLDS\", \"5\"))\n epochs = 12\n batch_studies = 8\n lr_backbone = 8e-6 # the encoder is being adapted, not retrained\n lr_head = 1e-3\n weight_decay = 0.02\n unfreeze_last = 6 # trainable transformer blocks, from the output end\n eval_batch = 12\n time_budget = float(os.environ.get(\"RSNA_TIME_BUDGET\", 8.0 * 3600))\n\n w_text = 0.35 # weight of the distilled text model in the label blend\n w_protocol = 0.10 # weight of the protocol-only model in the final fusion\n gold_weight = 3.0\n\n\n# Six slots: three planes crossed with the acquisition axes. The fat-suppressed\n# fluid-sensitive series exist for nearly every study; the T1 and the non-suppressed\n# fluid-sensitive series are scarcer, which is what the presence mask is for.\nSLOTS_RECOVERED = [\n (\"SAG_FLUID_FS\", \"Sagittal\", True, True),\n (\"COR_FLUID_FS\", \"Coronal\", True, True),\n (\"AX_FLUID_FS\", \"Axial\", True, True),\n (\"SAG_FLUID_NOFS\", \"Sagittal\", True, False),\n (\"COR_T1\", \"Coronal\", False, False),\n (\"SAG_T1\", \"Sagittal\", False, False),\n]\n\n# The alternative: plane x the provided flag, ignoring the recovered weighting.\n#\n# Worth stating why the recovered scheme is the default. The EDA notebook found that\n# `Fluid_Sensitive` and `Fat_Suppression` are byte-identical across all 24,371 series\n# rows. They name two physically independent properties - contrast weighting, set by\n# TR/TE, and a fat-suppression preparation applied on top of any weighting - so two\n# identical columns carry one bit between them, not two. Recovering both axes from the\n# DICOM header is therefore not a refinement; it is the difference between six slots\n# and three.\nSLOTS_PUBLIC = [\n (\"SAG_FLUID\", \"Sagittal\", None, True),\n (\"COR_FLUID\", \"Coronal\", None, True),\n (\"AX_FLUID\", \"Axial\", None, True),\n (\"SAG_STRUCT\", \"Sagittal\", None, False),\n (\"COR_STRUCT\", \"Coronal\", None, False),\n (\"AX_STRUCT\", \"Axial\", None, False),\n]\n\nSLOT_SCHEME = os.environ.get(\"SLOT_SCHEME\", \"recovered\")\nSLOTS = SLOTS_PUBLIC if SLOT_SCHEME == \"public\" else SLOTS_RECOVERED\nN_SLOT = len(SLOTS)\n\nFATSAT_OPTS = {\"FS\", \"FATSAT\", \"FAT_SAT\", \"FSAT\"}\n_SEP = re.compile(r\"[_\\-.]\")\n_FATSAT_RX = re.compile(r\"\\bfs\\b|fatsat|fat sat|\\bstir\\b|\\bspair\\b|\\bspir\\b|\\bwe\\b|\"\n r\"water excit|\\btirm\\b|\\bsting\\b|\\bfatsup\\b\")\n_T1_RX = re.compile(r\"\\bt1\\b|\\bt1w\\b\")\n_T2_RX = re.compile(r\"\\bt2\\b|\\bt2w\\b\")\n_PD_RX = re.compile(r\"\\bpd\\b|\\bpdw\\b|proton|\\bdp\\b|dens\")\n\n\ndef seed_all(seed=CFG.seed):\n np.random.seed(seed)\n torch.manual_seed(seed)\n torch.cuda.manual_seed_all(seed)","id":"cell-10"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def find_root(explicit=None):\n if explicit is not None:\n p = Path(explicit)\n if (p / \"test.csv\").is_file():\n return p\n raise FileNotFoundError(f\"no test.csv under {p}\")\n for c in [Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\"),\n Path(\"/kaggle/input/rsna-knee-abnormality-detection\"),\n Path(\"data\"), Path(\".\")]:\n if (c / \"test.csv\").is_file() and (c / \"test_series\").is_dir():\n return c\n base = Path(\"/kaggle/input\")\n if base.is_dir():\n # last resort: two-level scan, because the mount is sometimes nested one deeper\n for depth1 in sorted(p for p in base.iterdir() if p.is_dir()):\n for cand in [depth1] + sorted(p for p in depth1.iterdir() if p.is_dir()):\n if (cand / \"test.csv\").is_file():\n return cand\n raise FileNotFoundError(\"competition mount not found\")\n\n\ndef find_dinov2(variant=\"small\"):\n \"\"\"Locate a mounted DINOv2 checkpoint directory by variant name.\"\"\"\n base = Path(\"/kaggle/input\")\n if not base.is_dir():\n return None\n hits = []\n for root, dirs, files in os.walk(base):\n dirs[:] = [d for d in dirs if d not in (\"train_series\", \"test_series\")]\n if \"config.json\" in files and \"dinov2\" in root.lower():\n hits.append(Path(root))\n for h in hits:\n if variant in str(h).lower():\n return h\n return hits[0] if hits else None","id":"cell-11"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"HDR_TAGS = [\"SeriesDescription\", \"SequenceName\", \"ScanOptions\", \"ScanningSequence\",\n \"RepetitionTime\", \"EchoTime\", \"Laterality\", \"PixelSpacing\", \"Rows\",\n \"Columns\", \"RescaleSlope\", \"RescaleIntercept\",\n # Laterality fallbacks. `Laterality` is absent on about half the series,\n # and a study with no side at all is never mirrored - so medial and lateral\n # sit on opposite sides of the image for left and right knees. Four of the\n # twelve targets are side-specific, so that is a third of the metric.\n \"ImageLaterality\", \"StudyDescription\", \"BodyPartExamined\",\n \"PatientPosition\", \"ImagePositionPatient\"]\n\n\ndef probe(item):\n split, study, series, path = item\n row = {\"split\": split, \"StudyInstanceUID\": study, \"SeriesInstanceUID\": series,\n \"dir\": path}\n try:\n files = sorted(e.name for e in os.scandir(path) if e.name.endswith(\".dcm\"))\n row[\"files\"] = files\n row[\"n_slices\"] = len(files)\n if not files:\n return row\n ds = pydicom.dcmread(os.path.join(path, files[len(files) // 2]),\n stop_before_pixels=True, force=True)\n for t in HDR_TAGS:\n v = getattr(ds, t, None)\n if v is None:\n row[t] = None\n elif isinstance(v, (list, tuple)) or type(v).__name__ == \"MultiValue\":\n row[t] = \"|\".join(str(x) for x in v)\n else:\n row[t] = str(v)\n except Exception as exc:\n row[\"err\"] = str(exc)[:120]\n return row\n\n\ndef walk(root, split):\n base = Path(root) / split\n items = []\n if not base.is_dir():\n return pd.DataFrame()\n for study in os.scandir(base):\n if study.is_dir():\n for series in os.scandir(study.path):\n if series.is_dir():\n items.append((split, study.name, series.name, series.path))\n with ThreadPoolExecutor(max_workers=CFG.hdr_threads) as pool:\n rows = list(pool.map(probe, items))\n return pd.DataFrame(rows)\n\n\ndef annotate(df):\n \"\"\"Recover fat suppression and pulse-sequence weighting from the header.\n\n Two string-matching cautions, both of which silently invert the answer if missed:\n underscore is a word character, so a token test for `we` (water excitation) never\n fires inside `t2_de3d_we_tra` unless separators are normalised first; and GE writes\n `SAT_GEMS` for *spatial* saturation, so ScanOptions must be matched as exact tokens\n or non-fat-suppressed series get marked as suppressed.\n \"\"\"\n if df.empty:\n return df\n for t in HDR_TAGS:\n if t not in df.columns:\n df[t] = None\n\n desc = (df[\"SeriesDescription\"].fillna(\"\") + \" \" + df[\"SequenceName\"].fillna(\"\"))\n desc = desc.str.lower().str.replace(_SEP, \" \", regex=True)\n\n opts = df[\"ScanOptions\"].fillna(\"\").str.upper().str.split(\"|\")\n opts_fs = opts.apply(lambda ts: any(t.strip() in FATSAT_OPTS for t in ts))\n df[\"fatsat\"] = desc.str.contains(_FATSAT_RX) | opts_fs\n\n tr = pd.to_numeric(df[\"RepetitionTime\"], errors=\"coerce\")\n te = pd.to_numeric(df[\"EchoTime\"], errors=\"coerce\")\n gre = df[\"ScanningSequence\"].fillna(\"\").str.upper().str.contains(\"GR\")\n t1 = desc.str.contains(_T1_RX)\n t2 = desc.str.contains(_T2_RX)\n pdw = desc.str.contains(_PD_RX)\n\n df[\"weight\"] = np.where(t1 & ~t2 & ~pdw, \"T1\",\n np.where(t2 & ~pdw, \"T2\",\n np.where(pdw, \"PD\",\n np.where(gre, \"GRE\",\n np.where(tr < 800, \"T1\",\n np.where(te > 60, \"T2\",\n np.where(tr >= 800, \"PD\", \"UNK\")))))))\n df[\"fluid\"] = np.isin(df[\"weight\"], [\"PD\", \"T2\"])\n df[\"px\"] = pd.to_numeric(\n df[\"PixelSpacing\"].fillna(\"\").str.split(\"|\").str[0].replace(\"\", np.nan),\n errors=\"coerce\")\n return df\n\n\ndef pick_slots(series_df, plane_map):\n \"\"\"One series per slot per study.\n\n Ties are broken toward the stack with the most slices: a thicker stack samples the\n joint more densely, and the three-slice sampler benefits from the margin.\n \"\"\"\n if series_df.empty:\n return {}\n series_df = series_df.copy()\n series_df[\"plane\"] = series_df[\"SeriesInstanceUID\"].map(plane_map)\n out = {}\n for study, g in series_df.groupby(\"StudyInstanceUID\"):\n chosen = {}\n for name, plane, fluid, fs in SLOTS:\n sel = (g[\"plane\"] == plane) & (g[\"fatsat\"] == fs)\n # fluid=None means \"do not condition on weighting\" - the public scheme,\n # where the single provided flag stands in for both axes at once.\n if fluid is not None:\n sel = sel & (g[\"fluid\"] == fluid)\n cand = g[sel]\n if len(cand) == 0 and fluid is False:\n # T1 slots are the scarcest; fall back to any non-fat-sat series in the\n # plane before giving up on the slot entirely.\n cand = g[(g[\"plane\"] == plane) & (~g[\"fatsat\"])]\n if len(cand):\n chosen[name] = cand.sort_values(\"n_slices\", ascending=False).iloc[0]\n out[study] = chosen\n return out","id":"cell-12"},{"cell_type":"markdown","metadata":{},"source":"## 4. Sampling at a fixed physical scale\n\nA slice of $N\\times N$ pixels at spacing $s$ mm/pixel covers $Ns$ millimetres. Both vary\nwidely here, so a fixed-pixel resize hands the encoder images whose physical scale\ndiffers by about $3\\times$, and a meniscus occupies a different number of pixels in\ndifferent studies for no anatomical reason.\n\nCropping to a constant physical extent $L$ first and resampling to $P\\times P$ leaves\n\n$$s_\\text{eff} = L/P \\ \\text{mm/pixel}$$\n\nindependent of the acquisition — $160\\,\\text{mm} / 224 = 0.71$ mm/pixel here.\n\nIntensity needs the same treatment: MR has no absolute scale, so each series is\nnormalised to its own 1st–99th percentile over the whole volume (not per slice, so\nslices keep their relative contrast), which removes an offset that would otherwise track\nthe site.\n\n### Laterality\n\nFour targets are medial/lateral pairs, and medial is defined against the body midline —\nso which side of the *image* it falls on depends on which knee was scanned. The\ncorrection differs by plane: coronal and axial mirror under a horizontal flip, but a\nsagittal stack does not, because there the medial–lateral direction is the *slice* axis\nand what differs is the order the stack traverses the joint. Where `Laterality` is\nabsent the volume is left alone — a wrong flip is worse than no flip, and the presence\nmask lets the head learn how much to trust each slot.","id":"cell-13"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"ORDER_TAGS = [\"ImagePositionPatient\", \"ImageOrientationPatient\", \"SliceLocation\",\n \"InstanceNumber\"]\n\n# Ordering happens on reader threads, so `d[k] += 1` - three bytecodes, not one - needs\n# the lock. An undercounted fallback reads as a clean run.\nimport threading\n\nORDER_STATS = {\"geometry\": 0, \"slice_location\": 0, \"instance_number\": 0, \"filename\": 0}\n_ORDER_LOCK = threading.Lock()\n\n# The pixel path had no instrumentation at all, so three separate silent failures were\n# invisible: a physical crop that skips itself, a slice that fails to decode and is\n# substituted with black, and an inverted photometric interpretation. Each one degrades\n# an image without raising anything. Counted here, they become measurements.\nPIX_STATS = {\"crop_applied\": 0, \"crop_padded\": 0, \"crop_skipped_short_fov\": 0,\n \"crop_no_spacing\": 0, \"decode_failed\": 0, \"mono1_inverted\": 0,\n \"slices_read\": 0}\n_PIX_LOCK = threading.Lock()\n\n\ndef _bump(key):\n with _ORDER_LOCK:\n ORDER_STATS[key] += 1\n\n\ndef _bump_pix(key, n=1):\n with _PIX_LOCK:\n PIX_STATS[key] += n\n\n\ndef order_slices(directory, files):\n \"\"\"Return `files` sorted by position along the stack, nearest-first.\n\n The filenames in this corpus are SOP Instance UIDs, which are random, so the sorted\n filename order is a random permutation of the anatomy. Measured on the competition\n data, filename order agrees with physical order on about 5% of slices.\n\n That is not a cosmetic problem. Two things downstream assume the list is ordered:\n the sampling window (\"the central portion of the stack\") and the grouping of\n adjacent slices into the encoder's three channels. Under a random permutation the\n window selects a uniform random sample of the whole series rather than its centre,\n and a group of three adjacent entries is three unrelated positions in the knee - so\n the 2.5D input carried no depth information at all.\n\n Fallback chain, because not every series records the same tags:\n 1. project ImagePositionPatient onto the normal of ImageOrientationPatient. This\n is the only one that is correct for an obliquely angled stack.\n 2. SliceLocation, which the scanner has already projected but which is signed\n inconsistently between vendors.\n 3. InstanceNumber, acquisition order - right for most stacks, wrong for\n interleaved acquisitions.\n 4. filename, i.e. give up, and say so in the counters.\n \"\"\"\n if len(files) < 2:\n return list(files)\n\n positions, locations, instances = [], [], []\n normal = None\n for name in files:\n try:\n ds = pydicom.dcmread(os.path.join(directory, name), stop_before_pixels=True,\n force=True, specific_tags=ORDER_TAGS)\n except Exception:\n positions.append(None)\n locations.append(None)\n instances.append(None)\n continue\n pos = getattr(ds, \"ImagePositionPatient\", None)\n orient = getattr(ds, \"ImageOrientationPatient\", None)\n if normal is None and orient is not None and len(orient) == 6:\n try:\n row_dir = np.array([float(v) for v in orient[:3]])\n col_dir = np.array([float(v) for v in orient[3:]])\n normal = np.cross(row_dir, col_dir)\n except Exception:\n normal = None\n try:\n positions.append(np.array([float(v) for v in pos]) if pos is not None else None)\n except Exception:\n positions.append(None)\n loc = getattr(ds, \"SliceLocation\", None)\n locations.append(float(loc) if loc is not None else None)\n num = getattr(ds, \"InstanceNumber\", None)\n instances.append(float(num) if num is not None else None)\n\n def sorted_by(values, stat):\n _bump(stat)\n return [f for _, f in sorted(zip(values, files), key=lambda pair: pair[0])]\n\n if normal is not None and all(p is not None for p in positions):\n return sorted_by([float(p @ normal) for p in positions], \"geometry\")\n if all(loc is not None for loc in locations):\n return sorted_by(locations, \"slice_location\")\n if all(num is not None for num in instances):\n return sorted_by(instances, \"instance_number\")\n _bump(\"filename\")\n return list(files)\n\n\n# Fraction of the ordered stack each plane is sampled across.\n#\n# These are wider than the 20-80% window this replaced, and the widening is part of the\n# ordering fix rather than a separate change. Before ordering, \"the central 60%\" of a\n# random permutation was a uniform sample of the *whole* series; applying that same\n# window to a correctly ordered stack would genuinely discard the outer 40% and quietly\n# bundle a coverage reduction into this change. These keep the physical span roughly as\n# it was, so the ordering is the only variable that moved.\n#\n# The per-plane split follows Will's (`wguesdon`) reasoning: the menisci sit at the\n# medial and lateral extremes of a sagittal stack and Baker's cysts sit posteriorly on\n# an axial one, so those planes need their edges; a coronal stack's useful\n# anterior-posterior range is narrower.\nPLANE_WINDOW = {\"Sagittal\": (0.10, 0.90), \"Axial\": (0.10, 0.90),\n \"Coronal\": (0.15, 0.85)}\n\n\ndef read_slot(rec, n_slice, out_size, plane=None, group=None):\n \"\"\"`n_slice` slices from one series, physically ordered, at `out_size` pixels.\n\n Returns uint8 [n_slice, out, out], normalised per-series to its 1st-99th percentile.\n Percentiles rather than min/max because MR intensity has no absolute scale and one\n bright vessel would otherwise compress the whole dynamic range.\n\n The slices come out as `n_slice // group` anchors spread across the sampling window,\n each anchor contributing `group` *physically adjacent* slices. That layout is what\n `take_group` slices back out into the encoder's three channels, so a channel triplet\n is a genuine depth neighbourhood - roughly 10 mm of knee at this corpus's slice gaps\n - rather than three arbitrary positions. Spreading all nine slices evenly instead\n would cover more of the joint but hand the encoder three views 20 mm apart, which is\n a slab, not a 2.5D triplet.\n\n A slice of N pixels at spacing s mm covers Ns millimetres. Both vary widely here, so\n a fixed-pixel resize hands the encoder images whose physical scale differs by about\n 3x. Cropping to a constant physical extent first fixes the effective scale at\n crop_mm / img mm per pixel regardless of acquisition.\n \"\"\"\n files, d, px = rec[\"files\"], rec[\"dir\"], rec[\"px\"]\n n = len(files)\n if n == 0:\n return None\n group = CFG.group if group is None else group\n\n files = order_slices(d, files)\n\n lo_f, hi_f = PLANE_WINDOW.get(plane, (0.10, 0.90))\n lo, hi = int(lo_f * (n - 1)), int(hi_f * (n - 1))\n if hi <= lo:\n lo, hi = 0, n - 1\n\n n_anchor = max(1, n_slice // group)\n anchors = (np.linspace(lo, hi, n_anchor).astype(int) if n_anchor > 1\n else np.array([(lo + hi) // 2]))\n\n idx = []\n for centre in anchors:\n # Adjacent slices around the anchor, clipped into the series rather than\n # wrapped: a wrap would put the far end of the knee in the same channel stack.\n start = int(np.clip(centre - group // 2, 0, max(0, n - group)))\n idx.extend(range(start, min(start + group, n)))\n while len(idx) < n_slice:\n idx.append(idx[-1])\n\n planes = []\n for i in idx[:n_slice]:\n try:\n ds = pydicom.dcmread(os.path.join(d, files[int(i)]), force=True)\n a = ds.pixel_array.astype(np.float32)\n sl = float(getattr(ds, \"RescaleSlope\", 1) or 1)\n ic = float(getattr(ds, \"RescaleIntercept\", 0) or 0)\n a = a * sl + ic\n # MONOCHROME1 means high value renders as black - the image is inverted\n # relative to MONOCHROME2. Left uncorrected, fluid reads dark in those\n # series while reading bright everywhere else, and Effusion and Synovitis\n # are fluid-bright findings. Invert into the MONOCHROME2 convention.\n if str(getattr(ds, \"PhotometricInterpretation\", \"\")).strip() == \"MONOCHROME1\":\n a = a.max() - a\n _bump_pix(\"mono1_inverted\")\n _bump_pix(\"slices_read\")\n except Exception:\n a = None\n _bump_pix(\"decode_failed\")\n planes.append(a)\n\n shp = next((p.shape for p in planes if p is not None), None)\n if shp is None:\n return None\n planes = [p if (p is not None and p.shape == shp) else np.zeros(shp, np.float32)\n for p in planes]\n vol = np.stack(planes)\n\n if px and np.isfinite(px) and px > 0:\n want = int(round(CFG.crop_mm / px))\n h, w = shp\n if 16 < want < min(h, w):\n cy, cx = h // 2, w // 2\n half = want // 2\n vol = vol[:, max(0, cy - half):cy + half, max(0, cx - half):cx + half]\n _bump_pix(\"crop_applied\")\n elif want >= min(h, w):\n # The series covers less than crop_mm, so there is nothing to crop away.\n # Resizing here would set mm/px from the acquisition instead of from\n # crop_mm, which is the one thing the constant crop exists to prevent.\n # Padding to the same physical extent keeps the scale fixed; the border is\n # black because no anatomy was acquired there.\n if CFG.pad_short_fov:\n py, pxd = max(0, want - h), max(0, want - w)\n vol = np.pad(vol, ((0, 0), (py // 2, py - py // 2),\n (pxd // 2, pxd - pxd // 2)))\n _bump_pix(\"crop_padded\")\n else:\n _bump_pix(\"crop_skipped_short_fov\")\n else:\n _bump_pix(\"crop_skipped_short_fov\")\n else:\n _bump_pix(\"crop_no_spacing\")\n\n lo_v, hi_v = np.percentile(vol, [1, 99])\n vol = np.clip((vol - lo_v) / max(hi_v - lo_v, 1e-6), 0, 1)\n\n t = torch.from_numpy(np.ascontiguousarray(vol)).unsqueeze(0)\n t = F.interpolate(t, size=(out_size, out_size), mode=\"bilinear\", align_corners=False)\n # uint8, not float32: intensity is already normalised into [0, 1] here, so eight bits\n # cost nothing a bilinear resize has not already cost, and the cache is a quarter the\n # size for it.\n return (t.squeeze(0) * 255).round().clamp(0, 255).to(torch.uint8)\n\n\n# Word-boundary side markers. Deliberately NOT including bare Turkish \"sag\" (right):\n# every sagittal protocol in this corpus is named `sag_...`, so that token would mark\n# most of the dataset as right knees. Same family of trap as GE's SAT_GEMS. The\n# diacritic form is unambiguous and is matched; the stripped form is not.\n_SIDE_RX = [\n (\"R\", re.compile(r\"\\b(right|rt|r_?knee|knee_?r|dexter|sağ|derech[ao]|rechts?|\"\n r\"droite?|δεξ\\w*)\\b\", re.I)),\n (\"L\", re.compile(r\"\\b(left|lt|l_?knee|knee_?l|sinister|sol|izquierd[ao]|links?|\"\n r\"gauche|αριστερ\\w*)\\b\", re.I)),\n]\n\nLAT_STATS = {\"tag\": 0, \"image_laterality\": 0, \"text\": 0, \"geometry\": 0, \"unknown\": 0,\n \"geom_agree\": 0, \"geom_disagree\": 0}\n\n\ndef _lat_from_text(*fields):\n \"\"\"Side from any free-text header field, or None if absent or contradictory.\"\"\"\n blob = \" \".join(str(f) for f in fields if f and str(f).lower() != \"nan\")\n blob = re.sub(r\"[^\\w\\s]\", \" \", blob)\n hits = {side for side, rx in _SIDE_RX if rx.search(blob)}\n return hits.pop() if len(hits) == 1 else None\n\n\ndef _lat_from_geometry(ipp):\n \"\"\"Side from the x coordinate of ImagePositionPatient.\n\n DICOM's patient coordinate system is LPS: +x runs toward the patient's LEFT. A knee\n is therefore centred at negative x when it is the right knee. This is independent of\n head-first/feet-first, because ImagePositionPatient is already expressed in patient\n coordinates rather than scanner coordinates.\n\n The threshold exists because a value near zero means the scan is near the midline,\n where the sign carries no information.\n \"\"\"\n try:\n x = float(str(ipp).split(\"|\")[0])\n except (TypeError, ValueError):\n return None\n if abs(x) < 20.0:\n return None\n return \"R\" if x < 0 else \"L\"\n\n\ndef resolve_laterality(g):\n \"\"\"Study -> 'L'/'R'/None, from the strongest available evidence.\n\n Order: the explicit tag, then ImageLaterality, then side words in any description\n field, then the scan geometry. Each source is counted so the log says how much of\n the corpus each one actually rescued, and the geometric rule is scored against the\n explicit tag wherever both exist - a fallback nobody has validated is a guess.\n \"\"\"\n vals = [str(x).strip().upper() for x in g.get(\"Laterality\", pd.Series(dtype=object))\n .dropna()]\n vals = [v[0] for v in vals if v and v[0] in (\"L\", \"R\")]\n\n geom = next((s for s in (_lat_from_geometry(v)\n for v in g.get(\"ImagePositionPatient\",\n pd.Series(dtype=object)).dropna()) if s), None)\n if vals and geom:\n _bump_lat(\"geom_agree\" if geom == vals[0] else \"geom_disagree\")\n if vals:\n _bump_lat(\"tag\")\n return vals[0]\n\n ivals = [str(x).strip().upper() for x in\n g.get(\"ImageLaterality\", pd.Series(dtype=object)).dropna()]\n ivals = [v[0] for v in ivals if v and v[0] in (\"L\", \"R\")]\n if ivals:\n _bump_lat(\"image_laterality\")\n return ivals[0]\n\n txt = _lat_from_text(*g.get(\"SeriesDescription\", pd.Series(dtype=object)).tolist(),\n *g.get(\"StudyDescription\", pd.Series(dtype=object)).tolist(),\n *g.get(\"BodyPartExamined\", pd.Series(dtype=object)).tolist())\n if txt:\n _bump_lat(\"text\")\n return txt\n\n if geom and CFG.lat_from_geometry:\n _bump_lat(\"geometry\")\n return geom\n\n _bump_lat(\"unknown\")\n return None\n\n\ndef _bump_lat(key):\n LAT_STATS[key] += 1\n\n\ndef normalise_laterality(img, plane, lat):\n \"\"\"Map every knee onto a left-knee convention.\n\n Four of the twelve targets are medial/lateral pairs, and medial is defined against\n the body midline, so which side of the *image* it falls on depends on which knee was\n scanned. Coronal and axial views mirror under a horizontal flip. Sagittal stacks do\n not - each slice is unchanged by mirroring, what differs is the direction the stack\n traverses the joint - so the slice order is reversed instead.\n\n Where `Laterality` is absent the volume is left alone: a wrong flip is worse than no\n flip, and the presence mask lets the head learn how much to trust each slot.\n \"\"\"\n if lat != \"R\":\n return img\n if plane in (\"Coronal\", \"Axial\"):\n return torch.flip(img, dims=[-1])\n return torch.flip(img, dims=[0])\n\n\ndef available_ram_gb():\n \"\"\"Physical RAM available to this process, in GB, or None if it cannot be read.\"\"\"\n try:\n import psutil\n return psutil.virtual_memory().available / 1024 ** 3\n except Exception:\n pass\n try: # Linux without psutil, which includes Kaggle\n with open(\"/proc/meminfo\") as fh:\n for line in fh:\n if line.startswith(\"MemAvailable:\"):\n return int(line.split()[1]) / 1024 ** 2\n except Exception:\n pass\n return None\n\n\ndef plan_cache(n_study):\n \"\"\"Choose how many slices per slot memory allows.\n\n The cache is n_study x n_slot x slices x img^2 bytes. It grows with the *square* of\n resolution and only linearly with slices, so coverage is the cheap axis and\n resolution the expensive one; when the budget binds it is the slice count that gives\n way. Deciding once, from the training corpus size, keeps train and test on the same\n group layout.\n\n The configured budget is a ceiling, not a promise. At full corpus size the cache is\n around 12 GB, which fits a 30 GB Kaggle instance and does not fit a 13 GB one - and\n the failure mode is a SIGKILL partway through the decode, which no `try` will catch\n and which costs the whole run. So the budget is also capped at a fraction of what\n the machine actually reports free, leaving room for the model, the reader queue and\n the test cache.\n \"\"\"\n budget = CFG.cache_budget_gb\n free = available_ram_gb()\n if free is not None:\n safe = CFG.ram_fraction * free\n if safe < budget:\n log(f\"cache budget trimmed {budget:.1f} -> {safe:.1f} GB \"\n f\"({free:.1f} GB free, keeping {1 - CFG.ram_fraction:.0%} headroom)\")\n budget = safe\n per_slice = n_study * N_SLOT * CFG.img * CFG.img\n afford = int(budget * 1024 ** 3 // max(per_slice, 1))\n groups = max(1, min(CFG.n_group_max, afford // CFG.group))\n if groups < CFG.n_group_max:\n log(f\"cache budget {budget:.1f} GB allows {groups} group(s) of \"\n f\"{CFG.group}, not {CFG.n_group_max}\")\n return groups\n\n\ndef build_cache(slot_map, plane_map, lat_map, tag, n_group):\n \"\"\"Decode every (study, slot) once into an in-memory uint8 array.\n\n Fine-tuning revisits the same pixels every epoch, and a study is on the order of a\n hundred and fifty files. Reading them from the mount each epoch would make the epoch\n count a function of I/O rather than of learning.\n\n Reads are issued in bounded chunks: submitting every job at once lets the reader\n threads run arbitrarily far ahead and the completed buffers accumulate without limit.\n \"\"\"\n cache_slices = CFG.group * n_group\n studies = sorted(slot_map)\n sidx = {s: i for i, s in enumerate(studies)}\n cache = np.zeros((len(studies), N_SLOT, cache_slices, CFG.img, CFG.img), np.uint8)\n mask = np.zeros((len(studies), N_SLOT), np.float32)\n log(f\"{tag}: cache {cache.shape} = {cache.nbytes / 1024 ** 3:.1f} GB\")\n\n jobs = [(st, k, plane, slot_map[st][name])\n for st in studies\n for k, (name, plane, _, _) in enumerate(SLOTS)\n if name in slot_map[st]]\n log(f\"{tag}: decoding {len(jobs)} slot-series\")\n\n # Both stat dicts are module-level, so without this the `test:` line reports\n # train+test cumulatively and its percentages describe neither pass.\n ORDER_STATS.update({k: 0 for k in ORDER_STATS})\n PIX_STATS.update({k: 0 for k in PIX_STATS})\n\n chunk = 512\n done = 0\n with ThreadPoolExecutor(max_workers=CFG.pix_threads) as pool:\n for c0 in range(0, len(jobs), chunk):\n block = jobs[c0:c0 + chunk]\n imgs = pool.map(lambda j: read_slot(j[3], cache_slices, CFG.img, j[2]), block)\n for (st, k, plane, _), img in zip(block, imgs):\n done += 1\n if img is None:\n continue\n cache[sidx[st], k] = normalise_laterality(img, plane,\n lat_map.get(st)).numpy()\n mask[sidx[st], k] = 1.0\n if done % 4096 < chunk:\n log(f\" {tag} {done}/{len(jobs)}\")\n if time.time() - T0 > CFG.time_budget:\n log(f\" {tag}: time budget reached during decode\")\n break\n total = sum(ORDER_STATS.values()) or 1\n log(f\"{tag}: slice ordering \" + \", \".join(\n f\"{k} {v / total:.1%}\" for k, v in ORDER_STATS.items() if v))\n if ORDER_STATS[\"filename\"] / total > 0.05:\n warnings.warn(\n f\"{ORDER_STATS['filename'] / total:.0%} of series fell back to filename \"\n \"order, which is SOP UID order and therefore random. Those slots carry no \"\n \"depth structure.\", RuntimeWarning, stacklevel=2)\n\n n_crop = sum(PIX_STATS[k] for k in\n (\"crop_applied\", \"crop_padded\", \"crop_skipped_short_fov\",\n \"crop_no_spacing\")) or 1\n log(f\"{tag}: physical scale \" + \", \".join(\n f\"{k.replace('crop_', '')} {PIX_STATS[k]} ({PIX_STATS[k] / n_crop:.1%})\"\n for k in (\"crop_applied\", \"crop_padded\", \"crop_skipped_short_fov\",\n \"crop_no_spacing\") if PIX_STATS[k]))\n log(f\"{tag}: slices read {PIX_STATS['slices_read']}, \"\n f\"decode failures {PIX_STATS['decode_failed']}, \"\n f\"MONOCHROME1 inverted {PIX_STATS['mono1_inverted']}\")\n skipped = PIX_STATS[\"crop_skipped_short_fov\"] / n_crop\n if skipped > 0.05:\n warnings.warn(\n f\"{skipped:.0%} of slots had a field of view at or below crop_mm=\"\n f\"{CFG.crop_mm:.0f}mm, so no physical normalisation was applied to them and \"\n \"their mm/px came from the acquisition. Set RSNA_PAD_SHORT_FOV=1 to pad \"\n \"instead.\", RuntimeWarning, stacklevel=2)\n gc.collect()\n return studies, cache, mask","id":"cell-14"},{"cell_type":"markdown","metadata":{},"source":"## 5. Aggregating slots into twelve decisions\n\nA study arrives as up to six slot embeddings $x_s$ with a presence mask $m_s$. Pooling\nthem identically would discard the reason the protocol has three planes at all. Give\nevery diagnosis $o$ its own query and let it attend over the slots, absent ones masked\nout of the softmax:\n\n$$h_s = \\phi(x_s) + e_s,\\qquad\n\\alpha_{o,s} = \\frac{\\exp(\\langle h_s, q_o\\rangle/\\sqrt{H})\\, m_s}\n{\\sum_{s'} \\exp(\\langle h_{s'}, q_o\\rangle/\\sqrt{H})\\, m_{s'}}$$\n\nThe masked softmax renormalises over whatever the study actually contains, so a missing\naxial series shifts a diagnosis's attention onto the sequences that are present rather\nthan feeding it a zero vector.\n\n**The head is deliberately this small.** The label is attached to the *study*, so\nnothing in the supervision says which part of it carries the finding. Extra attention\nparameters have no signal to learn that from and spend their capacity on noise. Where\nthe supervision is coarse, the aggregation should be too.\n\n**Why the encoder is trained rather than frozen.** A frozen self-supervised encoder is\nthe cheap option, and the way it stops paying is informative: resolution, encoder size,\nslice coverage and slot aggregation each move the score by less than validation noise.\nThat pattern is diagnostic — those four axes change how much the model looks and how it\nsummarises, but none changes the *vocabulary* it looks with. So the last blocks are\nadapted, at a learning rate two orders of magnitude below the head's, because the head\nis random at initialisation and the encoder starts from a good solution.","id":"cell-15"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"class SlotHead(nn.Module):\n \"\"\"Per-diagnosis attention over the slot embeddings of one study.\n\n Each finding is read on particular sequences - cruciates sagittally, collateral\n ligaments and the meniscal body coronally, patellar cartilage axially - so pooling\n the slots identically would dilute the one that carries the evidence with the rest.\n The softmax is masked, so a missing axial series shifts a diagnosis's attention onto\n the sequences that are present instead of feeding it a zero vector.\n\n The aggregation is deliberately this simple. The label is attached to the *study*, so\n nothing in the supervision says which part of a study carries the finding; extra\n attention parameters have no signal to learn that from and spend capacity on noise.\n Where the supervision is coarse, the aggregation should be too.\n \"\"\"\n\n def __init__(self, dim, n_slot, n_out, hidden=256, p=0.2):\n super().__init__()\n self.proj = nn.Sequential(nn.LayerNorm(dim), nn.Linear(dim, hidden), nn.GELU())\n self.slot_emb = nn.Parameter(torch.randn(n_slot, hidden) * 0.02)\n self.query = nn.Parameter(torch.randn(n_out, hidden) * 0.02)\n self.drop = nn.Dropout(p)\n self.out = nn.Linear(hidden, n_out)\n self.hidden = hidden\n\n def forward(self, x, mask):\n h = self.proj(x) + self.slot_emb\n att = torch.einsum(\"bsh,oh->bos\", h, self.query) / self.hidden ** 0.5\n att = att.masked_fill(mask.unsqueeze(1) < 0.5, -1e4).softmax(-1)\n ctx = self.drop(torch.einsum(\"bos,bsh->boh\", att, h))\n return (ctx * self.out.weight.unsqueeze(0)).sum(-1) + self.out.bias\n\n\nclass Encoder(nn.Module):\n \"\"\"Uniform [B,3,H,W] -> [B,dim] wrapper over either a HF or a timm backbone.\"\"\"\n\n def __init__(self, module, kind, dim):\n super().__init__()\n self.module = module\n self.kind = kind\n self.dim = dim\n\n def forward(self, x):\n if self.kind == \"hf\":\n out = self.module(pixel_values=x).last_hidden_state\n # CLS and mean-pooled patches: the CLS token carries the global summary and\n # the patch mean the spatially distributed evidence, and the findings here\n # need both.\n return torch.cat([out[:, 0], out[:, 1:].mean(1)], dim=1)\n return self.module(x)\n\n def blocks(self):\n if self.kind == \"hf\":\n return list(self.module.encoder.layer)\n for attr in (\"blocks\", \"layers\"):\n b = getattr(self.module, attr, None)\n if b is not None:\n return list(b)\n return []\n\n\nclass Model(nn.Module):\n \"\"\"Encoder plus head, trained end to end.\n\n A study arrives as a bag of slot images. The bag is flattened for the encoder and\n folded back before the head, so the encoder never sees the study structure and the\n head never sees pixels.\n \"\"\"\n\n def __init__(self, encoder):\n super().__init__()\n self.encoder = encoder\n self.head = SlotHead(encoder.dim, N_SLOT, len(TARGETS))\n self.register_buffer(\"mean\", torch.tensor([0.485, 0.456, 0.406]).view(1, 3, 1, 1))\n self.register_buffer(\"std\", torch.tensor([0.229, 0.224, 0.225]).view(1, 3, 1, 1))\n\n def forward(self, imgs, mask):\n b, s = imgs.shape[:2]\n x = imgs.reshape(b * s, *imgs.shape[2:]).float().div_(255.0)\n x = (x - self.mean) / self.std\n feat = self.encoder(x).reshape(b, s, -1)\n return self.head(feat, mask)\n\n\ndef build_encoder():\n \"\"\"Resolve a backbone, and refuse to train from random initialisation silently.\n\n Order: an explicit RSNA_BACKBONE timm name, then a mounted DINOv2 directory, then\n timm with whatever cached weights exist. A frozen self-supervised encoder saturates\n on this task - resolution, encoder size, slice coverage and slot aggregation all\n stop moving the score at the same place, which says the binding constraint is the\n representation, learned as it was on natural images where nothing resembles a torn\n meniscus on a proton-density sequence. So the last blocks are adapted.\n\n The random-initialisation path stays available but shouts. One of the two source\n notebooks builds `resnet18` with `pretrained=False` and blends 88% of its output\n into the submission, which is a randomly initialised network given two epochs on\n 4,407 studies; that is the failure this warning exists to make impossible to miss.\n \"\"\"\n name = os.environ.get(\"RSNA_BACKBONE\")\n if name:\n import timm\n m = timm.create_model(name, pretrained=True, num_classes=0, global_pool=\"avg\")\n log(f\"backbone: timm {name}, dim {m.num_features}\")\n enc = Encoder(m, \"timm\", m.num_features)\n else:\n p = find_dinov2(\"small\")\n if p is not None:\n from transformers import AutoModel\n bb = AutoModel.from_pretrained(str(p))\n enc = Encoder(bb, \"hf\", bb.config.hidden_size * 2)\n log(f\"backbone: DINOv2 from {p}, dim {enc.dim}\")\n else:\n import timm\n fallback = \"resnet18\"\n try:\n m = timm.create_model(fallback, pretrained=True, num_classes=0,\n global_pool=\"avg\")\n log(f\"backbone: timm {fallback} (pretrained), dim {m.num_features}\")\n except Exception as exc:\n m = timm.create_model(fallback, pretrained=False, num_classes=0,\n global_pool=\"avg\")\n warnings.warn(\n \"NO PRETRAINED WEIGHTS AVAILABLE - training from random \"\n f\"initialisation ({exc}). On 4,407 weakly labelled studies this \"\n \"will not learn anything useful. Attach a DINOv2 or timm weights \"\n \"dataset before trusting any score from this run.\",\n RuntimeWarning, stacklevel=2)\n log(\"!! backbone: RANDOM INIT - results are not meaningful\")\n enc = Encoder(m, \"timm\", m.num_features)\n\n for prm in enc.parameters():\n prm.requires_grad = False\n blocks = enc.blocks()\n if blocks:\n # The early blocks of a self-supervised transformer are generic edge and texture\n # filters. There is not enough supervision here to improve them and quite enough\n # to damage them, so only the last few move.\n for blk in blocks[max(0, len(blocks) - CFG.unfreeze_last):]:\n for prm in blk.parameters():\n prm.requires_grad = True\n if enc.kind == \"hf\":\n for prm in enc.module.layernorm.parameters():\n prm.requires_grad = True\n else:\n # A CNN exposes no block list to freeze by depth, and its early layers are far\n # cheaper to relearn than a transformer's, so the whole thing is trained.\n for prm in enc.parameters():\n prm.requires_grad = True\n n_tr = sum(p.numel() for p in enc.parameters() if p.requires_grad)\n where = (f\"last {min(CFG.unfreeze_last, len(blocks))} of {len(blocks)} blocks\"\n if blocks else \"all layers (no block structure to freeze by depth)\")\n log(f\" trainable: {where}, {n_tr / 1e6:.1f}M params\")\n return enc","id":"cell-16"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def augment(imgs):\n \"\"\"Small in-plane affine plus an intensity scale, applied to a whole bag at once.\n\n No horizontal flip: laterality was normalised onto a left-knee convention upstream,\n and a horizontal flip would reintroduce exactly the nuisance axis that removed.\n\n No vertical flip either, which is where this departs from the source notebook. A\n knee coronal or sagittal image has the femur above and the tibia below; flipping it\n top to bottom produces an anatomy that does not exist, and the three OA targets are\n compartment-specific, so the model is being asked to call a finding on a joint whose\n bones have swapped. A rotation of a few degrees is the augmentation that respects\n the acquisition - patient positioning really does vary by that much.\n \"\"\"\n b = imgs.shape[0]\n dev = imgs.device\n x = imgs.float()\n lead = x.shape[:-3]\n x = x.reshape(-1, *x.shape[-3:])\n\n ang = (torch.rand(b, device=dev) - 0.5) * (2 * 10 * np.pi / 180) # +/- 10 degrees\n scale = 1.0 + (torch.rand(b, device=dev) - 0.5) * 0.16 # +/- 8%\n tx = (torch.rand(b, device=dev) - 0.5) * 0.12 # +/- 6%\n ty = (torch.rand(b, device=dev) - 0.5) * 0.12\n cos, sin = torch.cos(ang) / scale, torch.sin(ang) / scale\n theta = torch.stack([torch.stack([cos, -sin, tx], 1),\n torch.stack([sin, cos, ty], 1)], 1) # [b,2,3]\n\n rep = x.shape[0] // b\n theta = theta.repeat_interleave(rep, 0)\n grid = F.affine_grid(theta, x.shape, align_corners=False)\n x = F.grid_sample(x, grid, mode=\"bilinear\", padding_mode=\"zeros\",\n align_corners=False)\n\n gain = 1.0 + (torch.rand(b, 1, 1, 1, device=dev) - 0.5) * 0.2\n x = x * gain.repeat_interleave(rep, 0)\n return x.clamp(0, 255).reshape(*lead, *x.shape[-3:]).to(imgs.dtype)","id":"cell-17"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def protocol_features(hdr, studies):\n \"\"\"Study-level acquisition features from the header pass.\n\n These are available at test time - `test_series.csv` and the test DICOM headers both\n exist at inference - so unlike the reports this is a legal feature, not a teacher.\n It carries real signal because protocol tracks indication: a knee scanned with an\n extra fat-suppressed axial series was scanned by someone looking for something.\n \"\"\"\n if hdr.empty:\n return pd.DataFrame(index=studies).astype(np.float32)\n h = hdr.copy()\n h[\"slot\"] = h[\"plane\"].astype(str).str[0] + \"_\" + h[\"fluid\"].astype(str) + \\\n \"_\" + h[\"fatsat\"].astype(str)\n f = pd.crosstab(h[\"StudyInstanceUID\"], h[\"slot\"])\n agg = h.groupby(\"StudyInstanceUID\").agg(\n n_series=(\"SeriesInstanceUID\", \"size\"),\n n_planes=(\"plane\", \"nunique\"),\n n_fatsat=(\"fatsat\", \"sum\"),\n n_fluid=(\"fluid\", \"sum\"),\n med_slices=(\"n_slices\", \"median\"),\n max_slices=(\"n_slices\", \"max\"),\n med_px=(\"px\", \"median\"),\n )\n f = f.join(agg)\n return f.reindex(studies).fillna(0.0).astype(np.float32)\n\n\ndef protocol_model(feat_tr, y_tr, groups, feat_te, n_folds=5):\n \"\"\"Ridge per target, out-of-fold on the training studies, mean over folds at test.\"\"\"\n from sklearn.linear_model import Ridge\n from sklearn.pipeline import make_pipeline\n from sklearn.preprocessing import StandardScaler\n\n cols = sorted(set(feat_tr.columns) & set(feat_te.columns))\n # A slot column present in no training study is a constant, and a constant column\n # makes the normal equations singular rather than merely uninformative.\n cols = [c for c in cols if feat_tr[c].std() > 1e-8]\n if not cols:\n return np.zeros_like(y_tr), np.full((len(feat_te), y_tr.shape[1]), 0.5, np.float32)\n X, Xt = feat_tr[cols].values, feat_te[cols].values\n oof = np.zeros_like(y_tr)\n pred = np.zeros((len(Xt), y_tr.shape[1]), np.float32)\n for f in range(n_folds):\n va = np.flatnonzero(groups == f)\n tr = np.flatnonzero(groups != f)\n if len(va) == 0 or len(tr) < 10:\n continue\n # SVD, not the default Cholesky: the slot-count columns sum exactly to\n # `n_series`, so the design is perfectly collinear by construction and the\n # normal equations are singular however large alpha is.\n m = make_pipeline(StandardScaler(), Ridge(alpha=8.0, solver=\"svd\"))\n m.fit(X[tr], y_tr[tr])\n oof[va] = m.predict(X[va])\n pred += m.predict(Xt) / n_folds\n return oof, pred","id":"cell-18"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def select_device():\n \"\"\"Pick a device, and prove it can run a kernel before committing the run to it.\n\n `torch.cuda.is_available()` reports only that a driver and a device exist, not that\n the installed build has code for it. Kaggle's P100 is compute capability 6.0 and the\n pinned PyTorch ships kernels for sm_70 and up, so availability returns True and the\n first real launch dies with `no kernel image is available for execution on the\n device`.\n\n Left unchecked that happens inside the first training step - which is after the\n header pass and the twenty-minute cache decode, so the run burns most of an hour of\n GPU quota to discover something knowable in five seconds. Hence a real launch here,\n at the top of the run, in the same autocast path training uses.\n \"\"\"\n if not torch.cuda.is_available():\n log(\"device: cpu (no CUDA device reported)\")\n return torch.device(\"cpu\")\n\n name = torch.cuda.get_device_name(0)\n major, minor = torch.cuda.get_device_capability(0)\n try:\n x = torch.zeros(8, 3, 32, 32, device=\"cuda\")\n w = torch.zeros(4, 3, 3, 3, device=\"cuda\")\n with torch.autocast(\"cuda\"):\n y = F.conv2d(x, w).flatten(1)\n y = y @ y.T\n y.sum().item()\n torch.cuda.synchronize()\n except Exception as exc:\n arches = \" \".join(torch.cuda.get_arch_list())\n raise RuntimeError(\n f\"CUDA device '{name}' (sm_{major}{minor}) reports available but cannot \"\n f\"execute a kernel.\\n underlying error: {exc}\\n\"\n f\" this torch build has kernels for: {arches}\\n\"\n \"This is the known Kaggle P100 mismatch - the default image's PyTorch \"\n \"carries no Pascal kernels. Fix: set \\\"machine_shape\\\": \\\"NvidiaTeslaT4\\\" \"\n \"in kernel-metadata.json, or choose GPU T4 x2 in the notebook's \"\n \"accelerator settings. Failing here rather than after the cache build.\"\n ) from exc\n\n log(f\"device: cuda ({name}, sm_{major}{minor}) - kernel launch verified\")\n return torch.device(\"cuda\")\n\n\ndef macro_auc(y, p):\n from sklearn.metrics import roc_auc_score\n return float(np.nanmean([roc_auc_score(y[:, j], p[:, j])\n if len(set(y[:, j].tolist())) > 1 else np.nan\n for j in range(y.shape[1])]))\n\n\ndef per_target_auc(y, p):\n from sklearn.metrics import roc_auc_score\n return {t: (roc_auc_score(y[:, j], p[:, j])\n if len(set(y[:, j].tolist())) > 1 else float(\"nan\"))\n for j, t in enumerate(TARGETS)}\n\n\n# Below this many positives (or negatives) the Hanley-McNeil normal approximation\n# stops describing the estimate, so the standard error is reported but not trusted.\nMIN_POS_TO_READ = 10\n\n\ndef auc_se(auc, n_pos, n_neg):\n \"\"\"Hanley-McNeil standard error of an AUC, from the counts alone.\n\n A low AUC on a rare label is ambiguous: the model may have learned nothing, or\n the estimate may simply be unreadable because there are nine positives. Those\n demand opposite responses, and the mean AUC cannot tell them apart. This makes\n the second case visible without a bootstrap.\n \"\"\"\n if n_pos < 1 or n_neg < 1 or not np.isfinite(auc):\n return float(\"nan\")\n q1 = auc / (2.0 - auc)\n q2 = 2.0 * auc ** 2 / (1.0 + auc)\n var = (auc * (1 - auc)\n + (n_pos - 1) * (q1 - auc ** 2)\n + (n_neg - 1) * (q2 - auc ** 2)) / (n_pos * n_neg)\n return float(np.sqrt(max(var, 0.0)))\n\n\ndef per_target_report(y, p):\n \"\"\"Per-label AUC with the positive count and standard error beside it.\n\n The competition metric is the unweighted mean of twelve AUCs, so every label is\n worth 1/12 regardless of how often it occurs. The mean alone hides which of the\n twelve has headroom and which is already unmeasurable.\n \"\"\"\n aucs = per_target_auc(y, p)\n out = []\n for j, t in enumerate(TARGETS):\n n_pos, n_neg = int(y[:, j].sum()), len(y) - int(y[:, j].sum())\n auc, se = aucs[t], auc_se(aucs[t], n_pos, len(y) - int(y[:, j].sum()))\n # Two separate reasons a number here may not be actionable, and they are not\n # the same reason. Too few positives means the standard error itself is\n # untrustworthy (the normal approximation needs a real sample, and at n_pos=1\n # it reports a confident interval around nothing). Enough positives but an\n # interval straddling 0.5 means the estimate is sound and the model has\n # genuinely learned nothing. Only the second is worth acting on.\n if n_pos < MIN_POS_TO_READ or n_neg < MIN_POS_TO_READ:\n verdict = f\"too few positives (<{MIN_POS_TO_READ}) - se is unreliable\"\n elif abs(auc - 0.5) <= 2 * se:\n verdict = \"at chance\"\n else:\n verdict = \"measured\"\n out.append({\"target\": t, \"n_pos\": n_pos, \"n_neg\": n_neg,\n \"auc\": auc, \"se\": se, \"verdict\": verdict})\n return pd.DataFrame(out)\n\n\n_NORM_PUNCT = re.compile(r\"[^\\w\\s]\")\n_NORM_NUM = re.compile(r\"\\d+\")\n_NORM_WS = re.compile(r\"\\s+\")\n\n\ndef normalise_report(s):\n \"\"\"Collapse a report to the form its derived labels actually depend on.\n\n Two reports differing only in punctuation, casing or measurement digits produce the\n same twelve labels, so for fold-assignment purposes they are the same document.\n Hashing the raw string treats them as unrelated and scatters them across folds.\n \"\"\"\n s = _NORM_PUNCT.sub(\" \", str(s).lower())\n return _NORM_WS.sub(\" \", _NORM_NUM.sub(\"#\", s)).strip()\n\n\ndef report_fold(report, n_folds):\n \"\"\"Deterministic fold from the normalised report - no state, no ordering.\"\"\"\n key = normalise_report(report) or str(report)\n return int(hashlib.md5(key.encode()).hexdigest()[:8], 16) % n_folds\n\n\ndef audit_folds(reports, folds, y, n_folds):\n \"\"\"Two questions the fold split has to answer before any score off it is believable.\n\n First: does any group of near-identical reports straddle a fold boundary? If so the\n distilled text model sees a document in training and its twin in validation, and the\n out-of-fold correlation that sets `w_text` is measuring memorisation.\n\n Second: does every label have enough positives *and* negatives inside every fold?\n AUC is undefined on a single class and unstable near it, and the competition metric\n weights all twelve equally - so one starved label quietly caps the achievable score.\n \"\"\"\n grp = pd.Series([normalise_report(r) for r in reports])\n span = pd.DataFrame({\"g\": grp, \"f\": folds}).groupby(\"g\")[\"f\"].nunique()\n split = int((span > 1).sum())\n if split:\n n_aff = int(grp.isin(span[span > 1].index).sum())\n log(f\" WARNING: {split} near-duplicate report groups span folds \"\n f\"({n_aff} studies) - text-model OOF is optimistic\")\n else:\n log(f\" duplicate check: every near-duplicate report group is inside one fold\")\n\n yb = (np.asarray(y) > 0.5).astype(int)\n starved = []\n for f in range(n_folds):\n va = np.flatnonzero(np.asarray(folds) == f)\n if len(va) == 0:\n continue\n for j, t in enumerate(TARGETS):\n npos = int(yb[va, j].sum())\n nneg = len(va) - npos\n if min(npos, nneg) < MIN_POS_TO_READ:\n starved.append((f, t, npos, nneg))\n if starved:\n log(f\" WARNING: {len(starved)} (fold, label) cells below \"\n f\"{MIN_POS_TO_READ} positives or negatives:\")\n for f, t, npos, nneg in starved[:12]:\n log(f\" fold {f} {t:<18} pos {npos:>4} neg {nneg:>4}\")\n else:\n log(f\" class balance: all {n_folds} x {len(TARGETS)} cells have \"\n f\">= {MIN_POS_TO_READ} of each class\")\n return {\"groups_split\": split, \"starved_cells\": len(starved)}\n\n\ndef take_group(rows, g):\n return rows[:, :, g * CFG.group:(g + 1) * CFG.group]\n\n\n@torch.no_grad()\ndef predict(model, cache, mask, idx, dev, n_group):\n \"\"\"Average the logits over the groups of each slot.\n\n Training sees one group at a time, which doubles as augmentation along the stack;\n inference averages over all of them.\n \"\"\"\n model.eval()\n out = []\n for b in range(0, len(idx), CFG.eval_batch):\n sel = idx[b:b + CFG.eval_batch]\n rows = torch.from_numpy(cache[sel]).to(dev)\n m = torch.from_numpy(mask[sel]).to(dev)\n acc = None\n for g in range(n_group):\n with torch.autocast(\"cuda\", enabled=dev.type == \"cuda\"):\n z = model(take_group(rows, g), m).float()\n acc = z if acc is None else acc + z\n out.append(torch.sigmoid(acc / n_group).cpu().numpy())\n return np.concatenate(out) if out else np.zeros((0, len(TARGETS)), np.float32)\n\n\ndef rank_mean(mats, weights=None):\n \"\"\"Average several prediction matrices in rank space.\n\n AUC is invariant under any strictly increasing map of a column, so calibration is\n worth nothing and only order is read. Averaging raw probabilities lets whichever\n model happens to be most confident dominate; averaging ranks combines the only\n information the metric reads.\n \"\"\"\n mats = [np.asarray(m, np.float64) for m in mats]\n weights = np.ones(len(mats)) if weights is None else np.asarray(weights, np.float64)\n weights = weights / weights.sum()\n acc = np.zeros_like(mats[0])\n for w, m in zip(weights, mats):\n acc += w * pd.DataFrame(m).rank(pct=True).values\n return acc","id":"cell-19"},{"cell_type":"markdown","metadata":{},"source":"## 6. Validating without fooling yourself\n\n**Shared reports.** Some reports are byte-identical across studies — a template read for\nan unremarkable knee — and every study in such a group gets the same derived target.\nSplit the group across folds and the model is scored on a target whose source it trained\non. Folds are assigned by a hash of the report text, which keeps duplicates whole.\n\n**Two references, two meanings.** Performance is reported against the derived targets,\nwhich cover every study and measure whether the imaging model learned what the text\nsays; and against the annotated studies, which are far fewer and measure agreement with\na reading of the *images*. The second resembles the test set. The first has enough\npositives per label to tell a real difference from noise.\n\n**The annotated reference has to be held out.** This is where the source notebook leaks:\nit scores the annotated AUC over *every* gold study, including the ones in the training\nsplit, and then uses that number to choose the checkpoint. Here gold studies take their\nfold assignment like everything else and only held-out gold is scored.\n\n**Choosing an epoch.** The obvious rule — best epoch on the larger reference — is wrong,\nbecause the two measure different things and an epoch that gains on one while losing on\nthe other has been shown to be *different*, not better. Ranking by\n\n$$\\text{score}(e) = \\min\\big(\\mathrm{AUC}^\\text{derived}_e,\\ \\mathrm{AUC}^\\text{annot}_e\\big)$$\n\nrefuses that trade, and makes training length a measured quantity rather than a guess.","id":"cell-20"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def train_fold(cache, mask, Y, W, tr, va, gold_idx, gold_y, n_group, dev, fold):\n \"\"\"Fit one fold and return (best model, val prediction, chosen-epoch diagnostics).\n\n Epoch selection reads two references by their *worse* value. They measure different\n things - the derived targets cover every study and say whether the imaging model\n learned what the text says, the annotated subset is far smaller and says whether it\n agrees with a reading of the images - so an epoch that gains on one while losing on\n the other has not been shown to be better, only different. Taking the minimum\n refuses that trade, and makes training length a measured quantity rather than a\n hyperparameter guessed in advance.\n \"\"\"\n seed_all(CFG.train_seed + fold)\n model = Model(build_encoder()).to(dev)\n enc_params = [p for p in model.encoder.parameters() if p.requires_grad]\n groups = [{\"params\": model.head.parameters(), \"lr\": CFG.lr_head}]\n max_lr = [CFG.lr_head]\n if enc_params:\n groups.insert(0, {\"params\": enc_params, \"lr\": CFG.lr_backbone})\n max_lr.insert(0, CFG.lr_backbone)\n opt = torch.optim.AdamW(groups, weight_decay=CFG.weight_decay)\n\n steps = max(CFG.epochs * (len(tr) // CFG.batch_studies), 1)\n sched = torch.optim.lr_scheduler.OneCycleLR(opt, max_lr=max_lr, total_steps=steps,\n pct_start=0.15)\n scaler = torch.amp.GradScaler(\"cuda\", enabled=dev.type == \"cuda\")\n\n yv = (Y[va] > 0.5).astype(int)\n best, best_state, best_ep = -1.0, None, -1\n stepped = 0\n\n for ep in range(CFG.epochs):\n model.train()\n perm = np.random.permutation(tr)\n tot, nstep = 0.0, 0\n for b in range(0, len(perm) - CFG.batch_studies + 1, CFG.batch_studies):\n sel = perm[b:b + CFG.batch_studies]\n rows = torch.from_numpy(cache[sel]).to(dev)\n g = int(torch.randint(n_group, (1,)).item())\n imgs = augment(take_group(rows, g))\n m = torch.from_numpy(mask[sel]).to(dev)\n y = torch.from_numpy(Y[sel]).to(dev)\n w = torch.from_numpy(W[sel]).to(dev)\n with torch.autocast(\"cuda\", enabled=dev.type == \"cuda\"):\n loss = (F.binary_cross_entropy_with_logits(\n model(imgs, m), y, reduction=\"none\") * w).mean()\n opt.zero_grad(set_to_none=True)\n scaler.scale(loss).backward()\n scaler.step(opt)\n scaler.update()\n if stepped < steps:\n sched.step()\n stepped += 1\n tot += loss.item()\n nstep += 1\n\n pv = predict(model, cache, mask, va, dev, n_group)\n d_auc = macro_auc(yv, pv)\n\n # The annotated reference is evaluated on the gold studies *in this fold's\n # validation split only*. Scoring it over every gold study - which is what the\n # source notebook does - measures the model on rows it just trained on, and that\n # number then chooses the checkpoint.\n g_auc = float(\"nan\")\n gv = np.intersect1d(gold_idx, va)\n if len(gv) >= 8:\n gy = gold_y[np.searchsorted(gold_idx, gv)]\n g_auc = macro_auc(gy, predict(model, cache, mask, gv, dev, n_group))\n\n score = d_auc if not np.isfinite(g_auc) else min(d_auc, g_auc)\n log(f\" fold {fold} epoch {ep + 1}/{CFG.epochs} loss {tot / max(nstep, 1):.4f}\"\n f\" derived {d_auc:.4f} annot {g_auc:.4f} worse-of-two {score:.4f}\")\n if score > best:\n best, best_ep = score, ep\n best_state = {k: v.detach().cpu().clone()\n for k, v in model.state_dict().items()}\n if time.time() - T0 > CFG.time_budget:\n log(\" time budget reached inside fold\")\n break\n\n if best_state is not None:\n model.load_state_dict(best_state)\n log(f\" fold {fold}: restored epoch {best_ep + 1} (worse-of-two {best:.4f})\")\n return model, predict(model, cache, mask, va, dev, n_group), best","id":"cell-21"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def write_submission(study_ids, pred, test_df, path=\"submission.csv\"):\n sub = pd.DataFrame(pred, columns=TARGETS)\n sub.insert(0, \"StudyInstanceUID\", study_ids)\n sub = test_df[[\"StudyInstanceUID\"]].merge(sub, on=\"StudyInstanceUID\", how=\"left\")\n sub[TARGETS] = sub[TARGETS].fillna(0.5)\n sub.to_csv(path, index=False)\n log(f\"{path} {sub.shape}; nulls {int(sub[TARGETS].isna().sum().sum())}\")\n return sub\n\n\ndef write_benchmark(root, path=\"submission.csv\"):\n \"\"\"Write the 0.5 benchmark file immediately.\n\n A submission that never writes scores nothing at all, which is strictly worse than\n scoring badly. A try/except covers exceptions, but a kill for memory is a SIGKILL and\n never reaches it, so a valid file has to exist from the first second.\n \"\"\"\n t = pd.read_csv(Path(root) / \"test.csv\")\n for c in TARGETS:\n t[c] = 0.5\n t.to_csv(path, index=False)","id":"cell-22"},{"cell_type":"markdown","metadata":{},"source":"## 7. Run\n\nA valid `submission.csv` is written before anything else and overwritten only once real\npredictions exist: a submission that never writes scores nothing at all, which is\nstrictly worse than scoring badly, and a kill for memory is a SIGKILL that no\n`try`/`except` will catch.\n\nKnobs worth touching, all via environment variables so nothing below needs editing:\n`RSNA_FOLDS` (default 5 — drop to 2 if the run is time-bound), `RSNA_TIME_BUDGET`\n(seconds, default 8 h), `RSNA_BACKBONE` (a timm name, overrides the DINOv2 search),\n`SLOT_SCHEME=public` (ablate the header recovery back to the provided flags).","id":"cell-23"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def main(root=None):\n root = find_root(root)\n log(f\"input root: {root}\")\n write_benchmark(root)\n seed_all()\n log(f\"seeds: data/folds {CFG.seed}, training {CFG.train_seed}\")\n # Before the header pass and the cache decode, not after: an unusable accelerator\n # is knowable now and costs most of an hour to discover later.\n dev = select_device()\n\n test_df = pd.read_csv(root / \"test.csv\")\n train_df = pd.read_csv(root / \"train.csv\")\n train_series = pd.read_csv(root / \"train_series.csv\")\n test_series = pd.read_csv(root / \"test_series.csv\")\n log(f\"train {train_df.shape} test {test_df.shape}\")\n\n both = pd.concat([train_series, test_series])\n plane_map = dict(zip(both[\"SeriesInstanceUID\"], both[\"Anatomical_Plane\"]))\n\n log(\"header pass: test\")\n hte = annotate(walk(root, \"test_series\"))\n log(f\" {len(hte)} test series\")\n log(\"header pass: train\")\n htr = annotate(walk(root, \"train_series\"))\n log(f\" {len(htr)} train series\")\n for h in (htr, hte):\n if not h.empty:\n h[\"plane\"] = h[\"SeriesInstanceUID\"].map(plane_map)\n\n def lat_of(h, tag=\"\"):\n \"\"\"Study -> 'L' / 'R' / None, with the evidence source counted per study.\"\"\"\n LAT_STATS.update({k: 0 for k in LAT_STATS})\n d = {}\n if h.empty:\n return d\n for st, g in h.groupby(\"StudyInstanceUID\"):\n d[st] = resolve_laterality(g)\n\n n = len(d) or 1\n log(f\"{tag}: laterality \" + \", \".join(\n f\"{k} {v} ({v / n:.1%})\" for k, v in LAT_STATS.items()\n if v and not k.startswith(\"geom_\")))\n both = LAT_STATS[\"geom_agree\"] + LAT_STATS[\"geom_disagree\"]\n if both:\n # A fallback validated against the source it is meant to replace. If this\n # agreement is not high, the geometric rule is wrong for this corpus and\n # RSNA_LAT_GEOMETRY must stay off.\n log(f\"{tag}: geometric rule agrees with the tag on \"\n f\"{LAT_STATS['geom_agree']}/{both} studies \"\n f\"({LAT_STATS['geom_agree'] / both:.1%})\")\n unknown = LAT_STATS[\"unknown\"] / n\n if unknown > 0.10:\n warnings.warn(\n f\"{unknown:.0%} of {tag} studies have no laterality, so they are never \"\n \"mirrored to a common convention. Medial and lateral then sit on \"\n \"opposite sides of the image for left and right knees, and four of the \"\n \"twelve targets are side-specific.\", RuntimeWarning, stacklevel=2)\n return d\n\n slots_tr, slots_te = pick_slots(htr, plane_map), pick_slots(hte, plane_map)\n cov = pd.Series([len(v) for v in slots_tr.values()] or [0]).describe()\n log(f\"train slots per study: mean {cov['mean']:.2f} min {cov['min']:.0f} \"\n f\"max {cov['max']:.0f}\")\n\n n_group = plan_cache(len(train_df))\n log(f\"cache layout: {n_group} groups x {CFG.group} slices\")\n lat_tr, lat_te = lat_of(htr, \"train\"), lat_of(hte, \"test\")\n st_tr, Ctr, Mtr = build_cache(slots_tr, plane_map, lat_tr, \"train\", n_group)\n st_te, Cte, Mte = build_cache(slots_te, plane_map, lat_te, \"test\", n_group)\n\n # ---- targets ----------------------------------------------------------- #\n t_lab = time.time()\n rules = pd.DataFrame([extract(r) for r in train_df[\"Report\"].fillna(\"\")])\n rules[\"StudyInstanceUID\"] = train_df[\"StudyInstanceUID\"].values\n rules = rules.set_index(\"StudyInstanceUID\")\n log(f\"rule labels for {len(rules)} studies in {time.time() - t_lab:.1f}s\")\n\n # Folds are assigned by a hash of the report text, so a template read for an\n # unremarkable knee stays whole inside one fold. Split such a group and the model is\n # scored on a target whose source it has trained on. Hashing the *normalised* text\n # rather than the raw bytes: measured on this corpus, 25 groups covering 125 studies\n # differ only in punctuation, casing or measurement digits, and a raw hash scattered\n # them across folds while their twelve derived labels were identical.\n rep = train_df.set_index(\"StudyInstanceUID\")[\"Report\"].fillna(\"\")\n fold_of = {s: report_fold(rep.get(s, s), CFG.n_folds) for s in rep.index}\n\n # The distilled text model shares those folds for the same reason.\n order = list(rules.index)\n rule_mat = rules[TARGETS].values.astype(np.float32)\n text_oof, text_corr = distill([rep.get(s, \"\") for s in order], rule_mat,\n np.array([fold_of.get(s, 0) for s in order]),\n n_splits=CFG.n_folds, seed=CFG.seed)\n # Weight the distilled model per target by how well it recovered the rules out of\n # fold. A target where it failed contributes nothing rather than noise.\n w_text = CFG.w_text * np.clip(text_corr, 0.0, 1.0)\n w_text[text_corr < 0.15] = 0.0\n log(f\"text-model weight per target: min {w_text.min():.2f} max {w_text.max():.2f}, \"\n f\"{int((w_text == 0).sum())}/{len(TARGETS)} gated out\")\n derived = blend(rule_mat, text_oof, w_text=w_text)\n derived = pd.DataFrame(derived, index=order, columns=TARGETS)\n conf = rules[[t + \"__conf\" for t in TARGETS]].values.astype(np.float32)\n conf = pd.DataFrame(conf, index=order, columns=TARGETS)\n\n gold = train_df.set_index(\"StudyInstanceUID\")[TARGETS]\n gold = gold[gold.notna().all(axis=1)]\n log(f\"annotated studies: {len(gold)}\")\n\n Y = np.zeros((len(st_tr), len(TARGETS)), np.float32)\n W = np.zeros_like(Y)\n for i, st in enumerate(st_tr):\n if st in gold.index:\n Y[i], W[i] = gold.loc[st].values, CFG.gold_weight\n elif st in derived.index:\n Y[i] = derived.loc[st].values\n # Silence on a finding is weak evidence, and should pull weakly.\n W[i] = 0.25 + 0.75 * conf.loc[st].values\n keep = np.where(W.sum(1) > 0)[0]\n folds = np.array([fold_of.get(s, 0) for s in st_tr])\n log(f\"supervised {len(keep)} of {len(st_tr)} studies\")\n\n log(\"fold audit:\")\n audit_folds([rep.get(s, \"\") for s in st_tr], folds, Y, CFG.n_folds)\n\n gold_pos = np.array(sorted(i for i, s in enumerate(st_tr) if s in gold.index))\n gold_y = (gold.loc[[st_tr[i] for i in gold_pos]].values.astype(int)\n if len(gold_pos) else np.zeros((0, len(TARGETS)), int))\n\n # ---- protocol-only model ------------------------------------------------ #\n feat_tr = protocol_features(htr, st_tr)\n feat_te = protocol_features(hte, st_te)\n proto_oof, proto_te = protocol_model(feat_tr, Y, folds, feat_te, CFG.n_folds)\n proto_auc = macro_auc((Y[keep] > 0.5).astype(int), proto_oof[keep])\n # Diversity is only worth having if it is diversity around something. Protocol\n # tracks indication - a knee given an extra fat-suppressed axial series was scanned\n # by someone looking for something - but that is an empirical claim, so it has to\n # clear chance out-of-fold before it is allowed into the fusion.\n w_protocol = CFG.w_protocol if proto_auc > 0.52 else 0.0\n log(f\"protocol model: {feat_tr.shape[1]} features, derived-OOF macro AUC \"\n f\"{proto_auc:.4f} -> fusion weight {w_protocol:.2f}\")\n\n # ---- image model, grouped K-fold ---------------------------------------- #\n oof = np.full((len(st_tr), len(TARGETS)), np.nan, np.float32)\n test_preds, fold_scores = [], []\n\n for f in range(min(CFG.max_folds_to_run, CFG.n_folds)):\n va = np.array([i for i in keep if folds[i] == f])\n tr = np.array([i for i in keep if folds[i] != f])\n if len(va) == 0 or len(tr) < CFG.batch_studies:\n log(f\"fold {f}: too small, skipped\")\n continue\n log(f\"fold {f}: train {len(tr)} / val {len(va)} \"\n f\"(annotated in val: {len(np.intersect1d(gold_pos, va))})\")\n model, pv, score = train_fold(Ctr, Mtr, Y, W, tr, va, gold_pos, gold_y,\n n_group, dev, f)\n oof[va] = pv\n fold_scores.append(score)\n test_preds.append(predict(model, Cte, Mte, np.arange(len(st_te)), dev, n_group))\n del model\n gc.collect()\n if dev.type == \"cuda\":\n torch.cuda.empty_cache()\n if time.time() - T0 > CFG.time_budget:\n log(\"time budget reached; stopping fold loop\")\n break\n\n if not test_preds:\n log(\"no fold completed - leaving the benchmark submission in place\")\n return None\n\n # ---- honest out-of-fold report ------------------------------------------ #\n done = np.flatnonzero(np.isfinite(oof).all(1))\n yd = (Y[done] > 0.5).astype(int)\n log(f\"OOF over {len(done)} studies\")\n log(f\" vs derived targets : {macro_auc(yd, oof[done]):.4f}\")\n gd = np.intersect1d(gold_pos, done)\n if len(gd) >= 8:\n gy = gold_y[np.searchsorted(gold_pos, gd)]\n log(f\" vs annotated ({len(gd)}) : {macro_auc(gy, oof[gd]):.4f}\")\n\n # Per-label, on the derived targets rather than the handful of annotated ones -\n # the metric is the mean of twelve, so the breakdown is what says where the\n # headroom is. `se` is the width of the estimate: a label whose AUC sits within\n # about two standard errors of 0.500 has not been measured, only sampled.\n rep = per_target_report(yd, oof[done])\n log(f\" per label (derived targets, n={len(done)}):\")\n log(f\" {'target':<18} {'n_pos':>6} {'auc':>7} {'se':>7} verdict\")\n for r in rep.itertuples():\n log(f\" {r.target:<18} {r.n_pos:>6} {r.auc:>7.3f} {r.se:>7.3f} {r.verdict}\")\n\n # Persist the raw out-of-fold matrix so the next question about these predictions\n # can be answered without spending another run producing them.\n oof_path = \"oof_predictions.csv\"\n pd.concat([pd.DataFrame({\"StudyInstanceUID\": np.asarray(st_tr)[done]}),\n pd.DataFrame(oof[done], columns=[f\"pred_{t}\" for t in TARGETS]),\n pd.DataFrame(yd, columns=[f\"target_{t}\" for t in TARGETS])],\n axis=1).to_csv(oof_path, index=False)\n rep.to_csv(\"oof_per_label.csv\", index=False)\n log(f\" wrote {oof_path} and oof_per_label.csv\")\n\n # ---- fuse and submit ----------------------------------------------------- #\n image_te = rank_mean(test_preds)\n final = (image_te if w_protocol == 0.0\n else rank_mean([image_te, proto_te],\n weights=[1.0 - w_protocol, w_protocol]))\n sub = write_submission(st_te, final, test_df)\n log(f\"folds run: {len(test_preds)}, mean worse-of-two {np.mean(fold_scores):.4f}\")\n return sub\n log(\"done\")","id":"cell-24"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"\nimport traceback\n\ntry:\n sub = main()\nexcept Exception:\n traceback.print_exc()\n # Fall back to the benchmark file rather than dying: an exception here still leaves\n # a scoreable submission behind.\n write_benchmark(find_root())\n sub = pd.read_csv(\"submission.csv\")\n print(\"wrote fallback submission.csv\")\nlog(\"done\")\nsub.head(20)","id":"cell-25"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11"}},"nbformat":4,"nbformat_minor":5}