MRG/ENH: apply_inverse_cov

doc/changes/latest.inc

@@ -17,6 +17,8 @@ Changelog
 
 - Improved :func:`mne.viz.plot_epochs` to label epoch counts starting from 0 `Sophie Herbst`_
 
+- Add :func:`minimum_norm.apply_inverse_cov` to compute static power by applying inverse solutions to a data covariance matrix by `Denis Engemann`_,`Luke Bloy`_, and `Eric Larson`_
+
 - Add :func:`mne.minimum_norm.resolution_metrics` to compute various resolution metrics for inverse solutions, by `Olaf Hauk`_
 
 - Add current source density :func:`mne.preprocessing.compute_current_source_density` to compute the surface Laplacian in order to reduce volume conduction in data by `Alex Rockhill`_

doc/glossary.rst

@@ -45,6 +45,11 @@ general neuroimaging concepts. If you think a term is missing, please consider
         a type, such as gradiometer, and a unit, such as Tesla/Meter that
         is used in the code base, e.g. for plotting.
 
+    DICS
+        Dynamic Imaging of Coherent Sources, a method for computing source
+        power in different frequency bands. see :ref:`ex-inverse-source-power`
+        and :func:`mne.beamformer.make_dics`.
+
     digitization
         Digitization is a procedure of recording the headshape of a subject and
         the fiducial coils (or :term:`HPI`) and/or eeg electrodes locations on
@@ -171,6 +176,11 @@ general neuroimaging concepts. If you think a term is missing, please consider
         diagrams of approximate sensor positions in top-down diagrams of the head,
         so-called topographies or topomaps).
 
+    LCMV beamformer
+        Linearly constrained minimum variance beamformer, which attempts to
+        estimate activity for a given source while suppressing cross-talk from
+        other regions, see :func:`mne.beamformer.make_lcmv`.
+
     minimum-norm estimation
         Minimum-norm estimation (abbr. ``MNE``) can be used to generate a distributed
         map of activation on a :term:`source space`, usually on a cortical surface.

doc/overview/roadmap.rst

@@ -137,6 +137,8 @@ several (mostly) separable steps:
 1. Refactor code to use traitlets
 2. GUI elements to use PyQt5 (rather than TraitsUI/pyface)
 3. 3D plotting to use our abstracted 3D viz functions rather than Mayavi
+4. Refactor distance/fitting classes to public ones to enable the example
+   from :gh:`6693`.
 
 Once this is done, we can effectively switch to a PyVista backend.
 
@@ -160,6 +162,7 @@ as well as:
 
 - `OpenNEURO <https://openneuro.org>`__
     "A free and open platform for sharing MRI, MEG, EEG, iEEG, and ECoG data."
+    See for example :gh:`6687`.
 - `Human Connectome Project Datasets <http://www.humanconnectome.org/data>`__
     Over a 3-year span (2012-2015), the Human Connectome Project (HCP) scanned
     1,200 healthy adult subjects. The available data includes MR structural

doc/python_reference.rst

@@ -651,6 +651,7 @@ Inverse Solutions
 
    InverseOperator
    apply_inverse
+   apply_inverse_cov
    apply_inverse_epochs
    apply_inverse_raw
    compute_source_psd

examples/inverse/plot_evoked_ers_source_power.py

