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

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

def test_volume_labels_morph(tmpdir):
    """Test generating a source space from volume label."""
    # see gh-5224
    evoked = mne.read_evokeds(fname_evoked)[0].crop(0, 0)
    evoked.pick_channels(evoked.ch_names[:306:8])
    evoked.info.normalize_proj()
    n_ch = len(evoked.ch_names)
    aseg_fname = op.join(subjects_dir, 'sample', 'mri', 'aseg.mgz')
    label_names = get_volume_labels_from_aseg(aseg_fname)
    src = setup_volume_source_space(
        'sample', subjects_dir=subjects_dir, volume_label=label_names[:2],
        mri=aseg_fname)
    assert len(src) == 2
    assert src.kind == 'volume'
    n_src = sum(s['nuse'] for s in src)
    sphere = make_sphere_model('auto', 'auto', evoked.info)
    fwd = make_forward_solution(evoked.info, fname_trans, src, sphere)
    assert fwd['sol']['data'].shape == (n_ch, n_src * 3)
    inv = make_inverse_operator(evoked.info, fwd, make_ad_hoc_cov(evoked.info),
                                loose=1.)
    stc = apply_inverse(evoked, inv)
    assert stc.data.shape == (n_src, 1)
    img = stc.as_volume(src, mri_resolution=True)
    n_on = np.array(img.dataobj).astype(bool).sum()
    assert n_on == 291  # was 291 on `master` before gh-5590
    img = stc.as_volume(src, mri_resolution=False)
    n_on = np.array(img.dataobj).astype(bool).sum()
    assert n_on == 44  # was 20 on `master` before gh-5590

def run_inverse(subject_id):
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)
    data_path = op.join(meg_dir, subject)

    fname_ave = op.join(data_path, '%s-ave.fif' % subject)
    fname_cov = op.join(data_path, '%s-cov.fif' % subject)
    fname_fwd = op.join(data_path, '%s-meg-%s-fwd.fif' % (subject, spacing))
    fname_inv = op.join(data_path, '%s-meg-%s-inv.fif' % (subject, spacing))

    evokeds = mne.read_evokeds(fname_ave, condition=[0, 1, 2, 3, 4, 5])
    cov = mne.read_cov(fname_cov)
    # cov = mne.cov.regularize(cov, evokeds[0].info,
    #                                mag=0.05, grad=0.05, eeg=0.1, proj=True)

    forward = mne.read_forward_solution(fname_fwd, surf_ori=True)
    # forward = mne.pick_types_forward(forward, meg=True, eeg=False)

    # make an M/EEG, MEG-only, and EEG-only inverse operators
    info = evokeds[0].info
    inverse_operator = make_inverse_operator(info, forward, cov,
                                             loose=0.2, depth=0.8)

    write_inverse_operator(fname_inv, inverse_operator)

    # Compute inverse solution
    snr = 3.0
    lambda2 = 1.0 / snr ** 2

    for evoked in evokeds:
        stc = apply_inverse(evoked, inverse_operator, lambda2, "dSPM",
                            pick_ori=None)

        stc.save(op.join(data_path, 'mne_dSPM_inverse-%s' % evoked.comment))

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

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

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

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

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

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

    assert_equal(len(epochs), 2)

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

    inv = read_inverse_operator(fname_inv)

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

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

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

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

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

