test_source_estimate.py

#
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import os
import re
from contextlib import nullcontext
from copy import deepcopy
from pathlib import Path
from shutil import copyfile

import numpy as np
import pytest
from numpy.fft import fft
from numpy.testing import (
    assert_allclose,
    assert_array_almost_equal,
    assert_array_equal,
    assert_array_less,
    assert_equal,
)
from scipy import sparse
from scipy.optimize import fmin_cobyla
from scipy.spatial.distance import cdist

import mne
from mne import (
    Epochs,
    EvokedArray,
    Label,
    MixedSourceEstimate,
    MixedVectorSourceEstimate,
    SourceEstimate,
    SourceSpaces,
    VectorSourceEstimate,
    VolSourceEstimate,
    VolVectorSourceEstimate,
    compute_source_morph,
    convert_forward_solution,
    extract_label_time_course,
    find_events,
    labels_to_stc,
    pick_info,
    pick_types,
    pick_types_forward,
    read_cov,
    read_evokeds,
    read_forward_solution,
    read_source_estimate,
    read_source_spaces,
    read_trans,
    scale_mri,
    setup_volume_source_space,
    spatial_inter_hemi_adjacency,
    spatial_src_adjacency,
    spatial_tris_adjacency,
    spatio_temporal_src_adjacency,
    spatio_temporal_tris_adjacency,
    stats,
    stc_near_sensors,
    write_source_spaces,
)
from mne._fiff.constants import FIFF
from mne.datasets import testing
from mne.fixes import _get_img_fdata
from mne.io import read_info, read_raw_fif
from mne.label import label_sign_flip, read_labels_from_annot
from mne.minimum_norm import (
    apply_inverse,
    apply_inverse_epochs,
    make_inverse_operator,
    read_inverse_operator,
)
from mne.morph_map import _make_morph_map_hemi
from mne.source_estimate import _get_vol_mask, _make_stc, grade_to_tris
from mne.source_space._source_space import _get_src_nn
from mne.transforms import apply_trans, invert_transform
from mne.utils import (
    _record_warnings,
    catch_logging,
)

data_path = testing.data_path(download=False)
subjects_dir = data_path / "subjects"
fname_inv = (
    data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif"
)
fname_inv_fixed = (
    data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-4-meg-fixed-inv.fif"
)
fname_fwd = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-4-fwd.fif"
fname_cov = data_path / "MEG" / "sample" / "sample_audvis_trunc-cov.fif"
fname_evoked = data_path / "MEG" / "sample" / "sample_audvis_trunc-ave.fif"
fname_raw = data_path / "MEG" / "sample" / "sample_audvis_trunc_raw.fif"
fname_t1 = data_path / "subjects" / "sample" / "mri" / "T1.mgz"
fname_fs_t1 = data_path / "subjects" / "fsaverage" / "mri" / "T1.mgz"
fname_aseg = data_path / "subjects" / "sample" / "mri" / "aseg.mgz"
fname_src = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-6-fwd.fif"
fname_src_fs = data_path / "subjects" / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
bem_path = data_path / "subjects" / "sample" / "bem"
fname_src_3 = bem_path / "sample-oct-4-src.fif"
fname_src_vol = bem_path / "sample-volume-7mm-src.fif"
fname_stc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg"
fname_vol = (
    data_path / "MEG" / "sample" / "sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w"
)
fname_vsrc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-vol-7-fwd.fif"
fname_inv_vol = (
    data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-vol-7-meg-inv.fif"
)
fname_nirx = data_path / "NIRx" / "nirscout" / "nirx_15_0_recording"
rng = np.random.RandomState(0)

pytest.importorskip("nibabel")


@testing.requires_testing_data
def test_stc_baseline_correction():
    """Test baseline correction for source estimate objects."""
    # test on different source estimates
    stcs = [read_source_estimate(fname_stc), read_source_estimate(fname_vol, "sample")]
    # test on different "baseline" intervals
    baselines = [(0.0, 0.1), (None, None)]

    for stc in stcs:
        times = stc.times

        for start, stop in baselines:
            # apply baseline correction, then check if it worked
            stc = stc.apply_baseline(baseline=(start, stop))

            t0 = start or stc.times[0]
            t1 = stop or stc.times[-1]
            # index for baseline interval (include boundary latencies)
            imin = np.abs(times - t0).argmin()
            imax = np.abs(times - t1).argmin() + 1
            # data matrix from baseline interval
            data_base = stc.data[:, imin:imax]
            mean_base = data_base.mean(axis=1)
            zero_array = np.zeros(mean_base.shape[0])
            # test if baseline properly subtracted (mean=zero for all sources)
            assert_array_almost_equal(mean_base, zero_array)


