{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","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":30648,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os, gc\nimport pandas as pd, numpy as np\nfrom glob import glob\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm\nimport pickle","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:27:31.181123Z","iopub.execute_input":"2024-03-04T06:27:31.181515Z","iopub.status.idle":"2024-03-04T06:27:32.395613Z","shell.execute_reply.started":"2024-03-04T06:27:31.181478Z","shell.execute_reply":"2024-03-04T06:27:32.393978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# intro\n* Hi! In this notebook, I will demonstrate how to extract signal features from EEG signals and perform basic exploratory data analysis (EDA).\n* The notebook focuses on obtaining 9 key signal features: 'peak_to_peak', 'excess_kurtosis', 'crest_factor', 'zero_cross_ratio', 'shape_factor', 'impulse_factor', 'clearance_factor', 'skewness', and 'rms'.\n* Additionally, basic EDA will be conducted, including the creation of a box plot for visualization.\n* If you're interested in exploring other signal features or enhancing your EDA, this notebook serves as a helpful guide.\n* I recommend creating multimodal analyses using these features with spectrograms for a more comprehensive understanding.","metadata":{}},{"cell_type":"markdown","source":"## Create dataset\n* We will load all unique EEG signals into memory, slice EEG sub-data, and calculate features for each EEG signal.\n* For this task, we will utilize the Numba JIT compiler, as calculating features for 100,000 EEG sub-signals, each with 20 signal channels and 9 features, can be time-consuming.\n* Numba significantly accelerates this process, making it 10 times faster.","metadata":{}},{"cell_type":"markdown","source":"#### load data","metadata":{}},{"cell_type":"code","source":"# load meta data\ndataset_path = '/kaggle/input/hms-harmful-brain-activity-classification'\ntraining_meta = pd.read_csv(dataset_path + '/train.csv')\ntraining_meta.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:27:32.397719Z","iopub.execute_input":"2024-03-04T06:27:32.398465Z","iopub.status.idle":"2024-03-04T06:27:32.728305Z","shell.execute_reply.started":"2024-03-04T06:27:32.398431Z","shell.execute_reply":"2024-03-04T06:27:32.727079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unique_eeg_id_list = training_meta.groupby(\"eeg_id\").sum().index.tolist() # get unique id","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:27:32.729557Z","iopub.execute_input":"2024-03-04T06:27:32.729881Z","iopub.status.idle":"2024-03-04T06:27:32.804243Z","shell.execute_reply.started":"2024-03-04T06:27:32.729855Z","shell.execute_reply":"2024-03-04T06:27:32.802859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# load all eeg signal into memory\neeg_arr = {}\nfor i,eeg_id in tqdm(enumerate(unique_eeg_id_list)):\n    egg_data = pd.read_parquet(dataset_path+f'/train_eegs/{eeg_id}.parquet')\n    eeg_arr[eeg_id] = egg_data\n    \n    ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-03-04T06:27:32.806939Z","iopub.execute_input":"2024-03-04T06:27:32.807290Z","iopub.status.idle":"2024-03-04T06:35:49.273155Z","shell.execute_reply.started":"2024-03-04T06:27:32.807261Z","shell.execute_reply":"2024-03-04T06:35:49.271585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### get features","metadata":{}},{"cell_type":"code","source":"feature_name = ['peak_to_peak', 'kurtosis', 'crest_factor', 'zero_cross_ratio', 'shape_factor', 'impulse_factor', 'clearance_factor', 'skewness', 'rms']\nfeatures_n = 9\nsignal_n = 20\nsampling_rate = 200","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:35:49.275458Z","iopub.execute_input":"2024-03-04T06:35:49.275811Z","iopub.status.idle":"2024-03-04T06:35:49.286735Z","shell.execute_reply.started":"2024-03-04T06:35:49.275783Z","shell.execute_reply":"2024-03-04T06:35:49.285286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from numba import jit\n\n@jit(nopython=True)\ndef get_features(signal: np.ndarray):\n    \n    n = len(signal)\n    mean = np.mean(signal)\n    std = np.std(signal)  \n    peak = np.max(np.abs(signal))\n    abs_signal = np.abs(signal)\n    squared_signal = signal**2\n    zero_crossings = np.where(np.diff(np.sign(signal)))[0]\n\n    # get features\n    rms = np.float32(np.sqrt(np.sum(squared_signal) / n))\n    skewness = (np.sum((signal - mean)**3) / n) / (std**3)\n    clearance_factor = peak/ np.float32((np.sum(np.sqrt(abs_signal)) / n)**2)\n    impulse_factor = np.float32(peak / (np.sum(abs_signal) / n))\n    shape_factor = np.float32(rms / (np.sum(abs_signal) / n))\n    zero_cross_ratio = np.float32(len(zero_crossings) / (n - 1))\n    crest_factor = np.float32(peak / rms)\n    kurtosis = np.float32(np.sum((signal - mean) ** 4) / (n * std ** 4))\n    peak_to_peak = np.float32(np.max(signal) - np.min(signal))\n    features = [peak_to_peak, kurtosis, crest_factor, zero_cross_ratio, shape_factor, impulse_factor, clearance_factor, skewness, rms]\n\n    return features","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:35:49.288435Z","iopub.execute_input":"2024-03-04T06:35:49.288830Z","iopub.status.idle":"2024-03-04T06:35:51.008988Z","shell.execute_reply.started":"2024-03-04T06:35:49.288797Z","shell.execute_reply":"2024-03-04T06:35:51.007735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfeatures = np.zeros([len(training_meta), features_n * signal_n])\nnoise = np.random.normal(0, 0.001, 10000) # for zero divide error\nfor ind in tqdm(training_meta.index):\n    \n    eeg_id = training_meta.loc[ind,'eeg_id']\n    eeg_label_offset_seconds = training_meta.loc[ind,'eeg_label_offset_seconds']\n    \n    egg_data = eeg_arr[eeg_id]\n\n    start_ind_sub_data = int(eeg_label_offset_seconds * sampling_rate)\n    end_ind_sub_data = int((eeg_label_offset_seconds + 50) * sampling_rate)\n    eeg_sub_data = egg_data[start_ind_sub_data: end_ind_sub_data]\n    \n    columns = eeg_sub_data.columns\n    eeg_signals = eeg_sub_data.values\n    \n    for i, colname in enumerate(columns):\n        feature = get_features(eeg_signals[:,i] + noise)\n        features[ind, i*features_n : (i+1)*features_n] = feature\n        \nfor i, colname in enumerate(columns):\n    feature_col_names = [feature + \"_\" + colname for feature in feature_name]\n    training_meta.loc[:,feature_col_names] = features[:, i*features_n : (i+1)*features_n]","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:35:51.010629Z","iopub.execute_input":"2024-03-04T06:35:51.011022Z","iopub.status.idle":"2024-03-04T06:47:34.604975Z","shell.execute_reply.started":"2024-03-04T06:35:51.010990Z","shell.execute_reply":"2024-03-04T06:47:34.603618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### basic eda using boxplot","metadata":{}},{"cell_type":"code","source":"# feature_name\nfeatures_col_list = training_meta.columns.tolist()[15:]\n\nplt.figure(figsize=[40,120])\nfor i in range(signal_n):\n    for j in range(features_n):\n        feature_name_i = features_col_list[j + features_n * i]\n        sub_data = training_meta.loc[:,['expert_consensus']+[feature_name_i]]\n        plt.subplot(signal_n,features_n,j + features_n * i + 1)\n#         plt.title(feature_name_i)\n        sns.boxplot(x='expert_consensus', y=feature_name_i, data=sub_data, showfliers=False) #do not plot outlier \n    ","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:47:34.606722Z","iopub.execute_input":"2024-03-04T06:47:34.607436Z","iopub.status.idle":"2024-03-04T06:48:39.968453Z","shell.execute_reply.started":"2024-03-04T06:47:34.607400Z","shell.execute_reply":"2024-03-04T06:48:39.967436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### train classifier","metadata":{}},{"cell_type":"code","source":"training_meta.loc[:,['expert_consensus']+features_col_list]","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:48:39.971753Z","iopub.execute_input":"2024-03-04T06:48:39.972234Z","iopub.status.idle":"2024-03-04T06:48:40.270498Z","shell.execute_reply.started":"2024-03-04T06:48:39.972195Z","shell.execute_reply":"2024-03-04T06:48:40.269158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def encode_label(target_arr):\n    encode_dic = {\n        'Seizure' : 0,\n        'LPD' : 1,\n        'GPD' : 2,\n        'LRDA' : 3,\n        'GRDA': 4,\n        'Other' : 5,\n    }\n    encoded_label_arr = np.array([encode_dic[label[0]] for label in target_arr])\n    return encoded_label_arr.reshape(-1,1)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:48:40.275124Z","iopub.execute_input":"2024-03-04T06:48:40.275546Z","iopub.status.idle":"2024-03-04T06:48:40.282966Z","shell.execute_reply.started":"2024-03-04T06:48:40.275513Z","shell.execute_reply":"2024-03-04T06:48:40.281464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train set\nfrom sklearn.model_selection import train_test_split\n\nX = training_meta.loc[:,features_col_list].values\ny = encode_label(training_meta.loc[:,['expert_consensus']].values)\nX_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, random_state=7)\nprint(X_train.shape, X_valid.shape, y_train.shape, y_valid.shape)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:48:40.284656Z","iopub.execute_input":"2024-03-04T06:48:40.285197Z","iopub.status.idle":"2024-03-04T06:48:41.548275Z","shell.execute_reply.started":"2024-03-04T06:48:40.285163Z","shell.execute_reply":"2024-03-04T06:48:41.547008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from xgboost import XGBClassifier\n\nmodel = XGBClassifier(objective='multi:softprob',eval_metric=['merror','mlogloss'])\nmodel.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=True)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:48:41.549640Z","iopub.execute_input":"2024-03-04T06:48:41.550996Z","iopub.status.idle":"2024-03-04T06:49:45.438743Z","shell.execute_reply.started":"2024-03-04T06:48:41.550952Z","shell.execute_reply":"2024-03-04T06:49:45.437392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"feature importance(SHAP)","metadata":{}},{"cell_type":"code","source":"import shap","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:49:45.440654Z","iopub.execute_input":"2024-03-04T06:49:45.441110Z","iopub.status.idle":"2024-03-04T06:49:51.718291Z","shell.execute_reply.started":"2024-03-04T06:49:45.441075Z","shell.execute_reply":"2024-03-04T06:49:51.716715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nexplainer = shap.TreeExplainer(model)\nshap_values = explainer(X_train)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T06:49:51.719944Z","iopub.execute_input":"2024-03-04T06:49:51.720606Z","iopub.status.idle":"2024-03-04T06:58:26.061329Z","shell.execute_reply.started":"2024-03-04T06:49:51.720570Z","shell.execute_reply":"2024-03-04T06:58:26.059946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'Seizure' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:12:25.024453Z","iopub.execute_input":"2024-03-04T07:12:25.024984Z","iopub.status.idle":"2024-03-04T07:12:46.008835Z","shell.execute_reply.started":"2024-03-04T07:12:25.024949Z","shell.execute_reply":"2024-03-04T07:12:46.007636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'LPD' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:12:50.165424Z","iopub.execute_input":"2024-03-04T07:12:50.166558Z","iopub.status.idle":"2024-03-04T07:13:11.306398Z","shell.execute_reply.started":"2024-03-04T07:12:50.166492Z","shell.execute_reply":"2024-03-04T07:13:11.304935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'GPD' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:13:11.308638Z","iopub.execute_input":"2024-03-04T07:13:11.309050Z","iopub.status.idle":"2024-03-04T07:13:32.106255Z","shell.execute_reply.started":"2024-03-04T07:13:11.309016Z","shell.execute_reply":"2024-03-04T07:13:32.104923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'LRDA' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:13:32.108139Z","iopub.execute_input":"2024-03-04T07:13:32.108692Z","iopub.status.idle":"2024-03-04T07:13:52.861552Z","shell.execute_reply.started":"2024-03-04T07:13:32.108644Z","shell.execute_reply":"2024-03-04T07:13:52.860204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'GRDA' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:13:52.864261Z","iopub.execute_input":"2024-03-04T07:13:52.864753Z","iopub.status.idle":"2024-03-04T07:14:13.674807Z","shell.execute_reply.started":"2024-03-04T07:13:52.864709Z","shell.execute_reply":"2024-03-04T07:14:13.673401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" #'Other' class feature importance\nshap.summary_plot(shap_values[:,:,0], X_train, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:14:13.676545Z","iopub.execute_input":"2024-03-04T07:14:13.676941Z","iopub.status.idle":"2024-03-04T07:14:34.416070Z","shell.execute_reply.started":"2024-03-04T07:14:13.676910Z","shell.execute_reply":"2024-03-04T07:14:34.414768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### prediction","metadata":{}},{"cell_type":"code","source":"# test set\ntest_meta = pd.read_csv(dataset_path + '/test.csv')","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:22.881310Z","iopub.execute_input":"2024-03-04T07:11:22.881776Z","iopub.status.idle":"2024-03-04T07:11:22.898349Z","shell.execute_reply.started":"2024-03-04T07:11:22.881744Z","shell.execute_reply":"2024-03-04T07:11:22.896903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:23.071605Z","iopub.execute_input":"2024-03-04T07:11:23.072383Z","iopub.status.idle":"2024-03-04T07:11:23.083621Z","shell.execute_reply.started":"2024-03-04T07:11:23.072347Z","shell.execute_reply":"2024-03-04T07:11:23.082367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get test features\nfeatures = np.zeros([len(test_meta), features_n * signal_n])\nnoise = np.random.normal(0, 0.001, 10000) # for zero divide error\nfor ind in tqdm(test_meta.index):\n    \n    eeg_id = test_meta.loc[ind,'eeg_id']\n    \n    \n    egg_data = pd.read_parquet(dataset_path+f'/test_eegs/{eeg_id}.parquet')\n    columns = egg_data.columns\n    eeg_signals = egg_data.values\n    \n    for i, colname in enumerate(columns):\n        feature = get_features(eeg_signals[:,i] + noise)\n        features[ind, i*features_n : (i+1)*features_n] = feature\n        \nfor i, colname in enumerate(columns):\n    feature_col_names = [feature + \"_\" + colname for feature in feature_name]\n    test_meta.loc[:,feature_col_names] = features[:, i*features_n : (i+1)*features_n]","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:23.289430Z","iopub.execute_input":"2024-03-04T07:11:23.289891Z","iopub.status.idle":"2024-03-04T07:11:23.528357Z","shell.execute_reply.started":"2024-03-04T07:11:23.289828Z","shell.execute_reply":"2024-03-04T07:11:23.527048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:23.617249Z","iopub.execute_input":"2024-03-04T07:11:23.617699Z","iopub.status.idle":"2024-03-04T07:11:23.649953Z","shell.execute_reply.started":"2024-03-04T07:11:23.617664Z","shell.execute_reply":"2024-03-04T07:11:23.649034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test = test_meta.loc[:,features_col_list].values\nprint(X_test.shape)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:24.018894Z","iopub.execute_input":"2024-03-04T07:11:24.020063Z","iopub.status.idle":"2024-03-04T07:11:24.037820Z","shell.execute_reply.started":"2024-03-04T07:11:24.020011Z","shell.execute_reply":"2024-03-04T07:11:24.035715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_pred_proba = np.float32(model.predict_proba(X_test))\nprint(y_pred_proba.shape)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:24.506007Z","iopub.execute_input":"2024-03-04T07:11:24.507177Z","iopub.status.idle":"2024-03-04T07:11:24.522474Z","shell.execute_reply.started":"2024-03-04T07:11:24.507135Z","shell.execute_reply":"2024-03-04T07:11:24.518925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']\nsub = pd.DataFrame({'eeg_id': test_meta.eeg_id.values})\nsub[TARGETS] = y_pred_proba\nsub.iloc[0,1:] = sub.iloc[0,1:].values\nsub.to_csv(f'submission.csv',index=False)\nprint(f'Submission shape: {sub.shape}')\nsub.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:24.925483Z","iopub.execute_input":"2024-03-04T07:11:24.926800Z","iopub.status.idle":"2024-03-04T07:11:24.956199Z","shell.execute_reply.started":"2024-03-04T07:11:24.926759Z","shell.execute_reply":"2024-03-04T07:11:24.954875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.sum(sub.iloc[0,1:].values)","metadata":{"execution":{"iopub.status.busy":"2024-03-04T07:11:27.332595Z","iopub.execute_input":"2024-03-04T07:11:27.333074Z","iopub.status.idle":"2024-03-04T07:11:27.343480Z","shell.execute_reply.started":"2024-03-04T07:11:27.333038Z","shell.execute_reply":"2024-03-04T07:11:27.341935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}