Skip to content

Extended Signal Subspace MUSIC

Solver ID: ExSo-MUSIC

Usage

from invert import Solver

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

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

Overview

ExSo-MUSIC extended-source localization via pseudo-disk scanning. Supports 2nd-order (q=1, covariance-based, default) and 4th-order (q=2, quadricovariance) subspace estimation. q=2 can separate correlated sources but needs many more time samples.

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

ExSo-MUSIC extended-source localization (q=1 or q=2).

Scans pseudo-disks on the source-space mesh, evaluating how well the compound steering vector h(theta)^{⊗q} (sum of leadfield columns for all vertices in the disk) projects onto the signal subspace estimated from the (2q)-order cumulant matrix. Disks with lowest metric are selected as active sources and reconstructed via weighted minimum-norm.

By default uses 2nd-order statistics (q=1, covariance-based). Set q=2 for the 4th-order (quadricovariance) variant, which can handle correlated sources but needs substantially more time samples for a reliable estimate.

Notes
  • This solver is scalar-orientation. For free-orientation forwards it reduces to scalar per-location leadfields (default: orientation='pca').
  • If noise_cov is provided, the data/leadfield are whitened first. For q=1 this means subtracting an identity noise covariance in the whitened space.
Source code in invert/solvers/music/exso_music.py 
 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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