def INVERSE(wdir, Subject, epoch_info, evokeds):

    # import parameters from configuration file
    from configuration import ( lambda2, method )

    # compute noise covariance from empty room data
    emptyroom_raw = mne.io.Raw(wdir + '/data/maxfilter/' + Subject + '/'+ Subject +'_empty_sss.fif')  
    noise_cov     = mne.compute_raw_data_covariance(emptyroom_raw)
    
    # compute dSPM solution
    fname_fwd     = wdir + '/data/forward/' + Subject + '/' + Subject + '_phase1_trans_sss_filt140_raw-ico5-fwd.fif'
    forward       = mne.read_forward_solution(fname_fwd, surf_ori=True)
    
    # create inverse operator
    inverse_operator = make_inverse_operator(epoch_info, forward, noise_cov, loose=0.4, depth=0.8)
    
    # Compute inverse solution
    stcs = []
    for evoked in evokeds:
        stcs.append(apply_inverse(evoked, inverse_operator, lambda2, method=method, pick_ori = None))
    
    # save a covariance picture for visual inspection
    mne.viz.plot_cov(noise_cov, epoch_info, colorbar=True, proj=True,show_svd=False,show=False)
    plt.savefig(wdir + "/plots/" + Subject + "_covmat")
    plt.close()
    
    return stcs

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 get_mne_stc(ndvar=False, vol=False, subject='sample'):
    """MNE-Python SourceEstimate

    Parameters
    ----------
    ndvar : bool
        Convert to NDVar (default False; src="ico-4" is false, but it works as
        long as the source space is not accessed).
    vol : bool
        Volume source estimate.
    """
    data_path = Path(mne.datasets.testing.data_path())
    meg_sdir = data_path / 'MEG/sample'
    subjects_dir = data_path / 'subjects'
    # scaled subject
    if subject == 'fsaverage_scaled':
        subject_dir = os.path.join(subjects_dir, subject)
        if not os.path.exists(subject_dir):
            mne.scale_mri('fsaverage', subject, .9, subjects_dir=subjects_dir, skip_fiducials=True, labels=False, annot=True)
        data_subject = 'fsaverage'
    else:
        data_subject = subject

    if vol:
        inv = mn.read_inverse_operator(str(meg_sdir / 'sample_audvis_trunc-meg-vol-7-meg-inv.fif'))
        evoked = mne.read_evokeds(str(meg_sdir / 'sample_audvis_trunc-ave.fif'), 'Left Auditory')
        stc = mn.apply_inverse(evoked, inv, method='MNE', pick_ori='vector')
        if data_subject == 'fsaverage':
            m = mne.compute_source_morph(stc, 'sample', data_subject, subjects_dir)
            stc = m.apply(stc)
            stc.subject = subject
        elif subject != 'sample':
            raise ValueError(f"subject={subject!r}")
        if ndvar:
            return load.fiff.stc_ndvar(stc, subject, 'vol-7', subjects_dir, 'MNE', sss_filename='{subject}-volume-7mm-src.fif')
        else:
            return stc
    stc_path = meg_sdir / f'{data_subject}_audvis_trunc-meg'
    if ndvar:
        return load.fiff.stc_ndvar(stc_path, subject, 'ico-5', subjects_dir)
    else:
        return mne.read_source_estimate(str(stc_path), subject)

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

def inverse_function(sub_id, session):
    """ Will calculate the inverse model based dSPM
    """
    data_path = "/media/mje/My_Book/Data/MEG/MEG_libet/sub_2_tests"
    fname = "sub_%d_%s_tsss_mc" % (sub_id, session)
    fname_epochs = data_path + fname + "_epochs.fif"
    fname_fwd_meg = data_path + fname + "_fwd.fif"
    fname_cov = data_path + fname + "_cov.fif"
    fname_inv = data_path + fname + "_inv.fif"
    fname_stcs = fname + "_mne_dSPM_inverse"

    epochs = mne.read_epochs(fname_epochs)
    evoked = epochs.average()

    snr = 3.0
    lambda2 = 1.0 / snr ** 2

    # Load data
    forward_meg = mne.read_forward_solution(fname_fwd_meg, surf_ori=True)
    noise_cov = mne.read_cov(fname_cov)

    # regularize noise covariance
    noise_cov = mne.cov.regularize(noise_cov, evoked.info,
                                   mag=0.05, grad=0.05, eeg=0.1, proj=True)

    # Restrict forward solution as necessary for MEG
    forward_meg = mne.fiff.pick_types_forward(forward_meg, meg=True, eeg=False)

    # make an M/EEG, MEG-only, and EEG-only inverse operators
    info = evoked.info
    inverse_operator_meg = make_inverse_operator(info, forward_meg, noise_cov,
                                                 loose=0.2, depth=0.8)

    write_inverse_operator(fname_inv, inverse_operator_meg)

    # Compute inverse solution
    stc = apply_inverse(evoked, inverse_operator_meg, lambda2, "dSPM",
                        pick_normal=False)

    # Save result in stc files
    stc.save(fname_stcs)

