def test_evoked_arithmetic():
    """Test evoked arithmetic
    """
    ev = read_evokeds(fname, condition=0)
    ev1 = EvokedArray(np.ones_like(ev.data), ev.info, ev.times[0], nave=20)
    ev2 = EvokedArray(-np.ones_like(ev.data), ev.info, ev.times[0], nave=10)

    # combine_evoked([ev1, ev2]) should be the same as ev1 + ev2:
    # data should be added according to their `nave` weights
    # nave = ev1.nave + ev2.nave
    ev = ev1 + ev2
    assert_equal(ev.nave, ev1.nave + ev2.nave)
    assert_allclose(ev.data, 1. / 3. * np.ones_like(ev.data))
    ev = ev1 - ev2
    assert_equal(ev.nave, ev1.nave + ev2.nave)
    assert_equal(ev.comment, ev1.comment + ' - ' + ev2.comment)
    assert_allclose(ev.data, np.ones_like(ev1.data))
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        ev = merge_evoked([ev1, ev2])
    assert_true(len(w) >= 1)
    assert_allclose(ev.data, 1. / 3. * np.ones_like(ev.data))

    # default comment behavior if evoked.comment is None
    old_comment1 = ev1.comment
    old_comment2 = ev2.comment
    ev1.comment = None
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        ev = ev1 - ev2
        assert_equal(ev.comment, 'unknown')
    ev1.comment = old_comment1
    ev2.comment = old_comment2

    # equal weighting
    ev = combine_evoked([ev1, ev2], weights='equal')
    assert_allclose(ev.data, np.zeros_like(ev1.data))

    # combine_evoked([ev1, ev2], weights=[1, 0]) should yield the same as ev1
    ev = combine_evoked([ev1, ev2], weights=[1, 0])
    assert_equal(ev.nave, ev1.nave)
    assert_allclose(ev.data, ev1.data)

    # simple subtraction (like in oddball)
    ev = combine_evoked([ev1, ev2], weights=[1, -1])
    assert_allclose(ev.data, 2 * np.ones_like(ev1.data))

    assert_raises(ValueError, combine_evoked, [ev1, ev2], weights='foo')
    assert_raises(ValueError, combine_evoked, [ev1, ev2], weights=[1])

    # grand average
    evoked1, evoked2 = read_evokeds(fname, condition=[0, 1], proj=True)
    ch_names = evoked1.ch_names[2:]
    evoked1.info['bads'] = ['EEG 008']  # test interpolation
    evoked1.drop_channels(evoked1.ch_names[:1])
    evoked2.drop_channels(evoked2.ch_names[1:2])
    gave = grand_average([evoked1, evoked2])
    assert_equal(gave.data.shape, [len(ch_names), evoked1.data.shape[1]])
    assert_equal(ch_names, gave.ch_names)
    assert_equal(gave.nave, 2)

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 test_make_inverse_operator_vector():
    """Test MNE inverse computation (vector result)."""
    fwd_surf = read_forward_solution_meg(fname_fwd, surf_ori=True)
    fwd_fixed = read_forward_solution_meg(fname_fwd, surf_ori=False)
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    # Make different version of the inverse operator
    inv_1 = make_inverse_operator(evoked.info, fwd_fixed, noise_cov, loose=1)
    inv_2 = make_inverse_operator(evoked.info, fwd_surf, noise_cov, depth=None,
                                  use_cps=True)
    inv_3 = make_inverse_operator(evoked.info, fwd_surf, noise_cov, fixed=True,
                                  use_cps=True)
    inv_4 = make_inverse_operator(evoked.info, fwd_fixed, noise_cov,
                                  loose=.2, depth=None)

    # Apply the inverse operators and check the result
    for ii, inv in enumerate((inv_1, inv_2, inv_4)):
        # Don't do eLORETA here as it will be quite slow
        methods = ['MNE', 'dSPM', 'sLORETA'] if ii < 2 else ['MNE']
        for method in methods:
            stc = apply_inverse(evoked, inv, method=method)
            stc_vec = apply_inverse(evoked, inv, pick_ori='vector',
                                    method=method)
            assert_allclose(stc.data, stc_vec.magnitude().data)

    # Vector estimates don't work when using fixed orientations
    pytest.raises(RuntimeError, apply_inverse, evoked, inv_3,
                  pick_ori='vector')

    # When computing with vector fields, computing the difference between two
    # evokeds and then performing the inverse should yield the same result as
    # computing the difference between the inverses.
    evoked0 = read_evokeds(fname_data, condition=0, baseline=(None, 0))
    evoked0.crop(0, 0.2)
    evoked1 = read_evokeds(fname_data, condition=1, baseline=(None, 0))
    evoked1.crop(0, 0.2)
    diff = combine_evoked((evoked0, evoked1), [1, -1])
    stc_diff = apply_inverse(diff, inv_1, method='MNE')
    stc_diff_vec = apply_inverse(diff, inv_1, method='MNE', pick_ori='vector')
    stc_vec0 = apply_inverse(evoked0, inv_1, method='MNE', pick_ori='vector')
    stc_vec1 = apply_inverse(evoked1, inv_1, method='MNE', pick_ori='vector')
    assert_allclose(stc_diff_vec.data, (stc_vec0 - stc_vec1).data,
                    atol=1e-20)
    assert_allclose(stc_diff.data, (stc_vec0 - stc_vec1).magnitude().data,
                    atol=1e-20)

