def test_calculate_chpi_positions():
    """Test calculation of cHPI positions
    """
    trans, rot, t = head_pos_to_trans_rot_t(read_head_pos(pos_fname))
    with warnings.catch_warnings(record=True):
        raw = Raw(chpi_fif_fname, allow_maxshield=True, preload=True)
    t -= raw.first_samp / raw.info['sfreq']
    quats = _calculate_chpi_positions(raw, verbose='debug')
    trans_est, rot_est, t_est = head_pos_to_trans_rot_t(quats)
    _compare_positions((trans, rot, t), (trans_est, rot_est, t_est), 0.003)

    # degenerate conditions
    raw_no_chpi = Raw(test_fif_fname)
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_no_chpi)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = 999
            break
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['r'] = np.ones(3)
    raw_bad.crop(0, 1., copy=False)
    with warnings.catch_warnings(record=True):  # bad pos
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, verbose=True)
    # ignore HPI info header and [done] footer
    for line in log_file.getvalue().strip().split('\n')[4:-1]:
        assert_true('0/5 good' in line)

def test_calculate_chpi_positions():
    """Test calculation of cHPI positions."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(pos_fname)
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True)
    # This is a little hack (aliasing while decimating) to make it much faster
    # for testing purposes only. We can relax this later if we find it breaks
    # something.
    raw_dec = _decimate_chpi(raw, 15)
    with catch_logging() as log:
        py_quats = _calculate_chpi_positions(raw_dec, t_step_max=1.,
                                             verbose='debug')
    assert_true(log.getvalue().startswith('HPIFIT'))
    _assert_quats(py_quats, mf_quats, dist_tol=0.004, angle_tol=2.5)

    # degenerate conditions
    raw_no_chpi = read_raw_fif(test_fif_fname)
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_no_chpi)
    raw_bad = raw.copy()
    del raw_bad.info['hpi_meas'][0]['hpi_coils'][0]['coil_freq']
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = FIFF.FIFFV_COORD_UNKNOWN
            break
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = FIFF.FIFFV_COORD_HEAD
            d['r'] = np.ones(3)
    raw_bad.crop(0, 1.)
    picks = np.concatenate([np.arange(306, len(raw_bad.ch_names)),
                            pick_types(raw_bad.info, meg=True)[::16]])
    raw_bad.pick_channels([raw_bad.ch_names[pick] for pick in picks])
    with warnings.catch_warnings(record=True):  # bad pos
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, t_step_min=1., verbose=True)
    # ignore HPI info header and [done] footer
    assert_true('0/5 good' in log_file.getvalue().strip().split('\n')[-2])

    # half the rate cuts off cHPI coils
    raw.info['lowpass'] /= 2.
    assert_raises_regex(RuntimeError, 'above the',
                        _calculate_chpi_positions, raw)

    # test on 5k artemis data
    raw = read_raw_artemis123(art_fname, preload=True)
    mf_quats = read_head_pos(art_mc_fname)
    with catch_logging() as log:
        py_quats = _calculate_chpi_positions(raw, t_step_min=2.,
                                             verbose='debug')
    _assert_quats(py_quats, mf_quats, dist_tol=0.004, angle_tol=2.5)

def test_calculate_chpi_positions():
    """Test calculation of cHPI positions."""
    trans, rot, t = head_pos_to_trans_rot_t(read_head_pos(pos_fname))
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True,
                       add_eeg_ref=False)
    t -= raw.first_samp / raw.info['sfreq']
    quats = _calculate_chpi_positions(raw, verbose='debug')
    trans_est, rot_est, t_est = head_pos_to_trans_rot_t(quats)
    _compare_positions((trans, rot, t), (trans_est, rot_est, t_est), 0.003)

    # degenerate conditions
    raw_no_chpi = read_raw_fif(test_fif_fname, add_eeg_ref=False)
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_no_chpi)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = 999
            break
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['r'] = np.ones(3)
    raw_bad.crop(0, 1., copy=False)
    with warnings.catch_warnings(record=True):  # bad pos
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, verbose=True)
    # ignore HPI info header and [done] footer
    assert_true('0/5 good' in log_file.getvalue().strip().split('\n')[-2])

    # half the rate cuts off cHPI coils
    with warnings.catch_warnings(record=True):  # uint cast suggestion
        raw.resample(300., npad='auto')
    assert_raises_regex(RuntimeError, 'above the',
                        _calculate_chpi_positions, raw)

def test_calculate_chpi_positions():
    """Test calculation of cHPI positions
    """
    trans, rot, t = head_pos_to_trans_rot_t(read_head_pos(pos_fname))
    raw = Raw(chpi_fif_fname, allow_maxshield="yes", preload=True)
    t -= raw.first_samp / raw.info["sfreq"]
    quats = _calculate_chpi_positions(raw, verbose="debug")
    trans_est, rot_est, t_est = head_pos_to_trans_rot_t(quats)
    _compare_positions((trans, rot, t), (trans_est, rot_est, t_est), 0.003)

    # degenerate conditions
    raw_no_chpi = Raw(test_fif_fname)
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_no_chpi)
    raw_bad = raw.copy()
    for d in raw_bad.info["dig"]:
        if d["kind"] == FIFF.FIFFV_POINT_HPI:
            d["coord_frame"] = 999
            break
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info["dig"]:
        if d["kind"] == FIFF.FIFFV_POINT_HPI:
            d["r"] = np.ones(3)
    raw_bad.crop(0, 1.0, copy=False)
    with warnings.catch_warnings(record=True):  # bad pos
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, verbose=True)
    # ignore HPI info header and [done] footer
    for line in log_file.getvalue().strip().split("\n")[4:-1]:
        assert_true("0/5 good" in line)

    # half the rate cuts off cHPI coils
    with warnings.catch_warnings(record=True):  # uint cast suggestion
        raw.resample(300.0, npad="auto")
    assert_raises_regex(RuntimeError, "above the", _calculate_chpi_positions, raw)

def test_read_write_head_pos(tmpdir):
    """Test reading and writing head position quaternion parameters."""
    temp_name = op.join(str(tmpdir), 'temp.pos')
    # This isn't a 100% valid quat matrix but it should be okay for tests
    head_pos_rand = np.random.RandomState(0).randn(20, 10)
    # This one is valid
    head_pos_read = read_head_pos(pos_fname)
    for head_pos_orig in (head_pos_rand, head_pos_read):
        write_head_pos(temp_name, head_pos_orig)
        head_pos = read_head_pos(temp_name)
        assert_allclose(head_pos_orig, head_pos, atol=1e-3)
    # Degenerate cases
    pytest.raises(TypeError, write_head_pos, 0, head_pos_read)  # not filename
    pytest.raises(ValueError, write_head_pos, temp_name, 'foo')  # not array
    pytest.raises(ValueError, write_head_pos, temp_name, head_pos_read[:, :9])
    pytest.raises(TypeError, read_head_pos, 0)
    pytest.raises(IOError, read_head_pos, temp_name + 'foo')

def test_calculate_head_pos_ctf():
    """Test extracting of cHPI positions from ctf data."""
    raw = read_raw_ctf(ctf_chpi_fname)
    quats = _calculate_head_pos_ctf(raw)
    mc_quats = read_head_pos(ctf_chpi_pos_fname)
    _assert_quats(quats, mc_quats, dist_tol=0.004, angle_tol=2.5)

    raw = read_raw_fif(ctf_fname)
    pytest.raises(RuntimeError, _calculate_head_pos_ctf, raw)

def test_calculate_chpi_positions_on_chpi5_in_shorter_steps():
    """Comparing estimated cHPI positions with MF results (smaller steps)."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(chpi5_pos_fname)
    raw = read_raw_fif(chpi5_fif_fname, allow_maxshield='yes')
    raw = _decimate_chpi(raw.crop(0., 15.).load_data(), decim=8)
    py_quats = _calculate_chpi_positions(raw, t_step_min=0.1, t_step_max=0.1,
                                         t_window=0.1, verbose='debug')
    # needs interpolation, tolerance must be increased
    _assert_quats(py_quats, mf_quats, dist_tol=0.001, angle_tol=0.6)

