def test_morph_data():
    """Test morphing of data
    """
    subject_from = 'sample'
    subject_to = 'fsaverage'
    fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
    stc_from = read_source_estimate(fname)
    stc_from.crop(0.09, 0.1)  # for faster computation
    # After running this:
    #    stc_from.save('%s_audvis-meg-cropped' % subject_from)
    # this was run from a command line:
    #    mne_make_movie --stcin sample_audvis-meg-cropped-lh.stc
    #        --subject sample --morph fsaverage --smooth 12 --morphgrade 3
    #        --stc fsaverage_audvis-meg-cropped
    # XXX These files should eventually be moved to the sample dataset and
    # removed from mne/fiff/tests/data/
    fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
                    'fsaverage_audvis-meg-cropped')
    stc_to = read_source_estimate(fname)
    stc_to1 = morph_data(subject_from, subject_to, stc_from,
                            grade=3, smooth=12, buffer_size=1000)
    stc_to1.save('%s_audvis-meg' % subject_to)
    stc_to2 = morph_data(subject_from, subject_to, stc_from,
                            grade=3, smooth=12, buffer_size=3)
    # indexing silliness here due to mne_make_movie's indexing oddities
    assert_array_almost_equal(stc_to.data, stc_to1.data[:, 0][:, None], 5)
    assert_array_almost_equal(stc_to1.data, stc_to2.data)
    # make sure precomputed morph matrices work
    vertices_to = grade_to_vertices(subject_to, grade=3)
    morph_mat = compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertno, vertices_to,
                                     smooth=12)
    stc_to3 = morph_data_precomputed(subject_from, subject_to,
                                     stc_from, vertices_to, morph_mat)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to1.data.mean(axis=0)
    assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)

    # test two types of morphing:
    # 1) make sure we can fill by morphing
    stc_to5 = morph_data(subject_from, subject_to, stc_from,
                            grade=None, smooth=12, buffer_size=3)
    assert_true(stc_to5.data.shape[0] == 163842 + 163842)

    # 2) make sure we can specify vertices
    vertices_to = [np.arange(10242), np.arange(10242)]
    stc_to3 = morph_data(subject_from, subject_to, stc_from,
                            grade=vertices_to, smooth=12, buffer_size=3)
    stc_to4 = morph_data(subject_from, subject_to, stc_from,
                            grade=5, smooth=12, buffer_size=3)
    assert_array_almost_equal(stc_to3.data, stc_to4.data)

