Skip to content

Coherent Maximum Entropy on the Mean

Solver ID: cMEM

Usage

from invert import Solver

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

solver = Solver("cMEM")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview

Maximum-entropy-on-the-mean approach using graphical models and parcel-wise optimization (data-driven parcellation) to estimate source activity.

References

  1. Amblard, C., Lapalme, E., & Bhatt, P. (2004). Biomagnetic source detection by maximum entropy and graphical models. IEEE Transactions on Biomedical Engineering, 51(3), 427-442.

API Reference

Bases: BaseSolver

Coherent Maximum Entropy on the Mean (cMEM) source localization solver.

Parameters:

Name Type Description Default
name str

Name of the solver.

 'cMEM'
num_parcels int

Number of parcels for data-driven parcellation.

 200
max_iter int

Maximum optimization iterations per time point.

 50
batch_size int

Batch size for time-point processing.

 100
References

Amblard, C., Lapalme, E., & Bhatt, P. (2004). Biomagnetic source detection by maximum entropy and graphical models. IEEE Transactions on Biomedical Engineering, 51(3), 427-442.

Source code in invert/solvers/bayesian/cmem.py 
 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
class SolverCMEM(BaseSolver):
    """Coherent Maximum Entropy on the Mean (cMEM) source localization solver.

    Parameters
    ----------
    name : str
        Name of the solver.
    num_parcels : int
        Number of parcels for data-driven parcellation.
    max_iter : int
        Maximum optimization iterations per time point.
    batch_size : int
        Batch size for time-point processing.

    References
    ----------
    Amblard, C., Lapalme, E., & Bhatt, P. (2004). Biomagnetic source
    detection by maximum entropy and graphical models. IEEE
    Transactions on Biomedical Engineering, 51(3), 427-442.
    """

    meta = SolverMeta(
        acronym="cMEM",
        full_name="Coherent Maximum Entropy on the Mean",
        category="Bayesian",
        description=(
            "Maximum-entropy-on-the-mean approach using graphical models and "
            "parcel-wise optimization (data-driven parcellation) to estimate "
            "source activity."
        ),
        references=[
            "Amblard, C., Lapalme, E., & Bhatt, P. (2004). Biomagnetic source detection by maximum entropy and graphical models. IEEE Transactions on Biomedical Engineering, 51(3), 427-442.",
        ],
    )

    def __init__(
        self,
        name="cMEM",
        num_parcels=200,
        max_iter=50,
        batch_size=100,
        n_outer=3,
        alpha_tol=5e-2,
        **kwargs,
    ):
        self.name = name
        self.num_parcels = num_parcels
        self.max_iter = max_iter
        self.batch_size = batch_size
        self.n_outer = n_outer
        self.alpha_tol = alpha_tol
        return super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj=None,
        *args,
        alpha="auto",
        noise_cov: mne.Covariance | None = None,
        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.
        adjacency : scipy.sparse matrix, optional
            Source adjacency matrix. Computed from forward if not provided.
        positions : numpy.ndarray, optional
            Source positions (n, 3).

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

        # Compute adjacency from forward model if not provided
        if adjacency is None:
            try:
                adjacency = mne.spatial_src_adjacency(self.forward["src"], verbose=0)
            except Exception:
                adjacency = (
                    None  # triggers KMeans fallback in _data_driven_parcellation
                )
        self._adjacency = adjacency

        J, parcels = _cmem(
            data,
            self.leadfield,
            A=adjacency,
            num_parcels=self.num_parcels,
            max_iter=self.max_iter,
            batch_size=self.batch_size,
            n_outer=self.n_outer,
            alpha_tol=self.alpha_tol,
        )

        # Cache the result for immediate apply() on the same object.
        # cMEM is iterative and expensive; notebooks often do:
        #   make_inverse_operator(fwd, evoked); apply_inverse_operator(evoked)
        # so returning the cached solution avoids re-running the full solver.
        self._cmem_cached_obj = mne_obj
        self._cmem_cached_J = J

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

    def apply_inverse_operator(self, mne_obj):
        """Apply the cMEM inverse operator.

        Since cMEM computes the full source time series during
        ``make_inverse_operator``, applying the operator re-runs the
        algorithm on the new data unless the same object was used during
        ``make_inverse_operator`` (cached).

        Parameters
        ----------
        mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
            The MNE data object.

        Return
        ------
        stc : mne.SourceEstimate
            The source estimate.
        """
        J_cached = getattr(self, "_cmem_cached_J", None)
        obj_cached = getattr(self, "_cmem_cached_obj", None)
        if J_cached is not None and obj_cached is mne_obj:
            return self.source_to_object(J_cached)

        data = self.unpack_data_obj(mne_obj)
        # Skip validate_operator_data_compatibility: cMEM stores the full
        # solution in InverseOperator (not an operator matrix), so the
        # dimension check would fail.  Re-running _cmem below is the
        # intended application path.
        data = self._sensor_transform @ data

        A = getattr(self, "_adjacency", None)

        J, self.parcels = _cmem(
            data,
            self.leadfield,
            A=A,
            num_parcels=self.num_parcels,
            max_iter=self.max_iter,
            batch_size=self.batch_size,
            n_outer=self.n_outer,
            alpha_tol=self.alpha_tol,
        )

        self._cmem_cached_obj = mne_obj
        self._cmem_cached_J = J

        stc = self.source_to_object(J)
        return stc

