def test_sensitivity_maps():
    """Test sensitivity map computation"""
    fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
    proj_eog = read_proj(eog_fname)
    decim = 6
    for ch_type in ['eeg', 'grad', 'mag']:
        w = read_source_estimate(sensmap_fname % (ch_type, 'lh')).data
        stc = sensitivity_map(fwd, projs=None, ch_type=ch_type,
                              mode='free', exclude='bads')
        assert_array_almost_equal(stc.data, w, decim)
        assert_true(stc.subject == 'sample')
        # let's just make sure the others run
        if ch_type == 'grad':
            # fixed (2)
            w = read_source_estimate(sensmap_fname % (ch_type, '2-lh')).data
            stc = sensitivity_map(fwd, projs=None, mode='fixed',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == 'mag':
            # ratio (3)
            w = read_source_estimate(sensmap_fname % (ch_type, '3-lh')).data
            stc = sensitivity_map(fwd, projs=None, mode='ratio',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == 'eeg':
            # radiality (4), angle (5), remaining (6), and  dampening (7)
            modes = ['radiality', 'angle', 'remaining', 'dampening']
            ends = ['4-lh', '5-lh', '6-lh', '7-lh']
            for mode, end in zip(modes, ends):
                w = read_source_estimate(sensmap_fname % (ch_type, end)).data
                stc = sensitivity_map(fwd, projs=proj_eog, mode=mode,
                                      ch_type=ch_type, exclude='bads')
                assert_array_almost_equal(stc.data, w, decim)

def plot_forward(fwd, sbj_id, sbj_dir):
    import mne
    import matplotlib.pyplot as plt
    
    leadfield = fwd['sol']['data']
    print 'Leadfield size : %d x %d' % leadfield.shape
    grad_map = mne.sensitivity_map(fwd, ch_type='grad', mode='fixed')
    mag_map  = mne.sensitivity_map(fwd, ch_type='mag', mode='fixed')
    picks_meg = mne.pick_types(fwd['info'], meg=True, eeg=False)
    
    fig, axes = plt.subplots(1, 1, figsize=(10, 8), sharex=True)
    fig.suptitle('Lead field matrix (500 dipoles only)', fontsize=14)
    im = axes.imshow(leadfield[picks_meg, :500], origin='lower', aspect='auto', cmap='RdBu_r')
    axes.set_title('meg'.upper())
    axes.set_xlabel('sources')
    axes.set_ylabel('sensors')
    plt.colorbar(im, ax=axes, cmap='RdBu_r')
    plt.show()
    plt.figure()
    plt.hist([grad_map.data.ravel(), mag_map.data.ravel()], bins=20, label=['Gradiometers', 'Magnetometers'], 
              color=['c', 'b'])
    plt.legend()
    plt.title('Normal orientation sensitivity')
    plt.xlabel('sensitivity')
    plt.ylabel('count')
    plt.show()
    
    grad_map.plot(subject=sbj_id, time_label='Gradiometer sensitivity', subjects_dir=sbj_dir, clim='auto')

def test_sensitivity_maps():
    """Test sensitivity map computation"""
    fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
    projs = None
    proj_eog = read_proj(eog_fname)
    decim = 6
    for ch_type in ['eeg', 'grad', 'mag']:
        w_lh = mne.read_w(sensmap_fname % (ch_type, 'lh'))
        w_rh = mne.read_w(sensmap_fname % (ch_type, 'rh'))
        w = np.r_[w_lh['data'], w_rh['data']]
        stc = sensitivity_map(fwd, projs=projs, ch_type=ch_type,
                              mode='free', exclude='bads')
        assert_array_almost_equal(stc.data.ravel(), w, decim)
        assert_true(stc.subject == 'sample')
        # let's just make sure the others run
        if ch_type == 'grad':
            # fixed
            w_lh = mne.read_w(sensmap_fname % (ch_type, '2-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '2-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=projs, mode='fixed',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data.ravel(), w, decim)
        if ch_type == 'mag':
            # ratio
            w_lh = mne.read_w(sensmap_fname % (ch_type, '3-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '3-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=projs, mode='ratio',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data.ravel(), w, decim)
        if ch_type == 'eeg':
            # radiality (4)
            w_lh = mne.read_w(sensmap_fname % (ch_type, '4-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '4-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=projs, mode='radiality',
                                  ch_type=ch_type, exclude='bads')
            # angle (5)
            w_lh = mne.read_w(sensmap_fname % (ch_type, '5-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '5-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=proj_eog, mode='angle',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data.ravel(), w, decim)
            # remaining (6)
            w_lh = mne.read_w(sensmap_fname % (ch_type, '6-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '6-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=proj_eog, mode='remaining',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data.ravel(), w, decim)
            # dampening (7)
            w_lh = mne.read_w(sensmap_fname % (ch_type, '7-lh'))
            w_rh = mne.read_w(sensmap_fname % (ch_type, '7-rh'))
            w = np.r_[w_lh['data'], w_rh['data']]
            stc = sensitivity_map(fwd, projs=proj_eog, mode='dampening',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data.ravel(), w, decim)

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 test_sensitivity_maps():
    """Test sensitivity map computation."""
    fwd = mne.read_forward_solution(fwd_fname)
    fwd = mne.convert_forward_solution(fwd, surf_ori=True)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        projs = read_proj(eog_fname)
        projs.extend(read_proj(ecg_fname))
    decim = 6
    for ch_type in ['eeg', 'grad', 'mag']:
        w = read_source_estimate(sensmap_fname % (ch_type, 'lh')).data
        stc = sensitivity_map(fwd, projs=None, ch_type=ch_type,
                              mode='free', exclude='bads')
        assert_array_almost_equal(stc.data, w, decim)
        assert_true(stc.subject == 'sample')
        # let's just make sure the others run
        if ch_type == 'grad':
            # fixed (2)
            w = read_source_estimate(sensmap_fname % (ch_type, '2-lh')).data
            stc = sensitivity_map(fwd, projs=None, mode='fixed',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == 'mag':
            # ratio (3)
            w = read_source_estimate(sensmap_fname % (ch_type, '3-lh')).data
            stc = sensitivity_map(fwd, projs=None, mode='ratio',
                                  ch_type=ch_type, exclude='bads')
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == 'eeg':
            # radiality (4), angle (5), remaining (6), and  dampening (7)
            modes = ['radiality', 'angle', 'remaining', 'dampening']
            ends = ['4-lh', '5-lh', '6-lh', '7-lh']
            for mode, end in zip(modes, ends):
                w = read_source_estimate(sensmap_fname % (ch_type, end)).data
                stc = sensitivity_map(fwd, projs=projs, mode=mode,
                                      ch_type=ch_type, exclude='bads')
                assert_array_almost_equal(stc.data, w, decim)

    # test corner case for EEG
    stc = sensitivity_map(fwd, projs=[make_eeg_average_ref_proj(fwd['info'])],
                          ch_type='eeg', exclude='bads')
    # test corner case for projs being passed but no valid ones (#3135)
    assert_raises(ValueError, sensitivity_map, fwd, projs=None, mode='angle')
    assert_raises(RuntimeError, sensitivity_map, fwd, projs=[], mode='angle')
    # test volume source space
    fname = op.join(sample_path, 'sample_audvis_trunc-meg-vol-7-fwd.fif')
    fwd = mne.read_forward_solution(fname)
    sensitivity_map(fwd)

def plot_sensitivity_map(fwd_sol,
                         subject,
                         fname_leadfield_plot,
                         fname_sensitvity_plot):
    """Estimates and plots sensitivity map of forward solution."""

    # estimate lead field
    leadfield = fwd_sol['sol']['data']

    pp = PdfPages(fname_leadfield_plot)

    # plot leadfield
    plt.matshow(leadfield[:, :500])
    plt.xlabel('sources')
    plt.ylabel('sensors')
    plt.title('Lead field matrix (500 dipoles only)')
    pp.savefig()
    plt.close()


    # estimate sensitivity map for magnetometer
    mag_map = mne.sensitivity_map(fwd_sol, ch_type='mag', mode='fixed')

    # plot histogram of sensitivity
    plt.hist(mag_map.data.ravel(),
                   bins=20,
                   label='Magnetometers')
    plt.legend()
    plt.title('Normal orientation sensitivity')
    pp.savefig()
    plt.close()
    pp.close()

    subjects_dir = os.environ.get('SUBJECTS_DIR')
    brain = mag_map.plot(subject=subject,
                         time_label='Magnetometer sensitivity',
                         subjects_dir=subjects_dir,
                         fmin=0.1,
                         fmid=0.5,
                         fmax=0.9,
                         smoothing_steps=7)

    brain.save_montage(fname_sensitvity_plot)
    brain.close()

def test_sensitivity_maps():
    """Test sensitivity map computation."""
    fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        projs = read_proj(eog_fname)
        projs.extend(read_proj(ecg_fname))
    decim = 6
    for ch_type in ["eeg", "grad", "mag"]:
        w = read_source_estimate(sensmap_fname % (ch_type, "lh")).data
        stc = sensitivity_map(fwd, projs=None, ch_type=ch_type, mode="free", exclude="bads")
        assert_array_almost_equal(stc.data, w, decim)
        assert_true(stc.subject == "sample")
        # let's just make sure the others run
        if ch_type == "grad":
            # fixed (2)
            w = read_source_estimate(sensmap_fname % (ch_type, "2-lh")).data
            stc = sensitivity_map(fwd, projs=None, mode="fixed", ch_type=ch_type, exclude="bads")
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == "mag":
            # ratio (3)
            w = read_source_estimate(sensmap_fname % (ch_type, "3-lh")).data
            stc = sensitivity_map(fwd, projs=None, mode="ratio", ch_type=ch_type, exclude="bads")
            assert_array_almost_equal(stc.data, w, decim)
        if ch_type == "eeg":
            # radiality (4), angle (5), remaining (6), and  dampening (7)
            modes = ["radiality", "angle", "remaining", "dampening"]
            ends = ["4-lh", "5-lh", "6-lh", "7-lh"]
            for mode, end in zip(modes, ends):
                w = read_source_estimate(sensmap_fname % (ch_type, end)).data
                stc = sensitivity_map(fwd, projs=projs, mode=mode, ch_type=ch_type, exclude="bads")
                assert_array_almost_equal(stc.data, w, decim)

    # test corner case for EEG
    stc = sensitivity_map(fwd, projs=[make_eeg_average_ref_proj(fwd["info"])], ch_type="eeg", exclude="bads")
    # test corner case for projs being passed but no valid ones (#3135)
    assert_raises(ValueError, sensitivity_map, fwd, projs=None, mode="angle")
    assert_raises(RuntimeError, sensitivity_map, fwd, projs=[], mode="angle")
    # test volume source space
    fname = op.join(sample_path, "sample_audvis_trunc-meg-vol-7-fwd.fif")
    fwd = mne.read_forward_solution(fname)
    sensitivity_map(fwd)

bem = subjects_dir + '0001/bem/0001-inner_skull-bem-sol.fif'


# Note that forward solutions can also be read with read_forward_solution
fwd = mne.make_forward_solution(raw_fname, trans, src, bem,
                                fname="0001-fwd.fif", meg=True, eeg=False,
                                mindist=5.0,
                                n_jobs=n_jobs, overwrite=True)

# convert to surface orientation for better visualization
fwd = mne.convert_forward_solution(fwd, surf_ori=True)
leadfield = fwd['sol']['data']

print("Leadfield size : %d x %d" % leadfield.shape)

grad_map = mne.sensitivity_map(fwd, ch_type='grad', mode='fixed')
mag_map = mne.sensitivity_map(fwd, ch_type='mag', mode='fixed')

###############################################################################
# Show gain matrix a.k.a. leadfield matrix with sensitivity map

picks_meg = mne.pick_types(fwd['info'], meg=True, eeg=False)

fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
fig.suptitle('Lead field matrix (500 dipoles only)', fontsize=14)
for ax, picks, ch_type in zip(axes, [picks_meg], ['meg']):
    im = ax.imshow(leadfield[picks, :], origin='lower', aspect='auto',
                   cmap='RdBu_r')
    ax.set_title(ch_type.upper())
    ax.set_xlabel('sources')
    ax.set_ylabel('sensors')

bems = mne.make_bem_model(subject, conductivity=(0.3,),
                          subjects_dir=subjects_dir,
                          ico=None)  # ico = None for morphed SP.
bem_sol = mne.make_bem_solution(bems)
bem_sol['surfs'][0]['coord_frame'] = 5

##############################################################################
# Now we can read the channels that we want to map to the cortical locations.
# Then we can compute the forward solution.

info = hcp.read_info(subject=subject, hcp_path=hcp_path, data_type='rest',
                     run_index=0)

picks = mne.pick_types(info, meg=True, ref_meg=False)
info = mne.pick_info(info, picks)

fwd = mne.make_forward_solution(info, trans=head_mri_t, bem=bem_sol,
                                src=src_subject)
mag_map = mne.sensitivity_map(
    fwd, projs=None, ch_type='mag', mode='fixed', exclude=[], verbose=None)

##############################################################################
# we display sensitivity map on the original surface with little smoothing
# and admire the expected curvature-driven sensitivity pattern.

mag_map = mag_map.to_original_src(src_fsaverage, subjects_dir=subjects_dir)
mag_map.plot(subject='fsaverage', subjects_dir=subjects_dir,
             clim=dict(kind='percent', lims=[0, 50, 99]),
             smoothing_steps=2)

print __doc__

import mne
from mne.datasets import sample

data_path = sample.data_path()

fname = data_path + "/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif"
subjects_dir = data_path + "/subjects"

fwd = mne.read_forward_solution(fname, surf_ori=True)
leadfield = fwd["sol"]["data"]

print "Leadfield size : %d x %d" % leadfield.shape

grad_map = mne.sensitivity_map(fwd, ch_type="grad", mode="fixed")
mag_map = mne.sensitivity_map(fwd, ch_type="mag", mode="fixed")
eeg_map = mne.sensitivity_map(fwd, ch_type="eeg", mode="fixed")

###############################################################################
# Show gain matrix a.k.a. leadfield matrix with sensitivy map

import pylab as pl

pl.matshow(leadfield[:, :500])
pl.xlabel("sources")
pl.ylabel("sensors")
pl.title("Lead field matrix (500 dipoles only)")

pl.figure()
pl.hist(

import matplotlib.pyplot as plt

from mne import read_forward_solution, read_proj, sensitivity_map
from mne.datasets import sample

print(__doc__)

data_path = sample.data_path()

subjects_dir = data_path + '/subjects'
fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
ecg_fname = data_path + '/MEG/sample/sample_audvis_ecg-proj.fif'

fwd = read_forward_solution(fname, surf_ori=True)
projs = read_proj(ecg_fname)
projs = projs[3:][::2]  # take only one projection per channel type

# Compute sensitivity map
ssp_ecg_map = sensitivity_map(fwd, ch_type='grad', projs=projs, mode='angle')

###############################################################################
# Show sensitivity map

plt.hist(ssp_ecg_map.data.ravel())
plt.show()

args = dict(clim=dict(kind='value', lims=(0.2, 0.6, 1.)), smoothing_steps=7,
            hemi='rh', subjects_dir=subjects_dir)
ssp_ecg_map.plot(subject='sample', time_label='ECG SSP sensitivity', **args)

# Get relevant channel information
info = mne.io.read_info(raw_fname)
info = mne.pick_info(info, mne.pick_types(info, meg=True, eeg=False,
                                          exclude=[]))

# Morph fsaverage's source space to sample
src_fs = mne.read_source_spaces(fname_src_fs)
src_morph = mne.morph_source_spaces(src_fs, subject_to='sample',
                                    subjects_dir=subjects_dir)

# Compute the forward with our morphed source space
fwd = mne.make_forward_solution(info, trans=fname_trans,
                                src=src_morph, bem=fname_bem)
# fwd = mne.convert_forward_solution(fwd, surf_ori=True, force_fixed=True)
mag_map = mne.sensitivity_map(fwd, ch_type='mag')

# Return this SourceEstimate (on sample's surfaces) to fsaverage's surfaces
mag_map_fs = mag_map.to_original_src(src_fs, subjects_dir=subjects_dir)

# Plot the result, which tracks the sulcal-gyral folding
# outliers may occur, we'll place the cutoff at 99 percent.
kwargs = dict(clim=dict(kind='percent', lims=[0, 50, 99]),
              # no smoothing, let's see the dipoles on the cortex.
              smoothing_steps=1, hemi='rh', views=['lat'])

# Now note that the dipoles on fsaverage are almost equidistant while
# morphing will distribute the dipoles unevenly across the given subject's
# cortical surface to achieve the closest approximation to the average brain.
# Our testing code suggests a correlation of higher than 0.99.

#.........这里部分代码省略.........
    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',
               'sample-5120-5120-5120-bem-sol.fif')
    fname = join(meg_dir, 'sample_audvis-eeg-oct-6-fwd.fif')
    fwd_eeg = mne.make_forward_solution(raw.info, trans, src, bem,
                                        fname=fname, meg=False, eeg=True,
                                        mindist=5.0, n_jobs=2, overwrite=True)

    # for both EEG and MEG
    fname = join(meg_dir, 'sample_audvis-meg-eeg-oct-6-fwd.fif')
    fwd = mne.make_forward_solution(raw.info, trans, src, bem,
                                    fname=fname, meg=True, eeg=True,
                                    mindist=5.0, n_jobs=2, overwrite=True)

    # Create various sensitivity maps
    grad_map = mne.sensitivity_map(fwd, ch_type='grad', mode='free')
    grad_map.save(join(meg_dir, 'sample_audvis-grad-oct-6-fwd-sensmap'),
                  ftype='w')
    mag_map = mne.sensitivity_map(fwd, ch_type='mag', mode='free')
    mag_map.save(join(meg_dir, 'sample_audvis-mag-oct-6-fwd-sensmap'),
                 ftype='w')
    eeg_map = mne.sensitivity_map(fwd, ch_type='eeg', mode='free')
    eeg_map.save(join(meg_dir, 'sample_audvis-eeg-oct-6-fwd-sensmap'),
                 ftype='w')
    grad_map2 = mne.sensitivity_map(fwd, ch_type='grad', mode='fixed')
    grad_map2.save(join(meg_dir, 'sample_audvis-grad-oct-6-fwd-sensmap-2'),
                   ftype='w')
    mag_map2 = mne.sensitivity_map(fwd, ch_type='mag', mode='ratio')
    mag_map2.save(join(meg_dir, 'sample_audvis-mag-oct-6-fwd-sensmap-3'),
                  ftype='w')

    # Compute some with the EOG + ECG projectors
    projs = ecg_proj + eog_proj + raw.info['projs']
    for map_type in ['radiality', 'angle', 'remaining', 'dampening']:
        eeg_map = mne.sensitivity_map(fwd, projs=projs, ch_type='eeg',
                                      mode=map_type)
        eeg_map.save(join(meg_dir,
                          'sample_audvis-eeg-oct-6-fwd-sensmap-' + map_type))

###############################################################################
    # Compute MNE inverse operators
    #
    # Note: The MEG/EEG forward solution could be used for all
    #
    inv_meg = make_inverse_operator(raw.info, fwd_meg, noise_cov, loose=0.2)
    fname = join(meg_dir, 'sample_audvis-meg-oct-6-meg-inv.fif')
    write_inverse_operator(fname, inv_meg)

raw = mne.io.read_raw_edf(raw_fname, preload=True)


# Clean channel names to be able to use a standard 1005 montage
ch_names = [c.replace('.', '') for c in raw.ch_names]
raw.rename_channels({old: new for old, new in zip(raw.ch_names, ch_names)})

# Read and set the EEG electrode locations
montage = mne.channels.read_montage('standard_1005', ch_names=raw.ch_names,
                                    transform=True)

raw.set_montage(montage)
raw.set_eeg_reference(projection=True)  # needed for inverse modeling

# Check that the locations of EEG electrodes is correct with respect to MRI
mne.viz.plot_alignment(
    raw.info, src=src, eeg=['original', 'projected'], trans=trans, dig=True)

##############################################################################
# Setup source space and compute forward
# --------------------------------------

fwd = mne.make_forward_solution(raw.info, trans=trans, src=src,
                                bem=bem, eeg=True, mindist=5.0, n_jobs=1)
print(fwd)

# for illustration purposes use fwd to compute the sensitivity map
eeg_map = mne.sensitivity_map(fwd, ch_type='eeg', mode='fixed')
eeg_map.plot(time_label='EEG sensitivity', subjects_dir=subjects_dir,
             clim=dict(lims=[5, 50, 100]))