diff --git a/moabb/analysis/plotting.py b/moabb/analysis/plotting.py
index 1e3e3f038..a97f02b00 100644
--- a/moabb/analysis/plotting.py
+++ b/moabb/analysis/plotting.py
@@ -46,7 +46,7 @@ def score_plot(data, pipelines=None):
     ax.axvline(0.5, linestyle='--', color='k', linewidth=2)
     ax.set_title('Scores per dataset and algorithm')
     handles, labels = ax.get_legend_handles_labels()
-    color_dict = {l: h.get_facecolor()[0] for l, h in zip(labels, handles)}
+    color_dict = {lb: h.get_facecolor()[0] for lb, h in zip(labels, handles)}
     plt.tight_layout()
     return fig, color_dict
 
@@ -105,8 +105,8 @@ def summary_plot(sig_df, effect_df, p_threshold=0.05, simplify=True):
                 fmt='', cmap=palette, linewidths=1,
                 linecolor='0.8', annot_kws={'size': 10}, cbar=False,
                 vmin=-np.log(0.05), vmax=-np.log(1e-100))
-    for l in ax.get_xticklabels():
-        l.set_rotation(45)
+    for lb in ax.get_xticklabels():
+        lb.set_rotation(45)
     ax.tick_params(axis='y', rotation=0.9)
     ax.set_title("Algorithm comparison")
     plt.tight_layout()
