# Author: Mark Wronkiewicz <wronk@uw.edu>
#
# License: BSD-3-Clause

from contextlib import contextmanager
import os.path as op
import pathlib
import re

import numpy as np
from numpy.testing import assert_allclose, assert_array_equal
import pytest
from scipy import sparse
from scipy.special import sph_harm

import mne
from mne import compute_raw_covariance, pick_types, concatenate_raws, pick_info
from mne.annotations import _annotations_starts_stops
from mne.chpi import read_head_pos, filter_chpi
from mne.forward import _prep_meg_channels
from mne.datasets import testing
from mne.forward import use_coil_def
from mne.io import (read_raw_fif, read_info, read_raw_bti, read_raw_kit,
                    BaseRaw, read_raw_ctf)
from mne.io.constants import FIFF
from mne.preprocessing import (maxwell_filter, find_bad_channels_maxwell,
                               annotate_amplitude, compute_maxwell_basis,
                               maxwell_filter_prepare_emptyroom,
                               annotate_movement)
from mne.preprocessing.maxwell import (
    _get_n_moments, _sss_basis_basic, _sh_complex_to_real,
    _sh_real_to_complex, _sh_negate, _bases_complex_to_real, _trans_sss_basis,
    _bases_real_to_complex, _prep_mf_coils)
from mne.rank import _get_rank_sss, _compute_rank_int, compute_rank
from mne.utils import (assert_meg_snr, catch_logging, _record_warnings,
                       object_diff, buggy_mkl_svd, use_log_level)

io_path = op.join(op.dirname(__file__), '..', '..', 'io', 'tests', 'data')
raw_small_fname = op.join(io_path, 'test_raw.fif')

data_path = testing.data_path(download=False)
sss_path = op.join(data_path, 'SSS')
pre = op.join(sss_path, 'test_move_anon_')
raw_fname = pre + 'raw.fif'
sss_std_fname = pre + 'stdOrigin_raw_sss.fif'
sss_nonstd_fname = pre + 'nonStdOrigin_raw_sss.fif'
sss_bad_recon_fname = pre + 'badRecon_raw_sss.fif'
sss_reg_in_fname = pre + 'regIn_raw_sss.fif'
sss_fine_cal_fname = pre + 'fineCal_raw_sss.fif'
sss_ctc_fname = pre + 'crossTalk_raw_sss.fif'
sss_trans_default_fname = pre + 'transDefault_raw_sss.fif'
sss_trans_sample_fname = pre + 'transSample_raw_sss.fif'
sss_st1FineCalCrossTalkRegIn_fname = \
    pre + 'st1FineCalCrossTalkRegIn_raw_sss.fif'
sss_st1FineCalCrossTalkRegInTransSample_fname = \
    pre + 'st1FineCalCrossTalkRegInTransSample_raw_sss.fif'
sss_movecomp_fname = pre + 'movecomp_raw_sss.fif'
sss_movecomp_reg_in_fname = pre + 'movecomp_regIn_raw_sss.fif'
sss_movecomp_reg_in_st4s_fname = pre + 'movecomp_regIn_st4s_raw_sss.fif'
skip_fname = op.join(data_path, 'misc', 'intervalrecording_raw.fif')

erm_fname = pre + 'erm_raw.fif'
sss_erm_std_fname = pre + 'erm_devOrigin_raw_sss.fif'
sss_erm_reg_in_fname = pre + 'erm_regIn_raw_sss.fif'
sss_erm_fine_cal_fname = pre + 'erm_fineCal_raw_sss.fif'
sss_erm_ctc_fname = pre + 'erm_crossTalk_raw_sss.fif'
sss_erm_st_fname = pre + 'erm_st1_raw_sss.fif'
sss_erm_st1FineCalCrossTalk_fname = pre + 'erm_st1FineCalCrossTalk_raw_sss.fif'
sss_erm_st1FineCalCrossTalkRegIn_fname = \
    pre + 'erm_st1FineCalCrossTalkRegIn_raw_sss.fif'

sample_fname = op.join(
    data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
sss_samp_reg_in_fname = op.join(data_path, 'SSS',
                                'sample_audvis_trunc_regIn_raw_sss.fif')
sss_samp_fname = op.join(data_path, 'SSS', 'sample_audvis_trunc_raw_sss.fif')

pos_fname = op.join(data_path, 'SSS', 'test_move_anon_raw.pos')

bases_fname = op.join(sss_path, 'sss_data.mat')
fine_cal_fname = op.join(sss_path, 'sss_cal_3053.dat')
fine_cal_fname_3d = op.join(sss_path, 'sss_cal_3053_3d.dat')
ctc_fname = op.join(sss_path, 'ct_sparse.fif')
fine_cal_mgh_fname = op.join(sss_path, 'sss_cal_mgh.dat')
ctc_mgh_fname = op.join(sss_path, 'ct_sparse_mgh.fif')

triux_path = op.join(data_path, 'SSS', 'TRIUX')
tri_fname = op.join(triux_path, 'triux_bmlhus_erm_raw.fif')
tri_sss_fname = op.join(triux_path, 'triux_bmlhus_erm_raw_sss.fif')
tri_sss_reg_fname = op.join(triux_path, 'triux_bmlhus_erm_regIn_raw_sss.fif')
tri_sss_st4_fname = op.join(triux_path, 'triux_bmlhus_erm_st4_raw_sss.fif')
tri_sss_ctc_fname = op.join(triux_path, 'triux_bmlhus_erm_ctc_raw_sss.fif')
tri_sss_cal_fname = op.join(triux_path, 'triux_bmlhus_erm_cal_raw_sss.fif')
tri_sss_ctc_cal_fname = op.join(
    triux_path, 'triux_bmlhus_erm_ctc_cal_raw_sss.fif')
tri_sss_ctc_cal_reg_in_fname = op.join(
    triux_path, 'triux_bmlhus_erm_ctc_cal_regIn_raw_sss.fif')
tri_ctc_fname = op.join(triux_path, 'ct_sparse_BMLHUS.fif')
tri_cal_fname = op.join(triux_path, 'sss_cal_BMLHUS.dat')

io_dir = op.join(op.dirname(__file__), '..', '..', 'io')
fname_ctf_raw = op.join(io_dir, 'tests', 'data', 'test_ctf_comp_raw.fif')
ctf_fname_continuous = op.join(data_path, 'CTF', 'testdata_ctf.ds')

# In some of the tests, use identical coil defs to what is used in
# MaxFilter
elekta_def_fname = op.join(op.dirname(mne.__file__), 'data',
                           'coil_def_Elekta.dat')

int_order, ext_order = 8, 3
mf_head_origin = (0., 0., 0.04)
mf_meg_origin = (0., 0.013, -0.006)

# 30 random bad MEG channels (20 grad, 10 mag) that were used in generation
bads = ['MEG0912', 'MEG1722', 'MEG2213', 'MEG0132', 'MEG1312', 'MEG0432',
        'MEG2433', 'MEG1022', 'MEG0442', 'MEG2332', 'MEG0633', 'MEG1043',
        'MEG1713', 'MEG0422', 'MEG0932', 'MEG1622', 'MEG1343', 'MEG0943',
        'MEG0643', 'MEG0143', 'MEG2142', 'MEG0813', 'MEG2143', 'MEG1323',
        'MEG0522', 'MEG1123', 'MEG0423', 'MEG2122', 'MEG2532', 'MEG0812']


def _assert_n_free(raw_sss, lower, upper=None):
    """Check the DOF."""
    upper = lower if upper is None else upper
    n_free = raw_sss.info['proc_history'][0]['max_info']['sss_info']['nfree']
    assert lower <= n_free <= upper, \
        'nfree fail: %s <= %s <= %s' % (lower, n_free, upper)


def _assert_mag_coil_type(info, coil_type):
    __tracebackhide__ = True
    picks = pick_types(info, meg='mag', exclude=())
    coil_types = set(info['chs'][pick]['coil_type'] for pick in picks)
    assert coil_types == {coil_type}


def read_crop(fname, lims=(0, None)):
    """Read and crop."""
    return read_raw_fif(fname, allow_maxshield='yes').crop(*lims)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_movement_compensation(tmp_path):
    """Test movement compensation."""
    temp_dir = str(tmp_path)
    lims = (0, 4)
    raw = read_crop(raw_fname, lims).load_data()
    head_pos = read_head_pos(pos_fname)

    #
    # Movement compensation, no regularization, no tSSS
    #
    _assert_mag_coil_type(raw.info, FIFF.FIFFV_COIL_VV_MAG_T3)
    assert_allclose(raw.info['chs'][2]['cal'], 4.14e-11, rtol=1e-6)
    raw.info['chs'][2]['coil_type'] = FIFF.FIFFV_COIL_VV_MAG_T2
    raw_sss = maxwell_filter(raw, head_pos=head_pos, origin=mf_head_origin,
                             regularize=None, bad_condition='ignore')
    _assert_mag_coil_type(raw_sss.info, FIFF.FIFFV_COIL_VV_MAG_T3)
    assert_meg_snr(raw_sss, read_crop(sss_movecomp_fname, lims),
                   4.6, 12.4, chpi_med_tol=58)
    # IO
    temp_fname = op.join(temp_dir, 'test_raw_sss.fif')
    raw_sss.save(temp_fname)
    raw_sss = read_crop(temp_fname)
    assert_meg_snr(raw_sss, read_crop(sss_movecomp_fname, lims),
                   4.6, 12.4, chpi_med_tol=58)

    #
    # Movement compensation,    regularization, no tSSS
    #
    raw_sss = maxwell_filter(raw, head_pos=head_pos, origin=mf_head_origin)
    assert_meg_snr(raw_sss, read_crop(sss_movecomp_reg_in_fname, lims),
                   0.5, 1.9, chpi_med_tol=121)

    #
    # Movement compensation,    regularization,    tSSS at the end
    #
    raw_nohpi = filter_chpi(raw.copy(), t_window=0.2)
    with pytest.warns(RuntimeWarning, match='untested'):
        raw_sss_mv = maxwell_filter(raw_nohpi, head_pos=head_pos,
                                    st_duration=4., origin=mf_head_origin,
                                    st_fixed=False)
    # Neither match is particularly good because our algorithm actually differs
    assert_meg_snr(raw_sss_mv, read_crop(sss_movecomp_reg_in_st4s_fname, lims),
                   0.6, 1.3)
    tSSS_fname = op.join(sss_path, 'test_move_anon_st4s_raw_sss.fif')
    assert_meg_snr(raw_sss_mv, read_crop(tSSS_fname, lims),
                   0.6, 1.0, chpi_med_tol=None)
    assert_meg_snr(read_crop(sss_movecomp_reg_in_st4s_fname),
                   read_crop(tSSS_fname), 0.8, 1.0, chpi_med_tol=None)

    #
    # Movement compensation,    regularization,    tSSS at the beginning
    #
    raw_sss_mc = maxwell_filter(raw_nohpi, head_pos=head_pos, st_duration=4.,
                                origin=mf_head_origin)
    assert_meg_snr(raw_sss_mc, read_crop(tSSS_fname, lims),
                   0.6, 1.0, chpi_med_tol=None)
    assert_meg_snr(raw_sss_mc, raw_sss_mv, 0.6, 1.4)

    # some degenerate cases
    raw_erm = read_crop(erm_fname)
    with pytest.raises(ValueError, match='positions can only be used'):
        maxwell_filter(raw_erm, coord_frame='meg',
                       head_pos=head_pos)
    with pytest.raises(ValueError, match=r'of shape \(N, 10\)'):
        maxwell_filter(raw, head_pos=head_pos[:, :9])
    with pytest.raises(TypeError, match='instance of ndarray'):
        maxwell_filter(raw, head_pos='foo')
    with pytest.raises(ValueError, match='ascending'):
        maxwell_filter(raw, head_pos=head_pos[::-1])
    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 0] = raw._first_time - 1e-2
    with pytest.raises(ValueError, match='greater than'):
        maxwell_filter(raw, head_pos=head_pos_bad)

    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 4] = 1.  # off by more than 1 m
    with pytest.warns(RuntimeWarning, match='greater than 1 m'):
        maxwell_filter(raw.copy().crop(0, 0.1), head_pos=head_pos_bad,
                       bad_condition='ignore')

    # make sure numerical error doesn't screw it up, though
    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 0] = raw._first_time - 5e-4
    raw_sss_tweak = maxwell_filter(
        raw.copy().crop(0, 0.05), head_pos=head_pos_bad, origin=mf_head_origin)
    assert_meg_snr(raw_sss_tweak, raw_sss.copy().crop(0, 0.05), 1.4, 8.,
                   chpi_med_tol=5)


