Python mne.convert_forward_solution函数代码示例

def _compare_forwards(fwd, fwd_py, n_sensors, n_src,
                      meg_rtol=1e-4, meg_atol=1e-9,
                      eeg_rtol=1e-3, eeg_atol=1e-3):
    """Helper to test forwards"""
    # check source spaces
    assert_equal(len(fwd['src']), len(fwd_py['src']))
    _compare_source_spaces(fwd['src'], fwd_py['src'], mode='approx')
    for surf_ori in [False, True]:
        if surf_ori:
            # use copy here to leave our originals unmodified
            fwd = convert_forward_solution(fwd, surf_ori, copy=True)
            fwd_py = convert_forward_solution(fwd, surf_ori, copy=True)

        for key in ['nchan', 'source_nn', 'source_rr', 'source_ori',
                    'surf_ori', 'coord_frame', 'nsource']:
            print(key)
            assert_allclose(fwd_py[key], fwd[key], rtol=1e-4, atol=1e-7)
        assert_allclose(fwd_py['mri_head_t']['trans'],
                        fwd['mri_head_t']['trans'], rtol=1e-5, atol=1e-8)

        assert_equal(fwd_py['sol']['data'].shape, (n_sensors, n_src))
        assert_equal(len(fwd['sol']['row_names']), n_sensors)
        assert_equal(len(fwd_py['sol']['row_names']), n_sensors)

        # check MEG
        assert_allclose(fwd['sol']['data'][:306],
                        fwd_py['sol']['data'][:306],
                        rtol=meg_rtol, atol=meg_atol,
                        err_msg='MEG mismatch')
        # check EEG
        if fwd['sol']['data'].shape[0] > 306:
            assert_allclose(fwd['sol']['data'][306:],
                            fwd_py['sol']['data'][306:],
                            rtol=eeg_rtol, atol=eeg_atol,
                            err_msg='EEG mismatch')

def _compare_forwards(fwd, fwd_py, n_sensors, n_src, meg_rtol=1e-4, meg_atol=1e-9):
    """Helper to test forwards"""
    # check source spaces
    assert_equal(len(fwd["src"]), len(fwd_py["src"]))
    _compare_source_spaces(fwd["src"], fwd_py["src"], mode="approx")
    for surf_ori in [False, True]:
        if surf_ori:
            # use copy here to leave our originals unmodified
            fwd = convert_forward_solution(fwd, surf_ori, copy=True)
            fwd_py = convert_forward_solution(fwd, surf_ori, copy=True)

        for key in ["nchan", "source_nn", "source_rr", "source_ori", "surf_ori", "coord_frame", "nsource"]:
            print(key)
            assert_allclose(fwd_py[key], fwd[key], rtol=1e-4, atol=1e-7)
        assert_allclose(fwd_py["mri_head_t"]["trans"], fwd["mri_head_t"]["trans"], rtol=1e-5, atol=1e-8)

        assert_equal(fwd_py["sol"]["data"].shape, (n_sensors, n_src))
        assert_equal(len(fwd["sol"]["row_names"]), n_sensors)
        assert_equal(len(fwd_py["sol"]["row_names"]), n_sensors)

        # check MEG
        print("check MEG")
        assert_allclose(fwd["sol"]["data"][:306], fwd_py["sol"]["data"][:306], rtol=meg_rtol, atol=meg_atol)
        # check EEG
        if fwd["sol"]["data"].shape[0] > 306:
            print("check EEG")
            assert_allclose(fwd["sol"]["data"][306:], fwd_py["sol"]["data"][306:], rtol=1e-3, atol=1e-3)

def test_convert_forward():
    """Test converting forward solution between different representations
    """
    fwd = read_forward_solution(fname_meeg)
    print(fwd)  # __repr__
    assert_true(isinstance(fwd, Forward))
    # look at surface orientation
    fwd_surf = convert_forward_solution(fwd, surf_ori=True)
    fwd_surf_io = read_forward_solution(fname_meeg, surf_ori=True)
    compare_forwards(fwd_surf, fwd_surf_io)
    # go back
    fwd_new = convert_forward_solution(fwd_surf, surf_ori=False)
    print(fwd_new)
    assert_true(isinstance(fwd, Forward))
    compare_forwards(fwd, fwd_new)
    # now go to fixed
    fwd_fixed = convert_forward_solution(fwd_surf, surf_ori=False,
                                         force_fixed=True)
    print(fwd_fixed)
    assert_true(isinstance(fwd_fixed, Forward))
    fwd_fixed_io = read_forward_solution(fname_meeg, surf_ori=False,
                                         force_fixed=True)
    compare_forwards(fwd_fixed, fwd_fixed_io)
    # now go back to cartesian (original condition)
    fwd_new = convert_forward_solution(fwd_fixed)
    print(fwd_new)
    assert_true(isinstance(fwd_new, Forward))
    compare_forwards(fwd, fwd_new)