tfce = dict(start=.2, step=.2)

t_obs, clusters, cluster_pv, h0 = spatio_temporal_cluster_test(
    X, tfce, n_permutations=100)  # a more standard number would be 1000+
significant_points = cluster_pv.reshape(t_obs.shape).T < .05
print(str(significant_points.sum()) + " points selected by TFCE ...")

##############################################################################
# The results of these mass univariate analyses can be visualised by plotting
# :class:`mne.Evoked` objects as images (via :class:`mne.Evoked.plot_image`)
# and masking points for significance.
# Here, we group channels by Regions of Interest to facilitate localising
# effects on the head.

# We need an evoked object to plot the image to be masked
evoked = mne.combine_evoked([long_words.average(), -short_words.average()],
                            weights='equal')  # calculate difference wave
time_unit = dict(time_unit="s")
evoked.plot_joint(title="Long vs. short words", ts_args=time_unit,
                  topomap_args=time_unit)  # show difference wave

# Create ROIs by checking channel labels
selections = make_1020_channel_selections(evoked.info, midline="12z")

# Visualize the results
fig, axes = plt.subplots(nrows=3, figsize=(8, 8))
axes = {sel: ax for sel, ax in zip(selections, axes.ravel())}
evoked.plot_image(axes=axes, group_by=selections, colorbar=False, show=False,
                  mask=significant_points, show_names="all", titles=None,
                  **time_unit)
plt.colorbar(axes["Left"].images[-1], ax=list(axes.values()), shrink=.3,
             label="uV")

