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 single_epoch_tf_source(epochs, inv, src, label):
    """Calculates single trail power

    Parameter
    ---------
    epochs : ???
        The subject number to use.
    inv = inverse_operator
        ...
    src : source space
        ...
    label : label
        ...
    """
    snr = 1.0
    lambda2 = 1.0 / snr ** 2
    method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

    frequencies = np.arange(8, 13, 1)
    stcs = apply_inverse_epochs(epochs, inv, lambda2=lambda2, method=method,
                                label=None, pick_ori=None)
    time_series = [stc.extract_label_time_course(labels=label, src=src,
                                                 mode="pca_flip")[0]
                   for stc in stcs]

    ts_signed = []
    for j in range(len(time_series)):
        tmp = time_series[j]
        tmp *= np.sign(tmp[np.argmax(np.abs(tmp))])
        ts_signed.append(tmp)

    tfr = cwt_morlet(np.asarray(ts_signed), epochs.info["sfreq"], frequencies,
                     use_fft=True, n_cycles=4)

    return tfr

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 _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 fiff_mne(ds, fwd='{fif}*fwd.fif', cov='{fif}*cov.fif', label=None, name=None,
             tstart= -0.1, tstop=0.6, baseline=(None, 0)):
    """
    adds data from one label as

    """
    if name is None:
        if label:
            _, lbl = os.path.split(label)
            lbl, _ = os.path.splitext(lbl)
            name = lbl.replace('-', '_')
        else:
            name = 'stc'

    info = ds.info['info']

    raw = ds.info['raw']
    fif_name = raw.info['filename']
    fif_name, _ = os.path.splitext(fif_name)
    if fif_name.endswith('raw'):
        fif_name = fif_name[:-3]

    fwd = fwd.format(fif=fif_name)
    if '*' in fwd:
        d, n = os.path.split(fwd)
        names = fnmatch.filter(os.listdir(d), n)
        if len(names) == 1:
            fwd = os.path.join(d, names[0])
        else:
            raise IOError("No unique fwd file matching %r" % fwd)

    cov = cov.format(fif=fif_name)
    if '*' in cov:
        d, n = os.path.split(cov)
        names = fnmatch.filter(os.listdir(d), n)
        if len(names) == 1:
            cov = os.path.join(d, names[0])
        else:
            raise IOError("No unique cov file matching %r" % cov)

    fwd = mne.read_forward_solution(fwd, force_fixed=False, surf_ori=True)
    cov = mne.Covariance(cov)
    inv = _mn.make_inverse_operator(info, fwd, cov, loose=0.2, depth=0.8)
    epochs = mne_Epochs(ds, tstart=tstart, tstop=tstop, baseline=baseline)

    # mne example:
    snr = 3.0
    lambda2 = 1.0 / snr ** 2

    if label is not None:
        label = mne.read_label(label)
    stcs = _mn.apply_inverse_epochs(epochs, inv, lambda2, dSPM=False, label=label)

    x = np.vstack(s.data.mean(0) for s in stcs)
    s = stcs[0]
    dims = ('case', var(s.times, 'time'),)
    ds[name] = ndvar(x, dims, properties=None, info='')

    return stcs

                       exclude='bads')

# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=(None, 0), preload=True,
                    reject=dict(grad=4000e-13, eog=150e-6),
                    decim=5)  # decimate to save memory and increase speed

###############################################################################
# Compute inverse solution
snr = 3.0
noise_cov = mne.read_cov(fname_cov)
inverse_operator = read_inverse_operator(fname_inv)

stcs = apply_inverse_epochs(epochs, inverse_operator,
                            lambda2=1.0 / snr ** 2, verbose=False,
                            method="dSPM", pick_ori="normal")

###############################################################################
# Decoding in sensor space using a logistic regression

# Retrieve source space data into an array
X = np.array([stc.lh_data for stc in stcs])  # only keep left hemisphere
y = epochs.events[:, 2]