def test_convert_forward():
    """Test converting forward solution between different representations
    """
    fwd = read_forward_solution(fname_meeg_grad)
    fwd_repr = repr(fwd)
    assert_true('306' in fwd_repr)
    assert_true('60' in fwd_repr)
    assert_true(fwd_repr)
    assert_true(isinstance(fwd, Forward))
    # look at surface orientation
    fwd_surf = convert_forward_solution(fwd, surf_ori=True)
    # go back
    fwd_new = convert_forward_solution(fwd_surf, surf_ori=False)
    assert_true(repr(fwd_new))
    assert_true(isinstance(fwd_new, Forward))
    compare_forwards(fwd, fwd_new)
    del fwd_new
    gc.collect()

    # now go to fixed
    fwd_fixed = convert_forward_solution(fwd_surf, surf_ori=True,
                                         force_fixed=True, use_cps=False)
    del fwd_surf
    gc.collect()
    assert_true(repr(fwd_fixed))
    assert_true(isinstance(fwd_fixed, Forward))
    assert_true(is_fixed_orient(fwd_fixed))
    # now go back to cartesian (original condition)
    fwd_new = convert_forward_solution(fwd_fixed, surf_ori=False,
                                       force_fixed=False)
    assert_true(repr(fwd_new))
    assert_true(isinstance(fwd_new, Forward))
    compare_forwards(fwd, fwd_new)
    del fwd, fwd_new, fwd_fixed
    gc.collect()

def test_priors():
    """Test prior computations."""
    # Depth prior
    fwd = read_forward_solution(fname_meeg)
    assert not is_fixed_orient(fwd)
    n_sources = fwd['nsource']
    info = read_info(fname_evoked)
    depth_prior = compute_depth_prior(fwd, info, exp=0.8)
    assert depth_prior.shape == (3 * n_sources,)
    depth_prior = compute_depth_prior(fwd, info, exp=0.)
    assert_array_equal(depth_prior, 1.)
    with pytest.raises(ValueError, match='must be "whiten"'):
        compute_depth_prior(fwd, info, limit_depth_chs='foo')
    with pytest.raises(ValueError, match='noise_cov must be a Covariance'):
        compute_depth_prior(fwd, info, limit_depth_chs='whiten')
    fwd_fixed = convert_forward_solution(fwd, force_fixed=True)
    with pytest.deprecated_call():
        depth_prior = compute_depth_prior(
            fwd_fixed['sol']['data'], info, is_fixed_ori=True)
    assert depth_prior.shape == (n_sources,)
    # Orientation prior
    orient_prior = compute_orient_prior(fwd, 1.)
    assert_array_equal(orient_prior, 1.)
    orient_prior = compute_orient_prior(fwd_fixed, 0.)
    assert_array_equal(orient_prior, 1.)
    with pytest.raises(ValueError, match='oriented in surface coordinates'):
        compute_orient_prior(fwd, 0.5)
    fwd_surf_ori = convert_forward_solution(fwd, surf_ori=True)
    orient_prior = compute_orient_prior(fwd_surf_ori, 0.5)
    assert all(np.in1d(orient_prior, (0.5, 1.)))
    with pytest.raises(ValueError, match='between 0 and 1'):
        compute_orient_prior(fwd_surf_ori, -0.5)
    with pytest.raises(ValueError, match='with fixed orientation'):
        compute_orient_prior(fwd_fixed, 0.5)