diff --git a/moabb/analysis/results.py b/moabb/analysis/results.py
index 7222d521f..e70ecd154 100644
--- a/moabb/analysis/results.py
+++ b/moabb/analysis/results.py
@@ -39,7 +39,7 @@ class Results:
     '''
 
     def __init__(self, evaluation_class, paradigm_class, suffix='',
-                 overwrite=False):
+                 overwrite=False, hdf5_path=None):
         """
         class that will abstract result storage
         """
@@ -49,8 +49,11 @@ class that will abstract result storage
         assert issubclass(evaluation_class, BaseEvaluation)
         assert issubclass(paradigm_class, BaseParadigm)
 
-        self.mod_dir = os.path.dirname(
-            os.path.abspath(inspect.getsourcefile(moabb)))
+        if hdf5_path is None:
+            self.mod_dir = os.path.dirname(
+                os.path.abspath(inspect.getsourcefile(moabb)))
+        else:
+            self.mod_dir = os.path.abspath(hdf5_path)
         self.filepath = os.path.join(self.mod_dir, 'results',
                                      paradigm_class.__name__,
                                      evaluation_class.__name__,
diff --git a/moabb/datasets/Weibo2014.py b/moabb/datasets/Weibo2014.py
index 14615ccca..4f974520c 100644
--- a/moabb/datasets/Weibo2014.py
+++ b/moabb/datasets/Weibo2014.py
@@ -127,7 +127,7 @@ def _get_single_subject_data(self, subject):
         ch_types[61] = 'misc'
         info = mne.create_info(ch_names=ch_names + ['STIM014'],
                                ch_types=ch_types + ['stim'],
-                               sfreq=200, montage=None)
+                               sfreq=200)
         # until we get the channel names montage is None
         event_ids = data['label'].ravel()
         raw_data = np.transpose(data['data'], axes=[2, 0, 1])
diff --git a/moabb/datasets/bbci_eeg_fnirs.py b/moabb/datasets/bbci_eeg_fnirs.py
index 5cfde5f24..20f37216e 100644
--- a/moabb/datasets/bbci_eeg_fnirs.py
+++ b/moabb/datasets/bbci_eeg_fnirs.py
@@ -129,8 +129,9 @@ def _convert_one_session(self, data, mrk, session, trig_offset=0):
 
         montage = make_standard_montage('standard_1005')
         info = create_info(ch_names=ch_names, ch_types=ch_types,
-                           sfreq=200., montage=montage)
+                           sfreq=200.)
         raw = RawArray(data=eeg, info=info, verbose=False)
+        raw.set_montage(montage)
         return {'run_0': raw}
 
     def data_path(self, subject, path=None, force_update=False,
diff --git a/moabb/datasets/bnci.py b/moabb/datasets/bnci.py
index ffd94bc88..105daf662 100644
--- a/moabb/datasets/bnci.py
+++ b/moabb/datasets/bnci.py
@@ -299,8 +299,7 @@ def _load_data_003_2015(subject,
     ch_types = ['eeg'] * 8 + ['stim'] * 2
     montage = make_standard_montage('standard_1005')
 
-    info = create_info(
-        ch_names=ch_names, ch_types=ch_types, sfreq=sfreq, montage=montage)
+    info = create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
 
     sessions = {}
     sessions['session_0'] = {}
@@ -325,6 +324,7 @@ def _load_data_003_2015(subject,
 
         eeg_data = np.r_[run[1:-2] * 1e-6, targets, flashs]
         raw = RawArray(data=eeg_data, info=info, verbose=verbose)
+        raw.set_montage(montage)
         sessions['session_0']['run_' + str(ri)] = raw
 
     return sessions
@@ -531,9 +531,9 @@ def _convert_run(run, ch_names=None, ch_types=None, verbose=None):
     ch_names = ch_names + ['stim']
     ch_types = ch_types + ['stim']
     event_id = {ev: (ii + 1) for ii, ev in enumerate(run.classes)}
-    info = create_info(
-        ch_names=ch_names, ch_types=ch_types, sfreq=sfreq, montage=montage)
+    info = create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
     raw = RawArray(data=eeg_data.T, info=info, verbose=verbose)
+    raw.set_montage(montage)
     return raw, event_id
 
 
@@ -551,9 +551,9 @@ def _convert_run_p300_sl(run, verbose=None):
     eeg_data = np.c_[eeg_data, run.y, flash_stim]
     event_id = {ev: (ii + 1) for ii, ev in enumerate(run.classes)}
     event_id.update({ev: (ii + 3) for ii, ev in enumerate(run.classes_stim)})
-    info = create_info(
-        ch_names=ch_names, ch_types=ch_types, sfreq=sfreq, montage=montage)
+    info = create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
     raw = RawArray(data=eeg_data.T, info=info, verbose=verbose)
+    raw.set_montage(montage)
     return raw, event_id
 
 
@@ -596,9 +596,9 @@ def _convert_run_bbci(run, ch_types, verbose=None):
     ch_names = ch_names + ['Target', 'Flash']
     ch_types = ch_types + ['stim'] * 2
 
-    info = create_info(
-        ch_names=ch_names, ch_types=ch_types, sfreq=sfreq, montage=montage)
+    info = create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
     raw = RawArray(data=eeg_data.T, info=info, verbose=verbose)
+    raw.set_montage(montage)
     return raw, event_id
 
 
@@ -628,9 +628,9 @@ def _convert_run_epfl(run, verbose=None):
     ch_types = ch_types + ['stim']
     event_id = {'correct': 1, 'error': 2}
 
-    info = create_info(
-        ch_names=ch_names, ch_types=ch_types, sfreq=sfreq, montage=montage)
+    info = create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
     raw = RawArray(data=eeg_data.T, info=info, verbose=verbose)
+    raw.set_montage(montage)
     return raw, event_id
 
 
diff --git a/moabb/datasets/braininvaders.py b/moabb/datasets/braininvaders.py
index 7e4956a10..2fb0fd7d5 100644
--- a/moabb/datasets/braininvaders.py
+++ b/moabb/datasets/braininvaders.py
@@ -169,7 +169,7 @@ def data_path(self, subject, path=None, force_update=False,
         meta_file = os.path.join('subject{:d}'.format(subject), 'meta.yml')
         meta_path = path_folder + meta_file
         with open(meta_path, 'r') as stream:
-            meta = yaml.load(stream)
+            meta = yaml.load(stream, Loader=yaml.FullLoader)
         conditions = []
         if self.adaptive:
             conditions = conditions + ['adaptive']
diff --git a/moabb/datasets/epfl.py b/moabb/datasets/epfl.py
index 1d2aff0bd..9b90dd560 100644
--- a/moabb/datasets/epfl.py
+++ b/moabb/datasets/epfl.py
@@ -5,6 +5,7 @@
 import datetime as dt
 from moabb.datasets.base import BaseDataset
 from moabb.datasets import download as dl
+from mne.channels import make_standard_montage
 from scipy.io import loadmat
 import zipfile
 
@@ -115,6 +116,13 @@ def _get_single_run_data(self, file_path):
             'MA2']
         ch_types = ['eeg'] * 32 + ['misc'] * 2
 
+        # The last X entries are 0 for all signals. This leads to
+        # artifacts when epoching and band-pass filtering the data.
+        # Correct the signals for this.
+        sig_i = np.where(
+            np.diff(np.all(signals == 0, axis=0).astype(int)) != 0)[0][0]
+        signals = signals[:, :sig_i]
+        signals *= 1e-6  # data is stored as uV, but MNE expects V
         # we have to re-reference the signals
         # the average signal on the mastoids electrodes is used as reference
         references = [32, 33]
@@ -149,11 +157,13 @@ def _get_single_run_data(self, file_path):
         signals = np.concatenate([signals, stim_channel[None, :]])
 
         # create info dictionary
-        info = mne.create_info(ch_names, sfreq, ch_types, montage='biosemi32')
+        info = mne.create_info(ch_names, sfreq, ch_types)
         info['description'] = 'EPFL P300 dataset'
 
         # create the Raw structure
         raw = mne.io.RawArray(signals, info, verbose=False)
+        montage = make_standard_montage('biosemi32')
+        raw.set_montage(montage)
 
         return raw
 
diff --git a/moabb/datasets/fake.py b/moabb/datasets/fake.py
index bd83805b0..38eb6ba3b 100644
--- a/moabb/datasets/fake.py
+++ b/moabb/datasets/fake.py
@@ -48,8 +48,9 @@ def _generate_raw(self):
         eeg_data = np.c_[eeg_data, y]
 
         info = create_info(ch_names=ch_names, ch_types=ch_types,
-                           sfreq=sfreq, montage=montage)
+                           sfreq=sfreq)
         raw = RawArray(data=eeg_data.T, info=info, verbose=False)
+        raw.set_montage(montage)
         return raw
 
     def data_path(self, subject, path=None, force_update=False,
diff --git a/moabb/datasets/gigadb.py b/moabb/datasets/gigadb.py
index b07c5b73b..7bf226bf1 100644
--- a/moabb/datasets/gigadb.py
+++ b/moabb/datasets/gigadb.py
@@ -105,8 +105,9 @@ def _get_single_subject_data(self, subject):
                     "continuous data -- edge effects present")
 
         info = create_info(ch_names=ch_names, ch_types=ch_types,
-                           sfreq=data.srate, montage=montage)
+                           sfreq=data.srate)
         raw = RawArray(data=eeg_data, info=info, verbose=False)
+        raw.set_montage(montage)
 
         return {'session_0': {'run_0': raw}}
 
diff --git a/moabb/datasets/schirrmeister2017.py b/moabb/datasets/schirrmeister2017.py
index 451ca3f56..40207e542 100644
--- a/moabb/datasets/schirrmeister2017.py
+++ b/moabb/datasets/schirrmeister2017.py
@@ -198,8 +198,9 @@ def get_all_sensors(filename, pattern=None):
         """
         with h5py.File(filename, 'r') as h5file:
             clab_set = h5file['nfo']['clab'][:].squeeze()
-            all_sensor_names = [''.join(chr(c) for c in h5file[obj_ref]) for
-                                obj_ref in clab_set]
+            all_sensor_names = [''.join(
+                chr(c.squeeze()) for c in h5file[obj_ref])
+                                for obj_ref in clab_set]
             if pattern is not None:
                 all_sensor_names = filter(
                     lambda sname: re.search(pattern, sname),
diff --git a/moabb/evaluations/base.py b/moabb/evaluations/base.py
index f690611a2..a3d41a7f7 100644
--- a/moabb/evaluations/base.py
+++ b/moabb/evaluations/base.py
@@ -33,10 +33,12 @@ class BaseEvaluation(ABC):
     '''
 
     def __init__(self, paradigm, datasets=None, random_state=None, n_jobs=1,
-                 overwrite=False, error_score='raise', suffix=''):
+                 overwrite=False, error_score='raise', suffix='',
+                 hdf5_path=None):
         self.random_state = random_state
         self.n_jobs = n_jobs
         self.error_score = error_score
