def test_dipole_fitting_fixed():
    """Test dipole fitting with a fixed position."""
    import matplotlib.pyplot as plt
    tpeak = 0.073
    sphere = make_sphere_model(head_radius=0.1)
    evoked = read_evokeds(fname_evo, baseline=(None, 0))[0]
    evoked.pick_types(meg=True)
    t_idx = np.argmin(np.abs(tpeak - evoked.times))
    evoked_crop = evoked.copy().crop(tpeak, tpeak)
    assert_equal(len(evoked_crop.times), 1)
    cov = read_cov(fname_cov)
    dip_seq, resid = fit_dipole(evoked_crop, cov, sphere)
    assert isinstance(dip_seq, Dipole)
    assert isinstance(resid, Evoked)
    assert_equal(len(dip_seq.times), 1)
    pos, ori, gof = dip_seq.pos[0], dip_seq.ori[0], dip_seq.gof[0]
    amp = dip_seq.amplitude[0]
    # Fix position, allow orientation to change
    dip_free, resid_free = fit_dipole(evoked, cov, sphere, pos=pos)
    assert isinstance(dip_free, Dipole)
    assert isinstance(resid_free, Evoked)
    assert_allclose(dip_free.times, evoked.times)
    assert_allclose(np.tile(pos[np.newaxis], (len(evoked.times), 1)),
                    dip_free.pos)
    assert_allclose(ori, dip_free.ori[t_idx])  # should find same ori
    assert (np.dot(dip_free.ori, ori).mean() < 0.9)  # but few the same
    assert_allclose(gof, dip_free.gof[t_idx])  # ... same gof
    assert_allclose(amp, dip_free.amplitude[t_idx])  # and same amp
    assert_allclose(resid.data, resid_free.data[:, [t_idx]])
    # Fix position and orientation
    dip_fixed, resid_fixed = fit_dipole(evoked, cov, sphere, pos=pos, ori=ori)
    assert (isinstance(dip_fixed, DipoleFixed))
    assert_allclose(dip_fixed.times, evoked.times)
    assert_allclose(dip_fixed.info['chs'][0]['loc'][:3], pos)
    assert_allclose(dip_fixed.info['chs'][0]['loc'][3:6], ori)
    assert_allclose(dip_fixed.data[1, t_idx], gof)
    assert_allclose(resid.data, resid_fixed.data[:, [t_idx]])
    _check_roundtrip_fixed(dip_fixed)
    # bad resetting
    evoked.info['bads'] = [evoked.ch_names[3]]
    dip_fixed, resid_fixed = fit_dipole(evoked, cov, sphere, pos=pos, ori=ori)
    # Degenerate conditions
    evoked_nan = evoked.copy().crop(0, 0)
    evoked_nan.data[0, 0] = None
    pytest.raises(ValueError, fit_dipole, evoked_nan, cov, sphere)
    pytest.raises(ValueError, fit_dipole, evoked, cov, sphere, ori=[1, 0, 0])
    pytest.raises(ValueError, fit_dipole, evoked, cov, sphere, pos=[0, 0, 0],
                  ori=[2, 0, 0])
    pytest.raises(ValueError, fit_dipole, evoked, cov, sphere, pos=[0.1, 0, 0])
    # copying
    dip_fixed_2 = dip_fixed.copy()
    dip_fixed_2.data[:] = 0.
    assert not np.isclose(dip_fixed.data, 0., atol=1e-20).any()
    # plotting
    plt.close('all')
    dip_fixed.plot()
    plt.close('all')

def test_dipole_fitting_ctf():
    """Test dipole fitting with CTF data."""
    raw_ctf = read_raw_ctf(fname_ctf).set_eeg_reference()
    events = make_fixed_length_events(raw_ctf, 1)
    evoked = Epochs(raw_ctf, events, 1, 0, 0, baseline=None).average()
    cov = make_ad_hoc_cov(evoked.info)
    sphere = make_sphere_model((0., 0., 0.))
    # XXX Eventually we should do some better checks about accuracy, but
    # for now our CTF phantom fitting tutorials will have to do
    # (otherwise we need to add that to the testing dataset, which is
    # a bit too big)
    fit_dipole(evoked, cov, sphere)

def test_min_distance_fit_dipole():
    """Test dipole min_dist to inner_skull."""
    subject = 'sample'
    raw = read_raw_fif(fname_raw, preload=True)

    # select eeg data
    picks = pick_types(raw.info, meg=False, eeg=True, exclude='bads')
    info = pick_info(raw.info, picks)

    # Let's use cov = Identity
    cov = read_cov(fname_cov)
    cov['data'] = np.eye(cov['data'].shape[0])

    # Simulated scal map
    simulated_scalp_map = np.zeros(picks.shape[0])
    simulated_scalp_map[27:34] = 1

    simulated_scalp_map = simulated_scalp_map[:, None]

    evoked = EvokedArray(simulated_scalp_map, info, tmin=0)

    min_dist = 5.  # distance in mm

    bem = read_bem_solution(fname_bem)
    dip, residual = fit_dipole(evoked, cov, bem, fname_trans,
                               min_dist=min_dist)

    dist = _compute_depth(dip, fname_bem, fname_trans, subject, subjects_dir)

    # Constraints are not exact, so bump the minimum slightly
    assert_true(min_dist - 0.1 < (dist[0] * 1000.) < (min_dist + 1.))

    assert_raises(ValueError, fit_dipole, evoked, cov, fname_bem, fname_trans,
                  -1.)