@@ -0,0 +1,130 @@
+"""
+====================================================================
+Compute evoked ERS source power using DICS, LCMV beamfomer, and dSPM
+====================================================================
+
+Here we examine 3 ways of localizing event-related synchronization (ERS) of
+beta band activity in this dataset: :ref:`somato-dataset` using
+:term:`DICS`, :term:`LCMV beamformer`, and :term:`dSPM` applied to active and
+baseline covariance matrices.
+"""
+# Authors: Luke Bloy <[email protected]>
+#          Eric Larson <[email protected]>
+#
+# License: BSD (3-clause)
+import os.path as op
+
+import numpy as np
+import mne
+from mne.cov import compute_covariance
+from mne.datasets import somato
+from mne.time_frequency import csd_morlet
+from mne.beamformer import (make_dics, apply_dics_csd, make_lcmv,
+                            apply_lcmv_cov)
+from mne.minimum_norm import (make_inverse_operator, apply_inverse_cov)
+
+print(__doc__)
+
+###############################################################################
+# Reading the raw data and creating epochs:
+data_path = somato.data_path()
+subject = '01'
+task = 'somato'
+raw_fname = op.join(data_path, 'sub-{}'.format(subject), 'meg',
+                    'sub-{}_task-{}_meg.fif'.format(subject, task))
+
+raw = mne.io.read_raw_fif(raw_fname)
+
+# We are interested in the beta band (12-30 Hz)
+raw.load_data().filter(12, 30)
+
+# The DICS beamformer currently only supports a single sensor type.
+# We'll use the gradiometers in this example.
+picks = mne.pick_types(raw.info, meg='grad', exclude='bads')
+
+# Read epochs
+events = mne.find_events(raw)
+epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, picks=picks,
+                    preload=True)
+
+# Read forward operator and point to freesurfer subject directory
+fname_fwd = op.join(data_path, 'derivatives', 'sub-{}'.format(subject),
+                    'sub-{}_task-{}-fwd.fif'.format(subject, task))
+subjects_dir = op.join(data_path, 'derivatives', 'freesurfer', 'subjects')
+
+fwd = mne.read_forward_solution(fname_fwd)
+
+###############################################################################
+# Compute covariances
+# -------------------
+# ERS activity starts at 0.5 seconds after stimulus onset.
+
+active_win = (0.5, 1.5)
+baseline_win = (-1, 0)
+
+baseline_cov = compute_covariance(epochs, tmin=baseline_win[0],
+                                  tmax=baseline_win[1], method='shrunk',
+                                  rank=None)
+active_cov = compute_covariance(epochs, tmin=active_win[0], tmax=active_win[1],
+                                method='shrunk', rank=None)
+
+# Weighted averaging is already in the addition of covariance objects.
+common_cov = baseline_cov + active_cov
+
+
+###############################################################################
+# Compute some source estimates
+# -----------------------------
+# Here we will use DICS, LCMV beamformer, and dSPM.
+#
+# See :ref:`ex-inverse-source-power` for more information about DICS.
+
+def _gen_dics(active_win, baseline_win, epochs):
+    freqs = np.logspace(np.log10(12), np.log10(30), 9)
+    csd = csd_morlet(epochs, freqs, tmin=-1, tmax=1.5, decim=20)
+    csd_baseline = csd_morlet(epochs, freqs, tmin=baseline_win[0],
+                              tmax=baseline_win[1], decim=20)
+    csd_ers = csd_morlet(epochs, freqs, tmin=active_win[0], tmax=active_win[1],
+                         decim=20)
+    filters = make_dics(epochs.info, fwd, csd.mean(), pick_ori='max-power')
+    stc_base, freqs = apply_dics_csd(csd_baseline.mean(), filters)
+    stc_act, freqs = apply_dics_csd(csd_ers.mean(), filters)
+    stc_act /= stc_base
+    return stc_act
+
+
+# generate lcmv source estimate
+def _gen_lcmv(active_cov, baseline_cov, common_cov):
+    filters = make_lcmv(epochs.info, fwd, common_cov, reg=0.05,
+                        noise_cov=None, pick_ori='max-power')
+    stc_base = apply_lcmv_cov(baseline_cov, filters)
+    stc_act = apply_lcmv_cov(active_cov, filters)
+    stc_act /= stc_base
+    return stc_act
+
+
+# generate mne/dSPM source estimate
+def _gen_mne(active_cov, baseline_cov, common_cov, fwd, info, method='dSPM'):
+    inverse_operator = make_inverse_operator(info, fwd, common_cov)
+    stc_act = apply_inverse_cov(active_cov, info, inverse_operator,
+                                method=method, verbose=True)
+    stc_base = apply_inverse_cov(baseline_cov, info, inverse_operator,
+                                 method=method, verbose=True)
+    stc_act /= stc_base
+    return stc_act
+
+
+# Compute source estimates
+stc_dics = _gen_dics(active_win, baseline_win, epochs)
+stc_lcmv = _gen_lcmv(active_cov, baseline_cov, common_cov)
+stc_dspm = _gen_mne(active_cov, baseline_cov, common_cov, fwd, epochs.info)
+
+###############################################################################
+# Plot source estimates
+# ---------------------
+
+for method, stc in zip(['DICS', 'LCMV', 'dSPM'],
+                       [stc_dics, stc_lcmv, stc_dspm]):
+    title = '%s source power in the 12-30 Hz frequency band' % method
+    brain = stc.plot(hemi='rh', subjects_dir=subjects_dir,
+                     subject=subject, time_label=title)

