Python Epochs.average方法代码示例

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_detrend():
    """Test detrending of epochs
    """
    # test first-order
    epochs_1 = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                      baseline=None, detrend=1)
    epochs_2 = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                      baseline=None, detrend=None)
    data_picks = pick_types(epochs_1.info, meg=True, eeg=True,
                                 exclude='bads')
    evoked_1 = epochs_1.average()
    evoked_2 = epochs_2.average()
    evoked_2.detrend(1)
    # Due to roundoff these won't be exactly equal, but they should be close
    assert_true(np.allclose(evoked_1.data, evoked_2.data,
                            rtol=1e-8, atol=1e-20))

    # test zeroth-order case
    for preload in [True, False]:
        epochs_1 = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                          baseline=(None, None), preload=preload)
        epochs_2 = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                          baseline=None, preload=preload, detrend=0)
        a = epochs_1.get_data()
        b = epochs_2.get_data()
        # All data channels should be almost equal
        assert_true(np.allclose(a[:, data_picks, :], b[:, data_picks, :],
                                rtol=1e-16, atol=1e-20))
        # There are non-M/EEG channels that should not be equal:
        assert_true(not np.allclose(a, b))

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_epochs_proj():
    """Test handling projection (apply proj in Raw or in Epochs)
    """
    exclude = raw.info["bads"] + ["MEG 2443", "EEG 053"]  # bads + 2 more
    this_picks = pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True, exclude=exclude)
    epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=this_picks, baseline=(None, 0), proj=True)
    assert_true(all(p["active"] is True for p in epochs.info["projs"]))
    evoked = epochs.average()
    assert_true(all(p["active"] is True for p in evoked.info["projs"]))
    data = epochs.get_data()

    raw_proj = io.Raw(raw_fname, proj=True)
    epochs_no_proj = Epochs(
        raw_proj, events[:4], event_id, tmin, tmax, picks=this_picks, baseline=(None, 0), proj=False
    )

    data_no_proj = epochs_no_proj.get_data()
    assert_true(all(p["active"] is True for p in epochs_no_proj.info["projs"]))
    evoked_no_proj = epochs_no_proj.average()
    assert_true(all(p["active"] is True for p in evoked_no_proj.info["projs"]))
    assert_true(epochs_no_proj.proj is True)  # as projs are active from Raw

    assert_array_almost_equal(data, data_no_proj, decimal=8)

    # make sure we can exclude avg ref
    this_picks = pick_types(raw.info, meg=True, eeg=True, stim=True, eog=True, exclude=exclude)
    epochs = Epochs(
        raw, events[:4], event_id, tmin, tmax, picks=this_picks, baseline=(None, 0), proj=True, add_eeg_ref=True
    )
    assert_true(_has_eeg_average_ref_proj(epochs.info["projs"]))
    epochs = Epochs(
        raw, events[:4], event_id, tmin, tmax, picks=this_picks, baseline=(None, 0), proj=True, add_eeg_ref=False
    )
    assert_true(not _has_eeg_average_ref_proj(epochs.info["projs"]))

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_evoked_io_from_epochs():
    """Test IO of evoked data made from epochs
    """
    # offset our tmin so we don't get exactly a zero value when decimating
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        epochs = Epochs(raw, events[:4], event_id, tmin + 0.011, tmax, picks=picks, baseline=(None, 0), decim=5)
    assert_true(len(w) == 1)
    evoked = epochs.average()
    evoked.save(op.join(tempdir, "evoked-ave.fif"))
    evoked2 = read_evokeds(op.join(tempdir, "evoked-ave.fif"))[0]
    assert_allclose(evoked.data, evoked2.data, rtol=1e-4, atol=1e-20)
    assert_allclose(evoked.times, evoked2.times, rtol=1e-4, atol=1 / evoked.info["sfreq"])

    # now let's do one with negative time
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        epochs = Epochs(raw, events[:4], event_id, 0.1, tmax, picks=picks, baseline=(0.1, 0.2), decim=5)
    evoked = epochs.average()
    evoked.save(op.join(tempdir, "evoked-ave.fif"))
    evoked2 = read_evokeds(op.join(tempdir, "evoked-ave.fif"))[0]
    assert_allclose(evoked.data, evoked2.data, rtol=1e-4, atol=1e-20)
    assert_allclose(evoked.times, evoked2.times, rtol=1e-4, atol=1e-20)

    # should be equivalent to a cropped original
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        epochs = Epochs(raw, events[:4], event_id, -0.2, tmax, picks=picks, baseline=(0.1, 0.2), decim=5)
    evoked = epochs.average()
    evoked.crop(0.099, None)
    assert_allclose(evoked.data, evoked2.data, rtol=1e-4, atol=1e-20)
    assert_allclose(evoked.times, evoked2.times, rtol=1e-4, atol=1e-20)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_read_epochs_bad_events():
    """Test epochs when events are at the beginning or the end of the file
    """
    # Event at the beginning
    epochs = Epochs(
        raw, np.array([[raw.first_samp, 0, event_id]]), event_id, tmin, tmax, picks=picks, baseline=(None, 0)
    )
    with warnings.catch_warnings(record=True):
        evoked = epochs.average()

    epochs = Epochs(
        raw, np.array([[raw.first_samp, 0, event_id]]), event_id, tmin, tmax, picks=picks, baseline=(None, 0)
    )
    epochs.drop_bad_epochs()
    with warnings.catch_warnings(record=True):
        evoked = epochs.average()

    # Event at the end
    epochs = Epochs(
        raw, np.array([[raw.last_samp, 0, event_id]]), event_id, tmin, tmax, picks=picks, baseline=(None, 0)
    )

    with warnings.catch_warnings(record=True):
        evoked = epochs.average()
        assert evoked
    warnings.resetwarnings()

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_subtract_evoked():
    """Test subtraction of Evoked from Epochs
    """
    epochs = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks, baseline=(None, 0))

    # make sure subraction fails if data channels are missing
    assert_raises(ValueError, epochs.subtract_evoked, epochs.average(picks[:5]))

    # do the subraction using the default argument
    epochs.subtract_evoked()

    # apply SSP now
    epochs.apply_proj()

    # use preloading and SSP from the start
    epochs2 = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True, proj=True)

    evoked = epochs2.average()
    epochs2.subtract_evoked(evoked)

    # this gives the same result
    assert_allclose(epochs.get_data(), epochs2.get_data())

    # if we compute the evoked response after subtracting it we get zero
    zero_evoked = epochs.average()
    data = zero_evoked.data
    assert_array_almost_equal(data, np.zeros_like(data), decimal=20)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_compute_proj_epochs():
    """Test SSP computation on epochs"""
    event_id, tmin, tmax = 1, -0.2, 0.3

    raw = Raw(raw_fname, preload=True)
    events = read_events(event_fname)
    bad_ch = 'MEG 2443'
    picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=False,
                       exclude=[])
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=None, proj=False)

    evoked = epochs.average()
    projs = compute_proj_epochs(epochs, n_grad=1, n_mag=1, n_eeg=0, n_jobs=1)
    write_proj(op.join(tempdir, 'proj.fif.gz'), projs)
    for p_fname in [proj_fname, proj_gz_fname,
                    op.join(tempdir, 'proj.fif.gz')]:
        projs2 = read_proj(p_fname)

        assert_true(len(projs) == len(projs2))

        for p1, p2 in zip(projs, projs2):
            assert_true(p1['desc'] == p2['desc'])
            assert_true(p1['data']['col_names'] == p2['data']['col_names'])
            assert_true(p1['active'] == p2['active'])
            # compare with sign invariance
            p1_data = p1['data']['data'] * np.sign(p1['data']['data'][0, 0])
            p2_data = p2['data']['data'] * np.sign(p2['data']['data'][0, 0])
            if bad_ch in p1['data']['col_names']:
                bad = p1['data']['col_names'].index('MEG 2443')
                mask = np.ones(p1_data.size, dtype=np.bool)
                mask[bad] = False
                p1_data = p1_data[:, mask]
                p2_data = p2_data[:, mask]
            corr = np.corrcoef(p1_data, p2_data)[0, 1]
            assert_array_almost_equal(corr, 1.0, 5)

    # test that you can compute the projection matrix
    projs = activate_proj(projs)
    proj, nproj, U = make_projector(projs, epochs.ch_names, bads=[])

    assert_true(nproj == 2)
    assert_true(U.shape[1] == 2)

    # test that you can save them
    epochs.info['projs'] += projs
    evoked = epochs.average()
    evoked.save(op.join(tempdir, 'foo.fif'))

    projs = read_proj(proj_fname)

    projs_evoked = compute_proj_evoked(evoked, n_grad=1, n_mag=1, n_eeg=0)
    assert_true(len(projs_evoked) == 2)
    # XXX : test something

    # test parallelization
    projs = compute_proj_epochs(epochs, n_grad=1, n_mag=1, n_eeg=0, n_jobs=2)
    projs = activate_proj(projs)
    proj_par, _, _ = make_projector(projs, epochs.ch_names, bads=[])
    assert_allclose(proj, proj_par, rtol=1e-8, atol=1e-16)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_ica_ctf():
    """Test run ICA computation on ctf data with/without compensation."""
    method = 'fastica'
    raw = read_raw_ctf(ctf_fname, preload=True)
    events = make_fixed_length_events(raw, 99999)
    for comp in [0, 1]:
        raw.apply_gradient_compensation(comp)
        epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
        evoked = epochs.average()

        # test fit
        for inst in [raw, epochs]:
            ica = ICA(n_components=2, random_state=0, max_iter=2,
                      method=method)
            with pytest.warns(UserWarning, match='did not converge'):
                ica.fit(inst)

        # test apply and get_sources
        for inst in [raw, epochs, evoked]:
            ica.apply(inst)
            ica.get_sources(inst)

    # test mixed compensation case
    raw.apply_gradient_compensation(0)
    ica = ICA(n_components=2, random_state=0, max_iter=2, method=method)
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw)
    raw.apply_gradient_compensation(1)
    epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
    evoked = epochs.average()
    for inst in [raw, epochs, evoked]:
        with pytest.raises(RuntimeError, match='Compensation grade of ICA'):
            ica.apply(inst)
        with pytest.raises(RuntimeError, match='Compensation grade of ICA'):
            ica.get_sources(inst)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_ica_eeg():
    """Test ICA on EEG."""
    method = 'fastica'
    raw_fif = read_raw_fif(fif_fname, preload=True)
    with pytest.warns(RuntimeWarning, match='events'):
        raw_eeglab = read_raw_eeglab(input_fname=eeglab_fname,
                                     montage=eeglab_montage, preload=True)
    for raw in [raw_fif, raw_eeglab]:
        events = make_fixed_length_events(raw, 99999, start=0, stop=0.3,
                                          duration=0.1)
        picks_meg = pick_types(raw.info, meg=True, eeg=False)[:2]
        picks_eeg = pick_types(raw.info, meg=False, eeg=True)[:2]
        picks_all = []
        picks_all.extend(picks_meg)
        picks_all.extend(picks_eeg)
        epochs = Epochs(raw, events, None, -0.1, 0.1, preload=True)
        evoked = epochs.average()

        for picks in [picks_meg, picks_eeg, picks_all]:
            if len(picks) == 0:
                continue
            # test fit
            for inst in [raw, epochs]:
                ica = ICA(n_components=2, random_state=0, max_iter=2,
                          method=method)
                with pytest.warns(None):
                    ica.fit(inst, picks=picks)

            # test apply and get_sources
            for inst in [raw, epochs, evoked]:
                ica.apply(inst)
                ica.get_sources(inst)

    with pytest.warns(RuntimeWarning, match='MISC channel'):
        raw = read_raw_ctf(ctf_fname2,  preload=True)
    events = make_fixed_length_events(raw, 99999, start=0, stop=0.2,
                                      duration=0.1)
    picks_meg = pick_types(raw.info, meg=True, eeg=False)[:2]
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)[:2]
    picks_all = picks_meg + picks_eeg
    for comp in [0, 1]:
        raw.apply_gradient_compensation(comp)
        epochs = Epochs(raw, events, None, -0.1, 0.1, preload=True)
        evoked = epochs.average()

        for picks in [picks_meg, picks_eeg, picks_all]:
            if len(picks) == 0:
                continue
            # test fit
            for inst in [raw, epochs]:
                ica = ICA(n_components=2, random_state=0, max_iter=2,
                          method=method)
                with pytest.warns(None):
                    ica.fit(inst)

            # test apply and get_sources
            for inst in [raw, epochs, evoked]:
                ica.apply(inst)
                ica.get_sources(inst)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_render_report():
    """Test rendering -*.fif files for mne report.
    """
    tempdir = _TempDir()
    raw_fname_new = op.join(tempdir, 'temp_raw.fif')
    event_fname_new = op.join(tempdir, 'temp_raw-eve.fif')
    cov_fname_new = op.join(tempdir, 'temp_raw-cov.fif')
    fwd_fname_new = op.join(tempdir, 'temp_raw-fwd.fif')
    inv_fname_new = op.join(tempdir, 'temp_raw-inv.fif')
    for a, b in [[raw_fname, raw_fname_new],
                 [event_fname, event_fname_new],
                 [cov_fname, cov_fname_new],
                 [fwd_fname, fwd_fname_new],
                 [inv_fname, inv_fname_new]]:
        shutil.copyfile(a, b)

    # create and add -epo.fif and -ave.fif files
    epochs_fname = op.join(tempdir, 'temp-epo.fif')
    evoked_fname = op.join(tempdir, 'temp-ave.fif')
    raw = Raw(raw_fname_new)
    picks = pick_types(raw.info, meg='mag', eeg=False)  # faster with one type
    epochs = Epochs(raw, read_events(event_fname), 1, -0.2, 0.2, picks=picks)
    epochs.save(epochs_fname)
    epochs.average().save(evoked_fname)

    report = Report(info_fname=raw_fname_new, subjects_dir=subjects_dir)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        report.parse_folder(data_path=tempdir)
    assert_true(len(w) >= 1)

    # Check correct paths and filenames
    fnames = glob.glob(op.join(tempdir, '*.fif'))
    for fname in fnames:
        assert_true(op.basename(fname) in
                    [op.basename(x) for x in report.fnames])
        assert_true(''.join(report.html).find(op.basename(fname)) != -1)

    assert_equal(len(report.fnames), len(fnames))
    assert_equal(len(report.html), len(report.fnames))

    # Check saving functionality
    report.data_path = tempdir
    report.save(fname=op.join(tempdir, 'report.html'), open_browser=False)
    assert_true(op.isfile(op.join(tempdir, 'report.html')))

    assert_equal(len(report.html), len(fnames))
    assert_equal(len(report.html), len(report.fnames))

    # Check saving same report to new filename
    report.save(fname=op.join(tempdir, 'report2.html'), open_browser=False)
    assert_true(op.isfile(op.join(tempdir, 'report2.html')))

    # Check overwriting file
    report.save(fname=op.join(tempdir, 'report.html'), open_browser=False,
                overwrite=True)
    assert_true(op.isfile(op.join(tempdir, 'report.html')))

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_compute_proj():
    """Test SSP computation"""
    event_id, tmin, tmax = 1, -0.2, 0.3

    raw = Raw(raw_fname)
    events = read_events(event_fname)
    exclude = []
    bad_ch = 'MEG 2443'
    picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=False,
                            exclude=exclude)
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                        baseline=None, proj=False)

    evoked = epochs.average()
    projs = compute_proj_epochs(epochs, n_grad=1, n_mag=1, n_eeg=0)

    projs2 = read_proj(proj_fname)

    assert_true(len(projs) == len(projs2))

    for p1, p2 in zip(projs, projs2):
        assert_true(p1['desc'] == p2['desc'])
        assert_true(p1['data']['col_names'] == p2['data']['col_names'])
        assert_true(p1['active'] == p2['active'])
        # compare with sign invariance
        p1_data = p1['data']['data'] * np.sign(p1['data']['data'][0, 0])
        p2_data = p2['data']['data'] * np.sign(p2['data']['data'][0, 0])
        if bad_ch in p1['data']['col_names']:
            bad = p1['data']['col_names'].index('MEG 2443')
            mask = np.ones(p1_data.size, dtype=np.bool)
            mask[bad] = False
            p1_data = p1_data[:, mask]
            p2_data = p2_data[:, mask]
        corr = np.corrcoef(p1_data, p2_data)[0, 1]
        assert_array_almost_equal(corr, 1.0, 7)

    # test that you can compute the projection matrix
    projs = activate_proj(projs)
    proj, nproj, U = make_projector(projs, epochs.ch_names, bads=[])

    assert_true(nproj == 2)
    assert_true(U.shape[1] == 2)

    # test that you can save them
    epochs.info['projs'] += projs
    evoked = epochs.average()
    evoked.save('foo.fif')

    projs = read_proj(proj_fname)

    projs_evoked = compute_proj_evoked(evoked, n_grad=1, n_mag=1, n_eeg=0)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_preload_epochs():
    """Test preload of epochs
    """
    epochs_preload = Epochs(
        raw, events[:16], event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True, reject=reject, flat=flat
    )
    data_preload = epochs_preload.get_data()

    epochs = Epochs(
        raw, events[:16], event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=False, reject=reject, flat=flat
    )
    data = epochs.get_data()
    assert_array_equal(data_preload, data)
    assert_array_almost_equal(epochs_preload.average().data, epochs.average().data, 18)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_evoked_arithmetic():
    """Test arithmetic of evoked data
    """
    epochs1 = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks, baseline=(None, 0))
    evoked1 = epochs1.average()
    epochs2 = Epochs(raw, events[4:8], event_id, tmin, tmax, picks=picks, baseline=(None, 0))
    evoked2 = epochs2.average()
    epochs = Epochs(raw, events[:8], event_id, tmin, tmax, picks=picks, baseline=(None, 0))
    evoked = epochs.average()
    evoked_sum = evoked1 + evoked2
    assert_array_equal(evoked.data, evoked_sum.data)
    assert_array_equal(evoked.times, evoked_sum.times)
    assert_true(evoked_sum.nave == (evoked1.nave + evoked2.nave))
    evoked_diff = evoked1 - evoked1
    assert_array_equal(np.zeros_like(evoked.data), evoked_diff.data)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_xdawn_apply_transform():
    """Test Xdawn apply and transform."""
    # Get data
    raw, events, picks = _get_data()
    raw.pick_types(eeg=True, meg=False)
    epochs = Epochs(raw, events, event_id, tmin, tmax, proj=False,
                    preload=True, baseline=None,
                    verbose=False)
    n_components = 2
    # Fit Xdawn
    xd = Xdawn(n_components=n_components, correct_overlap=False)
    xd.fit(epochs)

    # Apply on different types of instances
    for inst in [raw, epochs.average(), epochs]:
        denoise = xd.apply(inst)
    # Apply on other thing should raise an error
    assert_raises(ValueError, xd.apply, 42)

    # Transform on epochs
    xd.transform(epochs)
    # Transform on ndarray
    xd.transform(epochs._data)
    # Transform on someting else
    assert_raises(ValueError, xd.transform, 42)

    # Check numerical results with shuffled epochs
    np.random.seed(0)  # random makes unstable linalg
    idx = np.arange(len(epochs))
    np.random.shuffle(idx)
    xd.fit(epochs[idx])
    denoise_shfl = xd.apply(epochs)
    assert_array_almost_equal(denoise['cond2']._data,
                              denoise_shfl['cond2']._data)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_evoked_standard_error():
    """Test calculation and read/write of standard error
    """
    epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0))
    evoked = [epochs.average(), epochs.standard_error()]
    io.write_evokeds(op.join(tempdir, 'evoked.fif'), evoked)
    evoked2 = read_evokeds(op.join(tempdir, 'evoked.fif'), [0, 1])
    evoked3 = [read_evokeds(op.join(tempdir, 'evoked.fif'), 'Unknown'),
               read_evokeds(op.join(tempdir, 'evoked.fif'), 'Unknown',
                            kind='standard_error')]
    for evoked_new in [evoked2, evoked3]:
        assert_true(evoked_new[0]._aspect_kind ==
                    FIFF.FIFFV_ASPECT_AVERAGE)
        assert_true(evoked_new[0].kind == 'average')
        assert_true(evoked_new[1]._aspect_kind ==
                    FIFF.FIFFV_ASPECT_STD_ERR)
        assert_true(evoked_new[1].kind == 'standard_error')
        for ave, ave2 in zip(evoked, evoked_new):
            assert_array_almost_equal(ave.data, ave2.data)
            assert_array_almost_equal(ave.times, ave2.times)
            assert_equal(ave.nave, ave2.nave)
            assert_equal(ave._aspect_kind, ave2._aspect_kind)
            assert_equal(ave.kind, ave2.kind)
            assert_equal(ave.last, ave2.last)
            assert_equal(ave.first, ave2.first)

