{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\nModule 0: Data Engineering & Sanity Checks\nM3CV EEG Biometric Authentication Pipeline\n\nRun steps in order via the __main__ block at the bottom.\nRequires: pip install mne librosa scipy numpy pandas --break-system-packages\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport scipy.io as sio\nimport mne\nimport librosa\n\n# ----------------------------------------------------------------------\n# CONFIG — edit these paths to match your environment\n# ----------------------------------------------------------------------\nBASE_DIR = \"/kaggle/input/competitions/eeg-biometric-competition\"\n\nENROLLMENT_CSV = os.path.join(BASE_DIR, \"Enrollment_Info.csv\")\nCALIBRATION_CSV = os.path.join(BASE_DIR, \"Calibration_Info.csv\")\nTESTING_CSV = os.path.join(BASE_DIR, \"Testing_Info.csv\")\n\nENROLLMENT_DIR = os.path.join(BASE_DIR, \"Enrollment\")\nCALIBRATION_DIR = os.path.join(BASE_DIR, \"Calibration\")\nTESTING_DIR = os.path.join(BASE_DIR, \"Testing\")\n\nOUT_DIR = \"./module0_outputs\"\nCACHE_DIR = os.path.join(OUT_DIR, \"spectrogram_cache\")\nos.makedirs(OUT_DIR, exist_ok=True)\nos.makedirs(CACHE_DIR, exist_ok=True)\n\n# Sampling params (from base paper / M3CV documentation)\nORIG_SFREQ = 1000 # Hz, as recorded\nDOWNSAMPLE_SFREQ = 250 # Hz, after downsampling\nN_CHANNELS = 65 # electrodes in each epoch\nEPOCH_LEN = 1000 # samples per epoch pre-downsampling\nN_MELS = 128\n\n# Generic channel names, only needed for the optional MNE-based downsample step\n# in preprocess_epoch(). No montage/positions required since we no longer do\n# spatial interpolation or referencing (see Step 5 note below on why).\nCHANNEL_NAMES = [f\"CH{i+1}\" for i in range(N_CHANNELS)]\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Normalize CSV schemas\n# ----------------------------------------------------------------------\ndef normalize_csvs():\n \"\"\"Load the three info CSVs and rename columns to a shared schema:\n epoch_id, subject_id, session, task, usage. Saves normalized copies to OUT_DIR.\"\"\"\n\n enrollment = pd.read_csv(ENROLLMENT_CSV)\n enrollment = enrollment.rename(columns={\n \"EpochID\": \"epoch_id\",\n \"subject\": \"subject_id\",\n \"session\": \"session\",\n \"condition\": \"task\",\n \"usage\": \"usage\",\n })\n\n calibration = pd.read_csv(CALIBRATION_CSV)\n calibration = calibration.rename(columns={\n \"EpochID\": \"epoch_id\",\n \"SubjectID\": \"subject_id\",\n \"Session\": \"session\",\n \"Task\": \"task\",\n \"Usage\": \"usage\",\n })\n\n testing = pd.read_csv(TESTING_CSV)\n testing = testing.rename(columns={\n \"Row\": \"epoch_id\",\n \"SubjectID\": \"subject_id\",\n \"Session\": \"session\",\n \"Task\": \"task\",\n \"Usage\": \"usage\",\n })\n\n for name, df in [(\"enrollment\", enrollment), (\"calibration\", calibration), (\"testing\", testing)]:\n missing = {\"epoch_id\", \"subject_id\", \"session\", \"task\", \"usage\"} - set(df.columns)\n if missing:\n raise ValueError(f\"{name} CSV missing expected columns after rename: {missing}\")\n\n enrollment.to_csv(os.path.join(OUT_DIR, \"enrollment_normalized.csv\"), index=False)\n calibration.to_csv(os.path.join(OUT_DIR, \"calibration_normalized.csv\"), index=False)\n testing.to_csv(os.path.join(OUT_DIR, \"testing_normalized.csv\"), index=False)\n\n print(\"[normalize_csvs] Saved normalized CSVs to\", OUT_DIR)\n return enrollment, calibration, testing\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Verify task-code mapping\n# ----------------------------------------------------------------------\ndef verify_task_distribution(enrollment, calibration, testing):\n \"\"\"Print unique task codes and their counts per split, so you can cross-check\n against M3CV documentation before assigning task codes to cascade stages.\"\"\"\n\n for name, df in [(\"enrollment\", enrollment), (\"calibration\", calibration), (\"testing\", testing)]:\n print(f\"\\n[verify_task_distribution] {name} — unique task codes and counts:\")\n print(df[\"task\"].value_counts().sort_index())\n\n task_summary = {\n \"enrollment_tasks\": sorted(enrollment[\"task\"].unique().tolist()),\n \"calibration_tasks\": sorted(calibration[\"task\"].unique().tolist()),\n \"testing_tasks\": sorted(testing[\"task\"].unique().tolist()),\n }\n with open(os.path.join(OUT_DIR, \"task_distribution.json\"), \"w\") as f:\n json.dump(task_summary, f, indent=2)\n\n print(\"\\n[verify_task_distribution] Saved task_distribution.json — \"\n \"confirm task code → paradigm mapping against M3CV docs before Module 2.\")\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Identify the cross-session verification cohort\n# ----------------------------------------------------------------------\ndef get_verification_cohort(enrollment, calibration):\n \"\"\"Return the sorted list of subject IDs present in BOTH Enrollment (S1) and\n Calibration (S2) — this is the fixed cohort used for Module 1's cross-session\n verification protocol. Cache to disk so every downstream module reuses the\n same list instead of recomputing it.\"\"\"\n\n enrollment_subjects = set(enrollment[\"subject_id\"].unique())\n calibration_subjects = set(calibration[\"subject_id\"].unique())\n cohort = sorted(enrollment_subjects & calibration_subjects)\n\n print(f\"[get_verification_cohort] Enrollment subjects: {len(enrollment_subjects)}\")\n print(f\"[get_verification_cohort] Calibration subjects: {len(calibration_subjects)}\")\n print(f\"[get_verification_cohort] Overlap (verification cohort): {len(cohort)}\")\n\n with open(os.path.join(OUT_DIR, \"verification_cohort.json\"), \"w\") as f:\n json.dump(cohort, f, indent=2)\n\n return cohort\n\n\ndef cap_epochs_per_task(df, max_per_task=2000, random_state=42):\n \"\"\"Cap the number of epochs per task to max_per_task, sampled proportionally\n by subject_id within each task so no single subject dominates the capped\n set. Tasks already at or under max_per_task are left untouched. Applies to\n whatever dataframe is passed in — call this AFTER filtering to your\n verification cohort, so counts reflect what's actually being cached.\"\"\"\n\n capped_parts = []\n\n for task_val, task_df in df.groupby(\"task\"):\n if len(task_df) <= max_per_task:\n capped_parts.append(task_df)\n continue\n\n frac = max_per_task / len(task_df)\n\n def _sample_subject(g):\n n = max(1, round(len(g) * frac))\n n = min(n, len(g))\n return g.sample(n=n, random_state=random_state)\n\n sampled = task_df.groupby(\"subject_id\", group_keys=False).apply(_sample_subject)\n\n # Rounding can leave sampled slightly off max_per_task — trim or top up\n if len(sampled) > max_per_task:\n sampled = sampled.sample(n=max_per_task, random_state=random_state)\n elif len(sampled) < max_per_task:\n remaining = task_df.drop(sampled.index)\n extra_needed = max_per_task - len(sampled)\n if len(remaining) > 0:\n extra = remaining.sample(n=min(extra_needed, len(remaining)), random_state=random_state)\n sampled = pd.concat([sampled, extra])\n\n capped_parts.append(sampled)\n\n result = pd.concat(capped_parts).reset_index(drop=True)\n\n print(\"[cap_epochs_per_task] Before capping:\")\n print(df[\"task\"].value_counts().sort_index())\n print(\"[cap_epochs_per_task] After capping (max_per_task=\"\n f\"{max_per_task}), stratified by subject:\")\n print(result[\"task\"].value_counts().sort_index())\n\n return result\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Sanity-check .mat files\n# ----------------------------------------------------------------------\ndef load_epoch_mat(epoch_id, split_dir):\n \"\"\"Load a single epoch's .mat file and return the raw (65, 1000) array.\n Handles the common case where the array is stored under a non-obvious key.\"\"\"\n\n fname = f\"epoch{int(epoch_id):06d}.mat\" if str(epoch_id).isdigit() else f\"{epoch_id}.mat\"\n fpath = os.path.join(split_dir, fname)\n\n mat = sio.loadmat(fpath)\n data_keys = [k for k in mat.keys() if not k.startswith(\"__\")]\n if len(data_keys) != 1:\n raise ValueError(f\"Ambiguous keys in {fpath}: {data_keys}\")\n\n arr = mat[data_keys[0]]\n return np.asarray(arr)\n\n\ndef sanity_check_mats(enrollment, calibration, cohort, n_samples=5):\n \"\"\"Load a handful of .mat files from Enrollment and Calibration for cohort\n subjects, and confirm shape == (65, 1000) with no NaNs. Raises loudly if not,\n since a shape mismatch here silently breaks every downstream module.\"\"\"\n\n checks = [\n (\"enrollment\", enrollment, ENROLLMENT_DIR),\n (\"calibration\", calibration, CALIBRATION_DIR),\n ]\n\n for name, df, split_dir in checks:\n sub_df = df[df[\"subject_id\"].isin(cohort)]\n sample_ids = sub_df[\"epoch_id\"].sample(min(n_samples, len(sub_df)), random_state=0)\n\n for eid in sample_ids:\n arr = load_epoch_mat(eid, split_dir)\n assert arr.shape == (N_CHANNELS, EPOCH_LEN), (\n f\"[sanity_check_mats] {name} epoch {eid} has shape {arr.shape}, \"\n f\"expected ({N_CHANNELS}, {EPOCH_LEN})\"\n )\n assert not np.isnan(arr).any(), f\"[sanity_check_mats] {name} epoch {eid} contains NaNs\"\n\n print(f\"[sanity_check_mats] {name}: {len(sample_ids)} epochs OK, shape ({N_CHANNELS}, {EPOCH_LEN})\")\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Preprocessing pipeline\n#\n# IMPORTANT: M3CV's released .mat epochs are ALREADY bad-channel-repaired,\n# TP9/TP10-referenced, bandpass-filtered (0.01-200Hz), and ICA-cleaned by the\n# dataset creators (base paper, Section 4.2.2, describes this as work \"done\n# by them\" prior to release). Redoing these steps on isolated 1-second/\n# 1000-sample epochs is both unnecessary and statistically unsound (ICA and\n# a 0.01Hz highpass both assume much longer continuous recordings than a\n# single 1000-sample epoch provides). This pipeline therefore performs only\n# the steps that are actually novel to the base paper / your contributions:\n# optional downsampling, flattening, and mel-spectrogram conversion.\n#\n# If you later obtain RAW, unprocessed EEG (not M3CV's released epochs),\n# reinstate bad-channel repair / referencing / bandpass / ICA on the\n# CONTINUOUS recording before epoching -- not per-epoch.\n# ----------------------------------------------------------------------\ndef preprocess_epoch(raw_array, do_downsample=True):\n \"\"\"Minimal preprocessing for one (65, 1000) M3CV epoch array (already\n cleaned upstream). Optionally downsamples to DOWNSAMPLE_SFREQ. Returns a\n numpy array of shape (n_channels, new_len).\"\"\"\n\n if not do_downsample:\n return raw_array\n\n info = mne.create_info(ch_names=CHANNEL_NAMES, sfreq=ORIG_SFREQ, ch_types=\"eeg\")\n raw = mne.io.RawArray(raw_array, info, verbose=False)\n raw.resample(DOWNSAMPLE_SFREQ, verbose=False)\n\n return raw.get_data()\n\n\ndef signal_to_mel_spectrogram(signal_1d, sr=DOWNSAMPLE_SFREQ, n_mels=N_MELS):\n \"\"\"Convert a single flattened multi-channel signal into a mel spectrogram\n in dB, computed at the EEG signal's OWN native sample rate.\n\n NOTE: earlier versions of this function resampled to an audio-style rate\n (44100 Hz) before computing the spectrogram, following the base paper's\n librosa convention literally. That is wrong for this data: upsampling a\n 250Hz-derived signal to 44100Hz manufactures ~176x more samples with zero\n added information, producing multi-megabyte spectrograms per epoch and\n filling disk before the full cohort finishes caching. Computing the mel\n spectrogram directly at sr (the signal's real rate) avoids that entirely\n and produces a spectrogram at the actual information density of the data.\n\n n_fft/hop_length are sized down from librosa's audio-scale defaults\n (n_fft=2048, hop_length=512) to fit a ~16,000-sample epoch-length signal\n instead of a multi-second audio clip.\n \"\"\"\n\n n_fft = min(256, len(signal_1d))\n hop_length = max(n_fft // 4, 1)\n\n mel = librosa.feature.melspectrogram(\n y=signal_1d.astype(np.float32), sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels\n )\n mel_db = librosa.power_to_db(mel, ref=np.max)\n return mel_db.astype(np.float32)\n\n\ndef estimate_cache_size(enrollment, calibration, cohort, n_probe=10):\n \"\"\"Process a handful of real epochs end-to-end to measure actual on-disk\n spectrogram size, then extrapolate to the full cohort. Run this BEFORE\n process_and_cache_cohort so a disk-full crash doesn't happen mid-run.\"\"\"\n\n import shutil\n\n probe_rows = enrollment[enrollment[\"subject_id\"].isin(cohort)].sample(\n min(n_probe, len(enrollment)), random_state=RANDOM_STATE if \"RANDOM_STATE\" in globals() else 42\n )\n\n sizes = []\n for _, row in probe_rows.iterrows():\n raw_array = load_epoch_mat(row[\"epoch_id\"], ENROLLMENT_DIR)\n cleaned = preprocess_epoch(raw_array)\n flattened = cleaned.flatten()\n mel_db = signal_to_mel_spectrogram(flattened)\n sizes.append(mel_db.nbytes)\n\n avg_bytes = np.mean(sizes)\n total_epochs = len(enrollment[enrollment[\"subject_id\"].isin(cohort)]) + \\\n len(calibration[calibration[\"subject_id\"].isin(cohort)])\n projected_total_gb = (avg_bytes * total_epochs) / (1024 ** 3)\n\n free_bytes = shutil.disk_usage(OUT_DIR).free\n free_gb = free_bytes / (1024 ** 3)\n\n print(f\"[estimate_cache_size] Avg spectrogram size: {avg_bytes / 1024:.1f} KB\")\n print(f\"[estimate_cache_size] Projected total for {total_epochs} epochs: {projected_total_gb:.2f} GB\")\n print(f\"[estimate_cache_size] Free disk space available: {free_gb:.2f} GB\")\n\n if projected_total_gb > free_gb * 0.9:\n print(\"[estimate_cache_size] WARNING: projected cache size is close to or \"\n \"exceeds available disk space. Reduce N_MELS, increase hop_length \"\n \"in signal_to_mel_spectrogram, or free up disk before proceeding.\")\n else:\n print(\"[estimate_cache_size] OK: sufficient disk space for full cohort caching.\")\n\n return projected_total_gb, free_gb\n\n\ndef process_and_cache_cohort(enrollment, calibration, cohort):\n \"\"\"Run the full preprocessing + mel-spectrogram pipeline for every epoch\n belonging to the verification cohort, in both Enrollment and Calibration,\n and cache results as .npy files so no module reprocesses raw .mat files\n more than once. Wraps each epoch in try/except so a single failure (disk\n full, corrupt file) is logged and skipped rather than crashing the whole\n run and losing everything already cached.\"\"\"\n\n jobs = [\n (\"enrollment\", enrollment, ENROLLMENT_DIR),\n (\"calibration\", calibration, CALIBRATION_DIR),\n ]\n\n manifest = []\n failures = []\n\n for split_name, df, split_dir in jobs:\n sub_df = df[df[\"subject_id\"].isin(cohort)]\n\n for _, row in sub_df.iterrows():\n eid = row[\"epoch_id\"]\n out_path = os.path.join(CACHE_DIR, f\"{split_name}_{eid}.npy\")\n\n if os.path.exists(out_path):\n manifest.append({\"epoch_id\": eid, \"split\": split_name,\n \"subject_id\": row[\"subject_id\"], \"session\": row[\"session\"],\n \"task\": row[\"task\"], \"cache_path\": out_path})\n continue\n\n try:\n raw_array = load_epoch_mat(eid, split_dir)\n cleaned = preprocess_epoch(raw_array)\n flattened = cleaned.flatten()\n mel_db = signal_to_mel_spectrogram(flattened)\n\n np.save(out_path, mel_db)\n manifest.append({\"epoch_id\": eid, \"split\": split_name,\n \"subject_id\": row[\"subject_id\"], \"session\": row[\"session\"],\n \"task\": row[\"task\"], \"cache_path\": out_path})\n except OSError as e:\n print(f\"[process_and_cache_cohort] STOPPING: disk error on epoch {eid} ({split_name}): {e}\")\n print(f\"[process_and_cache_cohort] {len(manifest)} epochs cached successfully before failure. \"\n f\"Free disk space and rerun — already-cached epochs will be skipped automatically.\")\n manifest_df = pd.DataFrame(manifest)\n manifest_df.to_csv(os.path.join(OUT_DIR, \"spectrogram_manifest.csv\"), index=False)\n raise\n except Exception as e:\n print(f\"[process_and_cache_cohort] Skipping epoch {eid} ({split_name}) due to error: {e}\")\n failures.append({\"epoch_id\": eid, \"split\": split_name, \"error\": str(e)})\n continue\n\n print(f\"[process_and_cache_cohort] Finished {split_name}: {len(sub_df)} epochs\")\n\n manifest_df = pd.DataFrame(manifest)\n manifest_df.to_csv(os.path.join(OUT_DIR, \"spectrogram_manifest.csv\"), index=False)\n print(f\"[process_and_cache_cohort] Manifest saved with {len(manifest_df)} entries\")\n\n if failures:\n pd.DataFrame(failures).to_csv(os.path.join(OUT_DIR, \"cache_failures.csv\"), index=False)\n print(f\"[process_and_cache_cohort] {len(failures)} epochs failed (non-disk errors) — see cache_failures.csv\")\n\n return manifest_df\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 0: Data Engineering & Sanity Checks ===\\n\")\n\n enrollment, calibration, testing = normalize_csvs()\n verify_task_distribution(enrollment, calibration, testing)\n cohort = get_verification_cohort(enrollment, calibration)\n\n # Filter to cohort BEFORE capping — the printed task counts earlier were\n # over all 95 Enrollment subjects; capping needs to act on what's actually\n # going to be cached (the 20-subject cohort), not the full pool.\n enrollment_cohort = enrollment[enrollment[\"subject_id\"].isin(cohort)].reset_index(drop=True)\n calibration_cohort = calibration[calibration[\"subject_id\"].isin(cohort)].reset_index(drop=True)\n\n MAX_EPOCHS_PER_TASK = 2000\n print(\"\\n[main] Capping Enrollment epochs per task:\")\n enrollment_capped = cap_epochs_per_task(enrollment_cohort, max_per_task=MAX_EPOCHS_PER_TASK)\n print(\"\\n[main] Capping Calibration epochs per task:\")\n calibration_capped = cap_epochs_per_task(calibration_cohort, max_per_task=MAX_EPOCHS_PER_TASK)\n\n sanity_check_mats(enrollment_capped, calibration_capped, cohort, n_samples=5)\n\n # IMPORTANT: if you ran an earlier version of this script (before the\n # 44100Hz audio-rate resample bug was fixed), delete ./module0_outputs/\n # spectrogram_cache/*.npy first — those files are ~20-100x oversized and\n # mixing them with correctly-sized new spectrograms will produce\n # inconsistent shapes downstream in Module 1.\n estimate_cache_size(enrollment_capped, calibration_capped, cohort)\n\n manifest = process_and_cache_cohort(enrollment_capped, calibration_capped, cohort)\n\n print(\"\\n=== Module 0 complete ===\")\n print(f\"Outputs saved to: {OUT_DIR}\")\n print(f\"Spectrogram cache: {CACHE_DIR}\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# \"\"\"\n# Module 1: Cross-Session EEG Biometric Verification Protocol\n# M3CV EEG Biometric Authentication Pipeline\n\n# Depends on Module 0's outputs: ./module0_outputs/spectrogram_manifest.csv\n# and the cached .npy spectrograms in ./module0_outputs/spectrogram_cache/.\n\n# Run steps in order via the __main__ block at the bottom.\n# Requires: pip install tensorflow scikit-learn numpy pandas matplotlib --break-system-packages\n# \"\"\"\n\n# import os\n# import json\n# import numpy as np\n# import pandas as pd\n# import tensorflow as tf\n# from sklearn.model_selection import train_test_split\n# from sklearn.metrics import roc_curve, roc_auc_score\n# from sklearn.preprocessing import LabelEncoder\n# import matplotlib.pyplot as plt\n\n# # ----------------------------------------------------------------------\n# # CONFIG\n# # ----------------------------------------------------------------------\n# MODULE0_DIR = \"./module0_outputs\"\n# MANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\n\n# OUT_DIR = \"./module1_outputs\"\n# os.makedirs(OUT_DIR, exist_ok=True)\n\n# EMBEDDING_DIM = 64\n# BATCH_SIZE = 32\n# EPOCHS = 25\n# LEARNING_RATE = 0.001\n# RANDOM_STATE = 42\n\n# np.random.seed(RANDOM_STATE)\n# tf.random.set_seed(RANDOM_STATE)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 1: Load manifest and cached spectrograms\n# # ----------------------------------------------------------------------\n# def load_manifest():\n# \"\"\"Load Module 0's manifest, which lists every cached spectrogram with its\n# split (enrollment/calibration), subject_id, session, and task.\"\"\"\n\n# manifest = pd.read_csv(MANIFEST_PATH)\n# print(f\"[load_manifest] Loaded {len(manifest)} entries\")\n# print(manifest[\"split\"].value_counts())\n# return manifest\n\n\n# def load_spectrogram_array(cache_path, target_shape=None):\n# \"\"\"Load one cached spectrogram .npy and optionally pad/crop to a target\n# shape, since mel spectrogram time-frames can vary by a few frames\n# depending on epoch length after resampling.\"\"\"\n\n# arr = np.load(cache_path)\n\n# if target_shape is not None:\n# h, w = arr.shape\n# th, tw = target_shape\n# if h != th or w != tw:\n# padded = np.zeros((th, tw), dtype=arr.dtype)\n# h_use, w_use = min(h, th), min(w, tw)\n# padded[:h_use, :w_use] = arr[:h_use, :w_use]\n# arr = padded\n\n# return arr\n\n\n# def determine_target_shape(manifest, n_check=20):\n# \"\"\"Inspect a sample of cached spectrograms to find the most common shape,\n# which becomes the fixed input shape for the CNN. Flags shape mismatch\n# so silent padding/cropping doesn't hide a preprocessing bug.\"\"\"\n\n# sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n# shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n# shape_counts = pd.Series(shapes).value_counts()\n\n# print(\"[determine_target_shape] Shape distribution across sample:\")\n# print(shape_counts)\n\n# target_shape = shape_counts.idxmax()\n# if len(shape_counts) > 1:\n# print(f\"[determine_target_shape] WARNING: inconsistent shapes found — \"\n# f\"using most common {target_shape}, others will be padded/cropped\")\n\n# return target_shape\n\n\n# # ----------------------------------------------------------------------\n# # STEP 2: Build train/template dataset from Enrollment\n# # ----------------------------------------------------------------------\n# def prepare_enrollment_data(manifest, target_shape, val_split=0.1):\n# \"\"\"Split Enrollment (Session 1) epochs for the cohort into a train/val set\n# for the embedding-extractor CNN. All Enrollment epochs (train+val) are\n# later used to build per-subject templates.\"\"\"\n\n# enrollment = manifest[manifest[\"split\"] == \"enrollment\"].reset_index(drop=True)\n\n# le = LabelEncoder()\n# enrollment[\"label\"] = le.fit_transform(enrollment[\"subject_id\"])\n\n# train_df, val_df = train_test_split(\n# enrollment, test_size=val_split, stratify=enrollment[\"label\"], random_state=RANDOM_STATE\n# )\n\n# print(f\"[prepare_enrollment_data] Train: {len(train_df)}, Val: {len(val_df)}, \"\n# f\"Subjects: {enrollment['subject_id'].nunique()}\")\n\n# return enrollment, train_df, val_df, le\n\n\n# def make_dataset(df, target_shape, batch_size=BATCH_SIZE, shuffle=True):\n# \"\"\"Build a tf.data.Dataset that loads cached spectrograms on the fly and\n# yields (spectrogram, label) pairs for training.\"\"\"\n\n# paths = df[\"cache_path\"].values\n# labels = df[\"label\"].values\n\n# def gen():\n# for p, y in zip(paths, labels):\n# arr = load_spectrogram_array(p, target_shape)\n# arr = arr[..., np.newaxis].astype(np.float32) # add channel dim\n# yield arr, y\n\n# ds = tf.data.Dataset.from_generator(\n# gen,\n# output_signature=(\n# tf.TensorSpec(shape=(*target_shape, 1), dtype=tf.float32),\n# tf.TensorSpec(shape=(), dtype=tf.int64),\n# ),\n# )\n\n# if shuffle:\n# ds = ds.shuffle(buffer_size=len(df), seed=RANDOM_STATE)\n\n# return ds.batch(batch_size).prefetch(tf.data.AUTOTUNE)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 3: CNN embedding-extractor architecture\n# # ----------------------------------------------------------------------\n# def build_model(input_shape, n_classes, embedding_dim=EMBEDDING_DIM):\n# \"\"\"CNN architecture following the base paper's Conv-Pool-Conv-Pool-Dropout\n# -Dense structure (Figure 5), with an explicit embedding_dim-wide dense\n# layer whose output is later reused as the identity embedding for\n# verification (not just the final softmax classification).\"\"\"\n\n# inputs = tf.keras.Input(shape=input_shape)\n\n# x = tf.keras.layers.Conv2D(32, 3, activation=\"relu\")(inputs)\n# x = tf.keras.layers.MaxPooling2D(2)(x)\n# x = tf.keras.layers.Conv2D(64, 3, activation=\"relu\")(x)\n# x = tf.keras.layers.MaxPooling2D(2)(x)\n# x = tf.keras.layers.Dropout(0.2)(x)\n# x = tf.keras.layers.Flatten()(x)\n# x = tf.keras.layers.Dense(128, activation=\"relu\")(x)\n\n# embedding = tf.keras.layers.Dense(embedding_dim, activation=\"relu\", name=\"embedding\")(x)\n# outputs = tf.keras.layers.Dense(n_classes, activation=\"softmax\", name=\"classifier\")(embedding)\n\n# model = tf.keras.Model(inputs=inputs, outputs=outputs)\n# embedding_model = tf.keras.Model(inputs=inputs, outputs=embedding)\n\n# model.compile(\n# optimizer=tf.keras.optimizers.Adam(learning_rate=LEARNING_RATE),\n# loss=\"sparse_categorical_crossentropy\",\n# metrics=[\"accuracy\"],\n# )\n\n# return model, embedding_model\n\n\n# def train_embedding_extractor(train_df, val_df, target_shape, n_classes):\n# \"\"\"Train the CNN as a 20-way subject classifier purely to learn useful\n# features — the trained embedding sub-model (penultimate dense layer) is\n# what Module 1 actually uses, not the softmax output.\"\"\"\n\n# train_ds = make_dataset(train_df, target_shape, shuffle=True)\n# val_ds = make_dataset(val_df, target_shape, shuffle=False)\n\n# model, embedding_model = build_model((*target_shape, 1), n_classes)\n\n# history = model.fit(\n# train_ds,\n# validation_data=val_ds,\n# epochs=EPOCHS,\n# verbose=1,\n# )\n\n# model.save(os.path.join(OUT_DIR, \"full_classifier.keras\"))\n# embedding_model.save(os.path.join(OUT_DIR, \"embedding_model.keras\"))\n\n# with open(os.path.join(OUT_DIR, \"training_history.json\"), \"w\") as f:\n# json.dump({k: [float(v) for v in vals] for k, vals in history.history.items()}, f, indent=2)\n\n# return embedding_model\n\n\n# # ----------------------------------------------------------------------\n# # STEP 4: Extract embeddings and build per-subject templates\n# # ----------------------------------------------------------------------\n# def extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n# \"\"\"Run every epoch in df through the trained embedding model, returning\n# an (n_epochs, embedding_dim) array aligned with df's row order.\"\"\"\n\n# ds = make_dataset(df.assign(label=0), target_shape, batch_size=batch_size, shuffle=False)\n# # label is a dummy placeholder here since we only need embeddings, not loss\n\n# embeddings = []\n# for batch_x, _ in ds:\n# emb = embedding_model.predict(batch_x, verbose=0)\n# embeddings.append(emb)\n\n# return np.concatenate(embeddings, axis=0)\n\n\n# def build_subject_templates(embeddings, subject_ids):\n# \"\"\"Average all of a subject's Enrollment embeddings into one template\n# vector per subject. Returns a dict {subject_id: template_vector}.\"\"\"\n\n# templates = {}\n# for sid in np.unique(subject_ids):\n# mask = subject_ids == sid\n# templates[sid] = embeddings[mask].mean(axis=0)\n\n# print(f\"[build_subject_templates] Built {len(templates)} subject templates\")\n# return templates\n\n\n# # ----------------------------------------------------------------------\n# # STEP 5: Genuine/impostor pair generation and scoring\n# # ----------------------------------------------------------------------\n# def cosine_similarity(a, b):\n# return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# def generate_verification_scores(probe_embeddings, probe_subject_ids, templates):\n# \"\"\"For every probe epoch, compute similarity against its OWN template\n# (genuine trial) and against every OTHER subject's template (impostor\n# trials). Returns two arrays: genuine_scores, impostor_scores.\"\"\"\n\n# subject_list = list(templates.keys())\n# genuine_scores, impostor_scores = [], []\n\n# for emb, true_sid in zip(probe_embeddings, probe_subject_ids):\n# for sid in subject_list:\n# score = cosine_similarity(emb, templates[sid])\n# if sid == true_sid:\n# genuine_scores.append(score)\n# else:\n# impostor_scores.append(score)\n\n# return np.array(genuine_scores), np.array(impostor_scores)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 6: FAR / FRR / EER / ROC-AUC\n# # ----------------------------------------------------------------------\n# def compute_verification_metrics(genuine_scores, impostor_scores, label=\"\"):\n# \"\"\"Sweep thresholds over the combined score range to compute FAR/FRR at\n# each threshold, find EER (where FAR == FRR), and compute ROC-AUC using\n# genuine=1 / impostor=0 labels.\"\"\"\n\n# y_true = np.concatenate([np.ones_like(genuine_scores), np.zeros_like(impostor_scores)])\n# y_score = np.concatenate([genuine_scores, impostor_scores])\n\n# fpr, tpr, thresholds = roc_curve(y_true, y_score)\n# fnr = 1 - tpr # FRR\n# far = fpr\n\n# eer_idx = np.nanargmin(np.abs(far - fnr))\n# eer = (far[eer_idx] + fnr[eer_idx]) / 2\n# eer_threshold = thresholds[eer_idx]\n\n# auc = roc_auc_score(y_true, y_score)\n\n# print(f\"\\n[compute_verification_metrics] {label}\")\n# print(f\" EER: {eer:.4f} at threshold {eer_threshold:.4f}\")\n# print(f\" ROC-AUC: {auc:.4f}\")\n# print(f\" Genuine scores: mean={genuine_scores.mean():.4f}, std={genuine_scores.std():.4f}\")\n# print(f\" Impostor scores: mean={impostor_scores.mean():.4f}, std={impostor_scores.std():.4f}\")\n\n# return {\n# \"label\": label,\n# \"eer\": float(eer),\n# \"eer_threshold\": float(eer_threshold),\n# \"roc_auc\": float(auc),\n# \"fpr\": fpr.tolist(),\n# \"tpr\": tpr.tolist(),\n# \"genuine_mean\": float(genuine_scores.mean()),\n# \"impostor_mean\": float(impostor_scores.mean()),\n# }\n\n\n# def plot_roc_curves(results_list, save_path):\n# plt.figure(figsize=(7, 6))\n# for res in results_list:\n# plt.plot(res[\"fpr\"], res[\"tpr\"], label=f\"{res['label']} (AUC={res['roc_auc']:.3f})\")\n# plt.plot([0, 1], [0, 1], \"k--\", alpha=0.3)\n# plt.xlabel(\"False Positive Rate (FAR)\")\n# plt.ylabel(\"True Positive Rate (1 - FRR)\")\n# plt.title(\"ROC Curve: Cross-Session vs Same-Session Verification\")\n# plt.legend()\n# plt.savefig(save_path, dpi=150, bbox_inches=\"tight\")\n# plt.close()\n# print(f\"[plot_roc_curves] Saved to {save_path}\")\n\n\n# # ----------------------------------------------------------------------\n# # STEP 7: Same-session baseline (for comparison against cross-session result)\n# # ----------------------------------------------------------------------\n# def run_same_session_baseline(enrollment_df, embedding_model, target_shape):\n# \"\"\"Reproduce a same-session verification test: split Enrollment epochs\n# into a template half and a probe half (random split, same session),\n# exactly the kind of split the base paper used for its classification\n# task. This is the number your cross-session result should be compared\n# against.\"\"\"\n\n# template_df, probe_df = train_test_split(\n# enrollment_df, test_size=0.3, stratify=enrollment_df[\"subject_id\"], random_state=RANDOM_STATE\n# )\n\n# template_embeddings = extract_embeddings(embedding_model, template_df, target_shape)\n# templates = build_subject_templates(template_embeddings, template_df[\"subject_id\"].values)\n\n# probe_embeddings = extract_embeddings(embedding_model, probe_df, target_shape)\n# genuine, impostor = generate_verification_scores(\n# probe_embeddings, probe_df[\"subject_id\"].values, templates\n# )\n\n# return compute_verification_metrics(genuine, impostor, label=\"Same-session baseline\")\n\n\n# # ----------------------------------------------------------------------\n# # MAIN\n# # ----------------------------------------------------------------------\n# if __name__ == \"__main__\":\n# print(\"=== Module 1: Cross-Session EEG Verification ===\\n\")\n\n# manifest = load_manifest()\n# target_shape = determine_target_shape(manifest)\n\n# enrollment, train_df, val_df, label_encoder = prepare_enrollment_data(manifest, target_shape)\n# n_classes = enrollment[\"subject_id\"].nunique()\n\n# embedding_model = train_embedding_extractor(train_df, val_df, target_shape, n_classes)\n\n# # --- Cross-session protocol ---\n# enrollment_embeddings = extract_embeddings(embedding_model, enrollment, target_shape)\n# templates = build_subject_templates(enrollment_embeddings, enrollment[\"subject_id\"].values)\n\n# calibration = manifest[manifest[\"split\"] == \"calibration\"].reset_index(drop=True)\n# calibration_embeddings = extract_embeddings(embedding_model, calibration, target_shape)\n\n# genuine_scores, impostor_scores = generate_verification_scores(\n# calibration_embeddings, calibration[\"subject_id\"].values, templates\n# )\n# cross_session_result = compute_verification_metrics(\n# genuine_scores, impostor_scores, label=\"Cross-session (Enrollment->Calibration)\"\n# )\n\n# # --- Same-session baseline ---\n# baseline_result = run_same_session_baseline(enrollment, embedding_model, target_shape)\n\n# # --- Save results ---\n# all_results = [cross_session_result, baseline_result]\n# with open(os.path.join(OUT_DIR, \"verification_results.json\"), \"w\") as f:\n# json.dump(all_results, f, indent=2)\n\n# plot_roc_curves(all_results, os.path.join(OUT_DIR, \"roc_comparison.png\"))\n\n# comparison_df = pd.DataFrame([\n# {\"protocol\": r[\"label\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n# for r in all_results\n# ])\n# comparison_df.to_csv(os.path.join(OUT_DIR, \"comparison_table.csv\"), index=False)\n# print(\"\\n=== Comparison: Same-session baseline vs Cross-session ===\")\n# print(comparison_df.to_string(index=False))\n\n# print(f\"\\n=== Module 1 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# \"\"\"\n# Module 1: Cross-Session EEG Biometric Verification Protocol\n# M3CV EEG Biometric Authentication Pipeline\n\n# Depends on Module 0's outputs: ./module0_outputs/spectrogram_manifest.csv\n# and the cached .npy spectrograms in ./module0_outputs/spectrogram_cache/.\n\n# Run steps in order via the __main__ block at the bottom.\n# Requires: pip install tensorflow scikit-learn numpy pandas matplotlib --break-system-packages\n# \"\"\"\n\n# import os\n# import json\n# import numpy as np\n# import pandas as pd\n# import tensorflow as tf\n# from sklearn.model_selection import train_test_split\n# from sklearn.metrics import roc_curve, roc_auc_score\n# from sklearn.preprocessing import LabelEncoder\n# import matplotlib.pyplot as plt\n\n# # ----------------------------------------------------------------------\n# # CONFIG\n# # ----------------------------------------------------------------------\n# MODULE0_DIR = \"./module0_outputs\"\n# MANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\n\n# OUT_DIR = \"./module1_outputs\"\n# os.makedirs(OUT_DIR, exist_ok=True)\n\n# EMBEDDING_DIM = 64\n# BATCH_SIZE = 32\n# EPOCHS = 25\n# LEARNING_RATE = 0.001\n# RANDOM_STATE = 42\n\n# np.random.seed(RANDOM_STATE)\n# tf.random.set_seed(RANDOM_STATE)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 1: Load manifest and cached spectrograms\n# # ----------------------------------------------------------------------\n# def load_manifest():\n# \"\"\"Load Module 0's manifest, which lists every cached spectrogram with its\n# split (enrollment/calibration), subject_id, session, and task.\"\"\"\n\n# manifest = pd.read_csv(MANIFEST_PATH)\n# print(f\"[load_manifest] Loaded {len(manifest)} entries\")\n# print(manifest[\"split\"].value_counts())\n# return manifest\n\n\n# def load_spectrogram_array(cache_path, target_shape=None):\n# \"\"\"Load one cached spectrogram .npy and optionally pad/crop to a target\n# shape, since mel spectrogram time-frames can vary by a few frames\n# depending on epoch length after resampling.\"\"\"\n\n# arr = np.load(cache_path)\n\n# if target_shape is not None:\n# h, w = arr.shape\n# th, tw = target_shape\n# if h != th or w != tw:\n# padded = np.zeros((th, tw), dtype=arr.dtype)\n# h_use, w_use = min(h, th), min(w, tw)\n# padded[:h_use, :w_use] = arr[:h_use, :w_use]\n# arr = padded\n\n# return arr\n\n\n# def determine_target_shape(manifest, n_check=20):\n# \"\"\"Inspect a sample of cached spectrograms to find the most common shape,\n# which becomes the fixed input shape for the CNN. Flags shape mismatch\n# so silent padding/cropping doesn't hide a preprocessing bug.\"\"\"\n\n# sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n# shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n# shape_counts = pd.Series(shapes).value_counts()\n\n# print(\"[determine_target_shape] Shape distribution across sample:\")\n# print(shape_counts)\n\n# target_shape = shape_counts.idxmax()\n# if len(shape_counts) > 1:\n# print(f\"[determine_target_shape] WARNING: inconsistent shapes found — \"\n# f\"using most common {target_shape}, others will be padded/cropped\")\n\n# return target_shape\n\n\n# # ----------------------------------------------------------------------\n# # STEP 2: Build train/template dataset from Enrollment\n# # ----------------------------------------------------------------------\n# def prepare_enrollment_data(manifest, target_shape, val_split=0.1):\n# \"\"\"Split Enrollment (Session 1) epochs for the cohort into a train/val set\n# for the embedding-extractor CNN. All Enrollment epochs (train+val) are\n# later used to build per-subject templates.\"\"\"\n\n# enrollment = manifest[manifest[\"split\"] == \"enrollment\"].reset_index(drop=True)\n\n# le = LabelEncoder()\n# enrollment[\"label\"] = le.fit_transform(enrollment[\"subject_id\"])\n\n# train_df, val_df = train_test_split(\n# enrollment, test_size=val_split, stratify=enrollment[\"label\"], random_state=RANDOM_STATE\n# )\n\n# print(f\"[prepare_enrollment_data] Train: {len(train_df)}, Val: {len(val_df)}, \"\n# f\"Subjects: {enrollment['subject_id'].nunique()}\")\n\n# return enrollment, train_df, val_df, le\n\n\n# def make_dataset(df, target_shape, batch_size=BATCH_SIZE, shuffle=True):\n# \"\"\"Build a tf.data.Dataset that loads cached spectrograms on the fly and\n# yields (spectrogram, label) pairs for training.\"\"\"\n\n# paths = df[\"cache_path\"].values\n# labels = df[\"label\"].values\n\n# def gen():\n# for p, y in zip(paths, labels):\n# arr = load_spectrogram_array(p, target_shape)\n# arr = arr[..., np.newaxis].astype(np.float32) # add channel dim\n# yield arr, y\n\n# ds = tf.data.Dataset.from_generator(\n# gen,\n# output_signature=(\n# tf.TensorSpec(shape=(*target_shape, 1), dtype=tf.float32),\n# tf.TensorSpec(shape=(), dtype=tf.int64),\n# ),\n# )\n\n# if shuffle:\n# ds = ds.shuffle(buffer_size=len(df), seed=RANDOM_STATE)\n\n# return ds.batch(batch_size).prefetch(tf.data.AUTOTUNE)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 3: CNN embedding-extractor architecture\n# # ----------------------------------------------------------------------\n# def build_model(input_shape, n_classes, embedding_dim=EMBEDDING_DIM):\n# \"\"\"CNN architecture following the base paper's Conv-Pool-Conv-Pool-Dropout\n# -Dense structure (Figure 5), with an explicit embedding_dim-wide dense\n# layer whose output is later reused as the identity embedding for\n# verification (not just the final softmax classification).\"\"\"\n\n# inputs = tf.keras.Input(shape=input_shape)\n\n# x = tf.keras.layers.Conv2D(32, 3, activation=\"relu\")(inputs)\n# x = tf.keras.layers.MaxPooling2D(2)(x)\n# x = tf.keras.layers.Conv2D(64, 3, activation=\"relu\")(x)\n# x = tf.keras.layers.MaxPooling2D(2)(x)\n# x = tf.keras.layers.Dropout(0.2)(x)\n# x = tf.keras.layers.Flatten()(x)\n# x = tf.keras.layers.Dense(128, activation=\"relu\")(x)\n\n# embedding = tf.keras.layers.Dense(embedding_dim, activation=\"relu\", name=\"embedding\")(x)\n# outputs = tf.keras.layers.Dense(n_classes, activation=\"softmax\", name=\"classifier\")(embedding)\n\n# model = tf.keras.Model(inputs=inputs, outputs=outputs)\n# embedding_model = tf.keras.Model(inputs=inputs, outputs=embedding)\n\n# model.compile(\n# optimizer=tf.keras.optimizers.Adam(learning_rate=LEARNING_RATE),\n# loss=\"sparse_categorical_crossentropy\",\n# metrics=[\"accuracy\"],\n# )\n\n# return model, embedding_model\n\n\n# def train_embedding_extractor(train_df, val_df, target_shape, n_classes):\n# \"\"\"Train the CNN as a 20-way subject classifier purely to learn useful\n# features — the trained embedding sub-model (penultimate dense layer) is\n# what Module 1 actually uses, not the softmax output.\"\"\"\n\n# train_ds = make_dataset(train_df, target_shape, shuffle=True)\n# val_ds = make_dataset(val_df, target_shape, shuffle=False)\n\n# model, embedding_model = build_model((*target_shape, 1), n_classes)\n\n# history = model.fit(\n# train_ds,\n# validation_data=val_ds,\n# epochs=EPOCHS,\n# verbose=1,\n# )\n\n# model.save(os.path.join(OUT_DIR, \"full_classifier.keras\"))\n# embedding_model.save(os.path.join(OUT_DIR, \"embedding_model.keras\"))\n\n# with open(os.path.join(OUT_DIR, \"training_history.json\"), \"w\") as f:\n# json.dump({k: [float(v) for v in vals] for k, vals in history.history.items()}, f, indent=2)\n\n# return embedding_model\n\n\n# # ----------------------------------------------------------------------\n# # STEP 4: Extract embeddings and build per-subject templates\n# # ----------------------------------------------------------------------\n# def extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n# \"\"\"Run every epoch in df through the trained embedding model, returning\n# an (n_epochs, embedding_dim) array aligned with df's row order.\"\"\"\n\n# ds = make_dataset(df.assign(label=0), target_shape, batch_size=batch_size, shuffle=False)\n# # label is a dummy placeholder here since we only need embeddings, not loss\n\n# embeddings = []\n# for batch_x, _ in ds:\n# emb = embedding_model.predict(batch_x, verbose=0)\n# embeddings.append(emb)\n\n# return np.concatenate(embeddings, axis=0)\n\n\n# def build_subject_templates(embeddings, subject_ids):\n# \"\"\"Average all of a subject's Enrollment embeddings into one template\n# vector per subject. Returns a dict {subject_id: template_vector}.\"\"\"\n\n# templates = {}\n# for sid in np.unique(subject_ids):\n# mask = subject_ids == sid\n# templates[sid] = embeddings[mask].mean(axis=0)\n\n# print(f\"[build_subject_templates] Built {len(templates)} subject templates\")\n# return templates\n\n\n# # ----------------------------------------------------------------------\n# # STEP 5: Genuine/impostor pair generation and scoring\n# # ----------------------------------------------------------------------\n# def cosine_similarity(a, b):\n# return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# def generate_verification_scores(probe_embeddings, probe_subject_ids, templates):\n# \"\"\"For every probe epoch, compute similarity against its OWN template\n# (genuine trial) and against every OTHER subject's template (impostor\n# trials). Returns two arrays: genuine_scores, impostor_scores.\"\"\"\n\n# subject_list = list(templates.keys())\n# genuine_scores, impostor_scores = [], []\n\n# for emb, true_sid in zip(probe_embeddings, probe_subject_ids):\n# for sid in subject_list:\n# score = cosine_similarity(emb, templates[sid])\n# if sid == true_sid:\n# genuine_scores.append(score)\n# else:\n# impostor_scores.append(score)\n\n# return np.array(genuine_scores), np.array(impostor_scores)\n\n\n# # ----------------------------------------------------------------------\n# # STEP 6: FAR / FRR / EER / ROC-AUC\n# # ----------------------------------------------------------------------\n# def compute_verification_metrics(genuine_scores, impostor_scores, label=\"\"):\n# \"\"\"Sweep thresholds over the combined score range to compute FAR/FRR at\n# each threshold, find EER (where FAR == FRR), and compute ROC-AUC using\n# genuine=1 / impostor=0 labels.\"\"\"\n\n# y_true = np.concatenate([np.ones_like(genuine_scores), np.zeros_like(impostor_scores)])\n# y_score = np.concatenate([genuine_scores, impostor_scores])\n\n# fpr, tpr, thresholds = roc_curve(y_true, y_score)\n# fnr = 1 - tpr # FRR\n# far = fpr\n\n# eer_idx = np.nanargmin(np.abs(far - fnr))\n# eer = (far[eer_idx] + fnr[eer_idx]) / 2\n# eer_threshold = thresholds[eer_idx]\n\n# auc = roc_auc_score(y_true, y_score)\n\n# print(f\"\\n[compute_verification_metrics] {label}\")\n# print(f\" EER: {eer:.4f} at threshold {eer_threshold:.4f}\")\n# print(f\" ROC-AUC: {auc:.4f}\")\n# print(f\" Genuine scores: mean={genuine_scores.mean():.4f}, std={genuine_scores.std():.4f}\")\n# print(f\" Impostor scores: mean={impostor_scores.mean():.4f}, std={impostor_scores.std():.4f}\")\n\n# return {\n# \"label\": label,\n# \"eer\": float(eer),\n# \"eer_threshold\": float(eer_threshold),\n# \"roc_auc\": float(auc),\n# \"fpr\": fpr.tolist(),\n# \"tpr\": tpr.tolist(),\n# \"genuine_mean\": float(genuine_scores.mean()),\n# \"impostor_mean\": float(impostor_scores.mean()),\n# }\n\n\n# def plot_roc_curves(results_list, save_path):\n# plt.figure(figsize=(7, 6))\n# for res in results_list:\n# plt.plot(res[\"fpr\"], res[\"tpr\"], label=f\"{res['label']} (AUC={res['roc_auc']:.3f})\")\n# plt.plot([0, 1], [0, 1], \"k--\", alpha=0.3)\n# plt.xlabel(\"False Positive Rate (FAR)\")\n# plt.ylabel(\"True Positive Rate (1 - FRR)\")\n# plt.title(\"ROC Curve: Cross-Session vs Same-Session Verification\")\n# plt.legend()\n# plt.savefig(save_path, dpi=150, bbox_inches=\"tight\")\n# plt.close()\n# print(f\"[plot_roc_curves] Saved to {save_path}\")\n\n\n# # ----------------------------------------------------------------------\n# # STEP 7: Same-session baseline (for comparison against cross-session result)\n# # ----------------------------------------------------------------------\n# def run_same_session_baseline(enrollment_df, embedding_model, target_shape):\n# \"\"\"Reproduce a same-session verification test: split Enrollment epochs\n# into a template half and a probe half (random split, same session),\n# exactly the kind of split the base paper used for its classification\n# task. This is the number your cross-session result should be compared\n# against.\"\"\"\n\n# template_df, probe_df = train_test_split(\n# enrollment_df, test_size=0.3, stratify=enrollment_df[\"subject_id\"], random_state=RANDOM_STATE\n# )\n\n# template_embeddings = extract_embeddings(embedding_model, template_df, target_shape)\n# templates = build_subject_templates(template_embeddings, template_df[\"subject_id\"].values)\n\n# probe_embeddings = extract_embeddings(embedding_model, probe_df, target_shape)\n# genuine, impostor = generate_verification_scores(\n# probe_embeddings, probe_df[\"subject_id\"].values, templates\n# )\n\n# return compute_verification_metrics(genuine, impostor, label=\"Same-session baseline\")\n\n\n# # ----------------------------------------------------------------------\n# # STEP 8: Task-matched diagnostic (disentangle task-blending from session drift)\n# # ----------------------------------------------------------------------\n# def run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape):\n# \"\"\"The main cross-session result blends embeddings across all 13 task\n# types into one template per subject, which conflates two different\n# effects: genuine session-to-session EEG drift, and task-to-task\n# mismatch (comparing a motor-movement probe against a template partly\n# built from resting-state epochs). This function isolates session drift\n# alone by building templates and probes from the SAME task each time,\n# so template task == probe task, only session differs.\"\"\"\n\n# common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n# print(f\"\\n[run_task_matched_diagnostic] Testing {len(common_tasks)} tasks present in both splits\")\n\n# task_results = []\n\n# for task in common_tasks:\n# enrollment_task = enrollment[enrollment[\"task\"] == task]\n# calibration_task = calibration[calibration[\"task\"] == task]\n\n# # Need every cohort subject represented in both splits for this task,\n# # otherwise genuine pairs can't be formed for the missing subject(s)\n# subjects_in_both = set(enrollment_task[\"subject_id\"]) & set(calibration_task[\"subject_id\"])\n# if len(subjects_in_both) < 2:\n# print(f\"[run_task_matched_diagnostic] Task {task}: skipped, \"\n# f\"only {len(subjects_in_both)} subject(s) present in both splits\")\n# continue\n\n# enrollment_task = enrollment_task[enrollment_task[\"subject_id\"].isin(subjects_in_both)]\n# calibration_task = calibration_task[calibration_task[\"subject_id\"].isin(subjects_in_both)]\n\n# template_embeddings = extract_embeddings(embedding_model, enrollment_task, target_shape)\n# templates = build_subject_templates(template_embeddings, enrollment_task[\"subject_id\"].values)\n\n# probe_embeddings = extract_embeddings(embedding_model, calibration_task, target_shape)\n# genuine, impostor = generate_verification_scores(\n# probe_embeddings, calibration_task[\"subject_id\"].values, templates\n# )\n\n# if len(genuine) == 0 or len(impostor) == 0:\n# print(f\"[run_task_matched_diagnostic] Task {task}: skipped, insufficient pairs\")\n# continue\n\n# result = compute_verification_metrics(genuine, impostor, label=f\"Task {task} (cross-session, matched)\")\n# result[\"task\"] = int(task)\n# result[\"n_subjects\"] = len(subjects_in_both)\n# result[\"n_enrollment_epochs\"] = len(enrollment_task)\n# result[\"n_calibration_epochs\"] = len(calibration_task)\n# task_results.append(result)\n\n# task_df = pd.DataFrame([\n# {\"task\": r[\"task\"], \"n_subjects\": r[\"n_subjects\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n# for r in task_results\n# ]).sort_values(\"EER\")\n\n# task_df.to_csv(os.path.join(OUT_DIR, \"task_matched_diagnostic.csv\"), index=False)\n\n# print(\"\\n=== Task-matched cross-session EER (session drift isolated) ===\")\n# print(task_df.to_string(index=False))\n# print(f\"\\nBlended-template cross-session EER (for comparison): see verification_results.json\")\n\n# return task_df, task_results\n# if __name__ == \"__main__\":\n# print(\"=== Module 1: Cross-Session EEG Verification ===\\n\")\n\n# manifest = load_manifest()\n# target_shape = determine_target_shape(manifest)\n\n# enrollment, train_df, val_df, label_encoder = prepare_enrollment_data(manifest, target_shape)\n# n_classes = enrollment[\"subject_id\"].nunique()\n\n# embedding_model = train_embedding_extractor(train_df, val_df, target_shape, n_classes)\n\n# # --- Cross-session protocol ---\n# enrollment_embeddings = extract_embeddings(embedding_model, enrollment, target_shape)\n# templates = build_subject_templates(enrollment_embeddings, enrollment[\"subject_id\"].values)\n\n# calibration = manifest[manifest[\"split\"] == \"calibration\"].reset_index(drop=True)\n# calibration_embeddings = extract_embeddings(embedding_model, calibration, target_shape)\n\n# genuine_scores, impostor_scores = generate_verification_scores(\n# calibration_embeddings, calibration[\"subject_id\"].values, templates\n# )\n# cross_session_result = compute_verification_metrics(\n# genuine_scores, impostor_scores, label=\"Cross-session (Enrollment->Calibration)\"\n# )\n\n# # --- Same-session baseline ---\n# baseline_result = run_same_session_baseline(enrollment, embedding_model, target_shape)\n\n# # --- Save results ---\n# all_results = [cross_session_result, baseline_result]\n# with open(os.path.join(OUT_DIR, \"verification_results.json\"), \"w\") as f:\n# json.dump(all_results, f, indent=2)\n\n# plot_roc_curves(all_results, os.path.join(OUT_DIR, \"roc_comparison.png\"))\n\n# comparison_df = pd.DataFrame([\n# {\"protocol\": r[\"label\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n# for r in all_results\n# ])\n# comparison_df.to_csv(os.path.join(OUT_DIR, \"comparison_table.csv\"), index=False)\n# print(\"\\n=== Comparison: Same-session baseline vs Cross-session ===\")\n# print(comparison_df.to_string(index=False))\n\n# # --- Task-matched diagnostic: isolate session drift from task-blending ---\n# task_df, task_results = run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape)\n\n# print(f\"\\n=== Module 1 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nDiagnostic: check whether cached spectrograms are actually distinguishable\nacross subjects, or degenerate (near-identical regardless of input).\n\nRun this BEFORE retraining anything in Module 1.\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nMANIFEST_PATH = \"./module0_outputs/spectrogram_manifest.csv\"\n\nmanifest = pd.read_csv(MANIFEST_PATH)\nenrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n\n# Pick one epoch from 5 different subjects, same task, to isolate subject effect\ntask = enrollment[\"task\"].mode()[0]\nsame_task = enrollment[enrollment[\"task\"] == task]\nsample_subjects = same_task[\"subject_id\"].unique()[:5]\n\nspecs = []\nlabels = []\nfor sid in sample_subjects:\n row = same_task[same_task[\"subject_id\"] == sid].iloc[0]\n arr = np.load(row[\"cache_path\"])\n specs.append(arr)\n labels.append(f\"subj {sid}\")\n\n# 1. Visual check\nfig, axes = plt.subplots(1, len(specs), figsize=(4 * len(specs), 4))\nfor ax, spec, label in zip(axes, specs, labels):\n im = ax.imshow(spec, aspect=\"auto\", origin=\"lower\")\n ax.set_title(label)\nplt.tight_layout()\nplt.savefig(\"./spectrogram_diagnostic.png\", dpi=150)\nprint(\"Saved visual comparison to ./spectrogram_diagnostic.png — \"\n \"inspect it: if all 5 panels look visually identical (same blocky \"\n \"pattern, no subject-specific structure), the spectrograms are degenerate.\")\n\n# 2. Statistical check: how much of each spectrogram is at the floor value (near -80dB)?\nfor spec, label in zip(specs, labels):\n floor = spec.min()\n frac_at_floor = np.mean(np.isclose(spec, floor, atol=1.0))\n print(f\"{label}: min={floor:.2f} dB, {frac_at_floor*100:.1f}% of bins near floor, \"\n f\"unique values={len(np.unique(np.round(spec, 1)))}\")\n\n# 3. Pairwise correlation between different subjects' spectrograms (same task)\nprint(\"\\nPairwise correlation between subjects (same task) — should be well below \"\n \"0.99 if subjects are actually distinguishable:\")\nflat = [s.flatten() for s in specs]\nfor i in range(len(flat)):\n for j in range(i + 1, len(flat)):\n corr = np.corrcoef(flat[i], flat[j])[0, 1]\n print(f\" {labels[i]} vs {labels[j]}: corr={corr:.4f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 1: Cross-Session EEG Biometric Verification Protocol\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on Module 0's outputs: ./module0_outputs/spectrogram_manifest.csv\nand the cached .npy spectrograms in ./module0_outputs/spectrogram_cache/.\n\nRun steps in order via the __main__ block at the bottom.\nRequires: pip install tensorflow scikit-learn numpy pandas matplotlib --break-system-packages\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_curve, roc_auc_score\nfrom sklearn.preprocessing import LabelEncoder\nimport matplotlib.pyplot as plt\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\n\nOUT_DIR = \"./module1_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nEMBEDDING_DIM = 64\nBATCH_SIZE = 32\nEPOCHS = 25\nLEARNING_RATE = 0.001\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\ntf.random.set_seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load manifest and cached spectrograms\n# ----------------------------------------------------------------------\ndef load_manifest():\n \"\"\"Load Module 0's manifest, which lists every cached spectrogram with its\n split (enrollment/calibration), subject_id, session, and task.\"\"\"\n\n manifest = pd.read_csv(MANIFEST_PATH)\n print(f\"[load_manifest] Loaded {len(manifest)} entries\")\n print(manifest[\"split\"].value_counts())\n return manifest\n\n\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Load one cached spectrogram .npy, optionally pad/crop to a target\n shape, and apply per-sample z-score normalization.\n\n NORMALIZATION IS NOT OPTIONAL HERE: cached spectrograms are raw dB values\n (roughly -80 to 0, never centered near zero). Feeding that directly into\n Conv2D is a classic setup for dying ReLU -- large-magnitude, all-negative\n input can saturate an entire network to zero output on the first forward\n pass depending on how randomly-initialized kernel weights happen to\n interact with that scale, with no gradient path to recover. That failure\n mode produces exactly the symptom of a collapsed embedding model (every\n sample maps to the same constant output, cosine similarity ~1.0 for every\n pair) and is sensitive to random init -- meaning it can appear to work on\n one run and silently collapse on the next with no code change. Per-sample\n standardization removes the failure mode instead of hoping to avoid it.\n \"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n \"\"\"Inspect a sample of cached spectrograms to find the most common shape,\n which becomes the fixed input shape for the CNN. Flags shape mismatch\n so silent padding/cropping doesn't hide a preprocessing bug.\"\"\"\n\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n\n print(\"[determine_target_shape] Shape distribution across sample:\")\n print(shape_counts)\n\n target_shape = shape_counts.idxmax()\n if len(shape_counts) > 1:\n print(f\"[determine_target_shape] WARNING: inconsistent shapes found — \"\n f\"using most common {target_shape}, others will be padded/cropped\")\n\n return target_shape\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Build train/template dataset from Enrollment\n# ----------------------------------------------------------------------\ndef prepare_enrollment_data(manifest, target_shape, val_split=0.1):\n \"\"\"Split Enrollment (Session 1) epochs for the cohort into a train/val set\n for the embedding-extractor CNN. All Enrollment epochs (train+val) are\n later used to build per-subject templates.\"\"\"\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"].reset_index(drop=True)\n\n le = LabelEncoder()\n enrollment[\"label\"] = le.fit_transform(enrollment[\"subject_id\"])\n\n train_df, val_df = train_test_split(\n enrollment, test_size=val_split, stratify=enrollment[\"label\"], random_state=RANDOM_STATE\n )\n\n print(f\"[prepare_enrollment_data] Train: {len(train_df)}, Val: {len(val_df)}, \"\n f\"Subjects: {enrollment['subject_id'].nunique()}\")\n\n return enrollment, train_df, val_df, le\n\n\ndef make_dataset(df, target_shape, batch_size=BATCH_SIZE, shuffle=True):\n \"\"\"Build a tf.data.Dataset that loads cached spectrograms on the fly and\n yields (spectrogram, label) pairs for training.\"\"\"\n\n paths = df[\"cache_path\"].values\n labels = df[\"label\"].values\n\n def gen():\n for p, y in zip(paths, labels):\n arr = load_spectrogram_array(p, target_shape)\n arr = arr[..., np.newaxis].astype(np.float32) # add channel dim\n yield arr, y\n\n ds = tf.data.Dataset.from_generator(\n gen,\n output_signature=(\n tf.TensorSpec(shape=(*target_shape, 1), dtype=tf.float32),\n tf.TensorSpec(shape=(), dtype=tf.int64),\n ),\n )\n\n if shuffle:\n ds = ds.shuffle(buffer_size=len(df), seed=RANDOM_STATE)\n\n return ds.batch(batch_size).prefetch(tf.data.AUTOTUNE)\n\n\n# ----------------------------------------------------------------------\n# STEP 3: CNN embedding-extractor architecture\n# ----------------------------------------------------------------------\ndef build_model(input_shape, n_classes, embedding_dim=EMBEDDING_DIM):\n \"\"\"CNN architecture following the base paper's Conv-Pool-Conv-Pool-Dropout\n -Dense structure (Figure 5), with an explicit embedding_dim-wide dense\n layer whose output is later reused as the identity embedding for\n verification (not just the final softmax classification).\"\"\"\n\n inputs = tf.keras.Input(shape=input_shape)\n\n x = tf.keras.layers.Conv2D(32, 3, activation=\"relu\")(inputs)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Conv2D(64, 3, activation=\"relu\")(x)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Dropout(0.2)(x)\n x = tf.keras.layers.Flatten()(x)\n x = tf.keras.layers.Dense(128, activation=\"relu\")(x)\n\n embedding = tf.keras.layers.Dense(embedding_dim, activation=\"relu\", name=\"embedding\")(x)\n outputs = tf.keras.layers.Dense(n_classes, activation=\"softmax\", name=\"classifier\")(embedding)\n\n model = tf.keras.Model(inputs=inputs, outputs=outputs)\n embedding_model = tf.keras.Model(inputs=inputs, outputs=embedding)\n\n model.compile(\n optimizer=tf.keras.optimizers.Adam(learning_rate=LEARNING_RATE),\n loss=\"sparse_categorical_crossentropy\",\n metrics=[\"accuracy\"],\n )\n\n return model, embedding_model\n\n\ndef train_embedding_extractor(train_df, val_df, target_shape, n_classes):\n \"\"\"Train the CNN as a 20-way subject classifier purely to learn useful\n features — the trained embedding sub-model (penultimate dense layer) is\n what Module 1 actually uses, not the softmax output.\"\"\"\n\n train_ds = make_dataset(train_df, target_shape, shuffle=True)\n val_ds = make_dataset(val_df, target_shape, shuffle=False)\n\n model, embedding_model = build_model((*target_shape, 1), n_classes)\n\n history = model.fit(\n train_ds,\n validation_data=val_ds,\n epochs=EPOCHS,\n verbose=1,\n )\n\n model.save(os.path.join(OUT_DIR, \"full_classifier.keras\"))\n embedding_model.save(os.path.join(OUT_DIR, \"embedding_model.keras\"))\n\n with open(os.path.join(OUT_DIR, \"training_history.json\"), \"w\") as f:\n json.dump({k: [float(v) for v in vals] for k, vals in history.history.items()}, f, indent=2)\n\n return embedding_model\n\n\ndef check_embedding_diversity(embedding_model, df, target_shape, n_check=30):\n \"\"\"Extract embeddings for a handful of random samples and check they're\n not all collapsed to the same vector. Run this immediately after\n training, BEFORE spending time on the full verification pipeline --\n catching a collapsed model here costs seconds; catching it after a full\n genuine/impostor scoring run costs the whole runtime for a useless\n result.\"\"\"\n\n sample = df.sample(min(n_check, len(df)), random_state=RANDOM_STATE)\n embeddings = extract_embeddings(embedding_model, sample, target_shape)\n\n pairwise_std = embeddings.std(axis=0).mean()\n print(f\"\\n[check_embedding_diversity] Mean per-dimension std across {len(sample)} \"\n f\"random samples: {pairwise_std:.6f}\")\n\n if pairwise_std < 1e-3:\n raise RuntimeError(\n \"[check_embedding_diversity] EMBEDDING MODEL HAS COLLAPSED: \"\n f\"mean std across dimensions is {pairwise_std:.6f} (near zero), meaning \"\n \"the model produces nearly the same output for every input regardless \"\n \"of content. Do not proceed to verification scoring -- retrain (rerun \"\n \"train_embedding_extractor), and if this recurs, check that input \"\n \"normalization is actually being applied in load_spectrogram_array.\"\n )\n\n print(\"[check_embedding_diversity] OK: embeddings show real variation, proceeding.\")\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Extract embeddings and build per-subject templates\n# ----------------------------------------------------------------------\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Run every epoch in df through the trained embedding model, returning\n an (n_epochs, embedding_dim) array aligned with df's row order.\"\"\"\n\n ds = make_dataset(df.assign(label=0), target_shape, batch_size=batch_size, shuffle=False)\n # label is a dummy placeholder here since we only need embeddings, not loss\n\n embeddings = []\n for batch_x, _ in ds:\n emb = embedding_model.predict(batch_x, verbose=0)\n embeddings.append(emb)\n\n return np.concatenate(embeddings, axis=0)\n\n\ndef build_subject_templates(embeddings, subject_ids):\n \"\"\"Average all of a subject's Enrollment embeddings into one template\n vector per subject. Returns a dict {subject_id: template_vector}.\"\"\"\n\n templates = {}\n for sid in np.unique(subject_ids):\n mask = subject_ids == sid\n templates[sid] = embeddings[mask].mean(axis=0)\n\n print(f\"[build_subject_templates] Built {len(templates)} subject templates\")\n return templates\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Genuine/impostor pair generation and scoring\n# ----------------------------------------------------------------------\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\ndef generate_verification_scores(probe_embeddings, probe_subject_ids, templates):\n \"\"\"For every probe epoch, compute similarity against its OWN template\n (genuine trial) and against every OTHER subject's template (impostor\n trials). Returns two arrays: genuine_scores, impostor_scores.\"\"\"\n\n subject_list = list(templates.keys())\n genuine_scores, impostor_scores = [], []\n\n for emb, true_sid in zip(probe_embeddings, probe_subject_ids):\n for sid in subject_list:\n score = cosine_similarity(emb, templates[sid])\n if sid == true_sid:\n genuine_scores.append(score)\n else:\n impostor_scores.append(score)\n\n return np.array(genuine_scores), np.array(impostor_scores)\n\n\n# ----------------------------------------------------------------------\n# STEP 6: FAR / FRR / EER / ROC-AUC\n# ----------------------------------------------------------------------\ndef compute_verification_metrics(genuine_scores, impostor_scores, label=\"\"):\n \"\"\"Sweep thresholds over the combined score range to compute FAR/FRR at\n each threshold, find EER (where FAR == FRR), and compute ROC-AUC using\n genuine=1 / impostor=0 labels.\"\"\"\n\n y_true = np.concatenate([np.ones_like(genuine_scores), np.zeros_like(impostor_scores)])\n y_score = np.concatenate([genuine_scores, impostor_scores])\n\n fpr, tpr, thresholds = roc_curve(y_true, y_score)\n fnr = 1 - tpr # FRR\n far = fpr\n\n eer_idx = np.nanargmin(np.abs(far - fnr))\n eer = (far[eer_idx] + fnr[eer_idx]) / 2\n eer_threshold = thresholds[eer_idx]\n\n auc = roc_auc_score(y_true, y_score)\n\n print(f\"\\n[compute_verification_metrics] {label}\")\n print(f\" EER: {eer:.4f} at threshold {eer_threshold:.4f}\")\n print(f\" ROC-AUC: {auc:.4f}\")\n print(f\" Genuine scores: mean={genuine_scores.mean():.4f}, std={genuine_scores.std():.4f}\")\n print(f\" Impostor scores: mean={impostor_scores.mean():.4f}, std={impostor_scores.std():.4f}\")\n\n return {\n \"label\": label,\n \"eer\": float(eer),\n \"eer_threshold\": float(eer_threshold),\n \"roc_auc\": float(auc),\n \"fpr\": fpr.tolist(),\n \"tpr\": tpr.tolist(),\n \"genuine_mean\": float(genuine_scores.mean()),\n \"impostor_mean\": float(impostor_scores.mean()),\n }\n\n\ndef plot_roc_curves(results_list, save_path):\n plt.figure(figsize=(7, 6))\n for res in results_list:\n plt.plot(res[\"fpr\"], res[\"tpr\"], label=f\"{res['label']} (AUC={res['roc_auc']:.3f})\")\n plt.plot([0, 1], [0, 1], \"k--\", alpha=0.3)\n plt.xlabel(\"False Positive Rate (FAR)\")\n plt.ylabel(\"True Positive Rate (1 - FRR)\")\n plt.title(\"ROC Curve: Cross-Session vs Same-Session Verification\")\n plt.legend()\n plt.savefig(save_path, dpi=150, bbox_inches=\"tight\")\n plt.close()\n print(f\"[plot_roc_curves] Saved to {save_path}\")\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Same-session baseline (for comparison against cross-session result)\n# ----------------------------------------------------------------------\ndef run_same_session_baseline(enrollment_df, embedding_model, target_shape):\n \"\"\"Reproduce a same-session verification test: split Enrollment epochs\n into a template half and a probe half (random split, same session),\n exactly the kind of split the base paper used for its classification\n task. This is the number your cross-session result should be compared\n against.\"\"\"\n\n template_df, probe_df = train_test_split(\n enrollment_df, test_size=0.3, stratify=enrollment_df[\"subject_id\"], random_state=RANDOM_STATE\n )\n\n template_embeddings = extract_embeddings(embedding_model, template_df, target_shape)\n templates = build_subject_templates(template_embeddings, template_df[\"subject_id\"].values)\n\n probe_embeddings = extract_embeddings(embedding_model, probe_df, target_shape)\n genuine, impostor = generate_verification_scores(\n probe_embeddings, probe_df[\"subject_id\"].values, templates\n )\n\n return compute_verification_metrics(genuine, impostor, label=\"Same-session baseline\")\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Task-matched diagnostic (disentangle task-blending from session drift)\n# ----------------------------------------------------------------------\ndef run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape):\n \"\"\"The main cross-session result blends embeddings across all 13 task\n types into one template per subject, which conflates two different\n effects: genuine session-to-session EEG drift, and task-to-task\n mismatch (comparing a motor-movement probe against a template partly\n built from resting-state epochs). This function isolates session drift\n alone by building templates and probes from the SAME task each time,\n so template task == probe task, only session differs.\"\"\"\n\n common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n print(f\"\\n[run_task_matched_diagnostic] Testing {len(common_tasks)} tasks present in both splits\")\n\n task_results = []\n\n for task in common_tasks:\n enrollment_task = enrollment[enrollment[\"task\"] == task]\n calibration_task = calibration[calibration[\"task\"] == task]\n\n # Need every cohort subject represented in both splits for this task,\n # otherwise genuine pairs can't be formed for the missing subject(s)\n subjects_in_both = set(enrollment_task[\"subject_id\"]) & set(calibration_task[\"subject_id\"])\n if len(subjects_in_both) < 2:\n print(f\"[run_task_matched_diagnostic] Task {task}: skipped, \"\n f\"only {len(subjects_in_both)} subject(s) present in both splits\")\n continue\n\n enrollment_task = enrollment_task[enrollment_task[\"subject_id\"].isin(subjects_in_both)]\n calibration_task = calibration_task[calibration_task[\"subject_id\"].isin(subjects_in_both)]\n\n template_embeddings = extract_embeddings(embedding_model, enrollment_task, target_shape)\n templates = build_subject_templates(template_embeddings, enrollment_task[\"subject_id\"].values)\n\n probe_embeddings = extract_embeddings(embedding_model, calibration_task, target_shape)\n genuine, impostor = generate_verification_scores(\n probe_embeddings, calibration_task[\"subject_id\"].values, templates\n )\n\n if len(genuine) == 0 or len(impostor) == 0:\n print(f\"[run_task_matched_diagnostic] Task {task}: skipped, insufficient pairs\")\n continue\n\n result = compute_verification_metrics(genuine, impostor, label=f\"Task {task} (cross-session, matched)\")\n result[\"task\"] = int(task)\n result[\"n_subjects\"] = len(subjects_in_both)\n result[\"n_enrollment_epochs\"] = len(enrollment_task)\n result[\"n_calibration_epochs\"] = len(calibration_task)\n task_results.append(result)\n\n task_df = pd.DataFrame([\n {\"task\": r[\"task\"], \"n_subjects\": r[\"n_subjects\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n for r in task_results\n ]).sort_values(\"EER\")\n\n task_df.to_csv(os.path.join(OUT_DIR, \"task_matched_diagnostic.csv\"), index=False)\n\n print(\"\\n=== Task-matched cross-session EER (session drift isolated) ===\")\n print(task_df.to_string(index=False))\n print(f\"\\nBlended-template cross-session EER (for comparison): see verification_results.json\")\n\n return task_df, task_results\nif __name__ == \"__main__\":\n print(\"=== Module 1: Cross-Session EEG Verification ===\\n\")\n\n manifest = load_manifest()\n target_shape = determine_target_shape(manifest)\n\n enrollment, train_df, val_df, label_encoder = prepare_enrollment_data(manifest, target_shape)\n n_classes = enrollment[\"subject_id\"].nunique()\n\n embedding_model = train_embedding_extractor(train_df, val_df, target_shape, n_classes)\n\n # Fail fast if training collapsed -- costs seconds here vs. the full\n # verification runtime if caught later\n check_embedding_diversity(embedding_model, val_df, target_shape)\n\n # --- Cross-session protocol ---\n enrollment_embeddings = extract_embeddings(embedding_model, enrollment, target_shape)\n templates = build_subject_templates(enrollment_embeddings, enrollment[\"subject_id\"].values)\n\n calibration = manifest[manifest[\"split\"] == \"calibration\"].reset_index(drop=True)\n calibration_embeddings = extract_embeddings(embedding_model, calibration, target_shape)\n\n genuine_scores, impostor_scores = generate_verification_scores(\n calibration_embeddings, calibration[\"subject_id\"].values, templates\n )\n cross_session_result = compute_verification_metrics(\n genuine_scores, impostor_scores, label=\"Cross-session (Enrollment->Calibration)\"\n )\n\n # --- Same-session baseline ---\n baseline_result = run_same_session_baseline(enrollment, embedding_model, target_shape)\n\n # --- Save results ---\n all_results = [cross_session_result, baseline_result]\n with open(os.path.join(OUT_DIR, \"verification_results.json\"), \"w\") as f:\n json.dump(all_results, f, indent=2)\n\n plot_roc_curves(all_results, os.path.join(OUT_DIR, \"roc_comparison.png\"))\n\n comparison_df = pd.DataFrame([\n {\"protocol\": r[\"label\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n for r in all_results\n ])\n comparison_df.to_csv(os.path.join(OUT_DIR, \"comparison_table.csv\"), index=False)\n print(\"\\n=== Comparison: Same-session baseline vs Cross-session ===\")\n print(comparison_df.to_string(index=False))\n\n # --- Task-matched diagnostic: isolate session drift from task-blending ---\n task_df, task_results = run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape)\n\n print(f\"\\n=== Module 1 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 2: Confidence-Guided Sequential Multi-Task Authentication\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n ./module1_outputs/task_matched_diagnostic.csv\n\nRun steps in order via the __main__ block at the bottom.\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nTASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\n\nOUT_DIR = \"./module2_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nN_STAGES = 3\nTARGET_FAR = 0.05 # preferred operating point per stage; falls back to EER\n # threshold if this FAR isn't achievable on that stage's data\nBATCH_SIZE = 32\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n task_eer = pd.read_csv(TASK_EER_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n print(f\"[load_dependencies] Task EER table: {len(task_eer)} tasks\")\n print(\"[load_dependencies] Task reliability (ascending EER):\")\n print(task_eer.sort_values(\"EER\")[[\"task\", \"EER\", \"ROC_AUC\"]].to_string(index=False))\n\n return manifest, task_eer, embedding_model\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Assign cascade stages by measured reliability\n# ----------------------------------------------------------------------\ndef assign_cascade_stages(task_eer, n_stages=N_STAGES):\n \"\"\"Split tasks into n_stages tiers by ascending EER (most session-stable\n first), using roughly equal-sized groups. Stage 1 = most reliable tasks,\n later stages = progressively less reliable. This replaces any assumed\n paradigm-based ordering (EO/EC -> motor -> stimulus) with an ordering\n justified directly by Module 1's measured cross-session performance.\"\"\"\n\n sorted_tasks = task_eer.sort_values(\"EER\").reset_index(drop=True)\n stage_assignment = pd.qcut(sorted_tasks.index, n_stages, labels=[f\"stage_{i+1}\" for i in range(n_stages)])\n sorted_tasks[\"stage\"] = stage_assignment\n\n stage_tasks = {\n stage: sorted_tasks[sorted_tasks[\"stage\"] == stage][\"task\"].tolist()\n for stage in sorted_tasks[\"stage\"].unique()\n }\n\n print(\"\\n[assign_cascade_stages] Stage assignment:\")\n for stage, tasks in stage_tasks.items():\n stage_eers = sorted_tasks[sorted_tasks[\"task\"].isin(tasks)][\"EER\"]\n print(f\" {stage}: tasks {tasks}, EER range [{stage_eers.min():.3f}, {stage_eers.max():.3f}]\")\n\n with open(os.path.join(OUT_DIR, \"stage_assignment.json\"), \"w\") as f:\n json.dump({k: [int(t) for t in v] for k, v in stage_tasks.items()}, f, indent=2)\n\n return stage_tasks\n\n\ndef drop_worst_tier(stage_tasks, task_eer):\n \"\"\"Drop the highest-mean-EER stage entirely from the cascade. The full\n n-stage assignment is still computed and saved by assign_cascade_stages\n (for transparency / reporting all 13 tasks' characterization), but only\n the more reliable stages are actually used in the authentication\n pipeline. This is a deliberate scope simplification: the resulting\n cascade covers a subset of tasks, not the full 13 -- state this\n explicitly in the thesis writeup rather than presenting it as full\n multi-task coverage.\"\"\"\n\n stage_mean_eer = {\n stage: task_eer[task_eer[\"task\"].isin(tasks)][\"EER\"].mean()\n for stage, tasks in stage_tasks.items()\n }\n worst_stage = max(stage_mean_eer, key=stage_mean_eer.get)\n\n kept_stage_tasks = {s: t for s, t in stage_tasks.items() if s != worst_stage}\n\n print(f\"\\n[drop_worst_tier] Dropping {worst_stage} (mean EER={stage_mean_eer[worst_stage]:.3f}, \"\n f\"tasks={stage_tasks[worst_stage]})\")\n print(f\"[drop_worst_tier] Retained stages: {list(kept_stage_tasks.keys())}\")\n\n with open(os.path.join(OUT_DIR, \"dropped_tier.json\"), \"w\") as f:\n json.dump({\n \"dropped_stage\": worst_stage,\n \"dropped_tasks\": [int(t) for t in stage_tasks[worst_stage]],\n \"dropped_mean_eer\": float(stage_mean_eer[worst_stage]),\n \"retained_stages\": {s: [int(t) for t in tasks] for s, tasks in kept_stage_tasks.items()},\n }, f, indent=2)\n\n return kept_stage_tasks\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Embedding extraction helpers (mirrors Module 1's normalization)\n# ----------------------------------------------------------------------\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Same per-sample z-score normalization as Module 1 -- must match, or\n the embedding model (trained on normalized input) will see out-of-\n distribution input and produce meaningless embeddings.\"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Batch-extract embeddings for every row in df, returning an\n (n_rows, embedding_dim) array aligned with df's row order.\"\"\"\n\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n\n embeddings = embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n return embeddings\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Build per-stage, per-subject embeddings (enrollment templates +\n# calibration probe representations)\n# ----------------------------------------------------------------------\ndef build_stage_representations(manifest, embedding_model, target_shape, stage_tasks, split):\n \"\"\"For each subject and each stage, average that subject's embeddings\n across all epochs whose task falls in that stage's task list, within the\n given split (enrollment or calibration). Returns\n {stage: {subject_id: mean_embedding}}. A subject missing all epochs for\n a stage in this split will simply be absent from that stage's dict --\n callers must handle missing stages as \"unavailable\", not as an error.\"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_reps = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n reps = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid]\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n reps[sid] = embeddings.mean(axis=0)\n\n stage_reps[stage] = reps\n print(f\"[build_stage_representations] {split}/{stage}: {len(reps)} subjects represented \"\n f\"(from {len(stage_df)} epochs across tasks {tasks})\")\n\n return stage_reps\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Pairwise stage scores (all genuine + impostor claims)\n# ----------------------------------------------------------------------\ndef compute_pairwise_stage_scores(enrollment_reps, calibration_reps, cohort, stage_order):\n \"\"\"For every (true_subject, claimed_subject) pair -- true_subject's\n Calibration probe representation vs claimed_subject's Enrollment\n template -- compute a score at every stage where both are available.\n Returns a long-format DataFrame: one row per (true_subject,\n claimed_subject, stage, score).\"\"\"\n\n rows = []\n\n for true_sid in cohort:\n for claimed_sid in cohort:\n for stage in stage_order:\n probe_emb = calibration_reps.get(stage, {}).get(true_sid)\n template_emb = enrollment_reps.get(stage, {}).get(claimed_sid)\n\n if probe_emb is None or template_emb is None:\n continue\n\n score = cosine_similarity(probe_emb, template_emb)\n rows.append({\n \"true_subject\": true_sid,\n \"claimed_subject\": claimed_sid,\n \"stage\": stage,\n \"is_genuine\": true_sid == claimed_sid,\n \"score\": score,\n })\n\n df = pd.DataFrame(rows)\n print(f\"\\n[compute_pairwise_stage_scores] {len(df)} stage-level scores computed \"\n f\"({df['is_genuine'].sum()} genuine, {(~df['is_genuine']).sum()} impostor)\")\n return df\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Compute running fused scores (mean of first k stages) per pair\n# ----------------------------------------------------------------------\ndef compute_k_stage_fused_scores(scores_df, stage_order):\n \"\"\"For every (true_subject, claimed_subject) pair, compute the running\n mean score across the first k stages, for every k = 1..len(stage_order).\n fused_k1 == the raw stage-1 score; fused_k2 == mean(stage-1, stage-2);\n etc. This is the quantity the cascade actually reasons about at each\n decision point -- not independent per-stage scores.\"\"\"\n\n pivot = scores_df.pivot_table(\n index=[\"true_subject\", \"claimed_subject\", \"is_genuine\"], columns=\"stage\", values=\"score\"\n ).reset_index()\n\n for k in range(1, len(stage_order) + 1):\n cols = stage_order[:k]\n pivot[f\"fused_k{k}\"] = pivot[cols].mean(axis=1, skipna=True)\n\n return pivot\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Calibrate three-zone thresholds per escalation level\n# ----------------------------------------------------------------------\ndef calibrate_k_thresholds(fused_df, stage_order, target_far=TARGET_FAR, reject_percentile=5):\n \"\"\"For each k, calibrate:\n - accept_threshold[k]: score above which we confidently accept, using\n the SAME target-FAR ROC calibration as before, but computed on the\n actual fused_k score distribution -- so FAR is controlled directly\n on what the cascade will actually see at that decision point,\n instead of compounding independently-calibrated per-stage bars.\n - reject_threshold[k]: a low percentile of the GENUINE fused_k score\n distribution. Scores below this are rare for genuine users, so\n rejecting here immediately (no further escalation) costs little\n FRR while stopping impostors from getting further independent\n chances -- this is the fix for the OR-gate compounding problem.\n \"\"\"\n\n accept_thresholds, reject_thresholds = {}, {}\n\n for k in range(1, len(stage_order) + 1):\n col = f\"fused_k{k}\"\n genuine = fused_df[fused_df[\"is_genuine\"]][col].dropna().values\n impostor = fused_df[~fused_df[\"is_genuine\"]][col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, thresh = roc_curve(y_true, y_score)\n\n idx_candidates = np.where(fpr <= target_far)[0]\n if len(idx_candidates) > 0:\n idx = idx_candidates[np.argmax(tpr[idx_candidates])]\n accept_thr = float(thresh[idx])\n achieved_far = fpr[idx]\n else:\n eer_idx = np.nanargmin(np.abs(fpr - (1 - tpr)))\n accept_thr = float(thresh[eer_idx])\n achieved_far = fpr[eer_idx]\n\n reject_thr = float(np.percentile(genuine, reject_percentile))\n\n accept_thresholds[k] = accept_thr\n reject_thresholds[k] = reject_thr\n\n print(f\"[calibrate_k_thresholds] k={k}: accept>={accept_thr:.4f} (achieved FAR={achieved_far:.3f}), \"\n f\"reject<={reject_thr:.4f} (bottom {reject_percentile}% of genuine)\")\n\n return accept_thresholds, reject_thresholds\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Run the fused-confidence cascade\n# ----------------------------------------------------------------------\ndef run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds):\n \"\"\"Three-zone decision at each escalation level k:\n - fused_k >= accept_threshold[k] -> ACCEPT, stop (confident match)\n - fused_k <= reject_threshold[k] -> REJECT, stop (confidently not a\n match, escalating further is pointless -- this is what prevents\n impostors from getting repeated independent chances)\n - otherwise -> escalate to k+1 (genuinely\n ambiguous, more evidence justified)\n If every stage is exhausted still in the ambiguous zone, decide by the\n final k's accept threshold (last resort, no separate looser fallback).\"\"\"\n\n results = []\n\n for _, row in fused_df.iterrows():\n decision, resolved_k = None, None\n\n for k in range(1, len(stage_order) + 1):\n score = row.get(f\"fused_k{k}\")\n if pd.isna(score) or k not in accept_thresholds:\n continue\n\n if score >= accept_thresholds[k]:\n decision, resolved_k = \"accept\", k\n break\n elif score <= reject_thresholds[k]:\n decision, resolved_k = \"reject\", k\n break\n # else: ambiguous, continue to next k\n\n if decision is None:\n last_k = max(accept_thresholds.keys())\n last_score = row.get(f\"fused_k{last_k}\")\n if pd.notna(last_score) and last_score >= accept_thresholds[last_k]:\n decision, resolved_k = \"accept\", last_k\n else:\n decision, resolved_k = \"reject\", last_k\n\n results.append({\n \"true_subject\": row[\"true_subject\"],\n \"claimed_subject\": row[\"claimed_subject\"],\n \"is_genuine\": row[\"is_genuine\"],\n \"decision\": decision,\n \"resolved_k\": resolved_k,\n })\n\n return pd.DataFrame(results)\n\n\n# ----------------------------------------------------------------------\n# STEP 9: Evaluate cascade vs no-cascade (stage-1-only) baseline\n# ----------------------------------------------------------------------\ndef evaluate_cascade(cascade_results):\n genuine = cascade_results[cascade_results[\"is_genuine\"]]\n impostor = cascade_results[~cascade_results[\"is_genuine\"]]\n\n far = (impostor[\"decision\"] == \"accept\").mean()\n frr = (genuine[\"decision\"] == \"reject\").mean()\n\n print(f\"\\n[evaluate_cascade] Cascade FAR: {far:.4f}, FRR: {frr:.4f}\")\n\n resolution_dist = cascade_results.groupby([\"decision\"])[\"resolved_k\"].value_counts(normalize=False)\n print(\"[evaluate_cascade] Resolution counts by decision and escalation level (k):\")\n print(resolution_dist.to_string())\n\n genuine_accept_dist = genuine[genuine[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n print(\"[evaluate_cascade] Genuine users authenticated, by escalation level (k) reached:\")\n print(genuine_accept_dist.to_string())\n\n impostor_accept_dist = impostor[impostor[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n if len(impostor_accept_dist) > 0:\n print(\"[evaluate_cascade] Impostors wrongly accepted, by escalation level (k):\")\n print(impostor_accept_dist.to_string())\n\n return {\n \"far\": float(far),\n \"frr\": float(frr),\n \"genuine_resolution_by_k\": genuine_accept_dist.to_dict(),\n \"impostor_false_accept_by_k\": impostor_accept_dist.to_dict(),\n }\n\n\ndef evaluate_stage1_only_baseline(fused_df, accept_thresholds):\n \"\"\"No-cascade baseline: decide using ONLY the k=1 (first stage) fused\n score and its calibrated accept threshold, no escalation at all.\"\"\"\n\n genuine = fused_df[fused_df[\"is_genuine\"]]\n impostor = fused_df[~fused_df[\"is_genuine\"]]\n\n threshold = accept_thresholds[1]\n far = (impostor[\"fused_k1\"] >= threshold).mean()\n frr = (genuine[\"fused_k1\"] < threshold).mean()\n\n print(f\"\\n[evaluate_stage1_only_baseline] Stage-1-only FAR: {far:.4f}, FRR: {frr:.4f}\")\n return {\"far\": float(far), \"frr\": float(frr)}\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 2: Confidence-Guided Sequential Authentication ===\\n\")\n\n manifest, task_eer, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n stage_tasks_full = assign_cascade_stages(task_eer, n_stages=N_STAGES)\n stage_tasks = drop_worst_tier(stage_tasks_full, task_eer)\n stage_order = sorted(stage_tasks.keys()) # e.g. stage_1, stage_2 (worst tier dropped)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"\\n[main] Cohort: {len(cohort)} subjects\")\n\n enrollment_reps = build_stage_representations(manifest, embedding_model, target_shape, stage_tasks, \"enrollment\")\n calibration_reps = build_stage_representations(manifest, embedding_model, target_shape, stage_tasks, \"calibration\")\n\n scores_df = compute_pairwise_stage_scores(enrollment_reps, calibration_reps, cohort, stage_order)\n scores_df.to_csv(os.path.join(OUT_DIR, \"pairwise_stage_scores.csv\"), index=False)\n\n fused_df = compute_k_stage_fused_scores(scores_df, stage_order)\n fused_df.to_csv(os.path.join(OUT_DIR, \"fused_scores.csv\"), index=False)\n\n accept_thresholds, reject_thresholds = calibrate_k_thresholds(fused_df, stage_order)\n\n cascade_results = run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds)\n cascade_results.to_csv(os.path.join(OUT_DIR, \"cascade_decisions.csv\"), index=False)\n\n cascade_metrics = evaluate_cascade(cascade_results)\n baseline_metrics = evaluate_stage1_only_baseline(fused_df, accept_thresholds)\n\n comparison = pd.DataFrame([\n {\"protocol\": \"Stage-1-only (no cascade)\", **baseline_metrics},\n {\"protocol\": \"Full cascade (worst tier dropped)\",\n **{k: v for k, v in cascade_metrics.items() if k in (\"far\", \"frr\")}},\n ])\n comparison.to_csv(os.path.join(OUT_DIR, \"cascade_vs_baseline.csv\"), index=False)\n\n print(\"\\n=== Comparison: No-cascade baseline vs Full cascade ===\")\n print(comparison.to_string(index=False))\n\n with open(os.path.join(OUT_DIR, \"module2_summary.json\"), \"w\") as f:\n json.dump({\n \"stage_order\": stage_order,\n \"accept_thresholds\": accept_thresholds,\n \"reject_thresholds\": reject_thresholds,\n \"cascade_metrics\": cascade_metrics,\n \"baseline_metrics\": baseline_metrics,\n }, f, indent=2)\n\n print(f\"\\n=== Module 2 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 1: Cross-Session EEG Biometric Verification Protocol\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on Module 0's outputs: ./module0_outputs/spectrogram_manifest.csv\nand the cached .npy spectrograms in ./module0_outputs/spectrogram_cache/.\n\nRun steps in order via the __main__ block at the bottom.\nRequires: pip install tensorflow scikit-learn numpy pandas matplotlib --break-system-packages\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_curve, roc_auc_score, classification_report\nfrom sklearn.preprocessing import LabelEncoder\nimport matplotlib.pyplot as plt\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\n\nOUT_DIR = \"./module1_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nEMBEDDING_DIM = 64\nBATCH_SIZE = 32\nEPOCHS = 25\nLEARNING_RATE = 0.001\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\ntf.random.set_seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load manifest and cached spectrograms\n# ----------------------------------------------------------------------\ndef load_manifest():\n \"\"\"Load Module 0's manifest, which lists every cached spectrogram with its\n split (enrollment/calibration), subject_id, session, and task.\"\"\"\n\n manifest = pd.read_csv(MANIFEST_PATH)\n print(f\"[load_manifest] Loaded {len(manifest)} entries\")\n print(manifest[\"split\"].value_counts())\n return manifest\n\n\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Load one cached spectrogram .npy, optionally pad/crop to a target\n shape, and apply per-sample z-score normalization.\n\n NORMALIZATION IS NOT OPTIONAL HERE: cached spectrograms are raw dB values\n (roughly -80 to 0, never centered near zero). Feeding that directly into\n Conv2D is a classic setup for dying ReLU -- large-magnitude, all-negative\n input can saturate an entire network to zero output on the first forward\n pass depending on how randomly-initialized kernel weights happen to\n interact with that scale, with no gradient path to recover. That failure\n mode produces exactly the symptom of a collapsed embedding model (every\n sample maps to the same constant output, cosine similarity ~1.0 for every\n pair) and is sensitive to random init -- meaning it can appear to work on\n one run and silently collapse on the next with no code change. Per-sample\n standardization removes the failure mode instead of hoping to avoid it.\n \"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n \"\"\"Inspect a sample of cached spectrograms to find the most common shape,\n which becomes the fixed input shape for the CNN. Flags shape mismatch\n so silent padding/cropping doesn't hide a preprocessing bug.\"\"\"\n\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n\n print(\"[determine_target_shape] Shape distribution across sample:\")\n print(shape_counts)\n\n target_shape = shape_counts.idxmax()\n if len(shape_counts) > 1:\n print(f\"[determine_target_shape] WARNING: inconsistent shapes found — \"\n f\"using most common {target_shape}, others will be padded/cropped\")\n\n return target_shape\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Build train/template dataset from Enrollment\n# ----------------------------------------------------------------------\ndef prepare_enrollment_data(manifest, target_shape, val_split=0.1):\n \"\"\"Split Enrollment (Session 1) epochs for the cohort into a train/val set\n for the embedding-extractor CNN. All Enrollment epochs (train+val) are\n later used to build per-subject templates.\"\"\"\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"].reset_index(drop=True)\n\n le = LabelEncoder()\n enrollment[\"label\"] = le.fit_transform(enrollment[\"subject_id\"])\n\n train_df, val_df = train_test_split(\n enrollment, test_size=val_split, stratify=enrollment[\"label\"], random_state=RANDOM_STATE\n )\n\n print(f\"[prepare_enrollment_data] Train: {len(train_df)}, Val: {len(val_df)}, \"\n f\"Subjects: {enrollment['subject_id'].nunique()}\")\n\n return enrollment, train_df, val_df, le\n\n\ndef make_dataset(df, target_shape, batch_size=BATCH_SIZE, shuffle=True):\n \"\"\"Build a tf.data.Dataset that loads cached spectrograms on the fly and\n yields (spectrogram, label) pairs for training.\"\"\"\n\n paths = df[\"cache_path\"].values\n labels = df[\"label\"].values\n\n def gen():\n for p, y in zip(paths, labels):\n arr = load_spectrogram_array(p, target_shape)\n arr = arr[..., np.newaxis].astype(np.float32) # add channel dim\n yield arr, y\n\n ds = tf.data.Dataset.from_generator(\n gen,\n output_signature=(\n tf.TensorSpec(shape=(*target_shape, 1), dtype=tf.float32),\n tf.TensorSpec(shape=(), dtype=tf.int64),\n ),\n )\n\n if shuffle:\n ds = ds.shuffle(buffer_size=len(df), seed=RANDOM_STATE)\n\n return ds.batch(batch_size).prefetch(tf.data.AUTOTUNE)\n\n\n# ----------------------------------------------------------------------\n# STEP 3: CNN embedding-extractor architecture\n# ----------------------------------------------------------------------\ndef build_model(input_shape, n_classes, embedding_dim=EMBEDDING_DIM):\n \"\"\"CNN architecture following the base paper's Conv-Pool-Conv-Pool-Dropout\n -Dense structure (Figure 5), with an explicit embedding_dim-wide dense\n layer whose output is later reused as the identity embedding for\n verification (not just the final softmax classification).\"\"\"\n\n inputs = tf.keras.Input(shape=input_shape)\n\n x = tf.keras.layers.Conv2D(32, 3, activation=\"relu\")(inputs)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Conv2D(64, 3, activation=\"relu\")(x)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Dropout(0.2)(x)\n x = tf.keras.layers.Flatten()(x)\n x = tf.keras.layers.Dense(128, activation=\"relu\")(x)\n\n embedding = tf.keras.layers.Dense(embedding_dim, activation=\"relu\", name=\"embedding\")(x)\n outputs = tf.keras.layers.Dense(n_classes, activation=\"softmax\", name=\"classifier\")(embedding)\n\n model = tf.keras.Model(inputs=inputs, outputs=outputs)\n embedding_model = tf.keras.Model(inputs=inputs, outputs=embedding)\n\n model.compile(\n optimizer=tf.keras.optimizers.Adam(learning_rate=LEARNING_RATE),\n loss=\"sparse_categorical_crossentropy\",\n metrics=[\"accuracy\"],\n )\n\n return model, embedding_model\n\n\ndef train_embedding_extractor(train_df, val_df, target_shape, n_classes):\n \"\"\"Train the CNN as a 20-way subject classifier purely to learn useful\n features — the trained embedding sub-model (penultimate dense layer) is\n what Module 1 actually uses for verification, but the full classifier\n is also returned so precision/recall/F1 can be evaluated on the\n classification proxy task.\"\"\"\n\n train_ds = make_dataset(train_df, target_shape, shuffle=True)\n val_ds = make_dataset(val_df, target_shape, shuffle=False)\n\n model, embedding_model = build_model((*target_shape, 1), n_classes)\n\n history = model.fit(\n train_ds,\n validation_data=val_ds,\n epochs=EPOCHS,\n verbose=1,\n )\n\n model.save(os.path.join(OUT_DIR, \"full_classifier.keras\"))\n embedding_model.save(os.path.join(OUT_DIR, \"embedding_model.keras\"))\n\n with open(os.path.join(OUT_DIR, \"training_history.json\"), \"w\") as f:\n json.dump({k: [float(v) for v in vals] for k, vals in history.history.items()}, f, indent=2)\n\n return model, embedding_model\n\n\ndef check_embedding_diversity(embedding_model, df, target_shape, n_check=30):\n \"\"\"Extract embeddings for a handful of random samples and check they're\n not all collapsed to the same vector. Run this immediately after\n training, BEFORE spending time on the full verification pipeline --\n catching a collapsed model here costs seconds; catching it after a full\n genuine/impostor scoring run costs the whole runtime for a useless\n result.\"\"\"\n\n sample = df.sample(min(n_check, len(df)), random_state=RANDOM_STATE)\n embeddings = extract_embeddings(embedding_model, sample, target_shape)\n\n pairwise_std = embeddings.std(axis=0).mean()\n print(f\"\\n[check_embedding_diversity] Mean per-dimension std across {len(sample)} \"\n f\"random samples: {pairwise_std:.6f}\")\n\n if pairwise_std < 1e-3:\n raise RuntimeError(\n \"[check_embedding_diversity] EMBEDDING MODEL HAS COLLAPSED: \"\n f\"mean std across dimensions is {pairwise_std:.6f} (near zero), meaning \"\n \"the model produces nearly the same output for every input regardless \"\n \"of content. Do not proceed to verification scoring -- retrain (rerun \"\n \"train_embedding_extractor), and if this recurs, check that input \"\n \"normalization is actually being applied in load_spectrogram_array.\"\n )\n\n print(\"[check_embedding_diversity] OK: embeddings show real variation, proceeding.\")\n\n\ndef evaluate_classification_metrics(model, val_df, target_shape, label_encoder):\n \"\"\"Compute precision/recall/F1 on the 20-way subject classification proxy\n task (validation split), via sklearn.classification_report -- the same\n metrics reported in the base paper's Table 3, for direct comparison.\n\n IMPORTANT: this describes the closed-set classification task the CNN was\n trained on (can it tell 20 known subjects apart), NOT your actual\n verification protocol (FAR/FRR/EER, does it correctly accept/reject\n open-set claims). Do not present these numbers as verification results\n -- they answer a different, easier question, and reporting them as your\n headline result would put your framing back in the base paper's\n identification-only paradigm, which is the gap this thesis argues is\n insufficient in the first place.\n \"\"\"\n\n val_ds = make_dataset(val_df, target_shape, shuffle=False)\n\n y_true, y_pred = [], []\n for batch_x, batch_y in val_ds:\n preds = model.predict(batch_x, verbose=0)\n y_pred.extend(np.argmax(preds, axis=1))\n y_true.extend(batch_y.numpy())\n\n target_names = [str(c) for c in label_encoder.classes_]\n report = classification_report(y_true, y_pred, target_names=target_names,\n output_dict=True, zero_division=0)\n report_df = pd.DataFrame(report).transpose()\n report_df.to_csv(os.path.join(OUT_DIR, \"classification_report.csv\"))\n\n print(\"\\n[evaluate_classification_metrics] Precision/Recall/F1 -- \"\n \"20-way subject classification proxy task, validation split \"\n \"(NOT the verification protocol -- see docstring):\")\n print(report_df.to_string())\n\n return report\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Extract embeddings and build per-subject templates\n# ----------------------------------------------------------------------\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Run every epoch in df through the trained embedding model, returning\n an (n_epochs, embedding_dim) array aligned with df's row order.\"\"\"\n\n ds = make_dataset(df.assign(label=0), target_shape, batch_size=batch_size, shuffle=False)\n # label is a dummy placeholder here since we only need embeddings, not loss\n\n embeddings = []\n for batch_x, _ in ds:\n emb = embedding_model.predict(batch_x, verbose=0)\n embeddings.append(emb)\n\n return np.concatenate(embeddings, axis=0)\n\n\ndef build_subject_templates(embeddings, subject_ids):\n \"\"\"Average all of a subject's Enrollment embeddings into one template\n vector per subject. Returns a dict {subject_id: template_vector}.\"\"\"\n\n templates = {}\n for sid in np.unique(subject_ids):\n mask = subject_ids == sid\n templates[sid] = embeddings[mask].mean(axis=0)\n\n print(f\"[build_subject_templates] Built {len(templates)} subject templates\")\n return templates\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Genuine/impostor pair generation and scoring\n# ----------------------------------------------------------------------\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\ndef generate_verification_scores(probe_embeddings, probe_subject_ids, templates):\n \"\"\"For every probe epoch, compute similarity against its OWN template\n (genuine trial) and against every OTHER subject's template (impostor\n trials). Returns two arrays: genuine_scores, impostor_scores.\"\"\"\n\n subject_list = list(templates.keys())\n genuine_scores, impostor_scores = [], []\n\n for emb, true_sid in zip(probe_embeddings, probe_subject_ids):\n for sid in subject_list:\n score = cosine_similarity(emb, templates[sid])\n if sid == true_sid:\n genuine_scores.append(score)\n else:\n impostor_scores.append(score)\n\n return np.array(genuine_scores), np.array(impostor_scores)\n\n\n# ----------------------------------------------------------------------\n# STEP 6: FAR / FRR / EER / ROC-AUC\n# ----------------------------------------------------------------------\ndef compute_verification_metrics(genuine_scores, impostor_scores, label=\"\"):\n \"\"\"Sweep thresholds over the combined score range to compute FAR/FRR at\n each threshold, find EER (where FAR == FRR), and compute ROC-AUC using\n genuine=1 / impostor=0 labels.\"\"\"\n\n y_true = np.concatenate([np.ones_like(genuine_scores), np.zeros_like(impostor_scores)])\n y_score = np.concatenate([genuine_scores, impostor_scores])\n\n fpr, tpr, thresholds = roc_curve(y_true, y_score)\n fnr = 1 - tpr # FRR\n far = fpr\n\n eer_idx = np.nanargmin(np.abs(far - fnr))\n eer = (far[eer_idx] + fnr[eer_idx]) / 2\n eer_threshold = thresholds[eer_idx]\n\n auc = roc_auc_score(y_true, y_score)\n\n print(f\"\\n[compute_verification_metrics] {label}\")\n print(f\" EER: {eer:.4f} at threshold {eer_threshold:.4f}\")\n print(f\" ROC-AUC: {auc:.4f}\")\n print(f\" Genuine scores: mean={genuine_scores.mean():.4f}, std={genuine_scores.std():.4f}\")\n print(f\" Impostor scores: mean={impostor_scores.mean():.4f}, std={impostor_scores.std():.4f}\")\n\n return {\n \"label\": label,\n \"eer\": float(eer),\n \"eer_threshold\": float(eer_threshold),\n \"roc_auc\": float(auc),\n \"fpr\": fpr.tolist(),\n \"tpr\": tpr.tolist(),\n \"genuine_mean\": float(genuine_scores.mean()),\n \"impostor_mean\": float(impostor_scores.mean()),\n }\n\n\ndef plot_roc_curves(results_list, save_path):\n plt.figure(figsize=(7, 6))\n for res in results_list:\n plt.plot(res[\"fpr\"], res[\"tpr\"], label=f\"{res['label']} (AUC={res['roc_auc']:.3f})\")\n plt.plot([0, 1], [0, 1], \"k--\", alpha=0.3)\n plt.xlabel(\"False Positive Rate (FAR)\")\n plt.ylabel(\"True Positive Rate (1 - FRR)\")\n plt.title(\"ROC Curve: Cross-Session vs Same-Session Verification\")\n plt.legend()\n plt.savefig(save_path, dpi=150, bbox_inches=\"tight\")\n plt.close()\n print(f\"[plot_roc_curves] Saved to {save_path}\")\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Same-session baseline (for comparison against cross-session result)\n# ----------------------------------------------------------------------\ndef run_same_session_baseline(enrollment_df, embedding_model, target_shape):\n \"\"\"Reproduce a same-session verification test: split Enrollment epochs\n into a template half and a probe half (random split, same session),\n exactly the kind of split the base paper used for its classification\n task. This is the number your cross-session result should be compared\n against.\"\"\"\n\n template_df, probe_df = train_test_split(\n enrollment_df, test_size=0.3, stratify=enrollment_df[\"subject_id\"], random_state=RANDOM_STATE\n )\n\n template_embeddings = extract_embeddings(embedding_model, template_df, target_shape)\n templates = build_subject_templates(template_embeddings, template_df[\"subject_id\"].values)\n\n probe_embeddings = extract_embeddings(embedding_model, probe_df, target_shape)\n genuine, impostor = generate_verification_scores(\n probe_embeddings, probe_df[\"subject_id\"].values, templates\n )\n\n return compute_verification_metrics(genuine, impostor, label=\"Same-session baseline\")\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Task-matched diagnostic (disentangle task-blending from session drift)\n# ----------------------------------------------------------------------\ndef run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape):\n \"\"\"The main cross-session result blends embeddings across all 13 task\n types into one template per subject, which conflates two different\n effects: genuine session-to-session EEG drift, and task-to-task\n mismatch (comparing a motor-movement probe against a template partly\n built from resting-state epochs). This function isolates session drift\n alone by building templates and probes from the SAME task each time,\n so template task == probe task, only session differs.\"\"\"\n\n common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n print(f\"\\n[run_task_matched_diagnostic] Testing {len(common_tasks)} tasks present in both splits\")\n\n task_results = []\n\n for task in common_tasks:\n enrollment_task = enrollment[enrollment[\"task\"] == task]\n calibration_task = calibration[calibration[\"task\"] == task]\n\n # Need every cohort subject represented in both splits for this task,\n # otherwise genuine pairs can't be formed for the missing subject(s)\n subjects_in_both = set(enrollment_task[\"subject_id\"]) & set(calibration_task[\"subject_id\"])\n if len(subjects_in_both) < 2:\n print(f\"[run_task_matched_diagnostic] Task {task}: skipped, \"\n f\"only {len(subjects_in_both)} subject(s) present in both splits\")\n continue\n\n enrollment_task = enrollment_task[enrollment_task[\"subject_id\"].isin(subjects_in_both)]\n calibration_task = calibration_task[calibration_task[\"subject_id\"].isin(subjects_in_both)]\n\n template_embeddings = extract_embeddings(embedding_model, enrollment_task, target_shape)\n templates = build_subject_templates(template_embeddings, enrollment_task[\"subject_id\"].values)\n\n probe_embeddings = extract_embeddings(embedding_model, calibration_task, target_shape)\n genuine, impostor = generate_verification_scores(\n probe_embeddings, calibration_task[\"subject_id\"].values, templates\n )\n\n if len(genuine) == 0 or len(impostor) == 0:\n print(f\"[run_task_matched_diagnostic] Task {task}: skipped, insufficient pairs\")\n continue\n\n result = compute_verification_metrics(genuine, impostor, label=f\"Task {task} (cross-session, matched)\")\n result[\"task\"] = int(task)\n result[\"n_subjects\"] = len(subjects_in_both)\n result[\"n_enrollment_epochs\"] = len(enrollment_task)\n result[\"n_calibration_epochs\"] = len(calibration_task)\n task_results.append(result)\n\n task_df = pd.DataFrame([\n {\"task\": r[\"task\"], \"n_subjects\": r[\"n_subjects\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n for r in task_results\n ]).sort_values(\"EER\")\n\n task_df.to_csv(os.path.join(OUT_DIR, \"task_matched_diagnostic.csv\"), index=False)\n\n print(\"\\n=== Task-matched cross-session EER (session drift isolated) ===\")\n print(task_df.to_string(index=False))\n print(f\"\\nBlended-template cross-session EER (for comparison): see verification_results.json\")\n\n return task_df, task_results\nif __name__ == \"__main__\":\n print(\"=== Module 1: Cross-Session EEG Verification ===\\n\")\n\n manifest = load_manifest()\n target_shape = determine_target_shape(manifest)\n\n enrollment, train_df, val_df, label_encoder = prepare_enrollment_data(manifest, target_shape)\n n_classes = enrollment[\"subject_id\"].nunique()\n\n full_model, embedding_model = train_embedding_extractor(train_df, val_df, target_shape, n_classes)\n\n # Fail fast if training collapsed -- costs seconds here vs. the full\n # verification runtime if caught later\n check_embedding_diversity(embedding_model, val_df, target_shape)\n\n # Precision/recall/F1 on the classification proxy task -- see function\n # docstring for why these are NOT verification metrics\n evaluate_classification_metrics(full_model, val_df, target_shape, label_encoder)\n\n # --- Cross-session protocol ---\n enrollment_embeddings = extract_embeddings(embedding_model, enrollment, target_shape)\n templates = build_subject_templates(enrollment_embeddings, enrollment[\"subject_id\"].values)\n\n calibration = manifest[manifest[\"split\"] == \"calibration\"].reset_index(drop=True)\n calibration_embeddings = extract_embeddings(embedding_model, calibration, target_shape)\n\n genuine_scores, impostor_scores = generate_verification_scores(\n calibration_embeddings, calibration[\"subject_id\"].values, templates\n )\n cross_session_result = compute_verification_metrics(\n genuine_scores, impostor_scores, label=\"Cross-session (Enrollment->Calibration)\"\n )\n\n # --- Same-session baseline ---\n baseline_result = run_same_session_baseline(enrollment, embedding_model, target_shape)\n\n # --- Save results ---\n all_results = [cross_session_result, baseline_result]\n with open(os.path.join(OUT_DIR, \"verification_results.json\"), \"w\") as f:\n json.dump(all_results, f, indent=2)\n\n plot_roc_curves(all_results, os.path.join(OUT_DIR, \"roc_comparison.png\"))\n\n comparison_df = pd.DataFrame([\n {\"protocol\": r[\"label\"], \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]}\n for r in all_results\n ])\n comparison_df.to_csv(os.path.join(OUT_DIR, \"comparison_table.csv\"), index=False)\n print(\"\\n=== Comparison: Same-session baseline vs Cross-session ===\")\n print(comparison_df.to_string(index=False))\n\n # --- Task-matched diagnostic: isolate session drift from task-blending ---\n task_df, task_results = run_task_matched_diagnostic(enrollment, calibration, embedding_model, target_shape)\n\n print(f\"\\n=== Module 1 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 2: Confidence-Guided Sequential Multi-Task Authentication\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n ./module1_outputs/task_matched_diagnostic.csv\n\nRun steps in order via the __main__ block at the bottom.\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nTASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\n\nOUT_DIR = \"./module2_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nN_STAGES = 3\nTARGET_FAR = 0.05 # preferred operating point per stage; falls back to EER\n # threshold if this FAR isn't achievable on that stage's data\nEPOCHS_PER_TRIAL = 10 # calibration epochs are chunked into groups of this size;\n # each chunk becomes one independent probe \"trial\" instead\n # of averaging a subject's entire session into one probe.\n # ~150-200 calibration epochs/subject/stage -> ~15-20\n # trials/subject/stage, turning 20 genuine trials total\n # into several hundred.\nBATCH_SIZE = 32\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n task_eer = pd.read_csv(TASK_EER_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n print(f\"[load_dependencies] Task EER table: {len(task_eer)} tasks\")\n print(\"[load_dependencies] Task reliability (ascending EER):\")\n print(task_eer.sort_values(\"EER\")[[\"task\", \"EER\", \"ROC_AUC\"]].to_string(index=False))\n\n return manifest, task_eer, embedding_model\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Assign cascade stages by measured reliability\n# ----------------------------------------------------------------------\ndef assign_cascade_stages(task_eer, n_stages=N_STAGES):\n \"\"\"Split tasks into n_stages tiers by ascending EER (most session-stable\n first), using roughly equal-sized groups. Stage 1 = most reliable tasks,\n later stages = progressively less reliable. This replaces any assumed\n paradigm-based ordering (EO/EC -> motor -> stimulus) with an ordering\n justified directly by Module 1's measured cross-session performance.\"\"\"\n\n sorted_tasks = task_eer.sort_values(\"EER\").reset_index(drop=True)\n stage_assignment = pd.qcut(sorted_tasks.index, n_stages, labels=[f\"stage_{i+1}\" for i in range(n_stages)])\n sorted_tasks[\"stage\"] = stage_assignment\n\n stage_tasks = {\n stage: sorted_tasks[sorted_tasks[\"stage\"] == stage][\"task\"].tolist()\n for stage in sorted_tasks[\"stage\"].unique()\n }\n\n print(\"\\n[assign_cascade_stages] Stage assignment:\")\n for stage, tasks in stage_tasks.items():\n stage_eers = sorted_tasks[sorted_tasks[\"task\"].isin(tasks)][\"EER\"]\n print(f\" {stage}: tasks {tasks}, EER range [{stage_eers.min():.3f}, {stage_eers.max():.3f}]\")\n\n with open(os.path.join(OUT_DIR, \"stage_assignment.json\"), \"w\") as f:\n json.dump({k: [int(t) for t in v] for k, v in stage_tasks.items()}, f, indent=2)\n\n return stage_tasks\n\n\ndef drop_worst_tier(stage_tasks, task_eer):\n \"\"\"Drop the highest-mean-EER stage entirely from the cascade. The full\n n-stage assignment is still computed and saved by assign_cascade_stages\n (for transparency / reporting all 13 tasks' characterization), but only\n the more reliable stages are actually used in the authentication\n pipeline. This is a deliberate scope simplification: the resulting\n cascade covers a subset of tasks, not the full 13 -- state this\n explicitly in the thesis writeup rather than presenting it as full\n multi-task coverage.\"\"\"\n\n stage_mean_eer = {\n stage: task_eer[task_eer[\"task\"].isin(tasks)][\"EER\"].mean()\n for stage, tasks in stage_tasks.items()\n }\n worst_stage = max(stage_mean_eer, key=stage_mean_eer.get)\n\n kept_stage_tasks = {s: t for s, t in stage_tasks.items() if s != worst_stage}\n\n print(f\"\\n[drop_worst_tier] Dropping {worst_stage} (mean EER={stage_mean_eer[worst_stage]:.3f}, \"\n f\"tasks={stage_tasks[worst_stage]})\")\n print(f\"[drop_worst_tier] Retained stages: {list(kept_stage_tasks.keys())}\")\n\n with open(os.path.join(OUT_DIR, \"dropped_tier.json\"), \"w\") as f:\n json.dump({\n \"dropped_stage\": worst_stage,\n \"dropped_tasks\": [int(t) for t in stage_tasks[worst_stage]],\n \"dropped_mean_eer\": float(stage_mean_eer[worst_stage]),\n \"retained_stages\": {s: [int(t) for t in tasks] for s, tasks in kept_stage_tasks.items()},\n }, f, indent=2)\n\n return kept_stage_tasks\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Embedding extraction helpers (mirrors Module 1's normalization)\n# ----------------------------------------------------------------------\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Same per-sample z-score normalization as Module 1 -- must match, or\n the embedding model (trained on normalized input) will see out-of-\n distribution input and produce meaningless embeddings.\"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Batch-extract embeddings for every row in df, returning an\n (n_rows, embedding_dim) array aligned with df's row order.\"\"\"\n\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n\n embeddings = embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n return embeddings\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Build per-stage, per-subject embeddings (enrollment templates +\n# calibration probe representations)\n# ----------------------------------------------------------------------\ndef build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, split=\"enrollment\"):\n \"\"\"For each subject and stage, average ALL of that subject's epochs\n (within the given split) whose task falls in that stage. This is the\n correct approach for ENROLLMENT: a template should be one stable\n reference built from everything available, not chunked into trials.\n Returns {stage: {subject_id: mean_embedding}}.\"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_reps = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n reps = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid]\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n reps[sid] = embeddings.mean(axis=0)\n\n stage_reps[stage] = reps\n print(f\"[build_stage_templates] {split}/{stage}: {len(reps)} subjects, \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_reps\n\n\ndef build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, split=\"calibration\",\n epochs_per_trial=EPOCHS_PER_TRIAL):\n \"\"\"For each subject and stage, split that subject's epochs (within the\n given split) into chunks of epochs_per_trial, averaging each chunk into\n one independent probe \"trial\" embedding. This is the fix for the n=20\n genuine-trial statistical-power problem: instead of one averaged probe\n per subject per stage (whatever epochs exist -> one number), each\n subject now contributes multiple quasi-independent login-attempt trials.\n\n Returns {stage: {subject_id: [trial_embedding_0, trial_embedding_1, ...]}}\n \"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_trials = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n trials = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid].reset_index(drop=True)\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n\n n_chunks = max(1, len(embeddings) // epochs_per_trial)\n chunks = np.array_split(embeddings, n_chunks)\n trials[sid] = [c.mean(axis=0) for c in chunks]\n\n stage_trials[stage] = trials\n n_total_trials = sum(len(v) for v in trials.values())\n avg_trials = n_total_trials / max(len(trials), 1)\n print(f\"[build_stage_probe_trials] {split}/{stage}: {len(trials)} subjects, \"\n f\"{n_total_trials} total trials ({avg_trials:.1f} avg trials/subject), \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_trials\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Pairwise stage scores (all genuine + impostor claims)\n# ----------------------------------------------------------------------\ndef compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order):\n \"\"\"For every true_subject, every \"trial\" (simulated login attempt), and\n every claimed_subject, compute a score at every stage where both the\n trial's probe embedding and the claimed template are available.\n\n Trial alignment across stages: subject A's trial index i uses\n calibration_trials[stage][A][i] for whichever stage has that many\n trials; if stage_2 only has 12 trials but stage_1 has 18, trials 12-17\n simply have no stage_2 score (handled downstream as unavailable, not an\n error) -- this models a realistic scenario where not every login\n attempt necessarily has data for every stage.\n\n Returns a long-format DataFrame: one row per\n (trial_id, true_subject, claimed_subject, stage, score).\n \"\"\"\n\n rows = []\n\n for true_sid in cohort:\n n_trials = max(\n (len(calibration_trials.get(stage, {}).get(true_sid, [])) for stage in stage_order),\n default=0,\n )\n\n for trial_idx in range(n_trials):\n trial_id = f\"{true_sid}_trial{trial_idx}\"\n\n for stage in stage_order:\n subject_trials = calibration_trials.get(stage, {}).get(true_sid, [])\n if trial_idx >= len(subject_trials):\n continue # this stage has fewer trials for this subject\n\n probe_emb = subject_trials[trial_idx]\n\n for claimed_sid in cohort:\n template_emb = enrollment_templates.get(stage, {}).get(claimed_sid)\n if template_emb is None:\n continue\n\n score = cosine_similarity(probe_emb, template_emb)\n rows.append({\n \"trial_id\": trial_id,\n \"true_subject\": true_sid,\n \"claimed_subject\": claimed_sid,\n \"stage\": stage,\n \"is_genuine\": true_sid == claimed_sid,\n \"score\": score,\n })\n\n df = pd.DataFrame(rows)\n n_genuine_trials = df[df[\"is_genuine\"]][\"trial_id\"].nunique()\n n_impostor_pairs = df[~df[\"is_genuine\"]][[\"trial_id\", \"claimed_subject\"]].drop_duplicates().shape[0]\n print(f\"\\n[compute_pairwise_stage_scores] {len(df)} stage-level scores computed \"\n f\"({n_genuine_trials} genuine trials, {n_impostor_pairs} impostor trial-claim pairs)\")\n return df\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Compute running fused scores (mean of first k stages) per pair\n# ----------------------------------------------------------------------\ndef compute_k_stage_fused_scores(scores_df, stage_order):\n \"\"\"For every (true_subject, claimed_subject) pair, compute the running\n mean score across the first k stages, for every k = 1..len(stage_order).\n fused_k1 == the raw stage-1 score; fused_k2 == mean(stage-1, stage-2);\n etc. This is the quantity the cascade actually reasons about at each\n decision point -- not independent per-stage scores.\"\"\"\n\n pivot = scores_df.pivot_table(\n index=[\"trial_id\", \"true_subject\", \"claimed_subject\", \"is_genuine\"], columns=\"stage\", values=\"score\"\n ).reset_index()\n\n for k in range(1, len(stage_order) + 1):\n cols = stage_order[:k]\n pivot[f\"fused_k{k}\"] = pivot[cols].mean(axis=1, skipna=True)\n\n return pivot\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Calibrate three-zone thresholds per escalation level\n# ----------------------------------------------------------------------\ndef calibrate_k_thresholds(fused_df, stage_order, target_far=TARGET_FAR, reject_percentile=5):\n \"\"\"For each k, calibrate:\n - accept_threshold[k]: score above which we confidently accept, using\n the SAME target-FAR ROC calibration as before, but computed on the\n actual fused_k score distribution -- so FAR is controlled directly\n on what the cascade will actually see at that decision point,\n instead of compounding independently-calibrated per-stage bars.\n - reject_threshold[k]: a low percentile of the GENUINE fused_k score\n distribution. Scores below this are rare for genuine users, so\n rejecting here immediately (no further escalation) costs little\n FRR while stopping impostors from getting further independent\n chances -- this is the fix for the OR-gate compounding problem.\n \"\"\"\n\n accept_thresholds, reject_thresholds = {}, {}\n\n for k in range(1, len(stage_order) + 1):\n col = f\"fused_k{k}\"\n genuine = fused_df[fused_df[\"is_genuine\"]][col].dropna().values\n impostor = fused_df[~fused_df[\"is_genuine\"]][col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, thresh = roc_curve(y_true, y_score)\n\n idx_candidates = np.where(fpr <= target_far)[0]\n if len(idx_candidates) > 0:\n idx = idx_candidates[np.argmax(tpr[idx_candidates])]\n accept_thr = float(thresh[idx])\n achieved_far = fpr[idx]\n else:\n eer_idx = np.nanargmin(np.abs(fpr - (1 - tpr)))\n accept_thr = float(thresh[eer_idx])\n achieved_far = fpr[eer_idx]\n\n reject_thr = float(np.percentile(genuine, reject_percentile))\n\n accept_thresholds[k] = accept_thr\n reject_thresholds[k] = reject_thr\n\n print(f\"[calibrate_k_thresholds] k={k}: accept>={accept_thr:.4f} (achieved FAR={achieved_far:.3f}), \"\n f\"reject<={reject_thr:.4f} (bottom {reject_percentile}% of genuine)\")\n\n return accept_thresholds, reject_thresholds\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Run the fused-confidence cascade\n# ----------------------------------------------------------------------\ndef run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds):\n \"\"\"Three-zone decision at each escalation level k:\n - fused_k >= accept_threshold[k] -> ACCEPT, stop (confident match)\n - fused_k <= reject_threshold[k] -> REJECT, stop (confidently not a\n match, escalating further is pointless -- this is what prevents\n impostors from getting repeated independent chances)\n - otherwise -> escalate to k+1 (genuinely\n ambiguous, more evidence justified)\n If every stage is exhausted still in the ambiguous zone, decide by the\n final k's accept threshold (last resort, no separate looser fallback).\"\"\"\n\n results = []\n\n for _, row in fused_df.iterrows():\n decision, resolved_k = None, None\n\n for k in range(1, len(stage_order) + 1):\n score = row.get(f\"fused_k{k}\")\n if pd.isna(score) or k not in accept_thresholds:\n continue\n\n if score >= accept_thresholds[k]:\n decision, resolved_k = \"accept\", k\n break\n elif score <= reject_thresholds[k]:\n decision, resolved_k = \"reject\", k\n break\n # else: ambiguous, continue to next k\n\n if decision is None:\n last_k = max(accept_thresholds.keys())\n last_score = row.get(f\"fused_k{last_k}\")\n if pd.notna(last_score) and last_score >= accept_thresholds[last_k]:\n decision, resolved_k = \"accept\", last_k\n else:\n decision, resolved_k = \"reject\", last_k\n\n results.append({\n \"trial_id\": row[\"trial_id\"],\n \"true_subject\": row[\"true_subject\"],\n \"claimed_subject\": row[\"claimed_subject\"],\n \"is_genuine\": row[\"is_genuine\"],\n \"decision\": decision,\n \"resolved_k\": resolved_k,\n })\n\n return pd.DataFrame(results)\n\n\n# ----------------------------------------------------------------------\n# STEP 9: Evaluate cascade vs no-cascade (stage-1-only) baseline\n# ----------------------------------------------------------------------\ndef evaluate_cascade(cascade_results):\n genuine = cascade_results[cascade_results[\"is_genuine\"]]\n impostor = cascade_results[~cascade_results[\"is_genuine\"]]\n\n far = (impostor[\"decision\"] == \"accept\").mean()\n frr = (genuine[\"decision\"] == \"reject\").mean()\n\n print(f\"\\n[evaluate_cascade] Cascade FAR: {far:.4f}, FRR: {frr:.4f}\")\n\n resolution_dist = cascade_results.groupby([\"decision\"])[\"resolved_k\"].value_counts(normalize=False)\n print(\"[evaluate_cascade] Resolution counts by decision and escalation level (k):\")\n print(resolution_dist.to_string())\n\n genuine_accept_dist = genuine[genuine[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n print(\"[evaluate_cascade] Genuine users authenticated, by escalation level (k) reached:\")\n print(genuine_accept_dist.to_string())\n\n impostor_accept_dist = impostor[impostor[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n if len(impostor_accept_dist) > 0:\n print(\"[evaluate_cascade] Impostors wrongly accepted, by escalation level (k):\")\n print(impostor_accept_dist.to_string())\n\n return {\n \"far\": float(far),\n \"frr\": float(frr),\n \"genuine_resolution_by_k\": genuine_accept_dist.to_dict(),\n \"impostor_false_accept_by_k\": impostor_accept_dist.to_dict(),\n }\n\n\ndef evaluate_stage1_only_baseline(fused_df, accept_thresholds):\n \"\"\"No-cascade baseline: decide using ONLY the k=1 (first stage) fused\n score and its calibrated accept threshold, no escalation at all.\"\"\"\n\n genuine = fused_df[fused_df[\"is_genuine\"]]\n impostor = fused_df[~fused_df[\"is_genuine\"]]\n\n threshold = accept_thresholds[1]\n far = (impostor[\"fused_k1\"] >= threshold).mean()\n frr = (genuine[\"fused_k1\"] < threshold).mean()\n\n print(f\"\\n[evaluate_stage1_only_baseline] Stage-1-only FAR: {far:.4f}, FRR: {frr:.4f}\")\n return {\"far\": float(far), \"frr\": float(frr)}\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 2: Confidence-Guided Sequential Authentication ===\\n\")\n\n manifest, task_eer, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n stage_tasks_full = assign_cascade_stages(task_eer, n_stages=N_STAGES)\n stage_tasks = drop_worst_tier(stage_tasks_full, task_eer)\n stage_order = sorted(stage_tasks.keys()) # e.g. stage_1, stage_2 (worst tier dropped)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"\\n[main] Cohort: {len(cohort)} subjects\")\n\n enrollment_templates = build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, \"enrollment\")\n calibration_trials = build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, \"calibration\")\n\n scores_df = compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order)\n scores_df.to_csv(os.path.join(OUT_DIR, \"pairwise_stage_scores.csv\"), index=False)\n\n fused_df = compute_k_stage_fused_scores(scores_df, stage_order)\n fused_df.to_csv(os.path.join(OUT_DIR, \"fused_scores.csv\"), index=False)\n\n accept_thresholds, reject_thresholds = calibrate_k_thresholds(fused_df, stage_order)\n\n cascade_results = run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds)\n cascade_results.to_csv(os.path.join(OUT_DIR, \"cascade_decisions.csv\"), index=False)\n\n cascade_metrics = evaluate_cascade(cascade_results)\n baseline_metrics = evaluate_stage1_only_baseline(fused_df, accept_thresholds)\n\n comparison = pd.DataFrame([\n {\"protocol\": \"Stage-1-only (no cascade)\", **baseline_metrics},\n {\"protocol\": \"Full cascade (worst tier dropped)\",\n **{k: v for k, v in cascade_metrics.items() if k in (\"far\", \"frr\")}},\n ])\n comparison.to_csv(os.path.join(OUT_DIR, \"cascade_vs_baseline.csv\"), index=False)\n\n print(\"\\n=== Comparison: No-cascade baseline vs Full cascade ===\")\n print(comparison.to_string(index=False))\n\n with open(os.path.join(OUT_DIR, \"module2_summary.json\"), \"w\") as f:\n json.dump({\n \"stage_order\": stage_order,\n \"accept_thresholds\": accept_thresholds,\n \"reject_thresholds\": reject_thresholds,\n \"cascade_metrics\": cascade_metrics,\n \"baseline_metrics\": baseline_metrics,\n }, f, indent=2)\n\n print(f\"\\n=== Module 2 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 2: Confidence-Guided Sequential Multi-Task Authentication\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n ./module1_outputs/task_matched_diagnostic.csv\n\nRun steps in order via the __main__ block at the bottom.\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nTASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\n\nOUT_DIR = \"./module2_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nN_STAGES = 3\nTARGET_FAR = 0.05 # preferred operating point per stage; falls back to EER\n # threshold if this FAR isn't achievable on that stage's data\nEPOCHS_PER_TRIAL = 10 # calibration epochs are chunked into groups of this size;\n # each chunk becomes one independent probe \"trial\" instead\n # of averaging a subject's entire session into one probe.\n # ~150-200 calibration epochs/subject/stage -> ~15-20\n # trials/subject/stage, turning 20 genuine trials total\n # into several hundred.\nBATCH_SIZE = 32\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n task_eer = pd.read_csv(TASK_EER_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n print(f\"[load_dependencies] Task EER table: {len(task_eer)} tasks\")\n print(\"[load_dependencies] Task reliability (ascending EER):\")\n print(task_eer.sort_values(\"EER\")[[\"task\", \"EER\", \"ROC_AUC\"]].to_string(index=False))\n\n return manifest, task_eer, embedding_model\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Assign cascade stages by measured reliability\n# ----------------------------------------------------------------------\ndef assign_cascade_stages(task_eer, n_stages=N_STAGES):\n \"\"\"Split tasks into n_stages tiers by ascending EER (most session-stable\n first), using roughly equal-sized groups. Stage 1 = most reliable tasks,\n later stages = progressively less reliable. This replaces any assumed\n paradigm-based ordering (EO/EC -> motor -> stimulus) with an ordering\n justified directly by Module 1's measured cross-session performance.\"\"\"\n\n sorted_tasks = task_eer.sort_values(\"EER\").reset_index(drop=True)\n stage_assignment = pd.qcut(sorted_tasks.index, n_stages, labels=[f\"stage_{i+1}\" for i in range(n_stages)])\n sorted_tasks[\"stage\"] = stage_assignment\n\n stage_tasks = {\n stage: sorted_tasks[sorted_tasks[\"stage\"] == stage][\"task\"].tolist()\n for stage in sorted_tasks[\"stage\"].unique()\n }\n\n print(\"\\n[assign_cascade_stages] Stage assignment:\")\n for stage, tasks in stage_tasks.items():\n stage_eers = sorted_tasks[sorted_tasks[\"task\"].isin(tasks)][\"EER\"]\n print(f\" {stage}: tasks {tasks}, EER range [{stage_eers.min():.3f}, {stage_eers.max():.3f}]\")\n\n with open(os.path.join(OUT_DIR, \"stage_assignment.json\"), \"w\") as f:\n json.dump({k: [int(t) for t in v] for k, v in stage_tasks.items()}, f, indent=2)\n\n return stage_tasks\n\n\ndef drop_worst_tier(stage_tasks, task_eer):\n \"\"\"Drop the highest-mean-EER stage entirely from the cascade. The full\n n-stage assignment is still computed and saved by assign_cascade_stages\n (for transparency / reporting all 13 tasks' characterization), but only\n the more reliable stages are actually used in the authentication\n pipeline. This is a deliberate scope simplification: the resulting\n cascade covers a subset of tasks, not the full 13 -- state this\n explicitly in the thesis writeup rather than presenting it as full\n multi-task coverage.\"\"\"\n\n stage_mean_eer = {\n stage: task_eer[task_eer[\"task\"].isin(tasks)][\"EER\"].mean()\n for stage, tasks in stage_tasks.items()\n }\n worst_stage = max(stage_mean_eer, key=stage_mean_eer.get)\n\n kept_stage_tasks = {s: t for s, t in stage_tasks.items() if s != worst_stage}\n\n print(f\"\\n[drop_worst_tier] Dropping {worst_stage} (mean EER={stage_mean_eer[worst_stage]:.3f}, \"\n f\"tasks={stage_tasks[worst_stage]})\")\n print(f\"[drop_worst_tier] Retained stages: {list(kept_stage_tasks.keys())}\")\n\n with open(os.path.join(OUT_DIR, \"dropped_tier.json\"), \"w\") as f:\n json.dump({\n \"dropped_stage\": worst_stage,\n \"dropped_tasks\": [int(t) for t in stage_tasks[worst_stage]],\n \"dropped_mean_eer\": float(stage_mean_eer[worst_stage]),\n \"retained_stages\": {s: [int(t) for t in tasks] for s, tasks in kept_stage_tasks.items()},\n }, f, indent=2)\n\n return kept_stage_tasks\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Embedding extraction helpers (mirrors Module 1's normalization)\n# ----------------------------------------------------------------------\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Same per-sample z-score normalization as Module 1 -- must match, or\n the embedding model (trained on normalized input) will see out-of-\n distribution input and produce meaningless embeddings.\"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Batch-extract embeddings for every row in df, returning an\n (n_rows, embedding_dim) array aligned with df's row order.\"\"\"\n\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n\n embeddings = embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n return embeddings\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Build per-stage, per-subject embeddings (enrollment templates +\n# calibration probe representations)\n# ----------------------------------------------------------------------\ndef build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, split=\"enrollment\"):\n \"\"\"For each subject and stage, average ALL of that subject's epochs\n (within the given split) whose task falls in that stage. This is the\n correct approach for ENROLLMENT: a template should be one stable\n reference built from everything available, not chunked into trials.\n Returns {stage: {subject_id: mean_embedding}}.\"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_reps = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n reps = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid]\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n reps[sid] = embeddings.mean(axis=0)\n\n stage_reps[stage] = reps\n print(f\"[build_stage_templates] {split}/{stage}: {len(reps)} subjects, \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_reps\n\n\ndef build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, split=\"calibration\",\n epochs_per_trial=EPOCHS_PER_TRIAL):\n \"\"\"For each subject and stage, split that subject's epochs (within the\n given split) into chunks of epochs_per_trial, averaging each chunk into\n one independent probe \"trial\" embedding. This is the fix for the n=20\n genuine-trial statistical-power problem: instead of one averaged probe\n per subject per stage (whatever epochs exist -> one number), each\n subject now contributes multiple quasi-independent login-attempt trials.\n\n Returns {stage: {subject_id: [trial_embedding_0, trial_embedding_1, ...]}}\n \"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_trials = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n trials = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid].reset_index(drop=True)\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n\n n_chunks = max(1, len(embeddings) // epochs_per_trial)\n chunks = np.array_split(embeddings, n_chunks)\n trials[sid] = [c.mean(axis=0) for c in chunks]\n\n stage_trials[stage] = trials\n n_total_trials = sum(len(v) for v in trials.values())\n avg_trials = n_total_trials / max(len(trials), 1)\n print(f\"[build_stage_probe_trials] {split}/{stage}: {len(trials)} subjects, \"\n f\"{n_total_trials} total trials ({avg_trials:.1f} avg trials/subject), \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_trials\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Pairwise stage scores (all genuine + impostor claims)\n# ----------------------------------------------------------------------\ndef compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order):\n \"\"\"For every true_subject, every \"trial\" (simulated login attempt), and\n every claimed_subject, compute a score at every stage where both the\n trial's probe embedding and the claimed template are available.\n\n Trial alignment across stages: subject A's trial index i uses\n calibration_trials[stage][A][i] for whichever stage has that many\n trials; if stage_2 only has 12 trials but stage_1 has 18, trials 12-17\n simply have no stage_2 score (handled downstream as unavailable, not an\n error) -- this models a realistic scenario where not every login\n attempt necessarily has data for every stage.\n\n Returns a long-format DataFrame: one row per\n (trial_id, true_subject, claimed_subject, stage, score).\n \"\"\"\n\n rows = []\n\n for true_sid in cohort:\n n_trials = max(\n (len(calibration_trials.get(stage, {}).get(true_sid, [])) for stage in stage_order),\n default=0,\n )\n\n for trial_idx in range(n_trials):\n trial_id = f\"{true_sid}_trial{trial_idx}\"\n\n for stage in stage_order:\n subject_trials = calibration_trials.get(stage, {}).get(true_sid, [])\n if trial_idx >= len(subject_trials):\n continue # this stage has fewer trials for this subject\n\n probe_emb = subject_trials[trial_idx]\n\n for claimed_sid in cohort:\n template_emb = enrollment_templates.get(stage, {}).get(claimed_sid)\n if template_emb is None:\n continue\n\n score = cosine_similarity(probe_emb, template_emb)\n rows.append({\n \"trial_id\": trial_id,\n \"true_subject\": true_sid,\n \"claimed_subject\": claimed_sid,\n \"stage\": stage,\n \"is_genuine\": true_sid == claimed_sid,\n \"score\": score,\n })\n\n df = pd.DataFrame(rows)\n n_genuine_trials = df[df[\"is_genuine\"]][\"trial_id\"].nunique()\n n_impostor_pairs = df[~df[\"is_genuine\"]][[\"trial_id\", \"claimed_subject\"]].drop_duplicates().shape[0]\n print(f\"\\n[compute_pairwise_stage_scores] {len(df)} stage-level scores computed \"\n f\"({n_genuine_trials} genuine trials, {n_impostor_pairs} impostor trial-claim pairs)\")\n return df\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Compute running fused scores (mean of first k stages) per pair\n# ----------------------------------------------------------------------\ndef compute_k_stage_fused_scores(scores_df, stage_order):\n \"\"\"For every (true_subject, claimed_subject) pair, compute the running\n mean score across the first k stages, for every k = 1..len(stage_order).\n fused_k1 == the raw stage-1 score; fused_k2 == mean(stage-1, stage-2);\n etc. This is the quantity the cascade actually reasons about at each\n decision point -- not independent per-stage scores.\"\"\"\n\n pivot = scores_df.pivot_table(\n index=[\"trial_id\", \"true_subject\", \"claimed_subject\", \"is_genuine\"], columns=\"stage\", values=\"score\"\n ).reset_index()\n\n for k in range(1, len(stage_order) + 1):\n cols = stage_order[:k]\n pivot[f\"fused_k{k}\"] = pivot[cols].mean(axis=1, skipna=True)\n\n return pivot\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Calibrate three-zone thresholds per escalation level\n# ----------------------------------------------------------------------\ndef calibrate_k_thresholds(fused_df, stage_order, target_far=TARGET_FAR, reject_percentile=5):\n \"\"\"For each k, calibrate:\n - accept_threshold[k]: score above which we confidently accept, using\n the SAME target-FAR ROC calibration as before, computed on the\n actual fused_k score distribution.\n - reject_threshold[k]: a low percentile of the GENUINE fused_k score\n distribution. Scores below this are rare for genuine users, so\n rejecting here immediately (no further escalation) costs little\n FRR while stopping impostors from getting further independent\n chances.\n\n ORDERING GUARANTEE: reject_threshold must be strictly below\n accept_threshold, or the ambiguous escalation zone between them\n collapses to empty/inverted, silently disabling escalation entirely\n (this happened in an earlier run: good class separation pushed the\n genuine 5th-percentile ABOVE the impostor-based accept cutoff, making\n every non-accept trial automatically satisfy the reject condition too\n -- meaning the elif always fired and nothing ever reached an\n ambiguous, genuinely-escalated state). When the naturally-computed\n reject threshold would violate this ordering, it is clamped just\n below accept_threshold instead -- this is the CORRECT behavior for\n well-separated classes: a thin or near-zero ambiguous zone, not a\n broken one.\n \"\"\"\n\n accept_thresholds, reject_thresholds = {}, {}\n\n for k in range(1, len(stage_order) + 1):\n col = f\"fused_k{k}\"\n genuine = fused_df[fused_df[\"is_genuine\"]][col].dropna().values\n impostor = fused_df[~fused_df[\"is_genuine\"]][col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, thresh = roc_curve(y_true, y_score)\n\n idx_candidates = np.where(fpr <= target_far)[0]\n if len(idx_candidates) > 0:\n idx = idx_candidates[np.argmax(tpr[idx_candidates])]\n accept_thr = float(thresh[idx])\n achieved_far = fpr[idx]\n else:\n eer_idx = np.nanargmin(np.abs(fpr - (1 - tpr)))\n accept_thr = float(thresh[eer_idx])\n achieved_far = fpr[eer_idx]\n\n reject_thr_raw = float(np.percentile(genuine, reject_percentile))\n\n # Enforce ordering: reject must sit strictly below accept\n reject_thr = min(reject_thr_raw, accept_thr - 1e-4)\n was_clamped = reject_thr_raw >= accept_thr\n\n accept_thresholds[k] = accept_thr\n reject_thresholds[k] = reject_thr\n\n clamp_note = \" [CLAMPED -- classes well-separated, ambiguous zone near-empty]\" if was_clamped else \"\"\n print(f\"[calibrate_k_thresholds] k={k}: accept>={accept_thr:.4f} (achieved FAR={achieved_far:.3f}), \"\n f\"reject<={reject_thr:.4f} (raw genuine-5th-pct={reject_thr_raw:.4f}){clamp_note}\")\n\n return accept_thresholds, reject_thresholds\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Run the fused-confidence cascade\n# ----------------------------------------------------------------------\ndef run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds):\n \"\"\"Three-zone decision at each escalation level k:\n - fused_k >= accept_threshold[k] -> ACCEPT, stop (confident match)\n - fused_k <= reject_threshold[k] -> REJECT, stop (confidently not a\n match, escalating further is pointless -- this is what prevents\n impostors from getting repeated independent chances)\n - otherwise -> escalate to k+1 (genuinely\n ambiguous, more evidence justified)\n If every stage is exhausted still in the ambiguous zone, decide by the\n final k's accept threshold (last resort, no separate looser fallback).\"\"\"\n\n results = []\n\n for _, row in fused_df.iterrows():\n decision, resolved_k = None, None\n first_available_k = None\n\n for k in range(1, len(stage_order) + 1):\n score = row.get(f\"fused_k{k}\")\n if pd.isna(score) or k not in accept_thresholds:\n continue\n\n if first_available_k is None:\n first_available_k = k\n\n if score >= accept_thresholds[k]:\n decision, resolved_k = \"accept\", k\n break\n elif score <= reject_thresholds[k]:\n decision, resolved_k = \"reject\", k\n break\n # else: ambiguous, continue to next k\n\n if decision is None:\n last_k = max(accept_thresholds.keys())\n last_score = row.get(f\"fused_k{last_k}\")\n if pd.notna(last_score) and last_score >= accept_thresholds[last_k]:\n decision, resolved_k = \"accept\", last_k\n else:\n decision, resolved_k = \"reject\", last_k\n\n # Distinguish a genuine ambiguity-driven escalation (Stage-1 data\n # existed but was inconclusive) from a missing-data default (Stage-1\n # simply had no data for this trial, so resolution started later)\n was_true_escalation = (\n resolved_k is not None and first_available_k is not None\n and resolved_k > first_available_k and first_available_k == 1\n )\n\n results.append({\n \"trial_id\": row[\"trial_id\"],\n \"true_subject\": row[\"true_subject\"],\n \"claimed_subject\": row[\"claimed_subject\"],\n \"is_genuine\": row[\"is_genuine\"],\n \"decision\": decision,\n \"resolved_k\": resolved_k,\n \"first_available_k\": first_available_k,\n \"true_escalation\": was_true_escalation,\n })\n\n return pd.DataFrame(results)\n\n\n# ----------------------------------------------------------------------\n# STEP 9: Evaluate cascade vs no-cascade (stage-1-only) baseline\n# ----------------------------------------------------------------------\ndef evaluate_cascade(cascade_results):\n genuine = cascade_results[cascade_results[\"is_genuine\"]]\n impostor = cascade_results[~cascade_results[\"is_genuine\"]]\n\n far = (impostor[\"decision\"] == \"accept\").mean()\n frr = (genuine[\"decision\"] == \"reject\").mean()\n\n print(f\"\\n[evaluate_cascade] Cascade FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n\n resolution_dist = cascade_results.groupby([\"decision\"])[\"resolved_k\"].value_counts(normalize=False)\n print(\"[evaluate_cascade] Resolution counts by decision and escalation level (k):\")\n print(resolution_dist.to_string())\n\n n_missing_stage1 = (cascade_results[\"first_available_k\"] != 1).sum()\n n_true_escalations = cascade_results[\"true_escalation\"].sum()\n print(f\"[evaluate_cascade] Trials with no Stage-1 data at all: {n_missing_stage1}\")\n print(f\"[evaluate_cascade] TRUE ambiguity-driven escalations (Stage-1 present but \"\n f\"inconclusive): {n_true_escalations}\")\n\n genuine_accept_dist = genuine[genuine[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n print(\"[evaluate_cascade] Genuine users authenticated, by escalation level (k) reached \"\n \"(includes missing-data defaults, see true_escalation count above for the real signal):\")\n print(genuine_accept_dist.to_string())\n\n impostor_accept_dist = impostor[impostor[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n if len(impostor_accept_dist) > 0:\n print(\"[evaluate_cascade] Impostors wrongly accepted, by escalation level (k):\")\n print(impostor_accept_dist.to_string())\n\n return {\n \"far\": float(far),\n \"frr\": float(frr),\n \"n_genuine\": len(genuine),\n \"n_impostor\": len(impostor),\n \"genuine_resolution_by_k\": genuine_accept_dist.to_dict(),\n \"impostor_false_accept_by_k\": impostor_accept_dist.to_dict(),\n }\n\n\ndef evaluate_stage1_only_baseline(fused_df, accept_thresholds):\n \"\"\"No-cascade baseline: decide using ONLY the k=1 (first stage) fused\n score and its calibrated accept threshold, no escalation at all.\n\n IMPORTANT: a trial missing fused_k1 (no Stage-1 data available for that\n trial -- can happen since trial count is max() across stages, so some\n subjects have more Stage-2 chunks than Stage-1 chunks) cannot be\n decided by a stage-1-only system at all. A naive `NaN >= threshold`\n comparison silently evaluates False in pandas for BOTH >= and <\n checks, meaning such a trial is counted as neither accept nor reject\n in the numerator while still occupying a denominator slot -- this\n artificially DEFLATES both FAR and FRR, making the baseline look\n better than it should. These trials are dropped from the denominator\n explicitly instead, with the drop count reported so it's visible.\"\"\"\n\n valid = fused_df.dropna(subset=[\"fused_k1\"])\n n_dropped = len(fused_df) - len(valid)\n if n_dropped > 0:\n print(f\"[evaluate_stage1_only_baseline] Dropped {n_dropped} trial(s) with no \"\n f\"Stage-1 data (undecidable by a stage-1-only system) from the denominator\")\n\n genuine = valid[valid[\"is_genuine\"]]\n impostor = valid[~valid[\"is_genuine\"]]\n\n threshold = accept_thresholds[1]\n far = (impostor[\"fused_k1\"] >= threshold).mean()\n frr = (genuine[\"fused_k1\"] < threshold).mean()\n\n print(f\"\\n[evaluate_stage1_only_baseline] Stage-1-only FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n return {\"far\": float(far), \"frr\": float(frr), \"n_genuine\": len(genuine), \"n_impostor\": len(impostor)}\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 2: Confidence-Guided Sequential Authentication ===\\n\")\n\n manifest, task_eer, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n stage_tasks_full = assign_cascade_stages(task_eer, n_stages=N_STAGES)\n stage_tasks = drop_worst_tier(stage_tasks_full, task_eer)\n stage_order = sorted(stage_tasks.keys()) # e.g. stage_1, stage_2 (worst tier dropped)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"\\n[main] Cohort: {len(cohort)} subjects\")\n\n enrollment_templates = build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, \"enrollment\")\n calibration_trials = build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, \"calibration\")\n\n scores_df = compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order)\n scores_df.to_csv(os.path.join(OUT_DIR, \"pairwise_stage_scores.csv\"), index=False)\n\n fused_df = compute_k_stage_fused_scores(scores_df, stage_order)\n fused_df.to_csv(os.path.join(OUT_DIR, \"fused_scores.csv\"), index=False)\n\n accept_thresholds, reject_thresholds = calibrate_k_thresholds(fused_df, stage_order)\n\n cascade_results = run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds)\n cascade_results.to_csv(os.path.join(OUT_DIR, \"cascade_decisions.csv\"), index=False)\n\n cascade_metrics = evaluate_cascade(cascade_results)\n baseline_metrics = evaluate_stage1_only_baseline(fused_df, accept_thresholds)\n\n comparison = pd.DataFrame([\n {\"protocol\": \"Stage-1-only (no cascade)\", **baseline_metrics},\n {\"protocol\": \"Full cascade (worst tier dropped)\",\n **{k: v for k, v in cascade_metrics.items() if k in (\"far\", \"frr\", \"n_genuine\", \"n_impostor\")}},\n ])\n comparison.to_csv(os.path.join(OUT_DIR, \"cascade_vs_baseline.csv\"), index=False)\n\n print(\"\\n=== Comparison: No-cascade baseline vs Full cascade ===\")\n print(comparison.to_string(index=False))\n\n with open(os.path.join(OUT_DIR, \"module2_summary.json\"), \"w\") as f:\n json.dump({\n \"stage_order\": stage_order,\n \"accept_thresholds\": accept_thresholds,\n \"reject_thresholds\": reject_thresholds,\n \"cascade_metrics\": cascade_metrics,\n \"baseline_metrics\": baseline_metrics,\n }, f, indent=2)\n\n print(f\"\\n=== Module 2 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 2: Confidence-Guided Sequential Multi-Task Authentication\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n ./module1_outputs/task_matched_diagnostic.csv\n\nRun steps in order via the __main__ block at the bottom.\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nTASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\n\nOUT_DIR = \"./module2_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nN_STAGES = 3\nTARGET_FAR = 0.05 # preferred operating point per stage; falls back to EER\n # threshold if this FAR isn't achievable on that stage's data\nEPOCHS_PER_TRIAL = 3 # calibration epochs are chunked into groups of this size;\n # each chunk becomes one independent probe \"trial\". A\n # PRIOR VALUE of 10 over-smoothed each trial into an\n # artificially clean score (averaging 10 epochs cancels\n # out most of the noise a real single authentication\n # attempt would actually have), which produced classes\n # so well-separated that ZERO trials ever landed in the\n # ambiguous escalation zone -- meaning the cascade's\n # core mechanism was never actually exercised. A\n # smaller chunk size is both more realistic (a real\n # login attempt doesn't get 10 free epochs of\n # averaging) and gives the cascade genuine ambiguous\n # cases to resolve.\nBATCH_SIZE = 32\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n task_eer = pd.read_csv(TASK_EER_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n print(f\"[load_dependencies] Task EER table: {len(task_eer)} tasks\")\n print(\"[load_dependencies] Task reliability (ascending EER):\")\n print(task_eer.sort_values(\"EER\")[[\"task\", \"EER\", \"ROC_AUC\"]].to_string(index=False))\n\n return manifest, task_eer, embedding_model\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Assign cascade stages by measured reliability\n# ----------------------------------------------------------------------\ndef assign_cascade_stages(task_eer, n_stages=N_STAGES):\n \"\"\"Split tasks into n_stages tiers by ascending EER (most session-stable\n first), using roughly equal-sized groups. Stage 1 = most reliable tasks,\n later stages = progressively less reliable. This replaces any assumed\n paradigm-based ordering (EO/EC -> motor -> stimulus) with an ordering\n justified directly by Module 1's measured cross-session performance.\"\"\"\n\n sorted_tasks = task_eer.sort_values(\"EER\").reset_index(drop=True)\n stage_assignment = pd.qcut(sorted_tasks.index, n_stages, labels=[f\"stage_{i+1}\" for i in range(n_stages)])\n sorted_tasks[\"stage\"] = stage_assignment\n\n stage_tasks = {\n stage: sorted_tasks[sorted_tasks[\"stage\"] == stage][\"task\"].tolist()\n for stage in sorted_tasks[\"stage\"].unique()\n }\n\n print(\"\\n[assign_cascade_stages] Stage assignment:\")\n for stage, tasks in stage_tasks.items():\n stage_eers = sorted_tasks[sorted_tasks[\"task\"].isin(tasks)][\"EER\"]\n print(f\" {stage}: tasks {tasks}, EER range [{stage_eers.min():.3f}, {stage_eers.max():.3f}]\")\n\n with open(os.path.join(OUT_DIR, \"stage_assignment.json\"), \"w\") as f:\n json.dump({k: [int(t) for t in v] for k, v in stage_tasks.items()}, f, indent=2)\n\n return stage_tasks\n\n\ndef drop_worst_tier(stage_tasks, task_eer):\n \"\"\"Drop the highest-mean-EER stage entirely from the cascade. The full\n n-stage assignment is still computed and saved by assign_cascade_stages\n (for transparency / reporting all 13 tasks' characterization), but only\n the more reliable stages are actually used in the authentication\n pipeline. This is a deliberate scope simplification: the resulting\n cascade covers a subset of tasks, not the full 13 -- state this\n explicitly in the thesis writeup rather than presenting it as full\n multi-task coverage.\"\"\"\n\n stage_mean_eer = {\n stage: task_eer[task_eer[\"task\"].isin(tasks)][\"EER\"].mean()\n for stage, tasks in stage_tasks.items()\n }\n worst_stage = max(stage_mean_eer, key=stage_mean_eer.get)\n\n kept_stage_tasks = {s: t for s, t in stage_tasks.items() if s != worst_stage}\n\n print(f\"\\n[drop_worst_tier] Dropping {worst_stage} (mean EER={stage_mean_eer[worst_stage]:.3f}, \"\n f\"tasks={stage_tasks[worst_stage]})\")\n print(f\"[drop_worst_tier] Retained stages: {list(kept_stage_tasks.keys())}\")\n\n with open(os.path.join(OUT_DIR, \"dropped_tier.json\"), \"w\") as f:\n json.dump({\n \"dropped_stage\": worst_stage,\n \"dropped_tasks\": [int(t) for t in stage_tasks[worst_stage]],\n \"dropped_mean_eer\": float(stage_mean_eer[worst_stage]),\n \"retained_stages\": {s: [int(t) for t in tasks] for s, tasks in kept_stage_tasks.items()},\n }, f, indent=2)\n\n return kept_stage_tasks\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Embedding extraction helpers (mirrors Module 1's normalization)\n# ----------------------------------------------------------------------\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Same per-sample z-score normalization as Module 1 -- must match, or\n the embedding model (trained on normalized input) will see out-of-\n distribution input and produce meaningless embeddings.\"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n \"\"\"Batch-extract embeddings for every row in df, returning an\n (n_rows, embedding_dim) array aligned with df's row order.\"\"\"\n\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n\n embeddings = embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n return embeddings\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Build per-stage, per-subject embeddings (enrollment templates +\n# calibration probe representations)\n# ----------------------------------------------------------------------\ndef build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, split=\"enrollment\"):\n \"\"\"For each subject and stage, average ALL of that subject's epochs\n (within the given split) whose task falls in that stage. This is the\n correct approach for ENROLLMENT: a template should be one stable\n reference built from everything available, not chunked into trials.\n Returns {stage: {subject_id: mean_embedding}}.\"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_reps = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n reps = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid]\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n reps[sid] = embeddings.mean(axis=0)\n\n stage_reps[stage] = reps\n print(f\"[build_stage_templates] {split}/{stage}: {len(reps)} subjects, \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_reps\n\n\ndef build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, split=\"calibration\",\n epochs_per_trial=EPOCHS_PER_TRIAL):\n \"\"\"For each subject and stage, split that subject's epochs (within the\n given split) into chunks of epochs_per_trial, averaging each chunk into\n one independent probe \"trial\" embedding. This is the fix for the n=20\n genuine-trial statistical-power problem: instead of one averaged probe\n per subject per stage (whatever epochs exist -> one number), each\n subject now contributes multiple quasi-independent login-attempt trials.\n\n Returns {stage: {subject_id: [trial_embedding_0, trial_embedding_1, ...]}}\n \"\"\"\n\n split_df = manifest[manifest[\"split\"] == split]\n stage_trials = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n trials = {}\n\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid].reset_index(drop=True)\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n\n n_chunks = max(1, len(embeddings) // epochs_per_trial)\n chunks = np.array_split(embeddings, n_chunks)\n trials[sid] = [c.mean(axis=0) for c in chunks]\n\n stage_trials[stage] = trials\n n_total_trials = sum(len(v) for v in trials.values())\n avg_trials = n_total_trials / max(len(trials), 1)\n print(f\"[build_stage_probe_trials] {split}/{stage}: {len(trials)} subjects, \"\n f\"{n_total_trials} total trials ({avg_trials:.1f} avg trials/subject), \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_trials\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Pairwise stage scores (all genuine + impostor claims)\n# ----------------------------------------------------------------------\ndef compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order):\n \"\"\"For every true_subject, every \"trial\" (simulated login attempt), and\n every claimed_subject, compute a score at every stage where both the\n trial's probe embedding and the claimed template are available.\n\n Trial alignment across stages: subject A's trial index i uses\n calibration_trials[stage][A][i] for whichever stage has that many\n trials; if stage_2 only has 12 trials but stage_1 has 18, trials 12-17\n simply have no stage_2 score (handled downstream as unavailable, not an\n error) -- this models a realistic scenario where not every login\n attempt necessarily has data for every stage.\n\n Returns a long-format DataFrame: one row per\n (trial_id, true_subject, claimed_subject, stage, score).\n \"\"\"\n\n rows = []\n\n for true_sid in cohort:\n n_trials = max(\n (len(calibration_trials.get(stage, {}).get(true_sid, [])) for stage in stage_order),\n default=0,\n )\n\n for trial_idx in range(n_trials):\n trial_id = f\"{true_sid}_trial{trial_idx}\"\n\n for stage in stage_order:\n subject_trials = calibration_trials.get(stage, {}).get(true_sid, [])\n if trial_idx >= len(subject_trials):\n continue # this stage has fewer trials for this subject\n\n probe_emb = subject_trials[trial_idx]\n\n for claimed_sid in cohort:\n template_emb = enrollment_templates.get(stage, {}).get(claimed_sid)\n if template_emb is None:\n continue\n\n score = cosine_similarity(probe_emb, template_emb)\n rows.append({\n \"trial_id\": trial_id,\n \"true_subject\": true_sid,\n \"claimed_subject\": claimed_sid,\n \"stage\": stage,\n \"is_genuine\": true_sid == claimed_sid,\n \"score\": score,\n })\n\n df = pd.DataFrame(rows)\n n_genuine_trials = df[df[\"is_genuine\"]][\"trial_id\"].nunique()\n n_impostor_pairs = df[~df[\"is_genuine\"]][[\"trial_id\", \"claimed_subject\"]].drop_duplicates().shape[0]\n print(f\"\\n[compute_pairwise_stage_scores] {len(df)} stage-level scores computed \"\n f\"({n_genuine_trials} genuine trials, {n_impostor_pairs} impostor trial-claim pairs)\")\n return df\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Compute running fused scores (mean of first k stages) per pair\n# ----------------------------------------------------------------------\ndef compute_k_stage_fused_scores(scores_df, stage_order):\n \"\"\"For every (true_subject, claimed_subject) pair, compute the running\n mean score across the first k stages, for every k = 1..len(stage_order).\n fused_k1 == the raw stage-1 score; fused_k2 == mean(stage-1, stage-2);\n etc. This is the quantity the cascade actually reasons about at each\n decision point -- not independent per-stage scores.\"\"\"\n\n pivot = scores_df.pivot_table(\n index=[\"trial_id\", \"true_subject\", \"claimed_subject\", \"is_genuine\"], columns=\"stage\", values=\"score\"\n ).reset_index()\n\n for k in range(1, len(stage_order) + 1):\n cols = stage_order[:k]\n pivot[f\"fused_k{k}\"] = pivot[cols].mean(axis=1, skipna=True)\n\n return pivot\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Calibrate three-zone thresholds per escalation level\n# ----------------------------------------------------------------------\ndef calibrate_k_thresholds(fused_df, stage_order, target_far=TARGET_FAR, reject_percentile=5):\n \"\"\"For each k, calibrate:\n - accept_threshold[k]: score above which we confidently accept, using\n the SAME target-FAR ROC calibration as before, computed on the\n actual fused_k score distribution.\n - reject_threshold[k]: a low percentile of the GENUINE fused_k score\n distribution. Scores below this are rare for genuine users, so\n rejecting here immediately (no further escalation) costs little\n FRR while stopping impostors from getting further independent\n chances.\n\n ORDERING GUARANTEE: reject_threshold must be strictly below\n accept_threshold, or the ambiguous escalation zone between them\n collapses to empty/inverted, silently disabling escalation entirely\n (this happened in an earlier run: good class separation pushed the\n genuine 5th-percentile ABOVE the impostor-based accept cutoff, making\n every non-accept trial automatically satisfy the reject condition too\n -- meaning the elif always fired and nothing ever reached an\n ambiguous, genuinely-escalated state). When the naturally-computed\n reject threshold would violate this ordering, it is clamped just\n below accept_threshold instead -- this is the CORRECT behavior for\n well-separated classes: a thin or near-zero ambiguous zone, not a\n broken one.\n \"\"\"\n\n accept_thresholds, reject_thresholds = {}, {}\n\n for k in range(1, len(stage_order) + 1):\n col = f\"fused_k{k}\"\n genuine = fused_df[fused_df[\"is_genuine\"]][col].dropna().values\n impostor = fused_df[~fused_df[\"is_genuine\"]][col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, thresh = roc_curve(y_true, y_score)\n\n idx_candidates = np.where(fpr <= target_far)[0]\n if len(idx_candidates) > 0:\n idx = idx_candidates[np.argmax(tpr[idx_candidates])]\n accept_thr = float(thresh[idx])\n achieved_far = fpr[idx]\n else:\n eer_idx = np.nanargmin(np.abs(fpr - (1 - tpr)))\n accept_thr = float(thresh[eer_idx])\n achieved_far = fpr[eer_idx]\n\n reject_thr_raw = float(np.percentile(genuine, reject_percentile))\n\n # Enforce ordering: reject must sit strictly below accept\n reject_thr = min(reject_thr_raw, accept_thr - 1e-4)\n was_clamped = reject_thr_raw >= accept_thr\n\n accept_thresholds[k] = accept_thr\n reject_thresholds[k] = reject_thr\n\n clamp_note = \" [CLAMPED -- classes well-separated, ambiguous zone near-empty]\" if was_clamped else \"\"\n print(f\"[calibrate_k_thresholds] k={k}: accept>={accept_thr:.4f} (achieved FAR={achieved_far:.3f}), \"\n f\"reject<={reject_thr:.4f} (raw genuine-5th-pct={reject_thr_raw:.4f}){clamp_note}\")\n\n return accept_thresholds, reject_thresholds\n\n\n# ----------------------------------------------------------------------\n# STEP 8: Run the fused-confidence cascade\n# ----------------------------------------------------------------------\ndef run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds):\n \"\"\"Three-zone decision at each escalation level k:\n - fused_k >= accept_threshold[k] -> ACCEPT, stop (confident match)\n - fused_k <= reject_threshold[k] -> REJECT, stop (confidently not a\n match, escalating further is pointless -- this is what prevents\n impostors from getting repeated independent chances)\n - otherwise -> escalate to k+1 (genuinely\n ambiguous, more evidence justified)\n If every stage is exhausted still in the ambiguous zone, decide by the\n final k's accept threshold (last resort, no separate looser fallback).\"\"\"\n\n results = []\n\n for _, row in fused_df.iterrows():\n decision, resolved_k = None, None\n first_available_k = None\n\n for k in range(1, len(stage_order) + 1):\n score = row.get(f\"fused_k{k}\")\n if pd.isna(score) or k not in accept_thresholds:\n continue\n\n if first_available_k is None:\n first_available_k = k\n\n if score >= accept_thresholds[k]:\n decision, resolved_k = \"accept\", k\n break\n elif score <= reject_thresholds[k]:\n decision, resolved_k = \"reject\", k\n break\n # else: ambiguous, continue to next k\n\n if decision is None:\n last_k = max(accept_thresholds.keys())\n last_score = row.get(f\"fused_k{last_k}\")\n if pd.notna(last_score) and last_score >= accept_thresholds[last_k]:\n decision, resolved_k = \"accept\", last_k\n else:\n decision, resolved_k = \"reject\", last_k\n\n # Distinguish a genuine ambiguity-driven escalation (Stage-1 data\n # existed but was inconclusive) from a missing-data default (Stage-1\n # simply had no data for this trial, so resolution started later)\n was_true_escalation = (\n resolved_k is not None and first_available_k is not None\n and resolved_k > first_available_k and first_available_k == 1\n )\n\n results.append({\n \"trial_id\": row[\"trial_id\"],\n \"true_subject\": row[\"true_subject\"],\n \"claimed_subject\": row[\"claimed_subject\"],\n \"is_genuine\": row[\"is_genuine\"],\n \"decision\": decision,\n \"resolved_k\": resolved_k,\n \"first_available_k\": first_available_k,\n \"true_escalation\": was_true_escalation,\n })\n\n return pd.DataFrame(results)\n\n\ndef evaluate_cascade_by_data_availability(cascade_results):\n \"\"\"Split cascade results into two groups for a fair comparison:\n - matched: trials that HAD Stage-1 data -- directly comparable to the\n Stage-1-only baseline, since both are deciding the exact same\n trials. If true_escalation counts are near zero, expect this\n subgroup's FAR/FRR to closely match the baseline.\n - coverage_extension: trials with NO Stage-1 data at all -- these are\n trials the baseline literally cannot decide (dropped from its\n denominator). Cascade decides them via Stage-2 alone. Report this\n separately rather than blending it into one number, since blending\n makes cascade look worse than baseline for reasons that have\n nothing to do with escalation quality.\"\"\"\n\n matched = cascade_results[cascade_results[\"first_available_k\"] == 1]\n extension = cascade_results[cascade_results[\"first_available_k\"] != 1]\n\n def _far_frr(df):\n genuine = df[df[\"is_genuine\"]]\n impostor = df[~df[\"is_genuine\"]]\n far = (impostor[\"decision\"] == \"accept\").mean() if len(impostor) > 0 else float(\"nan\")\n frr = (genuine[\"decision\"] == \"reject\").mean() if len(genuine) > 0 else float(\"nan\")\n return far, frr, len(genuine), len(impostor)\n\n m_far, m_frr, m_ng, m_ni = _far_frr(matched)\n e_far, e_frr, e_ng, e_ni = _far_frr(extension)\n\n print(\"\\n[evaluate_cascade_by_data_availability] Matched population (has Stage-1 data, \"\n \"directly comparable to baseline):\")\n print(f\" FAR={m_far:.4f}, FRR={m_frr:.4f} (n_genuine={m_ng}, n_impostor={m_ni})\")\n print(\"[evaluate_cascade_by_data_availability] Coverage-extension population (no Stage-1 \"\n \"data, undecidable by baseline, decided via Stage-2 fallback):\")\n print(f\" FAR={e_far:.4f}, FRR={e_frr:.4f} (n_genuine={e_ng}, n_impostor={e_ni})\")\n\n return {\n \"matched\": {\"far\": float(m_far), \"frr\": float(m_frr), \"n_genuine\": m_ng, \"n_impostor\": m_ni},\n \"coverage_extension\": {\"far\": float(e_far), \"frr\": float(e_frr), \"n_genuine\": e_ng, \"n_impostor\": e_ni},\n }\n\n\n# ----------------------------------------------------------------------\n# STEP 10: Evaluate cascade vs no-cascade (stage-1-only) baseline\n# ----------------------------------------------------------------------\ndef evaluate_cascade(cascade_results):\n genuine = cascade_results[cascade_results[\"is_genuine\"]]\n impostor = cascade_results[~cascade_results[\"is_genuine\"]]\n\n far = (impostor[\"decision\"] == \"accept\").mean()\n frr = (genuine[\"decision\"] == \"reject\").mean()\n\n print(f\"\\n[evaluate_cascade] Cascade FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n\n resolution_dist = cascade_results.groupby([\"decision\"])[\"resolved_k\"].value_counts(normalize=False)\n print(\"[evaluate_cascade] Resolution counts by decision and escalation level (k):\")\n print(resolution_dist.to_string())\n\n n_missing_stage1 = (cascade_results[\"first_available_k\"] != 1).sum()\n n_true_escalations = cascade_results[\"true_escalation\"].sum()\n print(f\"[evaluate_cascade] Trials with no Stage-1 data at all: {n_missing_stage1}\")\n print(f\"[evaluate_cascade] TRUE ambiguity-driven escalations (Stage-1 present but \"\n f\"inconclusive): {n_true_escalations}\")\n\n genuine_accept_dist = genuine[genuine[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n print(\"[evaluate_cascade] Genuine users authenticated, by escalation level (k) reached \"\n \"(includes missing-data defaults, see true_escalation count above for the real signal):\")\n print(genuine_accept_dist.to_string())\n\n impostor_accept_dist = impostor[impostor[\"decision\"] == \"accept\"][\"resolved_k\"].value_counts(normalize=True)\n if len(impostor_accept_dist) > 0:\n print(\"[evaluate_cascade] Impostors wrongly accepted, by escalation level (k):\")\n print(impostor_accept_dist.to_string())\n\n return {\n \"far\": float(far),\n \"frr\": float(frr),\n \"n_genuine\": len(genuine),\n \"n_impostor\": len(impostor),\n \"genuine_resolution_by_k\": genuine_accept_dist.to_dict(),\n \"impostor_false_accept_by_k\": impostor_accept_dist.to_dict(),\n }\n\n\ndef evaluate_stage1_only_baseline(fused_df, accept_thresholds):\n \"\"\"No-cascade baseline: decide using ONLY the k=1 (first stage) fused\n score and its calibrated accept threshold, no escalation at all.\n\n IMPORTANT: a trial missing fused_k1 (no Stage-1 data available for that\n trial -- can happen since trial count is max() across stages, so some\n subjects have more Stage-2 chunks than Stage-1 chunks) cannot be\n decided by a stage-1-only system at all. A naive `NaN >= threshold`\n comparison silently evaluates False in pandas for BOTH >= and <\n checks, meaning such a trial is counted as neither accept nor reject\n in the numerator while still occupying a denominator slot -- this\n artificially DEFLATES both FAR and FRR, making the baseline look\n better than it should. These trials are dropped from the denominator\n explicitly instead, with the drop count reported so it's visible.\"\"\"\n\n valid = fused_df.dropna(subset=[\"fused_k1\"])\n n_dropped = len(fused_df) - len(valid)\n if n_dropped > 0:\n print(f\"[evaluate_stage1_only_baseline] Dropped {n_dropped} trial(s) with no \"\n f\"Stage-1 data (undecidable by a stage-1-only system) from the denominator\")\n\n genuine = valid[valid[\"is_genuine\"]]\n impostor = valid[~valid[\"is_genuine\"]]\n\n threshold = accept_thresholds[1]\n far = (impostor[\"fused_k1\"] >= threshold).mean()\n frr = (genuine[\"fused_k1\"] < threshold).mean()\n\n print(f\"\\n[evaluate_stage1_only_baseline] Stage-1-only FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n return {\"far\": float(far), \"frr\": float(frr), \"n_genuine\": len(genuine), \"n_impostor\": len(impostor)}\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 2: Confidence-Guided Sequential Authentication ===\\n\")\n\n manifest, task_eer, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n stage_tasks_full = assign_cascade_stages(task_eer, n_stages=N_STAGES)\n stage_tasks = drop_worst_tier(stage_tasks_full, task_eer)\n stage_order = sorted(stage_tasks.keys()) # e.g. stage_1, stage_2 (worst tier dropped)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"\\n[main] Cohort: {len(cohort)} subjects\")\n\n enrollment_templates = build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, \"enrollment\")\n calibration_trials = build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, \"calibration\")\n\n scores_df = compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order)\n scores_df.to_csv(os.path.join(OUT_DIR, \"pairwise_stage_scores.csv\"), index=False)\n\n fused_df = compute_k_stage_fused_scores(scores_df, stage_order)\n fused_df.to_csv(os.path.join(OUT_DIR, \"fused_scores.csv\"), index=False)\n\n accept_thresholds, reject_thresholds = calibrate_k_thresholds(fused_df, stage_order)\n\n cascade_results = run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds)\n cascade_results.to_csv(os.path.join(OUT_DIR, \"cascade_decisions.csv\"), index=False)\n\n cascade_metrics = evaluate_cascade(cascade_results)\n subgroup_metrics = evaluate_cascade_by_data_availability(cascade_results)\n baseline_metrics = evaluate_stage1_only_baseline(fused_df, accept_thresholds)\n\n comparison = pd.DataFrame([\n {\"protocol\": \"Stage-1-only (no cascade)\", **baseline_metrics},\n {\"protocol\": \"Full cascade (worst tier dropped)\",\n **{k: v for k, v in cascade_metrics.items() if k in (\"far\", \"frr\", \"n_genuine\", \"n_impostor\")}},\n ])\n comparison.to_csv(os.path.join(OUT_DIR, \"cascade_vs_baseline.csv\"), index=False)\n\n print(\"\\n=== Comparison: No-cascade baseline vs Full cascade ===\")\n print(comparison.to_string(index=False))\n\n with open(os.path.join(OUT_DIR, \"module2_summary.json\"), \"w\") as f:\n json.dump({\n \"stage_order\": stage_order,\n \"accept_thresholds\": accept_thresholds,\n \"reject_thresholds\": reject_thresholds,\n \"cascade_metrics\": cascade_metrics,\n \"subgroup_metrics\": subgroup_metrics,\n \"baseline_metrics\": baseline_metrics,\n }, f, indent=2)\n\n print(f\"\\n=== Module 2 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 3: Task-Aware Metric-Learning Model for Robust EEG Identity Representation\nM3CV EEG Biometric Authentication Pipeline\n\nDepends on Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\nCompares against Module 1 outputs: ./module1_outputs/task_matched_diagnostic.csv\n\nRun steps in order via the __main__ block at the bottom.\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import roc_curve, roc_auc_score\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nMODULE1_TASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\n\nOUT_DIR = \"./module3_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nEMBEDDING_DIM = 64\nTASK_EMBED_DIM = 8\nMARGIN = 0.3 # triplet loss margin, in squared-L2-distance units for\n # L2-normalized embeddings (range [0,4]; 0.3 is a\n # standard starting point, roughly a cosine-similarity\n # gap of ~0.15)\nP_SUBJECTS = 10 # subjects per training batch\nK_SAMPLES = 4 # epochs per subject per batch -> batch size = P*K = 40\nSTEPS_PER_EPOCH = 60\nEPOCHS = 30\nLEARNING_RATE = 0.0005\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\ntf.random.set_seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load manifest and spectrogram helpers (must match Module 1 exactly)\n# ----------------------------------------------------------------------\ndef load_manifest():\n manifest = pd.read_csv(MANIFEST_PATH)\n print(f\"[load_manifest] Loaded {len(manifest)} entries\")\n return manifest\n\n\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Identical to Module 1's normalization -- per-sample z-score. This is\n not optional (see Module 1's dead-ReLU collapse incident): unnormalized\n dB-range input reliably breaks CNN training regardless of loss function.\"\"\"\n\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# STEP 2: PK batch sampler (P subjects x K epochs per batch)\n# ----------------------------------------------------------------------\ndef build_pk_generator(df, target_shape, label_encoder, task_encoder,\n p_subjects=P_SUBJECTS, k_samples=K_SAMPLES, steps_per_epoch=STEPS_PER_EPOCH):\n \"\"\"Yields batches of (p_subjects * k_samples) spectrograms, each batch\n containing exactly k_samples epochs from each of p_subjects randomly\n chosen subjects -- required for batch-hard triplet loss to have valid\n positive AND negative pairs within every batch.\"\"\"\n\n subject_rows = {\n sid: list(zip(g[\"cache_path\"], g[\"task\"]))\n for sid, g in df.groupby(\"subject_id\")\n }\n subjects = list(subject_rows.keys())\n\n def gen():\n rng = np.random.RandomState(RANDOM_STATE)\n for _ in range(steps_per_epoch):\n chosen_subjects = rng.choice(subjects, size=min(p_subjects, len(subjects)), replace=False)\n\n specs, tasks, labels = [], [], []\n for sid in chosen_subjects:\n rows = subject_rows[sid]\n sampled_idx = rng.choice(len(rows), size=k_samples, replace=len(rows) < k_samples)\n for idx in sampled_idx:\n cache_path, task = rows[idx]\n specs.append(load_spectrogram_array(cache_path, target_shape))\n tasks.append(task)\n labels.append(sid)\n\n specs = np.stack(specs)[..., np.newaxis].astype(np.float32)\n tasks = np.array([task_encoder[t] for t in tasks], dtype=np.int32)\n labels = np.array(label_encoder.transform(labels), dtype=np.int32)\n\n yield (specs, tasks), labels\n\n batch_size = p_subjects * k_samples\n output_signature = (\n (\n tf.TensorSpec(shape=(batch_size, *target_shape, 1), dtype=tf.float32),\n tf.TensorSpec(shape=(batch_size,), dtype=tf.int32),\n ),\n tf.TensorSpec(shape=(batch_size,), dtype=tf.int32),\n )\n\n return tf.data.Dataset.from_generator(gen, output_signature=output_signature).prefetch(tf.data.AUTOTUNE)\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Batch-hard triplet loss\n# ----------------------------------------------------------------------\ndef pairwise_squared_distances(embeddings):\n \"\"\"Squared Euclidean distance matrix for L2-normalized embeddings.\n For unit vectors, ||a-b||^2 = 2 - 2*cos(a,b), so this is a monotonic\n transform of cosine similarity -- range [0, 4].\"\"\"\n\n dot = tf.matmul(embeddings, embeddings, transpose_b=True)\n sq_norms = tf.linalg.diag_part(dot)\n distances = tf.expand_dims(sq_norms, 1) - 2.0 * dot + tf.expand_dims(sq_norms, 0)\n return tf.maximum(distances, 0.0)\n\n\ndef batch_hard_triplet_loss(labels, embeddings, margin=MARGIN):\n \"\"\"For each anchor in the batch: find its HARDEST positive (same\n identity, max distance) and HARDEST negative (different identity, min\n distance). Loss = mean(relu(margin + hardest_positive - hardest_negative)).\n This is the standard Hermans et al. 2017 formulation -- focuses learning\n on the pairs that are actually difficult, rather than averaging over\n many trivially-easy triplets.\"\"\"\n\n labels = tf.reshape(labels, [-1, 1])\n pdist = pairwise_squared_distances(embeddings)\n\n adjacency = tf.equal(labels, tf.transpose(labels))\n batch_size = tf.shape(embeddings)[0]\n not_self = tf.logical_not(tf.eye(batch_size, dtype=tf.bool))\n\n mask_positive = tf.logical_and(adjacency, not_self)\n pdist_pos = tf.where(mask_positive, pdist, tf.zeros_like(pdist))\n hardest_positive = tf.reduce_max(pdist_pos, axis=1)\n\n mask_negative = tf.logical_not(adjacency)\n max_dist = tf.reduce_max(pdist)\n pdist_neg = tf.where(mask_negative, pdist, max_dist * tf.ones_like(pdist))\n hardest_negative = tf.reduce_min(pdist_neg, axis=1)\n\n loss = tf.maximum(hardest_positive - hardest_negative + margin, 0.0)\n return tf.reduce_mean(loss)\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Task-conditioned embedding model\n# ----------------------------------------------------------------------\ndef build_metric_model(input_shape, n_tasks, embedding_dim=EMBEDDING_DIM, task_embed_dim=TASK_EMBED_DIM):\n \"\"\"Same CNN backbone family as Module 1 (Conv-Pool-Conv-Pool-Dropout-\n Dense), but with a task-ID embedding concatenated in before the final\n embedding layer -- letting the network condition its identity\n representation on which paradigm (resting/motor/stimulus) produced the\n input, per Contribution 3's task-aware design. Output is L2-normalized\n (required for the triplet loss's distance formulation, and standard\n practice for metric-learning embeddings generally).\"\"\"\n\n spec_input = tf.keras.Input(shape=input_shape, name=\"spectrogram\")\n task_input = tf.keras.Input(shape=(), dtype=tf.int32, name=\"task_id\")\n\n x = tf.keras.layers.Conv2D(32, 3, activation=\"relu\")(spec_input)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Conv2D(64, 3, activation=\"relu\")(x)\n x = tf.keras.layers.MaxPooling2D(2)(x)\n x = tf.keras.layers.Dropout(0.2)(x)\n x = tf.keras.layers.Flatten()(x)\n x = tf.keras.layers.Dense(128, activation=\"relu\")(x)\n\n task_embed = tf.keras.layers.Embedding(input_dim=n_tasks + 1, output_dim=task_embed_dim)(task_input)\n task_embed = tf.keras.layers.Reshape((task_embed_dim,))(task_embed)\n\n combined = tf.keras.layers.Concatenate()([x, task_embed])\n embedding = tf.keras.layers.Dense(embedding_dim)(combined)\n embedding = tf.keras.layers.Lambda(\n lambda t: tf.math.l2_normalize(t, axis=1), name=\"l2_normalize\"\n )(embedding)\n\n model = tf.keras.Model(inputs=[spec_input, task_input], outputs=embedding)\n return model\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Train\n# ----------------------------------------------------------------------\ndef train_metric_model(train_df, val_df, target_shape, label_encoder, task_encoder, n_tasks):\n model = build_metric_model((*target_shape, 1), n_tasks)\n\n model.compile(\n optimizer=tf.keras.optimizers.Adam(learning_rate=LEARNING_RATE),\n loss=lambda y_true, y_pred: batch_hard_triplet_loss(y_true, y_pred, margin=MARGIN),\n )\n\n train_ds = build_pk_generator(train_df, target_shape, label_encoder, task_encoder)\n val_ds = build_pk_generator(val_df, target_shape, label_encoder, task_encoder,\n steps_per_epoch=max(10, STEPS_PER_EPOCH // 4))\n\n history = model.fit(\n train_ds,\n validation_data=val_ds,\n epochs=EPOCHS,\n steps_per_epoch=STEPS_PER_EPOCH,\n validation_steps=max(10, STEPS_PER_EPOCH // 4),\n verbose=1,\n )\n\n model.save(os.path.join(OUT_DIR, \"metric_embedding_model.keras\"))\n\n with open(os.path.join(OUT_DIR, \"training_history.json\"), \"w\") as f:\n json.dump({k: [float(v) for v in vals] for k, vals in history.history.items()}, f, indent=2)\n\n return model\n\n\ndef check_embedding_diversity(model, df, target_shape, task_encoder, n_check=30):\n \"\"\"Same fail-fast check as Module 1 -- catch a collapsed embedding\n model in seconds, not after a full verification run.\"\"\"\n\n sample = df.sample(min(n_check, len(df)), random_state=RANDOM_STATE)\n specs = np.stack([load_spectrogram_array(p, target_shape) for p in sample[\"cache_path\"]])[..., np.newaxis]\n tasks = np.array([task_encoder[t] for t in sample[\"task\"]], dtype=np.int32)\n\n embeddings = model.predict([specs.astype(np.float32), tasks], verbose=0)\n pairwise_std = embeddings.std(axis=0).mean()\n print(f\"\\n[check_embedding_diversity] Mean per-dimension std: {pairwise_std:.6f}\")\n\n if pairwise_std < 1e-3:\n raise RuntimeError(\n f\"[check_embedding_diversity] EMBEDDING MODEL HAS COLLAPSED (std={pairwise_std:.6f}). \"\n \"Do not proceed to verification -- check triplet loss is decreasing in training_history.json.\"\n )\n print(\"[check_embedding_diversity] OK: embeddings show real variation.\")\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Extract embeddings, build templates, run task-matched verification\n# (identical protocol to Module 1's task-matched diagnostic, for a\n# direct ablation comparison)\n# ----------------------------------------------------------------------\ndef extract_embeddings(model, df, target_shape, task_encoder, batch_size=64):\n specs = np.stack([load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]])[..., np.newaxis]\n tasks = np.array([task_encoder[t] for t in df[\"task\"]], dtype=np.int32)\n return model.predict([specs.astype(np.float32), tasks], batch_size=batch_size, verbose=0)\n\n\ndef build_subject_templates(embeddings, subject_ids):\n templates = {}\n for sid in np.unique(subject_ids):\n mask = subject_ids == sid\n templates[sid] = embeddings[mask].mean(axis=0)\n return templates\n\n\ndef compute_verification_metrics(genuine_scores, impostor_scores, label=\"\"):\n y_true = np.concatenate([np.ones_like(genuine_scores), np.zeros_like(impostor_scores)])\n y_score = np.concatenate([genuine_scores, impostor_scores])\n\n fpr, tpr, thresholds = roc_curve(y_true, y_score)\n fnr = 1 - tpr\n eer_idx = np.nanargmin(np.abs(fpr - fnr))\n eer = (fpr[eer_idx] + fnr[eer_idx]) / 2\n auc = roc_auc_score(y_true, y_score)\n\n return {\"label\": label, \"eer\": float(eer), \"roc_auc\": float(auc)}\n\n\ndef run_task_matched_verification(manifest, model, target_shape, task_encoder, cohort):\n \"\"\"Identical protocol to Module 1's run_task_matched_diagnostic: for\n each task present in both splits, build per-task templates from\n Enrollment and score Calibration probes against them, cross-session,\n same task both sides. Produces a directly comparable table to Module\n 1's task_matched_diagnostic.csv.\"\"\"\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n\n results = []\n\n for task in common_tasks:\n enrollment_task = enrollment[(enrollment[\"task\"] == task) & (enrollment[\"subject_id\"].isin(cohort))]\n calibration_task = calibration[(calibration[\"task\"] == task) & (calibration[\"subject_id\"].isin(cohort))]\n\n subjects_in_both = set(enrollment_task[\"subject_id\"]) & set(calibration_task[\"subject_id\"])\n if len(subjects_in_both) < 2:\n continue\n\n enrollment_task = enrollment_task[enrollment_task[\"subject_id\"].isin(subjects_in_both)]\n calibration_task = calibration_task[calibration_task[\"subject_id\"].isin(subjects_in_both)]\n\n template_embeddings = extract_embeddings(model, enrollment_task, target_shape, task_encoder)\n templates = build_subject_templates(template_embeddings, enrollment_task[\"subject_id\"].values)\n\n probe_embeddings = extract_embeddings(model, calibration_task, target_shape, task_encoder)\n\n genuine_scores, impostor_scores = [], []\n for emb, true_sid in zip(probe_embeddings, calibration_task[\"subject_id\"].values):\n for sid, template in templates.items():\n score = cosine_similarity(emb, template)\n (genuine_scores if sid == true_sid else impostor_scores).append(score)\n\n if len(genuine_scores) == 0 or len(impostor_scores) == 0:\n continue\n\n result = compute_verification_metrics(\n np.array(genuine_scores), np.array(impostor_scores), label=f\"Task {task}\"\n )\n result[\"task\"] = int(task)\n result[\"n_subjects\"] = len(subjects_in_both)\n results.append(result)\n\n task_df = pd.DataFrame([{\"task\": r[\"task\"], \"n_subjects\": r[\"n_subjects\"],\n \"EER\": r[\"eer\"], \"ROC_AUC\": r[\"roc_auc\"]} for r in results])\n task_df = task_df.sort_values(\"EER\").reset_index(drop=True)\n task_df.to_csv(os.path.join(OUT_DIR, \"task_matched_metric_learning.csv\"), index=False)\n\n print(\"\\n=== Task-matched cross-session EER (metric-learning embeddings) ===\")\n print(task_df.to_string(index=False))\n\n return task_df\n\n\n# ----------------------------------------------------------------------\n# STEP 7: Ablation comparison against Module 1's softmax-derived embeddings\n# ----------------------------------------------------------------------\ndef compare_against_module1(metric_task_df):\n if not os.path.exists(MODULE1_TASK_EER_PATH):\n print(f\"\\n[compare_against_module1] Module 1 comparison file not found at \"\n f\"{MODULE1_TASK_EER_PATH} -- skipping ablation comparison.\")\n return None\n\n module1_df = pd.read_csv(MODULE1_TASK_EER_PATH)[[\"task\", \"EER\", \"ROC_AUC\"]]\n module1_df = module1_df.rename(columns={\"EER\": \"EER_softmax\", \"ROC_AUC\": \"AUC_softmax\"})\n\n metric_df = metric_task_df[[\"task\", \"EER\", \"ROC_AUC\"]].rename(\n columns={\"EER\": \"EER_metric_learning\", \"ROC_AUC\": \"AUC_metric_learning\"}\n )\n\n merged = module1_df.merge(metric_df, on=\"task\", how=\"inner\")\n merged[\"EER_improvement\"] = merged[\"EER_softmax\"] - merged[\"EER_metric_learning\"]\n merged = merged.sort_values(\"task\").reset_index(drop=True)\n merged.to_csv(os.path.join(OUT_DIR, \"ablation_comparison.csv\"), index=False)\n\n print(\"\\n=== Ablation: softmax-derived vs metric-learning embeddings (per task) ===\")\n print(merged.to_string(index=False))\n print(f\"\\nMean EER (softmax-derived): {merged['EER_softmax'].mean():.4f}\")\n print(f\"Mean EER (metric-learning): {merged['EER_metric_learning'].mean():.4f}\")\n print(f\"Mean improvement: {merged['EER_improvement'].mean():.4f}\")\n\n return merged\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 3: Task-Aware Metric-Learning EEG Identity Representation ===\\n\")\n\n manifest = load_manifest()\n target_shape = determine_target_shape(manifest)\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"].reset_index(drop=True)\n calibration = manifest[manifest[\"split\"] == \"calibration\"].reset_index(drop=True)\n cohort = sorted(set(enrollment[\"subject_id\"].unique()) & set(calibration[\"subject_id\"].unique()))\n print(f\"[main] Cohort: {len(cohort)} subjects\")\n\n label_encoder = LabelEncoder()\n label_encoder.fit(cohort)\n\n all_tasks = sorted(manifest[\"task\"].unique())\n task_encoder = {t: i for i, t in enumerate(all_tasks)}\n n_tasks = len(all_tasks)\n\n enrollment_cohort = enrollment[enrollment[\"subject_id\"].isin(cohort)]\n train_df, val_df = train_test_split(\n enrollment_cohort, test_size=0.1, stratify=enrollment_cohort[\"subject_id\"], random_state=RANDOM_STATE\n )\n\n model = train_metric_model(train_df, val_df, target_shape, label_encoder, task_encoder, n_tasks)\n check_embedding_diversity(model, val_df, target_shape, task_encoder)\n\n metric_task_df = run_task_matched_verification(manifest, model, target_shape, task_encoder, cohort)\n compare_against_module1(metric_task_df)\n\n print(f\"\\n=== Module 3 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 2 variant: Top-5-Task, 3-Stage Cascade (exploratory, self-contained)\n\nTests whether narrowing the cascade to only the 5 most session-stable tasks\n(by Module 1's EER ranking), restructured as a proper 3-stage cascade,\noutperforms the earlier 8-task/2-stage design. Self-contained -- no\ndependency on module2_cascade.py being present in the same directory.\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n ./module1_outputs/task_matched_diagnostic.csv\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nTASK_EER_PATH = os.path.join(MODULE1_DIR, \"task_matched_diagnostic.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\n\nOUT_DIR = \"./module2_top5_outputs\"\nos.makedirs(OUT_DIR, exist_ok=True)\n\nTOP_N_TASKS = 5\nN_STAGES = 3\nTARGET_FAR = 0.05\nEPOCHS_PER_TRIAL = 3\nBATCH_SIZE = 32\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n task_eer = pd.read_csv(TASK_EER_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries, Task EER table: {len(task_eer)} tasks\")\n return manifest, task_eer, embedding_model\n\n\ndef select_top_n_tasks(task_eer, n=TOP_N_TASKS):\n top = task_eer.sort_values(\"EER\").head(n).reset_index(drop=True)\n print(f\"[select_top_n_tasks] Top {n} tasks by EER (Module 1, softmax-derived embeddings):\")\n print(top[[\"task\", \"EER\", \"ROC_AUC\"]].to_string(index=False))\n return top\n\n\ndef assign_stages_manual(top_tasks_df, n_stages=N_STAGES):\n sorted_tasks = top_tasks_df.sort_values(\"EER\").reset_index(drop=True)\n task_list = sorted_tasks[\"task\"].tolist()\n eer_list = sorted_tasks[\"EER\"].tolist()\n\n groups = np.array_split(np.arange(len(task_list)), n_stages)\n stage_tasks = {}\n\n print(f\"\\n[assign_stages_manual] Stage assignment (top-{len(task_list)} tasks, {n_stages} stages):\")\n for i, idx_group in enumerate(groups):\n if len(idx_group) == 0:\n continue\n stage_name = f\"stage_{i + 1}\"\n stage_tasks[stage_name] = [int(task_list[j]) for j in idx_group]\n eers = [eer_list[j] for j in idx_group]\n print(f\" {stage_name}: tasks {stage_tasks[stage_name]}, EER range [{min(eers):.3f}, {max(eers):.3f}]\")\n\n return stage_tasks\n\n\n# ----------------------------------------------------------------------\n# Embedding helpers (identical normalization to Module 1/2)\n# ----------------------------------------------------------------------\ndef load_spectrogram_array(cache_path, target_shape=None):\n arr = np.load(cache_path)\n\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n shape_counts = pd.Series(shapes).value_counts()\n return shape_counts.idxmax()\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=BATCH_SIZE):\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n return embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\n# ----------------------------------------------------------------------\n# Templates (enrollment) and probe trials (calibration)\n# ----------------------------------------------------------------------\ndef build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, split=\"enrollment\"):\n split_df = manifest[manifest[\"split\"] == split]\n stage_reps = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n reps = {}\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid]\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n reps[sid] = embeddings.mean(axis=0)\n stage_reps[stage] = reps\n print(f\"[build_stage_templates] {split}/{stage}: {len(reps)} subjects, \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_reps\n\n\ndef build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, split=\"calibration\",\n epochs_per_trial=EPOCHS_PER_TRIAL):\n split_df = manifest[manifest[\"split\"] == split]\n stage_trials = {}\n\n for stage, tasks in stage_tasks.items():\n stage_df = split_df[split_df[\"task\"].isin(tasks)]\n trials = {}\n for sid in stage_df[\"subject_id\"].unique():\n sub_df = stage_df[stage_df[\"subject_id\"] == sid].reset_index(drop=True)\n embeddings = extract_embeddings(embedding_model, sub_df, target_shape)\n n_chunks = max(1, len(embeddings) // epochs_per_trial)\n chunks = np.array_split(embeddings, n_chunks)\n trials[sid] = [c.mean(axis=0) for c in chunks]\n stage_trials[stage] = trials\n n_total_trials = sum(len(v) for v in trials.values())\n avg_trials = n_total_trials / max(len(trials), 1)\n print(f\"[build_stage_probe_trials] {split}/{stage}: {len(trials)} subjects, \"\n f\"{n_total_trials} total trials ({avg_trials:.1f} avg trials/subject), \"\n f\"from {len(stage_df)} epochs across tasks {tasks}\")\n\n return stage_trials\n\n\n# ----------------------------------------------------------------------\n# Pairwise scores, fused scores, threshold calibration, cascade\n# ----------------------------------------------------------------------\ndef compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order):\n rows = []\n\n for true_sid in cohort:\n n_trials = max(\n (len(calibration_trials.get(stage, {}).get(true_sid, [])) for stage in stage_order),\n default=0,\n )\n\n for trial_idx in range(n_trials):\n trial_id = f\"{true_sid}_trial{trial_idx}\"\n\n for stage in stage_order:\n subject_trials = calibration_trials.get(stage, {}).get(true_sid, [])\n if trial_idx >= len(subject_trials):\n continue\n\n probe_emb = subject_trials[trial_idx]\n\n for claimed_sid in cohort:\n template_emb = enrollment_templates.get(stage, {}).get(claimed_sid)\n if template_emb is None:\n continue\n\n score = cosine_similarity(probe_emb, template_emb)\n rows.append({\n \"trial_id\": trial_id,\n \"true_subject\": true_sid,\n \"claimed_subject\": claimed_sid,\n \"stage\": stage,\n \"is_genuine\": true_sid == claimed_sid,\n \"score\": score,\n })\n\n df = pd.DataFrame(rows)\n n_genuine_trials = df[df[\"is_genuine\"]][\"trial_id\"].nunique()\n n_impostor_pairs = df[~df[\"is_genuine\"]][[\"trial_id\", \"claimed_subject\"]].drop_duplicates().shape[0]\n print(f\"\\n[compute_pairwise_stage_scores] {len(df)} stage-level scores computed \"\n f\"({n_genuine_trials} genuine trials, {n_impostor_pairs} impostor trial-claim pairs)\")\n return df\n\n\ndef compute_k_stage_fused_scores(scores_df, stage_order):\n pivot = scores_df.pivot_table(\n index=[\"trial_id\", \"true_subject\", \"claimed_subject\", \"is_genuine\"], columns=\"stage\", values=\"score\"\n ).reset_index()\n\n for k in range(1, len(stage_order) + 1):\n cols = stage_order[:k]\n pivot[f\"fused_k{k}\"] = pivot[cols].mean(axis=1, skipna=True)\n\n return pivot\n\n\ndef calibrate_k_thresholds(fused_df, stage_order, target_far=TARGET_FAR, reject_percentile=5):\n accept_thresholds, reject_thresholds = {}, {}\n\n for k in range(1, len(stage_order) + 1):\n col = f\"fused_k{k}\"\n genuine = fused_df[fused_df[\"is_genuine\"]][col].dropna().values\n impostor = fused_df[~fused_df[\"is_genuine\"]][col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, thresh = roc_curve(y_true, y_score)\n\n idx_candidates = np.where(fpr <= target_far)[0]\n if len(idx_candidates) > 0:\n idx = idx_candidates[np.argmax(tpr[idx_candidates])]\n accept_thr = float(thresh[idx])\n achieved_far = fpr[idx]\n else:\n eer_idx = np.nanargmin(np.abs(fpr - (1 - tpr)))\n accept_thr = float(thresh[eer_idx])\n achieved_far = fpr[eer_idx]\n\n reject_thr_raw = float(np.percentile(genuine, reject_percentile))\n reject_thr = min(reject_thr_raw, accept_thr - 1e-4)\n was_clamped = reject_thr_raw >= accept_thr\n\n accept_thresholds[k] = accept_thr\n reject_thresholds[k] = reject_thr\n\n clamp_note = \" [CLAMPED]\" if was_clamped else \"\"\n print(f\"[calibrate_k_thresholds] k={k}: accept>={accept_thr:.4f} (achieved FAR={achieved_far:.3f}), \"\n f\"reject<={reject_thr:.4f} (raw genuine-5th-pct={reject_thr_raw:.4f}){clamp_note}\")\n\n return accept_thresholds, reject_thresholds\n\n\ndef run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds):\n results = []\n\n for _, row in fused_df.iterrows():\n decision, resolved_k, first_available_k = None, None, None\n\n for k in range(1, len(stage_order) + 1):\n score = row.get(f\"fused_k{k}\")\n if pd.isna(score) or k not in accept_thresholds:\n continue\n\n if first_available_k is None:\n first_available_k = k\n\n if score >= accept_thresholds[k]:\n decision, resolved_k = \"accept\", k\n break\n elif score <= reject_thresholds[k]:\n decision, resolved_k = \"reject\", k\n break\n\n if decision is None:\n last_k = max(accept_thresholds.keys())\n last_score = row.get(f\"fused_k{last_k}\")\n if pd.notna(last_score) and last_score >= accept_thresholds[last_k]:\n decision, resolved_k = \"accept\", last_k\n else:\n decision, resolved_k = \"reject\", last_k\n\n was_true_escalation = (\n resolved_k is not None and first_available_k is not None\n and resolved_k > first_available_k and first_available_k == 1\n )\n\n results.append({\n \"trial_id\": row[\"trial_id\"],\n \"true_subject\": row[\"true_subject\"],\n \"claimed_subject\": row[\"claimed_subject\"],\n \"is_genuine\": row[\"is_genuine\"],\n \"decision\": decision,\n \"resolved_k\": resolved_k,\n \"first_available_k\": first_available_k,\n \"true_escalation\": was_true_escalation,\n })\n\n return pd.DataFrame(results)\n\n\n# ----------------------------------------------------------------------\n# Evaluation\n# ----------------------------------------------------------------------\ndef evaluate_cascade(cascade_results):\n genuine = cascade_results[cascade_results[\"is_genuine\"]]\n impostor = cascade_results[~cascade_results[\"is_genuine\"]]\n\n far = (impostor[\"decision\"] == \"accept\").mean()\n frr = (genuine[\"decision\"] == \"reject\").mean()\n\n print(f\"\\n[evaluate_cascade] Cascade FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n\n n_missing_stage1 = (cascade_results[\"first_available_k\"] != 1).sum()\n n_true_escalations = cascade_results[\"true_escalation\"].sum()\n print(f\"[evaluate_cascade] Trials with no Stage-1 data at all: {n_missing_stage1}\")\n print(f\"[evaluate_cascade] TRUE ambiguity-driven escalations: {n_true_escalations}\")\n\n return {\"far\": float(far), \"frr\": float(frr), \"n_genuine\": len(genuine), \"n_impostor\": len(impostor)}\n\n\ndef evaluate_cascade_by_data_availability(cascade_results):\n matched = cascade_results[cascade_results[\"first_available_k\"] == 1]\n extension = cascade_results[cascade_results[\"first_available_k\"] != 1]\n\n def _far_frr(df):\n genuine = df[df[\"is_genuine\"]]\n impostor = df[~df[\"is_genuine\"]]\n far = (impostor[\"decision\"] == \"accept\").mean() if len(impostor) > 0 else float(\"nan\")\n frr = (genuine[\"decision\"] == \"reject\").mean() if len(genuine) > 0 else float(\"nan\")\n return far, frr, len(genuine), len(impostor)\n\n m_far, m_frr, m_ng, m_ni = _far_frr(matched)\n e_far, e_frr, e_ng, e_ni = _far_frr(extension)\n\n print(f\"\\n[evaluate_cascade_by_data_availability] Matched population: \"\n f\"FAR={m_far:.4f}, FRR={m_frr:.4f} (n_genuine={m_ng}, n_impostor={m_ni})\")\n print(f\"[evaluate_cascade_by_data_availability] Coverage-extension: \"\n f\"FAR={e_far:.4f}, FRR={e_frr:.4f} (n_genuine={e_ng}, n_impostor={e_ni})\")\n\n return {\"matched\": {\"far\": m_far, \"frr\": m_frr}, \"coverage_extension\": {\"far\": e_far, \"frr\": e_frr}}\n\n\ndef evaluate_stage1_only_baseline(fused_df, accept_thresholds):\n valid = fused_df.dropna(subset=[\"fused_k1\"])\n n_dropped = len(fused_df) - len(valid)\n if n_dropped > 0:\n print(f\"[evaluate_stage1_only_baseline] Dropped {n_dropped} trial(s) with no Stage-1 data\")\n\n genuine = valid[valid[\"is_genuine\"]]\n impostor = valid[~valid[\"is_genuine\"]]\n\n threshold = accept_thresholds[1]\n far = (impostor[\"fused_k1\"] >= threshold).mean()\n frr = (genuine[\"fused_k1\"] < threshold).mean()\n\n print(f\"\\n[evaluate_stage1_only_baseline] Stage-1-only FAR: {far:.4f}, FRR: {frr:.4f} \"\n f\"(n_genuine={len(genuine)}, n_impostor={len(impostor)})\")\n return {\"far\": float(far), \"frr\": float(frr), \"n_genuine\": len(genuine), \"n_impostor\": len(impostor)}\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 2 variant: Top-5-Task, 3-Stage Cascade ===\\n\")\n\n manifest, task_eer, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n top5_eer = select_top_n_tasks(task_eer, n=TOP_N_TASKS)\n stage_tasks = assign_stages_manual(top5_eer, n_stages=N_STAGES)\n stage_order = sorted(stage_tasks.keys())\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"\\n[main] Cohort: {len(cohort)} subjects\")\n\n enrollment_templates = build_stage_templates(manifest, embedding_model, target_shape, stage_tasks, \"enrollment\")\n calibration_trials = build_stage_probe_trials(manifest, embedding_model, target_shape, stage_tasks, \"calibration\")\n\n scores_df = compute_pairwise_stage_scores(enrollment_templates, calibration_trials, cohort, stage_order)\n scores_df.to_csv(os.path.join(OUT_DIR, \"pairwise_stage_scores.csv\"), index=False)\n\n fused_df = compute_k_stage_fused_scores(scores_df, stage_order)\n fused_df.to_csv(os.path.join(OUT_DIR, \"fused_scores.csv\"), index=False)\n\n accept_thresholds, reject_thresholds = calibrate_k_thresholds(fused_df, stage_order)\n\n cascade_results = run_cascade(fused_df, stage_order, accept_thresholds, reject_thresholds)\n cascade_results.to_csv(os.path.join(OUT_DIR, \"cascade_decisions.csv\"), index=False)\n\n cascade_metrics = evaluate_cascade(cascade_results)\n subgroup_metrics = evaluate_cascade_by_data_availability(cascade_results)\n baseline_metrics = evaluate_stage1_only_baseline(fused_df, accept_thresholds)\n\n comparison = pd.DataFrame([\n {\"protocol\": \"Stage-1-only (no cascade)\", **baseline_metrics},\n {\"protocol\": \"Top-5-task, 3-stage cascade\",\n **{k: v for k, v in cascade_metrics.items() if k in (\"far\", \"frr\", \"n_genuine\", \"n_impostor\")}},\n ])\n comparison.to_csv(os.path.join(OUT_DIR, \"cascade_vs_baseline.csv\"), index=False)\n\n print(\"\\n=== Comparison: No-cascade baseline vs Top-5-task 3-stage cascade ===\")\n print(comparison.to_string(index=False))\n\n with open(os.path.join(OUT_DIR, \"module2_top5_summary.json\"), \"w\") as f:\n json.dump({\n \"stage_tasks\": stage_tasks,\n \"accept_thresholds\": accept_thresholds,\n \"reject_thresholds\": reject_thresholds,\n \"cascade_metrics\": cascade_metrics,\n \"subgroup_metrics\": subgroup_metrics,\n \"baseline_metrics\": baseline_metrics,\n }, f, indent=2)\n\n print(f\"\\n=== Module 2 (top-5 variant) complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 4: Explainable and Responsible EEG Biometric Authentication Evaluation\nM3CV EEG Biometric Authentication Pipeline\n\nSCOPE NOTE: All 65 EEG channels are flattened into one signal before mel\nspectrogram conversion (Module 0, matching the base paper's approach).\nChannel identity is destroyed before the CNN sees the data -- so Grad-CAM\nhere attributes importance to FREQUENCY and TIME regions of the spectrogram,\nNOT individual electrodes. State this explicitly as a limitation in the\nthesis; do not describe results as \"channel-level\" explainability.\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n - Module 2 outputs: ./module2_outputs/fused_scores.csv (for subject-wise\n EER and calibration -- uses the 8-task/2-stage run)\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport scipy.io as sio\nimport tensorflow as tf\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMODULE2_DIR = \"./module2_outputs\"\n\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\nFUSED_SCORES_PATH = os.path.join(MODULE2_DIR, \"fused_scores.csv\")\n\nENROLLMENT_DIR = \"/kaggle/input/competitions/eeg-biometric-competition/Enrollment\"\nCALIBRATION_DIR = \"/kaggle/input/competitions/eeg-biometric-competition/Calibration\"\n\nOUT_DIR = \"./module4_outputs\"\nGRADCAM_DIR = os.path.join(OUT_DIR, \"gradcam_examples\")\nos.makedirs(OUT_DIR, exist_ok=True)\nos.makedirs(GRADCAM_DIR, exist_ok=True)\n\nN_CHANNELS = 65\nDOWNSAMPLE_SFREQ = 250\nORIG_SFREQ = 1000\nN_MELS = 40\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n return manifest, embedding_model\n\n\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Identical normalization to Modules 1-3.\"\"\"\n arr = np.load(cache_path)\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n return pd.Series(shapes).value_counts().idxmax()\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=32):\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n return embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n\n\ndef build_subject_templates(embeddings, subject_ids):\n templates = {}\n for sid in np.unique(subject_ids):\n mask = subject_ids == sid\n templates[sid] = embeddings[mask].mean(axis=0)\n return templates\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Grad-CAM tied to the verification decision (cosine similarity to\n# claimed template), not the classification proxy task\n# ----------------------------------------------------------------------\ndef get_last_conv_layer(model):\n for layer in reversed(model.layers):\n if isinstance(layer, tf.keras.layers.Conv2D):\n return layer\n raise ValueError(\"No Conv2D layer found in model\")\n\n\ndef build_gradcam_model(embedding_model):\n conv_layer = get_last_conv_layer(embedding_model)\n print(f\"[build_gradcam_model] Using conv layer: {conv_layer.name}, output shape {conv_layer.output.shape}\")\n return tf.keras.Model(inputs=embedding_model.input, outputs=[conv_layer.output, embedding_model.output])\n\n\ndef compute_gradcam(grad_model, spectrogram, template_embedding):\n \"\"\"Grad-CAM heatmap for the SCALAR target: cosine_similarity(embedding\n of this spectrogram, claimed template). This ties the attribution\n directly to the authentication decision, not to a classification\n logit. Returns a (freq, time) heatmap resized to the input spectrogram\n shape, plus the similarity score itself.\"\"\"\n\n spec_tensor = tf.convert_to_tensor(spectrogram[np.newaxis, ..., np.newaxis], dtype=tf.float32)\n template_tensor = tf.convert_to_tensor(template_embedding, dtype=tf.float32)\n\n with tf.GradientTape() as tape:\n conv_output, embedding_output = grad_model(spec_tensor)\n embedding_output = tf.squeeze(embedding_output, axis=0)\n similarity = tf.reduce_sum(embedding_output * template_tensor) / (\n tf.norm(embedding_output) * tf.norm(template_tensor) + 1e-8\n )\n\n grads = tape.gradient(similarity, conv_output)\n pooled_grads = tf.reduce_mean(grads, axis=(0, 1, 2))\n conv_output = conv_output[0]\n\n heatmap = tf.reduce_sum(conv_output * pooled_grads, axis=-1)\n heatmap = tf.maximum(heatmap, 0)\n heatmap = heatmap / (tf.reduce_max(heatmap) + 1e-8)\n\n heatmap_resized = tf.image.resize(\n heatmap[..., tf.newaxis], spectrogram.shape, method=\"bilinear\"\n ).numpy().squeeze()\n\n return heatmap_resized, float(similarity.numpy())\n\n\ndef save_gradcam_visualization(spectrogram, heatmap, similarity, title, save_path):\n fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n axes[0].imshow(spectrogram, aspect=\"auto\", origin=\"lower\", cmap=\"viridis\")\n axes[0].set_title(\"Input spectrogram (normalized)\")\n axes[1].imshow(heatmap, aspect=\"auto\", origin=\"lower\", cmap=\"jet\")\n axes[1].set_title(\"Grad-CAM (freq x time attribution)\")\n axes[2].imshow(spectrogram, aspect=\"auto\", origin=\"lower\", cmap=\"gray\")\n axes[2].imshow(heatmap, aspect=\"auto\", origin=\"lower\", cmap=\"jet\", alpha=0.5)\n axes[2].set_title(\"Overlay\")\n fig.suptitle(f\"{title} (similarity={similarity:.3f})\")\n plt.tight_layout()\n plt.savefig(save_path, dpi=120, bbox_inches=\"tight\")\n plt.close()\n\n\ndef generate_gradcam_examples(manifest, embedding_model, target_shape, cohort, n_examples=5):\n \"\"\"Generate and save Grad-CAM visualizations for a handful of genuine\n cross-session verification trials (Enrollment template vs Calibration\n probe, same subject).\"\"\"\n\n grad_model = build_gradcam_model(embedding_model)\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n example_subjects = np.random.RandomState(RANDOM_STATE).choice(cohort, size=min(n_examples, len(cohort)), replace=False)\n\n for sid in example_subjects:\n enrollment_sub = enrollment[enrollment[\"subject_id\"] == sid]\n calibration_sub = calibration[calibration[\"subject_id\"] == sid]\n if len(enrollment_sub) == 0 or len(calibration_sub) == 0:\n continue\n\n template_embeddings = extract_embeddings(embedding_model, enrollment_sub, target_shape)\n template = template_embeddings.mean(axis=0)\n\n probe_row = calibration_sub.sample(1, random_state=RANDOM_STATE).iloc[0]\n probe_spec = load_spectrogram_array(probe_row[\"cache_path\"], target_shape)\n\n heatmap, similarity = compute_gradcam(grad_model, probe_spec, template)\n save_gradcam_visualization(\n probe_spec, heatmap, similarity,\n title=f\"Subject {sid}, task {probe_row['task']}\",\n save_path=os.path.join(GRADCAM_DIR, f\"gradcam_{sid}.png\"),\n )\n\n print(f\"[generate_gradcam_examples] Saved {len(example_subjects)} Grad-CAM visualizations to {GRADCAM_DIR}\")\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Cross-session explanation consistency\n# ----------------------------------------------------------------------\ndef evaluate_explanation_consistency(manifest, embedding_model, target_shape, cohort, n_subjects=10):\n \"\"\"For each subject, compute Grad-CAM on an Enrollment (S1) epoch and a\n Calibration (S2) epoch of the SAME task, both scored against that\n subject's own template. Compare heatmaps via 2D correlation -- a\n consistently-attending model should highlight similar time-frequency\n regions across sessions, even though the raw signal differs.\"\"\"\n\n grad_model = build_gradcam_model(embedding_model)\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n subjects_to_test = np.random.RandomState(RANDOM_STATE).choice(cohort, size=min(n_subjects, len(cohort)), replace=False)\n\n consistency_scores = []\n\n for sid in subjects_to_test:\n for task in common_tasks[:1]: # one task is enough to demonstrate the check; extend if needed\n e_rows = enrollment[(enrollment[\"subject_id\"] == sid) & (enrollment[\"task\"] == task)]\n c_rows = calibration[(calibration[\"subject_id\"] == sid) & (calibration[\"task\"] == task)]\n if len(e_rows) < 2 or len(c_rows) < 1:\n continue\n\n template_embeddings = extract_embeddings(embedding_model, e_rows, target_shape)\n template = template_embeddings.mean(axis=0)\n\n s1_spec = load_spectrogram_array(e_rows.iloc[0][\"cache_path\"], target_shape)\n s2_spec = load_spectrogram_array(c_rows.iloc[0][\"cache_path\"], target_shape)\n\n heatmap_s1, _ = compute_gradcam(grad_model, s1_spec, template)\n heatmap_s2, _ = compute_gradcam(grad_model, s2_spec, template)\n\n corr = np.corrcoef(heatmap_s1.flatten(), heatmap_s2.flatten())[0, 1]\n consistency_scores.append({\"subject_id\": sid, \"task\": int(task), \"heatmap_correlation\": float(corr)})\n\n consistency_df = pd.DataFrame(consistency_scores)\n consistency_df.to_csv(os.path.join(OUT_DIR, \"explanation_consistency.csv\"), index=False)\n\n print(f\"\\n[evaluate_explanation_consistency] Mean heatmap correlation across \"\n f\"{len(consistency_df)} subject/task pairs: {consistency_df['heatmap_correlation'].mean():.4f}\")\n print(consistency_df.to_string(index=False))\n\n return consistency_df\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Subject-wise EER breakdown\n# ----------------------------------------------------------------------\ndef evaluate_subject_wise_eer(fused_scores_path=FUSED_SCORES_PATH, score_col=\"fused_k1\"):\n \"\"\"For each subject X, treat X as a fixed claimed identity: gather all\n genuine attempts (X claiming X) and all impostor attempts (anyone else\n claiming X), compute EER for that subject's identity specifically.\n Surfaces whether any individual subject is especially vulnerable\n (low EER for impostors = easy to spoof) or especially fragile\n (high FRR = frequently locked out).\"\"\"\n\n if not os.path.exists(fused_scores_path):\n print(f\"[evaluate_subject_wise_eer] {fused_scores_path} not found -- skipping.\")\n return None\n\n fused_df = pd.read_csv(fused_scores_path)\n results = []\n\n for claimed_sid in sorted(fused_df[\"claimed_subject\"].unique()):\n subject_rows = fused_df[fused_df[\"claimed_subject\"] == claimed_sid]\n genuine = subject_rows[subject_rows[\"is_genuine\"]][score_col].dropna().values\n impostor = subject_rows[~subject_rows[\"is_genuine\"]][score_col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, _ = roc_curve(y_true, y_score)\n fnr = 1 - tpr\n eer_idx = np.nanargmin(np.abs(fpr - fnr))\n eer = (fpr[eer_idx] + fnr[eer_idx]) / 2\n\n results.append({\n \"subject_id\": claimed_sid, \"EER\": float(eer),\n \"n_genuine\": len(genuine), \"n_impostor\": len(impostor),\n })\n\n subject_eer_df = pd.DataFrame(results).sort_values(\"EER\", ascending=False).reset_index(drop=True)\n subject_eer_df.to_csv(os.path.join(OUT_DIR, \"subject_wise_eer.csv\"), index=False)\n\n print(\"\\n=== Subject-wise EER (highest to lowest -- top rows are the most vulnerable identities) ===\")\n print(subject_eer_df.to_string(index=False))\n print(f\"\\nMean subject-wise EER: {subject_eer_df['EER'].mean():.4f}, \"\n f\"std: {subject_eer_df['EER'].std():.4f}\")\n\n return subject_eer_df\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Confidence calibration (reliability diagram)\n# ----------------------------------------------------------------------\ndef evaluate_calibration(fused_scores_path=FUSED_SCORES_PATH, score_col=\"fused_k1\", n_bins=10):\n \"\"\"Bin scores into n_bins buckets and compare each bucket's mean score\n (treated as a proxy 'confidence') against the empirical fraction of\n genuine trials actually falling in that bucket. A well-calibrated\n score should show these two roughly track each other -- systematic\n divergence means the raw cosine similarity should NOT be read as a\n literal probability without further calibration (e.g. Platt scaling).\"\"\"\n\n if not os.path.exists(fused_scores_path):\n print(f\"[evaluate_calibration] {fused_scores_path} not found -- skipping.\")\n return None\n\n fused_df = pd.read_csv(fused_scores_path).dropna(subset=[score_col])\n fused_df[\"bin\"] = pd.cut(fused_df[score_col], bins=n_bins)\n\n calibration_table = fused_df.groupby(\"bin\", observed=True).agg(\n mean_score=(score_col, \"mean\"),\n fraction_genuine=(\"is_genuine\", \"mean\"),\n n_samples=(score_col, \"count\"),\n ).reset_index()\n\n calibration_table.to_csv(os.path.join(OUT_DIR, \"calibration_table.csv\"), index=False)\n\n fig, ax = plt.subplots(figsize=(6, 6))\n ax.plot(calibration_table[\"mean_score\"], calibration_table[\"fraction_genuine\"], \"o-\", label=\"Observed\")\n ax.plot([0, 1], [0, 1], \"k--\", alpha=0.3, label=\"Perfect calibration (reference)\")\n ax.set_xlabel(\"Mean score in bin (proxy confidence)\")\n ax.set_ylabel(\"Fraction actually genuine\")\n ax.set_title(\"Reliability diagram (Stage-1 fused score)\")\n ax.legend()\n plt.savefig(os.path.join(OUT_DIR, \"reliability_diagram.png\"), dpi=120, bbox_inches=\"tight\")\n plt.close()\n\n print(\"\\n=== Calibration table ===\")\n print(calibration_table.to_string(index=False))\n\n return calibration_table\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Robustness under missing channels / noise\n# ----------------------------------------------------------------------\ndef load_epoch_mat(epoch_id, split_dir):\n fname = f\"epoch{int(epoch_id):06d}.mat\" if str(epoch_id).isdigit() else f\"{epoch_id}.mat\"\n fpath = os.path.join(split_dir, fname)\n mat = sio.loadmat(fpath)\n data_keys = [k for k in mat.keys() if not k.startswith(\"__\")]\n return np.asarray(mat[data_keys[0]])\n\n\ndef corrupt_and_spectrogram(raw_array, target_shape, channel_dropout_frac=0.0, noise_std=0.0):\n \"\"\"Simulate corruption at the RAW epoch level (1000Hz, before\n downsampling/flattening), then run the SAME preprocessing order as\n Module 0 (downsample -> flatten -> mel spectrogram -> pad/crop ->\n normalize) to produce a spectrogram shaped identically to what the\n model was trained on.\n\n A PRIOR VERSION of this function skipped the downsample-to-250Hz step\n and flattened the raw 1000Hz signal directly -- 4x longer than what\n Module 0 actually caches -- producing a differently-shaped spectrogram\n that crashed the model's first Dense layer (shape mismatch: the Dense\n layer's weight matrix is fixed at training time to match target_shape).\n \"\"\"\n\n import librosa\n import mne\n\n corrupted = raw_array.copy().astype(np.float32)\n\n if channel_dropout_frac > 0:\n n_drop = int(N_CHANNELS * channel_dropout_frac)\n drop_idx = np.random.RandomState(RANDOM_STATE).choice(N_CHANNELS, size=n_drop, replace=False)\n corrupted[drop_idx, :] = 0.0\n\n if noise_std > 0:\n corrupted += np.random.RandomState(RANDOM_STATE).normal(0, noise_std, corrupted.shape)\n\n # Downsample 1000Hz -> 250Hz, matching Module 0's preprocess_epoch order\n channel_names = [f\"CH{i+1}\" for i in range(N_CHANNELS)]\n info = mne.create_info(ch_names=channel_names, sfreq=ORIG_SFREQ, ch_types=\"eeg\")\n raw_mne = mne.io.RawArray(corrupted, info, verbose=False)\n raw_mne.resample(DOWNSAMPLE_SFREQ, verbose=False)\n downsampled = raw_mne.get_data()\n\n flattened = downsampled.flatten()\n\n n_fft = min(1024, len(flattened))\n hop_length = max(n_fft // 8, 1)\n mel = librosa.feature.melspectrogram(y=flattened, sr=DOWNSAMPLE_SFREQ, n_fft=n_fft,\n hop_length=hop_length, n_mels=N_MELS)\n mel_db = librosa.power_to_db(mel, ref=np.max).astype(np.float32)\n\n # Defensive pad/crop to target_shape -- same safety net used everywhere\n # else in the pipeline (load_spectrogram_array), so any future shape\n # drift is handled explicitly instead of crashing at the model\n h, w = mel_db.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=mel_db.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = mel_db[:h_use, :w_use]\n mel_db = padded\n\n mel_db = (mel_db - mel_db.mean()) / (mel_db.std() + 1e-6)\n return mel_db\n\n\ndef evaluate_robustness(manifest, embedding_model, target_shape, cohort, n_probes=100):\n \"\"\"Build clean templates from Enrollment, then score a sample of\n Calibration probes under increasing channel-dropout and noise\n corruption, reporting EER degradation at each severity level.\"\"\"\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n template_embeddings = extract_embeddings(embedding_model, enrollment[enrollment[\"subject_id\"].isin(cohort)], target_shape)\n templates = build_subject_templates(template_embeddings, enrollment[enrollment[\"subject_id\"].isin(cohort)][\"subject_id\"].values)\n\n probe_pool = calibration[calibration[\"subject_id\"].isin(cohort)].sample(\n min(n_probes, len(calibration)), random_state=RANDOM_STATE\n )\n\n conditions = [\n (\"clean\", 0.0, 0.0),\n (\"channel_dropout_10pct\", 0.10, 0.0),\n (\"channel_dropout_25pct\", 0.25, 0.0),\n (\"channel_dropout_50pct\", 0.50, 0.0),\n (\"noise_low\", 0.0, 5.0),\n (\"noise_high\", 0.0, 20.0),\n ]\n\n degradation_results = []\n\n for cond_name, drop_frac, noise_std in conditions:\n genuine_scores, impostor_scores = [], []\n\n for _, row in probe_pool.iterrows():\n split_dir = ENROLLMENT_DIR if row[\"split\"] == \"enrollment\" else CALIBRATION_DIR\n try:\n raw_array = load_epoch_mat(row[\"epoch_id\"], split_dir)\n except Exception as e:\n continue\n\n spec = corrupt_and_spectrogram(raw_array, target_shape, channel_dropout_frac=drop_frac, noise_std=noise_std)\n spec_input = spec[np.newaxis, ..., np.newaxis].astype(np.float32)\n embedding = embedding_model.predict(spec_input, verbose=0)[0]\n\n true_sid = row[\"subject_id\"]\n for sid, template in templates.items():\n score = cosine_similarity(embedding, template)\n (genuine_scores if sid == true_sid else impostor_scores).append(score)\n\n if len(genuine_scores) == 0 or len(impostor_scores) == 0:\n continue\n\n y_true = np.concatenate([np.ones(len(genuine_scores)), np.zeros(len(impostor_scores))])\n y_score = np.concatenate([genuine_scores, impostor_scores])\n fpr, tpr, _ = roc_curve(y_true, y_score)\n fnr = 1 - tpr\n eer_idx = np.nanargmin(np.abs(fpr - fnr))\n eer = (fpr[eer_idx] + fnr[eer_idx]) / 2\n\n degradation_results.append({\"condition\": cond_name, \"EER\": float(eer)})\n print(f\"[evaluate_robustness] {cond_name}: EER={eer:.4f}\")\n\n degradation_df = pd.DataFrame(degradation_results)\n degradation_df.to_csv(os.path.join(OUT_DIR, \"robustness_degradation.csv\"), index=False)\n\n print(\"\\n=== Robustness: EER under corruption ===\")\n print(degradation_df.to_string(index=False))\n\n return degradation_df\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 4: Explainable and Responsible EEG Biometric Authentication Evaluation ===\\n\")\n\n manifest, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"[main] Cohort: {len(cohort)} subjects\\n\")\n\n generate_gradcam_examples(manifest, embedding_model, target_shape, cohort, n_examples=5)\n evaluate_explanation_consistency(manifest, embedding_model, target_shape, cohort, n_subjects=10)\n evaluate_subject_wise_eer()\n evaluate_calibration()\n evaluate_robustness(manifest, embedding_model, target_shape, cohort, n_probes=100)\n\n print(f\"\\n=== Module 4 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nModule 4: Explainable and Responsible EEG Biometric Authentication Evaluation\nM3CV EEG Biometric Authentication Pipeline\n\nSCOPE NOTE: All 65 EEG channels are flattened into one signal before mel\nspectrogram conversion (Module 0, matching the base paper's approach).\nChannel identity is destroyed before the CNN sees the data -- so Score-CAM\nhere attributes importance to FREQUENCY and TIME regions of the spectrogram,\nNOT individual electrodes. State this explicitly as a limitation in the\nthesis; do not describe results as \"channel-level\" explainability.\n\nDepends on:\n - Module 0 outputs: ./module0_outputs/spectrogram_manifest.csv + cache\n - Module 1 outputs: ./module1_outputs/embedding_model.keras\n - Module 2 outputs: ./module2_outputs/fused_scores.csv (for subject-wise\n EER and calibration -- uses the 8-task/2-stage run)\n\"\"\"\n\nimport os\nimport json\nimport numpy as np\nimport pandas as pd\nimport scipy.io as sio\nimport tensorflow as tf\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import roc_curve\n\n# ----------------------------------------------------------------------\n# CONFIG\n# ----------------------------------------------------------------------\nMODULE0_DIR = \"./module0_outputs\"\nMODULE1_DIR = \"./module1_outputs\"\nMODULE2_DIR = \"./module2_outputs\"\n\nMANIFEST_PATH = os.path.join(MODULE0_DIR, \"spectrogram_manifest.csv\")\nEMBEDDING_MODEL_PATH = os.path.join(MODULE1_DIR, \"embedding_model.keras\")\nFUSED_SCORES_PATH = os.path.join(MODULE2_DIR, \"fused_scores.csv\")\n\nENROLLMENT_DIR = \"/kaggle/input/competitions/eeg-biometric-competition/Enrollment\"\nCALIBRATION_DIR = \"/kaggle/input/competitions/eeg-biometric-competition/Calibration\"\n\nOUT_DIR = \"./module4_outputs\"\nSCORECAM_DIR = os.path.join(OUT_DIR, \"scorecam_examples\")\nos.makedirs(OUT_DIR, exist_ok=True)\nos.makedirs(SCORECAM_DIR, exist_ok=True)\n\nN_CHANNELS = 65\nDOWNSAMPLE_SFREQ = 250\nORIG_SFREQ = 1000\nN_MELS = 40\nRANDOM_STATE = 42\n\nnp.random.seed(RANDOM_STATE)\n\n\n# ----------------------------------------------------------------------\n# STEP 1: Load dependencies\n# ----------------------------------------------------------------------\ndef load_dependencies():\n manifest = pd.read_csv(MANIFEST_PATH)\n embedding_model = tf.keras.models.load_model(EMBEDDING_MODEL_PATH)\n print(f\"[load_dependencies] Manifest: {len(manifest)} entries\")\n return manifest, embedding_model\n\n\ndef load_spectrogram_array(cache_path, target_shape=None):\n \"\"\"Identical normalization to Modules 1-3.\"\"\"\n arr = np.load(cache_path)\n if target_shape is not None:\n h, w = arr.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=arr.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = arr[:h_use, :w_use]\n arr = padded\n arr = (arr - arr.mean()) / (arr.std() + 1e-6)\n return arr\n\n\ndef determine_target_shape(manifest, n_check=20):\n sample = manifest.sample(min(n_check, len(manifest)), random_state=RANDOM_STATE)\n shapes = [np.load(p).shape for p in sample[\"cache_path\"]]\n return pd.Series(shapes).value_counts().idxmax()\n\n\ndef cosine_similarity(a, b):\n return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-8)\n\n\ndef extract_embeddings(embedding_model, df, target_shape, batch_size=32):\n arrays = np.stack([\n load_spectrogram_array(p, target_shape) for p in df[\"cache_path\"]\n ])[..., np.newaxis].astype(np.float32)\n return embedding_model.predict(arrays, batch_size=batch_size, verbose=0)\n\n\ndef build_subject_templates(embeddings, subject_ids):\n templates = {}\n for sid in np.unique(subject_ids):\n mask = subject_ids == sid\n templates[sid] = embeddings[mask].mean(axis=0)\n return templates\n\n\n# ----------------------------------------------------------------------\n# STEP 2: Score-CAM tied to the verification decision (cosine similarity to\n# claimed template), not the classification proxy task\n#\n# Score-CAM (Wang et al., 2020) is used instead of Grad-CAM (Selvaraju et\n# al., 2017) because it is gradient-free: rather than backpropagating\n# through the network (which can be noisy or saturate in deeper layers --\n# this project already hit one dead-ReLU collapse earlier in training),\n# Score-CAM masks the input using each conv channel's own activation\n# pattern, runs a genuine forward pass, and uses the resulting change in\n# the model's real output score as that channel's importance weight. This\n# grounds the attribution in actual forward-pass evidence rather than\n# gradient magnitude, and is widely reported to produce sharper, more\n# faithful heatmaps than Grad-CAM.\n# ----------------------------------------------------------------------\ndef get_last_conv_layer(model):\n for layer in reversed(model.layers):\n if isinstance(layer, tf.keras.layers.Conv2D):\n return layer\n raise ValueError(\"No Conv2D layer found in model\")\n\n\ndef build_scorecam_conv_model(embedding_model):\n \"\"\"Dual-output model used ONCE per example to extract the last conv\n layer's activation maps for the clean (unmasked) input -- these maps\n are what get upsampled and used as masks in compute_scorecam.\"\"\"\n\n conv_layer = get_last_conv_layer(embedding_model)\n print(f\"[build_scorecam_conv_model] Using conv layer: {conv_layer.name}, output shape {conv_layer.output.shape}\")\n return tf.keras.Model(inputs=embedding_model.input, outputs=[conv_layer.output, embedding_model.output])\n\n\ndef compute_scorecam(embedding_model, conv_model, spectrogram, template_embedding, batch_size=32):\n \"\"\"Score-CAM heatmap for the SCALAR target: cosine_similarity(embedding\n of a masked spectrogram, claimed template). Ties attribution directly\n to the authentication decision, not a classification logit. Returns a\n (freq, time) heatmap resized to the input spectrogram shape, plus the\n CLEAN (unmasked) similarity score for reference.\"\"\"\n\n spec_tensor = tf.convert_to_tensor(spectrogram[np.newaxis, ..., np.newaxis], dtype=tf.float32)\n template_vec = np.asarray(template_embedding, dtype=np.float32)\n\n # 1. Get conv activation maps for the clean input (single forward pass)\n conv_output, _ = conv_model(spec_tensor)\n conv_output = conv_output[0].numpy() # (h, w, k)\n n_channels = conv_output.shape[-1]\n h_in, w_in = spectrogram.shape\n\n # 2. Build one masked input per conv channel: upsample that channel's\n # activation map to input resolution, normalize to [0,1], multiply\n # element-wise with the original spectrogram\n masked_inputs = []\n for k in range(n_channels):\n channel_map = conv_output[:, :, k]\n upsampled = tf.image.resize(\n channel_map[np.newaxis, ..., np.newaxis], (h_in, w_in), method=\"bilinear\"\n ).numpy().squeeze()\n\n cmin, cmax = upsampled.min(), upsampled.max()\n norm_map = np.zeros_like(upsampled) if (cmax - cmin) < 1e-8 else (upsampled - cmin) / (cmax - cmin)\n\n masked_inputs.append(spectrogram * norm_map)\n\n masked_inputs = np.stack(masked_inputs)[..., np.newaxis].astype(np.float32) # (k, h, w, 1)\n\n # 3. Forward pass ALL masked inputs through the full embedding model --\n # this is the gradient-free step: importance comes from how much\n # each mask actually shifts the real similarity score\n masked_embeddings = embedding_model.predict(masked_inputs, batch_size=batch_size, verbose=0) # (k, dim)\n\n norms = np.linalg.norm(masked_embeddings, axis=1) * np.linalg.norm(template_vec) + 1e-8\n channel_scores = (masked_embeddings @ template_vec) / norms # (k,)\n\n # 4. Softmax the scores into channel weights (standard Score-CAM step),\n # then take the weighted sum of the ORIGINAL (unmasked) activation\n # maps to build the final heatmap\n weights = np.exp(channel_scores - np.max(channel_scores))\n weights = weights / (weights.sum() + 1e-8)\n\n heatmap = np.zeros(conv_output.shape[:2], dtype=np.float32)\n for k in range(n_channels):\n heatmap += weights[k] * conv_output[:, :, k]\n\n heatmap = np.maximum(heatmap, 0)\n heatmap = heatmap / (heatmap.max() + 1e-8)\n\n heatmap_resized = tf.image.resize(\n heatmap[..., np.newaxis], spectrogram.shape, method=\"bilinear\"\n ).numpy().squeeze()\n\n clean_embedding = embedding_model.predict(spec_tensor, verbose=0)[0]\n clean_similarity = float(\n np.dot(clean_embedding, template_vec) /\n (np.linalg.norm(clean_embedding) * np.linalg.norm(template_vec) + 1e-8)\n )\n\n return heatmap_resized, clean_similarity\n\n\ndef save_xai_visualization(spectrogram, heatmap, similarity, title, save_path):\n fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n axes[0].imshow(spectrogram, aspect=\"auto\", origin=\"lower\", cmap=\"viridis\")\n axes[0].set_title(\"Input spectrogram (normalized)\")\n axes[1].imshow(heatmap, aspect=\"auto\", origin=\"lower\", cmap=\"jet\")\n axes[1].set_title(\"Score-CAM (freq x time attribution)\")\n axes[2].imshow(spectrogram, aspect=\"auto\", origin=\"lower\", cmap=\"gray\")\n axes[2].imshow(heatmap, aspect=\"auto\", origin=\"lower\", cmap=\"jet\", alpha=0.5)\n axes[2].set_title(\"Overlay\")\n fig.suptitle(f\"{title} (similarity={similarity:.3f})\")\n plt.tight_layout()\n plt.savefig(save_path, dpi=120, bbox_inches=\"tight\")\n plt.show()\n plt.close()\n\n\ndef generate_scorecam_examples(manifest, embedding_model, target_shape, cohort, n_examples=5):\n \"\"\"Generate, save, AND DISPLAY Score-CAM visualizations for a handful\n of genuine cross-session verification trials (Enrollment template vs\n Calibration probe, same subject).\"\"\"\n\n conv_model = build_scorecam_conv_model(embedding_model)\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n example_subjects = np.random.RandomState(RANDOM_STATE).choice(cohort, size=min(n_examples, len(cohort)), replace=False)\n\n for sid in example_subjects:\n enrollment_sub = enrollment[enrollment[\"subject_id\"] == sid]\n calibration_sub = calibration[calibration[\"subject_id\"] == sid]\n if len(enrollment_sub) == 0 or len(calibration_sub) == 0:\n continue\n\n template_embeddings = extract_embeddings(embedding_model, enrollment_sub, target_shape)\n template = template_embeddings.mean(axis=0)\n\n probe_row = calibration_sub.sample(1, random_state=RANDOM_STATE).iloc[0]\n probe_spec = load_spectrogram_array(probe_row[\"cache_path\"], target_shape)\n\n heatmap, similarity = compute_scorecam(embedding_model, conv_model, probe_spec, template)\n save_xai_visualization(\n probe_spec, heatmap, similarity,\n title=f\"Subject {sid}, task {probe_row['task']}\",\n save_path=os.path.join(SCORECAM_DIR, f\"scorecam_{sid}.png\"),\n )\n\n print(f\"[generate_scorecam_examples] Saved {len(example_subjects)} Score-CAM visualizations to {SCORECAM_DIR}\")\n\n\n# ----------------------------------------------------------------------\n# STEP 3: Cross-session explanation consistency\n# ----------------------------------------------------------------------\ndef evaluate_explanation_consistency(manifest, embedding_model, target_shape, cohort, n_subjects=10):\n \"\"\"For each subject, compute Score-CAM on an Enrollment (S1) epoch and a\n Calibration (S2) epoch of the SAME task, both scored against that\n subject's own template. Compare heatmaps via 2D correlation -- a\n consistently-attending model should highlight similar time-frequency\n regions across sessions, even though the raw signal differs.\"\"\"\n\n conv_model = build_scorecam_conv_model(embedding_model)\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n common_tasks = sorted(set(enrollment[\"task\"].unique()) & set(calibration[\"task\"].unique()))\n subjects_to_test = np.random.RandomState(RANDOM_STATE).choice(cohort, size=min(n_subjects, len(cohort)), replace=False)\n\n consistency_scores = []\n\n for sid in subjects_to_test:\n for task in common_tasks[:1]: # one task is enough to demonstrate the check; extend if needed\n e_rows = enrollment[(enrollment[\"subject_id\"] == sid) & (enrollment[\"task\"] == task)]\n c_rows = calibration[(calibration[\"subject_id\"] == sid) & (calibration[\"task\"] == task)]\n if len(e_rows) < 2 or len(c_rows) < 1:\n continue\n\n template_embeddings = extract_embeddings(embedding_model, e_rows, target_shape)\n template = template_embeddings.mean(axis=0)\n\n s1_spec = load_spectrogram_array(e_rows.iloc[0][\"cache_path\"], target_shape)\n s2_spec = load_spectrogram_array(c_rows.iloc[0][\"cache_path\"], target_shape)\n\n heatmap_s1, _ = compute_scorecam(embedding_model, conv_model, s1_spec, template)\n heatmap_s2, _ = compute_scorecam(embedding_model, conv_model, s2_spec, template)\n\n corr = np.corrcoef(heatmap_s1.flatten(), heatmap_s2.flatten())[0, 1]\n consistency_scores.append({\"subject_id\": sid, \"task\": int(task), \"heatmap_correlation\": float(corr)})\n\n consistency_df = pd.DataFrame(consistency_scores)\n consistency_df.to_csv(os.path.join(OUT_DIR, \"explanation_consistency.csv\"), index=False)\n\n print(f\"\\n[evaluate_explanation_consistency] Mean heatmap correlation across \"\n f\"{len(consistency_df)} subject/task pairs: {consistency_df['heatmap_correlation'].mean():.4f}\")\n print(consistency_df.to_string(index=False))\n\n return consistency_df\n\n\n# ----------------------------------------------------------------------\n# STEP 4: Subject-wise EER breakdown\n# ----------------------------------------------------------------------\ndef evaluate_subject_wise_eer(fused_scores_path=FUSED_SCORES_PATH, score_col=\"fused_k1\"):\n \"\"\"For each subject X, treat X as a fixed claimed identity: gather all\n genuine attempts (X claiming X) and all impostor attempts (anyone else\n claiming X), compute EER for that subject's identity specifically.\n Surfaces whether any individual subject is especially vulnerable\n (low EER for impostors = easy to spoof) or especially fragile\n (high FRR = frequently locked out).\"\"\"\n\n if not os.path.exists(fused_scores_path):\n print(f\"[evaluate_subject_wise_eer] {fused_scores_path} not found -- skipping.\")\n return None\n\n fused_df = pd.read_csv(fused_scores_path)\n results = []\n\n for claimed_sid in sorted(fused_df[\"claimed_subject\"].unique()):\n subject_rows = fused_df[fused_df[\"claimed_subject\"] == claimed_sid]\n genuine = subject_rows[subject_rows[\"is_genuine\"]][score_col].dropna().values\n impostor = subject_rows[~subject_rows[\"is_genuine\"]][score_col].dropna().values\n\n if len(genuine) == 0 or len(impostor) == 0:\n continue\n\n y_true = np.concatenate([np.ones_like(genuine), np.zeros_like(impostor)])\n y_score = np.concatenate([genuine, impostor])\n fpr, tpr, _ = roc_curve(y_true, y_score)\n fnr = 1 - tpr\n eer_idx = np.nanargmin(np.abs(fpr - fnr))\n eer = (fpr[eer_idx] + fnr[eer_idx]) / 2\n\n results.append({\n \"subject_id\": claimed_sid, \"EER\": float(eer),\n \"n_genuine\": len(genuine), \"n_impostor\": len(impostor),\n })\n\n subject_eer_df = pd.DataFrame(results).sort_values(\"EER\", ascending=False).reset_index(drop=True)\n subject_eer_df.to_csv(os.path.join(OUT_DIR, \"subject_wise_eer.csv\"), index=False)\n\n print(\"\\n=== Subject-wise EER (highest to lowest -- top rows are the most vulnerable identities) ===\")\n print(subject_eer_df.to_string(index=False))\n print(f\"\\nMean subject-wise EER: {subject_eer_df['EER'].mean():.4f}, \"\n f\"std: {subject_eer_df['EER'].std():.4f}\")\n\n return subject_eer_df\n\n\n# ----------------------------------------------------------------------\n# STEP 5: Confidence calibration (reliability diagram)\n# ----------------------------------------------------------------------\ndef evaluate_calibration(fused_scores_path=FUSED_SCORES_PATH, score_col=\"fused_k1\", n_bins=10):\n \"\"\"Bin scores into n_bins buckets and compare each bucket's mean score\n (treated as a proxy 'confidence') against the empirical fraction of\n genuine trials actually falling in that bucket. A well-calibrated\n score should show these two roughly track each other -- systematic\n divergence means the raw cosine similarity should NOT be read as a\n literal probability without further calibration (e.g. Platt scaling).\"\"\"\n\n if not os.path.exists(fused_scores_path):\n print(f\"[evaluate_calibration] {fused_scores_path} not found -- skipping.\")\n return None\n\n fused_df = pd.read_csv(fused_scores_path).dropna(subset=[score_col])\n fused_df[\"bin\"] = pd.cut(fused_df[score_col], bins=n_bins)\n\n calibration_table = fused_df.groupby(\"bin\", observed=True).agg(\n mean_score=(score_col, \"mean\"),\n fraction_genuine=(\"is_genuine\", \"mean\"),\n n_samples=(score_col, \"count\"),\n ).reset_index()\n\n calibration_table.to_csv(os.path.join(OUT_DIR, \"calibration_table.csv\"), index=False)\n\n fig, ax = plt.subplots(figsize=(6, 6))\n ax.plot(calibration_table[\"mean_score\"], calibration_table[\"fraction_genuine\"], \"o-\", label=\"Observed\")\n ax.plot([0, 1], [0, 1], \"k--\", alpha=0.3, label=\"Perfect calibration (reference)\")\n ax.set_xlabel(\"Mean score in bin (proxy confidence)\")\n ax.set_ylabel(\"Fraction actually genuine\")\n ax.set_title(\"Reliability diagram (Stage-1 fused score)\")\n ax.legend()\n plt.savefig(os.path.join(OUT_DIR, \"reliability_diagram.png\"), dpi=120, bbox_inches=\"tight\")\n plt.close()\n\n print(\"\\n=== Calibration table ===\")\n print(calibration_table.to_string(index=False))\n\n return calibration_table\n\n\n# ----------------------------------------------------------------------\n# STEP 6: Robustness under missing channels / noise\n# ----------------------------------------------------------------------\ndef load_epoch_mat(epoch_id, split_dir):\n fname = f\"epoch{int(epoch_id):06d}.mat\" if str(epoch_id).isdigit() else f\"{epoch_id}.mat\"\n fpath = os.path.join(split_dir, fname)\n mat = sio.loadmat(fpath)\n data_keys = [k for k in mat.keys() if not k.startswith(\"__\")]\n return np.asarray(mat[data_keys[0]])\n\n\ndef corrupt_and_spectrogram(raw_array, target_shape, channel_dropout_frac=0.0, noise_std=0.0):\n \"\"\"Simulate corruption at the RAW epoch level (1000Hz, before\n downsampling/flattening), then run the SAME preprocessing order as\n Module 0 (downsample -> flatten -> mel spectrogram -> pad/crop ->\n normalize) to produce a spectrogram shaped identically to what the\n model was trained on.\n\n A PRIOR VERSION of this function skipped the downsample-to-250Hz step\n and flattened the raw 1000Hz signal directly -- 4x longer than what\n Module 0 actually caches -- producing a differently-shaped spectrogram\n that crashed the model's first Dense layer (shape mismatch: the Dense\n layer's weight matrix is fixed at training time to match target_shape).\n \"\"\"\n\n import librosa\n import mne\n\n corrupted = raw_array.copy().astype(np.float32)\n\n if channel_dropout_frac > 0:\n n_drop = int(N_CHANNELS * channel_dropout_frac)\n drop_idx = np.random.RandomState(RANDOM_STATE).choice(N_CHANNELS, size=n_drop, replace=False)\n corrupted[drop_idx, :] = 0.0\n\n if noise_std > 0:\n corrupted += np.random.RandomState(RANDOM_STATE).normal(0, noise_std, corrupted.shape)\n\n # Downsample 1000Hz -> 250Hz, matching Module 0's preprocess_epoch order\n channel_names = [f\"CH{i+1}\" for i in range(N_CHANNELS)]\n info = mne.create_info(ch_names=channel_names, sfreq=ORIG_SFREQ, ch_types=\"eeg\")\n raw_mne = mne.io.RawArray(corrupted, info, verbose=False)\n raw_mne.resample(DOWNSAMPLE_SFREQ, verbose=False)\n downsampled = raw_mne.get_data()\n\n flattened = downsampled.flatten()\n\n n_fft = min(1024, len(flattened))\n hop_length = max(n_fft // 8, 1)\n mel = librosa.feature.melspectrogram(y=flattened, sr=DOWNSAMPLE_SFREQ, n_fft=n_fft,\n hop_length=hop_length, n_mels=N_MELS)\n mel_db = librosa.power_to_db(mel, ref=np.max).astype(np.float32)\n\n # Defensive pad/crop to target_shape -- same safety net used everywhere\n # else in the pipeline (load_spectrogram_array), so any future shape\n # drift is handled explicitly instead of crashing at the model\n h, w = mel_db.shape\n th, tw = target_shape\n if h != th or w != tw:\n padded = np.zeros((th, tw), dtype=mel_db.dtype)\n h_use, w_use = min(h, th), min(w, tw)\n padded[:h_use, :w_use] = mel_db[:h_use, :w_use]\n mel_db = padded\n\n mel_db = (mel_db - mel_db.mean()) / (mel_db.std() + 1e-6)\n return mel_db\n\n\ndef evaluate_robustness(manifest, embedding_model, target_shape, cohort, n_probes=100):\n \"\"\"Build clean templates from Enrollment, then score a sample of\n Calibration probes under increasing channel-dropout and noise\n corruption, reporting EER degradation at each severity level.\"\"\"\n\n enrollment = manifest[manifest[\"split\"] == \"enrollment\"]\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n\n template_embeddings = extract_embeddings(embedding_model, enrollment[enrollment[\"subject_id\"].isin(cohort)], target_shape)\n templates = build_subject_templates(template_embeddings, enrollment[enrollment[\"subject_id\"].isin(cohort)][\"subject_id\"].values)\n\n probe_pool = calibration[calibration[\"subject_id\"].isin(cohort)].sample(\n min(n_probes, len(calibration)), random_state=RANDOM_STATE\n )\n\n conditions = [\n (\"clean\", 0.0, 0.0),\n (\"channel_dropout_10pct\", 0.10, 0.0),\n (\"channel_dropout_25pct\", 0.25, 0.0),\n (\"channel_dropout_50pct\", 0.50, 0.0),\n (\"noise_low\", 0.0, 5.0),\n (\"noise_high\", 0.0, 20.0),\n ]\n\n degradation_results = []\n\n for cond_name, drop_frac, noise_std in conditions:\n genuine_scores, impostor_scores = [], []\n\n for _, row in probe_pool.iterrows():\n split_dir = ENROLLMENT_DIR if row[\"split\"] == \"enrollment\" else CALIBRATION_DIR\n try:\n raw_array = load_epoch_mat(row[\"epoch_id\"], split_dir)\n except Exception as e:\n continue\n\n spec = corrupt_and_spectrogram(raw_array, target_shape, channel_dropout_frac=drop_frac, noise_std=noise_std)\n spec_input = spec[np.newaxis, ..., np.newaxis].astype(np.float32)\n embedding = embedding_model.predict(spec_input, verbose=0)[0]\n\n true_sid = row[\"subject_id\"]\n for sid, template in templates.items():\n score = cosine_similarity(embedding, template)\n (genuine_scores if sid == true_sid else impostor_scores).append(score)\n\n if len(genuine_scores) == 0 or len(impostor_scores) == 0:\n continue\n\n y_true = np.concatenate([np.ones(len(genuine_scores)), np.zeros(len(impostor_scores))])\n y_score = np.concatenate([genuine_scores, impostor_scores])\n fpr, tpr, _ = roc_curve(y_true, y_score)\n fnr = 1 - tpr\n eer_idx = np.nanargmin(np.abs(fpr - fnr))\n eer = (fpr[eer_idx] + fnr[eer_idx]) / 2\n\n degradation_results.append({\"condition\": cond_name, \"EER\": float(eer)})\n print(f\"[evaluate_robustness] {cond_name}: EER={eer:.4f}\")\n\n degradation_df = pd.DataFrame(degradation_results)\n degradation_df.to_csv(os.path.join(OUT_DIR, \"robustness_degradation.csv\"), index=False)\n\n print(\"\\n=== Robustness: EER under corruption ===\")\n print(degradation_df.to_string(index=False))\n\n return degradation_df\n\n\ndef debug_compare_cached_vs_fresh(manifest, target_shape, n_samples=5):\n \"\"\"Diagnostic: for a handful of calibration epochs, compare the CACHED\n spectrogram (used by every other module) against a FRESHLY regenerated\n one via corrupt_and_spectrogram with zero corruption -- these should be\n nearly identical if the fresh-regeneration path genuinely replicates\n Module 0's caching pipeline. Run this BEFORE trusting evaluate_robustness\n results; if clean EER is near-random (~0.5), this diagnostic should\n reveal why.\"\"\"\n\n calibration = manifest[manifest[\"split\"] == \"calibration\"]\n sample = calibration.sample(min(n_samples, len(calibration)), random_state=RANDOM_STATE)\n\n print(\"\\n[debug_compare_cached_vs_fresh] Comparing cached vs freshly-regenerated spectrograms:\")\n\n for _, row in sample.iterrows():\n cached = load_spectrogram_array(row[\"cache_path\"], target_shape)\n\n try:\n raw_array = load_epoch_mat(row[\"epoch_id\"], CALIBRATION_DIR)\n except Exception as e:\n print(f\" epoch {row['epoch_id']}: FAILED to load raw .mat -- {e}\")\n continue\n\n fresh = corrupt_and_spectrogram(raw_array, target_shape, channel_dropout_frac=0.0, noise_std=0.0)\n\n corr = np.corrcoef(cached.flatten(), fresh.flatten())[0, 1]\n max_diff = np.max(np.abs(cached - fresh))\n mean_diff = np.mean(np.abs(cached - fresh))\n\n print(f\" epoch {row['epoch_id']} (subject {row['subject_id']}, task {row['task']}): \"\n f\"correlation={corr:.4f}, max_abs_diff={max_diff:.3f}, mean_abs_diff={mean_diff:.3f}\")\n\n print(\"[debug_compare_cached_vs_fresh] Correlation should be well above 0.95 for these to be \"\n \"considered equivalent. If correlation is low, the fresh-regeneration path has a real \"\n \"mismatch against Module 0's actual caching pipeline.\")\n\n\n# ----------------------------------------------------------------------\n# MAIN\n# ----------------------------------------------------------------------\nif __name__ == \"__main__\":\n print(\"=== Module 4: Explainable and Responsible EEG Biometric Authentication Evaluation ===\\n\")\n\n manifest, embedding_model = load_dependencies()\n target_shape = determine_target_shape(manifest)\n\n cohort = sorted(\n set(manifest[manifest[\"split\"] == \"enrollment\"][\"subject_id\"].unique()) &\n set(manifest[manifest[\"split\"] == \"calibration\"][\"subject_id\"].unique())\n )\n print(f\"[main] Cohort: {len(cohort)} subjects\\n\")\n\n generate_scorecam_examples(manifest, embedding_model, target_shape, cohort, n_examples=5)\n evaluate_explanation_consistency(manifest, embedding_model, target_shape, cohort, n_subjects=10)\n evaluate_subject_wise_eer()\n evaluate_calibration()\n\n debug_compare_cached_vs_fresh(manifest, target_shape, n_samples=5)\n evaluate_robustness(manifest, embedding_model, target_shape, cohort, n_probes=100)\n\n print(f\"\\n=== Module 4 complete. Outputs saved to {OUT_DIR} ===\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}