diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..a05cee0 --- /dev/null +++ b/.dockerignore @@ -0,0 +1,68 @@ +# Git +.git/ +.gitignore + +# Python cache +__pycache__/ +*.py[cod] +*.pyo +*.pyd +*.so +*.egg +*.egg-info/ +dist/ +build/ + +# Virtual environments +.venv/ +venv/ +env/ +ENV/ + +# Jupyter +.ipynb_checkpoints/ +*.ipynb + +# ML artifacts +*.pkl +*.pickle +*.joblib +*.h5 +*.hdf5 +*.pt +*.pth +*.ckpt +*.safetensors +models/ +checkpoints/ +runs/ + +# Data (large files not needed at build time) +*.csv.gz +*.parquet +*.feather +data/ +datasets/ +raw/ + +# Outputs +output/ +outputs/ +results/ +logs/ +*.log + +# IDE / OS +.idea/ +.vscode/ +.DS_Store +*.swp +*.swo + +# Environment variables +.env +.env.* + +# Docs +README.md +LICENSE diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fa1523c --- /dev/null +++ b/.gitignore @@ -0,0 +1,103 @@ +# Python +__pycache__/ +*.py[cod] +*.pyo +*.pyd +*.so +*.egg +*.egg-info/ +dist/ +build/ +.eggs/ +wheels/ + +# Virtual environments +.venv/ +venv/ +env/ +ENV/ + +# Jupyter +.ipynb_checkpoints/ +*.ipynb + +# Local dev working dirs (build_dev_subset.py outputs) +dev_subset/ +dev_model/ +dev_outputs/ + +# ML artifacts +*.pkl +*.pickle +*.joblib +*.h5 +*.hdf5 +*.pt +*.pth +*.ckpt +*.safetensors +models/ +checkpoints/ +runs/ + +# Data +*.csv.gz +*.parquet +*.feather +data/ +datasets/ +raw/ + +# Outputs +output/ +outputs/ +results/ +logs/ +*.log + + +# macOS +.DS_Store +.AppleDouble +.LSOverride + +# IDE +.idea/ +.vscode/ +*.swp +*.swo + +# Environment variables +.env +.env.* + +# pytest +.pytest_cache/ +.coverage +htmlcov/ + + +#measurment artifacts +*.edf +*.sav +full_model/ +full_outputs/ +docker_model/ +dev_model/ +dev_outputs/ +dev_model_v2/ +dev_outputs_v2/ +features/FEATURES.md +CLAUDE.md +build_dev_subset.py +learning_log.md +learning_log_p2.md +eda_outputs/ +docker_outputs/ + +#files +phase1_eda.py +/loso_cv.py +archive_outputs/ +eda/ +/tools \ No newline at end of file diff --git a/extract.py b/extract.py new file mode 100644 index 0000000..55b2798 --- /dev/null +++ b/extract.py @@ -0,0 +1,51 @@ +import sys +import numpy as np + +# Adjust these paths to match your repo layout +sys.path.insert(0, 'tools') # so `from loso_cv import extract_all_features` resolves +sys.path.insert(0, '.') # repo root, for helper_code/features imports loso_cv.py needs + +from loso_cv import extract_all_features + +DATA_PATH = '/Users/briandillon/.cache/kagglehub/datasets/physionet/physionetchallenge2026data/versions/7' + +# Ground truth: loso_cv.py's own extraction (same function that produced every logged result) +X_truth, y_truth, ages_truth, sites_truth, ids_truth = extract_all_features(DATA_PATH, verbose=True) + +# Cache: build_features_cache.py's output +d = np.load('tools/small_features.npz', allow_pickle=True) +X_cache, y_cache, ids_cache = d['X'], d['y'], d['patient_id'] + +print(f"Ground truth: {X_truth.shape}, cache: {X_cache.shape}") + +# Align by patient_id — order may differ between the two extraction runs +truth_lookup = {pid: i for i, pid in enumerate(ids_truth)} +cache_lookup = {pid: i for i, pid in enumerate(ids_cache)} + +common_ids = sorted(set(truth_lookup) & set(cache_lookup)) +print(f"Common patient IDs: {len(common_ids)} (expect 1103)") + +max_diffs = [] +mismatched_patients = [] + +for pid in common_ids: + row_truth = X_truth[truth_lookup[pid]] + row_cache = X_cache[cache_lookup[pid]] + # NaN-safe comparison — both matrices use NaN for missing CAISR data + diff = np.abs(np.nan_to_num(row_truth, nan=-999) - np.nan_to_num(row_cache, nan=-999)) + max_diff = diff.max() + max_diffs.append(max_diff) + if max_diff > 1e-4: + mismatched_patients.append((pid, max_diff, np.argmax(diff))) + +max_diffs = np.array(max_diffs) +print(f"\nMax abs diff across all patients/features: {max_diffs.max():.6f}") +print(f"Patients with any mismatch (>1e-4): {len(mismatched_patients)}") + +if mismatched_patients: + print("\nFirst 5 mismatches (patient_id, max_diff, feature_column_index):") + for m in mismatched_patients[:5]: + print(" ", m) +else: + print("\nCLEAN — every patient's 48-length vector matches exactly between " + "loso_cv.py's extraction and build_features_cache.py's cache.") \ No newline at end of file diff --git a/features/__init__.py b/features/__init__.py new file mode 100644 index 0000000..9a2ca78 --- /dev/null +++ b/features/__init__.py @@ -0,0 +1,87 @@ +# features/__init__.py +# Shared constants for the Narnia_ML feature pipeline. +# Update these whenever a feature module adds or removes features. +# +# Entry 2 (2026-07-01): Dropped absolute Hjorth (site-confounded per LOSO). +# Added within-recording ratio features (site-stable by construction). +# Feature count: 82 → 48. +# +# Entry 4 (2026-07-02): Added 2 age-residualized features (CA_rate, +# EEG_var_REM_Wake) via features/age_residuals.py — see that module's +# header for the age-residualized EDA rationale. These are appended by an +# AgeResidualizer pipeline step AFTER extraction/hstack, not by a feature +# module called from team_code.py's extraction block, so they live at the +# end of the vector rather than inside the demo/base/enriched/ratio blocks. +# Feature count: 48 → 50. +# +# NEW candidate, not yet validated (2026-07-20): Added 3 stage-conditional +# limb movement features (Limb_REM, Limb_NREM, Limb_REM_NREM_ratio) via +# features/caisr_enriched.py — same enrichment pattern already proven for +# AHI (REM_AHI/NREM_AHI/REM_NREM_AHI_ratio), applied to Limb_idx instead. +# Appended at the END of the enriched block (after N3_entropy, before the +# ratio block starts) — inserted, not appended to the whole vector, so +# downstream index constants that are defined relative to IDX_RATIO_START +# (e.g. IDX_EEG_VAR_REM_WAKE below) auto-update correctly; nothing already +# shipped needs a manual index bump. Feature count: 50 → 53 (raw 48 → 51, +# plus the same 2 age-residual features as before). Requires regenerating +# any --features-cache built before this change — the raw hstack length +# changed from 48 to 51. + + +N_DEMOGRAPHIC_FEATURES = 10 # features/demographic.py +N_CAISR_BASE_FEATURES = 12 # features/caisr_base.py +N_CAISR_ENRICHED_FEATURES = 14 # features/caisr_enriched.py (was 11 — +3 limb features, 2026-07-20) +N_RATIO_FEATURES = 15 # features/physiological_ratios.py +N_AGE_RESIDUAL_FEATURES = 2 # features/age_residuals.py (Entry 4, pipeline-appended) + +N_PHYSIOLOGICAL_FEATURES = 0 +N_ALGORITHMIC_FEATURES = N_CAISR_BASE_FEATURES + N_CAISR_ENRICHED_FEATURES # 23 + +N_TOTAL_FEATURES = ( + N_DEMOGRAPHIC_FEATURES + + N_CAISR_BASE_FEATURES + + N_CAISR_ENRICHED_FEATURES + + N_RATIO_FEATURES + + N_AGE_RESIDUAL_FEATURES +) # 50 + +IDX_DEMO_START = 0 +IDX_DEMO_END = N_DEMOGRAPHIC_FEATURES # 10 +IDX_BASE_START = IDX_DEMO_END +IDX_BASE_END = IDX_BASE_START + N_CAISR_BASE_FEATURES # 22 +IDX_ENRICHED_START = IDX_BASE_END +IDX_ENRICHED_END = IDX_ENRICHED_START + N_CAISR_ENRICHED_FEATURES # 33 +IDX_RATIO_START = IDX_ENRICHED_END +IDX_RATIO_END = IDX_RATIO_START + N_RATIO_FEATURES # 48 +IDX_AGE_RESID_START = IDX_RATIO_END +IDX_AGE_RESID_END = IDX_AGE_RESID_START + N_AGE_RESIDUAL_FEATURES # 50 + +IDX_AGE = IDX_DEMO_START # Age is demo feature 0 +IDX_CA_RATE = IDX_ENRICHED_START + 1 # CA_rate is caisr_enriched feature 1 +IDX_EEG_VAR_REM_WAKE = IDX_RATIO_START + 11 # EEG_var_REM_Wake is ratio feature 11 + +# Human-readable names for the 48-length PRE-AgeResidualizer vector, index-aligned. +# Pulled directly from FEATURES.md's Entry 2 (48-feature) index table — real +# names, not placeholders. Used by reg_sweep.py to label coefficients / top +# features so sweep output is directly readable against your own registry. +FEATURE_NAMES_48 = [ + 'Age', 'Sex_F', 'Sex_M', 'Sex_Unk', 'Race_Asian', 'Race_Black', + 'Race_Other', 'Race_Unavailable', 'Race_White', 'BMI', + 'AHI_total', 'Arousal_idx', 'Limb_idx', 'Wake_pct', 'N1_pct', 'N2_pct', + 'N3_pct', 'REM_pct', 'Sleep_eff', 'Prob_W', 'Prob_N3', 'Prob_arous', + 'OA_rate', 'CA_rate', 'HY_rate', 'RERA_rate', 'CA_total_ratio', + 'REM_AHI', 'NREM_AHI', 'REM_NREM_AHI_ratio', 'N3_gradient', + 'Spont_arousal_idx', 'N3_entropy', + 'Limb_REM', 'Limb_NREM', 'Limb_REM_NREM_ratio', + 'EEG_std_N3_Wake', 'EEG_mav_N3_Wake', 'EEG_zcr_N3_Wake', 'EEG_rms_N3_Wake', + 'EEG_var_N3_Wake', 'EEG_mob_N3_Wake', 'EEG_cplx_N3_Wake', + 'EEG_std_REM_Wake', 'EEG_mav_REM_Wake', 'EEG_zcr_REM_Wake', 'EEG_rms_REM_Wake', + 'EEG_var_REM_Wake', 'EEG_mob_REM_Wake', 'EEG_cplx_REM_Wake', + 'Chin_mav_REM_NREM_atonia', +] +assert len(FEATURE_NAMES_48) == IDX_RATIO_END, \ + f"FEATURE_NAMES_48 length {len(FEATURE_NAMES_48)} != {IDX_RATIO_END}" +assert FEATURE_NAMES_48[IDX_CA_RATE] == 'CA_rate' +assert FEATURE_NAMES_48[IDX_EEG_VAR_REM_WAKE] == 'EEG_var_REM_Wake' + +FEATURE_NAMES_50 = FEATURE_NAMES_48 + ['CA_rate_age_residual', 'EEG_var_REM_Wake_age_residual'] \ No newline at end of file diff --git a/features/age_residuals.py b/features/age_residuals.py new file mode 100644 index 0000000..cd184c4 --- /dev/null +++ b/features/age_residuals.py @@ -0,0 +1,71 @@ +# features/age_residuals.py +# Age-residualized features — Entry 4 +# +# Origin: 2026-07-02 age-residualized EDA (see learning_log.md). CA_rate and +# EEG_var_REM_Wake were previously written off as uninformative in the raw +# Phase 1 EDA (AUROC ~0.52, BH-FDR p_adj > 0.7 both) — but regressing Age out +# of every feature and re-testing the residual showed both carry real signal +# that raw testing couldn't see: +# +# Feature Raw AUROC / p_adj Residual AUROC / p_adj +# CA_rate 0.520 / 0.805 0.622 / 0.0075 +# EEG_var_REM_Wake 0.520 / 0.808 0.618 / 0.0098 +# +# Age is available at real inference time (unlike the ICD-subtype features +# ruled out earlier for exactly that reason) — this is feature engineering, +# not leakage. +# +# Fitting discipline: identical to SimpleImputer. The linear age-regression +# coefficient for each source feature is fit ONCE on training data (inside +# AgeResidualizer.fit(), called by Pipeline.fit()) and frozen. transform() +# only ever applies the already-fitted coefficients — it must never re-fit, +# including at inference time on a single patient. +# +# Source features (already present earlier in the vector — see +# features/__init__.py for the index constants used to locate them): +# Age index 0 (features/demographic.py) +# CA_rate IDX_ENRICHED_START + 1 (features/caisr_enriched.py) +# EEG_var_REM_Wake IDX_RATIO_START + 11 (features/physiological_ratios.py) +# +# These 2 new features are APPENDED at the end of the vector (indices 48-49) +# rather than inserted mid-vector, so the existing demo/base/enriched/ratio +# index constants don't shift. Total feature count: 48 -> 50. + +import numpy as np +from sklearn.base import BaseEstimator, TransformerMixin + +N_AGE_RESIDUAL_FEATURES = 2 +_MIN_FIT_ROWS = 10 + + +class AgeResidualizer(BaseEstimator, TransformerMixin): + def __init__(self, age_idx, ca_rate_idx, eeg_var_rem_wake_idx): + self.age_idx = age_idx + self.ca_rate_idx = ca_rate_idx + self.eeg_var_rem_wake_idx = eeg_var_rem_wake_idx + + @staticmethod + def _fit_one(age, feature): + mask = ~(np.isnan(age) | np.isnan(feature)) + if int(mask.sum()) < _MIN_FIT_ROWS: + return 0.0, 0.0 + slope, intercept = np.polyfit(age[mask], feature[mask], 1) + return float(slope), float(intercept) + + def fit(self, X, y=None): + X = np.asarray(X, dtype=float) + age = X[:, self.age_idx] + self.ca_rate_slope_, self.ca_rate_intercept_ = self._fit_one( + age, X[:, self.ca_rate_idx]) + self.eeg_slope_, self.eeg_intercept_ = self._fit_one( + age, X[:, self.eeg_var_rem_wake_idx]) + return self + + def transform(self, X): + X = np.asarray(X, dtype=float) + age = X[:, self.age_idx] + ca_rate_resid = X[:, self.ca_rate_idx] - ( + self.ca_rate_intercept_ + self.ca_rate_slope_ * age) + eeg_resid = X[:, self.eeg_var_rem_wake_idx] - ( + self.eeg_intercept_ + self.eeg_slope_ * age) + return np.hstack([X, ca_rate_resid.reshape(-1, 1), eeg_resid.reshape(-1, 1)]) \ No newline at end of file diff --git a/features/caisr_base.py b/features/caisr_base.py new file mode 100644 index 0000000..3338cdd --- /dev/null +++ b/features/caisr_base.py @@ -0,0 +1,78 @@ +# features/caisr_base.py +# Baseline CAISR annotation features — features 59-70 +# +# Iteration history: +# v1 (2026-06-30): Initial — AHI, arousal idx, limb idx, +# stage percentages (W/N1/N2/N3/REM), +# sleep efficiency, mean P(Wake), mean P(N3), +# mean P(arousal) [12 features] + +import numpy as np + + +def extract_caisr_base_features(algo_data): + """ + Extracts baseline CAISR features from annotation channels. + + Returns np.ndarray of length 12: + [0] AHI total (events/hr) + [1] Arousal index (events/hr) + [2] Limb movement index (events/hr) + [3] Wake % + [4] N1 % + [5] N2 % + [6] N3 % + [7] REM % + [8] Sleep efficiency + [9] Mean CAISR P(Wake) + [10] Mean CAISR P(N3) + [11] Mean CAISR P(arousal) + """ + if not algo_data: + return np.full(12, float('nan')) + + resp = algo_data.get('resp_caisr', np.array([])) + arousal = algo_data.get('arousal_caisr', np.array([])) + limb = algo_data.get('limb_caisr', np.array([])) + stages_raw = algo_data.get('stage_caisr', np.array([])) + + total_hours_resp = len(resp) / 3600.0 if len(resp) > 0 else 0.0 + total_hours_arousal = len(arousal) / 7200.0 if len(arousal) > 0 else 0.0 + total_hours_limb = len(limb) / 3600.0 if len(limb) > 0 else 0.0 + + def count_events(sig, t_hours): + if len(sig) == 0 or t_hours <= 0: + return float('nan') + binary = (np.asarray(sig) > 0).astype(int) + edges = np.diff(binary, prepend=0) + return np.count_nonzero(edges == 1) / t_hours + + ahi_auto = count_events(resp, total_hours_resp) + arousal_idx = count_events(arousal, total_hours_arousal) + limb_idx = count_events(limb, total_hours_limb) + + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) > 0: + w_pct = float(np.mean(valid_stages == 5)) + n1_pct = float(np.mean(valid_stages == 3)) + n2_pct = float(np.mean(valid_stages == 2)) + n3_pct = float(np.mean(valid_stages == 1)) + rem_pct = float(np.mean(valid_stages == 4)) + efficiency = float(np.mean((valid_stages >= 1) & (valid_stages <= 4))) + else: + w_pct = n1_pct = n2_pct = n3_pct = rem_pct = efficiency = float('nan') + + prob_w = float(np.mean(algo_data.get('caisr_prob_w', [float('nan')]))) + prob_n3 = float(np.mean(algo_data.get('caisr_prob_n3', [float('nan')]))) + prob_arou = float(np.mean(algo_data.get('caisr_prob_arous', [float('nan')]))) + + prob_w = prob_w if prob_w <= 1.0 else float('nan') + prob_n3 = prob_n3 if prob_n3 <= 1.0 else float('nan') + prob_arou = prob_arou if prob_arou <= 1.0 else float('nan') + + return np.array([ + ahi_auto, arousal_idx, limb_idx, + w_pct, n1_pct, n2_pct, n3_pct, rem_pct, efficiency, + prob_w, prob_n3, prob_arou, + ], dtype=float) \ No newline at end of file diff --git a/features/caisr_enriched.py b/features/caisr_enriched.py new file mode 100644 index 0000000..c7bd02d --- /dev/null +++ b/features/caisr_enriched.py @@ -0,0 +1,833 @@ +# features/caisr_enriched.py +# Enriched CAISR annotation features — features 71-81 (v1), extended to +# 71-84 (v5, 2026-07-20) with stage-conditional limb movement features +# +# Iteration history: +# v1 (2026-06-30): First enrichment pass — respiratory breakdown, +# stage-conditional AHI, N3 gradient, spontaneous +# arousal index, N3 confidence entropy [11 features] +# v5 (2026-07-20): Added Limb_REM, Limb_NREM, Limb_REM_NREM_ratio — +# stage-conditional limb movement index, same pattern +# as v1's REM_AHI/NREM_AHI/REM_NREM_AHI_ratio, applied +# to Limb_idx (this project's most robust coefficient, +# never previously feature-engineered) [14 features] +# +# Planned (not yet implemented): +# - PLM periodicity index (PLM / isolated LM ratio) Week 3+ +# - CAISR no-arousal confidence entropy (caisr_prob_no-ar) Week 3+ +# - Stage transition rate Week 3+ +# - REM latency Week 3+ + +import numpy as np + + +def extract_caisr_enriched_features(algo_data): + """ + Enriched CAISR features requiring temporal cross-referencing + between annotation channels. + + Returns np.ndarray of length 14: + [0] OA rate (events/hr) — obstructive apnea + [1] CA rate (events/hr) — central apnea + [2] HY rate (events/hr) — hypopnea + [3] RERA rate (events/hr) — effort-related arousals + [4] CA / total AHI ratio — neurological apnea fraction + [5] REM-AHI (events/hr) — apneas during REM + [6] NREM-AHI (events/hr) — apneas during NREM + [7] REM-AHI / NREM-AHI ratio — REM-predominant OSA marker + [8] N3 first-half / second-half ratio — temporal N3 gradient + [9] Spontaneous arousal index — non-respiratory arousals/hr + [10] N3 confidence entropy — slow-wave staging ambiguity + [11] Limb_REM (events/hr) — limb movements during REM + [12] Limb_NREM (events/hr) — limb movements during NREM + [13] Limb_REM / Limb_NREM ratio — see precedent-risk note below + """ + if not algo_data: + return np.full(14, float('nan')) + + resp = algo_data.get('resp_caisr', np.array([])) + arousal = algo_data.get('arousal_caisr', np.array([])) + limb = algo_data.get('limb_caisr', np.array([])) + stages_raw = algo_data.get('stage_caisr', np.array([])) + + total_hours_resp = len(resp) / 3600.0 if len(resp) > 0 else 0.0 + valid_stages = (stages_raw[stages_raw < 9.0] + if len(stages_raw) > 0 else np.array([])) + + # ── Event-type counting helper ──────────────────────────────────────────── + def count_event_type(resp_sig, event_code, t_hours): + if len(resp_sig) == 0 or t_hours <= 0: + return float('nan') + binary = (np.asarray(resp_sig) == event_code).astype(int) + edges = np.diff(binary, prepend=0) + return np.count_nonzero(edges == 1) / t_hours + + def count_events(sig, t_hours): + if len(sig) == 0 or t_hours <= 0: + return float('nan') + binary = (np.asarray(sig) > 0).astype(int) + edges = np.diff(binary, prepend=0) + return np.count_nonzero(edges == 1) / t_hours + + # [0-3] Respiratory event breakdown + # resp_caisr values: 0=none 1=OA 2=CA 3=MA 4=HY 5=RERA + oa_rate = count_event_type(resp, 1, total_hours_resp) + ca_rate = count_event_type(resp, 2, total_hours_resp) + hy_rate = count_event_type(resp, 4, total_hours_resp) + rera_rate = count_event_type(resp, 5, total_hours_resp) + + # [4] CA / total AHI ratio + ahi_total = count_events(resp, total_hours_resp) + ca_total_ratio = (ca_rate / ahi_total + if (not np.isnan(ca_rate) and not np.isnan(ahi_total) + and ahi_total > 0) + else float('nan')) + + # [5-7] Stage-conditional AHI + # Requires upsampling stage_caisr (30s epochs → 1s) before cross-referencing + rem_ahi = nrem_ahi = rem_nrem_ratio = float('nan') + if len(valid_stages) > 0 and len(resp) > 0: + stage_1s = np.repeat(valid_stages, 30)[:len(resp)] + rem_mask = (stage_1s == 4) + nrem_mask = (stage_1s >= 1) & (stage_1s <= 3) + + rem_hours = float(rem_mask.sum()) / 3600.0 + nrem_hours = float(nrem_mask.sum()) / 3600.0 + + if rem_hours > 0: + resp_rem = np.where(rem_mask, resp[:len(rem_mask)], 0) + edges_rem = np.diff((resp_rem > 0).astype(int), prepend=0) + rem_ahi = np.count_nonzero(edges_rem == 1) / rem_hours + + if nrem_hours > 0: + resp_nrem = np.where(nrem_mask, resp[:len(nrem_mask)], 0) + edges_nrem = np.diff((resp_nrem > 0).astype(int), prepend=0) + nrem_ahi = np.count_nonzero(edges_nrem == 1) / nrem_hours + + if (not np.isnan(rem_ahi) and not np.isnan(nrem_ahi) + and nrem_ahi > 0): + rem_nrem_ratio = rem_ahi / nrem_ahi + + # [8] N3 temporal gradient + # Healthy sleep front-loads N3 (ratio >> 1). + # CI-risk pattern: ratio compressed toward 1. + # + # Bug fixed (2026-07-01): original code used 1e-6 as denominator guard, + # producing values of 20,000-150,000 when n3_second == 0 (common in + # elderly / CI+ patients who have no N3 in the second half of the night). + # Fix: return NaN when n3_second < 1% (degenerate case — not meaningful + # as a ratio), and cap at 10.0 for remaining edge cases. + if len(valid_stages) > 1: + mid = len(valid_stages) // 2 + n3_first = float(np.mean(valid_stages[:mid] == 1)) + n3_second = float(np.mean(valid_stages[mid:] == 1)) + if n3_second < 0.01: + # No meaningful N3 in second half — gradient is not defined as ratio. + # Use a signed indicator instead: positive = any first-half N3 exists. + n3_gradient = float('nan') + else: + n3_gradient = min(n3_first / n3_second, 10.0) # cap at 10 + else: + n3_gradient = float('nan') + + # [9] Spontaneous arousal index — REDESIGNED v2 (2026-07-01) + # Original: computed per total recording hour → confounded by REM amount + # (CI+ have less REM → fewer arousals regardless of biology). + # Fix: compute per NREM hour only — removes the REM confound. + # Spontaneous = arousals not coincident with a respiratory event (±15s). + spont_arousal_idx = float('nan') + if len(arousal) > 0 and len(resp) > 0: + n_1s = min(len(arousal) // 2, len(resp)) + if n_1s > 0: + ar_reshaped = arousal[:n_1s * 2].reshape(-1, 2) + arousal_1s = (ar_reshaped.max(axis=1) > 0).astype(int) + resp_1s = (resp[:n_1s] > 0).astype(int) + edges_ar = np.diff(arousal_1s, prepend=0) + ar_starts = np.where(edges_ar == 1)[0] + n_resp_coincident = 0 + for s in ar_starts: + window = resp_1s[max(0, s - 5):min(n_1s, s + 15)] + if len(window) > 0 and window.any(): + n_resp_coincident += 1 + n_spont = max(0, len(ar_starts) - n_resp_coincident) + + # Use NREM hours as denominator (not total recording hours) + # NREM mask from stage_caisr upsampled to 1s resolution + if len(valid_stages) > 0: + stage_1s_ar = np.repeat(valid_stages, 30)[:n_1s] + nrem_mask_ar = (stage_1s_ar >= 1) & (stage_1s_ar <= 3) + nrem_hours_ar = float(nrem_mask_ar.sum()) / 3600.0 + if nrem_hours_ar > 0: + spont_arousal_idx = n_spont / nrem_hours_ar + else: + # Fallback: use total recording hours if no stage data + t_hours_1s = n_1s / 3600.0 + if t_hours_1s > 0: + spont_arousal_idx = n_spont / t_hours_1s + + # [10] N3 confidence entropy + # Binary entropy on caisr_prob_n3: higher = more ambiguous staging. + # Lower entropy = CAISR is confident = cleaner slow waves. + n3_entropy = float('nan') + prob_n3_arr = algo_data.get('caisr_prob_n3', np.array([])) + if len(prob_n3_arr) > 0: + p = np.clip(np.asarray(prob_n3_arr, dtype=float), 1e-9, 1.0 - 1e-9) + h = -(p * np.log(p) + (1.0 - p) * np.log(1.0 - p)) + n3_entropy = float(np.mean(h)) + + # [11-13] Stage-conditional limb movement index — NEW (2026-07-20) + # Direct analog of [5-7]'s REM/NREM-AHI split, applied to limb_caisr + # instead of resp_caisr. Limb_idx (caisr_base.py) is currently a single + # whole-night scalar — the same shape AHI was in v1, before this exact + # enrichment was applied to it. Limb_idx has been the single most + # robust non-demographic coefficient across all three coefficient- + # stability analyses this project has run (50/39/33-feature models), + # yet has never received any feature-engineering attention (all of it + # went to Spont_arousal_idx, which failed on every angle — see + # FEATURES.md's "Closed candidate" section). This is the natural next + # move on a proven-robust feature, using an already-validated + # intervention pattern rather than a new hypothesis. + # + # Precedent risk, stated up front rather than discovered later: + # REM_NREM_AHI_ratio (the exact analogous ratio for AHI) was one of + # the 11 Stage 1 sign-flip drops — only the two ABSOLUTE stage- + # conditional rates (REM_AHI, NREM_AHI) survived. Computing + # Limb_REM_NREM_ratio here anyway (rather than skipping it) so the + # coefficient-stability analysis can make that call on real evidence, + # same "let LOSO/coefficient-stability decide, don't pre-guess" + # discipline as everything else in this file — but if it flips sign + # the way its AHI analog did, that would not be a surprise. + limb_rem = limb_nrem = limb_rem_nrem_ratio = float('nan') + if len(valid_stages) > 0 and len(limb) > 0: + stage_1s_limb = np.repeat(valid_stages, 30)[:len(limb)] + rem_mask_limb = (stage_1s_limb == 4) + nrem_mask_limb = (stage_1s_limb >= 1) & (stage_1s_limb <= 3) + + rem_hours_limb = float(rem_mask_limb.sum()) / 3600.0 + nrem_hours_limb = float(nrem_mask_limb.sum()) / 3600.0 + + if rem_hours_limb > 0: + limb_rem_sig = np.where(rem_mask_limb, limb[:len(rem_mask_limb)], 0) + edges_rem = np.diff((limb_rem_sig > 0).astype(int), prepend=0) + limb_rem = np.count_nonzero(edges_rem == 1) / rem_hours_limb + + if nrem_hours_limb > 0: + limb_nrem_sig = np.where(nrem_mask_limb, limb[:len(nrem_mask_limb)], 0) + edges_nrem = np.diff((limb_nrem_sig > 0).astype(int), prepend=0) + limb_nrem = np.count_nonzero(edges_nrem == 1) / nrem_hours_limb + + if (not np.isnan(limb_rem) and not np.isnan(limb_nrem) + and limb_nrem > 0): + limb_rem_nrem_ratio = limb_rem / limb_nrem + + return np.array([ + oa_rate, ca_rate, hy_rate, rera_rate, + ca_total_ratio, + rem_ahi, nrem_ahi, rem_nrem_ratio, + n3_gradient, + spont_arousal_idx, + n3_entropy, + limb_rem, limb_nrem, limb_rem_nrem_ratio, + ], dtype=float) + + +# ── Wake fragmentation — NOT wired into extract_caisr_enriched_features ── +# Standalone, ablation-testable-first, per the same discipline as +# features/spectral.py's candidates. Hypothesis (learning_log.md, +# 2026-07-12): total wake time (W_pct) and arousal rate (Arousal_idx, +# a brief-event construct) don't distinguish "many short wake bouts" +# from "few long wake bouts" at the same total wake time — this +# fragmentation-vs-duration distinction, especially jointly with REM%, +# is hypothesized to carry CI-risk signal not captured by either +# existing feature. +# +# CAISR-annotation-only — no raw physiological EDF needed for THIS +# feature specifically. (Note: a full pipeline ablation against the +# shipped 48-feature baseline still needs phys_data, since +# physiological_ratios.py's block does — the raw-EDF-free speedup only +# applies to computing this feature itself / a standalone univariate +# signal check, not to a full baseline-vs-with-feature model comparison.) + +N_WAKE_FRAGMENTATION_FEATURES = 2 + +# Minimum sleep-period length to trust the bout count/duration estimate — +# same spirit as physiological_ratios.py's _MIN_STAGE_SECONDS, expressed +# in epochs here since this function works entirely in 30s-epoch space +# (no upsampling to signal rate needed, unlike every spectral.py feature). +_MIN_SLEEP_EPOCHS = 4 # 2 minutes + + +def extract_wake_fragmentation_features(algo_data): + """ + Wake fragmentation within WASO (wake after sleep onset). + + WASO window: from the first sleep-stage epoch (any of N1/N2/N3/REM) + through the last sleep-stage epoch, inclusive. This excludes initial + sleep-onset latency and any trailing wake after the final sleep + epoch — standard WASO convention, chosen specifically to avoid + ambiguity about whether to count a long "I'm done sleeping" wake + tail at the end of the recording. + + Returns np.ndarray of length N_WAKE_FRAGMENTATION_FEATURES (2): + [0] wake_bout_count — number of discrete contiguous wake + episodes within the WASO window. Uses the same rising-edge + counting technique as count_events() above, applied to the + wake-stage binary mask instead of a respiratory/arousal + signal. + [1] mean_wake_bout_duration — total WASO seconds / wake_bout_count. + The feature the underlying hypothesis is actually about: + long total wake time with LOW bout count means few, long + consolidated wake episodes rather than fragmented brief + awakenings, at the same total wake time. + + NaN fallback: + - Missing CAISR annotations → both NaN + - No sleep-stage epochs found at all → both NaN + - WASO window shorter than _MIN_SLEEP_EPOCHS → both NaN + - Zero wake bouts in the WASO window (fully consolidated sleep, + a real and valid outcome) → wake_bout_count=0 (a real value, + not NaN), mean_wake_bout_duration=NaN (0/0 is genuinely + undefined, not zero) + """ + if not algo_data: + return np.full(N_WAKE_FRAGMENTATION_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) == 0: + return np.full(N_WAKE_FRAGMENTATION_FEATURES, float('nan')) + + sleep_mask = (valid_stages >= 1) & (valid_stages <= 4) + if not sleep_mask.any(): + return np.full(N_WAKE_FRAGMENTATION_FEATURES, float('nan')) + + sleep_indices = np.where(sleep_mask)[0] + first_sleep, last_sleep = sleep_indices[0], sleep_indices[-1] + + waso_window = valid_stages[first_sleep:last_sleep + 1] + if len(waso_window) < _MIN_SLEEP_EPOCHS: + return np.full(N_WAKE_FRAGMENTATION_FEATURES, float('nan')) + + wake_binary = (waso_window == 5).astype(int) + edges = np.diff(wake_binary, prepend=0) + wake_bout_count = int(np.count_nonzero(edges == 1)) + + total_waso_epochs = int(wake_binary.sum()) + total_waso_seconds = total_waso_epochs * 30.0 + + if wake_bout_count > 0: + mean_wake_bout_duration = total_waso_seconds / wake_bout_count + else: + mean_wake_bout_duration = float('nan') + + return np.array([wake_bout_count, mean_wake_bout_duration], dtype=float) + + +# ── REM latency — NOT wired into extract_caisr_enriched_features ───────── +# Planned since Week 2-3 (FEATURES.md, "Always" gate — not scale-gated +# like the spectral candidates), never implemented until now. Classic +# clinical sleep-architecture marker. Motivated by the 2026-07-11 MCI/age +# finding: MCI is the earlier, subtler-signal subtype, and REM latency +# shifts are a well-established early marker in the cognitive-decline +# literature — a stronger a priori case than any of the three raw-EDF +# spectral candidates tested 2026-07-12, none of which survived the full +# pipeline. + +N_REM_LATENCY_FEATURES = 1 + + +def extract_rem_latency_features(algo_data): + """ + REM latency: time from sleep onset (first epoch in any sleep stage, + N1/N2/N3/REM) to the first REM epoch, in minutes. + + Returns np.ndarray of length N_REM_LATENCY_FEATURES (1): + [0] REM latency (minutes) + + NaN fallback: + - Missing CAISR annotations → NaN + - No sleep at all → NaN (latency undefined without a sleep onset) + - No REM occurs anywhere in the recording → NaN (a real clinical + possibility — some CI+ patients may have zero REM — but + "latency to an event that never happens" is undefined, not + zero or infinite; do not silently encode it as either) + """ + if not algo_data: + return np.full(N_REM_LATENCY_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) == 0: + return np.full(N_REM_LATENCY_FEATURES, float('nan')) + + sleep_mask = (valid_stages >= 1) & (valid_stages <= 4) + if not sleep_mask.any(): + return np.full(N_REM_LATENCY_FEATURES, float('nan')) + + first_sleep_idx = np.where(sleep_mask)[0][0] + + rem_mask = (valid_stages == 4) + if not rem_mask.any(): + return np.full(N_REM_LATENCY_FEATURES, float('nan')) + + first_rem_idx = np.where(rem_mask)[0][0] + + if first_rem_idx < first_sleep_idx: + # Should not happen (REM before any sleep stage is detected is + # incoherent given REM is itself a sleep stage) — guard against + # a data/annotation anomaly rather than returning a negative + # latency. + return np.full(N_REM_LATENCY_FEATURES, float('nan')) + + latency_epochs = first_rem_idx - first_sleep_idx + latency_minutes = latency_epochs * 30.0 / 60.0 + + return np.array([latency_minutes], dtype=float) + + +# ── No-arousal confidence entropy — NOT wired in ────────────────────────── +# Planned since Week 2-3 (FEATURES.md, "Always" gate), never implemented. +# Mechanically identical to the existing, already-shipped N3_entropy +# feature (caisr_enriched.py index 10) — same binary-entropy calculation, +# applied to caisr_prob_no-ar instead of caisr_prob_n3. Cheapest possible +# thing to test since it reuses proven, working logic verbatim. + +N_NO_AROUSAL_ENTROPY_FEATURES = 1 + + +def extract_no_arousal_entropy_features(algo_data): + """ + Binary entropy on caisr_prob_no-ar: higher = more ambiguous + arousal/no-arousal staging confidence. Identical calculation to the + shipped N3_entropy feature, different source channel. + + Returns np.ndarray of length N_NO_AROUSAL_ENTROPY_FEATURES (1). + + NaN fallback: missing CAISR annotations, or the caisr_prob_no-ar + channel not present in this dataset (key name per FEATURES.md's + original plan — not yet confirmed against real data; if this key + doesn't match what's actually in the CAISR annotation EDFs, this + will silently return NaN for every patient, which the NaN-rate + check in the EDA script will surface immediately). + """ + if not algo_data: + return np.full(N_NO_AROUSAL_ENTROPY_FEATURES, float('nan')) + + prob_no_ar_arr = algo_data.get('caisr_prob_no-ar', np.array([])) + if len(prob_no_ar_arr) == 0: + return np.full(N_NO_AROUSAL_ENTROPY_FEATURES, float('nan')) + + p = np.clip(np.asarray(prob_no_ar_arr, dtype=float), 1e-9, 1.0 - 1e-9) + h = -(p * np.log(p) + (1.0 - p) * np.log(1.0 - p)) + + return np.array([float(np.mean(h))], dtype=float) + +# ── Gating design decision (2026-07-14) ─────────────────────────────────── +# The 2026-07-12 learning_log.md entry describes a duration guard +# (_MIN_RECORDING_HOURS / _recording_too_short()) as added to this file. +# Confirmed against the actual repo: it was never implemented — this was +# a genuinely undecided design question at the time (exclude short +# recordings outright, vs. keep specific components that may still carry +# usable signal even when the whole-night measure doesn't), not a lost +# commit. Resolved here, but NOT by resurrecting a blanket duration cutoff. +# +# For NF1 specifically, duration is the wrong thing to gate on. The real +# constraint is estimability of the transition matrix itself: a short +# recording (e.g. ~12 minutes, ~24 epochs) yields ~23 transitions to +# populate a 5x5 matrix — most cells land at 0 or 1 count. The stationary +# distribution, solved via eigendecomposition of that matrix, isn't just +# "less confident" at that point; a sparse/disconnected empirical chain +# can produce a numerically unstable or meaningless eigenvector at +# eigenvalue 1. That's a different failure mode than "weak signal," and +# gating on wall-clock duration is a proxy for it, not the thing itself — +# a technically-short recording with a clean, well-sampled stage sequence +# would be discarded for no real reason under a duration cutoff. +# +# Fix: gate NF1 on total valid TRANSITION COUNT directly (_MIN_TRANSITIONS +# below), the same way NF2 already gates on arousal EVENT COUNT +# (_MIN_AROUSAL_EVENTS) rather than duration. This targets whether the +# component is estimable, not how long the file happens to be. +# _MIN_TRANSITIONS=40 is a starting value, not empirically validated — +# worth a quick synthetic check (does the stationary-distribution +# eigensolve behave reasonably at exactly 40 transitions?) before +# treating it as final. +# +# Given the original duration audit found effect sizes already near zero +# and only ~1.3% of small-set patients affected either way, this is a +# correctness fix, not expected to move the ablation result — but it's +# the right way to gate these two specific features regardless. + +_MIN_TRANSITIONS = 40 + + +# ── NF1: Markov stage-transition features — NOT wired in ───────────────── +# Per TEST_PLAN.md's dependency graph: RB1 came back null (2026-07-13, +# rem_latency at C=1.0, max |z|=0.13sigma across folds), which rules out +# the regularization-budget explanation for the 3-for-3 CAISR null +# pattern. Team decision: pursue branch (a) — test NF1+NF2 as a COMBINED +# feature set rather than one-at-a-time, since single-feature ablation +# is structurally blind to interaction effects. This deliberately breaks +# the project's one-variable-per-test discipline, per TEST_PLAN.md's +# explicit allowance for this branch. +# +# CAISR-annotation-only. Stage codes per this project's established +# convention (see extract_rem_latency_features, extract_wake_fragmentation +# _features above): Wake=5, N1=3, N2=2, N3=1, REM=4. + +N_MARKOV_FEATURES = 9 + +# State code -> index mapping, fixed order: Wake, N1, N2, N3, REM. +_STATE_CODES = [5, 3, 2, 1, 4] +_STATE_LABELS = ['Wake', 'N1', 'N2', 'N3', 'REM'] +_STATE_IDX = {code: i for i, code in enumerate(_STATE_CODES)} + + +def extract_markov_transition_features(algo_data): + """ + Stage-transition dynamics via a first-order Markov chain fitted to + the epoch-level stage sequence. Static stage percentages (already + shipped: Wake_pct, N1_pct, N2_pct, N3_pct, REM_pct) don't capture + HOW a patient moves between stages — this does. + + Returns np.ndarray of length N_MARKOV_FEATURES (9): + [0] transition_entropy — Shannon entropy (nats) of the + flattened joint transition-probability distribution + P(s_t, s_t+1) over all 25 (5x5) state pairs. This is the + entropy of the WHOLE transition matrix, not a per-row + average — a deliberate choice: per-row entropy averaged by + visit frequency would double-count the marginal stage + distribution already captured by [4:9] below, whereas the + joint-pair entropy captures how predictable/chaotic the + SEQUENCE of transitions is, which is the actually-new + information this feature is meant to add. + [1] n3_to_wake_rate (events/hr) — reported empirically, NOT + assumed rare or pathological a priori (per TEST_PLAN.md's + explicit caution against over-framing this transition). + [2] rem_to_n3_rate (events/hr) — same caution applies. + [3] escalating_rate (events/hr) — rate of the exact 3-step + sequence N3->N2->N1->Wake (4 consecutive epochs in that + order), the "smooth de-escalation" pattern contrasted + against the abrupt transitions above. + [4:9] stationary distribution — the fitted Markov chain's + stationary distribution over [Wake, N1, N2, N3, REM], + summing to 1. Solved via the row-normalized transition + matrix's left eigenvector at eigenvalue 1 (pi P = pi), + NOT simply the raw marginal stage-epoch proportions — those + are only equal in the infinite-sample/perfectly-ergodic + limit. Falls back to raw marginal proportions if the + eigen-solve fails (e.g. an absorbing/disconnected + empirical chain from a short or unusual recording). + + REQUIRED before ablation-testing [4:9] as independent: check + correlation against the existing Wake_pct/N1_pct/N2_pct/N3_pct/ + REM_pct features on the small set (per TEST_PLAN.md's NF1 internal + dependency). Not computed here — this function only extracts the + feature; the correlation check belongs in the ablation/EDA script + so it's visible per-run, not buried in extraction code. + + NaN fallback (all 9 values): + - Missing CAISR annotations + - Fewer than 2 valid stage epochs (can't form even 1 transition) + - Fewer than _MIN_TRANSITIONS (40) total valid transitions — + gated on estimability of the transition matrix itself, not + wall-clock duration (see module-level "Gating design decision" + note above for why duration was rejected as the gate here). + """ + if not algo_data: + return np.full(N_MARKOV_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + mask = np.isin(valid_stages, _STATE_CODES) + seq = valid_stages[mask] + + if len(seq) < 2: + return np.full(N_MARKOV_FEATURES, float('nan')) + + n_transitions = len(seq) - 1 + if n_transitions < _MIN_TRANSITIONS: + return np.full(N_MARKOV_FEATURES, float('nan')) + + total_hours = len(seq) * 30.0 / 3600.0 + + cur_idx = np.array([_STATE_IDX[c] for c in seq[:-1]]) + nxt_idx = np.array([_STATE_IDX[c] for c in seq[1:]]) + + trans_counts = np.zeros((5, 5)) + for i in range(len(cur_idx)): + trans_counts[cur_idx[i], nxt_idx[i]] += 1 + + total_trans = trans_counts.sum() + if total_trans == 0: + return np.full(N_MARKOV_FEATURES, float('nan')) + + # [0] Transition entropy — joint distribution over all (i,j) pairs. + p_joint = trans_counts / total_trans + flat_nonzero = p_joint.flatten() + flat_nonzero = flat_nonzero[flat_nonzero > 0] + transition_entropy = float(-np.sum(flat_nonzero * np.log(flat_nonzero))) + + # [1-2] Abrupt transition rates — reported as-is, no a priori framing. + n3_to_wake_rate = trans_counts[_STATE_IDX[1], _STATE_IDX[5]] / total_hours + rem_to_n3_rate = trans_counts[_STATE_IDX[4], _STATE_IDX[1]] / total_hours + + # [3] Escalating rate — exact 4-epoch N3->N2->N1->Wake sequence. + esc_count = 0 + for i in range(len(seq) - 3): + if seq[i] == 1 and seq[i+1] == 2 and seq[i+2] == 3 and seq[i+3] == 5: + esc_count += 1 + escalating_rate = esc_count / total_hours + + # [4-8] Stationary distribution via left eigenvector at eigenvalue 1. + row_sums = trans_counts.sum(axis=1, keepdims=True) + row_sums_safe = np.where(row_sums == 0, 1, row_sums) + p_rownorm = trans_counts / row_sums_safe + + stationary = None + try: + eigvals, eigvecs = np.linalg.eig(p_rownorm.T) + best_idx = int(np.argmin(np.abs(eigvals - 1))) + vec = np.real(eigvecs[:, best_idx]) + vec = np.clip(vec, 0, None) + if vec.sum() > 1e-9: + stationary = vec / vec.sum() + except np.linalg.LinAlgError: + stationary = None + + if stationary is None: + # Fallback: raw marginal proportion, same state order. + marginal_counts = np.array( + [np.sum(seq == code) for code in _STATE_CODES], dtype=float) + stationary = marginal_counts / marginal_counts.sum() + + return np.array([ + transition_entropy, n3_to_wake_rate, rem_to_n3_rate, escalating_rate, + *stationary, + ], dtype=float) + + +# ── NF2: Micro-arousal temporal clustering — NOT wired in ──────────────── +# Same sequencing note as NF1: built for the combined-feature test +# (branch (a) after RB1's null), not a standalone single-feature +# ablation candidate. +# +# Gated on event count only (_MIN_AROUSAL_EVENTS), not duration — this +# was already the right design even before the 2026-07-14 gating +# correction above; CV and burst index need enough EVENTS to be +# meaningful, not enough wall-clock time. The duration-guard call this +# patch originally (incorrectly) included has been removed — the +# event-count gate below is sufficient on its own. +# +# ASSUMPTION FLAGGED: arousal_caisr's sample rate is inferred as 2Hz +# (0.5s/sample) from the existing spont_arousal_idx code's +# reshape(-1, 2) pairing against 1s-resolution resp, in +# extract_caisr_enriched_features above. This has NOT been independently +# confirmed against real data / a header spec. It affects burst_index's +# absolute "<3 minute" threshold; it does NOT affect the CV feature +# (unitless ratio, invariant to a constant sample-rate scaling error). +# Confirm before trusting burst_index specifically. + +N_AROUSAL_CLUSTERING_FEATURES = 2 +_MIN_AROUSAL_EVENTS = 5 +_AROUSAL_SAMPLE_HZ = 2.0 # SEE ASSUMPTION FLAG ABOVE — confirm before trusting. + + +def extract_arousal_clustering_features(algo_data): + """ + Temporal clustering of arousal events — a different statistical + property of the arousal-event sequence than the rate-only + Arousal_idx/Spontaneous_arousal_idx already shipped. Tightly + clustered arousals (thalamocortical gating failure, hypothesized) + vs. evenly-spaced arousals (benign) could carry signal a rate + measure can't see. + + Returns np.ndarray of length N_AROUSAL_CLUSTERING_FEATURES (2): + [0] cv_inter_arousal_interval — coefficient of variation + (std/mean) of the time gaps between consecutive arousal + event onsets. Unitless — robust to the sample-rate + assumption above. + [1] burst_index — fraction of inter-arousal intervals under + 3 minutes (180s). Depends on the sample-rate assumption + above for its absolute threshold. + + REQUIRED before ablation-testing as independent: check correlation + against Arousal_idx / Spontaneous_arousal_idx on the small set (per + TEST_PLAN.md's NF2 internal dependency) — not computed here, same + reasoning as NF1's collinearity note above. + + NaN fallback (both values): + - Missing CAISR annotations + - Fewer than _MIN_AROUSAL_EVENTS (5) arousal events detected — + the sole gate for this feature; CV and burst fraction aren't + meaningful from 1-4 events regardless of recording length. + """ + if not algo_data: + return np.full(N_AROUSAL_CLUSTERING_FEATURES, float('nan')) + + arousal = algo_data.get('arousal_caisr', np.array([])) + if len(arousal) == 0: + return np.full(N_AROUSAL_CLUSTERING_FEATURES, float('nan')) + + binary = (np.asarray(arousal) > 0).astype(int) + edges = np.diff(binary, prepend=0) + start_indices = np.where(edges == 1)[0] + + if len(start_indices) < _MIN_AROUSAL_EVENTS: + return np.full(N_AROUSAL_CLUSTERING_FEATURES, float('nan')) + + times_sec = start_indices / _AROUSAL_SAMPLE_HZ + intervals = np.diff(times_sec) + + if len(intervals) == 0: + return np.full(N_AROUSAL_CLUSTERING_FEATURES, float('nan')) + + mu = float(np.mean(intervals)) + sigma = float(np.std(intervals)) + cv = sigma / mu if mu > 0 else float('nan') + burst_index = float(np.mean(intervals < 180.0)) + + return np.array([cv, burst_index], dtype=float) + +# ── Backup tier rank 1: Time-of-night-conditional transitions ──────────── +# Per TESTING_STRATEGY_post_entry5.md: opened after the trimmed NF1+NF2 +# result (2026-07-15) resolved to the "branch (a) exhausted for this +# feature family" gate row (S0001 AUROC flat at +0.18σ, reward pattern +# shifted rather than strengthened). This is rank 1 of the backup tier — +# cheapest to build (reuses the exact night-segmenting logic already +# established by the shipped N3-gradient feature), tried before +# conceding the transition-dynamics direction entirely. +# +# NOT wired into extract_caisr_enriched_features — same discipline as +# every other candidate function in this file. +# +# Reuses _STATE_CODES / _STATE_IDX from the NF1 section above this one — +# do not redefine; this function assumes NF1's code block already exists +# earlier in this same file. + +N_TIME_OF_NIGHT_FEATURES = 2 + +# Minimum occurrences of the "from" stage WITHIN a given third-of-night +# segment to trust a self-transition probability computed from it. This +# is deliberately smaller than NF1's _MIN_TRANSITIONS (40, for a whole- +# night 5x5 matrix) because each of these two features only needs ONE +# specific stage's self-transition count within ONE third of the night, +# not a full matrix — a much lower bar for what counts as "enough data." +# Not empirically validated, same caveat as _MIN_TRANSITIONS. +_MIN_SEGMENT_FROM_COUNT = 3 + + +def _segment_stage_sequence(seq, n_segments=3): + """ + Split an already-filtered stage-code sequence into n_segments + contiguous, equal-length chunks by epoch COUNT (not by target + duration) — the same segmenting philosophy as the shipped + n3_gradient feature (which splits at len(valid_stages)//2), just + generalized to thirds instead of halves. + """ + n = len(seq) + edges = np.linspace(0, n, n_segments + 1).astype(int) + return [seq[edges[i]:edges[i + 1]] for i in range(n_segments)] + + +def _self_transition_prob(segment, stage_code, min_from_count=_MIN_SEGMENT_FROM_COUNT): + """ + P(stage_code -> stage_code) within this segment specifically — i.e. + given the patient was in `stage_code` at epoch t (within this + segment), what fraction of the time were they still in it at t+1. + NaN if the segment never visited this stage enough times to trust + the estimate (per min_from_count), not just if the segment is short + overall — a segment could be long but rarely visit N3, for instance. + """ + if len(segment) < 2: + return float('nan') + from_mask = segment[:-1] == stage_code + n_from = int(from_mask.sum()) + if n_from < min_from_count: + return float('nan') + stayed = int(np.sum(from_mask & (segment[1:] == stage_code))) + return stayed / n_from + + +def extract_time_of_night_transition_features(algo_data): + """ + Differential transition features comparing early-night vs. late- + night stage persistence — targets whether a patient's TRANSITION + BEHAVIOR itself shifts over the night, not just their overall stage + percentages (which the shipped n3_gradient feature already partly + captures, via raw N3 presence rather than transition probability). + + The night is split into three equal-epoch-count thirds (early, mid, + late); the middle third is computed but not used by either feature + below — reserved for future extensions, not dead code by mistake. + + Returns np.ndarray of length N_TIME_OF_NIGHT_FEATURES (2): + [0] homeostatic_decay_index = + P_early[N3->N3] - P_late[N3->N3] + Healthy sleep: N3 self-transition (staying in deep sleep + once there) should be HIGH early in the night and DECLINE + by the end, as homeostatic sleep pressure dissipates. A + positive value is the expected/healthy direction; a value + near zero or negative (N3 persistence NOT declining, or + even increasing, late in the night) is the hypothesized + CI-risk pattern — not assumed true, this is what the + ablation is meant to test. + [1] rem_escalation_delta = + P_late[REM->REM] - P_early[REM->REM] + Healthy sleep: REM self-transition (staying in REM once + there) should INCREASE over the night as REM-drive + dominates the later cycles. A positive value is the + expected/healthy direction. + + REQUIRED before ablation-testing [0] as independent: check + correlation against the existing n3_gradient feature (caisr_enriched + index 8, "N3 first-half/second-half ratio") — both features target + a similar early-vs-late N3 concept via different mechanisms (raw + presence ratio vs. self-transition probability), so this is a + genuine collinearity risk, not a formality. Not computed here — see + the ablation script. + + NaN fallback (either or both values independently): + - Missing CAISR annotations -> both NaN + - Fewer than 2 valid stage epochs -> both NaN + - The relevant stage's self-transition probability is undefined + in the relevant segment (fewer than _MIN_SEGMENT_FROM_COUNT + occurrences of that stage as a "from" state in that third of + the night) -> that specific value NaN, independently of the + other one (e.g. homeostatic_decay_index can be NaN while + rem_escalation_delta is a real number, if N3 was too sparse + early or late but REM wasn't) + """ + if not algo_data: + return np.full(N_TIME_OF_NIGHT_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + mask = np.isin(valid_stages, _STATE_CODES) + seq = valid_stages[mask] + + if len(seq) < 2: + return np.full(N_TIME_OF_NIGHT_FEATURES, float('nan')) + + early, mid, late = _segment_stage_sequence(seq, n_segments=3) + + p_n3_early = _self_transition_prob(early, 1) # N3 code = 1 + p_n3_late = _self_transition_prob(late, 1) + p_rem_early = _self_transition_prob(early, 4) # REM code = 4 + p_rem_late = _self_transition_prob(late, 4) + + homeostatic_decay = ( + p_n3_early - p_n3_late + if not (np.isnan(p_n3_early) or np.isnan(p_n3_late)) + else float('nan') + ) + rem_escalation = ( + p_rem_late - p_rem_early + if not (np.isnan(p_rem_late) or np.isnan(p_rem_early)) + else float('nan') + ) + + return np.array([homeostatic_decay, rem_escalation], dtype=float) \ No newline at end of file diff --git a/features/demographic.py b/features/demographic.py new file mode 100644 index 0000000..125688e --- /dev/null +++ b/features/demographic.py @@ -0,0 +1,49 @@ +# features/demographic.py +# Demographic feature extraction — features 0-9 +# +# Iteration history: +# v1 (2026-06-30): Initial — age, sex onehot, race onehot, BMI (10 features) +# v1.1 (2026-07-01): Docstring corrections only, no logic change — +# - Age is NOT capped at 90. helper_code.py's load_age() returns a raw +# float(age) with NaN fallback; no clipping happens anywhere in the +# pipeline. The prior docstring was aspirational/stale, not a bug in +# behavior, but it matters: loso_cv.py's compute_prevalence() keys on +# exact age values, and a phantom cap would have implied a mismatch +# between the age used as a model feature and the age used for reward/ +# prevalence lookups that doesn't actually exist. +# - BMI NaN rate corrected: Phase 1 EDA found 75.9% missing, not ~15%. + +import numpy as np +from helper_code import load_age, load_sex, load_bmi, load_race + + +def extract_demographic_features(data): + """ + Extracts and encodes demographic features from a metadata dictionary. + + Returns np.ndarray of length 10: + [0] Age (continuous — not capped; NaN if missing/unparseable) + [1-3] Sex one-hot (Female, Male, Unknown) + [4-8] Race one-hot (Asian, Black, Others, Unavailable, White) + [9] BMI (continuous — 75.9% NaN per Phase 1 EDA, imputed downstream) + """ + age = np.array([load_age(data)]) + + sex = load_sex(data, standardize=True) + sex_vec = np.zeros(3) + if sex == 'Female': sex_vec[0] = 1 + elif sex == 'Male': sex_vec[1] = 1 + else: sex_vec[2] = 1 + + race = load_race(data, standardize=True) + race_vec = np.zeros(5) + if race == 'Asian': race_vec[0] = 1 + elif race == 'Black': race_vec[1] = 1 + elif race == 'Others': race_vec[2] = 1 + elif race == 'Unavailable': race_vec[3] = 1 + elif race == 'White': race_vec[4] = 1 + else: race_vec[2] = 1 # default to Others + + bmi = np.array([load_bmi(data)]) + + return np.concatenate([age, sex_vec, race_vec, bmi]) \ No newline at end of file diff --git a/features/feature_selection.py b/features/feature_selection.py new file mode 100644 index 0000000..f3921f7 --- /dev/null +++ b/features/feature_selection.py @@ -0,0 +1,69 @@ +# features/feature_selection.py +# Team Narnia — PhysioNet Challenge 2026 +# +# Entry 8: drops 11 SIGN-FLIPPING features (coefficients that change +# direction depending on which 2 training sites fit the model — direct, +# measured evidence of site-specific overfitting rather than genuine +# biological signal; see the cross-fold coefficient analysis, 2026-07-19, +# and learning_log.md for the full derivation). +# +# This is "Stage 1" of the two-stage drop tested in tools/ablation_drop_ +# signflip_features.py. Validated: +8.2% relative reward, mean age- +# conditioned AUROC 0.6459 -> 0.6711 (+0.0252), ALL THREE LOSO folds +# improved simultaneously (not a mixed result), cross-fold spread shrank +# 0.1264 -> 0.1145 (more stable, not just better on average — the +# specific signature the sign-flip theory predicted). Reproduced bit- +# for-bit across an independent Kaggle kernel restart. +# +# Stage 2 (CA_rate_age_residual, Race_Black/White/Asian) is deliberately +# NOT included here. CA_rate_age_residual carries separate validation +# history (Entry 4, BH-FDR significance on residual AUROC — a different +# question than raw-coefficient sign stability). The Race_* features +# carry real, documented demographic signal (dementia diagnosis odds +# genuinely differ by race in the published literature) that a numerically +# unstable one-hot coefficient does not necessarily invalidate. Stage 2 +# showed a further, real improvement in LOSO (clears the pre-committed +# 1.0-sigma gate outright, vs. Stage 1's 0.97 just under it) but was held +# back given the extra scrutiny both warrant and the proximity of the +# deadline — not promoted to production as of Entry 8. Revisit only with +# a specific, separate justification for each, not just because Stage 1 +# worked. +# +# Promoted into this shared features/ module (rather than living only in +# the ablation script) so team_code.py and any future validation run +# share exactly the same drop list — same precedent as AgeResidualizer +# and compute_value_weighted_sample_weights. + +import numpy as np +from sklearn.base import BaseEstimator, TransformerMixin + +from features import FEATURE_NAMES_50 + +STAGE1_DROP_FEATURES = [ + 'BMI', 'EEG_mob_REM_Wake', 'RERA_rate', 'N1_pct', 'Spont_arousal_idx', + 'EEG_std_REM_Wake', 'EEG_cplx_REM_Wake', 'REM_NREM_AHI_ratio', + 'CA_rate', 'CA_total_ratio', 'EEG_zcr_REM_Wake', +] + + +class FeatureDropper(BaseEstimator, TransformerMixin): + """ + Drops named columns from the 50-length (post-AgeResidualizer) vector. + MUST be inserted immediately after the 'age_residual' pipeline step + and before 'imputer' — AgeResidualizer's own source features (e.g. + CA_rate, used to compute CA_rate_age_residual, which IS kept in + Stage 1) need to still be present when AgeResidualizer runs, even + though the raw CA_rate column itself is dropped from what the + classifier ultimately sees. Dropping it any earlier in the pipeline + would break that computation. + """ + def __init__(self, drop_feature_names): + self.drop_feature_names = list(drop_feature_names) + + def fit(self, X, y=None): + return self + + def transform(self, X): + drop_indices = [FEATURE_NAMES_50.index(name) for name in self.drop_feature_names] + keep = [i for i in range(np.asarray(X).shape[1]) if i not in drop_indices] + return np.asarray(X)[:, keep] \ No newline at end of file diff --git a/features/human.py b/features/human.py new file mode 100644 index 0000000..0f3cd88 --- /dev/null +++ b/features/human.py @@ -0,0 +1,65 @@ +# features/human.py +# Human annotation feature extraction — not used in the model +# +# Human annotations are ONLY available in the training set. +# They are intentionally excluded from the model feature vector — +# run_model() must work without them on the hidden validation/test sets. +# This module is retained for research and EDA purposes only. +# +# Iteration history: +# v1 (2026-06-30): Initial — AHI, arousal, limb, stage%, transitions, +# WASO, REM latency [12 features — not in model] + +import numpy as np + + +def extract_human_annotations_features(human_data): + """ + Extracts features from expert-scored human annotations. + NOT INCLUDED in the model feature vector — for research/EDA only. + Returns np.ndarray of length 12. + """ + if not human_data or 'resp_expert' not in human_data: + return np.full(12, float('nan')) + + features = [] + total_seconds = len(human_data.get('resp_expert', [])) + total_hours = total_seconds / 3600.0 + + def count_events(key): + if key not in human_data or total_hours <= 0: + return float('nan') + sig = (human_data[key] > 0).astype(int) + edges = np.diff(sig, prepend=0) + return np.count_nonzero(edges == 1) / total_hours + + features.extend([ + count_events('resp_expert'), + count_events('arousal_expert'), + count_events('limb_expert'), + ]) + + stages = human_data.get('stage_expert', np.array([])) + valid = stages[stages < 9.0] if len(stages) > 0 else np.array([]) + if len(valid) > 0: + features.extend([ + float(np.mean(valid == 5)), + float(np.mean(valid == 4)), + float(np.mean(valid == 3)), + float(np.mean(valid == 2)), + float(np.mean(valid == 1)), + float(np.mean(valid > 0)), + ]) + else: + features.extend([float('nan')] * 6) + + if len(valid) > 1: + features.extend([ + float(np.count_nonzero(np.diff(valid)) / total_hours), + float(np.count_nonzero(valid == 0) * 30 / 60.0), + float(np.where(valid == 4)[0][0]) if np.any(valid == 4) else float('nan'), + ]) + else: + features.extend([float('nan')] * 3) + + return np.array(features) \ No newline at end of file diff --git a/features/physiological.py b/features/physiological.py new file mode 100644 index 0000000..f3a9fdd --- /dev/null +++ b/features/physiological.py @@ -0,0 +1,129 @@ +# features/physiological.py +# Physiological signal feature extraction — features 10-58 +# +# Iteration history: +# v1 (2026-06-30): Initial — 7 Hjorth parameters × 7 lead types (49 features) +# + + +import numpy as np +import os +from helper_code import ( + load_rename_rules, standardize_channel_names_rename_only, + derive_bipolar_signal +) + +# Resolve channel_table.csv relative to this file's location +_FEATURES_DIR = os.path.dirname(os.path.abspath(__file__)) +_REPO_DIR = os.path.dirname(_FEATURES_DIR) +DEFAULT_CSV_PATH = os.path.join(_REPO_DIR, 'channel_table.csv') + + +def extract_physiological_features(physiological_data, physiological_fs, + csv_path=DEFAULT_CSV_PATH): + """ + Standardises channel names, derives bipolar signals, and computes + 7 Hjorth-style time-domain features for 7 lead types. + + Lead order: EEG, EOG, ChinEMG, LegEMG, ECG, Resp, SpO2 + Feature order per lead: std, MAV, ZCR, RMS, variance, mobility, complexity + + Returns np.ndarray of length 49 (7 features × 7 leads). + NaN padding used for unavailable leads. + """ + original_labels = list(physiological_data.keys()) + + rename_rules = load_rename_rules(os.path.abspath(csv_path)) + rename_map, cols_to_drop = standardize_channel_names_rename_only( + original_labels, rename_rules) + + processed_channels = {} + processed_fs = {} + for old_label, data in physiological_data.items(): + if old_label in cols_to_drop: + continue + new_label = rename_map.get(old_label, old_label.lower()) + processed_channels[new_label] = data + if old_label in physiological_fs: + processed_fs[new_label] = physiological_fs[old_label] + else: + raise KeyError(f'No sampling rate for channel: {old_label}') + + if 'physiological_data' in locals(): + del physiological_data + + # Bipolar derivations + bipolar_configs = [ + ('f3-m2', 'f3', ['m2']), + ('f4-m1', 'f4', ['m1']), + ('c3-m2', 'c3', ['m2']), + ('c4-m1', 'c4', ['m1']), + ('o1-m2', 'o1', ['m2']), + ('o2-m1', 'o2', ['m1']), + ('e1-m2', 'e1', ['m2']), + ('e2-m1', 'e2', ['m1']), + ('chin1-chin2', 'chin 1', ['chin 2']), + ('lat', 'lleg+', ['lleg-']), + ('rat', 'rleg+', ['rleg-']), + ] + + for target, pos, neg_list in bipolar_configs: + if target in processed_channels or pos not in processed_channels: + continue + if not all(n in processed_channels for n in neg_list): + continue + all_involved = [pos] + neg_list + fs_vals = [processed_fs[ch] for ch in all_involved] + if len(set(fs_vals)) > 1: + raise ValueError( + f'Sampling rate mismatch for {target}: ' + f'{dict(zip(all_involved, fs_vals))}') + ref = (processed_channels[neg_list[0]] if len(neg_list) == 1 + else tuple(processed_channels[n] for n in neg_list)) + derived = derive_bipolar_signal(processed_channels[pos], ref) + if derived is not None: + processed_channels[target] = derived + processed_fs[target] = processed_fs[pos] + + # Lead selection — first available channel wins per lead type + leads_to_check = { + 'eeg': ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1'], + 'eog': ['e1-m2', 'e2-m1'], + 'chin': ['chin1-chin2', 'chin'], + 'leg': ['lat', 'rat'], + 'ecg': ['ecg', 'ekg'], + 'resp': ['airflow', 'ptaf', 'abd', 'chest'], + 'spo2': ['spo2', 'sao2'], + } + + final_features = [] + for lead_type, candidates in leads_to_check.items(): + sig = None + for candidate in candidates: + if candidate in processed_channels and processed_channels[candidate] is not None: + sig = processed_channels[candidate] + break + + if sig is not None and len(sig) > 1: + std_val = np.std(sig) + mav_val = np.mean(np.abs(sig)) + zcr = np.mean(np.diff(np.sign(sig)) != 0) + rms = np.sqrt(np.mean(sig ** 2)) + activity = np.var(sig) + diff1 = np.diff(sig) + mobility = (np.sqrt(np.var(diff1) / activity) + if activity > 0 else 0.0) + diff2 = np.diff(diff1) + var_d1 = np.var(diff1) + var_d2 = np.var(diff2) + complexity = ((np.sqrt(var_d2 / var_d1) / mobility) + if (var_d1 > 0 and mobility > 0) else 0.0) + final_features.extend( + [std_val, mav_val, zcr, rms, activity, mobility, complexity]) + else: + final_features.extend([float('nan')] * 7) + + if 'processed_channels' in locals(): + del processed_channels + + return np.array(final_features) \ No newline at end of file diff --git a/features/physiological_ratios.py b/features/physiological_ratios.py new file mode 100644 index 0000000..e32961e --- /dev/null +++ b/features/physiological_ratios.py @@ -0,0 +1,220 @@ +# features/physiological_ratios.py +# Stage-conditional ratio features — replaces absolute Hjorth (entry 2+) +# +# Iteration history: +# v2 (2026-07-01): LOSO ablation showed absolute Hjorth features hurt +# cross-site AUROC by +0.09 (I0006) and +0.03 (S0001). +# Root cause: absolute signal amplitude is equipment- +# dependent. Within-recording ratios cancel this out. +# +# All 15 features are ratios within a single patient's recording. +# Equipment differences affect numerator and denominator equally → cancel. +# +# Index range: 33-47 in entry 2 feature vector. + +import numpy as np +import os +from helper_code import ( + load_rename_rules, standardize_channel_names_rename_only, + derive_bipolar_signal +) + +_FEATURES_DIR = os.path.dirname(os.path.abspath(__file__)) +_REPO_DIR = os.path.dirname(_FEATURES_DIR) +DEFAULT_CSV_PATH = os.path.join(_REPO_DIR, 'channel_table.csv') + +N_RATIO_FEATURES = 15 # EEG N3/Wake (7) + EEG REM/Wake (7) + Chin atonia (1) + +# Minimum recording seconds per stage to trust a ratio +_MIN_STAGE_SECONDS = 120 # 2 minutes + +# Cap on ratio values — prevents extreme outliers from dominating +_RATIO_CAP = 20.0 + + +# ── Hjorth parameter vector ─────────────────────────────────────────────────── +def _hjorth(sig): + """ + Compute 7 Hjorth-style features for a 1-D signal array. + Returns np.ndarray(7,) or NaN array if signal is too short. + """ + if sig is None or len(sig) < 10: + return np.full(7, float('nan')) + activity = float(np.var(sig)) + if activity < 1e-20: + return np.array([0., 0., 0., 0., 0., 0., 0.], dtype=float) + diff1 = np.diff(sig) + var_d1 = float(np.var(diff1)) + mobility = float(np.sqrt(var_d1 / activity)) if activity > 0 else 0. + diff2 = np.diff(diff1) + var_d2 = float(np.var(diff2)) + complexity = float(np.sqrt(var_d2 / var_d1) / mobility) \ + if (var_d1 > 0 and mobility > 0) else 0. + return np.array([ + float(np.std(sig)), + float(np.mean(np.abs(sig))), + float(np.mean(np.diff(np.sign(sig)) != 0)), + float(np.sqrt(np.mean(sig ** 2))), + activity, + mobility, + complexity, + ], dtype=float) + + +def _safe_ratio(num_vec, denom_vec): + """ + Element-wise ratio with NaN guards and cap. + Returns NaN wherever either input is NaN or denominator ≈ 0. + """ + out = np.full(len(num_vec), float('nan')) + for i, (n, d) in enumerate(zip(num_vec, denom_vec)): + if np.isnan(n) or np.isnan(d) or abs(d) < 1e-12: + continue + out[i] = float(np.clip(n / d, -_RATIO_CAP, _RATIO_CAP)) + return out + + +def _get_channel(channels, candidates): + """Return first available channel from candidates list.""" + for c in candidates: + if c in channels and channels[c] is not None: + return channels[c] + return None + + +def _get_fs(fs_map, candidates, default=200.0): + for c in candidates: + if c in fs_map: + return float(fs_map[c]) + return default + + +# ── Main extraction function ────────────────────────────────────────────────── +def extract_physiological_ratio_features(phys_data, phys_fs, algo_data, + csv_path=DEFAULT_CSV_PATH): + """ + Stage-conditional ratio features requiring both physiological EDF + (raw signal) and CAISR annotations (stage timing). + + All features are within-recording ratios — site-stable by construction + because equipment scale differences cancel in numerator and denominator. + + Returns np.ndarray of length 15: + [0-6] EEG Hjorth(N3) / Hjorth(Wake) — slow-wave quality ratio + std, MAV, ZCR, RMS, var, mobility, complexity + Interpretation: <1 means N3 signal weaker than wake (CI-risk) + >1 means N3 signal stronger (healthy) + + [7-13] EEG Hjorth(REM) / Hjorth(Wake) — REM EEG quality ratio + Same 7 parameters. REM theta/alpha signature vs wake. + + [14] Chin MAV(REM) / Chin MAV(NREM) — atonia quality ratio + Healthy: near 0 (atonia in REM). + CI-risk: elevated (REM Behaviour Disorder marker). + Requires physiological EDF + stage_caisr. + + NaN fallback: + - Missing physio EDF → all 15 NaN + - Missing CAISR → all 15 NaN + - Insufficient samples in a stage (<120s) → that ratio NaN + """ + if not phys_data or not algo_data: + return np.full(N_RATIO_FEATURES, float('nan')) + + # ── Stage annotations ───────────────────────────────────────────────────── + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) < 10: + return np.full(N_RATIO_FEATURES, float('nan')) + + # ── Standardise channel names ───────────────────────────────────────────── + original_labels = list(phys_data.keys()) + rename_rules = load_rename_rules(os.path.abspath(csv_path)) + rename_map, cols_to_drop = standardize_channel_names_rename_only( + original_labels, rename_rules) + + channels = {} + fs_map = {} + for old_label, data in phys_data.items(): + if old_label in cols_to_drop: + continue + new_label = rename_map.get(old_label, old_label.lower()) + channels[new_label] = data + if old_label in phys_fs: + fs_map[new_label] = phys_fs[old_label] + + # ── Bipolar derivations for EEG and Chin ────────────────────────────────── + for target, pos, neg_list in [ + ('f3-m2', 'f3', ['m2']), + ('f4-m1', 'f4', ['m1']), + ('c3-m2', 'c3', ['m2']), + ('c4-m1', 'c4', ['m1']), + ('chin1-chin2', 'chin 1', ['chin 2']), + ]: + if target in channels or pos not in channels: + continue + if not all(n in channels for n in neg_list): + continue + ref = (channels[neg_list[0]] if len(neg_list) == 1 + else tuple(channels[n] for n in neg_list)) + derived = derive_bipolar_signal(channels[pos], ref) + if derived is not None: + channels[target] = derived + fs_map[target] = fs_map.get(pos, 200.0) + + # ── EEG stage-conditional Hjorth ratios [0-13] ──────────────────────────── + result = np.full(N_RATIO_FEATURES, float('nan')) + + eeg_sig = _get_channel(channels, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + eeg_fs = _get_fs(fs_map, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + + if eeg_sig is not None and len(eeg_sig) > 0: + spe = int(round(30.0 * eeg_fs)) # samples per 30s epoch + min_smp = int(_MIN_STAGE_SECONDS * eeg_fs) + + # Upsample stage to EEG sample rate + stage_eeg = np.repeat(valid_stages, spe)[:len(eeg_sig)] + + wake_mask = (stage_eeg == 5) + n3_mask = (stage_eeg == 1) + rem_mask = (stage_eeg == 4) + + def stage_hjorth(mask): + if int(mask.sum()) < min_smp: + return np.full(7, float('nan')) + seg = eeg_sig[:len(mask)][mask] + return _hjorth(seg) + + h_wake = stage_hjorth(wake_mask) + h_n3 = stage_hjorth(n3_mask) + h_rem = stage_hjorth(rem_mask) + + # [0-6] N3 / Wake ratios + result[0:7] = _safe_ratio(h_n3, h_wake) + # [7-13] REM / Wake ratios + result[7:14] = _safe_ratio(h_rem, h_wake) + + # ── Chin atonia ratio [14] ──────────────────────────────────────────────── + chin_sig = _get_channel(channels, ['chin1-chin2', 'chin']) + chin_fs = _get_fs(fs_map, ['chin1-chin2', 'chin']) + + if chin_sig is not None and len(chin_sig) > 0: + spe_c = int(round(30.0 * chin_fs)) + min_c = int(_MIN_STAGE_SECONDS * chin_fs) + + stage_chin = np.repeat(valid_stages, spe_c)[:len(chin_sig)] + rem_mask_c = (stage_chin == 4) + nrem_mask_c = (stage_chin >= 1) & (stage_chin <= 3) + + if (int(rem_mask_c.sum()) >= min_c and + int(nrem_mask_c.sum()) >= min_c): + chin_rem = float(np.mean( + np.abs(chin_sig[:len(stage_chin)][rem_mask_c]))) + chin_nrem = float(np.mean( + np.abs(chin_sig[:len(stage_chin)][nrem_mask_c]))) + if chin_nrem > 1e-12: + result[14] = float(np.clip(chin_rem / chin_nrem, + 0., _RATIO_CAP)) + + return result \ No newline at end of file diff --git a/features/pipeline.py b/features/pipeline.py new file mode 100644 index 0000000..68c6566 --- /dev/null +++ b/features/pipeline.py @@ -0,0 +1,455 @@ +# features/pipeline.py +# Shared model pipeline construction — single source of truth for +# team_code.py (the shipped submission) AND loso_cv.py (the validation +# harness). +# +# Root cause this file fixes (2026-07-03): loso_cv.py had its own +# hand-copied Pipeline() construction, predating AgeResidualizer (Entry 4). +# When AgeResidualizer was added to team_code.py, loso_cv.py's copy was +# never updated — LOSO silently kept validating the OLD 48-feature +# pipeline while team_code.py had already moved to 50. Caught because a +# post-"Entry 4" LOSO run matched the logged Entry 3 numbers to 4 decimal +# places across every metric, including per-fold top-feature importance — +# if AgeResidualizer had actually been active, at least some numerical +# drift would be expected (different feature count into the imputer, +# different fitted trees). +# +# Fix: both team_code.py and loso_cv.py now import build_pipeline() from +# here instead of each constructing their own Pipeline inline. There is +# now exactly ONE place that defines "what the model is" — changing it +# here changes it everywhere that matters, and it is structurally +# impossible for the validation harness and the submission to silently +# diverge again the way they just did. +# +# Lives under features/ (not repo root) deliberately: check_submission_files.py +# only resolves `from features.X import ...` / `from features import ...` +# into required-file paths. A root-level shared module would instead be +# misclassified as a third-party pip package by that script's import parser +# and flagged as "missing from requirements.txt" — a confusing false +# warning for something that isn't a package at all. + +# features/pipeline.py +# Shared model pipeline construction — single source of truth for +# team_code.py AND loso_cv.py / reg_sweep.py. +# +# 2026-07-08 addition: build_logreg_pipeline(), for the Ridge/Lasso/ +# ElasticNet regularization sweep. Reuses the same AgeResidualizer -> +# SimpleImputer -> [CalibratedClassifierCV] skeleton as build_pipeline(), +# swapping XGBClassifier for LogisticRegression, with two changes that +# matter specifically for a regularization sweep: +# +# 1. Added StandardScaler between the imputer and the classifier. This +# was NOT verified to already exist in whatever build_logreg_pipeline() +# produced the 0.6002 large-set LOSO result (learning_log.md, +# 2026-07-07) — that code wasn't available when this was written. +# A C sweep is not meaningfully comparable across features without +# scaling first: Age (~0-100), BMI (~15-50), event-rate features +# (~0-30/hr), and ratio features (~0-3) sit on wildly different +# scales, so an unscaled C controls regularization strength +# inconsistently across coefficients. IF a scaler already existed in +# the original build_logreg_pipeline(), the 0.6002 baseline was +# already fit this way and nothing changes. If it did NOT exist, +# this sweep's results are not directly comparable to that 0.6002 +# number at face value — confirm which case you're in before treating +# a sweep win as a clean improvement over the existing baseline. +# 2. class_weight='balanced' instead of an XGBoost-style scale_pos_weight +# — sklearn's LogisticRegression equivalent, same "derive from actual +# data, don't hardcode" principle already established for XGBoost +# (learning_log.md, 2026-06-30). +# +# solver='saga' is required for L1 and ElasticNet penalties; used +# uniformly (including for L2) so penalty type alone varies across the +# sweep, not solver + penalty together. +# features/pipeline.py +# Shared model pipeline construction — single source of truth for +# team_code.py (the shipped submission) AND loso_cv.py (the validation +# harness). +# +# Root cause this file fixes (2026-07-03): loso_cv.py had its own +# hand-copied Pipeline() construction, predating AgeResidualizer (Entry 4). +# When AgeResidualizer was added to team_code.py, loso_cv.py's copy was +# never updated — LOSO silently kept validating the OLD 48-feature +# pipeline while team_code.py had already moved to 50. Caught because a +# post-"Entry 4" LOSO run matched the logged Entry 3 numbers to 4 decimal +# places across every metric, including per-fold top-feature importance — +# if AgeResidualizer had actually been active, at least some numerical +# drift would be expected (different feature count into the imputer, +# different fitted trees). +# +# Fix: both team_code.py and loso_cv.py now import build_pipeline() from +# here instead of each constructing their own Pipeline inline. There is +# now exactly ONE place that defines "what the model is" — changing it +# here changes it everywhere that matters, and it is structurally +# impossible for the validation harness and the submission to silently +# diverge again the way they just did. +# +# Lives under features/ (not repo root) deliberately: check_submission_files.py +# only resolves `from features.X import ...` / `from features import ...` +# into required-file paths. A root-level shared module would instead be +# misclassified as a third-party pip package by that script's import parser +# and flagged as "missing from requirements.txt" — a confusing false +# warning for something that isn't a package at all. + +import numpy as np +from sklearn.pipeline import Pipeline +from sklearn.impute import SimpleImputer +from sklearn.calibration import CalibratedClassifierCV +from sklearn.model_selection import train_test_split +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from xgboost import XGBClassifier + +# Compatibility shim: cv='prefit' was removed in scikit-learn 1.6+ in favor +# of wrapping the fitted estimator in FrozenEstimator. This repo's pinned +# requirements.txt uses scikit-learn==1.5.2 (pre-FrozenEstimator), but +# loso_cv.py is also run in environments with newer scikit-learn (e.g. the +# Kaggle Notebook used for the large-set run, which had 1.8.0 and rejected +# cv='prefit' outright). Try the modern import; fall back to the legacy +# string if it's unavailable, so PrefitCalibratedModel works either way +# without needing to know which environment it's running in. +try: + from sklearn.frozen import FrozenEstimator + _HAS_FROZEN_ESTIMATOR = True +except ImportError: + _HAS_FROZEN_ESTIMATOR = False + +from features.age_residuals import AgeResidualizer +from features import IDX_AGE, IDX_CA_RATE, IDX_EEG_VAR_REM_WAKE + +# XGBoost hyperparameters, frozen since Entry 2 (depth 4->3, added L1/L2 + +# min_child_weight, tightened subsample/colsample after LOSO ablation showed +# the Entry 1 config overfit to within-site patterns). Do not tune these +# without a LOSO regression check — see learning_log.md, 2026-07-01. +XGB_PARAMS = dict( + n_estimators=300, + max_depth=3, + learning_rate=0.05, + subsample=0.7, + colsample_bytree=0.6, + reg_alpha=0.1, + reg_lambda=2.0, + min_child_weight=5, + random_state=42, + eval_metric='auc', + verbosity=0, +) + + +def _scale_pos_weight(y_train): + y_train = np.asarray(y_train) + n_pos = int((y_train == 1).sum()) + n_neg = int((y_train == 0).sum()) + return n_neg / n_pos if n_pos > 0 else 1.0 + + +def build_pipeline(y_train, calibrated=True, use_age_residual=True, + calibration_ensemble=True, xgb_overrides=None): + """ + Builds the full Narnia_ML model pipeline: + + [AgeResidualizer ->] SimpleImputer(median) -> XGBoost [+ CalibratedClassifierCV] + + y_train is needed up front (not deferred to .fit() time) only to + compute scale_pos_weight, matching every prior entry's convention of + deriving it from the actual training split rather than hardcoding it. + + calibrated=True (default — matches team_code.py / Entry 3+): wraps + XGBoost in CalibratedClassifierCV(method='sigmoid', cv=5). + calibrated=False: uncalibrated single XGBoost. Used by loso_cv.py's + run_ablation() and any fast/uncalibrated fold metrics, where the + extra ~5x training cost of calibration buys nothing. + + calibration_ensemble=True (default — matches shipped team_code.py): + CalibratedClassifierCV's own `ensemble=True` behavior — 5 models + trained on different folds, their 5 calibrated outputs averaged + at predict time. + calibration_ensemble=False (2026-07-06): passes `ensemble=False` to + CalibratedClassifierCV. RULED OUT as a fix for the large-set + reward collapse (2026-07-06 entry — identical raw scores, no + measurable reward/AUROC change vs ensemble=True). Kept available + for completeness, not because it's still a live hypothesis. + + use_age_residual=True (default — matches team_code.py / Entry 4): + includes the AgeResidualizer step (2 extra features, 48 -> 50). + use_age_residual=False (2026-07-06): builds the Entry-3-EQUIVALENT + pipeline. Also largely ruled out as a driver of the large-set + reward collapse (2026-07-06 attribution entry: negligible effect + at both scales) — kept for completeness/future feature work, not + because it's still the leading hypothesis for that specific problem. + + xgb_overrides=None (default): XGBoost uses XGB_PARAMS unchanged, the + config validated via LOSO at Entry 2 (small-scale n). Pass a dict + to override specific keys, e.g. {'reg_lambda': 10.0, + 'min_child_weight': 25} — added 2026-07-06 for loso_cv.py's + --tune-hyperparams sweep, testing whether XGB_PARAMS (never + re-validated above ~900-patient folds) is under-regularized at + large-scale fold sizes (~1,000-5,500 patients). See + learning_log.md, 2026-07-06 (raw score inflation entry) for the + finding motivating this: the I0006 holdout fold's raw + (pre-calibration) XGBoost output shifted from median 0.12 (small + scale) to median 0.69 (large scale) — a base-model phenomenon, + confirmed independent of calibration architecture. team_code.py + never passes this — it always wants the validated XGB_PARAMS — + this exists for loso_cv.py's comparison runs only. + + Returns an UNFITTED sklearn Pipeline. Caller is responsible for + calling .fit(X_train, y_train) — this function never fits anything, + so the same discipline (fit on train, freeze for inference/LOSO-test) + that AgeResidualizer and SimpleImputer already enforce internally is + also enforced at the call-site level: nothing in this function can + accidentally see test data. + """ + spw = _scale_pos_weight(y_train) + params = {**XGB_PARAMS, **(xgb_overrides or {})} + xgb = XGBClassifier(scale_pos_weight=spw, **params) + + classifier = ( + CalibratedClassifierCV(xgb, method='sigmoid', cv=5, ensemble=calibration_ensemble) + if calibrated else xgb + ) + + steps = [] + if use_age_residual: + steps.append(('age_residual', AgeResidualizer( + age_idx=IDX_AGE, + ca_rate_idx=IDX_CA_RATE, + eeg_var_rem_wake_idx=IDX_EEG_VAR_REM_WAKE, + ))) + steps.append(('imputer', SimpleImputer(strategy='median'))) + steps.append(('classifier', classifier)) + + return Pipeline(steps) + + +# Logistic regression hyperparameters. class_weight='balanced' stands in +# for XGB_PARAMS' scale_pos_weight -- same purpose (correct for the ~7% +# positive rate), different mechanism; LogisticRegression has no +# scale_pos_weight argument. C is passed separately per-call (see +# build_logreg_pipeline) since it's the one parameter under active +# investigation, not frozen the way XGB_PARAMS is. +LOGREG_PARAMS = dict( + penalty='l2', + class_weight='balanced', + max_iter=2000, + random_state=42, +) + + +def build_logreg_pipeline(y_train, calibrated=True, use_age_residual=True, + calibration_ensemble=True, C=0.01, + penalty=None, l1_ratio=None, max_iter=None): + """ + Logistic regression alternative to build_pipeline()'s XGBoost model. + Added 2026-07-06 after the model-family diagnostic (loso_cv.py + --diagnostic-models) showed logistic regression beating XGBoost's + age-conditioned AUROC on BOTH large-set LOSO folds: + + Model I0006 AUROC S0001 AUROC Mean + xgboost_baseline 0.549 0.539 0.544 + logreg_strong_l2 0.615 0.585 0.600 (C=0.01) + + Mean 0.600 is the first result in the whole large-set investigation + to clear the "real, not noise" bar established from the leaderboard's + own noise-floor analysis (~north of 0.57-0.58). It also showed LESS + severe raw-score inflation on the I0006 fold specifically (median + 0.47 vs XGBoost's 0.69) — see learning_log.md, 2026-07-06 (model + diagnostic entry) for the full comparison. The hyperparameter sweep + that ran alongside this diagnostic ruled out "XGBoost just needs more + regularization" — every more-conservative XGBoost candidate made the + I0006 inflation WORSE, not better — so this isn't a same-model tweak, + it's a different model family that empirically behaves better on the + exact problem under investigation. + + Same [AgeResidualizer ->] Imputer -> ... -> Classifier [+ Calibration] + shape as build_pipeline(), with two necessary differences: + - Adds a StandardScaler step after the imputer. Unlike tree-based + XGBoost, a linear model is sensitive to the wildly different + scales across this feature vector (e.g. Age in years vs. Hjorth + ratios near 1.0) and won't behave sensibly without it. + - class_weight='balanced' (see LOGREG_PARAMS) stands in for + XGBoost's scale_pos_weight. + + C=0.01 (default) matches the diagnostic's best-performing candidate. + Exposed as a parameter for further tuning — not yet swept the way + XGB_PARAMS was; this is the first calibrated test of this model + family, not a converged config. + + penalty/l1_ratio/max_iter (added 2026-07-09, for the Ridge/Lasso/ + ElasticNet regularization sweep): NEW — LOGREG_PARAMS hardcodes + penalty='l2', so l1/elasticnet were never tried before this. Left as + None by default so an un-parameterized call reproduces the EXACT + original tested config (penalty='l2' from LOGREG_PARAMS, default + solver, max_iter=2000) byte-for-byte — passing any of these three + overrides that default. l1/elasticnet require solver='saga' + (LOGREG_PARAMS' implicit lbfgs-default doesn't support them), so the + solver is switched automatically only when a non-l2 penalty is + requested — the confirmed-working l2 path never changes solver. + + calibrated / use_age_residual / calibration_ensemble: same meaning + and defaults as build_pipeline(). Calibration mode was ruled out as + the driver of the ORIGINAL XGBoost reward problem, but that finding + doesn't automatically transfer to a different base model — worth + re-checking cv5 vs cv5-single-curve on THIS model too if it looks + promising, rather than assuming the earlier conclusion still applies. + + Returns an UNFITTED sklearn Pipeline. Caller is responsible for + calling .fit(X_train, y_train) — same contract as build_pipeline(). + team_code.py does not use this function; it remains on build_pipeline() + (XGBoost) unless/until a decision is logged in learning_log.md to + switch, per this project's standing discipline. + """ + steps = [] + if use_age_residual: + steps.append(('age_residual', AgeResidualizer( + age_idx=IDX_AGE, + ca_rate_idx=IDX_CA_RATE, + eeg_var_rem_wake_idx=IDX_EEG_VAR_REM_WAKE, + ))) + steps.append(('imputer', SimpleImputer(strategy='median'))) + steps.append(('scaler', StandardScaler())) + + lr_params = dict(LOGREG_PARAMS) + if max_iter is not None: + lr_params['max_iter'] = max_iter + if penalty is not None: + lr_params['penalty'] = penalty + if lr_params['penalty'] in ('l1', 'elasticnet'): + lr_params['solver'] = 'saga' # LOGREG_PARAMS' l2 default (lbfgs) can't do l1/elasticnet + if lr_params['penalty'] == 'elasticnet': + if l1_ratio is None: + raise ValueError("penalty='elasticnet' requires l1_ratio in [0, 1].") + lr_params['l1_ratio'] = l1_ratio + + logreg = LogisticRegression(C=C, **lr_params) + classifier = ( + CalibratedClassifierCV(logreg, method='sigmoid', cv=5, ensemble=calibration_ensemble) + if calibrated else logreg + ) + steps.append(('classifier', classifier)) + + return Pipeline(steps) + + +class PrefitCalibratedModel: + """ + Alternative to build_pipeline(calibrated=True)'s CalibratedClassifierCV + (..., cv=5) ensembling. Added 2026-07-06 after two independent pieces + of evidence that cv=5's per-fold calibration curves may not transfer + stably across training-set scale: + + 1. Real leaderboard, Entry 2 -> Entry 3 (calibration added, nothing + else changed): standard AUROC rose (0.746 -> 0.772) while + age-conditioned AUROC and reward both fell, and accuracy dropped + sharply (0.866 -> 0.803) with F-measure barely moving -- more + patients being called positive, at a cost. + 2. Large-set LOSO (Entry 4 config): reward at the identical + THRESHOLD=0.12 dropped from 0.1148 (small set) to 0.0726 (large + set), traced to pooled specificity collapsing from ~93% to ~68% + -- age-conditioned AUROC barely moved, so this is a calibrated- + PROBABILITY-magnitude problem, not a ranking problem. + + Mechanism under test: cv=5 (without cv='prefit') fits 5 separate + XGBoost models on different 80% training slices and averages their + independently-fit Platt curves. Each of those 5 curves is itself + sensitive to whatever slice it happened to see, on top of the whole + ensemble's sensitivity to the overall training set's size/composition + (911 patients at small scale vs several thousand at large scale, per + LOSO fold). THRESHOLD=0.12 is applied AFTER this mapping, so if the + mapping itself shifts between training runs, the same fixed number + stops meaning the same thing -- this is the leading explanation for + both observations above. + + This class fits ONE XGBoost model on a training split, then fits ONE + Platt curve (cv='prefit') on a SEPARATE held-out calibration split -- + removing the 5-curves-averaged-together source of instability + entirely, so a comparison against build_pipeline(calibrated=True) + isolates whether cv=5's ensembling specifically is the cause. + + NOT a sklearn Pipeline: cv='prefit' calibration inherently needs two + different data slices (one to fit the base estimator, a separate one + to fit the calibrator), which doesn't fit Pipeline's contract of one + X flowing through every step via a single .fit(X, y) call. Exposes + .fit(X, y) / .predict_proba(X) with the same signatures a fitted + Pipeline would, so it's a drop-in replacement at loso_cv.py's call + sites without those call sites needing to know which mode is active. + + Caveat: splitting off a separate calibration slice costs real training + data, which matters more at small scale. At the large training set's + per-fold sizes (~1,000-5,500), a calib_fraction=0.2 split still leaves + a healthy number of positives on both sides. At the SMALL set's sizes + (~250-900 per fold, 28-64 positives), the same split leaves considerably + less on each side -- if comparing this against cv=5 on the small set + specifically, read the result with that in mind rather than treating + it as a clean apples-to-apples comparison of calibration mode alone. + """ + + def __init__(self, use_age_residual=True, calib_fraction=0.2, random_state=42): + self.use_age_residual = use_age_residual + self.calib_fraction = calib_fraction + self.random_state = random_state + + def fit(self, X, y): + X = np.asarray(X, dtype=float) + y = np.asarray(y) + + idx_train, idx_calib = train_test_split( + np.arange(len(y)), test_size=self.calib_fraction, + stratify=y, random_state=self.random_state) + + X_train_raw, y_train = X[idx_train], y[idx_train] + X_calib_raw, y_calib = X[idx_calib], y[idx_calib] + + # Preprocessing (age residual + imputer) fit ONLY on the training + # slice -- same train-only-fit discipline as everywhere else in + # this pipeline. The calibration slice only ever gets .transform()'d. + steps = [] + if self.use_age_residual: + steps.append(('age_residual', AgeResidualizer( + age_idx=IDX_AGE, + ca_rate_idx=IDX_CA_RATE, + eeg_var_rem_wake_idx=IDX_EEG_VAR_REM_WAKE, + ))) + steps.append(('imputer', SimpleImputer(strategy='median'))) + self.preprocessing_ = Pipeline(steps) + + X_train = self.preprocessing_.fit_transform(X_train_raw) + X_calib = self.preprocessing_.transform(X_calib_raw) + + spw = _scale_pos_weight(y_train) + xgb = XGBClassifier(scale_pos_weight=spw, **XGB_PARAMS) + xgb.fit(X_train, y_train) + self.base_estimator_ = xgb # kept for feature_importances_ access + + if _HAS_FROZEN_ESTIMATOR: + self.calibrated_ = CalibratedClassifierCV(FrozenEstimator(xgb), method='sigmoid') + else: + self.calibrated_ = CalibratedClassifierCV(xgb, method='sigmoid', cv='prefit') + self.calibrated_.fit(X_calib, y_calib) + return self + + def predict_proba(self, X): + X_t = self.preprocessing_.transform(np.asarray(X, dtype=float)) + return self.calibrated_.predict_proba(X_t) + + @property + def feature_importances_(self): + return self.base_estimator_.feature_importances_ + + +def extract_fitted_coefficients(fitted_pipeline): + """ + ADDED 2026-07-09, not part of the original file — pulls linear + coefficients back out of a fitted build_logreg_pipeline() Pipeline, + for the coefficient-stability-across-folds check. Handles both + calibrated=True (CalibratedClassifierCV wraps N cloned+refit estimators + — averages their coefficients) and calibrated=False. + + Returns a 1D np.array in STANDARDIZED-feature space (coefficients on + scaled features) — comparable across folds/features for a stability + check, but not directly interpretable as raw-unit effect sizes. + """ + clf = fitted_pipeline.named_steps['classifier'] + if isinstance(clf, CalibratedClassifierCV): + coefs = np.stack([cc.estimator.coef_[0] for cc in clf.calibrated_classifiers_]) + return coefs.mean(axis=0) + return clf.coef_[0] \ No newline at end of file diff --git a/features/sample_weighting.py b/features/sample_weighting.py new file mode 100644 index 0000000..cf4d124 --- /dev/null +++ b/features/sample_weighting.py @@ -0,0 +1,86 @@ +# features/sample_weighting.py +# Team Narnia — PhysioNet Challenge 2026 +# +# Entry 6 (2026-07-16/19): compute_value_weighted_sample_weights(), promoted +# from tools/reg_sweep.py (where it was developed and validated on Kaggle — +# see learning_log.md / learning_log_3, 2026-07-16 entries) into a shared +# features/ module, same precedent as AgeResidualizer (features/age_residuals.py): +# one definition, imported by both team_code.py (the real submission path) +# and tools/reg_sweep.py (the validation harness it was proven in), so the +# two can never silently diverge the way loso_cv.py's hand-copied Pipeline +# once did (see features/pipeline.py header for that incident). +# +# Validated result this function is responsible for (alpha=1.0, on top of +# C=0.001): +12.9% relative reward vs. C=0.001 alone, AUROC flat (max 0.21σ +# across all 3 folds). Reproduced bit-for-bit on an independent Kaggle run +# (learning_log_3, 2026-07-16 final entry). Combined with the C=0.001 change: +# +52% relative reward vs. the original shipped Entry 5 baseline, AUROC flat +# throughout the entire chain. + +import numpy as np + +from evaluate_model import compute_prevalence + + +def _lookup_prevalence_train(age_to_prevalence, age): + """Same nearest-key fallback logic as tools/test_age_banded_threshold.py's + _lookup_prevalence — kept as a local duplicate rather than a cross-file + import, matching reg_sweep.py's original standalone-on-Kaggle discipline.""" + key = round(age) + if key in age_to_prevalence: + return age_to_prevalence[key] + nearest = min(age_to_prevalence.keys(), key=lambda k: abs(k - age)) + return age_to_prevalence[nearest] + + +def compute_value_weighted_sample_weights(y_train, age_train, alpha): + """ + Added 2026-07-16 (tools/reg_sweep.py), promoted here 2026-07-19 for + Entry 6. Builds a per-training-sample weight array that layers a + reward-VALUE-aware boost on top of whatever class_weight='balanced' + already does (LOGREG_PARAMS — confirmed already active in every logreg + config tested this project; this is an ADDITION, not a replacement). + sklearn composes class_weight and an explicit sample_weight array + multiplicatively, so both apply together. + + ONLY positive-class samples get boosted — negatives keep weight 1.0 + regardless of age, matching the actual intent (prioritize getting + high-value positives right, not reweighting the whole population by + age indiscriminately). + + CRITICAL LEAKAGE DISCIPLINE: age_to_prevalence here is fit using ONLY + y_train/age_train. At LOSO-validation time (reg_sweep.py) that means + only the patients in THIS fold's training set, never the held-out test + fold. At real submission time (team_code.py) there is no fold rotation + at all — train_model() trains once on whatever training set the + organizers provide, so passing that full training set here is itself + the correct, leakage-safe usage; there is no held-out slice to leak + from in production. This is deliberately different from how + compute_prevalence is used elsewhere in this project: reward SCORING + correctly uses the FULL population as its reference (confirmed + 2026-07-09 finding), because that's evaluating an already-fixed + decision against reality. This is different — it directly shapes what + the MODEL LEARNS from these exact training patients, so using anything + beyond the training set's own data here would leak test-correlated + information into training, the same class of mistake AgeResidualizer's + train-fold-only fitting discipline exists to prevent. + + alpha=0.0 reproduces IDENTICAL behavior to no weighting at all (returns + all-ones) — backward compatible, matches every previously logged result + exactly when alpha isn't explicitly set. + """ + if alpha == 0.0: + return np.ones(len(y_train)) + + age_to_prevalence_train = compute_prevalence(age_train, y_train, age_train, gap=2) + prevalences = np.array([_lookup_prevalence_train(age_to_prevalence_train, a) + for a in age_train]) + value = (1.0 / prevalences) - 1.0 + # Normalize using this training set's own min/max — not any global + # reference — same train-only discipline as the prevalence fit above. + value_norm = (value - value.min()) / (value.max() - value.min() + 1e-12) + + weights = np.ones(len(y_train)) + pos_mask = (y_train == 1) + weights[pos_mask] = 1.0 + alpha * value_norm[pos_mask] + return weights \ No newline at end of file diff --git a/features/spectral.py b/features/spectral.py new file mode 100644 index 0000000..62b486c --- /dev/null +++ b/features/spectral.py @@ -0,0 +1,522 @@ +# features/spectral.py +# Raw-EDF spectral features — gated feature category (FEATURES.md, +# "Week 5+: Raw EDF features", gate condition: CAISR-only AUROC < 0.62 +# after entry 3-4). Entry 5 landed at 0.6002, treated as satisfied-in- +# spirit given deadline pressure (see TEST_PLAN.md, 2026-07-11 K3 note). +# +# Unlike caisr_base/caisr_enriched (algorithmic annotation summaries) and +# physiological_ratios (Hjorth-derived time-domain features), this module +# computes actual frequency-domain content from the raw EEG signal. This +# is the first feature category that touches real signal content rather +# than CAISR-derived or time-domain-derived summary statistics. +# +# Priority order confirmed by the 2026-07-11 MCI/age cross-tab +# (learning_log.md): MCI is 44% of positives, significantly younger than +# other subtypes (68.4 vs 71.8 yr, p=0.00001), and has the lowest +# sensitivity (0.632 vs 0.766) — landing disproportionately in the +# highest-reward-value age band. This points toward early/subtle-stage +# slow-wave markers (this module) over autonomic/HRV markers (hrv.py, +# not yet built), since MCI is the earlier, subtler-signal stage of +# decline. +# +# Iteration history: +# v1 (this file, first feature only): N3 delta power. Deliberately +# scoped to ONE feature first, per the team's one-variable-per-test +# discipline (learning_log.md decision log) and to get a real LOSO +# read before committing to the full gated list (delta/theta ratio, +# spindle density — both planned, not yet implemented, see TODOs below). +# +# NOT YET WIRED IN: this module is standalone. Adding its output to the +# actual feature vector requires updating features/__init__.py (new +# IDX_SPECTRAL_START/END constants, N_TOTAL_FEATURES bump) and every +# hstack call site (team_code.py, tools/loso_cv.py, +# tools/build_features_cache.py) consistently. That's a deliberate +# follow-up step, not done here — do not wire this in without updating +# all three call sites together, per the same discipline that made +# build_features_cache.py's verification necessary in the first place. + +import numpy as np +import os +from scipy.signal import welch +from helper_code import ( + load_rename_rules, standardize_channel_names_rename_only, + derive_bipolar_signal +) + +_FEATURES_DIR = os.path.dirname(os.path.abspath(__file__)) +_REPO_DIR = os.path.dirname(_FEATURES_DIR) +DEFAULT_CSV_PATH = os.path.join(_REPO_DIR, 'channel_table.csv') + +N_SPECTRAL_FEATURES = 1 # N3 delta power only, for now. Planned: +2 + # (delta/theta ratio, N2 spindle density) once + # this feature has a real LOSO read. + +# Same threshold as physiological_ratios.py — insufficient stage duration +# means the estimate isn't trustworthy, return NaN rather than a noisy +# number. +_MIN_STAGE_SECONDS = 120 # 2 minutes + +# Standard delta band definition (AASM convention). +_DELTA_BAND = (0.5, 4.0) # Hz + +# Welch PSD window length. 4s windows give ~0.25Hz frequency resolution, +# fine enough to resolve the delta band's upper edge (4Hz) cleanly, short +# enough that a single N3 segment (often only a few minutes per patient) +# still yields several averaged windows rather than one noisy periodogram. +_WELCH_WINDOW_SECONDS = 4.0 + + +def _get_channel(channels, candidates): + """Return first available channel from candidates list.""" + for c in candidates: + if c in channels and channels[c] is not None: + return channels[c] + return None + + +def _get_fs(fs_map, candidates, default=200.0): + for c in candidates: + if c in fs_map: + return float(fs_map[c]) + return default + + +def _bandpower(sig, fs, band, window_seconds=_WELCH_WINDOW_SECONDS): + """ + Welch PSD bandpower for a single 1-D signal segment. + + Returns NaN if the segment is too short for even one Welch window, + or if the signal is degenerate (near-zero variance — e.g. a flat/ + disconnected channel). + """ + if sig is None or len(sig) < 10: + return float('nan') + + nperseg = int(round(window_seconds * fs)) + if nperseg < 8 or len(sig) < nperseg: + return float('nan') + + if float(np.var(sig)) < 1e-20: + return 0.0 + + freqs, psd = welch(sig, fs=fs, nperseg=nperseg) + + band_mask = (freqs >= band[0]) & (freqs <= band[1]) + if not np.any(band_mask): + return float('nan') + + # np.trapz was removed in NumPy 2.0+ (renamed np.trapezoid). Kaggle's + # image and local dev environments may not be on the same NumPy + # version, so don't hardcode either name — this bit a local sanity + # test during development (NumPy 2.4.4 has no np.trapz at all). + _trapz_fn = getattr(np, 'trapezoid', None) or np.trapz + return float(_trapz_fn(psd[band_mask], freqs[band_mask])) + + +def extract_spectral_features(phys_data, phys_fs, algo_data, + csv_path=DEFAULT_CSV_PATH): + """ + Raw-EDF spectral features requiring both physiological EDF (raw + signal) and CAISR annotations (stage timing) — same dual-source + requirement as physiological_ratios.py. + + Returns np.ndarray of length N_SPECTRAL_FEATURES (currently 1): + [0] EEG delta-band (0.5-4Hz) power during N3 sleep only. + Absolute power, NOT a ratio — unlike physiological_ratios.py's + features, this is not expected to be site-stable by + construction. Equipment scale could plausibly reintroduce + the same confound that made absolute Hjorth features fail + cross-site LOSO (2026-07-01 finding). This needs its OWN + LOSO ablation before being trusted as a keeper, not an + assumption that the delta-power literature backing + transfers automatically to this equipment-confound-prone + measurement style. Flagging explicitly so this isn't + silently assumed safe. + + NaN fallback: + - Missing physio EDF → NaN + - Missing CAISR → NaN + - Insufficient N3 duration (<120s) → NaN + - Signal too short for even one Welch window → NaN + """ + if not phys_data or not algo_data: + return np.full(N_SPECTRAL_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) < 10: + return np.full(N_SPECTRAL_FEATURES, float('nan')) + + # Same channel-standardization pattern as physiological_ratios.py — + # reused deliberately, not reimplemented, so both modules stay + # consistent if channel_table.csv or rename rules ever change. + original_labels = list(phys_data.keys()) + rename_rules = load_rename_rules(os.path.abspath(csv_path)) + rename_map, cols_to_drop = standardize_channel_names_rename_only( + original_labels, rename_rules) + + channels = {} + fs_map = {} + for old_label, data in phys_data.items(): + if old_label in cols_to_drop: + continue + new_label = rename_map.get(old_label, old_label.lower()) + channels[new_label] = data + if old_label in phys_fs: + fs_map[new_label] = phys_fs[old_label] + + # Same bipolar derivation candidates as physiological_ratios.py. + for target, pos, neg_list in [ + ('f3-m2', 'f3', ['m2']), + ('f4-m1', 'f4', ['m1']), + ('c3-m2', 'c3', ['m2']), + ('c4-m1', 'c4', ['m1']), + ]: + if target in channels or pos not in channels: + continue + if not all(n in channels for n in neg_list): + continue + ref = channels[neg_list[0]] + derived = derive_bipolar_signal(channels[pos], ref) + if derived is not None: + channels[target] = derived + fs_map[target] = fs_map.get(pos, 200.0) + + result = np.full(N_SPECTRAL_FEATURES, float('nan')) + + eeg_sig = _get_channel(channels, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + eeg_fs = _get_fs(fs_map, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + + if eeg_sig is None or len(eeg_sig) == 0: + return result + + spe = int(round(30.0 * eeg_fs)) # samples per 30s epoch + min_smp = int(_MIN_STAGE_SECONDS * eeg_fs) + + # Same upsample-stage-to-signal-rate pattern as physiological_ratios.py. + stage_eeg = np.repeat(valid_stages, spe)[:len(eeg_sig)] + n3_mask = (stage_eeg == 1) # N3 stage code — confirmed via + # caisr_base.py's stage convention + # (W=5, N1=3, N2=2, N3=1, REM=4) + + if int(n3_mask.sum()) < min_smp: + return result + + n3_segment = eeg_sig[:len(stage_eeg)][n3_mask] + result[0] = _bandpower(n3_segment, eeg_fs, _DELTA_BAND) + + return result + + +# Standard theta band definition (AASM convention: 4-8Hz). Some clinical +# references use 4-7Hz — 4-8Hz chosen here for a clean non-overlapping +# boundary with the delta band above (0.5-4Hz) and the alpha band below +# (typically 8-13Hz). Worth revisiting if delta/theta ratio results look +# sensitive to this exact boundary. +_THETA_BAND = (4.0, 8.0) # Hz + +N_DELTA_THETA_FEATURES = 1 +N_SPINDLE_FEATURES = 1 + +# Sigma band — standard sleep spindle frequency range (AASM convention). +_SIGMA_BAND = (11.0, 16.0) # Hz + +# Spindle duration constraints (AASM: spindles are 0.5-3s bursts). +_MIN_SPINDLE_DURATION_S = 0.5 +_MAX_SPINDLE_DURATION_S = 3.0 + +# Envelope smoothing window — short enough to preserve individual +# spindle onsets/offsets, long enough to avoid triggering on single- +# sample envelope noise. +_ENVELOPE_SMOOTH_S = 0.1 + +# Adaptive threshold: mean + K * std of the envelope, computed from the +# SAME N2 segment being scored. This is deliberately per-recording, not +# a fixed absolute value — if equipment scale multiplies the signal (and +# therefore the envelope) by some gain k, both the mean and std scale by +# k too, so "envelope > mean + K*std" is unaffected by k. Same site- +# stability reasoning as the delta/theta ratio, applied to event +# detection instead of a power ratio. +_THRESHOLD_STD_MULTIPLIER = 1.5 + +# Same cap convention as physiological_ratios.py — prevents a +# near-zero-theta patient from producing an extreme outlier ratio. +_RATIO_CAP = 20.0 + + +def _safe_ratio(num, denom, cap=_RATIO_CAP): + """Scalar-safe ratio with NaN/near-zero guards. Same discipline as + physiological_ratios.py's _safe_ratio, just for a single value + instead of a feature vector (this module only produces one ratio + per patient so far).""" + if num is None or denom is None or np.isnan(num) or np.isnan(denom): + return float('nan') + if abs(denom) < 1e-12: + return float('nan') + return float(np.clip(num / denom, 0.0, cap)) + + +def extract_delta_theta_ratio_features(phys_data, phys_fs, algo_data, + csv_path=DEFAULT_CSV_PATH): + """ + N3 delta/theta bandpower ratio — kept as a SEPARATE function from + extract_spectral_features (raw delta power) deliberately, so the two + can be ablation-tested independently. Raw delta power's small-set + ablation (learning_log.md, 2026-07-12) came back null-to-mildly- + negative on the best-powered fold (S0001) — this ratio version is a + distinct hypothesis, not an extension of that result, and should be + judged on its own LOSO read. + + Rationale for expecting better cross-site behavior than raw delta + power: if equipment introduces a multiplicative gain confound on the + EEG channel, it inflates delta and theta power equally, so it + cancels in the ratio — same mechanism that made the existing + within-recording ratio features (physiological_ratios.py) survive + cross-site LOSO where the original absolute Hjorth features didn't + (2026-07-01 finding). This is a hypothesis carried over by analogy, + not yet confirmed for THIS specific ratio — needs its own ablation + read, same as every other feature in this project. + + Returns np.ndarray of length N_DELTA_THETA_FEATURES (1): + [0] N3 delta-band power / N3 theta-band power. Both bands + computed from the SAME N3-masked segment via the SAME Welch + PSD call (not two separate extractions) — cheaper and + guarantees identical windowing between numerator and + denominator. + + NaN fallback: identical conditions to extract_spectral_features + (missing physio/CAISR, insufficient N3 duration, degenerate signal). + """ + if not phys_data or not algo_data: + return np.full(N_DELTA_THETA_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) < 10: + return np.full(N_DELTA_THETA_FEATURES, float('nan')) + + original_labels = list(phys_data.keys()) + rename_rules = load_rename_rules(os.path.abspath(csv_path)) + rename_map, cols_to_drop = standardize_channel_names_rename_only( + original_labels, rename_rules) + + channels = {} + fs_map = {} + for old_label, data in phys_data.items(): + if old_label in cols_to_drop: + continue + new_label = rename_map.get(old_label, old_label.lower()) + channels[new_label] = data + if old_label in phys_fs: + fs_map[new_label] = phys_fs[old_label] + + for target, pos, neg_list in [ + ('f3-m2', 'f3', ['m2']), + ('f4-m1', 'f4', ['m1']), + ('c3-m2', 'c3', ['m2']), + ('c4-m1', 'c4', ['m1']), + ]: + if target in channels or pos not in channels: + continue + if not all(n in channels for n in neg_list): + continue + ref = channels[neg_list[0]] + derived = derive_bipolar_signal(channels[pos], ref) + if derived is not None: + channels[target] = derived + fs_map[target] = fs_map.get(pos, 200.0) + + result = np.full(N_DELTA_THETA_FEATURES, float('nan')) + + eeg_sig = _get_channel(channels, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + eeg_fs = _get_fs(fs_map, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + + if eeg_sig is None or len(eeg_sig) == 0: + return result + + spe = int(round(30.0 * eeg_fs)) + min_smp = int(_MIN_STAGE_SECONDS * eeg_fs) + + stage_eeg = np.repeat(valid_stages, spe)[:len(eeg_sig)] + n3_mask = (stage_eeg == 1) + + if int(n3_mask.sum()) < min_smp: + return result + + n3_segment = eeg_sig[:len(stage_eeg)][n3_mask] + + delta_power = _bandpower(n3_segment, eeg_fs, _DELTA_BAND) + theta_power = _bandpower(n3_segment, eeg_fs, _THETA_BAND) + + result[0] = _safe_ratio(delta_power, theta_power) + + return result + + +def _detect_spindles(sig, fs, min_duration_s=_MIN_SPINDLE_DURATION_S, + max_duration_s=_MAX_SPINDLE_DURATION_S, + threshold_k=_THRESHOLD_STD_MULTIPLIER): + """ + Simplified amplitude-envelope spindle detector: + 1. Bandpass filter to sigma band (11-16Hz) + 2. Hilbert transform -> amplitude envelope + 3. Smooth envelope (short moving average) + 4. Threshold at mean + K*std of THIS segment's own envelope + (adaptive, not absolute — see module-level comment on why) + 5. Count contiguous above-threshold runs whose duration falls + within [min_duration_s, max_duration_s] + + This is a simplified single-channel amplitude detector, not a + validated clinical-grade algorithm (e.g. it does not implement + multi-channel consensus, doesn't model the characteristic waxing/ + waning spindle envelope shape beyond a simple threshold, and doesn't + merge near-adjacent events). Documented explicitly as a limitation — + do not present spindle counts from this function as clinically + validated without further work. It IS a legitimate first-pass + implementation of the standard "filter -> envelope -> threshold -> + duration-gate" approach used by simpler published detectors. + + Returns the count of qualifying events (int), or NaN if the signal + is too short to filter/analyze. + """ + from scipy.signal import butter, filtfilt, hilbert + + if sig is None or len(sig) < int(fs * 2): # need at least ~2s for a stable filter + return float('nan') + + nyq = fs / 2.0 + low = _SIGMA_BAND[0] / nyq + high = min(_SIGMA_BAND[1] / nyq, 0.99) # guard against fs too low + if low <= 0 or high <= low: + return float('nan') + + try: + b, a = butter(4, [low, high], btype='band') + filtered = filtfilt(b, a, sig) + except Exception: + return float('nan') + + envelope = np.abs(hilbert(filtered)) + + smooth_win = max(1, int(round(_ENVELOPE_SMOOTH_S * fs))) + if smooth_win > 1: + kernel = np.ones(smooth_win) / smooth_win + envelope = np.convolve(envelope, kernel, mode='same') + + env_mean = float(np.mean(envelope)) + env_std = float(np.std(envelope)) + if env_std < 1e-20: + return 0 # flat/degenerate signal — no spindles detectable, not an error + + threshold = env_mean + threshold_k * env_std + above = envelope > threshold + + # Count contiguous above-threshold runs meeting the duration gate. + min_samples = int(round(min_duration_s * fs)) + max_samples = int(round(max_duration_s * fs)) + + count = 0 + run_length = 0 + for val in above: + if val: + run_length += 1 + else: + if min_samples <= run_length <= max_samples: + count += 1 + run_length = 0 + if min_samples <= run_length <= max_samples: # trailing run at signal end + count += 1 + + return count + + +def extract_spindle_density_features(phys_data, phys_fs, algo_data, + csv_path=DEFAULT_CSV_PATH): + """ + N2 sleep spindle density — spindles per minute of N2 sleep. + + Unlike extract_spectral_features (bandpower) and + extract_delta_theta_ratio_features (power ratio), this is an + EVENT-RATE measure: it counts discrete oscillatory bursts via + _detect_spindles() rather than integrating power in a band. See that + function's docstring for the detection method and its limitations. + + Returns np.ndarray of length N_SPINDLE_FEATURES (1): + [0] Spindle count during N2 / N2 duration in minutes. + + NaN fallback: missing physio/CAISR, insufficient N2 duration + (<120s, same _MIN_STAGE_SECONDS threshold as every other feature in + this module and physiological_ratios.py), or signal too short to + filter. + """ + if not phys_data or not algo_data: + return np.full(N_SPINDLE_FEATURES, float('nan')) + + stages_raw = algo_data.get('stage_caisr', np.array([])) + valid_stages = stages_raw[stages_raw < 9.0] if len(stages_raw) > 0 else np.array([]) + + if len(valid_stages) < 10: + return np.full(N_SPINDLE_FEATURES, float('nan')) + + original_labels = list(phys_data.keys()) + rename_rules = load_rename_rules(os.path.abspath(csv_path)) + rename_map, cols_to_drop = standardize_channel_names_rename_only( + original_labels, rename_rules) + + channels = {} + fs_map = {} + for old_label, data in phys_data.items(): + if old_label in cols_to_drop: + continue + new_label = rename_map.get(old_label, old_label.lower()) + channels[new_label] = data + if old_label in phys_fs: + fs_map[new_label] = phys_fs[old_label] + + for target, pos, neg_list in [ + ('f3-m2', 'f3', ['m2']), + ('f4-m1', 'f4', ['m1']), + ('c3-m2', 'c3', ['m2']), + ('c4-m1', 'c4', ['m1']), + ]: + if target in channels or pos not in channels: + continue + if not all(n in channels for n in neg_list): + continue + ref = channels[neg_list[0]] + derived = derive_bipolar_signal(channels[pos], ref) + if derived is not None: + channels[target] = derived + fs_map[target] = fs_map.get(pos, 200.0) + + result = np.full(N_SPINDLE_FEATURES, float('nan')) + + eeg_sig = _get_channel(channels, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + eeg_fs = _get_fs(fs_map, ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1']) + + if eeg_sig is None or len(eeg_sig) == 0: + return result + + spe = int(round(30.0 * eeg_fs)) + min_smp = int(_MIN_STAGE_SECONDS * eeg_fs) + + stage_eeg = np.repeat(valid_stages, spe)[:len(eeg_sig)] + n2_mask = (stage_eeg == 2) # N2 stage code — confirmed via + # caisr_base.py (W=5, N1=3, N2=2, + # N3=1, REM=4) + + if int(n2_mask.sum()) < min_smp: + return result + + n2_segment = eeg_sig[:len(stage_eeg)][n2_mask] + n2_minutes = len(n2_segment) / eeg_fs / 60.0 + + count = _detect_spindles(n2_segment, eeg_fs) + if isinstance(count, float) and np.isnan(count): + return result + + result[0] = float(count) / n2_minutes if n2_minutes > 0 else float('nan') + + return result \ No newline at end of file diff --git a/features/subtype_weighting.py b/features/subtype_weighting.py new file mode 100644 index 0000000..7d8a83c --- /dev/null +++ b/features/subtype_weighting.py @@ -0,0 +1,175 @@ +# features/subtype_weighting.py +# Team Narnia — PhysioNet Challenge 2026 +# +# Entry 7: MCI-targeted sample weighting, layered on top of the Entry 6 +# age-based value weighting (features/sample_weighting.py). Promoted from +# eda/build_positive_patient_report.py's CODE_FAMILY/classify_family/ +# dominant_family (originally interpretation-only, never touched by +# team_code.py) into a shared features/ module, because this is now a +# TRAINING-TIME input to the real submission path — same "single source +# of truth, not duplicated" precedent as AgeResidualizer and +# compute_value_weighted_sample_weights. +# +# CODE_FAMILY here includes 3 codes added 2026-07-19, after auditing every +# unique code in the real ICD_codes_CI.csv against the official +# ICD_codes.csv reference table (moody-challenge.physionet.org/2026/ +# ICD_codes.csv): G31.85 (Corticobasal degeneration), 331.11 (Pick's +# disease, ICD-9), F10.27 (Alcohol-induced persisting dementia) — all +# three were present in the real data but missing from the original +# eda-only dict. Verified zero practical impact on any 2026-07-19 finding +# (0 of 498 patients had their dominant_family() classification changed +# by this gap — every affected row co-occurred with an already-correctly- +# classified code for the same patient), but added here for completeness +# given this is now a training-time dependency, not just interpretation. +# +# Note: the official reference table only recognizes 3 top-level +# categories (MCI, Alzheimer's disease, Dementia) — the finer split used +# here (Lewy body / vascular / frontotemporal / unspecified / etc., all +# "Dementia" officially) is this project's own clinically-motivated +# breakdown, not the challenge's. Only the "Mild cognitive impairment" +# category is actually load-bearing for compute_mci_boosted_weights below +# — the finer dementia sub-buckets exist for interpretation/EDA +# elsewhere (site_subtype_breakdown.py), not for anything in this module. +# +# Validated candidate (2026-07-19 ablation, tools/ablation_mci_sample_ +# weight.py): beta_mci=0.75 on top of alpha_age=1.0 — reward +5.6% +# relative vs. alpha_age=1.0 alone, S0001 AUROC flat (+0.01 sigma), MCI +# sensitivity 76.8% -> 85.9%. A fine sweep (0.60-0.90) confirmed 0.75 as a +# genuine local optimum — a smooth, broad plateau, not a single-point +# fluke. See learning_log.md, 2026-07-19, for full methodology. +# +# LEAKAGE DISCIPLINE: ICD/subtype data is used ONLY to shape a training- +# time sample weight, exactly like the age-based weighting it extends — +# never a model input, never touched by run_model()/inference. Official +# ICD_codes_CI.csv is guaranteed present in both training_set_small and +# training_set_large (confirmed against the official challenge data- +# access documentation, moody-challenge.physionet.org/2026, 2026-07-19), +# at the same directory level as demographics.csv. + +import os +from collections import defaultdict + +import numpy as np +import pandas as pd + +from features.sample_weighting import compute_value_weighted_sample_weights + +ICD_CODES_FILENAME = 'ICD_codes_CI.csv' + +CODE_FAMILY = { + # ICD-10 + "G30": "Alzheimer disease", + "G31.84": "Mild cognitive impairment", + "G31.83": "Dementia with Lewy bodies", + "G31.0": "Frontotemporal dementia", + "G31.85": "Dementia in other disease", # added 2026-07-19 (Corticobasal degeneration) + "F01": "Vascular dementia", + "F02": "Dementia in other disease", + "F03": "Unspecified dementia", + "F10.27": "Dementia in other disease", # added 2026-07-19 (alcohol-induced persisting dementia) + # ICD-9 + "331.0": "Alzheimer disease", + "331.83": "Mild cognitive impairment", + "331.82": "Dementia with Lewy bodies", + "331.19": "Frontotemporal dementia", + "331.11": "Dementia in other disease", # added 2026-07-19 (Pick's disease) + "290": "Senile/presenile dementia", + "294": "Dementia in conditions classified elsewhere", +} + +_MCI_LABEL = "Mild cognitive impairment" + +_PRIORITY = [ + "Dementia with Lewy bodies", "Frontotemporal dementia", "Alzheimer disease", + "Vascular dementia", _MCI_LABEL, "Dementia in other disease", + "Unspecified dementia", "Senile/presenile dementia", + "Dementia in conditions classified elsewhere", +] + + +def _clean_icd_value(val): + """CSV cells missing a value come back as float NaN even under + dtype=str (pandas' missing-value representation doesn't respect the + dtype hint for empty cells). NaN is truthy in Python, so a naive + `row.get('ICD10','') or row.get('ICD9','')` silently returns NaN + instead of '' whenever ICD10 is missing — classify_family() then + calls .startswith() on a float and crashes. This is the normal case, + not an edge case: most patients have only one of ICD10/ICD9 filled. + Same fix as eda/site_subtype_breakdown.py's identical bug, caught + 2026-07-19 testing this module against the real ICD_codes_CI.csv.""" + if pd.isna(val): + return '' + return str(val) + + +def classify_family(icd10, icd9): + code = icd10 or icd9 + if not code: + return "unknown" + for prefix, family in CODE_FAMILY.items(): + if code.startswith(prefix): + return family + return f"other ({code})" + + +def dominant_family(families): + for p in _PRIORITY: + if p in families: + return p + return sorted(families)[0] if families else "unknown" + + +def load_icd_subtype_lookup(data_folder, verbose=False): + """ + Loads ICD_codes_CI.csv from the training data folder and returns a + dict: BDSPPatientID (str) -> dominant subtype (str). + + Graceful fallback, not a crash: returns an empty dict if the file is + missing. compute_mci_boosted_weights() treats any patient absent from + this lookup as "not MCI" (weight multiplier 1.0) — same behavior as + if the file were present but that patient simply had no ICD entries. + This means a missing file silently degrades to alpha-only weighting + (Entry 6 behavior) rather than crashing train_model() — deliberate, + since the ICD file is training-time-only interpretive data, not core + to the pipeline the way demographics.csv is. + """ + icd_path = os.path.join(data_folder, ICD_CODES_FILENAME) + if not os.path.exists(icd_path): + if verbose: + print(f' ! {ICD_CODES_FILENAME} not found in {data_folder} — ' + f'MCI-boosted weighting will fall back to alpha-only (Entry 6 behavior).') + return {} + df = pd.read_csv(icd_path, dtype=str) + by_patient = defaultdict(set) + for _, row in df.iterrows(): + pid = row['BDSPPatientID'] + fam = classify_family(_clean_icd_value(row.get('ICD10')), _clean_icd_value(row.get('ICD9'))) + by_patient[pid].add(fam) + return {pid: dominant_family(fams) for pid, fams in by_patient.items()} + + +def compute_mci_boosted_weights(y_train, age_train, patient_id_train, subtype_lookup, + alpha_age=1.0, beta_mci=0.75): + """ + Combines the Entry 6 age-based value weighting with an MCI-specific + boost. beta_mci=0.0 reproduces compute_value_weighted_sample_weights + exactly, unchanged — this is a pure additive extension, not a + replacement. + + Only positive training samples are ever boosted; negatives always get + weight 1.0, same convention as the age-based weighting it extends. + patient_id_train entries absent from subtype_lookup (including the + case where subtype_lookup is empty, e.g. ICD_codes_CI.csv missing) + are treated as "not MCI" — no boost, silent fallback to alpha-only + behavior for those patients. + """ + w = compute_value_weighted_sample_weights(y_train, age_train, alpha=alpha_age) + if beta_mci == 0.0: + return w + w = w.copy() + for i in range(len(y_train)): + if y_train[i] == 1: + pid = str(patient_id_train[i]) + if subtype_lookup.get(pid) == _MCI_LABEL: + w[i] *= (1.0 + beta_mci) + return w \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 28b6875..0f5a7c6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,8 @@ edfio==0.4.10 joblib==1.4.2 numpy==2.0.2 pandas==2.2.2 -scikit-learn==1.6.0 +scikit-learn==1.5.2 scipy==1.13.1 -tqdm==4.67.1 \ No newline at end of file +tqdm==4.67.1 +xgboost==2.1.3 + \ No newline at end of file diff --git a/team_code.py b/team_code.py index 502b03d..7c9396b 100644 --- a/team_code.py +++ b/team_code.py @@ -1,533 +1,406 @@ #!/usr/bin/env python -# Edit this script to add your team's code. Some functions are *required*, but you can edit most parts of the required functions, -# change or remove non-required functions, and add your own functions. +# Team Narnia — PhysioNet Challenge 2026 +# Entry 3: Added Platt-scaled calibration (CalibratedClassifierCV) and +# explicit threshold (0.12, tuned via loso_cv.py pooled sweep) +# replacing model.predict()'s default 0.5 cutoff. +# Entry 4: Added AgeResidualizer pipeline step (2 new features: +# CA_rate_age_residual, EEG_var_REM_Wake_age_residual). Both were +# previously written off as uninformative but age-residualized EDA +# (2026-07-02) showed real signal masked by an age confound. Fit at +# train time only, applied identically at inference — see +# features/age_residuals.py. Feature count: 48 -> 50. +# Entry 4 (fix, 2026-07-03): moved model-pipeline construction (XGBoost +# params, calibration, AgeResidualizer wiring) out of this file and +# into features/pipeline.py's build_pipeline(), shared with +# loso_cv.py. This closes a real bug where loso_cv.py had its own +# hand-copied Pipeline() that silently never got the AgeResidualizer +# step added here — LOSO was validating a stale 48-feature model +# while this file had already moved to 50. See features/pipeline.py +# for the full incident writeup. +# Entry 5 (2026-07-10): FIRST LARGE-TRAINING-SET SUBMISSION. Switched +# model_family from XGBoost to logistic regression +# (build_logreg_pipeline, features/pipeline.py). Evidence: on the +# large training set, logreg's mean age-conditioned AUROC (0.6002, +# both folds independently) beats every XGBoost/CatBoost variant +# tested (all clustered ~0.544), a real (2.36 sigma) effect, not +# noise — see learning_log.md, 2026-07-06/07 model-family +# diagnostic entries. Reward is unaffected either way (all model +# families land in the same 0.072-0.076 band at large scale) — this +# switch costs nothing on the metric that matters most and gains +# real ground on the metric that was weaker. +# THRESHOLD changed 0.12 -> 0.10 to match: logreg's own pooled +# threshold sweep picked a DIFFERENT optimum than XGBoost's — do +# not carry XGBoost's tuned value over by habit. +# xgboost stays in requirements.txt: features/pipeline.py still +# imports XGBClassifier unconditionally at module level (build_pipeline() +# still exists, just unused by this file now), so the dependency +# doesn't go away just because this file stopped calling it. +# CAVEAT: this is the first-ever large-set submission. Every number +# above is a LOSO estimate — zero large-set leaderboard data points +# exist yet to confirm the LOSO->leaderboard transfer at this scale. +# Treat this submission itself as the calibration point, not a +# confirmed result. +# Entry 6 (2026-07-19): Two independently-validated, compounding levers on +# top of the Entry 5 logreg pipeline — branch (a) feature work +# (NF1/NF2, time-of-night) and the age-banded threshold both closed +# with null/falsified results and are NOT reopened here (see +# learning_log_3, 2026-07-16 entries). +# 1. C changed 0.01 -> 0.001 (build_logreg_pipeline's C= arg): +# reward +35% relative at 3-site LOSO, AUROC flat. Reproduced +# exactly across two independent Kaggle runs. +# 2. Value-weighted sample weighting added at .fit() time, alpha=1.0 +# (compute_value_weighted_sample_weights, features/ +# sample_weighting.py — promoted from tools/reg_sweep.py, same +# shared-module precedent as AgeResidualizer): additional +12.9% +# relative reward on top of #1, AUROC still flat (max 0.21σ +# across all 3 folds). Reproduced bit-for-bit on an independent +# Kaggle run. train_model() has no fold rotation, so the weights +# are computed on the full training set passed in — this is +# itself the leakage-safe usage (see that function's docstring). +# Combined: +52% relative reward vs. the original shipped Entry 5 +# baseline (0.1354 vs. 0.089 at 3-site LOSO), AUROC flat throughout +# the entire chain (0.6449 -> 0.6485 -> 0.6459, all within noise). +# THRESHOLD changed 0.10 -> 0.08 to match the new probability +# distribution's own pooled threshold sweep optimum. +# Entry 7 (2026-07-19, TESTED BUT NOT SHIPPED): MCI-targeted sample +# weighting (compute_mci_boosted_weights, features/subtype_ +# weighting.py), beta_mci=0.75 layered on top of Entry 6's +# alpha=1.0. Validated on its own: +5.6% relative reward, S0001 +# AUROC flat, MCI sensitivity 76.8% -> 85.9% (tools/ablation_mci_ +# sample_weight.py). Deliberately held back from Entry 8 below — +# both changes touch the sample-weighting layer, and shipping +# them together would make it impossible to attribute any +# outcome to either one individually (the exact ambiguity that +# made Entry 6's own regression so hard to diagnose). beta_mci +# =0.75 was also tuned against the 50-feature model, not Entry +# 8's 39-feature one — the value would need re-validating against +# the new feature set, not carried over unchanged. features/ +# subtype_weighting.py remains in the repo, unused by this file, +# for exactly that future re-test. See learning_log.md, 2026- +# 07-19/20. +# Entry 8 (2026-07-20): drops 11 SIGN-FLIPPING features — coefficients +# that changed direction depending on which 2 of the 3 training +# sites fit the model (features/feature_selection.py, +# STAGE1_DROP_FEATURES) — direct, measured evidence of site- +# specific overfitting, from the cross-fold coefficient analysis. +# Validated via tools/ablation_drop_signflip_features.py "Stage +# 1": +8.2% relative reward, mean age-conditioned AUROC 0.6459 -> +# 0.6711, ALL THREE LOSO folds improved simultaneously (not a +# mixed result), cross-fold spread shrank 0.1264 -> 0.1145 (more +# STABLE, not just better on average — the specific signature +# this theory predicted). Reproduced bit-for-bit across an +# independent Kaggle kernel restart. THRESHOLD changed 0.08 -> +# 0.10 to match this config's own pooled reward-sweep optimum +# (t=0.10 gave best reward 0.1434, vs. the prior t=0.08's 0.1325 +# on the unchanged 50-feature model) — same "re-tune the +# threshold for the new probability distribution, don't carry +# the old value over" discipline as every prior threshold change. +# Stage 2 of that same ablation (also dropping CA_rate_age_ +# residual and Race_Black/White/Asian) scored higher still and +# clears the pre-committed 1.0-sigma gate outright — deliberately +# NOT included here. CA_rate_age_residual carries its own Entry 4 +# validation history; the Race_* features carry real, documented +# demographic signal that coefficient instability alone doesn't +# invalidate. Both deserve separate, specific justification +# before being dropped — not bundled in just because Stage 1 +# worked. See features/feature_selection.py header and +# learning_log.md, 2026-07-20, for the full reasoning. +# +# See features/FEATURES.md and LEARNING_LOG.md for full rationale. ################################################################################ -# -# Optional libraries, functions, and variables. You can change or remove them. -# +# Libraries ################################################################################ import joblib import numpy as np import os -from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor -from sklearn.pipeline import Pipeline -from sklearn.impute import SimpleImputer -import sys from tqdm import tqdm from helper_code import * +# Feature modules +from features.demographic import extract_demographic_features +from features.caisr_base import extract_caisr_base_features +from features.caisr_enriched import extract_caisr_enriched_features +from features.physiological_ratios import (extract_physiological_ratio_features, + N_RATIO_FEATURES) +from features.human import extract_human_annotations_features +from features.pipeline import build_logreg_pipeline +from features.sample_weighting import compute_value_weighted_sample_weights +from features.feature_selection import FeatureDropper, STAGE1_DROP_FEATURES +from features import (N_CAISR_BASE_FEATURES, N_CAISR_ENRICHED_FEATURES, + N_DEMOGRAPHIC_FEATURES, IDX_AGE) + +from sklearn.pipeline import Pipeline + ################################################################################ -# Path & Constant Configuration (Added for Robustness) +# Configuration ################################################################################ -# Get the absolute directory where this script is located -SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) - -# Build the absolute path to the CSV file relative to the script location +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) DEFAULT_CSV_PATH = os.path.join(SCRIPT_DIR, 'channel_table.csv') +# Fallback sizes for missing files +_N_BASE = N_CAISR_BASE_FEATURES # 12 +_N_ENRICHED = N_CAISR_ENRICHED_FEATURES # 11 +_N_RATIO = N_RATIO_FEATURES # 15 + +# Entry 8 — new pooled reward-sweep optimum for the 39-feature (Stage 1 +# sign-flip features dropped) probability distribution. Confirmed via +# tools/ablation_drop_signflip_features.py: reward peaks at t=0.10 +# (0.1434) — a different optimum than the prior 50-feature model's +# t=0.08, because dropping features reshapes the calibrated output, not +# just the ranking. See learning_log.md, 2026-07-20. +#verified THRESHOLD +THRESHOLD = 0.10 + +# Entry 6 — value-weighted sample-weighting boost applied at .fit() time +# (compute_value_weighted_sample_weights, features/sample_weighting.py). +# alpha=1.0 is the validated, reproduced value — do not change without a +# new sweep + reproduction run, same discipline as C and THRESHOLD above. +SAMPLE_WEIGHT_ALPHA = 1.0 + ################################################################################ -# -# Required functions. Edit these functions to add your code, but do not change the arguments for the functions. -# +# Required functions ################################################################################ -# Train your models. This function is *required*. You should edit this function to add your code, but do *not* change the arguments -# of this function. If you do not train one of the models, then you can return None for the model. - -# Train your model. def train_model(data_folder, model_folder, verbose, csv_path=DEFAULT_CSV_PATH): - # Find the data files. if verbose: print('Finding the Challenge data...') - patient_data_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) + patient_data_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) patient_metadata_list = find_patients(patient_data_file) - num_records = len(patient_metadata_list) + num_records = len(patient_metadata_list) if num_records == 0: raise FileNotFoundError('No data were provided.') - # Extract the features and labels from the data. if verbose: print('Extracting features and labels from the data...') - # Iterate over the records to extract the features and labels. - features = list() - labels = list() - - pbar = tqdm(range(num_records), desc="Extracting Features", unit="record", disable=not verbose) + features = [] + labels = [] + + pbar = tqdm(range(num_records), desc='Extracting Features', + unit='record', disable=not verbose) for i in pbar: try: - # Extract identifiers for this specific record - record = patient_metadata_list[i] + record = patient_metadata_list[i] patient_id = record[HEADERS['bids_folder']] site_id = record[HEADERS['site_id']] session_id = record[HEADERS['session_id']] if verbose: - pbar.set_postfix({"patient": patient_id}) - - # Load the patient data. - patient_data_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) - patient_data = load_demographics(patient_data_file, patient_id, session_id) - demographic_features = extract_demographic_features(patient_data) - - # Load signal data. - - # Load the physiological signal. - physiological_data_file = os.path.join(data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, site_id, f"{patient_id}_ses-{session_id}.edf") - # --- Check if the file actually exists before proceeding --- - if not os.path.exists(physiological_data_file): + pbar.set_postfix({'patient': patient_id}) + + # ── Demographics ───────────────────────────────────────────────── + demo_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) + patient_data = load_demographics(demo_file, patient_id, session_id) + demo_f = extract_demographic_features(patient_data) + + # ── Physiological EDF ───────────────────────────────────────────── + phys_file = os.path.join( + data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, + site_id, f'{patient_id}_ses-{session_id}.edf') + if not os.path.exists(phys_file): if verbose: - print(f" ! Missing physiological data for {patient_id}. Skipping...") - continue # skip record - physiological_data, physiological_fs = load_signal_data(physiological_data_file) - physiological_features = extract_physiological_features(physiological_data, physiological_fs, csv_path=csv_path) # This function can rename, re-reference, resample, etc. the signal data. - - # Load the algorithmic annotations. - algorithmic_annotations_file = os.path.join(data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, site_id, f"{patient_id}_ses-{session_id}_caisr_annotations.edf") - algorithmic_annotations, algorithmic_fs = load_signal_data(algorithmic_annotations_file) - algorithmic_features = extract_algorithmic_annotations_features(algorithmic_annotations) - - # Load the human annotations; these data will not be available in the hidden validation and test sets. - human_annotations_file = os.path.join(data_folder, HUMAN_ANNOTATIONS_SUBFOLDER, site_id, f"{patient_id}_ses-{session_id}_expert_annotations.edf") - human_annotations, human_fs = load_signal_data(human_annotations_file) - human_features = extract_human_annotations_features(human_annotations) - - # Load the diagnoses; these data will not be available in the hidden validation and test sets. - diagnosis_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) - label = load_diagnoses(diagnosis_file, patient_id) - - # Store the features and labels, but the human annotations are not available on the hidden validation and test sets. + tqdm.write( + f' ! Missing physiological EDF for {patient_id} — skipping.') + continue + phys_data, phys_fs = load_signal_data(phys_file) + + # ── CAISR Annotations ───────────────────────────────────────────── + algo_file = os.path.join( + data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + site_id, f'{patient_id}_ses-{session_id}_caisr_annotations.edf') + if os.path.exists(algo_file): + algo_data, _ = load_signal_data(algo_file) + caisr_base = extract_caisr_base_features(algo_data) + caisr_enriched = extract_caisr_enriched_features(algo_data) + # Ratio features require BOTH phys + CAISR + ratio_f = extract_physiological_ratio_features( + phys_data, phys_fs, algo_data, csv_path=csv_path) + else: + tqdm.write( + f'Error loading EDF file: [Errno 2] No such file or directory: ' + f"'{algo_file}'") + caisr_base = np.full(_N_BASE, float('nan')) + caisr_enriched = np.full(_N_ENRICHED, float('nan')) + ratio_f = np.full(_N_RATIO, float('nan')) + + # ── Human Annotations (training only — not used in model) ───────── + human_file = os.path.join( + data_folder, HUMAN_ANNOTATIONS_SUBFOLDER, + site_id, f'{patient_id}_ses-{session_id}_expert_annotations.edf') + if os.path.exists(human_file): + human_data, _ = load_signal_data(human_file) + _ = extract_human_annotations_features(human_data) + + # ── Label ───────────────────────────────────────────────────────── + label = load_diagnoses( + os.path.join(data_folder, DEMOGRAPHICS_FILE), patient_id) + if label == 0 or label == 1: - features.append(np.hstack([demographic_features, physiological_features, algorithmic_features])) + # 48-feature extraction vector. AgeResidualizer (pipeline + # step, Entry 4) appends the 2 age-residualized features + # at model.fit() time — do not hstack them here. + features.append(np.hstack([ + demo_f, + caisr_base, + caisr_enriched, + ratio_f, + ])) labels.append(label) - if 'physiological_data' in locals(): - del physiological_data - if 'algorithmic_annotations' in locals(): - del algorithmic_annotations + if 'phys_data' in locals(): del phys_data + if 'algo_data' in locals(): del algo_data + if 'human_data' in locals(): del human_data except Exception as e: - # If an error occurs (e.g., a record is corrupted), log it and move to the next - tqdm.write(f" !!! Error processing record {i+1} ({patient_id}): {e}") + tqdm.write(f' !!! Error on record {i+1} ({patient_id}): {e}') continue pbar.close() features = np.asarray(features, dtype=np.float32) - labels = np.asarray(labels, dtype=bool) + labels = np.asarray(labels, dtype=bool) - # Train the models on the features. if verbose: - print('Training the model on the data...') - - # This very simple model trains a random forest model with very simple features. - - # Define the parameters for the random forest classifier and regressor. - n_estimators = 12 # Number of trees in the forest. - max_leaf_nodes = 34 # Maximum number of leaf nodes in each tree. - random_state = 56 # Random state; set for reproducibility. - - # Created a Pipeline wrapping SimpleImputer and RandomForestClassifier. - # This automatically injects median values for any missing data (NaN) during both fit() and predict() calls. - rf = RandomForestClassifier( - n_estimators=n_estimators, max_leaf_nodes=max_leaf_nodes, random_state=random_state) - - model = Pipeline([ - ('imputer', SimpleImputer(strategy='median')), - ('classifier', rf) - ]) - - # Fit the model. - model.fit(features, labels) - - # Create a folder for the model if it does not already exist. - os.makedirs(model_folder, exist_ok=True) + n_pos = int(labels.sum()) + n_neg = int((~labels).sum()) + print(f'Training on {len(labels)} patients ' + f'({n_pos} positive, {n_neg} negative)...') + print(f'Feature vector shape: {features.shape}') + + # ── Model pipeline (Entry 8: logreg, C=0.001, value-weighted sample + # weighting at alpha=1.0, calibrated, AgeResidualizer on, Stage 1 + # sign-flip features dropped) ─────────────────────────────────────── + # Built via features/pipeline.py's build_logreg_pipeline() — the same + # shared definition loso_cv.py/reg_sweep.py uses, so the validation + # harness and this submission can never silently diverge on what "the + # model" is (see features/pipeline.py header for the incident that + # motivated this discipline in the first place). C=0.001 and alpha=1.0 + # are both independently reproduced on Kaggle (see learning_log_3, + # 2026-07-16 entries) — do not change either without a fresh + # reproduction run. + # + # FeatureDropper (features/feature_selection.py) is inserted + # immediately after 'age_residual' and before 'imputer' — it MUST sit + # there, not earlier, so AgeResidualizer's own source features (e.g. + # CA_rate, used to compute CA_rate_age_residual, which is KEPT here) + # are still present when AgeResidualizer runs, even though the raw + # CA_rate column itself is one of the 11 dropped. Do not reorder these + # steps. + _base_model = build_logreg_pipeline(labels, calibrated=True, C=0.001) + model = Pipeline( + [_base_model.steps[0], ('feature_dropper', FeatureDropper(STAGE1_DROP_FEATURES))] + + _base_model.steps[1:] + ) + + # train_model() has no fold rotation — it trains once on whatever + # training set the organizers provide — so computing the weights on + # the full `labels`/age_train passed in here IS the leakage-safe usage + # (see compute_value_weighted_sample_weights' docstring for why this + # differs from reg_sweep.py's per-LOSO-fold call, which only ever sees + # that fold's training split). + age_train = features[:, IDX_AGE] + sample_weight = compute_value_weighted_sample_weights( + labels, age_train, alpha=SAMPLE_WEIGHT_ALPHA) + + # 'classifier' is build_logreg_pipeline's final step name — sklearn's + # Pipeline routes fit_params via the 'stepname__paramname' convention. + # CalibratedClassifierCV.fit() accepts and forwards sample_weight to + # the wrapped LogisticRegression. (2026-07-16 incident: a leftover + # duplicate unweighted .fit() call was caught and removed during + # implementation on Kaggle — it would have silently discarded the + # weighting entirely with zero errors raised. There is only ONE .fit() + # call below; keep it that way.) + model.fit(features, labels, classifier__sample_weight=sample_weight) - # Save the model. + os.makedirs(model_folder, exist_ok=True) save_model(model_folder, model) if verbose: print('Done.') print() -# Load your trained models. This function is *required*. You should edit this function to add your code, but do *not* change the -# arguments of this function. If you do not train one of the models, then you can return None for the model. + def load_model(model_folder, verbose): - model_filename = os.path.join(model_folder, 'model.sav') - model = joblib.load(model_filename) - return model + return joblib.load(os.path.join(model_folder, 'model.sav')) + -# Run your trained model. This function is *required*. You should edit this function to add your code, but do *not* change the -# arguments of this function. def run_model(model, record, data_folder, verbose): - # Load the model. model = model['model'] - # Extract identifiers from the record dictionary patient_id = record[HEADERS['bids_folder']] site_id = record[HEADERS['site_id']] session_id = record[HEADERS['session_id']] - # Load the patient data. - patient_data_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) - patient_data = load_demographics(patient_data_file, patient_id, session_id) - demographic_features = extract_demographic_features(patient_data) + # ── Demographics ────────────────────────────────────────────────────────── + demo_file = os.path.join(data_folder, DEMOGRAPHICS_FILE) + patient_data = load_demographics(demo_file, patient_id, session_id) + demo_f = extract_demographic_features(patient_data) - # Load the signal data. - phys_file = os.path.join(data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, site_id, f"{patient_id}_ses-{session_id}.edf") + # ── Physiological EDF ───────────────────────────────────────────────────── + phys_file = os.path.join( + data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, + site_id, f'{patient_id}_ses-{session_id}.edf') + phys_data = phys_fs = None if os.path.exists(phys_file): phys_data, phys_fs = load_signal_data(phys_file) - # Ensure csv_path is accessible or defined - physiological_features = extract_physiological_features(phys_data, phys_fs) - else: - physiological_features = np.full(49, float('nan')) # Fallback if signal data does not exist - # Load the algorithmic annotations. - algo_file = os.path.join(data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, site_id, f"{patient_id}_ses-{session_id}_caisr_annotations.edf") + # ── CAISR Annotations ───────────────────────────────────────────────────── + algo_file = os.path.join( + data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + site_id, f'{patient_id}_ses-{session_id}_caisr_annotations.edf') if os.path.exists(algo_file): algo_data, _ = load_signal_data(algo_file) - algorithmic_features = extract_algorithmic_annotations_features(algo_data) + caisr_base = extract_caisr_base_features(algo_data) + caisr_enriched = extract_caisr_enriched_features(algo_data) + if phys_data is not None: + ratio_f = extract_physiological_ratio_features( + phys_data, phys_fs, algo_data) + else: + ratio_f = np.full(_N_RATIO, float('nan')) else: - algorithmic_features = np.full(12, float('nan')) # Fallback if algorithmic annotations do not exist - - features = np.hstack([demographic_features, physiological_features, algorithmic_features]).reshape(1, -1) - - # Apply the model to the features. - binary_output = model.predict(features)[0] + caisr_base = np.full(_N_BASE, float('nan')) + caisr_enriched = np.full(_N_ENRICHED, float('nan')) + ratio_f = np.full(_N_RATIO, float('nan')) + + # NOTE: this is the 48-feature extraction vector (demo/base/enriched/ + # ratios). The loaded `model` pipeline's age_residual step appends the + # 2 Entry 4 age-residualized features automatically — do not hstack + # them here, and do not hardcode 50 anywhere in this function. + features = np.hstack([ + demo_f, + caisr_base, + caisr_enriched, + ratio_f, + ]).reshape(1, -1) + + # ── Entry 8: calibrated probability + explicit threshold ────────────────── + # model.predict_proba() is already calibrated — calibration is baked into + # the pipeline itself (see train_model()). model.predict() is NOT used + # here: its default 0.5 cutoff is far too conservative for the reward + # metric at low local prevalence (this was Entry 1/2's mistake — reward + # 0.011 despite reasonable AUROC). THRESHOLD is set from tools/ablation_ + # drop_signflip_features.py's pooled reward sweep for the 39-feature + # (Stage 1 dropped) probability distribution specifically (t=0.10) — not + # carried over from Entry 6's t=0.08, since dropping features reshapes + # the calibrated output, not just the ranking. probability_output = model.predict_proba(features)[0][1] + binary_output = int(probability_output > THRESHOLD) return binary_output, probability_output + ################################################################################ -# -# Optional functions. You can change or remove these functions and/or add new functions. -# +# Utilities ################################################################################ -def extract_demographic_features(data): - """ - Extracts and encodes demographic features from a metadata dictionary. - - Inputs: - data (dict): A dictionary containing patient metadata (e.g., from a CSV row). - - Returns: - np.array: A feature vector of length 11: - - [0]: Age (Continuous) - - [1:4]: Sex (One-hot: Female, Male, Other/Unknown) - - [4:9]: Race (One-hot: Asian, Black, Other, Unavailable, White) - - [9]: BMI (Continuous) - """ - # 1. Age - age = load_age(data) - age = np.array([age]) - - # 2. Sex feature (one-hot encoding for Female, Male, Other/Unknown) - # Uses lowercase prefix matching to handle variants like 'F', 'Female', 'M', or 'Male' - sex = load_sex(data, standardize=True) - sex_vec = np.zeros(3) - if sex == 'Female': - sex_vec[0] = 1 # Index 0: Female - elif sex == 'Male': - sex_vec[1] = 1 # Index 1: Male - else: - sex_vec[2] = 1 # Index 2: Other/Unknown - - # 3. Race One-Hot Encoding (5 dimensions) - # Standardizes the raw text into one of five categories using the helper function - race = load_race(data, standardize=True) - race_vec = np.zeros(5) - # Pre-defined mapping for index consistency - if race == 'Asian': - race_vec[0] = 1 - elif race == 'Black': - race_vec[1] = 1 - elif race == 'Others': - race_vec[2] = 1 - elif race == 'Unavailable': - race_vec[3] = 1 - elif race == 'White': - race_vec[4] = 1 - else: - race_vec[2] = 1 # Default to 'Others' for any unrecognized - - # 4. Body mass index (BMI) - bmi = load_bmi(data) - bmi = np.array([bmi]) - - # 5. Concatenate all components into a single vector (1 + 3 + 5 + 1 = 10) - - return np.concatenate([age, sex_vec, race_vec, bmi]) - - -def extract_physiological_features(physiological_data, physiological_fs, csv_path=DEFAULT_CSV_PATH): - """ - Standardizes channels and extracts statistical/spectral features. - """ - original_labels = list(physiological_data.keys()) - - # Step 1: Load rules and standardize names - # Note: Use script-relative path or absolute path for robustness - rename_rules = load_rename_rules(os.path.abspath(csv_path)) - rename_map, cols_to_drop = standardize_channel_names_rename_only(original_labels, rename_rules) - - # Step 2: Apply renaming to BOTH signals and their corresponding FS - processed_channels = {} - processed_fs = {} - for old_label, data in physiological_data.items(): - if old_label in cols_to_drop: - continue - new_label = rename_map.get(old_label, old_label.lower()) - processed_channels[new_label] = data - # Mapping the sampling rate to the new label - if old_label in physiological_fs: - processed_fs[new_label] = physiological_fs[old_label] - else: - # Report error and stop if no FS is found for a kept channel - raise KeyError(f"Sampling frequency (fs) not found for channel '{old_label}' ") - - if 'physiological_data' in locals(): del physiological_data - - # Step 3: Construct Bipolar Derivations - bipolar_configs = [ - ('f3-m2', 'f3', ['m2']), ('f4-m1', 'f4', ['m1']), - ('c3-m2', 'c3', ['m2']), ('c4-m1', 'c4', ['m1']), - ('o1-m2', 'o1', ['m2']), ('o2-m1', 'o2', ['m1']), - ('e1-m2', 'e1', ['m2']), ('e2-m1', 'e2', ['m1']), - ('chin1-chin2', 'chin 1', ['chin 2']), - ('lat', 'lleg+', ['lleg-']), ('rat', 'rleg+', ['rleg-']) - ] - - for target, pos, neg_list in bipolar_configs: - # 1. Skip if target already exists or pos channel missing - if target in processed_channels or pos not in processed_channels: - continue - - # 2. Check all neg channels exist - if not all(n in processed_channels for n in neg_list): - continue - - # 3. Check sampling rate consistency - all_involved = [pos] + neg_list - fs_values = [processed_fs[ch] for ch in all_involved] - - if len(set(fs_values)) > 1: - raise ValueError(f"Sampling rate mismatch for {target}: {dict(zip(all_involved, fs_values))}") - - # 4. Derive bipolar signal - ref_sig = processed_channels[neg_list[0]] if len(neg_list) == 1 else tuple(processed_channels[n] for n in neg_list) - - derived = derive_bipolar_signal(processed_channels[pos], ref_sig) - - if derived is not None: - processed_channels[target] = derived - processed_fs[target] = processed_fs[pos] - - leads_to_check = { - 'eeg': ['f3-m2', 'f4-m1', 'c3-m2', 'c4-m1'], - 'eog': ['e1-m2', 'e2-m1'], - 'chin': ['chin1-chin2', 'chin'], - 'leg': ['lat', 'rat'], - 'ecg': ['ecg', 'ekg'], - 'resp': ['airflow', 'ptaf', 'abd', 'chest'], - 'spo2': ['spo2', 'sao2'] # Added sao2 as fallback for spo2 - } - - final_features = [] - for lead_type, candidates in leads_to_check.items(): - sig = None - fs = None - - # Identify the first available candidate - for candidate in candidates: - if candidate in processed_channels and processed_channels[candidate] is not None: - sig = processed_channels[candidate] - fs = processed_fs.get(candidate) - break - - if sig is not None and len(sig) > 1: - # --- Time Domain Features (Very Fast) --- - std_val = np.std(sig) - mav_val = np.mean(np.abs(sig)) - - # Zero Crossing Rate (Proxy for frequency/slowing) - zcr = np.mean(np.diff(np.sign(sig)) != 0) - - # Root Mean Square - rms = np.sqrt(np.mean(sig**2)) - - # Signal Activity (Variance) - activity = np.var(sig) - - # Mobility (Hjorth Parameter) - Proxy for mean frequency - # sqrt(var(diff(sig)) / var(sig)) - diff_sig = np.diff(sig) - mobility = np.sqrt(np.var(diff_sig) / activity) if activity > 0 else 0.0 - - # Complexity (Hjorth Parameter) - Proxy for bandwidth - diff2_sig = np.diff(diff_sig) - var_d2 = np.var(diff2_sig) - var_d1 = np.var(diff_sig) - complexity = (np.sqrt(var_d2 / var_d1) / mobility) if (var_d1 > 0 and mobility > 0) else 0.0 - - final_features.extend([std_val, mav_val, zcr, rms, activity, mobility, complexity]) - - else: - # Padding: 7 features per lead type - final_features.extend([float('nan')] * 7) - - if 'processed_channels' in locals(): del processed_channels - - return np.array(final_features) - -def extract_algorithmic_annotations_features(algo_data): - """ - Extracts sleep architecture and event density features from CAISR outputs. - Output vector length: 12 - """ - if not algo_data: - return np.full(12, float('nan')) - - features = [] - - # --- 1. Respiratory & Arousal Event Densities --- - # Total duration in hours (assuming 1Hz for event traces) - # If the signal exists, we calculate events per hour (Index) - total_hours = len(algo_data.get('resp_caisr', [])) / 3600.0 - - def count_discrete_events(key): - if key not in algo_data or total_hours <= 0: - return float('nan') - - sig = algo_data[key].astype(float) - # Create a binary mask: 1 if there is an event, 0 if not - binary_sig = (sig > 0).astype(int) - - # Detect rising edges: 0 to 1 transition - # diff will be 1 at the start of an event, -1 at the end - diff = np.diff(binary_sig, prepend=0) - num_events = np.count_nonzero(diff == 1) - - return num_events / total_hours - - ahi_auto = count_discrete_events('resp_caisr') # Automated Apnea-Hypopnea Index - arousal_auto = count_discrete_events('arousal_caisr') # Automated Arousal Index - limb_auto = count_discrete_events('limb_caisr') # Automated Limb Movement Index - - features.extend([ahi_auto, arousal_auto, limb_auto]) - - # --- 2. Sleep Architecture (from stage_caisr) --- - # Standard labels: 5=W, 4=R, 3=N1, 2=N2, 1=N3 (or similar mapping) - stages = algo_data.get('stage_caisr', np.array([])) - # Filter out invalid/background values (like the 9.0 in your sample) - valid_stages = stages[stages < 9.0] - - if len(valid_stages) > 0: - total_epochs = len(valid_stages) - # Percentage of each stage - w_pct = np.mean(valid_stages == 5) - r_pct = np.mean(valid_stages == 4) - n1_pct = np.mean(valid_stages == 3) - n2_pct = np.mean(valid_stages == 2) - n3_pct = np.mean(valid_stages == 1) - - # Sleep Efficiency: (N1+N2+N3+R) / Total - efficiency = np.mean((valid_stages >= 1) & (valid_stages <= 4)) - else: - w_pct = n1_pct = n2_pct = n3_pct = r_pct = efficiency = float('nan') - - features.extend([w_pct, n1_pct, n2_pct, n3_pct, r_pct, efficiency]) - - # --- 3. Model Confidence / Uncertainty --- - # Mean probability of Wake and REM (indicators of sleep stability) - # We use the raw probability traces - prob_w = np.mean(algo_data.get('caisr_prob_w', [float('nan')])) - prob_n3 = np.mean(algo_data.get('caisr_prob_n3', [float('nan')])) - prob_arous = np.mean(algo_data.get('caisr_prob_arous', [float('nan')])) - - # Standardize '9.0' or other filler values to NaN - clean_prob = lambda x: x if x < 1.0 else float('nan') - features.extend([clean_prob(prob_w), clean_prob(prob_n3), clean_prob(prob_arous)]) - - return np.array(features) - -def extract_human_annotations_features(human_data): - """ - Extracts features from expert-scored human annotations. - Output vector length: 12 (to match algorithmic feature length) - """ - # If data is missing (common in hidden test sets), return a zero vector - if not human_data or 'resp_expert' not in human_data: - return np.full(12, float('nan')) - - features = [] - - # --- 1. Human Event Indices (Events per Hour) --- - # Total duration in hours based on 1Hz signal - total_seconds = len(human_data.get('resp_expert', [])) - total_hours = total_seconds / 3600.0 - - def count_discrete_events(key): - if key not in human_data or total_hours <= 0: - return float('nan') - sig = (human_data[key] > 0).astype(int) - # Identify the start of each continuous event block - diff = np.diff(sig, prepend=0) - return np.count_nonzero(diff == 1) / total_hours - - ahi_human = count_discrete_events('resp_expert') # Human AHI - arousal_human = count_discrete_events('arousal_expert') # Human Arousal Index - limb_human = count_discrete_events('limb_expert') # Human PLMI - - features.extend([ahi_human, arousal_human, limb_human]) - - # --- 2. Human Sleep Architecture --- - # Standard labels: 0=W, 1=N1, 2=N2, 3=N3, 4=R, 5=Unknown/Movement - stages = human_data.get('stage_expert', np.array([])) - - # Filter out label 5 (often used by experts for movement/unscored) - valid_mask = (stages < 9.0) - valid_stages = stages[valid_mask] - - if len(valid_stages) > 0: - w_pct = np.mean(valid_stages == 5) - r_pct = np.mean(valid_stages == 4) - n1_pct = np.mean(valid_stages == 3) - n2_pct = np.mean(valid_stages == 2) - n3_pct = np.mean(valid_stages == 1) - efficiency = np.mean(valid_stages > 0) - else: - w_pct = n1_pct = n2_pct = n3_pct = r_pct = efficiency = float('nan') - - features.extend([w_pct, n1_pct, n2_pct, n3_pct, r_pct, efficiency]) - - # --- 3. Fragmentation & Stability (Replacing Probabilities) --- - # These metrics quantify how "broken" the sleep is, which is a key marker. - if len(valid_stages) > 1: - # Number of stage transitions - transitions = np.count_nonzero(np.diff(valid_stages)) / total_hours - # Wake After Sleep Onset (WASO) proxy: non-zero stages followed by zero - waso_minutes = (np.count_nonzero(valid_stages == 0) * 30) / 60.0 - # REM Latency (epochs until first REM) - rem_indices = np.where(valid_stages == 4)[0] - rem_latency = rem_indices[0] if len(rem_indices) > 0 else float('nan') - else: - transitions = waso_minutes = rem_latency = float('nan') - - features.extend([transitions, waso_minutes, rem_latency]) - - return np.array(features) - - -# Save your trained model. def save_model(model_folder, model): - d = {'model': model} - filename = os.path.join(model_folder, 'model.sav') - joblib.dump(d, filename, protocol=0) \ No newline at end of file + joblib.dump({'model': model}, + os.path.join(model_folder, 'model.sav'), + protocol=0) \ No newline at end of file diff --git a/tools/check_submission_files.py b/tools/check_submission_files.py new file mode 100644 index 0000000..4082396 --- /dev/null +++ b/tools/check_submission_files.py @@ -0,0 +1,347 @@ +""" +check_submission_files.py +Team Narnia — PhysioNet Challenge 2026 + +Verifies every file the Docker submission actually needs is present, +before you commit/push/submit. Designed to be run from a cluttered repo +root (dev_model/, dev_model_v2/, docker_outputs/, dev_subset/, etc. all +present) without being confused by any of that — it only checks for the +specific files the submission depends on. + +Two categories checked: + 1. Fixed set — files the Challenge harness needs regardless of what + team_code.py does (locked files, Dockerfile, requirements.txt, + channel_table.csv). + 2. Dynamic set — parsed directly from team_code.py's own import + statements, so this list can't drift out of sync with the actual + code the way a hardcoded list could. If team_code.py changes what + it imports, this script's requirements update automatically next run. + +Also does two cheap sanity checks while it's at it: + - Warns if team_code.py imports anything from a known dev-only tool + (loso_cv, build_dev_subset, verify_label_integrity, phase1_eda, + age_residualized_eda) — that would mean the submission accidentally + depends on a script that isn't meant to ship. + - Cross-checks that third-party packages actually imported by + team_code.py / features/*.py are listed in requirements.txt. + +Optionally also runs the ratchet check (see ratchet_check.py): compares a +fresh loso_cv.py result against a hand-curated "best confirmed so far" +baseline (ratchet_baselines.json), using Hanley-McNeil SE so a noisy dip +doesn't get mistaken for a real regression. Only runs if --loso-results is +passed — if omitted, this is a visible WARNING, not a silent skip, since +submitting without checking against the ratchet is exactly the gap this +was built to close. + +Usage: + python check_submission_files.py [--repo-root .] + python check_submission_files.py --loso-results loso_results.csv \\ + --ratchet-baseline small_entry3 [--candidate-reward 0.11] + +Exit code: 0 if everything required is present AND (if run) the ratchet +check passes, 1 otherwise (so this can gate a pre-commit hook or a +submission checklist script if you want). +""" + +import argparse +import ast +import os +import re +import sys +from pathlib import Path + +from ratchet_check import check_ratchet, _print_report + +# Files the Challenge harness needs regardless of team_code.py's contents. +# These are the "do not edit" files from CLAUDE.md plus build/config files. +FIXED_REQUIRED = [ + "team_code.py", + "train_model.py", + "run_model.py", + "helper_code.py", + "evaluate_model.py", + "channel_table.csv", + "requirements.txt", + "Dockerfile", +] + +# Present in the documented file structure but not imported/called at +# runtime by train_model.py or run_model.py — flagged as recommended, +# not build-breaking if absent. +RECOMMENDED = [ + "create_labels.py", + ".dockerignore", +] + +# Dev-only tools that should never be an import dependency of team_code.py. +# Their presence as a FILE in the repo is fine; their presence as an +# IMPORT inside team_code.py would mean the submission depends on +# something not meant to ship. +DEV_ONLY_MODULES = { + "loso_cv", "build_dev_subset", "verify_label_integrity", + "phase1_eda", "age_residualized_eda", "entry3_calibration", +} + +# import-name -> requirements.txt package-name, where they differ. +IMPORT_TO_PACKAGE = { + "sklearn": "scikit-learn", + "yaml": "pyyaml", +} + +# Standard library modules — never expected in requirements.txt. +STDLIB_SKIP = { + "os", "sys", "re", "ast", "argparse", "json", "csv", "collections", + "datetime", "itertools", "functools", "pathlib", "typing", "math", + "time", "warnings", "copy", "io", "abc", +} + + +def _parse_imports(source: str, filename: str): + """ + Parses one file's AST for import statements. Returns + (feature_module_paths, third_party_packages, dev_only_imports) — kept + as a standalone function so it can be called once per file during + transitive BFS resolution below. + """ + tree = ast.parse(source, filename=filename) + + feature_module_paths = set() + third_party_packages = set() + dev_only_imports = set() + + for node in ast.walk(tree): + if isinstance(node, ast.ImportFrom): + module = node.module or "" + top_level = module.split(".")[0] + + if module == "features": + feature_module_paths.add("features/__init__.py") + elif module.startswith("features."): + submodule = module.split(".", 1)[1] + feature_module_paths.add(f"features/{submodule}.py") + feature_module_paths.add("features/__init__.py") + elif top_level in DEV_ONLY_MODULES: + dev_only_imports.add(top_level) + elif top_level not in STDLIB_SKIP and top_level not in ( + "helper_code", "team_code", "evaluate_model", "train_model", "run_model" + ): + third_party_packages.add(top_level) + + elif isinstance(node, ast.Import): + for alias in node.names: + top_level = alias.name.split(".")[0] + if top_level in DEV_ONLY_MODULES: + dev_only_imports.add(top_level) + elif top_level not in STDLIB_SKIP: + third_party_packages.add(top_level) + + return feature_module_paths, third_party_packages, dev_only_imports + + +def parse_team_code_imports(team_code_path: Path): + """ + Parses team_code.py's imports, then TRANSITIVELY parses every + features/*.py file discovered that way for further `features.X` + imports, repeating until no new files are found (BFS). This closes a + real gap (found 2026-07-08, on a live merged repo): team_code.py + importing `features.pipeline`, which itself imports + `features.age_residuals`, was going COMPLETELY UNDETECTED as a required + file by a single-file-only parse — that version would PASS a submission + missing a real runtime dependency, only failing at Docker build/import + time instead of at this pre-commit check. (This is the same failure + mode this file's own header comment already claimed was fixed — + it wasn't, in the code that was actually shipping. Fixed for real now.) + + Returns: + feature_module_paths: set of relative file paths under features/ + that must exist, resolved transitively. + third_party_packages: set of top-level third-party import names, + collected from team_code.py AND every transitively-discovered + features/*.py file (a feature module can import a package + team_code.py itself never mentions — e.g. features/pipeline.py + importing sklearn.linear_model, which team_code.py doesn't). + dev_only_imports: set of dev-only module names imported anywhere in + the transitive closure (should be empty). + """ + repo_root = team_code_path.parent + source = team_code_path.read_text() + + feature_module_paths, third_party_packages, dev_only_imports = _parse_imports( + source, str(team_code_path)) + + parsed = set() + to_parse = set(feature_module_paths) + while to_parse: + rel_path = to_parse.pop() + if rel_path in parsed: + continue + parsed.add(rel_path) + + full_path = repo_root / rel_path + if not full_path.exists(): + continue # reported as MISSING by the caller + + sub_source = full_path.read_text() + sub_features, sub_third_party, sub_dev_only = _parse_imports( + sub_source, str(full_path)) + + third_party_packages |= sub_third_party + dev_only_imports |= sub_dev_only + + newly_found = sub_features - feature_module_paths + feature_module_paths |= sub_features + to_parse |= newly_found + + return feature_module_paths, third_party_packages, dev_only_imports + + +def check_requirements_coverage(repo_root: Path, packages: set): + """ + Cross-checks that each third-party package imported by team_code.py + appears (case-insensitively, substring match) somewhere in + requirements.txt. Returns list of packages NOT found. + """ + req_path = repo_root / "requirements.txt" + if not req_path.exists(): + return sorted(packages) # everything is "missing" if the file itself is gone + + req_text = req_path.read_text().lower() + missing = [] + for pkg in sorted(packages): + pkg_name = IMPORT_TO_PACKAGE.get(pkg, pkg) + if pkg_name.lower() not in req_text: + missing.append(pkg) + return missing + + +def run(repo_root: str, loso_results: str = None, ratchet_baseline: str = None, + baselines_file: str = "ratchet_baselines.json", candidate_reward: float = None, + sigma_threshold: float = 1.0): + root = Path(repo_root).resolve() + print(f"Checking submission files in: {root}\n") + + problems = [] + warnings = [] + + # ── Fixed required files ───────────────────────────────────────────────── + print("Fixed required files:") + for rel_path in FIXED_REQUIRED: + full = root / rel_path + ok = full.exists() and full.is_file() + status = "OK" if ok else "MISSING" + print(f" [{status:^7}] {rel_path}") + if not ok: + problems.append(rel_path) + + print("\nRecommended (not build-breaking if absent):") + for rel_path in RECOMMENDED: + full = root / rel_path + ok = full.exists() and full.is_file() + status = "OK" if ok else "missing" + print(f" [{status:^7}] {rel_path}") + if not ok: + warnings.append(f"Recommended file missing: {rel_path}") + + # ── Dynamic: parse team_code.py's actual imports ───────────────────────── + team_code_path = root / "team_code.py" + if not team_code_path.exists(): + print("\nCannot parse team_code.py imports — file is missing (see above).") + feature_paths, third_party, dev_only = set(), set(), set() + else: + print("\nFeature modules required by team_code.py's actual imports:") + feature_paths, third_party, dev_only = parse_team_code_imports(team_code_path) + for rel_path in sorted(feature_paths): + full = root / rel_path + ok = full.exists() and full.is_file() + status = "OK" if ok else "MISSING" + print(f" [{status:^7}] {rel_path}") + if not ok: + problems.append(rel_path) + + # ── Dev-only import leakage check ──────────────────────────────────────── + if dev_only: + print(f"\n!! WARNING: team_code.py imports from dev-only tool(s): " + f"{', '.join(sorted(dev_only))}") + print(" These are not meant to ship as a submission dependency.") + warnings.append(f"team_code.py imports dev-only module(s): {sorted(dev_only)}") + else: + print("\nNo dev-only tool imports found in team_code.py — clean.") + + # ── requirements.txt coverage check ────────────────────────────────────── + if third_party: + print(f"\nThird-party packages imported by team_code.py: {sorted(third_party)}") + missing_reqs = check_requirements_coverage(root, third_party) + if missing_reqs: + print(f"!! WARNING: not found in requirements.txt: {missing_reqs}") + warnings.append(f"Packages imported but not in requirements.txt: {missing_reqs}") + else: + print("All imported packages found in requirements.txt — clean.") + + # ── Ratchet check (optional — see ratchet_check.py) ────────────────────── + if loso_results: + print(f"\nRatchet check ({loso_results} vs baseline '{ratchet_baseline}'):\n") + try: + passed, results, messages = check_ratchet( + loso_results, ratchet_baseline, baselines_file, + candidate_reward=candidate_reward, sigma_threshold=sigma_threshold, + ) + _print_report(ratchet_baseline, results, messages) + if not passed: + problems.append( + f"Ratchet check FAILED against baseline '{ratchet_baseline}' " + f"— see report above.") + except (FileNotFoundError, KeyError, ValueError) as e: + print(f"!! RATCHET CHECK ERROR: {e}") + problems.append(f"Ratchet check could not run: {e}") + else: + print("\n!! WARNING: no --loso-results passed — ratchet check SKIPPED.") + print(" This submission has not been checked against the best confirmed") + print(" LOSO baseline. Run with --loso-results and --ratchet-baseline") + print(" before submitting, or confirm you're deliberately skipping this.") + warnings.append("Ratchet check skipped — no --loso-results provided.") + + # ── Summary ─────────────────────────────────────────────────────────────── + print(f"\n{'='*60}") + if not problems: + print("PASS — all required files present.") + else: + print(f"FAIL — {len(problems)} required file(s) missing:") + for p in problems: + print(f" - {p}") + if warnings: + print(f"\n{len(warnings)} warning(s) (not build-breaking, worth a look):") + for w in warnings: + print(f" - {w}") + print(f"{'='*60}") + + return len(problems) == 0 + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--repo-root", default=".", help="Path to the repo root (default: current directory)") + parser.add_argument("--loso-results", default=None, + help="Path to loso_cv.py's per-fold results CSV. If provided, runs the " + "ratchet check (see ratchet_check.py) against --ratchet-baseline.") + parser.add_argument("--ratchet-baseline", default=None, + help="Key into ratchet_baselines.json (e.g. small_entry3). " + "Required if --loso-results is passed.") + parser.add_argument("--baselines-file", + default=os.path.join(os.path.dirname(os.path.abspath(__file__)), "ratchet_baselines.json")) + parser.add_argument("--candidate-reward", type=float, default=None, + help="Pooled reward at your chosen threshold, if you have one. Optional — " + "without it, only the AUROC ratchet is checked.") + parser.add_argument("--sigma-threshold", type=float, default=1.0, + help="AUROC regression must exceed this many pooled Hanley-McNeil SEs " + "to fail the ratchet. Default 1.0.") + args = parser.parse_args() + + if args.loso_results and not args.ratchet_baseline: + parser.error("--ratchet-baseline is required when --loso-results is passed.") + + success = run(args.repo_root, loso_results=args.loso_results, + ratchet_baseline=args.ratchet_baseline, + baselines_file=args.baselines_file, + candidate_reward=args.candidate_reward, + sigma_threshold=args.sigma_threshold) + sys.exit(0 if success else 1) \ No newline at end of file diff --git a/tools/confusion_at_threshold.py b/tools/confusion_at_threshold.py new file mode 100644 index 0000000..4bb153e --- /dev/null +++ b/tools/confusion_at_threshold.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +""" +confusion_at_threshold.py — Team Narnia + +Computes confusion matrices (TP/FP/FN/TN, sensitivity, specificity, PPV) +from a loso_probabilities.csv file at one or more chosen thresholds, using +calibrated_probability (not the uncalibrated model.predict() default that +loso_results.csv's TP/FP/FN/TN columns are built from). + +Built to answer a specific question: does logreg's small-set reward +advantage over XGBoost come from fewer false positives (specificity) or +better true-positive capture, and does it hold at a SHARED threshold or +only at each model's own tuned optimum? + +Usage: + python confusion_at_threshold.py --probs loso_outputs_xgb/loso_probabilities.csv \\ + --probs loso_outputs_logreg_l1/loso_probabilities.csv \\ + --thresholds 0.10,0.12,0.15,0.20 \\ + --per-site + + # If you only have one file and want to sweep it: + python confusion_at_threshold.py --probs loso_outputs_logreg_l1/loso_probabilities.csv \\ + --thresholds 0.05,0.08,0.10,0.12,0.15,0.20 +""" + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd + + +def confusion_at(labels, probs, threshold): + preds = (probs > threshold).astype(int) + tp = int(((preds == 1) & (labels == 1)).sum()) + fp = int(((preds == 1) & (labels == 0)).sum()) + fn = int(((preds == 0) & (labels == 1)).sum()) + tn = int(((preds == 0) & (labels == 0)).sum()) + sensitivity = tp / (tp + fn) if (tp + fn) > 0 else float('nan') + specificity = tn / (tn + fp) if (tn + fp) > 0 else float('nan') + ppv = tp / (tp + fp) if (tp + fp) > 0 else float('nan') + accuracy = (tp + tn) / max(tp + fp + fn + tn, 1) + return dict(TP=tp, FP=fp, FN=fn, TN=tn, + sensitivity=round(sensitivity, 3), specificity=round(specificity, 3), + ppv=round(ppv, 3), accuracy=round(accuracy, 3)) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--probs', action='append', required=True, + help='Path to a loso_probabilities.csv. Repeat for multiple runs to compare.') + parser.add_argument('--thresholds', default='0.05,0.08,0.10,0.12,0.15,0.20') + parser.add_argument('--per-site', action='store_true', + help='Also break out I0006 vs S0001 separately, not just pooled.') + parser.add_argument('--prob-col', default='calibrated_probability', + help='Column to threshold. Use "raw_probability" to compare uncalibrated instead.') + args = parser.parse_args() + + thresholds = [float(t) for t in args.thresholds.split(',')] + + for path in args.probs: + path = Path(path) + df = pd.read_csv(path) + label = f'{path.parent.name}/{path.name}' + print(f'\n{"="*70}\n{label} (n={len(df)}, {int(df.label.sum())} positive)\n{"="*70}') + + labels = df['label'].to_numpy() + probs = df[args.prob_col].to_numpy() + + print(f'\n POOLED (both folds combined):') + for t in thresholds: + c = confusion_at(labels, probs, t) + print(f' t={t:.2f} TP={c["TP"]:>3} FP={c["FP"]:>3} FN={c["FN"]:>3} TN={c["TN"]:>3} ' + f'sens={c["sensitivity"]:.3f} spec={c["specificity"]:.3f} ppv={c["ppv"]:.3f}') + + if args.per_site: + for site in sorted(df['SiteID'].unique()): + sub = df[df['SiteID'] == site] + print(f'\n {site} only (n={len(sub)}, {int(sub.label.sum())} positive):') + for t in thresholds: + c = confusion_at(sub['label'].to_numpy(), sub[args.prob_col].to_numpy(), t) + print(f' t={t:.2f} TP={c["TP"]:>3} FP={c["FP"]:>3} FN={c["FN"]:>3} TN={c["TN"]:>3} ' + f'sens={c["sensitivity"]:.3f} spec={c["specificity"]:.3f} ppv={c["ppv"]:.3f}') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/tools/loso_cv.py b/tools/loso_cv.py new file mode 100644 index 0000000..5734605 --- /dev/null +++ b/tools/loso_cv.py @@ -0,0 +1,806 @@ +#!/usr/bin/env python3 +""" +loso_cv.py — Team Narnia +Leave-One-Site-Out cross-validation for honest local performance estimation. + +RECONCILED 2026-07-09: the local copy of this file had fallen behind every +Kaggle dataset version actually used to produce real results (confirmed via +grep against the real local file -- it had none of --model-family, +--calibration, --tune-hyperparams, --diagnostic-models). This version +reconciles the most complete Kaggle snapshot seen (matching the CONFIRMED +real features/pipeline.py signature: build_logreg_pipeline(y_train, +calibrated, use_age_residual, calibration_ensemble, C, penalty, l1_ratio, +max_iter)) back into the local repo, PLUS adds --penalty/--l1-ratio CLI +passthrough (new, 2026-07-09) for the Ridge/Lasso/ElasticNet regularization +sweep. + +STILL MISSING, confirmed absent from every version seen so far, local or +Kaggle: --rotate-i0002 (3-site LOSO rotation). A run using 3-site rotation +DID happen (loso_summary.txt printed "LOSO rotation: I0006, S0001, I0002 +(3-site)") but the code that produced it hasn't surfaced in any file shared +so far -- holdout_sites is still hardcoded to 2-site here. +""" + +import argparse +import os +import sys +import warnings +import numpy as np +import pandas as pd +from pathlib import Path +from tqdm import tqdm + +warnings.filterwarnings('ignore') + +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +from helper_code import ( + find_patients, load_demographics, load_signal_data, + load_diagnoses, DEMOGRAPHICS_FILE, + PHYSIOLOGICAL_DATA_SUBFOLDER, + ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + HEADERS, +) +from features.demographic import extract_demographic_features +from features.caisr_base import extract_caisr_base_features +from features.caisr_enriched import extract_caisr_enriched_features +from features.physiological_ratios import (extract_physiological_ratio_features, + N_RATIO_FEATURES) +from features import N_CAISR_BASE_FEATURES, N_CAISR_ENRICHED_FEATURES + +_N_BASE = N_CAISR_BASE_FEATURES +_N_ENRICHED = N_CAISR_ENRICHED_FEATURES +_N_RATIO = N_RATIO_FEATURES + +from features.pipeline import build_pipeline, build_logreg_pipeline, PrefitCalibratedModel +from features.age_residuals import AgeResidualizer +from features import IDX_AGE, IDX_CA_RATE, IDX_EEG_VAR_REM_WAKE +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import Pipeline +from sklearn.impute import SimpleImputer +from sklearn.calibration import CalibratedClassifierCV + +from evaluate_model import compute_reward, compute_prevalence + +ENTRY3_THRESHOLDS = [0.05, 0.08, 0.10, 0.12, 0.15, 0.20] + +FEATURE_NAMES = [ + 'Age','Sex_F','Sex_M','Sex_U', + 'Race_As','Race_Bl','Race_Ot','Race_Un','Race_Wh','BMI', + 'AHI_total','Arousal_idx','Limb_idx', + 'W_pct','N1_pct','N2_pct','N3_pct','REM_pct','Sleep_eff', + 'Prob_W','Prob_N3','Prob_arous', + 'OA_rate','CA_rate','HY_rate','RERA_rate', + 'CA_total_ratio','REM_AHI','NREM_AHI','REM_NREM_ratio', + 'N3_gradient','Spont_arousal_idx','N3_entropy', + 'EEG_std_N3_Wake','EEG_mav_N3_Wake','EEG_zcr_N3_Wake', + 'EEG_rms_N3_Wake','EEG_var_N3_Wake','EEG_mob_N3_Wake','EEG_cplx_N3_Wake', + 'EEG_std_REM_Wake','EEG_mav_REM_Wake','EEG_zcr_REM_Wake', + 'EEG_rms_REM_Wake','EEG_var_REM_Wake','EEG_mob_REM_Wake','EEG_cplx_REM_Wake', + 'Chin_atonia_ratio', + 'CA_rate_age_resid','EEG_var_REM_Wake_age_resid', +] + + +def age_conditioned_auroc(labels, probs, ages, gap=2): + idx_pos = np.where(labels == 1)[0] + idx_neg = np.where(labels == 0)[0] + numer = 0.0 + denom = 0 + for i in idx_pos: + for j in idx_neg: + if abs(ages[i] - ages[j]) <= gap: + if probs[i] > probs[j]: + numer += 1.0 + elif probs[i] == probs[j]: + numer += 0.5 + denom += 1 + if denom == 0: + return float('nan'), 0 + return round(numer / denom, 4), denom + + +def standard_auroc(labels, probs): + from scipy.stats import mannwhitneyu + pos = probs[labels == 1] + neg = probs[labels == 0] + if len(pos) == 0 or len(neg) == 0: + return float('nan') + U, _ = mannwhitneyu(pos, neg, alternative='greater') + return round(U / (len(pos) * len(neg)), 4) + + +def extract_all_features(data_folder, verbose=True): + demo_path = os.path.join(data_folder, DEMOGRAPHICS_FILE) + patient_list = find_patients(demo_path) + + features_list = [] + labels_list = [] + ages_list = [] + sites_list = [] + patient_ids = [] + + pbar = tqdm(patient_list, desc='Extracting features', unit='patient', + disable=not verbose) + + for record in pbar: + try: + pid = record[HEADERS['bids_folder']] + site = record[HEADERS['site_id']] + sess = record[HEADERS['session_id']] + + if verbose: + pbar.set_postfix({'patient': pid}) + + pdata = load_demographics(demo_path, pid, sess) + demo_f = extract_demographic_features(pdata) + age = float(demo_f[0]) + bdsp_id = pdata.get('BDSPPatientID', pid) + + phys_file = os.path.join(data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, + site, f'{pid}_ses-{sess}.edf') + if not os.path.exists(phys_file): + continue + phys_data, phys_fs = load_signal_data(phys_file) + + algo_file = os.path.join(data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + site, f'{pid}_ses-{sess}_caisr_annotations.edf') + if os.path.exists(algo_file): + algo_data, _ = load_signal_data(algo_file) + base_f = extract_caisr_base_features(algo_data) + enriched_f = extract_caisr_enriched_features(algo_data) + ratio_f = extract_physiological_ratio_features( + phys_data, phys_fs, algo_data) + del algo_data + else: + base_f = np.full(_N_BASE, float('nan')) + enriched_f = np.full(_N_ENRICHED, float('nan')) + ratio_f = np.full(_N_RATIO, float('nan')) + + del phys_data + + label = load_diagnoses(demo_path, pid) + if label not in (0, 1): + continue + + features_list.append( + np.hstack([demo_f, base_f, enriched_f, ratio_f])) + labels_list.append(label) + ages_list.append(age) + sites_list.append(site) + patient_ids.append(str(bdsp_id)) + + del base_f, enriched_f, ratio_f + + except Exception as ex: + tqdm.write(f' ! Error on {pid}: {ex}') + continue + + pbar.close() + + return ( + np.asarray(features_list, dtype=np.float32), + np.asarray(labels_list, dtype=int), + np.asarray(ages_list, dtype=float), + np.asarray(sites_list), + np.asarray(patient_ids), + ) + + +def build_model(X_train, y_train, use_age_residual=True, model_family='xgboost', + C=0.01, penalty=None, l1_ratio=None): + """Uncalibrated model. + + C/penalty/l1_ratio (2026-07-09): only used when model_family='logreg'. + Default C=0.01, penalty=None reproduces the original tested config. + """ + if model_family == 'logreg': + model = build_logreg_pipeline(y_train, calibrated=False, use_age_residual=use_age_residual, + C=C, penalty=penalty, l1_ratio=l1_ratio) + elif model_family == 'catboost': + model = build_catboost_pipeline(y_train, calibrated=False, use_age_residual=use_age_residual) + else: + model = build_pipeline(y_train, calibrated=False, use_age_residual=use_age_residual) + model.fit(X_train, y_train) + return model + + +def build_calibrated_model(X_train, y_train, use_age_residual=True, + calibration_mode='cv5', model_family='xgboost', + C=0.01, penalty=None, l1_ratio=None): + if calibration_mode == 'prefit': + if model_family != 'xgboost': + raise ValueError('calibration_mode="prefit" is XGBoost-specific ' + '(PrefitCalibratedModel) — not supported with ' + f'model_family="{model_family}".') + model = PrefitCalibratedModel(use_age_residual=use_age_residual) + model.fit(X_train, y_train) + return model + + calibration_ensemble = (calibration_mode != 'cv5-single-curve') + if model_family == 'logreg': + model = build_logreg_pipeline(y_train, calibrated=True, use_age_residual=use_age_residual, + calibration_ensemble=calibration_ensemble, + C=C, penalty=penalty, l1_ratio=l1_ratio) + elif model_family == 'catboost': + model = build_catboost_pipeline(y_train, calibrated=True, use_age_residual=use_age_residual, + calibration_ensemble=calibration_ensemble) + else: + model = build_pipeline(y_train, calibrated=True, use_age_residual=use_age_residual, + calibration_ensemble=calibration_ensemble) + model.fit(X_train, y_train) + return model + + +def sweep_thresholds(probs, labels, ages, age_to_prevalence, + thresholds=ENTRY3_THRESHOLDS, label=''): + print(f'\n Threshold sweep {label}:') + rows = [] + for t in thresholds: + preds = (probs > t).astype(int) + reward = compute_reward(labels, preds, ages, age_to_prevalence) + auroc, _ = age_conditioned_auroc(labels, probs, ages) + rows.append({'threshold': t, 'reward': round(reward, 4), 'auroc': auroc}) + print(f' t={t:.2f} reward={reward:.4f} auroc={auroc:.4f}') + return rows + + +def run_fold(holdout_site, features, labels, ages, sites, patient_ids, + age_to_prevalence, fold_name, use_age_residual=True, + calibration_mode='cv5', model_family='xgboost', + C=0.01, penalty=None, l1_ratio=None): + train_mask = (sites != holdout_site) + test_mask = (sites == holdout_site) + + X_train = features[train_mask] + y_train = labels[train_mask] + X_test = features[test_mask] + y_test = labels[test_mask] + ages_test = ages[test_mask] + ids_test = patient_ids[test_mask] + sites_test = sites[test_mask] + + n_pos_train = int(y_train.sum()) + n_neg_train = int((y_train == 0).sum()) + n_pos_test = int(y_test.sum()) + n_neg_test = int((y_test == 0).sum()) + + print(f'\n Train: {len(y_train)} patients ' + f'({n_pos_train} pos, {n_neg_train} neg)') + print(f' Test: {len(y_test)} patients ' + f'({n_pos_test} pos, {n_neg_test} neg)') + + if n_pos_test < 3: + print(f' SKIP — fewer than 3 positives in holdout site.') + return None, [] + + print(' Fitting model (uncalibrated) ...') + model = build_model(X_train, y_train, use_age_residual=use_age_residual, + model_family=model_family, C=C, penalty=penalty, l1_ratio=l1_ratio) + + probs = model.predict_proba(X_test)[:, 1] + preds = model.predict(X_test) + + print(f' Fitting model (Platt-calibrated, mode={calibration_mode}, family={model_family}) ...') + calibrated_model = build_calibrated_model(X_train, y_train, + use_age_residual=use_age_residual, + calibration_mode=calibration_mode, + model_family=model_family, + C=C, penalty=penalty, l1_ratio=l1_ratio) + calibrated_probs = calibrated_model.predict_proba(X_test)[:, 1] + + age_auroc, n_pairs = age_conditioned_auroc(y_test, probs, ages_test) + std_auroc = standard_auroc(y_test, probs) + accuracy = float((preds == y_test).mean()) + tp = int(((preds == 1) & (y_test == 1)).sum()) + fp = int(((preds == 1) & (y_test == 0)).sum()) + fn = int(((preds == 0) & (y_test == 1)).sum()) + tn = int(((preds == 0) & (y_test == 0)).sum()) + sensitivity = tp / (tp + fn) if (tp + fn) > 0 else float('nan') + specificity = tn / (tn + fp) if (tn + fp) > 0 else float('nan') + ppv = tp / (tp + fp) if (tp + fp) > 0 else float('nan') + + sweep_thresholds(calibrated_probs, y_test, ages_test, age_to_prevalence, + label=f'(holdout={holdout_site}, calibrated)') + + classifier_step = model.named_steps['classifier'] + if hasattr(classifier_step, 'feature_importances_'): + imp = classifier_step.feature_importances_ + elif hasattr(classifier_step, 'coef_'): + imp = np.abs(classifier_step.coef_[0]) + else: + imp = np.zeros(features.shape[1]) + top5_idx = np.argsort(imp)[::-1][:5] + top5_feat = [(FEATURE_NAMES[i] if i < len(FEATURE_NAMES) else f'feat_{i}', + round(float(imp[i]), 4)) for i in top5_idx] + + result = { + 'fold': fold_name, + 'holdout_site': holdout_site, + 'n_train': len(y_train), + 'n_pos_train': n_pos_train, + 'n_neg_train': n_neg_train, + 'n_test': len(y_test), + 'n_pos_test': n_pos_test, + 'n_neg_test': n_neg_test, + 'age_cond_auroc': age_auroc, + 'n_age_pairs': n_pairs, + 'std_auroc': std_auroc, + 'accuracy': round(accuracy, 4), + 'sensitivity': round(sensitivity, 4), + 'specificity': round(specificity, 4), + 'ppv': round(ppv, 4), + 'TP': tp, 'FP': fp, 'FN': fn, 'TN': tn, + 'top5_features': str(top5_feat), + } + + print(f' Age-conditioned AUROC : {age_auroc:.4f} ({n_pairs} pairs)') + print(f' Standard AUROC : {std_auroc:.4f}') + print(f' Sensitivity : {sensitivity:.3f}') + print(f' Specificity : {specificity:.3f}') + print(f' Top features: {top5_feat[:3]}') + + per_patient_rows = [] + for i in range(len(ids_test)): + per_patient_rows.append({ + 'BDSPPatientID': ids_test[i], + 'SiteID': sites_test[i], + 'Age': ages_test[i], + 'label': int(y_test[i]), + 'raw_probability': round(float(probs[i]), 6), + 'calibrated_probability': round(float(calibrated_probs[i]), 6), + }) + + return result, per_patient_rows + + +def run_ablation(holdout_site, features, labels, ages, sites): + # Derived from features/__init__.py's IDX_* constants rather than + # hardcoded — this exact dict went stale silently when caisr_enriched.py + # grew from 11 to 14 features (2026-07-20 limb enrichment) and would + # have kept reporting a group ablation against the WRONG 11 columns + # with no error raised. Same class of bug pipeline.py's header + # documents for loso_cv.py's old hand-copied Pipeline() — one + # canonical source, not a second hand-maintained copy. + from features import (IDX_DEMO_START, IDX_DEMO_END, IDX_BASE_START, IDX_BASE_END, + IDX_ENRICHED_START, IDX_ENRICHED_END, IDX_RATIO_START, IDX_RATIO_END) + GROUPS = { + 'demo': list(range(IDX_DEMO_START, IDX_DEMO_END)), + 'caisr_base': list(range(IDX_BASE_START, IDX_BASE_END)), + 'caisr_enriched': list(range(IDX_ENRICHED_START, IDX_ENRICHED_END)), + 'ratios': list(range(IDX_RATIO_START, IDX_RATIO_END)), + } + + train_mask = sites != holdout_site + test_mask = sites == holdout_site + ages_test = ages[test_mask] + y_test = labels[test_mask] + + if y_test.sum() < 3: + return [] + + baseline_model = build_model(features[train_mask], labels[train_mask]) + baseline_probs = baseline_model.predict_proba(features[test_mask])[:, 1] + baseline_auroc, _ = age_conditioned_auroc(y_test, baseline_probs, ages_test) + + rows = [] + for group_name, indices in GROUPS.items(): + X_ablated = features.copy() + X_ablated[:, indices] = float('nan') + + ablated_model = build_model(X_ablated[train_mask], labels[train_mask]) + ablated_probs = ablated_model.predict_proba(X_ablated[test_mask])[:, 1] + ablated_auroc, _ = age_conditioned_auroc(y_test, ablated_probs, ages_test) + + delta = round(ablated_auroc - baseline_auroc, 4) + rows.append({ + 'holdout_site': holdout_site, + 'ablated_group': group_name, + 'n_features': len(indices), + 'baseline_auroc': baseline_auroc, + 'ablated_auroc': ablated_auroc, + 'delta': delta, + 'interpretation': ( + 'HURTS generalisation' if delta < -0.02 else + 'HELPS generalisation' if delta > 0.02 else + 'NEUTRAL' + ), + }) + print(f' Ablate {group_name:<18s}: AUROC {ablated_auroc:.4f} ' + f'(delta={delta:+.4f})') + + return rows + + +HYPERPARAM_CANDIDATES = { + 'baseline': {}, + 'more_l2': {'reg_lambda': 10.0}, + 'more_min_child_wt': {'min_child_weight': 25}, + 'shallower': {'max_depth': 2}, + 'fewer_trees': {'n_estimators': 100}, + 'conservative_combo': {'max_depth': 2, 'reg_lambda': 8.0, + 'min_child_weight': 25, 'n_estimators': 150}, +} + + +def run_hyperparameter_sweep(holdout_site, features, labels, ages, sites, use_age_residual=True): + train_mask = sites != holdout_site + test_mask = sites == holdout_site + X_train, y_train = features[train_mask], labels[train_mask] + X_test, y_test = features[test_mask], labels[test_mask] + ages_test = ages[test_mask] + + if int(y_test.sum()) < 3: + print(f' SKIP — fewer than 3 positives in holdout site {holdout_site}.') + return [] + + rows = [] + for name, overrides in HYPERPARAM_CANDIDATES.items(): + model = build_pipeline(y_train, calibrated=False, + use_age_residual=use_age_residual, + xgb_overrides=overrides) + model.fit(X_train, y_train) + raw_probs = model.predict_proba(X_test)[:, 1] + + auroc, n_pairs = age_conditioned_auroc(y_test, raw_probs, ages_test) + row = { + 'holdout_site': holdout_site, + 'candidate': name, + 'overrides': str(overrides), + 'raw_median': round(float(np.median(raw_probs)), 4), + 'raw_frac_above_0.5': round(float((raw_probs > 0.5).mean()), 4), + 'raw_frac_above_0.12': round(float((raw_probs > 0.12).mean()), 4), + 'age_cond_auroc': auroc, + 'n_age_pairs': n_pairs, + } + rows.append(row) + print(f" {name:20s} raw_median={row['raw_median']:.4f} " + f"frac>0.5={row['raw_frac_above_0.5']:.3f} " + f"AUROC={auroc}") + + return rows + + +DIAGNOSTIC_MODELS = { + 'xgboost_baseline': {'family': 'xgboost', 'params': {}}, + 'logreg_strong_l2': {'family': 'logreg', 'params': {'C': 0.01}}, + 'logreg_moderate_l2':{'family': 'logreg', 'params': {'C': 1.0}}, +} + + +def _build_diagnostic_pipeline(y_train, family, params, use_age_residual=True): + if family == 'xgboost': + return build_pipeline(y_train, calibrated=False, + use_age_residual=use_age_residual, + xgb_overrides=params) + + if family == 'logreg': + steps = [] + if use_age_residual: + steps.append(('age_residual', AgeResidualizer( + age_idx=IDX_AGE, ca_rate_idx=IDX_CA_RATE, + eeg_var_rem_wake_idx=IDX_EEG_VAR_REM_WAKE))) + steps.append(('imputer', SimpleImputer(strategy='median'))) + steps.append(('scaler', StandardScaler())) + steps.append(('classifier', LogisticRegression( + penalty='l2', class_weight='balanced', max_iter=2000, + random_state=42, **params))) + return Pipeline(steps) + + raise ValueError(f'Unknown model family: {family}') + + +def run_model_family_diagnostic(holdout_site, features, labels, ages, sites, + use_age_residual=True): + train_mask = sites != holdout_site + test_mask = sites == holdout_site + X_train, y_train = features[train_mask], labels[train_mask] + X_test, y_test = features[test_mask], labels[test_mask] + ages_test = ages[test_mask] + + if int(y_test.sum()) < 3: + print(f' SKIP — fewer than 3 positives in holdout site {holdout_site}.') + return [] + + rows = [] + for name, spec in DIAGNOSTIC_MODELS.items(): + model = _build_diagnostic_pipeline(y_train, spec['family'], spec['params'], + use_age_residual=use_age_residual) + model.fit(X_train, y_train) + raw_probs = model.predict_proba(X_test)[:, 1] + + auroc, n_pairs = age_conditioned_auroc(y_test, raw_probs, ages_test) + row = { + 'holdout_site': holdout_site, + 'model': name, + 'family': spec['family'], + 'params': str(spec['params']), + 'raw_median': round(float(np.median(raw_probs)), 4), + 'raw_frac_above_0.5': round(float((raw_probs > 0.5).mean()), 4), + 'raw_frac_above_0.12': round(float((raw_probs > 0.12).mean()), 4), + 'age_cond_auroc': auroc, + 'n_age_pairs': n_pairs, + } + rows.append(row) + print(f" {name:20s} raw_median={row['raw_median']:.4f} " + f"frac>0.5={row['raw_frac_above_0.5']:.3f} " + f"AUROC={auroc}") + + return rows + + +CATBOOST_PARAMS = dict( + iterations=300, + depth=3, + learning_rate=0.05, + l2_leaf_reg=2.0, + random_seed=42, + verbose=0, +) + + +def _catboost_scale_pos_weight(y_train): + y_train = np.asarray(y_train) + n_pos = int((y_train == 1).sum()) + n_neg = int((y_train == 0).sum()) + return n_neg / n_pos if n_pos > 0 else 1.0 + + +def build_catboost_pipeline(y_train, calibrated=True, use_age_residual=True, + calibration_ensemble=True): + try: + from catboost import CatBoostClassifier + except ImportError as e: + raise ImportError( + "model_family='catboost' requires the catboost package, which " + "isn't installed. Run: pip install catboost" + ) from e + + spw = _catboost_scale_pos_weight(y_train) + cb = CatBoostClassifier(scale_pos_weight=spw, **CATBOOST_PARAMS) + + classifier = ( + CalibratedClassifierCV(cb, method='sigmoid', cv=5, ensemble=calibration_ensemble) + if calibrated else cb + ) + + steps = [] + if use_age_residual: + steps.append(('age_residual', AgeResidualizer( + age_idx=IDX_AGE, + ca_rate_idx=IDX_CA_RATE, + eeg_var_rem_wake_idx=IDX_EEG_VAR_REM_WAKE, + ))) + steps.append(('imputer', SimpleImputer(strategy='median'))) + steps.append(('classifier', classifier)) + + return Pipeline(steps) + + +def parse_args(): + p = argparse.ArgumentParser(description='LOSO CV — Narnia_ML') + p.add_argument('--data', required=True, help='Training data folder') + p.add_argument('--out', required=True, help='Output directory') + p.add_argument('--ablation', action='store_true', + help='Run feature group ablation (adds ~2x runtime)') + p.add_argument('--no-age-residual', action='store_true', + help='Build the Entry-3-EQUIVALENT pipeline (no AgeResidualizer).') + p.add_argument('--calibration', choices=['cv5', 'cv5-single-curve', 'prefit'], default='cv5', + help='Calibration mode. See build_calibrated_model() docstring.') + p.add_argument('--tune-hyperparams', action='store_true', + help='Run the HYPERPARAM_CANDIDATES sweep instead of normal fold evaluation.') + p.add_argument('--diagnostic-models', action='store_true', + help='Run the DIAGNOSTIC_MODELS comparison instead of normal fold evaluation.') + p.add_argument('--model-family', choices=['xgboost', 'logreg', 'catboost'], default='xgboost', + help='Base model family.') + p.add_argument('--C', type=float, default=0.01, + help=('Inverse regularization strength for --model-family logreg. ' + 'Default 0.01 matches build_logreg_pipeline()\'s tested default. ' + 'Ignored for other model families. Added 2026-07-09.')) + p.add_argument('--penalty', choices=['l1', 'l2', 'elasticnet'], default=None, + help=('Regularization penalty for --model-family logreg. Default ' + '(unset) uses LOGREG_PARAMS\' l2 -- the exact original tested ' + 'config. l1/elasticnet are UNTESTED territory as of this flag ' + '(2026-07-09). Ignored for other model families.')) + p.add_argument('--l1-ratio', type=float, default=None, + help='Required if --penalty elasticnet. Ignored otherwise.') + p.add_argument('--holdout-sites', default='I0006,S0001', + help=('Comma-separated list of sites to rotate as LOSO holdout. ' + 'Default "I0006,S0001" matches all prior runs (I0002 always ' + 'in training). Added 2026-07-09 for cheap single-fold ' + 'backfills, e.g. --holdout-sites I0002 to run ONLY the ' + 'I0002-held-out fold (fast — no need to also re-run the ' + 'I0006/S0001 folds, which are unaffected by whether I0002 ' + 'itself is ever held out). Pass "I0006,S0001,I0002" for ' + 'full 3-site rotation.')) + p.add_argument('--verbose', action='store_true', default=True) + return p.parse_args() + + +def main(): + args = parse_args() + out_dir = Path(args.out) + out_dir.mkdir(parents=True, exist_ok=True) + use_age_residual = not args.no_age_residual + + if args.penalty and args.model_family != 'logreg': + print(f'!! WARNING: --penalty was passed but --model-family is ' + f'"{args.model_family}", not "logreg" -- --penalty/--C/--l1-ratio ' + f'are silently ignored for this model family.\n') + + print(f'\n{"="*60}') + print(' Leave-One-Site-Out CV — Team Narnia') + print(f'{"="*60}\n') + print(f' Pipeline mode: {"Entry 4 (WITH age-residual features, 50-wide)" if use_age_residual else "Entry 3-EQUIVALENT (NO age-residual features, 48-wide)"}') + print(f' Calibration mode: {args.calibration}') + print(f' Model family: {args.model_family}') + if args.model_family == 'logreg': + print(f' Logreg C={args.C} penalty={args.penalty or "l2 (default)"}' + f'{f" l1_ratio={args.l1_ratio}" if args.penalty == "elasticnet" else ""}') + print(f'{"="*60}\n') + + print('Step 1: Extracting features from all patients (~16-18 min) ...\n') + features, labels, ages, sites, patient_ids = extract_all_features( + args.data, verbose=args.verbose) + + print(f'\nFeature matrix: {features.shape}') + print(f'Labels: {labels.sum()} positive, {(labels==0).sum()} negative') + print(f'Sites: {dict(zip(*np.unique(sites, return_counts=True)))}') + + age_to_prevalence = compute_prevalence(ages, labels, ages, gap=2) + + print('\nSite breakdown:') + for site in np.unique(sites): + mask = sites == site + n_pos = int(labels[mask].sum()) + n_neg = int((labels[mask] == 0).sum()) + print(f' {site}: {mask.sum()} patients ' + f'({n_pos} pos / {n_neg} neg) ' + f'CI rate={n_pos/mask.sum():.1%}') + + if args.tune_hyperparams: + print(f'\n{"="*60}') + print(' Hyperparameter sweep mode (--tune-hyperparams)') + print(f'{"="*60}\n') + print(f' Candidates: {list(HYPERPARAM_CANDIDATES.keys())}\n') + + sweep_rows = [] + for holdout in ['I0006', 'S0001']: + print(f'\n{"─"*50}') + print(f' Sweep: hold out {holdout}') + print(f'{"─"*50}') + rows = run_hyperparameter_sweep(holdout, features, labels, ages, sites, + use_age_residual=use_age_residual) + sweep_rows.extend(rows) + + sweep_df = pd.DataFrame(sweep_rows) + sweep_df.to_csv(out_dir / 'loso_hyperparam_sweep.csv', index=False) + print(f'\nWrote {out_dir / "loso_hyperparam_sweep.csv"}') + return + + if args.diagnostic_models: + print(f'\n{"="*60}') + print(' Model-family diagnostic mode (--diagnostic-models)') + print(f'{"="*60}\n') + print(f' Models: {list(DIAGNOSTIC_MODELS.keys())}\n') + + diag_rows = [] + for holdout in ['I0006', 'S0001']: + print(f'\n{"─"*50}') + print(f' Diagnostic: hold out {holdout}') + print(f'{"─"*50}') + rows = run_model_family_diagnostic(holdout, features, labels, ages, sites, + use_age_residual=use_age_residual) + diag_rows.extend(rows) + + diag_df = pd.DataFrame(diag_rows) + diag_df.to_csv(out_dir / 'loso_model_diagnostic.csv', index=False) + print(f'\nWrote {out_dir / "loso_model_diagnostic.csv"}') + return + + # NOTE (2026-07-09): was hardcoded to ['I0006', 'S0001']. Now configurable + # via --holdout-sites, so a single-fold backfill (e.g. just I0002, to + # complete the model-family x site matrix without re-running the full + # 2-fold cycle) doesn't require a separate script. + holdout_sites = [s.strip() for s in args.holdout_sites.split(',') if s.strip()] + print(f' Holdout sites: {holdout_sites}') + + print('\nStep 2: Running LOSO folds ...') + fold_results = [] + ablation_rows = [] + all_patient_rows = [] + + for holdout in holdout_sites: + fold_name = f'holdout_{holdout}' + print(f'\n{"─"*50}') + print(f' Fold: hold out {holdout}') + print(f'{"─"*50}') + + result, patient_rows = run_fold(holdout, features, labels, ages, sites, + patient_ids, age_to_prevalence, fold_name, + use_age_residual=use_age_residual, + calibration_mode=args.calibration, + model_family=args.model_family, + C=args.C, penalty=args.penalty, l1_ratio=args.l1_ratio) + if result: + fold_results.append(result) + all_patient_rows.extend(patient_rows) + + if args.ablation: + print(f'\n Feature group ablation (holdout={holdout}):') + ab = run_ablation(holdout, features, labels, ages, sites) + ablation_rows.extend(ab) + + results_df = pd.DataFrame(fold_results) + results_df.to_csv(out_dir / 'loso_results.csv', index=False) + + if ablation_rows: + ablation_df = pd.DataFrame(ablation_rows) + ablation_df.to_csv(out_dir / 'loso_feature_impact.csv', index=False) + + pooled_sweep_rows = [] + if all_patient_rows: + patients_df = pd.DataFrame(all_patient_rows) + patients_df.to_csv(out_dir / 'loso_probabilities.csv', index=False) + print(f'\nWrote per-patient probabilities: ' + f'{out_dir / "loso_probabilities.csv"} ' + f'({len(patients_df)} patients)') + + pooled_probs = patients_df['calibrated_probability'].to_numpy() + pooled_labels = patients_df['label'].to_numpy() + pooled_ages = patients_df['Age'].to_numpy() + pooled_sweep_rows = sweep_thresholds( + pooled_probs, pooled_labels, pooled_ages, age_to_prevalence, + label='(POOLED — both folds combined, use this to pick submission threshold)') + pd.DataFrame(pooled_sweep_rows).to_csv( + out_dir / 'loso_threshold_sweep.csv', index=False) + + lines = [ + '=' * 60, + ' LOSO CV SUMMARY — Team Narnia', + '=' * 60, + '', + f' Pipeline mode: {"Entry 4 (WITH age-residual features, 50-wide)" if use_age_residual else "Entry 3-EQUIVALENT (NO age-residual features, 48-wide)"}', + f' Calibration mode: {args.calibration}', + f' Model family: {args.model_family}', + ] + if args.model_family == 'logreg': + lines.append(f' Logreg C={args.C} penalty={args.penalty or "l2 (default)"}') + lines.append('') + + if len(fold_results) > 0: + aurocs = [r['age_cond_auroc'] for r in fold_results + if not np.isnan(r['age_cond_auroc'])] + mean_auroc = np.mean(aurocs) if aurocs else float('nan') + + lines += ['── Per-fold results ──────────────────────────────────────', ''] + for r in fold_results: + lines.append( + f" Holdout {r['holdout_site']:<8s} " + f"Age-cond AUROC={r['age_cond_auroc']:.4f} " + f"Std AUROC={r['std_auroc']:.4f} " + f"Sens={r['sensitivity']:.3f} Spec={r['specificity']:.3f} " + f"n_pos={r['n_pos_test']}" + ) + lines += ['', f' Mean age-conditioned AUROC: {mean_auroc:.4f}', ''] + + if pooled_sweep_rows: + lines += ['── Pooled threshold sweep (both folds combined) ──', ''] + best_reward_row = max(pooled_sweep_rows, key=lambda r: r['reward']) + for row in pooled_sweep_rows: + marker = ' <-- highest reward' if row is best_reward_row else '' + lines.append( + f" t={row['threshold']:.2f} reward={row['reward']:.4f} " + f"auroc={row['auroc']:.4f}{marker}" + ) + lines.append('') + + summary = '\n'.join(lines) + print('\n' + summary) + (out_dir / 'loso_summary.txt').write_text(summary) + print(f'\nOutputs written to: {out_dir}') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/tools/ratchet_baselines.json b/tools/ratchet_baselines.json new file mode 100644 index 0000000..bb3853b --- /dev/null +++ b/tools/ratchet_baselines.json @@ -0,0 +1,88 @@ +{ + "_readme": "Hand-curated ratchet baselines. Do NOT let a script auto-overwrite these — same discipline as learning_log.md's decision log. Only promote a new run to a baseline here after you've deliberately decided to trust it, not automatically after every loso_cv.py run. Schema (2026-07-08, corrected): 'folds' is a list of per-holdout-site {site, age_cond_auroc, n_pos, n_neg}. age_cond_auroc/reward at the top level are DERIVED (simple mean across folds) and shown for convenience/display only — the ratchet check recomputes from 'folds' directly, matching loso_cv.py's own 'Mean age-conditioned AUROC' convention (a plain average across folds, NOT weighted by fold size). This was corrected 2026-07-08 after a size-weighted version gave a materially different, wrong number (0.5411 vs the confirmed-correct 0.5443) on a real CatBoost run.\n\n'leaderboard_confirmed' (added 2026-07-08): whether THIS baseline's config has ever actually been submitted to the real Challenge leaderboard, as opposed to only evaluated locally via LOSO. This matters because the LOSO->leaderboard offset has already been shown to be unstable at small scale (Entry 1: LOSO 0.499 -> leaderboard 0.624, delta +0.125; Entry 2: LOSO 0.540 -> leaderboard 0.616, delta +0.076 — different sign of surprise each time, offset retired as a predictor 2026-07-02) and there is currently ZERO leaderboard data point at large scale — no large-set config has been submitted yet. A ratchet PASS against a leaderboard_confirmed=false baseline means 'at least as good as our best LOSO estimate so far', which is a materially weaker claim than 'at least as good as something we know performs on the real leaderboard.' See ratchet_check.py's printed CAUTION for how this is surfaced.", + + "small_entry3": { + "description": "Entry 3 config: 48 features, XGBoost (depth=3, reg_alpha=0.1, reg_lambda=2.0, min_child_weight=5), CalibratedClassifierCV(sigmoid, cv=5), THRESHOLD=0.12. Small training set (1,103 patients).", + "scale": "small", + "model_family": "xgboost", + "feature_count": 48, + "folds": [ + {"site": "I0006", "age_cond_auroc": 0.5374, "n_pos": 20, "n_neg": 172}, + {"site": "S0001", "age_cond_auroc": 0.5425, "n_pos": 56, "n_neg": 801} + ], + "age_cond_auroc": 0.5399, + "reward": 0.1148, + "reward_threshold": 0.12, + "date_confirmed": "2026-07-01", + "source": "learning_log.md, 2026-07-01 (later still) — Entry 3 LOSO results", + "leaderboard_confirmed": true, + "leaderboard_auroc": 0.606, + "leaderboard_reward": null, + "leaderboard_standard_auroc": 0.772, + "leaderboard_accuracy": 0.803, + "leaderboard_note": "Entry 3 WAS submitted (3rd of 3 official entries: 0.624 -> 0.616 -> 0.606 age-cond AUROC). Exact reward figure still not on record, but per features/pipeline.py's PrefitCalibratedModel docstring (uploaded 2026-07-09): going Entry 2 -> Entry 3 (calibration added, nothing else changed), standard AUROC ROSE (0.746 -> 0.772) while age-conditioned AUROC AND REWARD BOTH FELL, and accuracy dropped sharply (0.866 -> 0.803) with F-measure barely moving — more patients called positive, at a cost. This corrects an earlier, more optimistic assumption in this project's history that calibration would likely also improve reward at Entry 3; it did not on the real leaderboard." + }, + + "large_logreg": { + "description": "Logistic regression, large training set (6,281 patients), cv5 calibration, age-residual features on. First result to clear the ~0.57-0.58 noise floor on BOTH folds independently. Fold n_pos/n_neg backfilled 2026-07-08 from the CatBoost run's loso_results.csv — same LOSO split, sizes are model-independent.", + "scale": "large", + "model_family": "logreg", + "feature_count": 50, + "folds": [ + {"site": "I0006", "age_cond_auroc": 0.6154, "n_pos": 112, "n_neg": 1030}, + {"site": "S0001", "age_cond_auroc": 0.5850, "n_pos": 334, "n_neg": 4805} + ], + "age_cond_auroc": 0.6002, + "reward": 0.0719, + "reward_threshold": 0.10, + "date_confirmed": "2026-07-07", + "source": "learning_log.md, 2026-07-07 — Calibrated logistic regression entry; fold sizes confirmed 2026-07-08", + "leaderboard_confirmed": false, + "leaderboard_auroc": null, + "leaderboard_reward": null, + "leaderboard_note": "NOT YET SUBMITTED. No large-training-set config of any kind has been submitted to the real leaderboard as of 2026-07-08 — this is a LOSO-only estimate. Treat any ratchet PASS against this baseline as 'not worse than our best local estimate', not as a leaderboard guarantee." + }, + + "entry5_3site": { + "description": "Entry 5 config as ACTUALLY SHIPPED: logreg, C=0.01 (l2, LOGREG_PARAMS default), CalibratedClassifierCV(sigmoid, cv=5), AgeResidualizer on, NO sample weighting, THRESHOLD=0.10. Large training set, 3-site LOSO rotation (I0006/S0001/I0002, via --rotate-i0002). This is the correct baseline for any Entry 6+ ratchet check — do NOT use 'large_logreg' for that purpose, which is a 2-site-rotation, pre-Entry-5 config and is not a fair comparison (different fold composition, independent of any model change).", + "scale": "large", + "model_family": "logreg", + "feature_count": 50, + "folds": [ + {"site": "I0006", "age_cond_auroc": 0.6184, "n_pos": 112, "n_neg": 1030}, + {"site": "S0001", "age_cond_auroc": 0.5905, "n_pos": 333, "n_neg": 4805}, + {"site": "I0002", "age_cond_auroc": 0.7259, "n_pos": 52, "n_neg": 267} + ], + "age_cond_auroc": 0.6449, + "reward": 0.0882, + "reward_threshold": 0.10, + "date_confirmed": "2026-07-15", + "source": "learning_log_p2.md, 2026-07-15 (L1/L2 grid reproducibility-floor entry) — this run's own fresh l2/C=0.01 re-run via reg_sweep.py, explicitly established there as the reference point for every subsequent branch (b) comparison in the C=0.001/alpha=1.0 chain (0.6449 -> 0.6485 -> 0.6459). Supersedes the original 2026-07-08 loso_cv.py numbers (I0006 0.6154/S0001 0.5850/I0002 0.7328, mean 0.6444) — same config, ~0.005-0.007 AUROC run-to-run noise floor established in that same entry, not a real difference. Fold n_pos/n_neg cross-checked against the uploaded alpha1_reproduction/loso_results.csv (2026-07-16 final reproduction run) — identical data split (post-denylist); note S0001 n_pos=333 here vs. 334 in the older 2026-07-08 log, consistent with the 1-patient denylist exclusion applied in later runs.", + "leaderboard_confirmed": true, + "leaderboard_auroc": 0.660, + "leaderboard_reward": 0.101, + "leaderboard_standard_auroc": 0.826, + "leaderboard_accuracy": 0.869, + "leaderboard_note": "Entry 5 WAS submitted (submission 2156, commit 8060971c9195bdaa2fbd4127b5a9c9a2be18e961): reward 0.101, age-conditioned AUROC 0.660, age-weighted AUROC 0.664, standard AUROC 0.826, AUPRC 0.303, accuracy 0.869, F-measure 0.372. Leaderboard age-cond AUROC (0.660) exceeds this LOSO estimate (0.6449) by +0.011 — consistent with, though smaller than, this project's established LOSO-always-underestimates-leaderboard pattern (prior gaps: +0.125 Entry 1, +0.076 Entry 2, +0.060 Entry 5 vs. an earlier LOSO estimate)." + }, + + "large_xgboost_age_resid": { + "description": "XGBoost, large training set, cv5 calibration, WITH age-residual features (Entry 4 config, 50 features). Reference point for the reward-collapse investigation — NOT a recommended submission candidate given logreg's AUROC advantage at no reward cost.", + "scale": "large", + "model_family": "xgboost", + "feature_count": 50, + "folds": [ + {"site": "I0006", "age_cond_auroc": 0.5485, "n_pos": 112, "n_neg": 1030}, + {"site": "S0001", "age_cond_auroc": 0.5392, "n_pos": 334, "n_neg": 4805} + ], + "age_cond_auroc": 0.5438, + "reward": 0.0726, + "reward_threshold": 0.12, + "date_confirmed": "2026-07-05", + "source": "learning_log.md, 2026-07-05 / 2026-07-06 — large-set XGBoost runs; fold sizes confirmed 2026-07-08", + "leaderboard_confirmed": false, + "leaderboard_auroc": null, + "leaderboard_reward": null, + "leaderboard_note": "NOT YET SUBMITTED. Same caveat as large_logreg — no large-set config has ever been submitted." + } +} \ No newline at end of file diff --git a/tools/ratchet_check.py b/tools/ratchet_check.py new file mode 100644 index 0000000..e702eaf --- /dev/null +++ b/tools/ratchet_check.py @@ -0,0 +1,376 @@ +#!/usr/bin/env python + +""" +ratchet_check.py +Team Narnia — PhysioNet Challenge 2026 + +Answers one question before you submit: "is this LOSO result actually worse +than what we've already confirmed, or is it inside the noise band we already +characterized?" Motivated directly by the Entry 1-3 leaderboard sequence +(0.624 -> 0.616 -> 0.606), which LOOKED like a monotonic decline but was +confirmed via Hanley-McNeil SE to be ~0.12-0.40 sigma per step — i.e. noise, +not regression. This script automates that same check instead of re-deriving +it by hand every time. + +One-directional ratchet: only ever fails on a REGRESSION beyond the sigma +threshold. An improvement of any size always passes — this is a gate against +backsliding, not a two-sided significance test. + +Baselines live in ratchet_baselines.json, hand-curated (see that file's +_readme). This script never writes to it — promoting a new result to a +baseline is a deliberate decision, not something a script should do for you. + +Usage (standalone): + python ratchet_check.py --loso-results loso_results.csv --baseline small_entry3 + + # If your loso_results.csv uses different column names than the + # defaults below, point at them explicitly: + python ratchet_check.py --loso-results loso_results.csv --baseline small_entry3 \\ + --auroc-col age_cond_auroc --n-pos-col n_pos_test --n-neg-col n_neg_test + +Also importable — see check_ratchet() for use from check_submission_files.py. +""" + +import argparse +import json +import os +import sys +from pathlib import Path + +_DEFAULT_BASELINES_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'ratchet_baselines.json') + +import numpy as np +import pandas as pd + +# Confirmed against a real loso_results.csv (2026-07-08, CatBoost run): +# columns are fold, holdout_site, n_train, n_pos_train, n_neg_train, n_test, +# n_pos_test, n_neg_test, age_cond_auroc, n_age_pairs, std_auroc, accuracy, +# sensitivity, specificity, ppv, TP, FP, FN, TN, top5_features. Candidates +# below are ordered with the confirmed real names first; older/alternate +# names kept as fallbacks in case a different loso_cv.py revision is used. +_AUROC_COL_CANDIDATES = ['age_cond_auroc', 'age_conditioned_auroc', 'auroc_age', 'auroc'] +_N_POS_COL_CANDIDATES = ['n_pos_test', 'n_pos', 'num_pos', 'n_positive'] +_N_NEG_COL_CANDIDATES = ['n_neg_test', 'n_neg', 'num_neg', 'n_negative'] +_SITE_COL_CANDIDATES = ['holdout_site', 'fold', 'site'] + + +def _resolve_column(df, candidates, explicit, quantity_name): + if explicit is not None: + if explicit not in df.columns: + raise ValueError( + f"--{quantity_name}-col '{explicit}' not found in loso_results.csv. " + f"Actual columns: {list(df.columns)}") + return explicit + for cand in candidates: + if cand in df.columns: + return cand + raise ValueError( + f"Could not find a column for '{quantity_name}' in loso_results.csv. " + f"Tried: {candidates}. Actual columns: {list(df.columns)}. " + f"Pass --{quantity_name}-col explicitly to point at the right one.") + + +def hanley_mcneil_se(auroc, n_pos, n_neg): + """ + Standard error of an AUROC estimate (Hanley & McNeil, 1982). + Same formula already used by hand in learning_log.md's 2026-07-01 + entry to confirm the Entry 1-3 leaderboard decline was noise. + """ + if n_pos <= 1 or n_neg <= 1: + raise ValueError( + f"Need n_pos > 1 and n_neg > 1 for a Hanley-McNeil SE " + f"(got n_pos={n_pos}, n_neg={n_neg}).") + + q1 = auroc / (2 - auroc) + q2 = (2 * auroc ** 2) / (1 + auroc) + + variance = ( + auroc * (1 - auroc) + + (n_pos - 1) * (q1 - auroc ** 2) + + (n_neg - 1) * (q2 - auroc ** 2) + ) / (n_pos * n_neg) + + return float(np.sqrt(max(variance, 0.0))) + + +def se_of_fold_mean(fold_aurocs_n_pos_n_neg): + """ + SE of a SIMPLE (unweighted) mean of per-fold AUROCs — matching + loso_cv.py's own "Mean age-conditioned AUROC" convention (confirmed + 2026-07-08 against a real loso_summary.txt: it's a plain average + across folds, e.g. (0.5493+0.5393)/2 = 0.5443 for the 2-fold I0006/S0001 + LOSO setup — NOT weighted by each fold's test-set size). + + For n independent fold estimates, Var(mean) = (1/n^2) * sum(Var_i), so + SE(mean) = sqrt(sum(SE_i^2)) / n. This is standard for an average of + independent estimates and differs from a single Hanley-McNeil SE + computed on pooled/summed n_pos+n_neg, which would implicitly (and + wrongly here) assume one AUROC was computed over the full pooled + sample rather than averaged from per-site folds. + """ + n = len(fold_aurocs_n_pos_n_neg) + if n == 0: + raise ValueError("Need at least one fold to compute a mean SE.") + variance_sum = sum( + hanley_mcneil_se(auroc, n_pos, n_neg) ** 2 + for auroc, n_pos, n_neg in fold_aurocs_n_pos_n_neg + ) + return float(np.sqrt(variance_sum) / n) + + +def compare_auroc(candidate_folds, baseline_folds, sigma_threshold=1.0): + """ + Compares a candidate's mean age-conditioned AUROC against a baseline's, + where "mean" means the same simple per-fold average loso_cv.py itself + reports (see se_of_fold_mean). Each *_folds argument is a list of + (auroc, n_pos, n_neg) tuples, one per LOSO holdout fold. + + One-directional ratchet: only fails on a regression beyond + sigma_threshold pooled SEs. An improving candidate always passes. + """ + candidate_auroc = float(np.mean([f[0] for f in candidate_folds])) + baseline_auroc = float(np.mean([f[0] for f in baseline_folds])) + + se_candidate = se_of_fold_mean(candidate_folds) + se_baseline = se_of_fold_mean(baseline_folds) + se_diff = float(np.sqrt(se_candidate ** 2 + se_baseline ** 2)) + + delta = candidate_auroc - baseline_auroc + sigmas = delta / se_diff if se_diff > 0 else float('inf') + + if delta >= 0: + verdict = 'PASS (improvement)' + is_regression = False + elif sigmas >= -sigma_threshold: + verdict = 'PASS (within noise band)' + is_regression = False + else: + verdict = 'FAIL (regression exceeds noise band)' + is_regression = True + + return { + 'metric': 'age_cond_auroc', + 'candidate': candidate_auroc, + 'baseline': baseline_auroc, + 'delta': delta, + 'se_diff': se_diff, + 'sigmas': sigmas, + 'sigma_threshold': sigma_threshold, + 'verdict': verdict, + 'is_regression': is_regression, + } + + +def compare_reward(candidate_reward, baseline, max_pct_drop=0.15): + """ + Reward has no clean analytic SE — it's a prevalence-weighted score per + patient (compute_reward in evaluate_model.py), not a simple binomial + rate, so Hanley-McNeil doesn't apply. This is a blunt percentage-drop + heuristic, explicitly informational rather than a rigorous statistical + gate. A proper version would bootstrap over loso_probabilities.csv + (resample patients with replacement, recompute reward each time) — not + implemented here; flag if you want that added. + """ + baseline_reward = baseline['reward'] + if baseline_reward == 0: + pct_drop = float('inf') if candidate_reward < 0 else 0.0 + else: + pct_drop = (baseline_reward - candidate_reward) / abs(baseline_reward) + + if candidate_reward >= baseline_reward: + verdict = 'PASS (improvement)' + is_regression = False + elif pct_drop <= max_pct_drop: + verdict = 'PASS (within tolerance, NOT statistically rigorous)' + is_regression = False + else: + verdict = 'WARN (drop exceeds tolerance, NOT statistically rigorous)' + is_regression = True + + return { + 'metric': 'reward', + 'candidate': candidate_reward, + 'baseline': baseline_reward, + 'pct_drop': pct_drop, + 'max_pct_drop': max_pct_drop, + 'verdict': verdict, + 'is_regression': is_regression, + 'note': 'Percentage-drop heuristic only — reward has no analytic SE. ' + 'Not as rigorous as the AUROC check.', + } + + +def load_baselines(path=_DEFAULT_BASELINES_PATH): + path = Path(path) + if not path.exists(): + raise FileNotFoundError( + f"Baseline file not found: {path}. Run from the repo root, or " + f"pass --baselines-file to point at it explicitly.") + with open(path) as f: + data = json.load(f) + data.pop('_readme', None) + return data + + +def summarize_loso_results(csv_path, auroc_col=None, n_pos_col=None, n_neg_col=None, + site_col=None): + """ + Parses per-fold LOSO results into the (auroc, n_pos, n_neg) tuple list + compare_auroc() needs. Does NOT collapse to a single pooled number here + — the simple-mean-of-folds convention (matching loso_cv.py's own + "Mean age-conditioned AUROC" line) is applied by the caller via + compare_auroc(), so this function stays a pure, honest parse of the CSV. + """ + df = pd.read_csv(csv_path) + + auroc_col = _resolve_column(df, _AUROC_COL_CANDIDATES, auroc_col, 'auroc') + n_pos_col = _resolve_column(df, _N_POS_COL_CANDIDATES, n_pos_col, 'n-pos') + n_neg_col = _resolve_column(df, _N_NEG_COL_CANDIDATES, n_neg_col, 'n-neg') + try: + site_col = _resolve_column(df, _SITE_COL_CANDIDATES, site_col, 'site') + sites = df[site_col].astype(str).tolist() + except ValueError: + sites = [f'fold_{i}' for i in range(len(df))] # site label is cosmetic only + + folds = list(zip(df[auroc_col].astype(float), df[n_pos_col].astype(int), + df[n_neg_col].astype(int))) + return folds, sites + + +def _folds_from_baseline(baseline): + """Builds the (auroc, n_pos, n_neg) tuple list from a baseline dict's + 'folds' list. Raises a clear error if the baseline predates this schema + (i.e. still uses the old pooled n_pos/n_neg format).""" + if 'folds' not in baseline: + raise ValueError( + "This baseline uses the old pooled-total schema (n_pos/n_neg at " + "the top level) and can't be used with the corrected per-fold SE " + "calculation. Update ratchet_baselines.json to include a 'folds' " + "list: [{'site':..., 'age_cond_auroc':..., 'n_pos':..., 'n_neg':...}, ...]") + return [(f['age_cond_auroc'], f['n_pos'], f['n_neg']) for f in baseline['folds']] + + +def check_ratchet(loso_results_csv, baseline_key, baselines_file=_DEFAULT_BASELINES_PATH, + candidate_reward=None, sigma_threshold=1.0, max_reward_pct_drop=0.15, + auroc_col=None, n_pos_col=None, n_neg_col=None, site_col=None): + """ + Main entry point — also called from check_submission_files.py. + Returns (passed: bool, results: list[dict], messages: list[str]). + """ + messages = [] + baselines = load_baselines(baselines_file) + + if baseline_key not in baselines: + raise KeyError( + f"Baseline key '{baseline_key}' not found in {baselines_file}. " + f"Available: {list(baselines.keys())}") + + baseline = baselines[baseline_key] + baseline_folds = _folds_from_baseline(baseline) + + if not baseline.get('leaderboard_confirmed', False): + messages.append( + f"CAUTION: baseline '{baseline_key}' is LOSO-only — it has never been " + f"submitted to the real leaderboard (leaderboard_confirmed=false). " + f"A PASS here means 'not worse than your best LOSO estimate so far', " + f"NOT 'confirmed to perform on the real leaderboard'. Your only two " + f"known LOSO->leaderboard deltas (both at small scale) had DIFFERENT " + f"signs of surprise (+0.125, then +0.076) — do not assume this offset " + f"transfers to large scale, where you have zero leaderboard data points " + f"so far." + ) + + candidate_folds, candidate_sites = summarize_loso_results( + loso_results_csv, auroc_col, n_pos_col, n_neg_col, site_col) + + if len(candidate_folds) != len(baseline_folds): + messages.append( + f"candidate has {len(candidate_folds)} fold(s), baseline has " + f"{len(baseline_folds)} — simple-mean comparison still works but isn't " + f"apples-to-apples if these represent different LOSO rotations " + f"(e.g. 2-site vs 3-site with I0002 included).") + + auroc_result = compare_auroc(candidate_folds, baseline_folds, sigma_threshold=sigma_threshold) + auroc_result['candidate_sites'] = candidate_sites + results = [auroc_result] + + if candidate_reward is not None: + results.append(compare_reward( + candidate_reward, baseline, max_pct_drop=max_reward_pct_drop)) + else: + messages.append( + "No --candidate-reward passed — reward ratchet check skipped. " + "AUROC-only gate is NOT a full picture; pass reward explicitly " + "once you have a pooled threshold sweep result.") + + passed = not any(r['is_regression'] for r in results) + return passed, results, messages + + +def _print_report(baseline_key, results, messages): + print(f"Ratchet check against baseline: {baseline_key}\n") + + cautions = [m for m in messages if m.startswith('CAUTION')] + other_messages = [m for m in messages if not m.startswith('CAUTION')] + for c in cautions: + print(f"!! {c}\n") + + for r in results: + print(f" [{r['metric']}]") + print(f" candidate: {r['candidate']:.4f} baseline: {r['baseline']:.4f}") + if r['metric'] == 'age_cond_auroc': + sites = r.get('candidate_sites') + if sites: + print(f" candidate folds: {sites}") + print(f" delta: {r['delta']:+.4f} SE_diff: {r['se_diff']:.4f} " + f"sigmas: {r['sigmas']:+.2f} (threshold: {r['sigma_threshold']})") + else: + print(f" pct_drop: {r['pct_drop']:.1%} tolerance: {r['max_pct_drop']:.1%}") + print(f" note: {r['note']}") + print(f" verdict: {r['verdict']}\n") + + for m in other_messages: + print(f"NOTE: {m}\n") + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--loso-results', required=True, + help='Path to loso_cv.py\'s per-fold results CSV.') + parser.add_argument('--baseline', required=True, + help='Key into ratchet_baselines.json (e.g. small_entry3).') + parser.add_argument('--baselines-file', default=_DEFAULT_BASELINES_PATH) + parser.add_argument('--candidate-reward', type=float, default=None, + help='Pooled reward at your chosen threshold, if you have one ' + '(from a separate threshold-sweep step). Optional.') + parser.add_argument('--sigma-threshold', type=float, default=1.0, + help='AUROC regression must exceed this many pooled SEs to fail. ' + 'Default 1.0 (matches the "noise floor" framing already used ' + 'in learning_log.md).') + parser.add_argument('--max-reward-pct-drop', type=float, default=0.15, + help='Reward regression tolerance as a fraction of baseline. ' + 'Default 0.15 (15%%). NOT statistically derived.') + parser.add_argument('--auroc-col', default=None) + parser.add_argument('--n-pos-col', default=None) + parser.add_argument('--n-neg-col', default=None) + parser.add_argument('--site-col', default=None, + help='Optional — used only for display, not the SE calculation.') + args = parser.parse_args() + + try: + passed, results, messages = check_ratchet( + args.loso_results, args.baseline, args.baselines_file, + candidate_reward=args.candidate_reward, + sigma_threshold=args.sigma_threshold, + max_reward_pct_drop=args.max_reward_pct_drop, + auroc_col=args.auroc_col, n_pos_col=args.n_pos_col, n_neg_col=args.n_neg_col, + site_col=args.site_col, + ) + except (FileNotFoundError, KeyError, ValueError) as e: + print(f"ERROR: {e}", file=sys.stderr) + sys.exit(2) + + _print_report(args.baseline, results, messages) + print('PASS' if passed else 'FAIL') + sys.exit(0 if passed else 1) \ No newline at end of file diff --git a/tools/reg_sweep.py b/tools/reg_sweep.py new file mode 100644 index 0000000..8c8eaa7 --- /dev/null +++ b/tools/reg_sweep.py @@ -0,0 +1,394 @@ +#!/usr/bin/env python + +""" +reg_sweep.py +Team Narnia — PhysioNet Challenge 2026 + +Ridge / Lasso / ElasticNet regularization sweep for the logreg model family, +plus a per-fold coefficient-stability check (does the same feature keep a +similar sign/magnitude across LOSO folds, or is the additive signal +fold-specific?). + +Built standalone rather than bolted onto loso_cv.py, because that file's +current CLI surface (hyperparameter flags, single-fold-only mode) wasn't +available to build against. This script: + - Uses evaluate_model.py's ACTUAL compute_auroc_age / compute_reward / + compute_prevalence functions directly (imported, not reimplemented) — + metric fidelity to the official scorer is guaranteed by construction. + - Writes loso_results.csv / loso_threshold_sweep.csv in the EXACT column + schema confirmed against a real loso_cv.py output (2026-07-08, CatBoost + run) — drops straight into ratchet_check.py with zero changes. + - Implements its own LOSO fold splitting (2-site or 3-site via + --rotate-i0002) rather than assuming loso_cv.py exposes this the same way. + +IMPORTANT — things this script assumes that you should confirm before +trusting cross-comparisons against existing baselines: + +1. REWARD PREVALENCE REFERENCE: CONFIRMED 2026-07-09 against two real runs + (exact reproduction on every threshold tested) — the prevalence + reference is the FULL labeled population across ALL sites, including + whichever site wasn't held out in a given LOSO run. This script now + passes the complete --features-cache population (age, y) as that + reference automatically, since your cache should already include every + site regardless of rotation scheme (I0002 is always at least in + training). No separate prevalence file needed as long as your cache is + built from the full dataset, not a pre-filtered subset. + +2. SCALING: build_logreg_pipeline() (features/pipeline.py) adds a + StandardScaler that was not confirmed to exist in whatever produced the + 0.6002 large-set baseline. See that file's docstring for why this + matters for a C sweep specifically. This is the one remaining + unverified assumption — resolving it needs the actual pipeline code + that generated that baseline, not just its numeric output. + +3. BINARY THRESHOLD for the per-fold sensitivity/specificity/ppv/accuracy + columns is fixed at 0.5 on the CALIBRATED probability — matching + evaluate_model.py's own binary_predictions convention. The + reward-optimal threshold is chosen separately via the pooled threshold + sweep (loso_threshold_sweep.csv), same two-stage structure as your + existing Entry 3 THRESHOLD convention. + +Usage: + python reg_sweep.py --features-cache small_features.npz \\ + --penalties l2,l1,elasticnet --C 0.01,0.1,1,10 --l1-ratios 0.3,0.5,0.7 \\ + --output-dir outputs/reg_sweep_small + + python reg_sweep.py --features-cache large_features.npz --rotate-i0002 \\ + --penalties l2 --C 0.1,1 --output-dir outputs/reg_sweep_large_top2 + +Expected --features-cache format: a .npz file with arrays: + X — (n_patients, 48) float, PRE-AgeResidualizer feature vector + (same 48-length hstack team_code.py/loso_cv.py produce) + y — (n_patients,) int/bool, Cognitive_Impairment label + site — (n_patients,) str, SiteID (e.g. 'S0001', 'I0006', 'I0002') + age — (n_patients,) float, Age + +This is NOT the same as any file already in your repo — you'll need to +build/cache this once from your existing extraction step. If you already +have loso_probabilities.csv-style per-patient data with features attached, +adapt build_dataset_npz() below instead of re-extracting from EDFs. +""" + +import argparse +import ast +import itertools +import os +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +from sklearn.metrics import roc_auc_score + +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +from features import FEATURE_NAMES_48, FEATURE_NAMES_50 +from features.pipeline import build_logreg_pipeline, extract_fitted_coefficients +from evaluate_model import compute_prevalence, compute_reward, compute_auroc as _compute_auroc + +# ── Fold definitions ───────────────────────────────────────────────────────── +ALWAYS_TRAIN_SITE = 'I0002' # only relevant when --rotate-i0002 is NOT passed +TWO_SITE_ROTATION = ['I0006', 'S0001'] +THREE_SITE_ROTATION = ['I0006', 'S0001', 'I0002'] + +REQUIRED_LOSO_RESULTS_COLUMNS = [ + 'fold', 'holdout_site', 'n_train', 'n_pos_train', 'n_neg_train', + 'n_test', 'n_pos_test', 'n_neg_test', 'age_cond_auroc', 'n_age_pairs', + 'std_auroc', 'accuracy', 'sensitivity', 'specificity', 'ppv', + 'TP', 'FP', 'FN', 'TN', 'top5_features', +] + + +def compute_auroc_age_with_pairs(labels, predictions, ages, gap=0): + """ + Mirrors evaluate_model.py's compute_auroc_age EXACTLY (same loop, same + tie-handling) but also returns the pair count (denom), matching + loso_results.csv's n_age_pairs column. Do not let this drift from + evaluate_model.py's own logic — it's a read-only duplicate purely to + expose one more number for logging, never edit evaluate_model.py itself. + """ + labels = np.asarray(labels) + predictions = np.asarray(predictions) + ages = np.asarray(ages, dtype=float) + + idx_pos = np.where(labels == 1)[0] + idx_neg = np.where(labels == 0)[0] + + numer = 0.0 + denom = 0 + for i in idx_pos: + for j in idx_neg: + if abs(ages[i] - ages[j]) <= gap: + if predictions[i] > predictions[j]: + numer += 1 + elif predictions[i] == predictions[j]: + numer += 0.5 + denom += 1 + auroc = numer / denom if denom > 0 else float('nan') + return auroc, denom + + +def load_feature_cache(path): + path = Path(path) + if not path.exists(): + raise FileNotFoundError(f"Feature cache not found: {path}") + data = np.load(path, allow_pickle=True) + for key in ('X', 'y', 'site', 'age'): + if key not in data: + raise ValueError( + f"'{key}' missing from {path}. Found keys: {list(data.keys())}. " + f"See this script's module docstring for the expected format.") + X = np.asarray(data['X'], dtype=float) + y = np.asarray(data['y']).astype(int) + site = np.asarray(data['site']).astype(str) + age = np.asarray(data['age'], dtype=float) + if X.shape[1] != len(FEATURE_NAMES_48): + raise ValueError( + f"X has {X.shape[1]} columns, expected {len(FEATURE_NAMES_48)} " + f"(the pre-AgeResidualizer 48-length vector). Check your cache " + f"was built from the same extraction order as team_code.py.") + n = len(y) + if not (len(site) == len(age) == n == X.shape[0]): + raise ValueError( + f"Mismatched lengths: X={X.shape[0]}, y={n}, site={len(site)}, age={len(age)}") + return X, y, site, age + + +def make_folds(site, rotate_i0002): + """Yields (holdout_site, train_idx, test_idx). I0002 is excluded from + the rotation (always in training) unless --rotate-i0002 is passed — + same logic/justification as loso_cv.py's own KNOWN_MISLABELED_PATIENT_IDS- + adjacent fold setup (learning_log.md, 2026-07-05: I0002 no longer + justified as a permanent-training-only site at large scale).""" + holdout_sites = THREE_SITE_ROTATION if rotate_i0002 else TWO_SITE_ROTATION + for holdout in holdout_sites: + test_idx = np.where(site == holdout)[0] + train_idx = np.where(site != holdout)[0] + if len(test_idx) == 0: + raise ValueError(f"No patients found for holdout site '{holdout}' — check site labels.") + yield holdout, train_idx, test_idx + + +def _binary_metrics(y_true, probs, threshold=0.5): + pred = (probs > threshold).astype(int) + tp = int(np.sum((y_true == 1) & (pred == 1))) + fp = int(np.sum((y_true == 0) & (pred == 1))) + fn = int(np.sum((y_true == 1) & (pred == 0))) + tn = int(np.sum((y_true == 0) & (pred == 0))) + accuracy = (tp + tn) / max(tp + fp + fn + tn, 1) + sensitivity = tp / max(tp + fn, 1) + specificity = tn / max(tn + fp, 1) + ppv = tp / max(tp + fp, 1) + return dict(TP=tp, FP=fp, FN=fn, TN=tn, accuracy=accuracy, + sensitivity=sensitivity, specificity=specificity, ppv=ppv) + + +def run_one_config(X, y, age, site, rotate_i0002, penalty, C, l1_ratio, + calibrated=True, threshold=0.5, config_label=None): + """ + Runs one LOSO cycle for one (penalty, C, l1_ratio) config. + Returns (results_rows, coef_rows, pooled_df) where pooled_df has one row + per patient across ALL test folds (site, age, label, probability) — used + both for the reward threshold sweep and as this config's own prevalence + reference (see module docstring, assumption #1). + """ + results_rows = [] + coef_rows = [] + pooled_records = [] + + for fold_i, (holdout, train_idx, test_idx) in enumerate(make_folds(site, rotate_i0002)): + X_train, y_train = X[train_idx], y[train_idx] + X_test, y_test, age_test = X[test_idx], y[test_idx], age[test_idx] + + pipeline = build_logreg_pipeline( + y_train, penalty=penalty, C=C, l1_ratio=l1_ratio, calibrated=calibrated) + pipeline.fit(X_train, y_train) + + probs = pipeline.predict_proba(X_test)[:, 1] + + age_cond_auroc, n_age_pairs = compute_auroc_age_with_pairs(y_test, probs, age_test, gap=2) + std_auroc = _compute_auroc(y_test, probs) + bm = _binary_metrics(y_test, probs, threshold=threshold) + + coefs = extract_fitted_coefficients(pipeline) + top5_idx = np.argsort(-np.abs(coefs))[:5] + top5 = [(FEATURE_NAMES_50[i], round(float(coefs[i]), 4)) for i in top5_idx] + + results_rows.append({ + 'fold': f'holdout_{holdout}', + 'holdout_site': holdout, + 'n_train': len(train_idx), + 'n_pos_train': int(y_train.sum()), + 'n_neg_train': int((~y_train.astype(bool)).sum()), + 'n_test': len(test_idx), + 'n_pos_test': int(y_test.sum()), + 'n_neg_test': int((~y_test.astype(bool)).sum()), + 'age_cond_auroc': round(age_cond_auroc, 4), + 'n_age_pairs': n_age_pairs, + 'std_auroc': round(std_auroc, 4), + 'accuracy': round(bm['accuracy'], 4), + 'sensitivity': round(bm['sensitivity'], 4), + 'specificity': round(bm['specificity'], 4), + 'ppv': round(bm['ppv'], 4), + 'TP': bm['TP'], 'FP': bm['FP'], 'FN': bm['FN'], 'TN': bm['TN'], + 'top5_features': str(top5), + }) + + for feat_name, coef_val in zip(FEATURE_NAMES_50, coefs): + coef_rows.append({'fold': f'holdout_{holdout}', 'feature': feat_name, + 'coefficient': float(coef_val)}) + + for i, p in zip(test_idx, probs): + pooled_records.append({'site': site[i], 'age': age[i], 'label': y[i], 'probability': p}) + + pooled_df = pd.DataFrame(pooled_records) + return results_rows, coef_rows, pooled_df + + +def threshold_sweep(pooled_df, thresholds, prevalence_ages=None, prevalence_labels=None): + """ + Pooled reward/AUROC threshold sweep, matching your existing + loso_threshold_sweep.csv columns (threshold, reward, auroc). + + CONFIRMED 2026-07-09 against two real runs (2-site logreg + 3-site + I0002-rotation logreg, same underlying data): the prevalence reference + for compute_prevalence is the FULL available labeled population — + ALL sites, INCLUDING whichever site wasn't held out / scored in a given + run. Verified by exact reproduction on every tested threshold: scoring + the 2-site (6281-patient) pool using only itself as the prevalence + reference gave numbers 20-40% too high; scoring it using the full + 3-site (6600-patient) pool as the reference reproduced the real + reported reward EXACTLY at every threshold tested. This makes sense as + a design choice — local age-prevalence should reflect the true + population base rate from every available label, not just whichever + patients happen to be in a particular fold's test set. + + prevalence_ages/prevalence_labels: pass the FULL dataset's age/label + arrays (not just the pooled test predictions) — i.e. every patient in + your --features-cache, train and test folds combined, regardless of + rotation scheme. If omitted, falls back to the pooled test set as its + own reference (NOT recommended — this was the old, disproven default; + kept only so the function still runs standalone without a full dataset + handy, e.g. for quick synthetic smoke tests). + + The 'auroc' column is pooled age-conditioned AUROC over the combined + scored population (gap=2) — confirmed exact match separately, see + compute_auroc_age_with_pairs usage below. + """ + labels = pooled_df['label'].to_numpy() + ages = pooled_df['age'].to_numpy(dtype=float) + probs = pooled_df['probability'].to_numpy() + + if prevalence_ages is not None and prevalence_labels is not None: + ref_ages, ref_labels = prevalence_ages, prevalence_labels + else: + ref_ages, ref_labels = ages, labels + + age_to_prevalence = compute_prevalence(ages, ref_labels, ref_ages, gap=2) + pooled_age_cond_auroc, _ = compute_auroc_age_with_pairs(labels, probs, ages, gap=2) + + rows = [] + for t in thresholds: + preds = (probs > t).astype(int) + reward = compute_reward(labels, preds, ages, age_to_prevalence) + rows.append({'threshold': t, 'reward': round(reward, 4), 'auroc': round(pooled_age_cond_auroc, 4)}) + return pd.DataFrame(rows) + + +def config_dirname(penalty, C, l1_ratio): + if penalty == 'elasticnet': + return f'logreg_elasticnet_C{C}_l1r{l1_ratio}' + return f'logreg_{penalty}_C{C}' + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--features-cache', required=True) + parser.add_argument('--rotate-i0002', action='store_true', + help='Use 3-site LOSO rotation (I0006, S0001, I0002) instead of the default 2-site.') + parser.add_argument('--penalties', default='l2', help='Comma-separated: l1,l2,elasticnet') + parser.add_argument('--C', default='1.0', help='Comma-separated C values, e.g. 0.01,0.1,1,10') + parser.add_argument('--l1-ratios', default='0.5', help='Comma-separated, only used for elasticnet') + parser.add_argument('--threshold', type=float, default=0.5, + help='Binary decision threshold for per-fold sens/spec/ppv/accuracy columns.') + parser.add_argument('--reward-thresholds', default='0.05,0.08,0.10,0.12,0.15,0.20', + help='Thresholds swept for the pooled reward table.') + parser.add_argument('--no-calibration', action='store_true', + help='Skip CalibratedClassifierCV wrapping (faster, for quick screening).') + parser.add_argument('--output-dir', required=True) + parser.add_argument('--max-configs', type=int, default=40, + help='Safety cap on total sweep size — refuses to run a grid larger than this ' + 'without --force, so a typo in --C doesn\'t accidentally launch 200 fits.') + parser.add_argument('--force', action='store_true') + args = parser.parse_args() + + X, y, site, age = load_feature_cache(args.features_cache) + print(f"Loaded {len(y)} patients ({int(y.sum())} positive) from {args.features_cache}") + print(f"Sites: {dict(zip(*np.unique(site, return_counts=True)))}") + + penalties = args.penalties.split(',') + C_values = [float(c) for c in args.C.split(',')] + l1_ratios = [float(r) for r in args.l1_ratios.split(',')] if 'elasticnet' in penalties else [None] + + configs = [] + for penalty in penalties: + if penalty == 'elasticnet': + for C, l1r in itertools.product(C_values, l1_ratios): + configs.append((penalty, C, l1r)) + else: + for C in C_values: + configs.append((penalty, C, None)) + + if len(configs) > args.max_configs and not args.force: + print(f"ERROR: sweep grid has {len(configs)} configs, exceeding --max-configs={args.max_configs}. " + f"Narrow your grid or pass --force if this is intentional.", file=sys.stderr) + sys.exit(2) + + reward_thresholds = [float(t) for t in args.reward_thresholds.split(',')] + + output_dir = Path(args.output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + summary_rows = [] + for penalty, C, l1_ratio in configs: + label = config_dirname(penalty, C, l1_ratio) + print(f"\n=== {label} ===") + config_dir = output_dir / label + config_dir.mkdir(parents=True, exist_ok=True) + + results_rows, coef_rows, pooled_df = run_one_config( + X, y, age, site, args.rotate_i0002, penalty, C, l1_ratio, + calibrated=not args.no_calibration, threshold=args.threshold) + + results_df = pd.DataFrame(results_rows)[REQUIRED_LOSO_RESULTS_COLUMNS] + results_df.to_csv(config_dir / 'loso_results.csv', index=False) + + pd.DataFrame(coef_rows).to_csv(config_dir / 'coefficients.csv', index=False) + + sweep_df = threshold_sweep(pooled_df, reward_thresholds, + prevalence_ages=age, prevalence_labels=y) + sweep_df.to_csv(config_dir / 'loso_threshold_sweep.csv', index=False) + + mean_auroc = results_df['age_cond_auroc'].mean() + best_reward_row = sweep_df.loc[sweep_df['reward'].idxmax()] + print(f" mean age_cond_auroc: {mean_auroc:.4f}") + print(f" best reward: {best_reward_row['reward']:.4f} @ t={best_reward_row['threshold']}") + + summary_rows.append({ + 'config': label, 'penalty': penalty, 'C': C, 'l1_ratio': l1_ratio, + 'mean_age_cond_auroc': round(mean_auroc, 4), + 'best_reward': round(best_reward_row['reward'], 4), + 'best_reward_threshold': best_reward_row['threshold'], + }) + + summary_df = pd.DataFrame(summary_rows).sort_values('mean_age_cond_auroc', ascending=False) + summary_df.to_csv(output_dir / 'sweep_summary.csv', index=False) + print(f"\n{'='*60}\nSweep summary (sorted by mean age_cond_auroc):\n{'='*60}") + print(summary_df.to_string(index=False)) + print(f"\nPer-config loso_results.csv files are ready for ratchet_check.py, e.g.:") + print(f" python check_submission_files.py --loso-results " + f"{output_dir}/{summary_df.iloc[0]['config']}/loso_results.csv " + f"--ratchet-baseline large_logreg") + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/tools/test_combined_features.py b/tools/test_combined_features.py new file mode 100644 index 0000000..dc80828 --- /dev/null +++ b/tools/test_combined_features.py @@ -0,0 +1,456 @@ +#!/usr/bin/env python3 +""" +test_combined_features.py — Team Narnia + +Combined-feature ablation for NF1 (Markov transitions, 9 features) and +NF2 (arousal clustering, 2 features) — tested TOGETHER as an 11-feature +block against the shipped 48-feature baseline, NOT one-at-a-time. + +This deliberately departs from the project's one-variable-per-test +discipline (learning_log.md, 2026-07-01 decision log). Justification, +per TEST_PLAN.md's dependency graph: RB1 (2026-07-13) ruled out the +regularization-budget explanation for the 3-for-3 CAISR null pattern +(mean_wake_bout_duration, rem_latency, no_arousal_entropy — all null at +both C=0.01 and C=1.0). With regularization budget eliminated, the +remaining live hypothesis for why real univariate EDA signal keeps +dying in the full pipeline is that single-feature ablation is +structurally blind to INTERACTION effects — a feature that only helps +conditional on another feature already being present. Testing NF1+NF2 +combined is the direct test of that hypothesis. This is branch (a) of +the RB1-null decision point; branch (b) (ceiling + calibration/threshold +tuning) was not chosen this round. + +REQUIRED collinearity checks (per TEST_PLAN.md's NF1/NF2 internal +dependencies) are computed and printed BEFORE the ablation runs, not +buried afterward — a null or positive ablation result is uninterpretable +without knowing whether the new features are redundant with what's +already shipped: + - NF1's stationary distribution vs. existing Wake_pct/N1_pct/N2_pct/ + N3_pct/REM_pct (caisr_base.py) + - NF2's CV/burst_index vs. existing Arousal_idx/Spontaneous_arousal_idx + (caisr_base.py / caisr_enriched.py) + +Same fixed config as test_spectral_ablation.py: logreg, l2, calibrated +(cv5), AgeResidualizer on. --C defaults to 0.01 (shipped value) since +RB1 already established C=1.0 doesn't help single-feature candidates — +no a priori reason a combined test needs different regularization, but +the flag is included for consistency/completeness, not because it's +expected to matter here. + +Usage: + python tools/test_combined_features.py \\ + --data /path/to/training_set \\ + --out outputs/nf1_nf2_combined \\ + --rotate-i0002 + +Outputs: + collinearity_report.txt — correlations, printed AND written, so + they aren't lost if the run is re-piped. + ablation_results_combined.csv + ablation_summary_combined.txt +""" + +import argparse +import os +import sys +import warnings + +import numpy as np +import pandas as pd +from pathlib import Path +from tqdm import tqdm + +warnings.filterwarnings('ignore') + +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +from helper_code import ( + find_patients, load_demographics, load_signal_data, + load_diagnoses, DEMOGRAPHICS_FILE, + PHYSIOLOGICAL_DATA_SUBFOLDER, + ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + HEADERS, +) +from features.demographic import extract_demographic_features +from features.caisr_base import extract_caisr_base_features +from features.caisr_enriched import ( + extract_caisr_enriched_features, + extract_markov_transition_features, N_MARKOV_FEATURES, + extract_arousal_clustering_features, N_AROUSAL_CLUSTERING_FEATURES, +) +from features.physiological_ratios import (extract_physiological_ratio_features, + N_RATIO_FEATURES) +from features import N_CAISR_BASE_FEATURES, N_CAISR_ENRICHED_FEATURES +from features.pipeline import build_logreg_pipeline +from evaluate_model import compute_reward, compute_prevalence + +_N_BASE = N_CAISR_BASE_FEATURES +_N_ENRICHED = N_CAISR_ENRICHED_FEATURES +_N_RATIO = N_RATIO_FEATURES +N_COMBINED = N_MARKOV_FEATURES + N_AROUSAL_CLUSTERING_FEATURES # 11 + +# Same denylist as every other large-set-comparable script. +KNOWN_MISLABELED_PATIENT_IDS = {'115257116'} + +# Indices into the shipped 48-vector for the collinearity check. Per +# FEATURES.md: demo(10) + caisr_base(12) + caisr_enriched(11) + ratios(15). +# caisr_base occupies [10:22]; within it (FEATURES.md rows 59-70): +# Wake_pct=idx13(local2), N1_pct=14(3), N2_pct=15(4), N3_pct=16(5), +# REM_pct=17(6), Arousal_idx=idx11(local1) — CONFIRM against your +# actual caisr_base.py return order before trusting this block; written +# from FEATURES.md's documented order, not verified against the live +# extraction function's actual array order. +_IDX_AROUSAL_IDX = 10 + 1 # caisr_base[1] = Arousal_idx — unchanged, was already correct +_IDX_WAKE_PCT = 10 + 3 # caisr_base[3] = Wake % (was 10+2 — WRONG, fixed) +_IDX_N1_PCT = 10 + 4 # caisr_base[4] = N1 % (was 10+3 — WRONG, fixed) +_IDX_N2_PCT = 10 + 5 # caisr_base[5] = N2 % (was 10+4 — WRONG, fixed) +_IDX_N3_PCT = 10 + 6 # caisr_base[6] = N3 % (was 10+5 — WRONG, fixed) +_IDX_REM_PCT = 10 + 7 # caisr_base[7] = REM % (was 10+6 — WRONG, fixed) +_IDX_SPONT_AROUSAL = 10 + 12 + 9 # caisr_enriched[9] = Spontaneous_arousal_idx — unchanged + + +def extract_all_features_combined(data_folder, verbose=True): + """ + Extracts the existing 48-feature vector AND the combined NF1+NF2 + 11-feature block in one pass. Same one-EDF-load-per-patient + discipline as test_spectral_ablation.py's extraction loop. + """ + demo_path = os.path.join(data_folder, DEMOGRAPHICS_FILE) + patient_list = find_patients(demo_path) + + features_48_list = [] + combined_list = [] + labels_list = [] + ages_list = [] + sites_list = [] + patient_ids = [] + skipped_denied = 0 + + pbar = tqdm(patient_list, desc='Extracting (48-feature + NF1/NF2 combined)', + unit='patient', disable=not verbose) + + for record in pbar: + try: + pid = record[HEADERS['bids_folder']] + site = record[HEADERS['site_id']] + sess = record[HEADERS['session_id']] + + if verbose: + pbar.set_postfix({'patient': pid}) + + pdata = load_demographics(demo_path, pid, sess) + bdsp_id = str(pdata.get('BDSPPatientID', pid)) + + if bdsp_id in KNOWN_MISLABELED_PATIENT_IDS: + skipped_denied += 1 + continue + + demo_f = extract_demographic_features(pdata) + age = float(demo_f[0]) + + phys_file = os.path.join(data_folder, PHYSIOLOGICAL_DATA_SUBFOLDER, + site, f'{pid}_ses-{sess}.edf') + if not os.path.exists(phys_file): + continue + phys_data, phys_fs = load_signal_data(phys_file) + + algo_file = os.path.join(data_folder, ALGORITHMIC_ANNOTATIONS_SUBFOLDER, + site, f'{pid}_ses-{sess}_caisr_annotations.edf') + if os.path.exists(algo_file): + algo_data, _ = load_signal_data(algo_file) + base_f = extract_caisr_base_features(algo_data) + enriched_f = extract_caisr_enriched_features(algo_data) + ratio_f = extract_physiological_ratio_features( + phys_data, phys_fs, algo_data) + markov_f = extract_markov_transition_features(algo_data) + arousal_f = extract_arousal_clustering_features(algo_data) + del algo_data + else: + base_f = np.full(_N_BASE, float('nan')) + enriched_f = np.full(_N_ENRICHED, float('nan')) + ratio_f = np.full(_N_RATIO, float('nan')) + markov_f = np.full(N_MARKOV_FEATURES, float('nan')) + arousal_f = np.full(N_AROUSAL_CLUSTERING_FEATURES, float('nan')) + + del phys_data + + label = load_diagnoses(demo_path, pid) + if label not in (0, 1): + continue + + features_48_list.append( + np.hstack([demo_f, base_f, enriched_f, ratio_f])) + combined_list.append(np.hstack([markov_f, arousal_f])) + labels_list.append(label) + ages_list.append(age) + sites_list.append(site) + patient_ids.append(bdsp_id) + + del base_f, enriched_f, ratio_f, markov_f, arousal_f + + except Exception as ex: + tqdm.write(f' ! Error on {pid}: {ex}') + continue + + pbar.close() + print(f'\nSkipped {skipped_denied} denylisted patient(s).') + + return ( + np.asarray(features_48_list, dtype=np.float32), + np.asarray(combined_list, dtype=np.float32), + np.asarray(labels_list, dtype=int), + np.asarray(ages_list, dtype=float), + np.asarray(sites_list), + np.asarray(patient_ids), + ) + + +def collinearity_report(X48, combined): + """ + Required-before-ablation check (TEST_PLAN.md, NF1/NF2 internal + dependencies). Pearson correlation, computed on complete-case rows + only (both sides non-NaN) per pair — sample size varies per pair + and is reported alongside r, since NF1/NF2's NaN rates differ from + the existing features'. + """ + lines = ['=' * 60, ' Collinearity check (required before ablation)', '=' * 60, ''] + + stationary_labels = ['Wake', 'N1', 'N2', 'N3', 'REM'] + stationary_cols = combined[:, 4:9] # NF1 output indices [4:9] + existing_pct_idx = [_IDX_WAKE_PCT, _IDX_N1_PCT, _IDX_N2_PCT, + _IDX_N3_PCT, _IDX_REM_PCT] + + lines.append('NF1 stationary distribution vs. existing stage %:') + for i, (label, exist_idx) in enumerate(zip(stationary_labels, existing_pct_idx)): + a = stationary_cols[:, i] + b = X48[:, exist_idx] + valid = ~np.isnan(a) & ~np.isnan(b) + n = int(valid.sum()) + if n > 2: + r = float(np.corrcoef(a[valid], b[valid])[0, 1]) + lines.append(f' stationary_{label:<5} vs {label}_pct: r={r:+.3f} (n={n})') + else: + lines.append(f' stationary_{label:<5} vs {label}_pct: n={n}, too few to compute') + + lines.append('') + lines.append('NF2 CV/burst_index vs. existing arousal-rate features:') + cv_col = combined[:, 9] + burst_col = combined[:, 10] + for name, col in [('cv_inter_arousal', cv_col), ('burst_index', burst_col)]: + for exist_name, exist_idx in [('Arousal_idx', _IDX_AROUSAL_IDX), + ('Spontaneous_arousal_idx', _IDX_SPONT_AROUSAL)]: + b = X48[:, exist_idx] + valid = ~np.isnan(col) & ~np.isnan(b) + n = int(valid.sum()) + if n > 2: + r = float(np.corrcoef(col[valid], b[valid])[0, 1]) + lines.append(f' {name} vs {exist_name}: r={r:+.3f} (n={n})') + else: + lines.append(f' {name} vs {exist_name}: n={n}, too few to compute') + + lines.append('') + lines.append('High |r| (>0.9) means that component is likely redundant with an') + lines.append('existing feature — a null OR positive ablation result for that') + lines.append('specific component should be read in that light, per TEST_PLAN.md.') + lines.append('Index assumptions (_IDX_* constants) are taken from FEATURES.md\'s') + lines.append('documented order, NOT independently re-verified against the live') + lines.append('caisr_base.py/caisr_enriched.py return arrays — confirm before') + lines.append('trusting these numbers if anything looks implausible (e.g. |r|') + lines.append('near 0 for stationary_Wake vs Wake_pct would be a red flag).') + + return '\n'.join(lines) + + +def age_conditioned_auroc(labels, probs, ages, gap=2): + idx_pos = np.where(labels == 1)[0] + idx_neg = np.where(labels == 0)[0] + numer = 0.0 + denom = 0 + for i in idx_pos: + for j in idx_neg: + if abs(ages[i] - ages[j]) <= gap: + if probs[i] > probs[j]: + numer += 1.0 + elif probs[i] == probs[j]: + numer += 0.5 + denom += 1 + if denom == 0: + return float('nan'), 0 + return round(numer / denom, 4), denom + + +def hanley_mcneil_se(auc, n_pos, n_neg): + Q1 = auc / (2 - auc) + Q2 = 2 * auc**2 / (1 + auc) + var = (auc*(1-auc) + (n_pos-1)*(Q1-auc**2) + (n_neg-1)*(Q2-auc**2)) / (n_pos*n_neg) + return np.sqrt(var) + + +def run_fold(holdout_site, X, y, ages, sites, age_to_prevalence, label, C=0.01): + train_mask = sites != holdout_site + test_mask = sites == holdout_site + + y_train = y[train_mask] + y_test = y[test_mask] + ages_test = ages[test_mask] + + if int(y_test.sum()) < 3: + print(f' SKIP {label} — fewer than 3 positives in {holdout_site}.') + return None + + model = build_logreg_pipeline(y_train, calibrated=True, + use_age_residual=True, + C=C, penalty=None) + model.fit(X[train_mask], y_train) + probs = model.predict_proba(X[test_mask])[:, 1] + + auroc, n_pairs = age_conditioned_auroc(y_test, probs, ages_test) + preds = (probs > 0.10).astype(int) # shipped Entry 5 threshold + reward = compute_reward(y_test, preds, ages_test, age_to_prevalence) + + n_pos = int(y_test.sum()) + n_neg = int((y_test == 0).sum()) + + return { + 'config': label, + 'holdout_site': holdout_site, + 'C': C, + 'n_pos_test': n_pos, + 'n_neg_test': n_neg, + 'age_cond_auroc': auroc, + 'n_age_pairs': n_pairs, + 'reward_at_t0.10': round(reward, 4), + } + + +def main(): + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + p.add_argument('--data', required=True) + p.add_argument('--out', required=True) + p.add_argument('--C', type=float, default=0.01, + help='Logistic regression C. Default 0.01 matches shipped ' + 'config; RB1 (2026-07-13) already showed C=1.0 does not ' + 'help single-feature candidates, so this is included for ' + 'completeness, not because it is expected to matter here.') + p.add_argument('--rotate-i0002', action='store_true') + p.add_argument('--trimmed', action='store_true', + help='Drop the 5 stationary-distribution features (r=0.998-1.000 ' + 'vs existing stage %% features — near-pure duplicates) and ' + 'burst_index (r=0.73 vs Arousal_idx — meaningfully redundant), ' + 'per the 2026-07-14 collinearity findings. Keeps only the 5 ' + 'components with low redundancy: transition_entropy, ' + 'n3_to_wake_rate, rem_to_n3_rate, escalating_rate, ' + 'cv_inter_arousal. Full collinearity report is still computed ' + 'and printed on all 11 components regardless of this flag, ' + 'for reference.') + p.add_argument('--verbose', action='store_true', default=True) + args = p.parse_args() + + out_dir = Path(args.out) + out_dir.mkdir(parents=True, exist_ok=True) + + print('Extracting 48-feature vector + NF1 (Markov) + NF2 (arousal clustering) ...') + X48, combined, y, ages, sites, patient_ids = extract_all_features_combined( + args.data, verbose=args.verbose) + + print(f'\nExtracted {len(y)} patients ({int(y.sum())} positive)') + print(f'X48 shape: {X48.shape}, combined (NF1+NF2) shape: {combined.shape}') + + nan_rates = np.isnan(combined).mean(axis=0) + labels = ['transition_entropy', 'n3_to_wake_rate', 'rem_to_n3_rate', + 'escalating_rate', 'stationary_Wake', 'stationary_N1', + 'stationary_N2', 'stationary_N3', 'stationary_REM', + 'cv_inter_arousal', 'burst_index'] + print('\nPer-feature NaN rate:') + for lbl, rate in zip(labels, nan_rates): + print(f' {lbl:<22}: {rate:.1%}') + + report = collinearity_report(X48, combined) + print('\n' + report) + out_dir.mkdir(parents=True, exist_ok=True) # defensive re-create — a prior + # run hit FileNotFoundError here despite the mkdir() call at the top of + # main(), implying something removed the directory mid-run. Re-asserting + # it immediately before every write is cheap insurance against that. + (out_dir / 'collinearity_report.txt').write_text(report) + + # Trimmed subset, per 2026-07-14 collinearity findings: stationary + # distribution (indices 4-8) is near-pure duplicate of existing + # stage %% features (r=0.998-1.000); burst_index (index 10) is + # meaningfully redundant with Arousal_idx (r=0.73). Kept indices: + # [0]=transition_entropy, [1]=n3_to_wake_rate, [2]=rem_to_n3_rate, + # [3]=escalating_rate, [9]=cv_inter_arousal. + _TRIMMED_IDX = [0, 1, 2, 3, 9] + _TRIMMED_LABELS = ['transition_entropy', 'n3_to_wake_rate', 'rem_to_n3_rate', + 'escalating_rate', 'cv_inter_arousal'] + + if args.trimmed: + combined_for_model = combined[:, _TRIMMED_IDX] + block_label = f'{len(_TRIMMED_IDX)}-feature trimmed NF1+NF2' + print(f'\n--trimmed: using only {_TRIMMED_LABELS} ' + f'({combined_for_model.shape[1]} of 11 original components)') + else: + combined_for_model = combined + block_label = '11-feature NF1+NF2 combined' + + X_with = np.hstack([X48, combined_for_model]) + n_total_features = X_with.shape[1] + + age_to_prevalence = compute_prevalence(ages, y, ages, gap=2) + + holdout_sites = ['I0006', 'S0001'] + if args.rotate_i0002: + holdout_sites.append('I0002') + + rows = [] + for holdout in holdout_sites: + print(f'\n--- Fold: hold out {holdout} ---') + r_without = run_fold(holdout, X48, y, ages, sites, age_to_prevalence, + label='WITHOUT NF1+NF2 (shipped 48-feature)', + C=args.C) + r_with = run_fold(holdout, X_with, y, ages, sites, age_to_prevalence, + label=f'WITH {block_label} ({n_total_features}-feature)', + C=args.C) + if r_without: rows.append(r_without) + if r_with: rows.append(r_with) + + if r_without and r_with: + se_a = hanley_mcneil_se(r_without['age_cond_auroc'], + r_without['n_pos_test'], r_without['n_neg_test']) + se_b = hanley_mcneil_se(r_with['age_cond_auroc'], + r_with['n_pos_test'], r_with['n_neg_test']) + se_diff = np.sqrt(se_a**2 + se_b**2) + z = (r_with['age_cond_auroc'] - r_without['age_cond_auroc']) / se_diff + print(f" AUROC: without={r_without['age_cond_auroc']:.4f} " + f"with={r_with['age_cond_auroc']:.4f} " + f"delta={r_with['age_cond_auroc']-r_without['age_cond_auroc']:+.4f} " + f"z={z:+.2f}sigma") + print(f" Reward@0.10: without={r_without['reward_at_t0.10']:.4f} " + f"with={r_with['reward_at_t0.10']:.4f}") + + results_df = pd.DataFrame(rows) + suffix = 'trimmed' if args.trimmed else 'combined' + out_dir.mkdir(parents=True, exist_ok=True) # defensive re-create, see note above + out_csv = out_dir / f'ablation_results_{suffix}.csv' + results_df.to_csv(out_csv, index=False) + print(f'\nWrote {out_csv}') + + summary_lines = [ + '=' * 60, + f' NF1 (Markov transitions) + NF2 (arousal clustering) — {block_label} Ablation', + '=' * 60, + '', + 'Fixed config: logreg, C={}, l2, calibrated (cv5), AgeResidualizer on'.format(args.C), + f'Only variable: presence/absence of the {block_label} block', + '(NOT one-at-a-time — this is branch (a) after RB1\'s null, testing', + 'for interaction effects invisible to single-feature ablation)', + '', + results_df.to_string(index=False), + ] + summary = '\n'.join(summary_lines) + (out_dir / f'ablation_summary_{suffix}.txt').write_text(summary) + print('\n' + summary) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/verify_label_integrity.py b/verify_label_integrity.py new file mode 100644 index 0000000..18d20bb --- /dev/null +++ b/verify_label_integrity.py @@ -0,0 +1,192 @@ +""" +verify_label_integrity.py +Team Narnia — PhysioNet Challenge 2026 + +Checks whether a training set's demographics.csv is pre-curated to contain +ONLY clean Positive/Negative patients, per the official 3-way definition +on the Challenge scoring page, or whether it silently mixes "Other" +(ambiguous/excluded) patients into the Cognitive_Impairment=False bucket. + +Background: the small training set was verified clean (every negative has +>=6yr Time_to_Last_Visit, no single-code or wrong-timing patients hiding +in the negatives). There is no guarantee the large training set received +the same curation before release — this script checks that assumption +explicitly rather than trusting it. + +Official definition (physionetchallenges.org/2026/#scoring): + Positive: >=2 CI-related ICD codes, first one 1-6 years after the PSG, + >=7 days between first and last code. + Negative: zero CI-related ICD codes ever, AND >=6 years of encounters + (Time_to_Last_Visit) after the PSG. + Other: everything else (first code <1yr or >6yr out, exactly one + code, or zero codes with <6yr follow-up). EXCLUDED from + official scoring. + +Usage: + python verify_label_integrity.py \ + --demographics /path/to/training_set_large/demographics.csv \ + --icd /path/to/training_set_large/ICD_codes_CI.csv \ + --out large_set_label_audit.csv + +Exit behavior: + Prints a summary to the console either way. If any "Other" patients + are found hiding in the current True/False label column, writes their + IDs to --out so you can decide whether to exclude them from training + (recommended, to match what real scoring will never include as a + negative) before committing to a full training run on the large set. +""" + +import argparse +import csv +from collections import defaultdict +from datetime import datetime + +LOWER_BOUND_DAYS = 365.25 * 1 # 1 year +UPPER_BOUND_DAYS = 365.25 * 6 # 6 years +MIN_SPAN_DAYS = 7 # first-to-last code confirmation +MIN_FOLLOWUP_DAYS = 365.25 * 6 # negative follow-up requirement + + +def parse_date(s): + return datetime.strptime(s[:10], "%Y-%m-%d") + + +def load_icd_dates(icd_path): + """Returns dict: BDSPPatientID -> sorted list of ICD code dates (str).""" + by_patient = defaultdict(list) + with open(icd_path) as f: + for row in csv.DictReader(f): + by_patient[row["BDSPPatientID"]].append(row["ICDDate"]) + return by_patient + + +def classify_official(psg_date, icd_dates, time_to_last_visit_days): + """ + Reconstructs the official Positive/Negative/Other category from raw + ICD dates and follow-up duration, independent of whatever label the + local create_labels.py may have already assigned. + + Returns (category, reason) where category is one of + 'Positive', 'Negative', 'Other', and reason explains Other cases. + """ + n_codes = len(icd_dates) + + if n_codes == 0: + if time_to_last_visit_days is not None and time_to_last_visit_days >= MIN_FOLLOWUP_DAYS: + return "Negative", None + return "Other", "zero_codes_insufficient_followup" + + if n_codes == 1: + return "Other", "exactly_one_code" + + first = parse_date(min(icd_dates)) + last = parse_date(max(icd_dates)) + span_days = (last - first).days + gap_days = (first - psg_date).days + + if span_days < MIN_SPAN_DAYS: + return "Other", "span_under_7days" + if gap_days < LOWER_BOUND_DAYS: + return "Other", "diagnosed_too_early" + if gap_days > UPPER_BOUND_DAYS: + return "Other", "diagnosed_too_late" + + return "Positive", None + + +def run(demographics_path, icd_path, out_path): + icd_by_patient = load_icd_dates(icd_path) + + total = 0 + official_counts = defaultdict(int) + reason_counts = defaultdict(int) + mismatches = [] # patients whose current True/False label disagrees with official category + + with open(demographics_path) as f: + reader = csv.DictReader(f) + for row in reader: + total += 1 + pid = row["BDSPPatientID"] + site = row.get("SiteID", "") + psg_date = parse_date(row["CreationTime"]) + icd_dates = icd_by_patient.get(pid, []) + + ttlv_raw = row.get("Time_to_Last_Visit", "") + try: + ttlv_days = float(ttlv_raw) if ttlv_raw not in ("", "nan") else None + except ValueError: + ttlv_days = None + + category, reason = classify_official(psg_date, icd_dates, ttlv_days) + official_counts[category] += 1 + if reason: + reason_counts[reason] += 1 + + current_label = row.get("Cognitive_Impairment", "") + current_is_true = current_label == "True" + + # Flag any disagreement between what create_labels.py assigned + # and what the official 3-way rule would assign. + expected_true = (category == "Positive") + if current_is_true != expected_true or category == "Other": + mismatches.append({ + "BDSPPatientID": pid, + "SiteID": site, + "current_label": current_label, + "official_category": category, + "reason": reason or "", + "n_icd_codes": len(icd_dates), + "time_to_last_visit_days": ttlv_days if ttlv_days is not None else "", + }) + + # ── Report ──────────────────────────────────────────────────────────── + print(f"Total patients checked: {total}") + print(f"\nOfficial category breakdown:") + for cat in ("Positive", "Negative", "Other"): + n = official_counts.get(cat, 0) + pct = 100 * n / total if total else 0 + print(f" {cat:<10} {n:>6} ({pct:.1f}%)") + + if reason_counts: + print(f"\n'Other' breakdown:") + for reason, n in sorted(reason_counts.items(), key=lambda kv: -kv[1]): + print(f" {reason:<35} {n}") + + current_true = sum(1 for _ in open(demographics_path)) - 1 # not used directly; see below + n_other = official_counts.get("Other", 0) + + print(f"\n{'='*60}") + if n_other == 0: + print("CLEAN — no 'Other' patients found. This dataset appears to be") + print("pre-curated the same way the small training set was. Safe to") + print("train on the existing True/False labels as-is.") + else: + pct_of_total = 100 * n_other / total if total else 0 + print(f"NOT CLEAN — {n_other} patients ({pct_of_total:.1f}% of the dataset)") + print("are 'Other' by the official definition but are currently") + print("labeled False (Negative) in this demographics.csv.") + print(f"\nRecommendation: exclude these {n_other} patients from training") + print("entirely (do not train on them as either class) so your training") + print("distribution matches what real scoring will actually contain.") + print(f"\nFull list written to: {out_path}") + print(f"{'='*60}") + + if mismatches: + fieldnames = ["BDSPPatientID", "SiteID", "current_label", "official_category", + "reason", "n_icd_codes", "time_to_last_visit_days"] + with open(out_path, "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames) + writer.writeheader() + writer.writerows(mismatches) + + return official_counts, mismatches + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--demographics", required=True) + parser.add_argument("--icd", required=True) + parser.add_argument("--out", default="label_audit_mismatches.csv") + args = parser.parse_args() + + run(args.demographics, args.icd, args.out) \ No newline at end of file