{"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":7392775,"sourceType":"datasetVersion","datasetId":4297782}],"dockerImageVersionId":30673,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport tensorflow as tf\nfrom matplotlib import pyplot as plt\nimport keras\nfrom keras import ops\nfrom sklearn.model_selection import KFold, GroupKFold\nfrom sklearn.preprocessing import StandardScaler\nfrom keras.models import Sequential, Model\nfrom keras.layers import Input, Dense, Activation, Dropout\nfrom keras.layers import BatchNormalization\nfrom keras.utils import to_categorical\nfrom keras import regularizers, layers\nfrom keras.callbacks import History, EarlyStopping\nfrom keras.optimizers import SGD\nimport sys\n#from tf.keras.optimizers import Adam\nimport os, gc\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# sample = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/sample_submission.csv\")\n# test_dataset = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/test.csv\")\n# seizure_vote = [1/6] * test_dataset.size\n# lpd_vote = [1/6] * test_dataset.size\n# gpd_vote = [1/6] * test_dataset.size\n# lrda_vote = [1/6] * test_dataset.size\n# grda_vote = [1/6] * test_dataset.size\n# other_vote = [1/6] * test_dataset.size\n\n\n# prediction_list = list(zip(test_dataset[\"eeg_id\"],seizure_vote,lpd_vote,gpd_vote,lrda_vote,grda_vote,other_vote))\n# submission_df = pd.DataFrame(prediction_list, columns = [sample.columns])\n\n# submission_df.to_csv('submission.csv', index=False)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-14T00:10:09.930986Z","iopub.execute_input":"2024-04-14T00:10:09.931415Z","iopub.status.idle":"2024-04-14T00:10:23.700396Z","shell.execute_reply.started":"2024-04-14T00:10:09.931388Z","shell.execute_reply":"2024-04-14T00:10:23.699590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load and Investigate Dataset","metadata":{}},{"cell_type":"code","source":"train_dataset = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/train.csv\")\ntargets = train_dataset.columns[-6:]\n\nprint(\"Training data shape:\", train_dataset.shape)\nprint(\"Training targets:\", targets)\nprint(train_dataset)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:10:23.701919Z","iopub.execute_input":"2024-04-14T00:10:23.702456Z","iopub.status.idle":"2024-04-14T00:10:23.965760Z","shell.execute_reply.started":"2024-04-14T00:10:23.702429Z","shell.execute_reply":"2024-04-14T00:10:23.964726Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The 'train.csv' file contains 106,800 entries with primary keys in eeg_id, spectrogram_id, label_id, and patient_id. Of note is that an EEG and its corresponding spectrogram are not unique to a sample, and multiple test inputs can come from the same eeg_id/spectrogram_id. In these cases, the eeg_label_offset_seconds column is used to isolate a segment of the entire EEG recording for processing. An isolated EEG subsample consists of the time series data in the time range [offset_seconds, offset_seconds+50] in seconds.","metadata":{}},{"cell_type":"code","source":"print(\"Number of unique EEG ids:\", len(pd.unique(train_dataset['eeg_id'])))\nprint(\"Number of unique Spectrogram ids:\", len(pd.unique(train_dataset['spectrogram_id'])))\nprint(\"Number of unique patients:\", len(pd.unique(train_dataset['patient_id'])))","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:10:23.967020Z","iopub.execute_input":"2024-04-14T00:10:23.967493Z","iopub.status.idle":"2024-04-14T00:10:23.979905Z","shell.execute_reply.started":"2024-04-14T00:10:23.967433Z","shell.execute_reply":"2024-04-14T00:10:23.978909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Factoring for entries coming from the same EEG, the training dataset consists of 17,089 unique EEGs and 11,138 unique spectrograms, from 1,950 different patients. The EEG subsample consists of the 50-second window beginning at the specified offset, while its accompanying spectrogram consists of the 600-second window beginning at the specified offset.","metadata":{}},{"cell_type":"code","source":"sample_eeg_parquet = pd.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/1003529080.parquet', engine='pyarrow')\nprint(sample_eeg_parquet)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:10:23.982284Z","iopub.execute_input":"2024-04-14T00:10:23.983073Z","iopub.status.idle":"2024-04-14T00:10:24.180676Z","shell.execute_reply.started":"2024-04-14T00:10:23.983043Z","shell.execute_reply":"2024-04-14T00:10:24.179680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_spectrogram_parquet = pd.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/1000086677.parquet', engine='pyarrow')\nprint(sample_spectrogram_parquet)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:10:24.181948Z","iopub.execute_input":"2024-04-14T00:10:24.182308Z","iopub.status.idle":"2024-04-14T00:10:24.242223Z","shell.execute_reply.started":"2024-04-14T00:10:24.182274Z","shell.execute_reply":"2024-04-14T00:10:24.241273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each eeg_id contains a corresponding parquet file with a matching name to its numeric id value. A full EEG parquet consists of one or more (potentially overlapping) 50-second time windows of data, sampled at a rate of 200 samples/second. Time series data for 19 electrode locations and heart data are provided in an EEG. \n\nA full spectrogram parquet likewise consists of one or more (potentially overlapping) 600-second time windows of data, sampled at a rate of 0.5 samples/second (2 second intervals). Time series data for 400 columns are included, among the 4 regions [LL, RL, LP, RP] and with 100 frequencies from each region.","metadata":{}},{"cell_type":"markdown","source":"## Loading the dataset\n\nAs there are 11,138 unique spectrograms, loading each parquet file one at a time takes significant time. A dataset has been aggregated and publicly provided by Chris Doette to only require reading a single file significantly quicker. The following code allows the use of this dataset for both EEG data and spectrogram data.","metadata":{}},{"cell_type":"code","source":"%%time\nREAD_SPEC_FILES = False #Required to be false to load the single npy file\nREAD_EEG_SPEC_FILES = False\n\n# READ ALL SPECTROGRAMS\nPATH = '/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/'\nfiles = os.listdir(PATH)\nprint(f'There are {len(files)} spectrogram parquets')\nspectrograms = np.load('/kaggle/input/brain-spectrograms/specs.npy',allow_pickle=True).item()","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:10:24.243632Z","iopub.execute_input":"2024-04-14T00:10:24.244005Z","iopub.status.idle":"2024-04-14T00:11:24.529681Z","shell.execute_reply.started":"2024-04-14T00:10:24.243971Z","shell.execute_reply":"2024-04-14T00:11:24.528629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train_dataset.groupby('eeg_id')[['spectrogram_id','spectrogram_label_offset_seconds']].agg(\n    {'spectrogram_id':'first','spectrogram_label_offset_seconds':'min'}) #Find minimum offset for each unique id\ntrain.columns = ['spec_id','min']\n\ntmp = train_dataset.groupby('eeg_id')[['spectrogram_id','spectrogram_label_offset_seconds']].agg(\n    {'spectrogram_label_offset_seconds':'max'}) #Find maximum offset for each unique id\ntrain['max'] = tmp\n\ntmp = train_dataset.groupby('eeg_id')[['patient_id']].agg('first') #Add patient id to dataframe\ntrain['patient_id'] = tmp\n\ntmp = train_dataset.groupby('eeg_id')[targets].agg('sum') #Add vote counts per id to dataframe\nfor target in targets:\n    train[target] = tmp[target].values\n\n#Normalize vote counts\ny_data = train[targets].values\ny_data = y_data / y_data.sum(axis=1,keepdims=True)\ntrain[targets] = y_data\n\ntmp = train_dataset.groupby('eeg_id')[['expert_consensus']].agg('first') #Add consensus vote target to dataframe\ntrain['target'] = tmp\n\ntrain = train.reset_index()\nprint('Train non-overlapp eeg_id shape:', train.shape )\n\nprint(train)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:11:24.530709Z","iopub.execute_input":"2024-04-14T00:11:24.530959Z","iopub.status.idle":"2024-04-14T00:11:24.614094Z","shell.execute_reply.started":"2024-04-14T00:11:24.530937Z","shell.execute_reply":"2024-04-14T00:11:24.613119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Once the specs.npy file is read, and a list of unique spectrograms generated, ","metadata":{}},{"cell_type":"code","source":"%time\n# ENGINEER FEATURES\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# FEATURE NAMES\nSPEC_COLS = pd.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/1000086677.parquet').columns[1:]\nFEATURES = [f'{c}_mean_10m' for c in SPEC_COLS]\nFEATURES += [f'{c}_min_10m' for c in SPEC_COLS]\nFEATURES += [f'{c}_mean_20s' for c in SPEC_COLS]\nFEATURES += [f'{c}_min_20s' for c in SPEC_COLS]\nprint(f'We are creating {len(FEATURES)} features for {len(train)} rows... ',end='')\n\n# print(SPEC_COLS)\ndata = np.zeros((len(train),len(FEATURES)))\nfor k in range(len(train)):\n    if k%100==0: print(k,', ',end='')\n    row = train.iloc[k]\n    r = int( (row['min'] + row['max'])//4 ) \n\n    # 10 MINUTE WINDOW FEATURES (MEANS and MINS)\n    x = np.nanmean(spectrograms[row.spec_id][r:r+300,:],axis=0)\n    data[k,:400] = x\n    x = np.nanmin(spectrograms[row.spec_id][r:r+300,:],axis=0)\n    data[k,400:800] = x\n\n    # 20 SECOND WINDOW FEATURES (MEANS and MINS)\n    x = np.nanmean(spectrograms[row.spec_id][r+145:r+155,:],axis=0)\n    data[k,800:1200] = x\n    x = np.nanmin(spectrograms[row.spec_id][r+145:r+155,:],axis=0)\n    data[k,1200:1600] = x\n\n\ntrain[FEATURES] = data\nprint(); print('New train shape:',train.shape)\n\n# FREE MEMORY\ndel spectrograms, data\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:11:24.615099Z","iopub.execute_input":"2024-04-14T00:11:24.615373Z","iopub.status.idle":"2024-04-14T00:11:41.784593Z","shell.execute_reply.started":"2024-04-14T00:11:24.615338Z","shell.execute_reply":"2024-04-14T00:11:41.783331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:11:41.785748Z","iopub.execute_input":"2024-04-14T00:11:41.786002Z","iopub.status.idle":"2024-04-14T00:11:41.813960Z","shell.execute_reply.started":"2024-04-14T00:11:41.785980Z","shell.execute_reply":"2024-04-14T00:11:41.813076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A feature set of 4 values is created along each of the 400 region-frequency pairs for a total of 1600 features. Those values are:\n\n1. mean of the 10 min window\n2. min of the 10 min window\n3. mean of the 20 second window\n4. min of the 20 second window\n\nThese 4 quantities are appended to the training data and can be used to train a classifier.","metadata":{}},{"cell_type":"markdown","source":"# Buliding a model\n\nEach sample consists of a spectrogram containing a majority vote target among the 6 classes as a target, and 1600 features consisting of the 4 mean and min over time values of the 400 combinations of region+frequency.","metadata":{}},{"cell_type":"markdown","source":"# K-fold validation\nWe use a 5-fold cross validation for this model. ","metadata":{}},{"cell_type":"code","source":"def create_model():\n    model = Sequential()\n    model.add(Input(shape=(1600,)))\n    model.add(Dense(2100, activation=\"sigmoid\", name=\"HiddenLayer1\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1600, activation=\"sigmoid\", name=\"HiddenLayer2\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1600, activation=\"sigmoid\", name=\"HiddenLayer3\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1600, activation=\"sigmoid\", name=\"HiddenLayer4\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1024, activation=\"sigmoid\", name=\"HiddenLayer5\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1024, activation=\"sigmoid\", name=\"HiddenLayer6\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(1024, activation=\"sigmoid\", name=\"HiddenLayer7\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(512, activation=\"sigmoid\", name=\"HiddenLayer8\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))    \n\n    model.add(Dense(256, activation=\"sigmoid\", name=\"HiddenLayer9\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))    \n\n    model.add(Dense(32, activation=\"sigmoid\", name=\"HiddenLayer10\"))\n    model.add(BatchNormalization())\n    model.add(Dropout(0.2))\n\n    model.add(Dense(6, activation=\"softmax\", name=\"Output\"))\n\n    opt = SGD(learning_rate=1e-3, momentum=0.9, clipvalue=1)\n    model.compile(loss='sparse_categorical_crossentropy', optimizer=opt, metrics=['sparse_categorical_accuracy'])\n    return model","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:11:41.816663Z","iopub.execute_input":"2024-04-14T00:11:41.816977Z","iopub.status.idle":"2024-04-14T00:11:41.829881Z","shell.execute_reply.started":"2024-04-14T00:11:41.816952Z","shell.execute_reply":"2024-04-14T00:11:41.829010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_oof = []\nall_true = []\ntarget_map = {'Seizure':0, 'LPD':1, 'GPD':2, 'LRDA':3, 'GRDA':4, 'Other':5}\nhistory = []\n\nmodel = create_model()\nmodel.summary()\nmodel.save(\"TestModel.keras\")\n\ntrain[FEATURES] = train[FEATURES].fillna(0)\ngroup_kfold = GroupKFold(n_splits=5)\nfor i, (train_index, valid_index) in enumerate(group_kfold.split(train[:], train[:].target, train[:].patient_id)):\n    print(f'### Fold {i+1}')\n    print(f'### train size {len(train_index)}, valid size {len(valid_index)}')\n    \n    #training data\n    data_train = train.loc[train_index,FEATURES] #extract feature data\n    label_train = train.loc[train_index,'target'].map(target_map) #extract label data\n    #label_train = to_categorical(label_train, 6)\n    \n    #validation data\n    data_valid = train.loc[valid_index,FEATURES]\n    label_valid = train.loc[valid_index,'target'].map(target_map)\n    #label_valid = to_categorical(label_valid, 6)\n    \n    #Normalize input features\n    scaler = StandardScaler().fit(data_train)\n    data_train = scaler.transform(data_train)\n    data_valid = scaler.transform(data_valid)  \n \n    #print(np.any(np.isnan(data_train)))\n\n    assert not np.any(np.isnan(data_train))\n    \n    \n    # Checkpoint Callback\n    checkpoint_filepath = 'epoch-{epoch:02d}-loss-{val_loss:.3f}.weights.h5'\n    checkpoint = keras.callbacks.ModelCheckpoint(\n        filepath=checkpoint_filepath,\n        save_weights_only=True,\n        monitor='val_loss',\n        mode='min',\n        save_best_only=True,\n        save_freq='epoch')\n    es = EarlyStopping(monitor='val_loss',mode='min', patience=15)\n    callbacks_list = [checkpoint, es]\n    fold_history = model.fit(x=data_train, y=label_train, verbose=1, epochs=50, validation_data=(data_valid,label_valid),callbacks=callbacks_list)\n    history.append(fold_history)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:11:41.831010Z","iopub.execute_input":"2024-04-14T00:11:41.831272Z","iopub.status.idle":"2024-04-14T00:18:32.339344Z","shell.execute_reply.started":"2024-04-14T00:11:41.831250Z","shell.execute_reply":"2024-04-14T00:18:32.338541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(history[4].history['sparse_categorical_accuracy'])\nplt.plot(history[4].history['val_sparse_categorical_accuracy'])\nplt.title('model accuracy')\nplt.ylabel('accuracy')\nplt.xlabel('epoch')\nplt.legend(['train', 'val'], loc='upper left')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:18:32.341160Z","iopub.execute_input":"2024-04-14T00:18:32.341447Z","iopub.status.idle":"2024-04-14T00:18:32.674626Z","shell.execute_reply.started":"2024-04-14T00:18:32.341421Z","shell.execute_reply":"2024-04-14T00:18:32.673679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"validation_predictions = model.predict(data_valid)\nprint(validation_predictions)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T00:21:06.669372Z","iopub.execute_input":"2024-04-14T00:21:06.669744Z","iopub.status.idle":"2024-04-14T00:21:06.989280Z","shell.execute_reply.started":"2024-04-14T00:21:06.669717Z","shell.execute_reply":"2024-04-14T00:21:06.988315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference on Test Dataset\n## Load the test dataset and apply preprocessing","metadata":{}},{"cell_type":"code","source":"test_dataset = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/test.csv\")\nSPEC_PATH = '/kaggle/input/hms-harmful-brain-activity-classification/test_spectrograms/'\ndata = np.zeros((len(test_dataset),len(FEATURES)))\n\nfor k in range(len(test_dataset)):\n    row = test_dataset.iloc[k]\n    s = int( row.spectrogram_id )\n    spec = pd.read_parquet(f'{SPEC_PATH}{s}.parquet')\n    \n    # 10 MINUTE WINDOW FEATURES\n    x = np.nanmean( spec.iloc[:,1:].values, axis=0)\n    data[k,:400] = x\n    x = np.nanmin( spec.iloc[:,1:].values, axis=0)\n    data[k,400:800] = x\n\n    # 20 SECOND WINDOW FEATURES\n    x = np.nanmean( spec.iloc[145:155,1:].values, axis=0)\n    data[k,800:1200] = x\n    x = np.nanmin( spec.iloc[145:155,1:].values, axis=0)\n    data[k,1200:1600] = x\n\ndata = scaler.transform(data)\ntest_dataset[FEATURES] = data\nprint('New test shape',test_dataset.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T01:33:25.340336Z","iopub.execute_input":"2024-04-14T01:33:25.341293Z","iopub.status.idle":"2024-04-14T01:33:26.270798Z","shell.execute_reply.started":"2024-04-14T01:33:25.341254Z","shell.execute_reply":"2024-04-14T01:33:26.269881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict on the test dataset and save to csv","metadata":{}},{"cell_type":"code","source":"predictions = model.predict(data)\n\neeg_df = pd.DataFrame(test_dataset[\"eeg_id\"], columns = [\"eeg_id\"])\nsubmission_df = pd.DataFrame(predictions, columns = [\"seizure_vote\",\"lpd_vote\",\"gpd_vote\",\"lrda_vote\",\"grda_vote\",\"other_vote\"])\n\nsubmission_df = pd.concat([eeg_df, submission_df], axis = 1)\nprint(submission_df)\nsubmission_df.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2024-04-14T02:17:17.910329Z","iopub.execute_input":"2024-04-14T02:17:17.911275Z","iopub.status.idle":"2024-04-14T02:17:17.980532Z","shell.execute_reply.started":"2024-04-14T02:17:17.911237Z","shell.execute_reply":"2024-04-14T02:17:17.979636Z"},"trusted":true},"execution_count":null,"outputs":[]}]}