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)
    raw_chpi = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=True)
    # XXX we need to test that the cHPI signals are actually in the correct
    # place, but that should be a subsequent enhancement (not trivial to do so)
    psd_sim, freqs_sim = compute_raw_psd(raw_sim)
    psd_chpi, freqs_chpi = compute_raw_psd(raw_chpi)
    assert_array_equal(freqs_sim, freqs_chpi)
    hpi_freqs = np.array([x['custom_ref'][0]
                          for x in raw.info['hpi_meas'][0]['hpi_coils']])
    freq_idx = np.sort([np.argmin(np.abs(freqs_sim - f)) for f in hpi_freqs])
    picks_meg = pick_types(raw.info, meg=True, eeg=False)
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)
    assert_allclose(psd_sim[picks_eeg], psd_chpi[picks_eeg])
    assert_true((psd_chpi[picks_meg][:, freq_idx] >
                 100 * psd_sim[picks_meg][:, freq_idx]).all())

def test_simulate_raw_bem():
    """Test simulation of raw data with BEM"""
    seed = 42
    raw, src, stc, trans, sphere = _get_data()
    raw_sim_sph = simulate_raw(raw, stc, trans, src, sphere, cov=None, ecg=True, blink=True, random_state=seed)
    raw_sim_bem = simulate_raw(
        raw, stc, trans, src, bem_fname, cov=None, ecg=True, blink=True, random_state=seed, n_jobs=2
    )
    # some components (especially radial) might not match that well,
    # so just make sure that most components have high correlation
    assert_array_equal(raw_sim_sph.ch_names, raw_sim_bem.ch_names)
    picks = pick_types(raw.info, meg=True, eeg=True)
    n_ch = len(picks)
    corr = np.corrcoef(raw_sim_sph[picks][0], raw_sim_bem[picks][0])
    assert_array_equal(corr.shape, (2 * n_ch, 2 * n_ch))
    assert_true(np.median(np.diag(corr[:n_ch, -n_ch:])) > 0.9)

def test_simulate_round_trip():
    """Test simulate_raw round trip calculations."""
    # Check a diagonal round-trip
    raw, src, stc, trans, sphere = _get_data()
    raw.pick_types(meg=True, stim=True)
    bem = read_bem_solution(bem_1_fname)
    old_bem = bem.copy()
    old_src = src.copy()
    old_trans = trans.copy()
    fwd = make_forward_solution(raw.info, trans, src, bem)
    # no omissions
    assert (sum(len(s['vertno']) for s in src) ==
            sum(len(s['vertno']) for s in fwd['src']) ==
            36)
    # make sure things were not modified
    assert (old_bem['surfs'][0]['coord_frame'] ==
            bem['surfs'][0]['coord_frame'])
    assert trans == old_trans
    _compare_source_spaces(src, old_src)
    data = np.eye(fwd['nsource'])
    raw.crop(0, (len(data) - 1) / raw.info['sfreq'])
    stc = SourceEstimate(data, [s['vertno'] for s in fwd['src']],
                         0, 1. / raw.info['sfreq'])
    for use_cps in (False, True):
        this_raw = simulate_raw(raw, stc, trans, src, bem, cov=None,
                                use_cps=use_cps)
        this_raw.pick_types(meg=True, eeg=True)
        assert (old_bem['surfs'][0]['coord_frame'] ==
                bem['surfs'][0]['coord_frame'])
        assert trans == old_trans
        _compare_source_spaces(src, old_src)
        this_fwd = convert_forward_solution(fwd, force_fixed=True,
                                            use_cps=use_cps)
        assert_allclose(this_raw[:][0], this_fwd['sol']['data'],
                        atol=1e-12, rtol=1e-6)

