def test_morphing():
    mne.set_log_level('warning')
    sss = datasets._mne_source_space('fsaverage', 'ico-4', subjects_dir)
    vertices_to = [sss[0]['vertno'], sss[1]['vertno']]
    ds = datasets.get_mne_sample(-0.1, 0.1, src='ico', sub='index==0', stc=True)
    stc = ds['stc', 0]
    morph_mat = mne.compute_morph_matrix('sample', 'fsaverage', stc.vertno,
                                         vertices_to, None, subjects_dir)
    ndvar = ds['src']

    morphed_ndvar = morph_source_space(ndvar, 'fsaverage')
    morphed_stc = mne.morph_data_precomputed('sample', 'fsaverage', stc,
                                             vertices_to, morph_mat)
    assert_array_equal(morphed_ndvar.x[0], morphed_stc.data)
    morphed_stc_ndvar = load.fiff.stc_ndvar([morphed_stc], 'fsaverage', 'ico-4',
                                            subjects_dir, 'src', parc=None)
    assert_dataobj_equal(morphed_ndvar, morphed_stc_ndvar)

    def make_stcs(self, snr = 3.0, method='dSPM', save_to_disk = False, morph=True):
        """docstring for make_stcs"""

        morph_status = 'no_morph'

        lambda2 = 1.0 / snr ** 2.0
        pick_ori = None

        stcs = []

        for evoked_object, cond in self.evoked_list:

            print "Making source estimates for %s." %cond
            stc = mne.minimum_norm.apply_inverse(evoked_object, self.inverse_solution, method = method, lambda2 = lambda2, pick_ori = pick_ori)

            if morph == True:

                morph_status = 'morphed'
                # create morph map
                # get vertices to morph to (we'll take the fsaverage vertices)
                subject_to = 'fsaverage'
                fs = mne.read_source_spaces(self.subjects_dir + '%s/bem/%s-ico-4-src.fif' % (subject_to, subject_to))
                vertices_to = [fs[0]['vertno'], fs[1]['vertno']]

                # info of the stc we're going to morph is just the present stc
                subject_from = self.subject
                stc_from = stc

                # use the morph function
                morph_mat = mne.compute_morph_matrix(subject_from, subject_to, vertices_from=stc_from.vertices, vertices_to=vertices_to, subjects_dir=self.subjects_dir)
                stc = mne.morph_data_precomputed(subject_from, subject_to, stc_from, vertices_to, morph_mat)


                # stc_to.save('%s_audvis-meg' % subject_from)
                # mne.compute_morph_matrix('2-COMB','3-COMB',stcs[0].vertices, stcs[5].vertices, subjects_dir=self.subjects_dir)

            if save_to_disk:
                fact_morph_cond_path = op.join(self.stc_fact, morph_status, cond)
                if not os.path.isdir(fact_morph_cond_path):
                    os.makedirs(fact_morph_cond_path)

                stc.save(op.join(fact_morph_cond_path, '%s_%s_%s' %(self.subject, cond, morph_status)))

                self.add_preprocessing_notes("Saved source estimates (%s) for %s." %(morph_status, cond))

    def make_epoch_stcs(self, epochs, snr = 2.0, method='dSPM', morph=True, save_to_disk = False):
        """Apply inverse operator to epochs to get source estimates of each item"""


        lambda2 = 1.0 / snr ** 2.0

        inverse = mne.minimum_norm.read_inverse_operator( self.processed_files + self.subject + '_inv.fif' )

        eps = mne.minimum_norm.apply_inverse_epochs(epochs=epochs,inverse_operator=inverse,lambda2=lambda2,method = method)

        if morph == True:
            eps_morphed = []
            counter = 1
            morph_status = 'morphed'
            # create morph map
            # get vertices to morph to (we'll take the fsaverage vertices)
            subject_to = 'fsaverage'
            fs = mne.read_source_spaces(self.subjects_dir + '%s/bem/%s-ico-4-src.fif' % (subject_to, subject_to))
            vertices_to = [fs[0]['vertno'], fs[1]['vertno']]
            subject_from = self.subject

            for stc_from in eps:
                print "Morphing source estimate for epoch %d" %counter
            # use the morph function
                morph_mat = mne.compute_morph_matrix(subject_from, subject_to, vertices_from=stc_from.vertices, vertices_to=vertices_to, subjects_dir=self.subjects_dir)
                stc = mne.morph_data_precomputed(subject_from, subject_to, stc_from, vertices_to, morph_mat)
                # stc = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=1,
                #                     grade=vertices_to, subjects_dir=self.subjects_dir)

                eps_morphed.append(stc)
                counter += 1

            eps = eps_morphed

        if save_to_disk:
            with open(op.join(self.stc_cont, '%s_stc_epochs.pickled' %self.subject), 'w') as fileout:
                pickle.dump(eps, fileout)
            #save.pickle(obj=eps,dest=self.stc_cont + self.subject + '_stc_epochs')

        return eps