def test_morph_data():
    """Test morphing of data
    """
    subject_from = 'sample'
    subject_to = 'fsaverage'
    fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
    stc_from = read_source_estimate(fname, subject='sample')
    fname = op.join(data_path, 'MEG', 'sample', 'fsaverage_audvis-meg')
    stc_to = read_source_estimate(fname)
    # make sure we can specify grade
    stc_from.crop(0.09, 0.1)  # for faster computation
    stc_to.crop(0.09, 0.1)  # for faster computation
    stc_to1 = stc_from.morph(subject_to, grade=3, smooth=12, buffer_size=1000,
                             subjects_dir=subjects_dir)
    stc_to1.save(op.join(tempdir, '%s_audvis-meg' % subject_to))
    # make sure we can specify vertices
    vertices_to = grade_to_vertices(subject_to, grade=3)
    stc_to2 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=1000,
                         subjects_dir=subjects_dir)
    # make sure we can use different buffer_size
    stc_to3 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=3,
                         subjects_dir=subjects_dir)

    assert_array_almost_equal(stc_to.data, stc_to1.data, 5)
    assert_array_almost_equal(stc_to1.data, stc_to2.data)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)
    # make sure precomputed morph matrices work
    morph_mat = compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertno, vertices_to,
                                     smooth=12, subjects_dir=subjects_dir)
    stc_to3 = stc_from.morph_precomputed(subject_to, vertices_to, morph_mat)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to1.data.mean(axis=0)
    assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)

    # make sure we can fill by morphing
    stc_to5 = morph_data(subject_from, subject_to, stc_from, grade=None,
                         smooth=12, buffer_size=3, subjects_dir=subjects_dir)
    assert_true(stc_to5.data.shape[0] == 163842 + 163842)

    # test morphing to the same subject
    stc_to6 = stc_from.morph(subject_from, grade=stc_from.vertno, smooth=1,
                             subjects_dir=subjects_dir)
    mask = np.ones(stc_from.data.shape[0], dtype=np.bool)
    # XXX: there is a bug somewhere that causes a difference at 2 vertices..
    mask[6799] = False
    mask[6800] = False
    assert_array_almost_equal(stc_from.data[mask], stc_to6.data[mask], 5)

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 apply_inverse(fnevo, method='dSPM', snr=3.0, event='LLst', 
                  baseline=False, btmin=-0.3, btmax=-0.1, min_subject='fsaverage'):
    '''  
        Parameter
        ---------
        fnevo: string or list
            The evoked file with ECG, EOG and environmental noise free.
        method: inverse method, 'MNE' or 'dSPM'
        event: string
            The event name related with epochs.
        min_subject: string
            The subject name as the common brain.
        snr: signal to noise ratio for inverse solution. 
    '''
    #Get the default subjects_dir
    from mne.minimum_norm import apply_inverse
    fnlist = get_files_from_list(fnevo)
    # loop across all filenames
    for fname in fnlist:
        fn_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        stc_name = name[:name.rfind('-ave.fif')] 
        subject = name.split('_')[0]
        subject_path = subjects_dir + '/%s' %subject
        min_dir = subjects_dir + '/%s' %min_subject
        fn_trans = fn_path + '/%s-trans.fif' % subject
        fn_cov = fn_path + '/%s_empty,nr-cov.fif' % subject
        fn_src = subject_path + '/bem/%s-ico-5-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        snr = snr
        lambda2 = 1.0 / snr ** 2 
        #noise_cov = mne.read_cov(fn_cov)
        [evoked] = mne.read_evokeds(fname)
        noise_cov = mne.read_cov(fn_cov)
        # this path used for ROI definition
        stc_path = min_dir + '/%s_ROIs/%s' %(method,subject)
        #fn_cov = meg_path + '/%s_empty,fibp1-45,nr-cov.fif' % subject
        set_directory(stc_path)
        noise_cov = mne.cov.regularize(noise_cov, evoked.info,
                                        mag=0.05, grad=0.05, proj=True)
        fwd_ev = mne.make_forward_solution(evoked.info, trans=fn_trans,
                                            src=fn_src, bem=fn_bem,
                                            fname=None, meg=True, eeg=False,
                                            mindist=5.0, n_jobs=2,
                                            overwrite=True)
        fwd_ev = mne.convert_forward_solution(fwd_ev, surf_ori=True)
        forward_meg_ev = mne.pick_types_forward(fwd_ev, meg=True, eeg=False)
        inverse_operator_ev = mne.minimum_norm.make_inverse_operator(
            evoked.info, forward_meg_ev, noise_cov,
            loose=0.2, depth=0.8)
        # Compute inverse solution
        stc = apply_inverse(evoked, inverse_operator_ev, lambda2, method,
                            pick_ori=None)
        # Morph STC
        stc_morph = mne.morph_data(subject, min_subject, stc, grade=5, smooth=5)
        stc_morph.save(stc_path + '/%s' % (stc_name), ftype='stc')
        if baseline == True:
            stc_base = stc_morph.crop(btmin, btmax)
            stc_base.save(stc_path + '/%s_%s_baseline' % (subject, event), ftype='stc')

def test_morph_data():
    """Test morphing of data
    """
    subject_from = 'sample'
    subject_to = 'fsaverage'
    fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
    stc_from = SourceEstimate(fname)
    stc_from.crop(0.09, 0.1)  # for faster computation
    stc_to = morph_data(subject_from, subject_to, stc_from,
                            grade=3, smooth=12, buffer_size=1000)
    stc_to.save('%s_audvis-meg' % subject_to)

    stc_to2 = morph_data(subject_from, subject_to, stc_from,
                            grade=3, smooth=12, buffer_size=3)
    assert_array_almost_equal(stc_to.data, stc_to2.data)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to.data.mean(axis=0)
    assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)

