{"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":7886478,"sourceType":"datasetVersion","datasetId":4629648},{"sourceId":7917940,"sourceType":"datasetVersion","datasetId":4652608},{"sourceId":7943199,"sourceType":"datasetVersion","datasetId":4670238}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Probabilistic Classification of Epileptic Patterns using the Discrete Wavelet Transform","metadata":{}},{"cell_type":"markdown","source":"Doctors interpret electroencephalogram (EEG) signals manually, for the diagnosis of harmful brain activity such as epilepsy. This is a laborious process and can lead to discrepancies in the diagnosis, even between experts.  This study is based on a Kaggle competition. It aims to estimate the probability that a 10-second EEG sample belongs to one of 6 epileptic patterns. Patterns of interest are labeled as seizure (**SZ**), generalized periodic discharges (**GPD**), lateralized periodic discharges (**LPD**), lateralized rhythmic delta activity (**LRDA**), generalized rhythmic delta activity (**GRDA**), and the category “**Other**”. The last one is a not-defined pattern.","metadata":{}},{"cell_type":"code","source":"import os\nimport os.path\n\nimport pandas as pd, numpy as np\nimport pyarrow.parquet as pq\nimport matplotlib.pyplot as plt\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import savgol_filter\nimport tensorflow as tf\nfrom pywt import wavedec\nimport concurrent\nimport pickle,gzip\nimport xgboost as xgb\nimport sklearn\nimport keras\nfrom keras import layers\nfrom keras import ops","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:19:54.483026Z","iopub.execute_input":"2025-04-07T23:19:54.483386Z","iopub.status.idle":"2025-04-07T23:20:13.366185Z","shell.execute_reply.started":"2025-04-07T23:19:54.483348Z","shell.execute_reply":"2025-04-07T23:20:13.365049Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_train=pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')\ndf_search=df_train[[\"eeg_id\",\"eeg_label_offset_seconds\"]]\ndf_search_np=df_search.to_numpy(dtype=int)\n\nEEG_Tag=df_train['eeg_id'].unique()\nEEG_Path='/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/'\nprint(len(df_train['eeg_id'].unique()))\nprint(len(df_train['patient_id'].unique()))","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:13.367051Z","iopub.execute_input":"2025-04-07T23:20:13.367622Z","iopub.status.idle":"2025-04-07T23:20:13.665766Z","shell.execute_reply.started":"2025-04-07T23:20:13.367587Z","shell.execute_reply":"2025-04-07T23:20:13.664591Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The metadata file includes the patient label, the electroencephalogram (EEG) file label, its corresponding spectrogram, the diagnostic consensus among experts, and the expert votes for each type of brain activity, among other data.\n\n* There are **17089** EEG files taken from **1950** patients. Each file contains a n number of subsamples of 10-second windows.","metadata":{}},{"cell_type":"code","source":"df_train.head(5)","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:13.666845Z","iopub.execute_input":"2025-04-07T23:20:13.667257Z","iopub.status.idle":"2025-04-07T23:20:13.696564Z","shell.execute_reply.started":"2025-04-07T23:20:13.667209Z","shell.execute_reply":"2025-04-07T23:20:13.695445Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"In this study, only EEG data will be analyzed. EEG files contain data sampled at a frequency of 200 Hz from 19 electrodes (Fp1, F3, C3, P3, F7, T3, T5, O1, Fz, Cz, Pz, Fp2, F4, C4, P4, F8, T4, T6, O2) located on the scalp and an electrode corresponding to an electrocardiogram signal, labeled EKG, which is not considered for the analyses.","metadata":{}},{"cell_type":"code","source":"eeg= pd.read_parquet('/kaggle/input/example/2289322082.parquet')\neeg.columns\n","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:13.697637Z","iopub.execute_input":"2025-04-07T23:20:13.697924Z","iopub.status.idle":"2025-04-07T23:20:13.867330Z","shell.execute_reply.started":"2025-04-07T23:20:13.697902Z","shell.execute_reply":"2025-04-07T23:20:13.866445Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Example of an electrode signal taken from one of the EEG files**","metadata":{}},{"cell_type":"code","source":"eeg_np=eeg.to_numpy()\nx = np.linspace(0,30, 6000)\nfig, axs = plt.subplots(1,1, figsize=(8,4))\naxs.plot(x,eeg_np[:,0][:6000])\naxs.set_title('Fp1')\naxs.axis([0, 30, -850, 120])\naxs.set_xlabel('Time [s]')\naxs.set_ylabel('Voltage [μV]')\n\nthreshold = -120\n\naxs.fill_between(x,25,-600,where=eeg_np[:,0][:6000]< threshold  , color='cyan', alpha=0.4, transform=axs.get_xaxis_transform())\n\n###Filtro","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:13.868397Z","iopub.execute_input":"2025-04-07T23:20:13.868697Z","iopub.status.idle":"2025-04-07T23:20:14.198703Z","shell.execute_reply.started":"2025-04-07T23:20:13.868659Z","shell.execute_reply":"2025-04-07T23:20:14.197625Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The previous graph shows the first 30 seconds of the signal taken from the Fp1 electrode. One of the difficulties in classifying EEG signals is that they can be distorted by noise and artifacts, one of the causes being sudden movements, dropped electrodes, and scalp sweat, among others. In the previous figure, a segment of the signal affected by artifacts is shown in cyan.","metadata":{}},{"cell_type":"code","source":"x = np.linspace(0,20, 4000)\nfig, axs = plt.subplots(1,1, figsize=(8,4))\nsos = signal.butter(2, 0.1, 'hp', fs=200, output='sos')\nfiltered = signal.sosfilt(sos, eeg_np[:,0][:4000])\n\n\nx_sg=savgol_filter(filtered, 17, 2)\naxs.plot(x, x_sg)\naxs.set_title('Filtered signal (Fp1)')\naxs.set_xlabel('Time [s]')\naxs.set_ylabel('Voltage [μV]')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:14.201170Z","iopub.execute_input":"2025-04-07T23:20:14.201461Z","iopub.status.idle":"2025-04-07T23:20:14.446288Z","shell.execute_reply.started":"2025-04-07T23:20:14.201439Z","shell.execute_reply":"2025-04-07T23:20:14.445085Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The figure shows the first 20 seconds of the signal after passing through one of the filtering stages.","metadata":{}},{"cell_type":"markdown","source":"# Method","metadata":{}},{"cell_type":"markdown","source":"Three main procedures were performed on the multichannel signal coming from the 19 electrodes, of each EEG sample. The first step was to filter the signal. To reduce the effect of the artifacts, a second-order Butterworth high-pass filter was used, with a cut-off frequency equal to 0.1 Hz. In addition, the Savitzky-Golay filter attenuated the noise, acting as a pass low-filter and is characterized by maintaining the spectral properties of the EEG signals, its parameters are the polynomial order (k) and the window length (n). Based on [1], this filter was designed with k=2 and n=17.\n\nThe next step decomposes the filtered signal to extract its most important features. EEG signals are non-stationary. Discrete Wavelet Transform (DWT) is one of the most efficient tools for signal analysis of this type [1].\n\nDWT decomposes a signal into approximation coefficients and detail coefficients. The approximation coefficients at each level are again decomposed into approximation and detail coefficients. Each level represents different frequency bands. The extraction of features from each frequency band can be used in automated seizure detection systems [2].\n\nThe number of decomposition levels and the wavelet type must be indicated to compute the DWT. Using the PyWavelet library [3], 4 levels were specified, and the Daubechies wavelet of order 4. According to [1], it is the most appropriate for  EEG signal analysis.\n\nFor applying the previously described algorithm to an electrode signal, it is represented by 5 frequency sub-bands given by their respective coefficients. According to the computational complexity and the type of classification task required, There are different approaches to determining the feature extraction method. In [4], an algorithm derived from DWT is used to classify 7 seizure signal patterns, and 6 statistical parameters are obtained for each sub-band.\n\nIt is proposed to mathematically represent each frequency band by 6 features, with [4] as a model, for the feature extraction procedure. With M being the length of the signal in each sub-band, the 6 features are given by:\n\n1. \tMean absolute values of the coefficients in each sub-band, μ.\n\n$$ μ={1 \\over M} \\sum\\limits_{j=1}^{M}|{y_i}|$$\n\n2. \tAverage power of the coefficients in each sub-band, λ. \n\n$$λ= \\sqrt{{1 \\over M} \\sum\\limits_{j=1}^{M} {y_i^2}}$$\n\n3. Standard deviation of the coefficients in each sub-band, σ.\n\n$$σ= \\sqrt{{1 \\over M} \\sum\\limits_{j=1}^{M} {(y_i-μ)^2}}$$\n\n4. \tSkewness (skew) of the coefficients in each sub-band, ϕ.\n\n$$ϕ= \\sqrt{{1 \\over M} {\\sum\\limits_{j=1}^{M} {{(y_i-μ)^3}\\over σ ^3}}}$$\n\n5. \tKurtosis of the coefficients in each sub-band, K.\n\n$$K= \\sqrt{{1 \\over M} {\\sum\\limits_{j=1}^{M} {{(y_i-μ)^4}\\over σ ^4}}}$$\n\n","metadata":{}},{"cell_type":"markdown","source":"After computing the DWT, only the last 4 sub-bands were chosen for the feature extraction, since the feature analysis showed that the first one was irrelevant. Each EEG signal is composed of 18 channels, so after this step, it is transformed into a matrix of dimensions 18x5x4, which is flattened to obtain a vector of dimensions 1x360for the classification task.","metadata":{}},{"cell_type":"markdown","source":"\n**Next, the code and graphs of the coefficients are shown, after applying the DWT with 4 levels of decomposition to the first 20 seconds of the signal from electrode Fp1.**","metadata":{}},{"cell_type":"code","source":"(cA4, cD4, cD3, cD2,cD1) = wavedec(x_sg, 'db4', level=4,mode='smooth')\nx1 = np.linspace(0,20, len(cD1))\nfig, axs= plt.subplots(5,1, figsize=(10,20))\naxs[0].plot(x1[:], cD1)\naxs[0].set_title('cD1 Coeficient')\n\nx2 = np.linspace(0,20, len(cD2))\naxs[1].plot(x2, cD2)\naxs[1].set_title('cD2 Coeficient')\n\nx3 = np.linspace(0,20, len(cD3))\naxs[2].plot(x3, cD3)\naxs[2].set_title('cD3 Coeficient')\n\nx4 = np.linspace(0,20, len(cD4))\naxs[3].plot(x4, cD4)\naxs[3].set_title('cD4 Coeficient')\n\nx5= np.linspace(0,20, len(cA4))\naxs[4].plot(x4, cA4)\naxs[4].set_title('cA4 Coeficient')\naxs[4].set_xlabel('Time [s]')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:14.448313Z","iopub.execute_input":"2025-04-07T23:20:14.448614Z","iopub.status.idle":"2025-04-07T23:20:15.459457Z","shell.execute_reply.started":"2025-04-07T23:20:14.448587Z","shell.execute_reply":"2025-04-07T23:20:15.458362Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The following code is used to filter an unsegmented EEG label, that is, first, the entire signal is filtered, and then other functions are applied to segment it and calculate the feature vector of the multi-channel EEG. Before filtering it, if there are nan values, they are replaced by zeros.","metadata":{}},{"cell_type":"code","source":"def Filter(x):\n    x_clean = np.where(np.isnan(x), np.median(x[~np.isnan(x)]), x)\n    threshold_high = np.percentile(x_clean, 99)\n    threshold_low = np.percentile(x_clean, 1)\n    x_clean = np.clip(x_clean, threshold_low, threshold_high)\n    \n    sos = signal.butter(2, 0.1, 'hp', fs=200, output='sos')\n    filt_1Hz = signal.sosfilt(sos, x_clean)\n    x_filtered = savgol_filter(filt_1Hz, 17, 2)\n    return x_filtered\n\n\ndef Filtering(M):\n    \n    eeg_np=M.fillna(0).to_numpy()\n    for i in range(19):\n        eeg_np[:,i]=Filter(eeg_np[:,i])\n    return eeg_np","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-07T23:20:15.460483Z","iopub.execute_input":"2025-04-07T23:20:15.460750Z","iopub.status.idle":"2025-04-07T23:20:15.467072Z","shell.execute_reply.started":"2025-04-07T23:20:15.460728Z","shell.execute_reply":"2025-04-07T23:20:15.465901Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The following code is used to apply the DWT to the filtered EEG signal and extract the feature vector.","metadata":{}},{"cell_type":"code","source":"def Stats_DWT(subband_x):\n      Mean=np.mean(np.absolute(subband_x))\n      AVP=scipy.stats.pmean(np.absolute(subband_x),2,weights=None)\n      SD=subband_x.std()\n      Skew=scipy.stats.skew(subband_x)\n      Kurt=scipy.stats.kurtosis(subband_x)\n\n      return Mean,AVP,SD, Skew,Kurt\n        \ndef DWT_Features(xf, include_cd1=False):\n    (cA4, cD4, cD3, cD2, cD1) = wavedec(xf, 'db4', level=4, mode='smooth')\n    f_cA4 = Stats_DWT(cA4)\n    f_cD4 = Stats_DWT(cD4)\n    f_cD3 = Stats_DWT(cD3)\n    f_cD2 = Stats_DWT(cD2)\n    features = np.concatenate((f_cA4, f_cD4, f_cD3, f_cD2), axis=None)\n    \n    if include_cd1:\n        f_cD1 = Stats_DWT(cD1)\n        features = np.concatenate((features, f_cD1), axis=None)\n    \n    return features","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.468165Z","iopub.execute_input":"2025-04-07T23:20:15.468541Z","iopub.status.idle":"2025-04-07T23:20:15.490556Z","shell.execute_reply.started":"2025-04-07T23:20:15.468505Z","shell.execute_reply":"2025-04-07T23:20:15.489541Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"After filtering each electrode signal, the Transverse Central Parietal (TCP) montage was used to accentuate spikes, a characteristic behavior of epileptic signals. Its setup consists of differencing the signals taken from 2 electrodes (e.g., Fp1-F7, F7-T3) [4].\n\nAfter this procedure, 18 differential signals are obtained, which will be used for feature extraction. The code below will be used for this step.","metadata":{}},{"cell_type":"code","source":"def single_channel(eeg_np):\n\n    \n    Fp1_f7=eeg_np[:,0]-eeg_np[:,4]  \n    F7_T3=eeg_np[:,4]-eeg_np[:,5]  \n    T3_T5=eeg_np[:,5]-eeg_np[:,6]\n    T5_O1=eeg_np[:,6]-eeg_np[:,7]\n\n    Fp2_F8=eeg_np[:,11]-eeg_np[:,15]  \n    F8_T4=eeg_np[:,15]-eeg_np[:,16]  \n    T4_T6=eeg_np[:,16]-eeg_np[:,17]\n    T6_O2=eeg_np[:,17]-eeg_np[:,18]\n\n    Fp1_F3=eeg_np[:,0]-eeg_np[:,1]\n    F3_C3=eeg_np[:,1]-eeg_np[:,2] \n    C3_P3=eeg_np[:,2]-eeg_np[:,3]  \n    P3_O1=eeg_np[:,3]-eeg_np[:,7]\n\n    Fp2_F4=eeg_np[:,11]-eeg_np[:,12]\n    F4_C4=eeg_np[:,12]-eeg_np[:,13]\n    C4_P4=eeg_np[:,13]-eeg_np[:,14] \n    P4_O2=eeg_np[:,14]-eeg_np[:,18]\n\n    Fz_Cz=eeg_np[:,8]-eeg_np[:,9] \n    Cz_Pz =eeg_np[:,9]-eeg_np[:,10]\n        \n    ch1=DWT_Features(Fp1_f7)\n    ch2=DWT_Features(F7_T3)\n    ch3=DWT_Features(T3_T5)\n    ch4=DWT_Features(T5_O1)\n        \n    ch5=DWT_Features(Fp2_F8)\n    ch6=DWT_Features(F8_T4)\n    ch7=DWT_Features(T4_T6)\n    ch8=DWT_Features(T6_O2)\n        \n    ch9=DWT_Features(Fp1_F3)\n    ch10=DWT_Features(F3_C3)\n    ch11=DWT_Features(C3_P3)\n    ch12=DWT_Features(P3_O1)\n        \n    ch13=DWT_Features(Fp2_F4)\n    ch14=DWT_Features(F4_C4)\n    ch15=DWT_Features(C4_P4)\n    ch16=DWT_Features(P4_O2)\n        \n    ch17=DWT_Features(Fz_Cz)\n    ch18=DWT_Features(Cz_Pz)\n    channel_features=np.concatenate((ch1,ch2,ch3,ch4,ch5,ch6,ch7,ch8,ch9,ch10,ch11,ch12,ch13,ch14,ch15,ch16,ch17,ch18), axis=None)\n   \n    return channel_features\n      ","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.491611Z","iopub.execute_input":"2025-04-07T23:20:15.491962Z","iopub.status.idle":"2025-04-07T23:20:15.508371Z","shell.execute_reply.started":"2025-04-07T23:20:15.491903Z","shell.execute_reply":"2025-04-07T23:20:15.507283Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Code to extract the feature vector from each EEG subsample with a 10 seconds observation window.","metadata":{}},{"cell_type":"code","source":"def shift_offset(Multichannel, file_tag):\n    \n    signal_features=np.array([])\n    aux=df_search_np[df_search_np==file_tag]\n    row=file_tag\n    offset_values=df_search_np[:len(aux)][aux==row, :][:,1]\n    Matrix_features=np.empty([len(aux),360])\n\n    Multichannel=Filtering(Multichannel)\n    \n    for i in range(len(aux)):\n        if len(np.argwhere(np.isnan(Multichannel).any(axis=1)))<=10:\n            eeg_np=Multichannel[(offset_values[i]+20)*200:(offset_values[i]+30)*200]\n            channel_features=single_channel(eeg_np)\n            Matrix_features[i:]=channel_features\n\n        else:\n            Matrix_features[:]=np.nan\n\n    return Matrix_features   \n\n","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.509359Z","iopub.execute_input":"2025-04-07T23:20:15.509641Z","iopub.status.idle":"2025-04-07T23:20:15.528513Z","shell.execute_reply.started":"2025-04-07T23:20:15.509617Z","shell.execute_reply":"2025-04-07T23:20:15.527092Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The following code is used to process each parquet file of the training EEGs.","metadata":{}},{"cell_type":"code","source":"def process_parquet_file(file_tag):\n    file_path = f'{EEG_Path}{file_tag}.parquet'\n    try:\n        if os.path.isfile(file_path):\n            eeg_np = pd.read_parquet(file_path)\n            result = shift_offset(eeg_np, file_tag)\n        else:\n            raise FileNotFoundError(f\"File {file_path} not found.\")\n    except Exception as e:\n        print(f\"Error processing file {file_tag}: {e}\")\n        rows = df_search_np[df_search_np[:, 0] == file_tag].shape[0]\n        result = np.full((rows, 360), np.nan)\n    return result","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.529696Z","iopub.execute_input":"2025-04-07T23:20:15.530145Z","iopub.status.idle":"2025-04-07T23:20:15.549252Z","shell.execute_reply.started":"2025-04-07T23:20:15.530106Z","shell.execute_reply":"2025-04-07T23:20:15.548120Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def process_iteration(EEG_Tag):\n    L=[]\n\n    with concurrent.futures.ProcessPoolExecutor() as executor:\n        Files_created=EEG_Tag\n        results=executor.map(process_parquet_file,Files_created)\n        for result in results:\n            L.append(result)\n           \n        Matrix_Features=np.vstack(L)\n    return Matrix_Features    ","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.550341Z","iopub.execute_input":"2025-04-07T23:20:15.550775Z","iopub.status.idle":"2025-04-07T23:20:15.567839Z","shell.execute_reply.started":"2025-04-07T23:20:15.550739Z","shell.execute_reply":"2025-04-07T23:20:15.566622Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Full=process_iteration(EEG_Tag)\nFull.shape","metadata":{"execution":{"iopub.status.busy":"2025-04-07T23:20:15.569095Z","iopub.execute_input":"2025-04-07T23:20:15.569489Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Once the feature matrix is obtained, the nan values must be removed and the labels converted into a format suitable for the classification algorithm.","metadata":{}},{"cell_type":"code","source":"def clean(df_train,Full):\n    y_labels=np.array(df_train[['seizure_vote','lpd_vote', 'gpd_vote','lrda_vote','grda_vote' , 'other_vote']])\n    k=1/y_labels.sum(axis=1)\n    t=k.reshape(-1,1)\n    y_prob=t*y_labels\n\n    index_nan=np.argwhere(np.asarray(np.isnan(Full),dtype=int).sum(axis=1)>0)\n    Full_clean=np.delete(Full, index_nan, axis=0)\n    y_labels_clean=np.delete(y_labels,index_nan,axis=0)\n    y_prob_clean=np.delete(y_prob,index_nan,axis=0)\n    return (Full_clean, y_labels_clean, y_prob_clean)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Full_clean, y_labels_clean, y_prob_clean=clean(df_train,Full)\ny_labels_clean","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The number of votes of the experts in each class (pattern) was transformed into a probability value, to convert the labels to a format suitable for classification. If a class receives all the votes, the probability that an EEG signal sample belongs to that class is 1, being 0 if it received no votes.\n\n","metadata":{}},{"cell_type":"code","source":"def discrep(y_labels_clean):\n    max=y_labels_clean.max(axis=1).reshape(-1, 1)\n    y_full_one_hot = np.where(y_labels_clean == max, 1, 0)\n    index_discrep=np.argwhere(y_full_one_hot.sum(axis=1)>1)\n    return y_full_one_hot,index_discrep","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"y_full_one_hot,index_discrep=discrep(y_labels_clean)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def consensus(Full_clean,y_full_one_hot,y_prob_clean,index_discrep):\n    X_consensus=np.delete(Full_clean, index_discrep, axis=0)\n    y_consensus=np.delete(y_full_one_hot,index_discrep,axis=0)\n    y_prob_cons=np.delete(y_prob_clean,index_discrep,axis=0)\n#\n    y_consensus=(y_consensus > 0).nonzero()[1]\n    sample_weights=y_prob_cons.max(axis=1)\n    return (X_consensus,y_consensus,y_prob_cons,sample_weights)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_consensus,y_consensus,y_prob_cons,sample_weights=consensus(Full_clean,y_full_one_hot,y_prob_clean,index_discrep)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def permutation (X_consensus,y_consensus,y_prob_cons,sample_weights):\n    permutation = np.array([i for i in range(len(X_consensus))])\n    np.random.seed(1)\n    np.random.shuffle(permutation)\n    X_perm = np.array([X_consensus[i] for i in permutation])\n    y_perm = np.array([y_consensus[i] for i in permutation])\n    sample_weights=np.array([sample_weights[i] for i in permutation])\n    y_prob_perm = np.array([y_prob_cons[i] for i in permutation])\n    factor=int((X_perm.shape[0])*0.8)\n    weights_train=sample_weights[:factor].flatten()\n    weights_val=sample_weights[factor:].flatten()\n    X_train=X_perm[:factor]\n    y_train=y_perm[:factor]\n    y_prob_train=y_prob_perm[:factor]\n    X_val=X_perm[factor:]\n    y_val=y_perm[factor:]\n    y_prob_val=y_prob_perm[factor:]\n\n\n###\n    X_discrep=Full_clean[index_discrep].reshape(-1,360)\n    y_prob_discrep=y_prob_clean[index_discrep].reshape(-1,6)\n    return (weights_train,weights_val,X_train,y_train,y_prob_train,X_val,y_val,y_prob_val,X_discrep,y_prob_discrep)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"weights_train,weights_val,X_train,y_train,y_prob_train,X_val,y_val,y_prob_val,X_discrep,y_prob_discrep=permutation(X_consensus,y_consensus,y_prob_cons,sample_weights)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"markdown","source":"Ensemble learning method was used with two base models: the XGBClassifier model from the XGBoost library and a neural network from the Keras library due to its versatility. With this method, it was observed that the predicted probability distribution is closer to the real probability distribution (expert votes), compared to the individual predictions of the models.\n\nThe first model was trained using all the samples with consensus among experts, considering their probability (sample_weight: model parameter). Although there is consensus, the probability that an EEG sample belongs to a label can be less than 1.  The discrepancy in votes between experts was used to evaluate its performance.\nThe neural network was designed with 5 hidden layers, and the Kullback Leibler Divergence (KLD) was set as its cost function. Consequently, in its training, the predicted probability distribution became closer to the observed probability distribution. 20% of the data with consensus was used as validation data.\n\nEach model was trained and evaluated with different data portions to reduce the bias that could be incurred in the predictions. Both models contained numerous parameters configured through trial and error until the results of their evaluation metrics were similar, ensuring no overfitting.\n\nOnce all predictors are trained, the ensemble can make a prediction for a new instance aggregating the predictions of the two models. The aggregation function is the average of the estimated probability distribution of each model.","metadata":{}},{"cell_type":"markdown","source":"**XGBoost Model**","metadata":{}},{"cell_type":"code","source":"X=np.concatenate([X_train,X_val],axis=0)\ny=np.concatenate([y_train,y_val],axis=0)\nw=np.concatenate([weights_train,weights_val],axis=0)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"xgb_model = xgb.XGBClassifier(objective='multi:softprob',n_classes=6, max_depth=8,subsample=0.8, colsample_bytree=0.8,tree_method='hist',learning_rate=0.04, device='cuda:1',n_estimators=60,verbosity=0)\nxgb_model.fit(X, y,sample_weight=[(w)])","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Kullback-Leibler Divergence**","metadata":{}},{"cell_type":"code","source":"kl = tf.keras.losses.KLDivergence()\nprint(kl(y_prob_discrep,xgb_model.predict_proba(X_discrep)).numpy())\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"****Keras model****","metadata":{}},{"cell_type":"code","source":"#X_traink=np.concatenate([X_train,X_discrep],axis=0)\nX_traink=X_train\n#y_prob_traink=np.concatenate([y_prob_train,y_prob_discrep],axis=0)\ny_prob_traink=y_prob_train\nX_valk=X_val  \ny_prob_valk=y_prob_val","metadata":{"_kg_hide-output":true,"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"inputs = keras.Input(shape=(360,), name=\"digits\")\nx = layers.Dense(300, activation=\"relu\", name=\"dense_1\")(inputs)\nx = layers.Dense(200, activation=\"relu\", name=\"dense_2\")(x)\nx = layers.Dense(100, activation=\"relu\", name=\"dense_3\")(x)\nx = layers.Dense(100, activation=\"relu\", name=\"dense_4\")(x)\nx = layers.Dense(100, activation=\"relu\", name=\"dense_5\")(x)\noutputs = layers.Dense(6, activation=\"softmax\", name=\"predictions\")(x)\nmodel_x= keras.Model(inputs=inputs, outputs=outputs)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"opt=keras.optimizers.Adam(\n    learning_rate=0.0005,\n    beta_1=0.9,\n    beta_2=0.999,\n    epsilon=1e-07,\n    amsgrad=False,\n    weight_decay=None,\n    clipnorm=None,\n    clipvalue=None,\n    global_clipnorm=None,\n    use_ema=False,\n    ema_momentum=0.99,\n    ema_overwrite_frequency=None,\n    name=\"adam\"\n)\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model_x.compile( optimizer=opt,\n    loss=keras.losses.KLDivergence(reduction=\"sum_over_batch_size\", name=\"kl_divergence\"),\n                )","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"np.random.seed(8)\nkeras.utils.set_random_seed(8)\nprint(\"Fit model on training data\")\n\nfrom tensorflow.keras.callbacks import EarlyStopping\n\nearly_stop = EarlyStopping(monitor='val_loss', patience=3, restore_best_weights=True)\n\nhistory = model_x.fit(\n    X_traink,\n    y_prob_traink,\n    batch_size=32,\n    epochs=50,\n    validation_data=(X_valk, y_prob_valk),\n    callbacks=[early_stop]\n)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pd.DataFrame(history.history).plot(figsize=(8, 5))\nplt.grid(True)\nplt.gca().set_ylim(0, 4)\nplt.gca().set_title('Loss Function (Keras Model)')\nplt.gca().set_xlabel('Epoch')\nplt.gca().set_ylabel('Kullback Leibler Divergence')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_discrep.shape","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"p1=xgb_model.predict_proba(X_discrep)\np2=model_x.predict(X_discrep)\np12=p1*0.5+p2*0.5\n\nkl1=kl(y_prob_discrep,p1).numpy()\nkl2=kl(y_prob_discrep,p2).numpy()\nkl12=kl(y_prob_discrep,p12).numpy()\nd = {'Model': ['xgb_model', 'model_x', 'Ensemble'], 'KLD': [kl1, kl2,kl12]}\ndf = pd.DataFrame(data=d)\ndf","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"A significant decrease in the KLD value can be observed just by taking the average value of the prediction of the two models (ensemble model).","metadata":{}},{"cell_type":"markdown","source":"# **Discussion**","metadata":{}},{"cell_type":"markdown","source":"Several studies were found whose objective was to classify the correct label to which a certain epilepsy signal belonged, but did not consider its probability,  so these studies cannot be used as a reference for comparison. Therefore, the result metric obtained will be compared with the best metric in this competition (March 28, 2024).\n\nThe ideal value of KLD is 0, which implies that the predicted probability distribution and the observed probability distribution are equal. This study's minimum value was between **0.67** and **0.78**. In the competition, the minimum KLD value is **0.22**.\n\nOne of the reasons for the low performance is that the method used with DWT for feature extraction may have omitted important information for classification. On the other hand, the competition-winning models use transfer learning techniques, employing pre-built convolutional neural network architectures such as EfficientNet and ResNet. These architectures are more complex and have many more parameters, but they are state-of-the-art for EEG signal classification. Another reason is that only the EEG data were used for the design of the model, the competition also included spectrogram data that could have been processed to create a third model and to include its predictions in the aggregation function with the two previous models. It would have improved the estimation of the probability distribution.\n\n","metadata":{}},{"cell_type":"markdown","source":"**Test data file preprocessing**","metadata":{}},{"cell_type":"code","source":"df_test=pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/test.csv')\ndf_test","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_search=df_train[[\"eeg_id\"]]\ntest_search_np=df_search.to_numpy(dtype=int)\n\nEEG_TEST=df_test['eeg_id'].unique()\ntest_Path='/kaggle/input/hms-harmful-brain-activity-classification/test_eegs/'","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def no_offset(Multichannel, file_tag):\n     \n    Multichannel=Filtering(Multichannel)\n    if Multichannel.shape[0]>5000:\n        eeg_np=Multichannel[20*200:30*200]\n        if len(np.argwhere(np.isnan(Multichannel).any(axis=1)))<=10:\n                \n            channel_features=single_channel(eeg_np)\n            Matrix_features=channel_features\n\n        else:\n            Matrix_features=np.empty(1,360)\n            Matrix_features[:]=np.nan\n   \n    return Matrix_features   ","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def process_test_file(file_tag):\n\n    if os.path.isfile(f'{test_Path}{file_tag}.parquet'):\n        eeg_np = pd.read_parquet(f'{test_Path}{file_tag}.parquet')\n        result = no_offset(eeg_np,file_tag) \n        \n    return result","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"L=[]\n\nwith concurrent.futures.ProcessPoolExecutor() as executor:\n    Files_created=EEG_TEST\n    results=executor.map(process_test_file,Files_created)\n    for result in results:\n        L.append(result)\n        print(\"\", len(L), end=\"\")\n    Test_Features=np.vstack(L)\nTest_Features.shape","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Inference**","metadata":{}},{"cell_type":"code","source":"pred1=xgb_model.predict_proba(Test_Features)\npred2=model_x.predict(Test_Features)\n\nprediccion=pred1*0.5+pred2*0.5","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Submission**","metadata":{}},{"cell_type":"code","source":"egg_t = df_test[[\"eeg_id\"]].copy()\ntarget =[ 'seizure_vote','lpd_vote', 'gpd_vote','lrda_vote','grda_vote' , 'other_vote']\negg_t[target] = prediccion.tolist()\n\nsub= pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/sample_submission.csv')\nsub= sub[[\"eeg_id\"]].copy()\nsub = sub.merge(egg_t, on=\"eeg_id\", how=\"left\")\nsub.to_csv(\"submission.csv\", index=False)\nsub.head()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# References","metadata":{}},{"cell_type":"markdown","source":"1. I. B. Slimen, L. Boubchir , Z. Mbarki , H. Seddik1, “EEG epileptic seizure detection and classification based on dual-tree complex wavelet transform and machine learning algorithms”, The Journal of Biomedical Research, 34(3): 151-161, 2020.\n\n2. O. Faust, U. Rajendra Acharya b, H. Adeli, A. Adeli, “Wavelet-based EEG processing for computer-aided seizure detection and epilepsy diagnosis”, vol. 26, pp. 56-64, 2015.\n\n3. G. R. Lee, R. Gommers, F. Wasilewski, K.Wohlfahrt, A O’Leary, “PyWavelets: A Python package for wavelet analysis”, Journal of Open Source Software, 2019.\n\n4. H. Albaqami, G. Mubashar Hassan, A. Dat, “Wavelet-based multiclass seizure type classification system”, Applied Sciences, 2022.","metadata":{}}]}