def test_dipole_fitting_fixed():
    """Test dipole fitting with a fixed position"""
    tpeak = 0.073
    sphere = make_sphere_model(head_radius=0.1)
    evoked = read_evokeds(fname_evo, baseline=(None, 0))[0]
    evoked.pick_types(meg=True)
    t_idx = np.argmin(np.abs(tpeak - evoked.times))
    evoked_crop = evoked.copy().crop(tpeak, tpeak)
    assert_equal(len(evoked_crop.times), 1)
    cov = read_cov(fname_cov)
    dip_seq, resid = fit_dipole(evoked_crop, cov, sphere)
    assert_true(isinstance(dip_seq, Dipole))
    assert_equal(len(dip_seq.times), 1)
    pos, ori, gof = dip_seq.pos[0], dip_seq.ori[0], dip_seq.gof[0]
    amp = dip_seq.amplitude[0]
    # Fix position, allow orientation to change
    dip_free, resid_free = fit_dipole(evoked, cov, sphere, pos=pos)
    assert_true(isinstance(dip_free, Dipole))
    assert_allclose(dip_free.times, evoked.times)
    assert_allclose(np.tile(pos[np.newaxis], (len(evoked.times), 1)),
                    dip_free.pos)
    assert_allclose(ori, dip_free.ori[t_idx])  # should find same ori
    assert_true(np.dot(dip_free.ori, ori).mean() < 0.9)  # but few the same
    assert_allclose(gof, dip_free.gof[t_idx])  # ... same gof
    assert_allclose(amp, dip_free.amplitude[t_idx])  # and same amp
    assert_allclose(resid, resid_free[:, [t_idx]])
    # Fix position and orientation
    dip_fixed, resid_fixed = fit_dipole(evoked, cov, sphere, pos=pos, ori=ori)
    assert_true(isinstance(dip_fixed, DipoleFixed))
    assert_allclose(dip_fixed.times, evoked.times)
    assert_allclose(dip_fixed.info['chs'][0]['loc'][:3], pos)
    assert_allclose(dip_fixed.info['chs'][0]['loc'][3:6], ori)
    assert_allclose(dip_fixed.data[1, t_idx], gof)
    assert_allclose(resid, resid_fixed[:, [t_idx]])
    _check_roundtrip_fixed(dip_fixed)
    # Degenerate conditions
    assert_raises(ValueError, fit_dipole, evoked, cov, sphere, pos=[0])
    assert_raises(ValueError, fit_dipole, evoked, cov, sphere, ori=[1, 0, 0])
    assert_raises(ValueError, fit_dipole, evoked, cov, sphere, pos=[0, 0, 0],
                  ori=[2, 0, 0])
    assert_raises(ValueError, fit_dipole, evoked, cov, sphere, pos=[0.1, 0, 0])

def test_gamma_map_vol_sphere():
    """Gamma MAP with a sphere forward and volumic source space."""
    evoked = read_evokeds(fname_evoked, condition=0, baseline=(None, 0),
                          proj=False)
    evoked.resample(50, npad=100)
    evoked.crop(tmin=0.1, tmax=0.16)  # crop to window around peak

    cov = read_cov(fname_cov)
    cov = regularize(cov, evoked.info, rank=None)

    info = evoked.info
    sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.080)
    src = mne.setup_volume_source_space(subject=None, pos=30., mri=None,
                                        sphere=(0.0, 0.0, 0.0, 80.0),
                                        bem=None, mindist=5.0,
                                        exclude=2.0)
    fwd = mne.make_forward_solution(info, trans=None, src=src, bem=sphere,
                                    eeg=False, meg=True)

    alpha = 0.5
    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
                  loose=0, return_residual=False)

    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
                  loose=0.2, return_residual=False)

    stc = gamma_map(evoked, fwd, cov, alpha, tol=1e-4,
                    xyz_same_gamma=False, update_mode=2,
                    return_residual=False)

    assert_array_almost_equal(stc.times, evoked.times, 5)

    # Compare orientation obtained using fit_dipole and gamma_map
    # for a simulated evoked containing a single dipole
    stc = mne.VolSourceEstimate(50e-9 * np.random.RandomState(42).randn(1, 4),
                                vertices=stc.vertices[:1],
                                tmin=stc.tmin,
                                tstep=stc.tstep)
    evoked_dip = mne.simulation.simulate_evoked(fwd, stc, info, cov, nave=1e9,
                                                use_cps=True)

    dip_gmap = gamma_map(evoked_dip, fwd, cov, 0.1, return_as_dipoles=True)

    amp_max = [np.max(d.amplitude) for d in dip_gmap]
    dip_gmap = dip_gmap[np.argmax(amp_max)]
    assert (dip_gmap[0].pos[0] in src[0]['rr'][stc.vertices])

    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
    assert (np.abs(np.dot(dip_fit.ori[0], dip_gmap.ori[0])) > 0.99)

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 compute_bias(raw):
    events = find_events(raw, 'STI201', verbose=False)
    events = events[1:]  # first one has an artifact
    tmin, tmax = -0.2, 0.1
    epochs = mne.Epochs(raw, events, dipole_number, tmin, tmax,
                        baseline=(None, -0.01), preload=True, verbose=False)
    sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=None,
                                   verbose=False)
    cov = mne.compute_covariance(epochs, tmax=0, method='oas',
                                 rank=None, verbose=False)
    idx = epochs.time_as_index(0.036)[0]
    data = epochs.get_data()[:, :, idx].T
    evoked = mne.EvokedArray(data, epochs.info, tmin=0.)
    dip = fit_dipole(evoked, cov, sphere, n_jobs=1, verbose=False)[0]
    actual_pos = mne.dipole.get_phantom_dipoles()[0][dipole_number - 1]
    misses = 1000 * np.linalg.norm(dip.pos - actual_pos, axis=-1)
    return misses

