Skip to content

Safe Flexible-Extent ESMV

Solver ID: SAFE_FLEX_ESMV

Usage

from invert import Solver

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

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

Overview

Flex-extent ESMV variant with a capped ('safe') adaptive contrast mechanism designed to avoid collapsing multi-source structure.

References

  1. Lukas Hecker (2025). Unpublished.
  2. Jonmohamadi, Y., Poudel, G., Innes, C., Weiss, D., Krueger, R., & Jones, R. (2014). Comparison of beamformers for EEG source signal reconstruction. Biomedical Signal Processing and Control, 14, 175-188.

API Reference

Bases: BaseSolver

FlexESMV6: FlexESMV2 + safe adaptive contrast.

Unlike FlexESMV5, this caps the contrast exponent per sample so that the K-th strongest peak (K estimated via MDL on sensor covariance) stays above the EMD support threshold. This aims to reduce dispersion without collapsing multi-source structure (which spikes EMD).

Source code in invert/solvers/beamformers/safe_flex_esmv.py 
 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
class SolverSafeFlexESMV(BaseSolver):
    """FlexESMV6: FlexESMV2 + *safe* adaptive contrast.

    Unlike FlexESMV5, this caps the contrast exponent per sample so that the
    K-th strongest peak (K estimated via MDL on sensor covariance) stays above
    the EMD support threshold. This aims to reduce dispersion without
    collapsing multi-source structure (which spikes EMD).
    """

    meta = SolverMeta(
        slug="safe_flex_esmv",
        full_name="Safe Flexible-Extent ESMV",
        category="Beamformers",
        description=(
            "Flex-extent ESMV variant with a capped ('safe') adaptive contrast "
            "mechanism designed to avoid collapsing multi-source structure."
        ),
        references=[
            "Lukas Hecker (2025). Unpublished.",
            "Jonmohamadi, Y., Poudel, G., Innes, C., Weiss, D., Krueger, R., & Jones, R. "
            "(2014). Comparison of beamformers for EEG source signal reconstruction. "
            "Biomedical Signal Processing and Control, 14, 175-188.",
        ],
    )

    def __init__(
        self,
        name: str = "SafeFlexESMV (FlexESMV2+SafeContrast) Beamformer",
        params: _SafeContrastParams | None = None,
        reduce_rank: bool = True,
        rank: str | int = "auto",
        **kwargs,
    ):
        if params is None:
            params = _SafeContrastParams()
        self.name = name
        self._scp = params
        self._inner_kwargs = dict(kwargs)
        self._inner_kwargs.pop("verbose", None)
        return super().__init__(reduce_rank=reduce_rank, rank=rank, **kwargs)

    def make_inverse_operator(self, forward, mne_obj, *args, alpha="auto", **kwargs):
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
        self._x_cached = self.unpack_data_obj(mne_obj)
        self.inverse_operators = []
        return self

    def apply_inverse_operator(self, mne_obj):  # type: ignore[override]
        base = SolverFlexESMV2(
            reduce_rank=self.reduce_rank,
            rank=self.rank,
            verbose=self.verbose,
            **self._inner_kwargs,
        )
        base.make_inverse_operator(self.forward, mne_obj, alpha="auto")
        stc = base.apply_inverse_operator(mne_obj)

        y = stc.data
        if y.ndim != 2 or y.size == 0:
            return stc

        scp = self._scp
        p = np.mean(np.abs(y), axis=1)
        pmax = float(np.max(p)) + scp.eps
        rel = np.clip(p / pmax, scp.eps, 1.0)

        patchiness = float(np.mean(rel > scp.patchiness_threshold))
        gamma_target = 1.0 + scp.gamma_base_max * float(
            np.exp(-patchiness / scp.gamma_fscale)
        )

        # Estimate K via MDL on sensor covariance and cap gamma so top-K stay above min_keep
        x = getattr(self, "_x_cached", None)
        k_hat = 1
        if x is not None:
            n_times = x.shape[1]
            Cx = (x @ x.T) / max(1, n_times)
            eigvals = np.linalg.eigvalsh(Cx)[::-1]
            k_hat = _estimate_rank_mdl(
                eigvals, n_samples=n_times, max_rank=scp.max_sources
            )

        k_hat = int(np.clip(k_hat, 1, scp.max_sources))
        r_k = (
            float(np.sort(rel)[::-1][k_hat - 1])
            if rel.size >= k_hat
            else float(np.max(rel))
        )
        if r_k < 1.0 - 1e-9:
            gamma_max_allowed = float(np.log(scp.min_keep) / np.log(r_k))
            gamma = float(
                np.clip(min(gamma_target, gamma_max_allowed), 1.0, gamma_target)
            )
        else:
            gamma = float(gamma_target)

        scale = rel ** (gamma - 1.0)
        stc.data[:] = y * scale[:, None]
        return stc