@pytest.mark.slowtest
def test_other_systems():
    """Test Maxwell filtering on KIT, BTI, and CTF files."""
    # KIT
    kit_dir = op.join(io_dir, 'kit', 'tests', 'data')
    sqd_path = op.join(kit_dir, 'test.sqd')
    mrk_path = op.join(kit_dir, 'test_mrk.sqd')
    elp_path = op.join(kit_dir, 'test_elp.txt')
    hsp_path = op.join(kit_dir, 'test_hsp.txt')
    raw_kit = read_raw_kit(sqd_path, mrk_path, elp_path, hsp_path)
    with pytest.warns(RuntimeWarning, match='fit'):
        pytest.raises(RuntimeError, maxwell_filter, raw_kit)
    with catch_logging() as log:
        raw_sss = maxwell_filter(raw_kit, origin=(0., 0., 0.04),
                                 ignore_ref=True, verbose=True)
    assert '12/15 out' in log.getvalue()  # homogeneous fields removed
    _assert_n_free(raw_sss, 65, 65)
    raw_sss_auto = maxwell_filter(raw_kit, origin=(0., 0., 0.04),
                                  ignore_ref=True, mag_scale='auto')
    assert_allclose(raw_sss._data, raw_sss_auto._data)
    # The KIT origin fit is terrible
    with pytest.warns(RuntimeWarning, match='more than 20 mm'):
        with catch_logging() as log:
            pytest.raises(RuntimeError, maxwell_filter, raw_kit,
                          ignore_ref=True, regularize=None)  # bad condition
            raw_sss = maxwell_filter(raw_kit, origin='auto',
                                     ignore_ref=True, bad_condition='info',
                                     verbose=True)
    log = log.getvalue()
    assert 'badly conditioned' in log
    assert 'more than 20 mm from' in log
    # fits can differ slightly based on scipy version, so be lenient here
    _assert_n_free(raw_sss, 28, 34)  # bad origin == brutal reg
    # Let's set the origin
    with catch_logging() as log:
        raw_sss = maxwell_filter(raw_kit, origin=(0., 0., 0.04),
                                 ignore_ref=True, bad_condition='info',
                                 regularize=None, verbose=True)
    log = log.getvalue()
    assert 'badly conditioned' in log
    assert '80/80 in, 12/15 out' in log
    _assert_n_free(raw_sss, 80)
    # Now with reg
    with catch_logging() as log:
        raw_sss = maxwell_filter(raw_kit, origin=(0., 0., 0.04),
                                 ignore_ref=True, verbose=True)
    log = log.getvalue()
    assert 'badly conditioned' not in log
    assert '12/15 out' in log
    _assert_n_free(raw_sss, 65)

    # BTi
    bti_dir = op.join(io_dir, 'bti', 'tests', 'data')
    bti_pdf = op.join(bti_dir, 'test_pdf_linux')
    bti_config = op.join(bti_dir, 'test_config_linux')
    bti_hs = op.join(bti_dir, 'test_hs_linux')
    raw_bti = read_raw_bti(bti_pdf, bti_config, bti_hs, preload=False)
    picks = pick_types(raw_bti.info, meg='mag', exclude=())
    power = np.sqrt(np.sum(raw_bti[picks][0] ** 2))
    raw_sss = maxwell_filter(raw_bti)
    _assert_n_free(raw_sss, 70)
    _assert_shielding(raw_sss, power, 0.5)
    raw_sss_auto = maxwell_filter(raw_bti, mag_scale='auto', verbose=True)
    _assert_shielding(raw_sss_auto, power, 0.7)

    # CTF
    raw_ctf = read_crop(fname_ctf_raw)
    assert raw_ctf.compensation_grade == 3
    with pytest.raises(RuntimeError, match='compensated'):
        maxwell_filter(raw_ctf)
    raw_ctf.apply_gradient_compensation(0)
    with pytest.raises(ValueError, match='digitization points'):
        maxwell_filter(raw_ctf)
    raw_sss = maxwell_filter(raw_ctf, origin=(0., 0., 0.04))
    _assert_n_free(raw_sss, 68)
    _assert_shielding(raw_sss, raw_ctf, 1.8)
    with catch_logging() as log:
        raw_sss = maxwell_filter(raw_ctf, origin=(0., 0., 0.04),
                                 ignore_ref=True, verbose=True)
    assert ', 12/15 out' in log.getvalue()  # homogeneous fields removed
    _assert_n_free(raw_sss, 70)
    _assert_shielding(raw_sss, raw_ctf, 12)
    raw_sss_auto = maxwell_filter(raw_ctf, origin=(0., 0., 0.04),
                                  ignore_ref=True, mag_scale='auto')
    assert_allclose(raw_sss._data, raw_sss_auto._data)
    with catch_logging() as log:
        maxwell_filter(raw_ctf, origin=(0., 0., 0.04), regularize=None,
                       ignore_ref=True, verbose=True)
    assert '80/80 in, 12/15 out' in log.getvalue()  # homogeneous fields


def test_spherical_conversions():
    """Test spherical harmonic conversions."""
    # Test our real<->complex conversion functions
    az, pol = np.meshgrid(np.linspace(0, 2 * np.pi, 30),
                          np.linspace(0, np.pi, 20))
    for degree in range(1, int_order):
        for order in range(0, degree + 1):
            sph = sph_harm(order, degree, az, pol)
            # ensure that we satisfy the conjugation property
            assert_allclose(_sh_negate(sph, order),
                            sph_harm(-order, degree, az, pol))
            # ensure our conversion functions work
            sph_real_pos = _sh_complex_to_real(sph, order)
            sph_real_neg = _sh_complex_to_real(sph, -order)
            sph_2 = _sh_real_to_complex([sph_real_pos, sph_real_neg], order)
            assert_allclose(sph, sph_2, atol=1e-7)


