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 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 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 fiff_mne(ds, fwd='{fif}*fwd.fif', cov='{fif}*cov.fif', label=None, name=None,
             tstart= -0.1, tstop=0.6, baseline=(None, 0)):
    """
    adds data from one label as

    """
    if name is None:
        if label:
            _, lbl = os.path.split(label)
            lbl, _ = os.path.splitext(lbl)
            name = lbl.replace('-', '_')
        else:
            name = 'stc'

    info = ds.info['info']

    raw = ds.info['raw']
    fif_name = raw.info['filename']
    fif_name, _ = os.path.splitext(fif_name)
    if fif_name.endswith('raw'):
        fif_name = fif_name[:-3]

    fwd = fwd.format(fif=fif_name)
    if '*' in fwd:
        d, n = os.path.split(fwd)
        names = fnmatch.filter(os.listdir(d), n)
        if len(names) == 1:
            fwd = os.path.join(d, names[0])
        else:
            raise IOError("No unique fwd file matching %r" % fwd)

    cov = cov.format(fif=fif_name)
    if '*' in cov:
        d, n = os.path.split(cov)
        names = fnmatch.filter(os.listdir(d), n)
        if len(names) == 1:
            cov = os.path.join(d, names[0])
        else:
            raise IOError("No unique cov file matching %r" % cov)

    fwd = mne.read_forward_solution(fwd, force_fixed=False, surf_ori=True)
    cov = mne.Covariance(cov)
    inv = _mn.make_inverse_operator(info, fwd, cov, loose=0.2, depth=0.8)
    epochs = mne_Epochs(ds, tstart=tstart, tstop=tstop, baseline=baseline)

    # mne example:
    snr = 3.0
    lambda2 = 1.0 / snr ** 2

    if label is not None:
        label = mne.read_label(label)
    stcs = _mn.apply_inverse_epochs(epochs, inv, lambda2, dSPM=False, label=label)

    x = np.vstack(s.data.mean(0) for s in stcs)
    s = stcs[0]
    dims = ('case', var(s.times, 'time'),)
    ds[name] = ndvar(x, dims, properties=None, info='')

    return stcs

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 run():

    args = sys.argv
    if len(args) <= 1:
        print 'Usage: run_anatomy_tutorial.sh <sample data directory>'
        return

    sample_dir = args[1]
    subjects_dir = join(sample_dir, 'subjects')
    meg_dir = join(sample_dir, 'MEG', 'sample')

    os.environ['SUBJECTS_DIR'] = subjects_dir
    os.environ['MEG_DIR'] = meg_dir

    subject = 'sample'

    bem = join(subjects_dir, subject, 'bem', 'sample-5120-bem-sol.fif')
    mri = join(subjects_dir, subject, 'mri', 'T1.mgz')
    fname = join(subjects_dir, subject, 'bem', 'volume-7mm-src.fif')
    src = setup_volume_source_space(subject, fname=fname, pos=7, mri=mri,
                                    bem=bem, overwrite=True,
                                    subjects_dir=subjects_dir)

###############################################################################
    # Compute forward solution a.k.a. lead field

    raw = mne.io.Raw(join(meg_dir, 'sample_audvis_raw.fif'))
    fwd_fname = join(meg_dir, 'sample_audvis-meg-vol-7-fwd.fif')
    trans = join(meg_dir, 'sample_audvis_raw-trans.fif')
    # for MEG only
    fwd = make_forward_solution(raw.info, trans=trans, src=src, bem=bem,
                                fname=fwd_fname, meg=True, eeg=False,
                                overwrite=True)

    # Make a sensitivity map
    smap = mne.sensitivity_map(fwd, ch_type='grad', mode='free')
    smap.save(join(meg_dir, 'sample_audvis-grad-vol-7-fwd-sensmap'), ftype='w')

###############################################################################
    # Compute MNE inverse operators
    #
    # Note: The MEG/EEG forward solution could be used for all
    #
    noise_cov = mne.read_cov(join(meg_dir, 'sample_audvis-cov.fif'))
    inv = make_inverse_operator(raw.info, fwd, noise_cov)
    fname = join(meg_dir, 'sample_audvis-meg-vol-7-meg-inv.fif')
    write_inverse_operator(fname, inv)