def test_calculate_chpi_positions_on_chpi5_in_one_second_steps():
    """Comparing estimated cHPI positions with MF results (one second)."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(chpi5_pos_fname)
    raw = read_raw_fif(chpi5_fif_fname, allow_maxshield='yes')
    # the last two seconds contain a maxfilter problem!
    # fiff file timing: 26. to 43. seconds
    # maxfilter estimates a wrong head position for interval 16: 41.-42. sec
    raw = _decimate_chpi(raw.crop(0., 15.).load_data(), decim=8)
    # needs no interpolation, because maxfilter pos files comes with 1 s steps
    py_quats = _calculate_chpi_positions(raw, t_step_min=1.0, t_step_max=1.0,
                                         t_window=1.0, verbose='debug')
    _assert_quats(py_quats, mf_quats, dist_tol=0.0008, angle_tol=.5)

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_equal(len(raw.info['projs']), len(raw_tsss.info['projs']))
    assert_equal(raw_tsss.estimate_rank(), 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_equal(raw_tsss.estimate_rank(), 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_equal(raw_tsss.estimate_rank(), 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_equal(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_equal(raw.info['bads'], bads)
    assert_equal(raw_tsss.info['bads'], bads)  # don't reset
    raw_tsss = maxwell_filter(raw_tsss, head_pos=head_pos)
    assert_equal(raw_tsss.info['bads'], [])  # do reset MEG bads
    raw_tsss_2 = maxwell_filter(raw, st_duration=tmax, head_pos=head_pos)
    assert_equal(raw_tsss_2.info['bads'], [])
    assert_meg_snr(raw_tsss, raw_tsss_2, 1e5)

def test_spatiotemporal_only():
    """Test tSSS-only processing"""
    # Load raw testing data
    raw = Raw(raw_fname,
              allow_maxshield='yes').crop(0, 2, copy=False).load_data()
    picks = pick_types(raw.info, meg='mag', exclude=())
    power = np.sqrt(np.sum(raw[picks][0] ** 2))
    # basics
    raw_tsss = maxwell_filter(raw, st_duration=1., st_only=True)
    assert_equal(raw_tsss.estimate_rank(), 366)
    _assert_shielding(raw_tsss, power, 10)
    # temporal proj will actually reduce spatial DOF with small windows!
    raw_tsss = maxwell_filter(raw, st_duration=0.1, st_only=True)
    assert_true(raw_tsss.estimate_rank() < 350)
    _assert_shielding(raw_tsss, power, 40)
    # with movement
    head_pos = read_head_pos(pos_fname)
    raw_tsss = maxwell_filter(raw, st_duration=1., st_only=True,
                              head_pos=head_pos)
    assert_equal(raw_tsss.estimate_rank(), 366)
    _assert_shielding(raw_tsss, power, 12)
    with warnings.catch_warnings(record=True):  # st_fixed False
        raw_tsss = maxwell_filter(raw, st_duration=1., st_only=True,
                                  head_pos=head_pos, st_fixed=False)
    assert_equal(raw_tsss.estimate_rank(), 366)
    _assert_shielding(raw_tsss, power, 12)
    # should do nothing
    raw_tsss = maxwell_filter(raw, st_duration=1., st_correlation=1.,
                              st_only=True)
    assert_allclose(raw[:][0], raw_tsss[:][0])
    # degenerate
    assert_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=1., st_only=True)
    raw_tsss = maxwell_filter(raw_tsss)
    raw_tsss_2 = maxwell_filter(raw, st_duration=1.)
    assert_meg_snr(raw_tsss, raw_tsss_2, 1e5)
    # now also with head movement, and a bad MEG channel
    assert_equal(len(raw.info['bads']), 0)
    raw.info['bads'] = ['EEG001', 'MEG2623']
    raw_tsss = maxwell_filter(raw, st_duration=1., st_only=True,
                              head_pos=head_pos)
    assert_equal(raw.info['bads'], ['EEG001', 'MEG2623'])
    assert_equal(raw_tsss.info['bads'], ['EEG001', 'MEG2623'])  # don't reset
    raw_tsss = maxwell_filter(raw_tsss, head_pos=head_pos)
    assert_equal(raw_tsss.info['bads'], ['EEG001'])  # do reset MEG bads
    raw_tsss_2 = maxwell_filter(raw, st_duration=1., head_pos=head_pos)
    assert_equal(raw_tsss_2.info['bads'], ['EEG001'])
    assert_meg_snr(raw_tsss, raw_tsss_2, 1e5)

def test_simulate_raw_chpi():
    """Test simulation of raw data with cHPI."""
    raw = read_raw_fif(raw_chpi_fname, allow_maxshield='yes')
    picks = np.arange(len(raw.ch_names))
    picks = np.setdiff1d(picks, pick_types(raw.info, meg=True, eeg=True)[::4])
    raw.load_data().pick_channels([raw.ch_names[pick] for pick in picks])
    raw.info.normalize_proj()
    sphere = make_sphere_model('auto', 'auto', raw.info)
    # make sparse spherical source space
    sphere_vol = tuple(sphere['r0'] * 1000.) + (sphere.radius * 1000.,)
    src = setup_volume_source_space(sphere=sphere_vol, pos=70.)
    stc = _make_stc(raw, src)
    # simulate data with cHPI on
    with pytest.deprecated_call():
        raw_sim = simulate_raw(raw, stc, None, src, sphere, cov=None,
                               head_pos=pos_fname, interp='zero')
    # need to trim extra samples off this one
    with pytest.deprecated_call():
        raw_chpi = simulate_raw(raw, stc, None, src, sphere, cov=None,
                                chpi=True, head_pos=pos_fname, interp='zero')
    # test cHPI indication
    hpi_freqs, hpi_pick, hpi_ons = _get_hpi_info(raw.info)
    assert_allclose(raw_sim[hpi_pick][0], 0.)
    assert_allclose(raw_chpi[hpi_pick][0], hpi_ons.sum())
    # test that the cHPI signals make some reasonable values
    picks_meg = pick_types(raw.info, meg=True, eeg=False)
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)

    for picks in [picks_meg[:3], picks_eeg[:3]]:
        psd_sim, freqs_sim = psd_welch(raw_sim, picks=picks)
        psd_chpi, freqs_chpi = psd_welch(raw_chpi, picks=picks)

        assert_array_equal(freqs_sim, freqs_chpi)
        freq_idx = np.sort([np.argmin(np.abs(freqs_sim - f))
                            for f in hpi_freqs])
        if picks is picks_meg:
            assert (psd_chpi[:, freq_idx] >
                    100 * psd_sim[:, freq_idx]).all()
        else:
            assert_allclose(psd_sim, psd_chpi, atol=1e-20)

    # test localization based on cHPI information
    quats_sim = _calculate_chpi_positions(raw_chpi, t_step_min=10.)
    quats = read_head_pos(pos_fname)
    _assert_quats(quats, quats_sim, dist_tol=5e-3, angle_tol=3.5)

def test_simulate_raw_chpi():
    """Test simulation of raw data with cHPI"""
    with warnings.catch_warnings(record=True):  # MaxShield
        raw = Raw(raw_chpi_fname, allow_maxshield=True)
    sphere = make_sphere_model('auto', 'auto', raw.info)
    # make sparse spherical source space
    sphere_vol = tuple(sphere['r0'] * 1000.) + (sphere.radius * 1000.,)
    src = setup_volume_source_space('sample', sphere=sphere_vol, pos=70.)
    stc = _make_stc(raw, src)
    # simulate data with cHPI on
    raw_sim = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=False)
    # need to trim extra samples off this one
    raw_chpi = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=True,
                            head_pos=pos_fname)
    # test cHPI indication
    hpi_freqs, _, hpi_pick, hpi_ons = _get_hpi_info(raw.info)[:4]
    assert_allclose(raw_sim[hpi_pick][0], 0.)
    assert_allclose(raw_chpi[hpi_pick][0], hpi_ons.sum())
    # test that the cHPI signals make some reasonable values
    picks_meg = pick_types(raw.info, meg=True, eeg=False)
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)

    for picks in [picks_meg, picks_eeg]:
        psd_sim, freqs_sim = psd_welch(raw_sim, picks=picks)
        psd_chpi, freqs_chpi = psd_welch(raw_chpi, picks=picks)

        assert_array_equal(freqs_sim, freqs_chpi)
        freq_idx = np.sort([np.argmin(np.abs(freqs_sim - f))
                           for f in hpi_freqs])
        if picks is picks_meg:
            assert_true((psd_chpi[:, freq_idx] >
                         100 * psd_sim[:, freq_idx]).all())
        else:
            assert_allclose(psd_sim, psd_chpi, atol=1e-20)

    # test localization based on cHPI information
    quats_sim = _calculate_chpi_positions(raw_chpi)
    trans_sim, rot_sim, t_sim = head_pos_to_trans_rot_t(quats_sim)
    trans, rot, t = head_pos_to_trans_rot_t(read_head_pos(pos_fname))
    t -= raw.first_samp / raw.info['sfreq']
    _compare_positions((trans, rot, t), (trans_sim, rot_sim, t_sim),
                       max_dist=0.005)

def test_simulate_raw_chpi():
    """Test simulation of raw data with cHPI."""
    raw = read_raw_fif(raw_chpi_fname, allow_maxshield='yes')
    sphere = make_sphere_model('auto', 'auto', raw.info)
    # make sparse spherical source space
    sphere_vol = tuple(sphere['r0'] * 1000.) + (sphere.radius * 1000.,)
    src = setup_volume_source_space('sample', sphere=sphere_vol, pos=70.)
    stc = _make_stc(raw, src)
    # simulate data with cHPI on
    raw_sim = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=False)
    # need to trim extra samples off this one
    raw_chpi = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=True,
                            head_pos=pos_fname)
    # test cHPI indication
    hpi_freqs, _, hpi_pick, hpi_ons = _get_hpi_info(raw.info)[:4]
    assert_allclose(raw_sim[hpi_pick][0], 0.)
    assert_allclose(raw_chpi[hpi_pick][0], hpi_ons.sum())
    # test that the cHPI signals make some reasonable values
    picks_meg = pick_types(raw.info, meg=True, eeg=False)
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)

    for picks in [picks_meg, picks_eeg]:
        psd_sim, freqs_sim = psd_welch(raw_sim, picks=picks)
        psd_chpi, freqs_chpi = psd_welch(raw_chpi, picks=picks)

        assert_array_equal(freqs_sim, freqs_chpi)
        freq_idx = np.sort([np.argmin(np.abs(freqs_sim - f))
                           for f in hpi_freqs])
        if picks is picks_meg:
            assert_true((psd_chpi[:, freq_idx] >
                         100 * psd_sim[:, freq_idx]).all())
        else:
            assert_allclose(psd_sim, psd_chpi, atol=1e-20)

    # test localization based on cHPI information
    quats_sim = _calculate_chpi_positions(raw_chpi)
    quats = read_head_pos(pos_fname)
    _assert_quats(quats, quats_sim, dist_tol=0.006, angle_tol=4)

def test_calculate_chpi_positions():
    """Test calculation of cHPI positions."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(pos_fname)
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True)
    # This is a little hack (aliasing while decimating) to make it much faster
    # for testing purposes only. We can relax this later if we find it breaks
    # something.
    raw_dec = _decimate_chpi(raw, 15)
    with catch_logging() as log:
        py_quats = _calculate_chpi_positions(raw_dec, verbose='debug')
    assert_true(log.getvalue().startswith('HPIFIT'))
    _assert_quats(py_quats, mf_quats, dist_tol=0.004, angle_tol=2.5)

    # degenerate conditions
    raw_no_chpi = read_raw_fif(test_fif_fname)
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_no_chpi)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = 999
            break
    assert_raises(RuntimeError, _calculate_chpi_positions, raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['r'] = np.ones(3)
    raw_bad.crop(0, 1.)
    with warnings.catch_warnings(record=True):  # bad pos
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, t_step_min=5., verbose=True)
    # ignore HPI info header and [done] footer
    assert_true('0/5 good' in log_file.getvalue().strip().split('\n')[-2])

    # half the rate cuts off cHPI coils
    raw.info['lowpass'] /= 2.
    assert_raises_regex(RuntimeError, 'above the',
                        _calculate_chpi_positions, raw)