@testing.requires_testing_data
def test_multipolar_bases():
    """Test multipolar moment basis calculation using sensor information."""
    from scipy.io import loadmat
    # Test our basis calculations
    info = read_info(raw_fname)
    with use_coil_def(elekta_def_fname):
        coils = _prep_meg_channels(info, do_es=True)['defs']
    # Check against a known benchmark
    sss_data = loadmat(bases_fname)
    exp = dict(int_order=int_order, ext_order=ext_order)
    for origin in ((0, 0, 0.04), (0, 0.02, 0.02)):
        o_str = ''.join('%d' % (1000 * n) for n in origin)
        exp.update(origin=origin)
        S_tot = _sss_basis_basic(exp, coils, method='alternative')
        # Test our real<->complex conversion functions
        S_tot_complex = _bases_real_to_complex(S_tot, int_order, ext_order)
        S_tot_round = _bases_complex_to_real(S_tot_complex,
                                             int_order, ext_order)
        assert_allclose(S_tot, S_tot_round, atol=1e-7)

        S_tot_mat = np.concatenate([sss_data['Sin' + o_str],
                                    sss_data['Sout' + o_str]], axis=1)
        S_tot_mat_real = _bases_complex_to_real(S_tot_mat,
                                                int_order, ext_order)
        S_tot_mat_round = _bases_real_to_complex(S_tot_mat_real,
                                                 int_order, ext_order)
        assert_allclose(S_tot_mat, S_tot_mat_round, atol=1e-7)
        assert_allclose(S_tot_complex, S_tot_mat, rtol=1e-4, atol=1e-8)
        assert_allclose(S_tot, S_tot_mat_real, rtol=1e-4, atol=1e-8)

        # Now normalize our columns
        S_tot /= np.sqrt(np.sum(S_tot * S_tot, axis=0))[np.newaxis]
        S_tot_complex /= np.sqrt(np.sum(
            (S_tot_complex * S_tot_complex.conj()).real, axis=0))[np.newaxis]
        # Check against a known benchmark
        S_tot_mat = np.concatenate([sss_data['SNin' + o_str],
                                    sss_data['SNout' + o_str]], axis=1)
        # Check this roundtrip
        S_tot_mat_real = _bases_complex_to_real(S_tot_mat,
                                                int_order, ext_order)
        S_tot_mat_round = _bases_real_to_complex(S_tot_mat_real,
                                                 int_order, ext_order)
        assert_allclose(S_tot_mat, S_tot_mat_round, atol=1e-7)
        assert_allclose(S_tot_complex, S_tot_mat, rtol=1e-4, atol=1e-8)

        # Now test our optimized version
        S_tot = _sss_basis_basic(exp, coils)
        with use_coil_def(elekta_def_fname):
            S_tot_fast = _trans_sss_basis(
                exp, all_coils=_prep_mf_coils(info), trans=info['dev_head_t'])
        # there are some sign differences for columns (order/degrees)
        # in here, likely due to Condon-Shortley. Here we use a
        # Magnetometer channel to figure out the flips because the
        # gradiometer channels have effectively zero values for first three
        # external components (i.e., S_tot[grad_picks, 80:83])
        flips = (np.sign(S_tot_fast[2]) != np.sign(S_tot[2]))
        flips = 1 - 2 * flips
        assert_allclose(S_tot, S_tot_fast * flips, atol=1e-16)


@testing.requires_testing_data
def test_basic():
    """Test Maxwell filter basic version."""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    raw = read_crop(raw_fname, (0., 1.))
    raw_err = read_crop(raw_fname).apply_proj()
    raw_erm = read_crop(erm_fname)
    with pytest.raises(RuntimeError, match='cannot be applied'):
        maxwell_filter(raw_err)
    with pytest.raises(TypeError, match='instance of BaseRaw'):
        maxwell_filter(1.)
    with pytest.raises(ValueError, match='Number of requested bases'):
        maxwell_filter(raw, int_order=20)  # too many

    n_int_bases = int_order ** 2 + 2 * int_order
    n_ext_bases = ext_order ** 2 + 2 * ext_order
    nbases = n_int_bases + n_ext_bases

    # Check number of bases computed correctly
    assert _get_n_moments([int_order, ext_order]).sum() == nbases

    # Test SSS computation at the standard head origin
    assert len(raw.info['projs']) == 12  # 11 MEG projs + 1 AVG EEG
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore')
    assert len(raw_sss.info['projs']) == 1  # avg EEG
    assert raw_sss.info['projs'][0]['desc'] == 'Average EEG reference'
    assert_meg_snr(raw_sss, read_crop(sss_std_fname), 200., 1000.)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert len(py_cal) == 0
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert len(py_ctc) == 0
    py_st = raw_sss.info['proc_history'][0]['max_info']['max_st']
    assert len(py_st) == 0
    with pytest.raises(RuntimeError, match='cannot reapply'):
        maxwell_filter(raw_sss)

    # Test SSS computation at non-standard head origin
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, origin=[0., 0.02, 0.02], regularize=None,
                                 bad_condition='ignore')
    assert_meg_snr(raw_sss, read_crop(sss_nonstd_fname), 250., 700.)

    # Test SSS computation at device origin
    sss_erm_std = read_crop(sss_erm_std_fname)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             origin=mf_meg_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_erm_std, 70., 260.)
    for key in ('job', 'frame'):
        vals = [x.info['proc_history'][0]['max_info']['sss_info'][key]
                for x in [raw_sss, sss_erm_std]]
        assert vals[0] == vals[1]

    # Two equivalent things: at device origin in device coords (0., 0., 0.)
    # and at device origin at head coords info['dev_head_t'][:3, 3]
    raw_sss_meg = maxwell_filter(
        raw, coord_frame='meg', origin=(0., 0., 0.))
    raw_sss_head = maxwell_filter(
        raw, origin=raw.info['dev_head_t']['trans'][:3, 3])
    assert_meg_snr(raw_sss_meg, raw_sss_head, 100., 900.)

    # Check against SSS functions from proc_history
    assert _get_n_moments(int_order) == _get_rank_sss(raw_sss)

    # Degenerate cases
    with pytest.raises(ValueError, match='Invalid value'):
        maxwell_filter(raw, coord_frame='foo')
    with pytest.raises(ValueError, match='numerical array'):
        maxwell_filter(raw, origin='foo')
    with pytest.raises(ValueError, match='3-element array'):
        maxwell_filter(raw, origin=[0] * 4)
    with pytest.raises(ValueError, match='must be a float'):
        maxwell_filter(raw, mag_scale='foo')
    raw_missing = raw.copy().load_data()
    raw_missing.info['bads'] = ['MEG0111']
    raw_missing.pick_types(meg=True)  # will be missing the bad
    maxwell_filter(raw_missing)
    with pytest.warns(RuntimeWarning, match='not in data'):
        maxwell_filter(raw_missing, calibration=fine_cal_fname)


@testing.requires_testing_data
def test_maxwell_filter_additional(tmp_path):
    """Test processing of Maxwell filtered data."""
    # TODO: Future tests integrate with mne/io/tests/test_proc_history

    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    file_name = 'test_move_anon'
    raw_fname = op.join(data_path, 'SSS', file_name + '_raw.fif')

    # Use 2.0 seconds of data to get stable cov. estimate
    raw = read_crop(raw_fname, (0., 2.))

    # Get MEG channels, compute Maxwell filtered data
    raw.load_data()
    raw.pick_types(meg=True, eeg=False)
    int_order = 8
    raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')

    # Test io on processed data
    tempdir = str(tmp_path)
    test_outname = op.join(tempdir, 'test_raw_sss.fif')
    raw_sss.save(test_outname)
    raw_sss_loaded = read_crop(test_outname).load_data()

    # Some numerical imprecision since save uses 'single' fmt
    assert_allclose(raw_sss_loaded[:][0], raw_sss[:][0],
                    rtol=1e-6, atol=1e-20)

    # Test rank of covariance matrices for raw and SSS processed data
    cov_raw = compute_raw_covariance(raw)
    cov_sss = compute_raw_covariance(raw_sss)

    scalings = None
    cov_raw_rank = _compute_rank_int(
        cov_raw, scalings=scalings, info=raw.info, proj=False)
    cov_sss_rank = _compute_rank_int(
        cov_sss, scalings=scalings, info=raw_sss.info, proj=False)

    assert cov_raw_rank == raw.info['nchan']
    assert cov_sss_rank == _get_n_moments(int_order)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_bads_reconstruction():
    """Test Maxwell filter reconstruction of bad channels."""
    raw = read_crop(raw_fname, (0., 1.))
    raw.info['bads'] = bads
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore')
    assert_meg_snr(raw_sss, read_crop(sss_bad_recon_fname), 300.)