@testing.requires_testing_data
def test_spatial_inter_hemi_adjacency():
    """Test spatial adjacency between hemispheres."""
    # trivial cases
    conn = spatial_inter_hemi_adjacency(fname_src_3, 5e-6)
    assert_equal(conn.data.size, 0)
    conn = spatial_inter_hemi_adjacency(fname_src_3, 5e6)
    assert_equal(conn.data.size, np.prod(conn.shape) // 2)
    # actually interesting case (1cm), should be between 2 and 10% of verts
    src = read_source_spaces(fname_src_3)
    conn = spatial_inter_hemi_adjacency(src, 10e-3)
    conn = conn.tocsr()
    n_src = conn.shape[0]
    assert n_src * 0.02 < conn.data.size < n_src * 0.10
    assert_equal(conn[: src[0]["nuse"], : src[0]["nuse"]].data.size, 0)
    assert_equal(conn[-src[1]["nuse"] :, -src[1]["nuse"] :].data.size, 0)
    c = (conn.T + conn) / 2.0 - conn
    c.eliminate_zeros()
    assert_equal(c.data.size, 0)
    # check locations
    upper_right = conn[: src[0]["nuse"], src[0]["nuse"] :].toarray()
    assert_equal(upper_right.sum(), conn.sum() // 2)
    good_labels = ["S_pericallosal", "Unknown", "G_and_S_cingul-Mid-Post", "G_cuneus"]
    for hi, hemi in enumerate(("lh", "rh")):
        has_neighbors = src[hi]["vertno"][np.where(np.any(upper_right, axis=1 - hi))[0]]
        labels = read_labels_from_annot(
            "sample", "aparc.a2009s", hemi, subjects_dir=subjects_dir
        )
        use_labels = [
            label.name[:-3]
            for label in labels
            if np.isin(label.vertices, has_neighbors).any()
        ]
        assert set(use_labels) - set(good_labels) == set()


@pytest.mark.slowtest
@testing.requires_testing_data
def test_volume_stc(tmp_path):
    """Test volume STCs."""
    h5io = pytest.importorskip("h5io")

    N = 100
    data = np.arange(N)[:, np.newaxis]
    datas = [data, data, np.arange(2)[:, np.newaxis], np.arange(6).reshape(2, 3, 1)]
    vertno = np.arange(N)
    vertnos = [vertno, vertno[:, np.newaxis], np.arange(2)[:, np.newaxis], np.arange(2)]
    vertno_reads = [vertno, vertno, np.arange(2), np.arange(2)]
    for data, vertno, vertno_read in zip(datas, vertnos, vertno_reads):
        if data.ndim in (1, 2):
            stc = VolSourceEstimate(data, [vertno], 0, 1)
            ext = "stc"
            klass = VolSourceEstimate
        else:
            assert data.ndim == 3
            stc = VolVectorSourceEstimate(data, [vertno], 0, 1)
            ext = "h5"
            klass = VolVectorSourceEstimate
        fname_temp = tmp_path / ("temp-vl." + ext)
        stc_new = stc
        n = 3 if ext == "h5" else 2
        for ii in range(n):
            if ii < 2:
                stc_new.save(fname_temp, overwrite=True)
            else:
                # Pass stc.vertices[0], an ndarray, to ensure support for
                # the way we used to write volume STCs
                h5io.write_hdf5(
                    str(fname_temp),
                    dict(
                        vertices=stc.vertices[0],
                        data=stc.data,
                        tmin=stc.tmin,
                        tstep=stc.tstep,
                        subject=stc.subject,
                        src_type=stc._src_type,
                    ),
                    title="mnepython",
                    overwrite=True,
                )
            stc_new = read_source_estimate(fname_temp)
            assert isinstance(stc_new, klass)
            assert_array_equal(vertno_read, stc_new.vertices[0])
            assert_array_almost_equal(stc.data, stc_new.data)

    # now let's actually read a MNE-C processed file
    stc = read_source_estimate(fname_vol, "sample")
    assert isinstance(stc, VolSourceEstimate)

    assert "sample" in repr(stc)
    assert " KiB" in repr(stc)

    stc_new = stc
    fname_temp = tmp_path / ("temp-vl.stc")
    with pytest.raises(ValueError, match="'ftype' parameter"):
        stc.save(fname_vol, ftype="whatever", overwrite=True)
    for ftype in ["w", "h5"]:
        for _ in range(2):
            fname_temp = tmp_path / f"temp-vol.{ftype}"
            stc_new.save(fname_temp, ftype=ftype, overwrite=True)
            stc_new = read_source_estimate(fname_temp)
            assert isinstance(stc_new, VolSourceEstimate)
            assert_array_equal(stc.vertices[0], stc_new.vertices[0])
            assert_array_almost_equal(stc.data, stc_new.data)


@testing.requires_testing_data
def test_save_stc_as_gifti(tmp_path):
    """Save the stc as a GIFTI file and export."""
    nib = pytest.importorskip("nibabel")
    surfpath_src = bem_path / "sample-oct-6-src.fif"
    surfpath_stc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg"
    src = read_source_spaces(surfpath_src)  # need source space
    stc = read_source_estimate(surfpath_stc)  # need stc
    assert isinstance(src, SourceSpaces)
    assert isinstance(stc, SourceEstimate)

    surf_fname = tmp_path / "stc_write"

    stc.save_as_surface(surf_fname, src)

    # did structural get written?
    img_lh = nib.load(f"{surf_fname}-lh.gii")
    img_rh = nib.load(f"{surf_fname}-rh.gii")
    assert isinstance(img_lh, nib.gifti.gifti.GiftiImage)
    assert isinstance(img_rh, nib.gifti.gifti.GiftiImage)

    # did time series get written?
    img_timelh = nib.load(f"{surf_fname}-lh.time.gii")
    img_timerh = nib.load(f"{surf_fname}-rh.time.gii")
    assert isinstance(img_timelh, nib.gifti.gifti.GiftiImage)
    assert isinstance(img_timerh, nib.gifti.gifti.GiftiImage)


@testing.requires_testing_data
def test_stc_as_volume():
    """Test previous volume source estimate morph."""
    nib = pytest.importorskip("nibabel")
    inverse_operator_vol = read_inverse_operator(fname_inv_vol)

    # Apply inverse operator
    stc_vol = read_source_estimate(fname_vol, "sample")

    img = stc_vol.as_volume(inverse_operator_vol["src"], mri_resolution=True, dest="42")
    t1_img = nib.load(fname_t1)
    # always assure nifti and dimensionality
    assert isinstance(img, nib.Nifti1Image)
    assert img.header.get_zooms()[:3] == t1_img.header.get_zooms()[:3]

    img = stc_vol.as_volume(inverse_operator_vol["src"], mri_resolution=False)

    assert isinstance(img, nib.Nifti1Image)
    assert img.shape[:3] == inverse_operator_vol["src"][0]["shape"][:3]

    with pytest.raises(ValueError, match="Invalid value.*output.*"):
        stc_vol.as_volume(inverse_operator_vol["src"], format="42")


@testing.requires_testing_data
def test_save_vol_stc_as_nifti(tmp_path):
    """Save the stc as a nifti file and export."""
    nib = pytest.importorskip("nibabel")
    src = read_source_spaces(fname_vsrc)
    vol_fname = tmp_path / "stc.nii.gz"

    # now let's actually read a MNE-C processed file
    stc = read_source_estimate(fname_vol, "sample")
    assert isinstance(stc, VolSourceEstimate)

    stc.save_as_volume(vol_fname, src, dest="surf", mri_resolution=False)
    with _record_warnings():  # nib<->numpy
        img = nib.load(str(vol_fname))
    assert img.shape == src[0]["shape"] + (len(stc.times),)

    with _record_warnings():  # nib<->numpy
        t1_img = nib.load(fname_t1)
    stc.save_as_volume(vol_fname, src, dest="mri", mri_resolution=True, overwrite=True)
    with _record_warnings():  # nib<->numpy
        img = nib.load(str(vol_fname))
    assert img.shape == t1_img.shape + (len(stc.times),)
    assert_allclose(img.affine, t1_img.affine, atol=1e-5)

    # export without saving
    img = stc.as_volume(src, dest="mri", mri_resolution=True)
    assert img.shape == t1_img.shape + (len(stc.times),)
    assert_allclose(img.affine, t1_img.affine, atol=1e-5)

    src = SourceSpaces([src[0], src[0]])
    stc = VolSourceEstimate(
        np.r_[stc.data, stc.data],
        [stc.vertices[0], stc.vertices[0]],
        tmin=stc.tmin,
        tstep=stc.tstep,
        subject="sample",
    )
    img = stc.as_volume(src, dest="mri", mri_resolution=False)
    assert img.shape == src[0]["shape"] + (len(stc.times),)


@testing.requires_testing_data
def test_expand():
    """Test stc expansion."""
    stc_ = read_source_estimate(fname_stc, "sample")
    vec_stc_ = VectorSourceEstimate(
        np.zeros((stc_.data.shape[0], 3, stc_.data.shape[1])),
        stc_.vertices,
        stc_.tmin,
        stc_.tstep,
        stc_.subject,
    )

    for stc in [stc_, vec_stc_]:
        assert "sample" in repr(stc)
        labels_lh = read_labels_from_annot(
            "sample", "aparc", "lh", subjects_dir=subjects_dir
        )
        new_label = labels_lh[0] + labels_lh[1]
        stc_limited = stc.in_label(new_label)
        stc_new = stc_limited.copy()
        stc_new.data.fill(0)
        for label in labels_lh[:2]:
            stc_new += stc.in_label(label).expand(stc_limited.vertices)
        pytest.raises(TypeError, stc_new.expand, stc_limited.vertices[0])
        pytest.raises(ValueError, stc_new.expand, [stc_limited.vertices[0]])
        # make sure we can't add unless vertno agree
        pytest.raises(ValueError, stc.__add__, stc.in_label(labels_lh[0]))


def _fake_stc(n_time=10, is_complex=False):
    np.random.seed(7)
    verts = [np.arange(10), np.arange(90)]
    data = np.random.rand(100, n_time)
    if is_complex:
        data.astype(complex)
    return SourceEstimate(data, verts, 0, 1e-1, "foo")


def _fake_vec_stc(n_time=10, is_complex=False):
    np.random.seed(7)
    verts = [np.arange(10), np.arange(90)]
    data = np.random.rand(100, 3, n_time)
    if is_complex:
        data.astype(complex)
    return VectorSourceEstimate(data, verts, 0, 1e-1, "foo")


@testing.requires_testing_data
def test_stc_snr():
    """Test computing SNR from a STC."""
    inv = read_inverse_operator(fname_inv_fixed)
    fwd = read_forward_solution(fname_fwd)
    cov = read_cov(fname_cov)
    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0].crop(0, 0.01)
    stc = apply_inverse(evoked, inv)
    assert (stc.data < 0).any()
    with pytest.warns(RuntimeWarning, match="nAm"):
        stc.estimate_snr(evoked.info, fwd, cov)  # dSPM
    with _record_warnings(), pytest.warns(RuntimeWarning, match="free ori"):
        abs(stc).estimate_snr(evoked.info, fwd, cov)
    stc = apply_inverse(evoked, inv, method="MNE")
    snr = stc.estimate_snr(evoked.info, fwd, cov)
    assert_allclose(snr.times, evoked.times)
    snr = snr.data
    assert snr.max() < -10
    assert snr.min() > -120


def test_stc_attributes():
    """Test STC attributes."""
    stc = _fake_stc(n_time=10)
    vec_stc = _fake_vec_stc(n_time=10)

    n_times = len(stc.times)
    assert_equal(stc._data.shape[-1], n_times)
    assert_array_equal(stc.times, stc.tmin + np.arange(n_times) * stc.tstep)
    assert_array_almost_equal(
        stc.times, [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
    )

    def attempt_times_mutation(stc):
        stc.times -= 1

    def attempt_assignment(stc, attr, val):
        setattr(stc, attr, val)

    # .times is read-only
    pytest.raises(ValueError, attempt_times_mutation, stc)
    pytest.raises(ValueError, attempt_assignment, stc, "times", [1])

    # Changing .tmin or .tstep re-computes .times
    stc.tmin = 1
    assert isinstance(stc.tmin, float)
    assert_array_almost_equal(
        stc.times, [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9]
    )

    stc.tstep = 1
    assert isinstance(stc.tstep, float)
    assert_array_almost_equal(
        stc.times, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    )

    # tstep <= 0 is not allowed
    pytest.raises(ValueError, attempt_assignment, stc, "tstep", 0)
    pytest.raises(ValueError, attempt_assignment, stc, "tstep", -1)

    # Changing .data re-computes .times
    stc.data = np.random.rand(100, 5)
    assert_array_almost_equal(stc.times, [1.0, 2.0, 3.0, 4.0, 5.0])

    # .data must match the number of vertices
    pytest.raises(ValueError, attempt_assignment, stc, "data", [[1]])
    pytest.raises(ValueError, attempt_assignment, stc, "data", None)

    # .data much match number of dimensions
    pytest.raises(ValueError, attempt_assignment, stc, "data", np.arange(100))
    pytest.raises(ValueError, attempt_assignment, vec_stc, "data", [np.arange(100)])
    pytest.raises(ValueError, attempt_assignment, vec_stc, "data", [[[np.arange(100)]]])

    # .shape attribute must also work when ._data is None
    stc._kernel = np.zeros((2, 2))
    stc._sens_data = np.zeros((2, 3))
    stc._data = None
    assert_equal(stc.shape, (2, 3))

    # bad size of data
    stc = _fake_stc()
    data = stc.data[:, np.newaxis, :]
    with pytest.raises(ValueError, match="2 dimensions for SourceEstimate"):
        SourceEstimate(data, stc.vertices, 0, 1)
    stc = SourceEstimate(data[:, 0, 0], stc.vertices, 0, 1)
    assert stc.data.shape == (len(data), 1)


def test_io_stc(tmp_path):
    """Test IO for STC files."""
    stc = _fake_stc()
    stc.save(tmp_path / "tmp.stc")
    stc2 = read_source_estimate(tmp_path / "tmp.stc")

    assert_array_almost_equal(stc.data, stc2.data)
    assert_array_almost_equal(stc.tmin, stc2.tmin)
    assert_equal(len(stc.vertices), len(stc2.vertices))
    for v1, v2 in zip(stc.vertices, stc2.vertices):
        assert_array_almost_equal(v1, v2)
    assert_array_almost_equal(stc.tstep, stc2.tstep)
    # test warning for complex data
    stc2.data = stc2.data.astype(np.complex128)
    with pytest.raises(ValueError, match="Cannot save complex-valued STC"):
        stc2.save(tmp_path / "complex.stc")


@pytest.mark.parametrize("is_complex", (True, False))
@pytest.mark.parametrize("vector", (True, False))
def test_io_stc_h5(tmp_path, is_complex, vector):
    """Test IO for STC files using HDF5."""
    pytest.importorskip("h5io")
    if vector:
        stc = _fake_vec_stc(is_complex=is_complex)
    else:
        stc = _fake_stc(is_complex=is_complex)
    match = "can only be written" if vector else "Invalid value for the 'ftype"
    with pytest.raises(ValueError, match=match):
        stc.save(tmp_path / "tmp.h5", ftype="foo")
    out_name = tmp_path / "tmp"
    stc.save(out_name, ftype="h5")
    # test overwrite
    assert out_name.with_name(out_name.name + "-stc.h5").is_file()
    with pytest.raises(FileExistsError, match="Destination file exists"):
        stc.save(out_name, ftype="h5")
    stc.save(out_name, ftype="h5", overwrite=True)
    stc3 = read_source_estimate(out_name)
    stc4 = read_source_estimate(out_name.with_name(out_name.name + "-stc"))
    stc5 = read_source_estimate(out_name.with_name(out_name.name + "-stc.h5"))
    pytest.raises(RuntimeError, read_source_estimate, out_name, subject="bar")
    for stc_new in stc3, stc4, stc5:
        assert_equal(stc_new.subject, stc.subject)
        assert_array_equal(stc_new.data, stc.data)
        assert_array_equal(stc_new.tmin, stc.tmin)
        assert_array_equal(stc_new.tstep, stc.tstep)
        assert_equal(len(stc_new.vertices), len(stc.vertices))
        for v1, v2 in zip(stc_new.vertices, stc.vertices):
            assert_array_equal(v1, v2)


def test_io_w(tmp_path):
    """Test IO for w files."""
    stc = _fake_stc(n_time=1)
    w_fname = tmp_path / "fake"
    stc.save(w_fname, ftype="w")
    src = read_source_estimate(w_fname)
    src.save(tmp_path / "tmp", ftype="w")
    src2 = read_source_estimate(tmp_path / "tmp-lh.w")
    assert_array_almost_equal(src.data, src2.data)
    assert_array_almost_equal(src.lh_vertno, src2.lh_vertno)
    assert_array_almost_equal(src.rh_vertno, src2.rh_vertno)


def test_stc_arithmetic():
    """Test arithmetic for STC files."""
    stc = _fake_stc()
    data = stc.data.copy()
    vec_stc = _fake_vec_stc()
    vec_data = vec_stc.data.copy()

    out = list()
    for a in [data, stc, vec_data, vec_stc]:
        a = a + a * 3 + 3 * a - a**2 / 2

        a += a
        a -= a
        with np.errstate(invalid="ignore"):
            a /= 2 * a
        a *= -a

        a += 2
        a -= 1
        a *= -1
        a /= 2
        b = 2 + a
        b = 2 - a
        b = +a
        assert_array_equal(b.data, a.data)
        with np.errstate(invalid="ignore"):
            a **= 3
        out.append(a)

    assert_array_equal(out[0], out[1].data)
    assert_array_equal(out[2], out[3].data)
    assert_array_equal(stc.sqrt().data, np.sqrt(stc.data))
    assert_array_equal(vec_stc.sqrt().data, np.sqrt(vec_stc.data))
    assert_array_equal(abs(stc).data, abs(stc.data))
    assert_array_equal(abs(vec_stc).data, abs(vec_stc.data))

    stc_sum = stc.sum()
    assert_array_equal(stc_sum.data, stc.data.sum(1, keepdims=True))
    stc_mean = stc.mean()
    assert_array_equal(stc_mean.data, stc.data.mean(1, keepdims=True))
    vec_stc_mean = vec_stc.mean()
    assert_array_equal(vec_stc_mean.data, vec_stc.data.mean(2, keepdims=True))


@pytest.mark.slowtest
@testing.requires_testing_data
@pytest.mark.parametrize("kind", ("scalar", "vector"))
@pytest.mark.parametrize("method", ("fft", "polyphase"))
def test_stc_methods(kind, method):
    """Test stc methods lh_data, rh_data, bin(), resample()."""
    stc = read_source_estimate(fname_stc)

    if kind == "vector":
        # Make a vector version of the above source estimate
        x = stc.data[:, np.newaxis, :]
        yz = np.zeros((x.shape[0], 2, x.shape[2]))
        stc = VectorSourceEstimate(
            np.concatenate((x, yz), 1),
            stc.vertices,
            stc.tmin,
            stc.tstep,
            stc.subject,
        )

    # lh_data / rh_data
    assert_array_equal(stc.lh_data, stc.data[: len(stc.lh_vertno)])
    assert_array_equal(stc.rh_data, stc.data[len(stc.lh_vertno) :])

    # bin
    binned = stc.bin(0.12)
    a = np.mean(stc.data[..., : np.searchsorted(stc.times, 0.12)], axis=-1)
    assert_array_equal(a, binned.data[..., 0])

    stc = read_source_estimate(fname_stc)
    stc.subject = "sample"
    label_lh = read_labels_from_annot(
        "sample", "aparc", "lh", subjects_dir=subjects_dir
    )[0]
    label_rh = read_labels_from_annot(
        "sample", "aparc", "rh", subjects_dir=subjects_dir
    )[0]
    label_both = label_lh + label_rh
    for label in (label_lh, label_rh, label_both):
        assert isinstance(stc.shape, tuple) and len(stc.shape) == 2
        stc_label = stc.in_label(label)
        if label.hemi != "both":
            if label.hemi == "lh":
                verts = stc_label.vertices[0]
            else:  # label.hemi == 'rh':
                verts = stc_label.vertices[1]
            n_vertices_used = len(label.get_vertices_used(verts))
            assert_equal(len(stc_label.data), n_vertices_used)
    stc_lh = stc.in_label(label_lh)
    pytest.raises(ValueError, stc_lh.in_label, label_rh)
    label_lh.subject = "foo"
    pytest.raises(RuntimeError, stc.in_label, label_lh)

    stc_new = deepcopy(stc)
    o_sfreq = 1.0 / stc.tstep
    # note that using no padding for this STC reduces edge ringing...
    stc_new.resample(2 * o_sfreq, npad=0, method=method)
    assert stc_new.data.shape[1] == 2 * stc.data.shape[1]
    assert stc_new.tstep == stc.tstep / 2
    stc_new.resample(o_sfreq, npad=0, method=method)
    assert stc_new.data.shape[1] == stc.data.shape[1]
    assert stc_new.tstep == stc.tstep
    if method == "fft":
        # no low-passing so survives round-trip
        assert_allclose(stc_new.data, stc.data, atol=1e-5)
    else:
        # low-passing means we need something more flexible
        corr = np.corrcoef(stc_new.data.ravel(), stc.data.ravel())[0, 1]
        assert 0.99 < corr < 1


@testing.requires_testing_data
def test_stc_resamp_noop():
    """Tests resampling doesn't affect data if sfreq is identical."""
    stc = read_source_estimate(fname_stc)
    data_before = stc.data
    data_after = stc.resample(sfreq=1.0 / stc.tstep).data
    assert_array_equal(data_before, data_after)


@testing.requires_testing_data
def test_center_of_mass():
    """Test computing the center of mass on an stc."""
    stc = read_source_estimate(fname_stc)
    pytest.raises(ValueError, stc.center_of_mass, "sample")
    stc.lh_data[:] = 0
    vertex, hemi, t = stc.center_of_mass("sample", subjects_dir=subjects_dir)
    assert hemi == 1
    # XXX Should design a fool-proof test case, but here were the
    # results:
    assert_equal(vertex, 124791)
    assert_equal(np.round(t, 2), 0.12)


@testing.requires_testing_data
@pytest.mark.parametrize("kind", ("surface", "mixed"))
@pytest.mark.parametrize("vector", (False, True))
def test_extract_label_time_course(kind, vector):
    """Test extraction of label time courses from (Mixed)SourceEstimate."""
    n_stcs = 3
    n_times = 50

    src = read_inverse_operator(fname_inv)["src"]
    if kind == "mixed":
        pytest.importorskip("nibabel")
        label_names = ("Left-Cerebellum-Cortex", "Right-Cerebellum-Cortex")
        src += setup_volume_source_space(
            "sample",
            pos=20.0,
            volume_label=label_names,
            subjects_dir=subjects_dir,
            add_interpolator=False,
        )
        klass = MixedVectorSourceEstimate
    else:
        klass = VectorSourceEstimate
    if not vector:
        klass = klass._scalar_class
    vertices = [s["vertno"] for s in src]
    n_verts = np.array([len(v) for v in vertices])
    vol_means = np.arange(-1, 1 - len(src), -1)
    vol_means_t = np.repeat(vol_means[:, np.newaxis], n_times, axis=1)

    # get some labels
    labels_lh = read_labels_from_annot("sample", hemi="lh", subjects_dir=subjects_dir)
    labels_rh = read_labels_from_annot("sample", hemi="rh", subjects_dir=subjects_dir)
    labels = list()
    labels.extend(labels_lh[:5])
    labels.extend(labels_rh[:4])

    n_labels = len(labels)

    label_tcs = dict(mean=np.arange(n_labels)[:, None] * np.ones((n_labels, n_times)))
    label_tcs["max"] = label_tcs["mean"]
    label_tcs[None] = label_tcs["mean"]

    # compute the mean with sign flip
    label_tcs["mean_flip"] = np.zeros_like(label_tcs["mean"])
    for i, label in enumerate(labels):
        label_tcs["mean_flip"][i] = i * np.mean(label_sign_flip(label, src[:2]))

    # compute pca_flip
    label_flip = []
    for i, label in enumerate(labels):
        this_flip = i * label_sign_flip(label, src[:2])
        label_flip.append(this_flip)
    # compute pca_flip
    label_tcs["pca_flip"] = np.zeros_like(label_tcs["mean"])
    for i, (label, flip) in enumerate(zip(labels, label_flip)):
        sign = np.sign(np.dot(np.full((flip.shape[0]), i), flip))
        label_tcs["pca_flip"][i] = sign * label_tcs["mean"][i]

    # generate some stc's with known data
    stcs = list()
    pad = (((0, 0), (2, 0), (0, 0)), "constant")
    for i in range(n_stcs):
        data = np.zeros((n_verts.sum(), n_times))
        # set the value of the stc within each label
        for j, label in enumerate(labels):
            if label.hemi == "lh":
                idx = np.intersect1d(vertices[0], label.vertices)
                idx = np.searchsorted(vertices[0], idx)
            elif label.hemi == "rh":
                idx = np.intersect1d(vertices[1], label.vertices)
                idx = len(vertices[0]) + np.searchsorted(vertices[1], idx)
            data[idx] = label_tcs["mean"][j]
        for j in range(len(vol_means)):
            offset = n_verts[: 2 + j].sum()
            data[offset : offset + n_verts[j]] = vol_means[j]

        if vector:
            # the values it on the Z axis
            data = np.pad(data[:, np.newaxis], *pad)
        this_stc = klass(data, vertices, 0, 1)
        stcs.append(this_stc)

    if vector:
        for key in label_tcs:
            label_tcs[key] = np.pad(label_tcs[key][:, np.newaxis], *pad)
        vol_means_t = np.pad(vol_means_t[:, np.newaxis], *pad)

    # test some invalid inputs
    with pytest.raises(ValueError, match="Invalid value for the 'mode'"):
        extract_label_time_course(stcs, labels, src, mode="notamode")

    # have an empty label
    empty_label = labels[0].copy()
    empty_label.vertices += 1000000
    with pytest.raises(ValueError, match="does not contain any vertices"):
        extract_label_time_course(stcs, empty_label, src)

    # but this works:
    with pytest.warns(RuntimeWarning, match="does not contain any vertices"):
        tc = extract_label_time_course(stcs, empty_label, src, allow_empty=True)
    end_shape = (3, n_times) if vector else (n_times,)
    for arr in tc:
        assert arr.shape == (1 + len(vol_means),) + end_shape
        assert_array_equal(arr[:1], np.zeros((1,) + end_shape))
        if len(vol_means):
            assert_array_equal(arr[1:], vol_means_t)

    # test the different modes
    modes = ["mean", "mean_flip", "pca_flip", "max", "auto", None]

    for mode in modes:
        if vector and mode not in ("mean", "max", "auto"):
            with pytest.raises(ValueError, match="when using a vector"):
                extract_label_time_course(stcs, labels, src, mode=mode)
            continue
        with _record_warnings():  # SVD convergence on arm64
            label_tc = extract_label_time_course(stcs, labels, src, mode=mode)
        label_tc_method = [
            stc.extract_label_time_course(labels, src, mode=mode) for stc in stcs
        ]
        assert len(label_tc) == n_stcs
        assert len(label_tc_method) == n_stcs
        for j, (tc1, tc2) in enumerate(zip(label_tc, label_tc_method)):
            if mode is None:
                assert all(arr.shape[1] == tc1[0].shape[1] for arr in tc1)
                assert all(arr.shape[1] == tc2[0].shape[1] for arr in tc2)
                assert (len(tc1), tc1[0].shape[1]) == (n_labels,) + end_shape
                assert (len(tc2), tc2[0].shape[1]) == (n_labels,) + end_shape
                for arr1, arr2 in zip(tc1, tc2):  # list of arrays
                    assert_allclose(arr1, arr2, rtol=1e-8, atol=1e-16)
            else:
                assert tc1.shape == (n_labels + len(vol_means),) + end_shape
                assert tc2.shape == (n_labels + len(vol_means),) + end_shape
                assert_allclose(tc1, tc2, rtol=1e-8, atol=1e-16)
            if mode == "auto":
                use_mode = "mean" if vector else "mean_flip"
            else:
                use_mode = mode
            if mode == "pca_flip":
                for arr1, arr2 in zip(tc1, label_tcs[use_mode]):
                    assert_array_almost_equal(arr1, arr2)
            elif use_mode is None:
                for arr1, arr2 in zip(
                    tc1[:n_labels], label_tcs[use_mode]
                ):  # list of arrays
                    assert_allclose(
                        arr1, np.tile(arr2, (arr1.shape[0], 1)), rtol=1e-8, atol=1e-16
                    )
            elif use_mode in ("mean", "max", "mean_flip"):
                assert_array_almost_equal(tc1[:n_labels], label_tcs[use_mode])
            if mode is not None:
                assert_array_almost_equal(tc1[n_labels:], vol_means_t)

    # test label with very few vertices (check SVD conditionals)
    label = Label(vertices=src[0]["vertno"][:2], hemi="lh")
    x = label_sign_flip(label, src[:2])
    assert len(x) == 2
    label = Label(vertices=[], hemi="lh")
    x = label_sign_flip(label, src[:2])
    assert x.size == 0


@testing.requires_testing_data
@pytest.mark.parametrize(
    "label_type, mri_res, vector, test_label, cf, call",
    [
        (str, False, False, False, "head", "meth"),  # head frame
        (str, False, False, str, "mri", "func"),  # fastest, default for testing
        (str, False, True, int, "mri", "func"),  # vector
        (str, True, False, False, "mri", "func"),  # mri_resolution
        (list, True, False, False, "mri", "func"),  # volume label as list
        (dict, True, False, False, "mri", "func"),  # volume label as dict
    ],
)
def test_extract_label_time_course_volume(
    src_volume_labels, label_type, mri_res, vector, test_label, cf, call
):
    """Test extraction of label time courses from Vol(Vector)SourceEstimate."""
    src_labels, volume_labels, lut = src_volume_labels
    n_tot = 46
    assert n_tot == len(src_labels)
    inv = read_inverse_operator(fname_inv_vol)
    if cf == "head":
        src = inv["src"]
        assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD
        rr = apply_trans(invert_transform(inv["mri_head_t"]), src[0]["rr"])
    else:
        assert cf == "mri"
        src = read_source_spaces(fname_src_vol)
        assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_MRI
        rr = src[0]["rr"]
    for s in src_labels:
        assert_allclose(s["rr"], rr, atol=1e-7)
    assert len(src) == 1 and src.kind == "volume"
    klass = VolVectorSourceEstimate
    if not vector:
        klass = klass._scalar_class
    vertices = [src[0]["vertno"]]
    n_verts = len(src[0]["vertno"])
    n_times = 50
    data = vertex_values = np.arange(1, n_verts + 1)
    end_shape = (n_times,)
    if vector:
        end_shape = (3,) + end_shape
        data = np.pad(data[:, np.newaxis], ((0, 0), (2, 0)), "constant")
    data = np.repeat(data[..., np.newaxis], n_times, -1)
    stcs = [klass(data.astype(float), vertices, 0, 1)]

    def eltc(*args, **kwargs):
        if call == "func":
            return extract_label_time_course(stcs, *args, **kwargs)
        else:
            assert call == "meth"
            return [stcs[0].extract_label_time_course(*args, **kwargs)]

    with pytest.raises(RuntimeError, match="atlas vox_mri_t does not match"):
        eltc(fname_fs_t1, src, mri_resolution=mri_res)
    assert len(src_labels) == 46  # includes unknown
    assert_array_equal(
        src[0]["vertno"],  # src includes some in "unknown" space
        np.sort(np.concatenate([s["vertno"] for s in src_labels])),
    )
    # spot check
    assert src_labels[-1]["seg_name"] == "CC_Anterior"
    assert src[0]["nuse"] == 4157
    assert len(src[0]["vertno"]) == 4157
    assert sum(s["nuse"] for s in src_labels) == 4157
    assert_array_equal(src_labels[-1]["vertno"], [8011, 8032, 8557])
    assert_array_equal(
        np.where(np.isin(src[0]["vertno"], [8011, 8032, 8557]))[0], [2672, 2688, 2995]
    )
    # triage "labels" argument
    if mri_res:
        # All should be there
        missing = []
    else:
        # Nearest misses these
        missing = [
            "Left-vessel",
            "Right-vessel",
            "5th-Ventricle",
            "non-WM-hypointensities",
        ]
    n_want = len(src_labels)
    if label_type is str:
        labels = fname_aseg
    elif label_type is list:
        labels = (fname_aseg, volume_labels)
    else:
        assert label_type is dict
        labels = (fname_aseg, {k: lut[k] for k in volume_labels})
        assert mri_res
        assert len(missing) == 0
        # we're going to add one that won't exist
        missing = ["intentionally_bad"]
        labels[1][missing[0]] = 10000
        n_want += 1
        n_tot += 1
    n_want -= len(missing)

    # actually do the testing
    if cf == "head" and not mri_res:  # some missing
        with pytest.warns(RuntimeWarning, match="any vertices"):
            eltc(labels, src, allow_empty=True, mri_resolution=mri_res)
    for mode in ("mean", "max"):
        with catch_logging() as log:
            label_tc = eltc(
                labels,
                src,
                mode=mode,
                allow_empty="ignore",
                mri_resolution=mri_res,
                verbose=True,
            )
        log = log.getvalue()
        assert re.search("^Reading atlas.*aseg\\.mgz\n", log) is not None
        if len(missing):
            # assert that the missing ones get logged
            assert "does not contain" in log
            assert repr(missing) in log
        else:
            assert "does not contain" not in log
        assert f"\n{n_want}/{n_tot} atlas regions had at least" in log
        assert len(label_tc) == 1
        label_tc = label_tc[0]
        assert label_tc.shape == (n_tot,) + end_shape
        if vector:
            assert_array_equal(label_tc[:, :2], 0.0)
            label_tc = label_tc[:, 2]
        assert label_tc.shape == (n_tot, n_times)
        # let's test some actual values by trusting the masks provided by
        # setup_volume_source_space. mri_resolution=True does some
        # interpolation so we should not expect equivalence, False does
        # nearest so we should.
        if mri_res:
            rtol = 0.2 if mode == "mean" else 0.8  # max much more sensitive
        else:
            rtol = 0.0
        for si, s in enumerate(src_labels):
            func = dict(mean=np.mean, max=np.max)[mode]
            these = vertex_values[np.isin(src[0]["vertno"], s["vertno"])]
            assert len(these) == s["nuse"]
            if si == 0 and s["seg_name"] == "Unknown":
                continue  # unknown is crappy
            if s["nuse"] == 0:
                want = 0.0
                if mri_res:
                    # this one is totally due to interpolation, so no easy
                    # test here
                    continue
            else:
                want = func(these)
            assert_allclose(label_tc[si], want, atol=1e-6, rtol=rtol)
            # compare with in_label, only on every fourth for speed
            if test_label is not False and si % 4 == 0:
                label = s["seg_name"]
                if test_label is int:
                    label = lut[label]
                in_label = stcs[0].in_label(label, fname_aseg, src).data
                assert in_label.shape == (s["nuse"],) + end_shape
                if vector:
                    assert_array_equal(in_label[:, :2], 0.0)
                    in_label = in_label[:, 2]
                if want == 0:
                    assert in_label.shape[0] == 0
                else:
                    in_label = func(in_label)
                    assert_allclose(in_label, want, atol=1e-6, rtol=rtol)
        if mode == "mean" and not vector:  # check the reverse
            if label_type is dict:
                ctx = pytest.warns(RuntimeWarning, match="does not contain")
            else:
                ctx = nullcontext()
            with ctx:
                stc_back = labels_to_stc(labels, label_tc, src=src)
            assert stc_back.data.shape == stcs[0].data.shape
            corr = np.corrcoef(stc_back.data.ravel(), stcs[0].data.ravel())[0, 1]
            assert 0.6 < corr < 0.63
            assert_allclose(_varexp(label_tc, label_tc), 1.0)
            ve = _varexp(stc_back.data, stcs[0].data)
            assert 0.83 < ve < 0.85
            with _record_warnings():  # ignore no output
                label_tc_rt = extract_label_time_course(
                    stc_back, labels, src=src, mri_resolution=mri_res, allow_empty=True
                )
            assert label_tc_rt.shape == label_tc.shape
            corr = np.corrcoef(label_tc.ravel(), label_tc_rt.ravel())[0, 1]
            lower, upper = (0.99, 0.999) if mri_res else (0.95, 0.97)
            assert lower < corr < upper
            ve = _varexp(label_tc_rt, label_tc)
            lower, upper = (0.99, 0.999) if mri_res else (0.97, 0.99)
            assert lower < ve < upper


def _varexp(got, want):
    return max(
        1 - np.linalg.norm(got.ravel() - want.ravel()) ** 2 / np.linalg.norm(want) ** 2,
        0.0,
    )


@testing.requires_testing_data
def test_extract_label_time_course_equiv():
    """Test extraction of label time courses from stc equivalences."""
    label = read_labels_from_annot(
        "sample", "aparc", "lh", regexp="transv", subjects_dir=subjects_dir
    )
    assert len(label) == 1
    label = label[0]
    inv = read_inverse_operator(fname_inv)
    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0].crop(0, 0.01)
    stc = apply_inverse(evoked, inv, pick_ori="normal", label=label)
    stc_full = apply_inverse(evoked, inv, pick_ori="normal")
    stc_in_label = stc_full.in_label(label)
    mean = stc.extract_label_time_course(label, inv["src"])
    mean_2 = stc_in_label.extract_label_time_course(label, inv["src"])
    assert_allclose(mean, mean_2)
    inv["src"][0]["vertno"] = np.array([], int)
    assert len(stc_in_label.vertices[0]) == 22
    with pytest.raises(ValueError, match="22/22 left hemisphere.*missing"):
        stc_in_label.extract_label_time_course(label, inv["src"])


def _my_trans(data):
    """FFT that adds an additional dimension by repeating result."""
    data_t = fft(data)
    data_t = np.concatenate([data_t[:, :, None], data_t[:, :, None]], axis=2)
    return data_t, None


def test_transform_data():
    """Test applying linear (time) transform to data."""
    # make up some data
    n_sensors, n_vertices, n_times = 10, 20, 4
    kernel = rng.randn(n_vertices, n_sensors)
    sens_data = rng.randn(n_sensors, n_times)

    vertices = [np.arange(n_vertices)]
    data = np.dot(kernel, sens_data)

    for idx, tmin_idx, tmax_idx in zip(
        [None, np.arange(n_vertices // 2, n_vertices)], [None, 1], [None, 3]
    ):
        if idx is None:
            idx_use = slice(None, None)
        else:
            idx_use = idx

        data_f, _ = _my_trans(data[idx_use, tmin_idx:tmax_idx])

        for stc_data in (data, (kernel, sens_data)):
            stc = VolSourceEstimate(stc_data, vertices=vertices, tmin=0.0, tstep=1.0)
            stc_data_t = stc.transform_data(
                _my_trans, idx=idx, tmin_idx=tmin_idx, tmax_idx=tmax_idx
            )
            assert_allclose(data_f, stc_data_t)
    # bad sens_data
    sens_data = sens_data[..., np.newaxis]
    with pytest.raises(ValueError, match="sensor data must have 2"):
        VolSourceEstimate((kernel, sens_data), vertices, 0, 1)


def test_transform():
    """Test applying linear (time) transform to data."""
    # make up some data
    n_verts_lh, n_verts_rh, n_times = 10, 10, 10
    vertices = [np.arange(n_verts_lh), n_verts_lh + np.arange(n_verts_rh)]
    data = rng.randn(n_verts_lh + n_verts_rh, n_times)
    stc = SourceEstimate(data, vertices=vertices, tmin=-0.1, tstep=0.1)

    # data_t.ndim > 2 & copy is True
    stcs_t = stc.transform(_my_trans, copy=True)
    assert isinstance(stcs_t, list)
    assert_array_equal(stc.times, stcs_t[0].times)
    assert_equal(stc.vertices, stcs_t[0].vertices)

    data = np.concatenate(
        (stcs_t[0].data[:, :, None], stcs_t[1].data[:, :, None]), axis=2
    )
    data_t = stc.transform_data(_my_trans)
    assert_array_equal(data, data_t)  # check against stc.transform_data()

    # data_t.ndim > 2 & copy is False
    pytest.raises(ValueError, stc.transform, _my_trans, copy=False)

    # data_t.ndim = 2 & copy is True
    tmp = deepcopy(stc)
    stc_t = stc.transform(np.abs, copy=True)
    assert isinstance(stc_t, SourceEstimate)
    assert_array_equal(stc.data, tmp.data)  # xfrm doesn't modify original?

    # data_t.ndim = 2 & copy is False
    times = np.round(1000 * stc.times)
    verts = np.arange(len(stc.lh_vertno), len(stc.lh_vertno) + len(stc.rh_vertno), 1)
    verts_rh = stc.rh_vertno
    tmin_idx = np.searchsorted(times, 0)
    tmax_idx = np.searchsorted(times, 501)  # Include 500ms in the range
    data_t = stc.transform_data(np.abs, idx=verts, tmin_idx=tmin_idx, tmax_idx=tmax_idx)
    stc.transform(np.abs, idx=verts, tmin=-50, tmax=500, copy=False)
    assert isinstance(stc, SourceEstimate)
    assert_equal(stc.tmin, 0.0)
    assert_equal(stc.times[-1], 0.5)
    assert_equal(len(stc.vertices[0]), 0)
    assert_equal(stc.vertices[1], verts_rh)
    assert_array_equal(stc.data, data_t)

    times = np.round(1000 * stc.times)
    tmin_idx, tmax_idx = np.searchsorted(times, 0), np.searchsorted(times, 250)
    data_t = stc.transform_data(np.abs, tmin_idx=tmin_idx, tmax_idx=tmax_idx)
    stc.transform(np.abs, tmin=0, tmax=250, copy=False)
    assert_equal(stc.tmin, 0.0)
    assert_equal(stc.times[-1], 0.2)
    assert_array_equal(stc.data, data_t)


def test_spatio_temporal_tris_adjacency():
    """Test spatio-temporal adjacency from triangles."""
    pytest.importorskip("sklearn")
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    adjacency = spatio_temporal_tris_adjacency(tris, 2)
    x = [1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
    components = stats.cluster_level._get_components(np.array(x), adjacency)
    # _get_components works differently now...
    old_fmt = [0, 0, -2, -2, -2, -2, 0, -2, -2, -2, -2, 1]
    new_fmt = np.array(old_fmt)
    new_fmt = [np.nonzero(new_fmt == v)[0] for v in np.unique(new_fmt[new_fmt >= 0])]
    assert len(new_fmt) == len(components)
    for c, n in zip(components, new_fmt):
        assert_array_equal(c, n)


@testing.requires_testing_data
def test_spatio_temporal_src_adjacency():
    """Test spatio-temporal adjacency from source spaces."""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    adjacency = spatio_temporal_tris_adjacency(tris, 2).todense()
    assert_allclose(np.diag(adjacency), 1.0)
    src[0]["use_tris"] = np.array([[0, 1, 2]])
    src[1]["use_tris"] = np.array([[0, 1, 2]])
    src[0]["vertno"] = np.array([0, 1, 2])
    src[1]["vertno"] = np.array([0, 1, 2])
    src[0]["type"] = "surf"
    src[1]["type"] = "surf"
    adjacency2 = spatio_temporal_src_adjacency(src, 2)
    assert_array_equal(adjacency2.todense(), adjacency)
    # add test for dist adjacency
    src[0]["dist"] = np.ones((3, 3)) - np.eye(3)
    src[1]["dist"] = np.ones((3, 3)) - np.eye(3)
    src[0]["vertno"] = [0, 1, 2]
    src[1]["vertno"] = [0, 1, 2]
    src[0]["type"] = "surf"
    src[1]["type"] = "surf"
    adjacency3 = spatio_temporal_src_adjacency(src, 2, dist=2)
    assert_array_equal(adjacency3.todense(), adjacency)
    # add test for source space adjacency with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    src_ = inverse_operator["src"]
    with pytest.warns(RuntimeWarning, match="will have holes"):
        adjacency = spatio_temporal_src_adjacency(src_, n_times=2)
    a = adjacency.shape[0] / 2
    b = sum([s["nuse"] for s in inverse_operator["src"]])
    assert a == b

    assert_equal(grade_to_tris(5).shape, [40960, 3])


def test_to_data_frame():
    """Test stc Pandas exporter."""
    pytest.importorskip("pandas")
    n_vert, n_times = 10, 5
    vertices = [np.arange(n_vert, dtype=np.int64), np.empty(0, dtype=np.int64)]
    data = rng.randn(n_vert, n_times)
    stc_surf = SourceEstimate(
        data, vertices=vertices, tmin=0, tstep=1, subject="sample"
    )
    stc_vol = VolSourceEstimate(
        data, vertices=vertices[:1], tmin=0, tstep=1, subject="sample"
    )
    for stc in [stc_surf, stc_vol]:
        df = stc.to_data_frame()
        # test data preservation (first 2 dataframe elements are subj & time)
        assert_array_equal(df.values.T[2:], stc.data)
        # test long format
        df_long = stc.to_data_frame(long_format=True)
        assert len(df_long) == stc.data.size
        expected = ("subject", "time", "source", "value")
        assert set(expected) == set(df_long.columns)


@pytest.mark.parametrize("index", ("time", ["time", "subject"], None))
def test_to_data_frame_index(index):
    """Test index creation in stc Pandas exporter."""
    pytest.importorskip("pandas")
    n_vert, n_times = 10, 5
    vertices = [np.arange(n_vert, dtype=np.int64), np.empty(0, dtype=np.int64)]
    data = rng.randn(n_vert, n_times)
    stc = SourceEstimate(data, vertices=vertices, tmin=0, tstep=1, subject="sample")
    df = stc.to_data_frame(index=index)
    # test index setting
    if not isinstance(index, list):
        index = [index]
    assert list(df.index.names) == index
    # test that non-indexed data were present as columns
    non_index = list(set(["time", "subject"]) - set(index))
    if len(non_index):
        assert all(np.isin(non_index, df.columns))


@pytest.mark.parametrize("kind", ("surface", "mixed", "volume"))
@pytest.mark.parametrize("vector", (False, True))
@pytest.mark.parametrize("n_times", (5, 1))
def test_get_peak(kind, vector, n_times):
    """Test peak getter."""
    n_vert = 10
    vertices = [np.arange(n_vert)]
    if kind == "surface":
        klass = VectorSourceEstimate
        vertices += [np.empty(0, int)]
    elif kind == "mixed":
        klass = MixedVectorSourceEstimate
        vertices += [np.empty(0, int), np.empty(0, int)]
    else:
        assert kind == "volume"
        klass = VolVectorSourceEstimate
    data = np.zeros((n_vert, n_times))
    data[1, -1] = 1
    if vector:
        data = np.repeat(data[:, np.newaxis], 3, 1)
    else:
        klass = klass._scalar_class
    stc = klass(data, vertices, 0, 1)

    with pytest.raises(ValueError, match="out of bounds"):
        stc.get_peak(tmin=-100)
    with pytest.raises(ValueError, match="out of bounds"):
        stc.get_peak(tmax=90)
    with pytest.raises(ValueError, match="must be <=" if n_times > 1 else "out of"):
        stc.get_peak(tmin=0.002, tmax=0.001)

    vert_idx, time_idx = stc.get_peak()
    vertno = np.concatenate(stc.vertices)
    assert vert_idx in vertno
    assert time_idx in stc.times
    data_idx, time_idx = stc.get_peak(vert_as_index=True, time_as_index=True)
    if vector:
        use_data = stc.magnitude().data
    else:
        use_data = stc.data
    assert data_idx == 1
    assert time_idx == n_times - 1
    assert data_idx == np.argmax(np.abs(use_data[:, time_idx]))
    assert time_idx == np.argmax(np.abs(use_data[data_idx, :]))
    if kind == "surface":
        data_idx_2, time_idx_2 = stc.get_peak(
            vert_as_index=True, time_as_index=True, hemi="lh"
        )
        assert data_idx_2 == data_idx
        assert time_idx_2 == time_idx
        with pytest.raises(RuntimeError, match="no vertices"):
            stc.get_peak(hemi="rh")


@testing.requires_testing_data
def test_mixed_stc(tmp_path):
    """Test source estimate from mixed source space."""
    pytest.importorskip("h5io")
    N = 90  # number of sources
    T = 2  # number of time points
    S = 3  # number of source spaces

    data = rng.randn(N, T)
    vertno = S * [np.arange(N // S)]

    # make sure error is raised if vertices are not a list of length >= 2
    pytest.raises(ValueError, MixedSourceEstimate, data=data, vertices=[np.arange(N)])

    stc = MixedSourceEstimate(data, vertno, 0, 1)

    # make sure error is raised for plotting surface with volume source
    fname = tmp_path / "mixed-stc.h5"
    stc.save(fname)
    stc_out = read_source_estimate(fname)
    assert_array_equal(stc_out.vertices, vertno)
    assert_array_equal(stc_out.data, data)
    assert stc_out.tmin == 0
    assert stc_out.tstep == 1
    assert isinstance(stc_out, MixedSourceEstimate)


@pytest.mark.parametrize(
    "klass, kind",
    [
        (VectorSourceEstimate, "surf"),
        (VolVectorSourceEstimate, "vol"),
        (VolVectorSourceEstimate, "discrete"),
        (MixedVectorSourceEstimate, "mixed"),
    ],
)
@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128])
def test_vec_stc_basic(tmp_path, klass, kind, dtype):
    """Test (vol)vector source estimate."""
    pytest.importorskip("h5io")
    nn = np.array(
        [
            [1, 0, 0],
            [0, 1, 0],
            [np.sqrt(1.0 / 2.0), 0, np.sqrt(1.0 / 2.0)],
            [np.sqrt(1 / 3.0)] * 3,
        ],
        np.float32,
    )

    data = np.array(
        [
            [1, 0, 0],
            [0, 2, 0],
            [-3, 0, 0],
            [1, 1, 1],
        ],
        dtype,
    )[:, :, np.newaxis]
    amplitudes = np.array([1, 2, 3, np.sqrt(3)], dtype)
    magnitudes = amplitudes.copy()
    normals = np.array([1, 2, -3.0 / np.sqrt(2), np.sqrt(3)], dtype)
    if dtype in (np.complex64, np.complex128):
        data *= 1j
        amplitudes *= 1j
        normals *= 1j
    directions = np.array([[1, 0, 0], [0, 1, 0], [-1, 0, 0], [1.0 / np.sqrt(3)] * 3])
    vol_kind = kind if kind in ("discrete", "vol") else "vol"
    vol_src = SourceSpaces([dict(nn=nn, type=vol_kind)])
    assert vol_src.kind == dict(vol="volume").get(vol_kind, vol_kind)
    vol_verts = [np.arange(4)]
    surf_src = SourceSpaces(
        [dict(nn=nn[:2], type="surf"), dict(nn=nn[2:], type="surf")]
    )
    assert surf_src.kind == "surface"
    surf_verts = [np.array([0, 1]), np.array([0, 1])]
    if klass is VolVectorSourceEstimate:
        src = vol_src
        verts = vol_verts
    elif klass is VectorSourceEstimate:
        src = surf_src
        verts = surf_verts
    if klass is MixedVectorSourceEstimate:
        src = surf_src + vol_src
        verts = surf_verts + vol_verts
        assert src.kind == "mixed"
        data = np.tile(data, (2, 1, 1))
        amplitudes = np.tile(amplitudes, 2)
        magnitudes = np.tile(magnitudes, 2)
        normals = np.tile(normals, 2)
        directions = np.tile(directions, (2, 1))
    stc = klass(data, verts, 0, 1, "foo")
    amplitudes = amplitudes[:, np.newaxis]
    magnitudes = magnitudes[:, np.newaxis]

    # Magnitude of the vectors
    assert_array_equal(stc.magnitude().data, magnitudes)

    # Vector components projected onto the vertex normals
    if kind in ("vol", "mixed"):
        with pytest.raises(RuntimeError, match="surface or discrete"):
            stc.project("normal", src)[0]
    else:
        normal = stc.project("normal", src)[0]
        assert_allclose(normal.data[:, 0], normals)

    # Maximal-variance component, either to keep amps pos or to align to src-nn
    projected, got_directions = stc.project("pca")
    assert_allclose(got_directions, directions)
    assert_allclose(projected.data, amplitudes)
    projected, got_directions = stc.project("pca", src)
    flips = np.array([[1], [1], [-1.0], [1]])
    if klass is MixedVectorSourceEstimate:
        flips = np.tile(flips, (2, 1))
    assert_allclose(got_directions, directions * flips)
    assert_allclose(projected.data, amplitudes * flips)

    out_name = tmp_path / "temp.h5"
    stc.save(out_name)
    stc_read = read_source_estimate(out_name)
    assert_allclose(stc.data, stc_read.data)
    assert len(stc.vertices) == len(stc_read.vertices)
    for v1, v2 in zip(stc.vertices, stc_read.vertices):
        assert_array_equal(v1, v2)

    stc = klass(data[:, :, 0], verts, 0, 1)  # upbroadcast
    assert stc.data.shape == (len(data), 3, 1)
    # Bad data
    with pytest.raises(ValueError, match="must have shape.*3"):
        klass(data[:, :2], verts, 0, 1)
    data = data[:, :, np.newaxis]
    with pytest.raises(ValueError, match="3 dimensions for .*VectorSource"):
        klass(data, verts, 0, 1)


@pytest.mark.parametrize("real", (True, False))
def test_source_estime_project(real):
    """Test projecting a source estimate onto direction of max power."""
    n_src, n_times = 4, 100
    rng = np.random.RandomState(0)
    data = rng.randn(n_src, 3, n_times)
    if not real:
        data = data + 1j * rng.randn(n_src, 3, n_times)
        assert data.dtype == np.complex128
    else:
        assert data.dtype == np.float64

    # Make sure that the normal we get maximizes the power
    # (i.e., minimizes the negative power)
    want_nn = np.empty((n_src, 3))
    for ii in range(n_src):
        x0 = np.ones(3)

        def objective(x):
            x = x / np.linalg.norm(x)
            return -np.linalg.norm(np.dot(x, data[ii]))

        want_nn[ii] = fmin_cobyla(objective, x0, (), rhobeg=0.1, rhoend=1e-6)
    want_nn /= np.linalg.norm(want_nn, axis=1, keepdims=True)

    stc = VolVectorSourceEstimate(data, [np.arange(n_src)], 0, 1)
    _, directions = stc.project("pca")
    flips = np.sign(np.sum(directions * want_nn, axis=1, keepdims=True))
    directions *= flips
    assert_allclose(directions, want_nn, atol=2e-6)


@testing.requires_testing_data
def test_source_estime_project_label():
    """Test projecting a source estimate onto direction of max power."""
    fwd = read_forward_solution(fname_fwd)
    fwd = pick_types_forward(fwd, meg=True, eeg=False)

    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0]
    noise_cov = read_cov(fname_cov)
    free = make_inverse_operator(evoked.info, fwd, noise_cov, loose=1.0)
    stc_free = apply_inverse(evoked, free, pick_ori="vector")

    stc_pca = stc_free.project("pca", fwd["src"])[0]

    labels_lh = read_labels_from_annot(
        "sample", "aparc", "lh", subjects_dir=subjects_dir
    )
    new_label = labels_lh[0] + labels_lh[1]

    stc_in_label = stc_free.in_label(new_label)
    stc_pca_in_label = stc_pca.in_label(new_label)

    stc_in_label_pca = stc_in_label.project("pca", fwd["src"])[0]
    assert_array_equal(stc_pca_in_label.data, stc_in_label_pca.data)


@pytest.fixture(scope="module", params=[testing._pytest_param()])
def invs():
    """Inverses of various amounts of loose."""
    fwd = read_forward_solution(fname_fwd)
    fwd = pick_types_forward(fwd, meg=True, eeg=False)
    fwd_surf = convert_forward_solution(fwd, surf_ori=True)
    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0]
    noise_cov = read_cov(fname_cov)
    free = make_inverse_operator(evoked.info, fwd, noise_cov, loose=1.0)
    free_surf = make_inverse_operator(evoked.info, fwd_surf, noise_cov, loose=1.0)
    freeish = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.9999)
    fixed = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.0)
    fixedish = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.0001)
    assert_allclose(
        free["source_nn"], np.kron(np.ones(fwd["nsource"]), np.eye(3)).T, atol=1e-7
    )
    # This is the one exception:
    assert not np.allclose(free["source_nn"], free_surf["source_nn"])
    assert_allclose(
        free["source_nn"], np.tile(np.eye(3), (free["nsource"], 1)), atol=1e-7
    )
    # All others are similar:
    for other in (freeish, fixedish):
        assert_allclose(free_surf["source_nn"], other["source_nn"], atol=1e-7)
    assert_allclose(free_surf["source_nn"][2::3], fixed["source_nn"], atol=1e-7)
    expected_nn = np.concatenate([_get_src_nn(s) for s in fwd["src"]])
    assert_allclose(fixed["source_nn"], expected_nn, atol=1e-7)
    return evoked, free, free_surf, freeish, fixed, fixedish