evoked = read_evokeds(fname_evoked, condition=0, baseline=(None, 0))
evoked.pick_types(meg=True, eeg=False)

###############################################################################
# Then, we can load the precomputed inverse operator from a file.
fname_inv = data_path + '/MEG/sample/sample_audvis-meg-vol-7-meg-inv.fif'
inv = read_inverse_operator(fname_inv)
src = inv['src']

###############################################################################
# The source estimate is computed using the inverse operator and the
# sensor-space data.
snr = 3.0
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
stc = apply_inverse(evoked, inv, lambda2, method)
stc.crop(0.0, 0.2)

###############################################################################
# This time, we have a different container
# (:class:`VolSourceEstimate <mne.VolSourceEstimate>`) for the source time
# course.
print(stc)

###############################################################################
# This too comes with a convenient plot method.

stc.plot(src, subject='sample', subjects_dir=subjects_dir)

###############################################################################
# For this visualization, ``nilearn`` must be installed.

# Show the dipoles as arrows pointing along the surface normal
normals = lh['nn'][lh['vertno']]
mlab.quiver3d(dip_pos[:, 0], dip_pos[:, 1], dip_pos[:, 2],
              normals[:, 0], normals[:, 1], normals[:, 2],
              color=red, scale_factor=1E-3)

mlab.view(azimuth=180, distance=0.1)

###############################################################################
# Restricting the dipole orientations in this manner leads to the following
# source estimate for the sample data:

# Compute the source estimate for the 'left - auditory' condition in the sample
# dataset.
inv = make_inverse_operator(left_auditory.info, fwd, noise_cov, fixed=True)
stc = apply_inverse(left_auditory, inv, pick_ori=None)

# Visualize it at the moment of peak activity.
_, time_max = stc.get_peak(hemi='lh')
brain = stc.plot(surface='white', subjects_dir=subjects_dir,
                 initial_time=time_max, time_unit='s', size=(600, 400))

###############################################################################
# The direction of the estimated current is now restricted to two directions:
# inward and outward. In the plot, blue areas indicate current flowing inwards
# and red areas indicate current flowing outwards. Given the curvature of the
# cortex, groups of dipoles tend to point in the same direction: the direction
# of the electromagnetic field picked up by the sensors.

###############################################################################
# Loose dipole orientations

