Python inverse.apply_inverse函数代码示例

def _compare_inverses_approx(inv_1, inv_2, evoked, rtol, atol,
                             check_depth=True):
    # depth prior
    if check_depth:
        if inv_1['depth_prior'] is not None:
            assert_array_almost_equal(inv_1['depth_prior']['data'],
                                      inv_2['depth_prior']['data'], 5)
        else:
            assert_true(inv_2['depth_prior'] is None)
    # orient prior
    if inv_1['orient_prior'] is not None:
        assert_array_almost_equal(inv_1['orient_prior']['data'],
                                  inv_2['orient_prior']['data'])
    else:
        assert_true(inv_2['orient_prior'] is None)
    # source cov
    assert_array_almost_equal(inv_1['source_cov']['data'],
                              inv_2['source_cov']['data'])

    # These are not as close as we'd like XXX
    assert_array_almost_equal(np.abs(inv_1['eigen_fields']['data']),
                              np.abs(inv_2['eigen_fields']['data']), 0)
    assert_array_almost_equal(np.abs(inv_1['eigen_leads']['data']),
                              np.abs(inv_2['eigen_leads']['data']), 0)

    stc_1 = apply_inverse(evoked, inv_1, lambda2, "dSPM")
    stc_2 = apply_inverse(evoked, inv_2, lambda2, "dSPM")

    assert_true(stc_1.subject == stc_2.subject)
    assert_equal(stc_1.times, stc_2.times)
    assert_allclose(stc_1.data, stc_2.data, rtol=rtol, atol=atol)
    assert_true(inv_1['units'] == inv_2['units'])

def _compare_inverses_approx(inv_1, inv_2, evoked, stc_decimals,
                             check_depth=True):
    if check_depth:
        if inv_1['depth_prior'] is not None:
            assert_array_almost_equal(inv_1['depth_prior']['data'],
                                      inv_2['depth_prior']['data'])
        else:
            assert_true(inv_2['depth_prior'] is None)
    if inv_1['orient_prior'] is not None:
        assert_array_almost_equal(inv_1['orient_prior']['data'],
                                  inv_2['orient_prior']['data'])
    else:
        assert_true(inv_2['orient_prior'] is None)
    assert_array_almost_equal(inv_1['source_cov']['data'],
                              inv_2['source_cov']['data'])

    # These are not as close as we'd like XXX
    assert_array_almost_equal(np.abs(inv_1['eigen_fields']['data']),
                              np.abs(inv_2['eigen_fields']['data']), 0)
    assert_array_almost_equal(np.abs(inv_1['eigen_leads']['data']),
                              np.abs(inv_2['eigen_leads']['data']), 0)

    stc_1 = apply_inverse(evoked, inv_1, lambda2, "dSPM")
    stc_2 = apply_inverse(evoked, inv_2, lambda2, "dSPM")

    assert_equal(stc_1.times, stc_2.times)
    assert_array_almost_equal(stc_1.data, stc_2.data, stc_decimals)

