Merge remote-tracking branch 'upstream/main' into openmeeg

* upstream/main:
  BUG: Spectrum deprecation cleanup [circle deploy] (mne-tools#11115)
  Add API entry list and map (mne-tools#10999)
  Add legacy decorator (mne-tools#11097)
  [ENH, MRG] Add time-frequency epoch source estimation (mne-tools#11095)

Author: larsoner

doc/changes/latest.inc

@@ -38,6 +38,7 @@ Enhancements
 - Allow :func:`mne.beamformer.make_dics` to take ``pick_ori='vector'`` to compute vector source estimates (:gh:`19080` by `Alex Rockhill`_)
 - Add ``units`` parameter to :func:`mne.io.read_raw_edf` in case units are missing from the file (:gh:`11099` by `Alex Gramfort`_)
 - Add ``on_missing`` functionality to all of our classes that have a ``drop_channels`` method, to control what happens when channel names are not in the object (:gh:`11077` by `Andrew Quinn`_)
+- Add :func:`mne.minimum_norm.apply_inverse_tfr_epochs` to apply inverse methods to time-frequency resolved epochs (:gh:`11095` by `Alex Rockhill`_)
 
 Bugs
 ~~~~

doc/conf.py

@@ -468,6 +468,37 @@ def __call__(self, gallery_conf, fname, when):
     'matplotlib_animations': True,
     'compress_images': compress_images,
     'filename_pattern': '^((?!sgskip).)*$',
+    'api_usage_ignore': (
+        '('
+        '.*__.*__|'  # built-ins
+        '.*Base.*|.*Array.*|mne.Vector.*|mne.Mixed.*|mne.Vol.*|'  # inherited
+        'mne.coreg.Coregistration.*|'  # GUI
+        # common
+        '.*utils.*|.*verbose()|.*copy()|.*update()|.*save()|'
+        '.*get_data()|'
+        # mixins
+        '.*add_channels()|.*add_reference_channels()|'
+        '.*anonymize()|.*apply_baseline()|.*apply_function()|'
+        '.*apply_hilbert()|.*as_type()|.*decimate()|'
+        '.*drop()|.*drop_channels()|.*drop_log_stats()|'
+        '.*export()|.*get_channel_types()|'
+        '.*get_montage()|.*interpolate_bads()|.*next()|'
+        '.*pick()|.*pick_channels()|.*pick_types()|'
+        '.*plot_sensors()|.*rename_channels()|'
+        '.*reorder_channels()|.*savgol_filter()|'
+        '.*set_eeg_reference()|.*set_channel_types()|'
+        '.*set_meas_date()|.*set_montage()|.*shift_time()|'
+        '.*time_as_index()|.*to_data_frame()|'
+        # dictionary inherited
+        '.*clear()|.*fromkeys()|.*get()|.*items()|'
+        '.*keys()|.*pop()|.*popitem()|.*setdefault()|'
+        '.*values()|'
+        # sklearn inherited
+        '.*apply()|.*decision_function()|.*fit()|'
+        '.*fit_transform()|.*get_params()|.*predict()|'
+        '.*predict_proba()|.*set_params()|.*transform()|'
+        # I/O, also related to mixins
+        '.*.remove.*|.*.write.*)')
 }
 # Files were renamed from plot_* with:
 # find . -type f -name 'plot_*.py' -exec sh -c 'x="{}"; xn=`basename "${x}"`; git mv "$x" `dirname "${x}"`/${xn:5}' \;  # noqa

doc/inverse.rst

@@ -18,6 +18,7 @@ Inverse Solutions
    apply_inverse_cov
    apply_inverse_epochs
    apply_inverse_raw
+   apply_inverse_tfr_epochs
    compute_source_psd
    compute_source_psd_epochs
    compute_rank_inverse

examples/visualization/topo_customized.py

@@ -26,7 +26,6 @@
 import mne
 from mne.viz import iter_topography
 from mne import io
-from mne.time_frequency import psd_welch
 from mne.datasets import sample
 
 print(__doc__)
@@ -42,8 +41,9 @@
 tmin, tmax = 0, 120  # use the first 120s of data
 fmin, fmax = 2, 20  # look at frequencies between 2 and 20Hz
 n_fft = 2048  # the FFT size (n_fft). Ideally a power of 2
-psds, freqs = psd_welch(raw, picks=picks, tmin=tmin, tmax=tmax,
-                        fmin=fmin, fmax=fmax)
+spectrum = raw.compute_psd(
+    picks=picks, tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax)
+psds, freqs = spectrum.get_data(exclude=(), return_freqs=True)
 psds = 20 * np.log10(psds)  # scale to dB
 
 

mne/minimum_norm/__init__.py

@@ -2,9 +2,10 @@
 
 from .inverse import (InverseOperator, read_inverse_operator, apply_inverse,
                       apply_inverse_raw, make_inverse_operator,
-                      apply_inverse_epochs, write_inverse_operator,
-                      compute_rank_inverse, prepare_inverse_operator,
-                      estimate_snr, apply_inverse_cov, INVERSE_METHODS)
+                      apply_inverse_epochs, apply_inverse_tfr_epochs,
+                      write_inverse_operator, compute_rank_inverse,
+                      prepare_inverse_operator, estimate_snr,
+                      apply_inverse_cov, INVERSE_METHODS)
 from .time_frequency import (source_band_induced_power, source_induced_power,
                              compute_source_psd, compute_source_psd_epochs)
 from .resolution_matrix import (make_inverse_resolution_matrix,

mne/minimum_norm/inverse.py

@@ -27,7 +27,7 @@
 from ..io.pick import channel_type, pick_info, pick_types, pick_channels
 from ..cov import (compute_whitener, _read_cov, _write_cov, Covariance,
                    prepare_noise_cov)
-from ..epochs import BaseEpochs
+from ..epochs import BaseEpochs, EpochsArray
 from ..evoked import EvokedArray, Evoked
 from ..forward import (compute_depth_prior, _read_forward_meas_info,
                        is_fixed_orient, compute_orient_prior,
@@ -883,6 +883,8 @@ def apply_inverse(evoked, inverse_operator, lambda2=1. / 9., method="dSPM",
     --------
     apply_inverse_raw : Apply inverse operator to raw object.
     apply_inverse_epochs : Apply inverse operator to epochs object.
+    apply_inverse_tfr_epochs : Apply inverse operator to epochs tfr object.
+    apply_inverse_cov : Apply inverse operator to covariance object.
 
     Notes
     -----
@@ -1067,8 +1069,10 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, method="dSPM",
 
     See Also
     --------
-    apply_inverse_epochs : Apply inverse operator to epochs object.
     apply_inverse : Apply inverse operator to evoked object.
+    apply_inverse_epochs : Apply inverse operator to epochs object.
+    apply_inverse_tfr_epochs : Apply inverse operator to epochs tfr object.
+    apply_inverse_cov : Apply inverse operator to covariance object.
     """
     _validate_type(raw, BaseRaw, 'raw')
     _check_reference(raw, inverse_operator['info']['ch_names'])
@@ -1257,13 +1261,15 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, method="dSPM",
 
     Returns
     -------
-    stc : list of (SourceEstimate | VectorSourceEstimate | VolSourceEstimate)
+    stcs : list of (SourceEstimate | VectorSourceEstimate | VolSourceEstimate)
         The source estimates for all epochs.
 
     See Also
     --------
     apply_inverse_raw : Apply inverse operator to raw object.
     apply_inverse : Apply inverse operator to evoked object.
+    apply_inverse_tfr_epochs : Apply inverse operator to epochs tfr object.
+    apply_inverse_cov : Apply inverse operator to a covariance object.
     """
     stcs = _apply_inverse_epochs_gen(
         epochs, inverse_operator, lambda2, method=method, label=label,
@@ -1277,6 +1283,85 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, method="dSPM",
     return stcs
 
 
+def _apply_inverse_tfr_epochs_gen(epochs_tfr, inverse_operator, lambda2,
+                                  method, label, nave, pick_ori, prepared,
+                                  method_params, use_cps):
+    for freq_idx in range(epochs_tfr.freqs.size):
+        epochs = EpochsArray(epochs_tfr.data[:, :, freq_idx, :],
+                             epochs_tfr.info, events=epochs_tfr.events,
+                             tmin=epochs_tfr.tmin)
+        this_inverse_operator = inverse_operator[freq_idx] if \
+            isinstance(inverse_operator, (list, tuple)) else inverse_operator
+        stcs = _apply_inverse_epochs_gen(
+            epochs, this_inverse_operator, lambda2, method=method,
+            label=label, nave=nave, pick_ori=pick_ori, prepared=prepared,
+            method_params=method_params, use_cps=use_cps)
+        yield stcs
+
+
+@verbose
+def apply_inverse_tfr_epochs(epochs_tfr, inverse_operator, lambda2,
+                             method="dSPM", label=None, nave=1, pick_ori=None,
+                             return_generator=False, prepared=False,
+                             method_params=None, use_cps=True, verbose=None):
+    """Apply inverse operator to EpochsTFR.
+
+    Parameters
+    ----------
+    epochs_tfr : EpochsTFR object
+        Single trial, phase-amplitude (complex-valued), time-frequency epochs.
+    inverse_operator : list of dict | dict
+        The inverse operator for each frequency or a single inverse operator
+        to use for all frequencies.
+    lambda2 : float
+        The regularization parameter.
+    method : "MNE" | "dSPM" | "sLORETA" | "eLORETA"
+        Use minimum norm, dSPM (default), sLORETA, or eLORETA.
+    label : Label | None
+        Restricts the source estimates to a given label. If None,
+        source estimates will be computed for the entire source space.
+    nave : int
+        Number of averages used to regularize the solution.
+        Set to 1 on single Epoch by default.
+    %(pick_ori)s
+    return_generator : bool
+        Return a generator object instead of a list. This allows iterating
+        over the stcs without having to keep them all in memory.
+    prepared : bool
+        If True, do not call :func:`prepare_inverse_operator`.
+    method_params : dict | None
+        Additional options for eLORETA. See Notes of :func:`apply_inverse`.
+    %(use_cps_restricted)s
+    %(verbose)s
+
+    Returns
+    -------
+    stcs : list of list of (SourceEstimate | VectorSourceEstimate | VolSourceEstimate)
+        The source estimates for all frequencies (outside list) and for
+        all epochs (inside list).
+
+    See Also
+    --------
+    apply_inverse_raw : Apply inverse operator to raw object.
+    apply_inverse : Apply inverse operator to evoked object.
+    apply_inverse_epochs : Apply inverse operator to epochs object.
+    apply_inverse_cov : Apply inverse operator to a covariance object.
+    """  # noqa E501
+    from ..time_frequency.tfr import _check_tfr_complex
+    _check_tfr_complex(epochs_tfr)
+    if isinstance(inverse_operator, (list, tuple)) and \
+            len(inverse_operator) != epochs_tfr.freqs.size:
+        raise ValueError(f'Expected {epochs_tfr.freqs.size} inverse '
+                         f'operators, got {len(inverse_operator)}')
+    stcs = _apply_inverse_tfr_epochs_gen(
+        epochs_tfr, inverse_operator, lambda2,
+        method, label, nave, pick_ori, prepared,
+        method_params, use_cps)
+    if not return_generator:
+        stcs = [[stc for stc in tfr_stcs] for tfr_stcs in stcs]
+    return stcs
+
+
 @verbose
 def apply_inverse_cov(cov, info, inverse_operator, nave=1, lambda2=1 / 9,
                       method="dSPM", pick_ori=None, prepared=False,
@@ -1319,6 +1404,7 @@ def apply_inverse_cov(cov, info, inverse_operator, nave=1, lambda2=1 / 9,
     apply_inverse : Apply inverse operator to evoked object.
     apply_inverse_raw : Apply inverse operator to raw object.
     apply_inverse_epochs : Apply inverse operator to epochs object.
+    apply_inverse_tfr_epochs : Apply inverse operator to epochs tfr object.
 
     Notes
     -----

mne/minimum_norm/tests/test_inverse.py

@@ -25,12 +25,14 @@
                  convert_forward_solution, Covariance, combine_evoked,
                  SourceEstimate, make_sphere_model, make_ad_hoc_cov,
                  pick_channels_forward, compute_raw_covariance)
-from mne.io import read_raw_fif
+from mne.io import read_raw_fif, read_info
 from mne.minimum_norm import (apply_inverse, read_inverse_operator,
                               apply_inverse_raw, apply_inverse_epochs,
+                              apply_inverse_tfr_epochs,
                               make_inverse_operator, apply_inverse_cov,
                               write_inverse_operator, prepare_inverse_operator,
                               compute_rank_inverse, INVERSE_METHODS)
+from mne.time_frequency import EpochsTFR
 from mne.utils import catch_logging, _record_warnings
 
 test_path = testing.data_path(download=False)
@@ -1097,6 +1099,51 @@ def test_apply_mne_inverse_epochs():
             inverse_operator, 1.)
 
 
+@pytest.mark.slowtest
+@testing.requires_testing_data
+@pytest.mark.parametrize('return_generator', (True, False))
+def test_apply_inverse_tfr(return_generator):
+    """Test applying an inverse to time-frequency data."""
+    rng = np.random.default_rng(11)
+    n_epochs = 4
+    info = read_info(fname_raw)
+    inverse_operator = read_inverse_operator(fname_full)
+    freqs = np.arange(8, 10)
+    sfreq = info['sfreq']
+    times = np.arange(sfreq) / sfreq  # make epochs 1s long
+    data = rng.random((n_epochs, len(info.ch_names), freqs.size, times.size))
+    data = data + 1j * data  # make complex to simulate amplitude + phase
+    epochs_tfr = EpochsTFR(info, data, times=times, freqs=freqs)
+    epochs_tfr.apply_baseline((0, 0.5))
+    pick_ori = 'vector'
+
+    with pytest.raises(ValueError, match='Expected 2 inverse operators, '
+                                         'got 3'):
+        apply_inverse_tfr_epochs(epochs_tfr, [inverse_operator] * 3, lambda2)
+
+    # test epochs
+    stcs = apply_inverse_tfr_epochs(
+        epochs_tfr, inverse_operator, lambda2, "dSPM", pick_ori=pick_ori,
+        return_generator=return_generator)
+
+    n_orient = 3 if pick_ori == 'vector' else 1
+    if return_generator:
+        stcs = [[s for s in these_stcs] for these_stcs in stcs]
+    assert_allclose(stcs[0][0].times, times)
+    assert len(stcs) == freqs.size
+    assert all([len(s) == len(epochs_tfr) for s in stcs])
+    assert all([s.data.shape == (inverse_operator['nsource'],
+                                 n_orient, times.size)
+                for these_stcs in stcs for s in these_stcs])
+
+    evoked = EvokedArray(data.mean(axis=(0, 2)), info, epochs_tfr.tmin)
+    stc = apply_inverse(
+        evoked, inverse_operator, lambda2, "dSPM", pick_ori=pick_ori)
+    tfr_stc_data = np.array([[stc.data for stc in tfr_stcs]
+                             for tfr_stcs in stcs])
+    assert_allclose(stc.data, tfr_stc_data.mean(axis=(0, 1)))
+
+
 def test_make_inverse_operator_bads(evoked, noise_cov):
     """Test MNE inverse computation given a mismatch of bad channels."""
     fwd_op = read_forward_solution_meg(fname_fwd, surf_ori=True)

mne/time_frequency/psd.py

@@ -254,9 +254,10 @@ def psd_welch(inst, fmin=0, fmax=np.inf, tmin=None, tmax=None, n_fft=256,
     .. versionadded:: 0.12.0
     """
     spectrum = inst.compute_psd(
-        'welch', fmin, fmax, tmin, tmax, picks, proj, reject_by_annotation,
-        n_jobs=n_jobs, verbose=verbose, n_fft=n_fft, n_overlap=n_overlap,
-        n_per_seg=n_per_seg, average=average, window=window)
+        'welch', fmin=fmin, fmax=fmax, tmin=tmin, tmax=tmax, picks=picks,
+        proj=proj, n_jobs=n_jobs, verbose=verbose, n_fft=n_fft,
+        n_overlap=n_overlap, n_per_seg=n_per_seg, average=average,
+        window=window)
     return spectrum.get_data(return_freqs=True)
 
 
@@ -323,8 +324,8 @@ def psd_multitaper(inst, fmin=0, fmax=np.inf, tmin=None, tmax=None,
     .. footbibliography::
     """
     spectrum = inst.compute_psd(
-        'multitaper', fmin, fmax, tmin, tmax, picks, proj,
-        reject_by_annotation, n_jobs=n_jobs, verbose=verbose,
-        bandwidth=bandwidth, adaptive=adaptive, low_bias=low_bias,
-        normalization=normalization)
+        'multitaper', fmin=fmin, fmax=fmax, tmin=tmin, tmax=tmax, picks=picks,
+        proj=proj, reject_by_annotation=reject_by_annotation, n_jobs=n_jobs,
+        verbose=verbose, bandwidth=bandwidth, adaptive=adaptive,
+        low_bias=low_bias, normalization=normalization)
     return spectrum.get_data(return_freqs=True)

mne/time_frequency/tfr.py

@@ -2658,3 +2658,9 @@ def _preproc_tfr_instance(tfr, picks, tmin, tmax, fmin, fmax, vmin, vmax, dB,
     tfr.data = data
 
     return tfr
+
+
+def _check_tfr_complex(tfr, reason='source space estimation'):
+    """Check that time-frequency epochs or average data is complex."""
+    if not np.iscomplexobj(tfr.data):
+        raise RuntimeError(f'Time-frequency data must be complex for {reason}')

mne/utils/__init__.py

@@ -29,7 +29,7 @@
                      sys_info, _get_extra_data_path, _get_root_dir,
                      _get_numpy_libs)
 from .docs import (copy_function_doc_to_method_doc, copy_doc, linkcode_resolve,
-                   open_docs, deprecated, fill_doc, deprecated_alias,
+                   open_docs, deprecated, fill_doc, deprecated_alias, legacy,
                    copy_base_doc_to_subclass_doc, docdict as _docdict)
 from .fetching import _url_to_local_path
 from ._logging import (verbose, logger, set_log_level, set_log_file,