examples/inverse/plot_mne_cov_power.py

@@ -0,0 +1,139 @@
+"""
+===================================================================
+Compute source power estimate by projecting the covariance with MNE
+===================================================================
+
+We can apply the MNE inverse operator to a covariance matrix to obtain
+an estimate of source power. This is computationally more efficient than first
+estimating the source timecourses and then computing their power.
+"""
+# Author: Denis A. Engemann <[email protected]>
+#         Luke Bloy <[email protected]>
+#
+# License: BSD (3-clause)
+
+import os.path as op
+import numpy as np
+
+import mne
+from mne.datasets import sample
+from mne.minimum_norm import make_inverse_operator, apply_inverse_cov
+
+data_path = sample.data_path()
+subjects_dir = data_path + '/subjects'
+raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
+raw = mne.io.read_raw_fif(raw_fname)
+
+###############################################################################
+# Compute empty-room covariance
+# -----------------------------
+# First we compute an empty-room covariance, which captures noise from the
+# sensors and environment.
+
+raw_empty_room_fname = op.join(
+    data_path, 'MEG', 'sample', 'ernoise_raw.fif')
+raw_empty_room = mne.io.read_raw_fif(raw_empty_room_fname)
+raw_empty_room.crop(0, 60)
+raw_empty_room.info['bads'] = ['MEG 2443']
+raw_empty_room.info['projs'] = raw.info['projs']
+noise_cov = mne.compute_raw_covariance(
+    raw_empty_room, method=['empirical', 'shrunk'])
+del raw_empty_room
+
+###############################################################################
+# Epoch the data
+# --------------
+
+raw.info['bads'] = ['MEG 2443', 'EEG 053']
+raw.load_data().filter(4, 12)
+events = mne.find_events(raw, stim_channel='STI 014')
+event_id = dict(aud_l=1, aud_r=2, vis_l=3, vis_r=4)
+tmin, tmax = -0.2, 0.5
+baseline = (None, 0)  # means from the first instant to t = 0
+reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)
+epochs = mne.Epochs(raw.copy().filter(4, 12), events, event_id, tmin, tmax,
+                    proj=True, picks=('meg', 'eog'), baseline=None,
+                    reject=reject, preload=True)
+del raw
+
+###############################################################################
+# Compute and plot covariances
+# ----------------------------
+# In addition to the empty-room covariance above, we compute two additional
+# covariances:
+#
+# 1. Baseline covariance, which captures signals not of interest in our
+#    analysis (e.g., sensor noise, environmental noise, physiological
+#    artifacts, and also resting-state-like brain activity / "noise").
+# 2. Data covariance, which captures our activation of interest (in addition
+#    to noise sources).
+
+base_cov = mne.compute_covariance(
+    epochs, tmin=-0.2, tmax=0, method=['shrunk', 'empirical'], rank=None,
+    verbose=True)
+data_cov = mne.compute_covariance(
+    epochs, tmin=0., tmax=0.2, method=['shrunk', 'empirical'], rank=None,
+    verbose=True)
+
+fig_noise_cov = mne.viz.plot_cov(noise_cov, epochs.info, show_svd=False)
+fig_base_cov = mne.viz.plot_cov(base_cov, epochs.info, show_svd=False)
+fig_data_cov = mne.viz.plot_cov(data_cov, epochs.info, show_svd=False)
+
+###############################################################################
+# We can also look at the covariances using topomaps, here we just show the
+# baseline and data covariances, followed by the data covariance whitened
+# by the baseline covariance:
+
+evoked = epochs.average().pick('meg')
+evoked.drop_channels(evoked.info['bads'])
+evoked.plot(time_unit='s')
+evoked.plot_topomap(times=np.linspace(0.05, 0.15, 5), ch_type='mag',
+                    time_unit='s')
+
+evoked_noise_cov = mne.EvokedArray(data=np.diag(noise_cov['data'])[:, None],
+                                   info=evoked.info)
+evoked_data_cov = mne.EvokedArray(data=np.diag(data_cov['data'])[:, None],
+                                  info=evoked.info)
+evoked_data_cov_white = mne.whiten_evoked(evoked_data_cov, noise_cov)
+
+
+def plot_cov_diag_topomap(evoked, ch_type='grad'):
+    evoked.plot_topomap(
+        ch_type=ch_type, times=[0],
+        vmin=np.min, vmax=np.max, cmap='viridis',
+        units=dict(mag='None', grad='None'),
+        scalings=dict(mag=1, grad=1),
+        cbar_fmt=None)
+
+
+plot_cov_diag_topomap(evoked_noise_cov, 'grad')
+plot_cov_diag_topomap(evoked_data_cov, 'grad')
+plot_cov_diag_topomap(evoked_data_cov_white, 'grad')
+
+###############################################################################
+# Apply inverse operator to covariance
+# ------------------------------------
+# Finally, we can construct an inverse using the empty-room noise covariance:
+
+# Read the forward solution and compute the inverse operator
+fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
+fwd = mne.read_forward_solution(fname_fwd)
+
+# make an MEG inverse operator
+info = evoked.info
+inverse_operator = make_inverse_operator(info, fwd, noise_cov,
+                                         loose=0.2, depth=0.8)
+
+###############################################################################
+# Project our data and baseline covariance to source space:
+
+stc_data = apply_inverse_cov(data_cov, evoked.info, inverse_operator,
+                             nave=len(epochs), method='dSPM', verbose=True)
+stc_base = apply_inverse_cov(base_cov, evoked.info, inverse_operator,
+                             nave=len(epochs), method='dSPM', verbose=True)
+
+###############################################################################
+# And visualize power is relative to the baseline:
+stc_data /= stc_base
+brain = stc_data.plot(subject='sample', subjects_dir=subjects_dir,
+                      clim=dict(kind='percent', lims=(50, 90, 98)))