def test_morph_data():
    """Test morphing of data
    """
    subject_from = 'sample'
    subject_to = 'fsaverage'
    fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
    stc_from = read_source_estimate(fname, subject='sample')
    fname = op.join(data_path, 'MEG', 'sample', 'fsaverage_audvis-meg')
    stc_to = read_source_estimate(fname)
    # make sure we can specify grade
    stc_from.crop(0.09, 0.1)  # for faster computation
    stc_to.crop(0.09, 0.1)  # for faster computation
    stc_to1 = stc_from.morph(subject_to, grade=3, smooth=12, buffer_size=1000,
                             subjects_dir=subjects_dir)
    stc_to1.save(op.join(tempdir, '%s_audvis-meg' % subject_to))
    # make sure we can specify vertices
    vertices_to = grade_to_vertices(subject_to, grade=3)
    stc_to2 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=1000,
                         subjects_dir=subjects_dir)
    # make sure we can use different buffer_size
    stc_to3 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=3,
                         subjects_dir=subjects_dir)
    # indexing silliness here due to mne_make_movie's indexing oddities
    assert_array_almost_equal(stc_to.data, stc_to1.data, 5)
    assert_array_almost_equal(stc_to1.data, stc_to2.data)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)
    # make sure precomputed morph matrices work
    morph_mat = compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertno, vertices_to,
                                     smooth=12, subjects_dir=subjects_dir)
    stc_to3 = stc_from.morph_precomputed(subject_to, vertices_to, morph_mat)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to1.data.mean(axis=0)
    assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)

    # make sure we can fill by morphing
    stc_to5 = morph_data(subject_from, subject_to, stc_from, grade=None,
                         smooth=12, buffer_size=3, subjects_dir=subjects_dir)
    assert_true(stc_to5.data.shape[0] == 163842 + 163842)

def _morphed_epochs_files(params):
    subject, cond_name, stc_file_name, inverse_method, subjects_dir = params
    print('morphing {}'.format(stc_file_name))
    epoch_id = utils.namebase(stc_file_name).split('_')[2]
    morphed_stc_file_name = op.join(LOCAL_ROOT_DIR, 'stc_epochs_morphed',  '{}_{}_{}_{}'.format(subject, cond_name, epoch_id, inverse_method))
    if not op.isfile('{}-stc.h5'.format(morphed_stc_file_name)):
        stc = mne.read_source_estimate(stc_file_name)
        stc_morphed = mne.morph_data(subject, 'fsaverage', stc, grade=5, smooth=20,
            subjects_dir=subjects_dir)
        stc_morphed.save(morphed_stc_file_name, ftype='h5')
    else:
        print('{} {} {} already morphed'.format(subject, cond_name, epoch_id))

def test_morph_data():
    """Test morphing of data
    """
    import mne
    subject_from = 'sample'
    subject_to = 'morph'
    fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
    stc_from = mne.SourceEstimate(fname)
    stc_to = mne.morph_data(subject_from, subject_to, stc_from, 3)

    stc_to.save('%s_audvis-meg' % subject_to)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to.data.mean(axis=0)
    assert np.corrcoef(mean_to, mean_from).min() > 0.99

def smooth_ttest_results(tmin, tstep, subjects_dir, inverse_method='dSPM', n_jobs=1):
    for cond_id, cond_name in enumerate(events_id.keys()):
        for patient in get_patients():
            results_file_name = op.join(LOCAL_ROOT_DIR, 'permutation_ttest_results', '{}_{}_{}_clusters.npy'.format(patient, cond_name, inverse_method))
            if op.isfile(results_file_name):
                data = np.load(results_file_name)
                print('smoothing {} {}'.format(patient, cond_name))
                fsave_vertices = [np.arange(10242), np.arange(10242)]
                stc = _make_stc(data, fsave_vertices, tmin=tmin, tstep=tstep, subject='fsaverage')
                vertices_to = mne.grade_to_vertices('fsaverage', grade=None, subjects_dir=subjects_dir)
                print(stc.data.shape, vertices_to[0].shape)
                stc_smooth = mne.morph_data('fsaverage', 'fsaverage', stc, n_jobs=n_jobs, grade=vertices_to, subjects_dir=subjects_dir)
                stc_smooth.save(op.join(LOCAL_ROOT_DIR, 'results_for_blender', '{}_{}_{}'.format(patient, cond_name, inverse_method)), ftype='h5')
            else:
                print('no results for {} {}'.format(patient, cond_name))