__init__ 

__init__(
    name="cMEM",
    num_parcels=200,
    max_iter=50,
    batch_size=100,
    n_outer=3,
    alpha_tol=0.05,
    **kwargs,
)
Source code in invert/solvers/bayesian/cmem.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def __init__(
    self,
    name="cMEM",
    num_parcels=200,
    max_iter=50,
    batch_size=100,
    n_outer=3,
    alpha_tol=5e-2,
    **kwargs,
):
    self.name = name
    self.num_parcels = num_parcels
    self.max_iter = max_iter
    self.batch_size = batch_size
    self.n_outer = n_outer
    self.alpha_tol = alpha_tol
    return super().__init__(**kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: Covariance | None = None,
    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.

 None
alpha float

The regularization parameter.

 'auto'
adjacency scipy.sparse matrix

Source adjacency matrix. Computed from forward if not provided.

 None
positions ndarray

Source positions (n, 3).

 None
Return

self : object returns itself for convenience

Source code in invert/solvers/bayesian/cmem.py 
 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
def make_inverse_operator(
    self,
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: mne.Covariance | None = None,
    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.
    adjacency : scipy.sparse matrix, optional
        Source adjacency matrix. Computed from forward if not provided.
    positions : numpy.ndarray, optional
        Source positions (n, 3).

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

    # Compute adjacency from forward model if not provided
    if adjacency is None:
        try:
            adjacency = mne.spatial_src_adjacency(self.forward["src"], verbose=0)
        except Exception:
            adjacency = (
                None  # triggers KMeans fallback in _data_driven_parcellation
            )
    self._adjacency = adjacency

    J, parcels = _cmem(
        data,
        self.leadfield,
        A=adjacency,
        num_parcels=self.num_parcels,
        max_iter=self.max_iter,
        batch_size=self.batch_size,
        n_outer=self.n_outer,
        alpha_tol=self.alpha_tol,
    )

    # Cache the result for immediate apply() on the same object.
    # cMEM is iterative and expensive; notebooks often do:
    #   make_inverse_operator(fwd, evoked); apply_inverse_operator(evoked)
    # so returning the cached solution avoids re-running the full solver.
    self._cmem_cached_obj = mne_obj
    self._cmem_cached_J = J

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

apply_inverse_operator 

apply_inverse_operator(mne_obj)

Apply the cMEM inverse operator.

Since cMEM computes the full source time series during make_inverse_operator, applying the operator re-runs the algorithm on the new data unless the same object was used during make_inverse_operator (cached).

Parameters:

Name Type Description Default
mne_obj [Evoked, Epochs, Raw]

The MNE data object.

 required
Return

stc : mne.SourceEstimate The source estimate.

Source code in invert/solvers/bayesian/cmem.py 
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
def apply_inverse_operator(self, mne_obj):
    """Apply the cMEM inverse operator.

    Since cMEM computes the full source time series during
    ``make_inverse_operator``, applying the operator re-runs the
    algorithm on the new data unless the same object was used during
    ``make_inverse_operator`` (cached).

    Parameters
    ----------
    mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
        The MNE data object.

    Return
    ------
    stc : mne.SourceEstimate
        The source estimate.
    """
    J_cached = getattr(self, "_cmem_cached_J", None)
    obj_cached = getattr(self, "_cmem_cached_obj", None)
    if J_cached is not None and obj_cached is mne_obj:
        return self.source_to_object(J_cached)

    data = self.unpack_data_obj(mne_obj)
    # Skip validate_operator_data_compatibility: cMEM stores the full
    # solution in InverseOperator (not an operator matrix), so the
    # dimension check would fail.  Re-running _cmem below is the
    # intended application path.
    data = self._sensor_transform @ data

    A = getattr(self, "_adjacency", None)

    J, self.parcels = _cmem(
        data,
        self.leadfield,
        A=A,
        num_parcels=self.num_parcels,
        max_iter=self.max_iter,
        batch_size=self.batch_size,
        n_outer=self.n_outer,
        alpha_tol=self.alpha_tol,
    )

    self._cmem_cached_obj = mne_obj
    self._cmem_cached_J = J

    stc = self.source_to_object(J)
    return stc