# prepare a series of classifier applied at each time sample
clf = make_pipeline(StandardScaler(),  # z-score normalization
                    SelectKBest(f_classif, k=500),  # select features for speed
                    LinearModel(LogisticRegression(C=1)))
time_decod = SlidingEstimator(clf, scoring='roc_auc')

events = mne.make_fixed_length_events(raw, 1, start=0, stop=None, 
		duration=4.)

# Set up pick list: (MEG minus bad channels)
include = []
exclude = raw.info['bads']
picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True, 
		include=include, exclude=exclude)

# Read epochs and remove bad epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, 
		picks=picks, baseline=(None, 0), preload=True, 
		reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))

# Compute the inverse solution
stc = apply_inverse_epochs(epochs, inverse_operator, lambda2, method)
#stc.save(data_path + subj + '/' + subj + '_rest_raw_sss-oct-6-inv.fif')

# define frequencies of interest
fmin, fmax = 0., 70.
bandwidth = 4.  # bandwidth of the windows in Hz

# compute source space psd in label

# Note: By using "return_generator=True" stcs will be a generator object
# instead of a list. This allows us so to iterate without having to
# keep everything in memory.

stcs = compute_source_psd_epochs(epochs, inverse_operator, lambda2=lambda2,
                                 method=method, fmin=fmin, fmax=fmax,
                                 bandwidth=bandwidth, label=label, return_generator=True)

# load labels
labels = mne.read_labels_from_annot('0001', parc='PALS_B12_Lobes',
                                    # regexp="Bro",
                                    subjects_dir=subjects_dir)

labels_occ = [labels[9], labels[10]]

# Load data
inverse_operator = read_inverse_operator(fname_inv)
epochs = mne.read_epochs(fname_epochs)

epochs.crop(0, 1.1)
epochs.resample(200)

stcs_ent_left = apply_inverse_epochs(epochs["ent_left"], inverse_operator,
                                     lambda2, method, pick_ori="normal")

stcs_ent_right = apply_inverse_epochs(epochs["ent_right"], inverse_operator,
                                     lambda2, method, pick_ori="normal")

stcs_ctl_left = apply_inverse_epochs(epochs["ctl_left"], inverse_operator,
                                     lambda2, method, pick_ori="normal")

#src_ctl_l = np.asarray([stc.data.reshape(-1) for stc in stcs_ctl_left])
#src_ent_l = np.asarray([stc.data.reshape(-1) for stc in stcs_ent_left])

#data_ctl_l = np.squeeze(np.asarray(
#    mne.extract_label_time_course(stcs_ctl_left,
#                                  labels_occ[1],
#                                  inverse_operator["src"],
#                                  mode="pca_flip")))

lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
n_times = len(epochs.times)
n_vertices = 3732
n_epochs = len(epochs.events)

# Load data and compute inverse solution and stcs for each epoch.

noise_cov = mne.read_cov(fname_cov)
inverse_operator = read_inverse_operator(fname_inv)
X = np.zeros([n_epochs, n_vertices, n_times])

# to save memory, we'll load and transform our epochs step by step.
for condition_count, ep in zip([0, n_epochs / 2], epochs_list):
    stcs = apply_inverse_epochs(ep, inverse_operator, lambda2,
                            method, pick_normal=True,  # this saves us memory
                            return_generator=True)
    for jj, stc in enumerate(stcs):
        X[condition_count + jj] = stc.lh_data

###############################################################################
# Decoding in sensor space using a linear SVM

# Make arrays X and y such that :
# X is 3d with X.shape[0] is the total number of epochs to classify
# y is filled with integers coding for the class to predict
# We must have X.shape[0] equal to y.shape[0]

# we know the first half belongs to the first class, the second one
y = np.repeat([0, 1], len(X) / 2)   # belongs to the second class
X = X.reshape(n_epochs, n_vertices * n_times)

                                                    eog=150e-6))