def _compare_inverses_approx(inv_1, inv_2, evoked, rtol, atol,
                             depth_atol=1e-6, ctol=0.999999,
                             check_nn=True, check_K=True):
    """Compare inverses."""
    # depth prior
    if inv_1['depth_prior'] is not None:
        assert_allclose(inv_1['depth_prior']['data'],
                        inv_2['depth_prior']['data'], atol=depth_atol)
    else:
        assert_true(inv_2['depth_prior'] is None)
    # orient prior
    if inv_1['orient_prior'] is not None:
        assert_allclose(inv_1['orient_prior']['data'],
                        inv_2['orient_prior']['data'], atol=1e-7)
    else:
        assert_true(inv_2['orient_prior'] is None)
    # source cov
    assert_allclose(inv_1['source_cov']['data'], inv_2['source_cov']['data'],
                    atol=1e-7)
    for key in ('units', 'eigen_leads_weighted', 'nsource', 'coord_frame'):
        assert_equal(inv_1[key], inv_2[key], err_msg=key)
    assert_equal(inv_1['eigen_leads']['ncol'], inv_2['eigen_leads']['ncol'])
    K_1 = np.dot(inv_1['eigen_leads']['data'] * inv_1['sing'].astype(float),
                 inv_1['eigen_fields']['data'])
    K_2 = np.dot(inv_2['eigen_leads']['data'] * inv_2['sing'].astype(float),
                 inv_2['eigen_fields']['data'])
    # for free + surf ori, we only care about the ::2
    # (the other two dimensions have arbitrary direction)
    if inv_1['nsource'] * 3 == inv_1['source_nn'].shape[0]:
        # Technically this undersamples the free-orientation, non-surf-ori
        # inverse, but it's probably okay
        sl = slice(2, None, 3)
    else:
        sl = slice(None)
    if check_nn:
        assert_allclose(inv_1['source_nn'][sl], inv_2['source_nn'][sl],
                        atol=1e-4)
    if check_K:
        assert_allclose(np.abs(K_1[sl]), np.abs(K_2[sl]), rtol=rtol, atol=atol)

    # Now let's do some practical tests, too
    evoked = EvokedArray(np.eye(len(evoked.ch_names)), evoked.info)
    for method in ('MNE', 'dSPM'):
        stc_1 = apply_inverse(evoked, inv_1, lambda2, method)
        stc_2 = apply_inverse(evoked, inv_2, lambda2, method)
        assert_equal(stc_1.subject, stc_2.subject)
        assert_equal(stc_1.times, stc_2.times)
        stc_1 = stc_1.data
        stc_2 = stc_2.data
        norms = np.max(stc_1, axis=-1, keepdims=True)
        stc_1 /= norms
        stc_2 /= norms
        corr = np.corrcoef(stc_1.ravel(), stc_2.ravel())[0, 1]
        assert_true(corr > ctol, msg='%s < %s' % (corr, ctol))
        assert_allclose(stc_1, stc_2, rtol=rtol, atol=atol,
                        err_msg='%s: %s' % (method, corr))

def test_apply_inverse_operator():
    """Test MNE inverse computation (precomputed and non-precomputed)
    """
    inverse_operator = read_inverse_operator(fname_inv)
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    # Test old version of inverse computation starting from forward operator
    fwd_op = read_forward_solution(fname_fwd, surf_ori=True)
    my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov,
                                      loose=0.2, depth=0.8,
                                      limit_depth_chs=False)
    _compare_io(my_inv_op)
    assert_true(inverse_operator['units'] == 'Am')
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked, 2,
                             check_depth=False)
    # Inverse has 306 channels - 4 proj = 302
    assert_true(compute_rank_inverse(inverse_operator) == 302)

    # Test MNE inverse computation starting from forward operator
    my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov,
                                      loose=0.2, depth=0.8)
    _compare_io(my_inv_op)
    _compare_inverses_approx(my_inv_op, inverse_operator, evoked, 2)
    # Inverse has 306 channels - 4 proj = 302
    assert_true(compute_rank_inverse(inverse_operator) == 302)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "MNE")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 10e-10)
    assert_true(stc.data.mean() > 1e-11)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "sLORETA")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 10.0)
    assert_true(stc.data.mean() > 0.1)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "dSPM")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 35)
    assert_true(stc.data.mean() > 0.1)

    my_stc = apply_inverse(evoked, my_inv_op, lambda2, "dSPM")

    assert_true('dev_head_t' in my_inv_op['info'])
    assert_true('mri_head_t' in my_inv_op)

    assert_true(my_stc.subject == 'sample')
    assert_equal(stc.times, my_stc.times)
    assert_array_almost_equal(stc.data, my_stc.data, 2)

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 test_apply_inverse_sphere():
    """Test applying an inverse with a sphere model (rank-deficient)."""
    evoked = _get_evoked()
    evoked.pick_channels(evoked.ch_names[:306:8])
    evoked.info['projs'] = []
    cov = make_ad_hoc_cov(evoked.info)
    sphere = make_sphere_model('auto', 'auto', evoked.info)
    fwd = read_forward_solution(fname_fwd)
    vertices = [fwd['src'][0]['vertno'][::5],
                fwd['src'][1]['vertno'][::5]]
    stc = SourceEstimate(np.zeros((sum(len(v) for v in vertices), 1)),
                         vertices, 0., 1.)
    fwd = restrict_forward_to_stc(fwd, stc)
    fwd = make_forward_solution(evoked.info, fwd['mri_head_t'], fwd['src'],
                                sphere, mindist=5.)
    evoked = EvokedArray(fwd['sol']['data'].copy(), evoked.info)
    assert fwd['sol']['nrow'] == 39
    assert fwd['nsource'] == 101
    assert fwd['sol']['ncol'] == 303
    tempdir = _TempDir()
    temp_fname = op.join(tempdir, 'temp-inv.fif')
    inv = make_inverse_operator(evoked.info, fwd, cov, loose=1.)
    # This forces everything to be float32
    write_inverse_operator(temp_fname, inv)
    inv = read_inverse_operator(temp_fname)
    stc = apply_inverse(evoked, inv, method='eLORETA',
                        method_params=dict(eps=1e-2))
    # assert zero localization bias
    assert_array_equal(np.argmax(stc.data, axis=0),
                       np.repeat(np.arange(101), 3))