class SolverExSoMUSIC(BaseSolver):
    """ExSo-MUSIC extended-source localization (q=1 or q=2).

    Scans pseudo-disks on the source-space mesh, evaluating how well the
    compound steering vector ``h(theta)^{⊗q}`` (sum of leadfield columns for
    all vertices in the disk) projects onto the signal subspace estimated from
    the (2q)-order cumulant matrix.  Disks with lowest metric are selected as
    active sources and reconstructed via weighted minimum-norm.

    By default uses 2nd-order statistics (q=1, covariance-based).  Set q=2
    for the 4th-order (quadricovariance) variant, which can handle correlated
    sources but needs substantially more time samples for a reliable estimate.

    Notes
    -----
    - This solver is scalar-orientation. For free-orientation forwards it
      reduces to scalar per-location leadfields (default: ``orientation='pca'``).
    - If ``noise_cov`` is provided, the data/leadfield are whitened first. For
      q=1 this means subtracting an identity noise covariance in the whitened
      space.
    """

    meta = SolverMeta(
        acronym="ExSo-MUSIC",
        full_name="Extended Signal Subspace MUSIC",
        category="Subspace Methods",
        description=(
            "ExSo-MUSIC extended-source localization via pseudo-disk scanning.  "
            "Supports 2nd-order (q=1, covariance-based, default) and 4th-order "
            "(q=2, quadricovariance) subspace estimation.  q=2 can separate "
            "correlated sources but needs many more time samples."
        ),
        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="ExSo-MUSIC",
        *,
        q: int | str = "auto",
        max_disk_size: int | str = "auto",
        lambda_threshold: float | str = "auto",
        disk_hops: list[int] | str = "auto",
        sensor_rank: int | str | None = "auto",
        **kwargs,
    ):
        """Initialise ExSo-MUSIC solver.

        Parameters
        ----------
        q : int (1 or 2), or "auto"
            Order of the cumulant tensor used for subspace estimation.
            q=1 uses the ordinary (2nd-order) covariance matrix C₂ and works
            well when sources are uncorrelated.  q=2 builds the 4th-order
            cumulant matrix C₄ (quadricovariance), which can separate
            correlated sources because Gaussian noise vanishes in 4th-order
            cumulants.  The trade-off is that C₄ has dimension m² × m² and
            requires substantially more time samples for a reliable estimate.
            ``"auto"`` (default) selects q based on the number of time samples
            and sensors: uses q=2 when ``t >= 8*m`` (enough data for a stable
            4th-order cumulant estimate), otherwise q=1.
        max_disk_size : int or "auto"
            Maximum number of vertices a pseudo-disk may contain during the
            graph-expansion scan.  At each hop the disk grows by one ring of
            neighbours on the source-space mesh; expansion stops early if the
            vertex count exceeds this limit.  Prevents run-away growth on dense
            meshes (e.g. ico-5 with ~10k vertices per hemisphere).
            ``"auto"`` (default) uses 15, which captures exactly 1 hop of
            neighbours on typical icosahedral meshes (~6 neighbours per
            vertex).  This is mesh-resolution independent.
        lambda_threshold : float or "auto"
            Acceptance threshold for the ExSo-MUSIC metric Υ(E, h^{⊗q}).
            The metric ranges from 0 (perfect match to signal subspace) to 1
            (orthogonal).  When a float, every disk whose best metric across
            hop levels is ≤ λ is marked active.  When ``"auto"`` (default),
            thresholding is skipped and the ``n`` disks with lowest metric are
            selected globally via a min-heap.
        disk_hops : list of int, or "auto"
            Graph-hop distances at which the metric is evaluated during disk
            expansion.  At hop 0 the disk is just the centre vertex; at hop k
            it includes all vertices reachable within k edges on the
            source-space mesh.  The compound steering vector h is the sum of
            leadfield columns for all vertices in the current disk, so larger
            hops model broader extended sources.
            ``"auto"`` defaults to ``[0, 1, 2, 3, 4, 5, 7, 10, 15, 20]``.
            Note that ``max_disk_size`` caps the actual expansion regardless
            of the hop schedule.
        sensor_rank : int, "auto", or None
            Dimensionality reduction applied to the sensor space *before*
            computing 4th-order cumulants (q=2 only).  Reduces the C₄ matrix
            from m² × m² to r² × r² where r < m.  ``"auto"`` sets
            r = min(m, max(16, 4·n_sources)).  ``None`` disables reduction.
            Ignored when q=1.
        """
        self.name = name
        self.q = q
        self.max_disk_size = max_disk_size
        self.lambda_threshold = lambda_threshold
        self.disk_hops = disk_hops
        self.sensor_rank = sensor_rank
        # ExSo steering vectors are scalar; use PCA reduction for free-orientation
        # volume/discrete forwards unless the user explicitly overrides.
        kwargs.setdefault("orientation", "pca")
        return super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj=None,
        *args,
        alpha="auto",
        noise_cov: mne.Covariance | None = None,
        n="auto",
        q: int | None = None,
        lambda_threshold: float | str | None = None,
        adjacency=None,
        positions=None,
        disk_hops: list[int] | str | None = None,
        sensor_rank: int | str | None = None,
        source_weights=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".
        q : int
            ExSo-MUSIC order parameter (q=1 for 2nd-order, q=2 for 4th-order).
        lambda_threshold : float or "auto"
            Threshold λ in Eq. (7). If "auto", selects the best ``n`` disks globally.
        adjacency : numpy.ndarray, optional
            Source adjacency matrix (n, n).

        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

        m, t = data.shape
        q_raw = self.q if q is None else q
        q_eff = _coerce_q(q_raw, m=m, t=t)
        if q_eff not in (1, 2):
            raise ValueError(f"q must be 1 or 2, got {q_eff}")

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

        if adjacency is None:
            try:
                adjacency = build_source_adjacency(forward["src"])
            except Exception as exc:
                logger.warning(
                    "Could not build source adjacency from forward['src'] (%s). "
                    "Falling back to index-neighborhood disks.",
                    exc,
                )
                adjacency = None
        else:
            # If a free-orientation leadfield slips through (e.g. volume forward with
            # orientation='auto'), reduce to scalar so adjacency/disk scanning work.
            try:
                n_adj = int(getattr(adjacency, "shape", (0, 0))[0])
            except Exception:
                n_adj = 0
            if n_adj > 0 and self.leadfield.shape[1] == 3 * n_adj:
                logger.info(
                    "%s: reducing free-orientation leadfield (m x 3n) -> (m x n) "
                    "using PCA orientations for ExSo-MUSIC compatibility.",
                    self.name,
                )
                q_loc = estimate_orientation_pca(
                    self.leadfield,
                    data,
                    reg=float(self.orientation_pca_reg),
                    deterministic_sign=bool(self.orientation_pca_deterministic_sign),
                )
                self.leadfield = reduce_free_to_scalar(self.leadfield, q_loc)

        lam_eff = (
            self.lambda_threshold if lambda_threshold is None else lambda_threshold
        )
        hops_eff = self.disk_hops if disk_hops is None else disk_hops
        rank_eff = self.sensor_rank if sensor_rank is None else sensor_rank

        L_scan = self.leadfield
        if source_weights is None:
            source_weights = _vertex_area_weights_from_forward(forward)
        if source_weights is not None:
            w = np.asarray(source_weights, dtype=float).reshape(-1)
            if w.shape[0] == L_scan.shape[1]:
                L_scan = L_scan * w[np.newaxis, :]
            else:
                logger.warning(
                    "%s: source_weights has shape %s, expected (%d,); ignoring.",
                    self.name,
                    w.shape,
                    int(L_scan.shape[1]),
                )

        max_disk_eff = _coerce_max_disk_size(self.max_disk_size)

        source_map, metric_map = _exso_music(
            data,
            L_scan,
            q=q_eff,
            num_sources=num_sources,
            max_disk_size=max_disk_eff,
            adjacency=adjacency,
            lambda_threshold=lam_eff,
            disk_hops=hops_eff,
            sensor_rank=rank_eff,
            noise_cov_identity=bool(noise_cov is not None),
        )

        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 @ wf.sensor_transform, self.name),
        ]
        return self

