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 run_evoked(subject_id):
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)

    data_path = op.join(meg_dir, subject)
    epochs = mne.read_epochs(op.join(data_path, '%s-epo.fif' % subject),
                             preload=False)

    evoked_famous = epochs['face/famous'].average()
    evoked_scrambled = epochs['scrambled'].average()
    evoked_unfamiliar = epochs['face/unfamiliar'].average()

    # Simplify comment
    evoked_famous.comment = 'famous'
    evoked_scrambled.comment = 'scrambled'
    evoked_unfamiliar.comment = 'unfamiliar'

    contrast = mne.combine_evoked([evoked_famous, evoked_unfamiliar,
                                   evoked_scrambled], weights=[0.5, 0.5, -1.])
    contrast.comment = 'contrast'
    faces = mne.combine_evoked([evoked_famous, evoked_unfamiliar], 'nave')
    faces.comment = 'faces'

    mne.evoked.write_evokeds(op.join(data_path, '%s-ave.fif' % subject),
                             [evoked_famous, evoked_scrambled,
                              evoked_unfamiliar, contrast, faces])

    # take care of noise cov
    cov = mne.compute_covariance(epochs, tmax=0, method='shrunk')
    cov.save(op.join(data_path, '%s-cov.fif' % subject))

def compute_and_save_psd(epochs_fname, fmin=0, fmax=120,
                         method='welch', n_fft=256, n_overlap=0, 
                         picks=None, proj=False, n_jobs=1, verbose=None):
    """
    Load epochs from file,
    compute psd and save the result in numpy arrays
    """
    import numpy as np
    import os
    from mne import read_epochs
    epochs = read_epochs(epochs_fname)
    epochs_meg = epochs.pick_types(meg=True, eeg=False, eog=False, ecg=False)
    if method == 'welch':
        from mne.time_frequency import psd_welch
        psds, freqs = psd_welch(epochs_meg)
    elif method == 'multitaper':
        from mne.time_frequency import psd_multitaper
        psds, freqs = psd_multitaper(epochs_meg)
    else:
        raise Exception('nonexistent method for psd computation')
    path, name = os.path.split(epochs_fname)
    base, ext = os.path.splitext(name)
    psds_fname = base + '-psds.npz'
    # freqs_fname = base + '-freqs.npy'
    psds_fname = os.path.abspath(psds_fname)
    # print(psds.shape)
    np.savez(psds_fname, psds=psds, freqs=freqs)
    # np.save(freqs_file, freqs)
    return psds_fname

def test_access_by_name():
    """Test accessing epochs by event name and on_missing for rare events
    """
    assert_raises(ValueError, Epochs, raw, events, {1: 42, 2: 42}, tmin,
                  tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
                  tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
                  tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, 'foo', tmin, tmax,
                  picks=picks)
    # Test accessing non-existent events (assumes 12345678 does not exist)
    event_id_illegal = dict(aud_l=1, does_not_exist=12345678)
    assert_raises(ValueError, Epochs, raw, events, event_id_illegal,
                  tmin, tmax)
    # Test on_missing
    assert_raises(ValueError, Epochs, raw, events, 1, tmin, tmax,
                  on_missing='foo')
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        Epochs(raw, events, event_id_illegal, tmin, tmax, on_missing='warning')
        nw = len(w)
        assert_true(1 <= nw <= 2)
        Epochs(raw, events, event_id_illegal, tmin, tmax, on_missing='ignore')
        assert_equal(len(w), nw)
    epochs = Epochs(raw, events, {'a': 1, 'b': 2}, tmin, tmax, picks=picks)
    assert_raises(KeyError, epochs.__getitem__, 'bar')

    data = epochs['a'].get_data()
    event_a = events[events[:, 2] == 1]
    assert_true(len(data) == len(event_a))

    epochs = Epochs(raw, events, {'a': 1, 'b': 2}, tmin, tmax, picks=picks,
                    preload=True)
    assert_raises(KeyError, epochs.__getitem__, 'bar')
    epochs.save(op.join(tempdir, 'test-epo.fif'))
    epochs2 = read_epochs(op.join(tempdir, 'test-epo.fif'))

    for ep in [epochs, epochs2]:
        data = ep['a'].get_data()
        event_a = events[events[:, 2] == 1]
        assert_true(len(data) == len(event_a))

    assert_array_equal(epochs2['a'].events, epochs['a'].events)

    epochs3 = Epochs(raw, events, {'a': 1, 'b': 2, 'c': 3, 'd': 4},
                     tmin, tmax, picks=picks, preload=True)
    assert_equal(list(sorted(epochs3[['a', 'b']].event_id.values())),
                 [1, 2])
    epochs4 = epochs['a']
    epochs5 = epochs3['a']
    assert_array_equal(epochs4.events, epochs5.events)
    # 20 is our tolerance because epochs are written out as floats
    assert_array_almost_equal(epochs4.get_data(), epochs5.get_data(), 20)
    epochs6 = epochs3[['a', 'b']]
    assert_true(all(np.logical_or(epochs6.events[:, 2] == 1,
                                  epochs6.events[:, 2] == 2)))
    assert_array_equal(epochs.events, epochs6.events)
    assert_array_almost_equal(epochs.get_data(), epochs6.get_data(), 20)