def test_simulate_raw_bem(raw_data):
    """Test simulation of raw data with BEM."""
    raw, src, stc, trans, sphere = raw_data
    src = setup_source_space('sample', 'oct1', subjects_dir=subjects_dir)
    for s in src:
        s['nuse'] = 3
        s['vertno'] = src[1]['vertno'][:3]
        s['inuse'].fill(0)
        s['inuse'][s['vertno']] = 1
    # use different / more complete STC here
    vertices = [s['vertno'] for s in src]
    stc = SourceEstimate(np.eye(sum(len(v) for v in vertices)), vertices,
                         0, 1. / raw.info['sfreq'])
    with pytest.deprecated_call():
        raw_sim_sph = simulate_raw(raw, stc, trans, src, sphere, cov=None,
                                   verbose=True)
    with pytest.deprecated_call():
        raw_sim_bem = simulate_raw(raw, stc, trans, src, bem_fname, cov=None,
                                   n_jobs=2)
    # some components (especially radial) might not match that well,
    # so just make sure that most components have high correlation
    assert_array_equal(raw_sim_sph.ch_names, raw_sim_bem.ch_names)
    picks = pick_types(raw.info, meg=True, eeg=True)
    n_ch = len(picks)
    corr = np.corrcoef(raw_sim_sph[picks][0], raw_sim_bem[picks][0])
    assert_array_equal(corr.shape, (2 * n_ch, 2 * n_ch))
    med_corr = np.median(np.diag(corr[:n_ch, -n_ch:]))
    assert med_corr > 0.65
    # do some round-trip localization
    for s in src:
        transform_surface_to(s, 'head', trans)
    locs = np.concatenate([s['rr'][s['vertno']] for s in src])
    tmax = (len(locs) - 1) / raw.info['sfreq']
    cov = make_ad_hoc_cov(raw.info)
    # The tolerance for the BEM is surprisingly high (28) but I get the same
    # result when using MNE-C and Xfit, even when using a proper 5120 BEM :(
    for use_raw, bem, tol in ((raw_sim_sph, sphere, 2),
                              (raw_sim_bem, bem_fname, 31)):
        events = find_events(use_raw, 'STI 014')
        assert len(locs) == 6
        evoked = Epochs(use_raw, events, 1, 0, tmax, baseline=None).average()
        assert len(evoked.times) == len(locs)
        fits = fit_dipole(evoked, cov, bem, trans, min_dist=1.)[0].pos
        diffs = np.sqrt(np.sum((locs - fits) ** 2, axis=-1)) * 1000
        med_diff = np.median(diffs)
        assert med_diff < tol, '%s: %s' % (bem, med_diff)

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."""
    raw = read_raw_fif(raw_chpi_fname, allow_maxshield='yes',
                       add_eeg_ref=False)
    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 gen_x_raw(n_trials, raw_template, stc_sim, s_dict):
    """Helper to simulate multiple trials of raw data"""

    scaled_cov = deepcopy(s_dict['inv']['noise_cov'])
    scaled_cov['data'] = scaled_cov['data'] * s_dict['noise_scale']

    # XXX: blink rate: 9 to 21.5 blinks/min (higher than found experimentally)
    return simulate_raw(raw_template, stc_sim, s_dict['inv']['mri_head_t'],
                        src=s_dict['inv']['src'], bem=s_dict['bem_fname'],
                        cov=scaled_cov, blink=True, n_jobs=n_jobs,
                        verbose=False)

def test_simulate_raw_sphere():
    """Test simulation of raw data with sphere model."""
    seed = 42
    raw, src, stc, trans, sphere = _get_data()
    assert len(pick_types(raw.info, meg=False, ecg=True)) == 1

    # head pos
    head_pos_sim = dict()
    # these will be at 1., 2., ... sec
    shifts = [[0.001, 0., -0.001], [-0.001, 0.001, 0.]]

    for time_key, shift in enumerate(shifts):
        # Create 4x4 matrix transform and normalize
        temp_trans = deepcopy(raw.info['dev_head_t'])
        temp_trans['trans'][:3, 3] += shift
        head_pos_sim[time_key + 1.] = temp_trans['trans']

    #
    # Test raw simulation with basic parameters
    #
    raw_sim = simulate_raw(raw, stc, trans, src, sphere, read_cov(cov_fname),
                           head_pos=head_pos_sim,
                           blink=True, ecg=True, random_state=seed)
    raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                             cov_fname, head_pos=head_pos_sim,
                             blink=True, ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    std = dict(grad=2e-13, mag=10e-15, eeg=0.1e-6)
    raw_sim = simulate_raw(raw, stc, trans, src, sphere,
                           make_ad_hoc_cov(raw.info, std=std),
                           head_pos=head_pos_sim, blink=True, ecg=True,
                           random_state=seed)
    raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                             cov=std, head_pos=head_pos_sim, blink=True,
                             ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    sphere_norad = make_sphere_model('auto', None, raw.info)
    raw_meg = raw.copy().pick_types()
    raw_sim = simulate_raw(raw_meg, stc, trans, src, sphere_norad,
                           make_ad_hoc_cov(raw.info, std=None),
                           head_pos=head_pos_sim, blink=True, ecg=True,
                           random_state=seed)
    raw_sim_2 = simulate_raw(raw_meg, stc, trans_fname, src_fname,
                             sphere_norad,
                             cov='simple', head_pos=head_pos_sim, blink=True,
                             ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    # Test IO on processed data
    tempdir = _TempDir()
    test_outname = op.join(tempdir, 'sim_test_raw.fif')
    raw_sim.save(test_outname)

    raw_sim_loaded = read_raw_fif(test_outname, preload=True)
    assert_allclose(raw_sim_loaded[:][0], raw_sim[:][0], rtol=1e-6, atol=1e-20)
    del raw_sim, raw_sim_2
    # with no cov (no noise) but with artifacts, most time periods should match
    # but the EOG/ECG channels should not
    for ecg, eog in ((True, False), (False, True), (True, True)):
        raw_sim_3 = simulate_raw(raw, stc, trans, src, sphere,
                                 cov=None, head_pos=head_pos_sim,
                                 blink=eog, ecg=ecg, random_state=seed)
        raw_sim_4 = simulate_raw(raw, stc, trans, src, sphere,
                                 cov=None, head_pos=head_pos_sim,
                                 blink=False, ecg=False, random_state=seed)
        picks = np.arange(len(raw.ch_names))
        diff_picks = pick_types(raw.info, meg=False, ecg=ecg, eog=eog)
        these_picks = np.setdiff1d(picks, diff_picks)
        close = np.isclose(raw_sim_3[these_picks][0],
                           raw_sim_4[these_picks][0], atol=1e-20)
        assert np.mean(close) > 0.7
        far = ~np.isclose(raw_sim_3[diff_picks][0],
                          raw_sim_4[diff_picks][0], atol=1e-20)
        assert np.mean(far) > 0.99
    del raw_sim_3, raw_sim_4

    # make sure it works with EEG-only and MEG-only
    raw_sim_meg = simulate_raw(raw.copy().pick_types(meg=True, eeg=False),
                               stc, trans, src, sphere, cov=None,
                               ecg=True, blink=True, random_state=seed)
    raw_sim_eeg = simulate_raw(raw.copy().pick_types(meg=False, eeg=True),
                               stc, trans, src, sphere, cov=None,
                               ecg=True, blink=True, random_state=seed)
    raw_sim_meeg = simulate_raw(raw.copy().pick_types(meg=True, eeg=True),
                                stc, trans, src, sphere, cov=None,
                                ecg=True, blink=True, random_state=seed)
    assert_allclose(np.concatenate((raw_sim_meg[:][0], raw_sim_eeg[:][0])),
                    raw_sim_meeg[:][0], rtol=1e-7, atol=1e-20)
    del raw_sim_meg, raw_sim_eeg, raw_sim_meeg

    # check that raw-as-info is supported
    raw_sim = simulate_raw(raw, stc, trans, src, sphere, cov=None)
    n_samp = int(round(raw.info['sfreq']))
    for use_raw in (raw, raw.info):
        raw_sim_2 = simulate_raw(use_raw, stc, trans, src, sphere, cov=None,
                                 duration=1.)
        assert len(raw_sim_2.times) == n_samp
        assert_allclose(raw_sim[:, :n_samp][0],
                        raw_sim_2[:, :n_samp][0], rtol=1e-5, atol=1e-30)
    del raw_sim, raw_sim_2

#.........这里部分代码省略.........

#.........这里部分代码省略.........
        print('Omitting brain noise in the simulation')
        cov = None
    else:
        with printer('Loading covariance file for brain noise'):
            cov = mne.read_cov(cov)

    with printer('Setting up spherical model'):
        bem = mne.bem.make_sphere_model('auto', 'auto', raw.info,
                                        verbose=False)
        # check that our sources are reasonable
        rad = bem['layers'][0]['rad']
        r0 = bem['r0']
        outside = np.sqrt(np.sum((rr - r0) ** 2, axis=1)) >= rad
        n_outside = outside.sum()
        if n_outside > 0:
            print('')
            raise ValueError(
                '%s dipole%s outside the spherical model, are your positions '
                'in mm?' % (n_outside, 's were' if n_outside != 1 else ' was'))

    with printer('Constructing source estimate'):
        tmids = t.mean(axis=1)
        t = np.round(t * raw.info['sfreq']).astype(int)
        t[:, 1] += 1  # make it inclusive
        n_samp = int(np.ceil(duration * raw.info['sfreq']))
        data = np.zeros((n_dipoles, n_samp))
        for di, (t_, amp_) in enumerate(zip(t, amp)):
            data[di, t_[0]:t_[1]] = amp_ * np.hanning(t_[1] - t_[0])
        stc = mne.VolSourceEstimate(data, np.arange(n_dipoles),
                                    0, 1. / raw.info['sfreq'])

    # do the simulation
    print('')
    raw_mv = simulate_raw(raw, stc, trans, src, bem, cov=cov, head_pos=pos,
                          chpi=True, n_jobs=n_jobs, verbose=True)
    print('')

    if raw_out:
        with printer('Saving data'):
            raw_mv.save(raw_out, overwrite=overwrite)

    # plot results -- must be *after* save because we low-pass filter
    if plot_dipoles:
        with printer('Plotting dipoles'):
            fig, axs = plt.subplots(1, 3, figsize=(10, 3), facecolor='w')
            fig.canvas.set_window_title('Dipoles')
            meg_info = mne.pick_info(raw.info,
                                     mne.pick_types(raw.info,
                                                    meg=True, eeg=False))
            helmet_rr = [ch['coil_trans'][:3, 3].copy()
                         for ch in meg_info['chs']]
            helmet_nn = np.zeros_like(helmet_rr)
            helmet_nn[:, 2] = 1.
            surf = dict(rr=helmet_rr, nn=helmet_nn,
                        coord_frame=FIFF.FIFFV_COORD_DEVICE)
            helmet_rr = mne.surface.transform_surface_to(
                surf, 'head', meg_info['dev_head_t'])['rr']
            p = np.linspace(0, 2 * np.pi, 40)
            x_sphere, y_sphere = rad * np.sin(p), rad * np.cos(p)
            for ai, ax in enumerate(axs):
                others = np.setdiff1d(np.arange(3), [ai])
                ax.plot(helmet_rr[:, others[0]], helmet_rr[:, others[1]],
                        marker='o', linestyle='none', alpha=0.1,
                        markeredgecolor='none', markerfacecolor='b', zorder=-2)
                ax.plot(x_sphere + r0[others[0]], y_sphere + r0[others[1]],
                        color='y', alpha=0.25, zorder=-1)

def test_simulate_raw_sphere():
    """Test simulation of raw data with sphere model."""
    seed = 42
    raw, src, stc, trans, sphere = _get_data()
    assert_true(len(pick_types(raw.info, meg=False, ecg=True)) == 1)

    # head pos
    head_pos_sim = dict()
    # these will be at 1., 2., ... sec
    shifts = [[0.001, 0., -0.001], [-0.001, 0.001, 0.]]

    for time_key, shift in enumerate(shifts):
        # Create 4x4 matrix transform and normalize
        temp_trans = deepcopy(raw.info['dev_head_t'])
        temp_trans['trans'][:3, 3] += shift
        head_pos_sim[time_key + 1.] = temp_trans['trans']

    #
    # Test raw simulation with basic parameters
    #
    raw_sim = simulate_raw(raw, stc, trans, src, sphere, read_cov(cov_fname),
                           head_pos=head_pos_sim,
                           blink=True, ecg=True, random_state=seed)
    raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                             cov_fname, head_pos=head_pos_sim,
                             blink=True, ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    # Test IO on processed data
    tempdir = _TempDir()
    test_outname = op.join(tempdir, 'sim_test_raw.fif')
    raw_sim.save(test_outname)

    raw_sim_loaded = read_raw_fif(test_outname, preload=True)
    assert_allclose(raw_sim_loaded[:][0], raw_sim[:][0], rtol=1e-6, atol=1e-20)
    del raw_sim, raw_sim_2
    # with no cov (no noise) but with artifacts, most time periods should match
    # but the EOG/ECG channels should not
    for ecg, eog in ((True, False), (False, True), (True, True)):
        raw_sim_3 = simulate_raw(raw, stc, trans, src, sphere,
                                 cov=None, head_pos=head_pos_sim,
                                 blink=eog, ecg=ecg, random_state=seed)
        raw_sim_4 = simulate_raw(raw, stc, trans, src, sphere,
                                 cov=None, head_pos=head_pos_sim,
                                 blink=False, ecg=False, random_state=seed)
        picks = np.arange(len(raw.ch_names))
        diff_picks = pick_types(raw.info, meg=False, ecg=ecg, eog=eog)
        these_picks = np.setdiff1d(picks, diff_picks)
        close = np.isclose(raw_sim_3[these_picks][0],
                           raw_sim_4[these_picks][0], atol=1e-20)
        assert_true(np.mean(close) > 0.7)
        far = ~np.isclose(raw_sim_3[diff_picks][0],
                          raw_sim_4[diff_picks][0], atol=1e-20)
        assert_true(np.mean(far) > 0.99)
    del raw_sim_3, raw_sim_4

    # make sure it works with EEG-only and MEG-only
    raw_sim_meg = simulate_raw(raw.copy().pick_types(meg=True, eeg=False),
                               stc, trans, src, sphere, cov=None,
                               ecg=True, blink=True, random_state=seed)
    raw_sim_eeg = simulate_raw(raw.copy().pick_types(meg=False, eeg=True),
                               stc, trans, src, sphere, cov=None,
                               ecg=True, blink=True, random_state=seed)
    raw_sim_meeg = simulate_raw(raw.copy().pick_types(meg=True, eeg=True),
                                stc, trans, src, sphere, cov=None,
                                ecg=True, blink=True, random_state=seed)
    assert_allclose(np.concatenate((raw_sim_meg[:][0], raw_sim_eeg[:][0])),
                    raw_sim_meeg[:][0], rtol=1e-7, atol=1e-20)
    del raw_sim_meg, raw_sim_eeg, raw_sim_meeg

    # check that different interpolations are similar given small movements
    raw_sim_cos = simulate_raw(raw, stc, trans, src, sphere,
                               head_pos=head_pos_sim,
                               random_state=seed)
    raw_sim_lin = simulate_raw(raw, stc, trans, src, sphere,
                               head_pos=head_pos_sim, interp='linear',
                               random_state=seed)
    assert_allclose(raw_sim_cos[:][0], raw_sim_lin[:][0],
                    rtol=1e-5, atol=1e-20)
    del raw_sim_cos, raw_sim_lin

    # Make impossible transform (translate up into helmet) and ensure failure
    head_pos_sim_err = deepcopy(head_pos_sim)
    head_pos_sim_err[1.][2, 3] -= 0.1  # z trans upward 10cm
    assert_raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
                  ecg=False, blink=False, head_pos=head_pos_sim_err)
    assert_raises(RuntimeError, simulate_raw, raw, stc, trans, src,
                  bem_fname, ecg=False, blink=False,
                  head_pos=head_pos_sim_err)
    # other degenerate conditions
    assert_raises(TypeError, simulate_raw, 'foo', stc, trans, src, sphere)
    assert_raises(TypeError, simulate_raw, raw, 'foo', trans, src, sphere)
    assert_raises(ValueError, simulate_raw, raw, stc.copy().crop(0, 0),
                  trans, src, sphere)
    stc_bad = stc.copy()
    stc_bad.tstep += 0.1
    assert_raises(ValueError, simulate_raw, raw, stc_bad, trans, src, sphere)
    assert_raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
                  chpi=True)  # no cHPI info
    assert_raises(ValueError, simulate_raw, raw, stc, trans, src, sphere,
                  interp='foo')
#.........这里部分代码省略.........

pos = np.zeros((len(rots), 10))
for ii in range(len(pos)):
    pos[ii] = np.concatenate([[ii], rot_to_quat(rots[ii]), center, [0] * 3])
pos[:, 0] += raw.first_samp / raw.info['sfreq']  # initial offset

# Let's activate a vertices bilateral auditory cortices
src = mne.read_source_spaces(src_fname)
labels = mne.read_labels_from_annot('sample', 'aparc.a2009s', 'both',
                                    regexp='G_temp_sup-Plan_tempo',
                                    subjects_dir=subjects_dir)
assert len(labels) == 2  # one left, one right
vertices = [np.intersect1d(l.vertices, s['vertno'])
            for l, s in zip(labels, src)]
data = np.zeros([sum(len(v) for v in vertices), int(raw.info['sfreq'])])
activation = np.hanning(int(raw.info['sfreq'] * 0.2)) * 1e-9  # nAm
t_offset = int(np.ceil(0.2 * raw.info['sfreq']))  # 200 ms in (after baseline)
data[:, t_offset:t_offset + len(activation)] = activation
stc = mne.SourceEstimate(data, vertices, tmin=-0.2,
                         tstep=1. / raw.info['sfreq'])

# Simulate the movement
raw = simulate_raw(raw, stc, trans_fname, src, bem_fname,
                   head_pos=pos, interp='zero', n_jobs=-1)
raw_stat = simulate_raw(raw, stc, trans_fname, src, bem_fname,
                        head_pos=None, n_jobs=-1)
# Save the results
raw.save('simulated_movement_raw.fif', buffer_size_sec=1.)
raw_stat.save('simulated_stationary_raw.fif', buffer_size_sec=1.)
mne.chpi.write_head_quats('simulated_quats.pos', pos)
stc.save('simulated_activation')

def test_iterable():
    """Test iterable support for simulate_raw."""
    raw = read_raw_fif(raw_fname_short).load_data()
    raw.pick_channels(raw.ch_names[:10] + ['STI 014'])
    src = setup_volume_source_space(
        pos=dict(rr=[[-0.05, 0, 0], [0.1, 0, 0]],
                 nn=[[0, 1., 0], [0, 1., 0]]))
    assert src.kind == 'discrete'
    trans = None
    sphere = make_sphere_model(head_radius=None, info=raw.info)
    tstep = 1. / raw.info['sfreq']
    rng = np.random.RandomState(0)
    vertices = np.array([1])
    data = rng.randn(1, 2)
    stc = VolSourceEstimate(data, vertices, 0, tstep)
    assert isinstance(stc.vertices, np.ndarray)
    with pytest.raises(ValueError, match='at least three time points'):
        simulate_raw(raw.info, stc, trans, src, sphere, None)
    data = rng.randn(1, 1000)
    n_events = (len(raw.times) - 1) // 1000 + 1
    stc = VolSourceEstimate(data, vertices, 0, tstep)
    assert isinstance(stc.vertices, np.ndarray)
    with catch_logging() as log:
        with pytest.deprecated_call():
            raw_sim = simulate_raw(raw, stc, trans, src, sphere, None,
                                   verbose=True)
    log = log.getvalue()
    assert 'Making 15 copies of STC' in log
    assert_allclose(raw.times, raw_sim.times)
    events = find_events(raw_sim, initial_event=True)
    assert len(events) == n_events
    assert_array_equal(events[:, 2], 1)

    # Degenerate STCs
    with pytest.raises(RuntimeError,
                       match=r'Iterable did not provide stc\[0\]'):
        simulate_raw(raw.info, [], trans, src, sphere, None)
    with pytest.raises(RuntimeError,
                       match=r'Iterable did not provide stc\[2\].*duration'):
        with pytest.deprecated_call():
            simulate_raw(raw, [stc, stc], trans, src, sphere, None)
    # tuple with ndarray
    event_data = np.zeros(len(stc.times), int)
    event_data[0] = 3
    raw_new = simulate_raw(raw.info, [(stc, event_data)] * 15,
                           trans, src, sphere, None, first_samp=raw.first_samp)
    assert raw_new.n_times == 15000
    raw_new.crop(0, raw_sim.times[-1])
    _assert_iter_sim(raw_sim, raw_new, 3)
    with pytest.raises(ValueError, match='event data had shape .* but need'):
        simulate_raw(raw.info, [(stc, event_data[:-1])], trans, src, sphere,
                     None)
    with pytest.raises(ValueError, match='stim_data in a stc tuple .* int'):
        simulate_raw(raw.info, [(stc, event_data * 1.)], trans, src, sphere,
                     None)

    # iterable
    def stc_iter():
        stim_data = np.zeros(len(stc.times), int)
        stim_data[0] = 4
        ii = 0
        while ii < 100:
            ii += 1
            yield (stc, stim_data)
    with pytest.deprecated_call():
        raw_new = simulate_raw(raw, stc_iter(), trans, src, sphere, None)
    _assert_iter_sim(raw_sim, raw_new, 4)

    def stc_iter_bad():
        ii = 0
        while ii < 100:
            ii += 1
            yield (stc, 4, 3)
    with pytest.raises(ValueError, match='stc, if tuple, must be length'):
        simulate_raw(raw.info, stc_iter_bad(), trans, src, sphere, None)
    _assert_iter_sim(raw_sim, raw_new, 4)

    def stc_iter_bad():
        ii = 0
        while ii < 100:
            ii += 1
            stc_new = stc.copy()
            stc_new.vertices = np.array([ii % 2])
            yield stc_new
    with pytest.raises(RuntimeError, match=r'Vertex mismatch for stc\[1\]'):
        simulate_raw(raw.info, stc_iter_bad(), trans, src, sphere, None)

    # Forward omission
    vertices = np.array([0, 1])
    data = rng.randn(2, 1000)
    stc = VolSourceEstimate(data, vertices, 0, tstep)
    assert isinstance(stc.vertices, np.ndarray)
    # XXX eventually we should support filtering based on sphere radius, too,
    # by refactoring the code in source_space.py that does it!
    surf = _get_ico_surface(3)
    surf['rr'] *= 60  # mm
    model = _surfaces_to_bem([surf], [FIFF.FIFFV_BEM_SURF_ID_BRAIN], [0.3])
    bem = make_bem_solution(model)
    with pytest.warns(RuntimeWarning,
                      match='1 of 2 SourceEstimate vertices'):
#.........这里部分代码省略.........

def test_simulate_round_trip(raw_data):
    """Test simulate_raw round trip calculations."""
    # Check a diagonal round-trip
    raw, src, stc, trans, sphere = raw_data
    raw.pick_types(meg=True, stim=True)
    bem = read_bem_solution(bem_1_fname)
    old_bem = bem.copy()
    old_src = src.copy()
    old_trans = trans.copy()
    fwd = make_forward_solution(raw.info, trans, src, bem)
    # no omissions
    assert (sum(len(s['vertno']) for s in src) ==
            sum(len(s['vertno']) for s in fwd['src']) ==
            36)
    # make sure things were not modified
    assert (old_bem['surfs'][0]['coord_frame'] ==
            bem['surfs'][0]['coord_frame'])
    assert trans == old_trans
    _compare_source_spaces(src, old_src)
    data = np.eye(fwd['nsource'])
    raw.crop(0, (len(data) - 1) / raw.info['sfreq'])
    stc = SourceEstimate(data, [s['vertno'] for s in fwd['src']],
                         0, 1. / raw.info['sfreq'])
    for use_fwd in (None, fwd):
        if use_fwd is None:
            use_trans, use_src, use_bem = trans, src, bem
        else:
            use_trans = use_src = use_bem = None
        with pytest.deprecated_call():
            this_raw = simulate_raw(raw, stc, use_trans, use_src, use_bem,
                                    cov=None, forward=use_fwd)
        this_raw.pick_types(meg=True, eeg=True)
        assert (old_bem['surfs'][0]['coord_frame'] ==
                bem['surfs'][0]['coord_frame'])
        assert trans == old_trans
        _compare_source_spaces(src, old_src)
        this_fwd = convert_forward_solution(fwd, force_fixed=True)
        assert_allclose(this_raw[:][0], this_fwd['sol']['data'],
                        atol=1e-12, rtol=1e-6)
    with pytest.raises(ValueError, match='If forward is not None then'):
        simulate_raw(raw.info, stc, trans, src, bem, forward=fwd)
    # Not iterable
    with pytest.raises(TypeError, match='SourceEstimate, tuple, or iterable'):
        simulate_raw(raw.info, 0., trans, src, bem, None)
    # STC with a source that `src` does not have
    assert 0 not in src[0]['vertno']
    vertices = [[0, fwd['src'][0]['vertno'][0]], []]
    stc_bad = SourceEstimate(data[:2], vertices, 0, 1. / raw.info['sfreq'])
    with pytest.warns(RuntimeWarning,
                      match='1 of 2 SourceEstimate vertices'):
        simulate_raw(raw.info, stc_bad, trans, src, bem)
    assert 0 not in fwd['src'][0]['vertno']
    with pytest.warns(RuntimeWarning,
                      match='1 of 2 SourceEstimate vertices'):
        simulate_raw(raw.info, stc_bad, None, None, None, forward=fwd)
    # dev_head_t mismatch
    fwd['info']['dev_head_t']['trans'][0, 0] = 1.
    with pytest.raises(ValueError, match='dev_head_t.*does not match'):
        simulate_raw(raw.info, stc, None, None, None, forward=fwd)

def test_degenerate(raw_data):
    """Test degenerate conditions."""
    raw, src, stc, trans, sphere = raw_data
    info = raw.info
    # Make impossible transform (translate up into helmet) and ensure failure
    hp_err = _get_head_pos_sim(raw)
    hp_err[1.][2, 3] -= 0.1  # z trans upward 10cm
    with pytest.raises(RuntimeError, match='collided with inner skull'):
        simulate_raw(info, stc, trans, src, sphere, cov=None,
                     head_pos=hp_err)
    # other degenerate conditions
    with pytest.raises(TypeError, match='info must be an instance of'):
        simulate_raw('foo', stc, trans, src, sphere)
    with pytest.raises(TypeError, match='stc must be an instance of'):
        simulate_raw(info, 'foo', trans, src, sphere)
    with pytest.raises(ValueError, match='stc must have at least three time'):
        simulate_raw(info, stc.copy().crop(0, 0), trans, src, sphere)
    with pytest.raises(TypeError, match='must be an instance of Raw or Info'):
        simulate_raw(0, stc, trans, src, sphere)
    stc_bad = stc.copy()
    stc_bad.tstep += 0.1
    with pytest.raises(ValueError, match='same sample rate'):
        simulate_raw(info, stc_bad, trans, src, sphere)
    with pytest.raises(TypeError, match='Covariance matrix type'):
        with pytest.deprecated_call():
            simulate_raw(info, stc, trans, src, sphere, cov=0)
    with pytest.raises(RuntimeError, match='cHPI information not found'):
        with pytest.deprecated_call():
            simulate_raw(info, stc, trans, src, sphere, chpi=True)
    with pytest.raises(ValueError, match='interp must be one of'):
        simulate_raw(info, stc, trans, src, sphere, interp='foo')
    with pytest.raises(TypeError, match='unknown head_pos type'):
        simulate_raw(info, stc, trans, src, sphere, head_pos=1.)
    with pytest.raises(RuntimeError, match='All position times'):
        with pytest.deprecated_call():
            simulate_raw(raw, stc, trans, src, sphere, head_pos=pos_fname)
    head_pos_sim_err = _get_head_pos_sim(raw)
    head_pos_sim_err[-1.] = head_pos_sim_err[1.]  # negative time
    with pytest.raises(RuntimeError, match='All position times'):
        simulate_raw(info, stc, trans, src, sphere,
                     head_pos=head_pos_sim_err)
    raw_bad = raw.copy()
    raw_bad.info['dig'] = None
    with pytest.raises(RuntimeError, match='Cannot fit headshape'):
        with pytest.deprecated_call():
            simulate_raw(raw_bad, stc, trans, src, sphere, blink=True)
    with pytest.raises(RuntimeError, match='Cannot fit headshape'):
        add_eog(raw_bad)

def test_simulate_raw_sphere(raw_data, tmpdir):
    """Test simulation of raw data with sphere model."""
    seed = 42
    raw, src, stc, trans, sphere = raw_data
    assert len(pick_types(raw.info, meg=False, ecg=True)) == 1
    tempdir = str(tmpdir)

    # head pos
    head_pos_sim = _get_head_pos_sim(raw)

    #
    # Test raw simulation with basic parameters
    #
    raw.info.normalize_proj()
    cov = read_cov(cov_fname)
    cov['projs'] = raw.info['projs']
    raw.info['bads'] = raw.ch_names[:1]
    with pytest.deprecated_call(match='cov is deprecated'):
        raw_sim = simulate_raw(raw, stc, trans, src, sphere, cov,
                               head_pos=head_pos_sim,
                               blink=True, ecg=True, random_state=seed,
                               verbose=True)
    with pytest.warns(RuntimeWarning, match='applying projector with'):
        raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                                 cov_fname, head_pos=head_pos_sim,
                                 blink=True, ecg=True, random_state=seed)
    with pytest.raises(RuntimeError, match='Maximum number of STC iterations'):
        simulate_raw(raw.info, [stc] * 5, trans_fname, src_fname, sphere,
                     cov=None, max_iter=1)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    std = dict(grad=2e-13, mag=10e-15, eeg=0.1e-6)
    with pytest.deprecated_call():
        raw_sim = simulate_raw(raw, stc, trans, src, sphere,
                               make_ad_hoc_cov(raw.info, std=std),
                               head_pos=head_pos_sim, blink=True, ecg=True,
                               random_state=seed)
    with pytest.deprecated_call():
        raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                                 cov=std, head_pos=head_pos_sim, blink=True,
                                 ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    sphere_norad = make_sphere_model('auto', None, raw.info)
    raw_meg = raw.copy().pick_types()
    with pytest.deprecated_call():
        raw_sim = simulate_raw(raw_meg, stc, trans, src, sphere_norad,
                               cov=None,
                               head_pos=head_pos_sim, blink=True, ecg=True,
                               random_state=seed)
    with pytest.deprecated_call():
        raw_sim_2 = simulate_raw(raw_meg, stc, trans_fname, src_fname,
                                 sphere_norad, cov=None, head_pos=head_pos_sim,
                                 blink=True, ecg=True, random_state=seed)
    assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
    # Test IO on processed data
    test_outname = op.join(tempdir, 'sim_test_raw.fif')
    raw_sim.save(test_outname)

    raw_sim_loaded = read_raw_fif(test_outname, preload=True)
    assert_allclose(raw_sim_loaded[:][0], raw_sim[:][0], rtol=1e-6, atol=1e-20)
    del raw_sim, raw_sim_2
    # with no cov (no noise) but with artifacts, most time periods should match
    # but the EOG/ECG channels should not
    for ecg, eog in ((True, False), (False, True), (True, True)):
        with pytest.deprecated_call():
            raw_sim_3 = simulate_raw(raw, stc, trans, src, sphere,
                                     cov=None, head_pos=head_pos_sim,
                                     blink=eog, ecg=ecg, random_state=seed)
        with pytest.deprecated_call():
            raw_sim_4 = simulate_raw(raw, stc, trans, src, sphere,
                                     cov=None, head_pos=head_pos_sim,
                                     blink=False, ecg=False, random_state=seed)
        picks = np.arange(len(raw.ch_names))
        diff_picks = pick_types(raw.info, meg=False, ecg=ecg, eog=eog)
        these_picks = np.setdiff1d(picks, diff_picks)
        close = np.isclose(raw_sim_3[these_picks][0],
                           raw_sim_4[these_picks][0], atol=1e-20)
        assert np.mean(close) > 0.7
        far = ~np.isclose(raw_sim_3[diff_picks][0],
                          raw_sim_4[diff_picks][0], atol=1e-20)
        assert np.mean(far) > 0.99
    del raw_sim_3, raw_sim_4

    # make sure it works with EEG-only and MEG-only
    with pytest.deprecated_call():
        raw_sim_meg = simulate_raw(raw.copy().pick_types(meg=True, eeg=False),
                                   stc, trans, src, sphere, cov=None)
        raw_sim_eeg = simulate_raw(raw.copy().pick_types(meg=False, eeg=True),
                                   stc, trans, src, sphere, cov=None)
        raw_sim_meeg = simulate_raw(raw.copy().pick_types(meg=True, eeg=True),
                                    stc, trans, src, sphere, cov=None)
    for this_raw in (raw_sim_meg, raw_sim_eeg, raw_sim_meeg):
        add_eog(this_raw, random_state=seed)
    for this_raw in (raw_sim_meg, raw_sim_meeg):
        add_ecg(this_raw, random_state=seed)
    with pytest.raises(RuntimeError, match='only add ECG artifacts if MEG'):
        add_ecg(raw_sim_eeg)
    assert_allclose(np.concatenate((raw_sim_meg[:][0], raw_sim_eeg[:][0])),
                    raw_sim_meeg[:][0], rtol=1e-7, atol=1e-20)
    del raw_sim_meg, raw_sim_eeg, raw_sim_meeg

#.........这里部分代码省略.........

info.update(sfreq=sfreq, bads=[])

# Only use gradiometers
picks = mne.pick_types(info, meg='grad', stim=True, exclude=())
mne.pick_info(info, picks, copy=False)

# This is the raw object that will be used as a template for the simulation.
raw = mne.io.RawArray(np.zeros((info['nchan'], len(stc.times))), info)

# Define a covariance matrix for the simulated noise. In this tutorial, we use
# a simple diagonal matrix.
cov = mne.cov.make_ad_hoc_cov(info)
cov['data'] *= (20. / snr) ** 2  # Scale the noise to achieve the desired SNR

# Simulate the raw data, with a lowpass filter on the noise
raw = simulate_raw(raw, stc, trans_fname, src_fname, bem_fname, cov=cov,
                   random_state=rand, iir_filter=[4, -4, 0.8])


#####################################