def test_apply_inverse_operator():
    """Test MNE inverse application
    """
    inverse_operator = read_inverse_operator(fname_full)
    evoked = _get_evoked()

    # Inverse has 306 channels - 4 proj = 302
    assert_true(compute_rank_inverse(inverse_operator) == 302)

    # Inverse has 306 channels - 4 proj = 302
    assert_true(compute_rank_inverse(inverse_operator) == 302)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "MNE")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 10e-9)
    assert_true(stc.data.mean() > 1e-11)

    # test if using prepared and not prepared inverse operator give the same
    # result
    inv_op = prepare_inverse_operator(inverse_operator, nave=evoked.nave,
                                      lambda2=lambda2, method="MNE")
    stc2 = apply_inverse(evoked, inv_op, lambda2, "MNE")
    assert_array_almost_equal(stc.data, stc2.data)
    assert_array_almost_equal(stc.times, stc2.times)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "sLORETA")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 10.0)
    assert_true(stc.data.mean() > 0.1)

    stc = apply_inverse(evoked, inverse_operator, lambda2, "dSPM")
    assert_true(stc.subject == 'sample')
    assert_true(stc.data.min() > 0)
    assert_true(stc.data.max() < 35)
    assert_true(stc.data.mean() > 0.1)

    # Test we get errors when using custom ref or no average proj is present
    evoked.info['custom_ref_applied'] = True
    assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE")
    evoked.info['custom_ref_applied'] = False
    evoked.info['projs'] = []  # remove EEG proj
    assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE")

    def test_apply_inverse_operator(self):
        """Test MNE inverse computation (precomputed and non-precomputed)
        """
        # Test old version of inverse computation starting from forward
        # operator
        my_inv_op = make_inverse_operator(evoked.info, self.fwd_op, noise_cov,
                                          loose=0.2, depth=0.8,
                                          limit_depth_chs=False)
        _compare_io(my_inv_op)
        _compare_inverses_approx(my_inv_op, self.inv_op, evoked, 2,
                                 check_depth=False)
        # Inverse has 306 channels - 4 proj = 302
        assert_true(compute_rank_inverse(self.inv_op) == 302)

        # Test MNE inverse computation starting from forward operator
        my_inv_op = make_inverse_operator(evoked.info, self.fwd_op, noise_cov,
                                          loose=0.2, depth=0.8)
        _compare_io(my_inv_op)
        _compare_inverses_approx(my_inv_op, self.inv_op, evoked, 2)
        # Inverse has 306 channels - 4 proj = 302
        assert_true(compute_rank_inverse(self.inv_op) == 302)

        stc = apply_inverse(evoked, self.inv_op, lambda2, "MNE")
        assert_true(stc.data.min() > 0)
        assert_true(stc.data.max() < 10e-10)
        assert_true(stc.data.mean() > 1e-11)

        stc = apply_inverse(evoked, self.inv_op, lambda2, "sLORETA")
        assert_true(stc.data.min() > 0)
        assert_true(stc.data.max() < 9.0)
        assert_true(stc.data.mean() > 0.1)

        stc = apply_inverse(evoked, self.inv_op, lambda2, "dSPM")
        assert_true(stc.data.min() > 0)
        assert_true(stc.data.max() < 35)
        assert_true(stc.data.mean() > 0.1)

        my_stc = apply_inverse(evoked, my_inv_op, lambda2, "dSPM")

        assert_true('dev_head_t' in my_inv_op['info'])
        assert_true('mri_head_t' in my_inv_op)

        assert_equal(stc.times, my_stc.times)
        assert_array_almost_equal(stc.data, my_stc.data, 2)

