Allow making inverse solutions from fixed-orientation discrete forward models, enabling multi-dipole modeling

examples/inverse/multi_dipole_model.py

@@ -0,0 +1,166 @@
+# -*- coding: utf-8 -*-
+"""
+=================================================================
+Computing source timecourses with an XFit-like multi-dipole model
+=================================================================
+
+MEGIN's XFit program offers a "guided ECD modeling" interface, where multiple
+dipoles can be fitted interactively. By manually selecting subsets of sensors
+and time ranges, dipoles can be fitted to specific signal components. Then,
+source timecourses can be computed using a multi-dipole model. The advantage of
+using a multi-dipole model over fitting each dipole in isolation, is that when
+multiple dipoles contribute to the same signal component, the model can make
+sure that activity assigned to one dipole is not also assigned to another. This
+example shows how to build a multi-dipole model for estimating source
+timecourses for evokeds or single epochs.
+
+The XFit program is the recommended approach for guided ECD modeling, because
+it offers a convenient graphical user interface for it. These dipoles can then
+be imported into MNE-Python by using the :func:`mne.read_dipole` function for
+building and applying the multi-dipole model. In addition, this example will
+also demonstrate how to perform guided ECD modeling using only MNE-Python
+functionality, which is less convenient than using XFit, but has the benefit of
+being reproducible.
+"""
+# Author: Marijn van Vliet <[email protected]>
+#
+# License: BSD-3-Clause
+
+###############################################################################
+# Importing everything and setting up the data paths for the MNE-Sample
+# dataset.
+import mne
+from mne.datasets import sample
+from mne.channels import read_vectorview_selection
+from mne.minimum_norm import (make_inverse_operator, apply_inverse,
+                              apply_inverse_epochs)
+import matplotlib.pyplot as plt
+import numpy as np
+
+data_path = sample.data_path()
+meg_path = data_path / 'MEG' / 'sample'
+raw_fname = meg_path / 'sample_audvis_raw.fif'
+cov_fname = meg_path / 'sample_audvis-shrunk-cov.fif'
+bem_dir = data_path / 'subjects' / 'sample' / 'bem'
+bem_fname =  bem_dir / 'sample-5120-5120-5120-bem-sol.fif'
+
+###############################################################################
+# Read the MEG data from the audvis experiment. Make epochs and evokeds for the
+# left and right auditory conditions.
+raw = mne.io.read_raw_fif(raw_fname)
+raw = raw.pick_types(meg=True, eog=True, stim=True)
+info = raw.info
+
+# Create epochs for auditory events
+events = mne.find_events(raw)
+event_id = dict(right=1, left=2)
+epochs = mne.Epochs(raw, events, event_id,
+                    tmin=-0.1, tmax=0.3, baseline=(None, 0),
+                    reject=dict(mag=4e-12, grad=4000e-13, eog=150e-6))
+
+# Create evokeds for left and right auditory stimulation
+evoked_left = epochs['left'].average()
+evoked_right = epochs['right'].average()
+
+###############################################################################
+# Guided dipole modeling, meaning fitting dipoles to a manually selected subset
+# of sensors as a manually chosen time, can now be performed in MEGINs XFit on
+# the evokeds we computed above. However, it is possible to do it completely
+# in MNE-Python.
+
+# Setup conductor model
+cov = mne.read_cov(cov_fname)
+bem = mne.read_bem_solution(bem_fname)
+
+# Fit two dipoles at t=80ms. The first dipole is fitted using only the sensors
+# on the left side of the helmet. The second dipole is fitted using only the
+# sensors on the right side of the helmet.
+picks_left = read_vectorview_selection('Left', info=info)
+evoked_fit_left = evoked_left.copy().crop(0.08, 0.08)
+evoked_fit_left.pick_channels(picks_left)
+cov_fit_left = cov.copy().pick_channels(picks_left)
+
+picks_right = read_vectorview_selection('Right', info=info)
+evoked_fit_right = evoked_right.copy().crop(0.08, 0.08)
+evoked_fit_right.pick_channels(picks_right)
+cov_fit_right = cov.copy().pick_channels(picks_right)
+
+# Any SSS projections that are active on this data need to be re-normalized
+# after picking channels.
+evoked_fit_left.info.normalize_proj()
+evoked_fit_right.info.normalize_proj()
+cov_fit_left['projs'] = evoked_fit_left.info['projs']
+cov_fit_right['projs'] = evoked_fit_right.info['projs']
+
+# Fit the dipoles with the subset of sensors.
+dip_left, _ = mne.fit_dipole(evoked_fit_left, cov_fit_left, bem)
+dip_right, _ = mne.fit_dipole(evoked_fit_right, cov_fit_right, bem)
+
+###############################################################################
+# Now that we have the location and orientations of the dipoles, compute the
+# full timecourses using MNE, assigning activity to both dipoles at the same
+# time while preventing leakage between the two. We use a very low ``lambda``
+# value to ensure both dipoles are fully used.
+
+fwd, _ = mne.make_forward_dipole([dip_left, dip_right], bem, info)
+
+# Apply MNE inverse
+inv = make_inverse_operator(info, fwd, cov, fixed=True, depth=0)
+stc_left = apply_inverse(evoked_left, inv, method='MNE', lambda2=1E-6)
+stc_right = apply_inverse(evoked_right, inv, method='MNE', lambda2=1E-6)
+
+# Plot the timecourses of the resulting source estimate
+fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)
+axes[0].plot(stc_left.times, stc_left.data.T)
+axes[0].set_title('Left auditory stimulation')
+axes[0].legend(['Dipole 1', 'Dipole 2'])
+axes[1].plot(stc_right.times, stc_right.data.T)
+axes[1].set_title('Right auditory stimulation')
+axes[1].set_xlabel('Time (s)')
+fig.supylabel('Dipole amplitude')
+
+###############################################################################
+# We can also fit the timecourses to single epochs. Here, we do it for each
+# experimental condition separately.
+
+stcs_left = apply_inverse_epochs(epochs['left'], inv, lambda2=1E-6,
+                                 method='MNE')
+stcs_right = apply_inverse_epochs(epochs['right'], inv, lambda2=1E-6,
+                                  method='MNE')
+
+###############################################################################
+# To summarize and visualize the single-epoch dipole amplitudes, we will create
+# a detailed plot of the mean amplitude of the dipoles during different
+# experimental conditions.
+
+# Summarize the single epoch timecourses by computing the mean amplitude from
+# 60-90ms.
+amplitudes_left = []
+amplitudes_right = []
+for stc in stcs_left:
+    amplitudes_left.append(stc.crop(0.06, 0.09).mean().data)
+for stc in stcs_right:
+    amplitudes_right.append(stc.crop(0.06, 0.09).mean().data)
+amplitudes = np.vstack([amplitudes_left, amplitudes_right])
+
+# Visualize the epoch-by-epoch dipole ampltudes in a detailed figure.
+n = len(amplitudes)
+n_left = len(amplitudes_left)
+mean_left = np.mean(amplitudes_left, axis=0)
+mean_right = np.mean(amplitudes_right, axis=0)
+
+fig, ax = plt.subplots(figsize=(8, 4))
+ax.scatter(np.arange(n), amplitudes[:, 0], label='Dipole 1')
+ax.scatter(np.arange(n), amplitudes[:, 1], label='Dipole 2')
+transition_point = n_left - 0.5
+ax.plot([0, transition_point], [mean_left[0], mean_left[0]], color='C0')
+ax.plot([0, transition_point], [mean_left[1], mean_left[1]], color='C1')
+ax.plot([transition_point, n], [mean_right[0], mean_right[0]], color='C0')
+ax.plot([transition_point, n], [mean_right[1], mean_right[1]], color='C1')
+ax.axvline(transition_point, color='black')
+ax.set_xlabel('Epochs')
+ax.set_ylabel('Dipole amplitude')
+ax.legend()
+fig.suptitle('Single epoch dipole amplitudes')
+fig.text(0.30, 0.9, 'Left auditory stimulation', ha='center')
+fig.text(0.70, 0.9, 'Right auditory stimulation', ha='center')

