def per2test(X1, X2, p_thr, p, tstep, n_per=8192, fn_clu_out=None):

    #    Note that X needs to be a multi-dimensional array of shape
    #    samples (subjects) x time x space, so we permute dimensions
    n_subjects1 = X1.shape[2]
    n_subjects2 = X2.shape[2]
    fsave_vertices = [np.arange(X1.shape[0]/2), np.arange(X1.shape[0]/2)]
    X1 = np.transpose(X1, [2, 1, 0])
    X2 = np.transpose(X2, [2, 1, 0])
    X = [X1, X2]
    #    Now let's actually do the clustering. This can take a long time...
    #    Here we set the threshold quite high to reduce computation.
    f_threshold = stats.distributions.f.ppf(1. - p_thr / 2., n_subjects1 - 1, n_subjects2 - 1)
    print('Clustering.')
    connectivity = spatial_tris_connectivity(grade_to_tris(5))
    T_obs, clusters, cluster_p_values, H0 = clu = \
        spatio_temporal_cluster_test(X, n_permutations=n_per, 
                                    connectivity=connectivity, n_jobs=1,
                                    threshold=f_threshold)
    #    Now select the clusters that are sig. at p < 0.05 (note that this value
    #    is multiple-comparisons corrected).
    good_cluster_inds = np.where(cluster_p_values < p)[0]
    print 'the amount of significant clusters are: %d' %good_cluster_inds.shape
    ###############################################################################
    # Save the clusters as stc file
    # ----------------------
    assert good_cluster_inds.shape != 0, ('Current p_threshold is %f %p_thr,\
                                 maybe you need to reset a lower p_threshold')
    np.savez(fn_clu_out, clu=clu, tstep=tstep, fsave_vertices=fsave_vertices)

def run_permutation_ttest(tmin=None, tmax=None, p_threshold = 0.05, n_permutations=1024, inverse_method='dSPM', n_jobs=1):
    for cond_id, cond_name in enumerate(events_id.keys()):
        #todo: calc the 36
        controls_data = get_morphed_epochs_stcs(tmin, tmax, cond_name, get_healthy_controls(),
            36, inverse_method)
        controls_data = abs(controls_data)
        for patient in get_patients():
            try:
                print(patient, cond_name)
                patient_data = get_morphed_epochs_stcs(tmin, tmax, cond_name, [patient], None, inverse_method)
                patient_data = abs(patient_data)
                print(patient_data.shape, controls_data.shape)
                data = controls_data - patient_data
                del patient_data
                gc.collect()
                data = np.transpose(data, [2, 1, 0])
                connectivity = spatial_tris_connectivity(grade_to_tris(5))
                t_threshold = -stats.distributions.t.ppf(p_threshold / 2., data.shape[0] - 1)
                T_obs, clusters, cluster_p_values, H0 = \
                    spatio_temporal_cluster_1samp_test(data, connectivity=connectivity, n_jobs=n_jobs,
                        threshold=t_threshold, n_permutations=n_permutations)
                results_file_name = op.join(LOCAL_ROOT_DIR, 'permutation_ttest_results', '{}_{}_{}'.format(patient, cond_name, inverse_method))
                np.savez(results_file_name, T_obs=T_obs, clusters=clusters, cluster_p_values=cluster_p_values, H0=H0)
                good_cluster_inds = np.where(cluster_p_values < 0.05)[0]
                print('good_cluster_inds: {}'.format(good_cluster_inds))
            except:
                print('bummer! {}, {}'.format(patient, cond_name))
                print(traceback.format_exc())