def test_min_distance_fit_dipole():
    """Test dipole min_dist to inner_skull"""
    data_path = testing.data_path()
    raw_fname = data_path + '/MEG/sample/sample_audvis_trunc_raw.fif'

    subjects_dir = op.join(data_path, 'subjects')
    fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis-cov.fif')
    fname_trans = op.join(data_path, 'MEG', 'sample',
                          'sample_audvis_trunc-trans.fif')
    fname_bem = op.join(subjects_dir, 'sample', 'bem',
                        'sample-1280-1280-1280-bem-sol.fif')

    subject = 'sample'

    raw = Raw(raw_fname, preload=True)

    # select eeg data
    picks = pick_types(raw.info, meg=False, eeg=True, exclude='bads')
    info = pick_info(raw.info, picks)

    # Let's use cov = Identity
    cov = read_cov(fname_cov)
    cov['data'] = np.eye(cov['data'].shape[0])

    # Simulated scal map
    simulated_scalp_map = np.zeros(picks.shape[0])
    simulated_scalp_map[27:34] = 1

    simulated_scalp_map = simulated_scalp_map[:, None]

    evoked = EvokedArray(simulated_scalp_map, info, tmin=0)

    min_dist = 5.  # distance in mm

    dip, residual = fit_dipole(evoked, cov, fname_bem, fname_trans,
                               min_dist=min_dist)

    dist = _compute_depth(dip, fname_bem, fname_trans, subject, subjects_dir)

    assert_true(min_dist < (dist[0] * 1000.) < (min_dist + 1.))

    assert_raises(ValueError, fit_dipole, evoked, cov, fname_bem, fname_trans,
                  -1.)

def test_accuracy():
    """Test dipole fitting to sub-mm accuracy."""
    evoked = read_evokeds(fname_evo)[0].crop(0., 0.,)
    evoked.pick_types(meg=True, eeg=False)
    evoked.pick_channels([c for c in evoked.ch_names[::4]])
    for rad, perc_90 in zip((0.09, None), (0.002, 0.004)):
        bem = make_sphere_model('auto', rad, evoked.info,
                                relative_radii=(0.999, 0.998, 0.997, 0.995))
        src = read_source_spaces(fname_src)

        fwd = make_forward_solution(evoked.info, None, src, bem)
        fwd = convert_forward_solution(fwd, force_fixed=True, use_cps=True)
        vertices = [src[0]['vertno'], src[1]['vertno']]
        n_vertices = sum(len(v) for v in vertices)
        amp = 10e-9
        data = np.eye(n_vertices + 1)[:n_vertices]
        data[-1, -1] = 1.
        data *= amp
        stc = SourceEstimate(data, vertices, 0., 1e-3, 'sample')
        evoked.info.normalize_proj()
        sim = simulate_evoked(fwd, stc, evoked.info, cov=None, nave=np.inf)

        cov = make_ad_hoc_cov(evoked.info)
        dip = fit_dipole(sim, cov, bem, min_dist=0.001)[0]

        ds = []
        for vi in range(n_vertices):
            if vi < len(vertices[0]):
                hi = 0
                vertno = vi
            else:
                hi = 1
                vertno = vi - len(vertices[0])
            vertno = src[hi]['vertno'][vertno]
            rr = src[hi]['rr'][vertno]
            d = np.sqrt(np.sum((rr - dip.pos[vi]) ** 2))
            ds.append(d)
        # make sure that our median is sub-mm and the large majority are very
        # close (we expect some to be off by a bit e.g. because they are
        # radial)
        assert_true((np.percentile(ds, [50, 90]) < [0.0005, perc_90]).all())