+        self.hdf5_path = hdf5_path
 
         # check paradigm
         if not isinstance(paradigm, BaseParadigm):
@@ -82,7 +84,8 @@ def __init__(self, paradigm, datasets=None, random_state=None, n_jobs=1,
         self.results = Results(type(self),
                                type(self.paradigm),
                                overwrite=overwrite,
-                               suffix=suffix)
+                               suffix=suffix,
+                               hdf5_path=self.hdf5_path)
 
     def process(self, pipelines):
         '''Runs all pipelines on all datasets.
diff --git a/moabb/paradigms/base.py b/moabb/paradigms/base.py
index 9292a1360..5a8f0eb47 100644
--- a/moabb/paradigms/base.py
+++ b/moabb/paradigms/base.py
@@ -62,7 +62,7 @@ def prepare_process(self, dataset):
         """
         pass
 
-    def process_raw(self, raw, dataset):
+    def process_raw(self, raw, dataset, return_epochs=False):
         """
         Process one raw data file.
 
@@ -83,10 +83,16 @@ def process_raw(self, raw, dataset):
             The dataset corresponding to the raw file. mainly use to access
             dataset specific information.
 
+        return_epochs: boolean
+            This flag specifies whether to return only the data array or the
+            complete processed mne.Epochs
+
         returns
         -------
-        X : np.ndarray
+        X : Union[np.ndarray, mne.Epochs]
             the data that will be used as features for the model
+            Note: if return_epochs=True,  this is mne.Epochs
+                  if return_epochs=False, this is np.ndarray
 
         labels: np.ndarray
             the labels for training / evaluating the model
@@ -95,22 +101,22 @@ def process_raw(self, raw, dataset):
             A dataframe containing the metadata
 
         """
+        # get events id
+        event_id = self.used_events(dataset)
+
         # find the events, first check stim_channels then annotations
-        stim_channels = mne.utils._get_stim_channel(
-            None, raw.info, raise_error=False)
+        stim_channels = mne.utils._get_stim_channel(None, raw.info,
+                                                    raise_error=False)
         if len(stim_channels) > 0:
             events = mne.find_events(raw, shortest_event=0, verbose=False)
         else:
             events, _ = mne.events_from_annotations(raw, verbose=False)
 
-        channels = () if self.channels is None else self.channels
-
         # picks channels
-        picks = mne.pick_types(raw.info, eeg=True, stim=False,
-                               include=channels)
-
-        # get events id
-        event_id = self.used_events(dataset)
+        if self.channels is None:
+            picks = mne.pick_types(raw.info, eeg=True, stim=False)
+        else:
+            picks = mne.pick_types(raw.info, stim=False, include=self.channels)
 
         # pick events, based on event_id
         try:
@@ -137,11 +143,15 @@ def process_raw(self, raw, dataset):
                                 tmin=tmin, tmax=tmax, proj=False,
                                 baseline=None, preload=True,
                                 verbose=False, picks=picks,
+                                event_repeated='drop',
                                 on_missing='ignore')
             if self.resample is not None:
                 epochs = epochs.resample(self.resample)
             # rescale to work with uV
-            X.append(dataset.unit_factor * epochs.get_data())
+            if return_epochs:
+                X.append(epochs)
+            else:
+                X.append(dataset.unit_factor * epochs.get_data())
 
         inv_events = {k: v for v, k in event_id.items()}
         labels = np.array([inv_events[e] for e in epochs.events[:, -1]])
@@ -155,7 +165,7 @@ def process_raw(self, raw, dataset):
         metadata = pd.DataFrame(index=range(len(labels)))
         return X, labels, metadata
 
-    def get_data(self, dataset, subjects=None):
+    def get_data(self, dataset, subjects=None, return_epochs=False):
         """
         Return the data for a list of subject.
 
@@ -197,7 +207,8 @@ def get_data(self, dataset, subjects=None):
         for subject, sessions in data.items():
             for session, runs in sessions.items():
                 for run, raw in runs.items():
-                    proc = self.process_raw(raw, dataset)
+                    proc = self.process_raw(raw, dataset,
+                                            return_epochs=return_epochs)
 
                     if proc is None:
                         # this mean the run did not contain any selected event
@@ -211,12 +222,15 @@ def get_data(self, dataset, subjects=None):
                     metadata.append(met)
 
                     # grow X and labels in a memory efficient way. can be slow
-                    if len(X) > 0:
-                        X = np.append(X, x, axis=0)
-                        labels = np.append(labels, lbs, axis=0)
+                    if not return_epochs:
+                        if len(X) > 0:
+                            X = np.append(X, x, axis=0)
+                            labels = np.append(labels, lbs, axis=0)
+                        else:
+                            X = x
+                            labels = lbs
                     else:
-                        X = x
-                        labels = lbs
+                        X.append(x)
 
         metadata = pd.concat(metadata, ignore_index=True)
         return X, labels, metadata
diff --git a/moabb/paradigms/p300.py b/moabb/paradigms/p300.py
index 4e21c91e3..e6bf29551 100644
--- a/moabb/paradigms/p300.py
+++ b/moabb/paradigms/p300.py
@@ -79,7 +79,7 @@ def is_valid(self, dataset):
     def used_events(self, dataset):
         pass
 
-    def process_raw(self, raw, dataset):
+    def process_raw(self, raw, dataset, return_epochs=False):
         # find the events, first check stim_channels then annotations
         stim_channels = mne.utils._get_stim_channel(
             None, raw.info, raise_error=False)
@@ -126,7 +126,10 @@ def process_raw(self, raw, dataset):
             if self.resample is not None:
                 epochs = epochs.resample(self.resample)
             # rescale to work with uV
-            X.append(dataset.unit_factor * epochs.get_data())
+            if return_epochs:
+                X.append(epochs)
+            else:
+                X.append(dataset.unit_factor * epochs.get_data())
 
         inv_events = {k: v for v, k in event_id.items()}
         labels = np.array([inv_events[e] for e in epochs.events[:, -1]])
diff --git a/moabb/run.py b/moabb/run.py
index 36b9949cc..b2c0f6a43 100755
--- a/moabb/run.py
+++ b/moabb/run.py
@@ -112,7 +112,7 @@ def parse_pipelines_from_directory(d):
             content = _file.read()
 
             # load config
-            config_dict = yaml.load(content)
+            config_dict = yaml.load(content, Loader=yaml.FullLoader)
             ppl = create_pipeline_from_config(config_dict['pipeline'])
             pipeline_configs.append({'paradigms': config_dict['paradigms'],
                                      'pipeline': ppl,
@@ -187,7 +187,7 @@ def generate_paradigms(pipeline_configs, context={}):
     context_params = {}
     if options.context is not None:
         with open(options.context, 'r') as cfile:
-            context_params = yaml.load(cfile.read())
+            context_params = yaml.load(cfile.read(), Loader=yaml.FullLoader)
 
     paradigms = generate_paradigms(pipeline_configs, context_params)
 
diff --git a/moabb/tests/download.py b/moabb/tests/download.py
index eab7956ab..1f7703c83 100644
--- a/moabb/tests/download.py
+++ b/moabb/tests/download.py
@@ -2,19 +2,19 @@
 Tests to ensure that datasets download correctly
 '''
 # from moabb.datasets.gigadb import Cho2017