def stat_clus(X, tstep, time_thre=0, n_per=8192, p_threshold=0.01, p=0.05, fn_clu_out=None):
    print('Computing connectivity.')
    connectivity = spatial_tris_connectivity(grade_to_tris(5))
    #    Note that X needs to be a multi-dimensional array of shape
    #    samples (subjects) x time x space, so we permute dimensions
    X = np.transpose(X, [2, 1, 0])
    n_subjects = X.shape[0]
    fsave_vertices = [np.arange(X.shape[-1]/2), np.arange(X.shape[-1]/2)]
    #    Now let's actually do the clustering. This can take a long time...
    #    Here we set the threshold quite high to reduce computation.
    t_threshold = -stats.distributions.t.ppf(p_threshold / 2., n_subjects - 1)
    print('Clustering.')
    max_step = int(time_thre * 0.001 / tstep) + 1
    T_obs, clusters, cluster_p_values, H0 = clu = \
        spatio_temporal_cluster_1samp_test(X, connectivity=connectivity, n_jobs=1, max_step=max_step,
                                        threshold=t_threshold, n_permutations=n_per)
    #    Now select the clusters that are sig. at p < 0.05 (note that this value
    #    is multiple-comparisons corrected).
    good_cluster_inds = np.where(cluster_p_values < p)[0]
    print 'the amount of significant clusters are: %d' %good_cluster_inds.shape
    ###############################################################################
    # Save the clusters as stc file
    # ----------------------
    assert good_cluster_inds.shape != 0, ('Current p_threshold is %f %p_thr,\
                                 maybe you need to reset a lower p_threshold')
    np.savez(fn_clu_out, clu=clu, tstep=tstep, fsave_vertices=fsave_vertices)

def test_stc_to_label():
    """Test stc_to_label
    """
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        src = read_source_spaces(fwd_fname)
    src_bad = read_source_spaces(src_bad_fname)
    stc = read_source_estimate(stc_fname, 'sample')
    os.environ['SUBJECTS_DIR'] = op.join(data_path, 'subjects')
    labels1 = _stc_to_label(stc, src='sample', smooth=3)
    labels2 = _stc_to_label(stc, src=src, smooth=3)
    assert_equal(len(labels1), len(labels2))
    for l1, l2 in zip(labels1, labels2):
        assert_labels_equal(l1, l2, decimal=4)

    with warnings.catch_warnings(record=True) as w:  # connectedness warning
        warnings.simplefilter('always')
        labels_lh, labels_rh = stc_to_label(stc, src=src, smooth=True,
                                            connected=True)

    assert_true(len(w) > 0)
    assert_raises(ValueError, stc_to_label, stc, 'sample', smooth=True,
                  connected=True)
    assert_raises(RuntimeError, stc_to_label, stc, smooth=True, src=src_bad,
                  connected=True)
    assert_equal(len(labels_lh), 1)
    assert_equal(len(labels_rh), 1)

    # test getting tris
    tris = labels_lh[0].get_tris(src[0]['use_tris'], vertices=stc.vertices[0])
    assert_raises(ValueError, spatial_tris_connectivity, tris,
                  remap_vertices=False)
    connectivity = spatial_tris_connectivity(tris, remap_vertices=True)
    assert_true(connectivity.shape[0] == len(stc.vertices[0]))

    # "src" as a subject name
    assert_raises(TypeError, stc_to_label, stc, src=1, smooth=False,
                  connected=False, subjects_dir=subjects_dir)
    assert_raises(ValueError, stc_to_label, stc, src=SourceSpaces([src[0]]),
                  smooth=False, connected=False, subjects_dir=subjects_dir)
    assert_raises(ValueError, stc_to_label, stc, src='sample', smooth=False,
                  connected=True, subjects_dir=subjects_dir)
    assert_raises(ValueError, stc_to_label, stc, src='sample', smooth=True,
                  connected=False, subjects_dir=subjects_dir)
    labels_lh, labels_rh = stc_to_label(stc, src='sample', smooth=False,
                                        connected=False,
                                        subjects_dir=subjects_dir)
    assert_true(len(labels_lh) > 1)
    assert_true(len(labels_rh) > 1)

    # with smooth='patch'
    with warnings.catch_warnings(record=True) as w:  # connectedness warning
        warnings.simplefilter('always')
        labels_patch = stc_to_label(stc, src=src, smooth=True)
    assert_equal(len(w), 1)
    assert_equal(len(labels_patch), len(labels1))
    for l1, l2 in zip(labels1, labels2):
        assert_labels_equal(l1, l2, decimal=4)

