{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install -q langdetect pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg pydicom iterative-stratification","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:10:28.8434Z","iopub.execute_input":"2026-08-25T17:10:28.843847Z","iopub.status.idle":"2026-08-25T17:10:35.575006Z","shell.execute_reply.started":"2026-08-25T17:10:28.843811Z","shell.execute_reply":"2026-08-25T17:10:35.574229Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile phase2_eda.py","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:11:53.458782Z","iopub.execute_input":"2026-08-25T17:11:53.459341Z","iopub.status.idle":"2026-08-25T17:11:53.46454Z","shell.execute_reply.started":"2026-08-25T17:11:53.459305Z","shell.execute_reply":"2026-08-25T17:11:53.463255Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile phase2_eda.py\n\"\"\"\nPhase 2 EDA — RSNA Knee Abnormality Detection\n=============================================\nProduces every figure and every [FILL: ...] value referenced in\nPhase2_Data_Collection_Preprocessing_EDA.md\n\nRun inside a Kaggle notebook with the competition dataset attached.\n\n    !pip install -q langdetect pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg\n    !python phase2_eda.py\n\nOutputs\n-------\nfigures/*.png        12 figures\neda_summary.json     every FILL value, keyed by its tag\nstudy_features.parquet   engineered study-level feature table (Phase 2 §2.2.5)\n\"\"\"\n\nfrom __future__ import annotations\n\nimport json\nimport os\nimport random\nimport re\nimport unicodedata\nimport warnings\nfrom collections import Counter\nfrom datetime import date\nfrom pathlib import Path\n\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\n\nwarnings.filterwarnings(\"ignore\")\n\n# --------------------------------------------------------------------------\n# Config\n# --------------------------------------------------------------------------\nSEED = 42\nrandom.seed(SEED)\nnp.random.seed(SEED)\n\nROOT = Path(\"/kaggle/input/rsna-knee-abnormality-detection\")\nOUT = Path(\"figures\")\nOUT.mkdir(exist_ok=True)\n\nN_HEADER_SAMPLE = 400      # series to read DICOM headers from\nN_PIXEL_SAMPLE = 40        # series to actually decode pixels for\nMIN_REPORT_CHARS = 20\nMIN_LANG_N = 50            # min studies for a language group to appear in Fig 10\nTARGET_SLICES = 32         # matches the preprocessing decision in §2.2.2\n\nLABELS = [\n    \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n    \"Medial OA\", \"Lateral OA\", \"PF OA\",\n    \"Effusion\", \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\",\n]\n\nSNAKE = {\n    \"ACL\": \"acl\", \"MCL\": \"mcl\",\n    \"Medial Meniscus\": \"medial_meniscus\", \"Lateral Meniscus\": \"lateral_meniscus\",\n    \"Medial OA\": \"medial_oa\", \"Lateral OA\": \"lateral_oa\", \"PF OA\": \"pf_oa\",\n    \"Effusion\": \"effusion\", \"Synovitis\": \"synovitis\", \"Baker's\": \"bakers\",\n    \"Contusion\": \"contusion\", \"Fracture\": \"fracture\",\n}\n\nFILL: dict[str, object] = {\"access_date\": str(date.today())}\n\ndef note(tag: str, value) -> None:\n    \"\"\"Record a [FILL: tag] value and echo it.\"\"\"\n    if isinstance(value, (np.integer,)):\n        value = int(value)\n    elif isinstance(value, (np.floating,)):\n        value = round(float(value), 4)\n    FILL[tag] = value\n    print(f\"  [FILL: {tag}] = {value}\")\n\ndef save(fig, name: str) -> None:\n    fig.tight_layout()\n    fig.savefig(OUT / f\"{name}.png\", dpi=150, bbox_inches=\"tight\")\n    plt.close(fig)\n    print(f\"  -> figures/{name}.png\")\n\n# --------------------------------------------------------------------------\n# 1. Load tabular data\n# --------------------------------------------------------------------------\nprint(\"\\n=== Loading tabular data ===\")\ntrain = pd.read_csv(ROOT / \"train.csv\")\nseries = pd.read_csv(ROOT / \"train_series.csv\")\n\nnote(\"n_studies\", len(train))\nnote(\"n_series\", len(series))\n\ndups = int(train[\"StudyInstanceUID\"].duplicated().sum())\nnote(\"n_dup_studies\", dups)\nnote(\"dup_action\", \"none required\" if dups == 0 else \"kept first occurrence, dropped the rest\")\nif dups:\n    train = train.drop_duplicates(\"StudyInstanceUID\", keep=\"first\")\n\n# Which studies are gold-annotated? A study is labeled if any label column is non-null.\nlabel_present = train[LABELS].notna()\ntrain[\"is_labeled\"] = label_present.all(axis=1)\nn_lab = int(train[\"is_labeled\"].sum())\nnote(\"n_labeled\", n_lab)\nnote(\"n_unlabeled\", len(train) - n_lab)\nnote(\"pct_labeled\", round(100 * n_lab / len(train), 1))\n\n# Is missingness whole-row or per-cell?\npartial = int((label_present.any(axis=1) & ~label_present.all(axis=1)).sum())\nnote(\"label_missing_note\",\n     \"missingness is whole-row (a study is either fully annotated or not at all)\" if partial == 0\n     else f\"{partial} studies have PARTIAL annotation; masked BCE is required, not just row filtering\")\n\nlab = train[train[\"is_labeled\"]].copy()\nuniq = sorted(pd.unique(lab[LABELS].values.ravel()))\nnote(\"label_unique_vals\", str([int(u) for u in uniq]))\nlab[LABELS] = lab[LABELS].astype(\"int8\")\n\n# Orphan / missing cross-references\nstudies_with_series = set(series[\"StudyInstanceUID\"])\nnote(\"n_studies_no_series\", int((~train[\"StudyInstanceUID\"].isin(studies_with_series)).sum()))\nnote(\"n_orphan_series\", int((~series[\"StudyInstanceUID\"].isin(set(train[\"StudyInstanceUID\"]))).sum()))\n\n# --------------------------------------------------------------------------\n# 2. Figures 1-3: labels\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figures 1-3: label structure ===\")\n\nprev = lab[LABELS].mean().sort_values(ascending=False)\nfig, ax = plt.subplots(figsize=(10, 5))\nbars = ax.bar(range(len(prev)), prev.values * 100, color=\"#4C72B0\")\nax.set_xticks(range(len(prev)))\nax.set_xticklabels(prev.index, rotation=45, ha=\"right\")\nax.set_ylabel(\"Prevalence (%)\")\nax.set_title(f\"Figure 1 — Label prevalence across {n_lab:,} annotated studies\")\nfor b, v, n in zip(bars, prev.values, lab[LABELS][prev.index].sum()):\n    ax.text(b.get_x() + b.get_width() / 2, v * 100 + 0.6, f\"{v*100:.1f}%\\n(n={int(n)})\",\n            ha=\"center\", fontsize=7)\nax.grid(axis=\"y\", alpha=0.3)\nsave(fig, \"fig01_label_prevalence\")\nnote(\"fig1_obs\",\n     f\"prevalence ranges from {prev.min()*100:.1f}% ({prev.idxmin()}, n={int(lab[prev.idxmin()].sum())}) \"\n     f\"to {prev.max()*100:.1f}% ({prev.idxmax()}); \"\n     f\"{int((prev < 0.05).sum())} of 12 labels fall below 5%\")\n\nnpos = lab[LABELS].sum(axis=1)\nfig, ax = plt.subplots(figsize=(9, 4.5))\nax.hist(npos, bins=np.arange(-0.5, 13.5, 1), color=\"#55A868\", edgecolor=\"white\")\nax.set_xlabel(\"Number of positive labels in a study\")\nax.set_ylabel(\"Studies\")\nax.set_title(\"Figure 2 — Multi-label density per study\")\nax.grid(axis=\"y\", alpha=0.3)\nsave(fig, \"fig02_labels_per_study\")\nnote(\"fig2_obs\",\n     f\"mean {npos.mean():.2f} positive labels per study (median {int(npos.median())}, mode {int(npos.mode()[0])}); \"\n     f\"{100*(npos == 0).mean():.1f}% of studies are all-negative and \"\n     f\"{100*(npos >= 2).mean():.1f}% carry two or more findings\")\n\n# Phi coefficient matrix\nM = lab[LABELS].values.astype(float)\nphi = np.corrcoef(M.T)                      # for binary vars, Pearson r == phi\njac = np.zeros_like(phi)\nfor i in range(len(LABELS)):\n    for j in range(len(LABELS)):\n        inter = np.logical_and(M[:, i], M[:, j]).sum()\n        union = np.logical_or(M[:, i], M[:, j]).sum()\n        jac[i, j] = inter / union if union else 0.0\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 6.5))\nfor ax, mat, ttl, cmap in [(axes[0], phi, \"phi coefficient\", \"RdBu_r\"),\n                           (axes[1], jac, \"Jaccard index\", \"viridis\")]:\n    im = ax.imshow(mat, cmap=cmap, vmin=-1 if cmap == \"RdBu_r\" else 0, vmax=1)\n    ax.set_xticks(range(len(LABELS)))\n    ax.set_xticklabels(LABELS, rotation=45, ha=\"right\", fontsize=8)\n    ax.set_yticks(range(len(LABELS)))\n    ax.set_yticklabels(LABELS, fontsize=8)\n    ax.set_title(ttl)\n    for i in range(len(LABELS)):\n        for j in range(len(LABELS)):\n            ax.text(j, i, f\"{mat[i, j]:.2f}\", ha=\"center\", va=\"center\", fontsize=6)\n    fig.colorbar(im, ax=ax, fraction=0.046)\nfig.suptitle(\"Figure 3 — Label co-occurrence structure\", y=1.01)\nsave(fig, \"fig03_cooccurrence\")\noff = phi.copy()\nnp.fill_diagonal(off, 0)\ni, j = np.unravel_index(np.argmax(off), off.shape)\nnote(\"fig3_obs\",\n     f\"strongest association is {LABELS[i]}-{LABELS[j]} (phi={off[i, j]:.2f}); \"\n     f\"{int((np.abs(np.triu(off, 1)) > 0.3).sum())} label pairs exceed |phi|>0.30, \"\n     f\"confirming the labels are not independent\")\npd.DataFrame(phi, index=LABELS, columns=LABELS).to_csv(\"label_phi_matrix.csv\")\n\n# --------------------------------------------------------------------------\n# 3. Figures 4-5: report text\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figures 4-5: report text ===\")\n\ntrain[\"Report\"] = train[\"Report\"].fillna(\"\").astype(str)\ntrain[\"report_clean\"] = train[\"Report\"].map(\n    lambda s: re.sub(r\"\\s+\", \" \", unicodedata.normalize(\"NFKC\", s)).strip()\n)\ntrain[\"report_len_chars\"] = train[\"report_clean\"].str.len()\ntrain[\"report_len_tokens\"] = train[\"report_clean\"].str.split().str.len()\n\nnote(\"min_report_chars\", MIN_REPORT_CHARS)\nnote(\"n_empty_reports\", int((train[\"report_len_chars\"] < MIN_REPORT_CHARS).sum()))\nnote(\"pct_reports_over_512\", round(100 * (train[\"report_len_tokens\"] > 512).mean(), 2))\nnote(\"trunc_note\",\n     \"negligible\" if (train[\"report_len_tokens\"] > 512).mean() < 0.02\n     else \"non-trivial; a sliding-window or longer-context encoder should be considered\")\n\ntry:\n    from langdetect import DetectorFactory, detect\n    DetectorFactory.seed = SEED\n\n    def lang_of(s: str) -> str:\n        if len(s) < MIN_REPORT_CHARS:\n            return \"unknown\"\n        try:\n            return detect(s)\n        except Exception:\n            return \"unknown\"\n\n    train[\"report_lang\"] = train[\"report_clean\"].map(lang_of)\nexcept ImportError:\n    print(\"  ! langdetect not installed - run: pip install langdetect\")\n    train[\"report_lang\"] = \"unknown\"\n\nlang_counts = train[\"report_lang\"].value_counts()\nfig, ax = plt.subplots(figsize=(10, 4.5))\ntop = lang_counts.head(15)\nax.bar(top.index, top.values, color=\"#C44E52\")\nax.set_ylabel(\"Studies\")\nax.set_xlabel(\"Detected language\")\nax.set_title(\"Figure 4 — Report language distribution (top 15)\")\nfor k, (x, v) in enumerate(top.items()):\n    ax.text(k, v, f\"{100*v/len(train):.1f}%\", ha=\"center\", va=\"bottom\", fontsize=7)\nax.grid(axis=\"y\", alpha=0.3)\nsave(fig, \"fig04_report_language\")\nnote(\"fig4_obs\",\n     f\"{len(lang_counts)} languages detected; the dominant language is '{lang_counts.index[0]}' at \"\n     f\"{100*lang_counts.iloc[0]/len(train):.1f}%, and \"\n     f\"{int((lang_counts / len(train) < 0.01).sum())} languages each account for under 1% of studies\")\nnote(\"pct_dominant_lang\", round(100 * lang_counts.iloc[0] / len(train), 1))\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4.5))\naxes[0].hist(train[\"report_len_tokens\"].clip(upper=1200), bins=60, color=\"#8172B2\", edgecolor=\"white\")\naxes[0].axvline(512, color=\"red\", ls=\"--\", label=\"512-token budget\")\naxes[0].set_xlabel(\"Report length (whitespace tokens)\")\naxes[0].set_ylabel(\"Studies\")\naxes[0].legend()\naxes[0].set_title(\"Overall\")\ntop_langs = lang_counts.head(6).index\ndata = [train.loc[train[\"report_lang\"] == L, \"report_len_tokens\"].clip(upper=1200) for L in top_langs]\naxes[1].boxplot(data, labels=list(top_langs), showfliers=False)\naxes[1].set_ylabel(\"Report length (tokens)\")\naxes[1].set_title(\"By language (top 6)\")\nfor a in axes:\n    a.grid(axis=\"y\", alpha=0.3)\nfig.suptitle(\"Figure 5 — Report length distribution\", y=1.02)\nsave(fig, \"fig05_report_length\")\nnote(\"fig5_obs\",\n     f\"median {int(train['report_len_tokens'].median())} tokens \"\n     f\"(IQR {int(train['report_len_tokens'].quantile(.25))}-{int(train['report_len_tokens'].quantile(.75))}, \"\n     f\"max {int(train['report_len_tokens'].max())}); length varies noticeably by language, \"\n     f\"indicating site-specific reporting styles\")\n\n# --------------------------------------------------------------------------\n# 4. Figures 6-8: series structure\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figures 6-8: series structure ===\")\n\nper_study = series.groupby(\"StudyInstanceUID\").agg(\n    n_series=(\"SeriesInstanceUID\", \"count\"),\n    n_fluid_sensitive=(\"Fluid_Sensitive\", \"sum\"),\n    n_fat_suppressed=(\"Fat_Suppression\", \"sum\"),\n)\nplane_ct = (series.pivot_table(index=\"StudyInstanceUID\", columns=\"Anatomical_Plane\",\n                              values=\"SeriesInstanceUID\", aggfunc=\"count\")\n            .fillna(0).astype(int))\nfor p in [\"Sagittal\", \"Coronal\", \"Axial\"]:\n    if p not in plane_ct.columns:\n        plane_ct[p] = 0\nplane_ct = plane_ct[[\"Sagittal\", \"Coronal\", \"Axial\"]]\nper_study = per_study.join(plane_ct)\nper_study.columns = [\"n_series\", \"n_fluid_sensitive\", \"n_fat_suppressed\",\n                     \"n_sagittal\", \"n_coronal\", \"n_axial\"]\nfor p in [\"sagittal\", \"coronal\", \"axial\"]:\n    per_study[f\"has_{p}\"] = (per_study[f\"n_{p}\"] > 0).astype(int)\n\nnote(\"mean_series_per_study\", round(per_study[\"n_series\"].mean(), 2))\ncombo = per_study.apply(\n    lambda r: \"+\".join([n for n, f in [(\"Sag\", r.has_sagittal), (\"Cor\", r.has_coronal),\n                                       (\"Ax\", r.has_axial)] if f]) or \"none\", axis=1)\ncombo_counts = combo.value_counts()\nall_three = int((per_study[[\"has_sagittal\", \"has_coronal\", \"has_axial\"]].sum(axis=1) == 3).sum())\nnote(\"pct_missing_plane\", round(100 * (1 - all_three / len(per_study)), 1))\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4.5))\naxes[0].hist(per_study[\"n_series\"].clip(upper=20), bins=np.arange(0.5, 21.5, 1),\n             color=\"#4C72B0\", edgecolor=\"white\")\naxes[0].set_xlabel(\"Series per study\")\naxes[0].set_ylabel(\"Studies\")\naxes[0].set_title(\"Series count\")\naxes[1].bar(combo_counts.index, combo_counts.values, color=\"#55A868\")\naxes[1].set_ylabel(\"Studies\")\naxes[1].set_title(\"Plane availability\")\naxes[1].tick_params(axis=\"x\", rotation=30)\nfor a in axes:\n    a.grid(axis=\"y\", alpha=0.3)\nfig.suptitle(\"Figure 6 — Series per study and plane availability\", y=1.02)\nsave(fig, \"fig06_series_and_planes\")\nnote(\"fig6_obs\",\n     f\"median {int(per_study['n_series'].median())} series per study; \"\n     f\"{100*all_three/len(per_study):.1f}% of studies contain all three planes, so \"\n     f\"{FILL['pct_missing_plane']}% are missing at least one and require the plane-presence mask\")\n\nct = pd.crosstab(series[\"Fluid_Sensitive\"], series[\"Fat_Suppression\"])\nphi_ff = np.corrcoef(series[\"Fluid_Sensitive\"], series[\"Fat_Suppression\"])[0, 1]\nfig, ax = plt.subplots(figsize=(6, 5))\nim = ax.imshow(ct.values, cmap=\"Blues\")\nax.set_xticks(range(ct.shape[1])); ax.set_xticklabels(ct.columns)\nax.set_yticks(range(ct.shape[0])); ax.set_yticklabels(ct.index)\nax.set_xlabel(\"Fat_Suppression\"); ax.set_ylabel(\"Fluid_Sensitive\")\nfor i in range(ct.shape[0]):\n    for j in range(ct.shape[1]):\n        ax.text(j, i, f\"{ct.values[i, j]:,}\\n({100*ct.values[i, j]/ct.values.sum():.1f}%)\",\n                ha=\"center\", va=\"center\",\n                color=\"white\" if ct.values[i, j] > ct.values.max() / 2 else \"black\")\nax.set_title(f\"Figure 8 — Fluid_Sensitive x Fat_Suppression (phi = {phi_ff:.2f})\")\nfig.colorbar(im, ax=ax, fraction=0.046)\nsave(fig, \"fig08_fluid_vs_fatsup\")\noff_diag = ct.values.sum() - np.trace(ct.values)\nnote(\"fig8_obs\",\n     f\"phi = {phi_ff:.2f}; {100*off_diag/ct.values.sum():.1f}% of series fall off the diagonal, \"\n     f\"confirming the organisers' note that the two flags are correlated but not equivalent - \"\n     f\"they are retained as two separate features\")\n\n# --------------------------------------------------------------------------\n# 5. Figures 7 & 9: DICOM sampling\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figures 7 & 9: DICOM header/pixel sampling ===\")\n\ntry:\n    import pydicom\n\n    sample = series.sample(min(N_HEADER_SAMPLE, len(series)), random_state=SEED)\n    rows = []\n    for _, r in sample.iterrows():\n        d = ROOT / \"train_series\" / r[\"StudyInstanceUID\"] / r[\"SeriesInstanceUID\"]\n        if not d.exists():\n            continue\n        files = sorted(d.glob(\"*.dcm\"))\n        if not files:\n            continue\n        try:\n            hdr = pydicom.dcmread(files[0], stop_before_pixels=True)\n            iop = getattr(hdr, \"ImageOrientationPatient\", None)\n            derived = None\n            if iop is not None and len(iop) == 6:\n                nrm = np.cross(np.array(iop[:3], float), np.array(iop[3:], float))\n                derived = [\"Sagittal\", \"Coronal\", \"Axial\"][int(np.argmax(np.abs(nrm)))]\n            ps = getattr(hdr, \"PixelSpacing\", [np.nan, np.nan])\n            rows.append({\n                \"plane_given\": r[\"Anatomical_Plane\"],\n                \"plane_derived\": derived,\n                \"n_slices\": len(files),\n                \"rows\": int(getattr(hdr, \"Rows\", 0)),\n                \"cols\": int(getattr(hdr, \"Columns\", 0)),\n                \"pixel_spacing\": float(ps[0]) if ps is not None else np.nan,\n                \"transfer_syntax\": str(getattr(hdr.file_meta, \"TransferSyntaxUID\", \"unknown\")),\n                \"path\": str(d),\n            })\n        except Exception as e:\n            print(f\"  ! header read failed: {e}\")\n\n    hdrs = pd.DataFrame(rows)\n    print(f\"  read {len(hdrs)} series headers\")\n\n    mism = hdrs.dropna(subset=[\"plane_derived\"])\n    n_mismatch = int((mism[\"plane_given\"] != mism[\"plane_derived\"]).sum())\n    note(\"n_plane_mismatch\",\n         f\"{n_mismatch} of {len(mism)} sampled series \"\n         f\"({100*n_mismatch/max(len(mism),1):.1f}%)\")\n\n    fig, ax = plt.subplots(figsize=(9, 4.5))\n    ax.hist(hdrs[\"n_slices\"], bins=np.logspace(np.log10(max(hdrs[\"n_slices\"].min(), 1)),\n                                               np.log10(hdrs[\"n_slices\"].max()), 40),\n            color=\"#CCB974\", edgecolor=\"white\")\n    ax.set_xscale(\"log\")\n    ax.axvline(TARGET_SLICES, color=\"red\", ls=\"--\", label=f\"{TARGET_SLICES}-slice sampling target\")\n    ax.set_xlabel(\"Slices per series (log scale)\")\n    ax.set_ylabel(\"Series\")\n    ax.legend()\n    ax.set_title(f\"Figure 7 — Slice count distribution ({len(hdrs)} sampled series)\")\n    ax.grid(axis=\"y\", alpha=0.3)\n    save(fig, \"fig07_slices_per_series\")\n    over = 100 * (hdrs[\"n_slices\"] > 2 * TARGET_SLICES).mean()\n    note(\"fig7_obs\",\n         f\"median {int(hdrs['n_slices'].median())} slices \"\n         f\"(IQR {int(hdrs['n_slices'].quantile(.25))}-{int(hdrs['n_slices'].quantile(.75))}, \"\n         f\"max {int(hdrs['n_slices'].max())}); {over:.1f}% of series exceed 2x the \"\n         f\"{TARGET_SLICES}-slice budget and are therefore materially subsampled\")\n\n    # Figure 9 — pixel-level heterogeneity\n    fig, axes = plt.subplots(1, 3, figsize=(18, 4.5))\n    mat = (hdrs[\"rows\"].astype(str) + \"x\" + hdrs[\"cols\"].astype(str)).value_counts().head(10)\n    axes[0].bar(mat.index, mat.values, color=\"#4C72B0\")\n    axes[0].tick_params(axis=\"x\", rotation=45)\n    axes[0].set_title(\"Native matrix size (top 10)\")\n    axes[0].set_ylabel(\"Series\")\n\n    axes[1].hist(hdrs[\"pixel_spacing\"].dropna(), bins=40, color=\"#55A868\", edgecolor=\"white\")\n    axes[1].set_xlabel(\"In-plane pixel spacing (mm)\")\n    axes[1].set_title(\"Pixel spacing\")\n\n    decoded = 0\n    for _, r in hdrs.sample(min(N_PIXEL_SAMPLE, len(hdrs)), random_state=SEED).iterrows():\n        try:\n            f = sorted(Path(r[\"path\"]).glob(\"*.dcm\"))[len(list(Path(r[\"path\"]).glob(\"*.dcm\"))) // 2]\n            arr = pydicom.dcmread(f).pixel_array.astype(float)\n            axes[2].hist(arr.ravel(), bins=60, histtype=\"step\", alpha=0.5, density=True)\n            decoded += 1\n        except Exception:\n            continue\n    axes[2].set_xlabel(\"Raw pixel intensity\")\n    axes[2].set_title(f\"Raw intensity histograms ({decoded} volumes, un-normalised)\")\n    for a in axes:\n        a.grid(axis=\"y\", alpha=0.3)\n    fig.suptitle(\"Figure 9 — Pixel-level heterogeneity\", y=1.03)\n    save(fig, \"fig09_pixel_heterogeneity\")\n\n    note(\"n_pixel_sample\", decoded)\n    note(\"n_sampled_series\", len(hdrs))\n    note(\"fig9_obs\",\n         f\"{hdrs['rows'].nunique()} distinct row counts and \"\n         f\"{int(mat.shape[0])} common matrix sizes; pixel spacing spans \"\n         f\"{hdrs['pixel_spacing'].min():.2f}-{hdrs['pixel_spacing'].max():.2f} mm; \"\n         f\"raw intensity histograms do not share a common range, which is the direct \"\n         f\"justification for per-volume percentile normalisation\")\n    note(\"transfer_syntaxes\", hdrs[\"transfer_syntax\"].value_counts().to_dict())\n\n    mean_sl = hdrs[\"n_slices\"].median()\n    est_gb = len(series) * TARGET_SLICES * 256 * 256 / 1e9\n    note(\"cache_size_gb\", round(est_gb, 1))\n    note(\"cache_hours\", round(len(series) * mean_sl / 60000, 1))\n\nexcept ImportError:\n    print(\"  ! pydicom not installed - skipping Figures 7 and 9\")\n    for t in [\"n_plane_mismatch\", \"fig7_obs\", \"fig9_obs\", \"n_pixel_sample\",\n              \"n_sampled_series\", \"cache_size_gb\", \"cache_hours\"]:\n        note(t, \"NOT COMPUTED - install pydicom\")\n\n# --------------------------------------------------------------------------\n# 6. Figures 10-11: confound audits\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figures 10-11: confound audits ===\")\n\nnote(\"min_lang_n\", MIN_LANG_N)\nlab2 = train[train[\"is_labeled\"]].copy()\nlab2[LABELS] = lab2[LABELS].astype(\"int8\")\nkeep = [L for L, c in lab2[\"report_lang\"].value_counts().items() if c >= MIN_LANG_N]\nsub = lab2[lab2[\"report_lang\"].isin(keep)]\n\nif len(keep) >= 2:\n    prev_lang = sub.groupby(\"report_lang\")[LABELS].mean() * 100\n    fig, ax = plt.subplots(figsize=(15, 5.5))\n    w = 0.8 / len(prev_lang)\n    x = np.arange(len(LABELS))\n    for k, (L, row) in enumerate(prev_lang.iterrows()):\n        ax.bar(x + k * w, row[LABELS].values, w,\n               label=f\"{L} (n={int((sub['report_lang'] == L).sum())})\")\n    ax.set_xticks(x + 0.4 - w / 2)\n    ax.set_xticklabels(LABELS, rotation=45, ha=\"right\")\n    ax.set_ylabel(\"Prevalence (%)\")\n    ax.legend(fontsize=8, ncol=2)\n    ax.set_title(\"Figure 10 — Label prevalence by report language (site-confound audit)\")\n    ax.grid(axis=\"y\", alpha=0.3)\n    save(fig, \"fig10_prevalence_by_language\")\n    spread = (prev_lang.max() - prev_lang.min())\n    note(\"fig10_obs\",\n         f\"across {len(keep)} language groups, per-label prevalence spreads by up to \"\n         f\"{spread.max():.1f} percentage points (largest for {spread.idxmax()}); \"\n         f\"language is therefore predictive of the label independently of anatomy, \"\n         f\"a shortcut both branches could exploit\")\nelse:\n    note(\"fig10_obs\", \"insufficient language diversity above the minimum group size to run this audit\")\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4.5))\nfor flag, name, c in [(True, \"labeled\", \"#4C72B0\"), (False, \"unlabeled\", \"#C44E52\")]:\n    d = train.loc[train[\"is_labeled\"] == flag, \"report_len_tokens\"].clip(upper=1000)\n    if len(d):\n        axes[0].hist(d, bins=50, alpha=0.55, label=f\"{name} (n={len(d):,})\",\n                     density=True, color=c)\naxes[0].set_xlabel(\"Report length (tokens)\")\naxes[0].set_ylabel(\"Density\")\naxes[0].legend()\naxes[0].set_title(\"Report length\")\n\nlp = (train.groupby(\"is_labeled\")[\"report_lang\"].value_counts(normalize=True)\n      .unstack(0).fillna(0).head(10) * 100)\nlp.plot(kind=\"bar\", ax=axes[1], color=[\"#C44E52\", \"#4C72B0\"])\naxes[1].set_ylabel(\"% of pool\")\naxes[1].set_title(\"Language mix\")\naxes[1].legend([\"unlabeled\", \"labeled\"])\nfor a in axes:\n    a.grid(axis=\"y\", alpha=0.3)\nfig.suptitle(\"Figure 11 — Labeled vs unlabeled pool comparison\", y=1.02)\nsave(fig, \"fig11_labeled_vs_unlabeled\")\n\nml = train.loc[train[\"is_labeled\"], \"report_len_tokens\"].median()\nmu = train.loc[~train[\"is_labeled\"], \"report_len_tokens\"].median()\nlang_l = set(train.loc[train[\"is_labeled\"], \"report_lang\"].value_counts().head(5).index)\nlang_u = set(train.loc[~train[\"is_labeled\"], \"report_lang\"].value_counts().head(5).index)\nnote(\"fig11_obs\",\n     f\"labeled pool n={int(train['is_labeled'].sum()):,}, unlabeled n={int((~train['is_labeled']).sum()):,}; \"\n     f\"median report length {ml:.0f} vs {mu:.0f} tokens; \"\n     f\"{len(lang_l & lang_u)}/5 of the top-5 languages are shared between the pools\")\nnote(\"mar_check_result\",\n     f\"median report length {ml:.0f} tokens (labeled) vs {mu:.0f} (unlabeled); \"\n     + (\"the two pools look comparable, weakly supporting a missing-at-random assumption\"\n        if abs(ml - mu) / max(ml, 1) < 0.15 else\n        \"the pools differ materially, so the annotation subset is NOT plausibly random \"\n        \"and pseudo-labelling may not transfer\"))\n\n# --------------------------------------------------------------------------\n# 7. Figure 12: rule extractor validation\n# --------------------------------------------------------------------------\nprint(\"\\n=== Figure 12: rule extractor validation ===\")\n\n# Multilingual keyword lexicons. EXTEND THESE after reading real reports.\nLEX = {\n    \"ACL\": [r\"anterior cruciate\", r\"\\bacl\\b\", r\"\\blca\\b\", r\"ligamento cruzado anterior\",\n            r\"ligament crois[ée] ant[ée]rieur\", r\"vorderes kreuzband\"],\n    \"MCL\": [r\"medial collateral\", r\"\\bmcl\\b\", r\"ligamento colateral medial\",\n            r\"ligament collat[ée]ral m[ée]dial\", r\"innenband\"],\n    \"Medial Meniscus\": [r\"medial meniscus\", r\"menisco medial\", r\"m[ée]nisque m[ée]dial\",\n                        r\"innenmeniskus\"],\n    \"Lateral Meniscus\": [r\"lateral meniscus\", r\"menisco lateral\", r\"m[ée]nisque lat[ée]ral\",\n                         r\"au[sß]enmeniskus\"],\n    \"Medial OA\": [r\"medial (compartment )?(osteoarthr|chondral|cartilage loss|degenerat)\",\n                  r\"artrosis medial\", r\"gonarthrose m[ée]diale\"],\n    \"Lateral OA\": [r\"lateral (compartment )?(osteoarthr|chondral|cartilage loss|degenerat)\",\n                   r\"artrosis lateral\", r\"gonarthrose lat[ée]rale\"],\n    \"PF OA\": [r\"patellofemoral\", r\"patell[oa]-?femoral\", r\"retropatellar\", r\"chondromalacia patell\"],\n    \"Effusion\": [r\"effusion\", r\"joint fluid\", r\"derrame\", r\"[ée]panchement\", r\"erguss\"],\n    \"Synovitis\": [r\"synovitis\", r\"synovial (thickening|proliferat)\", r\"sinovitis\", r\"synovite\"],\n    \"Baker's\": [r\"baker'?s? cyst\", r\"popliteal cyst\", r\"quiste de baker\", r\"kyste de baker\",\n                r\"bakerzyste\"],\n    \"Contusion\": [r\"contusion\", r\"bone (bruise|marrow o?edema)\", r\"contusi[óo]n\",\n                  r\"knochenmark[sö]dem\"],\n    \"Fracture\": [r\"fracture\", r\"fractura\", r\"\\bfx\\b\", r\"fraktur\"],\n}\n\nNEG_CUES = [r\"\\bno\\b\", r\"\\bnot\\b\", r\"\\bwithout\\b\", r\"absence of\", r\"negative for\", r\"free of\",\n            r\"\\bnegative\\b\", r\"\\bsin\\b\", r\"ausencia de\", r\"\\bpas de\\b\", r\"\\bsans\\b\",\n            r\"\\bkein\\b\", r\"\\bohne\\b\", r\"\\bnessun\", r\"\\bsem\\b\", r\"unremarkable\", r\"intact\"]\nNEG_SCOPE = 8   # tokens after a negation cue\n\ndef negation_mask(text: str) -> list[bool]:\n    \"\"\"True for tokens that fall inside a negation scope.\"\"\"\n    toks = text.split()\n    mask = [False] * len(toks)\n    for i, t in enumerate(toks):\n        if any(re.search(c, t) for c in NEG_CUES):\n            for j in range(i, min(i + NEG_SCOPE + 1, len(toks))):\n                mask[j] = True\n    return mask\n\ndef extract(text: str, use_negation: bool = True) -> dict[str, int]:\n    toks = text.split()\n    mask = negation_mask(text) if use_negation else [False] * len(toks)\n    visible = \" \".join(t for t, m in zip(toks, mask) if not m)\n    return {L: int(any(re.search(p, visible) for p in pats)) for L, pats in LEX.items()}\n\nfrom sklearn.metrics import precision_recall_fscore_support  # noqa: E402\n\nresults = {}\nfor use_neg in (True, False):\n    preds = pd.DataFrame([extract(t, use_neg) for t in lab2[\"report_clean\"]], index=lab2.index)\n    p, r, f, _ = precision_recall_fscore_support(\n        lab2[LABELS].values, preds[LABELS].values, average=None, zero_division=0)\n    results[use_neg] = pd.DataFrame({\"precision\": p, \"recall\": r, \"f1\": f}, index=LABELS)\n\nfig, ax = plt.subplots(figsize=(14, 5))\nx = np.arange(len(LABELS))\nax.bar(x - 0.32, results[True][\"precision\"], 0.16, label=\"precision (neg ON)\", color=\"#4C72B0\")\nax.bar(x - 0.16, results[True][\"recall\"], 0.16, label=\"recall (neg ON)\", color=\"#8FB3DD\")\nax.bar(x + 0.16, results[False][\"precision\"], 0.16, label=\"precision (neg OFF)\", color=\"#C44E52\")\nax.bar(x + 0.32, results[False][\"recall\"], 0.16, label=\"recall (neg OFF)\", color=\"#E8A0A2\")\nax.set_xticks(x); ax.set_xticklabels(LABELS, rotation=45, ha=\"right\")\nax.set_ylabel(\"Score\"); ax.set_ylim(0, 1)\nax.legend(fontsize=8, ncol=2)\nax.set_title(\"Figure 12 — Rule-based extractor vs gold labels, with and without negation handling\")\nax.grid(axis=\"y\", alpha=0.3)\nsave(fig, \"fig12_rule_extractor\")\n\nd_prec = (results[True][\"precision\"] - results[False][\"precision\"]).mean()\nnote(\"fig12_obs\",\n     f\"with negation handling, macro precision {results[True]['precision'].mean():.3f} / \"\n     f\"recall {results[True]['recall'].mean():.3f} / F1 {results[True]['f1'].mean():.3f}; \"\n     f\"disabling negation changes mean precision by {d_prec:+.3f}, \"\n     f\"quantifying the value of the negation module. Best label: \"\n     f\"{results[True]['f1'].idxmax()} (F1={results[True]['f1'].max():.2f}); worst: \"\n     f\"{results[True]['f1'].idxmin()} (F1={results[True]['f1'].min():.2f})\")\nresults[True].to_csv(\"rule_extractor_scores.csv\")\n\n# --------------------------------------------------------------------------\n# 8. Engineered feature table + summary\n# --------------------------------------------------------------------------\nprint(\"\\n=== Writing outputs ===\")\n\nfeats = train[[\"StudyInstanceUID\", \"is_labeled\", \"report_len_chars\",\n               \"report_len_tokens\", \"report_lang\"]].copy()\nrule_df = pd.DataFrame([extract(t) for t in train[\"report_clean\"]], index=train.index)\nrule_df.columns = [f\"rule_hit_{SNAKE[c]}\" for c in rule_df.columns]\nfeats = pd.concat([feats, rule_df], axis=1)\nfeats[\"n_affirmed_terms\"] = rule_df.sum(axis=1)\nfeats = feats.merge(per_study, left_on=\"StudyInstanceUID\", right_index=True, how=\"left\")\nfor L in LABELS:\n    feats[SNAKE[L]] = train[L].values\nfeats[\"n_labels_positive\"] = train[LABELS].sum(axis=1)\n\ntry:\n    feats.to_parquet(\"study_features.parquet\", index=False)\n    print(f\"  -> study_features.parquet  {feats.shape}\")\nexcept ImportError:\n    feats.to_csv(\"study_features.csv\", index=False)\n    print(f\"  -> study_features.csv  {feats.shape}  (pyarrow unavailable; wrote CSV)\")\n\nnote(\"n_final\", int(feats[\"StudyInstanceUID\"].isin(studies_with_series).sum()))\n\nwith open(\"eda_summary.json\", \"w\") as fh:\n    json.dump(FILL, fh, indent=2, default=str)\nprint(\"  -> eda_summary.json\")\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"DONE. Paste values from eda_summary.json into the [FILL: ...] slots\")\nprint(\"in Phase2_Data_Collection_Preprocessing_EDA.md and Phase3_Methodology.md\")\nprint(\"=\" * 70)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:21:29.178811Z","iopub.execute_input":"2026-08-25T17:21:29.179333Z","iopub.status.idle":"2026-08-25T17:21:29.198082Z","shell.execute_reply.started":"2026-08-25T17:21:29.179306Z","shell.execute_reply":"2026-08-25T17:21:29.197434Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!python phase2_eda.py","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:21:41.91535Z","iopub.execute_input":"2026-08-25T17:21:41.915796Z","iopub.status.idle":"2026-08-25T17:21:43.188424Z","shell.execute_reply.started":"2026-08-25T17:21:41.915766Z","shell.execute_reply":"2026-08-25T17:21:43.187738Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nprint(os.listdir(\"/kaggle/input/\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:22:45.990243Z","iopub.execute_input":"2026-08-25T17:22:45.990816Z","iopub.status.idle":"2026-08-25T17:22:45.995395Z","shell.execute_reply.started":"2026-08-25T17:22:45.990781Z","shell.execute_reply":"2026-08-25T17:22:45.994786Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ROOT = Path(\"/kaggle/input/competitions\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:23:40.916887Z","iopub.execute_input":"2026-08-25T17:23:40.917614Z","iopub.status.idle":"2026-08-25T17:23:40.933474Z","shell.execute_reply.started":"2026-08-25T17:23:40.917587Z","shell.execute_reply":"2026-08-25T17:23:40.932497Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T21:14:20.904246Z","iopub.execute_input":"2026-08-25T21:14:20.904722Z","iopub.status.idle":"2026-08-25T21:14:20.908598Z","shell.execute_reply.started":"2026-08-25T21:14:20.904691Z","shell.execute_reply":"2026-08-25T21:14:20.907688Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom sklearn.feature_extraction.text import TfidfVectorizer\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import roc_auc_score\n\n# Function to auto-detect where Kaggle mounted the competition dataset\ndef get_dataset_root():\n    candidate_paths = [\n        Path(\"/kaggle/input/rsna-knee-abnormality-detection\"),\n        Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\"),\n    ]\n    for path in candidate_paths:\n        if (path / \"train.csv\").is_file():\n            return path\n    \n    # Fallback search if mounted under a custom directory\n    base = Path(\"/kaggle/input\")\n    if base.is_dir():\n        for file_path in base.rglob(\"train.csv\"):\n            return file_path.parent\n            \n    raise FileNotFoundError(\"Competition dataset is not attached! Click '+ Add Input' on the right sidebar.\")\n\n# 1. Set ROOT automatically\nROOT = get_dataset_root()\nprint(f\"Dataset successfully found at: {ROOT}\")\n\n# 2. Print files inside dataset root\nprint(\"\\nFiles inside ROOT:\", os.listdir(ROOT))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T21:16:45.46316Z","iopub.execute_input":"2026-08-25T21:16:45.463729Z","iopub.status.idle":"2026-08-25T21:16:45.473347Z","shell.execute_reply.started":"2026-08-25T21:16:45.4637Z","shell.execute_reply":"2026-08-25T21:16:45.472469Z"}},"outputs":[],"execution_count":null}]}