mne/forward/_make_forward.py

@@ -685,6 +685,10 @@ def make_forward_dipole(dipole, bem, info, trans=None, n_jobs=1, verbose=None):
     -----
     .. versionadded:: 0.12.0
     """
+    if isinstance(dipole, list):
+        from ..dipole import _concatenate_dipoles  # To avoid circular import
+        dipole = _concatenate_dipoles(dipole)
+
     # Make copies to avoid mangling original dipole
     times = dipole.times.copy()
     pos = dipole.pos.copy()

mne/forward/forward.py

@@ -1025,10 +1025,10 @@ def _triage_loose(src, loose, fixed='auto'):
     del fixed
 
     for key, this_loose in loose.items():
-        if key != 'surface' and this_loose != 1:
+        if key not in ('surface', 'discrete') and this_loose != 1:
             raise ValueError(
-                'loose parameter has to be 1 or "auto" for non-surface '
-                f'source spaces, got loose["{key}"] = {this_loose}')
+                'loose parameter has to be 1 or "auto" for non-surface/'
+                f'discrete source spaces, got loose["{key}"] = {this_loose}')
         if not 0 <= this_loose <= 1:
             raise ValueError(
                 f'loose ({key}) must be between 0 and 1, got {this_loose}')

mne/forward/tests/test_make_forward.py

@@ -459,6 +459,12 @@ def test_make_forward_dipole(tmp_path):
     assert isinstance(stc, VolSourceEstimate)
     assert_allclose(stc.times, np.arange(0., 0.003, 0.001))
 
+    # Test passing a list of Dipoles instead of a single Dipole object
+    fwd2, stc2 = make_forward_dipole([dip_even_samp[0], dip_even_samp[1:]],
+                                     sphere, info, trans=fname_trans)
+    assert_array_equal(fwd['sol']['data'], fwd2['sol']['data'])
+    assert_array_equal(stc.data, stc2.data)
+
 
 @testing.requires_testing_data
 def test_make_forward_no_meg(tmp_path):

mne/inverse_sparse/tests/test_gamma_map.py

@@ -130,16 +130,14 @@ def test_gamma_map_vol_sphere():
                                     eeg=False, meg=True)
 
     alpha = 0.5
-    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
-                  loose=0, return_residual=False)
-
-    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
-                  loose=0.2, return_residual=False)
-
     stc = gamma_map(evoked, fwd, cov, alpha, tol=1e-4,
                     xyz_same_gamma=False, update_mode=2,
                     return_residual=False)
+    assert_array_almost_equal(stc.times, evoked.times, 5)
 
+    # Computing inverse with restricted orientations should also work, since
+    # we have a discrete source space.
+    stc = gamma_map(evoked, fwd, cov, alpha, loose=0.2, return_residual=False)
     assert_array_almost_equal(stc.times, evoked.times, 5)
 
     # Compare orientation obtained using fit_dipole and gamma_map

mne/inverse_sparse/tests/test_mxne_inverse.py

@@ -205,13 +205,13 @@ def test_mxne_vol_sphere():
                                     bem=sphere, eeg=False, meg=True)
 
     alpha = 80.
-    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
-                  loose=0.0, return_residual=False,
-                  maxit=3, tol=1e-8, active_set_size=10)
 
-    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
-                  loose=0.2, return_residual=False,
-                  maxit=3, tol=1e-8, active_set_size=10)
+    # Computing inverse with restricted orientations should also work, since
+    # we have a discrete source space.
+    stc = mixed_norm(evoked_l21, fwd, cov, alpha, loose=0.2,
+                     return_residual=False, maxit=3, tol=1e-8,
+                     active_set_size=10)
+    assert_array_almost_equal(stc.times, evoked_l21.times, 5)
 
     # irMxNE tests
     with catch_logging() as log:

mne/minimum_norm/inverse.py

@@ -1395,6 +1395,26 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
     loose = _triage_loose(forward['src'], loose, fixed)
     del fixed
 
+    # Figure out what kind of inverse is requested
+    fixed_inverse = all(v == 0. for v in loose.values())
+    constrained_inverse = any(v < 1. for v in loose.values())
+
+    # We only support fixed orientations for surface and discrete source
+    # spaces. Not volume or mixed.
+    if fixed_inverse:
+        if len(loose) > 1:  # Mixed source space
+            raise ValueError('Computing inverse solutions for mixed source '
+                             'spaces with fixed orientations is not '
+                             'supported.')
+        if 'volume' in loose:
+            raise ValueError('Computing inverse solutions for volume source '
+                             'spaces with fixed orientations is not '
+                             'supported.')
+    if loose.get('volume', 1) < 1:
+        raise ValueError('Computing inverse solutions with restricted '
+                         'orientations (loose < 1) is not supported for '
+                         'volume source spaces.')
+
     # Deal with "depth"
     if exp is not None:
         exp = float(exp)
@@ -1403,10 +1423,10 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
                              f'equal to 0, got {exp}')
         exp = exp or None  # alias 0. -> None
 
-    # put the forward solution in correct orientation
+    # Put the forward solution in correct orientation.
     # (delaying for the case of fixed ori with depth weighting if
     # allow_fixed_depth is True)
-    if loose.get('surface', 1.) == 0. and len(loose) == 1:
+    if fixed_inverse:
         if not is_fixed_orient(forward):
             if allow_fixed_depth:
                 # can convert now
@@ -1424,7 +1444,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
                 'Forward operator has fixed orientation and can only '
                 'be used to make a fixed-orientation inverse '
                 'operator.')
-        if loose.get('surface', 1.) < 1. and not forward['surf_ori']:
+        if constrained_inverse and not forward['surf_ori']:
             logger.info('Converting forward solution to surface orientation')
             convert_forward_solution(
                 forward, surf_ori=True, use_cps=use_cps, copy=False)
@@ -1442,7 +1462,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
             rank=rank)
 
     # Deal with fixed orientation forward / inverse
-    if loose.get('surface', 1.) == 0. and len(loose) == 1:
+    if fixed_inverse:
         orient_prior = None
         if not is_fixed_orient(forward):
             if depth_prior is not None:
@@ -1454,8 +1474,8 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
             convert_forward_solution(
                 forward, surf_ori=True, force_fixed=True,
                 use_cps=use_cps, copy=False)
-    else:
-        if loose.get('surface', 1.) < 1:
+    else:  # Free or loose orientation
+        if constrained_inverse:
             assert forward['surf_ori']
         # In theory we could have orient_prior=None for loose=1., but
         # the MNE-C code does not do this

mne/minimum_norm/tests/test_inverse.py

@@ -831,6 +831,25 @@ def test_inverse_operator_volume(evoked, tmp_path):
         apply_inverse(evoked, inv_vol, pick_ori='normal')
 
 
+def test_inverse_operator_discrete(evoked, tmp_path):
+    """Test MNE inverse computation on discrete source space."""
+    # Make discrete source space
+    src = mne.setup_volume_source_space(
+        pos=dict(rr=[[0, 0, 0.1], [0, -0.01, 0.05]],
+                 nn=[[0, 1, 0], [1, 0, 0]]),
+        bem=fname_bem)
+
+    # Perform inverse
+    fwd = mne.make_forward_solution(
+        evoked.info, mne.Transform('head', 'mri'), src, fname_bem)
+    inv = make_inverse_operator(
+        evoked.info, fwd, make_ad_hoc_cov(evoked.info), loose=0, fixed=True,
+        depth=0)
+    stc = apply_inverse(evoked, inv)
+    assert (isinstance(stc, VolSourceEstimate))
+    assert stc.data.shape == (2, len(evoked.times))
+
+
 @pytest.mark.slowtest
 @testing.requires_testing_data
 def test_io_inverse_operator(tmp_path):

mne/utils/docs.py

@@ -801,10 +801,14 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
 """
 
 docdict['dipole'] = """
-dipole : instance of Dipole
+dipole : instance of Dipole | list of Dipole
     Dipole object containing position, orientation and amplitude of
     one or more dipoles. Multiple simultaneous dipoles may be defined by
-    assigning them identical times.
+    assigning them identical times. Alternatively, multiple simultaneous
+    dipoles may also be specified as a list of Dipole objects.
+
+    .. versionchanged:: 1.1
+        Added support for a list of :class:`mne.Dipole` instances.
 """
 
 docdict['distance'] = """