def compute_ts_inv_sol(raw, fwd_filename, cov_fname, snr, inv_method, aseg):
    import os.path as op
    import numpy as np
    import mne
    from mne.minimum_norm import make_inverse_operator, apply_inverse_raw
    from nipype.utils.filemanip import split_filename as split_f

    print '***** READ FWD SOL %s *****' % fwd_filename
    forward = mne.read_forward_solution(fwd_filename)

    # Convert to surface orientation for cortically constrained
    # inverse modeling
    if not aseg:
        forward = mne.convert_forward_solution(forward, surf_ori=True)

    lambda2 = 1.0 / snr ** 2

    # compute inverse operator
    print '***** COMPUTE INV OP *****'
    inverse_operator = make_inverse_operator(raw.info, forward, cov_fname,
                                             loose=0.2, depth=0.8)

    # apply inverse operator to the time windows [t_start, t_stop]s
    # TEST
    t_start = 0  # sec
    t_stop = 3  # sec
    start, stop = raw.time_as_index([t_start, t_stop])
    print '***** APPLY INV OP ***** [%d %d]sec' % (t_start, t_stop)
    stc = apply_inverse_raw(raw, inverse_operator, lambda2, inv_method,
                            label=None,
                            start=start, stop=stop, pick_ori=None)

    print '***'
    print 'stc dim ' + str(stc.shape)
    print '***'

    subj_path, basename, ext = split_f(raw.info['filename'])
    data = stc.data

    print 'data dim ' + str(data.shape)

    # save results in .npy file that will be the input for spectral node
    print '***** SAVE SOL *****'
    ts_file = op.abspath(basename + '.npy')
    np.save(ts_file, data)

    return ts_file

def estimate_inverse_solution(info,
                              noise_cov_mat,
                              fwd_sol=None,
                              fname_fwd=None,
                              fname_inv=None):
    """"Estimates inverse solution for the given data set."""

    if fwd_sol is not None:
        pass
    elif fname_fwd is not None:
        fwd_sol = mne.read_forward_solution(fname_fwd,
                                            surf_ori=True)
    else:
        print "ERROR: Neither a forward solution given nor the filename of one!"
        sys.exit()

    # restrict forward solution as necessary for MEG
    fwd = mne.fiff.pick_types_forward(fwd_sol,
                                      meg=True,
                                      eeg=False)


    # # regularize noise covariance
    # # --> not necessary as the data set to estimate the
    # #     noise-covariance matrix is quiet long, i.e.
    # #     the noise-covariance matrix is robust
    # noise_cov_mat = mne.cov.regularize(noise_cov_mat,
    #                                    info,
    #                                    mag=0.1,
    #                                    proj=True,
    #                                    verbose=verbose)

    # create the MEG inverse operators
    print ">>>> estimate inverse operator..."
    inverse_operator = min_norm.make_inverse_operator(info,
                                                      fwd,
                                                      noise_cov_mat,
                                                      loose=0.2,
                                                      depth=0.8)

    if fname_inv is not None:
        min_norm.write_inverse_operator(fname_inv, inverse_operator)

    return inverse_operator

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 apply_inverse_oper(fnepo, tmin=-0.2, tmax=0.8, subjects_dir=None):
    '''
    Apply inverse operator
    Parameter
    ---------
    fnepo: string or list
        The epochs file with ECG, EOG and environmental noise free.
    tmax, tmax:float
        The time period (second) of each epoch.
    '''
    # Get the default subjects_dir
    from mne import make_forward_solution
    from mne.minimum_norm import make_inverse_operator, write_inverse_operator

    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]
        subject_path = subjects_dir + '/%s' % subject
        fn_trans = fn_path + '/%s-trans.fif' % subject
        fn_cov = fn_path + '/%s_empty-cov.fif' % subject
        fn_src = subject_path + '/bem/%s-oct-6-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        fn_inv = fn_path + '/%s_epo-inv.fif' % subject

        epochs = mne.read_epochs(fname)
        epochs.crop(tmin, tmax)
        epochs.pick_types(meg=True, ref_meg=False)
        noise_cov = mne.read_cov(fn_cov)
        fwd = make_forward_solution(epochs.info, fn_trans, fn_src, fn_bem)
        fwd['surf_ori'] = True
        inv = make_inverse_operator(epochs.info, fwd, noise_cov, loose=0.2,
                                    depth=0.8, limit_depth_chs=False)
        write_inverse_operator(fn_inv, inv)

