Python minimum_norm.read_inverse_operator函数代码示例

def calc_sub_cortical_activity(events_id, evoked_fn, inv_fn, sub_corticals_codes_file, baseline_min_t=None,
        baseline_max_t = 0, snr = 3.0, inverse_method='dSPM'):

    sub_corticals = read_sub_corticals_code_file(sub_corticals_codes_file)
    if len(sub_corticals) == 0:
        return

    lambda2 = 1.0 / snr ** 2
    lut = read_freesurfer_lookup_table(FREE_SURFER_HOME)
    for cond in events_id.keys():
        evo = evoked_fn.format(cond=cond)
        evoked = {event:mne.read_evokeds(evo, baseline=(baseline_min_t, baseline_max_t))[0] for event in [event]}
        inverse_operator = read_inverse_operator(inv_fn.format(cond=cond))
        stc = apply_inverse(evoked[cond], inverse_operator, lambda2, inverse_method)
        read_vertices_from = len(stc.vertices[0])+len(stc.vertices[1])
        sub_corticals_activity = {}
        for sub_cortical_ind, sub_cortical_code in enumerate(sub_corticals):
            # +2 becasue the first two are the hemispheres
            sub_corticals_activity[sub_cortical_code] = stc.data[
                read_vertices_from: read_vertices_from + len(stc.vertices[sub_cortical_ind + 2])]
            read_vertices_from += len(stc.vertices[sub_cortical_ind + 2])

        if not os.path.isdir:
            os.mkdir(os.path.join(SUBJECT_MEG_FOLDER, 'subcorticals'))
        for sub_cortical_code, activity in sub_corticals_activity.iteritems():
            sub_cortical, _ = get_numeric_index_to_label(sub_cortical_code, lut)
            np.save(os.path.join(SUBJECT_MEG_FOLDER, 'subcorticals', '{}-{}-{}'.format(cond, sub_cortical, inverse_method)), activity.mean(0))
            np.save(os.path.join(SUBJECT_MEG_FOLDER, 'subcorticals', '{}-{}-{}-all-vertices'.format(cond, sub_cortical, inverse_method)), activity)

def test_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces."""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    src_ = inverse_operator['src']
    with pytest.warns(RuntimeWarning, match='will have holes'):
        connectivity = spatio_temporal_src_connectivity(src_, n_times=2)
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert (a == b)

    assert_equal(grade_to_tris(5).shape, [40960, 3])

def test_epochs_vector_inverse():
    """Test vector inverse consistency between evoked and epochs."""
    raw = read_raw_fif(fname_raw)
    events = find_events(raw, stim_channel='STI 014')[:2]
    reject = dict(grad=2000e-13, mag=4e-12, eog=150e-6)

    epochs = Epochs(raw, events, None, 0, 0.01, baseline=None,
                    reject=reject, preload=True)

    assert_equal(len(epochs), 2)

    evoked = epochs.average(picks=range(len(epochs.ch_names)))

    inv = read_inverse_operator(fname_inv)

    method = "MNE"
    snr = 3.
    lambda2 = 1. / snr ** 2

    stcs_epo = apply_inverse_epochs(epochs, inv, lambda2, method=method,
                                    pick_ori='vector', return_generator=False)
    stc_epo = np.mean(stcs_epo)

    stc_evo = apply_inverse(evoked, inv, lambda2, method=method,
                            pick_ori='vector')

    assert_allclose(stc_epo.data, stc_evo.data, rtol=1e-9, atol=0)

def test_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces"""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    with warnings.catch_warnings(record=True) as w:
        connectivity = spatio_temporal_src_connectivity(
                                            inverse_operator['src'], n_times=2)
        assert len(w) == 1
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert_true(a == b)

