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_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces."""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    src_ = inverse_operator['src']
    with pytest.warns(RuntimeWarning, match='will have holes'):
        connectivity = spatio_temporal_src_connectivity(src_, n_times=2)
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert (a == b)

    assert_equal(grade_to_tris(5).shape, [40960, 3])

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 test_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces"""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    with warnings.catch_warnings(record=True) as w:
        connectivity = spatio_temporal_src_connectivity(
                                            inverse_operator['src'], n_times=2)
        assert len(w) == 1
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert_true(a == b)

def test_spatio_temporal_src_connectivity():
    """Test spatio-temporal connectivity from source spaces."""
    tris = np.array([[0, 1, 2], [3, 4, 5]])
    src = [dict(), dict()]
    connectivity = spatio_temporal_tris_connectivity(tris, 2)
    src[0]['use_tris'] = np.array([[0, 1, 2]])
    src[1]['use_tris'] = np.array([[0, 1, 2]])
    src[0]['vertno'] = np.array([0, 1, 2])
    src[1]['vertno'] = np.array([0, 1, 2])
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity2 = spatio_temporal_src_connectivity(src, 2)
    assert_array_equal(connectivity.todense(), connectivity2.todense())
    # add test for dist connectivity
    src[0]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[1]['dist'] = np.ones((3, 3)) - np.eye(3)
    src[0]['vertno'] = [0, 1, 2]
    src[1]['vertno'] = [0, 1, 2]
    src[0]['type'] = 'surf'
    src[1]['type'] = 'surf'
    connectivity3 = spatio_temporal_src_connectivity(src, 2, dist=2)
    assert_array_equal(connectivity.todense(), connectivity3.todense())
    # add test for source space connectivity with omitted vertices
    inverse_operator = read_inverse_operator(fname_inv)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        src_ = inverse_operator['src']
        connectivity = spatio_temporal_src_connectivity(src_, n_times=2)
    assert_equal(len(w), 1)
    a = connectivity.shape[0] / 2
    b = sum([s['nuse'] for s in inverse_operator['src']])
    assert_true(a == b)

    assert_equal(grade_to_tris(5).shape, [40960, 3])

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 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 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 test_plot_snr():
    """Test plotting SNR estimate."""
    import matplotlib.pyplot as plt
    inv = read_inverse_operator(inv_fname)
    evoked = read_evokeds(evoked_fname, baseline=(None, 0))[0]
    plot_snr_estimate(evoked, inv)
    plt.close('all')

def test_psf_ctf():
    """Test computation of PSFs and CTFs for linear estimators
    """

    inverse_operator = read_inverse_operator(fname_inv)
    forward = read_forward_solution(fname_fwd, force_fixed=False,
                                    surf_ori=True)
    forward = pick_types_forward(forward, meg=True, eeg=False)
    labels = [mne.read_label(ss) for ss in fname_label]

    method = 'MNE'
    n_svd_comp = 2

    # Test PSFs (then CTFs)
    for mode in ('sum', 'svd'):
        stc_psf, psf_ev = point_spread_function(inverse_operator,
                                                forward,
                                                method=method,
                                                labels=labels,
                                                lambda2=lambda2,
                                                pick_ori='normal',
                                                mode=mode,
                                                n_svd_comp=n_svd_comp)

        n_vert, n_samples = stc_psf.shape
        should_n_vert = (inverse_operator['src'][1]['vertno'].shape[0] +
                         inverse_operator['src'][0]['vertno'].shape[0])
        if mode == 'svd':
            should_n_samples = len(labels) * n_svd_comp + 1
        else:
            should_n_samples = len(labels) + 1

        assert_true(n_vert == should_n_vert)
        assert_true(n_samples == should_n_samples)

        n_chan, n_samples = psf_ev.data.shape
        assert_true(n_chan == forward['nchan'])

    forward = read_forward_solution(fname_fwd, force_fixed=True, surf_ori=True)
    forward = pick_types_forward(forward, meg=True, eeg=False)

    # Test CTFs
    for mode in ('sum', 'svd'):
        stc_ctf = cross_talk_function(inverse_operator, forward,
                                      labels, method=method,
                                      lambda2=lambda2,
                                      signed=False, mode=mode,
                                      n_svd_comp=n_svd_comp)

        n_vert, n_samples = stc_ctf.shape
        should_n_vert = (inverse_operator['src'][1]['vertno'].shape[0] +
                         inverse_operator['src'][0]['vertno'].shape[0])
        if mode == 'svd':
            should_n_samples = len(labels) * n_svd_comp + 1
        else:
            should_n_samples = len(labels) + 1

        assert_true(n_vert == should_n_vert)
        assert_true(n_samples == should_n_samples)