def test_movement_compensation():
    """Test movement compensation"""
    temp_dir = _TempDir()
    lims = (0, 4)
    raw = Raw(raw_fname, allow_maxshield='yes', preload=True).crop(*lims)
    head_pos = read_head_pos(pos_fname)

    #
    # Movement compensation, no regularization, no tSSS
    #
    raw_sss = maxwell_filter(raw, head_pos=head_pos, origin=mf_head_origin,
                             regularize=None, bad_condition='ignore')
    assert_meg_snr(raw_sss, Raw(sss_movecomp_fname).crop(*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 = Raw(temp_fname)
    assert_meg_snr(raw_sss, Raw(sss_movecomp_fname).crop(*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, Raw(sss_movecomp_reg_in_fname).crop(*lims),
                   0.5, 1.9, chpi_med_tol=121)

    #
    # Movement compensation,    regularization,    tSSS at the end
    #
    raw_nohpi = filter_chpi(raw.copy())
    with warnings.catch_warnings(record=True) as w:  # untested feature
        raw_sss_mv = maxwell_filter(raw_nohpi, head_pos=head_pos,
                                    st_duration=4., origin=mf_head_origin,
                                    st_fixed=False)
    assert_equal(len(w), 1)
    assert_true('is untested' in str(w[0].message))
    # Neither match is particularly good because our algorithm actually differs
    assert_meg_snr(raw_sss_mv, Raw(sss_movecomp_reg_in_st4s_fname).crop(*lims),
                   0.6, 1.3)
    tSSS_fname = op.join(sss_path, 'test_move_anon_st4s_raw_sss.fif')
    assert_meg_snr(raw_sss_mv, Raw(tSSS_fname).crop(*lims),
                   0.6, 1.0, chpi_med_tol=None)
    assert_meg_snr(Raw(sss_movecomp_reg_in_st4s_fname), Raw(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, Raw(tSSS_fname).crop(*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 = Raw(erm_fname, allow_maxshield='yes')
    assert_raises(ValueError, maxwell_filter, raw_erm, coord_frame='meg',
                  head_pos=head_pos)  # can't do ERM file
    assert_raises(ValueError, maxwell_filter, raw,
                  head_pos=head_pos[:, :9])  # bad shape
    assert_raises(TypeError, maxwell_filter, raw, head_pos='foo')  # bad type
    assert_raises(ValueError, maxwell_filter, raw, head_pos=head_pos[::-1])
    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 0] = raw.first_samp / raw.info['sfreq'] - 1e-2
    assert_raises(ValueError, maxwell_filter, raw, head_pos=head_pos_bad)
    # make sure numerical error doesn't screw it up, though
    head_pos_bad[0, 0] = raw.first_samp / raw.info['sfreq'] - 1e-4
    raw_sss_tweak = maxwell_filter(raw, head_pos=head_pos_bad,
                                   origin=mf_head_origin)
    assert_meg_snr(raw_sss_tweak, raw_sss, 2., 10.)

def test_movement_compensation():
    """Test movement compensation."""
    temp_dir = _TempDir()
    lims = (0, 4)
    raw = read_crop(raw_fname, lims).load_data()
    head_pos = read_head_pos(pos_fname)

    #
    # Movement compensation, no regularization, no tSSS
    #
    raw_sss = maxwell_filter(raw, head_pos=head_pos, origin=mf_head_origin,
                             regularize=None, bad_condition='ignore')
    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())
    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)
    pytest.raises(ValueError, maxwell_filter, raw_erm, coord_frame='meg',
                  head_pos=head_pos)  # can't do ERM file
    pytest.raises(ValueError, maxwell_filter, raw,
                  head_pos=head_pos[:, :9])  # bad shape
    pytest.raises(TypeError, maxwell_filter, raw, head_pos='foo')  # bad type
    pytest.raises(ValueError, maxwell_filter, raw, head_pos=head_pos[::-1])
    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 0] = raw._first_time - 1e-2
    pytest.raises(ValueError, 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)

def test_movement_compensation():
    """Test movement compensation"""
    lims = (0, 8)
    with warnings.catch_warnings(record=True):  # maxshield
        raw = Raw(raw_fname, allow_maxshield=True, preload=True).crop(*lims)
    head_pos = read_head_pos(pos_fname)

    #
    # Movement compensation, no regularization, no tSSS
    #
    raw_sss = maxwell_filter(raw, head_pos=head_pos, origin=mf_head_origin,
                             regularize=None, bad_condition='ignore')
    assert_meg_snr(raw_sss, Raw(sss_movecomp_fname).crop(*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, Raw(sss_movecomp_reg_in_fname).crop(*lims),
                   0.7, 1.9, chpi_med_tol=121)

    #
    # Movement compensation,    regularization,    tSSS at the end
    #
    raw_nohpi = filter_chpi(raw.copy())
    with warnings.catch_warnings(record=True) as w:  # untested feature
        raw_sss_mv = maxwell_filter(raw_nohpi, head_pos=head_pos,
                                    st_duration=4., origin=mf_head_origin,
                                    st_fixed=False)
    assert_equal(len(w), 1)
    assert_true('is untested' in str(w[0].message))
    # Neither match is particularly good because our algorithm actually differs
    assert_meg_snr(raw_sss_mv, Raw(sss_movecomp_reg_in_st4s_fname).crop(*lims),
                   0.6, 1.3)
    tSSS_fname = op.join(sss_path, 'test_move_anon_st4s_raw_sss.fif')
    assert_meg_snr(raw_sss_mv, Raw(tSSS_fname).crop(*lims),
                   0.6, 1.0, chpi_med_tol=None)
    assert_meg_snr(Raw(sss_movecomp_reg_in_st4s_fname), Raw(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, Raw(tSSS_fname).crop(*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
    with warnings.catch_warnings(record=True):  # maxshield
        raw_erm = Raw(erm_fname, allow_maxshield=True)
    assert_raises(ValueError, maxwell_filter, raw_erm, coord_frame='meg',
                  head_pos=head_pos)  # can't do ERM file
    head_pos_bad = head_pos[:, :9]
    assert_raises(ValueError, maxwell_filter, raw,
                  head_pos=head_pos_bad)  # bad shape
    head_pos_bad = 'foo'
    assert_raises(TypeError, maxwell_filter, raw,
                  head_pos=head_pos_bad)  # bad type
    head_pos_bad = head_pos[::-1]
    assert_raises(ValueError, maxwell_filter, raw,
                  head_pos=head_pos_bad)
    head_pos_bad = head_pos.copy()
    head_pos_bad[0, 0] = 1.  # bad time given the first_samp...
    assert_raises(ValueError, maxwell_filter, raw, head_pos=head_pos_bad)