{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":102783,"databundleVersionId":12408531,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# MEG decoding for multiple face features\n\n## Task\n\nIn this notebook we'll predict multiple face features from optically pumped magnetometer magnetoencephalography (OPM-MEG) data. For every trial, we have a time series of 1001 time steps with up to 64 channels of data. From this time series, we want to classify the data in 16 classes.\n\nTraining data comes from 15 subjects (10 GByte), test data from 6 subjects (4 GByte).\n\n## Method\n\n1. We use only the 20 channels which are available for every trial. The other 44 channels are ignored.\n2. As the training and test subjects are disjunct, we need to ensure that the model generalizes to unseen subjects by cross-validating with a `GroupKFold`.\n3. We average timesteps 250:550 of every channel to get 20 features for the classifier.\n4. We use `RidgeClassifier` for classification.\n\n## References\n\n- [Competition](https://www.kaggle.com/competitions/meg-decoding-for-multiple-face-features-opm-meg)\n- Xu, W., Lyu, B., Ru, X., Li, D., Gu, W., Ma, X., Zheng, F., Li, T., Liao, P., Cheng, H., Yang, R., Song, J., Jin, Z., Li, C., He, K., & Gao, J. H. (2024). Decoding the Temporal Structures and Interactions of Multiple Face Dimensions Using Optically Pumped Magnetometer Magnetoencephalography (OPM-MEG). The Journal of neuroscience : the official journal of the Society for Neuroscience, 44(47), e2237232024. https://doi.org/10.1523/JNEUROSCI.2237-23.2024\n- [Hierarchical Data Format](https://en.wikipedia.org/wiki/Hierarchical_Data_Format)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"%%time\n# Unzip the data files\nimport os\nif not os.path.exists('test_data.h5'):\n    print('Unzipping...')\n    !cat /kaggle/input/meg-decoding-for-multiple-face-features-opm-meg/train_data.zip.??? > train_data.zip\n    !unzip train_data.zip\n    !rm train_data.zip\n    !cat /kaggle/input/meg-decoding-for-multiple-face-features-opm-meg/test_data.zip.??? > test_data.zip\n    !unzip test_data.zip\n    !rm test_data.zip\n    !ls -al\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T12:04:32.575180Z","iopub.execute_input":"2025-05-30T12:04:32.575927Z","iopub.status.idle":"2025-05-30T12:04:32.642321Z","shell.execute_reply.started":"2025-05-30T12:04:32.575899Z","shell.execute_reply":"2025-05-30T12:04:32.641317Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport h5py\nimport pandas as pd\nimport numpy as np\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import RidgeClassifier\nfrom sklearn.model_selection import cross_val_predict, GroupKFold\nfrom sklearn.metrics import accuracy_score\n\ntrain_subjects = [f\"sub{i}\" for i in range(1, 16)]\ntest_subjects = [f\"sub{i}\" for i in range(16, 22)]\n","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-05-30T12:04:32.643880Z","iopub.execute_input":"2025-05-30T12:04:32.644236Z","iopub.status.idle":"2025-05-30T12:04:34.174850Z","shell.execute_reply.started":"2025-05-30T12:04:32.644206Z","shell.execute_reply":"2025-05-30T12:04:34.173945Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We determine the subset of channels which is always available:","metadata":{}},{"cell_type":"code","source":"def determine_channel_intersection():\n    \"\"\"Determine the channels which are available for all subjects\"\"\"\n    global channel_name, channel_number\n    channel_name = {f\"MEG{i+1:02d}\" for i in range(64)}\n    with h5py.File(\"train_data.h5\", \"r\") as f:\n        for subject in train_subjects:\n            channel_name_sub = f[subject][\"ch\"][0] # shape (n_channels, )\n            channel_name_sub = {c.decode() for c in channel_name_sub}\n            channel_name = channel_name.intersection(channel_name_sub)\n    with h5py.File(\"test_data.h5\", \"r\") as f:\n        for subject in test_subjects:\n            channel_name_sub = f[subject][\"ch\"][0] # shape (n_channels, )\n            channel_name_sub = {c.decode() for c in channel_name_sub}\n            channel_name = channel_name.intersection(channel_name_sub)\n    channel_name = sorted(list(channel_name))\n    channel_number = [int(cn[3:]) for cn in channel_name]\n    print(f\"We use these {len(channel_number)} channels:\", channel_number)\n    # 20 channels: [1, 4, 7, 10, 15, 21, 22, 23, 24, 26, 28, 31, 35, 36, 37, 44, 48, 53, 54, 59]\n\ndetermine_channel_intersection()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T12:04:34.175798Z","iopub.execute_input":"2025-05-30T12:04:34.176632Z","iopub.status.idle":"2025-05-30T12:04:34.202659Z","shell.execute_reply.started":"2025-05-30T12:04:34.176605Z","shell.execute_reply":"2025-05-30T12:04:34.201725Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We now read the training data, which has 59802 samples, 20 channels, and 1001 timesteps.","metadata":{}},{"cell_type":"code","source":"%%time\ndef read_dataset(traintest):\n    \"\"\"Read the train or test dataset (all subjects, but only the intersection of channels)\"\"\"\n    X, y, group = [], [], []\n    if traintest == 'train':\n        filename, subject_list = \"train_data.h5\", train_subjects\n    else:\n        filename, subject_list = \"test_data.h5\", test_subjects\n    with h5py.File(filename, \"r\") as f:\n        for subject in subject_list:\n            print('Reading', subject)\n            channel_name_sub = f[subject][\"ch\"][0] # shape (n_channels, )\n            channel_name_sub = [c.decode() for c in channel_name_sub]\n            cols = [channel_name_sub.index(c) for c in channel_name] # subset consisting of 20 columns\n            \n            X_sub = f[subject][\"X\"][0]  # shape (n_epochs, n_channels, 1001) \n            X_sub = X_sub[:, cols, :] # select channel subset\n            X_sub = X_sub - X_sub.mean(axis=0, keepdims=True) # center every timestep for every subject\n            X.append(X_sub)\n\n            if traintest == 'train':\n                y_sub = f[subject][\"y\"][0]  # shape (n_epochs, )\n                assert y_sub.shape[0] == X_sub.shape[0]\n                y.append(y_sub)\n    \n                group.append(np.full(len(y_sub), int(subject[3:])))\n\n    X = np.vstack(X) # shape (59802, 23, 1001) \n    print('Dataset shape: (n_epochs, n_channels, n_timesteps) =', X.shape)\n    if traintest == 'train':\n        y = np.hstack(y) # shape (59802, )\n        group = np.hstack(group) # shape (59802, )\n        return X, y, group\n    else:\n        return X\n\nX, y, group = read_dataset('train')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T12:04:34.204765Z","iopub.execute_input":"2025-05-30T12:04:34.205331Z","iopub.status.idle":"2025-05-30T12:06:28.691890Z","shell.execute_reply.started":"2025-05-30T12:04:34.205306Z","shell.execute_reply":"2025-05-30T12:06:28.690889Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Cross-validation\n\nWe determine the best range for the timesteps, i.e., the timestep range which gives the highest accuracy in five-fold cross-validation.","metadata":{}},{"cell_type":"code","source":"%%time\n# Cross-validation for multiclass classification\nprint(f\"n_samples: {len(X)}\")\nprint(f\"n_features: {X.shape[1]}\")\nmodel = make_pipeline(StandardScaler(), RidgeClassifier())\nbest_acc = 0\nfor t in range(200, 700, 50):\n    for l in range(150, min(400, 1001-t), 50):\n        oof = cross_val_predict(model, X[:,:,t:t+l].mean(axis=-1), y, cv=GroupKFold(), groups=group)\n        acc = accuracy_score(y, oof)\n        # print(f\"# {t:4} {l:4} {acc=:.4f}\")\n        if acc > best_acc:\n            best_acc, best_t, best_l = acc, t, l\n\nprint(f\"# {best_t=} {best_l=} {best_acc=:.4f}\")\n# best_t=250 best_l=300 best_acc=0.0743","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T12:06:28.692949Z","iopub.execute_input":"2025-05-30T12:06:28.693276Z","iopub.status.idle":"2025-05-30T12:07:04.687914Z","shell.execute_reply.started":"2025-05-30T12:06:28.693249Z","shell.execute_reply":"2025-05-30T12:07:04.686891Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submission\n\nWe refit the model to the whole training set and compute the test predictions.","metadata":{}},{"cell_type":"code","source":"# Refit\nmodel.fit(X[:,:,best_t:best_t+best_l].mean(axis=-1), y)\n\n# Read the test dataset\nX_te = read_dataset('test')\nprint()\n\n# Prepare and write the submission file\nsubmission = pd.read_csv('/kaggle/input/meg-decoding-for-multiple-face-features-opm-meg/sample_submission.csv')\nsubmission['prediction'] = model.predict(X_te[:,:,best_t:best_t+best_l].mean(axis=-1))\nsubmission.to_csv('submission.csv', index=False)\n!head submission.csv\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T12:07:04.688972Z","iopub.execute_input":"2025-05-30T12:07:04.689334Z","iopub.status.idle":"2025-05-30T12:07:46.225071Z","shell.execute_reply.started":"2025-05-30T12:07:04.689298Z","shell.execute_reply":"2025-05-30T12:07:46.223988Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}