{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":87955,"databundleVersionId":10054796,"sourceType":"competition"}],"dockerImageVersionId":30918,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n! pip install mne\n! pip install torch_xla\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-02-22T21:06:42.274165Z","iopub.execute_input":"2025-02-22T21:06:42.274472Z","iopub.status.idle":"2025-02-22T21:06:55.892363Z","shell.execute_reply.started":"2025-02-22T21:06:42.274448Z","shell.execute_reply":"2025-02-22T21:06:55.891401Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Setup\n## Utils.py for general Utility\n[https://pypi.org/project/mne/](http://)\n\n---\n\n### Event markers:\n\n- 'S' (Sleep stage) → Red ('r')\n\n- 'K' (K-complex, a waveform in sleep EEG) → Green ('g')\n\n\nA brief, high spike in brain waves during deep sleep, helping with memory and keeping you asleep.\n\nK-complexes happen during deep (NREM => Non-Rapid Eye Movement) sleep to protect sleep by blocking minor disturbances and helping with memory processing. \n\n- 'REM' (Rapid Eye Movement sleep) → Blue ('b')\n\n- 'Son' (Sleep onset, start of sleep) → Magenta ('m')\n\n- 'Soff' (Sleep offset, wake-up time) → Cyan ('c')\n\n- 'A' (Arousal, brief wakefulness during sleep) → Yellow ('y')\n\n- 'MS' (Microstate, short EEG pattern states) → Black ('k')\n\n\nTiny patterns in brain activity that last a fraction of a second and show how the brain processes information.\n\nMicrostates occur in both wakefulness and sleep, reflecting rapid shifts in brain activity and thought processes.\n\n---","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nfrom scipy.io import loadmat\nimport numpy as np\nimport pandas as pd\nimport mne\nimport matplotlib.pyplot as plt\nimport pprint\nfrom matplotlib.lines import Line2D\n\n\n\nmarkers = ['S', 'K', 'REM', 'Son', 'Soff', 'A', 'MS', ]\nmarker_colors = {\n    'S': 'r',\n    'K': 'g',\n    'REM': 'b',\n    'Son': 'm',\n    'Soff': 'c',\n    'A': 'y',\n    'MS': 'k',\n}\n\ndef load_eeg(data_dir:Path, subject: str) -> mne.io.Raw:\n    fname = f\"{subject}_eeg_raw.mat\"\n    fpath = data_dir / fname\n    if not fpath.exists():\n        raise FileNotFoundError(f\"File {fpath} does not exist\")\n    \n    mat = loadmat(fpath)\n\n    pprint.pprint(mat.keys(), depth=2)\n    \n    eeg_data = mat[\"EEG\"][0, 0]['data']   \n    \n    info = mne.create_info(ch_names=['EEG1'], sfreq=250, ch_types=['eeg'])\n    raw = mne.io.RawArray(eeg_data, info)\n    print(raw.info)\n\n    return raw\n\ndef load_labels(data_dir:Path, subject: str) -> pd.DataFrame:\n    fname = f\"{subject}_labels.csv\"\n    fpath = data_dir / fname\n    if not fpath.exists():\n        raise FileNotFoundError(f\"File {fpath} does not exist\")\n    \n    df = pd.read_csv(fpath)\n\n    # for each marker M, theres M0 and M1. Change marker names to just M for both\n    df['Marker'] = df['Marker'].str.replace('0', '').str.replace('1', '')\n\n    # # keep only rows in which Marker has the '1' suffix\n    # df = df[df['Marker'].str.endswith('1')]\n\n    # # remove the '1' suffix from the Marker column\n    # df['Marker'] = df['Marker'].str.replace('1', '')\n\n    return df\n\n\ndef get_event_epochs(raw: mne.io.Raw, events_labels: pd.DataFrame, event:str, tmin:int=0, tmax:int=30) -> mne.Epochs:\n    df = events_labels\n    event_times = df[df['Marker'] == event][\"Timestamp_samples\"].values\n    event_epochs = df[df['Marker'] == event][\"Epoch\"].values\n\n\n    # Create an events array for MNE\n    events = np.column_stack([event_times, np.zeros(len(event_times), dtype=int), event_epochs])\n    try:\n        epochs = mne.Epochs(raw, events=events, tmin=tmin, tmax=tmax, baseline=None, preload=True)\n    except ValueError as e:\n        print(f\"Error creating epochs for event {event}: {e}\")\n        return None\n\n    return epochs\n\n\n\ndef sample2min(sample: int, sfreq: int) -> float:\n    return sample / sfreq / 60\n\n\n\n# Define stage names for legend labels\nstage_names = {\n    'S': 'Sleep stage',\n    'K': 'K-complex',\n    'REM': 'Rapid Eye Movement sleep',\n    'Son': 'Sleep onset',\n    'Soff': 'Sleep offset',\n    'A': 'Arousal',\n    'MS': 'Microstate'\n}\n\n\ndef plot_eeg(raw: mne.io.Raw, subject: str, t0: float, t1: float, events_labels: pd.DataFrame=None):\n    time = np.arange(raw._data.shape[1]) / raw.info['sfreq'] / 60\n\n    f, ax = plt.subplots(figsize=(14, 5))\n    ax.plot(time, raw._data[0], alpha=0.5)\n\n    if events_labels is not None:\n        for marker in markers:\n            events_samples = events_labels[events_labels['Marker'] == marker][\"Timestamp_samples\"].values\n            events_time = sample2min(events_samples, raw.info['sfreq'])\n            ax.vlines(events_time, -400, 400, alpha=0.5, color=marker_colors[marker])\n        \n        # Create legend handles with the full stage names\n        unique_markers = events_labels[\"Marker\"].unique()\n        legend_handles = [\n            Line2D([0], [0], color=marker_colors[marker], lw=4, label=stage_names[marker])\n            for marker in markers if marker in unique_markers\n        ]\n        # Place the legend below the graph, centered and taking the full width.\n        ax.legend(handles=legend_handles, title=\"Event Markers\", \n                  loc='upper center', bbox_to_anchor=(0.5, -0.15), \n                  ncol=len(legend_handles))\n\n    ax.set(\n        xlabel='Time (min)',\n        ylabel='EEG Amplitude',\n        ylim=[-400, 400],\n        xlim=[t0, t1],\n        title=f'EEG data for subject {subject} between {t0:.2f} and {t1:.2f} minutes'\n    )\n    \n    # Adjust layout to provide space for the legend at the bottom\n    f.tight_layout(rect=[0, 0.1, 1, 1])\n    return f, ax\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:06:55.893945Z","iopub.execute_input":"2025-02-22T21:06:55.894870Z","iopub.status.idle":"2025-02-22T21:06:58.997668Z","shell.execute_reply.started":"2025-02-22T21:06:55.894839Z","shell.execute_reply":"2025-02-22T21:06:58.996905Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"INPUT_DATA_DIR_PATH = Path(\"/kaggle/input/bci-i-idun-eeg-analysis-challenge\")\nOUTPUT_DATA_DIR_PATH=Path(\"/kaggle/working/\")\nSUBJECT = \"S001\"\nSUBJECTS_TRAIN =[\"S001\",\"S002\",\"S003\"]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:11.124679Z","iopub.execute_input":"2025-02-22T21:07:11.125186Z","iopub.status.idle":"2025-02-22T21:07:11.129457Z","shell.execute_reply.started":"2025-02-22T21:07:11.125161Z","shell.execute_reply":"2025-02-22T21:07:11.128581Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# General Explanation of the Data and functions used & Data Preparation\n\n#### Key Features of mne.io.Raw\n\n\n* Stores continuous EEG/MEG data (before epoching or event extraction)\n* Supports loading data from various formats (e.g., .fif, .edf, .mat).\n* Allows filtering, resampling, and referencing.\n* Can be visualized using built-in plotting functions.\n","metadata":{}},{"cell_type":"code","source":"eeg = load_eeg(INPUT_DATA_DIR_PATH,SUBJECT)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T18:48:18.647493Z","iopub.execute_input":"2025-02-22T18:48:18.647818Z","iopub.status.idle":"2025-02-22T18:48:19.244272Z","shell.execute_reply.started":"2025-02-22T18:48:18.647792Z","shell.execute_reply":"2025-02-22T18:48:19.242398Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"- Sampling Frequency: The number of times per second the EEG signal is recorded\n\n- High-Pass filter: Removes low-frequency noise (e.g., drift, slow artifacts), here captures all signals (0.00 Hz)\n\n- Low-Pass filter: Removes high-frequency noise (e.g., muscle activity, power line interference), Only frequencies below 125 Hz are kept.\n\nSince the sampling rate is 250 Hz, the maximum frequency that can be analyzed (*Nyquist limit*) is 125 Hz.\n\n> ### **Nyquist Limit (Nyquist Theorem) Explained Simply**\n> \n> #### **What is the Nyquist Limit?**  \n> The **Nyquist limit** (or **Nyquist frequency**) is the **highest frequency that can be accurately captured when sampling a signal**. It is **half of the sampling rate**:\n> \n> \n> **Nyquist Frequency = Sampling Rate / 2**\n>\n> \n> #### **Why is it Important?**  \n> - If a signal contains frequencies **above the Nyquist limit**, those frequencies will be **misrepresented** (aliasing occurs).  \n> - To accurately capture a signal, the **sampling rate must be at least twice the highest frequency in the signal**.  \n> \n> #### **Example (EEG Data with 250 Hz Sampling Rate)**  \n> - **Sampling Rate = 250 Hz**  \n> - **Nyquist Limit = \\(250/2 = 125\\) Hz**  \n> - Any EEG signal above **125 Hz** will be **distorted or lost**.  \n> \n> #### **What Happens if We Don't Follow Nyquist? (Aliasing)**  \n> If we try to sample a signal with frequencies **above the Nyquist limit**, it will appear as a **false (lower) frequency** in the recorded data. This is called **aliasing** and can create misleading results.\n> \n> #### **How to Prevent Aliasing?**  \n> 1. **Use a low-pass filter** before sampling (e.g., **125 Hz filter for 250 Hz sampling rate**).  \n> 2. **Increase the sampling rate** if higher frequencies need to be captured (e.g., **1000 Hz for brain signals up to 500 Hz**).  \n","metadata":{}},{"cell_type":"code","source":"labels = load_labels(INPUT_DATA_DIR_PATH, SUBJECT)\n\nlabels.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T18:48:21.232075Z","iopub.execute_input":"2025-02-22T18:48:21.232487Z","iopub.status.idle":"2025-02-22T18:48:21.247384Z","shell.execute_reply.started":"2025-02-22T18:48:21.232455Z","shell.execute_reply":"2025-02-22T18:48:21.246397Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"epoch_eeg_data = []\nepoch_labels = []\n\nfor marker in markers:\n    epochs = get_event_epochs(eeg, labels, marker)\n    if epochs is None:\n        print(f\"!!! --- No epochs found for marker {marker}\")\n        continue\n\n    _data = epochs.get_data()[:, 0, :]\n    epoch_eeg_data.append(_data)\n    epoch_labels.append(np.full(len(_data), marker))\n\nassert len(epoch_eeg_data) == len(epoch_labels)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T18:48:23.011684Z","iopub.execute_input":"2025-02-22T18:48:23.012038Z","iopub.status.idle":"2025-02-22T18:48:23.213646Z","shell.execute_reply.started":"2025-02-22T18:48:23.012012Z","shell.execute_reply":"2025-02-22T18:48:23.212702Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X = np.concatenate(epoch_eeg_data, axis=0)\ny = np.concatenate(epoch_labels, axis=0)\n\nprint(X.shape, y.shape)\n\nnp.save(OUTPUT_DATA_DIR_PATH / f\"{SUBJECT}_X.npy\", X)\nnp.save(OUTPUT_DATA_DIR_PATH / f\"{SUBJECT}_y.npy\", y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T18:48:26.670106Z","iopub.execute_input":"2025-02-22T18:48:26.670483Z","iopub.status.idle":"2025-02-22T18:48:26.715657Z","shell.execute_reply.started":"2025-02-22T18:48:26.670453Z","shell.execute_reply":"2025-02-22T18:48:26.714742Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Creating the npy files and storing\nfor subject in SUBJECTS_TRAIN:\n    eeg = load_eeg(INPUT_DATA_DIR_PATH,subject)\n\n    labels = load_labels(INPUT_DATA_DIR_PATH, subject)\n\n    epoch_eeg_data = []\n    epoch_labels = []\n    \n    for marker in markers:\n        epochs = get_event_epochs(eeg, labels, marker)\n        if epochs is None:\n            print(f\"!!! --- No epochs found for marker {marker}\")\n            continue\n    \n        _data = epochs.get_data()[:, 0, :]\n        epoch_eeg_data.append(_data)\n        epoch_labels.append(np.full(len(_data), marker))\n    \n    assert len(epoch_eeg_data) == len(epoch_labels)\n\n    X = np.concatenate(epoch_eeg_data, axis=0)\n    y = np.concatenate(epoch_labels, axis=0)\n    \n    print(X.shape, y.shape)\n    \n    np.save(OUTPUT_DATA_DIR_PATH / f\"{subject}_X.npy\", X)\n    np.save(OUTPUT_DATA_DIR_PATH / f\"{subject}_y.npy\", y)\n\n    \n    \n    \n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:19.930725Z","iopub.execute_input":"2025-02-22T21:07:19.931122Z","iopub.status.idle":"2025-02-22T21:07:22.594085Z","shell.execute_reply.started":"2025-02-22T21:07:19.931092Z","shell.execute_reply":"2025-02-22T21:07:22.592817Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The Juicy part - Data Analysis \nWe will not work with the data prepared earlier and work right back fresh","metadata":{}},{"cell_type":"code","source":"INPUT_DATA_DIR_PATH = Path(\"/kaggle/input/bci-i-idun-eeg-analysis-challenge\")\nOUTPUT_DATA_DIR_PATH=Path(\"/kaggle/working/\")\nSUBJECT = \"S001\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T20:55:38.377374Z","iopub.execute_input":"2025-02-22T20:55:38.377832Z","iopub.status.idle":"2025-02-22T20:55:38.382484Z","shell.execute_reply.started":"2025-02-22T20:55:38.377789Z","shell.execute_reply":"2025-02-22T20:55:38.380955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport mne\nimport scipy\nimport torch\n\n# choose a subject [0-3 only, subject 4 is in the test set!]\nsubject = SUBJECT","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:32.722522Z","iopub.execute_input":"2025-02-22T21:07:32.722848Z","iopub.status.idle":"2025-02-22T21:07:37.549515Z","shell.execute_reply.started":"2025-02-22T21:07:32.722823Z","shell.execute_reply":"2025-02-22T21:07:37.548776Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# this function plots the EEG data. You can pass a start/end time in minutes to zoom in on a particular time window\nplot = plot_eeg(eeg, SUBJECT, 40, 45)[0]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T09:47:19.149747Z","iopub.execute_input":"2025-02-22T09:47:19.150098Z","iopub.status.idle":"2025-02-22T09:47:19.783791Z","shell.execute_reply.started":"2025-02-22T09:47:19.150073Z","shell.execute_reply":"2025-02-22T09:47:19.782657Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"labels = load_labels(INPUT_DATA_DIR_PATH, SUBJECT)\nlabels = labels[labels.Marker.isin(markers)]\nlabels.head(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:18:17.571850Z","iopub.execute_input":"2025-02-22T10:18:17.572221Z","iopub.status.idle":"2025-02-22T10:18:17.587692Z","shell.execute_reply.started":"2025-02-22T10:18:17.572187Z","shell.execute_reply":"2025-02-22T10:18:17.586607Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"As we can see the epochs are not equaly balanced out, which is important to understand the data imbalance.","metadata":{}},{"cell_type":"code","source":"print(\"\\nCounts of each marker:\")\nprint(labels.Marker.value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:18:04.732198Z","iopub.execute_input":"2025-02-22T10:18:04.732555Z","iopub.status.idle":"2025-02-22T10:18:04.743162Z","shell.execute_reply.started":"2025-02-22T10:18:04.732527Z","shell.execute_reply":"2025-02-22T10:18:04.741694Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"f, ax = plt.subplots(figsize=(10, 5))\nax.bar(labels.Marker.value_counts().index, labels.Marker.value_counts().values)\nax.set(title=\"Marker counts\", xlabel=\"Marker\", ylabel=\"Count\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:19:19.887811Z","iopub.execute_input":"2025-02-22T10:19:19.888197Z","iopub.status.idle":"2025-02-22T10:19:20.093894Z","shell.execute_reply.started":"2025-02-22T10:19:19.888163Z","shell.execute_reply":"2025-02-22T10:19:20.092864Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Later, when we train models we will see the effect of this inbalance. For now, we can overlay the labels timestamps over the EEG data to get an idea of the timing of each event and what the signal looks like around each event.","metadata":{}},{"cell_type":"code","source":"plot = plot_eeg(eeg, subject, 0, 50, labels)[0]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:27:57.414904Z","iopub.execute_input":"2025-02-22T10:27:57.415420Z","iopub.status.idle":"2025-02-22T10:27:58.800945Z","shell.execute_reply.started":"2025-02-22T10:27:57.415380Z","shell.execute_reply":"2025-02-22T10:27:58.799658Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n- 'S' (Sleep stage) → Red ('r')\n\n- 'K' (K-complex, a waveform in sleep EEG) → Green ('g')\n\n\nA brief, high spike in brain waves during deep sleep, helping with memory and keeping you asleep.\n\nK-complexes happen during deep (NREM => Non-Rapid Eye Movement) sleep to protect sleep by blocking minor disturbances and helping with memory processing. \n\n- 'REM' (Rapid Eye Movement sleep) → Blue ('b')\n\n- 'Son' (Sleep onset, start of sleep) → Magenta ('m')\n\n- 'Soff' (Sleep offset, wake-up time) → Cyan ('c')\n\n- 'A' (Arousal, brief wakefulness during sleep) → Yellow ('y')\n\n- 'MS' (Microstate, short EEG pattern states) → Black ('k')","metadata":{}},{"cell_type":"markdown","source":"Finally, we can create \"epochs\" around each event and then collect all epochs for a particular event to get an idea for what the EEG signal looks like on avarage around each event. Note that epochs are taken to be 30 seconds long and starting at the event onset. You can play around with the tmin and tmax parameters to see how the epochs change for the function below to see how the epochs change, but keep in mind that pre-processed data and evaluation data assume this length and onset.\n\nhere, REM = Rapid Eye Movement Sleep","metadata":{}},{"cell_type":"code","source":"epochs = get_event_epochs(eeg, labels, 'REM')\nepochs.plot_image()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:30:17.908698Z","iopub.execute_input":"2025-02-22T10:30:17.909057Z","iopub.status.idle":"2025-02-22T10:30:25.205287Z","shell.execute_reply.started":"2025-02-22T10:30:17.909029Z","shell.execute_reply":"2025-02-22T10:30:25.204241Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Loading pre-processed data\n\n_X.npy and _Y.npy files contain the pre-processed data and labels, respectively. When loaded, they should be arrays of shape n_epochs x n_samples and n_epochs respectively. n_samples is 7501 (30 seconds of EEG at 250 Hz).The number of epochs changes across subjects, but the number of samples per epoch is fixed.","metadata":{}},{"cell_type":"code","source":"X = np.load(INPUT_DATA_DIR_PATH / f\"{SUBJECT}_X.npy\") # If we concat data for all subjects we will create a much better model\nY = np.load(INPUT_DATA_DIR_PATH / f\"{SUBJECT}_y.npy\")\n\nprint(X.shape, Y.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:32:29.635514Z","iopub.execute_input":"2025-02-22T10:32:29.635874Z","iopub.status.idle":"2025-02-22T10:32:29.668232Z","shell.execute_reply.started":"2025-02-22T10:32:29.635847Z","shell.execute_reply":"2025-02-22T10:32:29.667055Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Y contains marker associated with each epoch","metadata":{}},{"cell_type":"code","source":"Y[:5]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:32:42.176302Z","iopub.execute_input":"2025-02-22T10:32:42.176734Z","iopub.status.idle":"2025-02-22T10:32:42.184885Z","shell.execute_reply.started":"2025-02-22T10:32:42.176700Z","shell.execute_reply":"2025-02-22T10:32:42.183717Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"To get the epochs associated to a particular event, we can simply index the Y array. The following plot should look similar to the one above.","metadata":{}},{"cell_type":"code","source":"idxs = np.where(Y == 'REM')[0]\nrem_X = X[idxs]\n\nrem_average = rem_X.mean(axis=0)\n\nf, ax = plt.subplots(figsize=(10, 5))\nax.plot(rem_average)\nax.set(title=\"Average REM signal\", ylim=[-50, 50])\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:33:21.472061Z","iopub.execute_input":"2025-02-22T10:33:21.472439Z","iopub.status.idle":"2025-02-22T10:33:21.666125Z","shell.execute_reply.started":"2025-02-22T10:33:21.472411Z","shell.execute_reply":"2025-02-22T10:33:21.664902Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"As we can see, the number of epochs per event is different across events.\n","metadata":{}},{"cell_type":"code","source":"counts = {m: len(np.where(Y == m)[0]) for m in markers}\npercents = {m: round(counts[m]/len(Y) * 100, 2) for m in markers}\nsummary = pd.DataFrame({'Count': counts, 'Percent': percents})\nsummary","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:33:46.014436Z","iopub.execute_input":"2025-02-22T10:33:46.014838Z","iopub.status.idle":"2025-02-22T10:33:46.029013Z","shell.execute_reply.started":"2025-02-22T10:33:46.014803Z","shell.execute_reply":"2025-02-22T10:33:46.027602Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Fitting models\nIn this section we show how to fit a few simple models using scikit-learn and a more advanced deep learning model (time-series convolutional neural network) with PyTorch.\n\nFor simplicity, we will split the data from one subject into training and test sets. For the competition, however, you will have access to all the training data from three subjects and you will be asked to make predictions on data from a fourth subject for whom you do not have labels.\n\nThe code implementation will be similar, however, so you can use these examples to get started","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LogisticRegression \nfrom sklearn.metrics import balanced_accuracy_score, f1_score\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:40.507708Z","iopub.execute_input":"2025-02-22T21:07:40.508249Z","iopub.status.idle":"2025-02-22T21:07:41.300735Z","shell.execute_reply.started":"2025-02-22T21:07:40.508219Z","shell.execute_reply":"2025-02-22T21:07:41.300004Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Train/test split\nWe need to \"stratify\" the split to ensure that the training and test sets have the same proportions of each event type as the whole dataset.","metadata":{}},{"cell_type":"code","source":"def print_class_proportions(Y, label:str):\n    counts = {m: len(np.where(Y == m)[0]) for m in markers}\n    percents = {m: round(counts[m]/len(Y) * 100, 2) for m in markers}\n    summary = pd.DataFrame({'Count': counts, 'Percent': percents})\n    # sort by count\n    summary = summary.sort_values(by='Count', ascending=False)\n    print(f\"\\n{label} dataset:\")\n    print(summary)\n\n\n# counts for the whole dataset\nprint_class_proportions(Y, \"For Whole Subject 1 \")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T20:59:36.856986Z","iopub.execute_input":"2025-02-22T20:59:36.857446Z","iopub.status.idle":"2025-02-22T20:59:36.870949Z","shell.execute_reply.started":"2025-02-22T20:59:36.857406Z","shell.execute_reply":"2025-02-22T20:59:36.869760Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Example split without stratification.","metadata":{}},{"cell_type":"code","source":"# splitting the data without stratification\nX_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)\nprint_class_proportions(Y_train, \"S1 Training set\")\nprint_class_proportions(Y_test, \"S1 Test set\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:37:36.368749Z","iopub.execute_input":"2025-02-22T10:37:36.369086Z","iopub.status.idle":"2025-02-22T10:37:36.395428Z","shell.execute_reply.started":"2025-02-22T10:37:36.369060Z","shell.execute_reply":"2025-02-22T10:37:36.394337Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# the classess are not balanced, so we need to stratify the split\nX_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, stratify=Y)\nprint_class_proportions(Y_train, \"Training set\")\nprint_class_proportions(Y_test, \"Test set\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:39:24.648712Z","iopub.execute_input":"2025-02-22T10:39:24.649108Z","iopub.status.idle":"2025-02-22T10:39:24.679901Z","shell.execute_reply.started":"2025-02-22T10:39:24.649079Z","shell.execute_reply":"2025-02-22T10:39:24.678672Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"You should see much better results with the stratified splits. Our final pre-processing step is to standardize the EEG data.\n> As you can see its more evenly spread out for training data compared to without stratifying it. If we did not do that than the model might have overfitted due to the percentage of some classes being higher.\n> Now that we stratified, its more balanced out. It might be still imbalanced but its a step forward","metadata":{}},{"cell_type":"code","source":"# standardize the data\nscaler = StandardScaler()\nX_train = scaler.fit_transform(X_train)\nX_test = scaler.transform(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:42:08.718734Z","iopub.execute_input":"2025-02-22T10:42:08.719089Z","iopub.status.idle":"2025-02-22T10:42:08.807795Z","shell.execute_reply.started":"2025-02-22T10:42:08.719062Z","shell.execute_reply":"2025-02-22T10:42:08.806502Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Fitting models.\nWe define a couple of helper functions to fit models and plot the confusion matrix.\n\nOur main metric will be balanced F1 score which is ideal for imbalanced datasets, striking a balance between precision and recall. We also show balanced_accuracy_score. We also plot the confusion matrix to see how often each event type is predicted as another.","metadata":{}},{"cell_type":"code","source":"def plot_confusion_matrix(conf_mat, markers, title):\n    f, ax = plt.subplots(figsize=(10, 10))\n    sns.heatmap(conf_mat, annot=True, cmap='Blues', xticklabels=markers, yticklabels=markers, ax=ax)\n\n    # highlight the diagonal elements\n    for i in range(conf_mat.shape[0]):\n        ax.add_patch(plt.Rectangle((i, i), 1, 1, fill=False, edgecolor='red', lw=1))\n\n    ax.set(\n           title=title,\n           xlabel='Predicted',\n           ylabel='True',\n           aspect='equal'\n           )\n    plt.show()\n\n\ndef fit_and_eval_sk_classifier(clf, X_train, Y_train, X_test, Y_test):\n    clf.fit(X_train, Y_train)\n\n    Y_pred_train = clf.predict(X_train)\n    Y_pred_test = clf.predict(X_test)\n\n    if set(Y_pred_train) != set(Y_train):\n        print(f\"Predicted classes {set(Y_pred_train)} do not match true classes {set(Y_train)}\")\n    if set(Y_pred_test) != set(Y_test):\n        print(f\"Predicted classes {set(Y_pred_test)} do not match true classes {set(Y_test)}\")\n\n    train_acc = balanced_accuracy_score(Y_train, Y_pred_train)\n    test_acc = balanced_accuracy_score(Y_test, Y_pred_test)\n    print(f\"Train accuracy: {train_acc:.4f}, Test accuracy: {test_acc:.4f}\")\n\n    train_f1 = f1_score(Y_train, Y_pred_train, average='weighted')\n    test_f1 = f1_score(Y_test, Y_pred_test, average='weighted')\n    print(f\"Train F1: {train_f1:.4f}, Test F1: {test_f1:.4f}\")\n    \n    conf_mat = confusion_matrix(Y_test, Y_pred_test, normalize='true')\n    title=f\"{clf.__class__.__name__} - Train: {train_f1:.4f}, Test: {test_f1:.4f}\"\n    plot_confusion_matrix(conf_mat, markers, title)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:46.616365Z","iopub.execute_input":"2025-02-22T21:07:46.616831Z","iopub.status.idle":"2025-02-22T21:07:46.624767Z","shell.execute_reply.started":"2025-02-22T21:07:46.616786Z","shell.execute_reply":"2025-02-22T21:07:46.623889Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# fit multi-class logistic regression\nclf = LogisticRegression(multi_class='ovr', random_state=42, max_iter=1000, class_weight=None)\nfit_and_eval_sk_classifier(clf, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:44:13.762366Z","iopub.execute_input":"2025-02-22T10:44:13.762749Z","iopub.status.idle":"2025-02-22T10:44:21.387261Z","shell.execute_reply.started":"2025-02-22T10:44:13.762722Z","shell.execute_reply":"2025-02-22T10:44:21.386160Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"It seems to be predicting most data as Son. This is not surprising given that these are the most frequent event types.\n\n> as we saw the vertical lines of Son is more bluish hence it signifies model tends to predict that more as the other training data is less","metadata":{}},{"cell_type":"code","source":"print_class_proportions(Y, \"Whole dataset\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:49:54.854008Z","iopub.execute_input":"2025-02-22T10:49:54.854378Z","iopub.status.idle":"2025-02-22T10:49:54.864760Z","shell.execute_reply.started":"2025-02-22T10:49:54.854345Z","shell.execute_reply":"2025-02-22T10:49:54.863671Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We can use class_weight=\"balanced\" to give more weight to less frequent classes. This should improve the performance of the model on the minority classes.\n>  This is a mehtod of handling data imbalance","metadata":{}},{"cell_type":"code","source":"# fit multi-class logistic regression with class weights\nclf = LogisticRegression(multi_class='ovr', random_state=42, max_iter=1000, class_weight=\"balanced\")    \nfit_and_eval_sk_classifier(clf, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:50:39.927019Z","iopub.execute_input":"2025-02-22T10:50:39.927417Z","iopub.status.idle":"2025-02-22T10:50:48.427386Z","shell.execute_reply.started":"2025-02-22T10:50:39.927375Z","shell.execute_reply":"2025-02-22T10:50:48.426141Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Well for the code here: [https://github.com/BCI-I/BCII-IDUN-Challenge-tutorials/blob/main/02_tutorial.ipynb](https://github.com/BCI-I/BCII-IDUN-Challenge-tutorials/blob/main/02_tutorial.ipynb)\n\n## But here we see the f1 score decrease?! Interesting","metadata":{}},{"cell_type":"code","source":"from sklearn.svm import SVC\n\nsvc_cls = SVC(kernel='rbf', random_state=42, class_weight=\"balanced\")  # class_weight \"balanced\"\nfit_and_eval_sk_classifier(svc_cls, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:53:06.202255Z","iopub.execute_input":"2025-02-22T10:53:06.202722Z","iopub.status.idle":"2025-02-22T10:53:13.326841Z","shell.execute_reply.started":"2025-02-22T10:53:06.202691Z","shell.execute_reply":"2025-02-22T10:53:13.325404Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.svm import SVC\n\nsvc_cls = SVC(kernel='rbf', random_state=42)  # class_weight None\nfit_and_eval_sk_classifier(svc_cls, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:53:39.340674Z","iopub.execute_input":"2025-02-22T10:53:39.341070Z","iopub.status.idle":"2025-02-22T10:53:45.887989Z","shell.execute_reply.started":"2025-02-22T10:53:39.341039Z","shell.execute_reply":"2025-02-22T10:53:45.886723Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"You can also define your own custom weights based on the class frequencies. The example below is a bad set of weights, but you can play around with weights to see how to get the right balance across classes.","metadata":{}},{"cell_type":"code","source":"weights = {\n    'S' : 1,\n    'A': 1,\n    'K': 1,\n    'MS': 1,\n    'REM': 1,\n    'Son': 6,\n    'Soff': 6,\n}\n\nsvc_cls = SVC(kernel='rbf', random_state=42, class_weight=weights)  # pass  your own weights dictionary\nfit_and_eval_sk_classifier(svc_cls, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:57:18.087092Z","iopub.execute_input":"2025-02-22T10:57:18.087482Z","iopub.status.idle":"2025-02-22T10:57:24.897645Z","shell.execute_reply.started":"2025-02-22T10:57:18.087454Z","shell.execute_reply":"2025-02-22T10:57:24.896462Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's try Random Forest next.\n\n","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\n\nrf_cls = RandomForestClassifier(n_estimators=5, random_state=42, class_weight=\"balanced\")\nfit_and_eval_sk_classifier(rf_cls, X_train, Y_train, X_test, Y_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T10:58:07.080162Z","iopub.execute_input":"2025-02-22T10:58:07.080533Z","iopub.status.idle":"2025-02-22T10:58:08.062154Z","shell.execute_reply.started":"2025-02-22T10:58:07.080501Z","shell.execute_reply":"2025-02-22T10:58:08.060921Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"All of this model showed different levels of performance and different kinds of errors. The simple logistic regression was perhaps just too simple to achieve good performance. On the other hand, the SVM and Random Forest models seem to be overfitting the data, achieving high accuracy on the training set but low accuracy on the test set.","metadata":{}},{"cell_type":"markdown","source":"## A deep learning model\nDeep learning models are a natural choice for time-series classification tasks. We'll define a simple time-series convolutional neural network (tCNN) and train it using PyTorch.\n\nWe'll build on the lessons learned from the previous models to improve the performance of the deep learning model. We want something that is powerful enough to learn the complex relationships between EEG signals and the different event types, but we'll be careful to avoid overfitting. We'll apply three standard techniques to avoid overfitting:\n\nAugmentation: we'll create perturbed versions of the data to make the model more robust.\nRegularization: we'll add a penalty to the loss function to discourage complex models.\nEarly stopping: we'll do a short training run (you might want to use the validation loss to stop training more precisely)\nThe following block defines the model architecture.","metadata":{}},{"cell_type":"code","source":"#  define a time convolutional feed forward neural network with pytoimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import DataLoader, TensorDataset\nimport numpy as np\nimport torch\n\nclass TimeSeriesCNN(nn.Module):\n    def __init__(self, input_length, num_classes=7):\n        super(TimeSeriesCNN, self).__init__()\n        \n        # Convolutional layers\n        self.conv_layers = nn.Sequential(\n            nn.Conv1d(1, 32, kernel_size=8, stride=1, padding='same'),\n            nn.ReLU(),\n            nn.Dropout(0.2),\n            nn.MaxPool1d(kernel_size=2, stride=2),\n            \n            nn.Conv1d(32, 64, kernel_size=5, stride=1, padding='same'),\n            nn.ReLU(),\n            nn.Dropout(0.2),\n            nn.MaxPool1d(kernel_size=2, stride=2),\n            \n            nn.Conv1d(64, 128, kernel_size=3, stride=1, padding='same'),\n            nn.ReLU(),\n            nn.Dropout(0.2),\n            nn.MaxPool1d(kernel_size=2, stride=2)\n        )\n        \n        # Calculate the size of flattened features\n        self.flatten_size = self._get_flatten_size(input_length)\n        \n        # Fully connected layers\n        self.fc_layers = nn.Sequential(\n            nn.Linear(self.flatten_size, 256),\n            nn.ReLU(),\n            nn.Dropout(0.5),\n            nn.Linear(256, num_classes)\n        )\n        \n    def _get_flatten_size(self, input_length):\n        # Helper function to calculate flattened size\n        x = torch.randn(1, 1, input_length)\n        x = self.conv_layers(x)\n        return x.shape[1] * x.shape[2]\n    \n    def forward(self, x):\n        # x shape: (batch_size, n_samples)\n        # Add channel dimension\n        x = x.unsqueeze(1)  # (batch_size, 1, n_samples)\n        \n        x = self.conv_layers(x)\n        x = x.view(x.size(0), -1)  # Flatten\n        x = self.fc_layers(x)\n        return x","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:55.325720Z","iopub.execute_input":"2025-02-22T21:07:55.326087Z","iopub.status.idle":"2025-02-22T21:07:55.333743Z","shell.execute_reply.started":"2025-02-22T21:07:55.326058Z","shell.execute_reply":"2025-02-22T21:07:55.333058Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Next, we define a training function that will handle the augmentation, regularization, and early stopping for us.\n\n","metadata":{}},{"cell_type":"code","source":"# Training function\ndef train_model(model, train_loader, val_loader, num_epochs=50, device='cuda', lr=0.001, class_weights=None):\n    criterion = nn.CrossEntropyLoss(weight=class_weights) \n    optimizer = optim.Adam(model.parameters(), lr=lr, weight_decay=1e-5)  # weight_decay is L2 regularization\n        \n    model = model.to(device)\n    best_val_loss = np.inf\n    best_model = None\n    \n    for epoch in range(num_epochs):\n        # Training phase\n        model.train()\n        train_loss = 0\n        correct = 0\n        total = 0\n        \n        for inputs, labels in train_loader:\n            inputs, labels = inputs.to(device), labels.to(device)\n            \n            optimizer.zero_grad()\n            outputs = model(inputs)\n            loss = criterion(outputs, labels)\n            loss.backward()\n            optimizer.step()\n            \n            train_loss += loss.item()\n            _, predicted = outputs.max(1)\n            total += labels.size(0)\n            correct += predicted.eq(labels).sum().item()\n        \n        train_acc = 100. * correct / total\n        \n        # Validation phase\n        model.eval()\n        val_loss = 0\n        correct = 0\n        total = 0\n        \n        with torch.no_grad():\n            for inputs, labels in val_loader:\n                inputs, labels = inputs.to(device), labels.to(device)\n                outputs = model(inputs)\n                loss = criterion(outputs, labels)\n                \n                val_loss += loss.item()\n                _, predicted = outputs.max(1)\n                total += labels.size(0)\n                correct += predicted.eq(labels).sum().item()\n\n                if val_loss < best_val_loss:\n                    best_val_loss = val_loss\n                    best_model = model.state_dict()\n        \n        val_acc = 100. * correct / total\n        \n        if epoch == 0 or (epoch+1) % 5 == 0:\n            print(f'Epoch [{epoch+1}/{num_epochs}]')\n            print(f'Train Loss: {train_loss/len(train_loader):.4f}, Train Acc: {train_acc:.2f}%')\n            print(f'Val Loss: {val_loss/len(val_loader):.4f}, Val Acc: {val_acc:.2f}%')\n            print('--------------------')\n\n    return best_model","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:07:59.538551Z","iopub.execute_input":"2025-02-22T21:07:59.538857Z","iopub.status.idle":"2025-02-22T21:07:59.547346Z","shell.execute_reply.started":"2025-02-22T21:07:59.538835Z","shell.execute_reply":"2025-02-22T21:07:59.546507Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Finally, we define some custom Datasets to handle the augmentation and data loading.\n\n","metadata":{}},{"cell_type":"code","source":"from torch.utils.data import Dataset\n\nclass TimeSeriesAugmentation:\n    def __init__(self, noise_level=0.05, shift_range=10):\n        self.noise_level = noise_level\n        self.shift_range = shift_range\n    \n    def add_noise(self, x):\n        noise = torch.randn_like(x) * self.noise_level\n        return x + noise\n    \n    def time_shift(self, x):\n        # Generate random shift value\n        shift = torch.randint(-self.shift_range, self.shift_range, (1,)).item()\n        # Roll the tensor along the time dimension\n        return torch.roll(x, shifts=shift, dims=0)\n    \n    def scale(self, x):\n        scale_factor = torch.rand(1) * 0.4 + 0.8  # Random scale between 0.8 and 1.2\n        return x * scale_factor\n\nclass AugmentedDataset(Dataset):\n    def __init__(self, X, y, augment=True):\n        self.X = X\n        self.y = y\n        self.augment = augment\n        self.aug = TimeSeriesAugmentation()\n    \n    def __getitem__(self, idx):\n        x = self.X[idx].clone()  # Create a copy to avoid modifying original data\n        if self.augment:\n            if torch.rand(1) < 0.5:\n                x = self.aug.add_noise(x)\n            if torch.rand(1) < 0.5:\n                x = self.aug.time_shift(x)\n            if torch.rand(1) < 0.5:\n                x = self.aug.scale(x)\n        return x, self.y[idx]\n    \n    def __len__(self):\n        return len(self.X)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:03.807722Z","iopub.execute_input":"2025-02-22T21:08:03.808174Z","iopub.status.idle":"2025-02-22T21:08:03.815715Z","shell.execute_reply.started":"2025-02-22T21:08:03.808143Z","shell.execute_reply":"2025-02-22T21:08:03.815080Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Training the DL model\nWe need to convert the labels to integers and then to PyTorch tensors and turn the data into PyTorch tensors. Then we prepare dataloaders using our custom datasets","metadata":{}},{"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n\n\n# convert Y from text labels to integer labels\n_markers = list(set(Y))\n_Y_train = [_markers.index(m) for m in Y_train]\n_Y_test = [_markers.index(m) for m in Y_test]\n\n# Convert to PyTorch tensors\n_X_train = torch.FloatTensor(X_train)\n_Y_train = torch.LongTensor(_Y_train)\n_X_test = torch.FloatTensor(X_test)\n_Y_test = torch.LongTensor(_Y_test)\n\n# Create data loaders\ntrain_dataset = AugmentedDataset(_X_train, _Y_train)\nval_dataset = AugmentedDataset(_X_test, _Y_test)\ntrain_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)\nval_loader = DataLoader(val_dataset, batch_size=32, shuffle=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:06.439715Z","iopub.execute_input":"2025-02-22T21:08:06.440103Z","iopub.status.idle":"2025-02-22T21:08:06.558092Z","shell.execute_reply.started":"2025-02-22T21:08:06.440073Z","shell.execute_reply":"2025-02-22T21:08:06.557013Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"device","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T20:55:27.124084Z","iopub.execute_input":"2025-02-22T20:55:27.124502Z","iopub.status.idle":"2025-02-22T20:55:27.150061Z","shell.execute_reply.started":"2025-02-22T20:55:27.124473Z","shell.execute_reply":"2025-02-22T20:55:27.148496Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We're finally ready to train the model.\n\n","metadata":{}},{"cell_type":"code","source":"# Initialize model\ninput_length = X_train.shape[1]  # number of time samples\nmodel = TimeSeriesCNN(input_length=input_length, num_classes=len(set(Y)))\n\nmodel = model.to(device)\n\n# Train model\nbest_model = train_model(model, train_loader, val_loader, num_epochs=30, device=device, lr=1e-4)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T11:02:49.108433Z","iopub.execute_input":"2025-02-22T11:02:49.108808Z","iopub.status.idle":"2025-02-22T11:10:34.738373Z","shell.execute_reply.started":"2025-02-22T11:02:49.108780Z","shell.execute_reply":"2025-02-22T11:10:34.736325Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"To evaluate the model, we can load the best model and compute the balanced accuracy and F1 score. We also plot the confusion matrix to see how often each event type is predicted as another.","metadata":{}},{"cell_type":"code","source":"# load best model\nmodel.load_state_dict(best_model)\n\n# compute metric\nY_pred = model(_X_test.to(device)).argmax(1).cpu().numpy()\nscore = balanced_accuracy_score(_Y_test, Y_pred)\nf1 = f1_score(_Y_test, Y_pred, average='weighted')\n\n# plot confusion matrix\nconf_mat = confusion_matrix(_Y_test, Y_pred, normalize='true')\nplot_confusion_matrix(conf_mat, markers, f\"tCNN - Accuracy: {score:.4f}, F1: {f1:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T11:13:08.094851Z","iopub.execute_input":"2025-02-22T11:13:08.096371Z","iopub.status.idle":"2025-02-22T11:13:10.812635Z","shell.execute_reply.started":"2025-02-22T11:13:08.096323Z","shell.execute_reply":"2025-02-22T11:13:10.811330Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Findings (with data of just subject 1):\n\n- Son doesnt predict any Son, same with Soff; instead the predictions are opposite Son predicts Soff and vice versa\n- REM and A are the most predicted then comes MS\n- Accurate predictions: REM, MS, A, S ; None, None; which is surprising the data which we need to predict is sleepiness and we cant predict the Son and Soff ; thats the irony we cant predict Sleep onset and Sleep offset","metadata":{}},{"cell_type":"markdown","source":"Similarly to what we did for the previous models, we can balance the classes by computing the class weights based on the frequency of each class in the training set. Here we compute the weights of each class in the training set and pass them to the CrossEntropyLoss function.","metadata":{}},{"cell_type":"code","source":"# Calculate class weights based on frequency\ndef compute_class_weights(y_train):\n    class_counts = torch.bincount(y_train)\n    total = len(y_train)\n    weights = total / (len(class_counts) * class_counts.float())\n    return weights\n\n\nweights = compute_class_weights(_Y_train).to(device)\nbalanced_model = train_model(model, train_loader, val_loader, num_epochs=30, device=device, lr=1e-4, class_weights=weights)\n\n# load best model\nmodel.load_state_dict(best_model)\n\n# compute metric\nY_pred = model(_X_test.to(device)).argmax(1).cpu().numpy()\nscore = balanced_accuracy_score(_Y_test, Y_pred)\nf1 = f1_score(_Y_test, Y_pred, average='weighted')\n\n# plot confusion matrix\nconf_mat = confusion_matrix(_Y_test, Y_pred, normalize='true')\nplot_confusion_matrix(conf_mat, markers, f\"tCNN - Accuracy: {score:.4f}, F1: {f1:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T11:13:11.138397Z","iopub.execute_input":"2025-02-22T11:13:11.138921Z","iopub.status.idle":"2025-02-22T11:21:23.385839Z","shell.execute_reply.started":"2025-02-22T11:13:11.138874Z","shell.execute_reply":"2025-02-22T11:21:23.384613Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Conclusions\nThe models shown here are intended simply as a demonstration for how you can load and process data before fitting and evaluating your own models. There's much more that can be done to improve performance! This is where you come in, experiment with different ideas and come up with a great machine learning model.","metadata":{}},{"cell_type":"markdown","source":"# Now with full data\n\n","metadata":{}},{"cell_type":"markdown","source":"Run directly from here","metadata":{}},{"cell_type":"code","source":"#imports\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# CONCAT 3 FILES\n# List of subject IDs (adjust these IDs as needed)\nsubjects = [\"S001\", \"S002\", \"S003\"]\n\n# Initialize empty lists to store data for each subject\nX_list = []\nY_list = []\n\n# Loop over each subject and load their data\nfor subject in subjects:\n    # /kaggle/working/S001_X.npy\n    print()\n    X_subject = np.load(OUTPUT_DATA_DIR_PATH / f\"{subject}_X.npy\")\n    Y_subject = np.load(OUTPUT_DATA_DIR_PATH / f\"{subject}_y.npy\")\n    X_list.append(X_subject)\n    Y_list.append(Y_subject)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:17.682577Z","iopub.execute_input":"2025-02-22T21:08:17.682895Z","iopub.status.idle":"2025-02-22T21:08:17.742967Z","shell.execute_reply.started":"2025-02-22T21:08:17.682873Z","shell.execute_reply":"2025-02-22T21:08:17.742116Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# Concatenate the data from all subjects along the first axis\nX = np.concatenate(X_list, axis=0)\nY = np.concatenate(Y_list, axis=0)\n\nprint(X.shape, Y.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:19.040241Z","iopub.execute_input":"2025-02-22T21:08:19.040529Z","iopub.status.idle":"2025-02-22T21:08:19.099483Z","shell.execute_reply.started":"2025-02-22T21:08:19.040507Z","shell.execute_reply":"2025-02-22T21:08:19.098740Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def print_class_proportions(Y, label:str):\n    counts = {m: len(np.where(Y == m)[0]) for m in markers}\n    percents = {m: round(counts[m]/len(Y) * 100, 2) for m in markers}\n    summary = pd.DataFrame({'Count': counts, 'Percent': percents})\n    # sort by count\n    summary = summary.sort_values(by='Count', ascending=False)\n    print(f\"\\n{label} dataset:\")\n    print(summary)\n\n\n# counts for the whole dataset\nprint_class_proportions(Y, \"For Whole \")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:53.001702Z","iopub.execute_input":"2025-02-22T21:08:53.002032Z","iopub.status.idle":"2025-02-22T21:08:53.011854Z","shell.execute_reply.started":"2025-02-22T21:08:53.001971Z","shell.execute_reply":"2025-02-22T21:08:53.011078Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# splitting the data without stratification\nX_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)\nprint_class_proportions(Y_train, \"Full Training set\")\nprint_class_proportions(Y_test, \"Full Test set\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:08:55.050358Z","iopub.execute_input":"2025-02-22T21:08:55.050673Z","iopub.status.idle":"2025-02-22T21:08:55.112700Z","shell.execute_reply.started":"2025-02-22T21:08:55.050647Z","shell.execute_reply":"2025-02-22T21:08:55.111894Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# TPU\n# import torch\n# import torch_xla\n# import torch_xla.core.xla_model as xm\n\n# # Get the TPU device\n# device = xm.xla_device()\n\n# # Check if connected\n# print(f\"Using device: {device}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:06:21.911142Z","iopub.execute_input":"2025-02-22T21:06:21.911372Z","iopub.status.idle":"2025-02-22T21:06:21.915087Z","shell.execute_reply.started":"2025-02-22T21:06:21.911349Z","shell.execute_reply":"2025-02-22T21:06:21.914270Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n\n\n# convert Y from text labels to integer labels\n_markers = list(set(Y))\n_Y_train = [_markers.index(m) for m in Y_train]\n_Y_test = [_markers.index(m) for m in Y_test]\n\n# Convert to PyTorch tensors\n_X_train = torch.FloatTensor(X_train)\n_Y_train = torch.LongTensor(_Y_train)\n_X_test = torch.FloatTensor(X_test)\n_Y_test = torch.LongTensor(_Y_test)\n\n# Create data loaders\ntrain_dataset = AugmentedDataset(_X_train, _Y_train)\nval_dataset = AugmentedDataset(_X_test, _Y_test)\ntrain_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)\nval_loader = DataLoader(val_dataset, batch_size=32, shuffle=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:09:33.651870Z","iopub.execute_input":"2025-02-22T21:09:33.652247Z","iopub.status.idle":"2025-02-22T21:09:33.741561Z","shell.execute_reply.started":"2025-02-22T21:09:33.652221Z","shell.execute_reply":"2025-02-22T21:09:33.740565Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"device","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:09:36.220389Z","iopub.execute_input":"2025-02-22T21:09:36.220740Z","iopub.status.idle":"2025-02-22T21:09:36.226963Z","shell.execute_reply.started":"2025-02-22T21:09:36.220709Z","shell.execute_reply":"2025-02-22T21:09:36.226098Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Initialize model\ninput_length = X_train.shape[1]  # number of time samples\nmodel = TimeSeriesCNN(input_length=input_length, num_classes=len(set(Y)))\n\nmodel = model.to(device)\n\n# Train model\nbest_model = train_model(model, train_loader, val_loader, num_epochs=30, device=device, lr=1e-4)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:09:38.604417Z","iopub.execute_input":"2025-02-22T21:09:38.604696Z","iopub.status.idle":"2025-02-22T21:10:49.760113Z","shell.execute_reply.started":"2025-02-22T21:09:38.604675Z","shell.execute_reply":"2025-02-22T21:10:49.759185Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Calculate class weights based on frequency\ndef compute_class_weights(y_train):\n    class_counts = torch.bincount(y_train)\n    total = len(y_train)\n    weights = total / (len(class_counts) * class_counts.float())\n    return weights\n\n\nweights = compute_class_weights(_Y_train).to(device)\nbalanced_model = train_model(model, train_loader, val_loader, num_epochs=100, device=device, lr=1e-4, class_weights=weights)\n\n# load best model\nmodel.load_state_dict(best_model)\n\n# compute metric\nY_pred = model(_X_test.to(device)).argmax(1).cpu().numpy()\nscore = balanced_accuracy_score(_Y_test, Y_pred)\nf1 = f1_score(_Y_test, Y_pred, average='weighted')\n\n# plot confusion matrix\nconf_mat = confusion_matrix(_Y_test, Y_pred, normalize='true')\nplot_confusion_matrix(conf_mat, markers, f\"tCNN - Accuracy: {score:.4f}, F1: {f1:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:16:57.325450Z","iopub.execute_input":"2025-02-22T21:16:57.325746Z","iopub.status.idle":"2025-02-22T21:20:49.947095Z","shell.execute_reply.started":"2025-02-22T21:16:57.325723Z","shell.execute_reply":"2025-02-22T21:20:49.946212Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"subject = \"S004\"\nX_test = np.load(INPUT_DATA_DIR_PATH / f\"{subject}_X.npy\")\n\nY_pred = model(torch.FloatTensor(X_test).to(device)).argmax(1).cpu().numpy()\nprint(Y_pred[:5])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:20:49.948349Z","iopub.execute_input":"2025-02-22T21:20:49.948634Z","iopub.status.idle":"2025-02-22T21:20:50.035118Z","shell.execute_reply.started":"2025-02-22T21:20:49.948613Z","shell.execute_reply":"2025-02-22T21:20:50.034458Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:23:04.345878Z","iopub.execute_input":"2025-02-22T21:23:04.346233Z","iopub.status.idle":"2025-02-22T21:23:04.351798Z","shell.execute_reply.started":"2025-02-22T21:23:04.346209Z","shell.execute_reply":"2025-02-22T21:23:04.350941Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:23:40.865498Z","iopub.execute_input":"2025-02-22T21:23:40.865853Z","iopub.status.idle":"2025-02-22T21:23:40.885577Z","shell.execute_reply.started":"2025-02-22T21:23:40.865824Z","shell.execute_reply":"2025-02-22T21:23:40.884853Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"_markers","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:30:44.451934Z","iopub.execute_input":"2025-02-22T21:30:44.452295Z","iopub.status.idle":"2025-02-22T21:30:44.457489Z","shell.execute_reply.started":"2025-02-22T21:30:44.452267Z","shell.execute_reply":"2025-02-22T21:30:44.456629Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df[\"Marker\"].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:29:49.014276Z","iopub.execute_input":"2025-02-22T21:29:49.014567Z","iopub.status.idle":"2025-02-22T21:29:49.041795Z","shell.execute_reply.started":"2025-02-22T21:29:49.014545Z","shell.execute_reply":"2025-02-22T21:29:49.041034Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df[\"Marker\"] = df[\"Marker\"].apply(lambda x: _markers[x - 1])  # Adjust for 0-based indexing","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:42:38.960421Z","iopub.execute_input":"2025-02-22T21:42:38.960724Z","iopub.status.idle":"2025-02-22T21:42:38.965571Z","shell.execute_reply.started":"2025-02-22T21:42:38.960700Z","shell.execute_reply":"2025-02-22T21:42:38.964789Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:42:41.258701Z","iopub.execute_input":"2025-02-22T21:42:41.259036Z","iopub.status.idle":"2025-02-22T21:42:41.268523Z","shell.execute_reply.started":"2025-02-22T21:42:41.259012Z","shell.execute_reply":"2025-02-22T21:42:41.267706Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"IDs = np.arange(len(Y_pred))\n\ndf = pd.DataFrame({\"ID\":IDs, \"Marker\": Y_pred})\ndf.head()\ndf.to_csv(\"example_submission.csv\", index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:41:31.511631Z","iopub.execute_input":"2025-02-22T21:41:31.512047Z","iopub.status.idle":"2025-02-22T21:41:31.518328Z","shell.execute_reply.started":"2025-02-22T21:41:31.512016Z","shell.execute_reply":"2025-02-22T21:41:31.517502Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-22T21:43:54.446997Z","iopub.execute_input":"2025-02-22T21:43:54.447403Z","iopub.status.idle":"2025-02-22T21:43:54.457028Z","shell.execute_reply.started":"2025-02-22T21:43:54.447375Z","shell.execute_reply":"2025-02-22T21:43:54.456180Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}