Trap music (#11679)

doc/changes/0.19.inc

@@ -37,7 +37,7 @@ Changelog
 
 - Allow :meth:`mne.Annotations.crop` to support negative ``tmin`` and ``tmax`` by `Joan Massich`_
 
-- Unknown events code in GDF are now visible in the ``event_id`` by `Theodore Papadopoulo`_
+- Unknown events code in GDF are now visible in the ``event_id`` by `Théodore Papadopoulo`_
 
 - Now :func:`mne.io.read_raw_ctf` populates ``raw.annotations`` with the markers in ``MarkerFile.mrk`` if any by `Joan Massich`_
 

doc/changes/latest.inc

@@ -48,6 +48,7 @@ Enhancements
 - Add support for eyetracking data using :func:`mne.io.read_raw_eyelink` (:gh:`11152` by `Dominik Welke`_ and `Scott Huberty`_)
 - :func:`mne.channels.make_1020_channel_selections` gained a new parameter, ``return_ch_names``, to allow for easy retrieval of EEG channel names corresponding to the left, right, and midline portions of the montage (:gh:`11632` by `Richard Höchenberger`_)
 - Methods for setting the sensor types of channels (e.g., for raw data, :meth:`mne.io.Raw.set_channel_types`) gained a new parameter, ``on_unit_change``, to control behavior (raise an exception, emit a warning, or do nothing) in case the measurement unit is adjusted automatically (:gh:`11668` by `Richard Höchenberger`_)
+- :func:`mne.beamformer.trap_music` implements the TRAP-MUSIC localisation algorithm with the same signature as :func:`mne.beamformer.rap_music` (:gh:`11679` by `Théodore Papadopoulo`_)
 
 Bugs
 ~~~~

doc/changes/names.inc

@@ -508,7 +508,7 @@
 
 .. _Teon Brooks: https://teonbrooks.com
 
-.. _Theodore Papadopoulo: https://github.com/papadop
+.. _Théodore Papadopoulo: https://github.com/papadop
 
 .. _Thomas Binns: https://github.com/tsbinns
 

doc/conf.py

@@ -1481,6 +1481,7 @@ def reset_warnings(gallery_conf, fname):
     "psf_ctf_vertices_lcmv.py",
     "publication_figure.py",
     "rap_music.py",
+    "trap_music.py",
     "read_inverse.py",
     "read_neo_format.py",
     "read_noise_covariance_matrix.py",

doc/inverse.rst

@@ -74,6 +74,7 @@ Inverse Solutions
    apply_dics_epochs
    apply_dics_tfr_epochs
    rap_music
+   trap_music
    make_lcmv_resolution_matrix
 
 .. currentmodule:: mne

doc/references.bib

@@ -1190,6 +1190,17 @@ @article{Makeig1993
   year = {1993}
 }
 
+@article{Makela2018,
+  author = {Mäkelä, Niko and Stenroos, Matti and Sarvas, Jukka and Ilmoniemi, Risto J.},
+  doi = {10.1016/j.neuroimage.2017.11.013},
+  journal = {Neuroimage},
+  number = {},
+  pages = {73-83},
+  title = {Truncated RAP-MUSIC (TRAP-MUSIC) for MEG and EEG source localization},
+  volume = {167},
+  year = {2018}
+}
+
 @article{MarisOostenveld2007,
   author = {Maris, Eric and Oostenveld, Robert},
   doi = {10.1016/j.jneumeth.2007.03.024},
@@ -1259,6 +1270,17 @@ @article{MosherLeahy1999
   year = {1999}
 }
 
+@inproceedings{MosherLeahy1996,
+	title = {{EEG} and {MEG} source localization using recursively applied ({RAP}) {MUSIC}},
+	doi = {10.1109/ACSSC.1996.599135},
+	booktitle = {Conference {Record} of {The} {Thirtieth} {Asilomar} {Conference} on {Signals}, {Systems} and {Computers}},
+	author = {Mosher, J.C. and Leahy, R.M.},
+	month = nov,
+	year = {1996},
+	note = {ISSN: 1058-6393},
+	pages = {1201--1207 vol.2}
+}
+
 @inproceedings{MoukademEtAl2014,
   address = {{Lisbon}},
   author = {Moukadem, Ali and Bouguila, Zied and Abdeslam, Djaffar Ould and Dieterlen, Alain},
@@ -2426,4 +2448,4 @@ @article{TierneyEtAl2022
 	issn = {1053-8119},
 	doi = {10.1016/j.neuroimage.2022.119338},
 	author = {Tierney, Tim M. and Mellor, Stephanie nd O'Neill, George C. and Timms, Ryan C. and Barnes, Gareth R.},
-}
+}

examples/inverse/trap_music.py

@@ -0,0 +1,61 @@
+"""
+.. _ex-trap-music:
+
+=================================
+Compute Trap-Music on evoked data
+=================================
+
+Compute a Truncated Recursively Applied and Projected MUltiple Signal Classification
+(TRAP-MUSIC) :footcite:`Makela2018` on evoked data.
+"""
+
+# Author: Théodore Papadopoulo <Theodore.Papadopoulo@inria.fr>
+#
+# License: BSD-3-Clause
+
+# %%
+
+import mne
+
+from mne.datasets import sample
+from mne.beamformer import trap_music
+from mne.viz import plot_dipole_locations, plot_dipole_amplitudes
+
+print(__doc__)
+
+data_path = sample.data_path()
+subjects_dir = data_path / "subjects"
+meg_path = data_path / "MEG" / "sample"
+fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
+evoked_fname = meg_path / "sample_audvis-ave.fif"
+cov_fname = meg_path / "sample_audvis-cov.fif"
+
+# Read the evoked response and crop it
+condition = "Right Auditory"
+evoked = mne.read_evokeds(evoked_fname, condition=condition, baseline=(None, 0))
+# select N100
+evoked.crop(tmin=0.05, tmax=0.15)
+
+evoked.pick_types(meg=True, eeg=False)
+
+# Read the forward solution
+forward = mne.read_forward_solution(fwd_fname)
+
+# Read noise covariance matrix
+noise_cov = mne.read_cov(cov_fname)
+
+dipoles, residual = trap_music(
+    evoked, forward, noise_cov, n_dipoles=2, return_residual=True, verbose=True
+)
+trans = forward["mri_head_t"]
+plot_dipole_locations(dipoles, trans, "sample", subjects_dir=subjects_dir)
+plot_dipole_amplitudes(dipoles)
+
+# Plot the evoked data and the residual.
+evoked.plot(ylim=dict(grad=[-300, 300], mag=[-800, 800], eeg=[-6, 8]), time_unit="s")
+residual.plot(ylim=dict(grad=[-300, 300], mag=[-800, 800], eeg=[-6, 8]), time_unit="s")
+
+# %%
+# References
+# ----------
+# .. footbibliography::

mne/beamformer/__init__.py

@@ -15,5 +15,6 @@
     apply_dics_csd,
 )
 from ._rap_music import rap_music
+from ._rap_music import trap_music
 from ._compute_beamformer import Beamformer, read_beamformer
 from .resolution_matrix import make_lcmv_resolution_matrix

mne/beamformer/_rap_music.py

@@ -12,13 +12,14 @@
 from ..inverse_sparse.mxne_inverse import _make_dipoles_sparse
 from ..minimum_norm.inverse import _log_exp_var
 from ..utils import logger, verbose, _check_info_inv, fill_doc
-from ..dipole import Dipole
 from ._compute_beamformer import _prepare_beamformer_input
 
 
 @fill_doc
-def _apply_rap_music(data, info, times, forward, noise_cov, n_dipoles=2, picks=None):
-    """RAP-MUSIC for evoked data.
+def _apply_rap_music(
+    data, info, times, forward, noise_cov, n_dipoles=2, picks=None, use_trap=False
+):
+    """RAP-MUSIC or TRAP-MUSIC for evoked data.
 
     Parameters
     ----------
@@ -35,6 +36,8 @@ def _apply_rap_music(data, info, times, forward, noise_cov, n_dipoles=2, picks=N
         The number of dipoles to estimate. The default value is 2.
     picks : list of int
         Caller ensures this is a list of int.
+    use_trap : bool
+        Use the TRAP-MUSIC variant if True (default False).
 
     Returns
     -------
@@ -43,7 +46,6 @@ def _apply_rap_music(data, info, times, forward, noise_cov, n_dipoles=2, picks=N
     explained_data : array | None
         Data explained by the dipoles using a least square fitting with the
         selected active dipoles and their estimated orientation.
-        Computed only if return_explained_data is True.
     """
     from scipy import linalg
 
@@ -114,6 +116,8 @@ def _apply_rap_music(data, info, times, forward, noise_cov, n_dipoles=2, picks=N
         projection = _compute_proj(A[:, : k + 1])
         G_proj = np.einsum("ab,bso->aso", projection, G)
         phi_sig_proj = np.dot(projection, phi_sig)
+        if use_trap:
+            phi_sig_proj = phi_sig_proj[:, -(n_dipoles - k) :]
     del G, G_proj
 
     sol = linalg.lstsq(A, M)[0]
@@ -149,39 +153,6 @@ def _apply_rap_music(data, info, times, forward, noise_cov, n_dipoles=2, picks=N
     return dipoles, explained_data
 
 
-def _make_dipoles(times, poss, oris, sol, gof):
-    """Instantiate a list of Dipoles.
-
-    Parameters
-    ----------
-    times : array, shape (n_times,)
-        The time instants.
-    poss : array, shape (n_dipoles, 3)
-        The dipoles' positions.
-    oris : array, shape (n_dipoles, 3)
-        The dipoles' orientations.
-    sol : array, shape (n_times,)
-        The dipoles' amplitudes over time.
-    gof : array, shape (n_times,)
-        The goodness of fit of the dipoles.
-        Shared between all dipoles.
-
-    Returns
-    -------
-    dipoles : list
-        The list of Dipole instances.
-    """
-    oris = np.array(oris)
-
-    dipoles = []
-    for i_dip in range(poss.shape[0]):
-        i_pos = poss[i_dip][np.newaxis, :].repeat(len(times), axis=0)
-        i_ori = oris[i_dip][np.newaxis, :].repeat(len(times), axis=0)
-        dipoles.append(Dipole(times, i_pos, sol[i_dip], i_ori, gof))
-
-    return dipoles
-
-
 def _compute_subcorr(G, phi_sig):
     """Compute the subspace correlation."""
     from scipy import linalg
@@ -210,14 +181,47 @@ def _compute_proj(A):
     return np.identity(A.shape[0]) - np.dot(U, U.T.conjugate())
 
 
+def _rap_music(evoked, forward, noise_cov, n_dipoles, return_residual, use_trap):
+    """RAP-/TRAP-MUSIC implementation."""
+    info = evoked.info
+    data = evoked.data
+    times = evoked.times
+
+    picks = _check_info_inv(info, forward, data_cov=None, noise_cov=noise_cov)
+
+    data = data[picks]
+
+    dipoles, explained_data = _apply_rap_music(
+        data, info, times, forward, noise_cov, n_dipoles, picks, use_trap
+    )
+
+    if return_residual:
+        residual = evoked.copy().pick([info["ch_names"][p] for p in picks])
+        residual.data -= explained_data
+        active_projs = [p for p in residual.info["projs"] if p["active"]]
+        for p in active_projs:
+            p["active"] = False
+        residual.add_proj(active_projs, remove_existing=True)
+        residual.apply_proj()
+        return dipoles, residual
+    else:
+        return dipoles
+
+
 @verbose
 def rap_music(
-    evoked, forward, noise_cov, n_dipoles=5, return_residual=False, verbose=None
+    evoked,
+    forward,
+    noise_cov,
+    n_dipoles=5,
+    return_residual=False,
+    *,
+    verbose=None,
 ):
     """RAP-MUSIC source localization method.
 
     Compute Recursively Applied and Projected MUltiple SIgnal Classification
-    (RAP-MUSIC) on evoked data.
+    (RAP-MUSIC) :footcite:`MosherLeahy1999,MosherLeahy1996` on evoked data.
 
     .. note:: The goodness of fit (GOF) of all the returned dipoles is the
               same and corresponds to the GOF of the full set of dipoles.
@@ -247,43 +251,70 @@ def rap_music(
     See Also
     --------
     mne.fit_dipole
+    mne.beamformer.trap_music
 
     Notes
     -----
-    The references are:
+    .. versionadded:: 0.9.0
 
-        J.C. Mosher and R.M. Leahy. 1999. Source localization using recursively
-        applied and projected (RAP) MUSIC. Signal Processing, IEEE Trans. 47, 2
-        (February 1999), 332-340.
-        DOI=10.1109/78.740118 https://doi.org/10.1109/78.740118
+    References
+    ----------
+    .. footbibliography::
+    """
+    return _rap_music(evoked, forward, noise_cov, n_dipoles, return_residual, False)
 
-        Mosher, J.C.; Leahy, R.M., EEG and MEG source localization using
-        recursively applied (RAP) MUSIC, Signals, Systems and Computers, 1996.
-        pp.1201,1207 vol.2, 3-6 Nov. 1996
-        doi: 10.1109/ACSSC.1996.599135
 
-    .. versionadded:: 0.9.0
-    """
-    info = evoked.info
-    data = evoked.data
-    times = evoked.times
+@verbose
+def trap_music(
+    evoked,
+    forward,
+    noise_cov,
+    n_dipoles=5,
+    return_residual=False,
+    *,
+    verbose=None,
+):
+    """TRAP-MUSIC source localization method.
 
-    picks = _check_info_inv(info, forward, data_cov=None, noise_cov=noise_cov)
+    Compute Truncated Recursively Applied and Projected MUltiple SIgnal Classification
+    (TRAP-MUSIC) :footcite:`Makela2018` on evoked data.
 
-    data = data[picks]
+    .. note:: The goodness of fit (GOF) of all the returned dipoles is the
+              same and corresponds to the GOF of the full set of dipoles.
 
-    dipoles, explained_data = _apply_rap_music(
-        data, info, times, forward, noise_cov, n_dipoles, picks
-    )
+    Parameters
+    ----------
+    evoked : instance of Evoked
+        Evoked data to localize.
+    forward : instance of Forward
+        Forward operator.
+    noise_cov : instance of Covariance
+        The noise covariance.
+    n_dipoles : int
+        The number of dipoles to look for. The default value is 5.
+    return_residual : bool
+        If True, the residual is returned as an Evoked instance.
+    %(verbose)s
 
-    if return_residual:
-        residual = evoked.copy().pick([info["ch_names"][p] for p in picks])
-        residual.data -= explained_data
-        active_projs = [p for p in residual.info["projs"] if p["active"]]
-        for p in active_projs:
-            p["active"] = False
-        residual.add_proj(active_projs, remove_existing=True)
-        residual.apply_proj()
-        return dipoles, residual
-    else:
-        return dipoles
+    Returns
+    -------
+    dipoles : list of instance of Dipole
+        The dipole fits.
+    residual : instance of Evoked
+        The residual a.k.a. data not explained by the dipoles.
+        Only returned if return_residual is True.
+
+    See Also
+    --------
+    mne.fit_dipole
+    mne.beamformer.rap_music
+
+    Notes
+    -----
+    .. versionadded:: 1.4
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    return _rap_music(evoked, forward, noise_cov, n_dipoles, return_residual, True)