|
| 1 | +""".. _bcic-iv-2a-moabb-trial: |
| 2 | +
|
| 3 | +Basic Brain Decoding on EEG Data |
| 4 | +======================================== |
| 5 | +
|
| 6 | +This tutorial shows you how to train and test deep learning models with |
| 7 | +Braindecode in a classical EEG setting: you have trials of data with |
| 8 | +labels (e.g., Right Hand, Left Hand, etc.). |
| 9 | +
|
| 10 | +.. contents:: This example covers: |
| 11 | + :local: |
| 12 | + :depth: 2 |
| 13 | +
|
| 14 | +""" |
| 15 | + |
| 16 | +###################################################################### |
| 17 | +# Loading and preparing the data |
| 18 | +# ------------------------------------- |
| 19 | +# |
| 20 | + |
| 21 | + |
| 22 | +###################################################################### |
| 23 | +# Loading the dataset |
| 24 | +# ~~~~~~~~~~~~~~~~~~~~~~~ |
| 25 | +# |
| 26 | + |
| 27 | + |
| 28 | +###################################################################### |
| 29 | +# First, we load the data. In this tutorial, we load the BCI Competition |
| 30 | +# IV 2a data [1]_ using braindecode's wrapper to load via |
| 31 | +# `MOABB library <moabb_>`_ [2]_. |
| 32 | +# |
| 33 | +# .. note:: |
| 34 | +# To load your own datasets either via mne or from |
| 35 | +# preprocessed X/y numpy arrays, see :ref:`mne-dataset-example` |
| 36 | +# and :ref:`custom-dataset-example`. |
| 37 | +# |
| 38 | + |
| 39 | +from braindecode.datasets import MOABBDataset |
| 40 | + |
| 41 | +subject_id = 3 |
| 42 | +dataset = MOABBDataset(dataset_name="BNCI2014_001", subject_ids=[subject_id]) |
| 43 | + |
| 44 | + |
| 45 | +###################################################################### |
| 46 | +# Preprocessing |
| 47 | +# ~~~~~~~~~~~~~ |
| 48 | +# |
| 49 | + |
| 50 | + |
| 51 | +###################################################################### |
| 52 | +# Now we apply preprocessing like bandpass filtering to our dataset. |
| 53 | +# You can either apply functions provided by :class:`mne.io.Raw` or |
| 54 | +# :class:`mne.Epochs` or apply your own functions, either to the |
| 55 | +# MNE object or the underlying numpy array. |
| 56 | +# |
| 57 | +# .. note:: |
| 58 | +# Generally, braindecode prepocessing is directly applied to the loaded |
| 59 | +# data, and not applied on-the-fly as transformations, such as in |
| 60 | +# PyTorch-libraries like `<torchvision_>`_. |
| 61 | +# |
| 62 | + |
| 63 | +from numpy import multiply |
| 64 | + |
| 65 | +from braindecode.preprocessing import ( |
| 66 | + Preprocessor, |
| 67 | + exponential_moving_standardize, |
| 68 | + preprocess, |
| 69 | +) |
| 70 | + |
| 71 | +low_cut_hz = 4.0 # low cut frequency for filtering |
| 72 | +high_cut_hz = 38.0 # high cut frequency for filtering |
| 73 | +# Parameters for exponential moving standardization |
| 74 | +factor_new = 1e-3 |
| 75 | +init_block_size = 1000 |
| 76 | +# Factor to convert from V to uV |
| 77 | +factor = 1e6 |
| 78 | + |
| 79 | +preprocessors = [ |
| 80 | + Preprocessor("pick_types", eeg=True, meg=False, stim=False), # Keep EEG sensors |
| 81 | + Preprocessor(lambda data: multiply(data, factor)), # Convert from V to uV |
| 82 | + Preprocessor("filter", l_freq=low_cut_hz, h_freq=high_cut_hz), # Bandpass filter |
| 83 | + Preprocessor( |
| 84 | + exponential_moving_standardize, # Exponential moving standardization |
| 85 | + factor_new=factor_new, |
| 86 | + init_block_size=init_block_size, |
| 87 | + ), |
| 88 | +] |
| 89 | + |
| 90 | +# Transform the data |
| 91 | +preprocess(dataset, preprocessors, n_jobs=-1) |
| 92 | + |
| 93 | + |
| 94 | +###################################################################### |
| 95 | +# Extracting Compute Windows |
| 96 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 97 | +# |
| 98 | + |
| 99 | + |
| 100 | +###################################################################### |
| 101 | +# Now we extract compute windows from the signals, these will be the inputs |
| 102 | +# to the deep networks during training. In the case of trialwise |
| 103 | +# decoding, we just have to decide if we want to include some part |
| 104 | +# before and/or after the trial. For our work with this dataset, |
| 105 | +# it was often beneficial to also include the 500 ms before the trial. |
| 106 | +# |
| 107 | + |
| 108 | +from braindecode.preprocessing import create_windows_from_events |
| 109 | + |
| 110 | +trial_start_offset_seconds = -0.5 |
| 111 | +# Extract sampling frequency, check that they are same in all datasets |
| 112 | +sfreq = dataset.datasets[0].raw.info["sfreq"] |
| 113 | +assert all([ds.raw.info["sfreq"] == sfreq for ds in dataset.datasets]) |
| 114 | +# Calculate the trial start offset in samples. |
| 115 | +trial_start_offset_samples = int(trial_start_offset_seconds * sfreq) |
| 116 | + |
| 117 | +# Create windows using braindecode function for this. It needs parameters to define how |
| 118 | +# trials should be used. |
| 119 | +windows_dataset = create_windows_from_events( |
| 120 | + dataset, |
| 121 | + trial_start_offset_samples=trial_start_offset_samples, |
| 122 | + trial_stop_offset_samples=0, |
| 123 | + preload=True, |
| 124 | +) |
| 125 | + |
| 126 | + |
| 127 | +###################################################################### |
| 128 | +# Splitting the dataset into training and validation sets |
| 129 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 130 | +# |
| 131 | + |
| 132 | + |
| 133 | +###################################################################### |
| 134 | +# We can easily split the dataset using additional info stored in the |
| 135 | +# description attribute, in this case ``session`` column. We select |
| 136 | +# ``0train`` for training and ``1test`` for validation. |
| 137 | +# |
| 138 | + |
| 139 | +splitted = windows_dataset.split("session") |
| 140 | +train_set = splitted["0train"] # Session train |
| 141 | +valid_set = splitted["1test"] # Session evaluation |
| 142 | + |
| 143 | + |
| 144 | +###################################################################### |
| 145 | +# Creating a model |
| 146 | +# ---------------- |
| 147 | +# |
| 148 | + |
| 149 | + |
| 150 | +###################################################################### |
| 151 | +# Now we create the deep learning model! Braindecode comes with some |
| 152 | +# predefined convolutional neural network architectures for raw |
| 153 | +# time-domain EEG. Here, we use the :class:`ShallowFBCSPNet |
| 154 | +# <braindecode.models.ShallowFBCSPNet>` model from [3]_. These models are |
| 155 | +# pure `PyTorch <pytorch_>`_ deep learning models, therefore |
| 156 | +# to use your own model, it just has to be a normal PyTorch |
| 157 | +# :class:`torch.nn.Module`. |
| 158 | +# |
| 159 | + |
| 160 | +import torch |
| 161 | + |
| 162 | +from braindecode.models import ShallowFBCSPNet |
| 163 | +from braindecode.util import set_random_seeds |
| 164 | + |
| 165 | +cuda = torch.cuda.is_available() # check if GPU is available, if True chooses to use it |
| 166 | +device = "cuda" if cuda else "cpu" |
| 167 | +if cuda: |
| 168 | + torch.backends.cudnn.benchmark = True |
| 169 | +# Set random seed to be able to roughly reproduce results |
| 170 | +# Note that with cudnn benchmark set to True, GPU indeterminism |
| 171 | +# may still make results substantially different between runs. |
| 172 | +# To obtain more consistent results at the cost of increased computation time, |
| 173 | +# you can set `cudnn_benchmark=False` in `set_random_seeds` |
| 174 | +# or remove `torch.backends.cudnn.benchmark = True` |
| 175 | +seed = 20200220 |
| 176 | +set_random_seeds(seed=seed, cuda=cuda) |
| 177 | + |
| 178 | +n_classes = 4 |
| 179 | +classes = list(range(n_classes)) |
| 180 | +# Extract number of chans and time steps from dataset |
| 181 | +n_chans = train_set[0][0].shape[0] |
| 182 | +n_times = train_set[0][0].shape[1] |
| 183 | + |
| 184 | +model = ShallowFBCSPNet( |
| 185 | + n_chans, |
| 186 | + n_classes, |
| 187 | + n_times=n_times, |
| 188 | + final_conv_length="auto", |
| 189 | +) |
| 190 | + |
| 191 | +# Display torchinfo table describing the model |
| 192 | +print(model) |
| 193 | + |
| 194 | +# Send model to GPU |
| 195 | +if cuda: |
| 196 | + model = model.cuda() |
| 197 | + |
| 198 | + |
| 199 | +###################################################################### |
| 200 | +# Model Training |
| 201 | +# -------------- |
| 202 | +# |
| 203 | + |
| 204 | + |
| 205 | +###################################################################### |
| 206 | +# Now we will train the network! :class:`EEGClassifier |
| 207 | +# <braindecode.classifier.EEGClassifier>` is a Braindecode object |
| 208 | +# responsible for managing the training of neural networks. |
| 209 | +# It inherits from :class:`skorch.classifier.NeuralNetClassifier`, |
| 210 | +# so the training logic is the same as in `<skorch_>`_. |
| 211 | +# |
| 212 | + |
| 213 | + |
| 214 | +###################################################################### |
| 215 | +# .. note:: |
| 216 | +# In this tutorial, we use some default parameters that we |
| 217 | +# have found to work well for motor decoding, however we strongly |
| 218 | +# encourage you to perform your own hyperparameter optimization using |
| 219 | +# cross validation on your training data. |
| 220 | +# |
| 221 | + |
| 222 | +from skorch.callbacks import LRScheduler |
| 223 | +from skorch.helper import predefined_split |
| 224 | + |
| 225 | +from braindecode import EEGClassifier |
| 226 | + |
| 227 | +# We found these values to be good for the shallow network: |
| 228 | +lr = 0.0625 * 0.01 |
| 229 | +weight_decay = 0 |
| 230 | + |
| 231 | +# For deep4 they should be: |
| 232 | +# lr = 1 * 0.01 |
| 233 | +# weight_decay = 0.5 * 0.001 |
| 234 | + |
| 235 | +batch_size = 64 |
| 236 | +n_epochs = 4 |
| 237 | + |
| 238 | +clf = EEGClassifier( |
| 239 | + model, |
| 240 | + criterion=torch.nn.CrossEntropyLoss, |
| 241 | + optimizer=torch.optim.AdamW, |
| 242 | + train_split=predefined_split(valid_set), # using valid_set for validation |
| 243 | + optimizer__lr=lr, |
| 244 | + optimizer__weight_decay=weight_decay, |
| 245 | + batch_size=batch_size, |
| 246 | + callbacks=[ |
| 247 | + "accuracy", |
| 248 | + ("lr_scheduler", LRScheduler("CosineAnnealingLR", T_max=n_epochs - 1)), |
| 249 | + ], |
| 250 | + device=device, |
| 251 | + classes=classes, |
| 252 | +) |
| 253 | +# Model training for the specified number of epochs. ``y`` is ``None`` as it is |
| 254 | +# already supplied in the dataset. |
| 255 | +_ = clf.fit(train_set, y=None, epochs=n_epochs) |
| 256 | + |
| 257 | + |
| 258 | +###################################################################### |
| 259 | +# Plotting Results |
| 260 | +# ---------------- |
| 261 | +# |
| 262 | + |
| 263 | + |
| 264 | +###################################################################### |
| 265 | +# Now we use the history stored by skorch throughout training to plot |
| 266 | +# accuracy and loss curves. |
| 267 | +# |
| 268 | + |
| 269 | +import matplotlib.pyplot as plt |
| 270 | +import pandas as pd |
| 271 | +from matplotlib.lines import Line2D |
| 272 | + |
| 273 | +# Extract loss and accuracy values for plotting from history object |
| 274 | +results_columns = ["train_loss", "valid_loss", "train_accuracy", "valid_accuracy"] |
| 275 | +df = pd.DataFrame( |
| 276 | + clf.history[:, results_columns], |
| 277 | + columns=results_columns, |
| 278 | + index=clf.history[:, "epoch"], |
| 279 | +) |
| 280 | + |
| 281 | +# get percent of misclass for better visual comparison to loss |
| 282 | +df = df.assign( |
| 283 | + train_misclass=100 - 100 * df.train_accuracy, |
| 284 | + valid_misclass=100 - 100 * df.valid_accuracy, |
| 285 | +) |
| 286 | + |
| 287 | +fig, ax1 = plt.subplots(figsize=(8, 3)) |
| 288 | +df.loc[:, ["train_loss", "valid_loss"]].plot( |
| 289 | + ax=ax1, style=["-", ":"], marker="o", color="tab:blue", legend=False, fontsize=14 |
| 290 | +) |
| 291 | + |
| 292 | +ax1.tick_params(axis="y", labelcolor="tab:blue", labelsize=14) |
| 293 | +ax1.set_ylabel("Loss", color="tab:blue", fontsize=14) |
| 294 | + |
| 295 | +ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis |
| 296 | + |
| 297 | +df.loc[:, ["train_misclass", "valid_misclass"]].plot( |
| 298 | + ax=ax2, style=["-", ":"], marker="o", color="tab:red", legend=False |
| 299 | +) |
| 300 | +ax2.tick_params(axis="y", labelcolor="tab:red", labelsize=14) |
| 301 | +ax2.set_ylabel("Misclassification Rate [%]", color="tab:red", fontsize=14) |
| 302 | +ax2.set_ylim(ax2.get_ylim()[0], 85) # make some room for legend |
| 303 | +ax1.set_xlabel("Epoch", fontsize=14) |
| 304 | + |
| 305 | +# where some data has already been plotted to ax |
| 306 | +handles = [] |
| 307 | +handles.append( |
| 308 | + Line2D([0], [0], color="black", linewidth=1, linestyle="-", label="Train") |
| 309 | +) |
| 310 | +handles.append( |
| 311 | + Line2D([0], [0], color="black", linewidth=1, linestyle=":", label="Valid") |
| 312 | +) |
| 313 | +plt.legend(handles, [h.get_label() for h in handles], fontsize=14) |
| 314 | +plt.tight_layout() |
| 315 | + |
| 316 | + |
| 317 | +###################################################################### |
| 318 | +# Plotting a Confusion Matrix |
| 319 | +# ---------------------------- |
| 320 | +# |
| 321 | + |
| 322 | + |
| 323 | +####################################################################### |
| 324 | +# Here we generate a confusion matrix as in [3]_. |
| 325 | +# |
| 326 | + |
| 327 | + |
| 328 | +from sklearn.metrics import confusion_matrix |
| 329 | + |
| 330 | +from braindecode.visualization import plot_confusion_matrix |
| 331 | + |
| 332 | +# generate confusion matrices |
| 333 | +# get the targets |
| 334 | +y_true = valid_set.get_metadata().target |
| 335 | +y_pred = clf.predict(valid_set) |
| 336 | + |
| 337 | +# generating confusion matrix |
| 338 | +confusion_mat = confusion_matrix(y_true, y_pred) |
| 339 | + |
| 340 | +# add class labels |
| 341 | +# label_dict is class_name : str -> i_class : int |
| 342 | +label_dict = windows_dataset.datasets[0].window_kwargs[0][1]["mapping"] |
| 343 | +# sort the labels by values (values are integer class labels) |
| 344 | +labels = [k for k, v in sorted(label_dict.items(), key=lambda kv: kv[1])] |
| 345 | + |
| 346 | +# plot the basic conf. matrix |
| 347 | +plot_confusion_matrix(confusion_mat, class_names=labels) |
| 348 | + |
| 349 | +############################################################# |
| 350 | +# |
| 351 | +# |
| 352 | +# References |
| 353 | +# ---------- |
| 354 | +# |
| 355 | +# .. [1] Tangermann, M., Müller, K.R., Aertsen, A., Birbaumer, N., Braun, C., |
| 356 | +# Brunner, C., Leeb, R., Mehring, C., Miller, K.J., Mueller-Putz, G. |
| 357 | +# and Nolte, G., 2012. Review of the BCI competition IV. |
| 358 | +# Frontiers in neuroscience, 6, p.55. |
| 359 | +# |
| 360 | +# .. [2] Jayaram, Vinay, and Alexandre Barachant. |
| 361 | +# "MOABB: trustworthy algorithm benchmarking for BCIs." |
| 362 | +# Journal of neural engineering 15.6 (2018): 066011. |
| 363 | +# |
| 364 | +# .. [3] Schirrmeister, R.T., Springenberg, J.T., Fiederer, L.D.J., Glasstetter, M., |
| 365 | +# Eggensperger, K., Tangermann, M., Hutter, F., Burgard, W. and Ball, T. (2017), |
| 366 | +# Deep learning with convolutional neural networks for EEG decoding and visualization. |
| 367 | +# Hum. Brain Mapping, 38: 5391-5420. https://doi.org/10.1002/hbm.23730. |
| 368 | +# |
| 369 | +# .. include:: /links.inc |
0 commit comments