-from moabb.datasets.alex_mi import AlexMI
-from moabb.datasets.physionet_mi import PhysionetMI
-from moabb.datasets.bnci import (BNCI2014001, BNCI2014002, BNCI2014004,
-                                 BNCI2014008, BNCI2014009, BNCI2015001,
-                                 BNCI2015003, BNCI2015004)
-from moabb.datasets.bbci_eeg_fnirs import Shin2017A, Shin2017B
-from moabb.datasets.upper_limb import Ofner2017
-from moabb.datasets.mpi_mi import MunichMI
-from moabb.datasets.schirrmeister2017 import Schirrmeister2017
-from moabb.datasets.Weibo2014 import Weibo2014
-from moabb.datasets.Zhou2016 import Zhou2016
-from moabb.datasets.ssvep_exo import SSVEPExo
-from moabb.datasets.braininvaders import bi2013a
+# from moabb.datasets.alex_mi import AlexMI
+# from moabb.datasets.physionet_mi import PhysionetMI
+# from moabb.datasets.bnci import (BNCI2014001, BNCI2014002, BNCI2014004,
+#                                  BNCI2014008, BNCI2014009, BNCI2015001,
+#                                  BNCI2015003, BNCI2015004)
+# from moabb.datasets.bbci_eeg_fnirs import Shin2017A, Shin2017B
+# from moabb.datasets.upper_limb import Ofner2017
+# from moabb.datasets.mpi_mi import MunichMI
+# from moabb.datasets.schirrmeister2017 import Schirrmeister2017
+# from moabb.datasets.Weibo2014 import Weibo2014
+# from moabb.datasets.Zhou2016 import Zhou2016
+# from moabb.datasets.ssvep_exo import SSVEPExo
+# from moabb.datasets.braininvaders import bi2013a
 import unittest
 import mne
 