def test_inverse_operator_volume(evoked):
    """Test MNE inverse computation on volume source space."""
    tempdir = _TempDir()
    inv_vol = read_inverse_operator(fname_vol_inv)
    assert (repr(inv_vol))
    stc = apply_inverse(evoked, inv_vol, lambda2, 'dSPM')
    assert (isinstance(stc, VolSourceEstimate))
    # volume inverses don't have associated subject IDs
    assert (stc.subject is None)
    stc.save(op.join(tempdir, 'tmp-vl.stc'))
    stc2 = read_source_estimate(op.join(tempdir, 'tmp-vl.stc'))
    assert (np.all(stc.data > 0))
    assert (np.all(stc.data < 35))
    assert_array_almost_equal(stc.data, stc2.data)
    assert_array_almost_equal(stc.times, stc2.times)
    # vector source estimate
    stc_vec = apply_inverse(evoked, inv_vol, lambda2, 'dSPM', 'vector')
    assert (repr(stc_vec))
    assert_allclose(np.linalg.norm(stc_vec.data, axis=1), stc.data)

def test_inverse_operator_volume():
    """Test MNE inverse computation on volume source space"""
    evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0))
    inverse_operator_vol = read_inverse_operator(fname_vol_inv)
    stc = apply_inverse(evoked, inverse_operator_vol, lambda2, "dSPM")
    stc.save('tmp-vl.stc')
    stc2 = read_source_estimate('tmp-vl.stc')
    assert_true(np.all(stc.data > 0))
    assert_true(np.all(stc.data < 35))
    assert_array_almost_equal(stc.data, stc2.data)
    assert_array_almost_equal(stc.times, stc2.times)

def test_localization_bias_fixed(bias_params_fixed, method, lower, upper,
                                 depth):
    """Test inverse localization bias for fixed minimum-norm solvers."""
    evoked, fwd, noise_cov, _, want = bias_params_fixed
    fwd_use = convert_forward_solution(fwd, force_fixed=False)
    inv_fixed = make_inverse_operator(evoked.info, fwd_use, noise_cov,
                                      loose=0., depth=depth)
    loc = np.abs(apply_inverse(evoked, inv_fixed, lambda2, method).data)
    # Compute the percentage of sources for which there is no loc bias:
    perc = (want == np.argmax(loc, axis=0)).mean() * 100
    assert lower <= perc <= upper, method