__init__ 

__init__(
    name: str = "SafeFlexESMV (FlexESMV2+SafeContrast) Beamformer",
    params: _SafeContrastParams | None = None,
    reduce_rank: bool = True,
    rank: str | int = "auto",
    **kwargs,
)
Source code in invert/solvers/beamformers/safe_flex_esmv.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def __init__(
    self,
    name: str = "SafeFlexESMV (FlexESMV2+SafeContrast) Beamformer",
    params: _SafeContrastParams | None = None,
    reduce_rank: bool = True,
    rank: str | int = "auto",
    **kwargs,
):
    if params is None:
        params = _SafeContrastParams()
    self.name = name
    self._scp = params
    self._inner_kwargs = dict(kwargs)
    self._inner_kwargs.pop("verbose", None)
    return super().__init__(reduce_rank=reduce_rank, rank=rank, **kwargs)

make_inverse_operator 

make_inverse_operator(
    forward, mne_obj, *args, alpha="auto", **kwargs
)
Source code in invert/solvers/beamformers/safe_flex_esmv.py 
65
66
67
68
69
def make_inverse_operator(self, forward, mne_obj, *args, alpha="auto", **kwargs):
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
    self._x_cached = self.unpack_data_obj(mne_obj)
    self.inverse_operators = []
    return self

apply_inverse_operator 

apply_inverse_operator(mne_obj)
Source code in invert/solvers/beamformers/safe_flex_esmv.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
def apply_inverse_operator(self, mne_obj):  # type: ignore[override]
    base = SolverFlexESMV2(
        reduce_rank=self.reduce_rank,
        rank=self.rank,
        verbose=self.verbose,
        **self._inner_kwargs,
    )
    base.make_inverse_operator(self.forward, mne_obj, alpha="auto")
    stc = base.apply_inverse_operator(mne_obj)

    y = stc.data
    if y.ndim != 2 or y.size == 0:
        return stc

    scp = self._scp
    p = np.mean(np.abs(y), axis=1)
    pmax = float(np.max(p)) + scp.eps
    rel = np.clip(p / pmax, scp.eps, 1.0)

    patchiness = float(np.mean(rel > scp.patchiness_threshold))
    gamma_target = 1.0 + scp.gamma_base_max * float(
        np.exp(-patchiness / scp.gamma_fscale)
    )

    # Estimate K via MDL on sensor covariance and cap gamma so top-K stay above min_keep
    x = getattr(self, "_x_cached", None)
    k_hat = 1
    if x is not None:
        n_times = x.shape[1]
        Cx = (x @ x.T) / max(1, n_times)
        eigvals = np.linalg.eigvalsh(Cx)[::-1]
        k_hat = _estimate_rank_mdl(
            eigvals, n_samples=n_times, max_rank=scp.max_sources
        )

    k_hat = int(np.clip(k_hat, 1, scp.max_sources))
    r_k = (
        float(np.sort(rel)[::-1][k_hat - 1])
        if rel.size >= k_hat
        else float(np.max(rel))
    )
    if r_k < 1.0 - 1e-9:
        gamma_max_allowed = float(np.log(scp.min_keep) / np.log(r_k))
        gamma = float(
            np.clip(min(gamma_target, gamma_max_allowed), 1.0, gamma_target)
        )
    else:
        gamma = float(gamma_target)

    scale = rel ** (gamma - 1.0)
    stc.data[:] = y * scale[:, None]
    return stc