def test_confidence(tmpdir):
    """Test confidence limits."""
    evoked = read_evokeds(fname_evo_full, 'Left Auditory', baseline=(None, 0))
    evoked.crop(0.08, 0.08).pick_types()  # MEG-only
    cov = make_ad_hoc_cov(evoked.info)
    sphere = make_sphere_model((0., 0., 0.04), 0.08)
    dip_py = fit_dipole(evoked, cov, sphere)[0]
    fname_test = op.join(str(tmpdir), 'temp-dip.txt')
    dip_py.save(fname_test)
    dip_read = read_dipole(fname_test)
    with pytest.warns(RuntimeWarning, match="'noise/ft/cm', 'prob'"):
        dip_xfit = read_dipole(fname_dip_xfit)
    for dip_check in (dip_py, dip_read):
        assert_allclose(dip_check.pos, dip_xfit.pos, atol=5e-4)  # < 0.5 mm
        assert_allclose(dip_check.gof, dip_xfit.gof, atol=5e-1)  # < 0.5%
        assert_array_equal(dip_check.nfree, dip_xfit.nfree)  # exact match
        assert_allclose(dip_check.khi2, dip_xfit.khi2, rtol=2e-2)  # 2% miss
        assert set(dip_check.conf.keys()) == set(dip_xfit.conf.keys())
        for key in sorted(dip_check.conf.keys()):
            assert_allclose(dip_check.conf[key], dip_xfit.conf[key],
                            rtol=1.5e-1, err_msg=key)

def test_rap_music_sphere():
    """Test RAP-MUSIC with real data, sphere model, MEG only."""
    evoked, noise_cov = _get_data(ch_decim=8)
    sphere = mne.make_sphere_model(r0=(0., 0., 0.04))
    src = mne.setup_volume_source_space(subject=None, pos=10.,
                                        sphere=(0.0, 0.0, 40, 65.0),
                                        mindist=5.0, exclude=0.0)
    forward = mne.make_forward_solution(evoked.info, trans=None, src=src,
                                        bem=sphere)

    dipoles = rap_music(evoked, forward, noise_cov, n_dipoles=2)
    # Test that there is one dipole on each hemisphere
    pos = np.array([dip.pos[0] for dip in dipoles])
    assert_equal(pos.shape, (2, 3))
    assert_equal((pos[:, 0] < 0).sum(), 1)
    assert_equal((pos[:, 0] > 0).sum(), 1)
    # Check the amplitude scale
    assert_true(1e-10 < dipoles[0].amplitude[0] < 1e-7)
    # Check the orientation
    dip_fit = mne.fit_dipole(evoked, noise_cov, sphere)[0]
    assert_true(np.max(np.abs(np.dot(dip_fit.ori, dipoles[0].ori[0]))) > 0.99)
    assert_true(np.max(np.abs(np.dot(dip_fit.ori, dipoles[1].ori[0]))) > 0.99)

def test_confidence():
    """Test confidence limits."""
    tempdir = _TempDir()
    evoked = read_evokeds(fname_evo_full, 'Left Auditory', baseline=(None, 0))
    evoked.crop(0.08, 0.08).pick_types()  # MEG-only
    cov = make_ad_hoc_cov(evoked.info)
    sphere = make_sphere_model((0., 0., 0.04), 0.08)
    dip_py = fit_dipole(evoked, cov, sphere)[0]
    fname_test = op.join(tempdir, 'temp-dip.txt')
    dip_py.save(fname_test)
    dip_read = read_dipole(fname_test)
    with warnings.catch_warnings(record=True) as w:
        dip_xfit = read_dipole(fname_dip_xfit)
    assert_equal(len(w), 1)
    assert_true("['noise/ft/cm', 'prob']" in str(w[0].message))
    for dip_check in (dip_py, dip_read):
        assert_allclose(dip_check.pos, dip_xfit.pos, atol=5e-4)  # < 0.5 mm
        assert_allclose(dip_check.gof, dip_xfit.gof, atol=5e-1)  # < 0.5%
        assert_array_equal(dip_check.nfree, dip_xfit.nfree)  # exact match
        assert_allclose(dip_check.khi2, dip_xfit.khi2, rtol=2e-2)  # 2% miss
        assert_equal(set(dip_check.conf.keys()), set(dip_xfit.conf.keys()))
        for key in sorted(dip_check.conf.keys()):
            assert_allclose(dip_check.conf[key], dip_xfit.conf[key],
                            rtol=1.5e-1, err_msg=key)