def test_inverse_operator_volume():
    """Test MNE inverse computation on volume source space
    """
    inverse_operator_vol = read_inverse_operator(fname_vol_inv)
    _compare_io(inverse_operator_vol)
    stc = apply_inverse(evoked, inverse_operator_vol, lambda2, "dSPM")
    stc.save(op.join(tempdir, 'tmp-vl.stc'))
    stc2 = read_source_estimate(op.join(tempdir, 'tmp-vl.stc'))
    assert_true(np.all(stc.data > 0))
    assert_true(np.all(stc.data < 35))
    assert_array_almost_equal(stc.data, stc2.data)
    assert_array_almost_equal(stc.times, stc2.times)

def test_localization_bias_free(bias_params_free, method, lower, upper,
                                kwargs, depth):
    """Test inverse localization bias for free minimum-norm solvers."""
    evoked, fwd, noise_cov, _, want = bias_params_free
    inv_free = make_inverse_operator(evoked.info, fwd, noise_cov, loose=1.,
                                     depth=depth)
    loc = apply_inverse(evoked, inv_free, lambda2, method,
                        pick_ori='vector', **kwargs).data
    loc = np.linalg.norm(loc, axis=1)
    # Compute the percentage of sources for which there is no loc bias:
    perc = (want == np.argmax(loc, axis=0)).mean() * 100
    assert lower <= perc <= upper, method

def test_inverse_residual():
    """Test MNE inverse application."""
    # use fname_inv as it will be faster than fname_full (fewer verts and chs)
    evoked = _get_evoked().pick_types()
    inv = read_inverse_operator(fname_inv_fixed_depth)
    fwd = read_forward_solution(fname_fwd)
    fwd = convert_forward_solution(fwd, force_fixed=True, surf_ori=True)
    fwd = pick_channels_forward(fwd, evoked.ch_names)
    matcher = re.compile(r'.* ([0-9]?[0-9]?[0-9]?\.[0-9])% variance.*')
    for method in ('MNE', 'dSPM', 'sLORETA'):
        with catch_logging() as log:
            stc, residual = apply_inverse(
                evoked, inv, method=method, return_residual=True, verbose=True)
        log = log.getvalue()
        match = matcher.match(log.replace('\n', ' '))
        assert match is not None
        match = float(match.group(1))
        assert 45 < match < 50
        if method == 'MNE':  # must be first!
            recon = apply_forward(fwd, stc, evoked.info)
            proj_op = make_projector(evoked.info['projs'], evoked.ch_names)[0]
            recon.data[:] = np.dot(proj_op, recon.data)
            residual_fwd = evoked.copy()
            residual_fwd.data -= recon.data
        corr = np.corrcoef(residual_fwd.data.ravel(),
                           residual.data.ravel())[0, 1]
        assert corr > 0.999
    with catch_logging() as log:
        _, residual = apply_inverse(
            evoked, inv, 0., 'MNE', return_residual=True, verbose=True)
    log = log.getvalue()
    match = matcher.match(log.replace('\n', ' '))
    assert match is not None
    match = float(match.group(1))
    assert match == 100.
    assert_array_less(np.abs(residual.data), 1e-15)

    # Degenerate: we don't have the right representation for eLORETA for this
    with pytest.raises(ValueError, match='eLORETA does not .* support .*'):
        apply_inverse(evoked, inv, method="eLORETA", return_residual=True)

def test_localization_bias_loose(bias_params_fixed, method, lower, upper,
                                 depth):
    """Test inverse localization bias for loose minimum-norm solvers."""
    evoked, fwd, noise_cov, _, want = bias_params_fixed
    fwd = convert_forward_solution(fwd, surf_ori=False, force_fixed=False)
    assert not is_fixed_orient(fwd)
    inv_loose = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.2,
                                      depth=depth)
    loc = apply_inverse(evoked, inv_loose, lambda2, method).data
    assert (loc >= 0).all()
    # Compute the percentage of sources for which there is no loc bias:
    perc = (want == np.argmax(loc, axis=0)).mean() * 100
    assert lower <= perc <= upper, method