@pytest.mark.slowtest
@buggy_mkl_svd
@testing.requires_testing_data
def test_spatiotemporal():
    """Test Maxwell filter (tSSS) spatiotemporal processing."""
    # Load raw testing data
    raw = read_crop(raw_fname)

    # Test that window is less than length of data
    with pytest.raises(ValueError, match='must be'):
        maxwell_filter(raw, st_duration=1000.)

    # We could check both 4 and 10 seconds because Elekta handles them
    # differently (to ensure that std/non-std tSSS windows are correctly
    # handled), but the 4-sec case should hopefully be sufficient.
    st_durations = [4.]  # , 10.]
    tols = [(80, 100)]  # , 200.]
    kwargs = dict(origin=mf_head_origin, regularize=None,
                  bad_condition='ignore')
    for st_duration, tol in zip(st_durations, tols):
        # Load tSSS data depending on st_duration and get data
        tSSS_fname = op.join(sss_path,
                             'test_move_anon_st%0ds_raw_sss.fif' % st_duration)
        tsss_bench = read_crop(tSSS_fname)
        # Because Elekta's tSSS sometimes(!) lumps the tail window of data
        # onto the previous buffer if it's shorter than st_duration, we have to
        # crop the data here to compensate for Elekta's tSSS behavior.
        # if st_duration == 10.:
        #     tsss_bench.crop(0, st_duration)
        #     raw.crop(0, st_duration)

        # Test sss computation at the standard head origin. Same cropping issue
        # as mentioned above.
        raw_tsss = maxwell_filter(
            raw, st_duration=st_duration, **kwargs)
        assert _compute_rank_int(raw_tsss, proj=False) == 140
        assert_meg_snr(raw_tsss, tsss_bench, *tol)
        py_st = raw_tsss.info['proc_history'][0]['max_info']['max_st']
        assert (len(py_st) > 0)
        assert py_st['buflen'] == st_duration
        assert py_st['subspcorr'] == 0.98

    # Degenerate cases
    with pytest.raises(ValueError, match='Need 0 < st_correlation'):
        maxwell_filter(raw, st_duration=10., st_correlation=0.)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_spatiotemporal_only():
    """Test tSSS-only processing."""
    # Load raw testing data
    tmax = 0.5
    raw = read_crop(raw_fname, (0, tmax)).load_data()
    picks = pick_types(raw.info, meg=True, exclude='bads')[::2]
    raw.pick_channels([raw.ch_names[pick] for pick in picks])
    mag_picks = pick_types(raw.info, meg='mag', exclude=())
    power = np.sqrt(np.sum(raw[mag_picks][0] ** 2))
    # basics
    raw_tsss = maxwell_filter(raw, st_duration=tmax / 2., st_only=True)
    assert len(raw.info['projs']) == len(raw_tsss.info['projs'])
    assert _compute_rank_int(raw_tsss, proj=False) == len(picks)
    _assert_shielding(raw_tsss, power, 9)
    # with movement
    head_pos = read_head_pos(pos_fname)
    raw_tsss = maxwell_filter(raw, st_duration=tmax / 2., st_only=True,
                              head_pos=head_pos)
    assert _compute_rank_int(raw_tsss, proj=False) == len(picks)
    _assert_shielding(raw_tsss, power, 9)
    with pytest.warns(RuntimeWarning, match='st_fixed'):
        raw_tsss = maxwell_filter(raw, st_duration=tmax / 2., st_only=True,
                                  head_pos=head_pos, st_fixed=False)
    assert _compute_rank_int(raw_tsss, proj=False) == len(picks)
    _assert_shielding(raw_tsss, power, 9)
    # should do nothing
    raw_tsss = maxwell_filter(raw, st_duration=tmax, st_correlation=1.,
                              st_only=True)
    assert_allclose(raw[:][0], raw_tsss[:][0])
    # degenerate
    pytest.raises(ValueError, maxwell_filter, raw, st_only=True)  # no ST
    # two-step process equivalent to single-step process
    raw_tsss = maxwell_filter(raw, st_duration=tmax, st_only=True)
    raw_tsss = maxwell_filter(raw_tsss)
    raw_tsss_2 = maxwell_filter(raw, st_duration=tmax)
    assert_meg_snr(raw_tsss, raw_tsss_2, 1e5)
    # now also with head movement, and a bad MEG channel
    assert len(raw.info['bads']) == 0
    bads = [raw.ch_names[0]]
    raw.info['bads'] = list(bads)
    raw_tsss = maxwell_filter(raw, st_duration=tmax, st_only=True,
                              head_pos=head_pos)
    assert raw.info['bads'] == bads
    assert raw_tsss.info['bads'] == bads  # don't reset
    raw_tsss = maxwell_filter(raw_tsss, head_pos=head_pos)
    assert raw_tsss.info['bads'] == []  # do reset MEG bads
    raw_tsss_2 = maxwell_filter(raw, st_duration=tmax, head_pos=head_pos)
    assert raw_tsss_2.info['bads'] == []
    assert_meg_snr(raw_tsss, raw_tsss_2, 1e5)


@testing.requires_testing_data
def test_fine_calibration():
    """Test Maxwell filter fine calibration."""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    raw = read_crop(raw_fname, (0., 1.))
    sss_fine_cal = read_crop(sss_fine_cal_fname)

    # Test 1D SSS fine calibration
    with use_coil_def(elekta_def_fname):
        with catch_logging() as log:
            raw_sss = maxwell_filter(raw, calibration=fine_cal_fname,
                                     origin=mf_head_origin, regularize=None,
                                     bad_condition='ignore', verbose=True)
    log = log.getvalue()
    assert 'Using fine calibration' in log
    assert op.basename(fine_cal_fname) in log
    assert_meg_snr(raw_sss, sss_fine_cal, 82, 611)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert (py_cal is not None)
    assert (len(py_cal) > 0)
    mf_cal = sss_fine_cal.info['proc_history'][0]['max_info']['sss_cal']
    # we identify these differently
    mf_cal['cal_chans'][mf_cal['cal_chans'][:, 1] == 3022, 1] = 3024
    assert_allclose(py_cal['cal_chans'], mf_cal['cal_chans'])
    assert_allclose(py_cal['cal_corrs'], mf_cal['cal_corrs'],
                    rtol=1e-3, atol=1e-3)
    # with missing channels
    raw_missing = raw.copy().load_data()
    raw_missing.info['bads'] = ['MEG0111', 'MEG0943']  # 1 mag, 1 grad
    raw_missing.info._check_consistency()
    raw_sss_bad = maxwell_filter(
        raw_missing, calibration=fine_cal_fname, origin=mf_head_origin,
        regularize=None, bad_condition='ignore')
    raw_missing.pick_types(meg=True)  # actually remove bads
    raw_sss_bad.pick_channels(raw_missing.ch_names)  # remove them here, too
    with pytest.warns(RuntimeWarning, match='cal channels not in data'):
        raw_sss_missing = maxwell_filter(
            raw_missing, calibration=fine_cal_fname, origin=mf_head_origin,
            regularize=None, bad_condition='ignore')
    assert_meg_snr(raw_sss_missing, raw_sss_bad, 1000., 10000.)

    # Test 3D SSS fine calibration (no equivalent func in MaxFilter yet!)
    # very low SNR as proc differs, eventually we should add a better test
    raw_sss_3D = maxwell_filter(raw, calibration=fine_cal_fname_3d,
                                origin=mf_head_origin, regularize=None,
                                bad_condition='ignore')
    assert_meg_snr(raw_sss_3D, sss_fine_cal, 1.0, 6.)
    raw_ctf = read_crop(fname_ctf_raw).apply_gradient_compensation(0)
    pytest.raises(RuntimeError, maxwell_filter, raw_ctf, origin=(0., 0., 0.04),
                  calibration=fine_cal_fname)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_regularization():
    """Test Maxwell filter regularization."""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    min_tols = (20., 2.6, 1.0)
    med_tols = (200., 21., 3.7)
    origins = ((0., 0., 0.04), (0.,) * 3, (0., 0.02, 0.02))
    coord_frames = ('head', 'meg', 'head')
    raw_fnames = (raw_fname, erm_fname, sample_fname)
    sss_fnames = (sss_reg_in_fname, sss_erm_reg_in_fname,
                  sss_samp_reg_in_fname)
    comp_tols = [0, 1, 4]
    for ii, rf in enumerate(raw_fnames):
        raw = read_crop(rf, (0., 1.))
        sss_reg_in = read_crop(sss_fnames[ii])

        # Test "in" regularization
        raw_sss = maxwell_filter(raw, coord_frame=coord_frames[ii],
                                 origin=origins[ii])
        assert_meg_snr(raw_sss, sss_reg_in, min_tols[ii], med_tols[ii], msg=rf)

        # check components match
        _check_reg_match(raw_sss, sss_reg_in, comp_tols[ii])


def _check_reg_match(sss_py, sss_mf, comp_tol):
    """Check regularization."""
    info_py = sss_py.info['proc_history'][0]['max_info']['sss_info']
    assert (info_py is not None)
    assert (len(info_py) > 0)
    info_mf = sss_mf.info['proc_history'][0]['max_info']['sss_info']
    n_in = None
    for inf in (info_py, info_mf):
        if n_in is None:
            n_in = _get_n_moments(inf['in_order'])
        else:
            assert n_in == _get_n_moments(inf['in_order'])
        assert inf['components'][:n_in].sum() == inf['nfree']
    assert_allclose(info_py['nfree'], info_mf['nfree'],
                    atol=comp_tol, err_msg=sss_py._filenames[0])