def per2test(X1, X2, p_thr, p, tstep, n_per=8192, fn_clu_out=None):
    '''
      Calculate significant clusters using 2 sample ftest.

      Parameter
      ---------
      X1, X2: array
        The shape of X should be (Vertices, timepoints, subjects)
      tstep: float
        The interval between timepoints.
      n_per: int
        The permutation for ttest.
      p_thr: float
        The significant p_values.
      p: float
        The corrected p_values for comparisons.
      fn_clu_out: string
        The fnname for saving clusters.
    '''
    #    Note that X needs to be a multi-dimensional array of shape
    #    samples (subjects) x time x space, so we permute dimensions
    n_subjects1 = X1.shape[2]
    n_subjects2 = X2.shape[2]
    fsave_vertices = [np.arange(X1.shape[0]/2), np.arange(X1.shape[0]/2)]
    X1 = np.transpose(X1, [2, 1, 0])
    X2 = np.transpose(X2, [2, 1, 0])
    X = [X1, X2]

    #    Now let's actually do the clustering. This can take a long time...
    #    Here we set the threshold quite high to reduce computation.
    f_threshold = stats.distributions.f.ppf(1. - p_thr / 2., n_subjects1 - 1,
                                            n_subjects2 - 1)
    # t_threshold = stats.distributions.t.ppf(1. - p_thr / 2., n_subjects1 - 1,
    #                                         n_subjects2 - 1)

    print('Clustering...')
    connectivity = spatial_tris_connectivity(grade_to_tris(5))
    T_obs, clusters, cluster_p_values, H0 = clu = \
        spatio_temporal_cluster_test(X, n_permutations=n_per, #step_down_p=0.001,
                                     connectivity=connectivity, n_jobs=1,
                                     # threshold=t_threshold, stat_fun=stats.ttest_ind)
                                     threshold=f_threshold)

    #    Now select the clusters that are sig. at p < 0.05 (note that this value
    #    is multiple-comparisons corrected).
    good_cluster_inds = np.where(cluster_p_values < p)[0]
    print 'the amount of significant clusters are: %d' % good_cluster_inds.shape

    # Save the clusters as stc file
    np.savez(fn_clu_out, clu=clu, tstep=tstep, fsave_vertices=fsave_vertices)
    assert good_cluster_inds.shape != 0, ('Current p_threshold is %f %p_thr,\
                                 maybe you need to reset a lower p_threshold')

def create_spatial_connectivity(subject):
    try:
        connectivity_per_hemi = {}
        for hemi in utils.HEMIS:
            d = np.load(op.join(SUBJECTS_DIR, subject, 'mmvt', '{}.pial.npz'.format(hemi)))
            connectivity_per_hemi[hemi] = mne.spatial_tris_connectivity(d['faces'])
        utils.save(connectivity_per_hemi, op.join(MMVT_DIR, subject, 'spatial_connectivity.pkl'))
        success = True
    except:
        print('Error in create_spatial_connectivity!')
        print(traceback.format_exc())
        success = False
    return success