def plot_inverse_operator(subject, fnout_img=None, verbose=None):
    """"Reads in and plots the inverse solution."""

    # get name of inverse solution-file
    subjects_dir = os.environ.get('SUBJECTS_DIR')
    fname_dir_inv = os.path.join(subjects_dir,
                                 subject)

    for files in os.listdir(fname_dir_inv):
        if files.endswith(",cleaned_epochs_avg-7-src-meg-inv.fif"):
            fname_inv = os.path.join(fname_dir_inv,
                                     files)
            break

    try:
        fname_inv
    except NameError:
        print "ERROR: No forward solution found!"
        sys.exit()


    # read inverse solution
    inv = min_norm.read_inverse_operator(fname_inv)

    # print some information if desired
    if verbose is not None:
        print "Method: %s" % inv['methods']
        print "fMRI prior: %s" % inv['fmri_prior']
        print "Number of sources: %s" % inv['nsource']
        print "Number of channels: %s" % inv['nchan']


    # show 3D source space
    lh_points = inv['src'][0]['rr']
    lh_faces = inv['src'][0]['tris']
    rh_points = inv['src'][1]['rr']
    rh_faces = inv['src'][1]['tris']

    # create figure and plot results
    fig_inverse = mlab.figure(size=(1200, 1200),
                              bgcolor=(0, 0, 0))

    mlab.triangular_mesh(lh_points[:, 0],
                         lh_points[:, 1],
                         lh_points[:, 2],
                         lh_faces)

    mlab.triangular_mesh(rh_points[:, 0],
                         rh_points[:, 1],
                         rh_points[:, 2],
                         rh_faces)

    # save result
    if fnout_img is not None:
        mlab.savefig(fnout_img,
                     figure=fig_inverse,
                     size=(1200, 1200))

    mlab.close(all=True)

def calc_inv_kernel(fn_inv, method="dSPM", nave=1, snr=6.,
                    pick_ori="normal", verbose=None):

    """
    Interface for preparing the kernel of the inverse
    estimation.

        Parameters
        ----------
        fn_inv : String containing the filename
            of the inverse operator (must be a fif-file)
        method : string
            which source localization method should be used?
            MNE-based: "MNE" | "dSPM" | "sLORETA"
            default: method="dSPM"
        nave : number of averages used to regularize the solution
            default: nave=1
        snr : signal-to-noise ratio
            default: snr = 3.
        verbose : bool, str, int, or None
            If not None, override default verbose level
            (see mne.verbose).
            default: verbose=None
    """

    # -------------------------------------------
    # import necessary modules
    # -------------------------------------------
    import mne.minimum_norm as min_norm
    from mne.minimum_norm.inverse import _assemble_kernel
    import numpy as np

    # -------------------------------------------
    # estimate inverse kernel
    # -------------------------------------------
    # load inverse solution
    inv_operator = min_norm.read_inverse_operator(fn_inv, verbose=verbose)

    # set up the inverse according to the parameters
    lambda2      = 1. / snr ** 2.   # the regularization parameter.
    inv_operator = min_norm.prepare_inverse_operator(inv_operator, nave, lambda2, method)

    # estimate inverse kernel and noise normalization coefficient
    kernel, noise_norm, vertno = _assemble_kernel(inv_operator, None, method, pick_ori)

    if method == "MNE":
        noise_norm = np.ones((kernel.shape[0]/3))
        noise_norm = noise_norm[:, np.newaxis]


    # -------------------------------------------
    # return results
    # -------------------------------------------
    return kernel, noise_norm, vertno