@pytest.mark.parametrize("pick_ori", [None, "normal", "vector"])
def test_vec_stc_inv_free(invs, pick_ori):
    """Test vector STC behavior with two free-orientation inverses."""
    evoked, free, free_surf, _, _, _ = invs
    stc_free = apply_inverse(evoked, free, pick_ori=pick_ori)
    stc_free_surf = apply_inverse(evoked, free_surf, pick_ori=pick_ori)
    assert_allclose(stc_free.data, stc_free_surf.data, atol=1e-5)


@pytest.mark.parametrize("pick_ori", [None, "normal", "vector"])
def test_vec_stc_inv_free_surf(invs, pick_ori):
    """Test vector STC behavior with free and free-ish orientation invs."""
    evoked, _, free_surf, freeish, _, _ = invs
    stc_free = apply_inverse(evoked, free_surf, pick_ori=pick_ori)
    stc_freeish = apply_inverse(evoked, freeish, pick_ori=pick_ori)
    assert_allclose(stc_free.data, stc_freeish.data, atol=1e-3)


@pytest.mark.parametrize("pick_ori", (None, "normal", "vector"))
def test_vec_stc_inv_fixed(invs, pick_ori):
    """Test vector STC behavior with fixed-orientation inverses."""
    evoked, _, _, _, fixed, fixedish = invs
    stc_fixed = apply_inverse(evoked, fixed)
    stc_fixed_vector = apply_inverse(evoked, fixed, pick_ori="vector")
    assert_allclose(
        stc_fixed.data, stc_fixed_vector.project("normal", fixed["src"])[0].data
    )
    stc_fixedish = apply_inverse(evoked, fixedish, pick_ori=pick_ori)
    if pick_ori == "vector":
        assert_allclose(stc_fixed_vector.data, stc_fixedish.data, atol=1e-2)
        # two ways here: with magnitude...
        assert_allclose(abs(stc_fixed).data, stc_fixedish.magnitude().data, atol=1e-2)
        # ... and when picking the normal (signed)
        stc_fixedish = stc_fixedish.project("normal", fixedish["src"])[0]
    elif pick_ori is None:
        stc_fixed = abs(stc_fixed)
    else:
        assert pick_ori == "normal"  # no need to modify
    assert_allclose(stc_fixed.data, stc_fixedish.data, atol=1e-2)