def test_source_estimate():
    "Test SourceSpace dimension"
    ds = datasets.get_mne_sample(src='ico')
    dsa = ds.aggregate('side')

    # test auto-conversion
    asndvar('epochs', ds=ds)
    asndvar('epochs', ds=dsa)
    asndvar(dsa['epochs'][0])

    # source space clustering
    res = testnd.ttest_ind('src', 'side', ds=ds, samples=0, pmin=0.05,
                           tstart=0.05, mintime=0.02, minsource=10)
    assert_equal(res.clusters.n_cases, 52)

    # test morphing
    dsa = ds.aggregate('side')
    ndvar = dsa['src']
    stc = mne.SourceEstimate(ndvar.x[0], ndvar.source.vertno,
                             ndvar.time.tmin, ndvar.time.tstep,
                             ndvar.source.subject)
    subjects_dir = ndvar.source.subjects_dir
    path = ndvar.source._src_pattern.format(subject='fsaverage',
                                            src=ndvar.source.src,
                                            subjects_dir=subjects_dir)
    if os.path.exists(path):
        src_to = mne.read_source_spaces(path)
    else:
        src_to = mne.setup_source_space('fsaverage', path, 'ico4',
                                        subjects_dir=subjects_dir)
    vertices_to = [src_to[0]['vertno'], src_to[1]['vertno']]
    mm = mne.compute_morph_matrix('sample', 'fsaverage', ndvar.source.vertno,
                                  vertices_to, None, subjects_dir)
    stc_to = mne.morph_data_precomputed('sample', 'fsaverage', stc,
                                        vertices_to, mm)

    ndvar_m = morph_source_space(ndvar, 'fsaverage')
    assert_array_equal(ndvar_m.x[0], stc_to.data)

def morph_source_space(ndvar, subject_to, vertices_to=None, morph_mat=None,
                       copy=False):
    """Morph source estimate to a different MRI subject

    Parameters
    ----------
    ndvar : NDVar
        NDVar with SourceSpace dimension.
    subject_to : string
        Name of the subject on which to morph.
    vertices_to : None | list of array of int
        The vertices on the destination subject's brain. If ndvar contains a
        whole source space, vertices_to can be automatically loaded, although
        providing them as argument can speed up processing by a second or two.
    morph_mat : None | sparse matrix
        The morphing matrix. If ndvar contains a whole source space, the morph
        matrix can be automatically loaded, although providing a cached matrix
        can speed up processing by a second or two.
    copy : bool
        Make sure that the data of ``morphed_ndvar`` is separate from
        ``ndvar`` (default False).

    Returns
    -------
    morphed_ndvar : NDVar
        NDVar morphed to the destination subject.

    Notes
    -----
    This function is used to make sure a number of different NDVars are defined
    on the same MRI subject and handles scaled MRIs efficiently. If the MRI
    subject on which ``ndvar`` is defined is a scaled copy of ``subject_to``,
    by default a shallow copy of ``ndvar`` is returned. That means that it is
    not safe to assume that ``morphed_ndvar`` can be modified in place without
    altering ``ndvar``. To make sure the date of the output is independent from
    the data of the input, set the argument ``copy=True``.
    """
    subjects_dir = ndvar.source.subjects_dir
    subject_from = ndvar.source.subject
    src = ndvar.source.src
    if vertices_to is None:
        path = SourceSpace._src_pattern.format(subjects_dir=subjects_dir,
                                               subject=subject_to, src=src)
        src_to = mne.read_source_spaces(path)
        vertices_to = [src_to[0]['vertno'], src_to[1]['vertno']]
    elif not isinstance(vertices_to, list) or not len(vertices_to) == 2:
        raise ValueError('vertices_to must be a list of length 2')

    if subject_from == subject_to and _vertices_equal(ndvar.source.vertno,
                                                      vertices_to):
        if copy:
            return ndvar.copy()
        else:
            return ndvar

    axis = ndvar.get_axis('source')
    x = ndvar.x

    # check whether it is a scaled brain
    do_morph = True
    cfg_path = os.path.join(subjects_dir, subject_from,
                            'MRI scaling parameters.cfg')
    if os.path.exists(cfg_path):
        cfg = mne.coreg.read_mri_cfg(subject_from, subjects_dir)
        subject_from = cfg['subject_from']
        if subject_to == subject_from and _vertices_equal(ndvar.source.vertno,
                                                          vertices_to):
            if copy:
                x_ = x.copy()
            else:
                x_ = x
            vertices_to = ndvar.source.vertno
            do_morph = False

    if do_morph:
        vertices_from = ndvar.source.vertno
        if morph_mat is None:
            morph_mat = mne.compute_morph_matrix(subject_from, subject_to,
                                                 vertices_from, vertices_to,
                                                 None, subjects_dir)
        elif not sp.sparse.issparse(morph_mat):
            raise ValueError('morph_mat must be a sparse matrix')
        elif not sum(len(v) for v in vertices_to) == morph_mat.shape[0]:
            raise ValueError('morph_mat.shape[0] must match number of '
                             'vertices in vertices_to')

        # flatten data
        if axis != 0:
            x = x.swapaxes(0, axis)
        n_sources = len(x)
        if not n_sources == morph_mat.shape[1]:
            raise ValueError('ndvar source dimension length must be the same '
                             'as morph_mat.shape[0]')
        if ndvar.ndim > 2:
            shape = x.shape
            x = x.reshape((n_sources, -1))

        # apply morph matrix
        x_ = morph_mat * x