def test_stc_to_label():
    """Test stc_to_label."""
    src = read_source_spaces(fwd_fname)
    src_bad = read_source_spaces(src_bad_fname)
    stc = read_source_estimate(stc_fname, 'sample')
    os.environ['SUBJECTS_DIR'] = op.join(data_path, 'subjects')
    labels1 = _stc_to_label(stc, src='sample', smooth=3)
    labels2 = _stc_to_label(stc, src=src, smooth=3)
    assert_equal(len(labels1), len(labels2))
    for l1, l2 in zip(labels1, labels2):
        assert_labels_equal(l1, l2, decimal=4)

    with pytest.warns(RuntimeWarning, match='have holes'):
        labels_lh, labels_rh = stc_to_label(stc, src=src, smooth=True,
                                            connected=True)

    pytest.raises(ValueError, stc_to_label, stc, 'sample', smooth=True,
                  connected=True)
    pytest.raises(RuntimeError, stc_to_label, stc, smooth=True, src=src_bad,
                  connected=True)
    assert_equal(len(labels_lh), 1)
    assert_equal(len(labels_rh), 1)

    # test getting tris
    tris = labels_lh[0].get_tris(src[0]['use_tris'], vertices=stc.vertices[0])
    pytest.raises(ValueError, spatial_tris_connectivity, tris,
                  remap_vertices=False)
    connectivity = spatial_tris_connectivity(tris, remap_vertices=True)
    assert (connectivity.shape[0] == len(stc.vertices[0]))

    # "src" as a subject name
    pytest.raises(TypeError, stc_to_label, stc, src=1, smooth=False,
                  connected=False, subjects_dir=subjects_dir)
    pytest.raises(ValueError, stc_to_label, stc, src=SourceSpaces([src[0]]),
                  smooth=False, connected=False, subjects_dir=subjects_dir)
    pytest.raises(ValueError, stc_to_label, stc, src='sample', smooth=False,
                  connected=True, subjects_dir=subjects_dir)
    pytest.raises(ValueError, stc_to_label, stc, src='sample', smooth=True,
                  connected=False, subjects_dir=subjects_dir)
    labels_lh, labels_rh = stc_to_label(stc, src='sample', smooth=False,
                                        connected=False,
                                        subjects_dir=subjects_dir)
    assert (len(labels_lh) > 1)
    assert (len(labels_rh) > 1)

    # with smooth='patch'
    with pytest.warns(RuntimeWarning, match='have holes'):
        labels_patch = stc_to_label(stc, src=src, smooth=True)
    assert len(labels_patch) == len(labels1)
    for l1, l2 in zip(labels1, labels2):
        assert_labels_equal(l1, l2, decimal=4)

def test_spatial_src_connectivity():
    """Test spatial connectivity functionality."""
    # oct
    src = read_source_spaces(fname_src)
    assert src[0]['dist'] is not None  # distance info
    with pytest.warns(RuntimeWarning, match='will have holes'):
        con = spatial_src_connectivity(src).toarray()
    con_dist = spatial_src_connectivity(src, dist=0.01).toarray()
    assert (con == con_dist).mean() > 0.75
    # ico
    src = read_source_spaces(fname_src_fs)
    con = spatial_src_connectivity(src).tocsr()
    con_tris = spatial_tris_connectivity(grade_to_tris(5)).tocsr()
    assert con.shape == con_tris.shape
    assert_array_equal(con.data, con_tris.data)
    assert_array_equal(con.indptr, con_tris.indptr)
    assert_array_equal(con.indices, con_tris.indices)
    # one hemi
    con_lh = spatial_src_connectivity(src[:1]).tocsr()
    con_lh_tris = spatial_tris_connectivity(grade_to_tris(5)).tocsr()
    con_lh_tris = con_lh_tris[:10242, :10242].tocsr()
    assert_array_equal(con_lh.data, con_lh_tris.data)
    assert_array_equal(con_lh.indptr, con_lh_tris.indptr)
    assert_array_equal(con_lh.indices, con_lh_tris.indices)