@@ -60,60 +60,60 @@ def _get_events(raw):
     # def test_cho2017(self):
     #     self.run_dataset(Cho2017)
 
-    def test_bnci_1401(self):
-        self.run_dataset(BNCI2014001)
+    # def test_bnci_1401(self):
+    #     self.run_dataset(BNCI2014001)
 
-    def test_bnci_1402(self):
-        self.run_dataset(BNCI2014002)
+    # def test_bnci_1402(self):
+    #     self.run_dataset(BNCI2014002)
 
-    def test_bnci_1404(self):
-        self.run_dataset(BNCI2014004)
+    # def test_bnci_1404(self):
+    #     self.run_dataset(BNCI2014004)
 
-    def test_bnci_1408(self):
-        self.run_dataset(BNCI2014008)
+    # def test_bnci_1408(self):
+    #     self.run_dataset(BNCI2014008)
 
-    def test_bnci_1409(self):
-        self.run_dataset(BNCI2014009)
+    # def test_bnci_1409(self):
+    #     self.run_dataset(BNCI2014009)
 
-    def test_bnci_1501(self):
-        self.run_dataset(BNCI2015001)
+    # def test_bnci_1501(self):
+    #     self.run_dataset(BNCI2015001)
 
-    def test_bnci_1503(self):
-        self.run_dataset(BNCI2015003)
+    # def test_bnci_1503(self):
+    #     self.run_dataset(BNCI2015003)
 