#.........这里部分代码省略.........

vertno_max, time_max = stc.get_peak(hemi='rh')

subjects_dir = data_path + '/subjects'
brain = stc.plot(surface='inflated', hemi='rh', subjects_dir=subjects_dir,
                 clim=dict(kind='value', lims=[8, 12, 15]),
                 initial_time=time_max, time_unit='s')
brain.add_foci(vertno_max, coords_as_verts=True, hemi='rh', color='blue',
               scale_factor=0.6)
brain.show_view('lateral')

###############################################################################
# Morph data to average brain
# ---------------------------

fs_vertices = [np.arange(10242)] * 2
morph_mat = mne.compute_morph_matrix('sample', 'fsaverage', stc.vertices,
                                     fs_vertices, smooth=None,
                                     subjects_dir=subjects_dir)
stc_fsaverage = stc.morph_precomputed('fsaverage', fs_vertices, morph_mat)
brain_fsaverage = stc_fsaverage.plot(surface='inflated', hemi='rh',
                                     subjects_dir=subjects_dir,
                                     clim=dict(kind='value', lims=[8, 12, 15]),
                                     initial_time=time_max, time_unit='s')
brain_fsaverage.show_view('lateral')

###############################################################################
# Exercise
# --------
#    - By changing the method parameter to 'sloreta' recompute the source
#      estimates using the sLORETA method.

            
            n_epochs2[n][iCond] = evoked.nave
            
            stc = mne.minimum_norm.apply_inverse(evoked,inv,lambda2, method=method, pick_ori="normal")
            
            stc.crop(-0.1, 0.9)
            tmin = stc.tmin
            tstep = stc.tstep
            times = stc.times
            # Average brain
            """
            One should check if morph map is current and correct. Otherwise, it will spit out and error.
            Check SUBJECTS_DIR/morph-maps
            """
            morph_mat = mne.compute_morph_matrix(subs2[n], 'fsaverage', stc.vertices,
                                                 fs_vertices, smooth=None,
                                                 subjects_dir=fs_dir)
            stc_fsaverage = stc.morph_precomputed('fsaverage', fs_vertices, morph_mat)
        #    morph_dat = mne.morph_data(subs[n], 'fsaverage', stc, n_jobs=16,
        #                        grade=fs_vertices, subjects_dir=fs_dir)
            X2[:,:,n,iCond] = stc_fsaverage.data

    os.chdir(raw_dir)
    np.save(fname_data, X2)
    np.save('session2_times.npy',times)
    np.save('session2_tstep.npy',tstep)
    np.save('session2_n_averages.npy',n_epochs2)
else:
    os.chdir(raw_dir)
    X2 = np.load(fname_data)
    times = np.load('session2_times.npy')

