Skip to content

Gamma MAP with MSP Priors

Solver ID: Gamma-MAP-MSP

Usage

from invert import Solver

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

solver = Solver("Gamma-MAP-MSP")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview

Gamma-MAP sparse Bayesian inverse approach combined with MSP-style spatial smoothness/patch priors.

References

  1. Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3), 947–966.

API Reference

Bases: BaseSolver

Class for the Gamma Maximum A Posteriori (Gamma-MAP) inverse solution using multiple sparse priors (MSP).

References

Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3), 947-966.

Source code in invert/solvers/bayesian/gamma_map_msp.py 
 10
 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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
class SolverGammaMAPMSP(BaseSolver):
    """Class for the Gamma Maximum A Posteriori (Gamma-MAP) inverse solution
    using multiple sparse priors (MSP).

    References
    ----------
    Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG
    source imaging. NeuroImage, 44(3), 947-966.

    """

    meta = SolverMeta(
        acronym="Gamma-MAP-MSP",
        full_name="Gamma MAP with MSP Priors",
        category="Bayesian",
        description=(
            "Gamma-MAP sparse Bayesian inverse approach combined with MSP-style "
            "spatial smoothness/patch priors."
        ),
        references=[
            "Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3), 947–966.",
        ],
    )

    def __init__(self, name="Gamma-MAP-MSP", **kwargs):
        self.name = name
        return super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj,
        *args,
        alpha="auto",
        max_iter=100,
        p=0.5,
        smoothness_order=1,
        verbose=0,
        **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.
        max_iter : int
            Maximum numbers of iterations to find the optimal hyperparameters.
            max_iter = 1 corresponds to sLORETA-like solution.
        p : 0 < p < 2
            Hyperparameter which controls sparsity. Default: p = 0
        smoothness_order : int
            Controls the smoothness prior. The higher this integer, the higher
            the pursued smoothness of the inverse solution.

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

        inverse_operators = []
        self.get_alphas(reference=data_cov)
        for alpha in self.alphas:
            inverse_operator = self.make_source_map_inverse_operator(
                data, alpha, max_iter=max_iter, p=p, smoothness_order=smoothness_order
            )
            inverse_operators.append(inverse_operator)

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

    def make_source_map_inverse_operator(
        self, B, alpha, max_iter=100, p=0.5, smoothness_order=1
    ):
        """Computes the source MAP inverse operator based on the M/EEG data.

        Parameters
        ----------
        B : numpy.ndarray
            The M/EEG data matrix (channels, time points).
        alpha : float
            The regularization parameter.
        max_iter : int
            Maximum numbers of iterations to find the optimal hyperparameters.
            max_iter = 1 corresponds to sLORETA.
        p : 0 < p < 2
            Hyperparameter which controls sparsity. Default: p = 0.5
        smoothness_order : int
            Controls the smoothness prior. The higher this integer, the higher
            the pursued smoothness of the inverse solution.

        Return
        ------
        inverse_operator : numpy.ndarray
            The inverse operator which can be used to compute inverse solutions from new data.

        """
        L = deepcopy(self.leadfield)
        db, n = B.shape
        ds = L.shape[1]

        # Ensure Common average reference
        B -= B.mean(axis=0)
        L -= L.mean(axis=0)

        # Data Covariance Matrix
        # Cb = B @ B.T
        L_smooth, gradient = self.get_smooth_prior_cov(L, smoothness_order)
        gammas = np.ones(ds)
        sigma_e = alpha * np.identity(db)
        sigma_s = np.identity(
            ds
        )  # identity leads to weighted minimum L2 Norm-type solution
        sigma_b = sigma_e + L_smooth @ sigma_s @ L_smooth.T
        sigma_b_inv = np.linalg.inv(sigma_b)

        for _k in range(max_iter):
            # print(k)
            old_gammas = deepcopy(gammas)

            # gammas = ((1/n) * np.sqrt(np.sum(( np.diag(gammas) @ L_smooth.T @ sigma_b_inv @ B )**2, axis=1)))**((2-p)/2)

            term_1 = (gammas / np.sqrt(n)) * np.sqrt(
                np.sum((L_smooth.T @ sigma_b_inv @ B) ** 2, axis=1)
            )
            term_2 = 1 / np.sqrt(np.diagonal(L_smooth.T @ sigma_b_inv @ L_smooth))
            gammas = term_1 * term_2

            if np.linalg.norm(gammas) == 0:
                gammas = old_gammas
                break
            # print(gammas.min(), gammas.max())
            # gammas /= np.linalg.norm(gammas)

        # Smooth gammas according to smooth priors
        gammas_final = abs(gammas @ gradient)
        gammas_final = gammas / gammas.max()

        sigma_s_hat = np.diag(gammas_final) @ sigma_s
        inverse_operator = (
            sigma_s_hat @ L.T @ np.linalg.inv(sigma_e + L @ sigma_s_hat @ L.T)
        )
        # S = inverse_operator @ B
        return inverse_operator

    @staticmethod
    def frob(x):
        if len(x.shape) == 1:
            x = x[:, np.newaxis]
        return np.sqrt(np.trace(x @ x.T))

    def get_smooth_prior_cov(self, L, smoothness_order):
        """Create a smooth prior on the covariance matrix.

        Parameters
        ----------
        L : numpy.ndarray
            Leadfield matrix (channels, dipoles)
        smoothness_order : int
            The higher the order, the smoother the prior.

        Return
        ------
        L : numpy.ndarray
            The smoothed Leadfield matrix (channels, dipoles)
        gradient : numpy.ndarray
            The smoothness gradient (laplacian matrix)

        """
        adjacency = mne.spatial_src_adjacency(self.forward["src"], verbose=0)
        gradient = laplacian(adjacency).toarray().astype(np.float32)

        for _ in range(smoothness_order):
            gradient = gradient @ gradient
        L = L @ abs(gradient)
        # L -= L.mean(axis=0)
        L /= np.linalg.norm(L, axis=0)
        return L, gradient