def morph_STC(fn_stc, grade, template='fsaverage', event='LLst', 
              baseline=True, btmin=-0.3, btmax=0.):
    from mne import read_source_estimate, morph_data
    fnlist = get_files_from_list(fn_stc) 
    for fname in fnlist:  
        name = os.path.basename(fname)
        subject = name.split('_')[0]
        stc_name = name[:name.rfind('-ave.fif')] 
        min_dir = subjects_dir + '/%s' %template
        # this path used for ROI definition
        stc_path = min_dir + '/dSPM_ROIs/%s' %(subject)
        #fn_cov = meg_path + '/%s_empty,fibp1-45,nr-cov.fif' % subject
        set_directory(stc_path) 
        # Morph STC
        stc = read_source_estimate(fname)
        stc_morph = morph_data(subject, template, stc, grade=grade)
        stc_morph.save(stc_path + '/%s' % (stc_name), ftype='stc')
        if baseline == True:
            stc_base = stc_morph.crop(btmin, btmax)
            stc_base.save(stc_path + '/%s_%s_baseline' % (subject, event[:2]), ftype='stc') 

def morph_STC(fn_stc, grade, subjects_dir, template='fsaverage', event='LLst',
              baseline=True, btmin=-0.3, btmax=0.):
    '''
        Morph individual STC into the common brain space.

        Parameter
        ------------------------------------
        fn_stc: string or list
            The path of the individual STC.
        subjects_dir: The total bath of all the subjects.
        template: string
            The subject name as the common brain.
        event: string
            The name of event
        baseline: bool
            If true, prestimulus segment from 'btmin' to 'btmax' will be saved,
            If false, no baseline segment is saved.
        btmin, btmax: float
            If 'baseline' is True, baseline is croped using this period.

    '''
    from mne import read_source_estimate, morph_data
    fnlist = get_files_from_list(fn_stc)
    for fname in fnlist:
        name = os.path.basename(fname)
        subject = name.split('_')[0]
        stc_name = name[:name.rfind('-lh.stc')]
        min_dir = subjects_dir + '/%s' % template
        # this path used for ROI definition
        stc_path = min_dir + '/dSPM_ROIs/%s' % (subject)
        # fn_cov = meg_path + '/%s_empty,fibp1-45,nr-cov.fif' % subject
        set_directory(stc_path)
        # Morph STC
        stc = read_source_estimate(fname)
        stc_morph = morph_data(subject, template, stc, grade=grade)
        stc_morph.save(stc_path + '/%s' % (stc_name), ftype='stc')
        if baseline is True:
            stc_base = stc_morph.crop(btmin, btmax)
            stc_base.save(stc_path + '/%s_%s_baseline' % (subject, event[:2]),
                          ftype='stc')

    def __call__(self, data):

        output = deepcopy(data)

        from_surf_file = dirname(dirname(self.from_surf.surf_file))
        SUBJECTS_DIR, from_surf_name = split(from_surf_file)

        if 'lh' in self.from_surf.surf_file:  # TODO: not good, we need to check
            vertices = [arange(data.data[0].shape[0]), arange(0)]
        else:
            vertices = [arange(0), arange(data.data[0].shape[0])]

        stc = SourceEstimate(atleast_2d(data.data[0]).T, vertices=vertices,
                             tstep=0, tmin=0)
        m = morph_data(from_surf_name, self.to_surf, stc,
                       subjects_dir=SUBJECTS_DIR, grade=None,
                       smooth=self.smooth, verbose=False)

        output.data[0] = squeeze(m.data, axis=1)
        output.axis['surf'][0] = arange(m.data.shape[0])

        return output

