{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Abdelrahman Yehia\n# Mahmoud Hossam","metadata":{}},{"cell_type":"code","source":"# Imports\nimport numpy as np  # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport os\nfrom scipy import signal \nfrom scipy.io import loadmat\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import classification_report\nfrom sklearn.svm import SVC\nfrom tensorflow.keras.models  import Sequential\nfrom tensorflow.keras.layers import Conv2D, Activation, Dropout, Dense, MaxPooling2D, Flatten\nfrom tensorflow.keras.optimizers import SGD, Adam, Adadelta\nimport tensorflow as tf","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:14:12.700178Z","iopub.execute_input":"2022-02-08T18:14:12.700530Z","iopub.status.idle":"2022-02-08T18:14:19.434171Z","shell.execute_reply.started":"2022-02-08T18:14:12.700424Z","shell.execute_reply":"2022-02-08T18:14:19.433361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extracting the MAT files\n!unzip \"/kaggle/input/decoding-the-human-brain/train_01_06.zip\"  -d ./train\n!unzip \"/kaggle/input/decoding-the-human-brain/train_07_12.zip\"  -d ./train\n!unzip \"/kaggle/input/decoding-the-human-brain/train_13_16.zip\"  -d ./train\n!unzip \"/kaggle/input/decoding-the-human-brain/test_17_23.zip\"  -d ./test","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:14:19.435691Z","iopub.execute_input":"2022-02-08T18:14:19.435910Z","iopub.status.idle":"2022-02-08T18:15:55.641875Z","shell.execute_reply.started":"2022-02-08T18:14:19.435879Z","shell.execute_reply":"2022-02-08T18:15:55.640868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A walkaround the data\nThe training data contains 9414 trials from 16 participants.","metadata":{}},{"cell_type":"code","source":"# Run only if features.csv is not found or when implementing new features\n# Read MAT files\ntrain_mats = []\nfor dirname, _, filenames in os.walk('/kaggle/working/train/data'):\n    for filename in filenames:\n        print(\"Loading: \", os.path.join(dirname, filename))\n        train_mats.append(loadmat(os.path.join(dirname, filename)))\n\n# Merging data\ntrain_data = train_mats[0]['X']\ntrain_labl = train_mats[0]['y']\nfor i in range(1, len(train_mats)):\n    train_data = np.concatenate((train_data, train_mats[i]['X']), axis=0)\n    train_labl = np.concatenate((train_labl, train_mats[i]['y']), axis=0)\n\ndel train_mats # Free RAM pls \nprint('Done')\nprint(train_data.shape, type(train_data))\nprint(train_labl.shape, type(train_labl))","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:15:55.646228Z","iopub.execute_input":"2022-02-08T18:15:55.646489Z","iopub.status.idle":"2022-02-08T18:16:31.734247Z","shell.execute_reply.started":"2022-02-08T18:15:55.646457Z","shell.execute_reply":"2022-02-08T18:16:31.733290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the signal from the first channel from the first trail\nplt.figure(figsize=(16,4))\nplt.plot(train_data[0][0])\nplt.xlabel('Time')\nplt.ylabel('Value')\n","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:31.736642Z","iopub.execute_input":"2022-02-08T18:16:31.736871Z","iopub.status.idle":"2022-02-08T18:16:32.061127Z","shell.execute_reply.started":"2022-02-08T18:16:31.736840Z","shell.execute_reply":"2022-02-08T18:16:32.060440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The power spectral density (PSD) of the signal\nf,psd = signal.welch(train_data[0][0], fs=250.0)\nplt.figure(figsize=(16,4))\nplt.semilogy(f, psd)\nplt.xlabel('frequency [Hz]')\nplt.ylabel('PSD [V**2/Hz]')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:32.062082Z","iopub.execute_input":"2022-02-08T18:16:32.063023Z","iopub.status.idle":"2022-02-08T18:16:32.764654Z","shell.execute_reply.started":"2022-02-08T18:16:32.062970Z","shell.execute_reply":"2022-02-08T18:16:32.763680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing\nThe data description mentions that all preprocessing tasks were applied using mne-tools. However, we will just make some simple preprocessing tasks from a machine learning prospective.","metadata":{}},{"cell_type":"code","source":"# Is there any Null Values?\nprint(np.isnan(np.sum(train_data)))\nprint(np.isnan(np.sum(train_labl)))","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:32.766312Z","iopub.execute_input":"2022-02-08T18:16:32.766845Z","iopub.status.idle":"2022-02-08T18:16:33.286625Z","shell.execute_reply.started":"2022-02-08T18:16:32.766799Z","shell.execute_reply":"2022-02-08T18:16:33.285587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Wrong Labels?\nlabels = []\nfor label in train_labl:\n    if label not in labels:\n        labels.append(label)\nprint(labels)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.288277Z","iopub.execute_input":"2022-02-08T18:16:33.288826Z","iopub.status.idle":"2022-02-08T18:16:33.311775Z","shell.execute_reply.started":"2022-02-08T18:16:33.288779Z","shell.execute_reply":"2022-02-08T18:16:33.310753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Filtering","metadata":{}},{"cell_type":"code","source":"samp_freq = 1000 # Sample Frequency(Hz)\nnotch_freq = 50 # Frequency to be removed from the signal (Hz)\nquality_factor = 20.0\n\n# Design a notch filter using signal.iirnotch\nb_notch, a_notch = signal.iirnotch(notch_freq, quality_factor, samp_freq)\n\n# Compute magnitude response of the designed filter\nfreq, h = signal.freqz(b_notch, a_notch, fs=samp_freq)\n  \nfig = plt.figure(figsize=(8, 6))\n  \n# Plot magnitude response of the filter\nplt.plot(freq*samp_freq/(2*np.pi), 20 * np.log10(abs(h)),\n         'r', label='Bandpass filter', linewidth='2')\nplt.xlabel('Frequency [Hz]', fontsize=20)\nplt.ylabel('Magnitude [dB]', fontsize=20)\nplt.title('Notch Filter', fontsize=20)\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.313694Z","iopub.execute_input":"2022-02-08T18:16:33.314047Z","iopub.status.idle":"2022-02-08T18:16:33.561277Z","shell.execute_reply.started":"2022-02-08T18:16:33.313999Z","shell.execute_reply":"2022-02-08T18:16:33.560370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.signal import butter, lfilter,iirnotch, medfilt\ndef apply_filters(data):\n    # Low-pass filter 100Hz\n    low = 0.8  # low = 100Hz / (250Hz/2)\n    b_low, a_low = butter(5, low, btype='low') # 5th Order filter\n    \n    # Notch filter 50Hz (Powerline freq)\n    b_not, a_not = iirnotch(50.0, 30.0, 250.0)\n    \n    filtered = lfilter(b_low, a_low, data)     # Applying Low-pass filter\n    filtered = lfilter(b_not, a_not, filtered) # Applying Notch filter\n    #iltered = medfilt(filtered[t][c])         # Applying Median filter\n    \n    return filtered","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.562480Z","iopub.execute_input":"2022-02-08T18:16:33.562773Z","iopub.status.idle":"2022-02-08T18:16:33.569483Z","shell.execute_reply.started":"2022-02-08T18:16:33.562741Z","shell.execute_reply":"2022-02-08T18:16:33.568575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting\nn = np.linspace(0, 1, 375)\nfig = plt.figure(figsize=(8, 6))\nplt.subplot(211)\nplt.plot(n,train_data[0][0], color='r', linewidth=2)\nplt.xlabel('Time', fontsize=20)\nplt.ylabel('Magnitude', fontsize=18)\nplt.title('Original Signal', fontsize=20)\n\n# Apply notch filter to the noisy signal using signal.filtfilt\noutputSignal = apply_filters(train_data[0][0])\n\n# Plot notch-filtered version of signal\nplt.subplot(212)\n  \n# Plot output signal of notch filter\nplt.plot(n, outputSignal)\nplt.xlabel('Time', fontsize=20)\nplt.ylabel('Magnitude', fontsize=18)\nplt.title('Filtered Signal', fontsize=20)\nplt.subplots_adjust(hspace=0.5)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.574085Z","iopub.execute_input":"2022-02-08T18:16:33.574736Z","iopub.status.idle":"2022-02-08T18:16:33.925450Z","shell.execute_reply.started":"2022-02-08T18:16:33.574695Z","shell.execute_reply":"2022-02-08T18:16:33.924489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Extraction","metadata":{}},{"cell_type":"code","source":"def create_features(XX, tmin, tmax, sfreq, tmin_original=-0.5):\n    \"\"\"Creation of the feature space:\n    - restricting the time window of MEG data to [tmin, tmax]sec.\n    - Concatenating the 306 timeseries of each trial in one long\n      vector.\n    - Normalizing each feature independently (z-scoring).\n    \"\"\"\n    nsamples, nx, ny = XX.shape\n\n    #Applying 50Hz notch filter on all channels\n    x1_filt50hz = np.empty((nsamples,nx-162, ny))\n    print(x1_filt50hz.shape)\n    for i in range(nsamples):\n        for j in range(nx-162):\n            x1_filt50hz[i][j] = apply_filters(XX[i][162+j])#signal.filtfilt(b_notch, a_notch, XX[i][j])\n    \n    print(\"Applying the desired time window.\")\n    \n    beginning = np.round((tmin - tmin_original) * sfreq).astype(np.int)\n    end = np.round((tmax - tmin_original) * sfreq).astype(np.int)\n    XX = x1_filt50hz[:,:, beginning:end].copy()\n\n    print(\"2D Reshaping: concatenating all 306 timeseries.\")\n    XX = XX.reshape(XX.shape[0], XX.shape[1] * XX.shape[2])\n\n    print(\"Features Normalization.\")\n    XX -= XX.mean(0)\n    XX = np.nan_to_num(XX / XX.std(0))\n\n    return XX","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.926837Z","iopub.execute_input":"2022-02-08T18:16:33.927113Z","iopub.status.idle":"2022-02-08T18:16:33.937366Z","shell.execute_reply.started":"2022-02-08T18:16:33.927078Z","shell.execute_reply":"2022-02-08T18:16:33.936713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"DecMeg2014: https://www.kaggle.com/c/decoding-the-human-brain\")\nsubjects_train = range(1, 17) # use range(1, 17) for all subjects\nprint(\"Training on subjects\", subjects_train)\n\n# We throw away all the MEG data outside the first 0.5sec from when\n# the visual stimulus start:\ntmin = 0.0\ntmax = 0.500\nprint(\"Restricting MEG data to the interval [%s, %s]sec.\" % (tmin, tmax))\n\nX_train = []\ny_train = []\nX_test = []\nids_test = []\n\nsubject_samples = []\nprint(\"Creating the trainset.\")\nfor subject in subjects_train:\n    filename = '/kaggle/working/train/data/train_subject%02d.mat' % subject\n    print(\"Loading\", filename)\n    data = loadmat(filename, squeeze_me=True)\n    XX = data['X']\n    yy = data['y']\n    sfreq = data['sfreq']\n    tmin_original = data['tmin']\n    print(\"Dataset summary:\")\n    print(\"XX:\", XX.shape)\n    print(\"yy:\", yy.shape)\n    print(\"sfreq:\", sfreq)\n\n    XX = create_features(XX, tmin, tmax, sfreq)\n\n    X_train.append(XX)\n    y_train.append(yy)\n\nX_train = np.vstack(X_train)\ny_train = np.concatenate(y_train)\nprint(\"Trainset:\", X_train.shape)\n\nprint(\"Creating the testset.\")\nsubjects_test = range(17, 24)\nfor subject in subjects_test:\n    filename = '/kaggle/working/test/data/test_subject%02d.mat' % subject\n    print(\"Loading\", filename)\n    data = loadmat(filename, squeeze_me=True)\n    XX = data['X']\n    ids = data['Id']\n    sfreq = data['sfreq']\n    tmin_original = data['tmin']\n    print(\"Dataset summary:\")\n    print(\"XX:\", XX.shape)\n    print(\"ids:\", ids.shape)\n    print(\"sfreq:\", sfreq)\n\n    XX = create_features(XX, tmin, tmax, sfreq)\n\n    X_test.append(XX)\n    ids_test.append(ids)\n\nX_test = np.vstack(X_test)\nids_test = np.concatenate(ids_test)\nprint(\"Testset:\", X_test.shape)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:16:33.938459Z","iopub.execute_input":"2022-02-08T18:16:33.939139Z","iopub.status.idle":"2022-02-08T18:31:18.405499Z","shell.execute_reply.started":"2022-02-08T18:16:33.939102Z","shell.execute_reply":"2022-02-08T18:31:18.404049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(X_train.shape)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:31:18.408326Z","iopub.execute_input":"2022-02-08T18:31:18.408647Z","iopub.status.idle":"2022-02-08T18:31:18.416658Z","shell.execute_reply.started":"2022-02-08T18:31:18.408605Z","shell.execute_reply":"2022-02-08T18:31:18.415055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(X_train.shape)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:31:18.418341Z","iopub.execute_input":"2022-02-08T18:31:18.418602Z","iopub.status.idle":"2022-02-08T18:31:18.431069Z","shell.execute_reply.started":"2022-02-08T18:31:18.418567Z","shell.execute_reply":"2022-02-08T18:31:18.429990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Training","metadata":{}},{"cell_type":"code","source":"clf = SVC()\nprint(\"Training.\")\nclf.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T18:31:18.432631Z","iopub.execute_input":"2022-02-08T18:31:18.432862Z","iopub.status.idle":"2022-02-08T19:01:17.890318Z","shell.execute_reply.started":"2022-02-08T18:31:18.432832Z","shell.execute_reply":"2022-02-08T19:01:17.889163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Predicting.\")\ny_pred = clf.predict(X_test)","metadata":{"execution":{"iopub.status.busy":"2022-02-08T19:01:17.892238Z","iopub.execute_input":"2022-02-08T19:01:17.892961Z","iopub.status.idle":"2022-02-08T19:14:09.789587Z","shell.execute_reply.started":"2022-02-08T19:01:17.892919Z","shell.execute_reply":"2022-02-08T19:14:09.787948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(classification_report(y_train,np.round(clf.predict(X_train))))\n\nfilename_submission = \"Submission.csv\"\nprint(filename_submission)\nprint(\"Creating submission file\", filename_submission) \nf = open(filename_submission, \"w\")\nf.write(\"Id,Prediction\\n\")\nfor i in range(len(y_pred)):\n    f.write(str(ids_test[i]) + \",\" + str(np.round(y_pred[i])) + \"\\n\")\n\nf.close()\n# print(\"Done.\")","metadata":{"execution":{"iopub.status.busy":"2022-02-08T19:14:09.794533Z","iopub.execute_input":"2022-02-08T19:14:09.794894Z","iopub.status.idle":"2022-02-08T19:43:20.214658Z","shell.execute_reply.started":"2022-02-08T19:14:09.794851Z","shell.execute_reply":"2022-02-08T19:43:20.213675Z"},"trusted":true},"execution_count":null,"outputs":[]}]}