def _compare_forwards(fwd, fwd_py, n_sensors, n_src):
    """Helper to test forwards"""
    # check source spaces
    assert_equal(len(fwd['src']), len(fwd_py['src']))
    _compare_source_spaces(fwd['src'], fwd_py['src'], mode='approx')
    for surf_ori in [False, True]:
        if surf_ori:
            fwd = convert_forward_solution(fwd, surf_ori, copy=False)
            fwd_py = convert_forward_solution(fwd, surf_ori, copy=False)

        for key in ['nchan', 'source_nn', 'source_rr', 'source_ori',
                    'surf_ori', 'coord_frame', 'nsource']:
            print key
            assert_allclose(fwd_py[key], fwd[key], rtol=1e-4, atol=1e-7)
        assert_allclose(fwd_py['mri_head_t']['trans'],
                        fwd['mri_head_t']['trans'], rtol=1e-5, atol=1e-8)

        # check MEG
        assert_allclose(fwd['sol']['data'][:306],
                        fwd_py['sol']['data'][:306],
                        rtol=1e-4, atol=1e-9)
        # check EEG
        if fwd['sol']['data'].shape[0] > 306:
            assert_allclose(fwd['sol']['data'][306:],
                            fwd_py['sol']['data'][306:],
                            rtol=1e-3, atol=1e-3)
        assert_equal(fwd_py['sol']['data'].shape, (n_sensors, n_src))
        assert_equal(len(fwd['sol']['row_names']), n_sensors)
        assert_equal(len(fwd_py['sol']['row_names']), n_sensors)

def test_restrict_forward_to_stc():
    """Test restriction of source space to source SourceEstimate
    """
    start = 0
    stop = 5
    n_times = stop - start - 1
    sfreq = 10.0
    t_start = 0.123

    fwd = read_forward_solution(fname_meeg)
    fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True,
                                   use_cps=True)
    fwd = pick_types_forward(fwd, meg=True)

    vertno = [fwd['src'][0]['vertno'][0:15], fwd['src'][1]['vertno'][0:5]]
    stc_data = np.ones((len(vertno[0]) + len(vertno[1]), n_times))
    stc = SourceEstimate(stc_data, vertno, tmin=t_start, tstep=1.0 / sfreq)

    fwd_out = restrict_forward_to_stc(fwd, stc)
    assert_true(isinstance(fwd_out, Forward))

    assert_equal(fwd_out['sol']['ncol'], 20)
    assert_equal(fwd_out['src'][0]['nuse'], 15)
    assert_equal(fwd_out['src'][1]['nuse'], 5)
    assert_equal(fwd_out['src'][0]['vertno'], fwd['src'][0]['vertno'][0:15])
    assert_equal(fwd_out['src'][1]['vertno'], fwd['src'][1]['vertno'][0:5])

    fwd = read_forward_solution(fname_meeg)
    fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=False)
    fwd = pick_types_forward(fwd, meg=True)

    vertno = [fwd['src'][0]['vertno'][0:15], fwd['src'][1]['vertno'][0:5]]
    stc_data = np.ones((len(vertno[0]) + len(vertno[1]), n_times))
    stc = SourceEstimate(stc_data, vertno, tmin=t_start, tstep=1.0 / sfreq)

    fwd_out = restrict_forward_to_stc(fwd, stc)

    assert_equal(fwd_out['sol']['ncol'], 60)
    assert_equal(fwd_out['src'][0]['nuse'], 15)
    assert_equal(fwd_out['src'][1]['nuse'], 5)
    assert_equal(fwd_out['src'][0]['vertno'], fwd['src'][0]['vertno'][0:15])
    assert_equal(fwd_out['src'][1]['vertno'], fwd['src'][1]['vertno'][0:5])

    # Test saving the restricted forward object. This only works if all fields
    # are properly accounted for.
    temp_dir = _TempDir()
    fname_copy = op.join(temp_dir, 'copy-fwd.fif')
    with warnings.catch_warnings(record=True):
        warnings.simplefilter('always')
        write_forward_solution(fname_copy, fwd_out, overwrite=True)
    fwd_out_read = read_forward_solution(fname_copy)
    fwd_out_read = convert_forward_solution(fwd_out_read, surf_ori=True,
                                            force_fixed=False)
    compare_forwards(fwd_out, fwd_out_read)

def test_make_forward_solution_discrete():
    """Test making and converting a forward solution with discrete src."""
    # smoke test for depth weighting and discrete source spaces
    src = read_source_spaces(fname_src)[0]
    src = SourceSpaces([src] + setup_volume_source_space(
        pos=dict(rr=src['rr'][src['vertno'][:3]].copy(),
                 nn=src['nn'][src['vertno'][:3]].copy())))
    sphere = make_sphere_model()
    fwd = make_forward_solution(fname_raw, fname_trans, src, sphere,
                                meg=True, eeg=False)
    convert_forward_solution(fwd, surf_ori=True)