def plot_artefact_overview(raw_orig, raw_clean, stim_event_ids=[1],
                           stim_ch='STI 014', resp_ch=None,
                           resp_event_ids=None,
                           ecg_ch='EEG 002',
                           eog1_ch='EEG 001', eog2_ch='EEG 003',
                           eog_tmin=-0.5, eog_tmax=0.5, eog_id=998,
                           eog_lfreq=8., eog_hfreq=20.,
                           ecg_tmin=-0.5, ecg_tmax=0.5, ecg_id=999,
                           ecg_lfreq=8., ecg_hfreq=20.,
                           stim_tmin=-0.2, stim_tmax=0.8,
                           resp_tmin=-0.6, resp_tmax=0.4,
                           eve_output='onset', overview_fname=None):
    '''
    Plot an overview of the artefact rejection with ECG, EOG vertical and EOG
    horizontal channels. Shows the data before and after cleaning along with a
    difference plot.

    raw_orig: instance of mne.io.Raw | str
        File name of raw object of the uncleaned data.
    raw_clean: instance of mne.io.Raw | str
        File name of raw object of the cleaned data.
    stim_event_ids: list
        List of stim or resp event ids. Defaults to [1].
    resp_event_ids: list
        List of stim or resp event ids. Defaults to None.
    eve_output: 'onset' | 'offset' | 'step'
        Whether to report when events start, when events end, or both.
    overview_fname: str | None
        Name to save the plot generated. (considers raw_clean.info['filename'])

    Notes: Time is always shown in milliseconds (1e3) and the MEG data from mag
        is always in femtoTesla (fT) (1e15)
    '''

    from mne.preprocessing import create_ecg_epochs, create_eog_epochs

    raw = check_read_raw(raw_orig, preload=True)
    raw_clean = check_read_raw(raw_clean, preload=True)

    if not overview_fname:
        try:
            overview_fname = raw_clean.info['filename'].rsplit('-raw.fif')[0] + ',overview-plot.png'
        except:
            overview_fname = 'overview-plot.png'

    # stim related events
    events = mne.find_events(raw, stim_channel=stim_ch, output='onset')
    events_clean = mne.find_events(raw_clean, stim_channel=stim_ch, output='onset')

    epochs = mne.Epochs(raw, events, event_id=stim_event_ids,
                        tmin=stim_tmin, tmax=stim_tmax,
                        picks=mne.pick_types(raw.info, meg=True, exclude='bads'))
    evoked = epochs.average()
    epochs_clean = mne.Epochs(raw_clean, events_clean, event_id=stim_event_ids,
                              tmin=stim_tmin, tmax=stim_tmax,
                              picks=mne.pick_types(raw_clean.info, meg=True, exclude='bads'))
    evoked_clean = epochs_clean.average()

    stim_diff_signal = mne.combine_evoked([evoked, evoked_clean],
                                          weights=[1, -1])

    if resp_ch:
        # stim related events
        resp_events = mne.find_events(raw, stim_channel=resp_ch, output='onset')
        resp_events_clean = mne.find_events(raw_clean, stim_channel=resp_ch, output='onset')

        resp_epochs = mne.Epochs(raw, resp_events, event_id=resp_event_ids,
                            tmin=resp_tmin, tmax=resp_tmax,
                            picks=mne.pick_types(raw.info, meg=True, exclude='bads'))
        resp_evoked = resp_epochs.average()
        resp_epochs_clean = mne.Epochs(raw_clean, resp_events_clean, event_id=resp_event_ids,
                                  tmin=resp_tmin, tmax=resp_tmax,
                                  picks=mne.pick_types(raw_clean.info, meg=True, exclude='bads'))
        resp_evoked_clean = resp_epochs_clean.average()

        resp_diff_signal = mne.combine_evoked([resp_evoked, resp_evoked_clean],
                                              weights=[1, -1])

    # MEG signal around ECG events
    ecg_epochs = create_ecg_epochs(raw, ch_name=ecg_ch, event_id=ecg_id,
                                   picks=mne.pick_types(raw.info, meg=True, ecg=True, exclude=[ecg_ch]),
                                   tmin=ecg_tmin, tmax=ecg_tmax,
                                   l_freq=ecg_lfreq, h_freq=ecg_hfreq,
                                   preload=True, keep_ecg=False, baseline=(None, None))

    ecg_clean_epochs = create_ecg_epochs(raw_clean, ch_name=ecg_ch, event_id=ecg_id,
                                         picks=mne.pick_types(raw.info, meg=True, ecg=True, exclude=[ecg_ch]),
                                         tmin=ecg_tmin, tmax=ecg_tmax,
                                         l_freq=ecg_lfreq, h_freq=ecg_hfreq,
                                         preload=True, keep_ecg=False, baseline=(None, None))

    stim_diff_ecg = mne.combine_evoked([ecg_epochs.average(), ecg_clean_epochs.average()],
                                       weights=[1, -1])

    # MEG signal around EOG1 events
    eog1_epochs = create_eog_epochs(raw, ch_name=eog1_ch, event_id=eog_id,
                                    picks=mne.pick_types(raw.info, meg=True, exclude='bads'),
                                    tmin=eog_tmin, tmax=eog_tmax,
                                    l_freq=eog_lfreq, h_freq=eog_hfreq,
                                    preload=True, baseline=(None, None))
#.........这里部分代码省略.........