def test_dipole_fitting():
    """Test dipole fitting"""
    amp = 10e-9
    tempdir = _TempDir()
    rng = np.random.RandomState(0)
    fname_dtemp = op.join(tempdir, 'test.dip')
    fname_sim = op.join(tempdir, 'test-ave.fif')
    fwd = convert_forward_solution(read_forward_solution(fname_fwd),
                                   surf_ori=False, force_fixed=True)
    evoked = read_evokeds(fname_evo)[0]
    cov = read_cov(fname_cov)
    n_per_hemi = 5
    vertices = [np.sort(rng.permutation(s['vertno'])[:n_per_hemi])
                for s in fwd['src']]
    nv = sum(len(v) for v in vertices)
    stc = SourceEstimate(amp * np.eye(nv), vertices, 0, 0.001)
    with warnings.catch_warnings(record=True):  # semi-def cov
        evoked = generate_evoked(fwd, stc, evoked, cov, snr=20,
                                 random_state=rng)
    # For speed, let's use a subset of channels (strange but works)
    picks = np.sort(np.concatenate([
        pick_types(evoked.info, meg=True, eeg=False)[::2],
        pick_types(evoked.info, meg=False, eeg=True)[::2]]))
    evoked.pick_channels([evoked.ch_names[p] for p in picks])
    evoked.add_proj(make_eeg_average_ref_proj(evoked.info))
    write_evokeds(fname_sim, evoked)

    # Run MNE-C version
    run_subprocess([
        'mne_dipole_fit', '--meas', fname_sim, '--meg', '--eeg',
        '--noise', fname_cov, '--dip', fname_dtemp,
        '--mri', fname_fwd, '--reg', '0', '--tmin', '0',
    ])
    dip_c = read_dipole(fname_dtemp)

    # Run mne-python version
    sphere = make_sphere_model(head_radius=0.1)
    dip, residuals = fit_dipole(evoked, fname_cov, sphere, fname_fwd)

    # Sanity check: do our residuals have less power than orig data?
    data_rms = np.sqrt(np.sum(evoked.data ** 2, axis=0))
    resi_rms = np.sqrt(np.sum(residuals ** 2, axis=0))
    assert_true((data_rms > resi_rms).all())

    # Compare to original points
    transform_surface_to(fwd['src'][0], 'head', fwd['mri_head_t'])
    transform_surface_to(fwd['src'][1], 'head', fwd['mri_head_t'])
    src_rr = np.concatenate([s['rr'][v] for s, v in zip(fwd['src'], vertices)],
                            axis=0)
    src_nn = np.concatenate([s['nn'][v] for s, v in zip(fwd['src'], vertices)],
                            axis=0)

    # MNE-C skips the last "time" point :(
    dip.crop(dip_c.times[0], dip_c.times[-1])
    src_rr, src_nn = src_rr[:-1], src_nn[:-1]

    # check that we did at least as well
    corrs, dists, gc_dists, amp_errs, gofs = [], [], [], [], []
    for d in (dip_c, dip):
        new = d.pos
        diffs = new - src_rr
        corrs += [np.corrcoef(src_rr.ravel(), new.ravel())[0, 1]]
        dists += [np.sqrt(np.mean(np.sum(diffs * diffs, axis=1)))]
        gc_dists += [180 / np.pi * np.mean(np.arccos(np.sum(src_nn * d.ori,
                                                     axis=1)))]
        amp_errs += [np.sqrt(np.mean((amp - d.amplitude) ** 2))]
        gofs += [np.mean(d.gof)]
    assert_true(dists[0] >= dists[1], 'dists: %s' % dists)
    assert_true(corrs[0] <= corrs[1], 'corrs: %s' % corrs)
    assert_true(gc_dists[0] >= gc_dists[1], 'gc-dists (ori): %s' % gc_dists)
    assert_true(amp_errs[0] >= amp_errs[1], 'amplitude errors: %s' % amp_errs)

def test_mxne_vol_sphere():
    """Test (TF-)MxNE with a sphere forward and volumic source space."""
    evoked = read_evokeds(fname_data, condition=0, baseline=(None, 0))
    evoked.crop(tmin=-0.05, tmax=0.2)
    cov = read_cov(fname_cov)

    evoked_l21 = evoked.copy()
    evoked_l21.crop(tmin=0.081, tmax=0.1)

    info = evoked.info
    sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.080)
    src = mne.setup_volume_source_space(subject=None, pos=15., mri=None,
                                        sphere=(0.0, 0.0, 0.0, 80.0),
                                        bem=None, mindist=5.0,
                                        exclude=2.0)
    fwd = mne.make_forward_solution(info, trans=None, src=src,
                                    bem=sphere, eeg=False, meg=True)

    alpha = 80.
    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
                  loose=0.0, return_residual=False,
                  maxit=3, tol=1e-8, active_set_size=10)

    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
                  loose=0.2, return_residual=False,
                  maxit=3, tol=1e-8, active_set_size=10)

    # irMxNE tests
    stc = mixed_norm(evoked_l21, fwd, cov, alpha,
                     n_mxne_iter=1, maxit=30, tol=1e-8,
                     active_set_size=10)
    assert isinstance(stc, VolSourceEstimate)
    assert_array_almost_equal(stc.times, evoked_l21.times, 5)

    # Compare orientation obtained using fit_dipole and gamma_map
    # for a simulated evoked containing a single dipole
    stc = mne.VolSourceEstimate(50e-9 * np.random.RandomState(42).randn(1, 4),
                                vertices=stc.vertices[:1],
                                tmin=stc.tmin,
                                tstep=stc.tstep)
    evoked_dip = mne.simulation.simulate_evoked(fwd, stc, info, cov, nave=1e9,
                                                use_cps=True)

    dip_mxne = mixed_norm(evoked_dip, fwd, cov, alpha=80,
                          n_mxne_iter=1, maxit=30, tol=1e-8,
                          active_set_size=10, return_as_dipoles=True)

    amp_max = [np.max(d.amplitude) for d in dip_mxne]
    dip_mxne = dip_mxne[np.argmax(amp_max)]
    assert dip_mxne.pos[0] in src[0]['rr'][stc.vertices]

    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
    assert np.abs(np.dot(dip_fit.ori[0], dip_mxne.ori[0])) > 0.99

    # Do with TF-MxNE for test memory savings
    alpha = 60.  # overall regularization parameter
    l1_ratio = 0.01  # temporal regularization proportion

    stc, _ = tf_mixed_norm(evoked, fwd, cov, maxit=3, tol=1e-4,
                           tstep=16, wsize=32, window=0.1, alpha=alpha,
                           l1_ratio=l1_ratio, return_residual=True)
    assert isinstance(stc, VolSourceEstimate)
    assert_array_almost_equal(stc.times, evoked.times, 5)

