def test_cov_estimation_on_raw_segment():
    """Test estimation from raw on continuous recordings (typically empty room)
    """
    tempdir = _TempDir()
    raw = Raw(raw_fname, preload=False)
    cov = compute_raw_data_covariance(raw)
    cov_mne = read_cov(erm_cov_fname)
    assert_true(cov_mne.ch_names == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data, ord='fro')
                / linalg.norm(cov.data, ord='fro') < 1e-4)

    # test IO when computation done in Python
    cov.save(op.join(tempdir, 'test-cov.fif'))  # test saving
    cov_read = read_cov(op.join(tempdir, 'test-cov.fif'))
    assert_true(cov_read.ch_names == cov.ch_names)
    assert_true(cov_read.nfree == cov.nfree)
    assert_array_almost_equal(cov.data, cov_read.data)

    # test with a subset of channels
    picks = pick_channels(raw.ch_names, include=raw.ch_names[:5])
    cov = compute_raw_data_covariance(raw, picks=picks)
    assert_true(cov_mne.ch_names[:5] == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data[picks][:, picks],
                ord='fro') / linalg.norm(cov.data, ord='fro') < 1e-4)
    # make sure we get a warning with too short a segment
    raw_2 = raw.crop(0, 1)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        cov = compute_raw_data_covariance(raw_2)
    assert_true(len(w) == 1)

