{"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":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":7422302,"sourceType":"datasetVersion","datasetId":4318509}],"dockerImageVersionId":30635,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Introduction\n\n**eeglib** is a library for converting EEG waveforms into feature sets:\n\n[github: eeglib](https://github.com/Xiul109/eeglib)\n[eeglib: A Python module for EEG feature extraction](https://www.sciencedirect.com/science/article/pii/S2352711021000753)\n\nThis notebook uses **eeglib** to extract various features from the training set and\nvisualise them using violinplots.\n\nI've developed two plotting functions, one for left-right regions of the brain\nand one covering all of the spatial areas in covered by the eeg sensors in a pair-wise\nleft-right-front-back manner. These plot eeglib features against my interpretation\nof the probabilities.\n\n## Things to note\n\n1. **eeglib** is not in the standad Kaggle notebook libraries image fo you have to\n   offline pip install it using an offline pip install dataset. I've already created this\n   dataset here: [kaggle dataset: hms-hbac-offline-libs](https://www.kaggle.com/datasets/andrewscholan/hms-hbac-offline-libs).\n\n## tldr\n\nThe take on this is that eeglib features might be a useful decomposition of the eeg\nfile data but the data probably needs better preprocessing that eeglib preprocessing\nprovides.\n\nAnyway, there are a few pretty plots to look at!!!!\n\n## And finally\n\nIf you've found this notebook useful, please upvote it on kaggle.","metadata":{}},{"cell_type":"code","source":"!pip install \\\n   --requirement /kaggle/input/hms-hbac-offline-libs/requirements.txt \\\n   --no-index \\\n   --find-links file:///kaggle/input/hms-hbac-offline-libs/wheels","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-02T17:30:51.524468Z","iopub.execute_input":"2024-02-02T17:30:51.524836Z","iopub.status.idle":"2024-02-02T17:32:32.973261Z","shell.execute_reply.started":"2024-02-02T17:30:51.524804Z","shell.execute_reply":"2024-02-02T17:32:32.971095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# All imports in this code block\n\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport webcolors as wc\nimport math\nimport glob\n\nfrom koilerplate import INPUT_ROOT, WORKING_ROOT, TEMP_ROOT\nfrom pathlib import Path\nfrom enum import Enum\nfrom typing import List, Tuple, Dict\nfrom dataclasses import dataclass\nfrom tqdm.notebook import tqdm_notebook, tqdm\nfrom copy import deepcopy\nfrom pyarrow.parquet import ParquetDataset\n\nfrom eeglib.helpers import Helper\nfrom eeglib.eeg import EEG\n","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:32.976589Z","iopub.execute_input":"2024-02-02T17:32:32.977076Z","iopub.status.idle":"2024-02-02T17:32:40.324448Z","shell.execute_reply.started":"2024-02-02T17:32:32.977020Z","shell.execute_reply":"2024-02-02T17:32:40.322753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set up some basic file paths\nINPUT_PATH = Path(INPUT_ROOT)\nCOMPETITION_DATA_PATH = INPUT_PATH / \"hms-harmful-brain-activity-classification\"\nCOMPETITION_DATA_PATH","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:40.325850Z","iopub.execute_input":"2024-02-02T17:32:40.326408Z","iopub.status.idle":"2024-02-02T17:32:40.335659Z","shell.execute_reply.started":"2024-02-02T17:32:40.326375Z","shell.execute_reply":"2024-02-02T17:32:40.334232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the training set CSV\ntrain_info = pd.read_csv(COMPETITION_DATA_PATH/\"train.csv\")\ntrain_info","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:40.338435Z","iopub.execute_input":"2024-02-02T17:32:40.339753Z","iopub.status.idle":"2024-02-02T17:32:40.598457Z","shell.execute_reply.started":"2024-02-02T17:32:40.339698Z","shell.execute_reply":"2024-02-02T17:32:40.597089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Basic data visualisation\n\nThere's not a lot to look at in the basic training csv file. Just look at the expert categories to see if they are roughly the same order of magnitude each.","metadata":{}},{"cell_type":"markdown","source":"### Expert consensus\n\nJust look at the raw data from the training file.","metadata":{}},{"cell_type":"code","source":"# Visualise the expert consensus category\nsns.histplot(x=\"expert_consensus\", data=train_info)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:40.600270Z","iopub.execute_input":"2024-02-02T17:32:40.601543Z","iopub.status.idle":"2024-02-02T17:32:41.095789Z","shell.execute_reply.started":"2024-02-02T17:32:40.601486Z","shell.execute_reply":"2024-02-02T17:32:41.094878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Are all time offsets even ???\n\nJust check to see if there are any odd time offsets...","metadata":{}},{"cell_type":"code","source":"eeg_time_offset_odd = train_info[\"eeg_label_offset_seconds\"].mod(2) != 0\nspectrogram_time_offset_odd = train_info[\"spectrogram_label_offset_seconds\"].mod(2) != 0\nprint(f\"All EEG start time offsets are multiples of 2 sec: {eeg_time_offset_odd.sum()==0}\")\nprint(f\"All spectrogram start time offsets are multiples of 2 sec: {eeg_time_offset_odd.sum()==0}\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.096946Z","iopub.execute_input":"2024-02-02T17:32:41.098103Z","iopub.status.idle":"2024-02-02T17:32:41.111374Z","shell.execute_reply.started":"2024-02-02T17:32:41.097999Z","shell.execute_reply":"2024-02-02T17:32:41.109730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Labels\n\nAccording to the description in the dataset:\n\n- experts were given a 50 second extract of EEG waveform\n  data and each label corresponds to the centre 10 seconds of the sample.\n\n- the 50 second sample comes from a longer EEG waveform which may have a number\n  of labelled sections in it, each extracted as a snippet of 50 seconds with the\n  label corresponding to the centre 10 seconds.\n  \n- snippet sections may overlap.\n\n- labels could overlap.","metadata":{}},{"cell_type":"code","source":"labels = train_info['label_id'].unique()\nif len(labels)==train_info.shape[0]:\n    print(\"All label IDs are unique\")\nelse:\n    print(\"WARNING: Some label IDs are duplicated\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.113761Z","iopub.execute_input":"2024-02-02T17:32:41.114155Z","iopub.status.idle":"2024-02-02T17:32:41.129941Z","shell.execute_reply.started":"2024-02-02T17:32:41.114131Z","shell.execute_reply":"2024-02-02T17:32:41.128336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"eeg_ids = train_info['eeg_id'].unique()\nprint(f\"There are {len(eeg_ids)} EEG ids\")\nprint(f\"There are approximately {len(labels)/len(eeg_ids):.2f} labels per EEG id\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.131739Z","iopub.execute_input":"2024-02-02T17:32:41.132396Z","iopub.status.idle":"2024-02-02T17:32:41.139609Z","shell.execute_reply.started":"2024-02-02T17:32:41.132367Z","shell.execute_reply":"2024-02-02T17:32:41.138287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Add label start time and end times to the train info\nEEG_SNAPSHOT_DURATION = 50.0    # Experts analysed 50 second snapshots for labelling\nEEG_LABEL_DURATION = 10.0       # Labels relate to the centre 10 seconds of each 50s\ntrain_info[\"eeg_label_start_time\"] = (\n    train_info[\"eeg_label_offset_seconds\"] + \n    (EEG_SNAPSHOT_DURATION - EEG_LABEL_DURATION) / 2\n)\ntrain_info[\"eeg_label_end_time\"] = train_info[\"eeg_label_start_time\"] + EEG_LABEL_DURATION\ntrain_info","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.140788Z","iopub.execute_input":"2024-02-02T17:32:41.142091Z","iopub.status.idle":"2024-02-02T17:32:41.176492Z","shell.execute_reply.started":"2024-02-02T17:32:41.142023Z","shell.execute_reply":"2024-02-02T17:32:41.174920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"VOTES_AND_CONSENSUS = [\n    \"seizure_vote\", \n    \"lpd_vote\", \"gpd_vote\", \n    \"lrda_vote\", \"grda_vote\", \n    \"other_vote\", \n    \"expert_consensus\"\n]\n# First column of flags where there is an overlap between adjacent labels\ntrain_info[\"label_overlaps_next\"] = (\n    (train_info[\"eeg_id\"] == train_info[\"eeg_id\"].shift(-1)) &\n    (train_info[\"eeg_label_end_time\"] > train_info[\"eeg_label_start_time\"].shift(-1))\n)\n# Work out the overlap duration in seconds\ntrain_info[\"label_overlap_next_sec\"] = (\n    train_info[\"eeg_label_end_time\"] - train_info[\"eeg_label_start_time\"].shift(-1)\n).where(train_info[\"label_overlaps_next\"], 0)\ntrain_info","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.181433Z","iopub.execute_input":"2024-02-02T17:32:41.181790Z","iopub.status.idle":"2024-02-02T17:32:41.219421Z","shell.execute_reply.started":"2024-02-02T17:32:41.181767Z","shell.execute_reply":"2024-02-02T17:32:41.217528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"this_votes = train_info[[\"label_overlaps_next\"] + VOTES_AND_CONSENSUS]\nnext_votes = train_info[VOTES_AND_CONSENSUS].shift(-1)\nvoting = this_votes.join(next_votes, lsuffix=\"_current\", rsuffix=\"_next\")\nvoting[\"join_vote_conflict\"] = (\n    voting[\"label_overlaps_next\"] & (\n        (voting[\"seizure_vote_current\"] != voting[\"seizure_vote_next\"])\n        | (voting[\"lpd_vote_current\"] != voting[\"lpd_vote_next\"])\n        | (voting[\"gpd_vote_current\"] != voting[\"gpd_vote_next\"])\n        | (voting[\"lrda_vote_current\"] != voting[\"lrda_vote_next\"])\n        | (voting[\"grda_vote_current\"] != voting[\"grda_vote_next\"])\n        | (voting[\"other_vote_current\"] != voting[\"other_vote_next\"])\n    )\n)\nvoting[\"join_consensus_conflict\"] = (\n    voting[\"label_overlaps_next\"]\n    & (voting[\"expert_consensus_current\"] != voting[\"expert_consensus_next\"])\n)\nnum_vote_conflicts = voting['join_vote_conflict'].sum()\nnum_consensus_conflicts = voting['join_consensus_conflict'].sum()\nnum_overlaps = voting['label_overlaps_next'].sum()\n\nprint(\n    f\"{num_vote_conflicts} conflicts in expert voting when labels overlap \"\n    f\"({num_vote_conflicts/num_overlaps :.3%}).\"\n)\nprint(\n    f\"{num_consensus_conflicts} conflicts in expert consensus when labels overlap \"\n    f\"({num_consensus_conflicts/num_overlaps :.3%}).\"\n)\nprint(f\"{num_overlaps} total number of overlapping labels.\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.221211Z","iopub.execute_input":"2024-02-02T17:32:41.221594Z","iopub.status.idle":"2024-02-02T17:32:41.273130Z","shell.execute_reply.started":"2024-02-02T17:32:41.221563Z","shell.execute_reply":"2024-02-02T17:32:41.271133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Discussion\n\nThere are various conflicts between labels when they are merged, buy not as many as\nwe might fear.\n\nInterestingly, if you assume that an expert focusses on a smaller central region other\nthan the full 10 seconds of the label the conflicts reduce significantly.\n\n| Label length | Voting Conflicts | Consensus Conflicts | Overlapping labels |\n| ------------ | ---------------- | ------------------- | ------------------ |\n|        4 sec |     306 (0.729%) |        109 (0.260%) |              41960 |\n|        6 sec |     767 (1.279%) |        255 (0.425%) |              59966 |\n|        8 sec |    1247 (1.795%) |        427 (0.615%) |              69479 |\n|       10 sec |    1739 (2.314%) |        594 (0.790%) |              71546 |\n\n> **Note:** Above table not coded in this notebook. Shows results of some coding experiments.","metadata":{}},{"cell_type":"markdown","source":"## Dropping conflicted data\n\nThe easiest way to deal with the conflicted labels in the training set is to\ndrop those labels where there are conflicts in the consensus (which is less than\n1% of the data) as we treat them as outliers in the dataset.","metadata":{}},{"cell_type":"code","source":"train_info_clean = train_info[voting[\"join_consensus_conflict\"] != True].copy()\ntrain_info_clean","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.274779Z","iopub.execute_input":"2024-02-02T17:32:41.275181Z","iopub.status.idle":"2024-02-02T17:32:41.331934Z","shell.execute_reply.started":"2024-02-02T17:32:41.275150Z","shell.execute_reply":"2024-02-02T17:32:41.329981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Votes cast\n\nWe only know that there is a panel of experts, some presumably will not have voted and some may\nvote more based on either the spectrogram or the raw eeg.","metadata":{}},{"cell_type":"code","source":"# Do we know how many experts looked at each eeg trace?\n# Start by simply counting the votes for each entry \ntrain_info_clean[\"vote_count\"] = (\n    train_info_clean[\"seizure_vote\"] +\n    train_info_clean[\"lpd_vote\"] +\n    train_info_clean[\"gpd_vote\"] +\n    train_info_clean[\"lrda_vote\"] +\n    train_info_clean[\"grda_vote\"] +\n    train_info_clean[\"other_vote\"]\n)\ntrain_info_clean","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.334165Z","iopub.execute_input":"2024-02-02T17:32:41.335351Z","iopub.status.idle":"2024-02-02T17:32:41.367702Z","shell.execute_reply.started":"2024-02-02T17:32:41.335304Z","shell.execute_reply":"2024-02-02T17:32:41.366319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now get the maximum vote count based on EEG ID and Spectrogram ID\ntrain_info_clean[\"max_vote_count_eeg\"] = train_info_clean.groupby(\"eeg_id\")[\"vote_count\"].transform(\"max\")\ntrain_info_clean[\"max_vote_count_spectrogram\"] = train_info_clean.groupby(\"spectrogram_id\")[\"vote_count\"].transform(\"max\")\ntrain_info_clean","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.368918Z","iopub.execute_input":"2024-02-02T17:32:41.369207Z","iopub.status.idle":"2024-02-02T17:32:41.445844Z","shell.execute_reply.started":"2024-02-02T17:32:41.369184Z","shell.execute_reply":"2024-02-02T17:32:41.444830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now test how many entries there are where the raw count is different from the max\ntotal_rows = train_info_clean.shape[0]\nfewer_eeg_votes = train_info_clean[train_info_clean[\"vote_count\"] < train_info_clean[\"max_vote_count_eeg\"]]\nfewer_eeg_votes_rows = fewer_eeg_votes.shape[0]\nfewer_spectrogram_votes = train_info_clean[train_info_clean[\"vote_count\"] < train_info_clean[\"max_vote_count_spectrogram\"]]\nfewer_spectrogram_votes_rows = fewer_spectrogram_votes.shape[0]\ninconsistent_max_votes = train_info_clean[train_info_clean[\"max_vote_count_eeg\"] != train_info_clean[\"max_vote_count_spectrogram\"]]\ninconsistent_max_votes_rows = inconsistent_max_votes.shape[0]\nprint(f\"{fewer_eeg_votes_rows=}/{total_rows}; {fewer_eeg_votes_rows/total_rows : .2%}\")\nprint(f\"{fewer_spectrogram_votes_rows=}/{total_rows}; {fewer_spectrogram_votes_rows/total_rows : .2%}\")\nprint(f\"{inconsistent_max_votes_rows=}/{total_rows}; {inconsistent_max_votes_rows/total_rows : .2%}\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.447271Z","iopub.execute_input":"2024-02-02T17:32:41.447569Z","iopub.status.idle":"2024-02-02T17:32:41.472867Z","shell.execute_reply.started":"2024-02-02T17:32:41.447539Z","shell.execute_reply":"2024-02-02T17:32:41.471005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Explore further the inconsistent max votes\ninconsistent_fewer_eeg_votes = inconsistent_max_votes[\n    inconsistent_max_votes[\"max_vote_count_eeg\"] < inconsistent_max_votes[\"max_vote_count_spectrogram\"]\n]\ninconsistent_fewer_eeg_votes_rows = inconsistent_fewer_eeg_votes.shape[0]\ninconsistent_fewer_spectrogram_votes = inconsistent_max_votes[\n    inconsistent_max_votes[\"max_vote_count_spectrogram\"] < inconsistent_max_votes[\"max_vote_count_eeg\"]\n]\ninconsistent_fewer_spectrogram_votes_rows = inconsistent_fewer_spectrogram_votes.shape[0]\nprint(\n    f\"{inconsistent_fewer_eeg_votes_rows=}/{inconsistent_max_votes_rows}; \"\n    f\"{inconsistent_fewer_eeg_votes_rows/inconsistent_max_votes_rows : .2%}\"\n)\nprint(\n    f\"{inconsistent_fewer_spectrogram_votes_rows=}/{inconsistent_max_votes_rows}; \"\n    f\"{inconsistent_fewer_spectrogram_votes_rows/inconsistent_max_votes_rows : .2%}\"\n)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.474561Z","iopub.execute_input":"2024-02-02T17:32:41.474955Z","iopub.status.idle":"2024-02-02T17:32:41.488711Z","shell.execute_reply.started":"2024-02-02T17:32:41.474923Z","shell.execute_reply":"2024-02-02T17:32:41.487111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Discussion\n\nSo, what do we know:\n1. Sometimes experts will not agree on a label and will not cast a vote.\n2. The spectrogram IDs have more votes associated with them than the eeg ids.\n   this is probably consistent as the spectrograms cover a longer period of time\n   than the EEG traces.\n3. This probably means that the size of the panel could be derived from the max\n   number of votes associated with the spectrogram ID rather than the eeg ID.\n   \nTherefore, to convert our votes to probabilities we use the max count of spectrogram ID votes.","metadata":{}},{"cell_type":"markdown","source":"### Probabilities\n\nWe'll divide the vote numbers by the maximum number of votes grouped by spectrogram IDs","metadata":{}},{"cell_type":"code","source":"VOTE_COLUMNS = [\"seizure_vote\", \"lpd_vote\", \"gpd_vote\", \"lrda_vote\", \"grda_vote\", \"other_vote\"]\nP_COLUMNS = [\"P_sz\", \"P_lpd\", \"P_gpd\", \"P_lrda\", \"P_grda\", \"P_other\"]\nfor p_col, v_col in zip(P_COLUMNS, VOTE_COLUMNS):\n    train_info_clean[p_col] = train_info_clean[v_col] / train_info_clean[\"max_vote_count_spectrogram\"]\ntrain_info_clean = train_info_clean.copy()\ntrain_info_clean","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.490916Z","iopub.execute_input":"2024-02-02T17:32:41.491439Z","iopub.status.idle":"2024-02-02T17:32:41.590811Z","shell.execute_reply.started":"2024-02-02T17:32:41.491399Z","shell.execute_reply":"2024-02-02T17:32:41.589775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Work out probabilities across the whole training set\nP_all = np.asarray(\n    [train_info_clean[p_col].sum() for p_col in P_COLUMNS]\n) / total_rows\nP_all","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.592588Z","iopub.execute_input":"2024-02-02T17:32:41.593201Z","iopub.status.idle":"2024-02-02T17:32:41.603701Z","shell.execute_reply.started":"2024-02-02T17:32:41.593163Z","shell.execute_reply":"2024-02-02T17:32:41.602492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualise as a table\nP_df = pd.DataFrame(data=P_all.reshape((1,6)), columns=[\"P_sz\", \"P_lpd\", \"P_gpd\", \"P_lrda\", \"P_grda\", \"P_other\"])\nP_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.604884Z","iopub.execute_input":"2024-02-02T17:32:41.605178Z","iopub.status.idle":"2024-02-02T17:32:41.621382Z","shell.execute_reply.started":"2024-02-02T17:32:41.605152Z","shell.execute_reply":"2024-02-02T17:32:41.619867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# What's the probability sum across everything, don't expect this to be 1.0, but something close-ish\nP_all.sum()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.623876Z","iopub.execute_input":"2024-02-02T17:32:41.624626Z","iopub.status.idle":"2024-02-02T17:32:41.635684Z","shell.execute_reply.started":"2024-02-02T17:32:41.624574Z","shell.execute_reply":"2024-02-02T17:32:41.633539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Level of agreement\n\nFrom the competition overview:\n\n- 'idealized': High level of expert agreement\n- 'proto patterns':  Cases where ~1/2 of experts give a label as “other” and ~1/2\n   give one of the remaining five labels.\n- 'edge cases': Where experts are approximately split between 2 of the 5 named patterns\n\n> Not that easy to program!","metadata":{}},{"cell_type":"code","source":"# My stab at categorizing the agreement based on the vague description\ndef compute_agreement(P_sz:float, P_lpd:float, P_gpd:float, P_lrda:float, P_grda:float, P_other:float) -> str:\n    agreement = \"none\"\n    # Because of the way we have computed the probabilities for a row,\n    # the parameters passed may not add to 1.0. Lets fix this first.\n    P_tot = P_sz + P_lpd + P_gpd + P_lrda + P_grda + P_other\n    P_sz = P_sz / P_tot\n    P_lpd = P_lpd / P_tot\n    P_gpd = P_gpd / P_tot\n    P_lrda = P_lrda / P_tot\n    P_grda = P_grda / P_tot\n    P_other = P_other / P_tot\n    # Now rank them by probability value\n    p_dict = { \"sz\": P_sz, \"lpd\": P_lpd, \"gpd\": P_gpd, \"lrda\": P_lrda, \"grda\": P_grda, \"other\": P_other}\n    p_list = list(p_dict.items())\n    p_list.sort(key=lambda kv: 1-kv[1])\n    # Ignore all but the three highest probabilities, get the top ranking categories\n    first = p_list[0][0]\n    second = p_list[1][0]\n    third = p_list[2][0]\n    # Re-scale by the ignored probabilities of the lower 3 categories\n    p_first = p_dict[first]\n    p_second = p_dict[second]\n    p_third = p_dict[third]\n    p_top_3 = p_first + p_second + p_third\n    p_first /= p_top_3\n    p_second /= p_top_3\n    p_third /= p_top_3\n    if p_first > 0.75:\n        # If not other then there is strong agreement, if is other then no agreement\n        if first != \"other\":\n            agreement = \"idealized\"\n    elif (p_first + p_second) > 0.75:\n        # First two categories combined are significant\n        if first == \"other\" or second == \"other\":\n            # We'll assume that this is a proto-pattern\n            agreement = \"proto-pattern\"\n        else:\n            # We have two of equal-ish ranking\n            agreement = \"edge-case\"\n    else:\n        # No agreement (already set)\n        ...\n    return agreement","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.639204Z","iopub.execute_input":"2024-02-02T17:32:41.639662Z","iopub.status.idle":"2024-02-02T17:32:41.651459Z","shell.execute_reply.started":"2024-02-02T17:32:41.639627Z","shell.execute_reply":"2024-02-02T17:32:41.649519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_info_clean[\"agreement\"] = train_info_clean.apply(\n    lambda row : compute_agreement(\n        row[\"P_sz\"], row[\"P_lpd\"], row[\"P_gpd\"], row[\"P_lrda\"], row[\"P_grda\"], row[\"P_other\"]\n    ),\n    axis=1\n)\ntrain_info_clean","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:41.653352Z","iopub.execute_input":"2024-02-02T17:32:41.653678Z","iopub.status.idle":"2024-02-02T17:32:44.293800Z","shell.execute_reply.started":"2024-02-02T17:32:41.653654Z","shell.execute_reply":"2024-02-02T17:32:44.292828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualise the expert consensus category\nsns.histplot(x=\"agreement\", data=train_info_clean)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:44.295208Z","iopub.execute_input":"2024-02-02T17:32:44.295492Z","iopub.status.idle":"2024-02-02T17:32:44.579097Z","shell.execute_reply.started":"2024-02-02T17:32:44.295466Z","shell.execute_reply":"2024-02-02T17:32:44.577848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## How long are the EEG files?\n\nThe competition says that some of the EEG files have been merged so it would be useful to\nknow the range of file sizes (which gives up the number of samples).","metadata":{}},{"cell_type":"code","source":"# The EEG sampling rate is given in the competition data set info\nEEG_SAMPLING_RATE = 200","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:44.580712Z","iopub.execute_input":"2024-02-02T17:32:44.581759Z","iopub.status.idle":"2024-02-02T17:32:44.587224Z","shell.execute_reply.started":"2024-02-02T17:32:44.581715Z","shell.execute_reply":"2024-02-02T17:32:44.585571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Reading in training EEGs as a ParquetDataset, this takes some time...\")\neeg_files_dataset = ParquetDataset(COMPETITION_DATA_PATH / \"train_eegs\")\nprint(\"Now reading the length of each parquet file, also takes time...\")\nfile_durations:Dict[float, int] = {}\nfor fragment in tqdm_notebook(eeg_files_dataset.fragments):\n    file_duration = fragment.count_rows() / EEG_SAMPLING_RATE\n    if file_duration in file_durations:\n        file_durations[file_duration] += 1\n    else:\n        file_durations[file_duration] = 1 \nfile_durations_df = pd.DataFrame(sorted([(k,v) for k,v in file_durations.items()]), columns=[\"duration\", \"file_count\"])\nfile_durations_df ","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:32:44.589291Z","iopub.execute_input":"2024-02-02T17:32:44.589836Z","iopub.status.idle":"2024-02-02T17:34:39.939084Z","shell.execute_reply.started":"2024-02-02T17:32:44.589789Z","shell.execute_reply":"2024-02-02T17:34:39.937823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We're going to categorise the file length into various buckets...\ndef categorise_duration(duration: float, sort_order:bool = False) -> str:\n    if duration == 50.0:\n        bucket = (0, \"50s\")\n    elif duration <= 60.0:\n        bucket = (1, \"52-60s\")\n    elif duration <= 90.0:\n        bucket = (2, \"60-90s\")\n    elif duration <= 120.0:\n        bucket = (3, \"90s-2m\")\n    elif duration <= 180.0:\n        bucket = (4, \"2-3m\")\n    elif duration <= 300.0:\n        bucket = (5, \"3-5m\")\n    elif duration <= 600.0:\n        bucket = (6, \"5-10m\")\n    elif duration <= 1200.0:\n        bucket = (7, \"10-20m\")\n    else:\n        bucket = (8, \"> 20m\")\n    return bucket[0 if sort_order else 1]","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:39.940437Z","iopub.execute_input":"2024-02-02T17:34:39.941493Z","iopub.status.idle":"2024-02-02T17:34:39.950157Z","shell.execute_reply.started":"2024-02-02T17:34:39.941455Z","shell.execute_reply":"2024-02-02T17:34:39.948566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_durations_df[\"duration_bucket\"] = file_durations_df.apply(\n    lambda row: categorise_duration(row[\"duration\"]), axis=1\n)\nfile_durations_df[\"bucket_order\"] = file_durations_df.apply(\n    lambda row: categorise_duration(row[\"duration\"], True), axis=1\n)\nfile_durations_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:39.951457Z","iopub.execute_input":"2024-02-02T17:34:39.951759Z","iopub.status.idle":"2024-02-02T17:34:39.975895Z","shell.execute_reply.started":"2024-02-02T17:34:39.951734Z","shell.execute_reply":"2024-02-02T17:34:39.975223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at the spread of file durations\naxis = sns.barplot(\n    data=file_durations_df, \n    x=\"duration_bucket\", \n    y=\"file_count\", \n    estimator=\"sum\",\n    errorbar=None\n)\naxis.bar_label(axis.containers[0])\nprint(f\"Total number of EEG files: {file_durations_df['file_count'].sum()}\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:39.984237Z","iopub.execute_input":"2024-02-02T17:34:39.985011Z","iopub.status.idle":"2024-02-02T17:34:40.251006Z","shell.execute_reply.started":"2024-02-02T17:34:39.984964Z","shell.execute_reply":"2024-02-02T17:34:40.249467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## EEG Data exploration\n\nLoad in an arbitrary sample and have a look at the data","metadata":{}},{"cell_type":"code","source":"# Just choose a arbitrary sample to look at\nSAMPLE = 7_824\nsample_info = train_info_clean.iloc[SAMPLE]\nsample_info","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.253228Z","iopub.execute_input":"2024-02-02T17:34:40.253627Z","iopub.status.idle":"2024-02-02T17:34:40.263508Z","shell.execute_reply.started":"2024-02-02T17:34:40.253593Z","shell.execute_reply":"2024-02-02T17:34:40.262535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TIME_OFFSET_COLUMN = \"time_offset\"\nEEG_SAMPLING_PERIOD = 1.0 / EEG_SAMPLING_RATE\n\n# Function to add a time-channel to a dataframe\ndef add_time_channel(df: pd.DataFrame, period: float=EEG_SAMPLING_PERIOD) -> None:\n    df[TIME_OFFSET_COLUMN] = df.index * period","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.264621Z","iopub.execute_input":"2024-02-02T17:34:40.264888Z","iopub.status.idle":"2024-02-02T17:34:40.276454Z","shell.execute_reply.started":"2024-02-02T17:34:40.264862Z","shell.execute_reply":"2024-02-02T17:34:40.274254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Per second probabilities\n\nWe have labels that overlap in any given EEG file and we know that there are\ninconsistencies still between labels, however, we have dropped the overlapping\nlabels where the expert consensus \"agreement\" category differs.\n\nWe still need to translate out tabular label probabilities into something we can\nuse for training a model. Let's give each second in an EEG file a label probability\nbased on whatever labels happen to overlap.","metadata":{}},{"cell_type":"code","source":"def combine_probabilities(info: pd.DataFrame, row: pd.Series) -> pd.Series:\n    time_offset = float(row.name)   # This is the index value for the row\n    applicable_probabilities = info[P_COLUMNS].where(\n        (info[\"eeg_label_start_time\"]<=time_offset) & \n        (time_offset<info[\"eeg_label_end_time\"])\n    )\n    # Note apply fillna **after** mean as we don't want NaN rows to count towards the mean\n    return applicable_probabilities[P_COLUMNS].mean().fillna(0.0)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.278102Z","iopub.execute_input":"2024-02-02T17:34:40.278458Z","iopub.status.idle":"2024-02-02T17:34:40.288838Z","shell.execute_reply.started":"2024-02-02T17:34:40.278425Z","shell.execute_reply":"2024-02-02T17:34:40.287931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def eeg_P_per_sec(train_info: pd.DataFrame, eeg_id: int):\n    eeg_info = train_info[train_info[\"eeg_id\"]==eeg_id].copy()\n    last_label_start_offset = eeg_info[\"eeg_label_offset_seconds\"].iloc[-1]\n    # Note as we know all label boundaries are every 2 secs we step index by 2\n    p_per_sec = pd.DataFrame(0.0, index=np.arange(0.0, last_label_start_offset+EEG_SNAPSHOT_DURATION, 2.0), columns=P_COLUMNS)\n    p_per_sec = p_per_sec.apply(lambda row: combine_probabilities(eeg_info, row), axis=1)\n    p_per_sec[\"eeg_id\"] = eeg_id\n    p_per_sec[\"is_scored\"] = p_per_sec[P_COLUMNS].sum(axis=1) > 0\n    p_per_sec.replace(0.0, math.nan, inplace=True)\n    return p_per_sec.copy()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.289802Z","iopub.execute_input":"2024-02-02T17:34:40.290141Z","iopub.status.idle":"2024-02-02T17:34:40.305202Z","shell.execute_reply.started":"2024-02-02T17:34:40.290107Z","shell.execute_reply":"2024-02-02T17:34:40.304087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_eeg_p_per_sec = eeg_P_per_sec(train_info_clean, sample_info[\"eeg_id\"])\nsample_eeg_p_per_sec","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.306714Z","iopub.execute_input":"2024-02-02T17:34:40.307059Z","iopub.status.idle":"2024-02-02T17:34:40.435905Z","shell.execute_reply.started":"2024-02-02T17:34:40.307015Z","shell.execute_reply":"2024-02-02T17:34:40.434480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.lineplot(data=sample_eeg_p_per_sec[P_COLUMNS])","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.437942Z","iopub.execute_input":"2024-02-02T17:34:40.438432Z","iopub.status.idle":"2024-02-02T17:34:40.780977Z","shell.execute_reply.started":"2024-02-02T17:34:40.438394Z","shell.execute_reply":"2024-02-02T17:34:40.779425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load a single 50s snapshot for the EEG we've selected\n\nThis is the 50 second slice that the experts have seen.","metadata":{}},{"cell_type":"code","source":"# Function to load the eeg data\ndef load_eeg(\n    eeg_id:int, \n    start_time:float=0, \n    duration:float|None=None\n) -> pd.DataFrame:\n    df = pd.read_parquet(COMPETITION_DATA_PATH/\"train_eegs\"/f\"{eeg_id}.parquet\")\n    # Interpolated over the dataframe as some of the eeg files have the \n    # odd row of NaNs\n    df = df.interpolate()\n    # And front and back data may contain NaNs that won't interpolate\n    df = df.fillna(0.0)\n    # Add the time channel\n    add_time_channel(df)\n    # Get the duration for the whole file\n    file_duration = df[TIME_OFFSET_COLUMN].iloc[-1] + EEG_SAMPLING_PERIOD\n    # Now we want to extract the data for period in question\n    end_time = (\n        start_time + \n        (duration if duration is not None else file_duration - start_time)\n    )\n    df = df[\n        (df[TIME_OFFSET_COLUMN] >= start_time) & \n        (df[TIME_OFFSET_COLUMN] < end_time)\n    ]\n    return df.copy().reset_index()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.782947Z","iopub.execute_input":"2024-02-02T17:34:40.783449Z","iopub.status.idle":"2024-02-02T17:34:40.792902Z","shell.execute_reply.started":"2024-02-02T17:34:40.783409Z","shell.execute_reply":"2024-02-02T17:34:40.791352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the first 50 seconds of the sample EEG data using function\nsnapshot_eeg = load_eeg(\n    sample_info['eeg_id'], \n    sample_info['eeg_label_offset_seconds'],\n    EEG_SNAPSHOT_DURATION\n)\nsnapshot_eeg","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.794281Z","iopub.execute_input":"2024-02-02T17:34:40.794600Z","iopub.status.idle":"2024-02-02T17:34:40.888590Z","shell.execute_reply.started":"2024-02-02T17:34:40.794574Z","shell.execute_reply":"2024-02-02T17:34:40.887096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Single ended vs differential\n\nLooks like the data is single-ended (which is good) but we probably need to make it differential to produce the same type of plots in the sample data.","metadata":{}},{"cell_type":"code","source":"# What are the column names?\nfor col in snapshot_eeg.columns:\n    print(col)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.889975Z","iopub.execute_input":"2024-02-02T17:34:40.890442Z","iopub.status.idle":"2024-02-02T17:34:40.896487Z","shell.execute_reply.started":"2024-02-02T17:34:40.890398Z","shell.execute_reply":"2024-02-02T17:34:40.895503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Left lateral (LL)\nLL_EEG_CHANNELS = [\"Fp1-F7\", \"F7-T3\", \"T3-T5\", \"T5-O1\"]\n# Left parasagittal (LP)\nLP_EEG_CHANNELS = [\"Fp1-F3\", \"F3-C3\", \"C3-P3\", \"P3-O1\"]\n# Central\nCC_EEG_CHANNELS = [\"Fz-Cz\", \"Cz-Pz\"]\n# Right parasagittal (RP)\nRP_EEG_CHANNELS = [\"Fp2-F4\", \"F4-C4\", \"C4-P4\", \"P4-O2\"]\n# Right lateral (RL)\nRL_EEG_CHANNELS = [\"Fp2-F8\", \"F8-T4\", \"T4-T6\", \"T6-O2\"]\n# Auxiliary information columns\nAUX_EEG_COLUMNS = [ \"EKG\", TIME_OFFSET_COLUMN ]\n\n# Define how we want our EEG channels to be constructed\nDIFF_EEG_COLUMNS = \\\n    LL_EEG_CHANNELS + \\\n    LP_EEG_CHANNELS + \\\n    CC_EEG_CHANNELS + \\\n    RP_EEG_CHANNELS + \\\n    RL_EEG_CHANNELS + \\\n    AUX_EEG_COLUMNS\nDIFF_EEG_COLUMNS","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.897496Z","iopub.execute_input":"2024-02-02T17:34:40.897867Z","iopub.status.idle":"2024-02-02T17:34:40.912871Z","shell.execute_reply.started":"2024-02-02T17:34:40.897830Z","shell.execute_reply":"2024-02-02T17:34:40.911824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make a new dataframe with differential channels, rather than single ended\ndef make_differential(df:pd.DataFrame, columns:List[str] = DIFF_EEG_COLUMNS):\n    # Take copy as df is potentially a slice\n    df = df.copy();\n    to_drop: List[str] = []\n    for column in columns:\n        single_ended_channels = column.split(\"-\")\n        if len(single_ended_channels) == 2:\n            # We need to make the differential channel\n            df[column] = df[single_ended_channels[0]] - df[single_ended_channels[1]]\n            if single_ended_channels[0] not in to_drop:\n                to_drop.append(single_ended_channels[0])\n            if single_ended_channels[1] not in to_drop:\n                to_drop.append(single_ended_channels[1])\n            # Standardize column value\n            df[column] = (df[column] - df[column].mean()) / df[column].std()\n    # Drop non-differential columns\n    df.drop(to_drop, axis=1, inplace=True)\n    # Now return the dataframe in the correct column order\n    return df[columns]    ","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.914267Z","iopub.execute_input":"2024-02-02T17:34:40.914994Z","iopub.status.idle":"2024-02-02T17:34:40.925964Z","shell.execute_reply.started":"2024-02-02T17:34:40.914956Z","shell.execute_reply":"2024-02-02T17:34:40.925096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply the differential function\nextracted_diff_eeg = make_differential(snapshot_eeg)\nextracted_diff_eeg","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.927538Z","iopub.execute_input":"2024-02-02T17:34:40.928414Z","iopub.status.idle":"2024-02-02T17:34:40.993782Z","shell.execute_reply.started":"2024-02-02T17:34:40.928385Z","shell.execute_reply":"2024-02-02T17:34:40.992510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plotting\n\nWe're going to define a plotting process that can plot the dataframe columns as a grouped, stacked line plot using a combination of seaborn and matplot libraries.","metadata":{}},{"cell_type":"code","source":"# Define s dataclass that defines a group of channels to be plotted\n@dataclass\nclass PlotGroup:\n    channels: List[str]\n    fg_color: str\n    bg_color: str","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:40.995221Z","iopub.execute_input":"2024-02-02T17:34:40.995526Z","iopub.status.idle":"2024-02-02T17:34:41.002836Z","shell.execute_reply.started":"2024-02-02T17:34:40.995492Z","shell.execute_reply":"2024-02-02T17:34:41.001277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This function will plot a group of channels from the dataframe\ndef eeg_group_plot(\n    df: pd.DataFrame, \n    x: str, \n    plot_group: PlotGroup, \n    snapshot_offset:float,\n    axes:List[plt.axis], \n    plot_span_secs:float,\n    snapshot_duration\n):\n    xlim_lower = snapshot_offset + (snapshot_duration / 2) - (plot_span_secs / 2)\n    xlim_upper = xlim_lower + plot_span_secs\n    for channel, ax in zip(plot_group.channels, axes):\n        sns.lineplot(data=df, x=x, y=channel, ax=ax, color=plot_group.fg_color)\n        ax.set_facecolor(plot_group.bg_color)\n        ax.set_xlim((xlim_lower, xlim_upper))\n        ax.set_ylabel(channel, rotation=0, fontsize=12, horizontalalignment='right', verticalalignment='center')","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:41.004698Z","iopub.execute_input":"2024-02-02T17:34:41.005163Z","iopub.status.idle":"2024-02-02T17:34:41.018339Z","shell.execute_reply.started":"2024-02-02T17:34:41.005129Z","shell.execute_reply":"2024-02-02T17:34:41.016583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This function will plot and stack multiple groups together\ndef eeg_plot(\n    df: pd.DataFrame, \n    x: str, \n    plot_groups: List[plt.axis], \n    snapshot_offset: float, \n    plot_span_secs:float = EEG_LABEL_DURATION,\n    snapshot_duration = EEG_SNAPSHOT_DURATION,\n    fig_height=12, \n    fig_width=12\n):\n    num_channels = 0;\n    for plot_group in plot_groups:\n        num_channels += len(plot_group.channels)\n    figure, axes = plt.subplots(num_channels, 1)\n    figure.subplots_adjust(hspace=0)\n    figure.set_figheight(fig_height)\n    figure.set_figwidth(fig_width)\n    axis_num = 0\n    for plot_group in plot_groups:\n        num_group_chans = len(plot_group.channels)\n        group_axes = axes[axis_num:axis_num+num_group_chans]\n        eeg_group_plot(\n            df, x, plot_group, snapshot_offset, group_axes,\n            plot_span_secs=plot_span_secs,\n            snapshot_duration=snapshot_duration\n        )\n        axis_num += num_group_chans\n    figure.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:41.019945Z","iopub.execute_input":"2024-02-02T17:34:41.020338Z","iopub.status.idle":"2024-02-02T17:34:41.037273Z","shell.execute_reply.started":"2024-02-02T17:34:41.020301Z","shell.execute_reply":"2024-02-02T17:34:41.034933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now define what the plot groupings are and what the colours are to use.\n# We'll plot from the left side of the head to the right side, front to back.\n# We'll use red for left size and green for right (as per international navigation lights!)\n# EKG will be last channel\nDIFF_PLOT_GROUPS = [\n    PlotGroup([\"Fp1-F7\", \"F7-T3\", \"T3-T5\", \"T5-O1\"], 'red', wc.CSS3_NAMES_TO_HEX[\"lightpink\"]),\n    PlotGroup([\"Fp1-F3\", \"F3-C3\", \"C3-P3\", \"P3-O1\"], 'red', wc.CSS3_NAMES_TO_HEX[\"lightsalmon\"]),\n    PlotGroup([\"Fz-Cz\", \"Cz-Pz\"], 'black', wc.CSS3_NAMES_TO_HEX[\"gainsboro\"]),\n    PlotGroup([\"Fp2-F4\", \"F4-C4\", \"C4-P4\", \"P4-O2\"], 'green', wc.CSS3_NAMES_TO_HEX[\"darkseagreen\"]),\n    PlotGroup([\"Fp2-F8\", \"F8-T4\", \"T4-T6\", \"T6-O2\"], 'green', wc.CSS3_NAMES_TO_HEX[\"lightseagreen\"]),\n    PlotGroup([\"EKG\"], 'blue', wc.CSS3_NAMES_TO_HEX[\"powderblue\"]),\n]","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:41.039731Z","iopub.execute_input":"2024-02-02T17:34:41.040132Z","iopub.status.idle":"2024-02-02T17:34:41.056965Z","shell.execute_reply.started":"2024-02-02T17:34:41.040100Z","shell.execute_reply":"2024-02-02T17:34:41.055888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now plot the unprocessed data\neeg_plot(extracted_diff_eeg, TIME_OFFSET_COLUMN, DIFF_PLOT_GROUPS, sample_info['eeg_label_offset_seconds'])","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:41.058219Z","iopub.execute_input":"2024-02-02T17:34:41.058514Z","iopub.status.idle":"2024-02-02T17:34:44.675254Z","shell.execute_reply.started":"2024-02-02T17:34:41.058486Z","shell.execute_reply":"2024-02-02T17:34:44.674161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## EEGLIB preprocessing\n\neeglib has sparse documentation and may/may not be useful. Let's just explore what can be done with the preprocessing functions it provides.","metadata":{}},{"cell_type":"markdown","source":"### eeglib Helper class\n\nLooks like data has to start in a `Helper` object and we need to create this from a \nnumpy array (as there's no way to directly import a parquet file).\n\n**However:** If you simply create a `Helper` object from a numpy array the `Helper` \nonly saves the reference to the data and then, potentially, modifies it during \npre-processing. This is not a problem in a forward-running pipeline but when messing\naround in the Notebook this can have unintended side-effects.\n\nTherefore we create a helper function that safely creates the helper from the\nDataFrame without danger of modifying the input dataframe.","metadata":{}},{"cell_type":"code","source":"@dataclass\nclass BandPass:\n    low_cutoff: float|None\n    high_cutoff: float|None\n        \nDEFAULT_BANDPASS = BandPass(1.0, 50.0)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:38:05.279731Z","iopub.execute_input":"2024-02-02T17:38:05.280370Z","iopub.status.idle":"2024-02-02T17:38:05.285176Z","shell.execute_reply.started":"2024-02-02T17:38:05.280344Z","shell.execute_reply":"2024-02-02T17:38:05.284223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Makes an eeglib helper without potential side effects on the input dataframe\n# Returns the helper and the data dropped when we made the helper\ndef df_to_eeglib_helper(\n    df: pd.DataFrame,\n    columns: List[str]=DIFF_EEG_COLUMNS,\n    drop_columns: List[str]=AUX_EEG_COLUMNS,\n    sample_rate:int=EEG_SAMPLING_RATE, \n    window_size:int|None=int(EEG_LABEL_DURATION * EEG_SAMPLING_RATE),\n    band_pass:BandPass|None=None,\n    normalize:bool=False,\n    ica:bool=False\n) -> Helper:\n    required_cols = [col for col in columns if col not in drop_columns]\n    # Here is all important copy, gives us a new data array to be mutated by eeglib\n    copy_df = df[required_cols].copy()\n    dropped_df = df[drop_columns].copy()\n    data = copy_df.to_numpy().transpose()\n    helper = Helper(\n        data, \n        sampleRate=sample_rate, \n        names=required_cols, \n        windowSize=window_size,\n        highpass=band_pass.low_cutoff if (band_pass and band_pass.low_cutoff) else None,\n        lowpass=band_pass.high_cutoff if (band_pass and band_pass.high_cutoff) else None,\n        normalize=normalize,\n        ICA=ica,\n    )\n    return helper, dropped_df;","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:38:05.303119Z","iopub.execute_input":"2024-02-02T17:38:05.303489Z","iopub.status.idle":"2024-02-02T17:38:05.320553Z","shell.execute_reply.started":"2024-02-02T17:38:05.303459Z","shell.execute_reply":"2024-02-02T17:38:05.318940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert a eeglib EEG object back into a standard pandas dataframe\ndef eeg_to_df(\n    eeg: EEG, \n    eeg_channels:List[str], \n    restore_df:pd.DataFrame|None, \n    columns:List[str]=DIFF_EEG_COLUMNS\n):\n    df = pd.DataFrame(eeg.window.window.transpose(), columns=eeg_channels, copy=True)\n    if (restore_df is not None):\n        df = df.join(restore_df.reset_index())\n    return df[columns]","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:38:05.321810Z","iopub.execute_input":"2024-02-02T17:38:05.322571Z","iopub.status.idle":"2024-02-02T17:38:05.336943Z","shell.execute_reply.started":"2024-02-02T17:38:05.322540Z","shell.execute_reply":"2024-02-02T17:38:05.336084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"helper, aux_df = df_to_eeglib_helper(\n    extracted_diff_eeg,\n    band_pass=DEFAULT_BANDPASS,\n    normalize=True,\n)\n# Now we need to get an EEG object\n# Note, looks like a bug in eeglib, can only iterate once so we'll collect the eegs in a list.\n# Must make deep copies or will simply use underlying data for last iterator.\neegs = [deepcopy(eeg) for eeg in helper]\neeg = eegs[2] # Take centre 10 second window\nstart_time = (\n    sample_info['eeg_label_offset_seconds'] + \n    (EEG_SNAPSHOT_DURATION - EEG_LABEL_DURATION) / 2.0\n)\nend_time = start_time + EEG_LABEL_DURATION\neeg_aux_df = aux_df[(aux_df[TIME_OFFSET_COLUMN] >= start_time) & (aux_df[TIME_OFFSET_COLUMN] < end_time)]\neeglib_modified_df = eeg_to_df(\n    eeg, \n    helper.names, \n    eeg_aux_df\n)\neeg_plot(eeglib_modified_df, TIME_OFFSET_COLUMN, DIFF_PLOT_GROUPS, sample_info['eeg_label_offset_seconds'])","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:38:05.337916Z","iopub.execute_input":"2024-02-02T17:38:05.338383Z","iopub.status.idle":"2024-02-02T17:38:08.735742Z","shell.execute_reply.started":"2024-02-02T17:38:05.338355Z","shell.execute_reply":"2024-02-02T17:38:08.734555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### eeglib data input summary\n\nSo, we've ingressed some of the spectrogram data and it's been pre-processed (basically filtered and normalized).\n\n> Note that I can't get the ICA option to work, ignoring this!\n\nIt looks reasonably sensible although there are definite edge effects (looking at the training edge caused by the\nband-pass filter.","metadata":{}},{"cell_type":"markdown","source":"## eeglib features\n\nSo, why might we use eeglib?? The main reason appears to be that it can decompose\nthe time sequence eeg into a set of features. \nThere is a list in the [eeglib features documentation](https://eeglib.readthedocs.io/en/latest/features.html)\n\nThere _may_ be some visual correlations discernable between the features it produces\nand the expert votes in the training data.\n\nSo how to do this:\n\n1. Read in an entire eeg file into an EEGlib helper.\n2. Window this for every two seconds in the file data\n3. Generate a feature set from the window, per channel.\n4. Combine the results to produce a timed-based set of features.\n","metadata":{}},{"cell_type":"code","source":"# Read in the entire eeg file\nsample_eeg_df = make_differential(load_eeg(sample_info['eeg_id']))\nsample_eeg_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:47.611274Z","iopub.execute_input":"2024-02-02T17:34:47.611778Z","iopub.status.idle":"2024-02-02T17:34:47.679859Z","shell.execute_reply.started":"2024-02-02T17:34:47.611746Z","shell.execute_reply":"2024-02-02T17:34:47.678528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature generation\n\nDon't really know much about any of these feature types but we'll generate the single channel\nfeatures available in eeglib for the entire sample sliced into one second chunks","metadata":{}},{"cell_type":"code","source":"def extract_eeg_features(\n    helper: Helper,\n    scored_entries: pd.Series,\n    window_duration:float=EEG_LABEL_DURATION\n) -> pd.DataFrame:\n    # First we'll just build the features int python lists\n    bp_alpha = []\n    bp_beta = []\n    bp_delta = []\n    bp_theta = []\n    pfd = []\n    hfd = []\n    hjorth_activity = []\n    hjorth_mobility = []\n    hjorth_complexity = []\n    samp_en = []\n    lzc = []\n    dfa = []\n    num_channels = len(helper.names)\n    for eeg, is_scored in zip(helper, scored_entries):\n        if is_scored:\n            bp = eeg.bandPower()\n            # Split band power and convert to db\n            bp_alpha.append([10.0 * math.log10(ch[\"alpha\"]) for ch in bp])\n            bp_beta.append([10.0 * math.log10(ch[\"beta\"]) for ch in bp])\n            bp_delta.append([10.0 * math.log10(ch[\"delta\"]) for ch in bp])\n            bp_theta.append([10.0 * math.log10(ch[\"theta\"]) for ch in bp])\n            pfd.append(eeg.PFD())\n            hfd.append(eeg.HFD())\n            hjorth_activity.append(eeg.hjorthActivity())\n            hjorth_mobility.append(eeg.hjorthMobility())\n            hjorth_complexity.append(eeg.hjorthComplexity())\n            samp_en.append(eeg.sampEn())\n            lzc.append(eeg.LZC())\n            dfa.append(eeg.DFA())\n        else:\n            bp_alpha.append([math.nan] * num_channels)\n            bp_beta.append([math.nan] * num_channels)\n            bp_delta.append([math.nan] * num_channels)\n            bp_theta.append([math.nan] * num_channels)\n            pfd.append([math.nan] * num_channels)\n            hfd.append([math.nan] * num_channels)\n            hjorth_activity.append([math.nan] * num_channels)\n            hjorth_mobility.append([math.nan] * num_channels)\n            hjorth_complexity.append([math.nan] * num_channels)\n            samp_en.append([math.nan] * num_channels)\n            lzc.append([math.nan] * num_channels)\n            dfa.append([math.nan] * num_channels)\n    # Now make this into a dataframe\n    df_as_dict = {}\n    #print(bp_alpha)\n    for chan_idx, col_name in enumerate(helper.names):\n        df_as_dict[f\"{col_name}.bp_alpha\"] = [feature_values[chan_idx] for feature_values in bp_alpha]\n        df_as_dict[f\"{col_name}.bp_beta\"] = [feature_values[chan_idx] for feature_values in bp_beta]\n        df_as_dict[f\"{col_name}.bp_delta\"] = [feature_values[chan_idx] for feature_values in bp_delta]\n        df_as_dict[f\"{col_name}.bp_theta\"] = [feature_values[chan_idx] for feature_values in bp_theta]\n        df_as_dict[f\"{col_name}.pfd\"] = [feature_values[chan_idx] for feature_values in pfd]\n        df_as_dict[f\"{col_name}.hfd\"] = [feature_values[chan_idx] for feature_values in hfd]\n        df_as_dict[f\"{col_name}.hjorth_activity\"] = [feature_values[chan_idx] for feature_values in hjorth_activity]\n        df_as_dict[f\"{col_name}.hjorth_mobility\"] = [feature_values[chan_idx] for feature_values in hjorth_mobility]\n        df_as_dict[f\"{col_name}.hjorth_complexity\"] = [feature_values[chan_idx] for feature_values in hjorth_complexity]\n        df_as_dict[f\"{col_name}.samp_en\"] = [feature_values[chan_idx] for feature_values in samp_en]\n        df_as_dict[f\"{col_name}.lzc\"] = [feature_values[chan_idx] for feature_values in lzc]\n        df_as_dict[f\"{col_name}.dfa\"] = [feature_values[chan_idx] for feature_values in dfa]\n    feature_df = pd.DataFrame.from_dict(df_as_dict)    \n    # Set index as time channel\n    feature_df.index *= window_duration\n    return feature_df.copy() # This is because pandas warns the dataframe is fragmented!\n","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:47.681496Z","iopub.execute_input":"2024-02-02T17:34:47.681965Z","iopub.status.idle":"2024-02-02T17:34:47.704033Z","shell.execute_reply.started":"2024-02-02T17:34:47.681921Z","shell.execute_reply":"2024-02-02T17:34:47.702383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Build the sample features dataframe using the eeglib helper object\nFEATURE_WINDOW_DURATION = 2.0\nfeature_helper, _ = df_to_eeglib_helper(\n    sample_eeg_df,\n    band_pass=DEFAULT_BANDPASS,\n    normalize=True,\n    window_size=int(FEATURE_WINDOW_DURATION * EEG_SAMPLING_RATE)\n)\nfeature_df = extract_eeg_features(\n    feature_helper, \n    sample_eeg_p_per_sec[\"is_scored\"], \n    FEATURE_WINDOW_DURATION\n)\nfeature_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:47.706128Z","iopub.execute_input":"2024-02-02T17:34:47.706569Z","iopub.status.idle":"2024-02-02T17:34:53.720400Z","shell.execute_reply.started":"2024-02-02T17:34:47.706526Z","shell.execute_reply":"2024-02-02T17:34:53.718760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# There's a lot of numbers here! Try boiling some of it down into mean values per region...\nEEG_REGIONS = [\n    (\"LL\", LL_EEG_CHANNELS),\n    (\"LP\", LP_EEG_CHANNELS),\n    (\"CC\", CC_EEG_CHANNELS),\n    (\"RP\", RP_EEG_CHANNELS),\n    (\"RL\", RL_EEG_CHANNELS),\n]\nFEATURE_NAMES = [\n    \"bp_alpha\", \n    \"bp_beta\", \n    \"bp_delta\", \n    \"bp_theta\", \n    \"pfd\", \n    \"hfd\", \n    \"hjorth_activity\",\n    \"hjorth_mobility\", \n    \"hjorth_complexity\", \n    \"samp_en\", \n    \"lzc\", \n    \"dfa\"\n]\n\n# Function to add statistics channels per region\ndef add_feature_stats_by_region(df: pd.DataFrame) -> pd.DataFrame:\n    with_stats_df = df.copy()\n    for region_name, region_chans in EEG_REGIONS:\n        with_stats_df = with_stats_df.copy() # prevent DF fragmented warning\n        for feature_name in FEATURE_NAMES:\n            region_feature_chans = [f\"{chan}.{feature_name}\" for chan in region_chans]\n            with_stats_df[f\"{region_name}.{feature_name}.mean\"] = with_stats_df[region_feature_chans].mean(axis=1)\n    return with_stats_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.721921Z","iopub.execute_input":"2024-02-02T17:34:53.722377Z","iopub.status.idle":"2024-02-02T17:34:53.731490Z","shell.execute_reply.started":"2024-02-02T17:34:53.722342Z","shell.execute_reply":"2024-02-02T17:34:53.729991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features_and_region_stats_df = add_feature_stats_by_region(feature_df)\nfeatures_and_region_stats_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.732773Z","iopub.execute_input":"2024-02-02T17:34:53.733667Z","iopub.status.idle":"2024-02-02T17:34:53.872067Z","shell.execute_reply.started":"2024-02-02T17:34:53.733630Z","shell.execute_reply":"2024-02-02T17:34:53.870509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Joining the features data on probabilities for the eeg\n\nIn order to do some visualisation of the features we need to join the\nfeatures data with the probabilities for the same file data along the\ntime axis.","metadata":{}},{"cell_type":"code","source":"# Create a named function that does the joining\ndef join_features_and_probabilities(\n    features_df:pd.DataFrame, probabilities_df:pd.DataFrame\n) -> pd.DataFrame:\n    joined_df = features_df.join(probabilities_df)\n    return joined_df[joined_df[\"is_scored\"]].copy()\n\nlabelled_features = join_features_and_probabilities(\n    features_and_region_stats_df, sample_eeg_p_per_sec\n)\nlabelled_features","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.874296Z","iopub.execute_input":"2024-02-02T17:34:53.874711Z","iopub.status.idle":"2024-02-02T17:34:53.908871Z","shell.execute_reply.started":"2024-02-02T17:34:53.874676Z","shell.execute_reply":"2024-02-02T17:34:53.907342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature visualisation by brain region\n\nWe're first going to simply look at a nplot of each of the features from a single\nEEG file by the region of the brain that it emanated from.","metadata":{}},{"cell_type":"code","source":"REGION_MNEMONICS = [mne for mne, _ in EEG_REGIONS]\nREGION_MNEMONICS","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.910210Z","iopub.execute_input":"2024-02-02T17:34:53.910695Z","iopub.status.idle":"2024-02-02T17:34:53.917548Z","shell.execute_reply.started":"2024-02-02T17:34:53.910666Z","shell.execute_reply":"2024-02-02T17:34:53.915963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combine plot of scatterplot and density plot\ndef scatter_density_plot(\n    data:pd.DataFrame, x:str, y:str, color: str, ax:plt.Axes, xlim:Tuple[float,float]=(-0.1, 1.1)\n) -> None:\n    scatter_plot = sns.scatterplot(\n        data=data,\n        x=x, \n        y=y, \n        ax=ax,\n        color=color\n    ).set(xlim=xlim)\n#     CANNOT GET THIS TO WORK WITHOUT OCCASSIONAL DATA RELATED FAILURES\n#     levels = 5\n#     if data[y].count() >= levels:\n#         sns.kdeplot(\n#             data=data, \n#             x=x, \n#             y=y, \n#             levels=levels, \n#             fill=True, \n#             alpha=0.25, \n#             cut=2,\n#             ax=ax, \n#             warn_singular=False,\n#             color=color\n#         ).set(xlim=xlim)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.919346Z","iopub.execute_input":"2024-02-02T17:34:53.919681Z","iopub.status.idle":"2024-02-02T17:34:53.929125Z","shell.execute_reply.started":"2024-02-02T17:34:53.919653Z","shell.execute_reply":"2024-02-02T17:34:53.927926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot a scatterplot with density for each brain region \n# for a given probability channel\ndef regional_feature_probability_scatterplot(\n    data: pd.DataFrame, \n    feature_name:str,\n    probability_channels:List[str]=P_COLUMNS,\n    regions:List[str]=REGION_MNEMONICS,\n    fig_height:float=2.5,\n    fig_width:float=12,\n    force_color:str|None=None    \n) -> None:\n    figure, axes = plt.subplots(1, len(regions), sharey=True, sharex=True)\n    figure.subplots_adjust(wspace=0)\n    figure.set_figheight(fig_height)\n    figure.set_figwidth(fig_width)\n    visible_yaxis = True\n    my_color_cycler = plt.cycler(\n        'color',\n        plt.rcParams['axes.prop_cycle'].by_key()['color']\n    )\n    for region, ax in tqdm_notebook(\n        zip(regions, axes), \n        desc=\"Regions to plot\",\n        total=len(regions)\n    ):\n        feature_chan = f\"{region}.{feature_name}\"\n        props_cycle = my_color_cycler()\n        for probability_chan, props  in zip(probability_channels, props_cycle):\n            scatter_density_plot(\n                data=data, \n                x=probability_chan,\n                y=feature_chan,\n                ax=ax,\n                color=props['color'] if force_color is None else force_color\n            )\n        ax.title.set_text(region)\n        ax.yaxis.set_visible(visible_yaxis)\n        ax.set_xlabel(\"probability\")\n        ax.set_ylabel(None)\n        visible_yaxis = False\n    figure.supylabel(feature_name)\n    props_cycle = my_color_cycler()\n    legend_handles = [\n        mpatches.Patch(\n            color=props['color'] if force_color is None else force_color,\n            label=probability_chan\n        )\n        for probability_chan, props  in zip(probability_channels, props_cycle)\n    ]\n    figure.legend(handles=legend_handles)\n    figure.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:09:36.333469Z","iopub.execute_input":"2024-02-02T18:09:36.333749Z","iopub.status.idle":"2024-02-02T18:09:36.351484Z","shell.execute_reply.started":"2024-02-02T18:09:36.333727Z","shell.execute_reply":"2024-02-02T18:09:36.350378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regional_feature_probability_scatterplot(labelled_features, \"bp_alpha.mean\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:53.949157Z","iopub.execute_input":"2024-02-02T17:34:53.949573Z","iopub.status.idle":"2024-02-02T17:34:54.895000Z","shell.execute_reply.started":"2024-02-02T17:34:53.949538Z","shell.execute_reply":"2024-02-02T17:34:54.893660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for feature_name in tqdm_notebook(FEATURE_NAMES):\n    regional_feature_probability_scatterplot(labelled_features, f\"{feature_name}.mean\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:34:54.896774Z","iopub.execute_input":"2024-02-02T17:34:54.897272Z","iopub.status.idle":"2024-02-02T17:35:06.697770Z","shell.execute_reply.started":"2024-02-02T17:34:54.897230Z","shell.execute_reply":"2024-02-02T17:35:06.696454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Spatial Visualisation\n\nLet's try and spatially visualise the feature data by producing a grid of 14 scatterplots ordered as the probes are placed on the patient's head.","metadata":{}},{"cell_type":"code","source":"CC_EEG_CHANNELS_PADDED = [None] + CC_EEG_CHANNELS + [None]\nSPATIAL_CHANNELS = [\n    [ll, lp, cc, rp, rl] \n    for ll, lp, cc, rp, rl in zip(\n        LL_EEG_CHANNELS, LP_EEG_CHANNELS, CC_EEG_CHANNELS_PADDED, RP_EEG_CHANNELS, RL_EEG_CHANNELS\n    )\n]\nSPATIAL_CHANNELS","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:35:06.700193Z","iopub.execute_input":"2024-02-02T17:35:06.700931Z","iopub.status.idle":"2024-02-02T17:35:06.711503Z","shell.execute_reply.started":"2024-02-02T17:35:06.700894Z","shell.execute_reply":"2024-02-02T17:35:06.709670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def spatial_feature_probability_scatterplot(\n    data: pd.DataFrame,\n    feature_name:str,\n    probability_channels:List[str]=P_COLUMNS,\n    spatial_channels:List[List[str]]=SPATIAL_CHANNELS,\n    regions:List[str]=REGION_MNEMONICS,\n    fig_height:float=8,\n    fig_width:float=12,\n    force_color:str|None=None\n):\n    common_xlim = (-0.1, 1.1)\n    subplot_rows = len(spatial_channels)\n    subplot_cols = max([len(row) for row in spatial_channels])\n    figure, axes_2d = plt.subplots(subplot_rows, subplot_cols, sharey=True, sharex=True)\n    figure.subplots_adjust(wspace=0, hspace=0)\n    figure.set_figheight(fig_height)\n    figure.set_figwidth(fig_width)\n    my_color_cycler = plt.cycler(\n        'color',\n        plt.rcParams['axes.prop_cycle'].by_key()['color']\n    )\n    for row_idx, (channels, axes) in tqdm_notebook(\n        enumerate(zip(spatial_channels, axes_2d)),\n        total = subplot_rows,\n        desc=\"Front to back\"\n    ):\n        num_chans = len(channels)\n        visible_yaxis = True\n        first_row = 0 == row_idx\n        last_row = subplot_rows == row_idx+1\n        for col_idx, (region, ax) in tqdm_notebook(\n            enumerate(zip(regions, axes)),\n            total=subplot_cols,\n            desc=\"Left to right\"\n        ):\n            channel = channels[col_idx] if col_idx < num_chans else None\n            if channel is not None:\n                props_cycle = my_color_cycler()\n                for probability_channel, props in zip(probability_channels, props_cycle):\n                    scatter_density_plot(\n                        data=data, \n                        x=probability_channel, \n                        y=f\"{channel}.{feature_name}\", \n                        ax=ax,\n                        xlim=common_xlim,\n                        color=props[\"color\"] if force_color is None else force_color\n                    )\n            else:\n                ax.xaxis.set_visible(last_row)\n                ax.set_facecolor(wc.CSS3_NAMES_TO_HEX[\"gainsboro\"])\n                ax.set(xlim=common_xlim)\n            if first_row:\n                ax.title.set_text(region)\n            ax.set_xlabel(None if not last_row else \"probability\")\n            ax.set_ylabel(None)\n            ax.yaxis.set_visible(visible_yaxis)\n            visible_yaxis = False\n    figure.supylabel(feature_name)\n    props_cycle = my_color_cycler()\n    legend_handles = [\n        mpatches.Patch(\n            color=props['color'] if force_color is None else force_color,\n            label=probability_chan\n        )\n        for probability_chan, props  in zip(probability_channels, props_cycle)\n    ]\n    figure.legend(handles=legend_handles)\n    figure.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:09:36.289344Z","iopub.execute_input":"2024-02-02T18:09:36.289788Z","iopub.status.idle":"2024-02-02T18:09:36.307890Z","shell.execute_reply.started":"2024-02-02T18:09:36.289752Z","shell.execute_reply":"2024-02-02T18:09:36.306465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spatial_feature_probability_scatterplot(labelled_features, \"bp_alpha\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:35:06.732064Z","iopub.execute_input":"2024-02-02T17:35:06.732395Z","iopub.status.idle":"2024-02-02T17:35:11.172351Z","shell.execute_reply.started":"2024-02-02T17:35:06.732362Z","shell.execute_reply":"2024-02-02T17:35:11.170863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for feature_name in tqdm_notebook(FEATURE_NAMES, desc=\"Feature loop\"):\n    spatial_feature_probability_scatterplot(labelled_features, feature_name)","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:35:11.173960Z","iopub.execute_input":"2024-02-02T17:35:11.174351Z","iopub.status.idle":"2024-02-02T17:36:06.313093Z","shell.execute_reply.started":"2024-02-02T17:35:11.174317Z","shell.execute_reply":"2024-02-02T17:36:06.311312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Discussion\n\nSo can we see any correlation?\n\nI'm not sure however what we are looking at is statistically insignificant - \nwe're only looking at a single EEG sequence from a single patient.\n\nThe next step would be to feature-ise a larger proportion of the training dataset\nto see if the eeglib features begin to show correlations in the visualisations.\n","metadata":{}},{"cell_type":"markdown","source":"# Building a statistically larger the features dataset\n\nSo lets go ahead and build a larger dataset using the eeglib feature extraction\nbased on the above code.\n","metadata":{}},{"cell_type":"code","source":"# Function that combines all of the previous sub-processes to buid a\n# labelled features dataframe for a single eeg file\ndef featurize_eeg(train_info: pd.DataFrame, eeg_id: int) -> pd.DataFrame:\n    # Get temporal probabilities for this eeg\n    probabilities_df = eeg_P_per_sec(train_info, eeg_id)\n    # Load the EEG waveform data\n    eeg_df = load_eeg(eeg_id)\n    # Add a time channel to it\n    add_time_channel(eeg_df)\n    # Build the differential waveform\n    eeg_df = make_differential(eeg_df)\n    # Make an eeglib helper object\n    feature_helper, _ = df_to_eeglib_helper(\n        eeg_df,\n        band_pass=DEFAULT_BANDPASS,\n        normalize=True,\n        window_size=int(FEATURE_WINDOW_DURATION * EEG_SAMPLING_RATE)\n    )\n    # Extract the raw channel features from the helper\n    features_df = extract_eeg_features(feature_helper, probabilities_df[\"is_scored\"], FEATURE_WINDOW_DURATION)\n    # Add region statistics as well\n    features_df = add_feature_stats_by_region(features_df)\n    # And join the features and probabilities\n    features_df = join_features_and_probabilities(features_df, probabilities_df)\n    return features_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:36:06.314450Z","iopub.execute_input":"2024-02-02T17:36:06.314840Z","iopub.status.idle":"2024-02-02T17:36:06.323064Z","shell.execute_reply.started":"2024-02-02T17:36:06.314809Z","shell.execute_reply":"2024-02-02T17:36:06.321474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Just test the function works\ntest_eeg_df = featurize_eeg(train_info_clean, sample_info.eeg_id)\nspatial_feature_probability_scatterplot(test_eeg_df, \"dfa\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:36:06.325166Z","iopub.execute_input":"2024-02-02T17:36:06.325568Z","iopub.status.idle":"2024-02-02T17:36:11.965749Z","shell.execute_reply.started":"2024-02-02T17:36:06.325529Z","shell.execute_reply":"2024-02-02T17:36:11.964727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_eeg_ids = train_info_clean[\"eeg_id\"].unique()\nprint(f\"Total number of EEG files: {len(all_eeg_ids)}\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T17:36:11.967293Z","iopub.execute_input":"2024-02-02T17:36:11.968612Z","iopub.status.idle":"2024-02-02T17:36:11.976872Z","shell.execute_reply.started":"2024-02-02T17:36:11.968565Z","shell.execute_reply":"2024-02-02T17:36:11.976031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set the random seed\nRANDOM_SEED = 0x55aa6699\n\n# Sample percentage of the dataset for visualisation\nSAMPLE_FRAC_PC = 5.0 ","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:20:33.435369Z","iopub.execute_input":"2024-02-02T18:20:33.435806Z","iopub.status.idle":"2024-02-02T18:20:33.441455Z","shell.execute_reply.started":"2024-02-02T18:20:33.435772Z","shell.execute_reply":"2024-02-02T18:20:33.440372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NUM_SAMPLES = int(len(all_eeg_ids) * SAMPLE_FRAC_PC / 100.0)\nrandom_generator = np.random.default_rng(seed=RANDOM_SEED)\nrandom_eeg_ids = random_generator.choice(all_eeg_ids, size=NUM_SAMPLES, replace=False).tolist()\nprint(f\"Total number of sampled EEG files: {len(random_eeg_ids)}\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:20:37.113459Z","iopub.execute_input":"2024-02-02T18:20:37.113881Z","iopub.status.idle":"2024-02-02T18:20:37.122981Z","shell.execute_reply.started":"2024-02-02T18:20:37.113847Z","shell.execute_reply":"2024-02-02T18:20:37.121590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Some flags to control dataset generation\nGENERATE_DATASET = True \nITS_A_TEST_RUN = False\nTEST_RUN_SLICE_SIZE = 20\n\nsampled_eegs = random_eeg_ids\nif not GENERATE_DATASET or (GENERATE_DATASET and ITS_A_TEST_RUN):\n    # Just get a small slice of them\n    sampled_eegs = random_eeg_ids[:TEST_RUN_SLICE_SIZE]\nif GENERATE_DATASET:\n    print(\n        f\"{'TEST MODE: ' if ITS_A_TEST_RUN else ''}\"\n        f\"Will build dataset dataframe using {len(sampled_eegs)} EEG files.\"\n    )\nelse:\n    print(\"Dataset generation is disabled.\")","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:20:43.449250Z","iopub.execute_input":"2024-02-02T18:20:43.449685Z","iopub.status.idle":"2024-02-02T18:20:43.456237Z","shell.execute_reply.started":"2024-02-02T18:20:43.449651Z","shell.execute_reply":"2024-02-02T18:20:43.455347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now generate the dataset dataframes\noutput_dataset_df = None\nfor dataset_eeg_id in tqdm_notebook(sampled_eegs, desc=\"Files processed\"):\n    tqdm.write(f\"EEG id: {dataset_eeg_id}\")\n    dataset_eeg_df = featurize_eeg(train_info_clean, dataset_eeg_id)\n    if output_dataset_df is None:\n        output_dataset_df = dataset_eeg_df\n    else:\n        output_dataset_df = pd.concat(\n            [output_dataset_df, dataset_eeg_df], ignore_index=True, axis=0\n        )\noutput_dataset_df","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:20:54.368201Z","iopub.execute_input":"2024-02-02T18:20:54.369091Z","iopub.status.idle":"2024-02-02T18:44:34.055558Z","shell.execute_reply.started":"2024-02-02T18:20:54.369023Z","shell.execute_reply":"2024-02-02T18:44:34.053704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"props_cycler = plt.cycler(\n    'color',\n    plt.rcParams['axes.prop_cycle'].by_key()['color']\n)\nfor feature_name in tqdm_notebook(FEATURE_NAMES, desc=\"Feature loop\"):\n    props_cycle = props_cycler();\n    for probability_channel, props in zip(P_COLUMNS, props_cycle):\n        spatial_feature_probability_scatterplot(\n            output_dataset_df, \n            feature_name, \n            [probability_channel],\n            force_color=props['color']\n        )","metadata":{"execution":{"iopub.status.busy":"2024-02-02T18:45:02.905283Z","iopub.execute_input":"2024-02-02T18:45:02.905923Z","iopub.status.idle":"2024-02-02T18:54:32.149295Z","shell.execute_reply.started":"2024-02-02T18:45:02.905888Z","shell.execute_reply":"2024-02-02T18:54:32.147300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusion\n\nSo, I'm struggling to see much of a visual correlation between the eeglib\ngenerated features and the various probability levels.\n\nWe're looking for, essentially diagonal relationships in these graphs and there\nare not a lot of these and they're a bit unclear.\n\nThere's also obvious issues in the features, presumably as a result of\ninadequate pre-processing or the standardization process applied when we made the\ndifferential channels. These are significantly noticible in the band power\nfeatures.\n\nBut... this might form a basis for input to train a model...\n\nMany thanks,\nAndrew.","metadata":{}}]}