def test_surface_vector_source_morph():
    """Test surface and vector source estimate morph."""
    tempdir = _TempDir()

    inverse_operator_surf = read_inverse_operator(fname_inv_surf)

    stc_surf = read_source_estimate(fname_smorph, subject='sample')
    stc_surf.crop(0.09, 0.1)  # for faster computation

    stc_vec = _real_vec_stc()

    source_morph_surf = compute_source_morph(
        inverse_operator_surf['src'], subjects_dir=subjects_dir,
        smooth=1, warn=False)  # smooth 1 for speed
    assert source_morph_surf.subject_from == 'sample'
    assert source_morph_surf.subject_to == 'fsaverage'
    assert source_morph_surf.kind == 'surface'
    assert isinstance(source_morph_surf.src_data, dict)
    assert isinstance(source_morph_surf.src_data['vertices_from'], list)
    assert isinstance(source_morph_surf, SourceMorph)
    stc_surf_morphed = source_morph_surf.apply(stc_surf)
    assert isinstance(stc_surf_morphed, SourceEstimate)
    stc_vec_morphed = source_morph_surf.apply(stc_vec)
    with pytest.raises(ValueError, match='Only volume source estimates'):
        source_morph_surf.apply(stc_surf, output='nifti1')

    # check if correct class after morphing
    assert isinstance(stc_surf_morphed, SourceEstimate)
    assert isinstance(stc_vec_morphed, VectorSourceEstimate)

    # check __repr__
    assert 'surface' in repr(source_morph_surf)

    # check loading and saving for surf
    source_morph_surf.save(op.join(tempdir, '42.h5'))

    source_morph_surf_r = read_source_morph(op.join(tempdir, '42.h5'))

    assert (all([read == saved for read, saved in
                 zip(sorted(source_morph_surf_r.__dict__),
                     sorted(source_morph_surf.__dict__))]))

    # check wrong subject correction
    stc_surf.subject = None
    assert isinstance(source_morph_surf.apply(stc_surf), SourceEstimate)

    # degenerate
    stc_vol = read_source_estimate(fname_vol, 'sample')
    with pytest.raises(ValueError, match='stc_from was type'):
        source_morph_surf.apply(stc_vol)

def test_snr():
    """Test SNR calculation"""
    tempdir = _TempDir()
    inv = read_inverse_operator(fname_inv)
    evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0]
    snr = estimate_snr(evoked, inv)[0]
    orig_dir = os.getcwd()
    os.chdir(tempdir)
    try:
        cmd = ['mne_compute_mne', '--inv', fname_inv, '--meas', fname_evoked,
               '--snronly', '--bmin', '-200', '--bmax', '0']
        run_subprocess(cmd)
    except Exception:
        pass  # this returns 1 for some reason
    finally:
        os.chdir(orig_dir)
    snr_c = np.loadtxt(op.join(tempdir, 'SNR'))[:, 1]
    assert_allclose(snr, snr_c, atol=1e-2, rtol=1e-2)

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 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 load(typ, subject='fsaverage', analysis='analysis', block=999,
         download=True, preload=False):
    """Auxiliary saving function."""
    # get file name
    fname = paths(typ, subject=subject, analysis=analysis, block=block)

    # check if file exists
    if not op.exists(fname) and download:
        client.download(fname)

    # different data format depending file type
    if typ == 'behavior':
        from base import read_events
        out = read_events(fname)
    elif typ == 'sss':
        out = Raw(fname, preload=preload)
    elif typ in ['epo_block', 'epochs', 'epochs_decim', 'epochs_vhp']:
        out = read_epochs(fname, preload=preload)
    elif typ in ['cov']:
        from mne.cov import read_cov
        out = read_cov(fname)
    elif typ in ['fwd']:
        from mne import read_forward_solution
        out = read_forward_solution(fname, surf_ori=True)
    elif typ in ['inv']:
        from mne.minimum_norm import read_inverse_operator
        out = read_inverse_operator(fname)
    elif typ in ['evoked', 'decod', 'decod_tfr', 'score', 'score_tfr',
                 'evoked_source']:
        with open(fname, 'rb') as f:
            out = pickle.load(f)
    elif typ == 'morph':
        from scipy.sparse import csr_matrix
        loader = np.load(fname)
        out = csr_matrix((loader['data'], loader['indices'], loader['indptr']),
                         shape=loader['shape'])
    elif typ in ['score_source', 'score_pval']:
        out = np.load(fname)
    else:
        raise NotImplementedError()
    return out