@testing.requires_testing_data
def test_cross_talk(tmp_path):
    """Test Maxwell filter cross-talk cancellation."""
    raw = read_crop(raw_fname, (0., 1.))
    raw.info['bads'] = bads
    sss_ctc = read_crop(sss_ctc_fname)
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, cross_talk=pathlib.Path(ctc_fname),
                                 origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_ctc, 275.)
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert (len(py_ctc) > 0)
    with pytest.raises(TypeError, match='path-like'):
        maxwell_filter(raw, cross_talk=raw)
    pytest.raises(ValueError, maxwell_filter, raw, cross_talk=raw_fname)
    mf_ctc = sss_ctc.info['proc_history'][0]['max_info']['sss_ctc']
    del mf_ctc['block_id']  # we don't write this
    assert isinstance(py_ctc['decoupler'], sparse.csc_matrix)
    assert isinstance(mf_ctc['decoupler'], sparse.csc_matrix)
    assert_array_equal(py_ctc['decoupler'].toarray(),
                       mf_ctc['decoupler'].toarray())
    # I/O roundtrip
    tempdir = str(tmp_path)
    fname = op.join(tempdir, 'test_sss_raw.fif')
    sss_ctc.save(fname)
    sss_ctc_read = read_raw_fif(fname)
    mf_ctc_read = sss_ctc_read.info['proc_history'][0]['max_info']['sss_ctc']
    assert isinstance(mf_ctc_read['decoupler'], sparse.csc_matrix)
    assert_array_equal(mf_ctc_read['decoupler'].toarray(),
                       mf_ctc['decoupler'].toarray())
    assert object_diff(py_ctc, mf_ctc) == ''
    raw_ctf = read_crop(fname_ctf_raw).apply_gradient_compensation(0)
    raw_sss = maxwell_filter(raw_ctf, origin=(0., 0., 0.04))
    _assert_n_free(raw_sss, 68)
    raw_sss = maxwell_filter(raw_ctf, origin=(0., 0., 0.04), ignore_ref=True)
    _assert_n_free(raw_sss, 70)
    raw_missing = raw.copy().crop(0, 0.1).load_data().pick_channels(
        [raw.ch_names[pi] for pi in pick_types(raw.info, meg=True,
                                               exclude=())[3:]])
    with pytest.warns(RuntimeWarning, match='Not all cross-talk channels'):
        maxwell_filter(raw_missing, cross_talk=ctc_fname)
    # MEG channels not in cross-talk
    pytest.raises(RuntimeError, maxwell_filter, raw_ctf, origin=(0., 0., 0.04),
                  cross_talk=ctc_fname)


@testing.requires_testing_data
def test_head_translation():
    """Test Maxwell filter head translation."""
    raw = read_crop(raw_fname, (0., 1.))
    # First try with an unchanged destination
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, destination=raw_fname,
                                 origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore')
    assert_meg_snr(raw_sss, read_crop(sss_std_fname, (0., 1.)), 200.)
    # Now with default
    with use_coil_def(elekta_def_fname):
        with pytest.warns(RuntimeWarning, match='over 25 mm'):
            raw_sss = maxwell_filter(raw, destination=mf_head_origin,
                                     origin=mf_head_origin, regularize=None,
                                     bad_condition='ignore', verbose=True)
    assert_meg_snr(raw_sss, read_crop(sss_trans_default_fname), 125.)
    destination = np.eye(4)
    destination[2, 3] = 0.04
    assert_allclose(raw_sss.info['dev_head_t']['trans'], destination)
    # Now to sample's head pos
    with pytest.warns(RuntimeWarning, match='= 25.6 mm'):
        raw_sss = maxwell_filter(raw, destination=sample_fname,
                                 origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore', verbose=True)
    assert_meg_snr(raw_sss, read_crop(sss_trans_sample_fname), 13., 100.)
    assert_allclose(raw_sss.info['dev_head_t']['trans'],
                    read_info(sample_fname)['dev_head_t']['trans'])
    # Degenerate cases
    pytest.raises(RuntimeError, maxwell_filter, raw,
                  destination=mf_head_origin, coord_frame='meg')
    pytest.raises(ValueError, maxwell_filter, raw, destination=[0.] * 4)


# TODO: Eventually add simulation tests mirroring Taulu's original paper
# that calculates the localization error:
# http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=1495874

def _assert_shielding(raw_sss, erm_power, min_factor, max_factor=np.inf,
                      meg='mag'):
    """Assert a minimum shielding factor using empty-room power."""
    __tracebackhide__ = True
    picks = pick_types(raw_sss.info, meg=meg, ref_meg=False)
    if isinstance(erm_power, BaseRaw):
        picks_erm = pick_types(raw_sss.info, meg=meg, ref_meg=False)
        assert_allclose(picks, picks_erm)
        erm_power = np.sqrt((erm_power[picks_erm][0] ** 2).sum())
    sss_power = raw_sss[picks][0].ravel()
    sss_power = np.sqrt(np.sum(sss_power * sss_power))
    factor = erm_power / sss_power
    assert min_factor <= factor < max_factor, (
        'Shielding factor not %0.3f <= %0.3f < %0.3f'
        % (min_factor, factor, max_factor))


@buggy_mkl_svd
@testing.requires_testing_data
@pytest.mark.parametrize('regularize', ('in', None))
@pytest.mark.parametrize('bads', ([], ['MEG0111']))
def test_esss(regularize, bads):
    """Test extended-basis SSS."""
    # Make some fake "projectors" that actually contain external SSS bases
    raw_erm = read_crop(erm_fname).load_data().pick_types(meg=True)
    raw_erm.info['bads'] = bads
    proj_sss = mne.compute_proj_raw(raw_erm, meg='combined', verbose='error',
                                    n_mag=15, n_grad=15)
    good_info = pick_info(raw_erm.info, pick_types(raw_erm.info, meg=True))
    S_tot = _trans_sss_basis(
        dict(int_order=0, ext_order=3, origin=(0., 0., 0.)),
        all_coils=_prep_mf_coils(good_info), coil_scale=1., trans=None)
    assert S_tot.shape[-1] == len(proj_sss)
    for a, b in zip(proj_sss, S_tot.T):
        a['data']['data'][:] = b
    with catch_logging() as log:
        raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                                 regularize=regularize, verbose=True)
    log = log.getvalue()
    assert 'xtend' not in log
    with catch_logging() as log:
        raw_sss_2 = maxwell_filter(raw_erm, coord_frame='meg',
                                   regularize=regularize, ext_order=0,
                                   extended_proj=proj_sss, verbose=True)
    log = log.getvalue()
    assert 'Extending external SSS basis using 15 projection' in log
    assert_allclose(raw_sss_2._data, raw_sss._data, atol=1e-20)

    # This should work, as the projectors should be a superset
    raw_erm.info['bads'] = raw_erm.info['bads'] + ['MEG0112']
    maxwell_filter(raw_erm, coord_frame='meg', extended_proj=proj_sss)

    # Degenerate condititons
    proj_sss = proj_sss[:2]
    proj_sss[0]['data']['col_names'] = proj_sss[0]['data']['col_names'][:-1]
    with pytest.raises(ValueError, match='were missing'):
        maxwell_filter(raw_erm, coord_frame='meg', extended_proj=proj_sss)
    proj_sss[0] = 1.
    with pytest.raises(TypeError, match=r'extended_proj\[0\] must be an inst'):
        maxwell_filter(raw_erm, coord_frame='meg', extended_proj=proj_sss)
    with pytest.raises(TypeError, match='extended_proj must be an inst'):
        maxwell_filter(raw_erm, coord_frame='meg', extended_proj=1.)


@contextmanager
def get_n_projected():
    """Get the number of projected tSSS components from the log."""
    count = list()
    with use_log_level(True):
        with catch_logging() as log:
            yield count
    log = log.getvalue()
    assert 'Processing data using tSSS' in log
    log = log.splitlines()
    reg = re.compile(r'\s+Projecting\s+([0-9])+\s+intersecting tSSS .*')
    for line in log:
        m = reg.match(line)
        if m:
            count.append(int(m.group(1)))