def permutation_test_on_source_data_with_spatio_temporal_clustering(controls_data, patient_data, patient, cond_name,
                tstep, n_permutations, inverse_method='dSPM', n_jobs=6):
    try:
        print('permutation_test: patient {}, cond {}'.format(patient, cond_name))
        connectivity = spatial_tris_connectivity(grade_to_tris(5))
        #    Note that X needs to be a list of multi-dimensional array of shape
        #    samples (subjects_k) x time x space, so we permute dimensions
        print(controls_data.shape, patient_data.shape)
        X = [controls_data, patient_data]

        p_threshold = 0.05
        f_threshold = stats.distributions.f.ppf(1. - p_threshold / 2.,
                                                controls_data.shape[0] - 1, 1)
        print('Clustering. thtreshold = {}'.format(f_threshold))
        T_obs, clusters, cluster_p_values, H0 = clu =\
            spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=n_jobs, threshold=10, n_permutations=n_permutations)

        results_file_name = op.join(LOCAL_ROOT_DIR, 'clusters_results', '{}_{}_{}'.format(patient, cond_name, inverse_method))
        np.savez(results_file_name, T_obs=T_obs, clusters=clusters, cluster_p_values=cluster_p_values, H0=H0)
        #    Now select the clusters that are sig. at p < 0.05 (note that this value
        #    is multiple-comparisons corrected).
        good_cluster_inds = np.where(cluster_p_values < 0.05)[0]

        ###############################################################################
        # Visualize the clusters

        print('Visualizing clusters.')

        #    Now let's build a convenient representation of each cluster, where each
        #    cluster becomes a "time point" in the SourceEstimate
        fsave_vertices = [np.arange(10242), np.arange(10242)]
        stc_all_cluster_vis = summarize_clusters_stc(clu, tstep=tstep,
                                                     vertices=fsave_vertices,
                                                     subject='fsaverage')
        stc_all_cluster_vis.save(op.join(LOCAL_ROOT_DIR, 'stc_clusters', '{}_{}_{}'.format(patient, cond_name, inverse_method)), ftype='h5')

        # #    Let's actually plot the first "time point" in the SourceEstimate, which
        # #    shows all the clusters, weighted by duration
        # # blue blobs are for condition A != condition B
        # brain = stc_all_cluster_vis.plot('fsaverage', 'inflated', 'both',
        #                                  subjects_dir=subjects_dir, clim='auto',
        #                                  time_label='Duration significant (ms)')
        # brain.set_data_time_index(0)
        # brain.show_view('lateral')
        # brain.save_image('clusters.png')
    except:
        print('bummer! {}, {}'.format(patient, cond_name))
        print(traceback.format_exc())

def stat_clus(X, tstep, n_per=8192, p_threshold=0.01, p=0.05, fn_clu_out=None):
    '''
      Calculate significant clusters using 1sample ttest.

      Parameter
      ---------
      X: array
        The shape of X should be (Vertices, timepoints, subjects)
      tstep: float
        The interval between timepoints.
      n_per: int
        The permutation for ttest.
      p_threshold: float
        The significant p_values.
      p: float
        The corrected p_values for comparisons.
      fn_clu_out: string
        The fnname for saving clusters.
    '''

    print('Computing connectivity.')
    connectivity = spatial_tris_connectivity(grade_to_tris(5))

    #    Note that X needs to be a multi-dimensional array of shape
    #    samples (subjects) x time x space, so we permute dimensions
    X = np.transpose(X, [2, 1, 0])
    n_subjects = X.shape[0]
    fsave_vertices = [np.arange(X.shape[-1]/2), np.arange(X.shape[-1]/2)]

    #    Now let's actually do the clustering. This can take a long time...
    #    Here we set the threshold quite high to reduce computation.
    t_threshold = -stats.distributions.t.ppf(p_threshold / 2., n_subjects - 1)
    print('Clustering.')
    T_obs, clusters, cluster_p_values, H0 = clu = \
        spatio_temporal_cluster_1samp_test(X, connectivity=connectivity,
                                           n_jobs=1, threshold=t_threshold,
                                           n_permutations=n_per)

    #    Now select the clusters that are sig. at p < 0.05 (note that this value
    #    is multiple-comparisons corrected).
    good_cluster_inds = np.where(cluster_p_values < p)[0]
    print 'the amount of significant clusters are: %d' %good_cluster_inds.shape

    # Save the clusters as stc file
    np.savez(fn_clu_out, clu=clu, tstep=tstep, fsave_vertices=fsave_vertices)
    assert good_cluster_inds.shape != 0, ('Current p_threshold is %f %p_thr,\
                                 maybe you need to reset a lower p_threshold')