__init__ 

__init__(
    name="ExSo-MUSIC",
    *,
    q: int | str = "auto",
    max_disk_size: int | str = "auto",
    lambda_threshold: float | str = "auto",
    disk_hops: list[int] | str = "auto",
    sensor_rank: int | str | None = "auto",
    **kwargs,
)

Initialise ExSo-MUSIC solver.

Parameters:

Name Type Description Default
q int (1 or 2), or "auto"

Order of the cumulant tensor used for subspace estimation. q=1 uses the ordinary (2nd-order) covariance matrix C₂ and works well when sources are uncorrelated. q=2 builds the 4th-order cumulant matrix C₄ (quadricovariance), which can separate correlated sources because Gaussian noise vanishes in 4th-order cumulants. The trade-off is that C₄ has dimension m² × m² and requires substantially more time samples for a reliable estimate. "auto" (default) selects q based on the number of time samples and sensors: uses q=2 when t >= 8*m (enough data for a stable 4th-order cumulant estimate), otherwise q=1.

 'auto'
max_disk_size int or auto

Maximum number of vertices a pseudo-disk may contain during the graph-expansion scan. At each hop the disk grows by one ring of neighbours on the source-space mesh; expansion stops early if the vertex count exceeds this limit. Prevents run-away growth on dense meshes (e.g. ico-5 with ~10k vertices per hemisphere). "auto" (default) uses 15, which captures exactly 1 hop of neighbours on typical icosahedral meshes (~6 neighbours per vertex). This is mesh-resolution independent.

 'auto'
lambda_threshold float or auto

Acceptance threshold for the ExSo-MUSIC metric Υ(E, h^{⊗q}). The metric ranges from 0 (perfect match to signal subspace) to 1 (orthogonal). When a float, every disk whose best metric across hop levels is ≤ λ is marked active. When "auto" (default), thresholding is skipped and the n disks with lowest metric are selected globally via a min-heap.

 'auto'
disk_hops list of int, or "auto"

Graph-hop distances at which the metric is evaluated during disk expansion. At hop 0 the disk is just the centre vertex; at hop k it includes all vertices reachable within k edges on the source-space mesh. The compound steering vector h is the sum of leadfield columns for all vertices in the current disk, so larger hops model broader extended sources. "auto" defaults to [0, 1, 2, 3, 4, 5, 7, 10, 15, 20]. Note that max_disk_size caps the actual expansion regardless of the hop schedule.

 'auto'