def test_cov_estimation_on_raw_segment():
    """Estimate raw on continuous recordings (typically empty room)
    """
    raw = Raw(raw_fname)
    cov = compute_raw_data_covariance(raw)
    cov_mne = read_cov(erm_cov_fname)
    assert_true(cov_mne.ch_names == cov.ch_names)
    print (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro'))
    assert_true(linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-6

    # test IO when computation done in Python
    cov.save('test-cov.fif')  # test saving
    cov_read = read_cov('test-cov.fif')
    assert_true(cov_read.ch_names == cov.ch_names)
    assert_true(cov_read.nfree == cov.nfree)
    assert_true((linalg.norm(cov.data - cov_read.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-5)

    # test with a subset of channels
    picks = pick_channels(raw.ch_names, include=raw.ch_names[:5])
    cov = compute_raw_data_covariance(raw, picks=picks)
    assert_true(cov_mne.ch_names[:5] == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data[picks][:, picks],
                ord='fro') / linalg.norm(cov.data, ord='fro')) < 1e-6

def get_epochs_and_cov(X, y, window=500):
    """return epochs from array."""
    raw_train = toMNE(X, y)
    picks = range(len(getChannelNames()))

    events = list()
    events_id = dict()
    for j, eid in enumerate(getEventNames()):
        tmp = find_events(raw_train, stim_channel=eid, verbose=False)
        tmp[:, -1] = j + 1
        events.append(tmp)
        events_id[eid] = j + 1

    # concatenate and sort events
    events = np.concatenate(events, axis=0)
    order_ev = np.argsort(events[:, 0])
    events = events[order_ev]

    epochs = Epochs(raw_train, events, events_id,
                    tmin=-(window / 500.0) + 1 / 500.0 + 0.150,
                    tmax=0.150, proj=False, picks=picks, baseline=None,
                    preload=True, add_eeg_ref=False, verbose=False)

    cov_signal = compute_raw_data_covariance(raw_train, verbose=False)
    return epochs, cov_signal

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_lcmv_raw():
    """Test LCMV with raw data
    """
    raw, _, _, _, noise_cov, label, forward, _, _, _ =\
        _get_data(all_forward=False, epochs=False, data_cov=False)

    tmin, tmax = 0, 20
    start, stop = raw.time_as_index([tmin, tmax])

    # use only the left-temporal MEG channels for LCMV
    left_temporal_channels = mne.read_selection('Left-temporal')
    picks = mne.fiff.pick_types(raw.info, meg=True, exclude='bads',
                                selection=left_temporal_channels)

    data_cov = mne.compute_raw_data_covariance(raw, tmin=tmin, tmax=tmax)

    stc = lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=label,
                   start=start, stop=stop, picks=picks)

    assert_array_almost_equal(np.array([tmin, tmax]),
                              np.array([stc.times[0], stc.times[-1]]),
                              decimal=2)

    # make sure we get an stc with vertices only in the lh
    vertno = [forward['src'][0]['vertno'], forward['src'][1]['vertno']]
    assert_true(len(stc.vertno[0]) == len(np.intersect1d(vertno[0],
                                                         label.vertices)))
    assert_true(len(stc.vertno[1]) == 0)

def test_maxwell_filter_additional():
    """Test processing of Maxwell filtered data"""

    # TODO: Future tests integrate with mne/io/tests/test_proc_history

    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    data_path = op.join(testing.data_path(download=False))

    file_name = 'test_move_anon'

    raw_fname = op.join(data_path, 'SSS', file_name + '_raw.fif')

    with warnings.catch_warnings(record=True):  # maxshield
        # Use 2.0 seconds of data to get stable cov. estimate
        raw = Raw(raw_fname, preload=False, proj=False,
                  allow_maxshield=True).crop(0., 2., False)

    # Get MEG channels, compute Maxwell filtered data
    raw.preload_data()
    raw.pick_types(meg=True, eeg=False)
    int_order, ext_order = 8, 3
    raw_sss = maxwell.maxwell_filter(raw, int_order=int_order,
                                     ext_order=ext_order)

    # Test io on processed data
    tempdir = _TempDir()
    test_outname = op.join(tempdir, 'test_raw_sss.fif')
    raw_sss.save(test_outname)
    raw_sss_loaded = Raw(test_outname, preload=True, proj=False,
                         allow_maxshield=True)

    # Some numerical imprecision since save uses 'single' fmt
    assert_allclose(raw_sss_loaded._data[:, :], raw_sss._data[:, :],
                    rtol=1e-6, atol=1e-20)

    # Test rank of covariance matrices for raw and SSS processed data
    cov_raw = compute_raw_data_covariance(raw)
    cov_sss = compute_raw_data_covariance(raw_sss)

    scalings = None
    cov_raw_rank = _estimate_rank_meeg_cov(cov_raw['data'], raw.info, scalings)
    cov_sss_rank = _estimate_rank_meeg_cov(cov_sss['data'], raw_sss.info,
                                           scalings)

    assert_equal(cov_raw_rank, raw.info['nchan'])
    assert_equal(cov_sss_rank, maxwell.get_num_moments(int_order, 0))

def test_cov_estimation_on_raw_segment():
    """Estimate raw on continuous recordings (typically empty room)
    """
    raw = Raw(raw_fname)
    cov = mne.compute_raw_data_covariance(raw)
    cov_mne = mne.Covariance(erm_cov_fname)
    assert cov_mne.ch_names == cov.ch_names
    print (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro'))
    assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-6

    def _fit_xdawn(self, raw, epochs=None):
        # Epochs can be computed from raw info, but might have been computed already
        if epochs is None:
            epochs = self._epochs(raw, self._events())

        sig_cov = mne.compute_raw_data_covariance(raw, picks=self._picks(raw), verbose=False).data
        evoked_cov = np.cov(epochs.average().data)
        evals, evecs = eigh(evoked_cov, sig_cov)
        # sort eigenvectors by their corresponding eigenvalues and normalize them
        evecs = evecs[:, np.argsort(evals)[::-1]]
        self._V = evecs / np.sqrt(np.sum(evecs ** 2, axis=0))
        self._A = inv(self._V.T)

def comp_noise(sub_id):
    """ This function compute and save the noise covariance matrix from
    the empty room recordings
    """

    fname = "sub_%d_empty_room.fif" % sub_id
    outname = "sub_%d_empty_room-cov.fif" % sub_id
    raw = mne.fiff.Raw(fname, preload=True)

    picks = mne.fiff.pick_types(raw.info, meg=True, eeg=True, stim=False,
                                eog=True, exclude='bads')
    cov = mne.compute_raw_data_covariance(raw, picks=picks, reject=None)

    cov.save(outname)

def localize_activity(subj, tmax=np.Inf, clean_only=True, reg=0):
    import find_good_segments as fgs

    data_path = env.data + '/MEG_data/fifs/'
    fwd_path = env.data + '/MEG_data/analysis/rest/'
    fwd_fname = fwd_path + subj + '_rest_LP100_HP0.6_CP3_DS300_raw-5-fwd.fif'

    # preloading makes computing the covariance a lot faster
    raw = mne.fiff.Raw(data_path + subj + '_rest_LP100_HP0.6_CP3_DS300_raw.fif', preload=True)

    picks = mne.fiff.pick_channels_regexp(raw.info['ch_names'], 'M..-*')

    # we don't need to choose picks because we only work with MEG channels, and
    # all channels are good
    fwd = mne.read_forward_solution(fwd_fname)

    # Load real data as templates
    if clean_only:
        start, end, num_chans = fgs.find_good_segments(subj, threshold=3500e-15)
        if start == 0:
            start = start + 3

        if end - start > tmax:
            start = end - tmax

        cov = mne.compute_raw_data_covariance(raw, tmin=start, tmax=end, picks=picks)
    else:
        cov = mne.compute_raw_data_covariance(raw, picks=picks)

    weights = calculate_weights(fwd, cov, reg=reg)

    data, times = raw[picks, raw.time_as_index(start):raw.time_as_index(end)]
    print 'Multiplying data by beamformer weights...'
    sol = np.dot(weights, data)
    src = mne.SourceEstimate(sol, [fwd['src'][0]['vertno'], fwd['src'][1]['vertno']], times[0], times[1] - times[0])
    return src

def test_xdawn_fit():
    """Test Xdawn fit."""
    # get data
    raw, events, picks = _get_data()
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    preload=True, baseline=None, verbose=False)
    # =========== Basic Fit test =================
    # test base xdawn
    xd = Xdawn(n_components=2, correct_overlap='auto',
               signal_cov=None, reg=None)
    xd.fit(epochs)
    # with this parameters, the overlapp correction must be False
    assert_equal(xd.correct_overlap, False)
    # no overlapp correction should give averaged evoked
    evoked = epochs['cond2'].average()
    assert_array_equal(evoked.data, xd.evokeds_['cond2'].data)

    # ========== with signal cov provided ====================
    # provide covariance object
    signal_cov = compute_raw_data_covariance(raw, picks=picks)
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    xd.fit(epochs)
    # provide ndarray
    signal_cov = np.eye(len(picks))
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    xd.fit(epochs)
    # provide ndarray of bad shape
    signal_cov = np.eye(len(picks) - 1)
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    assert_raises(ValueError, xd.fit, epochs)
    # provide another type
    signal_cov = 42
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    assert_raises(ValueError, xd.fit, epochs)
    # fit with baseline correction and ovverlapp correction should throw an
    # error
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    preload=True, baseline=(None, 0), verbose=False)

    xd = Xdawn(n_components=2, correct_overlap=True)
    assert_raises(ValueError, xd.fit, epochs)

def test_lcmv_raw():
    """Test LCMV with raw data
    """
    forward = mne.read_forward_solution(fname_fwd)
    label = mne.read_label(fname_label)
    noise_cov = mne.read_cov(fname_cov)
    raw = mne.fiff.Raw(fname_raw, preload=False)

    tmin, tmax = 0, 20
    # Setup for reading the raw data
    raw.info['bads'] = ['MEG 2443', 'EEG 053']  # 2 bads channels

    # Set up pick list: EEG + MEG - bad channels (modify to your needs)
    left_temporal_channels = mne.read_selection('Left-temporal')
    picks = mne.fiff.pick_types(raw.info, meg=True, eeg=False, stim=True,
                                eog=True, exclude='bads',
                                selection=left_temporal_channels)

    noise_cov = mne.read_cov(fname_cov)
    noise_cov = mne.cov.regularize(noise_cov, raw.info,
                                   mag=0.05, grad=0.05, eeg=0.1, proj=True)

    start, stop = raw.time_as_index([tmin, tmax])

    # use only the left-temporal MEG channels for LCMV
    picks = mne.fiff.pick_types(raw.info, meg=True, exclude='bads',
                                selection=left_temporal_channels)

    data_cov = mne.compute_raw_data_covariance(raw, tmin=tmin, tmax=tmax)

    stc = lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=label,
                   start=start, stop=stop, picks=picks)

    assert_array_almost_equal(np.array([tmin, tmax]),
                              np.array([stc.times[0], stc.times[-1]]),
                              decimal=2)

    # make sure we get an stc with vertices only in the lh
    vertno = [forward['src'][0]['vertno'], forward['src'][1]['vertno']]
    assert_true(len(stc.vertno[0]) == len(np.intersect1d(vertno[0],
                                                         label.vertices)))
    assert_true(len(stc.vertno[1]) == 0)