def test_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces."""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        src_ = inverse_operator['src']
        connectivity = spatio_temporal_src_connectivity(src_, n_times=2)
    assert_equal(len(w), 1)
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert_true(a == b)

    assert_equal(grade_to_tris(5).shape, [40960, 3])

def _create_stcs(params):
    subject, events_id, epochs, evoked, inv, inverse_method, baseline, apply_for_epochs, apply_SSP_projection_vectors, add_eeg_ref = params
    snr = 3.0
    lambda2 = 1.0 / snr ** 2

    local_inv_file_name = op.join(LOCAL_ROOT_DIR, 'inv', '{}_ecr_nTSSS_conflict-inv.fif'.format(subject))
    if inv is None and os.path.isfile(local_inv_file_name):
        inv = read_inverse_operator(local_inv_file_name)
    if inv is None:
        return
    print([s['vertno'] for s in inv['src']])
    for cond_name in events_id.keys():
        if not apply_for_epochs:
            local_stc_file_name = op.join(LOCAL_ROOT_DIR, 'stc', '{}_{}_{}'.format(subject, cond_name, inverse_method))
            if os.path.isfile('{}-lh.stc'.format(local_stc_file_name)):
                print('stc was already calculated for {}'.format(subject))
            else:
                if evoked is None:
                    evoked_cond = mne.read_evokeds(op.join(LOCAL_ROOT_DIR, 'evo', '{}_ecr_{}-ave.fif'.format(subject, cond_name)))
                else:
                    evoked_cond = evoked[cond_name]
                stcs = apply_inverse(evoked_cond, inv, lambda2, inverse_method, pick_ori=None)
                stcs.save(local_stc_file_name, ftype='h5')
        else:
            local_stc_template = op.join(LOCAL_ROOT_DIR, 'stc_epochs', '{}_{}_{}_{}'.format(subject, cond_name, '{epoch_ind}', inverse_method))
            if len(glob.glob(local_stc_template.format(epoch_ind='*') == 36)):
                print('stc was already calculated for {}'.format(subject))
            else:
                stcs = apply_inverse_epochs(epochs[cond_name], inv, lambda2, inverse_method, pick_ori=None, return_generator=True)
                for epoch_ind, stc in enumerate(stcs):
                    if not os.path.isfile(local_stc_template.format(epoch_ind=epoch_ind)):
                        stc.save(local_stc_template.format(epoch_ind=epoch_ind), ftype='h5')

def apply_STC_ave(fnevo, method='dSPM', snr=3.0):
    ''' Inverse evoked data into the source space.
        Parameter
        ---------
        fnevo: string or list
            The evoked file with ECG, EOG and environmental noise free.
        method:string
            Inverse method, 'MNE' or 'dSPM'
        snr: float
            Signal to noise ratio for inverse solution.
    '''
    #Get the default subjects_dir
    from mne.minimum_norm import apply_inverse, read_inverse_operator
    fnlist = get_files_from_list(fnevo)
    # loop across all filenames
    for fname in fnlist:
        name = os.path.basename(fname)
        fn_path = os.path.split(fname)[0]
        fn_stc = fname[:fname.rfind('-ave.fif')]
        # fn_inv = fname[:fname.rfind('-ave.fif')] + ',ave-inv.fif'
        subject = name.split('_')[0]
        fn_inv = fn_path + '/%s_fibp1-45,ave-inv.fif' % subject
        snr = snr
        lambda2 = 1.0 / snr ** 2
        # noise_cov = mne.read_cov(fn_cov)
        [evoked] = mne.read_evokeds(fname)
        evoked.pick_types(meg=True, ref_meg=False)
        inv = read_inverse_operator(fn_inv)
        stc = apply_inverse(evoked, inv, lambda2, method,
                            pick_ori='normal')
        stc.save(fn_stc)

def apply_STC_epo(fnepo, event, method='MNE', snr=1.0, min_subject='fsaverage',
                  subjects_dir=None):

    from mne import morph_data
    from mne.minimum_norm import read_inverse_operator, apply_inverse_epochs

    fnlist = get_files_from_list(fnepo)
    # loop across all filenames
    for fname in fnlist:
        fn_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        subject = name.split('_')[0]
        min_dir = subjects_dir + '/%s' %min_subject
        snr = snr
        lambda2 = 1.0 / snr ** 2
        stcs_path = min_dir + '/stcs/%s/%s/' % (subject,event)
        reset_directory(stcs_path)
        # fn_inv = fname[:fname.rfind('-ave.fif')] + ',ave-inv.fif'
        fn_inv = fn_path + '/%s_epo-inv.fif' %subject

        # noise_cov = mne.read_cov(fn_cov)
        epo = mne.read_epochs(fname)
        epo.pick_types(meg=True, ref_meg=False)
        inv = read_inverse_operator(fn_inv)
        stcs = apply_inverse_epochs(epo, inv, lambda2, method,
                            pick_ori='normal')
        s = 0
        while s < len(stcs):
            stc_morph = morph_data(subject, min_subject, stcs[s])
            stc_morph.save(stcs_path + '/trial%s_fsaverage'
                           % (str(s)), ftype='stc')
            s = s + 1

def test_plot_snr():
    """Test plotting SNR estimate."""
    import matplotlib.pyplot as plt
    inv = read_inverse_operator(inv_fname)
    evoked = read_evokeds(evoked_fname, baseline=(None, 0))[0]
    plot_snr_estimate(evoked, inv)
    plt.close('all')

def test_psf_ctf():
    """Test computation of PSFs and CTFs for linear estimators
    """

    inverse_operator = read_inverse_operator(fname_inv)
    forward = read_forward_solution(fname_fwd, force_fixed=False,
                                    surf_ori=True)
    forward = pick_types_forward(forward, meg=True, eeg=False)
    labels = [mne.read_label(ss) for ss in fname_label]

    method = 'MNE'
    n_svd_comp = 2

    # Test PSFs (then CTFs)
    for mode in ('sum', 'svd'):
        stc_psf, psf_ev = point_spread_function(inverse_operator,
                                                forward,
                                                method=method,
                                                labels=labels,
                                                lambda2=lambda2,
                                                pick_ori='normal',
                                                mode=mode,
                                                n_svd_comp=n_svd_comp)

        n_vert, n_samples = stc_psf.shape
        should_n_vert = (inverse_operator['src'][1]['vertno'].shape[0] +
                         inverse_operator['src'][0]['vertno'].shape[0])
        if mode == 'svd':
            should_n_samples = len(labels) * n_svd_comp + 1
        else:
            should_n_samples = len(labels) + 1

        assert_true(n_vert == should_n_vert)
        assert_true(n_samples == should_n_samples)

        n_chan, n_samples = psf_ev.data.shape
        assert_true(n_chan == forward['nchan'])

    forward = read_forward_solution(fname_fwd, force_fixed=True, surf_ori=True)
    forward = pick_types_forward(forward, meg=True, eeg=False)

    # Test CTFs
    for mode in ('sum', 'svd'):
        stc_ctf = cross_talk_function(inverse_operator, forward,
                                      labels, method=method,
                                      lambda2=lambda2,
                                      signed=False, mode=mode,
                                      n_svd_comp=n_svd_comp)

        n_vert, n_samples = stc_ctf.shape
        should_n_vert = (inverse_operator['src'][1]['vertno'].shape[0] +
                         inverse_operator['src'][0]['vertno'].shape[0])
        if mode == 'svd':
            should_n_samples = len(labels) * n_svd_comp + 1
        else:
            should_n_samples = len(labels) + 1

        assert_true(n_vert == should_n_vert)
        assert_true(n_samples == should_n_samples)

def plot_inverse_operator(subject, fnout_img=None, verbose=None):
    """"Reads in and plots the inverse solution."""

    # get name of inverse solution-file
    subjects_dir = os.environ.get('SUBJECTS_DIR')
    fname_dir_inv = os.path.join(subjects_dir,
                                 subject)

    for files in os.listdir(fname_dir_inv):
        if files.endswith(",cleaned_epochs_avg-7-src-meg-inv.fif"):
            fname_inv = os.path.join(fname_dir_inv,
                                     files)
            break

    try:
        fname_inv
    except NameError:
        print "ERROR: No forward solution found!"
        sys.exit()


    # read inverse solution
    inv = min_norm.read_inverse_operator(fname_inv)

    # print some information if desired
    if verbose is not None:
        print "Method: %s" % inv['methods']
        print "fMRI prior: %s" % inv['fmri_prior']
        print "Number of sources: %s" % inv['nsource']
        print "Number of channels: %s" % inv['nchan']


    # show 3D source space
    lh_points = inv['src'][0]['rr']
    lh_faces = inv['src'][0]['tris']
    rh_points = inv['src'][1]['rr']
    rh_faces = inv['src'][1]['tris']

    # create figure and plot results
    fig_inverse = mlab.figure(size=(1200, 1200),
                              bgcolor=(0, 0, 0))

    mlab.triangular_mesh(lh_points[:, 0],
                         lh_points[:, 1],
                         lh_points[:, 2],
                         lh_faces)

    mlab.triangular_mesh(rh_points[:, 0],
                         rh_points[:, 1],
                         rh_points[:, 2],
                         rh_faces)

    # save result
    if fnout_img is not None:
        mlab.savefig(fnout_img,
                     figure=fig_inverse,
                     size=(1200, 1200))

    mlab.close(all=True)