def stft_source_localization(data, fn_inv, method="dSPM",
                             morph2fsaverage=False,
                             nave=1, snr=3.,
                             pick_ori="normal", verbose=True):

    """
    Apply inverse operator to data. In general, the
    L2-norm inverse solution is computed.

        Parameters
        ----------
        data: array of MEG data
            (only the data, no MNE-data object!)
        kernel: kernel of the inverse operator.
            (could be estimated using
            mne.minimum_norm._assemble_kernel())
        noise_norm: noise normalization array.
            (could be estimated using
            mne.minimum_norm._assemble_kernel())

        Returns
        -------
        src_loc: SourceEstimate | VolSourceEstimate
            The source estimates
        estimation_time: time needed to perform source
            localization (for one time slice)
    """

    # -------------------------------------------
    # import necessary modules
    # -------------------------------------------
    import numpy as np
    import types

    # check if data should be morphed
    if morph2fsaverage:
        from mne import compute_morph_matrix, grade_to_vertices, morph_data_precomputed
        from mne.source_estimate import SourceEstimate
        from os.path import basename, dirname


    # -------------------------------------------
    # estimate inverse kernel
    # -------------------------------------------
    kernel, noise_norm, vertno = calc_inv_kernel(fn_inv, method=method,
                                                 nave=nave, snr=snr,
                                                 pick_ori=pick_ori)


    # -------------------------------------------
    # get some information from the
    # input data
    # -------------------------------------------
    nfreq, nepochs, nchan = data.shape
    nvoxel = noise_norm.shape[0]

    if isinstance(data[0, 0, 0], types.ComplexType):
        src_loc_data = np.zeros((nfreq, nepochs, nvoxel), dtype=np.complex)
    else:
        src_loc_data = np.zeros((nfreq, nepochs, nvoxel))


    # -------------------------------------------
    # read in morphing matrix
    # -------------------------------------------
    if morph2fsaverage:
        subject_id = basename(fn_inv)[:6]
        subjects_dir = dirname(dirname(fn_inv))
        vertices_to = grade_to_vertices('fsaverage', grade=4,
                                        subjects_dir=subjects_dir)

        morph_mat = compute_morph_matrix(subject_id, 'fsaverage',
                                         vertno, vertices_to,
                                         subjects_dir=subjects_dir)
        nvoxel_morph = 2 * len(vertices_to[0])

        if isinstance(data[0, 0, 0], types.ComplexType):
            morphed_src_loc_data = np.zeros((nfreq, nepochs, nvoxel_morph), dtype=np.complex)
        else:
            morphed_src_loc_data = np.zeros((nfreq, nepochs, nvoxel_morph), dtype=np.complex)


    # -------------------------------------------
    # apply inverse operator for each time slice
    # -------------------------------------------
    for iepoch in range(nepochs):

        if verbose:
            from sys import stdout
            info = "\r" if iepoch > 0 else ""
            info += "... --> Epoch %d of %d done" % (iepoch+1, nepochs)
            stdout.write(info)
            stdout.flush()

        for ifreq in range(0, nfreq):
            # multiply measured data with inverse kernel
            loc_tmp = np.dot(kernel, data[ifreq, iepoch, :])

            if pick_ori != "normal":

#.........这里部分代码省略.........

src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'

# Read input stc file
stc_from = mne.read_source_estimate(fname)
# Morph using one method (supplying the vertices in fsaverage's source
# space makes it faster). Note that for any generic subject, you could do:
#     vertices_to = mne.grade_to_vertices(subject_to, grade=5)
# But fsaverage's source space was set up so we can just do this:
vertices_to = [np.arange(10242), np.arange(10242)]
stc_to = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=1,
                        grade=vertices_to, subjects_dir=subjects_dir)
stc_to.save('%s_audvis-meg' % subject_to)

# Morph using another method -- useful if you're going to do a lot of the
# same inter-subject morphing operations; you could save and load morph_mat
morph_mat = mne.compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertices, vertices_to,
                                     subjects_dir=subjects_dir)
stc_to_2 = mne.morph_data_precomputed(subject_from, subject_to,
                                      stc_from, vertices_to, morph_mat)
stc_to_2.save('%s_audvis-meg_2' % subject_to)

# View source activations
plt.plot(stc_from.times, stc_from.data.mean(axis=0), 'r', label='from')
plt.plot(stc_to.times, stc_to.data.mean(axis=0), 'b', label='to')
plt.plot(stc_to_2.times, stc_to.data.mean(axis=0), 'g', label='to_2')
plt.xlabel('time (ms)')
plt.ylabel('Mean Source amplitude')
plt.legend()
plt.show()