sensor_rank int, "auto", or None

Dimensionality reduction applied to the sensor space before computing 4th-order cumulants (q=2 only). Reduces the C₄ matrix from m² × m² to r² × r² where r < m. "auto" sets r = min(m, max(16, 4·n_sources)). None disables reduction. Ignored when q=1.

 'auto'
Source code in invert/solvers/music/exso_music.py 
 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
def __init__(
    self,
    name="ExSo-MUSIC",
    *,
    q: int | str = "auto",
    max_disk_size: int | str = "auto",
    lambda_threshold: float | str = "auto",
    disk_hops: list[int] | str = "auto",
    sensor_rank: int | str | None = "auto",
    **kwargs,
):
    """Initialise ExSo-MUSIC solver.

    Parameters
    ----------
    q : int (1 or 2), or "auto"
        Order of the cumulant tensor used for subspace estimation.
        q=1 uses the ordinary (2nd-order) covariance matrix C₂ and works
        well when sources are uncorrelated.  q=2 builds the 4th-order
        cumulant matrix C₄ (quadricovariance), which can separate
        correlated sources because Gaussian noise vanishes in 4th-order
        cumulants.  The trade-off is that C₄ has dimension m² × m² and
        requires substantially more time samples for a reliable estimate.
        ``"auto"`` (default) selects q based on the number of time samples
        and sensors: uses q=2 when ``t >= 8*m`` (enough data for a stable
        4th-order cumulant estimate), otherwise q=1.
    max_disk_size : int or "auto"
        Maximum number of vertices a pseudo-disk may contain during the
        graph-expansion scan.  At each hop the disk grows by one ring of
        neighbours on the source-space mesh; expansion stops early if the
        vertex count exceeds this limit.  Prevents run-away growth on dense
        meshes (e.g. ico-5 with ~10k vertices per hemisphere).
        ``"auto"`` (default) uses 15, which captures exactly 1 hop of
        neighbours on typical icosahedral meshes (~6 neighbours per
        vertex).  This is mesh-resolution independent.
    lambda_threshold : float or "auto"
        Acceptance threshold for the ExSo-MUSIC metric Υ(E, h^{⊗q}).
        The metric ranges from 0 (perfect match to signal subspace) to 1
        (orthogonal).  When a float, every disk whose best metric across
        hop levels is ≤ λ is marked active.  When ``"auto"`` (default),
        thresholding is skipped and the ``n`` disks with lowest metric are
        selected globally via a min-heap.
    disk_hops : list of int, or "auto"
        Graph-hop distances at which the metric is evaluated during disk
        expansion.  At hop 0 the disk is just the centre vertex; at hop k
        it includes all vertices reachable within k edges on the
        source-space mesh.  The compound steering vector h is the sum of
        leadfield columns for all vertices in the current disk, so larger
        hops model broader extended sources.
        ``"auto"`` defaults to ``[0, 1, 2, 3, 4, 5, 7, 10, 15, 20]``.
        Note that ``max_disk_size`` caps the actual expansion regardless
        of the hop schedule.
    sensor_rank : int, "auto", or None
        Dimensionality reduction applied to the sensor space *before*
        computing 4th-order cumulants (q=2 only).  Reduces the C₄ matrix
        from m² × m² to r² × r² where r < m.  ``"auto"`` sets
        r = min(m, max(16, 4·n_sources)).  ``None`` disables reduction.
        Ignored when q=1.
    """
    self.name = name
    self.q = q
    self.max_disk_size = max_disk_size
    self.lambda_threshold = lambda_threshold
    self.disk_hops = disk_hops
    self.sensor_rank = sensor_rank
    # ExSo steering vectors are scalar; use PCA reduction for free-orientation
    # volume/discrete forwards unless the user explicitly overrides.
    kwargs.setdefault("orientation", "pca")
    return super().__init__(**kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: Covariance | None = None,
    n="auto",
    q: int | None = None,
    lambda_threshold: float | str | None = None,
    adjacency=None,
    positions=None,
    disk_hops: list[int] | str | None = None,
    sensor_rank: int | str | None = None,
    source_weights=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'
n int or str

Number of sources to estimate, or "auto".

 'auto'
q int

ExSo-MUSIC order parameter (q=1 for 2nd-order, q=2 for 4th-order).

 None
lambda_threshold float or auto

Threshold λ in Eq. (7). If "auto", selects the best n disks globally.

 None
adjacency ndarray

Source adjacency matrix (n, n).

 None
Return

self : object returns itself for convenience

Source code in invert/solvers/music/exso_music.py 
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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
def make_inverse_operator(
    self,
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: mne.Covariance | None = None,
    n="auto",
    q: int | None = None,
    lambda_threshold: float | str | None = None,
    adjacency=None,
    positions=None,
    disk_hops: list[int] | str | None = None,
    sensor_rank: int | str | None = None,
    source_weights=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".
    q : int
        ExSo-MUSIC order parameter (q=1 for 2nd-order, q=2 for 4th-order).
    lambda_threshold : float or "auto"
        Threshold λ in Eq. (7). If "auto", selects the best ``n`` disks globally.
    adjacency : numpy.ndarray, optional
        Source adjacency matrix (n, n).

    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

    m, t = data.shape
    q_raw = self.q if q is None else q
    q_eff = _coerce_q(q_raw, m=m, t=t)
    if q_eff not in (1, 2):
        raise ValueError(f"q must be 1 or 2, got {q_eff}")

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

    if adjacency is None:
        try:
            adjacency = build_source_adjacency(forward["src"])
        except Exception as exc:
            logger.warning(
                "Could not build source adjacency from forward['src'] (%s). "
                "Falling back to index-neighborhood disks.",
                exc,
            )
            adjacency = None
    else:
        # If a free-orientation leadfield slips through (e.g. volume forward with
        # orientation='auto'), reduce to scalar so adjacency/disk scanning work.
        try:
            n_adj = int(getattr(adjacency, "shape", (0, 0))[0])
        except Exception:
            n_adj = 0
        if n_adj > 0 and self.leadfield.shape[1] == 3 * n_adj:
            logger.info(
                "%s: reducing free-orientation leadfield (m x 3n) -> (m x n) "
                "using PCA orientations for ExSo-MUSIC compatibility.",
                self.name,
            )
            q_loc = estimate_orientation_pca(
                self.leadfield,
                data,
                reg=float(self.orientation_pca_reg),
                deterministic_sign=bool(self.orientation_pca_deterministic_sign),
            )
            self.leadfield = reduce_free_to_scalar(self.leadfield, q_loc)

    lam_eff = (
        self.lambda_threshold if lambda_threshold is None else lambda_threshold
    )
    hops_eff = self.disk_hops if disk_hops is None else disk_hops
    rank_eff = self.sensor_rank if sensor_rank is None else sensor_rank

    L_scan = self.leadfield
    if source_weights is None:
        source_weights = _vertex_area_weights_from_forward(forward)
    if source_weights is not None:
        w = np.asarray(source_weights, dtype=float).reshape(-1)
        if w.shape[0] == L_scan.shape[1]:
            L_scan = L_scan * w[np.newaxis, :]
        else:
            logger.warning(
                "%s: source_weights has shape %s, expected (%d,); ignoring.",
                self.name,
                w.shape,
                int(L_scan.shape[1]),
            )

    max_disk_eff = _coerce_max_disk_size(self.max_disk_size)

    source_map, metric_map = _exso_music(
        data,
        L_scan,
        q=q_eff,
        num_sources=num_sources,
        max_disk_size=max_disk_eff,
        adjacency=adjacency,
        lambda_threshold=lam_eff,
        disk_hops=hops_eff,
        sensor_rank=rank_eff,
        noise_cov_identity=bool(noise_cov is not None),
    )

    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 @ wf.sensor_transform, self.name),
    ]
    return self