mne/inverse_sparse/_gamma_map.py

@@ -217,7 +217,7 @@ def gamma_map(evoked, forward, noise_cov, alpha, loose="auto", depth=0.8,
     %(rank_None)s
 
         .. versionadded:: 0.18
-    %(pick_ori-vec)s
+    %(pick_ori)s
     %(verbose)s
 
     Returns

mne/inverse_sparse/mxne_inverse.py

@@ -319,7 +319,7 @@ def mixed_norm(evoked, forward, noise_cov, alpha, loose='auto', depth=0.8,
     %(rank_None)s
 
         .. versionadded:: 0.18
-    %(pick_ori-vec)s
+    %(pick_ori)s
     %(verbose)s
 
     Returns
@@ -557,7 +557,7 @@ def tf_mixed_norm(evoked, forward, noise_cov,
     %(rank_None)s
 
         .. versionadded:: 0.18
-    %(pick_ori-vec)s
+    %(pick_ori)s
     n_tfmxne_iter : int
         Number of TF-MxNE iterations. If > 1, iterative reweighting is applied.
     %(verbose)s

mne/minimum_norm/__init__.py

@@ -4,7 +4,7 @@
                       apply_inverse_raw, make_inverse_operator,
                       apply_inverse_epochs, write_inverse_operator,
                       compute_rank_inverse, prepare_inverse_operator,
-                      estimate_snr)
+                      estimate_snr, apply_inverse_cov)
 from .psf_ctf import point_spread_function, cross_talk_function
 from .time_frequency import (source_band_induced_power, source_induced_power,
                              compute_source_psd, compute_source_psd_epochs)