def create_spatial_connectivity(subject):
    try:
        verts_neighbors_fname = op.join(MMVT_DIR, subject, 'verts_neighbors_{hemi}.pkl')
        connectivity_fname = op.join(MMVT_DIR, subject, 'spatial_connectivity.pkl')
        if utils.both_hemi_files_exist(verts_neighbors_fname) and op.isfile(verts_neighbors_fname):
            return True
        connectivity_per_hemi = {}
        for hemi in utils.HEMIS:
            neighbors = defaultdict(list)
            d = np.load(op.join(MMVT_DIR, subject, 'surf', '{}.pial.npz'.format(hemi)))
            connectivity_per_hemi[hemi] = mne.spatial_tris_connectivity(d['faces'])
            rows, cols = connectivity_per_hemi[hemi].nonzero()
            for ind in range(len(rows)):
                neighbors[rows[ind]].append(cols[ind])
            utils.save(neighbors, verts_neighbors_fname.format(hemi=hemi))
        utils.save(connectivity_per_hemi, connectivity_fname)
        success = True
    except:
        print('Error in create_spatial_connectivity!')
        print(traceback.format_exc())
        success = False
    return success

                        effects=effects, return_pvals=return_pvals)[0]
    # get f-values only.
    # Note. for further details on this ANOVA function consider the
    # corresponding time frequency example.

###############################################################################
# Compute clustering statistic

#    To use an algorithm optimized for spatio-temporal clustering, we
#    just pass the spatial connectivity matrix (instead of spatio-temporal)

source_space = grade_to_tris(5)
# as we only have one hemisphere we need only need half the connectivity
lh_source_space = source_space[source_space[:, 0] < 10242]
print('Computing connectivity.')
connectivity = spatial_tris_connectivity(lh_source_space)

#    Now let's actually do the clustering. Please relax, on a small
#    notebook and one single thread only this will take a couple of minutes ...
pthresh = 0.0005
f_thresh = f_threshold_twoway_rm(n_subjects, factor_levels, effects, pthresh)

#    To speed things up a bit we will ...
n_permutations = 128  # ... run fewer permutations (reduces sensitivity)

print('Clustering.')
T_obs, clusters, cluster_p_values, H0 = clu = \
    spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=1,
                                 threshold=f_thresh, stat_fun=stat_fun,
                                 n_permutations=n_permutations,
                                 buffer_size=None)

X2 = np.random.randn(n_vertices_fsave, n_times, n_subjects2) * 10
X1[:, :, :] += stc.data[:, :, np.newaxis]
# make the activity bigger for the second set of subjects
X2[:, :, :] += 3 * stc.data[:, :, np.newaxis]

#    We want to compare the overall activity levels for each subject
X1 = np.abs(X1)  # only magnitude
X2 = np.abs(X2)  # only magnitude

###############################################################################
# Compute statistic

#    To use an algorithm optimized for spatio-temporal clustering, we
#    just pass the spatial connectivity matrix (instead of spatio-temporal)
print('Computing connectivity.')
connectivity = spatial_tris_connectivity(grade_to_tris(5))

#    Note that X needs to be a list of multi-dimensional array of shape
#    samples (subjects_k) x time x space, so we permute dimensions
X1 = np.transpose(X1, [2, 1, 0])
X2 = np.transpose(X2, [2, 1, 0])
X = [X1, X2]

#    Now let's actually do the clustering. This can take a long time...
#    Here we set the threshold quite high to reduce computation.
p_threshold = 0.0001
f_threshold = stats.distributions.f.ppf(1. - p_threshold / 2.,
                                        n_subjects1 - 1, n_subjects2 - 1)
print('Clustering.')
T_obs, clusters, cluster_p_values, H0 = clu =\
    spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=2,