@testing.requires_testing_data
def test_epochs_vector_inverse():
    """Test vector inverse consistency between evoked and epochs."""
    raw = read_raw_fif(fname_raw)
    events = find_events(raw, stim_channel="STI 014")[:2]
    reject = dict(grad=2000e-13, mag=4e-12, eog=150e-6)

    epochs = Epochs(
        raw, events, None, 0, 0.01, baseline=None, reject=reject, preload=True
    )

    assert_equal(len(epochs), 2)

    evoked = epochs.average(picks=range(len(epochs.ch_names)))

    inv = read_inverse_operator(fname_inv)

    method = "MNE"
    snr = 3.0
    lambda2 = 1.0 / snr**2

    stcs_epo = apply_inverse_epochs(
        epochs, inv, lambda2, method=method, pick_ori="vector", return_generator=False
    )
    stc_epo = np.mean(stcs_epo)

    stc_evo = apply_inverse(evoked, inv, lambda2, method=method, pick_ori="vector")

    assert_allclose(stc_epo.data, stc_evo.data, rtol=1e-9, atol=0)


@testing.requires_testing_data
def test_vol_adjacency():
    """Test volume adjacency."""
    pytest.importorskip("sklearn")
    vol = read_source_spaces(fname_vsrc)

    pytest.raises(ValueError, spatial_src_adjacency, vol, dist=1.0)

    adjacency = spatial_src_adjacency(vol)
    n_vertices = vol[0]["inuse"].sum()
    assert_equal(adjacency.shape, (n_vertices, n_vertices))
    assert np.all(adjacency.data == 1)
    assert isinstance(adjacency, sparse.coo_array)

    adjacency2 = spatio_temporal_src_adjacency(vol, n_times=2)
    assert_equal(adjacency2.shape, (2 * n_vertices, 2 * n_vertices))
    assert np.all(adjacency2.data == 1)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_spatial_src_adjacency():
    """Test spatial adjacency functionality."""
    # oct
    src = read_source_spaces(fname_src)
    assert src[0]["dist"] is not None  # distance info
    with pytest.warns(RuntimeWarning, match="will have holes"):
        con = spatial_src_adjacency(src).toarray()
    con_dist = spatial_src_adjacency(src, dist=0.01).toarray()
    assert (con == con_dist).mean() > 0.75
    # ico
    src = read_source_spaces(fname_src_fs)
    con = spatial_src_adjacency(src).tocsr()
    con_tris = spatial_tris_adjacency(grade_to_tris(5)).tocsr()
    assert con.shape == con_tris.shape
    assert_array_equal(con.data, con_tris.data)
    assert_array_equal(con.indptr, con_tris.indptr)
    assert_array_equal(con.indices, con_tris.indices)
    # one hemi
    con_lh = spatial_src_adjacency(src[:1]).tocsr()
    con_lh_tris = spatial_tris_adjacency(grade_to_tris(5)).tocsr()
    con_lh_tris = con_lh_tris[:10242, :10242].tocsr()
    assert_array_equal(con_lh.data, con_lh_tris.data)
    assert_array_equal(con_lh.indptr, con_lh_tris.indptr)
    assert_array_equal(con_lh.indices, con_lh_tris.indices)