def morph_stcs_to_fsaverage(events_id, stc_per_epoch=False, inverse_method='dSPM', subjects_dir='', n_jobs=1):
    if subjects_dir is '':
        subjects_dir = os.environ['SUBJECTS_DIR']
    for subject in get_subjects():
        for cond_name in events_id.keys():
            print('morphing {}, {}'.format(subject, cond_name))
            if not stc_per_epoch:
                morphed_stc_file_name = op.join(LOCAL_ROOT_DIR, 'stc_morphed', '{}_{}_morphed_{}'.format(subject, cond_name, inverse_method))
                if op.isfile('{}-stc.h5'.format(morphed_stc_file_name)):
                    print('{} {} already morphed'.format(subject, cond_name))
                else:
                    local_stc_file_name = op.join(LOCAL_ROOT_DIR, 'stc', '{}_{}_{}'.format(subject, cond_name, inverse_method))
                    if op.isfile('{}-stc.h5'.format(local_stc_file_name)):
                        stc = mne.read_source_estimate(local_stc_file_name)
                        stc_morphed = mne.morph_data(subject, 'fsaverage', stc, grade=5, smooth=20,
                            subjects_dir=subjects_dir)
                        stc_morphed.save(morphed_stc_file_name, ftype='h5')
                    else:
                        print("can't find stc file for {}, {}".format(subject, cond_name))
            else:
                stcs = glob.glob(op.join(LOCAL_ROOT_DIR, 'stc_epochs', '{}_{}_*_{}-stc.h5'.format(subject, cond_name, inverse_method)))
                params = [(subject, cond_name, stc_file_name, inverse_method, subjects_dir) for stc_file_name in stcs]
                utils.parallel_run(pool, _morphed_epochs_files, params, n_jobs)

def apply_inverse(fn_epo, event='LLst',ctmin=0.05, ctmax=0.25, nctmin=-0.2, nctmax=0,
                  fmin=4, fmax=8, min_subject='fsaverage', STCs=False):
    """
    Inverse evokes into source space using DICS method.
    ----------
    fn_epo : epochs of raw data.
    event_id: event id related with epochs.
    ctmin: the min time for computing CSD
    ctmax: the max time for computing CSD
    fmin: min value of the interest frequency band
    fmax: max value of the interest frequency band 
    min_subject: the subject for the common brain space.
    STCs: bool, make STCs of epochs.
    """
    from mne import Epochs, pick_types
    from mne.io import Raw
    from mne.event import make_fixed_length_events
    fnlist = get_files_from_list(fn_epo)
    # loop across all filenames
    for fname in fnlist:
        subjects_dir = os.environ['SUBJECTS_DIR']
        # extract the subject infromation from the file name
        meg_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        stc_name = name[:name.rfind('-epo.fif')] 
        subject = name.split('_')[0]
        subject_path = subjects_dir + '/%s' %subject
        min_dir = subjects_dir + '/%s' %min_subject
        fn_trans = meg_path + '/%s-trans.fif' % subject
        fn_src = subject_path + '/bem/%s-ico-4-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        # Make sure the target path is exist
        stc_path = min_dir + '/DICS_ROIs/%s' % subject
        set_directory(stc_path)
        # Read the MNI source space
        epochs = mne.read_epochs(fname)
        evoked = epochs.average()
        forward = mne.make_forward_solution(epochs.info, trans=fn_trans,
                                            src=fn_src, bem=fn_bem,
                                            fname=None, meg=True, eeg=False,
                                            mindist=5.0, n_jobs=2,
                                            overwrite=True)
        forward = mne.convert_forward_solution(forward, surf_ori=True)
        from mne.time_frequency import compute_epochs_csd
        from mne.beamformer import dics
        data_csd = compute_epochs_csd(epochs, mode='multitaper', tmin=ctmin, tmax=ctmax, 
                                      fmin=fmin, fmax=fmax)

        noise_csd = compute_epochs_csd(epochs, mode='multitaper', tmin=nctmin, tmax=nctmax,
                                           fmin=fmin, fmax=fmax)
                
        stc = dics(evoked, forward, noise_csd, data_csd)
        from mne import morph_data
        stc_morph = morph_data(subject, min_subject, stc, grade=4, smooth=4)
        stc_morph.save(stc_path + '/%s_%d_%d' % (stc_name, fmin, fmax), ftype='stc')
        if STCs == True:
            stcs_path = stc_path + '/STCs-%s/' %event
            reset_directory(stcs_path)
            stcs = dics(epochs, forward, noise_csd, data_csd)
            s = 0
            while s < len(stcs):
                stc_morph = mne.morph_data(subject, min_subject, stcs[s], grade=4, smooth=4)
                stc_morph.save(stcs_path + '/trial_%s'
                                % (str(s)), ftype='stc')
                s = s + 1

