diff --git a/.circleci/config.yml b/.circleci/config.yml index 869cd781b68..5c4323bdb3e 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -28,6 +28,7 @@ jobs: # Spin up Xvfb - run: /sbin/start-stop-daemon --start --quiet --pidfile /tmp/custom_xvfb_99.pid --make-pidfile --background --exec /usr/bin/Xvfb -- :99 -screen 0 1400x900x24 -ac +extension GLX +render -noreset; - run: echo "export DISPLAY=:99" >> $BASH_ENV; + - run: echo "export OPENBLAS_NUM_THREADS=1" >> $BASH_ENV; # Fix libgcc_s.so.1 pthread_cancel bug: # https://github.com/ContinuumIO/anaconda-issues/issues/9190#issuecomment-386508136 # https://github.com/golemfactory/golem/issues/1019 @@ -45,7 +46,7 @@ jobs: conda env create --quiet -f environment.yml; source activate mne; conda install sphinx; - pip install sphinx_fontawesome sphinx_bootstrap_theme "https://api.github.com/repos/sphinx-gallery/sphinx-gallery/zipball/master" + pip install sphinx_fontawesome sphinx_bootstrap_theme "https://api.github.com/repos/sphinx-gallery/sphinx-gallery/zipball/master" memory_profiler pip uninstall --yes mne; echo "source activate mne" >> $BASH_ENV; - save_cache: diff --git a/.mailmap b/.mailmap index 25ebcd9bc9c..0bce911ac65 100644 --- a/.mailmap +++ b/.mailmap @@ -112,3 +112,6 @@ Erik Hornberger er06645810 Maksymenko Kostiantyn Nathalie Gayraud Nathalie Dominik Krzemiński dokato +Ezequiel Mikulan <39155887+ezemikulan@users.noreply.github.com> ezemikulan <39155887+ezemikulan@users.noreply.github.com> +Hubert Banville hubertjb +Sara Sommariva sarasommariva diff --git a/doc/Makefile b/doc/Makefile index fc3cd566658..68151a7f150 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -2,7 +2,7 @@ # # You can set these variables from the command line. -SPHINXOPTS = +SPHINXOPTS = -n -W -T --keep-going SPHINXBUILD = sphinx-build PAPER = @@ -33,13 +33,7 @@ help: @echo " doctest to run all doctests embedded in the documentation (if enabled)" clean: - -rm -rf _build/* - -rm -rf auto_examples - -rm -rf auto_tutorials - -rm -rf generated - -rm -rf *.stc - -rm -rf *.fif - -rm -rf *.nii.gz + -rm -rf _build/* auto_examples auto_tutorials generated *.stc *.fif *.nii.gz html_stable: $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) _build/html_stable diff --git a/doc/carousel.inc b/doc/carousel.inc index 4aaa8447575..cab12669030 100644 --- a/doc/carousel.inc +++ b/doc/carousel.inc @@ -36,7 +36,7 @@
-
Decoding
+
Decoding

Machine learning