# Compute inverse solution and for each epoch
snr = 1.0           # use smaller SNR for raw data
inv_method = 'dSPM'  # sLORETA, MNE, dSPM
parc = 'aparc'      # the parcellation to use, e.g., 'aparc' 'aparc.a2009s'

lambda2 = 1.0 / snr ** 2

# Compute inverse operator
inverse_operator = make_inverse_operator(raw.info, fwd, noise_cov,
                                         loose=None, depth=None,
                                         fixed=False)


stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method,
                            pick_ori=None, return_generator=True)

# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi
labels_parc = mne.read_labels_from_annot(subject, parc=parc,
                                         subjects_dir=subjects_dir)

# Average the source estimates within each label of the cortical parcellation
# and each sub structures contained in the src space
# If mode = 'mean_flip' this option is used only for the cortical label
src = inverse_operator['src']
label_ts = mne.extract_label_time_course(stcs, labels_parc, src,
                                         mode='mean_flip',
                                         allow_empty=True,
                                         return_generator=False)

# We compute the connectivity in the alpha band and plot it using a circular

method = "dSPM"
snr = 1.
lambda2 = 1. / snr**2

labels = mne.read_labels_from_annot(
    subject=subject, parc="PALS_B12_Brodmann", regexp="Brodmann")

condition = "interupt"

inv = read_inverse_operator(mne_folder + "%s_%s-inv.fif" % (subject, condition
                                                            ))
epochs = mne.read_epochs(epochs_folder + "%s_%s-epo.fif" % (subject, condition
                                                            ))
# epochs.resample(500)

stcs = apply_inverse_epochs(
    epochs["press"], inv, lambda2, method=method, pick_ori=None)
ts = [
    mne.extract_label_time_course(
        stc, labels, inv["src"], mode="mean_flip") for stc in stcs
]

# for h, tc in enumerate(ts):
#     for j, t in enumerate(tc):
#         t *= np.sign(t[np.argmax(np.abs(t))])
#         tc[j, :] = t
#     ts[h] = tc

ts = np.asarray(ts)
stc.save(source_folder + "%s_%s_epo" % (subject, condition))
np.save(source_folder + "ave_ts/%s_%s_ts-epo.npy" % (subject, condition), ts)