@testing.requires_testing_data
def test_vol_mask():
    """Test extraction of volume mask."""
    pytest.importorskip("nibabel")
    pytest.importorskip("sklearn")
    src = read_source_spaces(fname_vsrc)
    mask = _get_vol_mask(src)
    # Let's use an alternative way that should be equivalent
    vertices = [src[0]["vertno"]]
    n_vertices = len(vertices[0])
    data = (1 + np.arange(n_vertices))[:, np.newaxis]
    stc_tmp = VolSourceEstimate(data, vertices, tmin=0.0, tstep=1.0)
    img = stc_tmp.as_volume(src, mri_resolution=False)
    img_data = _get_img_fdata(img)[:, :, :, 0].T
    mask_nib = img_data != 0
    assert_array_equal(img_data[mask_nib], data[:, 0])
    assert_array_equal(np.where(mask_nib.ravel())[0], src[0]["vertno"])
    assert_array_equal(mask, mask_nib)
    assert_array_equal(img_data.shape, mask.shape)


@testing.requires_testing_data
def test_stc_near_sensors(tmp_path):
    """Test stc_near_sensors."""
    info = read_info(fname_evoked)
    # pick the left EEG sensors
    picks = pick_types(info, meg=False, eeg=True, exclude=())
    picks = [pick for pick in picks if info["chs"][pick]["loc"][0] < 0]
    pick_info(info, picks, copy=False)
    with info._unlock():
        info["projs"] = []
    info["bads"] = []
    assert info["nchan"] == 33
    evoked = EvokedArray(np.eye(info["nchan"]), info)
    trans = read_trans(fname_fwd)
    assert trans["to"] == FIFF.FIFFV_COORD_HEAD
    # testing does not have pial, so fake it
    os.makedirs(tmp_path / "sample" / "surf")
    for hemi in ("lh", "rh"):
        copyfile(
            subjects_dir / "sample" / "surf" / f"{hemi}.white",
            tmp_path / "sample" / "surf" / f"{hemi}.pial",
        )
    # here we use a distance is smaller than the inter-sensor distance
    kwargs = dict(
        subject="sample",
        trans=trans,
        subjects_dir=tmp_path,
        verbose=True,
        distance=0.005,
    )
    with pytest.raises(ValueError, match="No appropriate channels"):
        stc_near_sensors(evoked, **kwargs)
    evoked.set_channel_types({ch_name: "ecog" for ch_name in evoked.ch_names})
    with catch_logging() as log:
        stc = stc_near_sensors(evoked, **kwargs)
    log = log.getvalue()
    assert "Minimum projected intra-sensor distance: 7." in log  # 7.4
    # this should be left-hemisphere dominant
    assert 5000 > len(stc.vertices[0]) > 4000
    assert 200 > len(stc.vertices[1]) > 100
    # and at least one vertex should have the channel values
    dists = cdist(stc.data, evoked.data)
    assert np.isclose(dists, 0.0, atol=1e-6).any(0).all()

    src = read_source_spaces(fname_src)  # uses "white" but should be okay
    src._transform_to("head", trans)
    assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD
    stc_src = stc_near_sensors(evoked, src=src, **kwargs)
    assert len(stc_src.data) == 7928
    with pytest.warns(RuntimeWarning, match="not included"):  # some removed
        stc_src_full = compute_source_morph(
            stc_src,
            "sample",
            "sample",
            smooth=5,
            spacing=None,
            subjects_dir=subjects_dir,
        ).apply(stc_src)
    lh_idx = np.searchsorted(stc_src_full.vertices[0], stc.vertices[0])
    rh_idx = np.searchsorted(stc_src_full.vertices[1], stc.vertices[1])
    rh_idx += len(stc_src_full.vertices[0])
    sub_data = stc_src_full.data[np.concatenate([lh_idx, rh_idx])]
    assert sub_data.shape == stc.data.shape
    corr = np.corrcoef(stc.data.ravel(), sub_data.ravel())[0, 1]
    assert 0.6 < corr < 0.7

    # now single-weighting mode
    stc_w = stc_near_sensors(evoked, mode="single", **kwargs)
    assert_array_less(stc_w.data, stc.data + 1e-3)  # some tol
    assert len(stc_w.data) == len(stc.data)
    # at least one for each sensor should have projected right on it
    dists = cdist(stc_w.data, evoked.data)
    assert np.isclose(dists, 0.0, atol=1e-6).any(0).all()

    # finally, nearest mode: all should match
    stc_n = stc_near_sensors(evoked, mode="nearest", **kwargs)
    assert len(stc_n.data) == len(stc.data)
    # at least one for each sensor should have projected right on it
    dists = cdist(stc_n.data, evoked.data)
    assert np.isclose(dists, 0.0, atol=1e-6).any(1).all()  # all vert eq some ch

    # these are EEG electrodes, so the distance 0.01 is too small for the
    # scalp+skull. Even at a distance of 33 mm EEG 060 is too far:
    with pytest.warns(RuntimeWarning, match="Channel missing in STC: EEG 060"):
        stc = stc_near_sensors(
            evoked,
            trans,
            "sample",
            subjects_dir=tmp_path,
            project=False,
            distance=0.033,
        )
    assert stc.data.any(0).sum() == len(evoked.ch_names) - 1

    # and now with volumetric projection
    src = read_source_spaces(fname_vsrc)
    with catch_logging() as log:
        stc_vol = stc_near_sensors(
            evoked,
            trans,
            "sample",
            src=src,
            surface=None,
            subjects_dir=subjects_dir,
            distance=0.033,
            verbose=True,
        )
    assert isinstance(stc_vol, VolSourceEstimate)
    log = log.getvalue()
    assert "4157 volume vertices" in log