def test_mxne_inverse():
    """Test (TF-)MxNE inverse computation."""
    # Read noise covariance matrix
    cov = read_cov(fname_cov)

    # Handling average file
    loose = 0.0
    depth = 0.9

    evoked = read_evokeds(fname_data, condition=0, baseline=(None, 0))
    evoked.crop(tmin=-0.05, tmax=0.2)

    evoked_l21 = evoked.copy()
    evoked_l21.crop(tmin=0.081, tmax=0.1)
    label = read_label(fname_label)

    forward = read_forward_solution(fname_fwd)
    forward = convert_forward_solution(forward, surf_ori=True)

    # Reduce source space to make test computation faster
    inverse_operator = make_inverse_operator(evoked_l21.info, forward, cov,
                                             loose=loose, depth=depth,
                                             fixed=True, use_cps=True)
    stc_dspm = apply_inverse(evoked_l21, inverse_operator, lambda2=1. / 9.,
                             method='dSPM')
    stc_dspm.data[np.abs(stc_dspm.data) < 12] = 0.0
    stc_dspm.data[np.abs(stc_dspm.data) >= 12] = 1.
    weights_min = 0.5

    # MxNE tests
    alpha = 70  # spatial regularization parameter

    stc_prox = mixed_norm(evoked_l21, forward, cov, alpha, loose=loose,
                          depth=depth, maxit=300, tol=1e-8,
                          active_set_size=10, weights=stc_dspm,
                          weights_min=weights_min, solver='prox')
    with pytest.warns(None):  # CD
        stc_cd = mixed_norm(evoked_l21, forward, cov, alpha, loose=loose,
                            depth=depth, maxit=300, tol=1e-8,
                            active_set_size=10, weights=stc_dspm,
                            weights_min=weights_min, solver='cd',
                            pca=False)  # pca=False deprecated, doesn't matter
    stc_bcd = mixed_norm(evoked_l21, forward, cov, alpha, loose=loose,
                         depth=depth, maxit=300, tol=1e-8, active_set_size=10,
                         weights=stc_dspm, weights_min=weights_min,
                         solver='bcd')
    assert_array_almost_equal(stc_prox.times, evoked_l21.times, 5)
    assert_array_almost_equal(stc_cd.times, evoked_l21.times, 5)
    assert_array_almost_equal(stc_bcd.times, evoked_l21.times, 5)
    assert_allclose(stc_prox.data, stc_cd.data, rtol=1e-3, atol=0.0)
    assert_allclose(stc_prox.data, stc_bcd.data, rtol=1e-3, atol=0.0)
    assert_allclose(stc_cd.data, stc_bcd.data, rtol=1e-3, atol=0.0)
    assert stc_prox.vertices[1][0] in label.vertices
    assert stc_cd.vertices[1][0] in label.vertices
    assert stc_bcd.vertices[1][0] in label.vertices

    with pytest.warns(None):  # CD
        dips = mixed_norm(evoked_l21, forward, cov, alpha, loose=loose,
                          depth=depth, maxit=300, tol=1e-8, active_set_size=10,
                          weights=stc_dspm, weights_min=weights_min,
                          solver='cd', return_as_dipoles=True)
    stc_dip = make_stc_from_dipoles(dips, forward['src'])
    assert isinstance(dips[0], Dipole)
    assert stc_dip.subject == "sample"
    _check_stcs(stc_cd, stc_dip)

    with pytest.warns(None):  # CD
        stc, _ = mixed_norm(evoked_l21, forward, cov, alpha, loose=loose,
                            depth=depth, maxit=300, tol=1e-8,
                            active_set_size=10, return_residual=True,
                            solver='cd')
    assert_array_almost_equal(stc.times, evoked_l21.times, 5)
    assert stc.vertices[1][0] in label.vertices

    # irMxNE tests
    with pytest.warns(None):  # CD
        stc = mixed_norm(evoked_l21, forward, cov, alpha,
                         n_mxne_iter=5, loose=loose, depth=depth,
                         maxit=300, tol=1e-8, active_set_size=10,
                         solver='cd')
    assert_array_almost_equal(stc.times, evoked_l21.times, 5)
    assert stc.vertices[1][0] in label.vertices
    assert stc.vertices == [[63152], [79017]]

    # Do with TF-MxNE for test memory savings
    alpha = 60.  # overall regularization parameter
    l1_ratio = 0.01  # temporal regularization proportion

    stc, _ = tf_mixed_norm(evoked, forward, cov,
                           loose=loose, depth=depth, maxit=100, tol=1e-4,
                           tstep=4, wsize=16, window=0.1, weights=stc_dspm,
                           weights_min=weights_min, return_residual=True,
                           alpha=alpha, l1_ratio=l1_ratio)
    assert_array_almost_equal(stc.times, evoked.times, 5)
    assert stc.vertices[1][0] in label.vertices

    pytest.raises(ValueError, tf_mixed_norm, evoked, forward, cov,
                  alpha=101, l1_ratio=0.03)
    pytest.raises(ValueError, tf_mixed_norm, evoked, forward, cov,
                  alpha=50., l1_ratio=1.01)

evoked2_fname = data_path + 'ave_projon/9367_s5_Noun_Place_All-ave.fif'

snr = 3.0
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
inverse_operator = read_inverse_operator(fname_inv)
sample_vertices = [s['vertno'] for s in inverse_operator['src']]

#    Let's average and compute inverse, resampling to speed things up
#evoked1 = epochs1.average()
evoked1 = mne.read_evokeds(evoked1_fname, condition = 'epochs_TaggedWord')
print evoked1
print 
print inverse_operator
#evoked1.resample(50)
condition1 = apply_inverse(evoked1, inverse_operator, lambda2, method)

#evoked2 = epochs2.average()
evoked2 = mne.read_evokeds(evoked2_fname, condition = 'epochs_TaggedWord')
print evoked1
print 
print inverse_operator
#evoked2.resample(50)
condition2 = apply_inverse(evoked2, inverse_operator, lambda2, method)