def test_dipole_fitting():
    """Test dipole fitting."""
    amp = 100e-9
    tempdir = _TempDir()
    rng = np.random.RandomState(0)
    fname_dtemp = op.join(tempdir, 'test.dip')
    fname_sim = op.join(tempdir, 'test-ave.fif')
    fwd = convert_forward_solution(read_forward_solution(fname_fwd),
                                   surf_ori=False, force_fixed=True,
                                   use_cps=True)
    evoked = read_evokeds(fname_evo)[0]
    cov = read_cov(fname_cov)
    n_per_hemi = 5
    vertices = [np.sort(rng.permutation(s['vertno'])[:n_per_hemi])
                for s in fwd['src']]
    nv = sum(len(v) for v in vertices)
    stc = SourceEstimate(amp * np.eye(nv), vertices, 0, 0.001)
    evoked = simulate_evoked(fwd, stc, evoked.info, cov, nave=evoked.nave,
                             random_state=rng)
    # For speed, let's use a subset of channels (strange but works)
    picks = np.sort(np.concatenate([
        pick_types(evoked.info, meg=True, eeg=False)[::2],
        pick_types(evoked.info, meg=False, eeg=True)[::2]]))
    evoked.pick_channels([evoked.ch_names[p] for p in picks])
    evoked.add_proj(make_eeg_average_ref_proj(evoked.info))
    write_evokeds(fname_sim, evoked)

    # Run MNE-C version
    run_subprocess([
        'mne_dipole_fit', '--meas', fname_sim, '--meg', '--eeg',
        '--noise', fname_cov, '--dip', fname_dtemp,
        '--mri', fname_fwd, '--reg', '0', '--tmin', '0',
    ])
    dip_c = read_dipole(fname_dtemp)

    # Run mne-python version
    sphere = make_sphere_model(head_radius=0.1)
    with pytest.warns(RuntimeWarning, match='projection'):
        dip, residual = fit_dipole(evoked, cov, sphere, fname_fwd)
    assert isinstance(residual, Evoked)

    # Sanity check: do our residuals have less power than orig data?
    data_rms = np.sqrt(np.sum(evoked.data ** 2, axis=0))
    resi_rms = np.sqrt(np.sum(residual.data ** 2, axis=0))
    assert (data_rms > resi_rms * 0.95).all(), \
        '%s (factor: %s)' % ((data_rms / resi_rms).min(), 0.95)

    # Compare to original points
    transform_surface_to(fwd['src'][0], 'head', fwd['mri_head_t'])
    transform_surface_to(fwd['src'][1], 'head', fwd['mri_head_t'])
    assert_equal(fwd['src'][0]['coord_frame'], FIFF.FIFFV_COORD_HEAD)
    src_rr = np.concatenate([s['rr'][v] for s, v in zip(fwd['src'], vertices)],
                            axis=0)
    src_nn = np.concatenate([s['nn'][v] for s, v in zip(fwd['src'], vertices)],
                            axis=0)

    # MNE-C skips the last "time" point :(
    out = dip.crop(dip_c.times[0], dip_c.times[-1])
    assert (dip is out)
    src_rr, src_nn = src_rr[:-1], src_nn[:-1]

    # check that we did about as well
    corrs, dists, gc_dists, amp_errs, gofs = [], [], [], [], []
    for d in (dip_c, dip):
        new = d.pos
        diffs = new - src_rr
        corrs += [np.corrcoef(src_rr.ravel(), new.ravel())[0, 1]]
        dists += [np.sqrt(np.mean(np.sum(diffs * diffs, axis=1)))]
        gc_dists += [180 / np.pi * np.mean(np.arccos(np.sum(src_nn * d.ori,
                                                     axis=1)))]
        amp_errs += [np.sqrt(np.mean((amp - d.amplitude) ** 2))]
        gofs += [np.mean(d.gof)]
    # XXX possibly some OpenBLAS numerical differences make
    # things slightly worse for us
    factor = 0.7
    assert dists[0] / factor >= dists[1], 'dists: %s' % dists
    assert corrs[0] * factor <= corrs[1], 'corrs: %s' % corrs
    assert gc_dists[0] / factor >= gc_dists[1] * 0.8, \
        'gc-dists (ori): %s' % gc_dists
    assert amp_errs[0] / factor >= amp_errs[1],\
        'amplitude errors: %s' % amp_errs
    # This one is weird because our cov/sim/picking is weird
    assert gofs[0] * factor <= gofs[1] * 2, 'gof: %s' % gofs