def apply_inverse_ave(fnevo, subjects_dir):
    ''' Make individual inverse operator.

        Parameter
        ---------
        fnevo: string or list
            The evoked file with ECG, EOG and environmental noise free.
        subjects_dir: The total bath of all the subjects.

    '''
    from mne import make_forward_solution
    from mne.minimum_norm import make_inverse_operator, write_inverse_operator
    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)
        subject = name.split('_')[0]
        fn_inv = fn_path + '/%s_fibp1-45,ave-inv.fif' % subject
        subject_path = subjects_dir + '/%s' % subject

        fn_trans = fn_path + '/%s-trans.fif' % subject
        fn_cov = fn_path + '/%s_empty,fibp1-45-cov.fif' % subject
        fn_src = subject_path + '/bem/%s-oct-6-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        [evoked] = mne.read_evokeds(fname)
        evoked.pick_types(meg=True, ref_meg=False)
        noise_cov = mne.read_cov(fn_cov)
        # noise_cov = mne.cov.regularize(noise_cov, evoked.info,
        #                               mag=0.05, grad=0.05, proj=True)
        fwd = make_forward_solution(evoked.info, fn_trans, fn_src, fn_bem)
        fwd['surf_ori'] = True
        inv = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.2,
                                    depth=0.8, limit_depth_chs=False)
        write_inverse_operator(fn_inv, inv)

evoked = Evoked(fname_evoked, setno=0, baseline=(None, 0))
forward_meeg = mne.read_forward_solution(fname_fwd_meeg, 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_meeg, meg=True, eeg=False)
# Alternatively, you can just load a forward solution that is restricted
forward_eeg = mne.read_forward_solution(fname_fwd_eeg, surf_ori=True)

# make an M/EEG, MEG-only, and EEG-only inverse operators
info = evoked.info
inverse_operator_meeg = make_inverse_operator(info, forward_meeg, noise_cov,
                                              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",

fname_fwd = fwd_dir + '%s_task-5-fwd.fif' % subj
forward = mne.read_forward_solution(fname_fwd)
evoked_fname = data_dir + '%s_stop_parsed_matched_clean_BP1-100_DS300-ave.fif' % subj
evoked = mne.read_evokeds(evoked_fname)
epochs_fname = data_dir + '%s_stop_parsed_matched_clean_BP1-100_DS300-epo.fif.gz' % subj
epochs = mne.read_epochs(epochs_fname)

# we're interested in time points every 50ms (for now). That's sfreq of 20Hz.
# the niquist there is 10Hz, so let's downsample our data that way.
epochs_ds = epochs.copy()
epochs_ds.resample(20)
evoked_ds = [epochs_ds[name].average() for name in ['STI-correct', 'STI-incorrect']]

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

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

# Define epochs for left-auditory condition
event_id, tmin, tmax = 1, -0.2, 0.5
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
                                                    eog=150e-6))

# Compute inverse solution and for each epoch
snr = 1.0           # use smaller SNR for raw data
inv_method = 'dSPM'  # sLORETA, MNE, dSPM
parc = 'aparc'      # the parcellation to use, e.g., 'aparc' 'aparc.a2009s'

lambda2 = 1.0 / snr ** 2

# Compute inverse operator
inverse_operator = make_inverse_operator(raw.info, fwd, noise_cov,
                                         loose=None, depth=None,
                                         fixed=False)


stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method,
                            pick_ori=None, return_generator=True)

# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi
labels_parc = mne.read_labels_from_annot(subject, parc=parc,
                                         subjects_dir=subjects_dir)