def test_inverse_operator_channel_ordering():
    """Test MNE inverse computation is immune to channel reorderings
    """
    # These are with original ordering
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    fwd_orig = make_forward_solution(evoked.info, fname_trans, src_fname,
                                     fname_bem, eeg=True, mindist=5.0)
    fwd_orig = convert_forward_solution(fwd_orig, surf_ori=True)
    inv_orig = make_inverse_operator(evoked.info, fwd_orig, noise_cov,
                                     loose=0.2, depth=0.8,
                                     limit_depth_chs=False)
    stc_1 = apply_inverse(evoked, inv_orig, lambda2, "dSPM")

    # Assume that a raw reordering applies to both evoked and noise_cov,
    # so we don't need to create those from scratch. Just reorder them,
    # then try to apply the original inverse operator
    new_order = np.arange(len(evoked.info['ch_names']))
    randomiser = np.random.RandomState(42)
    randomiser.shuffle(new_order)
    evoked.data = evoked.data[new_order]
    evoked.info['chs'] = [evoked.info['chs'][n] for n in new_order]
    evoked.info._update_redundant()
    evoked.info._check_consistency()

    cov_ch_reorder = [c for c in evoked.info['ch_names']
                      if (c in noise_cov.ch_names)]

    new_order_cov = [noise_cov.ch_names.index(name) for name in cov_ch_reorder]
    noise_cov['data'] = noise_cov.data[np.ix_(new_order_cov, new_order_cov)]
    noise_cov['names'] = [noise_cov['names'][idx] for idx in new_order_cov]

    fwd_reorder = make_forward_solution(evoked.info, fname_trans, src_fname,
                                        fname_bem, eeg=True, mindist=5.0)
    fwd_reorder = convert_forward_solution(fwd_reorder, surf_ori=True)
    inv_reorder = make_inverse_operator(evoked.info, fwd_reorder, noise_cov,
                                        loose=0.2, depth=0.8,
                                        limit_depth_chs=False)

    stc_2 = apply_inverse(evoked, inv_reorder, lambda2, "dSPM")

    assert_equal(stc_1.subject, stc_2.subject)
    assert_array_equal(stc_1.times, stc_2.times)
    assert_allclose(stc_1.data, stc_2.data, rtol=1e-5, atol=1e-5)
    assert_true(inv_orig['units'] == inv_reorder['units'])

    # Reload with original ordering & apply reordered inverse
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    stc_3 = apply_inverse(evoked, inv_reorder, lambda2, "dSPM")
    assert_allclose(stc_1.data, stc_3.data, rtol=1e-5, atol=1e-5)

def _load_forward():
    """Load forward models."""
    fwd_free = mne.read_forward_solution(fname_fwd)
    fwd_free = mne.pick_types_forward(fwd_free, meg=True, eeg=False)
    fwd_free = mne.convert_forward_solution(fwd_free, surf_ori=False)
    fwd_surf = mne.convert_forward_solution(fwd_free, surf_ori=True,
                                            use_cps=False)
    fwd_fixed = mne.convert_forward_solution(fwd_free, force_fixed=True,
                                             use_cps=False)
    fwd_vol = mne.read_forward_solution(fname_fwd_vol)
    label = mne.read_label(fname_label)

    return fwd_free, fwd_surf, fwd_fixed, fwd_vol, label

def test_make_inverse_operator_loose(evoked):
    """Test MNE inverse computation (precomputed and non-precomputed)."""
    # Test old version of inverse computation starting from forward operator
    noise_cov = read_cov(fname_cov)
    inverse_operator = read_inverse_operator(fname_inv)
    fwd_op = convert_forward_solution(read_forward_solution_meg(fname_fwd),
                                      surf_ori=True, copy=False)
    with catch_logging() as log:
        with pytest.deprecated_call():  # limit_depth_chs
            my_inv_op = make_inverse_operator(
                evoked.info, fwd_op, noise_cov, loose=0.2, depth=0.8,
                limit_depth_chs=False, verbose=True)
    log = log.getvalue()
    assert 'MEG: rank 302 computed' in log
    assert 'limit = 1/%d' % fwd_op['nsource'] in log
    _compare_io(my_inv_op)
    assert_equal(inverse_operator['units'], 'Am')
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked,
                             rtol=1e-2, atol=1e-5, depth_atol=1e-3)
    # Test MNE inverse computation starting from forward operator
    with catch_logging() as log:
        my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov,
                                          loose='auto', depth=0.8,
                                          fixed=False, verbose=True)
    log = log.getvalue()
    assert 'MEG: rank 302 computed from 305' in log
    _compare_io(my_inv_op)
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked,
                             rtol=1e-3, atol=1e-5)
    assert ('dev_head_t' in my_inv_op['info'])
    assert ('mri_head_t' in my_inv_op)