def _calc_inverse(params):
    subject, epochs, overwrite = params
    epo = op.join(REMOTE_ROOT_DIR, 'ave', '{}_ecr_nTSSS_conflict-epo.fif'.format(subject))
    fwd = op.join(REMOTE_ROOT_DIR, 'fwd', '{}_ecr-fwd.fif'.format(subject))
    local_inv_file_name = op.join(LOCAL_ROOT_DIR, 'inv', '{}_ecr_nTSSS_conflict-inv.fif'.format(subject))

    if os.path.isfile(local_inv_file_name) and not overwrite:
        inverse_operator = read_inverse_operator(local_inv_file_name)
        print('inv already calculated for {}'.format(subject))
    else:
        if epochs is None:
            epochs = mne.read_epochs(epo)
        noise_cov = mne.compute_covariance(epochs.crop(None, 0, copy=True))
        inverse_operator = None
        if not os.path.isfile(fwd):
            print('no fwd for {}'.format(subject))
        else:
            forward = mne.read_forward_solution(fwd)
            inverse_operator = make_inverse_operator(epochs.info, forward, noise_cov,
                loose=None, depth=None)
            write_inverse_operator(local_inv_file_name, inverse_operator)
    return inverse_operator

def test_stc_as_volume():
    """Test previous volume source estimate morph."""
    import nibabel as nib
    inverse_operator_vol = read_inverse_operator(fname_inv_vol)

    # Apply inverse operator
    stc_vol = read_source_estimate(fname_vol, 'sample')

    img = stc_vol.as_volume(inverse_operator_vol['src'], mri_resolution=True,
                            dest='42')
    t1_img = nib.load(fname_t1)
    # always assure nifti and dimensionality
    assert isinstance(img, nib.Nifti1Image)
    assert img.header.get_zooms()[:3] == t1_img.header.get_zooms()[:3]

    img = stc_vol.as_volume(inverse_operator_vol['src'], mri_resolution=False)

    assert isinstance(img, nib.Nifti1Image)
    assert img.shape[:3] == inverse_operator_vol['src'][0]['shape'][:3]

    with pytest.raises(ValueError, match='invalid output'):
        stc_vol.as_volume(inverse_operator_vol['src'], format='42')

from my_settings import (epochs_folder, tf_folder, subjects_dir, mne_folder)
import mne
import sys
import numpy as np
import pandas as pd

from mne.minimum_norm import read_inverse_operator, source_induced_power

subject = sys.argv[1]

epochs = mne.read_epochs(epochs_folder + "%s_target-epo.fif" % subject)

inv = read_inverse_operator(mne_folder + "%s-inv.fif" % subject)

labels = mne.read_labels_from_annot(
    subject,
    parc='PALS_B12_Lobes',
    # regexp="Bro",
    subjects_dir=subjects_dir)
labels_selc = labels[9], labels[10]

frequencies = np.arange(8, 13, 1)  # define frequencies of interest
n_cycles = frequencies / 3.  # different number of cycle per frequency
method = "dSPM"

sides = ["left", "right"]
conditions = ["ctl", "ent"]
cor = ["correct", "incorrect"]
phase = ["in_phase", "out_phase"]
congrunet = ["cong", "incong"]