def morph_stc(subject_from, subject_to, stc_from_file):
    stc_from = mne.read_source_estimate(stc_from_file)
    vertices_to = [np.arange(10242), np.arange(10242)]
    stc_to = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=4, grade=vertices_to)
    stc_to.save("{}_{}.stc".format(stc_from_file[:-4], subject_to))

from mne.stats import spatio_temporal_cluster_test, summarize_clusters_stc
from mne.datasets import sample

print(__doc__)

###############################################################################
# Set parameters
data_path = sample.data_path()
stc_fname = data_path + "/MEG/sample/sample_audvis-meg-lh.stc"
subjects_dir = data_path + "/subjects"

# Load stc to in common cortical space (fsaverage)
stc = mne.read_source_estimate(stc_fname)
stc.resample(50)

stc = mne.morph_data("sample", "fsaverage", stc, grade=5, smooth=20, subjects_dir=subjects_dir)
n_vertices_fsave, n_times = stc.data.shape
tstep = stc.tstep

n_subjects1, n_subjects2 = 7, 9
print("Simulating data for %d and %d subjects." % (n_subjects1, n_subjects2))

#    Let's make sure our results replicate, so set the seed.
np.random.seed(0)
X1 = np.random.randn(n_vertices_fsave, n_times, n_subjects1) * 10
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

def test_morph_data():
    """Test morphing of data
    """
    tempdir = _TempDir()
    subject_from = 'sample'
    subject_to = 'fsaverage'
    stc_from = read_source_estimate(fname_smorph, subject='sample')
    stc_to = read_source_estimate(fname_fmorph)
    # make sure we can specify grade
    stc_from.crop(0.09, 0.1)  # for faster computation
    stc_to.crop(0.09, 0.1)  # for faster computation
    assert_array_equal(stc_to.time_as_index([0.09, 0.1], use_rounding=True),
                       [0, len(stc_to.times) - 1])
    assert_raises(ValueError, stc_from.morph, subject_to, grade=3, smooth=-1,
                  subjects_dir=subjects_dir)
    stc_to1 = stc_from.morph(subject_to, grade=3, smooth=12, buffer_size=1000,
                             subjects_dir=subjects_dir)
    stc_to1.save(op.join(tempdir, '%s_audvis-meg' % subject_to))
    # Morphing to a density that is too high should raise an informative error
    # (here we need to push to grade=6, but for some subjects even grade=5
    # will break)
    assert_raises(ValueError, stc_to1.morph, subject_from, grade=6,
                  subjects_dir=subjects_dir)
    # make sure we can specify vertices
    vertices_to = grade_to_vertices(subject_to, grade=3,
                                    subjects_dir=subjects_dir)
    stc_to2 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=1000,
                         subjects_dir=subjects_dir)
    # make sure we can use different buffer_size
    stc_to3 = morph_data(subject_from, subject_to, stc_from,
                         grade=vertices_to, smooth=12, buffer_size=3,
                         subjects_dir=subjects_dir)
    # make sure we get a warning about # of steps
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        morph_data(subject_from, subject_to, stc_from,
                   grade=vertices_to, smooth=1, buffer_size=3,
                   subjects_dir=subjects_dir)
    assert_equal(len(w), 2)

    assert_array_almost_equal(stc_to.data, stc_to1.data, 5)
    assert_array_almost_equal(stc_to1.data, stc_to2.data)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)
    # make sure precomputed morph matrices work
    morph_mat = compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertices, vertices_to,
                                     smooth=12, subjects_dir=subjects_dir)
    stc_to3 = stc_from.morph_precomputed(subject_to, vertices_to, morph_mat)
    assert_array_almost_equal(stc_to1.data, stc_to3.data)
    assert_raises(ValueError, stc_from.morph_precomputed,
                  subject_to, vertices_to, 'foo')
    assert_raises(ValueError, stc_from.morph_precomputed,
                  subject_to, [vertices_to[0]], morph_mat)
    assert_raises(ValueError, stc_from.morph_precomputed,
                  subject_to, [vertices_to[0][:-1], vertices_to[1]], morph_mat)
    assert_raises(ValueError, stc_from.morph_precomputed, subject_to,
                  vertices_to, morph_mat, subject_from='foo')

    # steps warning
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        compute_morph_matrix(subject_from, subject_to,
                             stc_from.vertices, vertices_to,
                             smooth=1, subjects_dir=subjects_dir)
    assert_equal(len(w), 2)

    mean_from = stc_from.data.mean(axis=0)
    mean_to = stc_to1.data.mean(axis=0)
    assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)

    # make sure we can fill by morphing
    stc_to5 = morph_data(subject_from, subject_to, stc_from, grade=None,
                         smooth=12, buffer_size=3, subjects_dir=subjects_dir)
    assert_true(stc_to5.data.shape[0] == 163842 + 163842)

    # Morph sparse data
    # Make a sparse stc
    stc_from.vertices[0] = stc_from.vertices[0][[100, 500]]
    stc_from.vertices[1] = stc_from.vertices[1][[200]]
    stc_from._data = stc_from._data[:3]

    assert_raises(RuntimeError, stc_from.morph, subject_to, sparse=True,
                  grade=5, subjects_dir=subjects_dir)

    stc_to_sparse = stc_from.morph(subject_to, grade=None, sparse=True,
                                   subjects_dir=subjects_dir)
    assert_array_almost_equal(np.sort(stc_from.data.sum(axis=1)),
                              np.sort(stc_to_sparse.data.sum(axis=1)))
    assert_equal(len(stc_from.rh_vertno), len(stc_to_sparse.rh_vertno))
    assert_equal(len(stc_from.lh_vertno), len(stc_to_sparse.lh_vertno))
    assert_equal(stc_to_sparse.subject, subject_to)
    assert_equal(stc_from.tmin, stc_from.tmin)
    assert_equal(stc_from.tstep, stc_from.tstep)

    stc_from.vertices[0] = np.array([], dtype=np.int64)
    stc_from._data = stc_from._data[:1]

    stc_to_sparse = stc_from.morph(subject_to, grade=None, sparse=True,
                                   subjects_dir=subjects_dir)
