{"cells":[{"cell_type":"markdown","id":"a8fdf11a","metadata":{},"source":"# RSNA Knee: the data, then a simple baseline\n\nThis notebook reads the competition data, shows what is in it, trains one small image model,\nand writes a submission. Everything runs in this notebook on one T4 GPU with the internet\nswitched off.\n\nThe short version of what the data forces on you:\n\n1. There are 4,407 training studies and only 58 of them carry the twelve labels. You cannot\n train an image model on 58 studies, so the labels have to come from somewhere else.\n2. Every training study carries the radiology report. `test.csv` has no `Report` column. Text is\n available while you fit and gone when you predict, so a report can only ever be a training\n label, never an input.\n3. The score is the average of twelve ROC AUC values. Only the order of your predictions inside\n each column is read, so calibration and thresholds are worth nothing.\n\nThe model here is a ResNet18 over twelve sagittal slices per study. What it reaches:\n\n| Measurement | Average AUC |\n|---|---|\n| Model against the report labels, on a held out fifth of the studies | about 0.78 |\n| Model against a radiologist, on the 58 labelled studies, none of which it trained on | 0.76 to 0.79 |\n| Public leaderboard score of this notebook | 0.798 |\n| The report labels themselves against the radiologist, on the same 58 | 0.90 |\n\nThe second row is a range because that is what 58 studies buy. Rerunning the notebook with one\nstudy moved between the two sides changed it by 0.03, which is the size section 3 predicts from\nthe sample alone. The holdout number in the first row moved by 0.005 over the same change.\n\nThe third row is the reason the second row is still worth printing. The 58 studies estimated\n0.778 and the leaderboard came back 0.798. They cannot rank two models that are 0.02 apart, but\nthey did say this model works, and the estimate landed close.\n\nThe last row is the ceiling on the second. The model is learning from labels that are 0.90\nagainst a radiologist, so better labels move it further than a bigger network does. Section 12\nsays what to change first.\n\nThe top of the leaderboard is around 0.95, so this is a starting point rather than a contender.\nThe whole notebook, reading all 4,407 studies and training from scratch, takes about 15 minutes\non one T4."},{"cell_type":"code","id":"400f9127","metadata":{},"execution_count":null,"outputs":[],"source":"import gc, hashlib, pathlib, re, time, unicodedata, warnings\nfrom concurrent.futures import ThreadPoolExecutor\n\nimport numpy as np\nimport pandas as pd\n\nwarnings.filterwarnings(\"ignore\")\nT0 = time.time()\n\ndef log(msg):\n print(f\"[{time.time() - T0:7.1f}s] {msg}\", flush=True)\n\n# One session on this competition gets a fixed amount of wall clock. Every stage below checks\n# these numbers, so a slow read shrinks the training set instead of killing the run.\nHARD_DEADLINE = T0 + 8.0 * 3600\nPREP_BUDGET_S = 2.5 * 3600\n\nN_SLICES = 12 # slices kept per study\nSIZE = 224 # pixels per side after the resize\nEPOCHS = 8\nBATCH = 32\nHOLDOUT = 0.20\nSEED = 0\nREAD_THREADS = 12\n\nLABELS = [\"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\", \"Medial OA\", \"Lateral OA\",\n \"PF OA\", \"Effusion\", \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\"]\n\n# The competition mounts under /kaggle/input/competitions/ on this image, and under\n# /kaggle/input/ on older ones.\nCOMP = next(p for p in map(pathlib.Path, [\n \"/kaggle/input/competitions/rsna-knee-abnormality-detection\",\n \"/kaggle/input/rsna-knee-abnormality-detection\",\n \"/tmp/rsnaknee\"]) if (p / \"train.csv\").exists())\n\ntrain = pd.read_csv(COMP / \"train.csv\")\ntrain_series = pd.read_csv(COMP / \"train_series.csv\")\ntest = pd.read_csv(COMP / \"test.csv\")\ntest_series = pd.read_csv(COMP / \"test_series.csv\")\n\nreports = train.set_index(\"StudyInstanceUID\").Report.fillna(\"\")\ngold = train[train[LABELS].notna().all(axis=1)].set_index(\"StudyInstanceUID\")[LABELS].astype(int)\n\nlog(f\"data root {COMP}\")\nlog(f\"train studies {len(train):,} train series {len(train_series):,}\")\nlog(f\"test studies {len(test):,} test series {len(test_series):,}\")\nlog(f\"labelled studies {len(gold)} ({len(gold) / len(train):.2%} of training)\")"},{"cell_type":"markdown","id":"5d343329","metadata":{},"source":"## Write a valid submission before doing anything else\n\nThis is a code competition. When you submit, Kaggle runs this notebook again against a test set\nyou have never seen, and a run that fails part way through scores nothing. So the first thing\nthe notebook does is write a submission of all 0.5 values. Every later stage overwrites it. If\nthe model stage fails, you still get a scored run and you can read the log to see why."},{"cell_type":"code","id":"1c334811","metadata":{},"execution_count":null,"outputs":[],"source":"def write_submission(pred: pd.DataFrame, path=\"submission.csv\") -> pd.DataFrame:\n \"\"\"Rank inside each column, then write one row per test study.\n\n The score reads only the order of the values in a column, so replacing the numbers by their\n ranks throws nothing away, and it removes any question about how the probabilities are\n spread out.\n \"\"\"\n sub = pred.rank(pct=True) if len(pred) > 1 else pred.copy()\n sub.index.name = \"StudyInstanceUID\"\n sub = sub.reindex(test.StudyInstanceUID).fillna(0.5).reset_index()\n sub = sub[[\"StudyInstanceUID\"] + LABELS]\n assert len(sub) == len(test), f\"{len(sub)} rows against {len(test)} test studies\"\n assert sub[LABELS].notna().all().all(), \"a prediction is missing\"\n sub.to_csv(path, index=False)\n return sub\n\nwrite_submission(pd.DataFrame(0.5, index=test.StudyInstanceUID, columns=LABELS))\nlog(\"wrote a placeholder submission.csv\")"},{"cell_type":"code","id":"e20fdcff","metadata":{},"execution_count":null,"outputs":[],"source":"import matplotlib.pyplot as plt\nfrom matplotlib.colors import LinearSegmentedColormap\n\n# One fixed set of colours, used in the same order everywhere below.\nBLUE, ORANGE, AQUA, YELLOW, PINK = \"#2a78d6\", \"#eb6834\", \"#1baf7a\", \"#eda100\", \"#e87ba4\"\nINK, MUTED, GRID = \"#0b0b0b\", \"#898781\", \"#e1e0d9\"\nSEQ = LinearSegmentedColormap.from_list(\"seq\", [\"#e8f1fd\", \"#9ec5f4\", \"#2a78d6\", \"#0d366b\"])\nDIV = LinearSegmentedColormap.from_list(\"div\", [\"#2a78d6\", \"#f0efec\", \"#d03b3b\"])\n\nplt.rcParams.update({\n \"figure.dpi\": 120, \"figure.facecolor\": \"#fcfcfb\", \"axes.facecolor\": \"#fcfcfb\",\n \"font.size\": 8.5, \"axes.grid\": True, \"grid.color\": GRID, \"grid.linewidth\": 0.6,\n \"axes.spines.top\": False, \"axes.spines.right\": False, \"axes.edgecolor\": \"#c3c2b7\",\n \"text.color\": INK, \"axes.labelcolor\": INK, \"xtick.color\": MUTED, \"ytick.color\": MUTED,\n \"axes.titlesize\": 9.5, \"legend.frameon\": False,\n})\n\n# A printed DataFrame wraps at the cell width and becomes unreadable, so every table below\n# goes through display() as HTML instead.\nfrom html import escape as _esc\nfrom IPython.display import display, HTML, Markdown\n\nTABLE_CSS = [\n {\"selector\": \"caption\", \"props\": [(\"caption-side\", \"top\"), (\"text-align\", \"left\"),\n (\"font-weight\", \"600\"), (\"font-size\", \"0.95rem\"),\n (\"padding\", \"0 0 0.45rem 0\"), (\"color\", INK)]},\n {\"selector\": \"th\", \"props\": [(\"background-color\", \"#f2f1ec\"), (\"color\", INK),\n (\"font-weight\", \"600\"), (\"text-align\", \"right\"),\n (\"padding\", \"5px 11px\"), (\"border-bottom\", f\"1px solid #c3c2b7\")]},\n {\"selector\": \"th.row_heading\", \"props\": [(\"text-align\", \"left\"),\n (\"background-color\", \"#fcfcfb\")]},\n {\"selector\": \"td\", \"props\": [(\"padding\", \"5px 11px\"), (\"text-align\", \"right\"),\n (\"border-bottom\", f\"1px solid {GRID}\")]},\n {\"selector\": \"\", \"props\": [(\"border-collapse\", \"collapse\"), (\"font-size\", \"0.86rem\"),\n (\"font-variant-numeric\", \"tabular-nums\"), (\"margin\", \"0.3rem 0\")]},\n]\n\ndef _cell(v):\n \"\"\"One value as text. Integers get thousands separators, floats three decimals.\"\"\"\n if isinstance(v, (bool, np.bool_)):\n return \"yes\" if v else \"no\"\n if isinstance(v, (int, np.integer)):\n return f\"{v:,}\"\n if isinstance(v, (float, np.floating)):\n if not np.isfinite(v):\n return \"\"\n return f\"{v:,.3f}\" if abs(v) < 1e4 else f\"{v:,.0f}\"\n return str(v)\n\ndef show(df, caption=\"\", grad=None, bars=None, vmin=None, vmax=None, hide_index=False):\n \"\"\"Display a table as HTML, optionally shaded or with in cell bars.\"\"\"\n st = (pd.DataFrame(df).style\n .set_caption(caption).set_table_styles(TABLE_CSS)\n .format(_cell, na_rep=\"\"))\n if grad is not None:\n st = st.background_gradient(cmap=SEQ, subset=grad, vmin=vmin, vmax=vmax)\n if bars is not None:\n st = st.bar(subset=bars, color=\"#cde2fb\", vmin=vmin, vmax=vmax)\n if hide_index:\n st = st.hide(axis=\"index\")\n display(st)\n\ndef facts(pairs, caption=\"\"):\n \"\"\"A two column table for a handful of single numbers.\"\"\"\n show(pd.DataFrame({\"value\": [v for _, v in pairs]}, index=[k for k, _ in pairs]), caption)\n\ndef note(md):\n display(Markdown(md))\n\ndef quote(text, caption=\"\"):\n \"\"\"Show a report with its line breaks intact.\n\n A markdown blockquote only quotes up to the first newline, and these reports are several\n lines long, so the rest would render as ordinary text.\n \"\"\"\n display(HTML(\n f\"
\"\n f\"
{_esc(caption)}
\"\n f\"
{_esc(str(text))}
\"))\n\ndef short_uid(u, n=10):\n \"\"\"The last n characters of a study identifier, which is enough to tell rows apart.\"\"\"\n return \"..\" + str(u)[-n:]"},{"cell_type":"markdown","id":"36b1c793","metadata":{},"source":"## 1. The task and how it is scored\n\nEach study is one knee MRI session. A session holds several series, and a series is a stack of\nimages taken with one set of scanner settings. You have to give each study twelve numbers, one\nfor each finding:\n\n| Column | What it means |\n|---|---|\n| ACL, MCL | Tear of the anterior cruciate or medial collateral ligament |\n| Medial Meniscus, Lateral Meniscus | Tear of the meniscus on the inner or outer side |\n| Medial OA, Lateral OA, PF OA | Osteoarthritis in the inner, outer, or kneecap compartment |\n| Effusion | Extra fluid in the joint |\n| Synovitis | Inflamed joint lining |\n| Baker's | A fluid filled cyst behind the knee |\n| Contusion | A bone bruise |\n| Fracture | A break in the bone |\n\nThe score is the average of twelve ROC AUC values, one per column. ROC AUC is the chance that a\nstudy which has the finding gets a higher number than a study which does not. A score of 0.5 is\nthe same as guessing and 1.0 is perfect.\n\nTwo things follow from that, and both remove work rather than adding it.\n\nThe first is that only the order of your numbers inside a column is read. If you sort a column\nand replace each value by its position, the score does not change. So you do not need to\ncalibrate anything, you do not need to pick a threshold, and when you combine two models you\nshould average their ranks rather than their probabilities. Averaging probabilities lets\nwhichever model happens to be more confident decide the answer.\n\nThe second is that all twelve findings cost the same. Write M for the average AUC a good model\nreaches. A column left at 0.5 gives up (M minus 0.5) divided by 12, no matter how well the\nother eleven do. At M of 0.85 that is 0.029 of your final score, which is much larger than the\ngap between neighbouring places near the top of the board. So a rare finding deserves more\nattention than a common one, because a rare finding is where a model most easily ends up\nguessing."},{"cell_type":"markdown","id":"a9837a28","metadata":{},"source":"## 2. What is in the files\n\nThere are five CSV files and two folders of images. `train.csv` holds the report and the twelve\nlabel columns. `train_series.csv` describes every series. The images sit at\n`train_series///.dcm`, one DICOM file per slice."},{"cell_type":"code","id":"fc493658","metadata":{},"execution_count":null,"outputs":[],"source":"show(pd.DataFrame([\n {\"file\": f, \"rows\": len(d), \"columns\": len(d.columns),\n \"carries the report\": \"Report\" in d.columns,\n \"carries the twelve labels\": all(c in d.columns for c in LABELS)}\n for f, d in [(\"train.csv\", train), (\"train_series.csv\", train_series),\n (\"test.csv\", test), (\"test_series.csv\", test_series)]]),\n \"What each file holds\", hide_index=True)\n\nshow(pd.DataFrame({\"column in train.csv\": list(train.columns),\n \"also in test.csv\": [c in test.columns for c in train.columns],\n \"filled in\": [f\"{train[c].notna().mean():.1%}\" for c in train.columns]}),\n \"Every column of train.csv, and whether it survives into test.csv\", hide_index=True)\n\nquote(reports.iloc[0][:500], \"One report, the first 500 characters\")\n\nsample = pd.read_csv(COMP / \"sample_submission.csv\").head(3)\nsample.insert(0, \"study\", sample.StudyInstanceUID.map(short_uid))\nshow(sample.drop(columns=\"StudyInstanceUID\"),\n \"sample_submission.csv, with the study identifier cut to its last ten characters\",\n hide_index=True)"},{"cell_type":"markdown","id":"b62503f3","metadata":{},"source":"`test.csv` has one column and it is the study identifier. That single fact decides the shape of\nevery solution to this competition. The report is available while you fit and absent when you\npredict, so you cannot build a model that reads text and images together. The report can only\nbe used to make training labels."},{"cell_type":"code","id":"d0af93b2","metadata":{},"execution_count":null,"outputs":[],"source":"per_study = train_series.groupby(\"StudyInstanceUID\").size()\nplanes = train_series.Anatomical_Plane.value_counts()\nplane_per_study = train_series.pivot_table(index=\"StudyInstanceUID\", columns=\"Anatomical_Plane\",\n values=\"SeriesInstanceUID\", aggfunc=\"count\").fillna(0)\n\nfig, (a1, a2, a3) = plt.subplots(1, 3, figsize=(11, 2.9))\n\na1.hist(per_study, bins=range(2, 16), color=BLUE, rwidth=0.8)\na1.set_xlabel(\"series in one study\"); a1.set_ylabel(\"studies\")\na1.set_title(f\"A study holds {per_study.median():.0f} series in the middle case\")\n\na2.bar(planes.index, planes.values, color=[BLUE, ORANGE, AQUA])\nfor i, v in enumerate(planes.values):\n a2.text(i, v + 200, f\"{v:,}\", ha=\"center\", fontsize=8)\na2.set_ylim(0, planes.max() * 1.18); a2.set_ylabel(\"series\")\na2.set_title(\"Sagittal is the most common plane\")\n\nshare = (plane_per_study > 0).mean().sort_values(ascending=False)\na3.barh(share.index, share.values, color=BLUE)\nfor i, v in enumerate(share.values):\n a3.text(v - 0.06, i, f\"{v:.1%}\", va=\"center\", ha=\"right\", color=\"white\", fontsize=8)\na3.set_xlim(0, 1.05); a3.set_xlabel(\"share of studies with at least one\")\na3.set_title(\"Every study has a sagittal series\")\na3.grid(axis=\"y\", visible=False)\n\nplt.tight_layout(); plt.show()\n\nfacts([(\"fewest series in a study\", int(per_study.min())),\n (\"middle\", int(per_study.median())),\n (\"most\", int(per_study.max())),\n (\"studies with no sagittal series\", int((plane_per_study.get(\"Sagittal\", 0) == 0).sum()))],\n \"Series per study\")"},{"cell_type":"markdown","id":"a9954fab","metadata":{},"source":"### The two scanner flags carry one fact, not two\n\n`train_series.csv` gives each series two flags. `Fluid_Sensitive` says fluid appears bright in\nthe image, which is what makes an effusion easy to see. `Fat_Suppression` says the scanner was\nset up to darken fat, which is what makes swelling in bone marrow stand out. These are two\ndifferent scanner settings and in general they vary on their own.\n\nIn this file they do not. Check whether the two columns ever disagree."},{"cell_type":"code","id":"bd1f743f","metadata":{},"execution_count":null,"outputs":[],"source":"agree = (train_series.Fluid_Sensitive == train_series.Fat_Suppression).mean()\nnote(f\"The two flags hold the same value on **{agree:.2%}** of the {len(train_series):,} \"\n f\"series rows.\")\nshow(pd.crosstab(train_series.Fluid_Sensitive, train_series.Fat_Suppression),\n \"Fluid_Sensitive against Fat_Suppression, series counts\")\nshow(pd.crosstab(train_series.Anatomical_Plane, train_series.Fluid_Sensitive),\n \"Anatomical_Plane against Fluid_Sensitive, series counts\")"},{"cell_type":"markdown","id":"67daf381","metadata":{},"source":"The two columns are equal on every one of the 24,371 rows. So the file gives you one piece of\ninformation about the scanner settings, under two names. If you want the two properties apart\nyou have to read them out of the DICOM headers, which section 6 looks at.\n\nFor a first model this is enough. It says which series show fluid brightly, and that is the\nseries to feed the model, because five of the twelve findings are about fluid or swelling."},{"cell_type":"markdown","id":"6f5d27d0","metadata":{},"source":"## 3. The 58 studies that carry labels\n\n58 studies have all twelve labels filled in by a radiologist looking at the images. The other\n4,349 have every label blank. Look at what those 58 contain."},{"cell_type":"code","id":"69a70b92","metadata":{},"execution_count":null,"outputs":[],"source":"prev = gold.mean().sort_values(ascending=False)\nco = pd.DataFrame(index=LABELS, columns=LABELS, dtype=float)\nfor a in LABELS:\n for b in LABELS:\n co.loc[a, b] = gold.loc[gold[a] == 1, b].mean() if gold[a].sum() else np.nan\n\nfig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.6), gridspec_kw={\"width_ratios\": [1, 1.15]})\n\na1.barh(prev.index[::-1], prev.values[::-1], color=BLUE)\nfor i, (c, v) in enumerate(zip(prev.index[::-1], prev.values[::-1])):\n a1.text(v + 0.012, i, f\"{v:.0%} ({int(gold[c].sum())})\", va=\"center\", fontsize=7.5)\na1.set_xlim(0, 0.78); a1.set_xlabel(\"share of the 58 studies with the finding\")\na1.set_title(\"Every finding is common in the 58,\\nand the rarest has 9 positives\")\na1.grid(axis=\"y\", visible=False)\n\norder = prev.index.tolist()\nim = a2.imshow(co.loc[order, order].values.astype(float), cmap=SEQ, vmin=0, vmax=1)\na2.set_xticks(range(12)); a2.set_xticklabels(order, rotation=90, fontsize=7)\na2.set_yticks(range(12)); a2.set_yticklabels(order, fontsize=7)\na2.set_title(\"Given the finding in the row,\\nhow often the column is also present\")\na2.grid(False)\nplt.colorbar(im, ax=a2, shrink=0.85)\nplt.tight_layout(); plt.show()\n\nfacts([(\"labelled studies\", len(gold)),\n (\"studies with no finding at all\", int((gold.sum(axis=1) == 0).sum())),\n (\"findings per study, average\", f\"{gold.sum(axis=1).mean():.1f}\"),\n (\"findings per study, most\", int(gold.sum(axis=1).max())),\n (\"rarest finding\", f\"{gold.sum().idxmin()}, {int(gold.sum().min())} positives\")],\n \"The 58 labelled studies\")"},{"cell_type":"markdown","id":"b7546274","metadata":{},"source":"These 58 studies are sick knees. Every finding is present in at least 15% of them, and the\naverage study has four of the twelve findings. That is not what a random hospital week looks\nlike, so treat these prevalences as a property of this subset and not of the test set.\n\n### What 58 studies can and cannot measure\n\nIt is tempting to use these 58 studies to choose between two ideas. Measure how much a number\nfrom 58 studies can move on its own before you do. The check below resamples the 58 studies\nwith replacement 2,000 times and reads the AUC again each time."},{"cell_type":"code","id":"1828eecc","metadata":{},"execution_count":null,"outputs":[],"source":"def bootstrap_auc(y: pd.DataFrame, p: pd.DataFrame, n_boot=2000, seed=0):\n \"\"\"95% interval for each finding's AUC and for the average, from resampling the rows.\"\"\"\n from sklearn.metrics import roc_auc_score\n rng = np.random.default_rng(seed)\n per, mean = {c: [] for c in LABELS}, []\n for _ in range(n_boot):\n i = rng.integers(0, len(y), len(y))\n vals = []\n for c in LABELS:\n yc = y[c].values[i]\n if yc.min() == yc.max():\n continue\n a = roc_auc_score(yc, p[c].values[i])\n per[c].append(a); vals.append(a)\n if len(vals) == len(LABELS):\n mean.append(np.mean(vals))\n width = {c: float(np.diff(np.percentile(v, [2.5, 97.5]))[0]) for c, v in per.items()}\n return width, float(np.diff(np.percentile(mean, [2.5, 97.5]))[0])\n\n# A stand in for a real model: the annotation itself with noise added, so the check measures the\n# sample size rather than any particular prediction.\nrng = np.random.default_rng(SEED)\nnoisy = gold + rng.normal(0, 0.9, gold.shape)\nw_per, w_mean = bootstrap_auc(gold, pd.DataFrame(noisy, index=gold.index, columns=LABELS))\n\nshow(pd.DataFrame({\"95% interval width\": pd.Series(w_per)})\n .sort_values(\"95% interval width\", ascending=False),\n \"How far one finding's AUC moves on 58 studies, from resampling the rows alone\",\n bars=[\"95% interval width\"], vmin=0, vmax=0.4)\nfacts([(\"widest single finding\", round(max(w_per.values()), 3)),\n (\"the average over twelve findings\", round(w_mean, 3))],\n \"The same measurement, summarised\")"},{"cell_type":"markdown","id":"c32d4715","metadata":{},"source":"One finding's AUC on 58 studies moves by 0.2 to 0.3 from resampling alone. The average over\ntwelve findings is steadier, at about 0.07, because averaging cancels some of the noise.\n\nThe rule that follows is worth keeping. Use the 58 studies to check that a model is not broken,\nbecause a mean AUC near 0.5 there means something is wrong. Do not use them to choose between\ntwo models that are 0.02 apart, because 58 studies cannot see a gap that small. For choosing,\nyou need a target that covers all 4,407 studies, and section 5 builds one."},{"cell_type":"markdown","id":"85f71c1b","metadata":{},"source":"## 4. The reports\n\nEvery training study carries the report the radiologist wrote. This is the only description of\n4,349 studies that you have, so it is worth knowing what it looks like."},{"cell_type":"code","id":"02f8aa0c","metadata":{},"execution_count":null,"outputs":[],"source":"HEADERS = {\n \"en\": [\"findings\", \"impression\", \"technique\", \"clinical\", \"conclusion:\", \"comparison\", \"history\"],\n \"es\": [\"tecnica\", \"resultados\", \"impresion\", \"hallazgos\", \"antecedentes\", \"estudio\"],\n \"nl\": [\"bevindingen\", \"klinische\", \"inlichtingen\", \"conclusie\", \"verslag\", \"vraagstelling\"],\n \"de\": [\"befund\", \"beurteilung\", \"fragestellung\", \"anamnese\", \"technik\"],\n \"fr\": [\"indication\", \"resultats\", \"technique:\", \"conclusion :\", \"protocole\"],\n \"tr\": [\"bulgular\", \"sonuc\", \"inceleme\", \"tetkik\", \"klinik bilgi\", \"oyku\"],\n \"hr\": [\"nalaz\", \"zaključak\", \"misljenje\", \"mr nalaz\", \"opis\"],\n}\nWORDS = {\n \"en\": [\"the\", \"and\", \"with\", \"there\", \"is\", \"of\", \"no \", \"are\", \"was\"],\n \"es\": [\"de la\", \"del\", \"con\", \"sin\", \"en el\", \"se \", \"los\", \"las\", \"una\"],\n \"nl\": [\"van\", \"met\", \"het\", \"een\", \"geen\", \"is er\", \"bij\", \"naar\"],\n \"de\": [\"und\", \"mit\", \"des\", \"der\", \"kein\", \"eine\", \"im \", \"zeigt\"],\n \"fr\": [\"avec\", \"sans\", \"du \", \"des \", \"le \", \"la \", \"une\", \"est\"],\n \"tr\": [\"ve \", \"ile\", \"sag\", \"sol\", \"var\", \"izlen\", \"mm \", \"olan\"],\n \"hr\": [\"je \", \"sa \", \"uz \", \"nema\", \"vidljiv\", \"desno\", \"lijevo\", \"prikaz\"],\n}\nDIACRITICS = {\"tr\": \"ığşİĞŞ\", \"hr\": \"čćžšđČĆŽŠĐ\", \"es\": \"ñáíóúÁÍÓÚ\",\n \"de\": \"äöüßÄÖÜ\", \"fr\": \"éèêçàùÉÈÇ\", \"nl\": \"ijëï\"}\n\ndef fold(s: str) -> str:\n \"\"\"Lower case, strip accents, and map the Turkish dotless i onto a plain i.\n\n The dotless i has no accent to strip, so without this line every Turkish word containing it\n fails to match a pattern written with a plain i.\n \"\"\"\n s = str(s).replace(\"ı\", \"i\").replace(\"İ\", \"i\").replace(\"đ\", \"d\").replace(\"Đ\", \"d\")\n s = unicodedata.normalize(\"NFKD\", s)\n return \"\".join(c for c in s if not unicodedata.combining(c)).lower()\n\ndef detect_language(text: str) -> str:\n \"\"\"Greek and Bulgarian are decided by the alphabet. The rest are scored on section\n headers first, then on common words, then on which accented letters appear.\"\"\"\n raw = str(text)\n letters = sum(c.isalpha() for c in raw) or 1\n if sum(\"Ͱ\" <= c <= \"Ͽ\" or \"ἀ\" <= c <= \"῿\" for c in raw) / letters > 0.30:\n return \"el\"\n if sum(\"Ѐ\" <= c <= \"ӿ\" for c in raw) / letters > 0.30:\n return \"bg\"\n f = fold(raw)\n score = {}\n for lg in WORDS:\n s = 6.0 * sum(h in f for h in HEADERS.get(lg, ()))\n s += sum(f.count(w) for w in WORDS[lg])\n s += 2.0 * sum(ch in raw for ch in DIACRITICS.get(lg, \"\"))\n score[lg] = s\n best = max(score, key=score.get)\n return best if score[best] > 0 else \"unk\"\n\nNAMES = {\"en\": \"English\", \"es\": \"Spanish\", \"tr\": \"Turkish\", \"hr\": \"Croatian\", \"el\": \"Greek\",\n \"de\": \"German\", \"bg\": \"Bulgarian\", \"nl\": \"Dutch\", \"fr\": \"French\", \"unk\": \"unknown\"}\nlang = reports.map(detect_language)\ncounts = lang.value_counts()\n\nfacts([(\"reports\", len(reports)),\n (\"distinct texts\", int(reports.nunique())),\n (\"studies sharing a report with another\", int(len(reports) - reports.nunique())),\n (\"largest group sharing one report\", int(reports.value_counts().iloc[0])),\n (\"shortest report, characters\", int(reports.str.len().min())),\n (\"middle\", int(reports.str.len().median())),\n (\"longest\", int(reports.str.len().max())),\n (\"not written in English\", f\"{(lang != 'en').mean():.1%}\")],\n \"The reports\")"},{"cell_type":"code","id":"c67dc7dc","metadata":{},"execution_count":null,"outputs":[],"source":"fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.1), gridspec_kw={\"width_ratios\": [1.25, 1]})\n\nnames = [NAMES[k] for k in counts.index]\na1.bar(names, counts.values, color=[BLUE if k == \"en\" else ORANGE for k in counts.index])\nfor i, v in enumerate(counts.values):\n a1.text(i, v + 25, f\"{v:,}\", ha=\"center\", fontsize=7.5)\na1.set_ylim(0, counts.max() * 1.15); a1.set_ylabel(\"studies\")\na1.set_title(\"The reports are written in nine languages\")\na1.tick_params(axis=\"x\", rotation=45)\nfor t in a1.get_xticklabels():\n t.set_ha(\"right\")\n\na2.hist(reports.str.len(), bins=50, color=BLUE)\na2.axvline(reports.str.len().median(), color=INK, ls=\":\", lw=1.2)\na2.text(reports.str.len().median() + 90, a2.get_ylim()[1] * 0.9,\n f\"median {int(reports.str.len().median())}\", fontsize=8)\na2.set_xlabel(\"characters in the report\"); a2.set_ylabel(\"studies\")\na2.set_title(\"Most reports are a few short paragraphs\")\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"fbb97610","metadata":{},"source":"Two facts here change what you do next.\n\n60% of the reports are not in English, and two of the nine languages are not written in the\nLatin alphabet. Any word list you write for reading these reports has to cover all nine or it\nwill quietly return \"nothing found\" for whole languages. Section 5 measures that.\n\n131 studies share a report with another study. The largest group is 37 studies with one\nidentical text, which is a template a radiologist reuses for a knee with nothing wrong. If the\ntraining labels come from the report, then every study in such a group gets the same label\nvector. Splitting that group across a training and validation divide would let the model see\nthe answer for a validation study during training. So the split below is made on a hash of the\nreport text, which keeps each group whole."},{"cell_type":"markdown","id":"aff09d94","metadata":{},"source":"## 5. Where the training labels come from\n\n58 labelled studies are not enough to train an image model, and the reports describe all 4,407.\nSo the reports have to be turned into labels. There are two ways to do it and this section\nmeasures both, because the quality of these labels is the ceiling on everything after.\n\n### A word list, written here\n\nThe first way is to look for words. For each finding, list the words a radiologist uses for the\nbody part and the words used for the problem, in all nine languages. Split the report into\nclauses, and for each clause ask whether a body part word and a problem word appear close\ntogether. If they do, the finding is present. If the clause also contains a word of denial such\nas \"no\" or \"sin\" or \"geen\", the finding is absent.\n\nFive findings need no body part word, because the word for the problem names the location\nalready. There is only one place a Baker's cyst can be."},{"cell_type":"code","id":"f40c390e","metadata":{},"execution_count":null,"outputs":[],"source":"def rx(*pats): return re.compile(\"|\".join(pats))\n\nA_MED_MEN = rx(r\"menisc\\w* (?:interno|medial|mediale)\", r\"medial menisc\\w*\", r\"mediale menisc\\w*\",\n r\"innenmeniskus\", r\"menisque (?:interne|medial)\", r\"medial menisk\\w*\", r\"ic menisk\",\n r\"unutrasnj\\w* menisk\", r\"esw menisk\", r\"εσω μηνισκ\", r\"μηνισκ\\w* εσω\",\n r\"вътрешния менис\", r\"медиалния менис\")\nA_LAT_MEN = rx(r\"menisc\\w* (?:externo|lateral|laterale)\", r\"lateral menisc\\w*\", r\"laterale menisc\\w*\",\n r\"aussenmeniskus\", r\"menisque (?:externe|lateral)\", r\"lateral menisk\\w*\", r\"dis menisk\",\n r\"vanjsk\\w* menisk\", r\"εξω μηνισκ\", r\"μηνισκ\\w* εξω\",\n r\"външния менис\", r\"латералния менис\")\nA_ACL = rx(r\"\\bacl\\b\", r\"anterior cruciate\", r\"cruzado anterior\", r\"\\blca\\b\", r\"voorste kruisband\",\n r\"\\bvkb\\b\", r\"vorder\\w* kreuzband\", r\"croise anterieur\", r\"on capraz\", r\"\\bocb\\b\",\n r\"prednj\\w* ukrizen\", r\"προσθι\\w* χιαστ\", r\"предна кръстн\")\nA_MCL = rx(r\"\\bmcl\\b\", r\"medial collateral\", r\"colateral (?:medial|interno)\", r\"innenband\",\n r\"mediale collaterale\", r\"\\blcm\\b\", r\"ic yan bag\", r\"medial yan bag\",\n r\"ligament collateral (?:medial|interne)\", r\"medijaln\\w* kolateral\",\n r\"εσω πλαγι\", r\"вътрешния колатерал\", r\"медиалния колатерал\")\nA_MEDC = rx(r\"\\bmedial\\b\", r\"\\bmediale\\b\", r\"\\binterno\\b\", r\"\\binterna\\b\", r\"\\binnen\", r\"\\bic\\b\",\n r\"medijaln\", r\"unutrasnj\", r\"εσω\", r\"медиал\", r\"вътрешн\")\nA_LATC = rx(r\"\\blateral\\b\", r\"\\blaterale\\b\", r\"\\bexterno\\b\", r\"\\bexterna\\b\", r\"\\baussen\",\n r\"\\bdis\\b\", r\"lateraln\", r\"vanjsk\", r\"εξω\", r\"латерал\", r\"външн\")\nA_PF = rx(r\"patell?o-?femoral\", r\"femoro-?patell\", r\"femoropatelar\", r\"retropatell\",\n r\"femoro-?patellaire\", r\"patellar cartilag\", r\"cartilag\\w* (?:of the )?patell\",\n r\"trochlea\", r\"troklea\", r\"επιγονατιδομηριαι\", r\"пателофеморал\")\n\nP_TEAR = rx(r\"\\btear\", r\"\\btorn\\b\", r\"ruptur\", r\"rotura\", r\"\\bdesgarr\", r\"disrupt\", r\"scheur\",\n r\"\\briss\", r\"yirtik\", r\"yirtil\", r\"macerat\", r\"puknu\", r\"lezij\\w* menisk\",\n r\"ρηξη\", r\"ρηγμα\", r\"разкъс\", r\"руптур\")\nP_OA = rx(r\"osteoarthrit\", r\"arthros\", r\"artros\", r\"gonarthros\", r\"gonartros\", r\"degenerat\",\n r\"dejenerat\", r\"joint space narrowing\", r\"pincement\", r\"osteofit\", r\"osteophyt\",\n r\"chondral thinning\", r\"cartilage (?:loss|thinning)\", r\"kraakbeen\", r\"knorpel\",\n r\"cartilag\", r\"kikirdak\", r\"hrskavic\", r\"οστεοαρθρ\", r\"αρθριτιδ\", r\"χονδροπαθ\",\n r\"артроз\", r\"остеоартр\", r\"хондропат\")\nP_EFF = rx(r\"effusion\", r\"derrame\", r\"\\bhydrops\\b\", r\"gelenkerguss\", r\"\\berguss\\b\", r\"epanchement\",\n r\"versamento\", r\"joint fluid\", r\"efuzyon\", r\"eklem (?:ici )?sivi\", r\"artikuler sivi\",\n r\"vocht in het gewricht\", r\"intra-?articular fluid\", r\"izliv\",\n r\"συλλογη (?:υγρου|αρθρικου)\", r\"αρθρικη συλλογη\", r\"изли[вя]\", r\"хидропс\")\nP_SYN = rx(r\"synovit\", r\"sinovit\", r\"synovial (?:thickening|proliferation|hypertroph)\", r\"sinovyal\",\n r\"sinovijalni\", r\"υμενιτιδ\", r\"αρθρικου υμενα\", r\"синовит\")\nP_BAK = rx(r\"\\bbaker\", r\"popliteal cyst\", r\"poplitea\\w* cyst\", r\"bakerzyste\", r\"quiste de baker\",\n r\"kyste de baker\", r\"cisti di baker\", r\"poplitealna cist\", r\"popliteal kist\",\n r\"κυστη baker\", r\"киста на бейкър\", r\"поплитеална кист\")\nP_CON = rx(r\"contusion\", r\"contusio\", r\"bone (?:bruise|contusion)\", r\"edema (?:oseo|osseo|medular)\",\n r\"bone marrow edema\", r\"beenmerg\\w*oedeem\", r\"knochenmarkodem\", r\"kneuzing\",\n r\"oedeme osseux\", r\"kemik kontuzyon\", r\"kontuzyon\", r\"medullar\\w* odem\",\n r\"οιδημα (?:μυελου|οστου)\", r\"костномозъчен оток\", r\"едем на костния\")\nP_FRA = rx(r\"fractur\", r\"fractuur\", r\"fraktur\", r\"\\bfissur\", r\"frattura\", r\"kirik\", r\"prijelom\",\n r\"breuk\", r\"καταγμα\", r\"фрактур\", r\"счупван\")\n\nNEG = rx(r\"\\bno\\b\", r\"\\bnot\\b\", r\"\\bwithout\\b\", r\"\\bnegative\\b\", r\"\\bintact\\b\", r\"\\bnormal\\b\",\n r\"\\bunremarkable\\b\", r\"\\bpreserved\\b\", r\"\\bsin\\b\", r\"\\bausencia\\b\", r\"\\bausente\",\n r\"\\bgeen\\b\", r\"\\bkein\", r\"\\bpas de\\b\", r\"\\bsans\\b\", r\"\\bnessun\", r\"\\byok\\b\", r\"\\bnema\\b\",\n r\"\\bbez\\b\", r\"δεν \", r\"ουδεν\", r\"χωρις\", r\"без \", r\"липсва\")\nSEV_LOW = rx(r\"\\btrace\\b\", r\"\\bsmall\\b\", r\"\\bminimal\\b\", r\"\\bmild\\b\", r\"\\bslight\\b\", r\"\\bscant\\b\",\n r\"\\bleve\\b\", r\"\\bminim\", r\"\\bpeque\", r\"\\bhafif\\b\", r\"\\bgering\", r\"\\bdiscret\",\n r\"\\bblag\", r\"\\bmicro\")\nSEV_HIGH = rx(r\"\\bmoderate\\b\", r\"\\blarge\\b\", r\"\\bmassive\\b\", r\"\\bmarked\\b\", r\"\\bgross\\b\",\n r\"\\bsignificant\\b\", r\"\\bextensive\\b\", r\"\\bmoderado\\b\", r\"\\babundante\\b\",\n r\"\\bimportante\\b\", r\"\\bbelirgin\\b\", r\"\\bausgepragt\", r\"\\bveel\\b\", r\"\\bizrazit\")\n\nCLAUSE = re.compile(r\"[.;:\\n]+\")\nPAIRED = {\"ACL\": (A_ACL, P_TEAR), \"MCL\": (A_MCL, P_TEAR),\n \"Medial Meniscus\": (A_MED_MEN, P_TEAR), \"Lateral Meniscus\": (A_LAT_MEN, P_TEAR),\n \"Medial OA\": (A_MEDC, P_OA), \"Lateral OA\": (A_LATC, P_OA), \"PF OA\": (A_PF, P_OA)}\nSOLO = {\"Effusion\": P_EFF, \"Synovitis\": P_SYN, \"Baker's\": P_BAK,\n \"Contusion\": P_CON, \"Fracture\": P_FRA}\n\ndef read_report(text: str, window: int = 70):\n \"\"\"Return one score per finding, and whether the report mentioned it at all.\n\n A score of 0.9 means a clause asserts the finding, 0.65 that it calls it small, 0.95 that it\n calls it large, and 0.1 that a clause denies it. A finding no clause touched stays at 0.0,\n which is the weakest part of this method and section 9 comes back to it.\n \"\"\"\n f = fold(text)\n score = {k: 0.0 for k in LABELS}\n seen = {k: False for k in LABELS}\n for clause in CLAUSE.split(f):\n if not clause.strip():\n continue\n if NEG.search(clause):\n value = 0.10\n elif SEV_HIGH.search(clause):\n value = 0.95\n elif SEV_LOW.search(clause):\n value = 0.65\n else:\n value = 0.90\n for k, pat in SOLO.items():\n if pat.search(clause):\n seen[k] = True\n score[k] = max(score[k], value)\n for k, (anat, path) in PAIRED.items():\n for m in anat.finditer(clause):\n lo, hi = max(0, m.start() - window), min(len(clause), m.end() + window)\n if path.search(clause[lo:hi]):\n seen[k] = True\n score[k] = max(score[k], value)\n break\n return score, seen\n\nparsed = [read_report(t) for t in reports]\nwordlist = pd.DataFrame([s for s, _ in parsed], index=reports.index)[LABELS]\nmentioned = pd.DataFrame([m for _, m in parsed], index=reports.index)[LABELS]\nlog(f\"read {len(reports):,} reports with the word list\")"},{"cell_type":"code","id":"886e3d57","metadata":{},"execution_count":null,"outputs":[],"source":"from sklearn.metrics import roc_auc_score\n\ndef score_against_gold(pred: pd.DataFrame) -> pd.Series:\n p = pred.reindex(gold.index)\n return pd.Series({c: roc_auc_score(gold[c], p[c].fillna(0.0)) for c in LABELS})\n\nwl_auc = score_against_gold(wordlist)\nshow(pd.DataFrame({\"AUC on the 58\": wl_auc, \"share of reports that mentioned it\":\n mentioned.mean()}).sort_values(\"AUC on the 58\", ascending=False),\n f\"The word list against the radiologist, average {wl_auc.mean():.3f}\",\n bars=[\"AUC on the 58\"], vmin=0.4, vmax=1.0)\nnote(f\"The word list found nothing at all in **{int((~mentioned.any(axis=1)).sum()):,}** of \"\n f\"{len(mentioned):,} studies, which then look like healthy knees.\")"},{"cell_type":"code","id":"eb7f5f8e","metadata":{},"execution_count":null,"outputs":[],"source":"cov = (pd.DataFrame({\"lang\": lang, \"fires\": mentioned.any(axis=1)})\n .groupby(\"lang\").agg(n=(\"fires\", \"size\"), fires=(\"fires\", \"mean\"))\n .sort_values(\"fires\"))\nby_lang = mentioned.groupby(lang).mean().loc[cov.index]\n\nfig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.3), gridspec_kw={\"width_ratios\": [1, 1.35]})\n\na1.bar([NAMES[k] for k in cov.index], cov.fires.values,\n color=[ORANGE if v < 0.8 else BLUE for v in cov.fires])\nfor i, (v, n) in enumerate(zip(cov.fires, cov.n)):\n a1.text(i, v + 0.03, f\"{v:.0%}\\nn={n}\", ha=\"center\", fontsize=7)\na1.set_ylim(0, 1.2); a1.set_ylabel(\"reports where anything at all fired\")\na1.set_title(\"Asking only whether something fired\")\na1.tick_params(axis=\"x\", rotation=45)\nfor t in a1.get_xticklabels():\n t.set_ha(\"right\")\n\nim = a2.imshow(by_lang.values, cmap=SEQ, vmin=0, vmax=1, aspect=\"auto\")\na2.set_xticks(range(12)); a2.set_xticklabels(LABELS, rotation=90, fontsize=7)\na2.set_yticks(range(len(by_lang))); a2.set_yticklabels([NAMES[k] for k in by_lang.index], fontsize=8)\na2.set_title(\"Asking per finding, which is the useful question\")\na2.grid(False)\nplt.colorbar(im, ax=a2, shrink=0.9, label=\"share of reports where the finding fired\")\nplt.tight_layout(); plt.show()\n\ncov_named = cov.rename(index=NAMES).rename(\n columns={\"n\": \"reports\", \"fires\": \"share where anything fired\"})\nshow(cov_named, \"Coverage as one number per language\",\n bars=[\"share where anything fired\"], vmin=0, vmax=1)\nshow(by_lang.rename(index=NAMES),\n \"Coverage per finding, which is where the gaps actually are\",\n grad=list(by_lang.columns), vmin=0, vmax=1)"},{"cell_type":"markdown","id":"b9f015ba","metadata":{},"source":"The left chart is the measure most people reach for, and it is the wrong one. It says Bulgarian\nis fine at 97%. The right chart shows what happens in Bulgarian. One term fires, the one for an\neffusion, and the other eleven findings are near zero. So a Bulgarian report comes out as a knee\nwith fluid in it and nothing else, every time.\n\nTurkish is the weakest language by the overall measure, at 60%, and 217 Turkish studies get\nnothing at all. Turkish builds words by adding endings, so the word for a torn meniscus is one\nlong word rather than two, and a fixed list of patterns misses most of the forms it takes.\n\nThe same reading applies elsewhere in that chart. Croatian never fires on either meniscus and\nGerman never fires on the lateral one. Neither is a fact about those patients. Both are terms\nmissing from the list.\n\nA study where nothing fires gets twelve zeros, and the model then learns it as a healthy knee.\nThis is the part that makes a word list expensive. It does not fail loudly. It reports that a\nwhole language of patients has nothing wrong with them.\n\nNotice also that PF OA came out below 0.5, which is worse than guessing. Osteoarthritis behind\nthe kneecap is graded from how thin the cartilage looks, and a word list cannot grade anything.\nIt matches the word \"patellofemoral\" whether the sentence describes damage or only names the\nimages that were taken.\n\n### A language model reading the same reports\n\nThe second way is to have a language model read each report and answer the twelve questions.\nSeveral people have done this and published the result as a Kaggle dataset. This notebook\nattaches one of them, from `flight0234/rsna-knee-hybrid-report-labels`, and measures it the\nsame way."},{"cell_type":"code","id":"48510cba","metadata":{},"execution_count":null,"outputs":[],"source":"SEARCH = [p for p in map(pathlib.Path, [\"/kaggle/input/datasets\", \"/kaggle/input\", \"/tmp/rsnallm\"])\n if p.exists()]\n\ndef find_file(name: str):\n \"\"\"Find an attached dataset file without walking the competition folder.\n\n rglob under /kaggle/input would descend into the competition mount and its ~700,000 DICOM\n files, which costs minutes for every lookup.\n \"\"\"\n for root in SEARCH:\n for pat in (name, f\"*/{name}\", f\"*/*/{name}\", f\"*/*/*/{name}\"):\n for hit in root.glob(pat):\n if \"competitions\" not in hit.parts:\n return hit\n return None\n\nPUBLISHED = {\n \"language model\": \"report_labels_v4hybrid.csv\", # flight0234/rsna-knee-hybrid-report-labels\n \"four sets merged\": \"report_labels_v5.csv\", # yunusgmsoy/rsna-knee-llm-labels-4-source-merged\n}\npublished, loaded = {}, []\nfor name, fname in PUBLISHED.items():\n hit = find_file(fname)\n if hit is None:\n loaded.append({\"labels\": name, \"file\": fname, \"found\": False, \"studies\": 0})\n continue\n d = pd.read_csv(hit).drop_duplicates(\"StudyInstanceUID\").set_index(\"StudyInstanceUID\")\n ok = all(c in d.columns for c in LABELS)\n if ok:\n published[name] = d[LABELS].astype(float)\n loaded.append({\"labels\": name, \"file\": hit.name, \"found\": True,\n \"studies\": int(len(d)) if ok else 0})\nshow(pd.DataFrame(loaded), \"What is attached\", hide_index=True)\n\nrows = [{\"labels\": \"word list, written here\", \"mean AUC on the 58\": wl_auc.mean(),\n \"studies covered\": len(wordlist), \"distinct values\": int(wordlist.stack().nunique())}]\nfor name, d in published.items():\n rows.append({\"labels\": name, \"mean AUC on the 58\": score_against_gold(d).mean(),\n \"studies covered\": len(d), \"distinct values\": int(d.stack().nunique())})\nshow(pd.DataFrame(rows).sort_values(\"mean AUC on the 58\", ascending=False),\n \"Three ways to label 4,407 studies from 4,407 reports\",\n bars=[\"mean AUC on the 58\"], vmin=0.5, vmax=1.0, hide_index=True)"},{"cell_type":"markdown","id":"d1939d7f","metadata":{},"source":"### Check that a published label set does not contain the answers\n\nOne of those two sets scores exactly 1.000. A label set built only from report text has no way\nto reproduce a radiologist's reading of the images perfectly, so a perfect score means the 58\nanswers were copied into the file. That is a reasonable thing to do when you train, because\nreal labels are better than derived ones. It stops being reasonable the moment you use those\nsame 58 studies to measure anything, which is the only thing they are otherwise for.\n\nThe check needs two parts. On the 58 rows, is every value exactly 0.0 or 1.0 and equal to the\nannotation? On the other 4,349 rows, is that not true? The second part is what separates a\ncopied answer key from a label set that rounds all of its output to 0 and 1 everywhere."},{"cell_type":"code","id":"19b2227c","metadata":{},"execution_count":null,"outputs":[],"source":"def contains_the_answers(labels: pd.DataFrame) -> dict:\n on = labels.reindex(gold.index)\n off = labels.loc[~labels.index.isin(gold.index)]\n exact_on = float(np.isin(on.values, [0.0, 1.0]).mean())\n exact_off = float(np.isin(off.values, [0.0, 1.0]).mean())\n match = float((on.values == gold.values).mean())\n return {\"exactly 0 or 1 on the 58\": exact_on,\n \"exactly 0 or 1 elsewhere\": exact_off,\n \"equals the annotation\": match,\n \"distinct values elsewhere\": int(len(np.unique(off.values))),\n \"contains the answers\": bool(exact_on > 0.999 and match > 0.999 and exact_off < 0.01)}\n\ncheck = pd.DataFrame({name: contains_the_answers(d) for name, d in published.items()}).T\nshow(check, \"The two part check, run on every attached label set\")\n\nTARGET_NAME = None\nfor name, d in published.items():\n if not contains_the_answers(d)[\"contains the answers\"]:\n TARGET_NAME = name\n break\n\nif TARGET_NAME is None:\n target = wordlist.reindex(train.StudyInstanceUID)\n TARGET_NAME = \"word list, written here\"\n note(\"No clean published set is attached, so the word list supplies the targets instead.\")\nelse:\n target = published[TARGET_NAME].reindex(train.StudyInstanceUID)\n note(f\"Training target: **{TARGET_NAME}**, average AUC on the 58 \"\n f\"**{score_against_gold(published[TARGET_NAME]).mean():.3f}**.\")\n\ntarget = target.fillna(0.0)"},{"cell_type":"markdown","id":"0444b5cb","metadata":{},"source":"The merged set is exactly 0.0 or 1.0 on all 58 rows, matches every annotation, and is exactly\n0.0 or 1.0 on none of the other 4,349 rows, where it takes hundreds of distinct values. The\nanswers were written in. The other set is not like that, so this notebook trains on it.\n\nThe gap between the two honest options is large. A word list reaches about 0.73 on the 58 and a\nlanguage model reading the same text reaches about 0.90. That difference is the single biggest\nlever in this competition, and it is bigger than anything an image model of this size will add\non top."},{"cell_type":"markdown","id":"d3b66871","metadata":{},"source":"## 6. What is inside the DICOM files\n\nNow the images. A DICOM file holds one slice of pixels plus a header describing how it was\ntaken. Read the headers of one series per study for a sample of studies and look at what varies."},{"cell_type":"code","id":"8a319916","metadata":{},"execution_count":null,"outputs":[],"source":"import pydicom\nfrom pydicom.errors import InvalidDicomError\n\nHDR = [\"ImagePositionPatient\", \"ImageOrientationPatient\", \"InstanceNumber\", \"Rows\", \"Columns\",\n \"PixelSpacing\", \"SliceThickness\", \"Laterality\", \"Manufacturer\", \"SeriesDescription\",\n \"RepetitionTime\", \"EchoTime\", \"MagneticFieldStrength\"]\n\ndef read_header(path: pathlib.Path):\n try:\n d = pydicom.dcmread(str(path), stop_before_pixels=True, specific_tags=HDR)\n except (InvalidDicomError, OSError, AttributeError):\n return None\n return d\n\ndef pick_series(series_df: pd.DataFrame, folder: pathlib.Path) -> pd.Series:\n \"\"\"One series per study: sagittal first, fluid sensitive next, then the thickest stack.\n\n Sagittal because both cruciate ligaments and both menisci are seen end on in that plane,\n and fluid sensitive because five of the twelve findings are about fluid or swelling. Where a\n study has several equally good candidates, the one with the most slices wins, since a thicker\n stack samples the joint more finely.\n \"\"\"\n d = series_df.copy()\n d[\"rank\"] = (d.Anatomical_Plane == \"Sagittal\") * 2 + d.Fluid_Sensitive.fillna(0)\n best = d[d[\"rank\"] == d.groupby(\"StudyInstanceUID\")[\"rank\"].transform(\"max\")]\n\n def n_files(row):\n p = folder / row.StudyInstanceUID / row.SeriesInstanceUID\n try:\n return sum(1 for _ in p.iterdir())\n except OSError:\n return 0\n\n with ThreadPoolExecutor(READ_THREADS) as ex:\n best = best.assign(n=list(ex.map(n_files, [r for _, r in best.iterrows()])))\n return (best.sort_values(\"n\", ascending=False)\n .groupby(\"StudyInstanceUID\").SeriesInstanceUID.first())\n\nchosen_train = pick_series(train_series, COMP / \"train_series\")\nlog(f\"chose one series for {len(chosen_train):,} of {len(train):,} training studies\")\n\nsample_ids = list(np.random.default_rng(SEED).permutation(chosen_train.index)[:250])\nrows = []\nfor sid in sample_ids:\n folder = COMP / \"train_series\" / sid / chosen_train[sid]\n files = sorted(folder.glob(\"*.dcm\"))\n if not files:\n continue\n d = read_header(files[len(files) // 2])\n if d is None:\n continue\n ps = getattr(d, \"PixelSpacing\", None)\n lat = getattr(d, \"Laterality\", None)\n ipp = getattr(d, \"ImagePositionPatient\", None)\n rows.append({\n \"study\": sid, \"slices\": len(files),\n \"rows_px\": int(getattr(d, \"Rows\", 0)), \"cols_px\": int(getattr(d, \"Columns\", 0)),\n \"mm_per_pixel\": float(ps[0]) if ps else np.nan,\n \"slice_mm\": float(getattr(d, \"SliceThickness\", np.nan) or np.nan),\n \"laterality\": str(lat) if lat else \"\",\n \"manufacturer\": str(getattr(d, \"Manufacturer\", \"\") or \"\").split()[0][:12].upper(),\n \"tesla\": float(getattr(d, \"MagneticFieldStrength\", np.nan) or np.nan),\n \"x_mm\": float(ipp[0]) if ipp else np.nan,\n })\nhdr = pd.DataFrame(rows)\nlog(f\"read headers for {len(hdr)} studies\")\n\nshow(hdr[[\"slices\", \"rows_px\", \"cols_px\", \"mm_per_pixel\", \"slice_mm\", \"tesla\"]]\n .describe().loc[[\"min\", \"25%\", \"50%\", \"75%\", \"max\"]]\n .rename(index={\"25%\": \"lower quarter\", \"50%\": \"middle\", \"75%\": \"upper quarter\"}),\n f\"What differs between studies, over a random {len(hdr)}\")\nshow(hdr.manufacturer.value_counts().to_frame(\"studies\"), \"Who made the scanner\")\nshow(hdr.laterality.replace(\"\", \"not recorded\").value_counts().to_frame(\"studies\"),\n f\"The Laterality tag, recorded on {(hdr.laterality != '').mean():.0%} of studies\")"},{"cell_type":"code","id":"37151d84","metadata":{},"execution_count":null,"outputs":[],"source":"fig, (a1, a2, a3) = plt.subplots(1, 3, figsize=(11, 2.9))\n\na1.hist(hdr.slices, bins=30, color=BLUE)\na1.set_xlabel(\"slices in the chosen series\"); a1.set_ylabel(\"studies\")\na1.set_title(f\"Median {hdr.slices.median():.0f} slices per series\")\n\na2.hist(hdr.mm_per_pixel.dropna(), bins=30, color=BLUE)\na2.set_xlabel(\"millimetres per pixel\"); a2.set_ylabel(\"studies\")\na2.set_title(\"Physical scale differs between studies\")\n\nlat_by_maker = (hdr.assign(has=hdr.laterality != \"\")\n .groupby(\"manufacturer\").has.agg([\"mean\", \"size\"])\n .sort_values(\"mean\", ascending=False))\na3.barh(lat_by_maker.index, lat_by_maker[\"mean\"], color=BLUE)\nfor i, (v, n) in enumerate(zip(lat_by_maker[\"mean\"], lat_by_maker[\"size\"])):\n a3.text(min(v + 0.03, 0.72), i, f\"{v:.0%} n={n}\", va=\"center\", fontsize=7)\na3.set_xlim(0, 1.05); a3.set_xlabel(\"share with a Laterality tag\")\na3.set_title(\"Whether the side is recorded\\ndepends on the scanner maker\")\na3.grid(axis=\"y\", visible=False)\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"fb035893","metadata":{},"source":"Three things in that output change how the images have to be prepared.\n\nThe number of slices per series is not fixed, so a model needs a rule for picking a fixed\nnumber of them. The millimetres per pixel is not fixed either, so resizing every image to the\nsame number of pixels hands the model knees at different physical scales. A stronger model\nwould crop a fixed number of millimetres instead. This one resizes, and that is one of the\nreasons it stops where it does.\n\nThe `Laterality` tag says whether the left or the right knee was scanned. It is missing for a\nlarge share of studies, and whether it is missing depends on which company made the scanner.\nThat matters more than it sounds, and the next two cells explain why."},{"cell_type":"markdown","id":"c9e98120","metadata":{},"source":"### The file names are not in slice order\n\nA series is a folder of files. The obvious way to read it is to sort the file names. That is\nwrong here, and it fails without any error. The file name is a unique identifier assigned by\nthe scanner, not a position, so sorting by it gives a stack in random anatomical order.\n\nThe true order is in the header. `ImagePositionPatient` gives the position of each slice in\nmillimetres, in a coordinate system fixed to the patient, where x runs from the patient's right\nto the patient's left. `ImageOrientationPatient` gives the two directions of the image plane,\nand the cross product of those two is the direction the stack travels. Projecting the position\nonto that direction gives one number per slice that increases along the stack.\n\nMeasure how far off the file name order is."},{"cell_type":"code","id":"7ceeb23e","metadata":{},"execution_count":null,"outputs":[],"source":"def slice_geometry(path: pathlib.Path):\n \"\"\"Position along the stack, the stack direction, and the patient x coordinate.\"\"\"\n d = read_header(path)\n if d is None:\n return None\n ipp = getattr(d, \"ImagePositionPatient\", None)\n iop = getattr(d, \"ImageOrientationPatient\", None)\n inst = getattr(d, \"InstanceNumber\", None)\n if ipp is None or iop is None or len(iop) != 6:\n return {\"k\": float(inst) if inst is not None else np.nan, \"nx\": np.nan,\n \"x\": np.nan, \"from_geometry\": False}\n p = np.asarray(ipp, dtype=float)\n r = np.asarray(iop, dtype=float)\n n = np.cross(r[:3], r[3:])\n return {\"k\": float(p @ n), \"nx\": float(n[0]), \"x\": float(p[0]), \"from_geometry\": True}\n\nfrom scipy.stats import spearmanr\n\nprobe = []\nfor sid in sample_ids[:40]:\n files = sorted((COMP / \"train_series\" / sid / chosen_train[sid]).glob(\"*.dcm\"))\n if len(files) < 5:\n continue\n with ThreadPoolExecutor(READ_THREADS) as ex:\n geo = list(ex.map(slice_geometry, files))\n g = [x for x in geo if x is not None and np.isfinite(x[\"k\"])]\n if len(g) < 5:\n continue\n ks = [x[\"k\"] for x in g]\n rho = float(spearmanr(np.arange(len(ks)), ks)[0])\n probe.append({\"study\": sid, \"n\": len(ks), \"rho_filename_vs_position\": rho,\n \"geometry_present\": all(x[\"from_geometry\"] for x in g)})\nprobe = pd.DataFrame(probe)\nassert len(probe), \"no series could be read, so the order check cannot run\"\n\nshow(probe.rho_filename_vs_position.describe()\n .to_frame(\"rank correlation\")\n .rename(index={\"25%\": \"lower quarter\", \"50%\": \"middle\", \"75%\": \"upper quarter\"}),\n f\"File name order against position along the stack, over {len(probe)} studies\")\nnote(f\"The header carries both position and orientation on \"\n f\"**{probe.geometry_present.mean():.0%}** of those studies, so the true order is always \"\n f\"recoverable here.\")"},{"cell_type":"code","id":"f615f6ad","metadata":{},"execution_count":null,"outputs":[],"source":"fig, ax = plt.subplots(figsize=(6.4, 2.6))\nax.hist(probe.rho_filename_vs_position, bins=np.linspace(-1, 1, 41), color=ORANGE)\nax.axvline(0, color=INK, ls=\":\", lw=1.2)\nax.set_xlabel(\"rank correlation, file name order against position\")\nax.set_ylabel(\"studies\")\nax.set_title(\"Sorting a series by file name gives an order unrelated to the anatomy\")\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"c8bad89b","metadata":{},"source":"The correlation sits around zero. So sorting by file name and then taking the middle slices\ngives a random handful of cross sections rather than the middle of the joint.\n\n### Putting every knee the same way round\n\nFive of the twelve findings name a side of the knee. Medial means the side facing the other\nleg and lateral means the side facing away. Which side of the picture that falls on depends on\nwhether the left or the right knee was scanned. If you do not correct for it, those five labels\nare being learned from a direction the model cannot see, and a horizontal flip as data\naugmentation would turn a medial tear into a lateral one.\n\nThe correction can be read from the geometry, which every slice carries, rather than from the\n`Laterality` tag, which many studies do not. The knee's x position tells you which leg it is,\nbecause the left leg sits on the positive x side of the body's midline. The x part of the stack\ndirection tells you which way the stack travels. Put together, they say whether the slices run\nfrom medial to lateral or the other way, and the ones that run the wrong way get reversed.\n\nCheck that rule against the `Laterality` tag on the studies that have one."},{"cell_type":"code","id":"d4d91f36","metadata":{},"execution_count":null,"outputs":[],"source":"def series_order(files: list, threads: int = 0):\n \"\"\"Sort a series into medial to lateral order and say which knee it is.\n\n Returns the files in order, plus 'L' or 'R' or '' when the geometry is missing. `threads`\n is 0 when this runs inside a worker thread already, because nesting one pool inside\n another multiplies the thread count.\n \"\"\"\n if threads:\n with ThreadPoolExecutor(threads) as ex:\n geo = list(ex.map(slice_geometry, files))\n else:\n geo = [slice_geometry(f) for f in files]\n pairs = [(g, f) for g, f in zip(geo, files) if g is not None and np.isfinite(g[\"k\"])]\n if not pairs:\n return files, \"\"\n pairs.sort(key=lambda t: t[0][\"k\"])\n xs = [g[\"x\"] for g, _ in pairs if np.isfinite(g[\"x\"])]\n nxs = [g[\"nx\"] for g, _ in pairs if np.isfinite(g[\"nx\"])]\n ordered = [f for _, f in pairs]\n if not xs or not nxs:\n return ordered, \"\"\n side = 1.0 if np.mean(xs) > 0 else -1.0 # positive x is the patient's left\n step = 1.0 if np.mean(nxs) > 0 else -1.0 # which way x moves along the stack\n if side * step < 0:\n ordered = ordered[::-1]\n return ordered, (\"L\" if side > 0 else \"R\")\n\ntag_by_study = hdr.set_index(\"study\").laterality\nrows = []\nfor sid in sample_ids[:120]:\n files = sorted((COMP / \"train_series\" / sid / chosen_train[sid]).glob(\"*.dcm\"))\n if len(files) < 5:\n continue\n _, side = series_order(files, READ_THREADS)\n tag = tag_by_study.get(sid, \"\")\n rows.append({\"study\": sid, \"from_geometry\": side, \"from_tag\": tag})\nsides = pd.DataFrame(rows)\nboth = sides[(sides.from_geometry != \"\") & (sides.from_tag != \"\")]\n\nfacts([(\"studies where the geometry gives a side\", f\"{(sides.from_geometry != '').mean():.0%}\"),\n (\"studies where the tag gives a side\", f\"{(sides.from_tag != '').mean():.0%}\"),\n (\"the two agree\", f\"{(both.from_geometry == both.from_tag).mean():.1%} of {len(both)}\"\n if len(both) else \"no study carries both\")],\n f\"Two ways to find out which knee was scanned, over {len(sides)} studies\")\nshow(pd.crosstab(sides.from_geometry.replace(\"\", \"not found\"),\n sides.from_tag.replace(\"\", \"not recorded\")),\n \"The side from the geometry, in rows, against the side from the tag, in columns\")"},{"cell_type":"markdown","id":"5bb7cebf","metadata":{},"source":"### One study, seen as a stack\n\nBefore building the whole training set, look at one study in the order the code above produces."},{"cell_type":"code","id":"0b129852","metadata":{},"execution_count":null,"outputs":[],"source":"import cv2\n\ndef read_pixels(path, size=SIZE):\n \"\"\"One slice, clipped at the 0.5 and 99.5 percentiles and resized.\"\"\"\n try:\n d = pydicom.dcmread(str(path))\n a = d.pixel_array.astype(np.float32)\n except Exception:\n return np.zeros((size, size), np.uint8)\n slope = float(getattr(d, \"RescaleSlope\", 1) or 1)\n inter = float(getattr(d, \"RescaleIntercept\", 0) or 0)\n a = a * slope + inter\n lo, hi = np.percentile(a, [0.5, 99.5])\n a = np.clip((a - lo) / max(hi - lo, 1e-6), 0, 1)\n return (cv2.resize(a, (size, size), interpolation=cv2.INTER_AREA) * 255).astype(np.uint8)\n\ndemo, demo_files = None, []\nfor sid in sample_ids:\n files = sorted((COMP / \"train_series\" / sid / chosen_train[sid]).glob(\"*.dcm\"))\n if len(files) >= N_SLICES:\n demo, demo_files = sid, files\n break\nassert demo is not None, \"no study had enough slices to show\"\ndemo_ordered, demo_side = series_order(demo_files, READ_THREADS)\npick = np.linspace(len(demo_ordered) * 0.1, len(demo_ordered) * 0.9 - 1,\n N_SLICES).round().astype(int)\ndemo_stack = np.stack([read_pixels(demo_ordered[i]) for i in pick])\n\nfig, axes = plt.subplots(2, 6, figsize=(11, 4))\nfor j, ax in enumerate(axes.ravel()):\n ax.imshow(demo_stack[j], cmap=\"gray\")\n ax.set_title(f\"slice {pick[j]}\", fontsize=7)\n ax.axis(\"off\")\nfig.suptitle(f\"One sagittal series in medial to lateral order, {demo_side} knee, \"\n f\"{len(demo_files)} slices reduced to {N_SLICES}\", y=1.0, fontsize=9.5)\nplt.tight_layout(); plt.show()\n\nquote(reports[demo][:600], f\"The report for that study, {NAMES[detect_language(reports[demo])]}, first 600 characters\")"},{"cell_type":"markdown","id":"6887ea1f","metadata":{},"source":"## 7. Building one input per study\n\nThe model gets one array per study, of twelve slices at 224 by 224. Twelve slices are taken at\neven spacing from the middle 80% of the stack, because the first and last few slices of a knee\nseries are usually outside the joint.\n\nThe whole set is read once and held in memory as 8 bit integers, because reading it again every\nepoch would cost far more than the training itself. Twelve slices at 224 by 224 for 4,407\nstudies is about 2.7 GB, which fits.\n\nEvery stage checks the clock. If reading is slower than expected, the training set gets smaller\nrather than the run failing."},{"cell_type":"code","id":"15ad78ca","metadata":{},"execution_count":null,"outputs":[],"source":"def build_study(args):\n \"\"\"One study, as twelve slices in medial to lateral order.\"\"\"\n folder, sid, series = args\n out = np.zeros((N_SLICES, SIZE, SIZE), np.uint8)\n files = sorted((folder / sid / series).glob(\"*.dcm\"))\n if not files:\n return sid, out, False\n ordered, _ = series_order(files)\n lo, hi = len(ordered) * 0.1, len(ordered) * 0.9 - 1\n idx = np.linspace(max(lo, 0), max(hi, 0), N_SLICES).round().astype(int)\n idx = np.clip(idx, 0, len(ordered) - 1)\n for j, i in enumerate(idx):\n out[j] = read_pixels(ordered[i])\n return sid, out, True\n\ndef build_set(ids, chosen, folder, budget_s, tag):\n X = np.zeros((len(ids), N_SLICES, SIZE, SIZE), np.uint8)\n kept, failed = [], 0\n start = time.time()\n jobs = [(folder, sid, chosen[sid]) for sid in ids]\n with ThreadPoolExecutor(READ_THREADS) as ex:\n for n, (sid, vol, ok) in enumerate(ex.map(build_study, jobs)):\n if ok:\n X[len(kept)] = vol\n kept.append(sid)\n else:\n failed += 1\n if (n + 1) % 250 == 0:\n rate = (n + 1) / (time.time() - start)\n log(f\"{tag}: {n + 1:,}/{len(ids):,} studies, {rate:.1f}/s, \"\n f\"{failed} unreadable\")\n if time.time() - start > budget_s or time.time() > HARD_DEADLINE - 1800:\n log(f\"{tag}: stopping at {n + 1:,} studies to stay inside the time budget\")\n break\n log(f\"{tag}: kept {len(kept):,} studies, {failed} unreadable, \"\n f\"{time.time() - start:.0f}s\")\n return X[:len(kept)], kept\n\nrng = np.random.default_rng(SEED)\ntrain_ids = list(chosen_train.index)\nrng.shuffle(train_ids)\n# Keep every labelled study, so the check in section 9 is always possible.\ntrain_ids = [s for s in gold.index if s in chosen_train.index] + \\\n [s for s in train_ids if s not in set(gold.index)]\n\nXtr, kept_ids = build_set(train_ids, chosen_train, COMP / \"train_series\", PREP_BUDGET_S, \"train\")\ngc.collect()\nfacts([(\"studies read\", len(kept_ids)),\n (\"array shape\", \" by \".join(str(d) for d in Xtr.shape)),\n (\"memory held\", f\"{Xtr.nbytes / 1e9:.2f} GB\"),\n (\"labelled studies inside it\", f\"{len(set(kept_ids) & set(gold.index))} of {len(gold)}\")],\n \"The training array\")"},{"cell_type":"markdown","id":"559b20c8","metadata":{},"source":"### Splitting on the report text\n\nOne fifth of the studies are held out. The split is decided by a hash of the report text so\nthat studies sharing a report stay on the same side, for the reason section 4 gives.\n\nThe 58 labelled studies are held out as well, whichever bucket their hash falls in, and so is\nany study that shares a report with one of them. Without that, the model would be trained on\nthe images of the studies used to check it against the radiologist, and that number would say\nhow well the model remembers rather than how well it reads. Holding out the whole report group\ncosts about 60 training studies out of 4,407, and the count printed below is the check that no\ngroup ended up on both sides."},{"cell_type":"code","id":"56782b30","metadata":{},"execution_count":null,"outputs":[],"source":"def bucket(sid: str) -> int:\n h = hashlib.md5(fold(reports.get(sid, sid)).encode()).hexdigest()\n return int(h[:8], 16) % 100\n\ngold_text = set(reports.reindex(gold.index).map(fold))\nis_val = np.array([bucket(s) < HOLDOUT * 100 or fold(reports.get(s, \"\")) in gold_text\n for s in kept_ids])\nYtr = target.reindex(kept_ids).values.astype(np.float32)\n\ngroups = reports.reindex(kept_ids).value_counts()\nshared = groups[groups > 1].index\nleaked = sum(1 for t in shared\n if len(set(is_val[[i for i, s in enumerate(kept_ids) if reports[s] == t]])) > 1)\nfacts([(\"studies used for training\", int((~is_val).sum())),\n (\"studies held out\", int(is_val.sum())),\n (\"labelled studies used for training\", len(set(np.array(kept_ids)[~is_val]) & set(gold.index))),\n (\"labelled studies held out\", f\"{len(set(np.array(kept_ids)[is_val]) & set(gold.index))}\"\n f\" of {len(gold)}\"),\n (\"report groups split across the divide\", leaked)],\n \"The split, and the two checks on it\")"},{"cell_type":"code","id":"bd0f2243","metadata":{},"execution_count":null,"outputs":[],"source":"fig, ax = plt.subplots(figsize=(7.4, 2.6))\nm = Ytr[~is_val].mean(axis=0)\no = np.argsort(-m)\nax.bar(np.array(LABELS)[o], m[o], color=BLUE)\nfor i, v in enumerate(m[o]):\n ax.text(i, v + 0.012, f\"{v:.2f}\", ha=\"center\", fontsize=7)\nax.set_ylim(0, m.max() * 1.2); ax.set_ylabel(\"average target value\")\nax.set_title(f\"What the model is asked to predict, from the {TARGET_NAME} labels\")\nax.tick_params(axis=\"x\", rotation=45)\nfor t in ax.get_xticklabels():\n t.set_ha(\"right\")\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"e5800ee3","metadata":{},"source":"## 8. A small model\n\nA ResNet18 that was trained on ImageNet, with two changes.\n\nThe first layer normally takes three colour channels. Here it takes twelve, one per slice, so\nthe twelve values at one pixel are the brightness of that spot as the stack passes through the\njoint. The first layer's weights are averaged across the three colours, copied twelve times,\nand divided by four so that the numbers coming out of the layer keep roughly the size the rest\nof the network expects.\n\nThe last layer becomes twelve outputs instead of a thousand, one per finding.\n\nThe training targets are values between 0 and 1 rather than 0 or 1, because they come from a\nlanguage model that reports how sure it is. Binary cross entropy accepts that directly and\ntreats a target of 0.7 as seven parts positive to three parts negative.\n\nThere is no horizontal flip in the augmentation. Mirroring an image of a knee turns it into an\nimage of the other knee, which relabels medial as lateral, so five of the twelve findings would\nbe trained against the wrong answer. Shifting the image and changing its brightness are safe.\n\nThe predictions come from the three epochs that scored best on the held out fifth, combined by\naveraging their ranks. Section 1 said the score reads only the order, so ranks are the right\nthing to average, and three epochs of one run cost nothing extra to keep."},{"cell_type":"code","id":"55cb0899","metadata":{},"execution_count":null,"outputs":[],"source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torchvision\n\ntorch.manual_seed(SEED)\nnp.random.seed(SEED) # batches() draws from the global generator\ndev = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nif dev == \"cuda\":\n # A real matrix multiply, so a broken GPU image fails here rather than in an hour.\n _ = (torch.randn(64, 64, device=\"cuda\") @ torch.randn(64, 64, device=\"cuda\")).sum().item()\nlog(f\"device {dev} {torch.cuda.get_device_name(0) if dev == 'cuda' else ''}\")\n\ndef build_model(n_ch=N_SLICES, n_out=len(LABELS)):\n m = torchvision.models.resnet18(weights=None)\n # find_file is bounded on purpose. An rglob under /kaggle/input would walk the competition\n # mount and its ~700,000 DICOM files.\n w = find_file(\"resnet18.pth\") or find_file(\"resnet18-f37072fd.pth\")\n if w is not None:\n sd = torch.load(str(w), map_location=\"cpu\", weights_only=False)\n sd = sd.get(\"state_dict\", sd)\n m.load_state_dict(sd)\n note(f\"Loaded ImageNet weights from `{w}`.\")\n else:\n note(\"No ImageNet weights are attached, so the network starts from random values.\")\n old = m.conv1.weight.data\n m.conv1 = nn.Conv2d(n_ch, 64, 7, 2, 3, bias=False)\n m.conv1.weight.data = old.mean(1, keepdim=True).repeat(1, n_ch, 1, 1) * (3.0 / n_ch)\n m.fc = nn.Linear(512, n_out)\n return m\n\nmodel = build_model().to(dev)\nfacts([(\"architecture\", \"ResNet18\"),\n (\"input\", f\"{N_SLICES} slices at {SIZE} by {SIZE}\"),\n (\"outputs\", len(LABELS)),\n (\"parameters\", f\"{sum(p.numel() for p in model.parameters()) / 1e6:.1f} million\"),\n (\"device\", torch.cuda.get_device_name(0) if dev == \"cuda\" else \"cpu\")],\n \"The model\")"},{"cell_type":"code","id":"888937ff","metadata":{},"execution_count":null,"outputs":[],"source":"def augment(x: torch.Tensor) -> torch.Tensor:\n \"\"\"Shift by up to 8% of the image, scale the brightness, add a little noise.\"\"\"\n b = x.shape[0]\n dx, dy = (torch.rand(2, b, device=x.device) * 0.16 - 0.08)\n theta = torch.zeros(b, 2, 3, device=x.device)\n theta[:, 0, 0] = 1; theta[:, 1, 1] = 1\n theta[:, 0, 2] = dx; theta[:, 1, 2] = dy\n grid = F.affine_grid(theta, x.shape, align_corners=False)\n x = F.grid_sample(x, grid, align_corners=False, padding_mode=\"border\")\n x = x * (0.9 + 0.2 * torch.rand(b, 1, 1, 1, device=x.device))\n return x + 0.01 * torch.randn_like(x)\n\ndef batches(idx, size, shuffle):\n idx = np.array(idx)\n if shuffle:\n idx = idx[np.random.permutation(len(idx))]\n for i in range(0, len(idx), size):\n yield idx[i:i + size]\n\ndef to_gpu(X, rows):\n return torch.from_numpy(X[rows]).to(dev).float().div_(255).sub_(0.5).div_(0.25)\n\n@torch.no_grad()\ndef predict(X, rows, bs=64):\n model.eval()\n out = []\n for b in batches(rows, bs, False):\n with torch.autocast(\"cuda\", enabled=(dev == \"cuda\")):\n out.append(torch.sigmoid(model(to_gpu(X, b))).float().cpu().numpy())\n return np.concatenate(out) if out else np.zeros((0, len(LABELS)), np.float32)\n\ntr_rows = np.where(~is_val)[0]\nva_rows = np.where(is_val)[0]\nYt = torch.from_numpy(Ytr).to(dev)\n\nhead = [p for n, p in model.named_parameters() if n.startswith((\"fc\", \"conv1\"))]\nbody = [p for n, p in model.named_parameters() if not n.startswith((\"fc\", \"conv1\"))]\nopt = torch.optim.AdamW([{\"params\": head, \"lr\": 6e-4}, {\"params\": body, \"lr\": 1.5e-4}],\n weight_decay=1e-4)\nsteps = max(1, EPOCHS * int(np.ceil(len(tr_rows) / BATCH)))\nsched = torch.optim.lr_scheduler.OneCycleLR(opt, max_lr=[6e-4, 1.5e-4], total_steps=steps,\n pct_start=0.25)\nscaler = torch.amp.GradScaler(\"cuda\", enabled=(dev == \"cuda\"))\n\ndef safe_auc(y, p):\n \"\"\"AUC, or nan when the column has no positives or no negatives to compare.\"\"\"\n y = np.asarray(y)\n if len(y) == 0 or y.min() == y.max():\n return np.nan\n return float(roc_auc_score(y, p))\n\ndef mean_auc(y, p):\n vals = [safe_auc((y[:, j] > 0.5).astype(int), p[:, j]) for j in range(len(LABELS))]\n vals = [v for v in vals if np.isfinite(v)]\n return float(np.mean(vals)) if vals else np.nan\n\nhistory, snapshots = [], []\ndone = 0\nfor epoch in range(EPOCHS):\n model.train()\n total = 0.0\n for b in batches(tr_rows, BATCH, True):\n if done >= steps:\n break\n x = augment(to_gpu(Xtr, b))\n with torch.autocast(\"cuda\", enabled=(dev == \"cuda\")):\n loss = F.binary_cross_entropy_with_logits(model(x), Yt[b])\n opt.zero_grad(set_to_none=True)\n scaler.scale(loss).backward()\n scaler.step(opt); scaler.update(); sched.step()\n total += loss.item() * len(b); done += 1\n pv = predict(Xtr, va_rows)\n va = mean_auc(Ytr[va_rows], pv)\n history.append({\"epoch\": epoch + 1, \"loss\": total / len(tr_rows), \"holdout_auc\": va})\n log(f\"epoch {epoch + 1}/{EPOCHS} loss {total / len(tr_rows):.4f} holdout AUC {va:.4f}\")\n snapshots.append((va, {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}))\n if time.time() > HARD_DEADLINE - 2400:\n log(\"stopping training early to leave time for the test set\")\n break\n\nassert snapshots, \"training did not finish one epoch, so there is nothing to report\"\nhist = pd.DataFrame(history)\n\n# Section 1 said the score reads only the order, so several models should be combined by\n# averaging their ranks rather than their probabilities. The three epochs that scored best on\n# the held out fifth are combined that way here. Averaging probabilities instead would let\n# whichever epoch happens to be most confident decide the answer.\nKEEP = min(3, len(snapshots))\ntop = sorted(snapshots, key=lambda t: -t[0])[:KEEP]\n\ndef predict_top(X, rows):\n \"\"\"Rank average the predictions of the KEEP best epochs.\"\"\"\n acc = None\n for _, sd in top:\n model.load_state_dict(sd)\n p = pd.DataFrame(predict(X, rows), columns=LABELS)\n p = p.rank(pct=True) if len(p) > 1 else p\n acc = p if acc is None else acc + p\n return (acc / len(top)).values\n\nmodel.load_state_dict(top[0][1]) # the CAM section reads one model, so leave the best loaded\nkept_epochs = [h[\"epoch\"] for h in history if h[\"holdout_auc\"] in [t[0] for t in top]]\nshow(hist.assign(**{\"combined by rank average\": hist.epoch.isin(kept_epochs)})\n .rename(columns={\"loss\": \"training loss\", \"holdout_auc\": \"holdout AUC\"}),\n \"Every epoch, and which three were kept\",\n bars=[\"holdout AUC\"], vmin=0.5, vmax=0.85, hide_index=True)\nfacts([(\"best single epoch, holdout AUC\", f\"{top[0][0]:.4f}\"),\n (\"epochs combined by rank average\", \", \".join(str(e) for e in kept_epochs))],\n \"What is used for every prediction below\")"},{"cell_type":"code","id":"a1482390","metadata":{},"execution_count":null,"outputs":[],"source":"fig, (a1, a2) = plt.subplots(1, 2, figsize=(9.5, 2.8))\na1.plot(hist.epoch, hist.loss, color=BLUE, lw=2, marker=\"o\", ms=4)\na1.set_xlabel(\"epoch\"); a1.set_ylabel(\"training loss\")\na1.set_title(\"Training loss\")\na2.plot(hist.epoch, hist.holdout_auc, color=ORANGE, lw=2, marker=\"o\", ms=4)\na2.axhline(0.5, color=MUTED, ls=\":\", lw=1)\na2.text(hist.epoch.iloc[0], 0.508, \"guessing\", color=MUTED, fontsize=7.5)\na2.set_xlabel(\"epoch\"); a2.set_ylabel(\"average AUC on the held out fifth\")\na2.set_title(\"Held out agreement with the report labels\")\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"c0880847","metadata":{},"source":"## 9. What the model learned\n\nThere are two different questions here and they need two different measurements.\n\nThe first is whether the model agrees with the labels it was trained on. That is measured on\nthe held out fifth, it covers hundreds of studies, and it is the number to use when comparing\ntwo versions of the model.\n\nThe second is whether the model agrees with a radiologist looking at the images, which is what\nthe competition scores. That can only be measured on the 58 labelled studies, and none of them\nwere trained on, so it is a fair measurement of the wrong size. Section 3 showed a single\nfinding's AUC on 58 studies moves by 0.2 to 0.3 from resampling alone, so this number says\nwhether the model works at all and not much more."},{"cell_type":"code","id":"6798375b","metadata":{},"execution_count":null,"outputs":[],"source":"val_ids = np.array(kept_ids)[va_rows]\npred_val = pd.DataFrame(predict_top(Xtr, va_rows), index=val_ids, columns=LABELS)\nsingle = pd.DataFrame(predict(Xtr, va_rows), index=val_ids, columns=LABELS)\n\ngold_rows = [i for i, s in enumerate(kept_ids) if s in gold.index]\npred_gold = pd.DataFrame(predict_top(Xtr, np.array(gold_rows)),\n index=[kept_ids[i] for i in gold_rows], columns=LABELS)\ng58 = gold.reindex(pred_gold.index)\n\nholdout_auc = pd.Series({c: safe_auc((Ytr[va_rows][:, j] > 0.5).astype(int), pred_val[c])\n for j, c in enumerate(LABELS)})\ngold_auc = pd.Series({c: safe_auc(g58[c], pred_gold[c]) for c in LABELS})\nlabel_auc = score_against_gold(target)\n\ntable = pd.DataFrame({\n \"positives in holdout\": (Ytr[va_rows] > 0.5).sum(axis=0).astype(int),\n \"model vs report labels\": holdout_auc,\n \"model vs radiologist (58)\": gold_auc,\n \"report labels vs radiologist (58)\": label_auc,\n})\nshow(table.sort_values(\"model vs radiologist (58)\", ascending=False),\n \"Per finding, three AUC values that mean three different things\",\n grad=[\"model vs report labels\", \"model vs radiologist (58)\",\n \"report labels vs radiologist (58)\"], vmin=0.5, vmax=1.0)\nfacts([(\"best single epoch, against the report labels on the holdout\",\n round(mean_auc(Ytr[va_rows], single.values), 3)),\n (\"three epochs combined, same holdout\", round(holdout_auc.mean(), 3)),\n (\"three epochs combined, against the radiologist on the 58\", round(gold_auc.mean(), 3)),\n (\"the report labels themselves, against the radiologist on the 58\",\n round(label_auc.mean(), 3))],\n \"Averages over the twelve findings\")"},{"cell_type":"code","id":"ed7bff89","metadata":{},"execution_count":null,"outputs":[],"source":"o = gold_auc.sort_values(ascending=False).index\nx = np.arange(len(LABELS))\nfig, ax = plt.subplots(figsize=(9.5, 3.2))\nax.bar(x - 0.22, holdout_auc[o], 0.44, color=BLUE, label=\"model against the report labels, holdout\")\nax.bar(x + 0.22, gold_auc[o], 0.44, color=ORANGE, label=\"model against the radiologist, 58 studies\")\nax.plot(x, label_auc[o], \"o\", color=INK, ms=5, label=\"report labels against the radiologist\")\nax.axhline(0.5, color=MUTED, ls=\":\", lw=1)\nax.set_xticks(x); ax.set_xticklabels(o, rotation=45, ha=\"right\")\nax.set_ylim(0.3, 1.0); ax.set_ylabel(\"ROC AUC\")\nax.set_title(\"Where the model works, and what its labels allowed\")\nax.legend(fontsize=7.5, loc=\"lower left\")\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"11b8759a","metadata":{},"source":"Read the two bars together with the black dot. The blue bar says how well the model reproduces\nthe labels it was given. The black dot says how good those labels were in the first place. The\norange bar is the product of the two, and it cannot beat the dot by much for long.\n\nThe submitted version of this notebook scored 0.798 on the public leaderboard against the\n0.778 in the orange bars. The two agreeing is not a given, because the 58 studies and the test\nset are different samples of different sizes, and section 3 says how much room there is for\nthem to disagree.\n\nMCL and Fracture are the two orange bars to be careful with. Both have the fewest positives in\nthe holdout, 130 and 77, and both are the findings whose orange bar swung most between runs of\nthis notebook. Their black dots are high, at 0.968 and 0.870, so the labels are not the limit\nthere. The model is.\n\n### Does the model tell the findings apart, or does it learn one thing twelve times\n\nA study with one finding often has several, as section 3 showed. So a model can score above 0.5\non all twelve columns by learning a single idea of how damaged a knee looks. Compare how much\nthe predictions move together against how much the labels do."},{"cell_type":"code","id":"c89f36eb","metadata":{},"execution_count":null,"outputs":[],"source":"pc = pred_val.corr(method=\"spearman\")\nlc = pd.DataFrame(Ytr[va_rows], columns=LABELS).corr(method=\"spearman\")\n\nfig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4))\nfor ax, m, t in [(a1, lc, \"The training labels\"), (a2, pc, \"The model's predictions\")]:\n im = ax.imshow(m.loc[LABELS, LABELS].values, cmap=DIV, vmin=-1, vmax=1)\n ax.set_xticks(x); ax.set_xticklabels(LABELS, rotation=90, fontsize=7)\n ax.set_yticks(x); ax.set_yticklabels(LABELS, fontsize=7)\n ax.set_title(f\"{t}\\naverage off diagonal correlation \"\n f\"{(m.values.sum() - 12) / (144 - 12):.2f}\")\n ax.grid(False)\n plt.colorbar(im, ax=ax, shrink=0.85)\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"a12ab18f","metadata":{},"source":"## 10. What the model looks at\n\nA number saying the model works is not the same as knowing what it responded to. The check\nbelow is a class activation map. It takes the last block of convolutions, asks how much each of\nits 512 feature maps pushed one output up, and adds the maps together with those amounts as\nweights. The result is a rough picture of which parts of the slice raised that finding's score.\n\nRead it as a hint and not as evidence. The map has the resolution of the last block, which is\n7 by 7 for a 224 pixel input, so it can point at a region and no smaller."},{"cell_type":"code","id":"8eaa8d33","metadata":{},"execution_count":null,"outputs":[],"source":"def class_activation_map(row: int, finding: str):\n \"\"\"Where in the stack the last block of the network pushed one finding up.\"\"\"\n model.eval()\n feats, grads = {}, {}\n h1 = model.layer4.register_forward_hook(lambda m, i, o: feats.__setitem__(\"a\", o))\n h2 = model.layer4.register_full_backward_hook(\n lambda m, gi, go: grads.__setitem__(\"g\", go[0]))\n x = to_gpu(Xtr, np.array([row]))\n out = model(x)\n model.zero_grad()\n out[0, LABELS.index(finding)].backward()\n h1.remove(); h2.remove()\n a, g = feats[\"a\"][0], grads[\"g\"][0]\n cam = F.relu((a * g.mean(dim=(1, 2), keepdim=True)).sum(0))\n cam = cam / (cam.max() + 1e-6)\n cam = cv2.resize(cam.detach().cpu().numpy(), (SIZE, SIZE))\n return cam, float(torch.sigmoid(out[0, LABELS.index(finding)]))\n\nbest_finding = gold_auc.idxmax()\nscores = pred_gold[best_finding].sort_values()\npicks = [(scores.index[-1], \"highest score\"), (scores.index[-2], \"second highest\"),\n (scores.index[1], \"second lowest\"), (scores.index[0], \"lowest score\")]\n\nfig, axes = plt.subplots(2, 4, figsize=(11, 5.6))\nfor col, (sid, tag) in enumerate(picks):\n row = kept_ids.index(sid)\n cam, p = class_activation_map(row, best_finding)\n mid = Xtr[row, N_SLICES // 2]\n axes[0, col].imshow(mid, cmap=\"gray\")\n axes[0, col].set_title(f\"{tag}\\nmodel {p:.2f}, radiologist \"\n f\"{int(gold.loc[sid, best_finding])}\", fontsize=8)\n axes[1, col].imshow(mid, cmap=\"gray\")\n axes[1, col].imshow(cam, cmap=\"inferno\", alpha=0.45)\n axes[1, col].set_title(\"where the score came from\", fontsize=8)\n for r in (0, 1):\n axes[r, col].axis(\"off\")\nfig.suptitle(f\"Middle slice and activation map for {best_finding}, \"\n f\"the finding the model does best on\", y=1.0, fontsize=9.5)\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","id":"c8a1292a","metadata":{},"source":"### How the scores are spread out\n\nThe last check on the predictions is whether they separate at all. A model that returns almost\nthe same number for every study has an AUC near 0.5 no matter how the number was computed."},{"cell_type":"code","id":"7fe3f91f","metadata":{},"execution_count":null,"outputs":[],"source":"fig, axes = plt.subplots(3, 4, figsize=(11, 5.4), sharex=True)\nfor ax, c in zip(axes.ravel(), LABELS):\n y = (Ytr[va_rows][:, LABELS.index(c)] > 0.5)\n ax.hist(pred_val[c][~y], bins=30, color=BLUE, alpha=0.75, label=\"label says absent\")\n ax.hist(pred_val[c][y], bins=30, color=ORANGE, alpha=0.75, label=\"label says present\")\n ax.set_title(f\"{c} AUC {holdout_auc[c]:.2f}\", fontsize=8)\n ax.set_yticks([])\naxes[0, 0].legend(fontsize=6.5)\nfig.suptitle(\"Predicted score on the held out studies, split by the training label\", y=1.0)\nplt.tight_layout(); plt.show()\n\nshow(pred_val.describe().loc[[\"min\", \"50%\", \"max\", \"std\"]].T\n .rename(columns={\"50%\": \"middle\", \"std\": \"spread\"}),\n \"How far apart the predictions are on the held out studies\",\n bars=[\"spread\"], vmin=0, vmax=0.35)"},{"cell_type":"markdown","id":"42b1d1bc","metadata":{},"source":"## 11. Writing the submission\n\nThe test images are read the same way as the training images, the model predicts, and each\ncolumn is replaced by its ranks. In the public run there are only a handful of test studies. In\nthe scored run this cell reads a test set nobody has seen.\n\nAny study the reader could not open, or did not reach before the clock ran out, gets 0.5 in\nevery column. That is the value of a study the model knows nothing about, and it keeps the file\nvalid, which is what the run is graded on."},{"cell_type":"code","id":"1538cafe","metadata":{},"execution_count":null,"outputs":[],"source":"chosen_test = pick_series(test_series, COMP / \"test_series\")\ntest_ids = [s for s in test.StudyInstanceUID if s in chosen_test.index]\nlog(f\"chose a series for {len(test_ids)} of {len(test)} test studies\")\n\nXte, kept_test = build_set(test_ids, chosen_test, COMP / \"test_series\",\n max(600, HARD_DEADLINE - time.time() - 900), \"test\")\n\npred_test = pd.DataFrame(predict_top(Xte, np.arange(len(kept_test))),\n index=kept_test, columns=LABELS)\n\nshow(pred_test.describe().loc[[\"min\", \"50%\", \"max\"]].T.rename(columns={\"50%\": \"middle\"}),\n f\"The model's output on the {len(pred_test)} test studies it could read\")\n\nsub = write_submission(pred_test)\nlog(f\"wrote submission.csv with {len(sub)} rows\")\nshown = sub.head(8).copy()\nshown.insert(0, \"study\", shown.StudyInstanceUID.map(short_uid))\nshow(shown.drop(columns=\"StudyInstanceUID\"),\n \"The first rows of submission.csv, after the ranks are taken\",\n grad=LABELS, vmin=0, vmax=1, hide_index=True)"},{"cell_type":"markdown","id":"2116bfbc","metadata":{},"source":"## 12. What to try next\n\nThe measurements above point at four changes, in the order they are worth making.\n\nBetter labels come first. Section 5 measured a gap of about 0.17 AUC between a word list and a\nlanguage model reading the same reports, and section 9 showed the model cannot beat its labels\nby much. Reading the reports more carefully is worth more than a bigger network.\n\nThen the other two planes. This model reads one sagittal series. Osteoarthritis in the inner\nand outer compartments is read from the coronal images, and the cartilage behind the kneecap is\nread from the axial ones, so three of the twelve findings are being asked for from a view that\nbarely shows them.\n\nThen physical scale. Section 6 showed millimetres per pixel varies between studies. Cropping a\nfixed number of millimetres around the joint and resizing that, instead of resizing the whole\nimage, gives every study the same scale and puts a floor on how small a feature can be and\nstill survive the resize. A meniscal tear is one to three millimetres across.\n\nLast, train once and predict many times. Nothing requires the scored run to be the run that\nlearned the weights. Training in a separate notebook, saving the weights as a dataset, and\nattaching them here would free the whole session for reading and predicting, which is the part\nthat genuinely cannot be done in advance."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11.0"}},"nbformat":4,"nbformat_minor":5}