def calc_inverse_operator(events_id, epochs_fn, fwd_sub_fn, inv_fn, min_crop_t=None, max_crop_t=0):
    for cond in events_id.keys():
        epochs = mne.read_epochs(epochs_fn.format(cond=cond))
        noise_cov = mne.compute_covariance(epochs.crop(min_crop_t, max_crop_t, copy=True))
        forward_sub = mne.read_forward_solution(fwd_sub_fn.format(cond=cond))
        inverse_operator_sub = make_inverse_operator(epochs.info, forward_sub, noise_cov,
            loose=None, depth=None)
        write_inverse_operator(inv_fn.format(cond=cond), inverse_operator_sub)

def _calc_epoches(params):
    subject, events_id, tmin, tmax = params
    out_file = op.join(LOCAL_ROOT_DIR, 'epo', '{}_ecr_nTSSS_conflict-epo.fif'.format(subject))
    if not op.isfile(out_file):
        events = mne.read_events(op.join(REMOTE_ROOT_DIR, 'events', '{}_ecr_nTSSS_conflict-eve.fif'.format(subject)))
        raw = mne.io.Raw(op.join(REMOTE_ROOT_DIR, 'raw', '{}_ecr_nTSSS_raw.fif'.format(subject)), preload=False)
        picks = mne.pick_types(raw.info, meg=True)
        epochs = find_epoches(raw, picks, events, events_id, tmin=tmin, tmax=tmax)
        epochs.save(out_file)
    else:
        epochs = mne.read_epochs(out_file)
    return epochs

def test_array_epochs():
    """Test creating epochs from array
    """

    # creating
    rng = np.random.RandomState(42)
    data = rng.random_sample((10, 20, 300))
    sfreq = 1e3
    ch_names = ["EEG %03d" % (i + 1) for i in range(20)]
    types = ["eeg"] * 20
    info = create_info(ch_names, sfreq, types)
    events = np.c_[np.arange(1, 600, 60), np.zeros(10), [1, 2] * 5]
    event_id = {"a": 1, "b": 2}
    epochs = EpochsArray(data, info, events=events, event_id=event_id, tmin=-0.2)

    # saving
    temp_fname = op.join(tempdir, "test-epo.fif")
    epochs.save(temp_fname)
    epochs2 = read_epochs(temp_fname)
    data2 = epochs2.get_data()
    assert_allclose(data, data2)
    assert_allclose(epochs.times, epochs2.times)
    assert_equal(epochs.event_id, epochs2.event_id)
    assert_array_equal(epochs.events, epochs2.events)

    # plotting
    import matplotlib

    matplotlib.use("Agg")  # for testing don't use X server
    epochs[0].plot()

    # indexing
    assert_array_equal(np.unique(epochs["a"].events[:, 2]), np.array([1]))
    assert_equal(len(epochs[:2]), 2)
    data[0, 5, 150] = 3000
    data[1, :, :] = 0
    data[2, 5, 210] = 3000
    data[3, 5, 260] = 0
    epochs = EpochsArray(
        data,
        info,
        events=events,
        event_id=event_id,
        tmin=0,
        reject=dict(eeg=1000),
        flat=dict(eeg=1e-1),
        reject_tmin=0.1,
        reject_tmax=0.2,
    )
    assert_equal(len(epochs), len(events) - 2)
    assert_equal(epochs.drop_log[0], ["EEG 006"])
    assert_equal(len(events), len(epochs.selection))