# painting an area with clicking and holding the left mouse button.
evoked_std.plot(window_title='Standard', gfp=True, time_unit='s')
evoked_dev.plot(window_title='Deviant', gfp=True, time_unit='s')


###############################################################################
# Show activations as topography figures.
times = np.arange(0.05, 0.301, 0.025)
evoked_std.plot_topomap(times=times, title='Standard', time_unit='s')
evoked_dev.plot_topomap(times=times, title='Deviant', time_unit='s')

###############################################################################
# We can see the MMN effect more clearly by looking at the difference between
# the two conditions. P50 and N100 are no longer visible, but MMN/P200 and
# P300 are emphasised.
evoked_difference = combine_evoked([evoked_dev, -evoked_std], weights='equal')
evoked_difference.plot(window_title='Difference', gfp=True, time_unit='s')

###############################################################################
# Source estimation.
# We compute the noise covariance matrix from the empty room measurement
# and use it for the other runs.
reject = dict(mag=4e-12)
cov = mne.compute_raw_covariance(raw_erm, reject=reject)
cov.plot(raw_erm.info)
del raw_erm

###############################################################################
# The transformation is read from a file. More information about coregistering
# the data, see :ref:`ch_interactive_analysis` or
# :func:`mne.gui.coregistration`.

event_id = {"left/auditory": 1, "right/auditory": 2, "left/visual": 3, "right/visual": 4}
epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, reject=reject)
epochs = mne.Epochs(raw, **epochs_params)

print(epochs)

###############################################################################
# Next, we create averages of stimulation-left vs stimulation-right trials.
# We can use basic arithmetic to, for example, construct and plot
# difference ERPs.

left, right = epochs["left"].average(), epochs["right"].average()

# create and plot difference ERP
mne.combine_evoked([left, -right], weights="equal").plot_joint()

###############################################################################
# This is an equal-weighting difference. If you have imbalanced trial numbers,
# you could also consider either equalizing the number of events per
# condition (using
# :meth:`epochs.equalize_epochs_counts <mne.Epochs.equalize_event_counts`).
# As an example, first, we create individual ERPs for each condition.

aud_l = epochs["auditory", "left"].average()
aud_r = epochs["auditory", "right"].average()
vis_l = epochs["visual", "left"].average()
vis_r = epochs["visual", "right"].average()

all_evokeds = [aud_l, aud_r, vis_l, vis_r]
print(all_evokeds)

# Fit ICA, find and remove major artifacts
ica = ICA(n_components=0.95, random_state=0).fit(raw, decim=1, reject=reject)

# compute correlation scores, get bad indices sorted by score
eog_epochs = create_eog_epochs(raw, ch_name='MRT31-2908', reject=reject)
eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name='MRT31-2908')
ica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on
ica.plot_components(eog_inds)  # view topographic sensitivity of components
ica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar
ica.plot_overlay(eog_epochs.average())  # inspect artifact removal
ica.apply(epochs)  # clean data, default in place

evoked = [epochs[k].average() for k in event_ids]

contrast = combine_evoked(evoked, weights=[-1, 1])  # Faces - scrambled

evoked.append(contrast)

for e in evoked:
    e.plot(ylim=dict(mag=[-400, 400]))

plt.show()

# estimate noise covarariance
noise_cov = mne.compute_covariance(epochs, tmax=0, method='shrunk',
                                   rank=None)

###############################################################################
# Visualize fields on MEG helmet

# :meth:`~mne.Evoked.plot_topomap`. Here we'll examine just the EEG channels,
# and see the classic auditory evoked N100-P200 pattern over dorso-frontal
# electrodes, then plot scalp topographies at some additional arbitrary times:

# sphinx_gallery_thumbnail_number = 13
aud_evoked.plot_joint(picks='eeg')
aud_evoked.plot_topomap(times=[0., 0.08, 0.1, 0.12, 0.2], ch_type='eeg')

##############################################################################
# Evoked objects can also be combined to show contrasts between conditions,
# using the :func:`mne.combine_evoked` function. A simple difference can be
# generated by negating one of the :class:`~mne.Evoked` objects passed into the
# function. We'll then plot the difference wave at each sensor using
# :meth:`~mne.Evoked.plot_topo`:

