{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# XGBoost — hand-crafted time-series features (WORKSHOP VERSION, self-contained)\n#\n# Role in the workshop: this is the \"general / basic preprocessing\" track — plain\n# statistical + frequency-domain features computed directly on the raw 20-channel\n# EEG signal. No bipolar montage, no mu-law encoding, no other EEG-specific tricks.\n# (Compare against the 1D CNN notebook, which uses the same raw signal but adds\n# those domain-specific preprocessing steps — that's where \"domain knowledge\" is\n# meant to show up in this workshop.)\n#\n# Features per channel (7 total, x 19 channels = 133-dim feature vector):\n#   mean, variance, zero-crossing rate, delta/theta/alpha/beta band power\n#\n# No sample weighting: the workshop subset is already class-balanced\n# (100 rows/class train, 20 rows/class val), so there's no imbalance to correct\n# for. (Full-scale training on the real, imbalanced ~5,900-row clean dataset is\n# where sample weighting / 2-step training matter — see the canonical pipeline\n# in kaggle_upload/src for that version.)\n#\n# Fully self-contained — no external .py imports.","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, time, random\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom scipy.stats import skew, kurtosis\nimport xgboost as xgb\nfrom sklearn.metrics import f1_score\n\nIS_KAGGLE = os.path.exists('/kaggle')\nprint(f\"Environment: {'Kaggle' if IS_KAGGLE else 'Local'} | xgboost {xgb.__version__}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:26:57.425663Z","iopub.execute_input":"2026-07-30T15:26:57.426527Z","iopub.status.idle":"2026-07-30T15:27:00.781428Z","shell.execute_reply.started":"2026-07-30T15:26:57.426498Z","shell.execute_reply":"2026-07-30T15:27:00.780611Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Config ============\nSEED       = 42\nFS         = 200          # Hz\nWINDOW_LEN = 10_000        # samples = 50 s at 200 Hz\n\nXGB_PARAMS = {\n    \"objective\"        : \"multi:softprob\",\n    \"num_class\"        : 6,\n    \"eval_metric\"      : \"mlogloss\",\n    \"tree_method\"      : \"hist\",\n    \"max_depth\"        : 6,\n    \"learning_rate\"    : 0.05,\n    \"subsample\"        : 0.8,\n    \"colsample_bytree\" : 0.8,\n    \"min_child_weight\" : 5,\n    \"seed\"             : SEED,\n}\nNUM_BOOST_ROUND        = 500\nEARLY_STOPPING_ROUNDS  = 30\n\nCHANNELS    = ['Fp1','F3','C3','P3','F7','T3','T5','O1','Fz','Cz','Pz',\n               'Fp2','F4','C4','P4','F8','T4','T6','O2']  # 19 EEG channels (EKG excluded)\nTARGETS     = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']\nCLASS_NAMES = ['Seizure', 'LPD', 'GPD', 'LRDA', 'GRDA', 'Other']\n\nrandom.seed(SEED)\nnp.random.seed(SEED)\n\nif IS_KAGGLE:\n    DATA_ROOT       = '/kaggle/input/competitions/hms-harmful-brain-activity-classification'\n    RAW_TRAIN_PATH  = os.path.join(DATA_ROOT, 'train.csv')\n    EEG_DIR         = os.path.join(DATA_ROOT, 'train_eegs')\n    SAMPLE_IDS_PATH = '/kaggle/input/datasets/xiaosufrankhu/midas-summer-academy-wk3-eeg/workshop_sample_ids.csv'\n    EEG_CACHE_DIR   = '/kaggle/working/eeg_cache_xgb'\nelse:\n    DATA_ROOT       = os.path.abspath('../')\n    RAW_TRAIN_PATH  = os.path.abspath('../data_raw/train.csv')\n    EEG_DIR         = os.path.join(DATA_ROOT, 'train_eegs')\n    SAMPLE_IDS_PATH = os.path.abspath('../data_raw/workshop_sample_ids.csv')\n    EEG_CACHE_DIR   = os.path.abspath('../eeg_cache_xgb')\n\n# NOTE: if running on Kaggle and this path doesn't exist, run !ls /kaggle/input\n# and update DATA_ROOT / EEG_DIR to match how the competition data was attached.\nprint(f\"Competition data path exists: {os.path.exists(DATA_ROOT)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:32:24.049895Z","iopub.execute_input":"2026-07-30T15:32:24.050268Z","iopub.status.idle":"2026-07-30T15:32:24.060057Z","shell.execute_reply.started":"2026-07-30T15:32:24.050223Z","shell.execute_reply":"2026-07-30T15:32:24.059196Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Data loading (workshop subset) ============\nraw_df     = pd.read_csv(RAW_TRAIN_PATH)\nsample_ids = pd.read_csv(SAMPLE_IDS_PATH)\n\nmerged = raw_df.merge(\n    sample_ids[[\"eeg_id\", \"eeg_sub_id\", \"split\"]],\n    on=[\"eeg_id\", \"eeg_sub_id\"],\n    how=\"inner\",\n)\ntrain_df = merged[merged[\"split\"] == \"train\"].reset_index(drop=True)\nval_df   = merged[merged[\"split\"] == \"val\"].reset_index(drop=True)\n\nprint(f'Train: {len(train_df):,} rows')\nprint(f'Val   : {len(val_df):,} rows')\nprint(train_df['expert_consensus'].value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:33:10.239158Z","iopub.execute_input":"2026-07-30T15:33:10.239500Z","iopub.status.idle":"2026-07-30T15:33:10.467342Z","shell.execute_reply.started":"2026-07-30T15:33:10.239478Z","shell.execute_reply":"2026-07-30T15:33:10.466458Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Raw EEG caching ============\n# Several rows may share the same eeg_id but use a different label window\n# (eeg_label_offset_seconds), so we cache each eeg_id's full raw recording once.\n\nos.makedirs(EEG_CACHE_DIR, exist_ok=True)\n\ndef cache_one_eeg(eeg_id: int, eeg_dir: str, cache_dir: str) -> None:\n    dst = os.path.join(cache_dir, f\"{eeg_id}.npy\")\n    if os.path.exists(dst):\n        return\n    src = os.path.join(eeg_dir, f\"{eeg_id}.parquet\")\n    eeg = pd.read_parquet(src, columns=CHANNELS).to_numpy(dtype=np.float32)  # (T, 19)\n    np.save(dst, eeg)\n\nneeded_eeg_ids = set(train_df['eeg_id']).union(val_df['eeg_id'])\nfor eeg_id in needed_eeg_ids:\n    cache_one_eeg(int(eeg_id), EEG_DIR, EEG_CACHE_DIR)\n\nprint(f'Cache ready: {len(needed_eeg_ids)} raw EEG recordings for workshop subset')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:34:12.301381Z","iopub.execute_input":"2026-07-30T15:34:12.301791Z","iopub.status.idle":"2026-07-30T15:34:19.038177Z","shell.execute_reply.started":"2026-07-30T15:34:12.301758Z","shell.execute_reply":"2026-07-30T15:34:19.037108Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Feature extraction (general / basic — no montage, no mu-law) ============\n# 7 features x 19 raw channels = 133-dim feature vector per row.\n\n_FREQS = np.fft.rfftfreq(WINDOW_LEN, d=1.0 / FS)\n_DELTA = (_FREQS >= 0.5)  & (_FREQS <  4.0)\n_THETA = (_FREQS >= 4.0)  & (_FREQS <  8.0)\n_ALPHA = (_FREQS >= 8.0)  & (_FREQS < 13.0)\n_BETA  = (_FREQS >= 13.0) & (_FREQS < 30.0)\n\n\ndef zero_crossing_rate(sig: np.ndarray) -> np.ndarray:\n    \"\"\"sig: (T, C) -> (C,). Fraction of adjacent-sample sign changes.\"\"\"\n    signs = np.sign(sig)\n    signs[signs == 0] = 1  # treat exact zero as positive, avoids spurious crossings\n    crossings = (np.diff(signs, axis=0) != 0).sum(axis=0)\n    return crossings / (len(sig) - 1)\n\n\ndef extract_features(window: np.ndarray) -> np.ndarray:\n    \"\"\"window: (T, 19) raw, unmontaged EEG samples -> (133,) feature vector.\"\"\"\n    window = np.nan_to_num(window, nan=0.0, posinf=0.0, neginf=0.0)\n\n    feat_mean = window.mean(axis=0)                                   # (19,)\n    feat_var  = window.var(axis=0)                                    # (19,)\n    feat_zcr  = zero_crossing_rate(window)                             # (19,)\n\n    fft_power = (np.abs(np.fft.rfft(window, axis=0)) ** 2) / len(window)  # (F, 19)\n    feat_delta = fft_power[_DELTA].sum(axis=0)\n    feat_theta = fft_power[_THETA].sum(axis=0)\n    feat_alpha = fft_power[_ALPHA].sum(axis=0)\n    # (beta band folded in as a 4th band-power feature, matching \"band power\" in the design doc)\n    feat_beta  = fft_power[_BETA].sum(axis=0)\n\n    return np.concatenate([\n        feat_mean, feat_var, feat_zcr, feat_delta, feat_theta, feat_alpha, feat_beta\n    ])  # 7 x 19 = 133 features\n\n\ndef build_features(df: pd.DataFrame, cache_dir: str, window_len: int = WINDOW_LEN):\n    X, y_hard, y_soft = [], [], []\n    for _, row in df.iterrows():\n        eeg_id = int(row['eeg_id'])\n        offset = int(row['eeg_label_offset_seconds']) * FS\n\n        eeg = np.load(os.path.join(cache_dir, f\"{eeg_id}.npy\"))  # (T, 19)\n        window = eeg[offset:offset + window_len]\n        if len(window) < window_len:\n            pad_top = window_len - len(window)\n            window  = np.pad(window, ((0, pad_top), (0, 0)), mode=\"constant\")\n\n        feats = extract_features(window)\n        soft  = row[TARGETS].to_numpy(dtype=np.float32)\n        soft  = soft / soft.sum() if soft.sum() > 0 else np.full(6, 1 / 6)\n\n        X.append(feats)\n        y_hard.append(int(np.argmax(soft)))\n        y_soft.append(soft)\n\n    return (np.array(X, dtype=np.float32),\n            np.array(y_hard),\n            np.array(y_soft, dtype=np.float32))\n\n\nt0 = time.time()\nX_train, y_train, y_train_soft = build_features(train_df, EEG_CACHE_DIR)\nX_val,   y_val,   y_val_soft   = build_features(val_df,   EEG_CACHE_DIR)\nprint(f\"X_train {X_train.shape} | X_val {X_val.shape}  [{time.time()-t0:.0f}s]\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:40:31.750331Z","iopub.execute_input":"2026-07-30T15:40:31.751395Z","iopub.status.idle":"2026-07-30T15:40:33.897436Z","shell.execute_reply.started":"2026-07-30T15:40:31.751327Z","shell.execute_reply":"2026-07-30T15:40:33.896768Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Train XGBoost ============\ndtrain = xgb.DMatrix(X_train, label=y_train)\ndval   = xgb.DMatrix(X_val,   label=y_val)\n\nevals_result = {}\nmodel = xgb.train(\n    params                = XGB_PARAMS,\n    dtrain                = dtrain,\n    num_boost_round       = NUM_BOOST_ROUND,\n    evals                 = [(dtrain, \"train\"), (dval, \"val\")],\n    early_stopping_rounds = EARLY_STOPPING_ROUNDS,\n    evals_result          = evals_result,\n    verbose_eval          = 25,\n)\n\nn_rounds_trained = len(evals_result[\"val\"][\"mlogloss\"])\nprint(f\"\\nRounds trained         : {n_rounds_trained}\")\nprint(f\"Best round (mlogloss)  : {model.best_iteration}  (val mlogloss={model.best_score:.4f})\")\n\n# ---- Stage val KL / F1 every 5 rounds ----\n# mlogloss is a convenient proxy XGBoost can optimize round-by-round, but KL divergence\n# is the metric we actually care about. They usually track each other but are not\n# guaranteed to peak at the same round — staging KL/F1 lets us pick the round that is\n# best on the metric that matters, and shows learners when/why the two disagree.\nSTAGE_EVERY  = 5\nstage_rounds = list(range(STAGE_EVERY, n_rounds_trained + 1, STAGE_EVERY))\nif not stage_rounds or stage_rounds[-1] != n_rounds_trained:\n    stage_rounds.append(n_rounds_trained)\n\nstage_kl, stage_f1 = [], []\nfor r in stage_rounds:\n    prob_r = model.predict(dval, iteration_range=(0, r)).reshape(-1, 6)\n    kl_r   = (y_val_soft * np.log(np.clip(y_val_soft, 1e-7, 1) / np.clip(prob_r, 1e-7, 1))).sum(axis=1).mean()\n    f1_r   = f1_score(y_val, prob_r.argmax(axis=1), average=\"macro\", zero_division=0)\n    stage_kl.append(kl_r)\n    stage_f1.append(f1_r)\n\nbest_idx   = int(np.argmin(stage_kl))\nBEST_ROUND = stage_rounds[best_idx]  # <- used for all \"final\" metrics/plots below\nprint(f\"Best round (val KL)    : {BEST_ROUND}  (val KL={stage_kl[best_idx]:.4f})\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:41:08.858494Z","iopub.execute_input":"2026-07-30T15:41:08.859470Z","iopub.status.idle":"2026-07-30T15:41:10.693187Z","shell.execute_reply.started":"2026-07-30T15:41:08.859444Z","shell.execute_reply":"2026-07-30T15:41:10.692675Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============ Metrics (reported at the KL-optimal round) ============\nprob_val = model.predict(dval, iteration_range=(0, BEST_ROUND)).reshape(-1, 6)\npred_val = prob_val.argmax(axis=1)\n\nmacro_f1  = f1_score(y_val, pred_val, average=\"macro\", zero_division=0)\nper_class = f1_score(y_val, pred_val, average=None,    zero_division=0)\n\nkl_per_sample = (y_val_soft * np.log(np.clip(y_val_soft, 1e-7, 1) / np.clip(prob_val, 1e-7, 1))).sum(axis=1)\nval_kl = kl_per_sample.mean()\n\n# Also report metrics at the mlogloss-best round, so \"loss-based\" and \"KL-based\"\n# selection can be compared side by side (they don't have to agree).\nprob_val_loss = model.predict(dval, iteration_range=(0, model.best_iteration + 1)).reshape(-1, 6)\npred_val_loss = prob_val_loss.argmax(axis=1)\nloss_f1 = f1_score(y_val, pred_val_loss, average=\"macro\", zero_division=0)\nloss_kl = (y_val_soft * np.log(np.clip(y_val_soft, 1e-7, 1) / np.clip(prob_val_loss, 1e-7, 1))).sum(axis=1).mean()\n\nprint(f\"--- selected by mlogloss (round {model.best_iteration}) ---\")\nprint(f\"Val macro F1     : {loss_f1:.4f}\")\nprint(f\"Val KL divergence: {loss_kl:.4f}\")\nprint(f\"\\n--- selected by val KL (round {BEST_ROUND}) ---\")\nprint(f\"Val macro F1     : {macro_f1:.4f}\")\nprint(f\"Val KL divergence: {val_kl:.4f}\")\nprint(f\"(random-guess baseline for 6 balanced classes: macro F1 \\u2248 0.167)\")\nprint(\"\\nPer-class F1 (at the val-KL-selected round):\")\nfor name, f in zip(CLASS_NAMES, per_class):\n    print(f\"  {name:<10} {f:.4f}\")\n\nfig, axes = plt.subplots(2, 2, figsize=(13, 9))\n\nrounds_all = range(len(evals_result[\"train\"][\"mlogloss\"]))\naxes[0, 0].plot(rounds_all, evals_result[\"train\"][\"mlogloss\"], label=\"train\")\naxes[0, 0].plot(rounds_all, evals_result[\"val\"][\"mlogloss\"],   label=\"val\")\naxes[0, 0].axvline(model.best_iteration, color=\"gray\", linestyle=\"--\",\n                    label=f\"mlogloss-best={model.best_iteration}\")\naxes[0, 0].axvline(BEST_ROUND, color=\"red\", linestyle=\"--\", label=f\"KL-best={BEST_ROUND}\")\naxes[0, 0].set_xlabel(\"Round\"); axes[0, 0].set_ylabel(\"mlogloss\")\naxes[0, 0].set_title(\"XGBoost mlogloss\"); axes[0, 0].legend(fontsize=8)\n\naxes[0, 1].plot(stage_rounds, stage_kl, marker=\"o\", color=\"darkorange\")\naxes[0, 1].axvline(BEST_ROUND, color=\"red\", linestyle=\"--\", label=f\"best={BEST_ROUND}\")\naxes[0, 1].set_xlabel(\"Round\"); axes[0, 1].set_ylabel(\"Val KL divergence\")\naxes[0, 1].set_title(\"Val KL divergence (every 5 rounds)\"); axes[0, 1].legend(fontsize=8)\n\naxes[1, 0].plot(stage_rounds, stage_f1, marker=\"o\", color=\"seagreen\")\naxes[1, 0].axhline(1/6, color=\"gray\", linestyle=\"--\", linewidth=1, label=\"random guess\")\naxes[1, 0].axvline(BEST_ROUND, color=\"red\", linestyle=\"--\", label=f\"best={BEST_ROUND}\")\naxes[1, 0].set_xlabel(\"Round\"); axes[1, 0].set_ylabel(\"Val macro F1\")\naxes[1, 0].set_title(\"Val macro F1 (every 5 rounds)\"); axes[1, 0].legend(fontsize=8)\n\naxes[1, 1].bar(CLASS_NAMES, per_class, color=\"steelblue\")\naxes[1, 1].axhline(macro_f1, color=\"red\", linestyle=\"--\", label=f\"macro F1 = {macro_f1:.3f}\")\naxes[1, 1].set_ylim(0, 1); axes[1, 1].set_ylabel(\"F1\")\naxes[1, 1].set_title(f\"Per-class F1 (val, round {BEST_ROUND})\"); axes[1, 1].legend(fontsize=8)\n\nplt.tight_layout()\nplt.savefig(\"xgb_workshop_eval.png\", dpi=150)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-30T15:41:47.546664Z","iopub.execute_input":"2026-07-30T15:41:47.546951Z","iopub.status.idle":"2026-07-30T15:41:48.568411Z","shell.execute_reply.started":"2026-07-30T15:41:47.546931Z","shell.execute_reply":"2026-07-30T15:41:48.567520Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Tweak & Compare (~20 min)\n\nPick 1-2 changes below, rerun the cell, and log your result in the shared sheet.\n\n| Knob | Default | Try |\n|---|---|---|\n| `max_depth` | 6 | 3 / 9 |\n| `num_boost_round` | 500 | 100 / 600 |\n| feature subset | all bands | drop spectral power (mean/var/zero-crossing only) |\n\n**Log your run:** What you changed | Val F1 | Val KL | one-line observation\n","metadata":{}},{"cell_type":"code","source":"# ============ Tweak & Compare ============\nTWEAK_MAX_DEPTH       = 6      # <- change me (try 3 or 9)\nTWEAK_NUM_BOOST_ROUND = 500    # <- change me (try 100 or 600)\nTWEAK_DROP_SPECTRAL   = False  # <- change me (try True to drop spectral-power features)\n\n\ndef build_features_tweak(df, cache_dir, drop_spectral=False, window_len=WINDOW_LEN):\n    X, y_hard, y_soft = [], [], []\n    for _, row in df.iterrows():\n        eeg_id = int(row[\"eeg_id\"])\n        offset = int(row[\"eeg_label_offset_seconds\"]) * FS\n\n        eeg    = np.load(os.path.join(cache_dir, f\"{eeg_id}.npy\"))\n        window = eeg[offset:offset + window_len]\n        if len(window) < window_len:\n            window = np.pad(window, ((0, window_len - len(window)), (0, 0)), mode=\"constant\")\n        window = np.nan_to_num(window, nan=0.0, posinf=0.0, neginf=0.0)\n\n        feat_mean = window.mean(axis=0)\n        feat_var  = window.var(axis=0)\n        feat_zcr  = zero_crossing_rate(window)\n\n        if drop_spectral:\n            feats = np.concatenate([feat_mean, feat_var, feat_zcr])  # 3 x 19 = 57 features\n        else:\n            fft_power  = (np.abs(np.fft.rfft(window, axis=0)) ** 2) / len(window)\n            feat_delta = fft_power[_DELTA].sum(axis=0)\n            feat_theta = fft_power[_THETA].sum(axis=0)\n            feat_alpha = fft_power[_ALPHA].sum(axis=0)\n            feat_beta  = fft_power[_BETA].sum(axis=0)\n            feats = np.concatenate([feat_mean, feat_var, feat_zcr,\n                                     feat_delta, feat_theta, feat_alpha, feat_beta])  # 7 x 19 = 133\n\n        soft = row[TARGETS].to_numpy(dtype=np.float32)\n        soft = soft / soft.sum() if soft.sum() > 0 else np.full(6, 1 / 6)\n\n        X.append(feats)\n        y_hard.append(int(np.argmax(soft)))\n        y_soft.append(soft)\n\n    return (np.array(X, dtype=np.float32),\n            np.array(y_hard),\n            np.array(y_soft, dtype=np.float32))\n\n\nXt_train, yt_train, yt_train_soft = build_features_tweak(train_df, EEG_CACHE_DIR, TWEAK_DROP_SPECTRAL)\nXt_val,   yt_val,   yt_val_soft   = build_features_tweak(val_df,   EEG_CACHE_DIR, TWEAK_DROP_SPECTRAL)\n\ntweak_params = dict(XGB_PARAMS)\ntweak_params[\"max_depth\"] = TWEAK_MAX_DEPTH\n\ndtrain_t = xgb.DMatrix(Xt_train, label=yt_train)\ndval_t   = xgb.DMatrix(Xt_val,   label=yt_val)\n\nevals_result_t = {}\nmodel_t = xgb.train(\n    params                = tweak_params,\n    dtrain                = dtrain_t,\n    num_boost_round       = TWEAK_NUM_BOOST_ROUND,\n    evals                 = [(dtrain_t, \"train\"), (dval_t, \"val\")],\n    early_stopping_rounds = EARLY_STOPPING_ROUNDS,\n    evals_result          = evals_result_t,\n    verbose_eval          = False,\n)\n\n# stage val KL/F1 every 5 rounds, same as the baseline Metrics cell, so \"best round\"\n# is picked the same way — otherwise this comparison isn't apples-to-apples\nn_rounds_t   = len(evals_result_t[\"val\"][\"mlogloss\"])\nstage_every  = 5\nstage_rounds_t = list(range(stage_every, n_rounds_t + 1, stage_every))\nif not stage_rounds_t or stage_rounds_t[-1] != n_rounds_t:\n    stage_rounds_t.append(n_rounds_t)\n\nstage_kl_t = []\nfor r in stage_rounds_t:\n    p = model_t.predict(dval_t, iteration_range=(0, r)).reshape(-1, 6)\n    stage_kl_t.append((yt_val_soft * np.log(np.clip(yt_val_soft, 1e-7, 1) / np.clip(p, 1e-7, 1))).sum(axis=1).mean())\n\nkl_best_round_t = stage_rounds_t[int(np.argmin(stage_kl_t))]\n\n\ndef _f1_kl_at(model, dval, y_val, y_val_soft, iteration_range):\n    p = model.predict(dval, iteration_range=iteration_range).reshape(-1, 6)\n    f1 = f1_score(y_val, p.argmax(axis=1), average=\"macro\", zero_division=0)\n    kl = (y_val_soft * np.log(np.clip(y_val_soft, 1e-7, 1) / np.clip(p, 1e-7, 1))).sum(axis=1).mean()\n    return f1, kl\n\n\nloss_f1_t, loss_kl_t = _f1_kl_at(model_t, dval_t, yt_val, yt_val_soft, (0, model_t.best_iteration + 1))\nkl_f1_t,   kl_kl_t   = _f1_kl_at(model_t, dval_t, yt_val, yt_val_soft, (0, kl_best_round_t))\n\nprint(f\"max_depth={TWEAK_MAX_DEPTH} | num_boost_round={TWEAK_NUM_BOOST_ROUND} | drop_spectral={TWEAK_DROP_SPECTRAL}\")\nprint(f\"--- selected by mlogloss (round {model_t.best_iteration}) ---\")\nprint(f\"Val macro F1 : {loss_f1_t:.4f}\")\nprint(f\"Val KL       : {loss_kl_t:.4f}\")\nprint(f\"\\n--- selected by val KL (round {kl_best_round_t}) ---\")\nprint(f\"Val macro F1 : {kl_f1_t:.4f}\")\nprint(f\"Val KL       : {kl_kl_t:.4f}\")\nprint(f\"\\n(baseline for comparison: mlogloss-best F1={loss_f1:.4f}/KL={loss_kl:.4f} | \"\n      f\"val-KL-best F1={macro_f1:.4f}/KL={val_kl:.4f})\")\nprint(\"-> log this row in the shared sheet: what you changed / F1 / KL / one-line observation\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Done early? Extend (~10 min)\n\n**Task:** add ONE new engineered feature to the feature vector and see if it moves val F1 / KL.\nA good candidate: a per-channel **spectral peak frequency** (the frequency with the highest power\nin each channel) — different information from the band-power features we already compute.\n\n**Copy-paste prompt for Claude / ChatGPT:**\n\n> Here is my feature extraction function and the function that builds a feature matrix from a\n> dataframe: [paste both `extract_features` and `build_features` from earlier in this notebook].\n> Add one new feature: for each channel, the frequency (in Hz) with the highest power in the FFT,\n> using the `FS`, `WINDOW_LEN`, and `_FREQS` variables already defined in this notebook. Give me\n> updated versions of both functions — call them `extract_features_with_peak_freq` and\n> `build_features_with_peak_freq` — with the new feature concatenated onto the existing feature\n> vector and everything else unchanged.\n>\n> Then, using these new functions, build train/val feature matrices, and train an XGBoost model\n> with:\n> ```\n> XGB_PARAMS = {\n>     \"objective\": \"multi:softprob\", \"num_class\": 6, \"eval_metric\": \"mlogloss\",\n>     \"tree_method\": \"hist\", \"max_depth\": 6, \"learning_rate\": 0.05,\n>     \"subsample\": 0.8, \"colsample_bytree\": 0.8, \"min_child_weight\": 5, \"seed\": 42,\n> }\n> NUM_BOOST_ROUND = 500\n> EARLY_STOPPING_ROUNDS = 30\n> ```\n> Print the resulting Val F1 and Val KL so I can compare against the baseline.\n\nPaste the AI's code into the cell below, run it, then compare F1/KL against the baseline above.\n\n**Note**: AI-generated code from a prompt like this typically gets you 80-90% of the way there — it may not run as-is on the first try (missing an import, a shape mismatch, a variable name that doesn't quite match what's already defined in this notebook). That gap is expected, and closing it is part of the exercise: read the error message, check it against the code already defined above, and fix it yourself. If you get stuck for more than ~10 minutes, the reference solution cell below shows one way to close that gap.\n\n","metadata":{}},{"cell_type":"code","source":"# ============ Your extended feature function goes here ============\n# Paste the code Claude/ChatGPT gives you in response to the prompt above, e.g.:\n#\n# def extract_features_extended(window):\n#     ...\n#\n# then rerun build_features_tweak-style code using your new function and compare\n# F1/KL against the 133-feature baseline from the \"Metrics\" cell above.\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Reference solution (only look if you're stuck)\n","metadata":{}},{"cell_type":"code","source":"# ============ Reference solution ============\ndef extract_features_with_peak_freq(window: np.ndarray) -> np.ndarray:\n    \"\"\"Same as extract_features(), plus one new feature: per-channel peak frequency.\"\"\"\n    window = np.nan_to_num(window, nan=0.0, posinf=0.0, neginf=0.0)\n\n    feat_mean = window.mean(axis=0)\n    feat_var  = window.var(axis=0)\n    feat_zcr  = zero_crossing_rate(window)\n\n    fft_power  = (np.abs(np.fft.rfft(window, axis=0)) ** 2) / len(window)\n    feat_delta = fft_power[_DELTA].sum(axis=0)\n    feat_theta = fft_power[_THETA].sum(axis=0)\n    feat_alpha = fft_power[_ALPHA].sum(axis=0)\n    feat_beta  = fft_power[_BETA].sum(axis=0)\n\n    # NEW: peak frequency per channel\n    peak_idx  = fft_power.argmax(axis=0)\n    feat_peak = _FREQS[peak_idx]\n\n    return np.concatenate([feat_mean, feat_var, feat_zcr,\n                            feat_delta, feat_theta, feat_alpha, feat_beta,\n                            feat_peak])  # 8 x 19 = 152 features\n\n\ndef build_features_with_peak_freq(df, cache_dir, window_len=WINDOW_LEN):\n    X, y_hard, y_soft = [], [], []\n    for _, row in df.iterrows():\n        eeg_id = int(row[\"eeg_id\"])\n        offset = int(row[\"eeg_label_offset_seconds\"]) * FS\n        eeg    = np.load(os.path.join(cache_dir, f\"{eeg_id}.npy\"))\n        window = eeg[offset:offset + window_len]\n        if len(window) < window_len:\n            window = np.pad(window, ((0, window_len - len(window)), (0, 0)), mode=\"constant\")\n\n        feats = extract_features_with_peak_freq(window)\n        soft  = row[TARGETS].to_numpy(dtype=np.float32)\n        soft  = soft / soft.sum() if soft.sum() > 0 else np.full(6, 1 / 6)\n\n        X.append(feats); y_hard.append(int(np.argmax(soft))); y_soft.append(soft)\n    return (np.array(X, dtype=np.float32), np.array(y_hard), np.array(y_soft, dtype=np.float32))\n\n\nXp_train, yp_train, _ = build_features_with_peak_freq(train_df, EEG_CACHE_DIR)\nXp_val,   yp_val,   yp_val_soft = build_features_with_peak_freq(val_df, EEG_CACHE_DIR)\n\ndtrain_p = xgb.DMatrix(Xp_train, label=yp_train)\ndval_p   = xgb.DMatrix(Xp_val,   label=yp_val)\nmodel_p  = xgb.train(\n    params=XGB_PARAMS, dtrain=dtrain_p, num_boost_round=NUM_BOOST_ROUND,\n    evals=[(dtrain_p, \"train\"), (dval_p, \"val\")],\n    early_stopping_rounds=EARLY_STOPPING_ROUNDS, verbose_eval=False,\n)\n\n# evaluated at its own mlogloss-best round — same criterion as the baseline's loss_f1/loss_kl\n# numbers from the Metrics cell above, so this stays a fair, apples-to-apples comparison\nprob_p = model_p.predict(dval_p, iteration_range=(0, model_p.best_iteration + 1)).reshape(-1, 6)\nf1_p = f1_score(yp_val, prob_p.argmax(axis=1), average=\"macro\", zero_division=0)\nkl_p = (yp_val_soft * np.log(np.clip(yp_val_soft, 1e-7, 1) / np.clip(prob_p, 1e-7, 1))).sum(axis=1).mean()\n\nprint(\"133 baseline features + peak-freq (152 total), evaluated at its own mlogloss-best round:\")\nprint(f\"Val F1: {f1_p:.4f} | Val KL: {kl_p:.4f}\")\nprint(f\"(baseline at its own mlogloss-best, for a fair comparison: F1={loss_f1:.4f} | KL={loss_kl:.4f})\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}