for fname in files:
    # exclude due to psychotropic factors
    if fname.find('WCEYBWMO') < 0:
        stc = mne.read_source_estimate(fname[:-7])
        if after_zero:
            stc.crop(tmin=0)
        X.append(stc.data)
        subj = fname.split('/')[-1].split('_')[0]
        gf_loc.append(np.nonzero(gf.maskid == subj)[0][0])
if len(X) != len(gf_loc):
    print 'Mismatch between subject data and gf!'

# re-sort subject order in gf to match X
gf = gf.iloc[gf_loc]

connectivity = mne.spatial_tris_connectivity(mne.grade_to_tris(5))

for thresh in p_threshold:
    if my_test == 'nvVSper':
        g0 = [X[i].T for i, group in enumerate(gf.group) if group == 'NV']
        g1 = [X[i].T for i, group in enumerate(gf.group) if group == 'persistent']
        data = [np.array(g0), np.array(g1)]
        stat_obs, clusters, p_values, H0 = \
            mne.stats.spatio_temporal_cluster_test(data, n_jobs=njobs,
                                                   threshold=thresh,
                                                   connectivity=connectivity,
                                                   stat_fun=group_comparison,
                                                   tail=1,
                                                   n_permutations=nperms,
                                                   verbose=True)
    elif my_test == 'nvVSrem':

def apply_sigSTC(fn_list_v, vevent, mevent, method='dSPM', vtmin=0., vtmax=0.35,
                 mtmin=-0.3, mtmax=0.05, radius=10.0):
    from mne import spatial_tris_connectivity, grade_to_tris
    from mne.stats import spatio_temporal_cluster_test
    from scipy import stats as stats
    X1, X2 = [], []
    stcs_trial = []
    for fn_v in fn_list_v:
        fn_m = fn_v[: fn_v.rfind('evtW')] + 'evtW_%s_bc_norm_1-lh.stc' %mevent
        stc_v = mne.read_source_estimate(fn_v)
        stcs_trial.append(stc_v.copy())
        stc_m = mne.read_source_estimate(fn_m)
        stc_v.resample(200)
        stc_m.resample(200)
        X1.append(stc_v.copy().crop(vtmin, vtmax).data)
        X2.append(stc_m.copy().crop(mtmin, mtmax).data)
    stcs_path = subjects_dir+'/fsaverage/%s_ROIs/conditions/' %method
    reset_directory(stcs_path)
    fn_avg = stcs_path + '%s' %(vevent)
    stcs = np.array(stcs_trial)
    stc_avg = np.sum(stcs, axis=0)/stcs.shape[0]
    stc_avg.save(fn_avg, ftype='stc')    
    X1 = np.array(X1).transpose(0, 2, 1)
    X2 = np.array(X2).transpose(0, 2, 1)
    ###############################################################################
    # Compute statistic
    
    #    To use an algorithm optimized for spatio-temporal clustering, we
    #    just pass the spatial connectivity matrix (instead of spatio-temporal)
    print('Computing connectivity.')
    connectivity = spatial_tris_connectivity(grade_to_tris(5))
    
    #    Note that X needs to be a list of multi-dimensional array of shape
    #    samples (subjects_k) x time x space, so we permute dimensions
    X = [X1, X2]
    #    Now let's actually do the clustering. This can take a long time...
    #    Here we set the threshold quite high to reduce computation.
    p_threshold = 0.0001
    f_threshold = stats.distributions.f.ppf(1. - p_threshold / 2.,
                                        X1.shape[0] - 1, X1.shape[0] - 1)
    print('Clustering.')
   
    clu = spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=2,
                                    threshold=f_threshold)
    #    Now select the clusters that are sig. at p < 0.05 (note that this value
    #    is multiple-comparisons corrected).
    #fsave_vertices = [np.arange(10242), np.arange(10242)]
    tstep = stc_v.tstep
    #stc_all_cluster_vis = summarize_clusters_stc(clu, tstep=tstep,
    #                                            vertices=fsave_vertices,
    #                                            subject='fsaverage')
    #stc_sig = stc_all_cluster_vis.mean()
    #fn_sig = subjects_dir+'/fsaverage/%s_ROIs/%s' %(method,vevent)
    #stc_sig.save(fn_sig)
    tstep = stc_v.tstep
    T_obs, clusters, clu_pvals, _ = clu
    n_times, n_vertices = T_obs.shape
    good_cluster_inds = np.where(clu_pvals < 0.05)[0]
    seeds = []
    #  Build a convenient representation of each cluster, where each
    #  cluster becomes a "time point" in the SourceEstimate
    T_obs = abs(T_obs)
    if len(good_cluster_inds) > 0:
        data = np.zeros((n_vertices, n_times))
        for cluster_ind in good_cluster_inds:
            data.fill(0)
            v_inds = clusters[cluster_ind][1]
            t_inds = clusters[cluster_ind][0]
            data[v_inds, t_inds] = T_obs[t_inds, v_inds]
            # Store a nice visualization of the cluster by summing across time
            data = np.sign(data) * np.logical_not(data == 0) * tstep
            seed = np.argmax(data.sum(axis=-1))
            seeds.append(seed)
    min_subject = 'fsaverage'
    labels_path = subjects_dir + '/fsaverage/dSPM_ROIs/%s/ini' %vevent
    reset_directory(labels_path)
    seeds = np.array(seeds)
    non_index_lh = seeds[seeds < 10242]
    if non_index_lh.shape != []:    
        func_labels_lh = mne.grow_labels(min_subject, non_index_lh,
                                        extents=radius, hemis=0, 
                                        subjects_dir=subjects_dir, n_jobs=1)
        i = 0
        while i < len(func_labels_lh):
            func_label = func_labels_lh[i]
            func_label.save(labels_path + '/%s_%d' %(vevent, i))
            i = i + 1
            
    seeds_rh = seeds - 10242
    non_index_rh = seeds_rh[seeds_rh > 0]
    if non_index_rh.shape != []:
        func_labels_rh = mne.grow_labels(min_subject, non_index_rh,
                                        extents=radius, hemis=1,
                                        subjects_dir=subjects_dir, n_jobs=1)                                             
   
        # right hemisphere definition
        j = 0
        while j < len(func_labels_rh):
            func_label = func_labels_rh[j]
            func_label.save(labels_path + '/%s_%d' %(vevent, j))