evoked_diff = mne.combine_evoked([aud_evoked, -vis_evoked], weights='equal')
evoked_diff.pick_types('mag').plot_topo(color='r', legend=False)

##############################################################################
# Inverse modeling
# ^^^^^^^^^^^^^^^^
#
# Finally, we can estimate the origins of the evoked activity by projecting the
# sensor data into this subject's :term:`source space` (a set of points either
# on the cortical surface or within the cortical volume of that subject, as
# estimated by structural MRI scans). MNE-Python supports lots of ways of doing
# this (dynamic statistical parametric mapping, dipole fitting, beamformers,
# etc.); here we'll use minimum-norm estimation (MNE) to generate a continuous
# map of activation constrained to the cortical surface. MNE uses a linear
# :term:`inverse operator` to project EEG+MEG sensor measurements into the
# source space. The inverse operator is computed from the

evoked_condition_1 = epochs_condition_1.average()
evoked_condition_2 = epochs_condition_2.average()

plt.figure()
plt.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)

plt.subplot(2, 1, 1)
# Create new stats image with only significant clusters
T_obs_plot = np.nan * np.ones_like(T_obs)
for c, p_val in zip(clusters, cluster_p_values):
    if p_val <= 0.05:
        T_obs_plot[c] = T_obs[c]

plt.imshow(T_obs,
           extent=[times[0], times[-1], freqs[0], freqs[-1]],
           aspect='auto', origin='lower', cmap='gray')
plt.imshow(T_obs_plot,
           extent=[times[0], times[-1], freqs[0], freqs[-1]],
           aspect='auto', origin='lower', cmap='RdBu_r')

plt.xlabel('Time (ms)')
plt.ylabel('Frequency (Hz)')
plt.title('Induced power (%s)' % ch_name)

ax2 = plt.subplot(2, 1, 2)
evoked_contrast = mne.combine_evoked([evoked_condition_1, evoked_condition_2],
                                     weights=[1, -1])
evoked_contrast.plot(axes=ax2, time_unit='s')

plt.show()

evoked_hcp = None
hcp_evokeds = hcp.read_evokeds(onset='stim', **hcp_params)

for ev in hcp_evokeds:
    if not ev.comment == 'Wrkmem_LM-TIM-face_BT-diff_MODE-mag':
        continue

# Once more we add and interpolate missing channels
evoked_hcp = preproc.interpolate_missing(ev, **hcp_params)


##############################################################################
# Time to compare the outputs
#

evoked = mne.combine_evoked(evokeds, weights='equal')
evoked_from_epochs_hcp = mne.combine_evoked(
    evokeds_from_epochs_hcp, weights='equal')

fig1, axes = plt.subplots(3, 1, figsize=(12, 8))

evoked.plot(axes=axes[0], show=False)
axes[0].set_title('MNE-HCP')

evoked_from_epochs_hcp.plot(axes=axes[1], show=False)
axes[1].set_title('HCP epochs')

evoked_hcp.plot(axes=axes[2], show=False)
axes[2].set_title('HCP evoked')
fig1.canvas.draw()

# Set up events
events = mne.find_events(raw)
event_id = {'Aud/L': 1, 'Aud/R': 2}
tmin, tmax = -.1, .5

# regular epoching
picks = mne.pick_types(raw.info, meg=True)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, reject=None,
                    baseline=None, preload=True, verbose=False)

# rERF
evokeds = linear_regression_raw(raw, events=events, event_id=event_id,
                                reject=None, tmin=tmin, tmax=tmax)
# linear_regression_raw returns a dict of evokeds
# select conditions similarly to mne.Epochs objects

# plot both results, and their difference
cond = "Aud/L"
fig, (ax1, ax2, ax3) = plt.subplots(3, 1)
params = dict(spatial_colors=True, show=False, ylim=dict(grad=(-200, 200)),
              time_unit='s')
epochs[cond].average().plot(axes=ax1, **params)
evokeds[cond].plot(axes=ax2, **params)
contrast = mne.combine_evoked([evokeds[cond], -epochs[cond].average()],
                              weights='equal')