@testing.requires_testing_data
def test_stc_near_sensors_picks():
    """Test using picks with stc_near_sensors."""
    pytest.importorskip("pymatreader")
    info = mne.io.read_raw_nirx(fname_nirx).info
    evoked = mne.EvokedArray(np.ones((len(info["ch_names"]), 1)), info)
    src = mne.read_source_spaces(fname_src_fs)
    kwargs = dict(
        evoked=evoked,
        subject="fsaverage",
        trans="fsaverage",
        subjects_dir=subjects_dir,
        src=src,
        surface=None,
        project=True,
    )
    with pytest.raises(ValueError, match="No appropriate channels"):
        stc_near_sensors(**kwargs)
    picks = np.arange(len(info["ch_names"]))
    data = stc_near_sensors(picks=picks, **kwargs).data
    assert len(data) == 20484
    assert (data >= 0).all()
    data = data[data > 0]
    n_pts = len(data)
    assert 500 < n_pts < 600
    lo, hi = np.percentile(data, (5, 95))
    assert 0.01 < lo < 0.1
    assert 1.3 < hi < 1.7  # > 1
    data = stc_near_sensors(picks=picks, mode="weighted", **kwargs).data
    assert (data >= 0).all()
    data = data[data > 0]
    assert len(data) == n_pts
    assert_array_equal(data, 1.0)  # values preserved