def morph_source_space(ndvar, subject_to, morph_mat=None, vertices_to=None):
    """Morph source estimate between subjects using a precomputed matrix

    Parameters
    ----------
    ndvar : NDVar
        NDVar with SourceSpace dimension.
    subject_to : string
        Name of the subject on which to morph.
    morph_mat : None | sparse matrix
        The morphing matrix. If ndvar contains a whole source space, the morph
        matrix can be automatically loaded, although providing a cached matrix
        can speed up processing by a second or two.
    vertices_to : None | list of array of int
        The vertices on the destination subject's brain. If ndvar contains a
        whole source space, vertices_to can be automatically loaded, although
        providing them as argument can speed up processing by a second or two.

    Returns
    -------
    morphed_ndvar : NDVar
        NDVar morphed to the destination subject.
    """
    src = ndvar.source.src
    subject_from = ndvar.source.subject
    subjects_dir = ndvar.source.subjects_dir
    vertices_from = ndvar.source.vertno
    if vertices_to is None:
        path = SourceSpace._src_pattern.format(subjects_dir=subjects_dir,
                                               subject=subject_to, src=src)
        src_to = mne.read_source_spaces(path)
        vertices_to = [src_to[0]['vertno'], src_to[1]['vertno']]
    elif not isinstance(vertices_to, list) or not len(vertices_to) == 2:
        raise ValueError('vertices_to must be a list of length 2')

    if morph_mat is None:
        morph_mat = mne.compute_morph_matrix(subject_from, subject_to,
                                             vertices_from, vertices_to, None,
                                             subjects_dir)
    elif not sp.sparse.issparse(morph_mat):
        raise ValueError('morph_mat must be a sparse matrix')
    elif not sum(len(v) for v in vertices_to) == morph_mat.shape[0]:
        raise ValueError('morph_mat.shape[0] must match number of vertices in '
                         'vertices_to')

    # flatten data
    axis = ndvar.get_axis('source')
    x = ndvar.x
    if axis != 0:
        x = x.swapaxes(0, axis)
    n_sources = len(x)
    if not n_sources == morph_mat.shape[1]:
        raise ValueError('ndvar source dimension length must be the same as '
                         'morph_mat.shape[0]')
    if ndvar.ndim > 2:
        shape = x.shape
        x = x.reshape((n_sources, -1))

    # apply morph matrix
    x_ = morph_mat * x

    # restore data shape
    if ndvar.ndim > 2:
        shape_ = (len(x_),) + shape[1:]
        x_ = x_.reshape(shape_)
    if axis != 0:
        x_ = x_.swapaxes(axis, 0)

    # package output NDVar
    parc = ndvar.source.parc
    source = SourceSpace(vertices_to, subject_to, src, subjects_dir, parc)
    dims = ndvar.dims[:axis] + (source,) + ndvar.dims[axis + 1:]
    info = ndvar.info.copy()
    out = NDVar(x_, dims, info, ndvar.name)
    return out

            continue
        raw_fname = paths('sss', subject=meg_subject, block=1)
        inv_fname = paths('inv', subject=meg_subject)
        cov = load('cov', subject=meg_subject)
        fwd = load('fwd', subject=meg_subject)
        info = read_info(raw_fname)
        inv = make_inverse_operator(info, fwd, cov, loose=0.2, depth=0.8)
        save(inv, 'inv', subject=meg_subject, overwrite=True)

# Morph anatomy into common model ---------------------------------------------
from mne import EvokedArray
from mne import compute_morph_matrix
from mne.minimum_norm import apply_inverse
if True:
    for meg_subject, subject in zip(range(1, 21), subjects_id):
        if subject in bad_mri:
            continue
        raw_fname = paths('sss', subject=meg_subject, block=1)
        info = read_info(raw_fname)

        # precompute morphing matrix for faster processing
        inv = load('inv', subject=meg_subject)
        evoked = EvokedArray(np.zeros((len(info['chs']), 2)), info, 0)
        evoked.pick_types(eeg=False, meg=True)
        stc = apply_inverse(evoked, inv, lambda2=1.0 / (2 ** 3.0),
                            method='dSPM', pick_ori=None)
        morph = compute_morph_matrix(subject, 'fsaverage', stc.vertices,
                                     vertices_to=[np.arange(10242)] * 2,
                                     subjects_dir=subjects_dir)
        save(morph, 'morph', subject=meg_subject)

print("Simulating data for %d subjects." % n_subjects)

#    Let's make sure our results replicate, so set the seed.
np.random.seed(0)
X = randn(n_vertices_sample, n_times, n_subjects, 4) * 10
for ii, condition in enumerate(conditions):
    X[:, :, :, ii] += condition.lh_data[:, :, np.newaxis]

#    It's a good idea to spatially smooth the data, and for visualization
#    purposes, let's morph these to fsaverage, which is a grade 5 source space
#    with vertices 0:10242 for each hemisphere. Usually you'd have to morph
#    each subject's data separately (and you might want to use morph_data
#    instead), but here since all estimates are on 'sample' we can use one
#    morph matrix for all the heavy lifting.
fsave_vertices = [np.arange(10242), np.array([])]  # right hemisphere is empty
morph_mat = compute_morph_matrix("sample", "fsaverage", sample_vertices, fsave_vertices, 20, subjects_dir)
n_vertices_fsave = morph_mat.shape[0]

#    We have to change the shape for the dot() to work properly
X = X.reshape(n_vertices_sample, n_times * n_subjects * 4)
print("Morphing data.")
X = morph_mat.dot(X)  # morph_mat is a sparse matrix
X = X.reshape(n_vertices_fsave, n_times, n_subjects, 4)

#    Now we need to prepare the group matrix for the ANOVA statistic.
#    To make the clustering function work correctly with the
#    ANOVA function X needs to be a list of multi-dimensional arrays
#    (one per condition) of shape: samples (subjects) x time x space

X = np.transpose(X, [2, 1, 0, 3])  # First we permute dimensions
# finally we split the array into a list a list of conditions