@buggy_mkl_svd
@pytest.mark.slowtest
@testing.requires_testing_data
def test_shielding_factor(tmp_path):
    """Test Maxwell filter shielding factor using empty room."""
    raw_erm = read_crop(erm_fname).load_data().pick_types(meg=True)
    erm_power = raw_erm[pick_types(raw_erm.info, meg='mag')][0]
    erm_power = np.sqrt(np.sum(erm_power * erm_power))
    erm_power_grad = raw_erm[pick_types(raw_erm.info, meg='grad')][0]
    erm_power_grad = np.sqrt(np.sum(erm_power * erm_power))

    # Vanilla SSS (second value would be for meg=True instead of meg='mag')
    _assert_shielding(read_crop(sss_erm_std_fname), erm_power, 10)  # 1.5)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None)
    _assert_shielding(raw_sss, erm_power, 12, 13)  # 1.5)
    _assert_shielding(raw_sss, erm_power_grad, 0.45, 0.55, 'grad')  # 1.5)

    # No external basis
    raw_sss_0 = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                               ext_order=0)
    _assert_shielding(raw_sss_0, erm_power, 1.0, 1.1)
    del raw_sss_0

    # Regularization
    _assert_shielding(read_crop(sss_erm_std_fname), erm_power, 10)  # 1.5)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg')
    _assert_shielding(raw_sss, erm_power, 14.5, 15.5)

    #
    # Extended (eSSS)
    #

    # Show that using empty-room projectors increase shielding factor
    proj = mne.compute_proj_raw(raw_erm, meg='combined', verbose='error',
                                n_mag=15, n_grad=15)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 38, 39)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             extended_proj=proj)
    _assert_shielding(raw_sss, erm_power, 49, 51)
    # Now with reg
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 42, 44)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             extended_proj=proj)
    _assert_shielding(raw_sss, erm_power, 59, 61)

    #
    # Different mag_scale values
    #
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             mag_scale='auto')
    _assert_shielding(raw_sss, erm_power, 12, 13)
    _assert_shielding(raw_sss, erm_power_grad, 0.48, 0.58, 'grad')
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             mag_scale=1.)  # not a good choice
    _assert_shielding(raw_sss, erm_power, 7.3, 8.)
    _assert_shielding(raw_sss, erm_power_grad, 0.2, 0.3, 'grad')
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             mag_scale=1000., bad_condition='ignore')
    _assert_shielding(raw_sss, erm_power, 4.0, 5.0)
    _assert_shielding(raw_sss, erm_power_grad, 0.1, 0.2, 'grad')

    #
    # Fine cal
    #
    _assert_shielding(read_crop(sss_erm_fine_cal_fname), erm_power, 12)  # 2.0)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             origin=mf_meg_origin,
                             calibration=pathlib.Path(fine_cal_fname))
    _assert_shielding(raw_sss, erm_power, 12, 13)  # 2.0)

    #
    # Crosstalk
    #
    _assert_shielding(read_crop(sss_erm_ctc_fname), erm_power, 12)  # 2.1)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             origin=mf_meg_origin,
                             cross_talk=ctc_fname)
    _assert_shielding(raw_sss, erm_power, 12, 13)  # 2.1)

    # Fine cal + Crosstalk
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             calibration=fine_cal_fname,
                             origin=mf_meg_origin,
                             cross_talk=ctc_fname)
    _assert_shielding(raw_sss, erm_power, 13, 14)  # 2.2)
    # Fine cal + Crosstalk + Extended
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             calibration=fine_cal_fname,
                             origin=mf_meg_origin,
                             cross_talk=ctc_fname, extended_proj=proj)
    _assert_shielding(raw_sss, erm_power, 28, 30)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             calibration=fine_cal_fname,
                             origin=mf_meg_origin,
                             cross_talk=ctc_fname, extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 25, 27)

    # tSSS
    _assert_shielding(read_crop(sss_erm_st_fname), erm_power, 37)  # 5.8)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             origin=mf_meg_origin, st_duration=1.)
    _assert_shielding(raw_sss, erm_power, 37, 38)  # 5.8)

    # Crosstalk + tSSS
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                                 cross_talk=ctc_fname, origin=mf_meg_origin,
                                 st_duration=1.)
    _assert_shielding(raw_sss, erm_power, 38, 39)  # 5.91)
    assert counts[0] == 4

    # Fine cal + tSSS
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                                 calibration=fine_cal_fname,
                                 origin=mf_meg_origin, st_duration=1.)
    _assert_shielding(raw_sss, erm_power, 38, 39)  # 5.98)
    assert counts[0] == 4

    # Extended + tSSS
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                                 origin=mf_meg_origin, st_duration=1.,
                                 extended_proj=proj)
    _assert_shielding(raw_sss, erm_power, 40, 42)
    assert counts[0] == 0
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                                 origin=mf_meg_origin, st_duration=1.,
                                 extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 35, 37)
    assert counts[0] == 1

    # Fine cal + Crosstalk + tSSS
    _assert_shielding(read_crop(sss_erm_st1FineCalCrossTalk_fname),
                      erm_power, 39, 40)  # 6.07)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             calibration=fine_cal_fname, origin=mf_meg_origin,
                             cross_talk=ctc_fname, st_duration=1.)
    _assert_shielding(raw_sss, erm_power, 39, 40)  # 6.05)

    # Fine cal + Crosstalk + tSSS + Extended (a bit worse)
    _assert_shielding(read_crop(sss_erm_st1FineCalCrossTalk_fname),
                      erm_power, 39, 40)  # 6.07)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg', regularize=None,
                             calibration=fine_cal_fname, origin=mf_meg_origin,
                             cross_talk=ctc_fname, st_duration=1.,
                             extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 34, 36)

    # Fine cal + Crosstalk + tSSS + Reg-in
    _assert_shielding(read_crop(sss_erm_st1FineCalCrossTalkRegIn_fname),
                      erm_power, 57, 58)  # 6.97)
    raw_sss = maxwell_filter(raw_erm, calibration=fine_cal_fname,
                             cross_talk=ctc_fname, st_duration=1.,
                             origin=mf_meg_origin,
                             coord_frame='meg', regularize='in')
    _assert_shielding(raw_sss, erm_power, 53, 54)  # 6.64)
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, calibration=fine_cal_fname,
                                 cross_talk=ctc_fname, st_duration=1.,
                                 coord_frame='meg', regularize='in')
    _assert_shielding(raw_sss, erm_power, 58, 59)  # 7.0)
    _assert_shielding(raw_sss, erm_power_grad, 1.6, 1.7, 'grad')
    assert counts[0] == 4
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, calibration=fine_cal_fname,
                                 cross_talk=ctc_fname, st_duration=1.,
                                 coord_frame='meg', regularize='in',
                                 mag_scale='auto')
    _assert_shielding(raw_sss, erm_power, 51, 52)
    _assert_shielding(raw_sss, erm_power_grad, 1.5, 1.6, 'grad')
    assert counts[0] == 3
    with get_n_projected() as counts:
        with _record_warnings():  # SVD convergence on arm64
            raw_sss = maxwell_filter(raw_erm, calibration=fine_cal_fname_3d,
                                     cross_talk=ctc_fname, st_duration=1.,
                                     coord_frame='meg', regularize='in')
    # Our 3D cal has worse defaults for this ERM than the 1D file
    _assert_shielding(raw_sss, erm_power, 57, 58)
    assert counts[0] == 3
    # Show it by rewriting the 3D as 1D and testing it
    temp_dir = str(tmp_path)
    temp_fname = op.join(temp_dir, 'test_cal.dat')
    with open(fine_cal_fname_3d, 'r') as fid:
        with open(temp_fname, 'w') as fid_out:
            for line in fid:
                fid_out.write(' '.join(line.strip().split(' ')[:14]) + '\n')
    with get_n_projected() as counts:
        with _record_warnings():  # SVD convergence sometimes
            raw_sss = maxwell_filter(raw_erm, calibration=temp_fname,
                                     cross_talk=ctc_fname, st_duration=1.,
                                     coord_frame='meg', regularize='in')
    # Our 3D cal has worse defaults for this ERM than the 1D file
    _assert_shielding(raw_sss, erm_power, 44, 45)
    assert counts[0] == 3

    # Fine cal + Crosstalk + tSSS + Reg-in + Extended
    with get_n_projected() as counts:
        raw_sss = maxwell_filter(raw_erm, calibration=fine_cal_fname,
                                 cross_talk=ctc_fname, st_duration=1.,
                                 coord_frame='meg', regularize='in',
                                 extended_proj=proj[:3])
    _assert_shielding(raw_sss, erm_power, 48, 50)
    assert counts[0] == 1


@pytest.mark.slowtest
@testing.requires_testing_data
def test_all():
    """Test maxwell filter using all options."""
    raw_fnames = (raw_fname, raw_fname, erm_fname, sample_fname)
    sss_fnames = (sss_st1FineCalCrossTalkRegIn_fname,
                  sss_st1FineCalCrossTalkRegInTransSample_fname,
                  sss_erm_st1FineCalCrossTalkRegIn_fname,
                  sss_samp_fname)
    fine_cals = (fine_cal_fname,
                 fine_cal_fname,
                 fine_cal_fname,
                 fine_cal_mgh_fname)
    coord_frames = ('head', 'head', 'meg', 'head')
    ctcs = (ctc_fname, ctc_fname, ctc_fname, ctc_mgh_fname)
    mins = (3.5, 3.5, 1.2, 0.9)
    meds = (10.8, 10.4, 3.2, 6.)
    st_durs = (1., 1., 1., None)
    destinations = (None, sample_fname, None, None)
    origins = (mf_head_origin,
               mf_head_origin,
               mf_meg_origin,
               mf_head_origin)
    for ii, rf in enumerate(raw_fnames):
        raw = read_crop(rf, (0., 1.))
        with _record_warnings():  # sometimes the fit is bad
            sss_py = maxwell_filter(
                raw, calibration=fine_cals[ii], cross_talk=ctcs[ii],
                st_duration=st_durs[ii], coord_frame=coord_frames[ii],
                destination=destinations[ii], origin=origins[ii])
        sss_mf = read_crop(sss_fnames[ii])
        assert_meg_snr(sss_py, sss_mf, mins[ii], meds[ii], msg=rf)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_triux():
    """Test TRIUX system support."""
    raw = read_crop(tri_fname, (0, 0.999))
    _assert_mag_coil_type(raw.info, FIFF.FIFFV_COIL_VV_MAG_T1)
    assert_allclose(raw.info['chs'][2]['cal'], 1.33e-10, rtol=1e-6)
    # standard
    with use_coil_def(elekta_def_fname):
        sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None)
    _assert_mag_coil_type(sss_py.info, FIFF.FIFFV_COIL_VV_MAG_T3)
    assert_meg_snr(sss_py, read_crop(tri_sss_fname), 37, 700)
    # cross-talk
    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
                            cross_talk=tri_ctc_fname)
    assert_meg_snr(sss_py, read_crop(tri_sss_ctc_fname), 31, 250)
    # fine cal
    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
                            calibration=tri_cal_fname)
    assert_meg_snr(sss_py, read_crop(tri_sss_cal_fname), 22, 200)
    # ctc+cal
    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
                            calibration=tri_cal_fname,
                            cross_talk=tri_ctc_fname)
    assert_meg_snr(sss_py, read_crop(tri_sss_ctc_cal_fname), 28, 200)
    # regularization
    sss_py = maxwell_filter(raw, coord_frame='meg', regularize='in')
    sss_mf = read_crop(tri_sss_reg_fname)
    assert_meg_snr(sss_py, sss_mf, 0.6, 9)
    _check_reg_match(sss_py, sss_mf, 1)
    # all three
    sss_py = maxwell_filter(raw, coord_frame='meg', regularize='in',
                            calibration=tri_cal_fname,
                            cross_talk=tri_ctc_fname)
    sss_mf = read_crop(tri_sss_ctc_cal_reg_in_fname)
    assert_meg_snr(sss_py, sss_mf, 0.6, 9)
    _check_reg_match(sss_py, sss_mf, 1)
    # tSSS
    raw = read_crop(tri_fname).fix_mag_coil_types()
    with use_coil_def(elekta_def_fname):
        sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
                                st_duration=4., verbose=True)
    assert_meg_snr(sss_py, read_crop(tri_sss_st4_fname), 700., 1600)