def create_noise_covariance_matrix(fname_empty_room, fname_out=None, verbose=None):
    """Creates the noise covariance matrix from an empty room file"""

    print ">>>> estimate noise covariance matrix from empty room file..."

    # read in data
    raw_empty = mne.fiff.Raw(fname_empty_room,
                             verbose=verbose)
    # filter data

    # pick only MEG channels
    picks = mne.fiff.pick_types(raw_empty.info,
                                meg=True,
                                exclude='bads')

    # calculate noise-covariance matrix
    noise_cov_mat = mne.compute_raw_data_covariance(raw_empty,
                                                    picks=picks,
                                                    verbose=verbose)
    # write noise-covariance matrix to disk
    if fname_out is not None:
        mne.write_cov(fname_out, noise_cov_mat)

    return noise_cov_mat

Estimate covariance matrix from a raw FIF file
==============================================

"""
# Author: Alexandre Gramfort <[email protected]>
#
# License: BSD (3-clause)

print __doc__

import mne
from mne import fiff
from mne.datasets import sample

data_path = sample.data_path('.')
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

raw = fiff.Raw(fname)

# Set up pick list: MEG + STI 014 - bad channels
cov = mne.compute_raw_data_covariance(raw, reject=dict(eeg=40e-6, eog=150e-6))
print cov

###############################################################################
# Show covariance
import pylab as pl
pl.figure()
pl.imshow(cov.data, interpolation="nearest")
pl.title('Full covariance matrix')
pl.show()

# License: BSD (3-clause)

print __doc__

import mne
from mne import fiff
from mne.datasets import sample

data_path = sample.data_path()
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

raw = fiff.Raw(fname)

include = []  # or stim channels ['STI 014']
raw.info['bads'] += ['EEG 053']  # bads + 1 more

# pick EEG channels
picks = fiff.pick_types(raw.info, meg=True, eeg=True, stim=False, eog=True,
                                            include=include, exclude='bads')
# setup rejection
reject = dict(eeg=80e-6, eog=150e-6)

# Compute the covariance from the raw data
cov = mne.compute_raw_data_covariance(raw, picks=picks, reject=reject)
print cov

###############################################################################
# Show covariance
mne.viz.plot_cov(cov, raw.info, colorbar=True, proj=True)
# try setting proj to False to see the effect

Process the files: <subject>_emptyroom.fif
Compute the noise matrix
Author: Praveen, Qunxi
"""
import mne, sys, os
import pylab as pl
import numpy as np