#.........这里部分代码省略.........
                                        fname=fname, meg=True, eeg=False,
                                        mindist=5.0, n_jobs=2, overwrite=True)

    # for EEG only
    bem = join(subjects_dir, 'sample', 'bem',
               'sample-5120-5120-5120-bem-sol.fif')
    fname = join(meg_dir, 'sample_audvis-eeg-oct-6-fwd.fif')
    fwd_eeg = mne.make_forward_solution(raw.info, trans, src, bem,
                                        fname=fname, meg=False, eeg=True,
                                        mindist=5.0, n_jobs=2, overwrite=True)

    # for both EEG and MEG
    fname = join(meg_dir, 'sample_audvis-meg-eeg-oct-6-fwd.fif')
    fwd = mne.make_forward_solution(raw.info, trans, src, bem,
                                    fname=fname, meg=True, eeg=True,
                                    mindist=5.0, n_jobs=2, overwrite=True)

    # Create various sensitivity maps
    grad_map = mne.sensitivity_map(fwd, ch_type='grad', mode='free')
    grad_map.save(join(meg_dir, 'sample_audvis-grad-oct-6-fwd-sensmap'),
                  ftype='w')
    mag_map = mne.sensitivity_map(fwd, ch_type='mag', mode='free')
    mag_map.save(join(meg_dir, 'sample_audvis-mag-oct-6-fwd-sensmap'),
                 ftype='w')
    eeg_map = mne.sensitivity_map(fwd, ch_type='eeg', mode='free')
    eeg_map.save(join(meg_dir, 'sample_audvis-eeg-oct-6-fwd-sensmap'),
                 ftype='w')
    grad_map2 = mne.sensitivity_map(fwd, ch_type='grad', mode='fixed')
    grad_map2.save(join(meg_dir, 'sample_audvis-grad-oct-6-fwd-sensmap-2'),
                   ftype='w')
    mag_map2 = mne.sensitivity_map(fwd, ch_type='mag', mode='ratio')
    mag_map2.save(join(meg_dir, 'sample_audvis-mag-oct-6-fwd-sensmap-3'),
                  ftype='w')

    # Compute some with the EOG + ECG projectors
    projs = ecg_proj + eog_proj + raw.info['projs']
    for map_type in ['radiality', 'angle', 'remaining', 'dampening']:
        eeg_map = mne.sensitivity_map(fwd, projs=projs, ch_type='eeg',
                                      mode=map_type)
        eeg_map.save(join(meg_dir,
                          'sample_audvis-eeg-oct-6-fwd-sensmap-' + map_type))

###############################################################################
    # Compute MNE inverse operators
    #
    # Note: The MEG/EEG forward solution could be used for all
    #
    inv_meg = make_inverse_operator(raw.info, fwd_meg, noise_cov, loose=0.2)
    fname = join(meg_dir, 'sample_audvis-meg-oct-6-meg-inv.fif')
    write_inverse_operator(fname, inv_meg)

    inv_eeg = make_inverse_operator(raw.info, fwd_eeg, noise_cov, loose=0.2)
    fname = join(meg_dir, 'sample_audvis-eeg-oct-6-eeg-inv.fif')
    write_inverse_operator(fname, inv_eeg)

    inv = make_inverse_operator(raw.info, fwd, noise_cov, loose=0.2)
    fname = join(meg_dir, 'sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif')
    write_inverse_operator(fname, inv)

    # inverse operator with fixed orientation (for testing). Not implemented
    #inv_fixed = make_inverse_operator(raw.info, fwd_meg, noise_cov,
    #                                  depth=None, fixed=True)
    #fname = join(meg_dir, 'sample_audvis-meg-oct-6-meg-nodepth-fixed-inv.fif')
    #write_inverse_operator(fname, inv_fixed)

    # produce two with diagonal noise (for testing)
    diag = noise_cov.as_diag()
    inv_meg_diag = make_inverse_operator(raw.info, fwd_meg, diag, loose=0.2)
    fname = join(meg_dir, 'sample_audvis-meg-oct-6-meg-diagnoise-inv.fif')
    write_inverse_operator(fname, inv_meg_diag)

    inv_eeg_diag = make_inverse_operator(raw.info, fwd, diag, loose=0.2)
    fname = join(meg_dir,
                 'sample_audvis-meg-eeg-oct-6-meg-eeg-diagnoise-inv.fif')
    write_inverse_operator(fname, inv_eeg_diag)

    # Produce stc files
    evoked.crop(0, 0.25)
    stc_meg = apply_inverse(evoked, inv_meg, method='MNE')
    stc_meg.save(join(meg_dir, 'sample_audvis-meg'))
    stc_eeg = apply_inverse(evoked, inv_eeg, method='MNE')
    stc_eeg.save(join(meg_dir, 'sample_audvis-eeg'))
    stc = apply_inverse(evoked, inv, method='MNE')
    stc.save(join(meg_dir, 'sample_audvis-meg-eeg'))

    # let's also morph to fsaverage
    stc_to = stc_meg.morph('fsaverage', grade=3, smooth=12)
    stc_to.save(join(meg_dir, 'fsaverage_audvis-meg'))
    stc_to = stc_eeg.morph('fsaverage', grade=3, smooth=12)
    stc_to.save(join(meg_dir, 'fsaverage_audvis-eeg'))
    stc_to = stc.morph('fsaverage', grade=3, smooth=12)
    stc_to.save(join(meg_dir, 'fsaverage_audvis-meg-eeg'))