#.........这里部分代码省略.........

        # downsample to 1 Hz effective sampling resolution. Use this instead of raw.resample() to avoid further filtering of the signal
        data = data[:, np.arange(0,data.shape[1],int(raw.info['sfreq']))]
        print 'Multiplying data by beamformer weights...'
        # get the abs() of Hilbert transform (Hilbert envelope)
        sol = abs(np.dot(weights, data))
        stc = mne.SourceEstimate(sol, [forward['src'][0]['vertno'], forward['src'][1]['vertno']], 0, 1, subject=subj)
        stc.save(dir_out + 'lcmv-%dto%d-'%(l_freq,h_freq) + subj)

# morph all subjects
subject_to = 'fsaverage'
for l_freq, h_freq in bands:
    for subj in subjs:
        fname = dir_out + 'lcmv-%dto%d-'%(l_freq,h_freq) + subj
        stc_from = mne.read_source_estimate(fname)
        vertices_to = [np.arange(10242), np.arange(10242)]
        stc = mne.morph_data(subj, subject_to, stc_from, grade=vertices_to)
        stc.save(dir_out + 'morphed-lcmv-%dto%d-'%(l_freq,h_freq) + subj)

# concatenate all subjects and apply ICA
band_ICs = []
band_corr_ICs = []
for l_freq, h_freq in bands:
    print 'Concatenating sources in band %d to %d Hz'%(l_freq, h_freq)
    init_sources = 20500
    init_time = 38500
    # create huge array so we can add all the data and then resize it appropriately
    data = np.empty([init_sources, init_time])
    data[:] = np.nan
    cnt = 0
    for subj in subjs:
        fname = dir_out + 'morphed-lcmv-%dto%d-'%(l_freq,h_freq) + subj