@testing.requires_testing_data
def test_MGH_cross_talk():
    """Test cross-talk."""
    raw = read_crop(raw_fname, (0., 1.))
    raw_sss = maxwell_filter(raw, cross_talk=ctc_mgh_fname)
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert (len(py_ctc) > 0)


@testing.requires_testing_data
def test_mf_skips():
    """Test processing of data with skips."""
    raw = read_raw_fif(skip_fname, preload=True)
    raw.fix_mag_coil_types()
    raw.pick_channels(raw.ch_names[:50])  # fast and inaccurate
    kwargs = dict(st_only=True, coord_frame='meg', int_order=4, ext_order=3)
    # smoke test that this runs
    maxwell_filter(raw, st_duration=17., skip_by_annotation=(), **kwargs)
    # and this one, too, which will process some all-zero data
    maxwell_filter(raw, st_duration=2., skip_by_annotation=(), **kwargs)
    with pytest.raises(ValueError, match='duration'):
        # skips decrease acceptable duration
        maxwell_filter(raw, st_duration=17., **kwargs)
    onsets, ends = _annotations_starts_stops(
        raw, ('edge', 'bad_acq_skip'), invert=True)
    assert (ends - onsets).min() / raw.info['sfreq'] == 2.
    assert (ends - onsets).max() / raw.info['sfreq'] == 3.
    for st_duration in (2., 3.):
        raw_sss = maxwell_filter(raw, st_duration=st_duration, **kwargs)
        for start, stop in zip(onsets, ends):
            orig_data = raw[:, start:stop][0]
            new_data = raw_sss[:, start:stop][0]
            if (stop - start) / raw.info['sfreq'] >= st_duration:
                # Should be modified
                assert not np.allclose(new_data, orig_data, atol=1e-20)
            else:
                # Should not be modified
                assert_allclose(new_data, orig_data, atol=1e-20)
    # Processing an individual file and concat should be equivalent to
    # concat then process
    raw.crop(0, 1)
    raw_sss = maxwell_filter(raw, st_duration=1., **kwargs)
    raw_sss_concat = concatenate_raws([raw_sss, raw_sss.copy()])
    raw_concat = concatenate_raws([raw.copy(), raw.copy()])
    raw_concat_sss = maxwell_filter(raw_concat, st_duration=1., **kwargs)
    raw_concat_sss_bad = maxwell_filter(raw_concat, st_duration=1.,
                                        skip_by_annotation=(), **kwargs)
    data_c = raw_concat[:][0]
    data_sc = raw_sss_concat[:][0]
    data_cs = raw_concat_sss[:][0]
    data_csb = raw_concat_sss_bad[:][0]
    assert not np.allclose(data_cs, data_c, atol=1e-20)
    assert not np.allclose(data_cs, data_csb, atol=1e-20)
    assert_allclose(data_sc, data_cs, atol=1e-20)


@testing.requires_testing_data
@pytest.mark.parametrize(
    ('fname', 'bads', 'annot', 'add_ch', 'ignore_ref', 'want_bads',
     'return_scores', 'h_freq', 'meas_date'), [
        # Neuromag data tested against MF
        (sample_fname, [], False, False, False, ['MEG 2443'], False, None,
         'raw'),
        # add 0111 to test picking, add annot to test it, and prepend chs for
        # idx
        (sample_fname, ['MEG 0111'], True, True, False, ['MEG 2443'], False,
         None, 'raw'),
        # CTF data seems to be sensitive to linalg lib (?) because some
        # channels are very close to the limit, so we just check that one shows
        # up
        (ctf_fname_continuous, [], False, False, False, {'BR1-4304'}, False,
         None, 'raw'),
        # faked
        (ctf_fname_continuous, [], False, False, True, ['MLC24-4304'], False,
         None, 'raw'),
        # For `return_scores=True`
        (sample_fname, ['MEG 0111'], True, True, False, ['MEG 2443'], True,
         50, 'raw'),
        (sample_fname, ['MEG 0111'], True, True, False, ['MEG 2443'], True,
         50, None),
    ])
def test_find_bad_channels_maxwell(fname, bads, annot, add_ch, ignore_ref,
                                   want_bads, return_scores, h_freq,
                                   meas_date):
    """Test automatic bad channel detection."""
    if fname.endswith('.ds'):
        raw = read_raw_ctf(fname).load_data()
        flat_idx = 33
    else:
        raw = read_raw_fif(fname)
        raw.fix_mag_coil_types().load_data().pick_types(meg=True, exclude=())
        flat_idx = 1
    if meas_date is None:
        raw.set_meas_date(None)
    else:
        assert meas_date == 'raw'
    raw.filter(None, 40)
    raw.info['bads'] = bads
    raw._data[flat_idx] = 0  # MaxFilter didn't have this but doesn't affect it
    want_flats = [raw.ch_names[flat_idx]]
    raw.apply_gradient_compensation(0)

    min_count = 5

    if add_ch:
        raw_eeg = read_raw_fif(fname)
        raw_eeg.pick_types(meg=False, eeg=True, exclude=()).load_data()
        with raw_eeg.info._unlock():
            raw_eeg.info['lowpass'] = 40.
        raw = raw_eeg.add_channels([raw])  # prepend the EEG channels
        assert 0 in pick_types(raw.info, meg=False, eeg=True)
    if ignore_ref:
        # Fake a bad one, otherwise we don't find any
        assert 42 in pick_types(raw.info, meg=True, ref_meg=False)
        assert raw.ch_names[42:43] == want_bads
        raw._data[42] += np.random.RandomState(0).randn(len(raw.times))
    # maxfilter -autobad on -v -f test_raw.fif -force -cal off -ctc off -regularize off -list -o test_raw.fif -f ~/mne_data/MNE-testing-data/MEG/sample/sample_audvis_trunc_raw.fif  # noqa: E501
    if annot:
        # do a problematic one (gh-7741): exactly one "step" unit
        step = int(round(raw.info['sfreq'] * 5.))
        dt = 1. / raw.info['sfreq']
        assert step == 1502
        raw.annotations.append(step * dt + raw._first_time, dt, 'BAD')
    with catch_logging() as log:
        return_vals = find_bad_channels_maxwell(
            raw, origin=(0., 0., 0.04), regularize=None,
            bad_condition='ignore', skip_by_annotation='BAD', verbose=True,
            ignore_ref=ignore_ref, min_count=min_count,
            return_scores=return_scores, h_freq=h_freq)

    if return_scores:
        assert len(return_vals) == 3
        got_bads, got_flats, got_scores = return_vals
    else:
        assert len(return_vals) == 2
        got_bads, got_flats = return_vals

    if isinstance(want_bads, list):
        assert got_bads == want_bads  # from MaxFilter
    else:
        assert want_bads.intersection(set(got_bads))
    assert got_flats == want_flats
    log = log.getvalue()
    assert 'Interval   1:    0.00' in log
    assert 'Interval   2:    5.00' in log

    if h_freq is not None and h_freq > raw.info['lowpass']:
        assert 'data has already been low-pass filtered' in log

    if return_scores:
        meg_chs = raw.copy().pick_types(meg=True, exclude=[]).ch_names
        ch_types = raw.get_channel_types(meg_chs)

        assert list(got_scores['ch_names']) == meg_chs
        assert list(got_scores['ch_types']) == ch_types
        # Check that time is monotonically increasing.
        assert (np.diff(got_scores['bins'].flatten()) >= 0).all()

        assert (got_scores['scores_flat'].shape ==
                got_scores['scores_noisy'].shape ==
                (len(meg_chs), len(got_scores['bins'])))

        assert (got_scores['limits_flat'].shape ==
                got_scores['limits_noisy'].shape ==
                (len(meg_chs), 1))

        # Check "flat" scores.
        scores_flat = got_scores['scores_flat']
        limits_flat = got_scores['limits_flat']
        # Deal with NaN's in the scores (can't use np.less directly because of
        # https://github.com/numpy/numpy/issues/17198)
        scores_flat[np.isnan(scores_flat)] = np.inf
        limits_flat[np.isnan(limits_flat)] = -np.inf
        n_segments_below_limit = (scores_flat < limits_flat).sum(-1)
        ch_idx = np.where(n_segments_below_limit >=
                          min(min_count, len(got_scores['bins'])))
        flats = set(got_scores['ch_names'][ch_idx])
        assert flats == set(want_flats)

        # Check "noisy" scores.
        scores_noisy = got_scores['scores_noisy']
        limits_noisy = got_scores['limits_noisy']
        scores_noisy[np.isnan(scores_noisy)] = -np.inf
        limits_noisy[np.isnan(limits_noisy)] = np.inf
        n_segments_beyond_limit = (scores_noisy > limits_noisy).sum(-1)
        ch_idx = np.where(n_segments_beyond_limit >=
                          min(min_count, len(got_scores['bins'])))
        bads = set(got_scores['ch_names'][ch_idx])
        assert bads == set(want_bads)