# Average the source estimates within each label of the cortical parcellation
# and each sub structures contained in the src space
# If mode = 'mean_flip' this option is used only for the cortical label
src = inverse_operator['src']
label_ts = mne.extract_label_time_course(stcs, labels_parc, src,

# 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.

###############################################################################

def compute_ROIs_inv_sol(raw, sbj_id, sbj_dir, fwd_filename, cov_fname, snr,
                         inv_method, parc, aseg, aseg_labels):
    import os.path as op
    import numpy as np
    import mne
    from mne.minimum_norm import make_inverse_operator, apply_inverse_raw
    from nipype.utils.filemanip import split_filename as split_f
    
    from neuropype_ephy.compute_inv_problem import get_aseg_labels

    print '***** READ noise covariance %s *****' % cov_fname
    noise_cov = mne.read_cov(cov_fname)

    print '***** READ FWD SOL %s *****' % fwd_filename
    forward = mne.read_forward_solution(fwd_filename)

    if not aseg:
        forward = mne.convert_forward_solution(forward, surf_ori=True,
                                               force_fixed=False)

    lambda2 = 1.0 / snr ** 2

    # compute inverse operator
    print '***** COMPUTE INV OP *****'
    if not aseg:
        loose = 0.2
        depth = 0.8
    else:
        loose = None
        depth = None

    inverse_operator = make_inverse_operator(raw.info, forward, noise_cov,
                                             loose=loose, depth=depth,
                                             fixed=False)

    # apply inverse operator to the time windows [t_start, t_stop]s
    print '***** APPLY INV OP *****'
    stc = apply_inverse_raw(raw, inverse_operator, lambda2, inv_method,
                            label=None,
                            start=None, stop=None,
                            buffer_size=1000,
                            pick_ori=None)  # None 'normal'

    print '***'
    print 'stc dim ' + str(stc.shape)
    print '***'

    labels_cortex = mne.read_labels_from_annot(sbj_id, parc=parc,
                                               subjects_dir=sbj_dir)

    src = inverse_operator['src']

    # allow_empty : bool -> Instead of emitting an error, return all-zero time
    # courses for labels that do not have any vertices in the source estimate
    # TODO cosa accade se la uso con solo la cortex? -> OK!!!
    label_ts = mne.extract_label_time_course_AP(stc, labels_cortex, src,
                                                mode='mean_flip',
                                                allow_empty=True,
                                                return_generator=False)

    # save results in .npy file that will be the input for spectral node
    print '***** SAVE SOL *****'
    subj_path, basename, ext = split_f(raw.info['filename'])
    ts_file = op.abspath(basename + '.npy')
    np.save(ts_file, label_ts)

    if aseg:
        labels_aseg = get_aseg_labels(src, sbj_dir, sbj_id, aseg_labels)
        labels = labels_cortex + labels_aseg
    else:
        labels = labels_cortex

    return ts_file, labels

from my_settings import (mne_folder, epochs_folder)
import sys
import mne
from mne.minimum_norm import make_inverse_operator

subject = sys.argv[1]

fwd = mne.read_forward_solution(mne_folder + "%s-fwd.fif" % subject,
                                surf_ori=False)
fwd = mne.pick_types_forward(fwd, meg="grad", eeg=False)
cov = mne.read_cov(mne_folder + "%s-cov.fif" % subject)
epochs = mne.read_epochs(epochs_folder +
                         "%s_trial_start-epo.fif" % subject,
                         preload=False)


inv = make_inverse_operator(epochs.info, fwd, cov,
                            loose=0.2, depth=0.8)

mne.minimum_norm.write_inverse_operator(mne_folder +
                                        "%s_grad-inv.fif" % subject,
                                        inv)

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)