def xhemi(ndvar, mask=None, hemi='lh', parc=True):
    """Project data from both hemispheres to ``hemi`` of fsaverage_sym

    Project data from both hemispheres to the same hemisphere for
    interhemisphere comparisons. The fsaverage_sym brain is a symmetric
    version of fsaverage to facilitate interhemisphere comparisons. It is
    included with FreeSurfer > 5.1 and can be obtained as described `here
    <http://surfer.nmr.mgh.harvard.edu/fswiki/Xhemi>`_. For statistical
    comparisons between hemispheres, use of the symmetric ``fsaverage_sym``
    model is recommended to minimize bias [1]_.

    Parameters
    ----------
    ndvar : NDVar
        NDVar with SourceSpace dimension.
    mask : bool
        Remove sources in "unknown-" labels (default is True unless ``ndvar``
        contains sources with "unknown-" label).
    hemi : 'lh' | 'rh'
        Hemisphere onto which to morph the data.
    parc : bool | str
        Parcellation for target source space; True to use same as in ``ndvar``
        (default).

    Returns
    -------
    lh : NDVAr
        Data from the left hemisphere on ``hemi`` of ``fsaverage_sym``.
    rh : NDVar
        Data from the right hemisphere on ``hemi`` of ``fsaverage_sym``.

    Notes
    -----
    Only symmetric volume source spaces are currently supported.

    References
    ----------
    .. [1] Greve D. N., Van der Haegen L., Cai Q., Stufflebeam S., Sabuncu M.
           R., Fischl B., Brysbaert M.
           A Surface-based Analysis of Language Lateralization and Cortical
           Asymmetry. Journal of Cognitive Neuroscience 25(9), 1477-1492, 2013.
    """
    source = ndvar.get_dim('source')
    other_hemi = 'rh' if hemi == 'lh' else 'lh'

    if isinstance(source, VolumeSourceSpace):
        ax = ndvar.get_axis('source')
        # extract hemi
        is_in_hemi = source.hemi == hemi
        source_out = source[is_in_hemi]
        # map other_hemi into hemi
        is_in_other = source.hemi == other_hemi
        coord_map = {tuple(source.coordinates[i]): i for i in np.flatnonzero(is_in_other)}
        other_source = [coord_map[-x, y, z] for x, y, z in source_out.coordinates]
        # combine data
        x_hemi = ndvar.x[index(is_in_hemi, at=ax)]
        x_other = ndvar.x[index(other_source, at=ax)]
        dims = list(ndvar.dims)
        dims[ax] = source_out
        out_same = NDVar(x_hemi, dims, ndvar.info, ndvar.name)
        out_other = NDVar(x_other, dims, ndvar.info, ndvar.name)
    else:
        if source.subject == 'fsaverage_sym':
            ndvar_sym = ndvar
        else:
            ndvar_sym = morph_source_space(ndvar, 'fsaverage_sym', parc=parc, mask=mask)

        vert_lh, vert_rh = ndvar_sym.source.vertices
        vert_from = [[], vert_rh] if hemi == 'lh' else [vert_lh, []]
        vert_to = [vert_lh, []] if hemi == 'lh' else [[], vert_rh]
        with warnings.catch_warnings():
            warnings.filterwarnings('ignore', r'\d+/\d+ vertices not included in smoothing', module='mne')
            morph_mat = mne.compute_morph_matrix(
                'fsaverage_sym', 'fsaverage_sym', vert_from, vert_to,
                subjects_dir=ndvar.source.subjects_dir, xhemi=True)

        out_same = ndvar_sym.sub(source=hemi)
        out_other = morph_source_space(
            ndvar_sym.sub(source=other_hemi), 'fsaverage_sym',
            out_same.source.vertices, morph_mat, parc=parc, xhemi=True, mask=mask)

    if hemi == 'lh':
        return out_same, out_other
    else:
        return out_other, out_same