def test_find_bads_maxwell_flat():
    """Test find_bads_maxwell when there are flat channels."""
    # See gh-9479
    raw = mne.io.read_raw_fif(raw_small_fname).load_data()
    assert_allclose(raw.times[-1], 23.97, atol=1e-2)
    noisy, flat = find_bad_channels_maxwell(raw, min_count=1)
    assert noisy == ['MEG 1032', 'MEG 2313', 'MEG 2443']
    assert flat == []
    n = int(round(raw.info['sfreq'] * 10))
    assert (len(raw.times) - n) / raw.info['sfreq'] > 10  # at least 10 sec
    with catch_logging() as log:
        want_noisy, want_flat = find_bad_channels_maxwell(
            raw.copy().crop(n / raw.info['sfreq'], None), min_count=1,
            verbose='debug')
    log = log.getvalue()
    assert 'in 2 intervals ' in log
    assert want_noisy == ['MEG 2313', 'MEG 2443']
    assert want_flat == []
    raw._data[:, :n] = 0
    with pytest.warns(RuntimeWarning, match='All-flat segment detected'):
        with catch_logging() as log:
            noisy, flat = find_bad_channels_maxwell(
                raw, min_count=1, verbose='debug')
    log = log.getvalue()
    assert ' in 4 intervals ' in log
    assert flat == raw.ch_names[:306]
    assert noisy == []  # none found because all flat
    # now do what we suggest in the warning
    annot, _ = annotate_amplitude(raw, flat=0., bad_percent=100,
                                  min_duration=1.)
    assert_allclose(annot.duration, 10., atol=1e-2)  # not even divisor sfreq
    raw.info['bads'] = []
    raw.set_annotations(annot)
    data_good = raw.get_data(reject_by_annotation='omit')
    assert data_good.shape[1] / raw.info['sfreq'] / 5. > 2  # at least 10 sec
    with catch_logging() as log:
        noisy, flat = find_bad_channels_maxwell(
            raw, min_count=1, skip_by_annotation='bad_flat', verbose='debug')
    log = log.getvalue()
    assert ' in 2 intervals ' in log, log
    assert flat == want_flat
    assert noisy == want_noisy


@pytest.mark.parametrize('regularize, n', [
    (None, 80),
    ('in', 71),
])
def test_compute_maxwell_basis(regularize, n):
    """Test compute_maxwell_basis."""
    raw = read_raw_fif(raw_small_fname).crop(0, 2)
    assert raw.info['bads'] == []
    raw.del_proj()
    rank = compute_rank(raw)['meg']
    assert rank == 306
    raw.info['bads'] = ['MEG 2443']
    kwargs = dict(regularize=regularize, verbose=True)
    raw_sss = maxwell_filter(raw, **kwargs)
    want = raw_sss.get_data('meg')
    rank = compute_rank(raw_sss)['meg']
    assert rank == n
    S, pS, reg_moments, n_use_in = compute_maxwell_basis(raw.info, **kwargs)
    assert n_use_in == n
    assert n_use_in == len(reg_moments) - 15  # no externals removed
    xform = S[:, :n_use_in] @ pS[:n_use_in]
    got = xform @ raw.pick_types(meg=True, exclude='bads').get_data()
    assert_allclose(got, want)


@testing.requires_testing_data
@pytest.mark.parametrize('bads', ('from_raw', 'union', 'keep'))
def test_prepare_emptyroom_bads(bads):
    """Test prepare_emptyroom."""
    raw = read_raw_fif(raw_fname, allow_maxshield='yes', verbose=False)
    names = [name for name in raw.ch_names if 'EEG' not in name]
    raw.pick_channels(names)
    raw_er = read_raw_fif(erm_fname, allow_maxshield='yes', verbose=False)
    raw_er.pick_channels(names)
    assert raw.ch_names == raw_er.ch_names
    assert raw_er.info['dev_head_t'] is None
    assert raw.info['dev_head_t'] is not None
    raw_er.set_montage(None)

    if bads == 'from_raw':
        raw_bads_orig = ['MEG0113', 'MEG2313']
        raw_er_bads_orig = []
    elif bads == 'union':
        raw_bads_orig = ['MEG0113']
        raw_er_bads_orig = ['MEG2313']
    elif bads == 'keep':
        raw_bads_orig = []
        raw_er_bads_orig = ['MEG0113', 'MEG2313']

    raw.info['bads'] = raw_bads_orig
    raw_er.info['bads'] = raw_er_bads_orig

    raw_er_prepared = maxwell_filter_prepare_emptyroom(
        raw_er=raw_er,
        raw=raw,
        bads=bads
    )
    assert raw_er_prepared.info['bads'] == ['MEG0113', 'MEG2313']
    assert raw_er_prepared.info['dev_head_t'] == raw.info['dev_head_t']

    montage_expected = raw.copy().pick_types(meg=True).get_montage()
    assert raw_er_prepared.get_montage() == montage_expected

    # Ensure the originals were not modified
    assert raw.info['bads'] == raw_bads_orig
    assert raw_er.info['bads'] == raw_er_bads_orig
    assert raw_er.info['dev_head_t'] is None
    assert raw_er.get_montage() is None


@testing.requires_testing_data
@pytest.mark.slowtest  # lots of params
@pytest.mark.parametrize('set_annot_when', ('before', 'after'))
@pytest.mark.parametrize('raw_meas_date', ('orig', None))
@pytest.mark.parametrize('raw_er_meas_date', ('orig', None))
def test_prepare_emptyroom_annot_first_samp(set_annot_when, raw_meas_date,
                                            raw_er_meas_date):
    """Test prepare_emptyroom."""
    raw = read_raw_fif(raw_fname, allow_maxshield='yes', verbose=False)
    raw_er = read_raw_fif(erm_fname, allow_maxshield='yes', verbose=False)
    names = raw.ch_names[:3]  # make it faster
    raw.pick_channels(names)
    raw_er.pick_channels(names)
    assert raw.ch_names == raw_er.ch_names
    assert raw.info['meas_date'] != raw_er.info['meas_date']
    if raw_meas_date is None:
        raw.set_meas_date(None)
    if raw_er_meas_date is None:
        raw_er.set_meas_date(None)
    # to make life easier, make it the same duration
    n_rep = max(int(np.ceil(len(raw.times) / len(raw_er.times))), 1)
    raw_er = mne.concatenate_raws([raw_er] * n_rep).crop(0, raw.times[-1])
    assert_allclose(raw.times, raw_er.times)
    raw_er_first_samp_orig = raw_er.first_samp
    assert len(raw.annotations) == 0
    pos = mne.chpi.read_head_pos(pos_fname)
    annot, _ = annotate_movement(raw, pos, 1.)
    # Add an annotation right at the beginning and end to make sure nothing
    # gets cropped
    onset = raw.times[[0, -1]]
    duration = 1. / raw.info['sfreq']
    annot.append(onset + raw.first_time * (raw.info['meas_date'] is not None),
                 duration, ['BAD_CUSTOM'])
    want_annot = 7  # 5 from annotate_movement plus our first and last samps
    if set_annot_when == 'before':
        raw.set_annotations(annot)
        meas_date = 'keep'
        want_date = raw_er.info['meas_date']
    else:
        assert set_annot_when == 'after'
        meas_date = 'from_raw'
        want_date = raw.info['meas_date']
    raw_er_prepared = maxwell_filter_prepare_emptyroom(
        raw_er=raw_er, raw=raw, meas_date=meas_date, emit_warning=True)
    assert raw_er.first_samp == raw_er_first_samp_orig
    assert raw_er_prepared.info['meas_date'] == want_date
    assert raw_er_prepared.first_samp == raw.first_samp

    # Ensure (movement) annotations carry over regardless of whether they're
    # set before or after preparation
    assert len(annot) == want_annot
    if set_annot_when == 'after':
        raw.set_annotations(annot)
        raw_er_prepared.set_annotations(annot)
    assert len(raw.annotations) == want_annot
    prop_bad = np.isnan(
        raw.get_data([0], reject_by_annotation='nan')).mean()
    assert 0.3 < prop_bad < 0.4
    assert len(raw_er_prepared.annotations) == want_annot
    prop_bad_er = np.isnan(
        raw_er_prepared.get_data([0], reject_by_annotation='nan')).mean()
    assert_allclose(prop_bad, prop_bad_er)