__init__ 

__init__(name='Gamma-MAP-MSP', **kwargs)
Source code in invert/solvers/bayesian/gamma_map_msp.py 
34
35
36
def __init__(self, name="Gamma-MAP-MSP", **kwargs):
    self.name = name
    return super().__init__(**kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj,
    *args,
    alpha="auto",
    max_iter=100,
    p=0.5,
    smoothness_order=1,
    verbose=0,
    **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'
max_iter int

Maximum numbers of iterations to find the optimal hyperparameters. max_iter = 1 corresponds to sLORETA-like solution.

 100
p 0 < p < 2

Hyperparameter which controls sparsity. Default: p = 0

 0.5
smoothness_order int

Controls the smoothness prior. The higher this integer, the higher the pursued smoothness of the inverse solution.

 1
Return

self : object returns itself for convenience

Source code in invert/solvers/bayesian/gamma_map_msp.py 
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
def make_inverse_operator(
    self,
    forward,
    mne_obj,
    *args,
    alpha="auto",
    max_iter=100,
    p=0.5,
    smoothness_order=1,
    verbose=0,
    **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.
    max_iter : int
        Maximum numbers of iterations to find the optimal hyperparameters.
        max_iter = 1 corresponds to sLORETA-like solution.
    p : 0 < p < 2
        Hyperparameter which controls sparsity. Default: p = 0
    smoothness_order : int
        Controls the smoothness prior. The higher this integer, the higher
        the pursued smoothness of the inverse solution.

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

    inverse_operators = []
    self.get_alphas(reference=data_cov)
    for alpha in self.alphas:
        inverse_operator = self.make_source_map_inverse_operator(
            data, alpha, max_iter=max_iter, p=p, smoothness_order=smoothness_order
        )
        inverse_operators.append(inverse_operator)

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