flow, fhigh = 1.0, 45.0
filter_type = 'butter'
filter_order = 4
njobs = 4

try:
    subject = sys.argv[1]#Get the subject
except:
    print "Please run with input file provided. Exiting"
    sys.exit()
subjects_dir = '/home/qdong/freesurfer/subjects/'
subject_path = subjects_dir + subject#Set the data path of the subject

raw_empty_fname = subject_path + '/MEG/%s_emptyroom.fif' %subject
raw_empty = mne.fiff.Raw(raw_empty_fname, preload=True)
#Filter the empty room data
picks_empty = mne.fiff.pick_types(raw_empty.info, meg=True, eeg=False, eog=True, ecg=True, stim=False, exclude='bads') 
raw_empty.filter(flow, fhigh, picks=picks_empty, n_jobs=njobs, method='iir', 
           iir_params={'ftype': filter_type, 'order': filter_order})          
#Get the basename
raw_empty_basename = os.path.splitext(os.path.basename(raw_empty_fname))[0]
cov = mne.compute_raw_data_covariance(raw_empty, picks=picks_empty)
mne.write_cov(subject_path+'/MEG/%s_cov.fif' %(raw_empty_basename), cov)

stc_data[1] = np.roll(stc_data[1], 80)
"""

stc = generate_sparse_stc(fwd["src"], labels, stc_data, tmin, tstep)

###############################################################################
# Generate noisy evoked data
picks = mne.fiff.pick_types(raw.info, meg=True)
iir_filter = iir_filter_raw(raw, order=5, picks=picks, tmin=60, tmax=180)
evoked = generate_evoked(fwd, stc, evoked_template, cov, snr)  # , iir_filter=iir_filter)
raw[:306, :ntimes] = evoked.data
raw = raw.crop(tmax=evoked.times[-1])
raw.save("t2-raw.fif")

raw2 = mne.fiff.Raw("t2-raw.fif")
cov = mne.compute_raw_data_covariance(raw2)

dsim = stc.in_label(labels[0])
pl.plot(dsim.data.T)

weights = ve.calculate_weights(fwd, cov, reg=0)

d, _ = raw2[:306, :]
sol = np.dot(weights, d)
src = mne.SourceEstimate(sol, (fwd["src"][0]["vertno"], fwd["src"][1]["vertno"]), tmin, tstep)

"""
# plotting the localized source
const = 1*10**-18
good_data = np.nonzero(abs(src.data[:,150])>const)
lv = np.nonzero(abs(src.lh_data[:,150])>const)

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

