{"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":"import os; \nprint(os.listdir(\"/kaggle/input/osic-pulmonary-fibrosis-progression\"))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:35:48.100108Z","iopub.execute_input":"2026-09-07T19:35:48.100449Z","iopub.status.idle":"2026-09-07T19:35:48.110678Z","shell.execute_reply.started":"2026-09-07T19:35:48.100418Z","shell.execute_reply":"2026-09-07T19:35:48.107775Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os; \nprint(os.listdir(\"/kaggle/input/osic-pulmonary-fibrosis-progression\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:39:54.349934Z","iopub.execute_input":"2026-09-07T19:39:54.350281Z","iopub.status.idle":"2026-09-07T19:39:54.364653Z","shell.execute_reply.started":"2026-09-07T19:39:54.350248Z","shell.execute_reply":"2026-09-07T19:39:54.362948Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nfor root, dirs, files in os.walk(\"/kaggle/input\"):\n    if \"train.csv\" in files:\n        print(\"DATA_DIR =\", root)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:40:05.650136Z","iopub.execute_input":"2026-09-07T19:40:05.650487Z","iopub.status.idle":"2026-09-07T19:40:55.035152Z","shell.execute_reply.started":"2026-09-07T19:40:05.650458Z","shell.execute_reply":"2026-09-07T19:40:55.034157Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os; os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:42:28.035215Z","iopub.execute_input":"2026-09-07T19:42:28.035587Z","iopub.status.idle":"2026-09-07T19:42:28.041537Z","shell.execute_reply.started":"2026-09-07T19:42:28.035559Z","shell.execute_reply":"2026-09-07T19:42:28.040643Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install shap lightgbm xgboost pydicom -q","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:43:09.674609Z","iopub.execute_input":"2026-09-07T19:43:09.674914Z","iopub.status.idle":"2026-09-07T19:43:15.241377Z","shell.execute_reply.started":"2026-09-07T19:43:09.674888Z","shell.execute_reply":"2026-09-07T19:43:15.239907Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os; os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\"\"\"\nSTEP 1 - CLINICAL PIPELINE (tabular data only)\n------------------------------------------------\nWhat it does:\n  1. Loads train.csv from the OSIC dataset\n  2. Builds (baseline visit -> future visit) training pairs per patient\n  3. Engineers time + history features\n  4. Runs 5-fold PATIENT-GROUPED cross-validation for 5 models\n  5. Trains quantile models (10/50/90 %) for prediction intervals\n  6. Computes MAE, RMSE, and the OSIC modified Laplace log-likelihood\n  7. Trains a rapid-progression classifier (>=10% FVC drop within 52 wk)\n  8. Saves results tables, SHAP plots and the final model to OUT_DIR\n\nHow to run on Kaggle:  paste this whole file into one notebook cell and run it.\n\"\"\"\n\nimport os, warnings, json\nimport numpy as np\nimport pandas as pd\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.linear_model import Ridge, LogisticRegression\nfrom sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error, roc_auc_score, confusion_matrix\nimport lightgbm as lgb\nimport xgboost as xgb\nimport joblib\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nCT_FEATURES_CSV = os.environ.get(\"CT_FEATURES_CSV\", \"\")   # leave empty for clinical-only run\nN_FOLDS = 5\nSEED = 42\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\nnp.random.seed(SEED)\n\nCAT_COLS = [\"Sex\", \"SmokingStatus\"]\n\n\n# ----------------------------------------------------------------\n# 1. LOAD + CLEAN\n# ----------------------------------------------------------------\ndef load_clinical():\n    df = pd.read_csv(os.path.join(DATA_DIR, \"train.csv\"))\n    df = df.drop_duplicates(subset=[\"Patient\", \"Weeks\"]).copy()\n    df = df.sort_values([\"Patient\", \"Weeks\"]).reset_index(drop=True)\n    print(f\"Loaded {len(df)} rows, {df.Patient.nunique()} patients\")\n    return df\n\n\n# ----------------------------------------------------------------\n# 2. BUILD TRAINING PAIRS  (baseline -> every later visit)\n# ----------------------------------------------------------------\ndef build_pairs(df):\n    rows = []\n    for pid, g in df.groupby(\"Patient\"):\n        g = g.sort_values(\"Weeks\").reset_index(drop=True)\n        base = g.iloc[0]                        # earliest visit = baseline\n        for i in range(1, len(g)):\n            tgt = g.iloc[i]\n            hist = g.iloc[:i]                   # visits known BEFORE the target\n            # patient-specific slope so far (needs >=2 points)\n            if len(hist) >= 2:\n                slope = np.polyfit(hist.Weeks, hist.FVC, 1)[0]\n            else:\n                slope = 0.0\n            rows.append(dict(\n                Patient=pid,\n                Age=base.Age, Sex=base.Sex, SmokingStatus=base.SmokingStatus,\n                base_fvc=base.FVC, base_pct=base.Percent, base_week=base.Weeks,\n                target_week=tgt.Weeks,\n                week_diff=tgt.Weeks - base.Weeks,\n                # ---- history features (Experiment B) ----\n                last_fvc=hist.FVC.iloc[-1],\n                last_week=hist.Weeks.iloc[-1],\n                weeks_since_last=tgt.Weeks - hist.Weeks.iloc[-1],\n                n_prev_visits=len(hist),\n                prev_fvc_mean=hist.FVC.mean(),\n                prev_fvc_std=hist.FVC.std() if len(hist) > 1 else 0.0,\n                slope_so_far=slope,\n                decline_so_far=hist.FVC.iloc[-1] - base.FVC,\n                # ---- targets ----\n                target_fvc=tgt.FVC,\n                target_delta=tgt.FVC - base.FVC,\n            ))\n    pairs = pd.DataFrame(rows)\n    pairs[\"horizon_bucket\"] = pd.cut(pairs.week_diff, [-1, 12, 26, 52, 999],\n                                     labels=[\"0-12\", \"13-26\", \"27-52\", \"52+\"])\n    print(f\"Built {len(pairs)} training pairs\")\n    return pairs\n\n\nFEATS_A = [\"Age\", \"Sex\", \"SmokingStatus\", \"base_fvc\", \"base_pct\", \"base_week\",\n           \"target_week\", \"week_diff\"]                                   # baseline-only\nFEATS_B = FEATS_A + [\"last_fvc\", \"last_week\", \"weeks_since_last\", \"n_prev_visits\",\n                     \"prev_fvc_mean\", \"prev_fvc_std\", \"slope_so_far\", \"decline_so_far\"]\n\n\ndef encode(X):\n    \"\"\"one-hot the two categorical columns\"\"\"\n    X = X.copy()\n    X[\"Sex\"] = (X[\"Sex\"] == \"Male\").astype(int)\n    for s in [\"Ex-smoker\", \"Never smoked\", \"Currently smokes\"]:\n        X[\"smk_\" + s.replace(\" \", \"_\")] = (X[\"SmokingStatus\"] == s).astype(int)\n    return X.drop(columns=[\"SmokingStatus\"])\n\n\n# ----------------------------------------------------------------\n# 3. METRIC  (OSIC modified Laplace log-likelihood)\n# ----------------------------------------------------------------\ndef laplace_ll(y_true, y_pred, sigma):\n    sigma = np.maximum(sigma, 70)\n    delta = np.minimum(np.abs(y_true - y_pred), 1000)\n    return np.mean(-np.sqrt(2) * delta / sigma - np.log(np.sqrt(2) * sigma))\n\n\n# ----------------------------------------------------------------\n# 4. MODELS\n# ----------------------------------------------------------------\ndef get_models():\n    return {\n        \"Ridge\": make_pipeline(StandardScaler(), Ridge(alpha=1.0)),\n        \"RandomForest\": RandomForestRegressor(n_estimators=400, min_samples_leaf=5,\n                                              random_state=SEED, n_jobs=-1),\n        \"XGBoost\": xgb.XGBRegressor(n_estimators=400, max_depth=4, learning_rate=0.03,\n                                    subsample=0.8, colsample_bytree=0.8, random_state=SEED),\n        \"LightGBM\": lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03,\n                                      min_child_samples=10, random_state=SEED, verbose=-1),\n    }\n\n\ndef quantile_models(alpha_list=(0.1, 0.5, 0.9)):\n    return {a: lgb.LGBMRegressor(objective=\"quantile\", alpha=a, n_estimators=400,\n                                 num_leaves=15, learning_rate=0.03, min_child_samples=10,\n                                 random_state=SEED, verbose=-1) for a in alpha_list}\n\n\n# ----------------------------------------------------------------\n# 5. CROSS-VALIDATION\n# ----------------------------------------------------------------\ndef run_cv(pairs, feats, label):\n    X = encode(pairs[feats])\n    y = pairs[\"target_fvc\"].values\n    groups = pairs[\"Patient\"].values\n    gkf = GroupKFold(n_splits=N_FOLDS)\n\n    oof = {name: np.zeros(len(y)) for name in get_models()}\n    oof_q = {a: np.zeros(len(y)) for a in (0.1, 0.5, 0.9)}\n\n    for fold, (tr, te) in enumerate(gkf.split(X, y, groups)):\n        for name, model in get_models().items():\n            model.fit(X.iloc[tr], y[tr])\n            oof[name][te] = model.predict(X.iloc[te])\n        for a, qm in quantile_models().items():\n            qm.fit(X.iloc[tr], y[tr])\n            oof_q[a][te] = qm.predict(X.iloc[te])\n        print(f\"  [{label}] fold {fold + 1}/{N_FOLDS} done\")\n\n    # patient-mean baseline (predicts baseline FVC for every horizon)\n    oof[\"Baseline(base_fvc)\"] = pairs[\"base_fvc\"].values\n\n    # Ensemble: simple average of tree models\n    oof[\"Ensemble\"] = np.mean([oof[\"RandomForest\"], oof[\"XGBoost\"], oof[\"LightGBM\"]], axis=0)\n\n    sigma = np.maximum((oof_q[0.9] - oof_q[0.1]) / 2.0, 70)   # interval half-width -> sigma\n    results = []\n    for name, p in oof.items():\n        results.append(dict(experiment=label, model=name,\n                            MAE=mean_absolute_error(y, p),\n                            RMSE=np.sqrt(mean_squared_error(y, p)),\n                            Laplace=laplace_ll(y, p, sigma)))\n    res = pd.DataFrame(results).sort_values(\"MAE\")\n    # coverage of the 80% interval (using quantile-median as point)\n    cover = np.mean((y >= oof_q[0.1]) & (y <= oof_q[0.9]))\n    res.attrs[\"coverage80\"] = cover\n\n    # per-horizon breakdown for the ensemble\n    hz = pd.DataFrame({\"bucket\": pairs.horizon_bucket, \"err\": np.abs(y - oof[\"Ensemble\"]),\n                       \"width\": oof_q[0.9] - oof_q[0.1]})\n    hz_tab = hz.groupby(\"bucket\").agg(n=(\"err\", \"size\"), MAE=(\"err\", \"mean\"),\n                                       interval_width=(\"width\", \"mean\")).reset_index()\n    hz_tab[\"experiment\"] = label\n    return res, hz_tab, oof, oof_q, sigma\n\n\n# ----------------------------------------------------------------\n# 6. RAPID PROGRESSION CLASSIFIER  (one row per patient)\n# ----------------------------------------------------------------\ndef rapid_progression(df):\n    rows = []\n    for pid, g in df.groupby(\"Patient\"):\n        g = g.sort_values(\"Weeks\")\n        base = g.iloc[0]\n        within = g[(g.Weeks > base.Weeks) & (g.Weeks - base.Weeks <= 52)]\n        if len(within) == 0:\n            continue\n        label = int((within.FVC <= 0.90 * base.FVC).any())\n        rows.append(dict(Patient=pid, Age=base.Age, Sex=base.Sex, SmokingStatus=base.SmokingStatus,\n                         base_fvc=base.FVC, base_pct=base.Percent, base_week=base.Weeks, label=label))\n    d = pd.DataFrame(rows)\n    X = encode(d[[\"Age\", \"Sex\", \"SmokingStatus\", \"base_fvc\", \"base_pct\", \"base_week\"]])\n    y = d.label.values\n    print(f\"Rapid-progression patients: {y.sum()} / {len(y)}\")\n    clf = RandomForestClassifier(n_estimators=400, class_weight=\"balanced\",\n                                 min_samples_leaf=3, random_state=SEED)\n    from sklearn.model_selection import StratifiedKFold, cross_val_predict\n    proba = cross_val_predict(clf, X, y, cv=StratifiedKFold(5, shuffle=True, random_state=SEED),\n                              method=\"predict_proba\")[:, 1]\n    pred = (proba >= 0.35).astype(int)      # low threshold -> favour sensitivity\n    tn, fp, fn, tp = confusion_matrix(y, pred).ravel()\n    out = dict(AUC=roc_auc_score(y, proba), sensitivity=tp / (tp + fn), specificity=tn / (tn + fp),\n               precision=tp / max(tp + fp, 1), n=len(y), positives=int(y.sum()))\n    clf.fit(X, y)\n    joblib.dump(clf, os.path.join(OUT_DIR, \"rapid_clf.joblib\"))\n    return out\n\n\n# ----------------------------------------------------------------\n# 7. SHAP + PLOTS\n# ----------------------------------------------------------------\ndef shap_plots(model, X, tag):\n    try:\n        import shap\n        expl = shap.TreeExplainer(model)\n        sv = expl.shap_values(X)\n        plt.figure()\n        shap.summary_plot(sv, X, show=False, max_display=15)\n        plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"shap_summary_{tag}.png\"), dpi=150); plt.close()\n        plt.figure()\n        shap.plots._waterfall.waterfall_legacy(expl.expected_value, sv[0], X.iloc[0], show=False, max_display=12)\n        plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"shap_waterfall_patient0_{tag}.png\"), dpi=150); plt.close()\n        print(f\"  SHAP plots saved ({tag})\")\n    except Exception as e:\n        print(\"  SHAP skipped:\", e)\n\n\ndef trajectory_plot(df, pairs, oof, oof_q, n=4):\n    pids = pairs.Patient.unique()[:n]\n    fig, axes = plt.subplots(1, n, figsize=(4 * n, 3.5))\n    for ax, pid in zip(np.atleast_1d(axes), pids):\n        g = df[df.Patient == pid]\n        m = pairs.Patient == pid\n        ax.plot(g.Weeks, g.FVC, \"ko-\", label=\"actual\")\n        ax.plot(pairs.target_week[m], oof[\"Ensemble\"][m], \"r--\", label=\"predicted\")\n        ax.fill_between(pairs.target_week[m], oof_q[0.1][m], oof_q[0.9][m], color=\"r\", alpha=.15, label=\"80% PI\")\n        ax.set_title(pid[:12]); ax.set_xlabel(\"week\"); ax.set_ylabel(\"FVC (mL)\")\n    axes[0].legend(fontsize=7) if n > 1 else axes.legend(fontsize=7)\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"trajectories.png\"), dpi=150); plt.close()\n\n\n# ----------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------\nif __name__ == \"__main__\":\n    df = load_clinical()\n    pairs = build_pairs(df)\n    pairs.to_csv(os.path.join(OUT_DIR, \"training_pairs.csv\"), index=False)\n\n    feats_A, feats_B = FEATS_A, FEATS_B\n    if CT_FEATURES_CSV and os.path.exists(CT_FEATURES_CSV):\n        ct = pd.read_csv(CT_FEATURES_CSV)\n        ct_cols = [c for c in ct.columns if c not in (\"Patient\", \"manufacturer\", \"slice_thickness\", \"n_slices\")]\n        pairs = pairs.merge(ct, on=\"Patient\", how=\"left\")\n        for c in ct_cols:\n            pairs[c] = pairs[c].fillna(pairs[c].median())\n        feats_B = feats_B + ct_cols\n        print(f\"CT features merged: {ct_cols}\")\n\n    print(\"\\n=== Experiment A: baseline-only ===\")\n    resA, hzA, _, _, _ = run_cv(pairs, feats_A, \"A_baseline_only\")\n    print(\"\\n=== Experiment B: history-aware ===\")\n    resB, hzB, oofB, oofqB, sigmaB = run_cv(pairs, feats_B, \"B_history\")\n\n    results = pd.concat([resA, resB])\n    horizon = pd.concat([hzA, hzB])\n    results.to_csv(os.path.join(OUT_DIR, \"cv_results.csv\"), index=False)\n    horizon.to_csv(os.path.join(OUT_DIR, \"horizon_results.csv\"), index=False)\n    print(\"\\n================ CV RESULTS ================\")\n    print(results.round(3).to_string(index=False))\n    print(f\"\\n80% interval coverage (Exp B): {resB.attrs['coverage80']:.3f}\")\n    print(\"\\n============ PER-HORIZON (Ensemble) ============\")\n    print(horizon.round(1).to_string(index=False))\n\n    print(\"\\n=== Rapid progression classifier ===\")\n    rp = rapid_progression(df)\n    print({k: round(v, 3) if isinstance(v, float) else v for k, v in rp.items()})\n    json.dump(rp, open(os.path.join(OUT_DIR, \"rapid_progression_metrics.json\"), \"w\"), indent=2)\n\n    # Fit final models on ALL data (for dashboard) + SHAP\n    Xall = encode(pairs[feats_B]); yall = pairs.target_fvc.values\n    final = lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03,\n                              min_child_samples=10, random_state=SEED, verbose=-1).fit(Xall, yall)\n    qfinal = {a: m.fit(Xall, yall) for a, m in quantile_models().items()}\n    joblib.dump(dict(model=final, quantiles=qfinal, columns=list(Xall.columns), feats=feats_B),\n                os.path.join(OUT_DIR, \"final_model.joblib\"))\n    shap_plots(final, Xall, \"expB\")\n    trajectory_plot(df, pairs, oofB, oofqB)\n    print(f\"\\nAll outputs saved in {OUT_DIR}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:45:18.612631Z","iopub.execute_input":"2026-09-07T19:45:18.613522Z","iopub.status.idle":"2026-09-07T19:46:03.035873Z","shell.execute_reply.started":"2026-09-07T19:45:18.613472Z","shell.execute_reply":"2026-09-07T19:46:03.034764Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"  # conformal correction: widen intervals so 80% of out-of-fold points are covered\n  import numpy as np\n  y = pairs.target_fvc.values\n  err = np.abs(y - oofqB[0.5])\n  half = (oofqB[0.9] - oofqB[0.1]) / 2\n  scale = np.quantile(err / half, 0.80)\n  print(f\"conformal scale factor = {scale:.2f}\")\n  cov = np.mean((y >= oofqB[0.5] - scale*half) & (y <= oofqB[0.5] + scale*half))\n  print(f\"coverage after conformal correction = {cov:.3f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:48:22.179117Z","iopub.execute_input":"2026-09-07T19:48:22.179848Z","iopub.status.idle":"2026-09-07T19:48:22.188533Z","shell.execute_reply.started":"2026-09-07T19:48:22.179814Z","shell.execute_reply":"2026-09-07T19:48:22.187511Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\nos.environ[\"MAX_PATIENTS\"] = \"10\"\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        hu = arr * slope + inter\n        imgs.append(hu)\n    vol = np.stack(imgs)\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm\n    try:\n        px = [float(x) for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [0.7, 0.7]\n    try:\n        st = float(slices[0].SliceThickness)\n        if st <= 0: raise ValueError\n    except Exception:\n        try:\n            st = abs(float(slices[1].ImagePositionPatient[2]) - float(slices[0].ImagePositionPatient[2]))\n        except Exception:\n            st = 1.0\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol):\n    \"\"\"threshold -> remove air touching the border -> keep 2 largest blobs -> fill holes\"\"\"\n    mask = vol < -320                                   # air-like voxels\n    mask = np.stack([segmentation.clear_border(m) for m in mask])   # per-slice: remove air outside the body\n    mask = morphology.binary_opening(mask, morphology.ball(1))\n    labels = measure.label(mask)\n    if labels.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(labels.ravel()); sizes[0] = 0\n    keep = np.argsort(sizes)[::-1][:2]                  # two largest = left + right lung\n    keep = [k for k in keep if sizes[k] > 0.2 * sizes.max()]\n    lung = np.isin(labels, keep)\n    lung = ndimage.binary_closing(lung, iterations=2)\n    lung = ndimage.binary_fill_holes(lung)\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            lung = segment_lungs(vol)\n            feats = extract_features(vol, lung, spacing)\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:50:36.204029Z","iopub.execute_input":"2026-09-07T19:50:36.204407Z","iopub.status.idle":"2026-09-07T19:51:48.926885Z","shell.execute_reply.started":"2026-09-07T19:50:36.204376Z","shell.execute_reply":"2026-09-07T19:51:48.925892Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        hu = arr * slope + inter\n        imgs.append(hu)\n    vol = np.stack(imgs)\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm\n    try:\n        px = [float(x) for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [0.7, 0.7]\n    try:\n        st = float(slices[0].SliceThickness)\n        if st <= 0: raise ValueError\n    except Exception:\n        try:\n            st = abs(float(slices[1].ImagePositionPatient[2]) - float(slices[0].ImagePositionPatient[2]))\n        except Exception:\n            st = 1.0\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol):\n    \"\"\"threshold -> remove air touching the border -> keep 2 largest blobs -> fill holes\"\"\"\n    mask = vol < -320                                   # air-like voxels\n    mask = np.stack([segmentation.clear_border(m) for m in mask])   # per-slice: remove air outside the body\n    mask = morphology.binary_opening(mask, morphology.ball(1))\n    labels = measure.label(mask)\n    if labels.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(labels.ravel()); sizes[0] = 0\n    keep = np.argsort(sizes)[::-1][:2]                  # two largest = left + right lung\n    keep = [k for k in keep if sizes[k] > 0.2 * sizes.max()]\n    lung = np.isin(labels, keep)\n    lung = ndimage.binary_closing(lung, iterations=2)\n    lung = ndimage.binary_fill_holes(lung)\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            lung = segment_lungs(vol)\n            feats = extract_features(vol, lung, spacing)\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:54:36.292331Z","iopub.execute_input":"2026-09-07T19:54:36.292683Z","iopub.status.idle":"2026-09-07T19:55:34.301775Z","shell.execute_reply.started":"2026-09-07T19:54:36.292655Z","shell.execute_reply":"2026-09-07T19:55:34.300652Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\"\"\"\nSTEP 3 - FUSION + ABLATION   (this is your headline experiment)\n----------------------------------------------------------------\nArm A : clinical features only\nArm B : clinical + interpretable CT features\nArm C : CT features only  (+ target week)\nSame patient-grouped folds for all arms, so the comparison is fair.\nOutputs: ablation_results.csv, ablation_horizon.csv, bootstrap test A vs B,\n         shap_summary_armB.png (shows which CT features matter), ct_feature_boxplots.png\nRequires: outputs/training_pairs.csv (from step 1) and outputs/ct_features.csv (from step 2)\n\"\"\"\n\nimport os, warnings\nimport numpy as np, pandas as pd\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom sklearn.ensemble import RandomForestRegressor\nimport joblib\n\nwarnings.filterwarnings(\"ignore\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nSEED, N_FOLDS = 42, 5\nnp.random.seed(SEED)\n\nCLIN = [\"Age\", \"Sex\", \"SmokingStatus\", \"base_fvc\", \"base_pct\", \"base_week\", \"target_week\", \"week_diff\",\n        \"last_fvc\", \"last_week\", \"weeks_since_last\", \"n_prev_visits\", \"prev_fvc_mean\", \"prev_fvc_std\",\n        \"slope_so_far\", \"decline_so_far\"]\nCT = [\"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n      \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\"]\n\n\ndef encode(X):\n    X = X.copy()\n    if \"Sex\" in X: X[\"Sex\"] = (X[\"Sex\"] == \"Male\").astype(int)\n    if \"SmokingStatus\" in X:\n        for s in [\"Ex-smoker\", \"Never smoked\", \"Currently smokes\"]:\n            X[\"smk_\" + s.replace(\" \", \"_\")] = (X[\"SmokingStatus\"] == s).astype(int)\n        X = X.drop(columns=[\"SmokingStatus\"])\n    return X\n\n\ndef laplace_ll(y, p, sigma):\n    sigma = np.maximum(sigma, 70); delta = np.minimum(np.abs(y - p), 1000)\n    return np.mean(-np.sqrt(2) * delta / sigma - np.log(np.sqrt(2) * sigma))\n\n\ndef models():\n    return {\n        \"RandomForest\": RandomForestRegressor(n_estimators=400, min_samples_leaf=5, random_state=SEED, n_jobs=-1),\n        \"XGBoost\": xgb.XGBRegressor(n_estimators=400, max_depth=4, learning_rate=0.03, subsample=0.8,\n                                    colsample_bytree=0.8, random_state=SEED),\n        \"LightGBM\": lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03,\n                                      min_child_samples=10, random_state=SEED, verbose=-1),\n    }\n\n\ndef run_arm(pairs, feats, arm, folds):\n    X = encode(pairs[feats]); y = pairs.target_fvc.values\n    oof = {n: np.zeros(len(y)) for n in models()}\n    q10, q90 = np.zeros(len(y)), np.zeros(len(y))\n    for tr, te in folds:\n        for n, m in models().items():\n            m.fit(X.iloc[tr], y[tr]); oof[n][te] = m.predict(X.iloc[te])\n        for a, arr in ((0.1, q10), (0.9, q90)):\n            qm = lgb.LGBMRegressor(objective=\"quantile\", alpha=a, n_estimators=400, num_leaves=15,\n                                   learning_rate=0.03, min_child_samples=10, random_state=SEED, verbose=-1)\n            qm.fit(X.iloc[tr], y[tr]); arr[te] = qm.predict(X.iloc[te])\n    oof[\"Ensemble\"] = np.mean(list(oof.values()), axis=0)\n    sigma = np.maximum((q90 - q10) / 2, 70)\n    res = [dict(arm=arm, model=n, MAE=mean_absolute_error(y, p), RMSE=np.sqrt(mean_squared_error(y, p)),\n                Laplace=laplace_ll(y, p, sigma), coverage80=np.mean((y >= q10) & (y <= q90)))\n           for n, p in oof.items()]\n    hz = pd.DataFrame({\"arm\": arm, \"bucket\": pairs.horizon_bucket, \"err\": np.abs(y - oof[\"Ensemble\"])})\n    hz = hz.groupby([\"arm\", \"bucket\"]).agg(n=(\"err\", \"size\"), MAE=(\"err\", \"mean\")).reset_index()\n    print(f\"  arm {arm} done\")\n    return pd.DataFrame(res), hz, oof[\"Ensemble\"], X\n\n\ndef paired_bootstrap(pairs, errA, errB, n_boot=2000):\n    \"\"\"bootstrap PATIENTS (not rows) - returns mean MAE difference and 95% CI\"\"\"\n    pid = pairs.Patient.values; uniq = np.unique(pid)\n    idx = {p: np.where(pid == p)[0] for p in uniq}\n    diffs = []\n    rng = np.random.default_rng(SEED)\n    for _ in range(n_boot):\n        samp = rng.choice(uniq, len(uniq), replace=True)\n        rows = np.concatenate([idx[p] for p in samp])\n        diffs.append(errA[rows].mean() - errB[rows].mean())\n    diffs = np.array(diffs)\n    return diffs.mean(), np.percentile(diffs, 2.5), np.percentile(diffs, 97.5), np.mean(diffs <= 0)\n\n\nif __name__ == \"__main__\":\n    pairs = pd.read_csv(os.path.join(OUT_DIR, \"training_pairs.csv\"))\n    ct = pd.read_csv(os.path.join(OUT_DIR, \"ct_features.csv\"))\n    n_before = pairs.Patient.nunique()\n    pairs = pairs.merge(ct[[\"Patient\"] + CT], on=\"Patient\", how=\"inner\")   # only patients with CT\n    print(f\"Patients with clinical data: {n_before}, with CT features too: {pairs.Patient.nunique()}\")\n\n    folds = list(GroupKFold(N_FOLDS).split(pairs, pairs.target_fvc, pairs.Patient))\n    print(\"Running arms...\")\n    resA, hzA, ensA, _ = run_arm(pairs, CLIN, \"A_clinical\", folds)\n    resB, hzB, ensB, XB = run_arm(pairs, CLIN + CT, \"B_clinical+CT\", folds)\n    resC, hzC, ensC, _ = run_arm(pairs, CT + [\"Age\", \"Sex\", \"target_week\", \"week_diff\"], \"C_CT_only\", folds)\n\n    res = pd.concat([resA, resB, resC]); hz = pd.concat([hzA, hzB, hzC])\n    res.to_csv(os.path.join(OUT_DIR, \"ablation_results.csv\"), index=False)\n    hz.to_csv(os.path.join(OUT_DIR, \"ablation_horizon.csv\"), index=False)\n    print(\"\\n================ ABLATION RESULTS ================\")\n    print(res.round(3).to_string(index=False))\n    print(\"\\n============ PER-HORIZON MAE (Ensemble) ============\")\n    print(hz.pivot(index=\"bucket\", columns=\"arm\", values=\"MAE\").round(1))\n\n    y = pairs.target_fvc.values\n    d, lo, hi, p_worse = paired_bootstrap(pairs, np.abs(y - ensA), np.abs(y - ensB))\n    print(f\"\\nPaired bootstrap (patients resampled): MAE(A) - MAE(B) = {d:.1f} mL, 95% CI [{lo:.1f}, {hi:.1f}]\")\n    print(\"  -> positive = CT features help.  CI containing 0 = no significant difference.\")\n    with open(os.path.join(OUT_DIR, \"bootstrap_A_vs_B.txt\"), \"w\") as f:\n        f.write(f\"MAE(A)-MAE(B)={d:.2f} mL, 95% CI [{lo:.2f},{hi:.2f}], P(B not better)={p_worse:.3f}\\n\")\n\n    # SHAP for arm B (which CT features matter?)\n    try:\n        import shap\n        m = lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03, min_child_samples=10,\n                              random_state=SEED, verbose=-1).fit(XB, y)\n        sv = shap.TreeExplainer(m).shap_values(XB)\n        plt.figure(); shap.summary_plot(sv, XB, show=False, max_display=20)\n        plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"shap_summary_armB.png\"), dpi=150); plt.close()\n        imp = pd.Series(np.abs(sv).mean(0), index=XB.columns).sort_values(ascending=False)\n        imp.to_csv(os.path.join(OUT_DIR, \"shap_importance_armB.csv\"))\n        print(\"\\nTop features by mean |SHAP| (arm B):\"); print(imp.head(12).round(2))\n        joblib.dump(dict(model=m, columns=list(XB.columns), feats=CLIN + CT), os.path.join(OUT_DIR, \"final_model_armB.joblib\"))\n    except Exception as e:\n        print(\"SHAP skipped:\", e)\n\n    # CT feature vs decline slope plot (does fibrosis proxy relate to progression?)\n    per_pat = pairs.groupby(\"Patient\").agg(slope=(\"slope_so_far\", \"last\"), haa=(\"pct_high_att\", \"first\"),\n                                           kurt=(\"hu_kurtosis\", \"first\"), vol=(\"lung_volume_ml\", \"first\"))\n    fig, ax = plt.subplots(1, 3, figsize=(11, 3.3))\n    for a, col, lab in zip(ax, [\"haa\", \"kurt\", \"vol\"], [\"% high attenuation (fibrosis proxy)\", \"HU kurtosis\", \"lung volume (mL)\"]):\n        a.scatter(per_pat[col], per_pat.slope, s=12, alpha=.7); a.set_xlabel(lab); a.set_ylabel(\"FVC slope (mL/week)\")\n        r = np.corrcoef(per_pat[col], per_pat.slope)[0, 1]; a.set_title(f\"r = {r:.2f}\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"ct_feature_boxplots.png\"), dpi=150); plt.close()\n    print(f\"\\nSaved ablation outputs to {OUT_DIR}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T19:58:26.951983Z","iopub.execute_input":"2026-09-07T19:58:26.952324Z","iopub.status.idle":"2026-09-07T19:58:45.337604Z","shell.execute_reply.started":"2026-09-07T19:58:26.952298Z","shell.execute_reply":"2026-09-07T19:58:45.336549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   !pip install pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg -q","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:01:51.696339Z","iopub.execute_input":"2026-09-07T20:01:51.696667Z","iopub.status.idle":"2026-09-07T20:01:55.938731Z","shell.execute_reply.started":"2026-09-07T20:01:51.696631Z","shell.execute_reply":"2026-09-07T20:01:55.937365Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   import os\n   os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n   os.environ[\"MAX_PATIENTS\"] = \"10\"\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: if there is no air anywhere, the intercept was missing/wrong -> shift to HU\n    if np.percentile(vol, 5) > -300:\n        vol = vol - 1024\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > -500                                     # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < -320) & body                             # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            lung = segment_lungs(vol)\n            feats = extract_features(vol, lung, spacing)\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:02:40.056881Z","iopub.execute_input":"2026-09-07T20:02:40.057283Z","iopub.status.idle":"2026-09-07T20:03:02.023496Z","shell.execute_reply.started":"2026-09-07T20:02:40.057244Z","shell.execute_reply":"2026-09-07T20:03:02.022226Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   import os\n   os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n   os.environ[\"MAX_PATIENTS\"] = \"10\"\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, -1024, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:07:12.510316Z","iopub.execute_input":"2026-09-07T20:07:12.510643Z","iopub.status.idle":"2026-09-07T20:07:52.305445Z","shell.execute_reply.started":"2026-09-07T20:07:12.510615Z","shell.execute_reply":"2026-09-07T20:07:52.304403Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   !pip install pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg python-gdcm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:11:42.219611Z","iopub.execute_input":"2026-09-07T20:11:42.219938Z","iopub.status.idle":"2026-09-07T20:11:49.587412Z","shell.execute_reply.started":"2026-09-07T20:11:42.219906Z","shell.execute_reply":"2026-09-07T20:11:49.585761Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   import os\n   os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n   os.environ[\"MAX_PATIENTS\"] = \"10\"\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, -1024, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:12:26.330257Z","iopub.execute_input":"2026-09-07T20:12:26.330631Z","iopub.status.idle":"2026-09-07T20:12:49.49239Z","shell.execute_reply.started":"2026-09-07T20:12:26.330595Z","shell.execute_reply":"2026-09-07T20:12:49.491056Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"   import os\n   os.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, -1024, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:13:47.134819Z","iopub.execute_input":"2026-09-07T20:13:47.135645Z","iopub.status.idle":"2026-09-07T20:14:10.688894Z","shell.execute_reply.started":"2026-09-07T20:13:47.135608Z","shell.execute_reply":"2026-09-07T20:14:10.686922Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg python-gdcm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:17:13.682223Z","iopub.execute_input":"2026-09-07T20:17:13.682576Z","iopub.status.idle":"2026-09-07T20:17:20.920561Z","shell.execute_reply.started":"2026-09-07T20:17:13.682539Z","shell.execute_reply":"2026-09-07T20:17:20.9193Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\nos.environ[\"MAX_PATIENTS\"] = \"10\"\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, 1024, -1024, 2048, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:18:36.065458Z","iopub.execute_input":"2026-09-07T20:18:36.066028Z","iopub.status.idle":"2026-09-07T20:19:22.793523Z","shell.execute_reply.started":"2026-09-07T20:18:36.06599Z","shell.execute_reply":"2026-09-07T20:19:22.792161Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, 1024, -1024, 2048, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:22:21.556256Z","iopub.execute_input":"2026-09-07T20:22:21.556997Z","iopub.status.idle":"2026-09-07T20:22:46.752825Z","shell.execute_reply.started":"2026-09-07T20:22:21.55696Z","shell.execute_reply":"2026-09-07T20:22:46.751751Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\nos.environ.pop(\"MAX_PATIENTS\", None)\n\n\"\"\"\nSTEP 2 - CT FEATURE EXTRACTION  (interpretable, no deep learning)\n------------------------------------------------------------------\nFor every patient folder of DICOM slices:\n  1. read slices, sort by position, convert to Hounsfield Units (HU)\n  2. segment the lungs with a simple threshold + connected components\n  3. compute 9 interpretable numbers (lung volume, HU stats, fibrosis proxies)\n  4. save one row per patient into ct_features.csv\n  5. save a few lung-mask pictures for your slides\n\nRun time: ~3-10 s per patient on Kaggle CPU. 176 patients ≈ 10-25 min.\nSet MAX_PATIENTS = 10 first to check it works, then set to None for all.\n\"\"\"\n\nimport os, glob, warnings, traceback\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy import ndimage\nfrom scipy.stats import skew, kurtosis\nfrom skimage import measure, morphology, segmentation\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings(\"ignore\")\n\n# ============================ CONFIG ============================\nDATA_DIR = os.environ.get(\"DATA_DIR\", \"/kaggle/input/osic-pulmonary-fibrosis-progression\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nMAX_PATIENTS = int(os.environ[\"MAX_PATIENTS\"]) if os.environ.get(\"MAX_PATIENTS\") else None\nSLICE_STEP = 2            # use every 2nd slice to halve run-time (set 1 for all slices)\nN_EXAMPLE_IMAGES = 4      # how many mask pictures to save for slides\n# ================================================================\nos.makedirs(OUT_DIR, exist_ok=True)\n\n\ndef load_scan(folder):\n    files = glob.glob(os.path.join(folder, \"*.dcm\"))\n    slices = []\n    for f in files:\n        try:\n            ds = pydicom.dcmread(f, force=True)\n            if not hasattr(ds, \"PixelData\"):\n                continue\n            slices.append(ds)\n        except Exception:\n            pass\n    if len(slices) < 5:\n        raise ValueError(\"too few readable slices\")\n    try:\n        slices.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    except Exception:\n        slices.sort(key=lambda s: int(getattr(s, \"InstanceNumber\", 0)))\n    return slices\n\n\ndef to_hu(slices):\n    \"\"\"stack slices into a 3D array of Hounsfield units (downsampled 2x in-plane for speed)\"\"\"\n    imgs = []\n    shape = slices[0].pixel_array.shape\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        if arr.shape != shape:                 # skip odd-sized slices\n            continue\n        slope = float(getattr(s, \"RescaleSlope\", 1) or 1)\n        inter = float(getattr(s, \"RescaleIntercept\", -1024) or -1024)\n        imgs.append((arr * slope + inter)[::2, ::2])\n    vol = np.stack(imgs)\n    # sanity: some OSIC scans store raw values (air ~ 0) with a wrong intercept.\n    # pick the shift (0 / -1024 / -2048) that gives the most \"air-like\" voxels, ignoring padding.\n    def hu_score(v):\n        valid = v[v > -1400]\n        if valid.size == 0: return 0\n        air = np.mean((valid > -1100) & (valid < -700))      # lungs / air\n        tissue = np.mean((valid > -150) & (valid < 300))     # soft tissue\n        return min(air, tissue)                              # a real HU scan has BOTH peaks\n    best = max([0, 1024, -1024, 2048, -2048], key=lambda sh: hu_score(vol + sh))\n    if best != 0:\n        vol = vol + best\n    vol = np.clip(vol, -1500, 1500)\n    # voxel size in mm (in-plane x2 because of downsampling, z x SLICE_STEP)\n    try:\n        px = [float(x) * 2 for x in slices[0].PixelSpacing]\n    except Exception:\n        px = [1.4, 1.4]\n    try:\n        zs = [float(s.ImagePositionPatient[2]) for s in slices]\n        st = float(np.median(np.abs(np.diff(zs))))\n        if not (0.1 < st < 50): raise ValueError\n    except Exception:\n        try:\n            st = float(slices[0].SliceThickness) * SLICE_STEP\n        except Exception:\n            st = 1.0 * SLICE_STEP\n    return vol, (st, px[0], px[1])\n\n\ndef segment_lungs(vol, body_thr=-500, lung_thr=-320):\n    \"\"\"body mask first (ignores scanner padding), then air INSIDE the body = lungs\"\"\"\n    lung = np.zeros(vol.shape, bool)\n    for i in range(vol.shape[0]):\n        sl = vol[i]\n        body = sl > body_thr                                 # tissue + bone\n        lab = measure.label(body)\n        if lab.max() == 0:\n            continue\n        sizes = np.bincount(lab.ravel()); sizes[0] = 0\n        body = lab == sizes.argmax()                          # largest blob = patient body\n        body = ndimage.binary_fill_holes(body)                # fill lungs into the body outline\n        cand = (sl < lung_thr) & body                         # air inside the body\n        cand = morphology.binary_opening(cand, morphology.disk(2))\n        lung[i] = cand\n    lab3 = measure.label(lung)\n    if lab3.max() == 0:\n        raise ValueError(\"no lung candidate found\")\n    sizes = np.bincount(lab3.ravel()); sizes[0] = 0\n    order = np.argsort(sizes)[::-1]\n    keep = [order[0]]                                         # largest component\n    if len(order) > 1 and sizes[order[1]] > 0.15 * sizes[order[0]]:\n        keep.append(order[1])                                 # second lung if separate\n    lung = np.isin(lab3, keep)\n    for i in range(lung.shape[0]):                            # fill vessels inside lungs\n        lung[i] = ndimage.binary_fill_holes(lung[i])\n    return lung\n\n\ndef extract_features(vol, lung, spacing):\n    hu = vol[lung]\n    if hu.size < 1000:\n        raise ValueError(\"lung mask too small\")\n    voxel_ml = spacing[0] * spacing[1] * spacing[2] / 1000.0\n    vol_ml = hu.size * voxel_ml\n    if not (800 < vol_ml < 9000):\n        raise ValueError(f\"implausible lung volume {vol_ml:.0f} mL (segmentation/spacing problem)\")\n    return dict(\n        lung_volume_ml=hu.size * voxel_ml,\n        hu_mean=float(hu.mean()),\n        hu_std=float(hu.std()),\n        hu_skew=float(skew(hu)),\n        hu_kurtosis=float(kurtosis(hu)),\n        pct_high_att=float(np.mean((hu >= -600) & (hu <= -250)) * 100),   # fibrosis proxy (HAA%)\n        pct_low_att=float(np.mean(hu < -950) * 100),                       # emphysema proxy\n        pct_dense=float(np.mean(hu > -500) * 100),                         # consolidation / dense tissue\n        pct_normal=float(np.mean((hu >= -950) & (hu < -700)) * 100),       # normally aerated lung\n    )\n\n\ndef save_example(vol, lung, pid, spacing):\n    z = vol.shape[0] // 2\n    fig, ax = plt.subplots(1, 2, figsize=(7, 3.5))\n    ax[0].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400); ax[0].set_title(f\"{pid[:14]} - CT (HU)\")\n    ax[1].imshow(vol[z], cmap=\"gray\", vmin=-1000, vmax=400)\n    ax[1].imshow(np.ma.masked_where(~lung[z], lung[z]), cmap=\"autumn\", alpha=.45); ax[1].set_title(\"lung mask\")\n    for a in ax: a.axis(\"off\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, f\"ctmask_{pid}.png\"), dpi=130); plt.close()\n\n\nif __name__ == \"__main__\":\n    folders = sorted(glob.glob(os.path.join(DATA_DIR, \"train\", \"*\")))\n    if MAX_PATIENTS: folders = folders[:MAX_PATIENTS]\n    rows, failed, n_img = [], [], 0\n    for i, folder in enumerate(folders):\n        pid = os.path.basename(folder)\n        try:\n            slices = load_scan(folder)[::SLICE_STEP]\n            vol, spacing = to_hu(slices)\n            feats, last_err = None, None\n            for body_thr, lung_thr in [(-500, -320), (-700, -450), (-300, -200), (-800, -600)]:\n                try:\n                    lung = segment_lungs(vol, body_thr, lung_thr)\n                    feats = extract_features(vol, lung, spacing)\n                    break\n                except Exception as e:\n                    last_err = e\n            if feats is None:\n                pcts = np.percentile(vol, [1, 25, 50, 75, 99]).round(0)\n                raise ValueError(f\"{last_err} | HU percentiles {pcts.tolist()}\")\n            feats.update(Patient=pid, n_slices=len(slices) * SLICE_STEP,\n                         slice_thickness=spacing[0],\n                         manufacturer=str(getattr(slices[0], \"Manufacturer\", \"unknown\")))\n            rows.append(feats)\n            if n_img < N_EXAMPLE_IMAGES:\n                save_example(vol, lung, pid, spacing); n_img += 1\n            print(f\"[{i + 1}/{len(folders)}] {pid}  vol={feats['lung_volume_ml']:.0f} mL  \"\n                  f\"HAA={feats['pct_high_att']:.1f}%  kurt={feats['hu_kurtosis']:.2f}\")\n        except Exception as e:\n            failed.append((pid, str(e)))\n            print(f\"[{i + 1}/{len(folders)}] {pid}  FAILED: {e}\")\n\n    df = pd.DataFrame(rows)\n    if len(df) == 0:\n        raise SystemExit(\"No patient succeeded - check DATA_DIR and ct_failed reasons above\")\n    df = df[[\"Patient\", \"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n             \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\",\n             \"n_slices\", \"slice_thickness\", \"manufacturer\"]]\n    df.to_csv(os.path.join(OUT_DIR, \"ct_features.csv\"), index=False)\n    pd.DataFrame(failed, columns=[\"Patient\", \"reason\"]).to_csv(os.path.join(OUT_DIR, \"ct_failed.csv\"), index=False)\n    print(f\"\\nDone. {len(df)} patients extracted, {len(failed)} failed. Saved to {OUT_DIR}/ct_features.csv\")\n    print(df.describe().round(2).T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:23:59.741482Z","iopub.execute_input":"2026-09-07T20:23:59.741815Z","iopub.status.idle":"2026-09-07T20:36:02.18384Z","shell.execute_reply.started":"2026-09-07T20:23:59.741789Z","shell.execute_reply":"2026-09-07T20:36:02.182628Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\"\"\"\nSTEP 3 - FUSION + ABLATION   (this is your headline experiment)\n----------------------------------------------------------------\nArm A : clinical features only\nArm B : clinical + interpretable CT features\nArm C : CT features only  (+ target week)\nSame patient-grouped folds for all arms, so the comparison is fair.\nOutputs: ablation_results.csv, ablation_horizon.csv, bootstrap test A vs B,\n         shap_summary_armB.png (shows which CT features matter), ct_feature_boxplots.png\nRequires: outputs/training_pairs.csv (from step 1) and outputs/ct_features.csv (from step 2)\n\"\"\"\n\nimport os, warnings\nimport numpy as np, pandas as pd\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom sklearn.ensemble import RandomForestRegressor\nimport joblib\n\nwarnings.filterwarnings(\"ignore\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nSEED, N_FOLDS = 42, 5\nnp.random.seed(SEED)\n\nCLIN = [\"Age\", \"Sex\", \"SmokingStatus\", \"base_fvc\", \"base_pct\", \"base_week\", \"target_week\", \"week_diff\",\n        \"last_fvc\", \"last_week\", \"weeks_since_last\", \"n_prev_visits\", \"prev_fvc_mean\", \"prev_fvc_std\",\n        \"slope_so_far\", \"decline_so_far\"]\nCT = [\"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n      \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\"]\n\n\ndef encode(X):\n    X = X.copy()\n    if \"Sex\" in X: X[\"Sex\"] = (X[\"Sex\"] == \"Male\").astype(int)\n    if \"SmokingStatus\" in X:\n        for s in [\"Ex-smoker\", \"Never smoked\", \"Currently smokes\"]:\n            X[\"smk_\" + s.replace(\" \", \"_\")] = (X[\"SmokingStatus\"] == s).astype(int)\n        X = X.drop(columns=[\"SmokingStatus\"])\n    return X\n\n\ndef laplace_ll(y, p, sigma):\n    sigma = np.maximum(sigma, 70); delta = np.minimum(np.abs(y - p), 1000)\n    return np.mean(-np.sqrt(2) * delta / sigma - np.log(np.sqrt(2) * sigma))\n\n\ndef models():\n    return {\n        \"RandomForest\": RandomForestRegressor(n_estimators=400, min_samples_leaf=5, random_state=SEED, n_jobs=-1),\n        \"XGBoost\": xgb.XGBRegressor(n_estimators=400, max_depth=4, learning_rate=0.03, subsample=0.8,\n                                    colsample_bytree=0.8, random_state=SEED),\n        \"LightGBM\": lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03,\n                                      min_child_samples=10, random_state=SEED, verbose=-1),\n    }\n\n\ndef run_arm(pairs, feats, arm, folds):\n    X = encode(pairs[feats]); y = pairs.target_fvc.values\n    oof = {n: np.zeros(len(y)) for n in models()}\n    q10, q90 = np.zeros(len(y)), np.zeros(len(y))\n    for tr, te in folds:\n        for n, m in models().items():\n            m.fit(X.iloc[tr], y[tr]); oof[n][te] = m.predict(X.iloc[te])\n        for a, arr in ((0.1, q10), (0.9, q90)):\n            qm = lgb.LGBMRegressor(objective=\"quantile\", alpha=a, n_estimators=400, num_leaves=15,\n                                   learning_rate=0.03, min_child_samples=10, random_state=SEED, verbose=-1)\n            qm.fit(X.iloc[tr], y[tr]); arr[te] = qm.predict(X.iloc[te])\n    oof[\"Ensemble\"] = np.mean(list(oof.values()), axis=0)\n    sigma = np.maximum((q90 - q10) / 2, 70)\n    res = [dict(arm=arm, model=n, MAE=mean_absolute_error(y, p), RMSE=np.sqrt(mean_squared_error(y, p)),\n                Laplace=laplace_ll(y, p, sigma), coverage80=np.mean((y >= q10) & (y <= q90)))\n           for n, p in oof.items()]\n    hz = pd.DataFrame({\"arm\": arm, \"bucket\": pairs.horizon_bucket, \"err\": np.abs(y - oof[\"Ensemble\"])})\n    hz = hz.groupby([\"arm\", \"bucket\"]).agg(n=(\"err\", \"size\"), MAE=(\"err\", \"mean\")).reset_index()\n    print(f\"  arm {arm} done\")\n    return pd.DataFrame(res), hz, oof[\"Ensemble\"], X\n\n\ndef paired_bootstrap(pairs, errA, errB, n_boot=2000):\n    \"\"\"bootstrap PATIENTS (not rows) - returns mean MAE difference and 95% CI\"\"\"\n    pid = pairs.Patient.values; uniq = np.unique(pid)\n    idx = {p: np.where(pid == p)[0] for p in uniq}\n    diffs = []\n    rng = np.random.default_rng(SEED)\n    for _ in range(n_boot):\n        samp = rng.choice(uniq, len(uniq), replace=True)\n        rows = np.concatenate([idx[p] for p in samp])\n        diffs.append(errA[rows].mean() - errB[rows].mean())\n    diffs = np.array(diffs)\n    return diffs.mean(), np.percentile(diffs, 2.5), np.percentile(diffs, 97.5), np.mean(diffs <= 0)\n\n\nif __name__ == \"__main__\":\n    pairs = pd.read_csv(os.path.join(OUT_DIR, \"training_pairs.csv\"))\n    ct = pd.read_csv(os.path.join(OUT_DIR, \"ct_features.csv\"))\n    n_before = pairs.Patient.nunique()\n    pairs = pairs.merge(ct[[\"Patient\"] + CT], on=\"Patient\", how=\"inner\")   # only patients with CT\n    print(f\"Patients with clinical data: {n_before}, with CT features too: {pairs.Patient.nunique()}\")\n\n    folds = list(GroupKFold(N_FOLDS).split(pairs, pairs.target_fvc, pairs.Patient))\n    print(\"Running arms...\")\n    resA, hzA, ensA, _ = run_arm(pairs, CLIN, \"A_clinical\", folds)\n    resB, hzB, ensB, XB = run_arm(pairs, CLIN + CT, \"B_clinical+CT\", folds)\n    resC, hzC, ensC, _ = run_arm(pairs, CT + [\"Age\", \"Sex\", \"target_week\", \"week_diff\"], \"C_CT_only\", folds)\n\n    res = pd.concat([resA, resB, resC]); hz = pd.concat([hzA, hzB, hzC])\n    res.to_csv(os.path.join(OUT_DIR, \"ablation_results.csv\"), index=False)\n    hz.to_csv(os.path.join(OUT_DIR, \"ablation_horizon.csv\"), index=False)\n    print(\"\\n================ ABLATION RESULTS ================\")\n    print(res.round(3).to_string(index=False))\n    print(\"\\n============ PER-HORIZON MAE (Ensemble) ============\")\n    print(hz.pivot(index=\"bucket\", columns=\"arm\", values=\"MAE\").round(1))\n\n    y = pairs.target_fvc.values\n    d, lo, hi, p_worse = paired_bootstrap(pairs, np.abs(y - ensA), np.abs(y - ensB))\n    print(f\"\\nPaired bootstrap (patients resampled): MAE(A) - MAE(B) = {d:.1f} mL, 95% CI [{lo:.1f}, {hi:.1f}]\")\n    print(\"  -> positive = CT features help.  CI containing 0 = no significant difference.\")\n    with open(os.path.join(OUT_DIR, \"bootstrap_A_vs_B.txt\"), \"w\") as f:\n        f.write(f\"MAE(A)-MAE(B)={d:.2f} mL, 95% CI [{lo:.2f},{hi:.2f}], P(B not better)={p_worse:.3f}\\n\")\n\n    # SHAP for arm B (which CT features matter?)\n    try:\n        import shap\n        m = lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03, min_child_samples=10,\n                              random_state=SEED, verbose=-1).fit(XB, y)\n        sv = shap.TreeExplainer(m).shap_values(XB)\n        plt.figure(); shap.summary_plot(sv, XB, show=False, max_display=20)\n        plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"shap_summary_armB.png\"), dpi=150); plt.close()\n        imp = pd.Series(np.abs(sv).mean(0), index=XB.columns).sort_values(ascending=False)\n        imp.to_csv(os.path.join(OUT_DIR, \"shap_importance_armB.csv\"))\n        print(\"\\nTop features by mean |SHAP| (arm B):\"); print(imp.head(12).round(2))\n        joblib.dump(dict(model=m, columns=list(XB.columns), feats=CLIN + CT), os.path.join(OUT_DIR, \"final_model_armB.joblib\"))\n    except Exception as e:\n        print(\"SHAP skipped:\", e)\n\n    # CT feature vs decline slope plot (does fibrosis proxy relate to progression?)\n    per_pat = pairs.groupby(\"Patient\").agg(slope=(\"slope_so_far\", \"last\"), haa=(\"pct_high_att\", \"first\"),\n                                           kurt=(\"hu_kurtosis\", \"first\"), vol=(\"lung_volume_ml\", \"first\"))\n    fig, ax = plt.subplots(1, 3, figsize=(11, 3.3))\n    for a, col, lab in zip(ax, [\"haa\", \"kurt\", \"vol\"], [\"% high attenuation (fibrosis proxy)\", \"HU kurtosis\", \"lung volume (mL)\"]):\n        a.scatter(per_pat[col], per_pat.slope, s=12, alpha=.7); a.set_xlabel(lab); a.set_ylabel(\"FVC slope (mL/week)\")\n        r = np.corrcoef(per_pat[col], per_pat.slope)[0, 1]; a.set_title(f\"r = {r:.2f}\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"ct_feature_boxplots.png\"), dpi=150); plt.close()\n    print(f\"\\nSaved ablation outputs to {OUT_DIR}\")\n","metadata":{"trusted":true,"execution":{"execution_failed":"2026-09-07T20:45:14.937Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"DATA_DIR\"] = \"/kaggle/input/competitions/osic-pulmonary-fibrosis-progression\"\n\"\"\"\nSTEP 3 - FUSION + ABLATION   (this is your headline experiment)\n----------------------------------------------------------------\nArm A : clinical features only\nArm B : clinical + interpretable CT features\nArm C : CT features only  (+ target week)\nSame patient-grouped folds for all arms, so the comparison is fair.\nOutputs: ablation_results.csv, ablation_horizon.csv, bootstrap test A vs B,\n         shap_summary_armB.png (shows which CT features matter), ct_feature_boxplots.png\nRequires: outputs/training_pairs.csv (from step 1) and outputs/ct_features.csv (from step 2)\n\"\"\"\n\nimport os, warnings\nimport numpy as np, pandas as pd\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom sklearn.ensemble import RandomForestRegressor\nimport joblib\n\nwarnings.filterwarnings(\"ignore\")\nOUT_DIR = os.environ.get(\"OUT_DIR\", \"/kaggle/working/outputs\")\nSEED, N_FOLDS = 42, 5\nnp.random.seed(SEED)\n\nCLIN = [\"Age\", \"Sex\", \"SmokingStatus\", \"base_fvc\", \"base_pct\", \"base_week\", \"target_week\", \"week_diff\",\n        \"last_fvc\", \"last_week\", \"weeks_since_last\", \"n_prev_visits\", \"prev_fvc_mean\", \"prev_fvc_std\",\n        \"slope_so_far\", \"decline_so_far\"]\nCT = [\"lung_volume_ml\", \"hu_mean\", \"hu_std\", \"hu_skew\", \"hu_kurtosis\",\n      \"pct_high_att\", \"pct_low_att\", \"pct_dense\", \"pct_normal\"]\n\n\ndef encode(X):\n    X = X.copy()\n    if \"Sex\" in X: X[\"Sex\"] = (X[\"Sex\"] == \"Male\").astype(int)\n    if \"SmokingStatus\" in X:\n        for s in [\"Ex-smoker\", \"Never smoked\", \"Currently smokes\"]:\n            X[\"smk_\" + s.replace(\" \", \"_\")] = (X[\"SmokingStatus\"] == s).astype(int)\n        X = X.drop(columns=[\"SmokingStatus\"])\n    return X\n\n\ndef laplace_ll(y, p, sigma):\n    sigma = np.maximum(sigma, 70); delta = np.minimum(np.abs(y - p), 1000)\n    return np.mean(-np.sqrt(2) * delta / sigma - np.log(np.sqrt(2) * sigma))\n\n\ndef models():\n    return {\n        \"RandomForest\": RandomForestRegressor(n_estimators=400, min_samples_leaf=5, random_state=SEED, n_jobs=-1),\n        \"XGBoost\": xgb.XGBRegressor(n_estimators=400, max_depth=4, learning_rate=0.03, subsample=0.8,\n                                    colsample_bytree=0.8, random_state=SEED),\n        \"LightGBM\": lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03,\n                                      min_child_samples=10, random_state=SEED, verbose=-1),\n    }\n\n\ndef run_arm(pairs, feats, arm, folds):\n    X = encode(pairs[feats]); y = pairs.target_fvc.values\n    oof = {n: np.zeros(len(y)) for n in models()}\n    q10, q90 = np.zeros(len(y)), np.zeros(len(y))\n    for tr, te in folds:\n        for n, m in models().items():\n            m.fit(X.iloc[tr], y[tr]); oof[n][te] = m.predict(X.iloc[te])\n        for a, arr in ((0.1, q10), (0.9, q90)):\n            qm = lgb.LGBMRegressor(objective=\"quantile\", alpha=a, n_estimators=400, num_leaves=15,\n                                   learning_rate=0.03, min_child_samples=10, random_state=SEED, verbose=-1)\n            qm.fit(X.iloc[tr], y[tr]); arr[te] = qm.predict(X.iloc[te])\n    oof[\"Ensemble\"] = np.mean(list(oof.values()), axis=0)\n    sigma = np.maximum((q90 - q10) / 2, 70)\n    res = [dict(arm=arm, model=n, MAE=mean_absolute_error(y, p), RMSE=np.sqrt(mean_squared_error(y, p)),\n                Laplace=laplace_ll(y, p, sigma), coverage80=np.mean((y >= q10) & (y <= q90)))\n           for n, p in oof.items()]\n    hz = pd.DataFrame({\"arm\": arm, \"bucket\": pairs.horizon_bucket, \"err\": np.abs(y - oof[\"Ensemble\"])})\n    hz = hz.groupby([\"arm\", \"bucket\"]).agg(n=(\"err\", \"size\"), MAE=(\"err\", \"mean\")).reset_index()\n    print(f\"  arm {arm} done\")\n    return pd.DataFrame(res), hz, oof[\"Ensemble\"], X\n\n\ndef paired_bootstrap(pairs, errA, errB, n_boot=2000):\n    \"\"\"bootstrap PATIENTS (not rows) - returns mean MAE difference and 95% CI\"\"\"\n    pid = pairs.Patient.values; uniq = np.unique(pid)\n    idx = {p: np.where(pid == p)[0] for p in uniq}\n    diffs = []\n    rng = np.random.default_rng(SEED)\n    for _ in range(n_boot):\n        samp = rng.choice(uniq, len(uniq), replace=True)\n        rows = np.concatenate([idx[p] for p in samp])\n        diffs.append(errA[rows].mean() - errB[rows].mean())\n    diffs = np.array(diffs)\n    return diffs.mean(), np.percentile(diffs, 2.5), np.percentile(diffs, 97.5), np.mean(diffs <= 0)\n\n\nif __name__ == \"__main__\":\n    pairs = pd.read_csv(os.path.join(OUT_DIR, \"training_pairs.csv\"))\n    ct = pd.read_csv(os.path.join(OUT_DIR, \"ct_features.csv\"))\n    n_before = pairs.Patient.nunique()\n    pairs = pairs.merge(ct[[\"Patient\"] + CT], on=\"Patient\", how=\"inner\")   # only patients with CT\n    print(f\"Patients with clinical data: {n_before}, with CT features too: {pairs.Patient.nunique()}\")\n\n    folds = list(GroupKFold(N_FOLDS).split(pairs, pairs.target_fvc, pairs.Patient))\n    print(\"Running arms...\")\n    resA, hzA, ensA, _ = run_arm(pairs, CLIN, \"A_clinical\", folds)\n    resB, hzB, ensB, XB = run_arm(pairs, CLIN + CT, \"B_clinical+CT\", folds)\n    resC, hzC, ensC, _ = run_arm(pairs, CT + [\"Age\", \"Sex\", \"target_week\", \"week_diff\"], \"C_CT_only\", folds)\n\n    res = pd.concat([resA, resB, resC]); hz = pd.concat([hzA, hzB, hzC])\n    res.to_csv(os.path.join(OUT_DIR, \"ablation_results.csv\"), index=False)\n    hz.to_csv(os.path.join(OUT_DIR, \"ablation_horizon.csv\"), index=False)\n    print(\"\\n================ ABLATION RESULTS ================\")\n    print(res.round(3).to_string(index=False))\n    print(\"\\n============ PER-HORIZON MAE (Ensemble) ============\")\n    print(hz.pivot(index=\"bucket\", columns=\"arm\", values=\"MAE\").round(1))\n\n    y = pairs.target_fvc.values\n    d, lo, hi, p_worse = paired_bootstrap(pairs, np.abs(y - ensA), np.abs(y - ensB))\n    print(f\"\\nPaired bootstrap (patients resampled): MAE(A) - MAE(B) = {d:.1f} mL, 95% CI [{lo:.1f}, {hi:.1f}]\")\n    print(\"  -> positive = CT features help.  CI containing 0 = no significant difference.\")\n    with open(os.path.join(OUT_DIR, \"bootstrap_A_vs_B.txt\"), \"w\") as f:\n        f.write(f\"MAE(A)-MAE(B)={d:.2f} mL, 95% CI [{lo:.2f},{hi:.2f}], P(B not better)={p_worse:.3f}\\n\")\n\n    # SHAP for arm B (which CT features matter?)\n    try:\n        import shap\n        m = lgb.LGBMRegressor(n_estimators=400, num_leaves=15, learning_rate=0.03, min_child_samples=10,\n                              random_state=SEED, verbose=-1).fit(XB, y)\n        sv = shap.TreeExplainer(m).shap_values(XB)\n        plt.figure(); shap.summary_plot(sv, XB, show=False, max_display=20)\n        plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"shap_summary_armB.png\"), dpi=150); plt.close()\n        imp = pd.Series(np.abs(sv).mean(0), index=XB.columns).sort_values(ascending=False)\n        imp.to_csv(os.path.join(OUT_DIR, \"shap_importance_armB.csv\"))\n        print(\"\\nTop features by mean |SHAP| (arm B):\"); print(imp.head(12).round(2))\n        joblib.dump(dict(model=m, columns=list(XB.columns), feats=CLIN + CT), os.path.join(OUT_DIR, \"final_model_armB.joblib\"))\n    except Exception as e:\n        print(\"SHAP skipped:\", e)\n\n    # CT feature vs decline slope plot (does fibrosis proxy relate to progression?)\n    per_pat = pairs.groupby(\"Patient\").agg(slope=(\"slope_so_far\", \"last\"), haa=(\"pct_high_att\", \"first\"),\n                                           kurt=(\"hu_kurtosis\", \"first\"), vol=(\"lung_volume_ml\", \"first\"))\n    fig, ax = plt.subplots(1, 3, figsize=(11, 3.3))\n    for a, col, lab in zip(ax, [\"haa\", \"kurt\", \"vol\"], [\"% high attenuation (fibrosis proxy)\", \"HU kurtosis\", \"lung volume (mL)\"]):\n        a.scatter(per_pat[col], per_pat.slope, s=12, alpha=.7); a.set_xlabel(lab); a.set_ylabel(\"FVC slope (mL/week)\")\n        r = np.corrcoef(per_pat[col], per_pat.slope)[0, 1]; a.set_title(f\"r = {r:.2f}\")\n    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, \"ct_feature_boxplots.png\"), dpi=150); plt.close()\n    print(f\"\\nSaved ablation outputs to {OUT_DIR}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-07T20:46:21.633216Z","iopub.execute_input":"2026-09-07T20:46:21.633919Z","iopub.status.idle":"2026-09-07T20:47:14.736023Z","shell.execute_reply.started":"2026-09-07T20:46:21.63388Z","shell.execute_reply":"2026-09-07T20:47:14.734978Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}