def _call_base_epochs_public_api(epochs, tmpdir):
    """Call all public API methods of an (non-empty) epochs object."""
    # make sure saving and loading returns the same data
    orig_data = epochs.get_data()
    export_file = tmpdir.join('test_rt-epo.fif')
    epochs.save(str(export_file))
    loaded_epochs = read_epochs(str(export_file))
    loaded_data = loaded_epochs.get_data()
    assert orig_data.shape == loaded_data.shape
    assert_allclose(loaded_data, orig_data)

    # decimation
    epochs_copy = epochs.copy()
    epochs_copy.decimate(1)
    assert epochs_copy.get_data().shape == orig_data.shape
    epochs_copy.info['lowpass'] = 10  # avoid warning
    epochs_copy.decimate(10)
    assert np.abs(10.0 - orig_data.shape[2] /
                  epochs_copy.get_data().shape[2]) <= 1

    # check that methods that require preloaded data fail
    with pytest.raises(RuntimeError):
        epochs.crop(tmin=epochs.tmin,
                    tmax=(epochs.tmin + (epochs.tmax - epochs.tmin) / 2))
    with pytest.raises(RuntimeError):
        epochs.drop_channels(epochs.ch_names[0:1])
    with pytest.raises(RuntimeError):
        epochs.resample(epochs.info['sfreq'] / 10)

    # smoke test
    epochs.standard_error()
    avg_evoked = epochs.average()
    epochs.subtract_evoked(avg_evoked)
    epochs.metadata
    epochs.events
    epochs.ch_names
    epochs.tmin
    epochs.tmax
    epochs.filename
    repr(epochs)
    epochs.plot(show=False)
    # save time by not calling all plot functions
    # epochs.plot_psd(show=False)
    # epochs.plot_drop_log(show=False)
    # epochs.plot_topo_image()
    # epochs.plot_psd_topomap()
    # epochs.plot_image()
    epochs.drop_bad()
    epochs_copy.apply_baseline()
    # do not call since we don't want to make assumptions about events
    # epochs_copy.equalize_event_counts(epochs.event_id.keys())
    epochs_copy.drop([0])

def load(subject, event_type):
    with open(op.join(path_data, '%s/behavior_%s.pkl') % (subject, event_type), 'rb') as f:
        events = pickle.load(f)
        # add explicit conditions
        events = complete_behavior(events)
    epochs = mne.read_epochs(op.join(path_data, '%s/epochs_%s.fif') % (subject, event_type))
    if event_type == 'Target':
        epochs.crop(0, .600)
    elif event_type == 'Cue':
        epochs.crop(0, .900)
    elif event_type == 'Probe':
        epochs.crop(0, .600)
    return epochs, events

def evok_epochs(sub_id, session):
    """ load a epoched file and average it and save the evoked file
    """
    fname = "sub_%d_%s" % (sub_id, session)
    f_load = fname + "_tsss_mc_epochs.fif"
    f_save = fname + "_tsss_mc_evk.fif"
    epochs = mne.read_epochs(f_load)

    evoked = epochs.average()
    evoked.comment = session
    evoked = mne.fiff.pick_types_evoked(evoked, meg='grad', exclude='bads')
    evoked.resample(sfreq=250)
    evoked.save(f_save)

