{"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"}],"dockerImageVersionId":30635,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Analysis Goal\nIn this competition, our aim is to classify brain activity using EEG and spectrogram recordings. The train.csv file includes segment IDs for both EEG and spectrogram data, along with expert consensus on brain states during specific periods. This notebook focuses on transforming EEG signals into images by extracting features like PSD, band power, spectral centroid, and a custom spectrogram (distinct from the dataset's provided spectrogram). We plan to train a CNN model on these EEG-derived images for classification. Additionally, we will train a separate CNN solely on the dataset's spectrograms to create a composite voting classifier for the final submission.\n\n# Load Data\nHere, we provide a brief overview of the data structure.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport os\nimport librosa\nimport librosa.display\nimport matplotlib.pyplot as plt\nimport warnings\nwarnings.filterwarnings('ignore')\nimport numpy as np\nimport pickle\nimport joblib\nfrom sklearn.model_selection import train_test_split\nfrom tensorflow.keras import models, layers\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential, Model\nfrom tensorflow.keras.layers import Dense, Conv2D, Flatten, MaxPooling2D, Dropout, GlobalAveragePooling2D\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder, normalize, MinMaxScaler\nfrom tensorflow.keras.applications import InceptionResNetV2\nfrom sklearn.metrics import confusion_matrix, accuracy_score\nfrom tensorflow.keras.utils import Sequence\nimport scipy.signal\nfrom scipy.signal import coherence, cwt, ricker\n\n# Load the dataset\nroot = '/kaggle/input/hms-harmful-brain-activity-classification/'\ndata = pd.read_csv(root + 'train.csv')\ndata","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-01-13T11:45:26.932261Z","iopub.execute_input":"2024-01-13T11:45:26.933715Z","iopub.status.idle":"2024-01-13T11:45:27.182992Z","shell.execute_reply.started":"2024-01-13T11:45:26.933682Z","shell.execute_reply":"2024-01-13T11:45:27.181003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_train_eeg = pd.read_parquet(\"/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/1000913311.parquet\")\nsample_train_eeg","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:27.188308Z","iopub.execute_input":"2024-01-13T11:45:27.188575Z","iopub.status.idle":"2024-01-13T11:45:27.283828Z","shell.execute_reply.started":"2024-01-13T11:45:27.188550Z","shell.execute_reply":"2024-01-13T11:45:27.278759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_train_spec = pd.read_parquet(\"/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/1000086677.parquet\")\nsample_train_spec","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:27.288078Z","iopub.execute_input":"2024-01-13T11:45:27.288347Z","iopub.status.idle":"2024-01-13T11:45:27.398326Z","shell.execute_reply.started":"2024-01-13T11:45:27.288323Z","shell.execute_reply":"2024-01-13T11:45:27.393855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Engineering\nTo create manageable-sized images from EEG time series for CNN training, we have developed several lightweight feature extraction functions. You can find more details in this [tutorial](https://towardsdatascience.com/cnns-for-audio-classification-6244954665ab).","metadata":{}},{"cell_type":"code","source":"fs = 200  # Sampling frequency\n\ndef calculate_psd(eeg_signal, fs=200, nfft=260):\n    f, Pxx_den = scipy.signal.periodogram(eeg_signal, fs=fs, nfft=nfft)\n    return Pxx_den[:130]\n\ndef calculate_band_power(eeg_signal, band, fs=200, window_size=200, step_size=77):\n    band_powers = []\n    for start in range(0, len(eeg_signal) - window_size + 1, step_size):\n        segment = eeg_signal[start:start + window_size]\n        f, Pxx = scipy.signal.welch(segment, fs, nperseg=window_size)\n        band_power = np.trapz(Pxx[np.logical_and(f >= band[0], f <= band[1])], f[np.logical_and(f >= band[0], f <= band[1])])\n        band_powers.append(band_power)\n    return np.array(band_powers)\n\ndef calculate_spectral_centroid(eeg_signal, fs=200, window_size=200, step_size=77):\n    spectral_centroids = []\n    for start in range(0, len(eeg_signal) - window_size + 1, step_size):\n        segment = eeg_signal[start:start + window_size]\n        f, Pxx_den = scipy.signal.periodogram(segment, fs)\n        spectral_centroid = (f * Pxx_den).sum() / Pxx_den.sum() if Pxx_den.sum() != 0 else 0\n        spectral_centroids.append(spectral_centroid)\n    return np.array(spectral_centroids)\n\ndef calculate_eeg_spectrogram(eeg_signal, fs=200):\n    nperseg = 256 # Window size for each FFT\n    noverlap = nperseg // 2  # 50% overlap\n    nfft = nperseg # Number of FFT points, usually the same as nperseg for simplicity\n    f, t, Sxx = scipy.signal.spectrogram(eeg_signal, fs=fs, nperseg=nperseg, noverlap=noverlap, nfft=nfft)\n    return Sxx\n\n# Function to pad the feature matrix to a fixed size\ndef padding(array, xx, yy):\n    h = array.shape[0]\n    w = array.shape[1]\n    a = max((xx - h) // 2, 0)\n    aa = max(0, xx - a - h)\n    b = max(0, (yy - w) // 2)\n    bb = max(yy - b - w, 0)\n    return np.pad(array, pad_width=((a, aa), (b, bb)), mode='constant')","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:27.401890Z","iopub.execute_input":"2024-01-13T11:45:27.402205Z","iopub.status.idle":"2024-01-13T11:45:27.423235Z","shell.execute_reply.started":"2024-01-13T11:45:27.402174Z","shell.execute_reply":"2024-01-13T11:45:27.421085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Generator\nSince EEG data is not directly available in `train.csv`, a data generator is defined to fetch EEG data corresponding to `train.csv` entries. This generator also performs feature engineering, ultimately producing images for training. The generated data image will be a tensor with depth 20 as the EEG signals contains 19 channels in different location, as well as a EKG channel, which may contain useful cardiac information.","metadata":{}},{"cell_type":"code","source":"class EEGDataGenerator(Sequence):\n    def __init__(self, df, batch_size, root, max_size=130, scaler=None):\n        self.df = df\n        self.batch_size = batch_size\n        self.root = root\n        self.max_size = max_size\n        self.scaler = scaler if scaler else StandardScaler()\n\n    def generate_eeg_features(self, eeg_data):\n        channel_features_list = []\n        max_feature_length = self.max_size\n        scaler = MinMaxScaler()\n\n        for channel in eeg_data.columns:\n            channel_data = eeg_data[channel]\n            # Calculate features\n            psd = calculate_psd(channel_data)\n            band_power_ts = calculate_band_power(channel_data, (8, 12), fs=200, window_size=256, step_size=128)\n            spectral_centroid_ts = calculate_spectral_centroid(channel_data, fs=200, window_size=256, step_size=128)\n            spectrogram = calculate_eeg_spectrogram(channel_data)\n\n            # Normalize features\n            norm_psd = scaler.fit_transform(psd.reshape(-1, 1))\n            norm_band_power = scaler.fit_transform(band_power_ts.reshape(-1, 1))\n            norm_spectral_centroid = scaler.fit_transform(spectral_centroid_ts.reshape(-1, 1))\n            norm_spectrogram = scaler.fit_transform(spectrogram)\n            \n            # Pad features for this channel\n            padded_psd = padding(norm_psd, max_feature_length, 1)\n            padded_band_power = padding(norm_band_power, max_feature_length, 1)\n            padded_spectral_centroid = padding(norm_spectral_centroid, max_feature_length, 1)\n            padded_spectrogram = padding(norm_spectrogram, max_feature_length, norm_spectrogram.shape[1])\n            \n            # Stack features for this channel\n            channel_features = np.hstack([\n                padded_psd, padded_band_power, \n                padded_spectral_centroid, padded_spectrogram\n            ])\n\n            channel_features_list.append(channel_features)\n\n        # Stack all channel features depth-wise to create a 3D tensor\n        eeg_tensor = np.stack(channel_features_list, axis=-1)\n        return eeg_tensor\n\n    def __len__(self):\n        return int(np.ceil(len(self.df) / self.batch_size))\n\n    def __getitem__(self, idx):\n        batch_df = self.df.iloc[idx * self.batch_size:(idx + 1) * self.batch_size]\n        X, y = [], []\n\n        for _, row in batch_df.iterrows():\n            eeg_id = row['eeg_id']\n            eeg_offset_seconds = row['eeg_label_offset_seconds']\n            eeg_data = pd.read_parquet(self.root + 'train_eegs/' + str(eeg_id) + '.parquet')\n\n            # Segment the EEG data\n            start_ind = int(eeg_samples_per_second * eeg_offset_seconds)\n            end_ind = start_ind + eeg_samples_per_second * eeg_seconds_per_subsample\n            subsample_eeg_data = eeg_data.iloc[start_ind:end_ind]\n\n            # Check if the subsample is complete\n            if len(subsample_eeg_data) == eeg_samples_per_second * eeg_seconds_per_subsample:\n                feature_image = self.generate_eeg_features(subsample_eeg_data)\n                X.append(feature_image)\n\n                # Prepare labels\n                total_votes = row[vote_list].sum()\n                normalized_votes = [vote / total_votes for vote in row[vote_list]]\n                y.append(normalized_votes)\n\n        return np.array(X), np.array(y)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:29.564065Z","iopub.execute_input":"2024-01-13T11:45:29.564426Z","iopub.status.idle":"2024-01-13T11:45:29.593814Z","shell.execute_reply.started":"2024-01-13T11:45:29.564401Z","shell.execute_reply":"2024-01-13T11:45:29.591920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dimensions of images generated by the data generator are verified by a dummy dataframe with the same sample number. One channel has one image, and all images are stacked along the depth direction, forming a 3D tensor for CNN training.","metadata":{}},{"cell_type":"code","source":"# Create a dummy instance of EEGDataGenerator\ndummy_df = pd.DataFrame()  # Empty DataFrame, as we won't use it here\ndummy_batch_size = 1       # Dummy batch size, as we only need to generate a single image\neeg_generator = EEGDataGenerator(dummy_df, dummy_batch_size, '')  # Assuming root is not needed for generating a single image\n\n# Sample EEG data segment\nsample_eeg_data = pd.DataFrame(np.random.randn(10000, 20), columns=[f'Channel_{i}' for i in range(1, 21)])\n\n# Generate EEG image\neeg_image = eeg_generator.generate_eeg_features(sample_eeg_data)\n\n# Print the shape of the generated image\nprint(\"Shape of the generated EEG image:\", eeg_image.shape)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:29.725840Z","iopub.execute_input":"2024-01-13T11:45:29.726208Z","iopub.status.idle":"2024-01-13T11:45:31.869566Z","shell.execute_reply.started":"2024-01-13T11:45:29.726180Z","shell.execute_reply":"2024-01-13T11:45:31.863963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Training\nWe start by defining constants for the data generator, setting the stage for training.","metadata":{}},{"cell_type":"code","source":"# Constants\neeg_samples_per_second = 200\neeg_seconds_per_subsample = 50\n    \n# Load the labels\nvote_list = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']\n\n# Create data generator\nbatch_size = 64\ndata_generator = EEGDataGenerator(data, batch_size, root)\n\n# Split data into training and validation sets\ntrain_df, val_df = train_test_split(data, test_size=0.2, random_state=42)\ntrain_generator = EEGDataGenerator(train_df, batch_size, root)\nval_generator = EEGDataGenerator(val_df, batch_size, root)\n\nprint(train_df.shape, val_df.shape)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:31.875454Z","iopub.execute_input":"2024-01-13T11:45:31.875801Z","iopub.status.idle":"2024-01-13T11:45:31.934935Z","shell.execute_reply.started":"2024-01-13T11:45:31.875754Z","shell.execute_reply":"2024-01-13T11:45:31.933101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training a pre-trained InceptionResNetV2","metadata":{}},{"cell_type":"code","source":"# # Define the CNN model\n# def build_cnn_model(input_shape, num_classes):\n#     model = Sequential()\n#     model.add(layers.Conv2D(32, (3, 3), strides=(2, 2), activation='relu', input_shape=input_shape))\n#     model.add(layers.MaxPooling2D((2, 2)))\n#     model.add(layers.Dropout(0.2))\n#     model.add(layers.Conv2D(64, (3, 3), strides=(2, 2), activation='relu'))\n#     model.add(layers.Flatten())\n#     model.add(layers.Dense(64, activation='relu'))\n#     model.add(layers.Dropout(0.2))\n#     model.add(Dense(num_classes, activation='sigmoid'))  # Sigmoid for multi-label classification\n#     model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy']) \n#     return model\n\n# # Model training\n# input_shape = eeg_image.shape  # Shape of the input data\n# num_classes = len(vote_list)  # Number of classes\n# model = build_cnn_model(input_shape, num_classes)\n# history = model.fit(train_generator, epochs=10, validation_data=val_generator)\n\n# # Evaluate the model\n# predictions = model.predict(X_val)\n\n# # Plot accuracy\n# plt.plot(history.history['accuracy'])\n# plt.plot(history.history['val_accuracy'])\n# plt.title('Model Accuracy')\n# plt.ylabel('Accuracy')\n# plt.xlabel('Epoch')\n# plt.legend(['Train', 'Validation'], loc='upper left')\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:31.941103Z","iopub.execute_input":"2024-01-13T11:45:31.941563Z","iopub.status.idle":"2024-01-13T11:45:31.966166Z","shell.execute_reply.started":"2024-01-13T11:45:31.941535Z","shell.execute_reply":"2024-01-13T11:45:31.962837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the pre-trained InceptionResNetV2 model\nbase_model = InceptionResNetV2(include_top=False, weights='imagenet', input_shape=input_shape)\n\n# Freeze the layers of the base model\nfor layer in base_model.layers:\n    layer.trainable = False\n\n# Add custom layers\nx = base_model.output\nx = GlobalAveragePooling2D()(x)\nx = Dense(1024, activation='relu')(x)  # Example: one dense layer\npredictions = Dense(num_classes, activation='sigmoid')(x)  # Final layer with sigmoid activation\n\n# Create the final model\nmodel = Model(inputs=base_model.input, outputs=predictions)\n\n# Compile the model\nmodel.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n\n# Train the model\nhistory = model.fit(train_generator, epochs=10, validation_data=val_generator)\n\n# Plot training results\nplt.plot(history.history['accuracy'])\nplt.plot(history.history['val_accuracy'])\nplt.title('Model Accuracy')\nplt.ylabel('Accuracy')\nplt.xlabel('Epoch')\nplt.legend(['Train', 'Validation'], loc='upper left')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:31.970344Z","iopub.execute_input":"2024-01-13T11:45:31.971264Z","iopub.status.idle":"2024-01-13T11:45:32.246476Z","shell.execute_reply.started":"2024-01-13T11:45:31.971229Z","shell.execute_reply":"2024-01-13T11:45:32.237883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"test_eeg_data_list = []\nfor eeg_id in test_df['eeg_id']:\n    eeg_data = pd.read_parquet(os.path.join(test_eeg_dir, str(eeg_id) + '.parquet'))\n    eeg_offset_seconds = 0\n    start_ind = int(eeg_samples_per_second * eeg_offset_seconds)\n    subsample_eeg_data = eeg_data[start_ind:start_ind + eeg_samples_per_second * eeg_seconds_per_subsample]\n\n    # Generate features for this subsample\n    feature_image = eeg_generator.generate_eeg_features(subsample_eeg_data)\n\n    # Reshape feature image for model prediction\n    feature_image_reshaped = feature_image.reshape((1,) + feature_image.shape)  # Adding batch dimension\n\n    test_eeg_data_list.append(feature_image_reshaped)\n\n# Predict using the CNN model\nsolution_list = []\nfor i, test_eeg_data in enumerate(test_eeg_data_list):\n    predictions = model.predict(test_eeg_data)\n    predictions = predictions.flatten()  # Flatten to get a 1D array of predictions\n\n    test_results_dict = {'eeg_id': test_df['eeg_id'].iloc[i]}\n    for ind, vote in enumerate(vote_list):\n        test_results_dict[vote] = predictions[ind]\n    solution_list.append(test_results_dict)\n\nsolution_df = pd.DataFrame(solution_list)\nsolution_df.to_csv(\"submission.csv\", index=False)\nprint(solution_df)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T11:45:32.257606Z","iopub.status.idle":"2024-01-13T11:45:32.258063Z","shell.execute_reply.started":"2024-01-13T11:45:32.257873Z","shell.execute_reply":"2024-01-13T11:45:32.257893Z"},"trusted":true},"execution_count":null,"outputs":[]}]}