{"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":"gpu","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":7465251,"sourceType":"datasetVersion","datasetId":4317718}],"dockerImageVersionId":30648,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd, numpy as np, os\nimport matplotlib.pyplot as plt\nimport gc\nfrom tqdm import tqdm\nfrom tensorflow.keras import layers, models, optimizers, losses, activations\nfrom keras.callbacks import ModelCheckpoint\nfrom sklearn.linear_model import LogisticRegression, SGDClassifier\nfrom sklearn.naive_bayes import GaussianNB, MultinomialNB\nfrom sklearn.ensemble import StackingClassifier, VotingClassifier\nfrom tensorflow.keras.utils import plot_model, to_categorical\nfrom sklearn.metrics import RocCurveDisplay\nfrom sklearn.model_selection import StratifiedKFold\nfrom scipy.signal import butter, lfilter\nimport tensorflow as tf\nfrom keras.saving import load_model\n\nfrom sklearn.model_selection import KFold, GroupKFold\nimport tensorflow.keras.backend as K, gc","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-05T12:45:03.365136Z","iopub.execute_input":"2024-02-05T12:45:03.365465Z","iopub.status.idle":"2024-02-05T12:45:17.186472Z","shell.execute_reply.started":"2024-02-05T12:45:03.365427Z","shell.execute_reply":"2024-02-05T12:45:17.185663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using [@cdeotte's](https://www.kaggle.com/cdeotte) data generation","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 = True","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:45:17.187943Z","iopub.execute_input":"2024-02-05T12:45:17.188454Z","iopub.status.idle":"2024-02-05T12:45:17.454093Z","shell.execute_reply.started":"2024-02-05T12:45:17.188428Z","shell.execute_reply":"2024-02-05T12:45:17.453171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-02-05T12:45:17.455173Z","iopub.execute_input":"2024-02-05T12:45:17.455457Z","iopub.status.idle":"2024-02-05T12:45:17.637544Z","shell.execute_reply.started":"2024-02-05T12:45:17.455432Z","shell.execute_reply":"2024-02-05T12:45:17.636663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('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-02-05T12:45:17.640277Z","iopub.execute_input":"2024-02-05T12:45:17.640739Z","iopub.status.idle":"2024-02-05T12:45:17.645920Z","shell.execute_reply.started":"2024-02-05T12:45:17.640705Z","shell.execute_reply":"2024-02-05T12:45:17.645106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\n\n","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:45:17.647452Z","iopub.execute_input":"2024-02-05T12:45:17.647754Z","iopub.status.idle":"2024-02-05T12:45:17.657126Z","shell.execute_reply.started":"2024-02-05T12:45:17.647719Z","shell.execute_reply":"2024-02-05T12:45:17.656194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nall_eegs = np.load('/kaggle/input/brain-eegs/eegs.npy',allow_pickle=True).item()","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:45:17.658292Z","iopub.execute_input":"2024-02-05T12:45:17.658648Z","iopub.status.idle":"2024-02-05T12:46:50.934652Z","shell.execute_reply.started":"2024-02-05T12:45:17.658616Z","shell.execute_reply":"2024-02-05T12:46:50.933656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# LOAD TRAIN\nEEG_IDS = train.eeg_id.unique()\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-02-05T12:46:50.935990Z","iopub.execute_input":"2024-02-05T12:46:50.936299Z","iopub.status.idle":"2024-02-05T12:46:51.228423Z","shell.execute_reply.started":"2024-02-05T12:46:50.936273Z","shell.execute_reply":"2024-02-05T12:46:51.227421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.signal import butter, lfilter\ndef butter_lowpass_filter(data, cutoff_freq=20, sampling_rate=200, order=7): \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-02-05T12:46:51.229938Z","iopub.execute_input":"2024-02-05T12:46:51.230365Z","iopub.status.idle":"2024-02-05T12:46:51.238731Z","shell.execute_reply.started":"2024-02-05T12:46:51.230327Z","shell.execute_reply":"2024-02-05T12:46:51.237651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class 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\n        return X[:,::self.downsample,:], y\n        #return decimate(X, self.downsample, axis=1), 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\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\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        return X,y","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:46:51.240210Z","iopub.execute_input":"2024-02-05T12:46:51.240530Z","iopub.status.idle":"2024-02-05T12:46:51.259748Z","shell.execute_reply.started":"2024-02-05T12:46:51.240475Z","shell.execute_reply":"2024-02-05T12:46:51.258662Z"},"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-02-05T12:46:51.262450Z","iopub.execute_input":"2024-02-05T12:46:51.262779Z","iopub.status.idle":"2024-02-05T12:46:51.274510Z","shell.execute_reply.started":"2024-02-05T12:46:51.262754Z","shell.execute_reply":"2024-02-05T12:46:51.273597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Implementation of ChronoNet architecture","metadata":{}},{"cell_type":"code","source":"def chronoNet(input_shape=(2000,8), num_classes=6):\n    print(f'Making ChronoNet with input {input_shape} and {num_classes} classes')\n    inp = layers.Input(shape=input_shape)\n    \n    # Inception Style Convolutional layers\n    cat = inp\n    for i in range(3):\n        conv1 = layers.Conv1D(128, 2, strides=2, padding='same')(cat)\n        conv2 = layers.Conv1D(128, 4, strides=2, padding='same')(cat)\n        conv3 = layers.Conv1D(128, 8, strides=2, padding='same')(cat)\n        cat = layers.concatenate([conv1, conv2, conv3])\n    \n    # DenseNet style RNN layers\n    rnn1 = layers.GRU(128, return_sequences=True)(cat)\n    rnn2 = layers.GRU(128, return_sequences=True)(rnn1)\n    cat = layers.concatenate([rnn1, rnn2])\n    for i in range(3):\n        rnn = layers.GRU(128, return_sequences=True)(cat)\n        cat = layers.concatenate([cat, rnn])\n\n    rnn = layers.GRU(128)(cat)\n\n    x = layers.Dense(6, activation='softmax', dtype='float32')(rnn)\n\n    # COMPILE MODEL\n    model = tf.keras.Model(inputs=inp, outputs=x)\n    loss = tf.keras.losses.KLDivergence()\n    opt = tf.keras.optimizers.AdamW(learning_rate = 1e-3, clipvalue=50.0)\n    model.compile(loss=loss, \n                  optimizer = opt, \n                  metrics=['accuracy', tf.keras.metrics.AUC()],)\n    return model\n    \nmodel = chronoNet()\nmodel.summary()\nplot_model(model)","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:46:51.275965Z","iopub.execute_input":"2024-02-05T12:46:51.276241Z","iopub.status.idle":"2024-02-05T12:46:54.220345Z","shell.execute_reply.started":"2024-02-05T12:46:51.276216Z","shell.execute_reply":"2024-02-05T12:46:54.219460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\nLR_START = 1e-6\nLR_MAX = 1e-3\nLR_MIN = 1e-6\nLR_RAMPUP_EPOCHS = 0\nLR_SUSTAIN_EPOCHS = 0\nEPOCHS2 = 30\n\ndef lrfn(epoch):\n    if epoch < LR_RAMPUP_EPOCHS:\n        lr = (LR_MAX - LR_START) / LR_RAMPUP_EPOCHS * epoch + LR_START\n    elif epoch < LR_RAMPUP_EPOCHS + LR_SUSTAIN_EPOCHS:\n        lr = LR_MAX\n    else:\n        decay_total_epochs = EPOCHS2 - LR_RAMPUP_EPOCHS - LR_SUSTAIN_EPOCHS - 1\n        decay_epoch_index = epoch - LR_RAMPUP_EPOCHS - LR_SUSTAIN_EPOCHS\n        phase = math.pi * decay_epoch_index / decay_total_epochs\n        cosine_decay = 0.5 * (1 + math.cos(phase))\n        lr = (LR_MAX - LR_MIN) * cosine_decay + LR_MIN\n    return lr\n\nrng = [i for i in range(EPOCHS2)]\nlr_y = [lrfn(x) for x in rng]\nplt.figure(figsize=(10, 4))\nplt.plot(rng, lr_y, '-o')\nplt.xlabel('epoch',size=14); plt.ylabel('learning rate',size=14)\nplt.title('Cosine Training Schedule',size=16); plt.show()\n\nLR2 = tf.keras.callbacks.LearningRateScheduler(lrfn, verbose = True)","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:46:54.221736Z","iopub.execute_input":"2024-02-05T12:46:54.222089Z","iopub.status.idle":"2024-02-05T12:46:54.445049Z","shell.execute_reply.started":"2024-02-05T12:46:54.222057Z","shell.execute_reply":"2024-02-05T12:46:54.444131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gkf = 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    mix_gen = DataGenerator(train.iloc[train_index], shuffle=True, batch_size=32)\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    model = chronoNet()\n    checkpoint_callback = ModelCheckpoint(\n            f'best_1D-Chrono_f{i}.h5',\n            monitor='val_loss',  # Choose the metric to monitor for saving the best model\n            save_best_only=True,  # Save only the best model\n            mode='min',  # 'min' if monitoring loss, 'max' if monitoring accuracy\n            verbose=1  # Display messages about the saving process\n    )\n\n    model.fit(train_gen, verbose=1,\n              validation_data = valid_gen,\n              epochs=30,\n              callbacks=[checkpoint_callback],\n                  )","metadata":{"execution":{"iopub.status.busy":"2024-02-05T12:46:54.446120Z","iopub.execute_input":"2024-02-05T12:46:54.446373Z"},"trusted":true},"execution_count":null,"outputs":[]}]}