###############################################################################
    # Do one dipole fitting
    evoked = evoked.pick_types(meg=True, eeg=False)
    evoked.crop(0.04, 0.095)
    bem = join(subjects_dir, 'sample', 'bem', 'sample-5120-bem-sol.fif')
    dip, _ = mne.fit_dipole(evoked, noise_cov, bem, trans)
    dip.save(join(meg_dir, 'sample_audvis_set1.dip'))

###############################################################################
# Let's do some dipole fits. We first compute the noise covariance,
# then do the fits for each event_id taking the time instant that maximizes
# the global field power.

cov = mne.compute_covariance(epochs, tmax=0)
data = []
for ii in event_id:
    evoked = epochs[str(ii)].average()
    idx_peak = np.argmax(evoked.copy().pick_types(meg='grad').data.std(axis=0))
    t_peak = evoked.times[idx_peak]
    evoked.crop(t_peak, t_peak)
    data.append(evoked.data[:, 0])
evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.)
del epochs, raw
dip = fit_dipole(evoked, cov, sphere, n_jobs=1)[0]

###############################################################################
# Now we can compare to the actual locations, taking the difference in mm:

actual_pos, actual_ori = mne.dipole.get_phantom_dipoles()
actual_amp = 100.  # nAm

fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, figsize=(6, 7))

diffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))
print('mean(position error) = %s' % (np.mean(diffs),))
ax1.bar(event_id, diffs)
ax1.set_xlabel('Dipole index')
ax1.set_ylabel('Loc. error (mm)')

# here we can get away with using method='oas' for speed (faster than "shrunk")
# but in general "shrunk" is usually better
cov = mne.compute_covariance(
    epochs, tmax=bmax, method='oas', rank=None)
mne.viz.plot_evoked_white(epochs['1'].average(), cov)

data = []
t_peak = 0.036  # true for Elekta phantom
for ii in event_id:
    # Avoid the first and last trials -- can contain dipole-switching artifacts
    evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak)
    data.append(evoked.data[:, 0])
evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.)
del epochs
dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1)

###############################################################################
# Do a quick visualization of how much variance we explained, putting the
# data and residuals on the same scale (here the "time points" are the
# 32 dipole peak values that we fit):

fig, axes = plt.subplots(2, 1)
evoked.plot(axes=axes)
for ax in axes:
    ax.texts = []
    for line in ax.lines:
        line.set_color('#98df81')
residual.plot(axes=axes)

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

fwd = mne.read_forward_solution(fwd_fname)

with mne.forward.use_coil_def(coil_def_fname):
    mne.viz.plot_alignment(
        raw.info, trans, subject, subjects_dir, meg='sensors', bem=bem,
        surfaces=('head', 'pial'))
mlab.view(45, 60, distance=0.4, focalpoint=(0.02, 0, 0.04))

###############################################################################
# Perform dipole fitting
# ----------------------

# Fit dipoles on a subset of time points
with mne.forward.use_coil_def(coil_def_fname):
    dip_opm, _ = mne.fit_dipole(evoked.copy().crop(0.015, 0.080),
                                cov, bem, trans, verbose=True)
idx = np.argmax(dip_opm.gof)
print('Best dipole at t=%0.1f ms with %0.1f%% GOF'
      % (1000 * dip_opm.times[idx], dip_opm.gof[idx]))

# Plot N20m dipole as an example
dip_opm.plot_locations(trans, subject, subjects_dir,
                       mode='orthoview', idx=idx)

###############################################################################
# Perform minimum-norm localization
# ---------------------------------
# Due to the small number of sensors, there will be some leakage of activity
# to areas with low/no sensitivity. Constraining the source space to
# areas we are sensitive to might be a good idea.