Skip to content

Fourth-Order Extended Signal Subspace MUSIC

Solver ID: 4-ExSo-MUSIC

Usage

from invert import Solver

# fwd = ...    (mne.Forward object)
# evoked = ... (mne.Evoked object)

solver = Solver("4-ExSo-MUSIC")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview

Fourth-order statistics (quadricovariance) variant of ExSo-MUSIC for localizing sources using higher-order cumulants, useful for correlated sources.

References

  1. Albera, L., Ferréol, A., Cosandier-Rimélé, D., Merlet, I., & Wendling, F. (2008). Brain source localization using a fourth-order deflation scheme. IEEE Transactions on Biomedical Engineering, 55(2), 490–501.

API Reference

Bases: BaseSolver

4-ExSo-MUSIC source localization solver.

Extended Signal Subspace MUSIC using fourth-order statistics (quadricovariance) for improved source localization, especially with correlated sources.

Parameters:

Name Type Description Default
name str

Name of the solver.

 '4-ExSo-MUSIC'
max_disk_size int

Maximum disk size for source patches.

 500
References

Albera, L., Ferréol, A., Cosandier-Rimélé, D., Merlet, I., & Wendling, F. (2008). Brain source localization using a fourth-order deflation scheme. IEEE Transactions on Biomedical Engineering, 55(2), 490-501.

Source code in invert/solvers/music/exso_music.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
class SolverExSoMUSIC(BaseSolver):
    """4-ExSo-MUSIC source localization solver.

    Extended Signal Subspace MUSIC using fourth-order statistics
    (quadricovariance) for improved source localization, especially
    with correlated sources.

    Parameters
    ----------
    name : str
        Name of the solver.
    max_disk_size : int
        Maximum disk size for source patches.

    References
    ----------
    Albera, L., Ferréol, A., Cosandier-Rimélé, D., Merlet, I., &
    Wendling, F. (2008). Brain source localization using a
    fourth-order deflation scheme. IEEE Transactions on Biomedical
    Engineering, 55(2), 490-501.
    """

    meta = SolverMeta(
        acronym="4-ExSo-MUSIC",
        full_name="Fourth-Order Extended Signal Subspace MUSIC",
        category="Subspace Methods",
        description=(
            "Fourth-order statistics (quadricovariance) variant of ExSo-MUSIC for "
            "localizing sources using higher-order cumulants, useful for correlated "
            "sources."
        ),
        references=[
            "Albera, L., Ferréol, A., Cosandier-Rimélé, D., Merlet, I., & Wendling, F. (2008). Brain source localization using a fourth-order deflation scheme. IEEE Transactions on Biomedical Engineering, 55(2), 490–501.",
        ],
    )

    def __init__(self, name="4-ExSo-MUSIC", max_disk_size=500, **kwargs):
        self.name = name
        self.max_disk_size = max_disk_size
        return super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj,
        *args,
        alpha="auto",
        n="auto",
        adjacency=None,
        positions=None,
        **kwargs,
    ):
        """Calculate inverse operator.

        Parameters
        ----------
        forward : mne.Forward
            The mne-python Forward model instance.
        mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
            The MNE data object.
        alpha : float
            The regularization parameter.
        n : int or str
            Number of sources to estimate, or "auto".
        adjacency : numpy.ndarray, optional
            Source adjacency matrix (n, n).
        positions : numpy.ndarray, optional
            Source positions (n, 3).

        Return
        ------
        self : object returns itself for convenience
        """
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
        data = self.unpack_data_obj(mne_obj)

        if not isinstance(n, int):
            num_sources = self.estimate_n_sources(data, method=n)
        else:
            num_sources = n

        source_map, metric_map = _exso_music(
            data,
            self.leadfield,
            A=adjacency,
            P=positions,
            num_sources=num_sources,
            max_disk_size=self.max_disk_size,
        )

        self.source_map = source_map
        self.metric_map = metric_map

        # Build a WMNE-based inverse operator on the selected dipoles
        dipole_idc = np.where(source_map > 0)[0]
        n_dipoles, n_chans = self.leadfield.shape[1], self.leadfield.shape[0]
        inverse_operator = np.zeros((n_dipoles, n_chans))

        if len(dipole_idc) > 0:
            L_sel = self.leadfield[:, dipole_idc]
            W = np.diag(np.linalg.norm(L_sel, axis=0))
            inverse_operator[dipole_idc, :] = (
                np.linalg.inv(L_sel.T @ L_sel + W.T @ W) @ L_sel.T
            )

        self.inverse_operators = [
            InverseOperator(inverse_operator, self.name),
        ]
        return self

__init__ 

__init__(name='4-ExSo-MUSIC', max_disk_size=500, **kwargs)
Source code in invert/solvers/music/exso_music.py 
47
48
49
50
def __init__(self, name="4-ExSo-MUSIC", max_disk_size=500, **kwargs):
    self.name = name
    self.max_disk_size = max_disk_size
    return super().__init__(**kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj,
    *args,
    alpha="auto",
    n="auto",
    adjacency=None,
    positions=None,
    **kwargs,
)

Calculate inverse operator.

Parameters:

Name Type Description Default
forward Forward

The mne-python Forward model instance.

 required
mne_obj [Evoked, Epochs, Raw]

The MNE data object.

 required
alpha float

The regularization parameter.

 'auto'
n int or str

Number of sources to estimate, or "auto".

 'auto'
adjacency ndarray

Source adjacency matrix (n, n).

 None
positions ndarray

Source positions (n, 3).

 None
Return

self : object returns itself for convenience

Source code in invert/solvers/music/exso_music.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def make_inverse_operator(
    self,
    forward,
    mne_obj,
    *args,
    alpha="auto",
    n="auto",
    adjacency=None,
    positions=None,
    **kwargs,
):
    """Calculate inverse operator.

    Parameters
    ----------
    forward : mne.Forward
        The mne-python Forward model instance.
    mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
        The MNE data object.
    alpha : float
        The regularization parameter.
    n : int or str
        Number of sources to estimate, or "auto".
    adjacency : numpy.ndarray, optional
        Source adjacency matrix (n, n).
    positions : numpy.ndarray, optional
        Source positions (n, 3).

    Return
    ------
    self : object returns itself for convenience
    """
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
    data = self.unpack_data_obj(mne_obj)

    if not isinstance(n, int):
        num_sources = self.estimate_n_sources(data, method=n)
    else:
        num_sources = n

    source_map, metric_map = _exso_music(
        data,
        self.leadfield,
        A=adjacency,
        P=positions,
        num_sources=num_sources,
        max_disk_size=self.max_disk_size,
    )

    self.source_map = source_map
    self.metric_map = metric_map

    # Build a WMNE-based inverse operator on the selected dipoles
    dipole_idc = np.where(source_map > 0)[0]
    n_dipoles, n_chans = self.leadfield.shape[1], self.leadfield.shape[0]
    inverse_operator = np.zeros((n_dipoles, n_chans))

    if len(dipole_idc) > 0:
        L_sel = self.leadfield[:, dipole_idc]
        W = np.diag(np.linalg.norm(L_sel, axis=0))
        inverse_operator[dipole_idc, :] = (
            np.linalg.inv(L_sel.T @ L_sel + W.T @ W) @ L_sel.T
        )

    self.inverse_operators = [
        InverseOperator(inverse_operator, self.name),
    ]
    return self