def test_rank():
    """Test cov rank estimation"""
    raw_sample = Raw(raw_fname)

    raw_sss = Raw(hp_fif_fname)
    raw_sss.add_proj(compute_proj_raw(raw_sss))

    cov_sample = compute_raw_data_covariance(raw_sample)
    cov_sample_proj = compute_raw_data_covariance(
        raw_sample.copy().apply_proj())

    cov_sss = compute_raw_data_covariance(raw_sss)
    cov_sss_proj = compute_raw_data_covariance(
        raw_sss.copy().apply_proj())

    picks_all_sample = pick_types(raw_sample.info, meg=True, eeg=True)
    picks_all_sss = pick_types(raw_sss.info, meg=True, eeg=True)

    info_sample = pick_info(raw_sample.info, picks_all_sample)
    picks_stack_sample = [('eeg', pick_types(info_sample, meg=False,
                                             eeg=True))]
    picks_stack_sample += [('meg', pick_types(info_sample, meg=True))]
    picks_stack_sample += [('all',
                            pick_types(info_sample, meg=True, eeg=True))]

    info_sss = pick_info(raw_sss.info, picks_all_sss)
    picks_stack_somato = [('eeg', pick_types(info_sss, meg=False, eeg=True))]
    picks_stack_somato += [('meg', pick_types(info_sss, meg=True))]
    picks_stack_somato += [('all',
                            pick_types(info_sss, meg=True, eeg=True))]

    iter_tests = list(itt.product(
        [(cov_sample, picks_stack_sample, info_sample),
         (cov_sample_proj, picks_stack_sample, info_sample),
         (cov_sss, picks_stack_somato, info_sss),
         (cov_sss_proj, picks_stack_somato, info_sss)],  # sss
        [dict(mag=1e15, grad=1e13, eeg=1e6)]
    ))

    for (cov, picks_list, this_info), scalings in iter_tests:
        for ch_type, picks in picks_list:

            this_very_info = pick_info(this_info, picks)

            # compute subset of projs
            this_projs = [c['active'] and
                          len(set(c['data']['col_names'])
                              .intersection(set(this_very_info['ch_names']))) >
                          0 for c in cov['projs']]
            n_projs = sum(this_projs)

            # count channel types
            ch_types = [channel_type(this_very_info, idx)
                        for idx in range(len(picks))]
            n_eeg, n_mag, n_grad = [ch_types.count(k) for k in
                                    ['eeg', 'mag', 'grad']]
            n_meg = n_mag + n_grad
            if ch_type in ('all', 'eeg'):
                n_projs_eeg = 1
            else:
                n_projs_eeg = 0

            # check sss
            if 'proc_history' in this_very_info:
                mf = this_very_info['proc_history'][0]['max_info']
                n_free = _get_sss_rank(mf)
                if 'mag' not in ch_types and 'grad' not in ch_types:
                    n_free = 0
                # - n_projs XXX clarify
                expected_rank = n_free + n_eeg
                if n_projs > 0 and ch_type in ('all', 'eeg'):
                    expected_rank -= n_projs_eeg
            else:
                expected_rank = n_meg + n_eeg - n_projs

            C = cov['data'][np.ix_(picks, picks)]
            est_rank = _estimate_rank_meeg_cov(C, this_very_info,
                                               scalings=scalings)

            assert_equal(expected_rank, est_rank)