#.........这里部分代码省略.........
            x = ndvar.x.copy() if copy else ndvar.x
            parc_arg = None if parc is True else parc
            if subject_is_scaled or parc_arg is not None:
                source_to = source._copy(subject_to, parc=parc_arg)
            else:
                source_to = source

            dims = (*ndvar.dims[:axis], source_to, *ndvar.dims[axis + 1:])
            return NDVar(x, dims, ndvar.info, ndvar.name)

    has_lh_out = bool(source.rh_n if xhemi else source.lh_n)
    has_rh_out = bool(source.lh_n if xhemi else source.rh_n)
    if vertices_to in (None, 'lh', 'rh'):
        default_vertices = source_space_vertices(source.kind, source.grade, subject_to, subjects_dir)
        lh_out = vertices_to == 'lh' or (vertices_to is None and has_lh_out)
        rh_out = vertices_to == 'rh' or (vertices_to is None and has_rh_out)
        vertices_to = [default_vertices[0] if lh_out else np.empty(0, int),
                       default_vertices[1] if rh_out else np.empty(0, int)]
        if mask is None:
            if source.parc is None:
                mask = False
            else:
                mask = not np.any(source.parc.startswith('unknown-'))
    elif not isinstance(vertices_to, list) or not len(vertices_to) == 2:
        raise ValueError(f"vertices_to={vertices_to!r}: must be a list of length 2")

    # check that requested data is available
    n_to_lh = len(vertices_to[0])
    n_to_rh = len(vertices_to[1])
    if n_to_lh and not has_lh_out:
        raise ValueError("Data on the left hemisphere was requested in vertices_to but is not available in ndvar")
    elif n_to_rh and not has_rh_out:
        raise ValueError("Data on the right hemisphere was requested in vertices_to but is not available in ndvar")
    elif n_to_lh == 0 and n_to_rh == 0:
        raise ValueError("No target vertices")

    # parc for new source space
    if parc is True:
        parc_to = None if source.parc is None else source.parc.name
    else:
        parc_to = parc
    if mask and parc_to is None:
        raise ValueError("Can't mask source space without parcellation...")
    # check that annot files are available
    if parc_to:
        fname = SourceSpace._ANNOT_PATH.format(
            subjects_dir=subjects_dir, subject=subject_to, hemi='%s',
            parc=parc_to)
        fnames = tuple(fname % hemi for hemi in ('lh', 'rh'))
        missing = tuple(fname for fname in fnames if not os.path.exists(fname))
        if missing:
            missing = '\n'.join(missing)
            raise IOError(f"Annotation files are missing for parc={parc_to!r}, subject={subject_to!r}. Use the parc parameter when morphing to set a different parcellation. The following files are missing:\n{missing}")
    # find target source space
    source_to = SourceSpace(vertices_to, subject_to, source.src, subjects_dir, parc_to)
    if mask is True:
        if parc is True:
            keep_labels = source.parc.cells
            if xhemi:
                keep_labels = [switch_hemi_tag(label) for label in keep_labels]
            index = source_to.parc.isin(keep_labels)
        else:
            index = source_to.parc.isnotin(('unknown-lh', 'unknown-rh'))
        source_to = source_to[index]
    elif mask not in (None, False):
        raise TypeError(f"mask={mask!r}")

    if morph_mat is None:
        with warnings.catch_warnings():
            warnings.filterwarnings('ignore', r'\d+/\d+ vertices not included in smoothing', module='mne')
            morph_mat = mne.compute_morph_matrix(subject_from, subject_to, source.vertices, source_to.vertices, None, subjects_dir, xhemi=xhemi)
    elif not sp.sparse.issparse(morph_mat):
        raise ValueError('morph_mat must be a sparse matrix')
    elif not sum(len(v) for v in source_to.vertices) == morph_mat.shape[0]:
        raise ValueError('morph_mat.shape[0] must match number of vertices in vertices_to')

    # flatten data
    x = ndvar.x
    if axis != 0:
        x = x.swapaxes(0, axis)
    n_sources = len(x)
    if not n_sources == morph_mat.shape[1]:
        raise ValueError('ndvar source dimension length must be the same as morph_mat.shape[0]')
    if ndvar.ndim > 2:
        shape = x.shape
        x = x.reshape((n_sources, -1))

    # apply morph matrix
    x_m = morph_mat * x

    # restore data shape
    if ndvar.ndim > 2:
        shape_ = (len(x_m),) + shape[1:]
        x_m = x_m.reshape(shape_)
    if axis != 0:
        x_m = x_m.swapaxes(axis, 0)

    # package output NDVar
    dims = (*ndvar.dims[:axis], source_to, *ndvar.dims[axis + 1:])
    return NDVar(x_m, dims, ndvar.info, ndvar.name)

src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'

# Read input stc file
stc_from = mne.read_source_estimate(fname)
# Morph using one method (supplying the vertices in fsaverage's source
# space makes it faster). Note that for any generic subject, you could do:
#     vertices_to = mne.grade_to_vertices(subject_to, grade=5)
# But fsaverage's source space was set up so we can just do this:
vertices_to = [np.arange(10242), np.arange(10242)]
stc_to = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=1,
                        grade=vertices_to)
stc_to.save('%s_audvis-meg' % subject_to)

# Morph using another method -- useful if you're going to do a lot of the
# same inter-subject morphing operations; you could save and load morph_mat
morph_mat = mne.compute_morph_matrix(subject_from, subject_to, stc_from.vertno,
                                     vertices_to)
stc_to_2 = mne.morph_data_precomputed(subject_from, subject_to,
                                      stc_from, vertices_to, morph_mat)
stc_to_2.save('%s_audvis-meg_2' % subject_to)

# View source activations
import pylab as pl
pl.plot(stc_from.times, stc_from.data.mean(axis=0), 'r', label='from')
pl.plot(stc_to.times, stc_to.data.mean(axis=0), 'b', label='to')
pl.plot(stc_to_2.times, stc_to.data.mean(axis=0), 'g', label='to_2')
pl.xlabel('time (ms)')
pl.ylabel('Mean Source amplitude')
pl.legend()
pl.show()

epochs_ds.resample(20)
evoked_ds = [epochs_ds[name].average() for name in ['STI-correct', 'STI-incorrect']]