def calc_inv_kernel(fn_inv, method="dSPM", nave=1, snr=6.,
                    pick_ori="normal", verbose=None):

    """
    Interface for preparing the kernel of the inverse
    estimation.

        Parameters
        ----------
        fn_inv : String containing the filename
            of the inverse operator (must be a fif-file)
        method : string
            which source localization method should be used?
            MNE-based: "MNE" | "dSPM" | "sLORETA"
            default: method="dSPM"
        nave : number of averages used to regularize the solution
            default: nave=1
        snr : signal-to-noise ratio
            default: snr = 3.
        verbose : bool, str, int, or None
            If not None, override default verbose level
            (see mne.verbose).
            default: verbose=None
    """

    # -------------------------------------------
    # import necessary modules
    # -------------------------------------------
    import mne.minimum_norm as min_norm
    from mne.minimum_norm.inverse import _assemble_kernel
    import numpy as np

    # -------------------------------------------
    # estimate inverse kernel
    # -------------------------------------------
    # load inverse solution
    inv_operator = min_norm.read_inverse_operator(fn_inv, verbose=verbose)

    # set up the inverse according to the parameters
    lambda2      = 1. / snr ** 2.   # the regularization parameter.
    inv_operator = min_norm.prepare_inverse_operator(inv_operator, nave, lambda2, method)

    # estimate inverse kernel and noise normalization coefficient
    kernel, noise_norm, vertno = _assemble_kernel(inv_operator, None, method, pick_ori)

    if method == "MNE":
        noise_norm = np.ones((kernel.shape[0]/3))
        noise_norm = noise_norm[:, np.newaxis]


    # -------------------------------------------
    # return results
    # -------------------------------------------
    return kernel, noise_norm, vertno

