{"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"}],"dockerImageVersionId":30635,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport pathlib\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\ndset_path = pathlib.PurePath(\"/kaggle/input/hms-harmful-brain-activity-classification/\")\ntrain_eegs = dset_path/\"train_eegs\"\ntrain_specs = dset_path/\"train_spectrograms\"","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-01-19T19:01:56.472572Z","iopub.execute_input":"2024-01-19T19:01:56.473286Z","iopub.status.idle":"2024-01-19T19:01:56.864942Z","shell.execute_reply.started":"2024-01-19T19:01:56.473249Z","shell.execute_reply":"2024-01-19T19:01:56.864165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta = pd.read_csv(dset_path/\"train.csv\", dtype={\"eeg_label_offset_seconds\": \"Int64\",\n                                                       \"spectrogram_label_offset_seconds\": \"Int64\",\n                                                       \"expert_consensus\": \"category\"})\ntrain_meta.info()","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:56.866420Z","iopub.execute_input":"2024-01-19T19:01:56.866989Z","iopub.status.idle":"2024-01-19T19:01:57.281220Z","shell.execute_reply.started":"2024-01-19T19:01:56.866960Z","shell.execute_reply":"2024-01-19T19:01:57.279997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def paired_eeg(eeg_id, offset):\n    consolidated_eeg = pd.read_parquet(train_eegs/f\"{eeg_id}.parquet\")\n    # 200 rows = 1 second\n    # Want 50 seconds starting at the offset\n    start = offset * 200\n    end = start + (200 * 50)\n    return consolidated_eeg.iloc[start:end,]\n\ndef paired_spectrogram(spec_id, offset):\n    consolidated_spec = pd.read_parquet(train_specs/f\"{spec_id}.parquet\")\n    start = offset\n    end = offset + 600\n    return consolidated_spec[(consolidated_spec[\"time\"] <= end) & (consolidated_spec[\"time\"] >= start)]","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:57.282960Z","iopub.execute_input":"2024-01-19T19:01:57.283422Z","iopub.status.idle":"2024-01-19T19:01:57.291292Z","shell.execute_reply.started":"2024-01-19T19:01:57.283383Z","shell.execute_reply":"2024-01-19T19:01:57.290180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"row = train_meta.sample(n=1).iloc[0,]\nrow","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:57.294757Z","iopub.execute_input":"2024-01-19T19:01:57.295581Z","iopub.status.idle":"2024-01-19T19:01:57.313171Z","shell.execute_reply.started":"2024-01-19T19:01:57.295540Z","shell.execute_reply":"2024-01-19T19:01:57.312005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting\n## EEG","metadata":{}},{"cell_type":"code","source":"eeg_df = paired_eeg(row[\"eeg_id\"], row[\"eeg_label_offset_seconds\"])\neeg_df","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:57.314352Z","iopub.execute_input":"2024-01-19T19:01:57.314648Z","iopub.status.idle":"2024-01-19T19:01:57.509777Z","shell.execute_reply.started":"2024-01-19T19:01:57.314621Z","shell.execute_reply":"2024-01-19T19:01:57.508763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The eeg dataframe contains 50 seconds, to somewhat replicate the plots in example_figures we should subset to the central 10 seconds.","metadata":{}},{"cell_type":"code","source":"def central_window_eeg(df, window_seconds=10):\n    start = df.index[0]\n    end = df.index[-1]\n    mid = (start + end) / 2\n    new_start = int(mid - window_seconds/2 * 200) - start + 1\n    new_end = int(mid + window_seconds/2 * 200) - start + 1\n    return df.iloc[new_start:new_end]\n    \ncentral_window_eeg(eeg_df, 10)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:57.511252Z","iopub.execute_input":"2024-01-19T19:01:57.511649Z","iopub.status.idle":"2024-01-19T19:01:57.557658Z","shell.execute_reply.started":"2024-01-19T19:01:57.511611Z","shell.execute_reply":"2024-01-19T19:01:57.556640Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_eeg(df, moving_avg=1):\n    fig, axs = plt.subplots(20, 1, figsize=(15, 10), sharex=True)\n    for i, ax in enumerate(axs):\n        ax.plot(df.iloc[:,i], color=\"black\")\n        ax.set_ylabel(df.columns[i], rotation=0)\n        ax.set_yticklabels([])\n        ax.set_yticks([])\n        ax.set_xticks([])\n        ax.spines[[\"top\", \"bottom\", \"left\", \"right\"]].set_visible(False)\n\nplot_eeg(central_window_eeg(eeg_df, 10))","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:57.558771Z","iopub.execute_input":"2024-01-19T19:01:57.559057Z","iopub.status.idle":"2024-01-19T19:01:58.921911Z","shell.execute_reply.started":"2024-01-19T19:01:57.559024Z","shell.execute_reply":"2024-01-19T19:01:58.921177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is a lot of noise in this plot, even when only plotting 10 seconds, especially for EKG.","metadata":{}},{"cell_type":"code","source":"_ = plt.specgram(eeg_df[\"EKG\"], Fs=200)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:58.922828Z","iopub.execute_input":"2024-01-19T19:01:58.923132Z","iopub.status.idle":"2024-01-19T19:01:59.254155Z","shell.execute_reply.started":"2024-01-19T19:01:58.923091Z","shell.execute_reply":"2024-01-19T19:01:59.253140Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looking at a spectrogram of the data the noise seems to be mostly at 60Hz and 0Hz. https://www.gehealthcare.co.uk/insights/article/a-guide-to-ecg-signal-filtering confirms that this is common for this kind of data.\n\nIt seems to be relatively common to filter eeg data using a bandpass filter between 1-40 or 1-30 Hz.","metadata":{}},{"cell_type":"code","source":"from scipy import signal\n\nfs = 200\n# b_notch, a_notch = signal.iirnotch(60, 30.0, fs)\n# b_lowpass, a_lowpass = signal.butter(4, 10, btype=\"lowpass\", fs=200)\nsos = signal.butter(4, [1, 40], btype=\"band\", fs=200, output=\"sos\")\n\ndef filt(data):\n#     data = signal.filtfilt(b_notch, a_notch, data)\n#     return signal.filtfilt(b_lowpass, a_lowpass, data)\n    return signal.sosfiltfilt(sos, data)\n\nto_filter = eeg_df[\"EKG\"]\nfig, axs = plt.subplots(1, 2, figsize=(15, 8), sharey=True)\naxs[0].plot(to_filter)\naxs[0].set_title(\"Original signal\")\naxs[1].plot(filt(to_filter))\naxs[1].set_title(\"Filtered\")","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:01:59.255531Z","iopub.execute_input":"2024-01-19T19:01:59.255838Z","iopub.status.idle":"2024-01-19T19:02:00.964778Z","shell.execute_reply.started":"2024-01-19T19:01:59.255810Z","shell.execute_reply":"2024-01-19T19:02:00.963637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_eeg(df):\n    fig, axs = plt.subplots(20, 1, figsize=(15, 10), sharex=True)\n    for i, ax in enumerate(axs):\n        ax.plot(filt(df.iloc[:,i]), color=\"black\")\n        ax.set_ylabel(df.columns[i], rotation=0)\n        ax.set_yticklabels([])\n        ax.set_yticks([])\n        ax.set_xticks([])\n        ax.spines[[\"top\", \"bottom\", \"left\", \"right\"]].set_visible(False)\n\nplot_eeg(central_window_eeg(eeg_df, 10))","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:00.968970Z","iopub.execute_input":"2024-01-19T19:02:00.969306Z","iopub.status.idle":"2024-01-19T19:02:02.322135Z","shell.execute_reply.started":"2024-01-19T19:02:00.969276Z","shell.execute_reply":"2024-01-19T19:02:02.321039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Spectrogram","metadata":{}},{"cell_type":"code","source":"spec = paired_spectrogram(row[\"spectrogram_id\"], row[\"spectrogram_label_offset_seconds\"])\nspec","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:02.323902Z","iopub.execute_input":"2024-01-19T19:02:02.324423Z","iopub.status.idle":"2024-01-19T19:02:02.404047Z","shell.execute_reply.started":"2024-01-19T19:02:02.324382Z","shell.execute_reply":"2024-01-19T19:02:02.403008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"~~This dataframe contains 10 minutes of data, for the purpose of replicating the plots we only need the central 10 seconds, though for making decisions we probably want to use a wider window.~~\n\nNevermind, I misread the plot, the plot in the examples is using the whole 10 minutes, but the central window function may still be useful.","metadata":{}},{"cell_type":"code","source":"def central_window_spec(df, window_seconds=10):\n    start = df.iloc[0, 0]\n    end = df.iloc[-1, 0]\n    mid = (start + end) / 2\n    new_start = mid - window_seconds/2\n    new_end = mid + window_seconds/2\n    return df[(df[\"time\"] >= new_start)].iloc[:window_seconds//2]\n\ncentral_window_spec(spec, 30)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:02.405549Z","iopub.execute_input":"2024-01-19T19:02:02.406207Z","iopub.status.idle":"2024-01-19T19:02:02.452646Z","shell.execute_reply.started":"2024-01-19T19:02:02.406165Z","shell.execute_reply":"2024-01-19T19:02:02.451758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy import signal\n\ndef plot_spec(df):\n    fig, axs = plt.subplots(4, 1, figsize=(15, 10), sharey=True)\n\n    prefixes = ['LL', 'RL', 'LP', 'RP']\n\n    for ax, prefix in zip(axs, prefixes):\n        cols = df.filter(regex=f\"^{prefix}_\").columns\n        ax.imshow(spec[cols].T, origin=\"lower\", norm=\"log\", cmap=\"plasma\", interpolation=\"none\")\n        ax.set_title(prefix)\n        ax.set_yticks(np.arange(0, 101, 25.))\n        ax.set_yticklabels([0, 5, 10, 15, 20])\n        ax.set_ylabel(\"Freq\")\n        ax.set_xticks(np.arange(0, 301, 75.))\n        ax.set_xticklabels(range(spec.iloc[0, 0], spec.iloc[-1, 0], (spec.iloc[-1, 0]-spec.iloc[0, 0])//4))\n        ax.set_xlabel(\"Seconds\")\n\n    plt.tight_layout()\n    \nplot_spec(spec)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:02.454150Z","iopub.execute_input":"2024-01-19T19:02:02.454831Z","iopub.status.idle":"2024-01-19T19:02:03.299161Z","shell.execute_reply.started":"2024-01-19T19:02:02.454793Z","shell.execute_reply":"2024-01-19T19:02:03.297965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some EDA","metadata":{}},{"cell_type":"code","source":"train_meta[\"expert_consensus\"].value_counts().plot(kind=\"bar\")","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:03.300755Z","iopub.execute_input":"2024-01-19T19:02:03.301439Z","iopub.status.idle":"2024-01-19T19:02:03.567275Z","shell.execute_reply.started":"2024-01-19T19:02:03.301380Z","shell.execute_reply":"2024-01-19T19:02:03.566136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta[\"patient_id\"].value_counts().head(30).plot(kind=\"bar\")","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:03.568373Z","iopub.execute_input":"2024-01-19T19:02:03.568667Z","iopub.status.idle":"2024-01-19T19:02:03.973748Z","shell.execute_reply.started":"2024-01-19T19:02:03.568642Z","shell.execute_reply":"2024-01-19T19:02:03.972595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_meta[\"patient_id\"].value_counts().mean())\nprint(train_meta[\"patient_id\"].value_counts().median())","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:03.975073Z","iopub.execute_input":"2024-01-19T19:02:03.975420Z","iopub.status.idle":"2024-01-19T19:02:03.988551Z","shell.execute_reply.started":"2024-01-19T19:02:03.975391Z","shell.execute_reply":"2024-01-19T19:02:03.987493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The number of examples for each patient seems to vary a lot. Potentially when learning should use less of the most frequent patients. I am assuming that the test/leaderboard data will include new patients, so trying to improve performance on unseen patients could be useful.","metadata":{}},{"cell_type":"code","source":"top_ids = train_meta[\"patient_id\"].value_counts().head(50).index\nmask = ~train_meta[\"patient_id\"].isin(top_ids)\nfiltered_train_meta = train_meta[mask]\nfiltered_train_meta[\"expert_consensus\"].value_counts().plot(kind=\"bar\")","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:04:15.595877Z","iopub.execute_input":"2024-01-19T19:04:15.596287Z","iopub.status.idle":"2024-01-19T19:04:15.843153Z","shell.execute_reply.started":"2024-01-19T19:04:15.596254Z","shell.execute_reply":"2024-01-19T19:04:15.842190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"However simply removing the most frequently occuring patient IDs (or downsampling) makes the imbalance in expert concensus a bigger issue.","metadata":{}},{"cell_type":"markdown","source":"## Spectrogram sums","metadata":{}},{"cell_type":"code","source":"sample = 500\ndef plot_summed_specs(df, sample):\n    acc = paired_spectrogram(df.iloc[0][\"spectrogram_id\"], df.iloc[0][\"spectrogram_label_offset_seconds\"])\n    for col in acc.columns:\n        acc[col].values[:] = 0\n    for i, row in df.sample(n=sample).iterrows():\n        spec = paired_spectrogram(row[\"spectrogram_id\"], row[\"spectrogram_label_offset_seconds\"]).fillna(0)\n#         normalised_spec = (spec - spec.min()) / (spec.max() - spec.min())\n        acc = acc.add(spec, fill_value=0)\n    plot_spec(acc)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:03.990883Z","iopub.execute_input":"2024-01-19T19:02:03.991301Z","iopub.status.idle":"2024-01-19T19:02:04.000119Z","shell.execute_reply.started":"2024-01-19T19:02:03.991264Z","shell.execute_reply":"2024-01-19T19:02:03.999060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"Seizure\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:04.001252Z","iopub.execute_input":"2024-01-19T19:02:04.001762Z","iopub.status.idle":"2024-01-19T19:02:31.275570Z","shell.execute_reply.started":"2024-01-19T19:02:04.001728Z","shell.execute_reply":"2024-01-19T19:02:31.274367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"LPD\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:02:31.277043Z","iopub.execute_input":"2024-01-19T19:02:31.277950Z","iopub.status.idle":"2024-01-19T19:03:02.880070Z","shell.execute_reply.started":"2024-01-19T19:02:31.277904Z","shell.execute_reply":"2024-01-19T19:03:02.878950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"GPD\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:03:02.881532Z","iopub.execute_input":"2024-01-19T19:03:02.882504Z","iopub.status.idle":"2024-01-19T19:03:15.082016Z","shell.execute_reply.started":"2024-01-19T19:03:02.882461Z","shell.execute_reply":"2024-01-19T19:03:15.080508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"LRDA\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:03:15.083251Z","iopub.status.idle":"2024-01-19T19:03:15.083894Z","shell.execute_reply.started":"2024-01-19T19:03:15.083699Z","shell.execute_reply":"2024-01-19T19:03:15.083720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"GRDA\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:03:15.085034Z","iopub.status.idle":"2024-01-19T19:03:15.085642Z","shell.execute_reply.started":"2024-01-19T19:03:15.085458Z","shell.execute_reply":"2024-01-19T19:03:15.085478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected = train_meta[train_meta[\"expert_consensus\"] == \"Other\"]\nplot_summed_specs(selected, sample)","metadata":{"execution":{"iopub.status.busy":"2024-01-19T19:03:15.086709Z","iopub.status.idle":"2024-01-19T19:03:15.087304Z","shell.execute_reply.started":"2024-01-19T19:03:15.087120Z","shell.execute_reply":"2024-01-19T19:03:15.087140Z"},"trusted":true},"execution_count":null,"outputs":[]}]}