def test_make_inverse_operator():
    """Test MNE inverse computation (precomputed and non-precomputed)."""
    # Test old version of inverse computation starting from forward operator
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)
    inverse_operator = read_inverse_operator(fname_inv)
    fwd_op = convert_forward_solution(read_forward_solution_meg(fname_fwd),
                                      surf_ori=True, copy=False)
    with catch_logging() as log:
        my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov,
                                          loose=0.2, depth=0.8,
                                          limit_depth_chs=False, verbose=True)
    log = log.getvalue()
    assert 'rank 302 (3 small eigenvalues omitted)' in log
    _compare_io(my_inv_op)
    assert_equal(inverse_operator['units'], 'Am')
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked,
                             rtol=1e-2, atol=1e-5, depth_atol=1e-3)
    # Test MNE inverse computation starting from forward operator
    with catch_logging() as log:
        my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov,
                                          loose=0.2, depth=0.8, verbose=True)
    log = log.getvalue()
    assert 'rank 302 (3 small eigenvalues omitted)' in log
    _compare_io(my_inv_op)
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked,
                             rtol=1e-3, atol=1e-5)
    assert ('dev_head_t' in my_inv_op['info'])
    assert ('mri_head_t' in my_inv_op)

def test_gamma_map():
    """Test Gamma MAP inverse"""
    forward = read_forward_solution(fname_fwd)
    forward = convert_forward_solution(forward, surf_ori=True)

    forward = pick_types_forward(forward, meg=False, eeg=True)
    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)

    alpha = 0.5
    stc = gamma_map(evoked, forward, cov, alpha, tol=1e-4,
                    xyz_same_gamma=True, update_mode=1)
    _check_stc(stc, evoked, 68477)

    stc = gamma_map(evoked, forward, cov, alpha, tol=1e-4,
                    xyz_same_gamma=False, update_mode=1)
    _check_stc(stc, evoked, 82010)

    dips = gamma_map(evoked, forward, cov, alpha, tol=1e-4,
                     xyz_same_gamma=False, update_mode=1,
                     return_as_dipoles=True)
    assert_true(isinstance(dips[0], Dipole))
    stc_dip = make_stc_from_dipoles(dips, forward['src'])
    _check_stcs(stc, stc_dip)

    # force fixed orientation
    stc = gamma_map(evoked, forward, cov, alpha, tol=1e-4,
                    xyz_same_gamma=False, update_mode=2,
                    loose=0, return_residual=False)
    _check_stc(stc, evoked, 85739, 20)

def _get_fwd_labels():
    fwd = read_forward_solution(fname_fwd)
    fwd = convert_forward_solution(fwd, force_fixed=True, use_cps=True)
    fwd = pick_types_forward(fwd, meg=True, eeg=False)
    labels = [read_label(op.join(data_path, 'MEG', 'sample', 'labels',
                         '%s.label' % label)) for label in label_names]
    return fwd, labels

def run_forward(subject_id):
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)
    data_path = op.join(meg_dir, subject)

    fname_ave = op.join(data_path, '%s-ave.fif' % subject)
    fname_fwd = op.join(data_path, '%s-meg-%s-fwd.fif' % (subject, spacing))
    fname_trans = op.join(study_path, 'ds117', subject, 'MEG', '%s-trans.fif' % subject)

    src = mne.setup_source_space(subject, spacing=spacing,
                                 subjects_dir=subjects_dir, overwrite=True,
                                 n_jobs=1, add_dist=False)

    src_fname = op.join(subjects_dir, subject, '%s-src.fif' % spacing)
    mne.write_source_spaces(src_fname, src)

    bem_model = mne.make_bem_model(subject, ico=4, subjects_dir=subjects_dir,
                                   conductivity=(0.3,))
    bem = mne.make_bem_solution(bem_model)
    info = mne.read_evokeds(fname_ave, condition=0).info
    fwd = mne.make_forward_solution(info, trans=fname_trans, src=src, bem=bem,
                                    fname=None, meg=True, eeg=False,
                                    mindist=mindist, n_jobs=1, overwrite=True)
    fwd = mne.convert_forward_solution(fwd, surf_ori=True)
    mne.write_forward_solution(fname_fwd, fwd, overwrite=True)