def test_compute_LF_matrix():
    import os
    import os.path as op
    import nipype.pipeline.engine as pe
    from nipype.interfaces.mne import WatershedBEM
    import mne
    import mne.io as io
    from mne.minimum_norm import make_inverse_operator, apply_inverse_raw
    from mne.report import Report
    from nipype.utils.filemanip import split_filename as split_f
    main_path = '/home/karim/Documents/pasca/data/resting_state/'
    sbj_id = 'K0002'
    sbj_dir = op.join(main_path, 'FSF')
    bem_dir = op.join(sbj_dir, sbj_id, 'bem')
    surface_dir = op.join(sbj_dir, sbj_id, 'bem/watershed')
    data_dir = op.join(main_path, 'MEG')
    raw_fname = op.join(data_dir, '%s/%s_rest_tsss_mc.fif' % (sbj_id, sbj_id))
    raw = io.Raw(raw_fname, preload=True)
    picks = mne.pick_types(raw.info, meg=True, ref_meg=False, exclude='bads')
    raw.filter(l_freq=0.1, h_freq=300, picks=picks, method='iir', n_jobs=2)
    raw.resample(sfreq=300, npad=0)
    report = Report()
    surfaces = [sbj_id + '_brain_surface',
     sbj_id + '_inner_skull_surface',
     sbj_id + '_outer_skull_surface',
     sbj_id + '_outer_skin_surface']
    new_surfaces = ['brain.surf',
     'inner_skull.surf',
     'outer_skull.surf',
     'outer_skin.surf']
    sbj_inner_skull_fname = op.join(bem_dir, sbj_id + '-' + new_surfaces[1])
    inner_skull_fname = op.join(bem_dir, new_surfaces[1])
    if not (op.isfile(sbj_inner_skull_fname) or op.isfile(inner_skull_fname)):
        bem_IF = WatershedBEM()
        bem_IF.inputs.subject_id = sbj_id
        bem_IF.inputs.subjects_dir = sbj_dir
        bem_IF.inputs.atlas_mode = True
        bem_IF.run()
        for i in range(len(surfaces)):
            os.system('cp %s %s' % (op.join(surface_dir, surfaces[i]), op.join(bem_dir, sbj_id + '-' + new_surfaces[i])))

    else:
        print '*** inner skull surface exists!!!'
    bem = op.join(bem_dir, '%s-5120-bem-sol.fif' % sbj_id)
    if not op.isfile(bem):
        os.system('$MNE_ROOT/bin/mne_setup_forward_model --subject ' + sbj_id + ' --homog --surf --ico 4')
    else:
        print '*** BEM solution file exists!!!'
    src_fname = op.join(bem_dir, '%s-ico-5-src.fif' % sbj_id)
    if not op.isfile(src_fname):
        src = mne.setup_source_space(sbj_id, fname=True, spacing='ico5', subjects_dir=sbj_dir, overwrite=True, n_jobs=2)
    else:
        print '*** source space file exists!!!'
        src = mne.read_source_spaces(src_fname)
    trans_fname = op.join(data_dir, '%s/%s-trans.fif' % (sbj_id, sbj_id))
    data_path, basename, ext = split_f(raw_fname)
    fwd_filename = op.join(data_path, '%s-fwd.fif' % basename)
    forward = mne.make_forward_solution(raw_fname, trans_fname, src, bem, fwd_filename, n_jobs=2, overwrite=True)
    forward = mne.convert_forward_solution(forward, surf_ori=True)
    snr = 1.0
    lambda2 = 1.0 / snr ** 2
    method = 'MNE'
    reject = dict(mag=4e-12, grad=4e-10, eog=0.00025)
    noise_cov = mne.compute_raw_data_covariance(raw, picks=picks, reject=reject)
    inverse_operator = make_inverse_operator(raw.info, forward, noise_cov, loose=0.2, depth=0.8)
    start, stop = raw.time_as_index([0, 30])
    stc = apply_inverse_raw(raw, inverse_operator, lambda2, method, label=None, start=start, stop=stop, pick_ori=None)
    print '***'
    stc.shape
    print '***'
    subj_path, basename, ext = split_f(raw_fname)
    stc_filename = op.join(subj_path, basename)
    stc.save(stc_filename)
    report_filename = op.join(subj_path, basename + '-BEM-report.html')
    print report_filename
    report.save(report_filename, open_browser=False, overwrite=True)
    return