# contruct two types of inverse solution: one based on baseline data (before the red square appears), and one based on blank data
cov_blank = mne.compute_covariance(epochs_ds['STB'], tmin=0, tmax=None, method='auto')
inv_blank = make_inverse_operator(epochs_ds.info, forward, cov_blank,
                                  loose=0.2, depth=0.8)
blank_idx = np.nonzero(epochs_ds.events[:, 2] == 15)[0]
epochs_ds.drop_epochs(blank_idx)
cov_base = mne.compute_covariance(epochs_ds, tmin=None, tmax=0, method='auto')
inv_base = make_inverse_operator(epochs_ds.info, forward, cov_base,
                                 loose=0.2, depth=0.8)

vertices_to = [np.arange(10242), np.arange(10242)]
vertices_from = [forward['src'][0]['vertno'], forward['src'][1]['vertno']]
morph_mat = mne.compute_morph_matrix(subj, 'fsaverage', vertices_from, vertices_to)
for c in range(len(conds)):
    # start with the simplest method, MNE + dSPM
    stc = apply_inverse(evoked_ds[c], inv_base, lambda2, method)
    stc = mne.morph_data_precomputed(subj, 'fsaverage', stc,
                                     vertices_to, morph_mat)
    fname = out_dir + '%s_%s_dSPM_base_clean' % (subj, conds[c])
    stc.save(fname)
    stc = apply_inverse(evoked_ds[c], inv_blank, lambda2, method)
    stc = mne.morph_data_precomputed(subj, 'fsaverage', stc,
                                     vertices_to, morph_mat)
    fname = out_dir + '%s_%s_dSPM_blank_clean' % (subj, conds[c])
    stc.save(fname)

    # the next estimate is LCMV beamformer in time
    data_cov = mne.compute_covariance(epochs_ds[conds[c]], tmin=0, tmax=None,

        inv = op.join(inv_dir, '%s-%d-sss-%s-inv.fif' % (subj, p.lp_cut, inv_type))
        inv = read_inverse_operator(inv)
        fname = op.join(inv_dir, '%s_%d-sss_eq_%s-ave.fif'
                        % (p.analyses[0], p.lp_cut, subj))
        aves = [mne.Evoked(fname, cond, baseline=(None, 0), proj=True,
                                kind='average') for cond in conditions]
        nave = np.unique([a.nave for a in aves])
        assert len(nave) == 1
        for ave, cond in zip(aves, conditions):
                assert ave.comment == cond
        naves[si] = nave[0]

        # apply inverse, bin, morph
        stcs = [apply_inverse(ave, inv, lambda2, 'dSPM') for ave in aves]
        stcs = [stc.bin(0.005) for stc in stcs]
        m = mne.compute_morph_matrix(struc, 'fsaverage', stcs[0].vertices,
                                     fs_verts, n_smooth)
        stcs = [stc.morph_precomputed('fsaverage', fs_verts, m)
                for stc in stcs]

        # put in big matrix
        if subj == p.subjects[0]:
            data = np.empty((len(stcs), len(p.subjects), stcs[0].shape[0],
                             stcs[0].shape[1]))
        for di, stc in enumerate(stcs):
            data[di, si, :, :] = stc.data
            times = stc.times
        print('Writing data...')
        np.savez_compressed(fname_data, data=data, times=times, naves=naves)
else:
    print('Loading saved data...')
    data = np.load(fname_data)

print('Simulating data for %d subjects.' % n_subjects)

#    Let's make sure our results replicate, so set the seed.
np.random.seed(0)
X = randn(n_vertices_sample, n_times, n_subjects, 4) * 10
for ii, condition in enumerate(conditions):
    X[:, :, :, ii] += condition.lh_data[:, :, np.newaxis]

#    It's a good idea to spatially smooth the data, and for visualization
#    purposes, let's morph these to fsaverage, which is a grade 5 source space
#    with vertices 0:10242 for each hemisphere. Usually you'd have to morph
#    each subject's data separately (and you might want to use morph_data
#    instead), but here since all estimates are on 'sample' we can use one
#    morph matrix for all the heavy lifting.
fsave_vertices = [np.arange(10242), np.array([], int)]  # right hemi is empty
morph_mat = compute_morph_matrix('sample', 'fsaverage', sample_vertices,
                                 fsave_vertices, 20, subjects_dir)
n_vertices_fsave = morph_mat.shape[0]

#    We have to change the shape for the dot() to work properly
X = X.reshape(n_vertices_sample, n_times * n_subjects * 4)
print('Morphing data.')
X = morph_mat.dot(X)  # morph_mat is a sparse matrix
X = X.reshape(n_vertices_fsave, n_times, n_subjects, 4)

#    Now we need to prepare the group matrix for the ANOVA statistic.
#    To make the clustering function work correctly with the
#    ANOVA function X needs to be a list of multi-dimensional arrays
#    (one per condition) of shape: samples (subjects) x time x space