#    Let's only deal with t > 0, cropping to reduce multiple comparisons
condition1.crop(0, None)
condition2.crop(0, None)
tmin = condition1.tmin
tstep = condition1.tstep
#

def test_mxne_inverse():
    """Test (TF-)MxNE inverse computation"""
    # Handling forward solution
    evoked = fiff.Evoked(fname_data, setno=1, baseline=(None, 0))

    # Read noise covariance matrix
    cov = read_cov(fname_cov)

    # Handling average file
    setno = 0
    loose = None
    depth = 0.9

    evoked = fiff.read_evoked(fname_data, setno=setno, baseline=(None, 0))
    evoked.crop(tmin=-0.1, tmax=0.4)

    evoked_l21 = copy.deepcopy(evoked)
    evoked_l21.crop(tmin=0.08, tmax=0.1)
    label = read_label(fname_label)
    weights_min = 0.5
    forward = read_forward_solution(fname_fwd, force_fixed=False,
                                    surf_ori=True)

    # Reduce source space to make test computation faster
    inverse_operator = make_inverse_operator(evoked.info, forward, cov,
                                             loose=loose, depth=depth,
                                             fixed=True)
    stc_dspm = apply_inverse(evoked_l21, inverse_operator, lambda2=1. / 9.,
                             method='dSPM')
    stc_dspm.data[np.abs(stc_dspm.data) < 12] = 0.0
    stc_dspm.data[np.abs(stc_dspm.data) >= 12] = 1.

    # MxNE tests
    alpha = 60  # spatial regularization parameter

    stc_prox = mixed_norm(evoked_l21, forward, cov, alpha, loose=None,
                          depth=0.9, maxit=1000, tol=1e-8, active_set_size=10,
                          solver='prox')
    stc_cd = mixed_norm(evoked_l21, forward, cov, alpha, loose=None,
                        depth=0.9, maxit=1000, tol=1e-8, active_set_size=10,
                        solver='cd')
    assert_array_almost_equal(stc_prox.times, evoked_l21.times, 5)
    assert_array_almost_equal(stc_cd.times, evoked_l21.times, 5)
    assert_array_almost_equal(stc_prox.data, stc_cd.data, 5)
    assert_true(stc_prox.vertno[1][0] in label.vertices)
    assert_true(stc_cd.vertno[1][0] in label.vertices)

    stc, _ = mixed_norm(evoked_l21, forward, cov, alpha, loose=None,
                        depth=depth, maxit=500, tol=1e-4, active_set_size=10,
                        weights=stc_dspm, weights_min=weights_min,
                        return_residual=True)

    assert_array_almost_equal(stc.times, evoked_l21.times, 5)
    assert_true(stc.vertno[1][0] in label.vertices)

    # Do with TF-MxNE for test memory savings
    alpha_space = 60.  # spatial regularization parameter
    alpha_time = 1.  # temporal regularization parameter

    stc, _ = tf_mixed_norm(evoked, forward, cov, alpha_space, alpha_time,
                           loose=loose, depth=depth, maxit=100, tol=1e-4,
                           tstep=4, wsize=16, window=0.1, weights=stc_dspm,
                           weights_min=weights_min, return_residual=True)

    assert_array_almost_equal(stc.times, evoked.times, 5)
    assert_true(stc.vertno[1][0] in label.vertices)

brain_gen = stc_gen.plot(clim=clim, figure=figs, **kwargs)

###############################################################################
# Simulate sensor-space signals
# -----------------------------
#
# Use the forward solution and add Gaussian noise to simulate sensor-space
# (evoked) data from the known source-space signals. The amount of noise is
# controlled by `nave` (higher values imply less noise).
#
evoked_gen = simulate_evoked(fwd, stc_gen, evoked.info, cov, nave,
                             random_state=seed)

# Map the simulated sensor-space data to source-space using the inverse
# operator.
stc_inv = apply_inverse(evoked_gen, inv_op, lambda2, method=method)