contrast.plot(axes=ax3, **params)
ax1.set_title("Traditional averaging")
ax2.set_title("rERF")
ax3.set_title("Difference")
plt.show()

    ##############################################################################################
    ####                                    By Condition                                      ####
    ##############################################################################################


    ####################################### Tactile/Tactile ######################################

    evoked_left_cued_TT = get_topo( raw, {'LT/LT': 25}, np.arange(.060, .120, .010), 0.005, tmin, tmax, reject_num = 100e-6 )
    evoked_right_cued_TT = get_topo(raw, {'RT/RT': 35}, np.arange(.060, .120, .010), 0.005, tmin, tmax, reject_num=100e-6)

    evoked_left_uncued_TT = get_topo( raw, {'RT/LT': 33}, np.arange(.060, .120, .010), 0.005, tmin, tmax, reject_num = 100e-6 )
    evoked_right_uncued_TT = get_topo( raw, {'LT/RT': 27}, np.arange(.060, .120, .010), 0.005, tmin, tmax, reject_num = 100e-6 )

    # get individual effect plots
    ind_left_red_TT = mne.combine_evoked([evoked_left_cued_TT, evoked_left_uncued_TT], weights=[1, -1])  # subtraction
    ind_left_red_TT.plot_topomap([.090], ch_type='eeg', show_names=True, average=0.030)
    fig = plt.gcf()
    fig.set_size_inches(18.5, 10.5)
    plt.savefig('%s/P_topo/%s_left_red_TT_topomap.png' % (directory, participant))
    if dont_plot:
        plt.close()

    # get individual effect plots
    ind_right_red_TT = mne.combine_evoked([evoked_right_cued_TT, evoked_right_uncued_TT],
                                         weights=[1, -1])  # subtraction
    ind_right_red_TT.plot_topomap([.090], ch_type='eeg', show_names=True, average=0.030)
    fig = plt.gcf()
    fig.set_size_inches(18.5, 10.5)
    plt.savefig('%s/P_topo/%s_right_red_TT_topomap.png' % (directory, participant))
    if dont_plot:

                               stim=True, exclude=bads)

        # create the real-time epochs object and start acquisition
        rt_epochs = RtEpochs(rt_client, event_id, tmin, tmax,
                             stim_channel='STI 014', picks=picks,
                             reject=dict(grad=4000e-13, eog=150e-6),
                             decim=1, isi_max=2.0, proj=None)
        rt_epochs.start()
        for ii, ev in enumerate(rt_epochs.iter_evoked()):
            print("Just got epoch %d" % (ii + 1))

            ev.pick_types(meg=True, eog=False)
            if ii == 0:
                evoked = ev
            else:
                evoked = mne.combine_evoked([evoked, ev], weights='nave')

            ax[0].cla()
            ax[1].cla()  # clear axis

            plot_events(rt_epochs.events[-5:], sfreq=ev.info['sfreq'],
                        first_samp=-rt_client.tmin_samp, axes=ax[0])

            # plot on second subplot
            evoked.plot(axes=ax[1], selectable=False, time_unit='s')
            ax[1].set_title('Evoked response for gradiometer channels'
                            '(event_id = %d)' % event_id)

            plt.pause(0.05 / speedup)
            plt.draw()
        rt_epochs.stop()

# condition (using ``Epochs.equalize_event_counts``), or the ``combine_evoked``
# function.
# As an example, first, we create individual ERPs for each condition.

aud_l = epochs["auditory", "left"].average()
aud_r = epochs["auditory", "right"].average()
vis_l = epochs["visual", "left"].average()
vis_r = epochs["visual", "right"].average()

all_evokeds = [aud_l, aud_r, vis_l, vis_r]

# This could have been much simplified with a list comprehension:
# all_evokeds = [epochs[cond] for cond in event_id]

# Then, we construct and plot an unweighted average of left vs. right trials.
mne.combine_evoked(all_evokeds, weights=(1, -1, 1, -1)).plot_joint()

###############################################################################
# Often, it makes sense to store Evoked objects in a dictionary or a list -
# either different conditions, or different subjects.

# If they are stored in a list, they can be easily averaged, for example,
# for a grand average across subjects (or conditions).
grand_average = mne.grand_average(all_evokeds)
mne.write_evokeds('/tmp/tmp-ave.fif', all_evokeds)