-    def test_bnci_1504(self):
-        self.run_dataset(BNCI2015004)
+    # def test_bnci_1504(self):
+    #     self.run_dataset(BNCI2015004)
 
-    def test_alexmi(self):
-        self.run_dataset(AlexMI)
+    # def test_alexmi(self):
+    #     self.run_dataset(AlexMI)
 
-    def test_physionet(self):
-        self.run_dataset(PhysionetMI)
+    # def test_physionet(self):
+    #     self.run_dataset(PhysionetMI)
 
-    def test_eegfnirs(self):
-        self.run_dataset(Shin2017A)
-        self.run_dataset(Shin2017B)
+    # def test_eegfnirs(self):
+    #     self.run_dataset(Shin2017A)
+    #     self.run_dataset(Shin2017B)
 
-    def test_upper_limb(self):
-        self.run_dataset(Ofner2017)
+    # def test_upper_limb(self):
+    #     self.run_dataset(Ofner2017)
 
-    def test_mpi_mi(self):
-        self.run_dataset(MunichMI)
+    # def test_mpi_mi(self):
+    #     self.run_dataset(MunichMI)
 
-    def test_schirrmeister2017(self):
-        self.run_dataset(Schirrmeister2017, subj=(0, 1))
+    # def test_schirrmeister2017(self):
+    #     self.run_dataset(Schirrmeister2017, subj=(0, 1))
 
-    def test_Weibo2014(self):
-        self.run_dataset(Weibo2014)
+    # def test_Weibo2014(self):
+    #     self.run_dataset(Weibo2014)
 
-    def test_Zhou2016(self):
-        self.run_dataset(Zhou2016)
+    # def test_Zhou2016(self):
+    #     self.run_dataset(Zhou2016)
 
-    def test_ssvep_exo(self):
-        self.run_dataset(SSVEPExo)
+    # def test_ssvep_exo(self):
+    #     self.run_dataset(SSVEPExo)
 
-    def test_bi2013a(self):
-        self.run_dataset(bi2013a)
+    # def test_bi2013a(self):
+    #     self.run_dataset(bi2013a)
 
 
 if __name__ == '__main__':