def _make_morph_map_hemi_same(subject_from, subject_to, subjects_dir, reg_from, reg_to):
    return _make_morph_map_hemi(
        subject_from, subject_from, subjects_dir, reg_from, reg_from
    )


@testing.requires_testing_data
@pytest.mark.parametrize(
    "kind",
    (
        pytest.param("volume", marks=[pytest.mark.slowtest]),
        "surface",
    ),
)
@pytest.mark.parametrize("scale", ((1.0, 0.8, 1.2), 1.0, 0.9))
def test_scale_morph_labels(kind, scale, monkeypatch, tmp_path):
    """Test label extraction, morphing, and MRI scaling relationships."""
    pytest.importorskip("nibabel")
    if kind == "volume":
        pytest.importorskip("dipy")
    subject_from = "sample"
    subject_to = "small"
    testing_dir = subjects_dir / subject_from
    from_dir = tmp_path / subject_from
    for root in ("mri", "surf", "label", "bem"):
        os.makedirs(from_dir / root, exist_ok=True)
    for hemi in ("lh", "rh"):
        for root, fname in (
            ("surf", "sphere"),
            ("surf", "white"),
            ("surf", "sphere.reg"),
            ("label", "aparc.annot"),
        ):
            use_fname = Path(root) / f"{hemi}.{fname}"
            copyfile(testing_dir / use_fname, from_dir / use_fname)

    for root, fname in (("mri", "aseg.mgz"), ("mri", "brain.mgz")):
        use_fname = Path(root) / fname
        copyfile(testing_dir / use_fname, from_dir / use_fname)
    del testing_dir
    if kind == "surface":
        src_from = read_source_spaces(fname_src_3)
        assert src_from[0]["dist"] is None
        assert src_from[0]["nearest"] is not None
        # avoid patch calc
        src_from[0]["nearest"] = src_from[1]["nearest"] = None
        assert len(src_from) == 2
        assert src_from[0]["nuse"] == src_from[1]["nuse"] == 258
        klass = SourceEstimate
        labels_from = read_labels_from_annot(subject_from, subjects_dir=tmp_path)
        n_labels = len(labels_from)
        write_source_spaces(
            tmp_path / subject_from / "bem" / f"{subject_from}-oct-4-src.fif",
            src_from,
        )
    else:
        assert kind == "volume"
        pytest.importorskip("dipy")
        src_from = read_source_spaces(fname_src_vol)
        src_from[0]["subject_his_id"] = subject_from
        labels_from = tmp_path / subject_from / "mri" / "aseg.mgz"
        n_labels = 46
        assert labels_from.is_file()
        klass = VolSourceEstimate
        assert len(src_from) == 1
        assert src_from[0]["nuse"] == 4157
        write_source_spaces(from_dir / "bem" / "sample-vol20-src.fif", src_from)
    scale_mri(
        subject_from,
        subject_to,
        scale,
        subjects_dir=tmp_path,
        annot=True,
        skip_fiducials=True,
        verbose=True,
        overwrite=True,
    )
    if kind == "surface":
        src_to = read_source_spaces(
            tmp_path / subject_to / "bem" / f"{subject_to}-oct-4-src.fif"
        )
        labels_to = read_labels_from_annot(subject_to, subjects_dir=tmp_path)
        # Save time since we know these subjects are identical
        monkeypatch.setattr(
            mne.morph_map, "_make_morph_map_hemi", _make_morph_map_hemi_same
        )
    else:
        src_to = read_source_spaces(
            tmp_path / subject_to / "bem" / f"{subject_to}-vol20-src.fif"
        )
        labels_to = tmp_path / subject_to / "mri" / "aseg.mgz"
    # 1. Label->STC->Label for the given subject should be identity
    #    (for surfaces at least; for volumes it's not as clean as this
    #     due to interpolation)
    n_times = 50
    rng = np.random.RandomState(0)
    label_tc = rng.randn(n_labels, n_times)
    # check that a random permutation of our labels yields a terrible
    # correlation
    corr = np.corrcoef(label_tc.ravel(), rng.permutation(label_tc).ravel())[0, 1]
    assert -0.06 < corr < 0.06
    # project label activations to full source space
    with pytest.raises(ValueError, match="subject"):
        labels_to_stc(labels_from, label_tc, src=src_from, subject="foo")
    stc = labels_to_stc(labels_from, label_tc, src=src_from)
    assert stc.subject == "sample"
    assert isinstance(stc, klass)
    label_tc_from = extract_label_time_course(stc, labels_from, src_from, mode="mean")
    if kind == "surface":
        assert_allclose(label_tc, label_tc_from, rtol=1e-12, atol=1e-12)
    else:
        corr = np.corrcoef(label_tc.ravel(), label_tc_from.ravel())[0, 1]
        assert 0.93 < corr < 0.95

    #
    # 2. Changing STC subject to the surrogate and then extracting
    #
    stc.subject = subject_to
    label_tc_to = extract_label_time_course(stc, labels_to, src_to, mode="mean")
    assert_allclose(label_tc_from, label_tc_to, rtol=1e-12, atol=1e-12)
    stc.subject = subject_from

    #
    # 3. Morphing STC to new subject then extracting
    #
    if isinstance(scale, tuple) and kind == "volume":
        ctx = nullcontext()
        test_morph = True
    elif kind == "surface":
        ctx = pytest.warns(RuntimeWarning, match="not included")
        test_morph = True
    else:
        ctx = nullcontext()
        test_morph = True
    with ctx:  # vertices not included
        morph = compute_source_morph(
            src_from,
            subject_to=subject_to,
            src_to=src_to,
            subjects_dir=tmp_path,
            niter_sdr=(),
            smooth=1,
            zooms=14.0,
            verbose=True,
        )  # speed up with higher zooms
    if kind == "volume":
        got_affine = morph.pre_affine.affine
        want_affine = np.eye(4)
        want_affine.ravel()[::5][:3] = 1.0 / np.array(scale, float)
        # just a scaling (to within 1% if zooms=None, 20% with zooms=10)
        assert_allclose(want_affine[:, :3], got_affine[:, :3], atol=0.4)
        assert got_affine[3, 3] == 1.0
        # little translation (to within `limit` mm)
        move = np.linalg.norm(got_affine[:3, 3])
        limit = 2.0 if scale == 1.0 else 12
        assert move < limit, scale
    if test_morph:
        stc_to = morph.apply(stc)
        label_tc_to_morph = extract_label_time_course(
            stc_to, labels_to, src_to, mode="mean"
        )
        if kind == "volume":
            corr = np.corrcoef(label_tc.ravel(), label_tc_to_morph.ravel())[0, 1]
            if isinstance(scale, tuple):
                # some other fixed constant
                # min_, max_ = 0.84, 0.855  # zooms='auto' values
                min_, max_ = 0.55, 0.67
            elif scale == 1:
                # min_, max_ = 0.85, 0.875  # zooms='auto' values
                min_, max_ = 0.72, 0.76
            else:
                # min_, max_ = 0.84, 0.855  # zooms='auto' values
                min_, max_ = 0.44, 0.63
            assert min_ < corr <= max_, scale
        else:
            assert_allclose(label_tc, label_tc_to_morph, atol=1e-12, rtol=1e-12)

    #
    # 4. The same round trip from (1) but in the warped space
    #
    stc = labels_to_stc(labels_to, label_tc, src=src_to)
    assert isinstance(stc, klass)
    label_tc_to = extract_label_time_course(stc, labels_to, src_to, mode="mean")
    if kind == "surface":
        assert_allclose(label_tc, label_tc_to, rtol=1e-12, atol=1e-12)
    else:
        corr = np.corrcoef(label_tc.ravel(), label_tc_to.ravel())[0, 1]
        assert 0.93 < corr < 0.96, scale