# If Evokeds objects are stored in a dictionary, they can be retrieved by name.
all_evokeds = dict((cond, epochs[cond].average()) for cond in event_id)
print(all_evokeds['left/auditory'])

# Besides for explicit access, this can be used for example to set titles.

    evoked_fname, condition='Left Auditory', baseline=(None, 0), proj=True)

##############################################################################
# Or another one stored in the same file:

evoked2 = mne.read_evokeds(
    evoked_fname, condition='Right Auditory', baseline=(None, 0), proj=True)

##############################################################################
# Two evoked objects can be contrasted using :func:`mne.combine_evoked`.
# This function can use ``weights='equal'``, which provides a simple
# element-by-element subtraction (and sets the
# ``mne.Evoked.nave`` attribute properly based on the underlying number
# of trials) using either equivalent call:

contrast = mne.combine_evoked([evoked1, evoked2], weights=[0.5, -0.5])
contrast = mne.combine_evoked([evoked1, -evoked2], weights='equal')
print(contrast)

##############################################################################
# To do a weighted sum based on the number of averages, which will give
# you what you would have gotten from pooling all trials together in
# :class:`mne.Epochs` before creating the :class:`mne.Evoked` instance,
# you can use ``weights='nave'``:

average = mne.combine_evoked([evoked1, evoked2], weights='nave')
print(contrast)

##############################################################################
# Instead of dealing with mismatches in the number of averages, we can use
# trial-count equalization before computing a contrast, which can have some

contrasts = list()

for subject_id in range(1, 20):
    if subject_id in exclude:
        continue
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)
    data_path = op.join(meg_dir, subject)
    contrast = mne.read_evokeds(op.join(data_path, '%s-ave.fif' % subject),
                                condition='contrast')
    contrast.pick_types(meg=False, eeg=True)
    contrast.apply_baseline((-0.2, 0.0))
    contrasts.append(contrast)

contrast = mne.combine_evoked(contrasts, 'equal')

channel = 'EEG070'
idx = contrast.ch_names.index(channel)
mne.viz.plot_compare_evokeds(contrast, [idx])

##############################################################################
# Assemble the data and run the cluster stats on channel data

data = np.array([c.data[idx] for c in contrasts])

n_permutations = 1000  # number of permutations to run

# set family-wise p-value
p_accept = 0.01

                     reject=reject)
epochs = mne.Epochs(raw, **epochs_params)

print(epochs)

###############################################################################
# Next, we create averages of stimulation-left vs stimulation-right trials.
# We can use basic arithmetic to, for example, construct and plot
# difference ERPs.

left, right = epochs["left"].average(), epochs["right"].average()

# create and plot difference ERP
joint_kwargs = dict(ts_args=dict(time_unit='s'),
                    topomap_args=dict(time_unit='s'))
mne.combine_evoked([left, -right], weights='equal').plot_joint(**joint_kwargs)

###############################################################################
# This is an equal-weighting difference. If you have imbalanced trial numbers,
# you could also consider either equalizing the number of events per
# condition (using
# :meth:`epochs.equalize_event_counts <mne.Epochs.equalize_event_counts>`).
# As an example, first, we create individual ERPs for each condition.

aud_l = epochs["auditory", "left"].average()
aud_r = epochs["auditory", "right"].average()
vis_l = epochs["visual", "left"].average()
vis_r = epochs["visual", "right"].average()

all_evokeds = [aud_l, aud_r, vis_l, vis_r]
print(all_evokeds)