def run():

    args = sys.argv
    if len(args) <= 1:
        print 'Usage: run_meg_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'

    src = setup_source_space(subject, fname=True, spacing='oct6', n_jobs=2,
                             overwrite=True)

    # If one wanted to use other source spaces, these types of options are
    # available
    src_fsaverage = setup_source_space('fsaverage', fname=True, spacing='ico5',
                                       n_jobs=2, overwrite=True,
                                       add_dist=False)
    morph_source_spaces(src_fsaverage, subject_to='sample')
    setup_source_space(subject, fname=True, spacing='all', overwrite=True,
                       n_jobs=2, add_dist=False)

    # Add distances to source space (if desired, takes a long time)
    bem_dir = join(subjects_dir, join('sample', 'bem'))
    os.rename(join(bem_dir, 'sample-oct-6-src.fif'),
              join(bem_dir, 'sample-oct-6-orig-src.fif'))
    new_src = add_source_space_distances(src, dist_limit=0.007)
    new_src.save(join(bem_dir, 'sample-oct-6-src.fif'))

    # Preprocessing
    raw = mne.io.Raw(join(meg_dir, 'sample_audvis_raw.fif'), preload=True)
    raw.info['bads'] = ['MEG 2443', 'EEG 053']
    reject = dict(grad=3000e-13, mag=4000e-15, eeg=100e-6)
    ecg_proj, _ = mne.preprocessing.compute_proj_ecg(raw, l_freq=1, h_freq=100,
                                                     ch_name='MEG 1531',
                                                     reject=reject)

    eog_proj, _ = mne.preprocessing.compute_proj_eog(raw, l_freq=1, h_freq=35,
                                                     reject=reject,
                                                     no_proj=True)

    events = mne.find_events(raw)
    mne.write_events(join(meg_dir, 'sample_audvis_raw-eve.fif'), events)
    event_id = [1, 2, 3, 4]
    tmin, tmax = -0.2, 0.5
    picks = mne.pick_types(raw.info, meg=True, eeg=True, stim=True, eog=True)

    # Average with no filter
    epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks)
    evoked = epochs.average()
    evoked.save(join(meg_dir, 'sample_audvis-no-filter-ave.fif'))

    raw.filter(l_freq=None, h_freq=40)
    raw_resampled = raw.resample(150)
    raw_resampled.save(join(meg_dir, 'sample_audvis_filt-0-40_raw.fif'),
                       overwrite=True)

    raw.add_proj(ecg_proj)
    raw.add_proj(eog_proj)

    resampled_events = mne.find_events(raw_resampled)
    mne.write_events(join(meg_dir, 'sample_audvis_filt-0-40_raw-eve.fif'),
                     resampled_events)

    # Average with filter
    epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks)
    evoked = epochs.average()
    evoked.save(join(meg_dir, 'sample_audvis-ave.fif'))

    # Compute the noise covariance matrix
    noise_cov = mne.compute_raw_data_covariance(raw, picks=picks)
    noise_cov.save(join(meg_dir, 'audvis.cov'))

    # Compute the empty-room noise covariance matrix
    ernoise_raw = mne.io.Raw(join(meg_dir, 'ernoise_raw.fif'), preload=True)
    ernoise_raw.info['bads'] = ['MEG 2443']
    ernoise_raw.filter(l_freq=None, h_freq=40)
    picks = mne.pick_types(ernoise_raw.info, meg=True, eeg=True, stim=True,
                           eog=True)
    ernoise_cov = mne.compute_raw_data_covariance(ernoise_raw, picks=picks)
    ernoise_cov.save(join(meg_dir, 'ernoise.cov'))

###############################################################################
    # Compute forward solution a.k.a. lead field
    trans = join(meg_dir, 'sample_audvis_raw-trans.fif')
    bem = join(subjects_dir, 'sample', 'bem', 'sample-5120-bem-sol.fif')
    # for MEG only
    fname = join(meg_dir, 'sample_audvis-meg-oct-6-fwd.fif')
    fwd_meg = mne.make_forward_solution(raw.info, trans, src, bem,
                                        fname=fname, meg=True, eeg=False,
                                        mindist=5.0, n_jobs=2, overwrite=True)

    # for EEG only
    bem = join(subjects_dir, 'sample', 'bem',
#.........这里部分代码省略.........