# 需要导入模块: from mne import Epochs [as 别名]
# 或者: from mne.Epochs import average [as 别名]
def test_acqparser_averaging():
    """Test averaging with AcqParserFIF vs. Elekta software."""
    raw = read_raw_fif(fname_raw_elekta, preload=True)
    acqp = AcqParserFIF(raw.info)
    for cat in acqp.categories:
        # XXX datasets match only when baseline is applied to both,
        # not sure where relative dc shift comes from
        cond = acqp.get_condition(raw, cat)
        eps = Epochs(raw, baseline=(-.05, 0), **cond)
        ev = eps.average()
        ev_ref = read_evokeds(fname_ave_elekta, cat['comment'],
                              baseline=(-.05, 0), proj=False)
        ev_mag = ev.copy()
        ev_mag.pick_channels(['MEG0111'])
        ev_grad = ev.copy()
        ev_grad.pick_channels(['MEG2643', 'MEG1622'])
        ev_ref_mag = ev_ref.copy()
        ev_ref_mag.pick_channels(['MEG0111'])
        ev_ref_grad = ev_ref.copy()
        ev_ref_grad.pick_channels(['MEG2643', 'MEG1622'])
        assert_allclose(ev_mag.data, ev_ref_mag.data,
                        rtol=0, atol=1e-15)  # tol = 1 fT
        # Elekta put these in a different order
        assert ev_grad.ch_names[::-1] == ev_ref_grad.ch_names
        assert_allclose(ev_grad.data[::-1], ev_ref_grad.data,
                        rtol=0, atol=1e-13)  # tol = 1 fT/cm