diff --git a/moabb/tests/paradigms.py b/moabb/tests/paradigms.py
index a7ccb27c9..480b13fc0 100644
--- a/moabb/tests/paradigms.py
+++ b/moabb/tests/paradigms.py
@@ -196,9 +196,9 @@ def test_P300_paradigm(self):
         dataset = FakeDataset(event_list=['Target', 'NonTarget'],
                               paradigm='p300')
         X, labels, metadata = paradigm.get_data(dataset, subjects=[1])
-        self.assertEquals(len(np.unique(labels)), 2)
-        self.assertEquals(list(np.unique(labels)),
-                          sorted(['Target', 'NonTarget']))
+        self.assertEqual(len(np.unique(labels)), 2)
+        self.assertEqual(list(np.unique(labels)),
+                         sorted(['Target', 'NonTarget']))
 
     def test_BaseImagery_noevent(self):
         # Assert error if events from paradigm and dataset dont overlap
diff --git a/tutorials/select_electrodes_resample.py b/tutorials/select_electrodes_resample.py
new file mode 100644
index 000000000..3ec8fa4bb
--- /dev/null
+++ b/tutorials/select_electrodes_resample.py
@@ -0,0 +1,88 @@
+"""
+================================
+Select electrodes and resampling
+================================
+
+Within paradigm, it is possible to restrict analysis only to a subset of
+electrodes and to resample to a specific sampling rate. There is also a
+utility function to select common electrodes shared between datasets.
+This tutorial demonstrates how to use this functionality.
+"""
+# Authors: Sylvain Chevallier <sylvain.chevallier@uvsq.fr>
+#
+# License: BSD (3-clause)
+from moabb.datasets import BNCI2014001, Zhou2016
+from moabb.paradigms import LeftRightImagery
+from moabb.evaluations import WithinSessionEvaluation
+from moabb.datasets.utils import find_intersecting_channels
+
+from sklearn.pipeline import make_pipeline
+from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
+from sklearn.linear_model import LogisticRegression as LR
+
+from mne.decoding import CSP
+from pyriemann.estimation import Covariances
+from pyriemann.tangentspace import TangentSpace
+
+import matplotlib.pyplot as plt
+import moabb.analysis.plotting as moabb_plt
+
+##############################################################################
+# Datasets
+# --------
+#
+# Select datasets for motor imagery
+
+datasets = [Zhou2016(), BNCI2014001()]
+
+##############################################################################
+# Paradigm
+# --------
+#
+# Restrict further analysis to specified channels, here C3, C4, and Cz.
+# Also, use a specific resampling. In this example, all datasets are
+# set to 200 Hz.
+
+paradigm = LeftRightImagery(channels=['C3', 'C4', 'Cz'], resample=200.)
+
+##############################################################################
+# Evaluation
+# ----------
+#
+# The evaluation is conducted on with CSP+LDA, only on the 3 electrodes, with
+# a sampling rate of 200 Hz.
+
+evaluation = WithinSessionEvaluation(paradigm=paradigm,
+                                     datasets=datasets)
+csp_lda = make_pipeline(CSP(n_components=2), LDA())
+ts_lr = make_pipeline(Covariances(estimator='oas'),
+                      TangentSpace(metric='riemann'),
+                      LR(C=1.0))
+results = evaluation.process({'csp+lda': csp_lda, 'ts+lr': ts_lr})
+print(results.head())
+
+##############################################################################
+# Electrode selection
+# -------------------
+#
+# It is possible to select the electrodes that are shared by all datasets
+# using the `find_intersecting_channels` function. Datasets that have 0
+# overlap with others are discarded. It returns the set of common channels,
+# as well as the list of datasets with valid channels.
+
+electrodes, datasets = find_intersecting_channels(datasets)
+evaluation = WithinSessionEvaluation(paradigm=paradigm,
+                                     datasets=datasets,
+                                     overwrite=True)
+results = evaluation.process({'csp+lda': csp_lda, 'ts+lr': ts_lr})
+print(results.head())
+
+##############################################################################
+# Plot results
+# ------------
+#
+# Compare the obtained results with the two pipelines, CSP+LDA and logistic
+# regression computed in the tangent space of the covariance matrices.
+
+fig = moabb_plt.paired_plot(results, 'csp+lda', 'ts+lr')
+plt.show()