{"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":"from pathlib import Path\nimport pandas as pd\nimport pydicom\nimport os\nimport numpy as np\nimport sys\nfrom datetime import datetime\nfrom scipy.ndimage import zoom, sobel\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.model_selection import KFold\nfrom xgboost import XGBClassifier\n\n# ==============================================================================\n# PIPELINE OVERVIEW & DOCUMENTATION\n# ==============================================================================\nprint(\"\"\"\n================================================================================\nRSNA KNEE ABNORMALITY DETECTION: TIMED EXTRACTION & MODEL OPTIMIZATION\nCompetition Link: https://www.kaggle.com/competitions/rsna-knee-abnormality-detection\n\nPIPELINE WORKFLOW & RELEVANT STEPS:\n1. Load Training Labels Data (LLM Report Labels) & Series Metadata\n2. Inspect Label Distribution, Binarize Probabilistic Scores (`y > 0.5`), \n   and Report Post-Conversion Patient Statistics (All 0s vs All 1s)\n3. Execute & Time High-Dimensional Feature Extraction (~5,120 Features)\n4. Initialize 10-Fold Cross-Validation Splitter \n5. Train XGBoost with Tuned Regularization & Subsampling\n6. Evaluate via Macro-Averaged ROC AUC on Test Folds \n7. Output Timing Performance & Top 15 Most Important Features\n================================================================================\n\"\"\")\n\n# ==========================================\n# Configuration Constants & File Paths\n# ==========================================\nTRY_TO_LOAD_DATA_FROM_DISK = False\nFEATURES_FILE_PATH = Path(\"extracted_features_5k.csv\")\n\nCOMP = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")\nLLM_LABELS_PATH = Path(\"/kaggle/input/datasets/pilkwang/rsna-knee-llm-labels/report_labels_v2.csv\")\nTRAIN_SERIES_PATH = COMP / \"train_series.csv\"\n\nLABELS = [\n    \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\", \"Medial OA\",\n    \"Lateral OA\", \"PF OA\", \"Effusion\", \"Synovitis\", \"Baker's\",\n    \"Contusion\", \"Fracture\",\n]\n\nID_COL = \"StudyInstanceUID\"\n\n\n# ==========================================\n# 1. Data Loading & Validation Functions\n# ==========================================\n\ndef load_datasets():\n    train_df = pd.read_csv(LLM_LABELS_PATH).drop(\n        columns=[c for c in pd.read_csv(LLM_LABELS_PATH, nrows=1).columns if \"__conf\" in c or \"__verdict\" in c]\n    )\n    train_series_df = pd.read_csv(TRAIN_SERIES_PATH)\n    return train_df, train_series_df\n\n\ndef validate_datasets(train_df, train_series_df):\n    missing_studies = set(train_df[ID_COL]) - set(train_series_df[ID_COL])\n    print(f\"1. All train_df StudyInstanceUID accounted for: {len(missing_studies) == 0}\")\n    print(f\"2. Train shape: {train_df.shape} | Train Series shape: {train_series_df.shape}\\n\")\n\n\ndef inspect_and_summarize_labels(train_df):\n    print(\"The raw LLM report labels contain probabilistic/continuous scores ranging from\")\n    print(\"0.08 to 0.94 representing model certainty or soft likelihoods. However, ranking\")\n    print(\"metrics like ROC AUC and standard classifiers require strict binary indicators\")\n    print(\"for training and evaluation. By applying a threshold of > 0.5, we convert these\")\n    print(\"soft scores into clear binary targets (0 = normal/unlikely, 1 = abnormal/likely).\\n\")\n    \n    print(\"================ LABELS DISTRIBUTION & NaN INSPECTION ================\")\n    summary_data = []\n    for label in LABELS:\n        nan_count = train_df[label].isna().sum()\n        value_counts = train_df[label].value_counts().to_dict()\n        summary_data.append({\n            \"Label\": label,\n            \"NaN_Count\": nan_count,\n            \"Value_Counts\": value_counts\n        })\n    \n    summary_df = pd.DataFrame(summary_data)\n    for _, row in summary_df.iterrows():\n        print(f\"{row['Label']:>16} {row['NaN_Count']:10} {row['Value_Counts']}\")\n    print(\"=====================================================================\\n\")\n\n\ndef inspect_binary_labels_and_patient_stats(train_df):\n    print(\"================ POST-CONVERSION (0/1) LABEL & PATIENT STATS ================\")\n    \n    # Create a binary DataFrame for the labels\n    df_binary = train_df[LABELS].applymap(lambda x: 1 if pd.notnull(x) and x > 0.5 else 0)\n    \n    summary_binary_data = []\n    for label in LABELS:\n        val_counts = df_binary[label].value_counts().to_dict()\n        summary_binary_data.append({\n            \"Label\": label,\n            \"Binary_Counts (0 vs 1)\": val_counts\n        })\n        \n    for item in summary_binary_data:\n        print(f\"{item['Label']:>16} -> {item['Binary_Counts (0 vs 1)']}\")\n        \n    # Calculate patients with all 0s vs all 1s across the 12 clinical targets\n    row_sums = df_binary.sum(axis=1)\n    all_zeros_count = (row_sums == 0).sum()\n    all_ones_count = (row_sums == len(LABELS)).sum()\n    \n    print(\"---------------------------------------------------------------------\")\n    print(f\"Total Patients (Studies) Evaluated: {len(df_binary)}\")\n    print(f\"Patients with ALL 0s (completely normal across all 12 targets): {all_zeros_count}\")\n    print(f\"Patients with ALL 1s (positive across all 12 targets): {all_ones_count}\")\n    print(\"=====================================================================\\n\")\n\n\n# ==========================================\n# 2. 5K Feature Extraction with Timing\n# ==========================================\n\nBASE_COLS = [\n    'img_mean', 'img_std', 'img_median', 'img_iqr', 'img_cv', \n    'img_p1', 'img_p5', 'img_p10', 'img_p25', 'img_p75', 'img_p90', \n    'img_p95', 'img_p99', 'img_max', 'img_min', 'edge_mean', 'edge_std', \n    'edge_max', 'sobel_mean', 'sobel_max'\n]\nfor r in range(10):\n    for c in range(10):\n        BASE_COLS.extend([f'z_{r}_{c}_mean', f'z_{r}_{c}_std', f'z_{r}_{c}_max'])\n\n\ndef extract_series_pixel_features(series_dir):\n    if not os.path.exists(series_dir):\n        return None\n        \n    dcm_files = [os.path.join(series_dir, f) for f in os.listdir(series_dir) if f.lower().endswith('.dcm')]\n    if not dcm_files:\n        return None\n    \n    dcm_files.sort()\n    step = max(1, len(dcm_files) // 16)\n    sampled_files = dcm_files[::step][:16]\n    \n    slice_features = []\n    \n    for f_path in sampled_files:\n        try:\n            ds = pydicom.dcmread(f_path)\n            if not hasattr(ds, 'pixel_array'):\n                continue\n            img = ds.pixel_array.astype(np.float32)\n            \n            if img.max() > img.min():\n                img = (img - img.min()) / (img.max() - img.min())\n            \n            if img.shape[0] != 120 or img.shape[1] != 120:\n                zoom_factors = (120 / img.shape[0], 120 / img.shape[1])\n                img = zoom(img, zoom_factors, order=1)\n\n            feat_dict = {\n                'img_mean': np.mean(img), 'img_std': np.std(img),\n                'img_median': np.median(img),\n                'img_iqr': np.percentile(img, 75) - np.percentile(img, 25),\n                'img_cv': np.std(img) / (np.mean(img) + 1e-5),\n                'img_p1': np.percentile(img, 1), 'img_p5': np.percentile(img, 5),\n                'img_p10': np.percentile(img, 10), 'img_p25': np.percentile(img, 25),\n                'img_p75': np.percentile(img, 75), 'img_p90': np.percentile(img, 90),\n                'img_p95': np.percentile(img, 95), 'img_p99': np.percentile(img, 99),\n                'img_max': np.max(img), 'img_min': np.min(img)\n            }\n            \n            gy, gx = np.gradient(img)\n            grad_mag = np.sqrt(gx**2 + gy**2)\n            sobel_mag = np.sqrt(sobel(img, axis=1)**2 + sobel(img, axis=0)**2)\n            \n            feat_dict.update({\n                'edge_mean': np.mean(grad_mag), 'edge_std': np.std(grad_mag), 'edge_max': np.max(grad_mag),\n                'sobel_mean': np.mean(sobel_mag), 'sobel_max': np.max(sobel_mag)\n            })\n            \n            blocks = img.reshape(10, 12, 10, 12).transpose(0, 2, 1, 3)\n            block_means = blocks.mean(axis=(2, 3)).flatten()\n            block_stds = blocks.std(axis=(2, 3)).flatten()\n            block_maxs = blocks.max(axis=(2, 3)).flatten()\n            \n            for i in range(100):\n                r, c = divmod(i, 10)\n                feat_dict[f'z_{r}_{c}_mean'] = block_means[i]\n                feat_dict[f'z_{r}_{c}_std'] = block_stds[i]\n                feat_dict[f'z_{r}_{c}_max'] = block_maxs[i]\n            \n            slice_features.append(feat_dict)\n        except Exception:\n            continue\n            \n    if not slice_features:\n        return None\n        \n    df_slice = pd.DataFrame(slice_features)\n    series_agg = {}\n    for col in BASE_COLS:\n        if col in df_slice.columns:\n            series_agg[f'{col}_mean'] = df_slice[col].mean()\n            series_agg[f'{col}_max'] = df_slice[col].max()\n            series_agg[f'{col}_std'] = df_slice[col].std()\n            series_agg[f'{col}_median'] = df_slice[col].median()\n        \n    return series_agg\n\n\ndef get_or_create_features(train_df, train_series_df):\n    if TRY_TO_LOAD_DATA_FROM_DISK and FEATURES_FILE_PATH.exists():\n        df_train_features = pd.read_csv(FEATURES_FILE_PATH).set_index(ID_COL)\n        print(f\"Loaded 5K features from disk. Shape: {df_train_features.shape}\")\n        return df_train_features, 0.0\n\n    print(f\"\\nExtracting ~5,120 features via Vectorized 10x10 spatial grids...\")\n    extraction_start = datetime.now()\n    \n    study_ids = train_df[ID_COL].unique()\n    total_cases = len(study_ids)\n    feature_records = []\n    \n    for i, study_id in enumerate(study_ids, 1):\n        if i % 50 == 0 or i == total_cases:\n            elapsed_so_far = (datetime.now() - extraction_start).total_seconds()\n            rate = i / max(1, elapsed_so_far)\n            eta = (total_cases - i) / max(0.1, rate)\n            sys.stdout.write(f\"\\rProcessing case {i}/{total_cases} | Elapsed: {elapsed_so_far:.1f}s | ETA: {eta:.1f}s\")\n            sys.stdout.flush()\n        \n        case_series = train_series_df[train_series_df[ID_COL] == study_id]\n        series_features_list = []\n        \n        for _, row in case_series.iterrows():\n            series_id = row.get(\"SeriesInstanceUID\", None)\n            if series_id:\n                series_dir = COMP / \"train\" / str(study_id) / str(series_id)\n                feat = extract_series_pixel_features(str(series_dir))\n                if feat is not None:\n                    series_features_list.append(feat)\n                        \n        study_feat = {ID_COL: study_id}\n        expected_series_cols = [f\"{c}_{agg}\" for c in BASE_COLS for agg in ['mean', 'max', 'std', 'median']]\n        \n        if series_features_list:\n            df_s = pd.DataFrame(series_features_list)\n            for col in expected_series_cols:\n                if col in df_s.columns:\n                    vals = df_s[col].astype(float)\n                    study_feat[f'{col}_study_mean'] = vals.mean()\n                    study_feat[f'{col}_study_max'] = vals.max()\n                    study_feat[f'{col}_study_std'] = vals.std()\n                    study_feat[f'{col}_study_median'] = vals.median()\n                else:\n                    study_feat[f'{col}_study_mean'] = 0.0\n                    study_feat[f'{col}_study_max'] = 0.0\n                    study_feat[f'{col}_study_std'] = 0.0\n                    study_feat[f'{col}_study_median'] = 0.0\n        else:\n            for col in expected_series_cols:\n                study_feat[f'{col}_study_mean'] = 0.0\n                study_feat[f'{col}_study_max'] = 0.0\n                study_feat[f'{col}_study_std'] = 0.0\n                study_feat[f'{col}_study_median'] = 0.0\n                \n        feature_records.append(study_feat)\n                \n    extraction_end = datetime.now()\n    extraction_duration = (extraction_end - extraction_start).total_seconds()\n    print(f\"\\nExtraction complete in {extraction_duration:.2f} seconds ({extraction_duration / 60:.2f} minutes).\")\n    \n    df_train_features = pd.DataFrame(feature_records).set_index(ID_COL).fillna(0)\n    df_train_features.to_csv(FEATURES_FILE_PATH)\n    return df_train_features, extraction_duration\n\n\n# ==========================================\n# 3. XGBoost Full-Feature Evaluation\n# ==========================================\n\ndef evaluate_xgboost_all_features(df_train_features, df_train_labels):\n    print(f\"\\n================ MODEL TRAINING (ALL {df_train_features.shape[1]} FEATURES) ================\")\n    kf = KFold(n_splits=10, shuffle=True, random_state=42)\n    X = df_train_features.values\n    feature_names = df_train_features.columns\n    global_importances = np.zeros(len(feature_names))\n    \n    mean_auc_list = []\n    \n    for label in LABELS:\n        valid_idx = df_train_labels[label].dropna().index\n        y = (df_train_labels.loc[valid_idx, label].values > 0.5).astype(int)\n        \n        if len(np.unique(y)) < 2:\n            continue\n            \n        sub_indices = [df_train_features.index.get_loc(idx_val) for idx_val in valid_idx]\n        X_label = X[sub_indices]\n        \n        fold_preds = np.zeros(len(y))\n        fold_trues = np.zeros(len(y))\n        label_importances = np.zeros(len(feature_names))\n        \n        for train_idx, test_idx in kf.split(X_label):\n            X_train, X_test = X_label[train_idx], X_label[test_idx]\n            y_train, y_test = y[train_idx], y[test_idx]\n            \n            clf = XGBClassifier(\n                n_estimators=150,\n                learning_rate=0.03,\n                max_depth=4,\n                colsample_bytree=0.15,\n                subsample=0.8,\n                reg_alpha=1.0,\n                reg_lambda=1.0,\n                random_state=42,\n                eval_metric='logloss',\n                n_jobs=-1\n            )\n            clf.fit(X_train, y_train)\n            \n            probas = clf.predict_proba(X_test)[:, 1]\n            fold_preds[test_idx] = probas\n            fold_trues[test_idx] = y_test\n            \n            label_importances += clf.feature_importances_ / 10.0\n            \n        global_importances += label_importances / len(LABELS)\n        \n        auc = roc_auc_score(fold_trues, fold_preds)\n        mean_auc_list.append(auc)\n        print(f\"[{label.rjust(16)}] 10-Fold AUC: {auc:.4f}\")\n        \n    overall_mean_auc = np.mean(mean_auc_list)\n    print(f\"--------------------------------------------------\")\n    print(f\"Overall Macro-Averaged AUC: {overall_mean_auc:.4f}\")\n    \n    top_idx = np.argsort(global_importances)[::-1][:15]\n    print(\"\\nTop 15 Most Predictive Features:\")\n    print(f\"{'Rank':<4} {'Feature Name':<45} {'Importance Score'}\")\n    print(\"-\" * 65)\n    for rank, i in enumerate(top_idx, 1):\n        print(f\"{rank:<4} {str(feature_names[i]):<45} {global_importances[i]:.5f}\")\n        \n    return overall_mean_auc\n\n\n# ==========================================\n# Main Execution Entry Point\n# ==========================================\nif __name__ == \"__main__\":\n    script_start_time = datetime.now()\n    train_df, train_series_df = load_datasets()\n    validate_datasets(train_df, train_series_df)\n    \n    # 1. Inspect initial label distribution and probabilistic score thresholding context\n    inspect_and_summarize_labels(train_df)\n    \n    # 2. Inspect post-conversion 0/1 label distribution and full-target patient sums (all 0s / all 1s)\n    inspect_binary_labels_and_patient_stats(train_df)\n\n    # 3. Feature Extraction with elapsed timing measurement\n    df_train_features, extraction_duration = get_or_create_features(train_df, train_series_df)\n\n    df_train_labels = train_df.set_index(ID_COL)[LABELS]\n    df_train_labels = df_train_labels.loc[df_train_features.index]\n\n    # 4. XGBoost Evaluation\n    overall_mean_auc = evaluate_xgboost_all_features(df_train_features, df_train_labels)\n\n    script_end_time = datetime.now()\n    total_elapsed = (script_end_time - script_start_time).total_seconds()\n\n    print(\"\\n================ FINAL SUMMARY ================\")\n    print(f\"1. Feature Extraction Time: {extraction_duration:.2f} seconds ({extraction_duration / 60:.2f} minutes)\")\n    print(f\"2. Overall Mean AUC: {overall_mean_auc:.4f}\")\n    print(f\"3. Total Pipeline Time: {total_elapsed:.2f} seconds ({total_elapsed / 60:.2f} minutes)\")\n    print(\"=============================================\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-08-15T21:27:25.515607Z","iopub.execute_input":"2026-08-15T21:27:25.516128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}