Python common.assert_meg_snr函数代码示例

def test_bads_reconstruction():
    """Test Maxwell filter reconstruction of bad channels."""
    raw = read_crop(raw_fname, (0., 1.))
    raw.info['bads'] = bads
    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.)

def test_chpi_subtraction():
    """Test subtraction of cHPI signals."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True)
    raw.info['bads'] = ['MEG0111']
    raw.del_proj()
    with catch_logging() as log:
        filter_chpi(raw, include_line=False, verbose=True)
    assert 'No average EEG' not in log.getvalue()
    assert '5 cHPI' in log.getvalue()
    # MaxFilter doesn't do quite as well as our algorithm with the last bit
    raw.crop(0, 16)
    # remove cHPI status chans
    raw_c = read_raw_fif(sss_hpisubt_fname).crop(0, 16).load_data()
    raw_c.pick_types(
        meg=True, eeg=True, eog=True, ecg=True, stim=True, misc=True)
    assert_meg_snr(raw, raw_c, 143, 624)

    # Degenerate cases
    raw_nohpi = read_raw_fif(test_fif_fname, preload=True)
    pytest.raises(RuntimeError, filter_chpi, raw_nohpi)

    # When MaxFliter downsamples, like::
    #     $ maxfilter -nosss -ds 2 -f test_move_anon_raw.fif \
    #           -o test_move_anon_ds2_raw.fif
    # it can strip out some values of info, which we emulate here:
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    raw = raw.crop(0, 1).load_data().resample(600., npad='auto')
    raw.info['lowpass'] = 200.
    del raw.info['maxshield']
    del raw.info['hpi_results'][0]['moments']
    del raw.info['hpi_subsystem']['event_channel']
    with catch_logging() as log:
        filter_chpi(raw, verbose=True)
    pytest.raises(ValueError, filter_chpi, raw, t_window=-1)
    assert '2 cHPI' in log.getvalue()

def test_regularization():
    """Test Maxwell filter regularization"""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    min_tols = (100., 2.6, 1.0)
    med_tols = (1000., 21.4, 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):
        with warnings.catch_warnings(record=True):  # maxshield
            raw = Raw(rf, allow_maxshield=True).crop(0., 1., False)
        sss_reg_in = Raw(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
        py_info = raw_sss.info['proc_history'][0]['max_info']['sss_info']
        assert_true(py_info is not None)
        assert_true(len(py_info) > 0)
        mf_info = sss_reg_in.info['proc_history'][0]['max_info']['sss_info']
        n_in = None
        for inf in py_info, mf_info:
            if n_in is None:
                n_in = _get_n_moments(inf['in_order'])
            else:
                assert_equal(n_in, _get_n_moments(inf['in_order']))
            assert_equal(inf['components'][:n_in].sum(), inf['nfree'])
        assert_allclose(py_info['nfree'], mf_info['nfree'],
                        atol=comp_tols[ii], err_msg=rf)

def test_cross_talk():
    """Test Maxwell filter cross-talk cancellation"""
    raw = Raw(raw_fname, allow_maxshield='yes').crop(0., 1., False)
    raw.info['bads'] = bads
    sss_ctc = Raw(sss_ctc_fname)
    raw_sss = maxwell_filter(raw, cross_talk=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_true(len(py_ctc) > 0)
    assert_raises(ValueError, maxwell_filter, raw, cross_talk=raw)
    assert_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_equal(object_diff(py_ctc, mf_ctc), '')
    raw_ctf = Raw(fname_ctf_raw)
    assert_raises(ValueError, maxwell_filter, raw_ctf)  # cannot fit headshape
    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.crop(0, 0.1, copy=True).load_data().pick_channels(
        [raw.ch_names[pi] for pi in pick_types(raw.info, meg=True,
                                               exclude=())[3:]])
    with warnings.catch_warnings(record=True) as w:
        maxwell_filter(raw_missing, cross_talk=ctc_fname)
    assert_equal(len(w), 1)
    assert_true('Not all cross-talk channels in raw' in str(w[0].message))
    # MEG channels not in cross-talk
    assert_raises(RuntimeError, maxwell_filter, raw_ctf, origin=(0., 0., 0.04),
                  cross_talk=ctc_fname)

def test_chpi_subtraction():
    """Test subtraction of cHPI signals"""
    raw = Raw(chpi_fif_fname, allow_maxshield="yes", preload=True)
    with catch_logging() as log:
        filter_chpi(raw, include_line=False, verbose=True)
    assert_true("5 cHPI" in log.getvalue())
    # MaxFilter doesn't do quite as well as our algorithm with the last bit
    raw.crop(0, 16, copy=False)
    # remove cHPI status chans
    raw_c = Raw(sss_hpisubt_fname).crop(0, 16, copy=False).load_data()
    raw_c.pick_types(meg=True, eeg=True, eog=True, ecg=True, stim=True, misc=True)
    assert_meg_snr(raw, raw_c, 143, 624)

    # Degenerate cases
    raw_nohpi = Raw(test_fif_fname, preload=True)
    assert_raises(RuntimeError, filter_chpi, raw_nohpi)

    # When MaxFliter downsamples, like::
    #     $ maxfilter -nosss -ds 2 -f test_move_anon_raw.fif \
    #           -o test_move_anon_ds2_raw.fif
    # it can strip out some values of info, which we emulate here:
    raw = Raw(chpi_fif_fname, allow_maxshield="yes")
    with warnings.catch_warnings(record=True):  # uint cast suggestion
        raw = raw.crop(0, 1).load_data().resample(600.0, npad="auto")
    raw.info["buffer_size_sec"] = np.float64(2.0)
    raw.info["lowpass"] = 200.0
    del raw.info["maxshield"]
    del raw.info["hpi_results"][0]["moments"]
    del raw.info["hpi_subsystem"]["event_channel"]
    with catch_logging() as log:
        filter_chpi(raw, verbose=True)
    assert_true("2 cHPI" in log.getvalue())

def test_bads_reconstruction():
    """Test Maxwell filter reconstruction of bad channels"""
    raw = Raw(raw_fname, allow_maxshield='yes').crop(0., 1., False)
    raw.info['bads'] = bads
    raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, Raw(sss_bad_recon_fname), 300.)

def test_bads_reconstruction():
    """Test Maxwell filter reconstruction of bad channels"""
    with warnings.catch_warnings(record=True):  # maxshield
        raw = Raw(raw_fname, allow_maxshield=True).crop(0.0, 1.0, False)
    raw.info["bads"] = bads
    raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None, bad_condition="ignore")
    assert_meg_snr(raw_sss, Raw(sss_bad_recon_fname), 300.0)

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.9, 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):
        with warnings.catch_warnings(record=True):  # maxshield
            raw = Raw(rf, allow_maxshield=True).crop(0., 1., copy=False)
        with warnings.catch_warnings(record=True):  # head fit off-center
            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 = Raw(sss_fnames[ii])
        assert_meg_snr(sss_py, sss_mf, mins[ii], meds[ii], msg=rf)

def test_fine_calibration():
    """Test Maxwell filter fine calibration"""

    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    with warnings.catch_warnings(record=True):  # maxshield
        raw = Raw(raw_fname, allow_maxshield=True).crop(0., 1., False)
    sss_fine_cal = Raw(sss_fine_cal_fname)

    # Test 1D SSS fine calibration
    raw_sss = maxwell_filter(raw, calibration=fine_cal_fname,
                             origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_fine_cal, 82, 611)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert_true(py_cal is not None)
    assert_true(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)

    # 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.)

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
    raw_sss = maxwell_filter(raw, calibration=fine_cal_fname,
                             origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_fine_cal, 82, 611)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert_true(py_cal is not None)
    assert_true(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)

    # 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)
    assert_raises(RuntimeError, maxwell_filter, raw_ctf, origin=(0., 0., 0.04),
                  calibration=fine_cal_fname)

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 pytest.warns(None):  # 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)

def test_chpi_subtraction():
    """Test subtraction of cHPI signals"""
    raw = Raw(chpi_fif_fname, allow_maxshield='yes', preload=True)
    with catch_logging() as log:
        filter_chpi(raw, include_line=False, verbose=True)
    assert_true('5 cHPI' in log.getvalue())
    # MaxFilter doesn't do quite as well as our algorithm with the last bit
    raw.crop(0, 16, copy=False)
    raw_c = Raw(sss_hpisubt_fname, preload=True).crop(0, 16, copy=False)
    raw_c.pick_types(meg=True, eeg=True, eog=True, ecg=True, stim=True,
                     misc=True, copy=False)  # remove cHPI status chans
    assert_meg_snr(raw, raw_c, 143, 624)

    # Degenerate cases
    raw_nohpi = Raw(test_fif_fname, preload=True)
    assert_raises(RuntimeError, filter_chpi, raw_nohpi)

    # When MaxFliter downsamples, like::
    #     $ maxfilter -nosss -ds 2 -f test_move_anon_raw.fif \
    #           -o test_move_anon_ds2_raw.fif
    # it can strip out some values of info, which we emulate here:
    raw = Raw(chpi_fif_fname, allow_maxshield='yes')
    raw.crop(0, 1, copy=False).load_data()
    raw.resample(600., npad='auto')
    raw.info['buffer_size_sec'] = np.float64(2.)
    raw.info['lowpass'] = 200.
    del raw.info['maxshield']
    del raw.info['hpi_results'][0]['moments']
    del raw.info['hpi_subsystem']['event_channel']
    with catch_logging() as log:
        filter_chpi(raw, verbose=True)
    assert_true('2 cHPI' in log.getvalue())

def test_spatiotemporal_maxwell():
    """Test Maxwell filter (tSSS) spatiotemporal processing"""
    # Load raw testing data
    raw = Raw(raw_fname, allow_maxshield='yes')

    # Test that window is less than length of data
    assert_raises(ValueError, maxwell_filter, raw, st_duration=1000.)

    # Check both 4 and 10 seconds because Elekta handles them differently
    # This is to ensure that std/non-std tSSS windows are correctly handled
    st_durations = [4., 10.]
    tols = [325., 200.]
    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 = Raw(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, copy=False)

        # Test sss computation at the standard head origin. Same cropping issue
        # as mentioned above.
        if st_duration == 10.:
            raw_tsss = maxwell_filter(raw.crop(0, st_duration),
                                      origin=mf_head_origin,
                                      st_duration=st_duration, regularize=None,
                                      bad_condition='ignore')
        else:
            raw_tsss = maxwell_filter(raw, st_duration=st_duration,
                                      origin=mf_head_origin, regularize=None,
                                      bad_condition='ignore', verbose=True)
            raw_tsss_2 = maxwell_filter(raw, st_duration=st_duration,
                                        origin=mf_head_origin, regularize=None,
                                        bad_condition='ignore', st_fixed=False,
                                        verbose=True)
            assert_meg_snr(raw_tsss, raw_tsss_2, 100., 1000.)
            assert_equal(raw_tsss.estimate_rank(), 140)
            assert_equal(raw_tsss_2.estimate_rank(), 140)
        assert_meg_snr(raw_tsss, tsss_bench, tol)
        py_st = raw_tsss.info['proc_history'][0]['max_info']['max_st']
        assert_true(len(py_st) > 0)
        assert_equal(py_st['buflen'], st_duration)
        assert_equal(py_st['subspcorr'], 0.98)

    # Degenerate cases
    assert_raises(ValueError, maxwell_filter, raw, st_duration=10.,
                  st_correlation=0.)

def test_chpi_subtraction():
    """Test subtraction of cHPI signals"""
    raw = Raw(chpi_fif_fname, allow_maxshield='yes', preload=True)
    filter_chpi(raw, include_line=False)
    # MaxFilter doesn't do quite as well as our algorithm with the last bit
    raw.crop(0, 16, copy=False)
    raw_c = Raw(sss_hpisubt_fname, preload=True).crop(0, 16, copy=False)
    raw_c.pick_types(meg=True, eeg=True, eog=True, ecg=True, stim=True,
                     misc=True, copy=False)  # remove cHPI status chans
    assert_meg_snr(raw, raw_c, 143, 624)

    # Degenerate cases
    raw_nohpi = Raw(test_fif_fname, preload=True)
    assert_raises(RuntimeError, filter_chpi, raw_nohpi)

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)