{"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":"markdown","source":"# Introduction of features commonly used for EEG signal analysis","metadata":{}},{"cell_type":"markdown","source":"In this notebook, I will introduce the feature for analysis of EEG signal.\n\nThe features I will introduce are:\n\n1. **Permutation entropy**\n2. **Spectral entropy**\n3. **Singular value decomposition entropy**\n4. **Hjorth mobility and complexity**\n5. **Number of zero-crossings**\n6. **Petrosian fractal dimension**\n7. **Katz fractal dimension**\n8. **Higuchi fractal dimension**\n9. **Detrended fluctuation analysis**\n\n\nIn the last cell I visualized these features using UMAP\n\nIf you want more information, please visit [documentation of antropy package](https://raphaelvallat.com/antropy/build/html/index.html)","metadata":{}},{"cell_type":"markdown","source":"### Please Upvote if you Find this Useful :)","metadata":{}},{"cell_type":"markdown","source":"# antropy is used for feature calculation","metadata":{}},{"cell_type":"code","source":"!pip install antropy","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-01-16T14:22:42.585109Z","iopub.execute_input":"2024-01-16T14:22:42.586014Z","iopub.status.idle":"2024-01-16T14:23:16.666083Z","shell.execute_reply.started":"2024-01-16T14:22:42.585968Z","shell.execute_reply":"2024-01-16T14:23:16.664896Z"},"_kg_hide-input":false,"jupyter":{"outputs_hidden":true},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport antropy as ant\nimport pywt\nfrom tqdm.notebook import tqdm\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:16.668538Z","iopub.execute_input":"2024-01-16T14:23:16.668955Z","iopub.status.idle":"2024-01-16T14:23:28.168831Z","shell.execute_reply.started":"2024-01-16T14:23:16.668918Z","shell.execute_reply":"2024-01-16T14:23:28.167562Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/train.csv\")\nTARGETS = df.columns[-6:]\nprint('Train shape:', df.shape )\nprint('Targets', list(TARGETS))\ndf.head()\ntrain = df.groupby('eeg_id')[['spectrogram_id','spectrogram_label_offset_seconds']].agg(\n    {'spectrogram_id':'first','spectrogram_label_offset_seconds':'min'})\ntrain.columns = ['spec_id','spectrogram_min']\n\ntmp = df.groupby('eeg_id')[['spectrogram_id','spectrogram_label_offset_seconds']].agg(\n    {'spectrogram_label_offset_seconds':'max'})\ntrain['spectrogram_max'] = tmp\n\ntmp = df.groupby(\"eeg_id\")[\"eeg_label_offset_seconds\"].agg(\"min\")\ntrain[\"eeg_min\"] = tmp\ntmp = df.groupby(\"eeg_id\")[\"eeg_label_offset_seconds\"].agg(\"max\")\ntrain[\"eeg_max\"] = tmp\n\n\ntmp = df.groupby('eeg_id')[['patient_id']].agg('first')\ntrain['patient_id'] = tmp\n\ntmp = df.groupby('eeg_id')[TARGETS].agg('sum')\nfor t in TARGETS:\n    train[t] = tmp[t].values\n\ny_data = train[TARGETS].values\ny_data = y_data / y_data.sum(axis=1,keepdims=True)\ntrain[TARGETS] = y_data\n\ntmp = df.groupby('eeg_id')[['expert_consensus']].agg('first')\ntrain['target'] = tmp\n\ntrain = train.reset_index()\nprint('Train non-overlapp eeg_id shape:', train.shape )\ntrain","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:28.171350Z","iopub.execute_input":"2024-01-16T14:23:28.172481Z","iopub.status.idle":"2024-01-16T14:23:28.721630Z","shell.execute_reply.started":"2024-01-16T14:23:28.172421Z","shell.execute_reply":"2024-01-16T14:23:28.720368Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true,"outputs_hidden":true},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train[\"ideal\"] = False\ntrain[\"ideal\"] = train.apply((lambda row: row[str.lower(row[\"target\"])+\"_vote\"]==1.0), axis=1)\ntrain[\"ideal\"].value_counts()\ntrain_ideal = train.query(\"ideal==1\")\n\ntarget = [\"seizure_vote\", \"lpd_vote\", \"gpd_vote\", \"lrda_vote\", \"grda_vote\", \"other_vote\"]\n\neegs = {key:[] for key in target}\n\nfor t in target:\n    \n    train_ideal_target = train_ideal[train_ideal[t]==1.0]\n    for j in range(10):\n        eegs[t] = train_ideal_target[\"eeg_id\"].sample(n=50, random_state=42).to_list()\n        \ntrain = train.set_index(\"eeg_id\")","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:28.724672Z","iopub.execute_input":"2024-01-16T14:23:28.725106Z","iopub.status.idle":"2024-01-16T14:23:29.092552Z","shell.execute_reply.started":"2024-01-16T14:23:28.725071Z","shell.execute_reply":"2024-01-16T14:23:29.091264Z"},"jupyter":{"source_hidden":true,"outputs_hidden":true},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def maddest(d, axis=None):\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef denoise(x, wavelet='haar', level=1):\n    ret = {key:[] for key in x.columns}\n    \n    for pos in x.columns:\n        coeff = pywt.wavedec(x[pos], wavelet, mode=\"per\")\n        sigma = (1/0.6745) * maddest(coeff[-level])\n\n        uthresh = sigma * np.sqrt(2*np.log(len(x)))\n        coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])\n\n        ret[pos]=pywt.waverec(coeff, wavelet, mode='per')\n    \n    return pd.DataFrame(ret)\n\nBR_ORDER = [\"Fp1-F7\", \"F7-T3\",   \"T3-T5\",   \"T5-O1\",\n            \"Fp1-F3\",  \"F3-C3\",   \"C3-P3\",   \"P3-O1\",\n            \"Fp2-F4\",  \"F4-C4\",  \"C4-P4\",   \"P4-O2\",\n            \"Fp2-F8\",  \"F8-T4\" ,  \"T4-T6\",   \"T6-O2\",\n            \"Fz-Cz\",   \"Cz-Pz\"]\nTARGET = [\"seizure_vote\", \"lpd_vote\", \"gpd_vote\", \"lrda_vote\", \"grda_vote\", \"other_vote\"]\n\ndef get_bipolar_referenced_eeg(eeg):\n\n    ret = {key:[] for key in BR_ORDER}\n\n    for bipos in BR_ORDER:\n        pos_a, pos_b = bipos.split(\"-\")[0], bipos.split(\"-\")[1]\n        ret[bipos] = eeg[pos_a] - eeg[pos_b]\n    \n    return pd.DataFrame(ret)\n\ndef get_denoised_eeg_signal(eeg_id):\n    path = \"/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/\"\n    eeg = pd.read_parquet(path + str(eeg_id) + \".parquet\")\n    eeg = denoise(eeg)\n    eeg = get_bipolar_referenced_eeg(eeg)\n    return eeg\n\ndef apply_function_to_eeg(eeg, features, feat_name, func, args):\n\n    for pos in eeg.columns:\n        try:\n            feat = func(eeg[pos], **args)\n        except:\n            feat = np.nan\n            \n        \n        if feat_name==\"hjorth_entropy\":\n            features[feat_name+\"_\"+\"m\"+\"_\"+pos].append(feat[0])\n            features[feat_name+\"_\"+\"c\"+\"_\"+pos].append(feat[1])\n        else:\n            features[feat_name+\"_\"+pos].append(feat)\n    return features\n\ndef get_feature(eegs, df, feat_name, func, args):\n    \n    if feat_name == \"hjorth_entropy\":\n        feat_brorder = [feat_name+\"_\"+f+\"_\"+pos for f in [\"m\",\"c\"] for pos in BR_ORDER]\n    else:\n        feat_brorder = [feat_name+\"_\"+pos for pos in BR_ORDER]\n    \n    features = {key:[] for key in feat_brorder}\n    features[\"eeg_id\"] = []\n    features[\"target\"] = []\n    \n    for t in target:\n        eeg_list = eegs[t]\n        print(f\"Get {t} feature\")\n        for eeg_id in tqdm(eeg_list):\n            features[\"eeg_id\"].append(eeg_id)\n            features[\"target\"].append(t)\n            \n            eeg = get_denoised_eeg_signal(eeg_id)\n            features = apply_function_to_eeg(eeg, features, feat_name, func, args)\n    \n    features = pd.DataFrame(features)\n    \n    features = features.fillna(features.median(numeric_only=True))\n    \n    return pd.DataFrame(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T15:30:07.749165Z","iopub.execute_input":"2024-01-16T15:30:07.749626Z","iopub.status.idle":"2024-01-16T15:30:07.775092Z","shell.execute_reply.started":"2024-01-16T15:30:07.749594Z","shell.execute_reply":"2024-01-16T15:30:07.773445Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_features = []","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:29.119605Z","iopub.execute_input":"2024-01-16T14:23:29.120692Z","iopub.status.idle":"2024-01-16T14:23:29.131547Z","shell.execute_reply.started":"2024-01-16T14:23:29.120648Z","shell.execute_reply":"2024-01-16T14:23:29.130536Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Permutaion entropy\n ","metadata":{}},{"cell_type":"code","source":"# Permutation entropy\nfeatures =  get_feature(eegs, df=train, feat_name=\"Permutaion_entropy\", func=ant.perm_entropy, args={\"normalize\":1})\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:29.133199Z","iopub.execute_input":"2024-01-16T14:23:29.134450Z","iopub.status.idle":"2024-01-16T14:23:55.664457Z","shell.execute_reply.started":"2024-01-16T14:23:29.134399Z","shell.execute_reply":"2024-01-16T14:23:55.663281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Spectral entropy","metadata":{}},{"cell_type":"code","source":"# Spectral entropy\nfeatures =  get_feature(eegs, df=train, feat_name=\"Spectral_entropy\", func=ant.spectral_entropy, args=dict(sf=100, method='welch', normalize=True))\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:23:55.666308Z","iopub.execute_input":"2024-01-16T14:23:55.666925Z","iopub.status.idle":"2024-01-16T14:24:13.973894Z","shell.execute_reply.started":"2024-01-16T14:23:55.666867Z","shell.execute_reply":"2024-01-16T14:24:13.972595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Singular value decomposition entropy","metadata":{}},{"cell_type":"code","source":"# Singular value decomposition entropy\nfeatures =  get_feature(eegs, df=train, feat_name=\"SVD_entropy\", func=ant.svd_entropy, args=dict(normalize=True))\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:24:13.975424Z","iopub.execute_input":"2024-01-16T14:24:13.975815Z","iopub.status.idle":"2024-01-16T14:24:30.968316Z","shell.execute_reply.started":"2024-01-16T14:24:13.975768Z","shell.execute_reply":"2024-01-16T14:24:30.966707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Hjorth mobility and complexity","metadata":{}},{"cell_type":"code","source":"# Hjorth mobility and complexity\nfeatures =  get_feature(eegs, df=train, feat_name=\"hjorth_entropy\", func=ant.hjorth_params, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T15:30:13.402808Z","iopub.execute_input":"2024-01-16T15:30:13.403327Z","iopub.status.idle":"2024-01-16T15:30:28.296776Z","shell.execute_reply.started":"2024-01-16T15:30:13.403286Z","shell.execute_reply":"2024-01-16T15:30:28.295563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Number of zero-crossings","metadata":{}},{"cell_type":"code","source":"# Number of zero-crossings\nfeatures =  get_feature(eegs, df=train, feat_name=\"zerocross\", func=ant.num_zerocross, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:24:30.981128Z","iopub.execute_input":"2024-01-16T14:24:30.982353Z","iopub.status.idle":"2024-01-16T14:24:44.825687Z","shell.execute_reply.started":"2024-01-16T14:24:30.982307Z","shell.execute_reply":"2024-01-16T14:24:44.824600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Petrosian fractal dimension","metadata":{}},{"cell_type":"code","source":"# Petrosian fractal dimension\nfeatures =  get_feature(eegs, df=train, feat_name=\"petrosian_fd\", func=ant.petrosian_fd, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:24:44.827770Z","iopub.execute_input":"2024-01-16T14:24:44.828481Z","iopub.status.idle":"2024-01-16T14:24:58.485418Z","shell.execute_reply.started":"2024-01-16T14:24:44.828420Z","shell.execute_reply":"2024-01-16T14:24:58.484200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Katz fractal dimension","metadata":{}},{"cell_type":"code","source":"# Katz fractal dimension\nfeatures =  get_feature(eegs, df=train, feat_name=\"katz_fd\", func=ant.katz_fd, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:24:58.486992Z","iopub.execute_input":"2024-01-16T14:24:58.487510Z","iopub.status.idle":"2024-01-16T14:25:12.848456Z","shell.execute_reply.started":"2024-01-16T14:24:58.487464Z","shell.execute_reply":"2024-01-16T14:25:12.847473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Higuchi fractal dimension","metadata":{}},{"cell_type":"code","source":"# Higuchi fractal dimension\nfeatures =  get_feature(eegs, df=train, feat_name=\"higuchi_fd\", func=ant.higuchi_fd, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:25:12.849960Z","iopub.execute_input":"2024-01-16T14:25:12.851221Z","iopub.status.idle":"2024-01-16T14:25:27.540861Z","shell.execute_reply.started":"2024-01-16T14:25:12.851171Z","shell.execute_reply":"2024-01-16T14:25:27.539415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"# Detrended fluctuation analysis","metadata":{}},{"cell_type":"code","source":"# Detrended fluctuation analysis\nfeatures =  get_feature(eegs, df=train, feat_name=\"detrended_fluction\", func=ant.detrended_fluctuation, args=dict())\nall_features.append(features)","metadata":{"execution":{"iopub.status.busy":"2024-01-16T14:25:27.542828Z","iopub.execute_input":"2024-01-16T14:25:27.543985Z","iopub.status.idle":"2024-01-16T14:26:11.807480Z","shell.execute_reply.started":"2024-01-16T14:25:27.543935Z","shell.execute_reply":"2024-01-16T14:26:11.806120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization using UMAP","metadata":{}},{"cell_type":"code","source":"all_feature = pd.concat(all_features, axis=1).drop([\"eeg_id\", \"target\"], axis=1)\nall_feature.head()","metadata":{"execution":{"iopub.status.busy":"2024-01-16T15:30:38.491221Z","iopub.execute_input":"2024-01-16T15:30:38.492216Z","iopub.status.idle":"2024-01-16T15:30:38.533004Z","shell.execute_reply.started":"2024-01-16T15:30:38.492162Z","shell.execute_reply":"2024-01-16T15:30:38.531763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import umap\n\n\nmapper = umap.UMAP(random_state=42,\n                   n_neighbors=128,\n                   min_dist=0.99,\n                   metric=\"mahalanobis\")\nembedding = mapper.fit_transform(all_feature)\n\nembedding_x = embedding[:,0]\nembedding_y = embedding[:,1]\n\nc_list = [\"r\", \"b\", \"y\", \"g\", \"m\", \"c\"]\nfor i, t in enumerate(TARGET):\n    plt.scatter(embedding_x[50*i:50*(i+1)], embedding_y[50*i:50*(i+1)], c=c_list[i], label=t)\nplt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-16T15:49:56.854615Z","iopub.execute_input":"2024-01-16T15:49:56.855148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}