def test_arithmetic():
    """Test evoked arithmetic."""
    ev = read_evokeds(fname, condition=0)
    ev1 = EvokedArray(np.ones_like(ev.data), ev.info, ev.times[0], nave=20)
    ev2 = EvokedArray(-np.ones_like(ev.data), ev.info, ev.times[0], nave=10)

    # combine_evoked([ev1, ev2]) should be the same as ev1 + ev2:
    # data should be added according to their `nave` weights
    # nave = ev1.nave + ev2.nave
    ev = combine_evoked([ev1, ev2], weights='nave')
    assert_equal(ev.nave, ev1.nave + ev2.nave)
    assert_allclose(ev.data, 1. / 3. * np.ones_like(ev.data))

    # with same trial counts, a bunch of things should be equivalent
    for weights in ('nave', 'equal', [0.5, 0.5]):
        ev = combine_evoked([ev1, ev1], weights=weights)
        assert_allclose(ev.data, ev1.data)
        assert_equal(ev.nave, 2 * ev1.nave)
        ev = combine_evoked([ev1, -ev1], weights=weights)
        assert_allclose(ev.data, 0., atol=1e-20)
        assert_equal(ev.nave, 2 * ev1.nave)
    ev = combine_evoked([ev1, -ev1], weights='equal')
    assert_allclose(ev.data, 0., atol=1e-20)
    assert_equal(ev.nave, 2 * ev1.nave)
    ev = combine_evoked([ev1, -ev2], weights='equal')
    expected = int(round(1. / (0.25 / ev1.nave + 0.25 / ev2.nave)))
    assert_equal(expected, 27)  # this is reasonable
    assert_equal(ev.nave, expected)

    # default comment behavior if evoked.comment is None
    old_comment1 = ev1.comment
    old_comment2 = ev2.comment
    ev1.comment = None
    ev = combine_evoked([ev1, -ev2], weights=[1, -1])
    assert_equal(ev.comment.count('unknown'), 2)
    assert_true('-unknown' in ev.comment)
    assert_true(' + ' in ev.comment)
    ev1.comment = old_comment1
    ev2.comment = old_comment2

    # equal weighting
    ev = combine_evoked([ev1, ev2], weights='equal')
    assert_allclose(ev.data, np.zeros_like(ev1.data))

    # combine_evoked([ev1, ev2], weights=[1, 0]) should yield the same as ev1
    ev = combine_evoked([ev1, ev2], weights=[1, 0])
    assert_equal(ev.nave, ev1.nave)
    assert_allclose(ev.data, ev1.data)

    # simple subtraction (like in oddball)
    ev = combine_evoked([ev1, ev2], weights=[1, -1])
    assert_allclose(ev.data, 2 * np.ones_like(ev1.data))

    assert_raises(ValueError, combine_evoked, [ev1, ev2], weights='foo')
    assert_raises(ValueError, combine_evoked, [ev1, ev2], weights=[1])

    # grand average
    evoked1, evoked2 = read_evokeds(fname, condition=[0, 1], proj=True)
    ch_names = evoked1.ch_names[2:]
    evoked1.info['bads'] = ['EEG 008']  # test interpolation
    evoked1.drop_channels(evoked1.ch_names[:1])
    evoked2.drop_channels(evoked2.ch_names[1:2])
    gave = grand_average([evoked1, evoked2])
    assert_equal(gave.data.shape, [len(ch_names), evoked1.data.shape[1]])
    assert_equal(ch_names, gave.ch_names)
    assert_equal(gave.nave, 2)
    assert_raises(ValueError, grand_average, [1, evoked1])

which can be useful for reproducing research results.

The MEGSIM files will be dowloaded automatically.

The datasets are documented in:
Aine CJ, Sanfratello L, Ranken D, Best E, MacArthur JA, Wallace T,
Gilliam K, Donahue CH, Montano R, Bryant JE, Scott A, Stephen JM
(2012) MEG-SIM: A Web Portal for Testing MEG Analysis Methods using
Realistic Simulated and Empirical Data. Neuroinformatics 10:141-158
"""

from mne import read_evokeds, combine_evoked
from mne.datasets.megsim import load_data

print(__doc__)

condition = 'visual'  # or 'auditory' or 'somatosensory'

# Load experimental RAW files for the visual condition
epochs_fnames = load_data(condition=condition, data_format='single-trial',
                          data_type='simulation', verbose=True)

# Take only 10 trials from the same simulation setup.
epochs_fnames = [f for f in epochs_fnames if 'sim6_trial_' in f][:10]

evokeds = [read_evokeds(f)[0] for f in epochs_fnames]
mean_evoked = combine_evoked(evokeds, weights='nave')

# Visualize the average
mean_evoked.plot()