diff --git a/doc/conf.py b/doc/conf.py index e6393d3b332..fdb2bf58083 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -12,18 +12,24 @@ # All configuration values have a default; values that are commented out # serve to show the default. -import inspect +from datetime import date +from distutils.version import LooseVersion import os -from os.path import relpath, dirname import sys -from datetime import date -import sphinx_gallery # noqa +import warnings + +import sphinx_gallery +from sphinx_gallery.sorting import FileNameSortKey import sphinx_bootstrap_theme from numpydoc import numpydoc, docscrape # noqa import sphinx_fontawesome import mne from mne.utils import linkcode_resolve # noqa, analysis:ignore +if LooseVersion(sphinx_gallery.__version__) < LooseVersion('0.2'): + raise ImportError('Must have at least version 0.2 of sphinx-gallery, got ' + '%s' % (sphinx_gallery.__version__,)) + # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. @@ -37,7 +43,7 @@ # -- General configuration ------------------------------------------------ # If your documentation needs a minimal Sphinx version, state it here. -#needs_sphinx = '1.0' +needs_sphinx = '1.8' # XXX This hack defines what extra methods numpydoc will document docscrape.ClassDoc.extra_public_methods = mne.utils._doc_special_members @@ -61,7 +67,7 @@ ] autosummary_generate = True -autodoc_default_flags = ['inherited-members'] +autodoc_default_options = {'inherited-members': None} # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -285,7 +291,7 @@ 'numpy': ('https://www.numpy.org/devdocs', None), 'scipy': ('https://scipy.github.io/devdocs', None), 'matplotlib': ('https://matplotlib.org', None), - 'sklearn': ('http://scikit-learn.org/stable', None), + 'sklearn': ('https://scikit-learn.org/stable', None), 'mayavi': ('http://docs.enthought.com/mayavi/mayavi', None), 'nibabel': ('http://nipy.org/nibabel', None), 'nilearn': ('http://nilearn.github.io', None), @@ -300,21 +306,63 @@ try: mlab = mne.utils._import_mlab() - find_mayavi_figures = True # Do not pop up any mayavi windows while running the # examples. These are very annoying since they steal the focus. mlab.options.offscreen = True + scrapers = ('matplotlib', 'mayavi') except Exception: - find_mayavi_figures = False - + scrapers = ('matplotlib',) +else: + # Let's do the same thing we do in tests: reraise traits exceptions + from traits.api import push_exception_handler + push_exception_handler(reraise_exceptions=True) + + +def reset_warnings(gallery_conf, fname): + """Ensure we are future compatible and ignore silly warnings.""" + # In principle, our examples should produce no warnings. + # Here we cause warnings to become errors, with a few exceptions. + # This list should be considered alongside + # setup.cfg -> [tool:pytest] -> filterwarnings + + # remove tweaks from other module imports or example runs + warnings.resetwarnings() + # restrict + warnings.filterwarnings('error') + # allow these, but show them + warnings.filterwarnings('always', '.*cannot make axes width small.*') + warnings.filterwarnings('always', '.*Axes that are not compatible.*') + warnings.filterwarnings('always', '.*FastICA did not converge.*') + warnings.filterwarnings( # xhemi morph (should probably update sample) + 'always', '.*does not exist, creating it and saving it.*') + warnings.filterwarnings('default', module='sphinx') # internal warnings + warnings.filterwarnings( + 'always', '.*converting a masked element to nan.*') # matplotlib? + # allow these warnings, but don't show them + warnings.filterwarnings('ignore', '.*is currently using agg.*') + warnings.filterwarnings( # SciPy-related warning (maybe 1.2.0 will fix it) + 'ignore', '.*the matrix subclass is not the recommended.*') + warnings.filterwarnings( # some joblib warning + 'ignore', '.*semaphore_tracker: process died unexpectedly.*') + warnings.filterwarnings( # needed until SciPy 1.2.0 is released + 'ignore', '.*will be interpreted as an array index.*', module='scipy') + for key in ('HasTraits', r'numpy\.testing', 'importlib', r'np\.loads', + r"it will be an error for 'np\.bool_'", # ndimage + "'U' mode is deprecated", # sphinx io + ): + warnings.filterwarnings( # deal with other modules having bad imports + 'ignore', message=".*%s.*" % key, category=DeprecationWarning) + # allow this ImportWarning, but don't show it + warnings.filterwarnings( + 'ignore', message="can't resolve package from", category=ImportWarning) + + +reset_warnings(None, None) sphinx_gallery_conf = { 'doc_module': ('mne',), - 'reference_url': { - 'mne': None, - }, + 'reference_url': dict(mne=None), 'examples_dirs': examples_dirs, 'gallery_dirs': gallery_dirs, - 'find_mayavi_figures': find_mayavi_figures, 'default_thumb_file': os.path.join('_static', 'mne_helmet.png'), 'backreferences_dir': 'generated', 'plot_gallery': 'True', # Avoid annoying Unicode/bool default warning @@ -322,6 +370,11 @@ 'thumbnail_size': (160, 112), 'min_reported_time': 1., 'abort_on_example_error': False, + 'reset_modules': ('matplotlib', reset_warnings), # called w/each script + 'image_scrapers': scrapers, + 'show_memory': True, + 'line_numbers': True, + 'within_subsection_order': FileNameSortKey, } numpydoc_class_members_toctree = False diff --git a/doc/contributing.rst b/doc/contributing.rst index 516edf20bd7..1b023ad35a6 100644 --- a/doc/contributing.rst +++ b/doc/contributing.rst @@ -120,9 +120,14 @@ see `Sphinx documentation`_ to learn more about editing them. Our code follows the `NumPy docstring standard`_. Documentation is automatically built remotely during pull requests. If -you want to also test documentation locally, you will need to install -``sphinx sphinx-gallery sphinx_bootstrap_theme sphinx_fontawesome``, and then -within the ``mne/doc`` directory do: +you want to also test documentation locally, you will need to install (using +``conda`` or ``pip`` as preferred): + +.. code-block:: console + + $ pip install sphinx sphinx-gallery sphinx_bootstrap_theme sphinx_fontawesome memory_profiler + +and then within the ``mne/doc`` directory do: .. code-block:: console diff --git a/doc/documentation.rst b/doc/documentation.rst index 8011d771194..0ca5b0b3463 100644 --- a/doc/documentation.rst +++ b/doc/documentation.rst @@ -553,7 +553,6 @@ There are also **examples**, which contain a short use-case to highlight MNE-fun .. toctree:: :maxdepth: 1 - manual/decoding.rst auto_tutorials/plot_sensors_decoding.rst .. raw:: html @@ -566,7 +565,6 @@ There are also **examples**, which contain a short use-case to highlight MNE-fun auto_examples/decoding/decoding_rsa.rst auto_examples/decoding/plot_decoding_csp_eeg.rst - auto_examples/decoding/plot_decoding_csp_space.rst auto_examples/decoding/plot_decoding_csp_timefreq.rst auto_examples/decoding/plot_decoding_spatio_temporal_source.rst auto_examples/decoding/plot_decoding_spoc_CMC.rst diff --git a/doc/manual/channel_interpolation.rst b/doc/manual/channel_interpolation.rst index a7a64316592..3730bbda0ed 100644 --- a/doc/manual/channel_interpolation.rst +++ b/doc/manual/channel_interpolation.rst @@ -60,11 +60,7 @@ To interpolate bad channels, one can simply do: >>> evoked.interpolate_bads(reset_bads=False) # doctest: +SKIP -and the bad channel will be fixed - -.. image:: ../../_images/sphx_glr_plot_interpolate_bad_channels_002.png - :align: center - :height: 300 px +and the bad channel will be fixed. .. topic:: Examples: diff --git a/doc/manual/decoding.rst b/doc/manual/decoding.rst deleted file mode 100644 index 2483b45a49e..00000000000 --- a/doc/manual/decoding.rst +++ /dev/null @@ -1,200 +0,0 @@ -.. include:: ../git_links.inc - -.. contents:: Contents - :local: - :depth: 3 - -.. _decoding: - -Decoding in MNE -############### - -For maximal compatibility with the Scikit-learn package, we follow the same API. Each estimator implements a ``fit``, a ``transform`` and a ``fit_transform`` method. In some cases, they also implement an ``inverse_transform`` method. For more details, visit the Scikit-learn page. - -For ease of comprehension, we will denote instantiations of the class using the same name as the class but in small caps instead of camel cases. - -Basic Estimators -================ - -Scaler -^^^^^^ -The :class:`mne.decoding.Scaler` will standardize the data based on channel scales. In the simplest modes ``scalings=None`` or ``scalings=dict(...)``, each data channel type (e.g., mag, grad, eeg) is treated separately and scaled by a constant. This is the approach used by e.g., :func:`mne.compute_covariance` to standardize channel scales. - -If ``scalings='mean'`` or ``scalings='median'``, each channel is scaled using empirical measures. Each channel is scaled independently by the mean and standand deviation, or median and interquartile range, respectively, across all epochs and time points during :class:`mne.decoding.Scaler.fit` (during training). The :meth:`mne.decoding.Scaler.transform` method is called to transform data (training or test set) by scaling all time points and epochs on a channel-by-channel basis. To perform both the ``fit`` and ``transform`` operations in a single call, the :meth:`mne.decoding.Scaler.fit_transform` method may be used. To invert the transform, :meth:`mne.decoding.Scaler.inverse_transform` can be used. For ``scalings='median'``, scikit-learn_ version 0.17+ is required. - -.. note:: This is different from directly applying :class:`sklearn.preprocessing.StandardScaler` or :class:`sklearn.preprocessing.RobustScaler` offered by scikit-learn_. The ``StandardScaler`` and ``RobustScaler`` scale each *classification feature*, e.g. each time point for each channel, with mean and standard deviation computed across epochs, whereas ``Scaler`` scales each *channel* using mean and standard deviation computed across all of its time points and epochs. - -Vectorizer -^^^^^^^^^^ -Scikit-learn API provides functionality to chain transformers and estimators by using :class:`sklearn.pipeline.Pipeline`. We can construct decoding pipelines and perform cross-validation and grid-search. However scikit-learn transformers and estimators generally expect 2D data (n_samples * n_features), whereas MNE transformers typically output data with a higher dimensionality (e.g. n_samples * n_channels * n_frequencies * n_times). A Vectorizer therefore needs to be applied between the MNE and the scikit-learn steps: e.g: make_pipeline(Xdawn(), Vectorizer(), LogisticRegression()) - -PSDEstimator -^^^^^^^^^^^^ -This estimator computes the power spectral density (PSD) using the multitaper method. It takes a 3D array as input, it into 2D and computes the PSD. - -FilterEstimator -^^^^^^^^^^^^^^^ -This estimator filters the 3D epochs data. - -.. warning:: This is meant for use in conjunction with ``RtEpochs``. It is not recommended in a normal processing pipeline as it may result in edge artifacts. - -Spatial filters -=============== - -Just like temporal filters, spatial filters provide weights to modify the data along the sensor dimension. They are popular in the BCI community because of their simplicity and ability to distinguish spatially-separated neural activity. - -Common Spatial Pattern -^^^^^^^^^^^^^^^^^^^^^^ - -This is a technique to analyze multichannel data based on recordings from two classes. Let :math:`X \in R^{C\times T}` be a segment of data with :math:`C` channels and :math:`T` time points. The data at a single time point is denoted by :math:`x(t)` such that :math:`X=[x(t), x(t+1), ..., x(t+T-1)]`. Common Spatial Pattern (CSP) finds a decomposition that projects the signal in the original sensor space to CSP space using the following transformation: - -.. math:: x_{CSP}(t) = W^{T}x(t) - :label: csp - -where each column of :math:`W \in R^{C\times C}` is a spatial filter and each row of :math:`x_{CSP}` is a CSP component. The matrix :math:`W` is also called the de-mixing matrix in other contexts. Let :math:`\Sigma^{+} \in R^{C\times C}` and :math:`\Sigma^{-} \in R^{C\times C}` be the estimates of the covariance matrices of the two conditions. -CSP analysis is given by the simultaneous diagonalization of the two covariance matrices - -.. math:: W^{T}\Sigma^{+}W = \lambda^{+} - :label: diagonalize_p -.. math:: W^{T}\Sigma^{-}W = \lambda^{-} - :label: diagonalize_n - -where :math:`\lambda^{C}` is a diagonal matrix whose entries are the eigenvalues of the following generalized eigenvalue problem - -.. math:: \Sigma^{+}w = \lambda \Sigma^{-}w - :label: eigen_problem - -Large entries in the diagonal matrix corresponds to a spatial filter which gives high variance in one class but low variance in the other. Thus, the filter facilitates discrimination between the two classes. - -.. topic:: Examples: - - * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_csp_eeg.py` - * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_csp_space.py` - -.. topic:: Spotlight: - - The winning entry of the Grasp-and-lift EEG competition in Kaggle uses the CSP implementation in MNE. It was featured as a `script of the week`_. - - -Source Power Comodulation (SPoC) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Source Power Comodulation (SPoC) [1]_ allows to identify the composition of orthogonal spatial filters that maximally correlate with a continuous target. - -SPoC can be seen as an extension of the CSP where the target is driven by a continuous variable rather than a discrete variable. Typical applications include extraction of motor patterns using EMG power or audio patterns using sound envelope. - -.. topic:: Examples - - * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_spoc_CMC.py` - -xDAWN -^^^^^ -Xdawn is a spatial filtering method designed to improve the signal to signal + noise ratio (SSNR) of the ERP responses. Xdawn was originally designed for P300 evoked potential by enhancing the target response with respect to the non-target response. The implementation in MNE-Python is a generalization to any type of ERP. - -.. topic:: Examples: - - * :ref:`sphx_glr_auto_examples_preprocessing_plot_xdawn_denoising.py` - * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_xdawn_eeg.py` - -Effect-matched spatial filtering -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The result is a spatial filter at each time point and a corresponding time course. Intuitively, the result gives the similarity between the filter at each time point and the data vector (sensors) at that time point. - -.. topic:: Examples - - * :ref:`sphx_glr_auto_examples_decoding_plot_ems_filtering.py` - -Patterns vs. filters -^^^^^^^^^^^^^^^^^^^^ - -When interpreting the components of the CSP, it is often more intuitive to think about how :math:`x(t)` is composed of the different CSP components :math:`x_{CSP}(t)`. In other words, we can rewrite Equation :eq:`csp` as follows: - -.. math:: x(t) = (W^{-1})^{T}x_{CSP}(t) - :label: patterns - -The columns of the matrix :math:`(W^{-1})^T` are called spatial patterns. This is also called the mixing matrix. The example :ref:`sphx_glr_auto_examples_decoding_plot_linear_model_patterns.py` demonstrates the difference between patterns and filters. - -Plotting a pattern is as simple as doing:: - - >>> info = epochs.info # doctest: +SKIP - >>> model.plot_patterns(info) # model is an instantiation of an estimator described in this section # doctest: +SKIP - -.. image:: ../../_images/sphx_glr_plot_linear_model_patterns_001.png - :align: center - :height: 100 px - -To plot the corresponding filter, you can do:: - - >>> model.plot_filters(info) # doctest: +SKIP - -.. image:: ../../_images/sphx_glr_plot_linear_model_patterns_002.png - :align: center - :height: 100 px - -Sensor-space decoding -===================== - -Decoding over time -^^^^^^^^^^^^^^^^^^ - -This strategy consists in fitting a multivariate predictive model on each -time instant and evaluating its performance at the same instant on new -epochs. The :class:`mne.decoding.SlidingEstimator` will take as input a -pair of features :math:`X` and targets :math:`y`, where :math:`X` has -more than 2 dimensions. For decoding over time the data :math:`X` -is the epochs data of shape n_epochs x n_channels x n_times. As the -last dimension of :math:`X` is the time an estimator will be fit -on every time instant. - -This approach is analogous to SlidingEstimator-based approaches in fMRI, -where here we are interested in when one can discriminate experimental -conditions and therefore figure out when the effect of interest happens. - -When working with linear models as estimators, this approach boils -down to estimating a discriminative spatial filter for each time instant. - -.. image:: ../../_images/sphx_glr_plot_sensors_decoding_001.png - :align: center - :width: 400px - -To generate this plot see our tutorial :ref:`sphx_glr_auto_tutorials_plot_sensors_decoding.py`. - -Temporal Generalization -^^^^^^^^^^^^^^^^^^^^^^^ - -Temporal Generalization is an extension of the decoding over time approach. -It consists in evaluating whether the model estimated at a particular -time instant accurately predicts any other time instant. It is analogous to -transferring a trained model to a distinct learning problem, where the problems -correspond to decoding the patterns of brain activity recorded at distinct time -instants. - -The object to for Temporal Generalization is -:class:`mne.decoding.GeneralizingEstimator`. It expects as input :math:`X` and -:math:`y` (similarly to :class:`mne.decoding.SlidingEstimator`) but, when generate -predictions from each model for all time instants. The class -:class:`mne.decoding.GeneralizingEstimator` is generic and will treat the last -dimension as the one to be used for generalization testing. For convenience, -here, we refer to it different tasks. If :math:`X` corresponds to epochs data -then the last dimension is time. - -.. image:: ../../_images/sphx_glr_plot_sensors_decoding_004.png - :align: center - :width: 400px - -To generate this plot see our tutorial :ref:`sphx_glr_auto_tutorials_plot_sensors_decoding.py`. - -Source-space decoding -===================== - -Source space decoding is also possible, but because the number of features can be much larger than in the sensor space, univariate feature selection using ANOVA f-test (or some other metric) can be done to reduce the feature dimension. Interpreting decoding results might be easier in source space as compared to sensor space. - -.. topic:: Examples: - - * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_spatio_temporal_source.py` - -.. _script of the week: http://blog.kaggle.com/2015/08/12/july-2015-scripts-of-the-week/ - -References -========== - -.. [1] Dahne, S., Meinecke, F. C., Haufe, S., Hohne, J., Tangermann, M., Muller, K. R., & Nikulin, V. V. (2014). SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters. NeuroImage, 86, 111-122. diff --git a/doc/python_reference.rst b/doc/python_reference.rst index 3f82bf141eb..967691389f3 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -56,7 +56,6 @@ Reading raw data anonymize_info find_edf_events - read_events_eeglab read_raw_artemis123 read_raw_bti read_raw_cnt diff --git a/doc/whats_new.rst b/doc/whats_new.rst index ba34c20e343..ee32eb9d7bd 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -111,7 +111,7 @@ Changelog - Add parameter ``rank=None`` to :func:`mne.beamformer.make_dics` by `Marijn van Vliet`_ -- Add parameter ``rank='full'`` to :func:`mne.beamformer.make_lcmv``, which can be set to ``None`` to auto-compute the rank of the covariance matrix before regularization by `Marijn van Vliet`_ +- Add parameter ``rank='full'`` to :func:`mne.beamformer.make_lcmv`, which can be set to ``None`` to auto-compute the rank of the covariance matrix before regularization by `Marijn van Vliet`_ Bug @@ -147,6 +147,8 @@ Bug - Fix bug with channel names in ``mgh70`` montage in :func:`mne.channels.read_montage` by `Eric Larson`_ +- Fix duplication of ``info['hpi_meas']`` and ``info['hpi_results']`` by `Sara Sommariva`_ + - Fix bug in :func:`mne.io.read_raw_edf` when reading large files on Windows by `Marcin Koculak`_ - Fix check in :func:`mne.viz.plot_sensors` for invalid channel locations by `Eric Larson`_ @@ -215,6 +217,8 @@ Bug - Fix normalization error in :func:`mne.beamformer.make_lcmv` when ``pick_ori='normal', weight_norm='unit_noise_gain'`` by `Marijn van Vliet`_ +- Update MNE-C installation instructions by `buildqa`_ + - Fix computation of max-power orientation in :func:`mne.beamformer.make_dics` when ``pick_ori='max-power', weight_norm='unit_noise_gain'`` by `Marijn van Vliet`_ API @@ -748,7 +752,7 @@ API - Allow passing a list of channel names as ``show_names`` in function :func:`mne.viz.plot_sensors` and methods :meth:`mne.Evoked.plot_sensors`, :meth:`mne.Epochs.plot_sensors` and :meth:`mne.io.Raw.plot_sensors` to show only a subset of channel names by `Jaakko Leppakangas`_ -- Make function `mne.io.eeglab.read_events_eeglab` public to allow loading overlapping events from EEGLAB files, by `Jona Sassenhagen`_. +- Make function ``mne.io.eeglab.read_events_eeglab`` public to allow loading overlapping events from EEGLAB files, by `Jona Sassenhagen`_. - :func:`mne.find_events` ``mask_type`` parameter will change from ``'not_and'`` to ``'and'`` in 0.16. @@ -1723,6 +1727,8 @@ Changelog - Adapt ``corrmap`` function (Viola et al. 2009) to semi-automatically detect similar ICs across data sets by `Jona Sassenhagen`_ and `Denis Engemann`_ and `Eric Larson`_ +- Clarify docstring for :class:`mne.preprocessing.ICA` by `jeythekey`_ + - New ``mne flash_bem`` command to compute BEM surfaces from Flash MRI images by `Lorenzo Desantis`_, `Alex Gramfort`_ and `Eric Larson`_. See :func:`mne.bem.make_flash_bem`. - New gfp parameter in :func:`mne.Evoked.plot` method to display Global Field Power (GFP) by `Eric Larson`_. @@ -2996,3 +3002,11 @@ of commits): .. _Marcin Koculak: https://github.com/mkoculak .. _David Sabbagh: https://github.com/DavidSabbagh + +.. _Hubert Banville: https://github.com/hubertjb + +.. _buildqa: https://github.com/buildqa + +.. _jeythekey: https://github.com/jeythekey + +.. _Sara Sommariva: http://www.dima.unige.it/~sommariva/ diff --git a/examples/connectivity/plot_mixed_source_space_connectivity.py b/examples/connectivity/plot_mixed_source_space_connectivity.py index db3c12fc6f8..2171eafa9ed 100644 --- a/examples/connectivity/plot_mixed_source_space_connectivity.py +++ b/examples/connectivity/plot_mixed_source_space_connectivity.py @@ -53,25 +53,24 @@ 'Right-Thalamus-Proper', 'Right-Cerebellum-Cortex'] -# Setup a surface-based source space +# Setup a surface-based source space, oct5 is not very dense (just used +# to speed up this example) src = setup_source_space(subject, subjects_dir=subjects_dir, - spacing='oct6', add_dist=False) + spacing='oct5', add_dist=False) # Setup a volume source space -# set pos=7.0 for speed issue -vol_src = setup_volume_source_space(subject, mri=fname_aseg, - pos=7.0, - bem=fname_model, - volume_label=labels_vol, - subjects_dir=subjects_dir) +# set pos=10.0 for speed, not very accurate +vol_src = setup_volume_source_space( + subject, mri=fname_aseg, pos=10.0, bem=fname_model, + add_interpolator=False, # just for speed, usually use True + volume_label=labels_vol, subjects_dir=subjects_dir) # Generate the mixed source space src += vol_src # compute the fwd matrix fwd = make_forward_solution(fname_raw, fname_trans, src, fname_bem, mindist=5.0, # ignore sources<=5mm from innerskull - meg=True, eeg=False, - n_jobs=1) + meg=True, eeg=False, n_jobs=1) # Load data raw = read_raw_fif(fname_raw, preload=True) diff --git a/examples/connectivity/plot_sensor_connectivity.py b/examples/connectivity/plot_sensor_connectivity.py index e48a307448e..1e81ee71999 100644 --- a/examples/connectivity/plot_sensor_connectivity.py +++ b/examples/connectivity/plot_sensor_connectivity.py @@ -102,7 +102,7 @@ points.module_manager.scalar_lut_manager.reverse_lut = True -mlab.scalarbar(title='Phase Lag Index (PLI)', nb_labels=4) +mlab.scalarbar(points, title='Phase Lag Index (PLI)', nb_labels=4) # Add the sensor names for the connections shown nodes_shown = list(set([n[0] for n in con_nodes] + diff --git a/examples/datasets/plot_megsim_data.py b/examples/datasets/plot_megsim_data.py index 12fd1956667..26e964d6cb4 100644 --- a/examples/datasets/plot_megsim_data.py +++ b/examples/datasets/plot_megsim_data.py @@ -48,5 +48,6 @@ evoked.plot(time_unit='s') # Compare to the simulated data (use verbose='error' b/c of naming) -evoked_sim = read_evokeds(evoked_fnames[0], condition=0, verbose='error') +evoked_sim = read_evokeds(evoked_fnames[0], condition=0, verbose='error', + baseline=(None, 0)) evoked_sim.plot(time_unit='s') diff --git a/examples/decoding/plot_decoding_csp_space.py b/examples/decoding/plot_decoding_csp_space.py deleted file mode 100644 index 02daf1f29d4..00000000000 --- a/examples/decoding/plot_decoding_csp_space.py +++ /dev/null @@ -1,113 +0,0 @@ -""" -==================================================================== -Decoding in sensor space data using the Common Spatial Pattern (CSP) -==================================================================== - -Decoding applied to MEG data in sensor space decomposed using CSP. -Here the classifier is applied to features extracted on CSP filtered signals. - -See http://en.wikipedia.org/wiki/Common_spatial_pattern and [1]_. - -References ----------- - -.. [1] Zoltan J. Koles. The quantitative extraction and topographic mapping - of the abnormal components in the clinical EEG. Electroencephalography - and Clinical Neurophysiology, 79(6):440--447, December 1991. -""" -# Authors: Alexandre Gramfort -# Romain Trachel -# -# License: BSD (3-clause) - -import numpy as np -import matplotlib.pyplot as plt - -import mne -from mne import io -from mne.datasets import sample - -print(__doc__) - -data_path = sample.data_path() - -############################################################################### -# Set parameters and read data -raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' -event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' -tmin, tmax = -0.2, 0.5 -event_id = dict(aud_l=1, vis_l=3) - -# Setup for reading the raw data -raw = io.read_raw_fif(raw_fname, preload=True) -raw.filter(2, None, fir_design='firwin') # replace baselining with high-pass -events = mne.read_events(event_fname) - -raw.info['bads'] = ['MEG 2443'] # set bad channels -picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=False, - exclude='bads') - -# Read epochs -epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, - picks=picks, baseline=None, preload=True) - -labels = epochs.events[:, -1] -evoked = epochs.average() - -############################################################################### -# Decoding in sensor space using a linear SVM - -from sklearn.svm import SVC # noqa -from sklearn.model_selection import ShuffleSplit # noqa -from mne.decoding import CSP # noqa - -n_components = 3 # pick some components -svc = SVC(C=1, kernel='linear') -csp = CSP(n_components=n_components, norm_trace=False) - -# Define a monte-carlo cross-validation generator (reduce variance): -cv = ShuffleSplit(n_splits=10, test_size=0.2, random_state=42) -scores = [] -epochs_data = epochs.get_data() - -for train_idx, test_idx in cv.split(labels): - y_train, y_test = labels[train_idx], labels[test_idx] - - X_train = csp.fit_transform(epochs_data[train_idx], y_train) - X_test = csp.transform(epochs_data[test_idx]) - - # fit classifier - svc.fit(X_train, y_train) - - scores.append(svc.score(X_test, y_test)) - -# Printing the results -class_balance = np.mean(labels == labels[0]) -class_balance = max(class_balance, 1. - class_balance) -print("Classification accuracy: %f / Chance level: %f" % (np.mean(scores), - class_balance)) - -# Or use much more convenient scikit-learn cross_val_score function using -# a Pipeline -from sklearn.pipeline import Pipeline # noqa -from sklearn.model_selection import cross_val_score # noqa -cv = ShuffleSplit(n_splits=10, test_size=0.2, random_state=42) -clf = Pipeline([('CSP', csp), ('SVC', svc)]) -scores = cross_val_score(clf, epochs_data, labels, cv=cv, n_jobs=1) -print(scores.mean()) # should match results above - -# And using reuglarized csp with Ledoit-Wolf estimator -csp = CSP(n_components=n_components, reg='ledoit_wolf', norm_trace=False) -clf = Pipeline([('CSP', csp), ('SVC', svc)]) -scores = cross_val_score(clf, epochs_data, labels, cv=cv, n_jobs=1) -print(scores.mean()) # should get better results than above - -# plot CSP patterns estimated on full data for visualization -csp.fit_transform(epochs_data, labels) -data = csp.patterns_ -fig, axes = plt.subplots(1, 4) -for idx in range(4): - mne.viz.plot_topomap(data[idx], evoked.info, axes=axes[idx], show=False) -fig.suptitle('CSP patterns') -fig.tight_layout() -mne.viz.utils.plt_show() diff --git a/examples/decoding/plot_decoding_spatio_temporal_source.py b/examples/decoding/plot_decoding_spatio_temporal_source.py index f9a3dbc4f84..f3390de98cf 100644 --- a/examples/decoding/plot_decoding_spatio_temporal_source.py +++ b/examples/decoding/plot_decoding_spatio_temporal_source.py @@ -1,4 +1,6 @@ """ +.. _tut_dec_st_source: + ========================== Decoding source space data ========================== @@ -8,47 +10,47 @@ classification to a small number of potentially relevant features. The classifier then is trained to selected features of epochs in source space. """ +# sphinx_gallery_thumbnail_number = 2 + # Author: Denis A. Engemann # Alexandre Gramfort # Jean-Remi King +# Eric Larson # # License: BSD (3-clause) -import os + import numpy as np import matplotlib.pyplot as plt + from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler from sklearn.feature_selection import SelectKBest, f_classif from sklearn.linear_model import LogisticRegression + import mne -from mne import io -from mne.datasets import sample from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator from mne.decoding import (cross_val_multiscore, LinearModel, SlidingEstimator, get_coef) print(__doc__) -data_path = sample.data_path() +data_path = mne.datasets.sample.data_path() fname_fwd = data_path + 'MEG/sample/sample_audvis-meg-oct-6-fwd.fif' fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif' subjects_dir = data_path + '/subjects' -subject = os.environ['SUBJECT'] = subjects_dir + '/sample' -os.environ['SUBJECTS_DIR'] = subjects_dir ############################################################################### # Set parameters raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif' -label_names = 'Aud-rh', 'Vis-rh' fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' tmin, tmax = -0.2, 0.8 event_id = dict(aud_r=2, vis_r=4) # load contra-lateral conditions # Setup for reading the raw data -raw = io.read_raw_fif(raw_fname, preload=True) +raw = mne.io.read_raw_fif(raw_fname, preload=True) raw.filter(None, 10., fir_design='firwin') events = mne.read_events(event_fname) @@ -83,7 +85,7 @@ # prepare a series of classifier applied at each time sample clf = make_pipeline(StandardScaler(), # z-score normalization SelectKBest(f_classif, k=500), # select features for speed - LinearModel(LogisticRegression(C=1))) + LinearModel(LogisticRegression(C=1, solver='liblinear'))) time_decod = SlidingEstimator(clf, scoring='roc_auc') # Run cross-validated decoding analyses: @@ -112,4 +114,5 @@ tmin=stc.tmin, tstep=stc.tstep, subject='sample') brain = stc_feat.plot(views=['lat'], transparent=True, - initial_time=0.1, time_unit='s') + initial_time=0.1, time_unit='s', + subjects_dir=subjects_dir) diff --git a/examples/decoding/plot_decoding_spoc_CMC.py b/examples/decoding/plot_decoding_spoc_CMC.py index 564a779d066..8469c1d1761 100644 --- a/examples/decoding/plot_decoding_spoc_CMC.py +++ b/examples/decoding/plot_decoding_spoc_CMC.py @@ -62,7 +62,8 @@ y = emg_epochs.get_data().var(axis=2)[:, 0] # target is EMG power # Classification pipeline with SPoC spatial filtering and Ridge Regression -clf = make_pipeline(SPoC(n_components=2, log=True, reg='oas'), Ridge()) +clf = make_pipeline(SPoC(n_components=2, log=True, reg='oas', + rank='full'), Ridge()) # Define a two fold cross-validation cv = KFold(n_splits=2, shuffle=False) diff --git a/examples/decoding/plot_decoding_xdawn_eeg.py b/examples/decoding/plot_decoding_xdawn_eeg.py index 5d03b677247..3cb6a158639 100644 --- a/examples/decoding/plot_decoding_xdawn_eeg.py +++ b/examples/decoding/plot_decoding_xdawn_eeg.py @@ -66,7 +66,8 @@ clf = make_pipeline(Xdawn(n_components=3), Vectorizer(), MinMaxScaler(), - LogisticRegression(penalty='l1')) + LogisticRegression(penalty='l1', solver='liblinear', + multi_class='auto')) # Get the labels labels = epochs.events[:, -1] diff --git a/examples/decoding/plot_ems_filtering.py b/examples/decoding/plot_ems_filtering.py index 5c6bf5b3c96..39fce25b4fa 100644 --- a/examples/decoding/plot_ems_filtering.py +++ b/examples/decoding/plot_ems_filtering.py @@ -79,7 +79,7 @@ # to overfit and cannot be used to estimate the variance of the # prediction within a given fold. -for train, test in StratifiedKFold().split(X, y): +for train, test in StratifiedKFold(n_splits=5).split(X, y): # In the original paper, the z-scoring is applied outside the CV. # However, we recommend to apply this preprocessing inside the CV. # Note that such scaling should be done separately for each channels if the @@ -121,7 +121,7 @@ # Visualize spatial filters across time evoked = EvokedArray(filters, epochs.info, tmin=epochs.tmin) -evoked.plot_topomap(time_unit='s') +evoked.plot_topomap(time_unit='s', scalings=1) ############################################################################# # Note that a similar transformation can be applied with `compute_ems` diff --git a/examples/inverse/plot_covariance_whitening_dspm.py b/examples/inverse/plot_covariance_whitening_dspm.py index 27d51936c31..0d280d60ea2 100644 --- a/examples/inverse/plot_covariance_whitening_dspm.py +++ b/examples/inverse/plot_covariance_whitening_dspm.py @@ -11,12 +11,12 @@ This example makes use of the previous, non-optimized code path that was used before implementing the suggestions presented in [1]_. +This example does quite a bit of processing, so even on a +fast machine it can take a couple of minutes to complete. + .. warning:: Please do not copy the patterns presented here for your own analysis, this is example is purely illustrative. -.. note:: This example does quite a bit of processing, so even on a - fast machine it can take a couple of minutes to complete. - References ---------- .. [1] Engemann D. and Gramfort A. (2015) Automated model selection in @@ -53,20 +53,18 @@ raw = raw.crop(0, 150.).load_data() picks = mne.pick_types(raw.info, meg=True, exclude='bads') -raw.filter(1, 20., n_jobs=1, fir_design='firwin') +raw.filter(None, 20.) events = mne.find_events(raw, stim_channel='UPPT001') event_ids = {"faces": 1, "scrambled": 2} tmin, tmax = -0.2, 0.5 -baseline = None # no baseline as high-pass is applied +baseline = (None, 0) reject = dict(mag=3e-12) -# Make source space - +# Make forward trans = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_raw-trans.fif' -src = mne.setup_source_space('spm', spacing='oct6', subjects_dir=subjects_dir, - add_dist=False) +src = data_path + '/subjects/spm/bem/spm-oct-6-src.fif' bem = data_path + '/subjects/spm/bem/spm-5120-5120-5120-bem-sol.fif' forward = mne.make_forward_solution(raw.info, trans, src, bem) del src @@ -75,7 +73,6 @@ conditions = 'faces', 'scrambled' snr = 3.0 lambda2 = 1.0 / snr ** 2 -method = 'dSPM' clim = dict(kind='value', lims=[0, 2.5, 5]) ############################################################################### @@ -173,4 +170,3 @@ ax_dynamics.legend(loc='upper left', fontsize=10) fig.subplots_adjust(hspace=0.2, left=0.01, right=0.99, wspace=0.03) -mne.viz.utils.plt_show() diff --git a/examples/inverse/plot_dics_source_power.py b/examples/inverse/plot_dics_source_power.py index ae5f9c420ad..994fbf32a1f 100644 --- a/examples/inverse/plot_dics_source_power.py +++ b/examples/inverse/plot_dics_source_power.py @@ -54,11 +54,16 @@ # Computing the cross-spectral density matrix at 4 evenly spaced frequencies # from 6 to 10 Hz. We use a decim value of 20 to speed up the computation in # this example at the loss of accuracy. +# +# .. warning:: The use of several sensor types with the DICS beamformer is +# not heavily tested yet. Here we use verbose='error' to +# suppress a warning along these lines. csd = csd_morlet(epochs, tmin=0, tmax=0.5, decim=20, - frequencies=np.linspace(6, 10, 4)) + frequencies=np.linspace(6, 10, 4), + n_cycles=2.5) # short signals, must live with few cycles # Compute DICS spatial filter and estimate source power. -filters = make_dics(epochs.info, forward, csd, reg=0.5) +filters = make_dics(epochs.info, forward, csd, reg=0.5, verbose='error') print(filters) stc, freqs = apply_dics_csd(csd, filters) diff --git a/examples/inverse/plot_gamma_map_inverse.py b/examples/inverse/plot_gamma_map_inverse.py index fafec55cb1c..21545814120 100644 --- a/examples/inverse/plot_gamma_map_inverse.py +++ b/examples/inverse/plot_gamma_map_inverse.py @@ -44,7 +44,7 @@ # Read noise noise covariance matrix and regularize it cov = mne.read_cov(cov_fname) -cov = mne.cov.regularize(cov, evoked.info) +cov = mne.cov.regularize(cov, evoked.info, rank=None) # Run the Gamma-MAP method with dipole output alpha = 0.5 diff --git a/examples/inverse/plot_lcmv_beamformer_volume.py b/examples/inverse/plot_lcmv_beamformer_volume.py index 263a6451f01..686b38feef1 100644 --- a/examples/inverse/plot_lcmv_beamformer_volume.py +++ b/examples/inverse/plot_lcmv_beamformer_volume.py @@ -63,9 +63,10 @@ # Compute covariance matrices, fit and apply spatial filter. # Read regularized noise covariance and compute regularized data covariance -noise_cov = mne.compute_covariance(epochs, tmin=tmin, tmax=0, method='shrunk') +noise_cov = mne.compute_covariance(epochs, tmin=tmin, tmax=0, method='shrunk', + rank=None) data_cov = mne.compute_covariance(epochs, tmin=0.04, tmax=0.15, - method='shrunk') + method='shrunk', rank=None) # Compute weights of free orientation (vector) beamformer with weight # normalization (neural activity index, NAI). Providing a noise covariance diff --git a/examples/inverse/plot_mixed_source_space_inverse.py b/examples/inverse/plot_mixed_source_space_inverse.py index 374ce5ca807..db406cbcf59 100644 --- a/examples/inverse/plot_mixed_source_space_inverse.py +++ b/examples/inverse/plot_mixed_source_space_inverse.py @@ -11,17 +11,14 @@ import os.path as op import matplotlib.pyplot as plt -import mne - -from mne.datasets import sample -from mne import setup_volume_source_space -from mne import make_forward_solution -from mne.minimum_norm import make_inverse_operator, apply_inverse from nilearn import plotting +import mne +from mne.minimum_norm import make_inverse_operator, apply_inverse + # Set dir -data_path = sample.data_path() +data_path = mne.datasets.sample.data_path() subject = 'sample' data_dir = op.join(data_path, 'MEG', subject) subjects_dir = op.join(data_path, 'subjects') @@ -52,22 +49,19 @@ 'Right-Thalamus-Proper', 'Right-Cerebellum-Cortex'] -# Get a surface-based source space. We could set one up like this:: -# -# >>> src = setup_source_space(subject, fname=None, spacing='oct6', -# add_dist=False, subjects_dir=subjects_dir) -# -# But we already have one saved: - -src = mne.read_source_spaces(op.join(bem_dir, 'sample-oct-6-src.fif')) +# Get a surface-based source space, here with few source points for speed: +src = mne.setup_source_space(subject, spacing='oct5', + add_dist=False, subjects_dir=subjects_dir) # Now we create a mixed src space by adding the volume regions specified in the # list labels_vol. First, read the aseg file and the source space bounds # using the inner skull surface (here using 10mm spacing to save time): -vol_src = setup_volume_source_space( - subject, mri=fname_aseg, pos=7.0, bem=fname_model, - volume_label=labels_vol, subjects_dir=subjects_dir, verbose=True) +vol_src = mne.setup_volume_source_space( + subject, mri=fname_aseg, pos=10.0, bem=fname_model, + volume_label=labels_vol, subjects_dir=subjects_dir, + add_interpolator=False, # just for speed, usually this should be True + verbose=True) # Generate the mixed source space src += vol_src @@ -78,23 +72,23 @@ n = sum(src[i]['nuse'] for i in range(len(src))) print('the src space contains %d spaces and %d points' % (len(src), n)) +############################################################################### # We could write the mixed source space with:: # # >>> write_source_spaces(fname_mixed_src, src, overwrite=True) # +# We can also export source positions to nift file and visualize it again: -############################################################################### -# Export source positions to nift file: nii_fname = op.join(bem_dir, '%s-mixed-src.nii' % subject) src.export_volume(nii_fname, mri_resolution=True) plotting.plot_img(nii_fname, cmap='nipy_spectral') -plt.show() # Compute the fwd matrix -fwd = make_forward_solution(fname_evoked, fname_trans, src, fname_bem, - mindist=5.0, # ignore sources<=5mm from innerskull - meg=True, eeg=False, n_jobs=1) +fwd = mne.make_forward_solution( + fname_evoked, fname_trans, src, fname_bem, + mindist=5.0, # ignore sources<=5mm from innerskull + meg=True, eeg=False, n_jobs=1) leadfield = fwd['sol']['data'] print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape) @@ -110,9 +104,9 @@ noise_cov = mne.read_cov(fname_cov) # Compute inverse solution and for each epoch -snr = 3.0 # use smaller SNR for raw data -inv_method = 'MNE' # sLORETA, MNE, dSPM -parc = 'aparc' # the parcellation to use, e.g., 'aparc' 'aparc.a2009s' +snr = 3.0 # use smaller SNR for raw data +inv_method = 'dSPM' # sLORETA, MNE, dSPM +parc = 'aparc' # the parcellation to use, e.g., 'aparc' 'aparc.a2009s' lambda2 = 1.0 / snr ** 2 @@ -120,29 +114,27 @@ inverse_operator = make_inverse_operator(evoked.info, fwd, noise_cov, depth=None, fixed=False) -stcs = apply_inverse(evoked, inverse_operator, lambda2, inv_method, - pick_ori=None) +stc = apply_inverse(evoked, inverse_operator, lambda2, inv_method, + pick_ori=None) # Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi -labels_parc = mne.read_labels_from_annot(subject, parc=parc, - subjects_dir=subjects_dir) +labels_parc = mne.read_labels_from_annot( + subject, parc=parc, subjects_dir=subjects_dir) +############################################################################### # Average the source estimates within each label of the cortical parcellation # and each sub structure contained in the src space -# If mode = 'mean_flip' this option is used only for the surface cortical label + src = inverse_operator['src'] -label_ts = mne.extract_label_time_course([stcs], labels_parc, src, - mode='mean', - allow_empty=True, - return_generator=False) +label_ts = mne.extract_label_time_course( + [stc], labels_parc, src, mode='mean', allow_empty=True) # plot the times series of 2 labels fig, axes = plt.subplots(1) -axes.plot(1e3 * stcs.times, label_ts[0][0, :], 'k', label='bankssts-lh') -axes.plot(1e3 * stcs.times, label_ts[0][71, :].T, 'r', +axes.plot(1e3 * stc.times, label_ts[0][0, :], 'k', label='bankssts-lh') +axes.plot(1e3 * stc.times, label_ts[0][71, :].T, 'r', label='Brain-stem') axes.set(xlabel='Time (ms)', ylabel='MNE current (nAm)') axes.legend() mne.viz.tight_layout() -plt.show() diff --git a/examples/inverse/plot_tf_lcmv.py b/examples/inverse/plot_tf_lcmv.py index b65f15a0a2d..75754d7d82e 100644 --- a/examples/inverse/plot_tf_lcmv.py +++ b/examples/inverse/plot_tf_lcmv.py @@ -126,7 +126,7 @@ tmin=tmin_plot, tmax=tmax_plot, baseline=None, proj=True) - noise_cov = compute_covariance(epochs_band, method='shrunk') + noise_cov = compute_covariance(epochs_band, method='shrunk', rank=None) noise_covs.append(noise_cov) del raw_band # to save memory diff --git a/examples/io/plot_read_and_write_raw_data.py b/examples/io/plot_read_and_write_raw_data.py index a1709be8a5a..4d4146906b6 100644 --- a/examples/io/plot_read_and_write_raw_data.py +++ b/examples/io/plot_read_and_write_raw_data.py @@ -37,7 +37,7 @@ data, times = raw[some_picks, start:(stop + 1)] # save 150s of MEG data in FIF file -raw.save('sample_audvis_meg_raw.fif', tmin=0, tmax=150, picks=picks, +raw.save('sample_audvis_meg_trunc_raw.fif', tmin=0, tmax=150, picks=picks, overwrite=True) ############################################################################### diff --git a/examples/preprocessing/plot_otp.py b/examples/preprocessing/plot_otp.py index b549b0ba6b9..84dc46778b5 100644 --- a/examples/preprocessing/plot_otp.py +++ b/examples/preprocessing/plot_otp.py @@ -58,8 +58,8 @@ def compute_bias(raw): baseline=(None, -0.01), preload=True, verbose=False) sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=None, verbose=False) - cov = mne.compute_covariance(epochs, tmax=0, method='shrunk', - verbose=False) + cov = mne.compute_covariance(epochs, tmax=0, method='oas', + rank=None, verbose=False) idx = epochs.time_as_index(0.036)[0] data = epochs.get_data()[:, :, idx].T evoked = mne.EvokedArray(data, epochs.info, tmin=0.) diff --git a/examples/realtime/plot_compute_rt_decoder.py b/examples/realtime/plot_compute_rt_decoder.py index 4ba4b916840..f6ef66db7c0 100644 --- a/examples/realtime/plot_compute_rt_decoder.py +++ b/examples/realtime/plot_compute_rt_decoder.py @@ -5,7 +5,7 @@ Supervised machine learning applied to MEG data in sensor space. Here the classifier is updated every 5 trials and the decoding -accuracy is plotted +accuracy is plotted. """ # Authors: Mainak Jas # @@ -14,6 +14,12 @@ import numpy as np import matplotlib.pyplot as plt +from sklearn import preprocessing +from sklearn.linear_model import LogisticRegression +from sklearn.pipeline import Pipeline +from sklearn.model_selection import cross_val_score, ShuffleSplit +from mne.decoding import Vectorizer, FilterEstimator + import mne from mne.realtime import MockRtClient, RtEpochs from mne.datasets import sample @@ -52,13 +58,6 @@ # Decoding in sensor space using a linear SVM n_times = len(rt_epochs.times) -from sklearn import preprocessing # noqa -from sklearn.svm import SVC # noqa -from sklearn.pipeline import Pipeline # noqa -from sklearn.model_selection import cross_val_score, ShuffleSplit # noqa -from mne.decoding import Vectorizer, FilterEstimator # noqa - - scores_x, scores, std_scores = [], [], [] # don't highpass filter because it's epoched data and the signal length @@ -66,7 +65,7 @@ filt = FilterEstimator(rt_epochs.info, None, 40, fir_design='firwin') scaler = preprocessing.StandardScaler() vectorizer = Vectorizer() -clf = SVC(C=1, kernel='linear') +clf = LogisticRegression(solver='lbfgs') concat_classifier = Pipeline([('filter', filt), ('vector', vectorizer), ('scaler', scaler), ('svm', clf)]) @@ -83,19 +82,20 @@ plt.show(block=False) for ev_num, ev in enumerate(rt_epochs.iter_evoked()): + if ev_num >= 50: # stop at 50 + break print("Just got epoch %d" % (ev_num + 1)) if ev_num == 0: - X = ev.data[None, data_picks, :] + X = ev.data[np.newaxis, data_picks, :] y = int(ev.comment) # the comment attribute contains the event_id else: - X = np.concatenate((X, ev.data[None, data_picks, :]), axis=0) + X = np.concatenate((X, ev.data[np.newaxis, data_picks, :]), axis=0) y = np.append(y, int(ev.comment)) - if ev_num >= min_trials: - - cv = ShuffleSplit(5, test_size=0.2, random_state=42) + if ev_num >= min_trials and ev_num % 5 == 0: + cv = ShuffleSplit(5, test_size=0.2, random_state=42) # 3 for speed scores_t = cross_val_score(concat_classifier, X, y, cv=cv, n_jobs=1) * 100 diff --git a/examples/visualization/plot_xhemi.py b/examples/visualization/plot_xhemi.py index 2fc3fca930e..6d6ce1bec52 100644 --- a/examples/visualization/plot_xhemi.py +++ b/examples/visualization/plot_xhemi.py @@ -31,7 +31,8 @@ # and vice-versa. morph = mne.compute_source_morph(stc, 'fsaverage_sym', 'fsaverage_sym', spacing=stc.vertices, warn=False, - subjects_dir=subjects_dir, xhemi=True) + subjects_dir=subjects_dir, xhemi=True, + verbose='error') # creating morph map stc_xhemi = morph.apply(stc) # Now we can subtract them and plot the result: diff --git a/ignore_words.txt b/ignore_words.txt index 44e255aa050..6597c4b6283 100644 --- a/ignore_words.txt +++ b/ignore_words.txt @@ -16,3 +16,4 @@ delimeters delimeter nd cas +thes diff --git a/mne/annotations.py b/mne/annotations.py index e002c2ba92e..73502ae67f3 100644 --- a/mne/annotations.py +++ b/mne/annotations.py @@ -66,8 +66,8 @@ class Annotations(object): >>> raw.set_annotations(annotations) # doctest: +SKIP >>> epochs = mne.Epochs(raw, events, event_id, tmin, tmax) # doctest: +SKIP - orig_time - ^^^^^^^^^ + **orig_time** + If ``orig_time`` is None, the annotations are synced to the start of the data (0 seconds). Otherwise the annotations are synced to sample 0 and ``raw.first_samp`` is taken into account the same way as with events. diff --git a/mne/connectivity/effective.py b/mne/connectivity/effective.py index 7b51918075b..3a603ddc3e9 100644 --- a/mne/connectivity/effective.py +++ b/mne/connectivity/effective.py @@ -153,10 +153,10 @@ def phase_slope_index(data, indices=None, sfreq=2 * np.pi, for fi, fj in zip(freq_idx, freq_idx[1:]): idx_fi[freq_dim] = fi idx_fj[freq_dim] = fj - acc += np.conj(cohy[idx_fi]) * cohy[idx_fj] + acc += np.conj(cohy[tuple(idx_fi)]) * cohy[tuple(idx_fj)] idx_fi[freq_dim] = band_idx - psi[idx_fi] = np.imag(acc) + psi[tuple(idx_fi)] = np.imag(acc) logger.info('[PSI Estimation Done]') return psi, freqs, times, n_epochs, n_tapers diff --git a/mne/decoding/csp.py b/mne/decoding/csp.py index 759f293edfb..96a00fcf450 100644 --- a/mne/decoding/csp.py +++ b/mne/decoding/csp.py @@ -706,6 +706,10 @@ class SPoC(CSP): Parameters to pass to :func:`mne.compute_covariance`. .. versionadded:: 0.16 + rank : None | int | dict | 'full' + See :func:`mne.compute_covariance`. + + .. versionadded:: 0.17 Attributes ---------- @@ -731,11 +735,12 @@ class SPoC(CSP): """ def __init__(self, n_components=4, reg=None, log=None, - transform_into='average_power', cov_method_params=None): + transform_into='average_power', cov_method_params=None, + rank=''): """Init of SPoC.""" super(SPoC, self).__init__(n_components=n_components, reg=reg, log=log, cov_est="epoch", norm_trace=False, - transform_into=transform_into, + transform_into=transform_into, rank=rank, cov_method_params=cov_method_params) # Covariance estimation have to be done on the single epoch level, # unlike CSP where covariance estimation can also be achieved through @@ -779,7 +784,8 @@ def fit(self, X, y): covs = np.empty((n_epochs, n_channels, n_channels)) for ii, epoch in enumerate(X): covs[ii] = _regularized_covariance( - epoch, reg=self.reg, method_params=self.cov_method_params) + epoch, reg=self.reg, method_params=self.cov_method_params, + rank=self.rank) C = covs.mean(0) Cz = np.mean(covs * target[:, np.newaxis, np.newaxis], axis=0) diff --git a/mne/decoding/transformer.py b/mne/decoding/transformer.py index 3e32558729d..2788dd827d0 100644 --- a/mne/decoding/transformer.py +++ b/mne/decoding/transformer.py @@ -480,6 +480,13 @@ class FilterEstimator(TransformerMixin): See Also -------- TemporalFilter + + Notes + ----- + This is primarily meant for use in conjunction with + :class:`mne.realtime.RtEpochs`. In general it is not recommended in a + normal processing pipeline as it may result in edge artifacts. Use with + caution. """ def __init__(self, info, l_freq, h_freq, picks=None, filter_length='auto', diff --git a/mne/gui/_kit2fiff_gui.py b/mne/gui/_kit2fiff_gui.py index ed5c9dba963..4da6d21ac12 100644 --- a/mne/gui/_kit2fiff_gui.py +++ b/mne/gui/_kit2fiff_gui.py @@ -7,7 +7,6 @@ from collections import Counter import os import sys -from warnings import warn import numpy as np from scipy.linalg import inv @@ -15,7 +14,7 @@ from ..externals.six.moves import queue from ..io.meas_info import _read_dig_points, _make_dig_points -from ..utils import get_config, set_config, logger +from ..utils import get_config, set_config, logger, warn from mayavi.core.ui.mayavi_scene import MayaviScene from mayavi.tools.mlab_scene_model import MlabSceneModel diff --git a/mne/io/proj.py b/mne/io/proj.py index 3cee1dcfaa9..f43c2f42e39 100644 --- a/mne/io/proj.py +++ b/mne/io/proj.py @@ -8,6 +8,7 @@ from copy import deepcopy from itertools import count from math import sqrt +import warnings import numpy as np from scipy import linalg @@ -108,9 +109,10 @@ def plot_topomap(self, layout=None, cmap=None, sensors=True, .. versionadded:: 0.15.0 """ # noqa: E501 from ..viz.topomap import plot_projs_topomap - return plot_projs_topomap([self], layout, cmap, sensors, colorbar, - res, size, show, outlines, - contours, image_interp, axes, info) + with warnings.catch_warnings(record=True): # tight_layout fails + return plot_projs_topomap([self], layout, cmap, sensors, colorbar, + res, size, show, outlines, + contours, image_interp, axes, info) class ProjMixin(object): diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py index 275a4f5e6a0..55a8f73860f 100644 --- a/mne/time_frequency/tfr.py +++ b/mne/time_frequency/tfr.py @@ -12,7 +12,6 @@ from copy import deepcopy from functools import partial from math import sqrt -from warnings import warn import numpy as np from scipy import linalg @@ -28,7 +27,7 @@ from ..io.meas_info import Info from ..utils import SizeMixin, _is_numeric from .multitaper import dpss_windows -from ..viz.utils import (figure_nobar, plt_show, _setup_cmap, +from ..viz.utils import (figure_nobar, plt_show, _setup_cmap, warn, _connection_line, _prepare_joint_axes, _setup_vmin_vmax, _set_title_multiple_electrodes) from ..externals.h5io import write_hdf5, read_hdf5 @@ -217,7 +216,7 @@ def _cwt(X, Ws, mode="same", decim=1, use_fft=True): msg = ('At least one of the wavelets is longer than the signal. ' 'Consider padding the signal or using shorter wavelets.') if use_fft: - warn(msg) + warn(msg, UserWarning) warn_me = False # Suppress further warnings else: raise ValueError(msg) diff --git a/mne/utils.py b/mne/utils.py index 7694c4f114e..bc31a0b3f29 100644 --- a/mne/utils.py +++ b/mne/utils.py @@ -371,7 +371,7 @@ def sum_squared(X): return np.dot(X_flat, X_flat) -def warn(message, category=RuntimeWarning): +def warn(message, category=RuntimeWarning, module='mne'): """Emit a warning with trace outside the mne namespace. This function takes arguments like warnings.warn, and sends messages @@ -386,6 +386,8 @@ def warn(message, category=RuntimeWarning): Warning message. category : instance of Warning The warning class. Defaults to ``RuntimeWarning``. + module : str + The name of the module emitting the warning. """ root_dir = op.dirname(__file__) frame = None @@ -412,8 +414,9 @@ def warn(message, category=RuntimeWarning): # We need to use this instead of warn(message, category, stacklevel) # because we move out of the MNE stack, so warnings won't properly # recognize the module name (and our warnings.simplefilter will fail) - warnings.warn_explicit(message, category, fname, lineno, - 'mne', globals().get('__warningregistry__', {})) + warnings.warn_explicit( + message, category, fname, lineno, module, + globals().get('__warningregistry__', {})) # To avoid a duplicate warning print, we only emit the logger.warning if # one of the handlers is a FileHandler. See gh-5592 if any(isinstance(h, logging.FileHandler) or getattr(h, '_mne_file_like', @@ -885,7 +888,13 @@ def _update_doc(self, olddoc): if self.extra: newdoc = "%s: %s" % (newdoc, self.extra) if olddoc: - newdoc = "%s\n\n %s" % (newdoc, olddoc) + # Get the spacing right to avoid sphinx warnings + n_space = 4 + for li, line in enumerate(olddoc.split('\n')): + if li > 0 and len(line.strip()): + n_space = len(line) - len(line.lstrip()) + break + newdoc = "%s\n\n%s%s" % (newdoc, ' ' * n_space, olddoc) return newdoc diff --git a/mne/viz/raw.py b/mne/viz/raw.py index a1ab96662a3..ad8f4faef28 100644 --- a/mne/viz/raw.py +++ b/mne/viz/raw.py @@ -8,7 +8,6 @@ import copy from functools import partial -from warnings import warn import numpy as np @@ -29,7 +28,7 @@ _radio_clicked, _set_radio_button, _handle_topomap_bads, _change_channel_group, _plot_annotations, _setup_butterfly, _handle_decim, _setup_plot_projector, _check_cov, - _set_ax_label_style, _draw_vert_line) + _set_ax_label_style, _draw_vert_line, warn) from .evoked import _plot_lines @@ -632,7 +631,7 @@ def _convert_psds(psds, dB, estimate, scaling, unit, ch_names): else: msg = "Zero value in PSD for channel(s) %s. " \ "These channels might be dead." % dead_ch - warn(msg) + warn(msg, UserWarning) if estimate == 'auto': estimate = 'power' if dB else 'amplitude' diff --git a/mne/viz/utils.py b/mne/viz/utils.py index 00a0b390879..a8642decee6 100644 --- a/mne/viz/utils.py +++ b/mne/viz/utils.py @@ -20,7 +20,7 @@ from copy import deepcopy from distutils.version import LooseVersion from itertools import cycle -from warnings import catch_warnings +import warnings from ..channels.layout import _auto_topomap_coords from ..channels.channels import _contains_ch_type @@ -129,20 +129,20 @@ def tight_layout(pad=1.2, h_pad=None, w_pad=None, fig=None): fig.canvas.draw() try: # see https://github.com/matplotlib/matplotlib/issues/2654 - with catch_warnings(record=True) as ws: + with warnings.catch_warnings(record=True) as ws: fig.tight_layout(pad=pad, h_pad=h_pad, w_pad=w_pad) except Exception: try: - with catch_warnings(record=True) as ws: + with warnings.catch_warnings(record=True) as ws: fig.set_tight_layout(dict(pad=pad, h_pad=h_pad, w_pad=w_pad)) except Exception: warn('Matplotlib function "tight_layout" is not supported.' ' Skipping subplot adjustment.') return for w in ws: - w = str(w.message) if hasattr(w, 'message') else w.get_message() - if not w.startswith('This figure includes Axes'): - warn(w) + w_msg = str(w.message) if hasattr(w, 'message') else w.get_message() + if not w_msg.startswith('This figure includes Axes'): + warn(w_msg, w.category, 'matplotlib') def _check_delayed_ssp(container): @@ -170,7 +170,7 @@ def _validate_if_list_of_axes(axes, obligatory_len=None): 'dimensions however. Try using ravel or flatten ' 'method of the array.' % axes.ndim) is_correct_type = np.array([isinstance(x, mpl.axes.Axes) - for x in axes]) + for x in axes]) if not np.all(is_correct_type): first_bad = np.where(np.logical_not(is_correct_type))[0][0] raise ValueError('axes must be a list or numpy array of matplotlib ' @@ -2666,10 +2666,10 @@ def _plot_masked_image(ax, data, times, mask=None, picks=None, yvals=None, # compute bounds between time samples time_diff = np.diff(times) / 2. if len(times) > 1 else [0.0005] time_lims = np.concatenate([[times[0] - time_diff[0]], times[:-1] + - time_diff, [times[-1] + time_diff[-1]]]) + time_diff, [times[-1] + time_diff[-1]]]) log_yvals = np.concatenate([[yvals[0] / ratio[0]], yvals, - [yvals[-1] * ratio[0]]]) + [yvals[-1] * ratio[0]]]) yval_lims = np.sqrt(log_yvals[:-1] * log_yvals[1:]) # construct a time-yvaluency bounds grid diff --git a/setup.cfg b/setup.cfg index a19582d0a79..a55bbdac016 100644 --- a/setup.cfg +++ b/setup.cfg @@ -18,6 +18,7 @@ doc-files = doc addopts = --showlocals --durations=20 --doctest-modules -rs --cov-report= --doctest-ignore-import-errors # Once SciPy updates not to have non-integer and non-tuple errors (1.2.0) we # should remove them from here. +# This list should also be considered alongside reset_warnings in doc/conf.py filterwarnings = error:: ignore::ImportWarning diff --git a/tutorials/plot_ica_from_raw.py b/tutorials/plot_ica_from_raw.py index a4041fcd05d..1e36e756086 100644 --- a/tutorials/plot_ica_from_raw.py +++ b/tutorials/plot_ica_from_raw.py @@ -55,8 +55,8 @@ picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False, stim=False, exclude='bads') -ica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13), - verbose='warning') # low iterations -> does not fully converge +# low iterations -> does not fully converge +ica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13)) # maximum number of components to reject n_max_ecg, n_max_eog = 3, 1 # here we don't expect horizontal EOG components diff --git a/tutorials/plot_mne_dspm_source_localization.py b/tutorials/plot_mne_dspm_source_localization.py index 2dfe01d14b7..b4dc1222b53 100644 --- a/tutorials/plot_mne_dspm_source_localization.py +++ b/tutorials/plot_mne_dspm_source_localization.py @@ -43,7 +43,7 @@ # For more details see :ref:`tut_compute_covariance`. noise_cov = mne.compute_covariance( - epochs, tmax=0., method=['shrunk', 'empirical'], verbose=True) + epochs, tmax=0., method=['shrunk', 'empirical'], rank=None, verbose=True) fig_cov, fig_spectra = mne.viz.plot_cov(noise_cov, raw.info) diff --git a/tutorials/plot_morph_stc.py b/tutorials/plot_morph_stc.py index 32ed7e42fdc..044bb9506b6 100644 --- a/tutorials/plot_morph_stc.py +++ b/tutorials/plot_morph_stc.py @@ -100,12 +100,6 @@ - Save a SourceMorph object to disk -Here is an example result for morphing **volume** source estimates: - -:ref:`sphx_glr_auto_examples_inverse_plot_morph_volume_stc.py` - -.. image:: ../../_images/sphx_glr_plot_morph_volume_stc_001.png - Morphing **surface** source estimates ===================================== @@ -236,12 +230,6 @@ Please see also Gramfort *et al.* (2013) [4]_. -Example result for morphing **surface** source estimates: - -:ref:`sphx_glr_auto_examples_inverse_plot_morph_surface_stc.py` - -.. image:: ../../_images/sphx_glr_plot_morph_surface_stc_001.png - References ========== .. [1] Greve D. N., Van der Haegen L., Cai Q., Stufflebeam S., Sabuncu M. diff --git a/tutorials/plot_receptive_field.py b/tutorials/plot_receptive_field.py index 9f8257ae830..2863ecc55af 100644 --- a/tutorials/plot_receptive_field.py +++ b/tutorials/plot_receptive_field.py @@ -150,7 +150,7 @@ elif ix_delay < 0: take[-1] = slice(-ix_delay, None) put[-1] = slice(None, ix_delay) - X_del[ii][put] = X[take] + X_del[ii][tuple(put)] = X[tuple(take)] # Now set the delayed axis to the 2nd dimension X_del = np.rollaxis(X_del, 0, 3) diff --git a/tutorials/plot_sensors_decoding.py b/tutorials/plot_sensors_decoding.py index a15930c6851..f9f89b121f5 100644 --- a/tutorials/plot_sensors_decoding.py +++ b/tutorials/plot_sensors_decoding.py @@ -1,12 +1,30 @@ -""" -================================= -Decoding sensor space data (MVPA) -================================= +r""" +=============== +Decoding (MVPA) +=============== + +.. contents:: Contents + :local: + :depth: 3 + +.. include:: ../git_links.inc + +Design philosophy +================= +Decoding (a.k.a. MVPA) in MNE largely follows the machine +learning API of the scikit-learn package. +Each estimator implements ``fit``, ``transform``, ``fit_transform``, and +(optionally) ``inverse_transform`` methods. For more details on this design, +visit scikit-learn_. For additional theoretical insights into the decoding +framework in MNE, see [1]_. + +For ease of comprehension, we will denote instantiations of the class using +the same name as the class but in small caps instead of camel cases. -Decoding, a.k.a MVPA or supervised machine learning, is applied to MEG -data in sensor space. Here the classifier is applied to every time -point. +Let's start by loading data for a simple two-class problem: """ +# sphinx_gallery_thumbnail_number = 6 + import numpy as np import matplotlib.pyplot as plt @@ -16,29 +34,22 @@ import mne from mne.datasets import sample -from mne.decoding import (SlidingEstimator, GeneralizingEstimator, - cross_val_multiscore, LinearModel, get_coef) +from mne.decoding import (SlidingEstimator, GeneralizingEstimator, Scaler, + cross_val_multiscore, LinearModel, get_coef, + Vectorizer, CSP) data_path = sample.data_path() -plt.close('all') - -# sphinx_gallery_thumbnail_number = 4 - -############################################################################### -# Set parameters raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' tmin, tmax = -0.200, 0.500 -event_id = dict(audio_left=1, visual_left=3) - -# Setup for reading the raw data +event_id = {'Auditory/Left': 1, 'Visual/Left': 3} # just use two raw = mne.io.read_raw_fif(raw_fname, preload=True) # The subsequent decoding analyses only capture evoked responses, so we can # low-pass the MEG data. Usually a value more like 40 Hz would be used, # but here low-pass at 20 so we can more heavily decimate, and allow -# the examlpe to run faster. -raw.filter(None, 20., fir_design='firwin') +# the examlpe to run faster. The 2 Hz high-pass helps improve CSP. +raw.filter(2, 20) events = mne.find_events(raw, 'STI 014') # Set up pick list: EEG + MEG - bad channels (modify to your needs) @@ -50,20 +61,240 @@ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks, baseline=(None, 0.), preload=True, reject=dict(grad=4000e-13, eog=150e-6), decim=10) -epochs.pick_types(meg=True, exclude='bads') +epochs.pick_types(meg=True, exclude='bads') # remove stim and EOG + +X = epochs.get_data() # MEG signals: n_epochs, n_meg_channels, n_times +y = epochs.events[:, 2] # target: Audio left or right ############################################################################### +# Transformation classes +# ====================== +# +# Scaler +# ^^^^^^ +# The :class:`mne.decoding.Scaler` will standardize the data based on channel +# scales. In the simplest modes ``scalings=None`` or ``scalings=dict(...)``, +# each data channel type (e.g., mag, grad, eeg) is treated separately and +# scaled by a constant. This is the approach used by e.g., +# :func:`mne.compute_covariance` to standardize channel scales. +# +# If ``scalings='mean'`` or ``scalings='median'``, each channel is scaled using +# empirical measures. Each channel is scaled independently by the mean and +# standand deviation, or median and interquartile range, respectively, across +# all epochs and time points during :class:`~mne.decoding.Scaler.fit` +# (during training). The :meth:`~mne.decoding.Scaler.transform` method is +# called to transform data (training or test set) by scaling all time points +# and epochs on a channel-by-channel basis. To perform both the ``fit`` and +# ``transform`` operations in a single call, the +# :meth:`~mne.decoding.Scaler.fit_transform` method may be used. To invert the +# transform, :meth:`~mne.decoding.Scaler.inverse_transform` can be used. For +# ``scalings='median'``, scikit-learn_ version 0.17+ is required. +# +# .. note:: Using this class is different from directly applying +# :class:`sklearn.preprocessing.StandardScaler` or +# :class:`sklearn.preprocessing.RobustScaler` offered by +# scikit-learn_. These scale each *classification feature*, e.g. +# each time point for each channel, with mean and standard +# deviation computed across epochs, whereas +# :class:`mne.decoding.Scaler` scales each *channel* using mean and +# standard deviation computed across all of its time points +# and epochs. +# +# Vectorizer +# ^^^^^^^^^^ +# Scikit-learn API provides functionality to chain transformers and estimators +# by using :class:`sklearn.pipeline.Pipeline`. We can construct decoding +# pipelines and perform cross-validation and grid-search. However scikit-learn +# transformers and estimators generally expect 2D data +# (n_samples * n_features), whereas MNE transformers typically output data +# with a higher dimensionality +# (e.g. n_samples * n_channels * n_frequencies * n_times). A Vectorizer +# therefore needs to be applied between the MNE and the scikit-learn steps +# like: + +# Uses all MEG sensors and time points as separate classification +# features, so the resulting filters used are spatio-temporal +clf = make_pipeline(Scaler(epochs.info), + Vectorizer(), + LogisticRegression(solver='lbfgs')) + +scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=1) + +# Mean scores across cross-validation splits +score = np.mean(scores, axis=0) +print('Spatio-temporal: %0.1f%%' % (100 * score,)) + +############################################################################### +# PSDEstimator +# ^^^^^^^^^^^^ +# The :class:`mne.decoding.PSDEstimator` +# computes the power spectral density (PSD) using the multitaper +# method. It takes a 3D array as input, converts it into 2D and computes the +# PSD. +# +# FilterEstimator +# ^^^^^^^^^^^^^^^ +# The :class:`mne.decoding.FilterEstimator` filters the 3D epochs data. +# +# Spatial filters +# =============== +# +# Just like temporal filters, spatial filters provide weights to modify the +# data along the sensor dimension. They are popular in the BCI community +# because of their simplicity and ability to distinguish spatially-separated +# neural activity. +# +# Common spatial pattern +# ^^^^^^^^^^^^^^^^^^^^^^ +# +# :class:`mne.decoding.CSP` is a technique to analyze multichannel data based +# on recordings from two classes [2]_ (see also +# http://en.wikipedia.org/wiki/Common_spatial_pattern). +# +# Let :math:`X \in R^{C\times T}` be a segment of data with +# :math:`C` channels and :math:`T` time points. The data at a single time point +# is denoted by :math:`x(t)` such that :math:`X=[x(t), x(t+1), ..., x(t+T-1)]`. +# Common spatial pattern (CSP) finds a decomposition that projects the signal +# in the original sensor space to CSP space using the following transformation: +# +# .. math:: x_{CSP}(t) = W^{T}x(t) +# :label: csp +# +# where each column of :math:`W \in R^{C\times C}` is a spatial filter and each +# row of :math:`x_{CSP}` is a CSP component. The matrix :math:`W` is also +# called the de-mixing matrix in other contexts. Let +# :math:`\Sigma^{+} \in R^{C\times C}` and :math:`\Sigma^{-} \in R^{C\times C}` +# be the estimates of the covariance matrices of the two conditions. +# CSP analysis is given by the simultaneous diagonalization of the two +# covariance matrices +# +# .. math:: W^{T}\Sigma^{+}W = \lambda^{+} +# :label: diagonalize_p +# .. math:: W^{T}\Sigma^{-}W = \lambda^{-} +# :label: diagonalize_n +# +# where :math:`\lambda^{C}` is a diagonal matrix whose entries are the +# eigenvalues of the following generalized eigenvalue problem +# +# .. math:: \Sigma^{+}w = \lambda \Sigma^{-}w +# :label: eigen_problem +# +# Large entries in the diagonal matrix corresponds to a spatial filter which +# gives high variance in one class but low variance in the other. Thus, the +# filter facilitates discrimination between the two classes. +# +# .. topic:: Examples +# +# * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_csp_eeg.py` +# * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_csp_timefreq.py` +# +# .. note:: +# +# The winning entry of the Grasp-and-lift EEG competition in Kaggle used +# the :class:`~mne.decoding.CSP` implementation in MNE and was featured as +# a `script of the week`_. +# +# .. _script of the week: http://blog.kaggle.com/2015/08/12/july-2015-scripts-of-the-week/ # noqa +# +# We can use CSP with these data with: + +csp = CSP(n_components=3, norm_trace=False) +clf = make_pipeline(csp, LogisticRegression(solver='lbfgs')) +scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=1) +print('CSP: %0.1f%%' % (100 * scores.mean(),)) + +############################################################################### +# Source power comodulation (SPoC) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Source Power Comodulation (:class:`mne.decoding.SPoC`) [3]_ +# identifies the composition of +# orthogonal spatial filters that maximally correlate with a continuous target. +# +# SPoC can be seen as an extension of the CSP where the target is driven by a +# continuous variable rather than a discrete variable. Typical applications +# include extraction of motor patterns using EMG power or audio patterns using +# sound envelope. +# +# .. topic:: Examples +# +# * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_spoc_CMC.py` +# +# xDAWN +# ^^^^^ +# :class:`mne.preprocessing.Xdawn` is a spatial filtering method designed to +# improve the signal to signal + noise ratio (SSNR) of the ERP responses [4]_. +# Xdawn was originally +# designed for P300 evoked potential by enhancing the target response with +# respect to the non-target response. The implementation in MNE-Python is a +# generalization to any type of ERP. +# +# .. topic:: Examples +# +# * :ref:`sphx_glr_auto_examples_preprocessing_plot_xdawn_denoising.py` +# * :ref:`sphx_glr_auto_examples_decoding_plot_decoding_xdawn_eeg.py` +# +# Effect-matched spatial filtering +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The result of :class:`mne.decoding.EMS` is a spatial filter at each time +# point and a corresponding time course [5]_. +# Intuitively, the result gives the similarity between the filter at +# each time point and the data vector (sensors) at that time point. +# +# .. topic:: Examples +# +# * :ref:`sphx_glr_auto_examples_decoding_plot_ems_filtering.py` +# +# Patterns vs. filters +# ^^^^^^^^^^^^^^^^^^^^ +# +# When interpreting the components of the CSP (or spatial filters in general), +# it is often more intuitive to think about how :math:`x(t)` is composed of +# the different CSP components :math:`x_{CSP}(t)`. In other words, we can +# rewrite Equation :eq:`csp` as follows: +# +# .. math:: x(t) = (W^{-1})^{T}x_{CSP}(t) +# :label: patterns +# +# The columns of the matrix :math:`(W^{-1})^T` are called spatial patterns. +# This is also called the mixing matrix. The example +# :ref:`sphx_glr_auto_examples_decoding_plot_linear_model_patterns.py` +# discusses the difference between patterns and filters. +# +# These can be plotted with: + +# Fit CSP on full data and plot +csp.fit(X, y) +csp.plot_patterns(epochs.info) +csp.plot_filters(epochs.info, scalings=1e-9) + +############################################################################### +# Decoding over time +# ================== +# +# This strategy consists in fitting a multivariate predictive model on each +# time instant and evaluating its performance at the same instant on new +# epochs. The :class:`mne.decoding.SlidingEstimator` will take as input a +# pair of features :math:`X` and targets :math:`y`, where :math:`X` has +# more than 2 dimensions. For decoding over time the data :math:`X` +# is the epochs data of shape n_epochs x n_channels x n_times. As the +# last dimension of :math:`X` is the time, an estimator will be fit +# on every time instant. +# +# This approach is analogous to SlidingEstimator-based approaches in fMRI, +# where here we are interested in when one can discriminate experimental +# conditions and therefore figure out when the effect of interest happens. +# +# When working with linear models as estimators, this approach boils +# down to estimating a discriminative spatial filter for each time instant. +# # Temporal decoding -# ----------------- +# ^^^^^^^^^^^^^^^^^ # # We'll use a Logistic Regression for a binary classification as machine # learning model. # We will train the classifier on all left visual vs auditory trials on MEG -X = epochs.get_data() # MEG signals: n_epochs, n_channels, n_times -y = epochs.events[:, 2] # target: Audio left or right - clf = make_pipeline(StandardScaler(), LogisticRegression(solver='lbfgs')) time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc', verbose=True) @@ -81,7 +312,6 @@ ax.legend() ax.axvline(.0, color='k', linestyle='-') ax.set_title('Sensor space decoding') -plt.show() ############################################################################### # You can retrieve the spatial filters and spatial patterns if you explicitly @@ -99,15 +329,28 @@ **joint_kwargs) ############################################################################### -# Temporal Generalization -# ----------------------- +# Temporal generalization +# ^^^^^^^^^^^^^^^^^^^^^^^ # -# This runs the analysis used in [1]_ and further detailed in [2]_ +# Temporal generalization is an extension of the decoding over time approach. +# It consists in evaluating whether the model estimated at a particular +# time instant accurately predicts any other time instant. It is analogous to +# transferring a trained model to a distinct learning problem, where the +# problems correspond to decoding the patterns of brain activity recorded at +# distinct time instants. # -# The idea is to fit the models on each time instant and see how it -# generalizes to any other time point. +# The object to for Temporal generalization is +# :class:`mne.decoding.GeneralizingEstimator`. It expects as input :math:`X` +# and :math:`y` (similarly to :class:`~mne.decoding.SlidingEstimator`) but +# generates predictions from each model for all time instants. The class +# :class:`~mne.decoding.GeneralizingEstimator` is generic and will treat the +# last dimension as the one to be used for generalization testing. For +# convenience, here, we refer to it as different tasks. If :math:`X` +# corresponds to epochs data then the last dimension is time. +# +# This runs the analysis used in [6]_ and further detailed in [7]_: -# define the Temporal Generalization object +# define the Temporal generalization object time_gen = GeneralizingEstimator(clf, n_jobs=1, scoring='roc_auc', verbose=True) @@ -125,40 +368,64 @@ ax.legend() ax.axvline(.0, color='k', linestyle='-') ax.set_title('Decoding MEG sensors over time') -plt.show() -# Plot the full matrix +############################################################################### +# Plot the full (generalization) matrix: + fig, ax = plt.subplots(1, 1) im = ax.imshow(scores, interpolation='lanczos', origin='lower', cmap='RdBu_r', extent=epochs.times[[0, -1, 0, -1]], vmin=0., vmax=1.) ax.set_xlabel('Testing Time (s)') ax.set_ylabel('Training Time (s)') -ax.set_title('Temporal Generalization') +ax.set_title('Temporal generalization') ax.axvline(0, color='k') ax.axhline(0, color='k') plt.colorbar(im, ax=ax) -plt.show() ############################################################################### +# Source-space decoding +# ===================== +# +# Source space decoding is also possible, but because the number of features +# can be much larger than in the sensor space, univariate feature selection +# using ANOVA f-test (or some other metric) can be done to reduce the feature +# dimension. Interpreting decoding results might be easier in source space as +# compared to sensor space. +# +# .. topic:: Examples +# +# * :ref:`tut_dec_st_source` +# # Exercise -# -------- -# - Can you improve the performance using full epochs and a common spatial -# pattern (CSP) used by most BCI systems? +# ======== +# # - Explore other datasets from MNE (e.g. Face dataset from SPM to predict # Face vs. Scrambled) # -# Have a look at the example -# :ref:`sphx_glr_auto_examples_decoding_plot_decoding_csp_space.py` -# # References # ========== -# -# .. [1] Jean-Remi King, Alexandre Gramfort, Aaron Schurger, Lionel Naccache +# .. [1] Jean-Rémi King et al. (2018) "Encoding and Decoding Neuronal Dynamics: +# Methodological Framework to Uncover the Algorithms of Cognition", +# in press. https://hal.archives-ouvertes.fr/hal-01848442/ +# .. [2] Zoltan J. Koles. The quantitative extraction and topographic mapping +# of the abnormal components in the clinical EEG. Electroencephalography +# and Clinical Neurophysiology, 79(6):440--447, December 1991. +# .. [3] Dahne, S., Meinecke, F. C., Haufe, S., Hohne, J., Tangermann, M., +# Muller, K. R., & Nikulin, V. V. (2014). SPoC: a novel framework for +# relating the amplitude of neuronal oscillations to behaviorally +# relevant parameters. NeuroImage, 86, 111-122. +# .. [4] Rivet, B., Souloumiac, A., Attina, V., & Gibert, G. (2009). xDAWN +# algorithm to enhance evoked potentials: application to +# brain-computer interface. Biomedical Engineering, IEEE Transactions +# on, 56(8), 2035-2043. +# .. [5] Aaron Schurger, Sebastien Marti, and Stanislas Dehaene, "Reducing +# multi-sensor data to a single time course that reveals experimental +# effects", BMC Neuroscience 2013, 14:122 +# .. [6] Jean-Remi King, Alexandre Gramfort, Aaron Schurger, Lionel Naccache # and Stanislas Dehaene, "Two distinct dynamic modes subtend the # detection of unexpected sounds", PLOS ONE, 2013, # http://www.ncbi.nlm.nih.gov/pubmed/24475052 -# -# .. [2] King & Dehaene (2014) 'Characterizing the dynamics of mental +# .. [7] King & Dehaene (2014) 'Characterizing the dynamics of mental # representations: the temporal generalization method', Trends In # Cognitive Sciences, 18(4), 203-210. # http://www.ncbi.nlm.nih.gov/pubmed/24593982 diff --git a/tutorials/plot_stats_cluster_spatio_temporal_2samp.py b/tutorials/plot_stats_cluster_spatio_temporal_2samp.py index 06524d33737..bc7ba59ace6 100644 --- a/tutorials/plot_stats_cluster_spatio_temporal_2samp.py +++ b/tutorials/plot_stats_cluster_spatio_temporal_2samp.py @@ -85,7 +85,7 @@ print('Clustering.') T_obs, clusters, cluster_p_values, H0 = clu =\ spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=1, - threshold=f_threshold) + threshold=f_threshold, buffer_size=None) # Now select the clusters that are sig. at p < 0.05 (note that this value # is multiple-comparisons corrected). good_cluster_inds = np.where(cluster_p_values < 0.05)[0] diff --git a/tutorials/plot_stats_spatio_temporal_cluster_sensors.py b/tutorials/plot_stats_spatio_temporal_cluster_sensors.py index f0a610f704f..2a26a443f3b 100644 --- a/tutorials/plot_stats_spatio_temporal_cluster_sensors.py +++ b/tutorials/plot_stats_spatio_temporal_cluster_sensors.py @@ -101,7 +101,7 @@ cluster_stats = spatio_temporal_cluster_test(X, n_permutations=1000, threshold=threshold, tail=1, - n_jobs=1, + n_jobs=1, buffer_size=None, connectivity=connectivity) T_obs, clusters, p_values, _ = cluster_stats diff --git a/tutorials/plot_visualize_stc.py b/tutorials/plot_visualize_stc.py index 1b230b2374b..df73766e4a2 100644 --- a/tutorials/plot_visualize_stc.py +++ b/tutorials/plot_visualize_stc.py @@ -53,9 +53,10 @@ # time using ``time_viewer=True``. # # In case ``mayavi`` is not available, we also offer a ``matplotlib`` -# backend. +# backend. Here we use verbose='error' to ignore a warning that not all +# vertices were used in plotting. stc.plot(subjects_dir=subjects_dir, initial_time=initial_time, - backend='matplotlib') + backend='matplotlib', verbose='error') ############################################################################### #