def test_surface_vector_source_morph():
    """Test surface and vector source estimate morph."""
    tempdir = _TempDir()

    inverse_operator_surf = read_inverse_operator(fname_inv_surf)

    stc_surf = read_source_estimate(fname_smorph, subject='sample')
    stc_surf.crop(0.09, 0.1)  # for faster computation

    stc_vec = _real_vec_stc()

    source_morph_surf = compute_source_morph(
        inverse_operator_surf['src'], subjects_dir=subjects_dir,
        smooth=1, warn=False)  # smooth 1 for speed
    assert source_morph_surf.subject_from == 'sample'
    assert source_morph_surf.subject_to == 'fsaverage'
    assert source_morph_surf.kind == 'surface'
    assert isinstance(source_morph_surf.src_data, dict)
    assert isinstance(source_morph_surf.src_data['vertices_from'], list)
    assert isinstance(source_morph_surf, SourceMorph)
    stc_surf_morphed = source_morph_surf.apply(stc_surf)
    assert isinstance(stc_surf_morphed, SourceEstimate)
    stc_vec_morphed = source_morph_surf.apply(stc_vec)
    with pytest.raises(ValueError, match='Only volume source estimates'):
        source_morph_surf.apply(stc_surf, output='nifti1')

    # check if correct class after morphing
    assert isinstance(stc_surf_morphed, SourceEstimate)
    assert isinstance(stc_vec_morphed, VectorSourceEstimate)

    # check __repr__
    assert 'surface' in repr(source_morph_surf)

    # check loading and saving for surf
    source_morph_surf.save(op.join(tempdir, '42.h5'))

    source_morph_surf_r = read_source_morph(op.join(tempdir, '42.h5'))

    assert (all([read == saved for read, saved in
                 zip(sorted(source_morph_surf_r.__dict__),
                     sorted(source_morph_surf.__dict__))]))

    # check wrong subject correction
    stc_surf.subject = None
    assert isinstance(source_morph_surf.apply(stc_surf), SourceEstimate)

    # degenerate
    stc_vol = read_source_estimate(fname_vol, 'sample')
    with pytest.raises(ValueError, match='stc_from was type'):
        source_morph_surf.apply(stc_vol)

