{"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":30823,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2024-12-19T12:49:28.177565Z","iopub.execute_input":"2024-12-19T12:49:28.177912Z","iopub.status.idle":"2024-12-19T12:49:28.18192Z","shell.execute_reply.started":"2024-12-19T12:49:28.177885Z","shell.execute_reply":"2024-12-19T12:49:28.180834Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-19T12:49:28.183036Z","iopub.execute_input":"2024-12-19T12:49:28.183278Z","iopub.status.idle":"2024-12-19T12:49:28.31977Z","shell.execute_reply.started":"2024-12-19T12:49:28.183256Z","shell.execute_reply":"2024-12-19T12:49:28.318648Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df1 = df.sample(10000, random_state=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-19T12:49:28.321389Z","iopub.execute_input":"2024-12-19T12:49:28.321711Z","iopub.status.idle":"2024-12-19T12:49:28.330731Z","shell.execute_reply.started":"2024-12-19T12:49:28.321675Z","shell.execute_reply":"2024-12-19T12:49:28.329789Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.fft import fft, ifft, fftshift, ifftshift\nfrom scipy.fft import dct, idct\n\ndef fdm(X, fs, fc, data_type='columns', filter_type='dct', sort_fc='descend', remove_mean=False, plot_subbands=True):\n# Take care of Vector inputs for einsum()\n# ValueError: einstein sum subscripts string contains too many subscripts for operand 0\n    if data_type == 'rows':\n        # axis = 1\n        X = np.transpose(X)\n    else:\n        axis = 0\n    if remove_mean:\n        X = X - np.mean(X, axis=0)\n    N = X.shape[0]\n    # print(N)\n    # fc=np.asarray([0, fs/64, fs/32, fs/16, fs/8, fs/4, fs/2])\n    fc = np.sort(fc)  # always sort fc in ascending order\n    # print(fc)\n    if fc[0] != 0:\n        fc = np.hstack((0, fc))\n    if fc[-1] != fs/2:\n        fc = np.hstack((fc, fs/2))\n    \n    if sort_fc == 'descend':\n        fc[::-1].sort()\n    # print(fc)\n    # Padding - Append data as symmetric rows (mirroring) at the beginning and end of the matrix\n    # if filter_type != 'dct':\n        # pad_rows = int(N * append_ratio);  # one-side length of appended rows\n        # # X_padded = np.pad(X, pad_width=((pad_rows,pad_rows), (0,0)), mode=\"symmetric\", reflect_type='even', )\n        # X_padded = np.pad(X, pad_width=((pad_rows,pad_rows), (0,0)), mode=\"reflect\", reflect_type='even')\n        # appended_length = X_padded.shape[0]  # Increased number of rows before and after appending data, twice of 'pad_rows'\n\n    # if filter_type == 'dct' or filter_type == 'dft':\n    if filter_type == 'dct':\n        dct_type = 2\n        K = np.round(2 * N * fc/fs).astype(int) # Vector operation - convert fc to k (digital freq)\n        no_of_subbands = K.shape[0] - 1         # For DCT implementation\n        Hk = np.zeros((N, 1, no_of_subbands))   # Create 3-D Matrix with third dimension as Fc filters\n\n    if filter_type == 'dft':\n        append_ratio = 0.02  # 2 %\n        pad_rows = int(N * append_ratio)  # one-side length of appended rows\n        X_padded = np.pad(X, pad_width=((pad_rows,pad_rows), (0,0)), mode=\"symmetric\", reflect_type='even')\n        # X_padded = np.pad(X, pad_width=((pad_rows,pad_rows), (0,0)), mode=\"reflect\", reflect_type='even')\n        appended_length = X_padded.shape[0]  # Increased number of rows before and after appending data, twice of 'pad_rows'\n        L = appended_length\n        N_fft = 2 * np.ceil(L/2).astype(int)        # make it even\n        K = np.round(N_fft * fc/fs).astype(int)     # Vector operation - convert fc to k (digital freq)\n        # print(K)\n        no_of_subbands = K.shape[0] - 1             # For FFT implementation\n        Hk = np.zeros((N_fft, 1, no_of_subbands))   # Create 3-D Matrix with third dimension as Fc filters\n\n    for i in range(no_of_subbands):\n        if filter_type == 'dct':\n            if sort_fc == 'ascend':\n                Hk[K[i] : K[i+1], :, i] = 1  # 0:48 slice means 0 to 47 only\n            if sort_fc == 'descend':\n                Hk[K[i+1] : K[i], :, i] = 1\n        if filter_type == 'dft':\n            if sort_fc == 'ascend':\n                Hk[K[i] : K[i+1], :, i] = 1\n                Hk[N_fft-K[i+1] : N_fft-K[i], :, i] = 1\n            if sort_fc == 'descend':\n                Hk[K[i+1] : K[i], :, i] = 1\n                Hk[N_fft-K[i] : N_fft-K[i+1], :, i] = 1\n    \n    if filter_type == 'dct':\n        Xk = dct(X, type=dct_type, n=N, axis=axis, norm='ortho', overwrite_x=False, workers=None)\n        # if X.shape[1] == 1:\n        #     Yk = Xk * Hk\n        # else:\n        #     Yk = np.einsum('ij,ijk->ijk', Xk, Hk)\n        Yk = np.einsum('ij,ijk->ijk', Xk, Hk)  # Python Broadcasting, Order of Multiplication of dimensions\n        Y = idct(Yk, type=dct_type, n=N, axis=axis, norm='ortho', overwrite_x=False, workers=None)\n    \n    if filter_type == 'dft':\n        Xk = 1/L * fft(X_padded, n=N_fft, axis=axis, norm=None)  # Divide by L because of scaling in FFT & IFFT formulae\n        # if X.shape[1] == 1:\n        #     Yk = Xk * Hk\n        # else:\n        #     Yk = np.einsum('ij,ijk->ijk', Xk, Hk)  # Python Broadcasting\n        Yk = np.einsum('ij,ijk->ijk', Xk, Hk)  # Python Broadcasting\n        Y = L * ifft(Yk, n=N_fft, axis=axis, norm=None, overwrite_x=False, workers=None)\n        Y = np.real(Y)  # Ignore the imaginary part\n        Y = Y[pad_rows+1 : appended_length-pad_rows+1, :, :]  # Remove padded data i.e., revert back to orginal size of data\n        \n    if X.shape[1] == 1:  # Output is returned as a 2D matrix if input to the function is a column vector or a row vector\n        FIBFs = np.squeeze(Y)\n    else:\n        FIBFs = Y  # Return output in 3D Matrix with the third dimension corresponding to different subbands    \n\n    if plot_subbands:\n        # plt.show()\n        t = np.arange(0, X.shape[0], 1) / fs\n        fc = np.sort(fc)  # always sort fc in ascending order\n        if fc[0] != 0:\n            fc = [0, fc]\n        if fc[-1] != fs/2:\n            fc = [fc, fs/2]\n        \n        if sort_fc == 'descend':\n            fc[::-1].sort()\n\n        no_of_subbands = len(fc) - 1\n        if no_of_subbands >= 6:\n            m = np.floor(no_of_subbands/2) + 1\n            n = 2  # Change to two-column tiled layout\n        else:\n            m = no_of_subbands + 1\n            n = 1\n\n        # First create a grid of plots\n        # ax will be an array of two Axes objects\n        # fig, ax = plt.subplots(m, n)\n        # fig, ax = plt.subplots(m, n, figsize = (15, 10))\n        # Call plot() method on the appropriate object\n        # ax[0].plot(x, np.sin(x))\n        # ax[1].plot(x, np.cos(x));\n        if n == 2:\n            plt.figure(figsize=(16, 16))\n        else:\n            plt.figure(figsize=(8, 16))\n        # print(m)\n        # print(n)\n        \n        for k in range(no_of_subbands+1):\n            plt.subplot(m, n, k+1)\n            if k == 0:\n                plt.plot(t, X[:, 0])\n                plt.title('Signal')\n            else:\n                if X.shape[1] == 1:  # Check if column vector or matrix\n                    plt.plot(t, FIBFs[:, k-1])  # Plotting just one channel of the multi channel data\n                else:\n                    plt.plot(t, FIBFs[:, 0, k-1])  # Plotting just one channel of the multi channel data\n                \n                if sort_fc == 'ascend':\n                    plt.title('FIBF {0}: Between {1} Hz to {2} Hz'.format(k, fc[k-1], fc[k]))\n                elif sort_fc == 'descend':\n                    plt.title('FIBF {0}: Between {1} Hz to {2} Hz'.format(k, fc[k], fc[k-1]))\n        # plt.title('Fourier Decomposition Method')\n\n    if data_type == 'rows':\n        if X.shape[1] == 1:  # Output is returned as a 2D matrix if input to the function is a column vector or a row vector\n            return FIBFs.T\n        else:\n            return FIBFs.transpose(1, 0, 2)  # The final 3D Matrix is a collection of multiple 2D matrices\n    return FIBFs","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-19T12:49:28.331898Z","iopub.execute_input":"2024-12-19T12:49:28.332216Z","iopub.status.idle":"2024-12-19T12:49:28.349825Z","shell.execute_reply.started":"2024-12-19T12:49:28.332183Z","shell.execute_reply":"2024-12-19T12:49:28.348964Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport tensorflow as tf\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import OneHotEncoder\nfrom scipy.signal import butter, filtfilt\nfrom scipy.fft import fft, ifft\nfrom scipy.fftpack import dct, idct\nimport h5py\nfrom tqdm import tqdm\n\ndef setup_eeg_files(base_path):\n    eeg_files = {}\n    for file in os.listdir(base_path):\n        if file.endswith('.parquet'):\n            try:\n                eeg_id = os.path.splitext(file)[0]\n                if eeg_id.isdigit():\n                    eeg_files[eeg_id] = os.path.join(base_path, file)\n            except Exception as e:\n                print(f\"Error processing file {file}: {e}\")\n    return eeg_files\n\ndef modified_fdm(X, fs, fc, data_type='columns'):\n    if data_type == 'rows':\n        X = np.transpose(X)\n    \n    X = X - np.mean(X, axis=0)\n    N = X.shape[0]\n    \n    fc = np.sort(fc)\n    if fc[0] != 0:\n        fc = np.concatenate(([0], fc))\n    if fc[-1] != fs/2:\n        fc = np.concatenate((fc, [fs/2]))\n    \n    Xk = dct(X, type=2, axis=0, norm='ortho')\n    no_of_subbands = len(fc) - 1\n    decomposed_signals = np.zeros((N, no_of_subbands, X.shape[1]))\n    \n    for i in range(no_of_subbands):\n        Hk = np.zeros_like(Xk)\n        k_low = int(np.round(2 * N * fc[i]/fs))\n        k_high = int(np.round(2 * N * fc[i+1]/fs))\n        Hk[k_low:k_high, :] = 1\n        \n        Yk = Xk * Hk\n        decomposed_signals[:, i, :] = idct(Yk, type=2, axis=0, norm='ortho')\n    \n    return decomposed_signals\n\ndef extract_features(segment):\n    \"\"\"Extract features maintaining correct dimensions.\"\"\"\n    # Time domain features\n    mean = np.mean(segment, axis=0)  # (bands, channels)\n    std = np.std(segment, axis=0)\n    var = np.var(segment, axis=0)\n    rms = np.sqrt(np.mean(segment**2, axis=0))\n    \n    # Frequency domain features\n    fft_vals = np.abs(np.fft.rfft(segment, axis=0))\n    freq_mean = np.mean(fft_vals, axis=0)\n    freq_std = np.std(fft_vals, axis=0)\n    \n    # Stack features and expand dimensions to match signal\n    features = np.stack([mean, std, var, rms, freq_mean, freq_std])  # (6, bands, channels)\n    return features\n\ndef preprocess_eeg(eeg_id, eeg_sub_id, label_offset, eeg_files):\n    filepath = eeg_files.get(str(eeg_id))\n    if not filepath:\n        return None\n    \n    try:\n        eeg_data = pd.read_parquet(filepath).values\n        start_index = int(label_offset * 100)\n        end_index = start_index + 5000\n        \n        if end_index > len(eeg_data):\n            return None\n        \n        segment = eeg_data[start_index:end_index]\n        \n        # Apply bandpass filter\n        b, a = butter(5, [0.5/(100/2), 40/(100/2)], btype='band')\n        filtered_segment = filtfilt(b, a, segment, axis=0)\n        \n        # Normalize\n        mean = np.mean(filtered_segment, axis=0)\n        std = np.std(filtered_segment, axis=0)\n        normalized_segment = (filtered_segment - mean) / (std + 1e-10)\n        \n        # Decompose signal\n        fc = np.array([0.5, 4, 8, 13, 30, 40])\n        decomposed = modified_fdm(normalized_segment, fs=100, fc=fc)  # (time, bands, channels)\n        \n        # Extract features\n        features = extract_features(decomposed)  # (features, bands, channels)\n        \n        # Reshape decomposed signal to (time, features/bands, channels)\n        # Add small time dimension to features to match decomposed signal\n        features_expanded = np.expand_dims(features, axis=0)  # (1, features, bands, channels)\n        features_tiled = np.tile(features_expanded, (decomposed.shape[0], 1, 1, 1))  # (time, features, bands, channels)\n        \n        # Combine decomposed signal and features\n        # Reshape decomposed signal to match features dimensions\n        decomposed_expanded = np.expand_dims(decomposed, axis=1)  # (time, 1, bands, channels)\n        \n        # Concatenate along the features dimension\n        final_data = np.concatenate([decomposed_expanded, features_tiled], axis=1)\n        \n        return final_data\n        \n    except Exception as e:\n        print(f\"Error processing EEG ID {eeg_id}: {str(e)}\")\n        return None\n\ndef process_batch(batch_indices, labels_df, eeg_files, output_file):\n    \"\"\"Process a batch of EEG signals and save to H5 file\"\"\"\n    with h5py.File(output_file, 'a') as h5f:\n        for idx in tqdm(batch_indices):\n            row = labels_df.iloc[idx]\n            signal = preprocess_eeg(row['eeg_id'], row['eeg_sub_id'], \n                                  row['eeg_label_offset_seconds'], eeg_files)\n            \n            if signal is not None:\n                # Create datasets if they don't exist\n                if 'X' not in h5f:\n                    h5f.create_dataset('X', data=[signal], maxshape=(None, *signal.shape), \n                                     chunks=True, compression='gzip')\n                    h5f.create_dataset('y', data=[row['expert_consensus']], \n                                     maxshape=(None,), chunks=True)\n                else:\n                    # Append new data\n                    h5f['X'].resize((h5f['X'].shape[0] + 1), axis=0)\n                    h5f['X'][-1] = signal\n                    h5f['y'].resize((h5f['y'].shape[0] + 1), axis=0)\n                    h5f['y'][-1] = row['expert_consensus']\n\ndef build_decomposition_cnn(input_shape, num_classes):\n    inputs = tf.keras.Input(shape=input_shape)\n    \n    x = tf.keras.layers.Conv2D(64, (7, 3), padding='same')(inputs)\n    x = tf.keras.layers.BatchNormalization()(x)\n    x = tf.keras.layers.Activation('relu')(x)\n    x = tf.keras.layers.MaxPooling2D((2, 1))(x)\n    \n    x = tf.keras.layers.Conv2D(128, (5, 3), padding='same')(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n    x = tf.keras.layers.Activation('relu')(x)\n    x = tf.keras.layers.MaxPooling2D((2, 1))(x)\n    \n    x = tf.keras.layers.Conv2D(256, (3, 3), padding='same')(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n    x = tf.keras.layers.Activation('relu')(x)\n    x = tf.keras.layers.GlobalAveragePooling2D()(x)\n    \n    x = tf.keras.layers.Dense(512, activation='relu')(x)\n    x = tf.keras.layers.Dropout(0.5)(x)\n    \n    x = tf.keras.layers.Dense(256, activation='relu')(x)\n    x = tf.keras.layers.Dropout(0.3)(x)\n    \n    outputs = tf.keras.layers.Dense(num_classes, activation='softmax')(x)\n    \n    return tf.keras.Model(inputs=inputs, outputs=outputs)\n\nclass EEGDataGenerator(tf.keras.utils.Sequence):\n    \"\"\"Memory-efficient data generator\"\"\"\n    def __init__(self, h5_file, batch_size=32, mode='train', validation_split=0.2):\n        self.h5_file = h5_file\n        self.batch_size = batch_size\n        self.mode = mode\n        \n        with h5py.File(self.h5_file, 'r') as h5f:\n            self.n_samples = len(h5f['y'])\n            \n        # Split indices for training/validation\n        indices = np.arange(self.n_samples)\n        split_idx = int(self.n_samples * (1 - validation_split))\n        self.indices = indices[:split_idx] if mode == 'train' else indices[split_idx:]\n        \n    def __len__(self):\n        return int(np.ceil(len(self.indices) / self.batch_size))\n    \n    def __getitem__(self, idx):\n        batch_indices = self.indices[idx * self.batch_size:(idx + 1) * self.batch_size]\n        \n        with h5py.File(self.h5_file, 'r') as h5f:\n            batch_x = h5f['X'][batch_indices]\n            batch_y = h5f['y'][batch_indices]\n            \n            # One-hot encode labels\n            encoder = OneHotEncoder(sparse_output=False)\n            batch_y = encoder.fit_transform(batch_y.reshape(-1, 1))\n            \n        return batch_x, batch_y\n\ndef main():\n    # Setup paths and load data\n    base_path = r'/kaggle/input/hms-harmful-brain-activity-classification/train_eegs'\n    eeg_files = setup_eeg_files(base_path)\n    labels_df = df1  # Assuming df1 is your labels dataframe\n    \n    # Process data in batches and save to H5 file\n    output_file = 'processed_eeg_data.h5'\n    batch_size = 100  # Adjust based on available memory\n    \n    total_samples = len(labels_df)\n    num_batches = (total_samples + batch_size - 1) // batch_size\n    \n    for i in range(num_batches):\n        start_idx = i * batch_size\n        end_idx = min((i + 1) * batch_size, total_samples)\n        batch_indices = range(start_idx, end_idx)\n        process_batch(batch_indices, labels_df, eeg_files, output_file)\n    \n    # Create data generators\n    train_gen = EEGDataGenerator(output_file, batch_size=32, mode='train')\n    val_gen = EEGDataGenerator(output_file, batch_size=32, mode='val')\n    \n    # Create and compile model\n    with h5py.File(output_file, 'r') as h5f:\n        input_shape = h5f['X'].shape[1:]\n    \n    model = build_decomposition_cnn(input_shape, num_classes=6)  # Adjust num_classes as needed\n    \n    model.compile(\n        optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),\n        loss='categorical_crossentropy',\n        metrics=['accuracy']\n    )\n    \n    callbacks = [\n        tf.keras.callbacks.EarlyStopping(\n            monitor='val_accuracy',\n            patience=15,\n            restore_best_weights=True\n        ),\n        tf.keras.callbacks.ReduceLROnPlateau(\n            monitor='val_loss',\n            factor=0.5,\n            patience=5,\n            min_lr=1e-6\n        ),\n        tf.keras.callbacks.ModelCheckpoint(\n            'best_decomposed_model.keras',\n            monitor='val_accuracy',\n            save_best_only=True\n        )\n    ]\n    \n    print(\"\\nTraining model...\")\n    history = model.fit(\n        train_gen,\n        validation_data=val_gen,\n        epochs=50,\n        callbacks=callbacks,\n        verbose=1\n    )\n    \n    model.save('eeg_decomposed_cnn_model_final.keras')\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-19T12:49:28.36626Z","iopub.execute_input":"2024-12-19T12:49:28.366459Z","execution_failed":"2024-12-19T12:56:52.9Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}