###############################################################################
# Plot the point-spread of corrupted signal
# -----------------------------------------
#
# Notice that after applying the forward- and inverse-operators to the known
# point sources that the point sources have spread across the source-space.
# This spread is due to the minimum norm solution so that the signal leaks to
# nearby vertices with similar orientations so that signal ends up crossing the
# sulci and gyri.
figs = [mlab.figure(5), mlab.figure(6), mlab.figure(7), mlab.figure(8)]
brain_inv = stc_inv.plot(figure=figs, **kwargs)

###############################################################################
# Exercises

inv = mne.minimum_norm.make_inverse_operator(evoked_std.info, fwd, cov)
snr = 3.0
lambda2 = 1.0 / snr ** 2
del fwd

###############################################################################
# The sources are computed using dSPM method and plotted on an inflated brain
# surface. For interactive controls over the image, use keyword
# ``time_viewer=True``.
# Standard condition.
stc_standard = mne.minimum_norm.apply_inverse(evoked_std, inv, lambda2, 'dSPM')
brain = stc_standard.plot(subjects_dir=subjects_dir, subject=subject,
                          surface='inflated', time_viewer=False, hemi='lh',
                          initial_time=0.1, time_unit='s')
del stc_standard, brain

###############################################################################
# Deviant condition.
stc_deviant = mne.minimum_norm.apply_inverse(evoked_dev, inv, lambda2, 'dSPM')
brain = stc_deviant.plot(subjects_dir=subjects_dir, subject=subject,
                         surface='inflated', time_viewer=False, hemi='lh',
                         initial_time=0.1, time_unit='s')
del stc_deviant, brain

###############################################################################
# Difference.
stc_difference = apply_inverse(evoked_difference, inv, lambda2, 'dSPM')
brain = stc_difference.plot(subjects_dir=subjects_dir, subject=subject,
                            surface='inflated', time_viewer=False, hemi='lh',
                            initial_time=0.15, time_unit='s')

                                              loose=0.2, depth=0.8)
inverse_operator_meg = make_inverse_operator(info, forward_meg, noise_cov,
                                              loose=0.2, depth=0.8)
inverse_operator_eeg = make_inverse_operator(info, forward_eeg, noise_cov,
                                              loose=0.2, depth=0.8)

write_inverse_operator('sample_audvis-meeg-oct-6-inv.fif',
                       inverse_operator_meeg)
write_inverse_operator('sample_audvis-meg-oct-6-inv.fif',
                       inverse_operator_meg)
write_inverse_operator('sample_audvis-eeg-oct-6-inv.fif',
                       inverse_operator_eeg)

# Compute inverse solution
stcs = dict()
stcs['meeg'] = apply_inverse(evoked, inverse_operator_meeg, lambda2, "dSPM",
                        pick_normal=False)
stcs['meg'] = apply_inverse(evoked, inverse_operator_meg, lambda2, "dSPM",
                        pick_normal=False)
stcs['eeg'] = apply_inverse(evoked, inverse_operator_eeg, lambda2, "dSPM",
                        pick_normal=False)

# Save result in stc files
names = ['meeg', 'meg', 'eeg']
for name in names:
    stcs[name].save('mne_dSPM_inverse-%s' % name)

###############################################################################
# View activation time-series
pl.close('all')
pl.figure(figsize=(8, 6))
for ii in range(len(stcs)):

def apply_inverse(fnepo, method='dSPM', event='LLst', min_subject='fsaverage', STC_US='ROI', 
                  snr=5.0):
    '''  
        Parameter
        ---------
        fnepo: string or list
            The epochs 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.
        STC_US: string
            The using of the inversion for further analysis.
            'ROI' stands for ROIs definition, 'CAU' stands for causality analysis.
        snr: signal to noise ratio for inverse solution. 
    '''
    #Get the default subjects_dir
    from mne.minimum_norm import (apply_inverse, apply_inverse_epochs)
    subjects_dir = os.environ['SUBJECTS_DIR']
    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)
        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 = fn_path + '/%s-trans.fif' % subject
        fn_cov = fn_path + '/%s_empty,nr-cov.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
        snr = snr
        lambda2 = 1.0 / snr ** 2 
        #noise_cov = mne.read_cov(fn_cov)
        epochs = mne.read_epochs(fname)
        noise_cov = mne.read_cov(fn_cov)
        if STC_US == 'ROI':
            # 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
            evoked = epochs.average()
            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,