def test_snr():
    """Test SNR calculation"""
    tempdir = _TempDir()
    inv = read_inverse_operator(fname_inv)
    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0]
    snr = estimate_snr(evoked, inv)[0]
    orig_dir = os.getcwd()
    os.chdir(tempdir)
    try:
        cmd = ['mne_compute_mne', '--inv', fname_inv, '--meas', fname_evoked,
               '--snronly', '--bmin', '-200', '--bmax', '0']
        run_subprocess(cmd)
    except Exception:
        pass  # this returns 1 for some reason
    finally:
        os.chdir(orig_dir)
    snr_c = np.loadtxt(op.join(tempdir, 'SNR'))[:, 1]
    assert_allclose(snr, snr_c, atol=1e-2, rtol=1e-2)

def apply_STC_ave(fnevo, method='dSPM', snr=3.0, min_subject='fsaverage'):
    ''' Inverse evoked data into the source space. 
        Parameter
        ---------
        fnevo: string or list
            The evoked file with ECG, EOG and environmental noise free.
        method:string
            Inverse method, 'MNE' or 'dSPM'
        snr: float
            Signal to noise ratio for inverse solution.
        event: string
            The event name related with evoked data.
        baseline: bool
            If true, prestimulus segment from 'btmin' to 'btmax' will be saved, 
            If false, no baseline segment is saved.
        btmin: float
            The start time point (second) of baseline.
        btmax: float
            The end time point(second) of baseline.
        min_subject: string
            The subject name as the common brain.
    '''
    #Get the default subjects_dir
    from mne.minimum_norm import apply_inverse, read_inverse_operator
    fnlist = get_files_from_list(fnevo)
    # loop across all filenames
    for fname in fnlist:
        name = os.path.basename(fname)
        fn_path = os.path.split(fname)[0]
        fn_stc = fname[:fname.rfind('-ave.fif')] 
        #fn_inv = fname[:fname.rfind('-ave.fif')] + ',ave-inv.fif' 
        subject = name.split('_')[0]
        fn_inv = fn_path + '/%s_fibp1-45,ave-inv.fif' %subject
        snr = snr
        lambda2 = 1.0 / snr ** 2 
        #noise_cov = mne.read_cov(fn_cov)
        [evoked] = mne.read_evokeds(fname)
        evoked.pick_types(meg=True, ref_meg=False)
        inv = read_inverse_operator(fn_inv)
        stc = apply_inverse(evoked, inv, lambda2, method,
                            pick_ori='normal')
        stc.save(fn_stc)