#.........这里部分代码省略.........

surf_connmat = normalize(surf_connmat, norm='l2')

# read surface seed roi mask
seedroi_gii = ng.read(seedroi_gii_path)
seedroi_data = seedroi_gii.darrays[darrays_int].data     # the values of seed vertex
print seedroi_data.shape
surfmask_inds = np.flatnonzero(seedroi_data)    # return the indices that are non-zeros in seedroi_data
print 'seed roi vertices ' + str(len(surfmask_inds))

# perform a hierarchical clustering considering spatial neighborhood (ward)
print 'ward: number of clusters:{}'.format(nb_clusters)
g = ng.read(mesh_path)
triangles = g.darrays[1].data

# compute the spatial neighbordhood over the seed surface
adjacency = sparse.coo_matrix(sparse.csc_matrix(sparse.csr_matrix(spatial_tris_connectivity(triangles))[surfmask_inds, :])
                              [:, surfmask_inds])
print " adjacency matrix: {}".format(adjacency.shape)

# ================================== ward clustering ====================================================
# ward = WardAgglomeration(n_clusters = nb_clusters, connectivity=connectivity, memory = 'nilearn_cache')
ward = AgglomerativeClustering(n_clusters=nb_clusters, linkage='ward', connectivity=adjacency)
ward.fit(surf_connmat)
labelsf = ward.labels_

for i in range(len(np.unique(labelsf))):
    print 'label %d: %d' % (i, len(labelsf[labelsf == i]))
# write the final gifti parcellation
print 'write parcellation.gii'

ii = 0