def calc_evoked(indices, epochs_fname, overwrite_epochs=False, overwrite_evoked=False):
    epochs = mne.read_epochs(epochs_fname, preload=False)
    print(epochs.events.shape)
    for event_name, event_indices in indices.items():
        evoked_event_fname = meg.get_cond_fname(meg.EVO, event_name)
        epochs_event_fname = meg.get_cond_fname(meg.EPO, event_name)
        if not op.isfile(epochs_event_fname) or overwrite_epochs:
            print('Saving {} epochs to {}, events num: {}'.format(event_name, epochs_event_fname, len(event_indices)))
            event_epochs = epochs[event_indices]
            event_epochs.save(epochs_event_fname)
        if not op.isfile(evoked_event_fname) or overwrite_evoked:
            print('Saving {} evoked to {}'.format(event_name, evoked_event_fname))
            mne.write_evokeds(evoked_event_fname, event_epochs.average())

def load(subject, event_type):
    fname = op.join(path_data, '%s/behavior_%s.hdf5') % (subject, event_type)
    events = read_hdf5(fname)
    # add explicit conditions
    events = complete_behavior(events)
    fname = op.join(path_data, '%s/epochs_%s.fif') % (subject, event_type)
    epochs = mne.read_epochs(fname)
    if event_type == 'Target':
        epochs.crop(0, .600)
    elif event_type == 'Cue':
        epochs.crop(0, .900)
    elif event_type == 'Probe':
        epochs.crop(0, .600)
    return epochs, events

def test_access_by_name():
    """Test accessing epochs by event name and on_missing for rare events
    """
    assert_raises(ValueError, Epochs, raw, events, {1: 42, 2: 42}, tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {"a": "spam", 2: "eggs"}, tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {"a": "spam", 2: "eggs"}, tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, "foo", tmin, tmax, picks=picks)
    # Test accessing non-existent events (assumes 12345678 does not exist)
    event_id_illegal = dict(aud_l=1, does_not_exist=12345678)
    assert_raises(ValueError, Epochs, raw, events, event_id_illegal, tmin, tmax)
    # Test on_missing
    assert_raises(ValueError, Epochs, raw, events, 1, tmin, tmax, on_missing="foo")
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        Epochs(raw, events, event_id_illegal, tmin, tmax, on_missing="warning")
        nw = len(w)
        assert_true(1 <= nw <= 2)
        Epochs(raw, events, event_id_illegal, tmin, tmax, on_missing="ignore")
        assert_equal(len(w), nw)
    epochs = Epochs(raw, events, {"a": 1, "b": 2}, tmin, tmax, picks=picks)
    assert_raises(KeyError, epochs.__getitem__, "bar")

    data = epochs["a"].get_data()
    event_a = events[events[:, 2] == 1]
    assert_true(len(data) == len(event_a))

    epochs = Epochs(raw, events, {"a": 1, "b": 2}, tmin, tmax, picks=picks, preload=True)
    assert_raises(KeyError, epochs.__getitem__, "bar")
    epochs.save(op.join(tempdir, "test-epo.fif"))
    epochs2 = read_epochs(op.join(tempdir, "test-epo.fif"))

    for ep in [epochs, epochs2]:
        data = ep["a"].get_data()
        event_a = events[events[:, 2] == 1]
        assert_true(len(data) == len(event_a))

    assert_array_equal(epochs2["a"].events, epochs["a"].events)

    epochs3 = Epochs(raw, events, {"a": 1, "b": 2, "c": 3, "d": 4}, tmin, tmax, picks=picks, preload=True)
    assert_equal(list(sorted(epochs3[["a", "b"]].event_id.values())), [1, 2])
    epochs4 = epochs["a"]
    epochs5 = epochs3["a"]
    assert_array_equal(epochs4.events, epochs5.events)
    # 20 is our tolerance because epochs are written out as floats
    assert_array_almost_equal(epochs4.get_data(), epochs5.get_data(), 20)
    epochs6 = epochs3[["a", "b"]]
    assert_true(all(np.logical_or(epochs6.events[:, 2] == 1, epochs6.events[:, 2] == 2)))
    assert_array_equal(epochs.events, epochs6.events)
    assert_array_almost_equal(epochs.get_data(), epochs6.get_data(), 20)

def test_read_epochs():
    event_id = 1
    tmin = -0.2
    tmax = 0.5

    #   Setup for reading the raw data
    raw = fiff.setup_read_raw(raw_fname)
    events = mne.read_events(event_name)

    # Set up pick list: MEG + STI 014 - bad channels (modify to your needs)
    include = ["STI 014"]
    want_meg = True
    want_eeg = False
    want_stim = False
    picks = fiff.pick_types(raw["info"], want_meg, want_eeg, want_stim, include, raw["info"]["bads"])

    data, times, channel_names = mne.read_epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0))

