{"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":"gpu","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":7392733,"sourceType":"datasetVersion","datasetId":4297749},{"sourceId":7465251,"sourceType":"datasetVersion","datasetId":4317718}],"dockerImageVersionId":30636,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# WaveNet Starter using RAW EEG Features!\nThis notebook is a WaveNet starter for Kaggle's Brain comp. It achieves `CV 0.91` and `LB 0.66`. Note that submitting train means achieves `CV 1.26` and `LB 0.97` [[here][1]]. So this notebook's WaveNet is successfully learning to predict brain events from raw EEG waveforms!\n\nThis model only uses two features. We can engineer more features and/or modify the model architecture to improve CV score and LB score. Furthermore we can build 1 model which inputs both spectrogram images and eeg waveforms. The two EEG features in this notebook are:\n* feature 1 : `Fp1 minus O1`\n* feature 2 : `Fp2 minus O2`\n\nFeature 1 is the beginning of the montage chains `LL` and `LP` minus the ending of montage `LL` and `LP`. And feature 2 is the beginning of the montage chains `RL` and `RP` minus the ending of montage `RL` and `RP`.\n![](https://raw.githubusercontent.com/cdeotte/Kaggle_Images/main/Jan-2024/montage.png)\n\n# UPDATE\nIn version 7 and 8, we add more features and update the model architecture to evaluate each montage chain separately and then concatenate the features. This new architecture is motivated by the discovery of a better formula to utilize EEG explained in discussion [here][2]\n\n* Version 5,6: Use 2 features - CV 0.91 LB 0.66\n* Version 7,8: Use 8 features grouped as 4 chains. Downsample time 5x - **CV 0.81 LB 0.53**, wow!\n\nWe train our new model in version 7, then save model weights. Then load model into version 8 to submit to LB.\n\n![](https://raw.githubusercontent.com/cdeotte/Kaggle_Images/main/Jan-2024/wave-model.png)\n\n# Train and Infer Tricks\nWe train the fold models in version 7 of the notebook and submit to Kaggle LB in version 8 of the notebook. This makes submission faster because we train the fold models for 30 minutes in version 7 then save them. In version 8, we just load the models without needing to retrain models during Kaggle submit. (And we train our old model in 5 annd submit in 6).\n\nVersion 4 uses `1xP100` GPU with full precision and takes 1 hour to train 5 folds 5 epochs of WaveNet. Version 5 uses `2xT4` GPU with mixed precision and takes 30 minutes to train 5 folds 5 epochs of WaveNet. \n\n[1]: https://www.kaggle.com/code/seshurajup/eda-train-csv\n[2]: https://www.kaggle.com/competitions/hms-harmful-brain-activity-classification/discussion/469760","metadata":{}},{"cell_type":"markdown","source":"# Load Train Data","metadata":{}},{"cell_type":"code","source":"import pandas as pd, numpy as np, os\nimport matplotlib.pyplot as plt\n\ntrain = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')\nprint( train.shape )\ndisplay( train.head() )\n\n# CHOICE TO CREATE OR LOAD EEGS FROM NOTEBOOK VERSION 1\nCREATE_EEGS = False\nTRAIN_MODEL = False","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:50:00.030732Z","iopub.execute_input":"2024-03-30T12:50:00.031141Z","iopub.status.idle":"2024-03-30T12:50:00.314493Z","shell.execute_reply.started":"2024-03-30T12:50:00.031112Z","shell.execute_reply":"2024-03-30T12:50:00.313548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Raw EEG Features","metadata":{}},{"cell_type":"code","source":"df = pd.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/1000913311.parquet')\nFEATS = df.columns\nprint(f'There are {len(FEATS)} raw eeg features')\nprint( list(FEATS) )","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:50:02.757804Z","iopub.execute_input":"2024-03-30T12:50:02.758656Z","iopub.status.idle":"2024-03-30T12:50:02.779127Z","shell.execute_reply.started":"2024-03-30T12:50:02.758626Z","shell.execute_reply":"2024-03-30T12:50:02.778147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- range(len(FEATS)): range(len(FEATS)) generates a sequence of numbers from 0 to the length of FEATS minus 1. This effectively creates a sequence of indices corresponding to the elements in FEATS.\n- zip(FEATS, range(len(FEATS))): zip combines the elements of FEATS with the corresponding indices generated by range(len(FEATS)). This pairs each element of FEATS with its index.\n- {x:y for x,y in ...}: This is a dictionary comprehension. iterate over the pairs generated by zip and creates a dictionary where the elements of FEATS are the keys (x) and their corresponding indices are the values (y).","metadata":{}},{"cell_type":"code","source":"# Feats = dataframe columns\nprint('We will use the following subset of raw EEG features:')\nFEATS = ['Fp1','T3','C3','O1','Fp2','C4','T4','O2']\nFEAT2IDX = {x:y for x,y in zip(FEATS,range(len(FEATS)))}\nprint( list(FEATS) )","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:50:03.257896Z","iopub.execute_input":"2024-03-30T12:50:03.258272Z","iopub.status.idle":"2024-03-30T12:50:03.264210Z","shell.execute_reply.started":"2024-03-30T12:50:03.258245Z","shell.execute_reply":"2024-03-30T12:50:03.263263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Read EEG data from Parquet file, process it, and optionally displays it as a plot.handle NaN values and extract the middle 50 seconds of the EEG data.\n\n\n- Reading Data: EEG data from a Parquet file specified by parquet_path using pd.read_parquet(). The columns to be read are specified by FEATS.\n\n- Extracting Middle 50 Seconds: It then extracts the middle 10,000 rows from the EEG data. The calculation to determine the offset is based on ensuring that the middle 50 seconds of data are selected.\n\n- Optional Display: If display parameter is set to True, plot the EEG data using matplotlib. create a plot for each feature (column) in the EEG data, with the x-axis representing the time and the y-axis representing the EEG values.\n\n- Converting to Numpy Array: convert the EEG data to a numpy array format. Before conversion, handle NaN (Not a Number) values by replacing them with the mean of the non-NaN values within each column.\n\n- Plotting (if display=True): It plots each EEG feature separately on the same plot, with each feature potentially having a different vertical offset for clarity.\n\n- Returning Data: returns EEG data in the form of a numpy array.","metadata":{}},{"cell_type":"code","source":"def eeg_from_parquet(parquet_path, display=False):\n    \n    # EXTRACT MIDDLE 50 SECONDS\n    eeg = pd.read_parquet(parquet_path, columns=FEATS)\n    rows = len(eeg)\n    offset = (rows-10_000)//2\n    eeg = eeg.iloc[offset:offset+10_000]\n    \n    if display: \n        plt.figure(figsize=(10,5))\n        offset = 0\n    \n    # CONVERT TO NUMPY\n    data = np.zeros((10_000,len(FEATS)))\n    for j,col in enumerate(FEATS):\n        \n        # FILL NAN\n        x = eeg[col].values.astype('float32')\n        m = np.nanmean(x)\n        if np.isnan(x).mean()<1: x = np.nan_to_num(x,nan=m)\n        else: x[:] = 0\n            \n        data[:,j] = x\n        \n        if display: \n            if j!=0: offset += x.max()\n            plt.plot(range(10_000),x-offset,label=col)\n            offset -= x.min()\n            \n    if display:\n        plt.legend()\n        name = parquet_path.split('/')[-1]\n        name = name.split('.')[0]\n        plt.title(f'EEG {name}',size=16)\n        plt.show()\n        \n    return data","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:50:05.418247Z","iopub.execute_input":"2024-03-30T12:50:05.418647Z","iopub.status.idle":"2024-03-30T12:50:05.430313Z","shell.execute_reply.started":"2024-03-30T12:50:05.418618Z","shell.execute_reply":"2024-03-30T12:50:05.429258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nall_eegs = {}\nDISPLAY = 4\nEEG_IDS = train.eeg_id.unique()\nPATH = '/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/'\n\nfor i,eeg_id in enumerate(EEG_IDS):\n    if (i%100==0)&(i!=0): print(i,', ',end='') \n    \n    # SAVE EEG TO PYTHON DICTIONARY OF NUMPY ARRAYS\n    data = eeg_from_parquet(f'{PATH}{eeg_id}.parquet', display=i<DISPLAY)              \n    all_eegs[eeg_id] = data\n    \n    if i==DISPLAY:\n        if CREATE_EEGS:\n            print(f'Processing {train.eeg_id.nunique()} eeg parquets... ',end='')\n        else:\n            print(f'Reading {len(EEG_IDS)} eeg NumPys from disk.')\n            break\n            \nif CREATE_EEGS: \n    np.save('eegs',all_eegs)\nelse:\n    all_eegs = np.load('/kaggle/input/brain-eegs/eegs.npy',allow_pickle=True).item()","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:50:05.678039Z","iopub.execute_input":"2024-03-30T12:50:05.678406Z","iopub.status.idle":"2024-03-30T12:51:37.031245Z","shell.execute_reply.started":"2024-03-30T12:50:05.678378Z","shell.execute_reply":"2024-03-30T12:51:37.029518Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Deduplicate Train EEG Id\n\n\n- Defining Targets:\n     - identify the target columns in the DataFrame df by selecting the last six columns and assigns them to the variable TARGETS.\n     - create a dictionary TARS that maps each target class to a numerical label.\n     - create a reverse dictionary TARS2 that maps numerical labels back to their corresponding target class names.\n\n- Preparing Training Data:\n     - groups data by the column 'eeg_id' and selects the first occurrence of 'patient_id' for each group, storing this in the train DataFrame.\n     - It aggregates the sum of target columns for each 'eeg_id' group and stores the result in the DataFrame tmp.\n     - normalize the target values by dividing each row by the sum of its elements.\n     - The normalized target values are assigned back to the train DataFrame.\n     - select the first occurrence of 'expert_consensus' for each 'eeg_id' group and stores it in tmp.\n     - The column 'target' is created in the train DataFrame and populated with the values from tmp.\n     \n- Filtering Data:\n     - reset the index of the train DataFrame.\n     - filter the data to include only rows where the 'eeg_id' is in the list EEG_IDS. This list is assumed to be defined elsewhere in the code.\n     - print the shape of the resulting training data DataFrame train.\n     - display the first few rows of the train DataFrame.","metadata":{}},{"cell_type":"code","source":"# LOAD TRAIN \ndf = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')\nTARGETS = df.columns[-6:]\nTARS = {'Seizure':0, 'LPD':1, 'GPD':2, 'LRDA':3, 'GRDA':4, 'Other':5}\nTARS2 = {x:y for y,x in TARS.items()}\n\ntrain = df.groupby('eeg_id')[['patient_id']].agg('first')\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()\ntrain = train.loc[train.eeg_id.isin(EEG_IDS)]\nprint('Train Data with unique eeg_id shape:', train.shape )\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:33.994346Z","iopub.execute_input":"2024-03-30T12:53:33.995093Z","iopub.status.idle":"2024-03-30T12:53:34.358264Z","shell.execute_reply.started":"2024-03-30T12:53:33.995061Z","shell.execute_reply":"2024-03-30T12:53:34.357345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Butterworth low-pass filter Example:\n- Butterworth filter involves several steps, including designing the filter coefficients, implementing the filter function, and applying it to the signal. ","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport scipy.signal\nimport matplotlib.pyplot as plt\n\ndef butter_lowpass(cutoff_freq, fs, order=5):\n    nyquist_freq = 0.5 * fs\n    normal_cutoff = cutoff_freq / nyquist_freq\n    b, a = scipy.signal.butter(order, normal_cutoff, btype='low', analog=False)\n    return b, a\n\ndef butter_lowpass_filter(data, cutoff_freq, fs, order=5):\n    b, a = butter_lowpass(cutoff_freq, fs, order=order)\n    y = scipy.signal.filtfilt(b, a, data)\n    return y\n\n# sample data\nfs = 1000.0  # Sample rate (Hz)\nt = np.arange(0, 10, 1/fs)  # 10 seconds of data\nf1 = 5.0  # Frequency of signal (Hz)\nf2 = 50.0  # Frequency of noise (Hz)\nsignal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)\nnoise = 0.5 * np.random.randn(len(t))  # Gaussian noise\ndata = signal + noise\n\n# filter parameters\ncutoff_freq = 30.0  # Cutoff frequency (Hz)\norder = 6  # Filter order\n\n# Applying Butterworth filter\nfiltered_data = butter_lowpass_filter(data, cutoff_freq, fs, order)\n\nplt.figure(figsize=(10, 6))\nplt.plot(t, data, 'b-', label='Original Signal')\nplt.plot(t, filtered_data, 'r-', linewidth=2, label='Filtered Signal')\nplt.xlabel('Time [s]')\nplt.ylabel('Amplitude')\nplt.title('Butterworth Low-pass Filter')\nplt.legend()\nplt.grid(True)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:36.893157Z","iopub.execute_input":"2024-03-30T12:53:36.894084Z","iopub.status.idle":"2024-03-30T12:53:37.746081Z","shell.execute_reply.started":"2024-03-30T12:53:36.894050Z","shell.execute_reply":"2024-03-30T12:53:37.744972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Butter Low-Pass Filter\n\n- four parameters:\n  - data: The input data to be filtered.\n  - cutoff_freq: The cutoff frequency of the low-pass filter. Frequencies higher than this will be attenuated.\n  - sampling_rate: The rate at which the data is sampled.\n  - order: The order of the Butterworth filter, determining the filter's roll-off.\n  \n-  0.5 * sampling_rate - The Nyquist frequency is half of the sampling rate. It represents the maximum frequency that can be represented accurately in the sampled data. This is important for filtering as it helps to define the normalized cutoff frequency.\n\n- normal_cutoff = cutoff_freq / nyquist - The cutoff frequency is normalized with respect to the Nyquist frequency. This normalization ensures that the cutoff frequency is expressed relative to the maximum representable frequency in the sampled data.\n\n- b, a = butter(order, normal_cutoff, btype='low', analog=False)  - butter function from the scipy.signal module is used to design a Butterworth filter. It returns the filter coefficients b (numerator) and a (denominator) for the difference equation representing the filter. The filter is specified as a low-pass filter (btype='low') and analog=False specifies that it's a digital filter.\n\n- filtered_data = lfilter(b, a, data, axis=0) - lfilter function from the scipy.signal module is used to apply the filter to the input data. It takes the filter coefficients b and a, along with the input data, and performs the filtering operation. The axis=0 argument specifies that the filtering should be done along the first axis of the data array.","metadata":{}},{"cell_type":"code","source":"from scipy.signal import butter, lfilter\n\ndef butter_lowpass_filter(data, cutoff_freq=20, sampling_rate=200, order=4):\n    nyquist = 0.5 * sampling_rate\n    normal_cutoff = cutoff_freq / nyquist\n    b, a = butter(order, normal_cutoff, btype='low', analog=False)\n    filtered_data = lfilter(b, a, data, axis=0)\n    return filtered_data","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:39.530671Z","iopub.execute_input":"2024-03-30T12:53:39.531498Z","iopub.status.idle":"2024-03-30T12:53:39.538757Z","shell.execute_reply.started":"2024-03-30T12:53:39.531458Z","shell.execute_reply":"2024-03-30T12:53:39.537672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"FREQS = [1,2,4,8,16][::-1]\nx = [all_eegs[EEG_IDS[0]][:,0]]\nfor k in FREQS:\n    x.append( butter_lowpass_filter(x[0], cutoff_freq=k) )\n\nplt.figure(figsize=(20,20))\nplt.plot(range(10_000),x[0], label='without filter')\nfor k in range(1,len(x)):\n    plt.plot(range(10_000),x[k]-k*(x[0].max()-x[0].min()), label=f'with filter {FREQS[k-1]}Hz')\nplt.legend()\nplt.title('Butter Low-Pass Filter Examples',size=18)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:41.788401Z","iopub.execute_input":"2024-03-30T12:53:41.788742Z","iopub.status.idle":"2024-03-30T12:53:42.677628Z","shell.execute_reply.started":"2024-03-30T12:53:41.788717Z","shell.execute_reply":"2024-03-30T12:53:42.676511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Loader with Butter Low-Pass Filter\n\n__data_generation Method:\n - Generate data for the given batch of samples.\n - Initialize arrays X and y to store input features and labels.\n - Loop over the indexes of the batch, retrieve corresponding data samples, and perform feature engineering on EEG data.\n - Standardize the feature values and apply Butterworth low-pass filter to the data.\n - Storeproce ssed samples in the X array and their corresponding labels in the y array.","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\n\nclass DataGenerator(tf.keras.utils.Sequence):\n    'Generates data for Keras'\n    def __init__(self, data, batch_size=32, shuffle=False, eegs=all_eegs, mode='train',\n                 downsample=5): \n\n        self.data = data\n        self.batch_size = batch_size\n        self.shuffle = shuffle\n        self.eegs = eegs\n        self.mode = mode\n        self.downsample = downsample\n        self.on_epoch_end()\n        \n    def __len__(self):\n        'Denotes the number of batches per epoch'\n        ct = int( np.ceil( len(self.data) / self.batch_size ) )\n        return ct\n\n    def __getitem__(self, index):\n        'Generate one batch of data'\n        indexes = self.indexes[index*self.batch_size:(index+1)*self.batch_size]\n        X, y = self.__data_generation(indexes)\n        return X[:,::self.downsample,:], y\n\n    def on_epoch_end(self):\n        'Updates indexes after each epoch'\n        self.indexes = np.arange( len(self.data) )\n        if self.shuffle: np.random.shuffle(self.indexes)\n                        \n    def __data_generation(self, indexes):\n        'Generates data containing batch_size samples' \n    \n        X = np.zeros((len(indexes),10_000,8),dtype='float32')\n        y = np.zeros((len(indexes),6),dtype='float32')\n        \n        sample = np.zeros((10_000,X.shape[-1]))\n        for j,i in enumerate(indexes):\n            row = self.data.iloc[i]      \n            data = self.eegs[row.eeg_id]\n            \n            # FEATURE ENGINEER\n            sample[:,0] = data[:,FEAT2IDX['Fp1']] - data[:,FEAT2IDX['T3']]\n            sample[:,1] = data[:,FEAT2IDX['T3']] - data[:,FEAT2IDX['O1']]\n            \n            sample[:,2] = data[:,FEAT2IDX['Fp1']] - data[:,FEAT2IDX['C3']]\n            sample[:,3] = data[:,FEAT2IDX['C3']] - data[:,FEAT2IDX['O1']]\n            \n            sample[:,4] = data[:,FEAT2IDX['Fp2']] - data[:,FEAT2IDX['C4']]\n            sample[:,5] = data[:,FEAT2IDX['C4']] - data[:,FEAT2IDX['O2']]\n            \n            sample[:,6] = data[:,FEAT2IDX['Fp2']] - data[:,FEAT2IDX['T4']]\n            sample[:,7] = data[:,FEAT2IDX['T4']] - data[:,FEAT2IDX['O2']]\n            \n            # STANDARDIZE\n            sample = np.clip(sample,-1024,1024)\n            sample = np.nan_to_num(sample, nan=0) / 32.0\n            \n            # BUTTER LOW-PASS FILTER\n            sample = butter_lowpass_filter(sample)\n            \n            X[j,] = sample\n            if self.mode!='test':\n                y[j] = row[TARGETS]\n            \n        return X,y","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:45.746026Z","iopub.execute_input":"2024-03-30T12:53:45.746655Z","iopub.status.idle":"2024-03-30T12:53:45.766912Z","shell.execute_reply.started":"2024-03-30T12:53:45.746626Z","shell.execute_reply":"2024-03-30T12:53:45.765965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Display Data Loader\n\n- for k in range(4):\n     - Loop through the first four samples in the batch.","metadata":{}},{"cell_type":"code","source":"gen = DataGenerator(train, shuffle=False)\n\nfor x,y in gen:\n    for k in range(4):\n        plt.figure(figsize=(20,4))\n        offset = 0\n        for j in range(x.shape[-1]):\n            if j!=0: offset -= x[k,:,j].min()\n            plt.plot(range(2_000),x[k,:,j]+offset,label=f'feature {j+1}')\n            offset += x[k,:,j].max()\n        tt = f'{y[k][0]:0.1f}'\n        for t in y[k][1:]:\n            tt += f', {t:0.1f}'\n        plt.title(f'EEG_Id = {EEG_IDS[k]}\\nTarget = {tt}',size=14)\n        plt.legend()\n        plt.show()\n    break","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:47.888351Z","iopub.execute_input":"2024-03-30T12:53:47.889300Z","iopub.status.idle":"2024-03-30T12:53:50.015948Z","shell.execute_reply.started":"2024-03-30T12:53:47.889259Z","shell.execute_reply":"2024-03-30T12:53:50.015075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Initialize GPUs","metadata":{}},{"cell_type":"code","source":"import os\nos.environ[\"CUDA_VISIBLE_DEVICES\"]=\"0,1\"\nimport tensorflow as tf\nprint('TensorFlow version =',tf.__version__)\n\n# USE MULTIPLE GPUS\ngpus = tf.config.list_physical_devices('GPU')\nif len(gpus)<=1: \n    strategy = tf.distribute.OneDeviceStrategy(device=\"/gpu:0\")\n    print(f'Using {len(gpus)} GPU')\nelse: \n    strategy = tf.distribute.MirroredStrategy()\n    print(f'Using {len(gpus)} GPUs')","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:52.398762Z","iopub.execute_input":"2024-03-30T12:53:52.399632Z","iopub.status.idle":"2024-03-30T12:53:52.410651Z","shell.execute_reply.started":"2024-03-30T12:53:52.399601Z","shell.execute_reply":"2024-03-30T12:53:52.409689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# USE MIXED PRECISION\nMIX = True\nif MIX:\n    tf.config.optimizer.set_experimental_options({\"auto_mixed_precision\": True})\n    print('Mixed precision enabled')\nelse:\n    print('Using full precision')","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:52.687854Z","iopub.execute_input":"2024-03-30T12:53:52.688196Z","iopub.status.idle":"2024-03-30T12:53:52.693798Z","shell.execute_reply.started":"2024-03-30T12:53:52.688172Z","shell.execute_reply":"2024-03-30T12:53:52.692639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Build WaveNet Model","metadata":{}},{"cell_type":"code","source":"# TRAIN SCHEDULE\ndef lrfn(epoch):\n        return [1e-3,1e-3,1e-4,1e-4,1e-5][epoch]\nLR = tf.keras.callbacks.LearningRateScheduler(lrfn, verbose = True)\nEPOCHS = 5","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:55.628393Z","iopub.execute_input":"2024-03-30T12:53:55.629717Z","iopub.status.idle":"2024-03-30T12:53:55.634362Z","shell.execute_reply.started":"2024-03-30T12:53:55.629681Z","shell.execute_reply":"2024-03-30T12:53:55.633483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n- dilation_rates = [2**i for i in range(n)]: create a list of dilation rates exponentially increasing from 1 to n with a base of 2. Dilation rate controls how far apart each unit in the convolutional layer is from its neighbors.\n\n- x = Conv1D(filters = filters, kernel_size = 1, padding = 'same')(x): apply a 1D convolutional layer with 1x1 kernel size to the input tensor x. This operation aims to transform the number of channels (filters) while preserving the input shape.\n\n- res_x = x: create a copy of the current state of x for later use.\n\n- Inside the loop:\n   - Two convolutional layers are applied to x. One with a 'tanh' activation function and the other with a 'sigmoid' activation function. These activations are then multiplied element-wise.\n   - The result is passed through another 1D convolutional layer.\n   - The original input x is added to the result of the convolution operation, forming a residual connection.\n","metadata":{}},{"cell_type":"code","source":"# Original code\n# from tensorflow.keras.layers import Input, Dense, Multiply, Add, Conv1D, Concatenate\n\n# def wave_block(x, filters, kernel_size, n):\n#     dilation_rates = [2**i for i in range(n)]\n#     x = Conv1D(filters = filters,\n#                kernel_size = 1,\n#                padding = 'same')(x)\n#     res_x = x\n#     for dilation_rate in dilation_rates:\n#         tanh_out = Conv1D(filters = filters,\n#                           kernel_size = kernel_size,\n#                           padding = 'same', \n#                           activation = 'tanh', \n#                           dilation_rate = dilation_rate)(x)\n#         sigm_out = Conv1D(filters = filters,\n#                           kernel_size = kernel_size,\n#                           padding = 'same',\n#                           activation = 'sigmoid', \n#                           dilation_rate = dilation_rate)(x)\n#         x = Multiply()([tanh_out, sigm_out])\n#         x = Conv1D(filters = filters,\n#                    kernel_size = 1,\n#                    padding = 'same')(x)\n#         res_x = Add()([res_x, x])\n#     return res_x","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:21:14.202946Z","iopub.execute_input":"2024-03-30T12:21:14.203375Z","iopub.status.idle":"2024-03-30T12:21:14.218091Z","shell.execute_reply.started":"2024-03-30T12:21:14.203344Z","shell.execute_reply":"2024-03-30T12:21:14.216853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from tensorflow.keras.layers import Conv1D, Multiply, Add\n# from tensorflow.keras import backend as K\n\n# def wave_block(x, filters, kernel_size, n):\n#     dilation_rates = [2**i for i in range(n)]\n#     x = Conv1D(filters=filters, kernel_size=1, padding='same')(x)\n#     res_x = x\n#     for dilation_rate in dilation_rates:\n#         tanh_out = Conv1D(filters=filters, kernel_size=kernel_size, padding='same',\n#                           activation='tanh', dilation_rate=dilation_rate)(x)\n#         sigm_out = Conv1D(filters=filters, kernel_size=kernel_size, padding='same',\n#                           activation='sigmoid', dilation_rate=dilation_rate)(x)\n#         merged = Multiply()([tanh_out, sigm_out])\n#         merged = Conv1D(filters=filters, kernel_size=1, padding='same')(merged)\n#         res_x = Add()([res_x, merged])\n#         x = merged\n#     return res_x","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:53:59.161616Z","iopub.execute_input":"2024-03-30T12:53:59.162328Z","iopub.status.idle":"2024-03-30T12:53:59.170383Z","shell.execute_reply.started":"2024-03-30T12:53:59.162297Z","shell.execute_reply":"2024-03-30T12:53:59.169229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"![](https://raw.githubusercontent.com/cdeotte/Kaggle_Images/main/Jan-2024/wave-model.png)","metadata":{}},{"cell_type":"code","source":"# def build_model():\n        \n#     # INPUT \n#     inp = tf.keras.Input(shape=(2_000,8))\n    \n#     ############\n#     # FEATURE EXTRACTION SUB MODEL\n#     inp2 = tf.keras.Input(shape=(2_000,1))\n#     x = wave_block(inp2, 8, 3, 12)\n#     x = wave_block(x, 16, 3, 8)\n#     x = wave_block(x, 32, 3, 4)\n#     x = wave_block(x, 64, 3, 1)\n#     model2 = tf.keras.Model(inputs=inp2, outputs=x)\n#     ###########\n    \n#     # LEFT TEMPORAL CHAIN\n#     x1 = model2(inp[:,:,0:1])\n#     x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)\n#     x2 = model2(inp[:,:,1:2])\n#     x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)\n#     z1 = tf.keras.layers.Average()([x1,x2])\n    \n#     # LEFT PARASAGITTAL CHAIN\n#     x1 = model2(inp[:,:,2:3])\n#     x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)\n#     x2 = model2(inp[:,:,3:4])\n#     x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)\n#     z2 = tf.keras.layers.Average()([x1,x2])\n    \n#     # RIGHT PARASAGITTAL CHAIN\n#     x1 = model2(inp[:,:,4:5])\n#     x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)\n#     x2 = model2(inp[:,:,5:6])\n#     x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)\n#     z3 = tf.keras.layers.Average()([x1,x2])\n    \n#     # RIGHT TEMPORAL CHAIN\n#     x1 = model2(inp[:,:,6:7])\n#     x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)\n#     x2 = model2(inp[:,:,7:8])\n#     x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)\n#     z4 = tf.keras.layers.Average()([x1,x2])\n    \n#     # COMBINE CHAINS\n#     y = tf.keras.layers.Concatenate()([z1,z2,z3,z4])\n#     y = tf.keras.layers.Dense(64, activation='relu')(y)\n#     y = tf.keras.layers.Dense(6,activation='softmax', dtype='float32')(y)\n    \n#     # COMPILE MODEL\n#     model = tf.keras.Model(inputs=inp, outputs=y)\n#     opt = tf.keras.optimizers.Adam(learning_rate = 1e-3)\n#     loss = tf.keras.losses.KLDivergence()\n#     model.compile(loss=loss, optimizer = opt)\n    \n#     return model","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:39:23.143222Z","iopub.execute_input":"2024-03-30T12:39:23.144290Z","iopub.status.idle":"2024-03-30T12:39:23.158339Z","shell.execute_reply.started":"2024-03-30T12:39:23.144257Z","shell.execute_reply":"2024-03-30T12:39:23.157354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow.keras.layers import Conv1D, Multiply, Add, GlobalAveragePooling1D, Average, Concatenate, Dense\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.losses import KLDivergence\nfrom tensorflow.keras import backend as K\n\ndef wave_block(x, filters, kernel_size, n, dropout_rate=0.2):\n    dilation_rates = [2**i for i in range(n)]\n    x = Conv1D(filters=filters, kernel_size=1, padding='same')(x)\n    res_x = x\n    for dilation_rate in dilation_rates:\n        tanh_out = Conv1D(filters=filters, kernel_size=kernel_size, padding='same',\n                          activation='tanh', dilation_rate=dilation_rate)(x)\n        sigm_out = Conv1D(filters=filters, kernel_size=kernel_size, padding='same',\n                          activation='sigmoid', dilation_rate=dilation_rate)(x)\n        merged = Multiply()([tanh_out, sigm_out])\n        merged = Conv1D(filters=filters, kernel_size=1, padding='same')(merged)\n        merged = tf.keras.layers.Dropout(dropout_rate)(merged)\n        res_x = Add()([res_x, merged])\n        x = merged\n    return res_x\n\ndef build_model(filters_list=(8, 16, 32, 64), kernel_size=3, dilation_rates=(12, 8, 4, 1), dropout_rate=0.2):\n    inp = tf.keras.Input(shape=(2000, 8))\n    \n    # Feature Extraction Sub Model\n    inp2 = tf.keras.Input(shape=(2000, 1))\n    x = inp2\n    for filters, dilation_rate in zip(filters_list, dilation_rates):\n        x = wave_block(x, filters, kernel_size, dilation_rate, dropout_rate=dropout_rate)\n    model2 = tf.keras.Model(inputs=inp2, outputs=x)\n    \n    # Chain Processing\n    chains = []\n    for i in range(4):\n        x1 = model2(inp[:,:,i*2:i*2+1])\n        x1 = GlobalAveragePooling1D()(x1)\n        x2 = model2(inp[:,:,i*2+1:i*2+2])\n        x2 = GlobalAveragePooling1D()(x2)\n        chains.append(Average()([x1, x2]))\n\n    # Combine Chains\n    y = Concatenate()(chains)\n    y = Dense(64, activation='relu')(y)\n    y = Dense(6, activation='softmax', dtype='float32')(y)\n    \n    # Compile Model\n    model = Model(inputs=inp, outputs=y)\n    opt = Adam(learning_rate=1e-3)\n    loss = KLDivergence()\n    model.compile(loss=loss, optimizer=opt)\n    \n    return model\n","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:55:19.020549Z","iopub.execute_input":"2024-03-30T12:55:19.021332Z","iopub.status.idle":"2024-03-30T12:55:19.036697Z","shell.execute_reply.started":"2024-03-30T12:55:19.021298Z","shell.execute_reply":"2024-03-30T12:55:19.035716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train Group KFold","metadata":{}},{"cell_type":"code","source":"VERBOSE = 1\nFOLDS_TO_TRAIN = 5\nif not os.path.exists('WaveNet_Model'):\n    os.makedirs('WaveNet_Model')\n\nfrom sklearn.model_selection import KFold, GroupKFold\nimport tensorflow.keras.backend as K, gc\n\nall_oof = []; all_oof2 = []; all_true = []\ngkf = GroupKFold(n_splits=5)\nfor i, (train_index, valid_index) in enumerate(gkf.split(train, train.target, train.patient_id)):   \n    \n    print('#'*25)\n    print(f'### Fold {i+1}')\n    train_gen = DataGenerator(train.iloc[train_index], shuffle=True, batch_size=32)\n    valid_gen = DataGenerator(train.iloc[valid_index], shuffle=False, batch_size=64, mode='valid')\n    print(f'### train size {len(train_index)}, valid size {len(valid_index)}')\n    print('#'*25)\n    \n    # TRAIN MODEL\n    K.clear_session()\n    with strategy.scope():\n        model = build_model()\n    if TRAIN_MODEL:\n        model.fit(train_gen, verbose=VERBOSE,\n              validation_data = valid_gen,\n              epochs=EPOCHS, callbacks = [LR])\n        model.save_weights(f'WaveNet_Model/WaveNet_fold{i}.h5')\n    else:\n        model.load_weights(f'/kaggle/input/brain-eegs/WaveNet_Model/WaveNet_fold{i}.h5')\n    \n    # WAVENET OOF\n    oof = model.predict(valid_gen, verbose=VERBOSE)\n    all_oof.append(oof)\n    all_true.append(train.iloc[valid_index][TARGETS].values)\n    \n    # TRAIN MEAN OOF\n    y_train = train.iloc[train_index][TARGETS].values\n    y_valid = train.iloc[valid_index][TARGETS].values\n    oof = y_valid.copy()\n    for j in range(6):\n        oof[:,j] = y_train[:,j].mean()\n    oof = oof / oof.sum(axis=1,keepdims=True)\n    all_oof2.append(oof)\n    \n    del model, oof, y_train, y_valid\n    gc.collect()\n    \n    if i==FOLDS_TO_TRAIN-1: break\n    \nall_oof = np.concatenate(all_oof)\nall_oof2 = np.concatenate(all_oof2)\nall_true = np.concatenate(all_true)","metadata":{"execution":{"iopub.status.busy":"2024-03-30T12:55:21.167879Z","iopub.execute_input":"2024-03-30T12:55:21.168253Z","iopub.status.idle":"2024-03-30T12:59:24.741977Z","shell.execute_reply.started":"2024-03-30T12:55:21.168224Z","shell.execute_reply":"2024-03-30T12:59:24.741143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CV Score for WaveNet","metadata":{}},{"cell_type":"code","source":"import sys\nsys.path.append('/kaggle/input/kaggle-kl-div')\nfrom kaggle_kl_div import score\n\noof = pd.DataFrame(all_oof.copy())\noof['id'] = np.arange(len(oof))\n\ntrue = pd.DataFrame(all_true.copy())\ntrue['id'] = np.arange(len(true))\n\ncv = score(solution=true, submission=oof, row_id_column_name='id')\nprint('CV Score with WaveNet Raw EEG =',cv)","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:02:25.428468Z","iopub.execute_input":"2024-03-30T13:02:25.428872Z","iopub.status.idle":"2024-03-30T13:02:25.493638Z","shell.execute_reply.started":"2024-03-30T13:02:25.428841Z","shell.execute_reply":"2024-03-30T13:02:25.492731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CV Score using Train Means","metadata":{}},{"cell_type":"code","source":"oof = pd.DataFrame(all_oof2.copy())\noof['id'] = np.arange(len(oof))\n\ntrue = pd.DataFrame(all_true.copy())\ntrue['id'] = np.arange(len(true))\n\ncv = score(solution=true, submission=oof, row_id_column_name='id')\nprint('CV Score with Train Means =',cv)","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:03:18.929775Z","iopub.execute_input":"2024-03-30T13:03:18.930160Z","iopub.status.idle":"2024-03-30T13:03:18.994835Z","shell.execute_reply.started":"2024-03-30T13:03:18.930132Z","shell.execute_reply":"2024-03-30T13:03:18.993868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submit to Kaggle LB","metadata":{}},{"cell_type":"code","source":"del all_eegs, train; gc.collect()\ntest = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/test.csv')\nprint('Test shape:',test.shape)\ntest.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:03:23.032692Z","iopub.execute_input":"2024-03-30T13:03:23.033093Z","iopub.status.idle":"2024-03-30T13:03:23.444408Z","shell.execute_reply.started":"2024-03-30T13:03:23.033062Z","shell.execute_reply":"2024-03-30T13:03:23.443356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_eegs2 = {}\nDISPLAY = 1\nEEG_IDS2 = test.eeg_id.unique()\nPATH2 = '/kaggle/input/hms-harmful-brain-activity-classification/test_eegs/'\n\nprint('Processing Test EEG parquets...'); print()\nfor i,eeg_id in enumerate(EEG_IDS2):\n        \n    # SAVE EEG TO PYTHON DICTIONARY OF NUMPY ARRAYS\n    data = eeg_from_parquet(f'{PATH2}{eeg_id}.parquet', i<DISPLAY)\n    all_eegs2[eeg_id] = data","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:03:26.164244Z","iopub.execute_input":"2024-03-30T13:03:26.164987Z","iopub.status.idle":"2024-03-30T13:03:27.024681Z","shell.execute_reply.started":"2024-03-30T13:03:26.164954Z","shell.execute_reply":"2024-03-30T13:03:27.023749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# INFER MLP ON TEST\npreds = []\nmodel = build_model()\ntest_gen = DataGenerator(test, shuffle=False, batch_size=64, eegs=all_eegs2, mode='test')\n\nprint('Inferring test... ',end='')\nfor i in range(FOLDS_TO_TRAIN):\n    print(f'fold {i+1}, ',end='')\n    if TRAIN_MODEL:\n        model.load_weights(f'WaveNet_Model/WaveNet_fold{i}.h5')\n    else:\n        model.load_weights(f'/kaggle/input/brain-eegs/WaveNet_Model/WaveNet_fold{i}.h5')\n    pred = model.predict(test_gen, verbose=0)\n    preds.append(pred)\npred = np.mean(preds,axis=0)\nprint()\nprint('Test preds shape',pred.shape)","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:03:29.703067Z","iopub.execute_input":"2024-03-30T13:03:29.703784Z","iopub.status.idle":"2024-03-30T13:04:05.754033Z","shell.execute_reply.started":"2024-03-30T13:03:29.703750Z","shell.execute_reply":"2024-03-30T13:04:05.752905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CREATE SUBMISSION.CSV\nfrom IPython.display import display\n\nsub = pd.DataFrame({'eeg_id':test.eeg_id.values})\nsub[TARGETS] = pred\nsub.to_csv('submission.csv',index=False)\nprint('Submission shape',sub.shape)\ndisplay( sub.head() )\n\n# SANITY CHECK TO CONFIRM PREDICTIONS SUM TO ONE\nprint('Sub row 0 sums to:',sub.iloc[0,-6:].sum())","metadata":{"execution":{"iopub.status.busy":"2024-03-30T13:04:18.578639Z","iopub.execute_input":"2024-03-30T13:04:18.579413Z","iopub.status.idle":"2024-03-30T13:04:18.598682Z","shell.execute_reply.started":"2024-03-30T13:04:18.579379Z","shell.execute_reply":"2024-03-30T13:04:18.597672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}