@testing.requires_testing_data
@pytest.mark.parametrize(
    "kind",
    [
        "surface",
        pytest.param("volume", marks=[pytest.mark.slowtest]),
    ],
)
def test_label_extraction_subject(kind):
    """Test that label extraction subject is treated properly."""
    if kind == "surface":
        inv = read_inverse_operator(fname_inv)
        labels = read_labels_from_annot("sample", subjects_dir=subjects_dir)
        labels_fs = read_labels_from_annot("fsaverage", subjects_dir=subjects_dir)
        labels_fs = [
            label for label in labels_fs if not label.name.startswith("unknown")
        ]
        assert all(label.subject == "sample" for label in labels)
        assert all(label.subject == "fsaverage" for label in labels_fs)
        assert len(labels) == len(labels_fs) == 68
        n_labels = 68
    else:
        assert kind == "volume"
        pytest.importorskip("nibabel")
        inv = read_inverse_operator(fname_inv_vol)
        inv["src"][0]["subject_his_id"] = "sample"  # modernize
        labels = subjects_dir / "sample" / "mri" / "aseg.mgz"
        labels_fs = subjects_dir / "fsaverage" / "mri" / "aseg.mgz"
        n_labels = 46
    src = inv["src"]
    assert src.kind == kind
    assert src._subject == "sample"
    ave = read_evokeds(fname_evoked)[0].apply_baseline((None, 0)).crop(0, 0.01)
    assert len(ave.times) == 4
    stc = apply_inverse(ave, inv)
    assert stc.subject == "sample"
    ltc = extract_label_time_course(stc, labels, src)
    stc.subject = "fsaverage"
    with pytest.raises(ValueError, match=r"source spac.*not match.* stc\.sub"):
        extract_label_time_course(stc, labels, src)
    stc.subject = "sample"
    assert ltc.shape == (n_labels, 4)
    if kind == "volume":
        with pytest.raises(RuntimeError, match="atlas.*not match.*source spa"):
            extract_label_time_course(stc, labels_fs, src)
    else:
        with pytest.raises(ValueError, match=r"label\.sub.*not match.* stc\."):
            extract_label_time_course(stc, labels_fs, src)
        stc.subject = None
        with pytest.raises(ValueError, match=r"label\.sub.*not match.* sour"):
            extract_label_time_course(stc, labels_fs, src)


def test_apply_function_stc():
    """Check the apply_function method for source estimate data."""
    # Create a sample _BaseSourceEstimate object
    n_vertices = 100
    n_times = 200
    vertices = [np.array(np.arange(50)), np.array(np.arange(50, 100))]
    tmin = 0.0
    tstep = 0.001
    data = np.random.default_rng(0).normal(size=(n_vertices, n_times))

    stc = _make_stc(data, vertices, tmin=tmin, tstep=tstep, src_type="surface")

    # A sample function to apply to the data
    def fun(data_row, **kwargs):
        return 2 * data_row

    # Test applying the function to all vertices without parallelization
    stc_copy = stc.copy()
    stc.apply_function(fun)
    for idx in range(n_vertices):
        assert_allclose(stc.data[idx, :], 2 * stc_copy.data[idx, :])

    # Test applying the function with parallelization
    stc.apply_function(fun, n_jobs=2)
    for idx in range(n_vertices):
        assert_allclose(stc.data[idx, :], 4 * stc_copy.data[idx, :])