def DICS_inverse(fn_epo, event_id=1,event='LLst', ctmin=0.05, ctmax=0.25, fmin=4, fmax=8, 
                  min_subject='fsaverage'):
    """
    Inverse evokes into source space using DICS method.
    ----------
    fn_epo : epochs of raw data.
    event_id: event id related with epochs.
    ctmin: the min time for computing CSD
    ctmax: the max time for computing CSD
    fmin: min value of the interest frequency band
    fmax: max value of the interest frequency band 
    min_subject: the subject for the common brain space.
    save_forward: Whether save the forward solution or not.
    """
    from mne import Epochs, pick_types
    from mne.io import Raw
    from mne.event import make_fixed_length_events
    fnlist = get_files_from_list(fn_epo)
    # loop across all filenames
    for fname in fnlist:
        meg_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        stc_name = name[:name.rfind('-epo.fif')] 
        subject = name.split('_')[0]
        subject_path = subjects_dir + '/%s' %subject
        min_dir = subjects_dir + '/%s' %min_subject
        fn_trans = meg_path + '/%s-trans.fif' % subject
        fn_src = subject_path + '/bem/%s-ico-5-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        # Make sure the target path is exist
        stc_path = min_dir + '/DICS_ROIs/%s' % subject
        set_directory(stc_path)
        # Read the MNI source space
        epochs = mne.read_epochs(fname)
        tmin = epochs.times.min()
        tmax = epochs.times.max()
        fn_empty = meg_path + '/%s_empty,nr-raw.fif' % subject
        raw_noise = Raw(fn_empty, preload=True)
        epochs.info['bads'] = raw_noise.info['bads']
        picks_noise = pick_types(raw_noise.info, meg='mag', exclude='bads')
        events_noise = make_fixed_length_events(raw_noise, event_id, duration=1.)
        epochs_noise = Epochs(raw_noise, events_noise, event_id, tmin,
                                tmax, proj=True, picks=picks_noise,
                                baseline=None, preload=True, reject=None)
        # Make sure the number of noise epochs is the same as data epochs
        epochs_noise = epochs_noise[:len(epochs.events)]
        evoked = epochs.average()
        forward = mne.make_forward_solution(epochs.info, trans=fn_trans,
                                            src=fn_src, bem=fn_bem,
                                            fname=None, meg=True, eeg=False,
                                            mindist=5.0, n_jobs=2,
                                            overwrite=True)
        forward = mne.convert_forward_solution(forward, surf_ori=True)
        from mne.time_frequency import compute_epochs_csd
        from mne.beamformer import dics
        data_csd = compute_epochs_csd(epochs, mode='multitaper', tmin=ctmin, tmax=ctmax, 
                                      fmin=fmin, fmax=fmax)

        noise_csd = compute_epochs_csd(epochs_noise, mode='multitaper', tmin=ctmin, tmax=ctmax,
                                           fmin=fmin, fmax=fmax)
                
        stc = dics(evoked, forward, noise_csd, data_csd)
        from mne import morph_data
        stc_morph = morph_data(subject, min_subject, stc, grade=5, smooth=5)
        stc_morph.save(stc_path + '/%s_%d_%d' % (event, fmin, fmax), ftype='stc')

raw.save(fname[:-4] + '_beta.fif')

# extract epochs
picks = mne.fiff.pick_types(raw.info, meg=True, eeg=True, eog=True,
                            exclude=raw.info['bads'])
events = mne.find_events(raw, stim_channel='STI 014')
epochs = mne.Epochs(raw, events, event_id=1, tmin=-0.2, tmax=0.5, proj=True,
                    picks=picks, baseline=(None, 0), preload=True,
                    reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))

# compute evoked response and noise covariance
evoked = epochs.average()
cov = mne.compute_covariance(epochs, tmax=0)

# Plot evoked
from mne.viz import plot_evoked
plot_evoked(evoked)

# Inverse modeling

# compute inverse operator
fwd_fname = 'sample_audvis-meg-eeg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
inv = mne.minimum_norm.make_inverse_operator(raw.info, fwd, cov, loose=0.2)

# compute inverse solution
stc = mne.minimum_norm.apply_inverse(evoked, inv, lambda2=1 / 3.0 ** 2, method='dSPM')

# morph it to average brain for group study
stc_avg = mne.morph_data('sample', 'fsaverage', stc, 5, smooth=5)