#.........这里部分代码省略.........
        Add sensor space data as NDVar as ``ds['meg']`` (default ``False``).
    src : False | 'ico' | 'vol'
        Add source space data as NDVar as ``ds['src']`` (default ``False``).
    sub : str | list | None
        Expresion for subset of events to load. For a very small dataset use e.g.
        ``[0,1]``.
    fixed : bool
        MNE inverse parameter.
    snr : scalar
        MNE inverse parameter.
    method : str
        MNE inverse parameter.
    rm : bool
        Pretend to be a repeated measures dataset (adds 'subject' variable).
    stc : bool
        Add mne SourceEstimate for source space data as ``ds['stc']`` (default
        ``False``).

    Returns
    -------
    ds : Dataset
        Dataset with epochs from the MNE sample dataset in ``ds['epochs']``.
    """
    data_dir = mne.datasets.sample.data_path()
    meg_dir = os.path.join(data_dir, "MEG", "sample")
    raw_file = os.path.join(meg_dir, "sample_audvis_filt-0-40_raw.fif")
    event_file = os.path.join(meg_dir, "sample_audvis_filt-0-40-eve.fif")
    subjects_dir = os.path.join(data_dir, "subjects")
    subject = "sample"
    label_path = os.path.join(subjects_dir, subject, "label", "%s.label")

    if not os.path.exists(event_file):
        raw = mne.io.Raw(raw_file)
        events = mne.find_events(raw, stim_channel="STI 014")
        mne.write_events(event_file, events)
    ds = load.fiff.events(raw_file, events=event_file)
    ds.index()
    ds.info["subjects_dir"] = subjects_dir
    ds.info["subject"] = subject
    ds.info["label"] = label_path

    # get the trigger variable form the dataset for eaier access
    trigger = ds["trigger"]

    # use trigger to add various labels to the dataset
    ds["condition"] = Factor(trigger, labels={1: "LA", 2: "RA", 3: "LV", 4: "RV", 5: "smiley", 32: "button"})
    ds["side"] = Factor(trigger, labels={1: "L", 2: "R", 3: "L", 4: "R", 5: "None", 32: "None"})
    ds["modality"] = Factor(trigger, labels={1: "A", 2: "A", 3: "V", 4: "V", 5: "None", 32: "None"})

    if rm:
        ds = ds.sub("trigger < 5")
        ds = ds.equalize_counts("side % modality")
        subject_f = ds.eval("side % modality").enumerate_cells()
        ds["subject"] = subject_f.as_factor("s%r", random=True)

    if sub:
        ds = ds.sub(sub)

    load.fiff.add_mne_epochs(ds, tmin, tmax, baseline)
    if sns:
        ds["meg"] = load.fiff.epochs_ndvar(ds["epochs"], data="mag", sysname="neuromag306mag")

    if not src:
        return ds
    elif src == "ico":
        src_tag = "ico-4"
    elif src == "vol":
        src_tag = "vol-10"
    else:
        raise ValueError("src = %r" % src)
    epochs = ds["epochs"]

    # get inverse operator
    inv_file = os.path.join(meg_dir, "sample_eelbrain_%s-inv.fif" % src_tag)
    if os.path.exists(inv_file):
        inv = mne.minimum_norm.read_inverse_operator(inv_file)
    else:
        fwd_file = os.path.join(meg_dir, "sample-%s-fwd.fif" % src_tag)
        bem_dir = os.path.join(subjects_dir, subject, "bem")
        bem_file = os.path.join(bem_dir, "sample-5120-5120-5120-bem-sol.fif")
        trans_file = os.path.join(meg_dir, "sample_audvis_raw-trans.fif")

        if os.path.exists(fwd_file):
            fwd = mne.read_forward_solution(fwd_file)
        else:
            src_ = _mne_source_space(subject, src_tag, subjects_dir)
            fwd = mne.make_forward_solution(epochs.info, trans_file, src_, bem_file, fwd_file)

        cov_file = os.path.join(meg_dir, "sample_audvis-cov.fif")
        cov = mne.read_cov(cov_file)
        inv = mn.make_inverse_operator(epochs.info, fwd, cov, None, None, fixed)
        mne.minimum_norm.write_inverse_operator(inv_file, inv)
    ds.info["inv"] = inv

    stcs = mn.apply_inverse_epochs(epochs, inv, 1.0 / (snr ** 2), method)
    ds["src"] = load.fiff.stc_ndvar(stcs, subject, src_tag, subjects_dir, method, fixed)
    if stc:
        ds["stc"] = stcs

    return ds

# Set up pick list: (MEG minus bad channels)
include = []
exclude = raw.info['bads']
picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True, 
		include=include, exclude=exclude)

# Read epochs and remove bad epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, 
		picks=picks, baseline=(None, 0), preload=True, 
		reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))

# Pull data for averaging later
epc_array = epochs.get_data()

# Compute the inverse solution
inv = apply_inverse_epochs(epochs, inverse_operator, lambda2, method, label=label)

#Need to add a line here to automatically create stc directory within subj
#

epoch_num = 1
epoch_num_str = str(epoch_num)
for i in inv:
	i.save(data_path + subj + '/tmp/' + label_name[3:] + '_rest_raw_sss-oct-6-inv' + epoch_num_str)
	epoch_num = epoch_num + 1
	epoch_num_str = str(epoch_num)

# The following is used to remove the empty opposing hemisphere files
# and then move the files to save into the appropriate directory

if hemi == 'left':

method = "MNE"
n_jobs = 1

# Load data
inverse_nrm = read_inverse_operator(inverse_fnrm)
inverse_hyp = read_inverse_operator(inverse_fhyp)

epochs_nrm = mne.read_epochs(epochs_fnrm)
epochs_hyp = mne.read_epochs(epochs_fhyp)

epochs_nrm = epochs_nrm["Tone"]
epochs_hyp = epochs_hyp["Tone"]


#
stcs_nrm = apply_inverse_epochs(epochs_nrm, inverse_nrm, lambda2, method, pick_ori="normal", return_generator=False)
stcs_hyp = apply_inverse_epochs(epochs_hyp, inverse_hyp, lambda2, method, pick_ori="normal", return_generator=False)


# resample
[stc.resample(300) for stc in stcs_nrm]
[stc.resample(300) for stc in stcs_hyp]

# Get labels from FreeSurfer cortical parcellation
labels = mne.read_labels_from_annot("subject_1", parc="PALS_B12_Brodmann", regexp="Brodmann", subjects_dir=subjects_dir)

# Average the source estimates within eachh label using sign-flips to reduce
# signal cancellations, also here we return a generator
src_nrm = inverse_nrm["src"]
label_ts_nrm = mne.extract_label_time_course(stcs_nrm, labels, src_nrm, mode="mean_flip", return_generator=False)

def compute_ROIs_inv_sol(raw_filename, sbj_id, sbj_dir, fwd_filename,
                         cov_fname, is_epoched=False, event_id=None,
                         t_min=None, t_max=None,
                         is_evoked=False, events_id=[],
                         snr=1.0, inv_method='MNE',
                         parc='aparc', aseg=False, aseg_labels=[],
                         is_blind=False, labels_removed=[], save_stc=False):
    import os
    import os.path as op
    import numpy as np
    import mne
    import pickle

    from mne.io import read_raw_fif
    from mne import read_epochs
    from mne.minimum_norm import make_inverse_operator, apply_inverse_raw
    from mne.minimum_norm import apply_inverse_epochs, apply_inverse
    from mne import get_volume_labels_from_src

    from nipype.utils.filemanip import split_filename as split_f

    from neuropype_ephy.preproc import create_reject_dict

    try:
        traits.undefined(event_id)
    except NameError:
        event_id = None

    print '\n*** READ raw filename %s ***\n' % raw_filename
    if is_epoched and event_id is None:
        epochs = read_epochs(raw_filename)
        info = epochs.info
    else:
        raw = read_raw_fif(raw_filename)
        info = raw.info

    subj_path, basename, ext = split_f(info['filename'])

    print '\n*** READ noise covariance %s ***\n' % cov_fname
    noise_cov = mne.read_cov(cov_fname)

    print '\n*** READ FWD SOL %s ***\n' % fwd_filename
    forward = mne.read_forward_solution(fwd_filename)

    if not aseg:
        forward = mne.convert_forward_solution(forward, surf_ori=True,
                                               force_fixed=False)

    lambda2 = 1.0 / snr ** 2

    # compute inverse operator
    print '\n*** COMPUTE INV OP ***\n'
    if not aseg:
        loose = 0.2
        depth = 0.8
    else:
        loose = None
        depth = None

    inverse_operator = make_inverse_operator(info, forward, noise_cov,
                                             loose=loose, depth=depth,
                                             fixed=False)

    # apply inverse operator to the time windows [t_start, t_stop]s
    print '\n*** APPLY INV OP ***\n'
    if is_epoched and event_id is not None:
        events = mne.find_events(raw)
        picks = mne.pick_types(info, meg=True, eog=True, exclude='bads')
        reject = create_reject_dict(info)

        if is_evoked:
            epochs = mne.Epochs(raw, events, events_id, t_min, t_max,
                                picks=picks, baseline=(None, 0), reject=reject)
            evoked = [epochs[k].average() for k in events_id]
            snr = 3.0
            lambda2 = 1.0 / snr ** 2

            ev_list = events_id.items()
            for k in range(len(events_id)):
                stc = apply_inverse(evoked[k], inverse_operator, lambda2,
                                    inv_method, pick_ori=None)

                print '\n*** STC for event %s ***\n' % ev_list[k][0]
                stc_file = op.abspath(basename + '_' + ev_list[k][0])

                print '***'
                print 'stc dim ' + str(stc.shape)
                print '***'

                if not aseg:
                    stc.save(stc_file)

        else:
            epochs = mne.Epochs(raw, events, event_id, t_min, t_max,
                                picks=picks, baseline=(None, 0), reject=reject)
            stc = apply_inverse_epochs(epochs, inverse_operator, lambda2,
                                       inv_method, pick_ori=None)

            print '***'
            print 'len stc %d' % len(stc)
#.........这里部分代码省略.........

events = mne.read_events(fname_event)

# Set up pick list
include = []
exclude = raw.info['bads'] + ['EEG 053']  # bads + 1 more

# pick MEG channels
picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                                            include=include, exclude=exclude)
# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
                                                    eog=150e-6))

# Compute inverse solution and stcs for each epoch
stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method, label,
                            pick_normal=True)

data = sum(stc.data for stc in stcs) / len(stcs)

# compute sign flip to avoid signal cancelation when averaging signed values
flip = mne.label_sign_flip(label, inverse_operator['src'])

label_mean = np.mean(data, axis=0)
label_mean_flip = np.mean(flip[:, np.newaxis] * data, axis=0)

###############################################################################
# View activation time-series
pl.figure()
h0 = pl.plot(1e3 * stcs[0].times, data.T, 'k')
h1 = pl.plot(1e3 * stcs[0].times, label_mean, 'r', linewidth=3)
h2 = pl.plot(1e3 * stcs[0].times, label_mean_flip, 'g', linewidth=3)

# pick MEG channels
picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                   exclude='bads')

# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
                                                    eog=150e-6))

# Compute inverse solution and for each epoch. Note that since we are passing
# the output to both extract_label_time_course and the phase_slope_index
# functions, we have to use "return_generator=False", since it is only possible
# to iterate over generators once.
snr = 1.0  # use lower SNR for single epochs
lambda2 = 1.0 / snr ** 2
stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method,
                            pick_normal=True, return_generator=False)

# Now, we generate seed time series by averaging the activity in the left
# visual corex
label = mne.read_label(fname_label)
src = inverse_operator['src']  # the source space used
seed_ts = mne.extract_label_time_course(stcs, label, src, mode='mean_flip')

# Combine the seed time course with the source estimates. There will be a total
# of 7500 signals:
# index 0: time course extracted from label
# index 1..7499: dSPM source space time courses
comb_ts = zip(seed_ts, stcs)

# Construct indices to estimate connectivity between the label time course
# and all source space time courses

def get_mne_sample(tmin=-0.1, tmax=0.4, baseline=(None, 0), sns=False,
                   src=None, sub="modality=='A'", fixed=False, snr=2,
                   method='dSPM'):
    """Load events and epochs from the MNE sample data

    Parameters
    ----------
    tmin, tmax baseline :
        Epoch parameters.
    sns : bool
        Add sensor space data as NDVar as ``ds['sns']``.
    src : None | 'ico' | 'vol'
        Add source space data as NDVar as ``ds['src']``.
    sub : str | None
        Expresion for subset of events to load.
    fixed : bool
        MNE inverse parameter.
    snr : scalar
        MNE inverse parameter.
    method : str
        MNE inverse parameter.

    Returns
    -------
    ds : Dataset
        Dataset with epochs from the MNE sample dataset.
    """
    data_dir = mne.datasets.sample.data_path()
    meg_dir = os.path.join(data_dir, 'MEG', 'sample')
    raw_file = os.path.join(meg_dir, 'sample_audvis_filt-0-40_raw.fif')
    subjects_dir = os.path.join(data_dir, 'subjects')
    subject = 'sample'
    label_path = os.path.join(subjects_dir, subject, 'label', '%s.label')

    ds = load.fiff.events(raw_file, stim_channel='STI 014')
    ds.info['subjects_dir'] = subjects_dir
    ds.info['subject'] = subject
    ds.info['label'] = label_path

    # get the trigger variable form the dataset for eaier access
    trigger = ds['trigger']

    # use trigger to add various labels to the dataset
    ds['condition'] = Factor(trigger, labels={1:'LA', 2:'RA', 3:'LV', 4:'RV',
                                              5:'smiley', 32:'button'})
    ds['side'] = Factor(trigger, labels={1: 'L', 2:'R', 3:'L', 4:'R',
                                         5:'None', 32:'None'})
    ds['modality'] = Factor(trigger, labels={1: 'A', 2:'A', 3:'V', 4:'V',
                                             5:'None', 32:'None'})

    if sub:
        ds = ds.sub(sub)

    load.fiff.add_mne_epochs(ds, tmin, tmax, baseline)
    if sns:
        ds['sns'] = load.fiff.epochs_ndvar(ds['epochs'])

    if not src:
        return ds

    bem_dir = os.path.join(subjects_dir, subject, 'bem')
    bem_file = os.path.join(bem_dir, 'sample-5120-5120-5120-bem-sol.fif')
    trans_file = os.path.join(meg_dir, 'sample_audvis_raw-trans.fif')
    epochs = ds['epochs']
    if src == 'ico':
        src_tag = 'ico-4'
    elif src == 'vol':
        src_tag = 'vol-10'
    else:
        raise ValueError("src = %r" % src)

    fwd_file = os.path.join(meg_dir, 'sample-%s-fwd.fif' % src_tag)
    if os.path.exists(fwd_file):
        fwd = mne.read_forward_solution(fwd_file)
    else:
        src_file = os.path.join(bem_dir, 'sample-%s-src.fif' % src_tag)
        if os.path.exists(src_file):
            src_ = src_file
        elif src == 'ico':
            src_ = mne.setup_source_space(subject, src_file, 'ico4',
                                          subjects_dir=subjects_dir)
        elif src == 'vol':
            mri_file = os.path.join(subjects_dir, subject, 'mri', 'orig.mgz')
            src_ = mne.setup_volume_source_space(subject, src_file, pos=10.,
                                                 mri=mri_file, bem=bem_file,
                                                 mindist=0., exclude=0.,
                                                 subjects_dir=subjects_dir)
        fwd = mne.make_forward_solution(epochs.info, trans_file, src_,
                                        bem_file, fwd_file)

    cov_file = os.path.join(meg_dir, 'sample_audvis-cov.fif')
    cov = mne.read_cov(cov_file)
    inv = mn.make_inverse_operator(epochs.info, fwd, cov, None, None,
                                   fixed)
    ds.info['inv'] = inv

    stcs = mn.apply_inverse_epochs(epochs, inv, 1. / (snr ** 2), method)
    ds['src'] = load.fiff.stc_ndvar(stcs, subject, src_tag, subjects_dir)

    return ds

#.........这里部分代码省略.........
        ``[0,1]``.
    fixed : bool
        MNE inverse parameter.
    snr : scalar
        MNE inverse parameter.
    method : str
        MNE inverse parameter.
    rm : bool
        Pretend to be a repeated measures dataset (adds 'subject' variable).
    stc : bool
        Add mne SourceEstimate for source space data as ``ds['stc']``.

    Returns
    -------
    ds : Dataset
        Dataset with epochs from the MNE sample dataset.
    """
    data_dir = mne.datasets.sample.data_path()
    meg_dir = os.path.join(data_dir, 'MEG', 'sample')
    raw_file = os.path.join(meg_dir, 'sample_audvis_filt-0-40_raw.fif')
    event_file = os.path.join(meg_dir, 'sample_audvis_filt-0-40-eve.fif')
    subjects_dir = os.path.join(data_dir, 'subjects')
    subject = 'sample'
    label_path = os.path.join(subjects_dir, subject, 'label', '%s.label')

    if not os.path.exists(event_file):
        raw = mne.io.Raw(raw_file)
        events = mne.find_events(raw, stim_channel='STI 014')
        mne.write_events(event_file, events)
    ds = load.fiff.events(raw_file, events=event_file)
    ds.index()
    ds.info['subjects_dir'] = subjects_dir
    ds.info['subject'] = subject
    ds.info['label'] = label_path

    # get the trigger variable form the dataset for eaier access
    trigger = ds['trigger']

    # use trigger to add various labels to the dataset
    ds['condition'] = Factor(trigger, labels={1:'LA', 2:'RA', 3:'LV', 4:'RV',
                                              5:'smiley', 32:'button'})
    ds['side'] = Factor(trigger, labels={1: 'L', 2:'R', 3:'L', 4:'R',
                                         5:'None', 32:'None'})
    ds['modality'] = Factor(trigger, labels={1: 'A', 2:'A', 3:'V', 4:'V',
                                             5:'None', 32:'None'})

    if rm:
        ds = ds.sub('trigger < 5')
        ds = ds.equalize_counts('side % modality')
        subject_f = ds.eval('side % modality').enumerate_cells()
        ds['subject'] = subject_f.as_factor('s%r', random=True)

    if sub:
        ds = ds.sub(sub)

    load.fiff.add_mne_epochs(ds, tmin, tmax, baseline)
    if sns:
        ds['sns'] = load.fiff.epochs_ndvar(ds['epochs'], data='mag')

    if src is None:
        return ds
    elif src == 'ico':
        src_tag = 'ico-4'
    elif src == 'vol':
        src_tag = 'vol-10'
    else:
        raise ValueError("src = %r" % src)
    epochs = ds['epochs']

    # get inverse operator
    inv_file = os.path.join(meg_dir, 'sample_eelbrain_%s-inv.fif' % src_tag)
    if os.path.exists(inv_file):
        inv = mne.minimum_norm.read_inverse_operator(inv_file)
    else:
        fwd_file = os.path.join(meg_dir, 'sample-%s-fwd.fif' % src_tag)
        bem_dir = os.path.join(subjects_dir, subject, 'bem')
        bem_file = os.path.join(bem_dir, 'sample-5120-5120-5120-bem-sol.fif')
        trans_file = os.path.join(meg_dir, 'sample_audvis_raw-trans.fif')

        if os.path.exists(fwd_file):
            fwd = mne.read_forward_solution(fwd_file)
        else:
            src_ = _mne_source_space(subject, src_tag, subjects_dir)
            fwd = mne.make_forward_solution(epochs.info, trans_file, src_,
                                            bem_file, fwd_file)

        cov_file = os.path.join(meg_dir, 'sample_audvis-cov.fif')
        cov = mne.read_cov(cov_file)
        inv = mn.make_inverse_operator(epochs.info, fwd, cov, None, None,
                                       fixed)
        mne.minimum_norm.write_inverse_operator(inv_file, inv)
    ds.info['inv'] = inv

    stcs = mn.apply_inverse_epochs(epochs, inv, 1. / (snr ** 2), method)
    ds['src'] = load.fiff.stc_ndvar(stcs, subject, src_tag, subjects_dir,
                                    method, fixed)
    if stc:
        ds['stc'] = stcs

    return ds

# Set up pick list
include = []

# Add a bad channel
raw.info['bads'] += ['EEG 053']  # bads + 1 more

# pick MEG channels
picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                   include=include, exclude='bads')
# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
                                                    eog=150e-6))

# Compute inverse solution and stcs for each epoch
stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method, label,
                            pick_ori="normal")

mean_stc = sum(stcs) /