def test_access_by_name():
    """Test accessing epochs by event name
    """
    assert_raises(ValueError, Epochs, raw, events, {1: 42, 2: 42}, tmin,
                  tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
                  tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
                  tmin, tmax, picks=picks)
    assert_raises(ValueError, Epochs, raw, events, 'foo', tmin, tmax,
                  picks=picks)
    epochs = Epochs(raw, events, {'a': 1, 'b': 2}, tmin, tmax, picks=picks)
    assert_raises(KeyError, epochs.__getitem__, 'bar')

    data = epochs['a'].get_data()
    event_a = events[events[:, 2] == 1]
    assert_true(len(data) == len(event_a))

    epochs = Epochs(raw, events, {'a': 1, 'b': 2}, tmin, tmax, picks=picks,
                    preload=True)
    assert_raises(KeyError, epochs.__getitem__, 'bar')
    epochs.save(op.join(tempdir, 'test-epo.fif'))
    epochs2 = read_epochs(op.join(tempdir, 'test-epo.fif'))

    for ep in [epochs, epochs2]:
        data = ep['a'].get_data()
        event_a = events[events[:, 2] == 1]
        assert_true(len(data) == len(event_a))

    assert_array_equal(epochs2['a'].events, epochs['a'].events)

    epochs3 = Epochs(raw, events, {'a': 1, 'b': 2, 'c': 3, 'd': 4},
                     tmin, tmax, picks=picks, preload=True)
    assert_equal(list(sorted(epochs3[['a', 'b']].event_id.values())),
                 [1, 2])
    epochs4 = epochs['a']
    epochs5 = epochs3['a']
    assert_array_equal(epochs4.events, epochs5.events)
    # 20 is our tolerance because epochs are written out as floats
    assert_array_almost_equal(epochs4.get_data(), epochs5.get_data(), 20)
    epochs6 = epochs3[['a', 'b']]
    assert_true(all(np.logical_or(epochs6.events[:, 2] == 1,
                                  epochs6.events[:, 2] == 2)))
    assert_array_equal(epochs.events, epochs6.events)
    assert_array_almost_equal(epochs.get_data(), epochs6.get_data(), 20)

def ep2ts(fif_file):
    """Read fif file with raw data or epochs and save
    timeseries to .npy
    """
    from mne import read_epochs
    from mne import pick_types

    from numpy import save
    import os.path as op

    with nostdout():
        epochs = read_epochs(fif_file)

    epochs_meg = epochs.pick_types(meg=True, eeg=False, eog=False, ecg=False)

    data = epochs_meg.get_data()
    save_path = op.abspath('ts_epochs.npy')
    save(save_path, data)
    return save_path

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 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)

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 run_time_frequency(subject_id):
    print("processing subject: %s" % subject_id)
    subject = "sub%03d" % subject_id
    data_path = op.join(meg_dir, subject)
    epochs = mne.read_epochs(op.join(data_path, '%s-epo.fif' % subject))

    faces = epochs['face']
    idx = [faces.ch_names.index('EEG070')]
    power_faces, itc_faces = mne.time_frequency.tfr_morlet(
        faces, freqs=freqs, return_itc=True, n_cycles=n_cycles, picks=idx)
    power_scrambled, itc_scrambled = mne.time_frequency.tfr_morlet(
        epochs['scrambled'], freqs=freqs, return_itc=True, n_cycles=n_cycles,
        picks=idx)

    power_faces.save(op.join(data_path, '%s-faces-tfr.h5' % subject),
                     overwrite=True)
    itc_faces.save(op.join(data_path, '%s-itc_faces-tfr.h5' % subject),
                   overwrite=True)

    power_scrambled.save(op.join(data_path, '%s-scrambled-tfr.h5' % subject),
                         overwrite=True)
    itc_scrambled.save(op.join(data_path, '%s-itc_scrambled-tfr.h5' % subject),
                       overwrite=True)