<jupyter_start><jupyter_text>Import the packages we will need for the solution.<jupyter_code>import os
import numpy as np
import mne
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Load the sample<jupyter_code>sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_filt-0-40_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file)<jupyter_output><empty_output><jupyter_text>PreprocessingNote that the data we are loading is filtered from 0.1 to 40 Hz so we do not need to filter it (for now).<jupyter_code># set up and fit the ICA
ica = mne.preprocessing.ICA(n_components=20, random_state=97, max_iter=800)
ica.fit(raw)
ica.exclude = [1, 2]  # details on how we picked these are omitted here
ica.plot_properties(raw, picks=ica.exclude);
orig_raw = raw.copy()
raw.load_data()
ica.apply(raw)

# show some frontal channels to clearly illustrate the artifact removal
chs = ['MEG 0111', 'MEG 0121', 'MEG 0131', 'MEG 0211', 'MEG 0221', 'MEG 0231',
       'MEG 0311', 'MEG 0321', 'MEG 0331', 'MEG 1511', 'MEG 1521', 'MEG 1531',
       'EEG 001', 'EEG 002', 'EEG 003', 'EEG 004', 'EEG 005', 'EEG 006',
       'EEG 007', 'EEG 008']
chan_idxs = [raw.ch_names.index(ch) for ch in chs]<jupyter_output><empty_output><jupyter_text>Detecting experimental eventsThe sample dataset includes several :term:`"STIM" channels `that recorded electricalsignals sent from the stimulus delivery computer (as brief DC shifts /squarewave pulses). These pulses (often called "triggers") are used in thisdataset to mark experimental events: stimulus onset, stimulus type, andparticipant response (button press). The individual STIM channels arecombined onto a single channel, in such a way that voltagelevels on that channel can be unambiguously decoded as a particular eventtype. On older Neuromag systems (such as that used to record the sample data)this summation channel was called ``STI 014``, so we can pass that channelname to the :func:`mne.find_events` function to recover the timing andidentity of the stimulus events.<jupyter_code>events = mne.find_events(raw, stim_channel='STI 014')
print(events[:5])  # show the first 5
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4, 'smiley': 5, 'buttonpress': 32}<jupyter_output><empty_output><jupyter_text>For paradigms that are not event-related (e.g., analysis of resting-statedata), you can extract regularly spaced (possibly overlapping) spans of databy creating events using :func:`mne.make_fixed_length_events` and thenproceeding with epoching as described in the next section. Epoching continuous dataThe :class:`~mne.io.Raw` object and the events array are the bare minimumneeded to create an :class:`~mne.Epochs` object, which we create with the:class:`~mne.Epochs` class constructor. Here we'll also specify some dataquality constraints: we'll reject any epoch where peak-to-peak signalamplitude is beyond reasonable limits for that channel type. This is donewith a *rejection dictionary*; you may include or omit thresholds for any ofthe channel types present in your data. The values given here are reasonablefor this particular dataset, but may need to be adapted for differenthardware or recording conditions. For a more automated approach, considerusing the `autoreject package`_.<jupyter_code>reject_criteria = dict(mag=4000e-15,     # 4000 fT
                       grad=4000e-13,    # 4000 fT/cm
                       eeg=150e-6,       # 150 µV
                       eog=250e-6)       # 250 µV<jupyter_output><empty_output><jupyter_text>We'll also pass the event dictionary as the ``event_id`` parameter (so we canwork with easy-to-pool event labels instead of the integer event IDs), andspecify ``tmin`` and ``tmax`` (the time relative to each event at which tostart and end each epoch). As mentioned above, by default:class:`~mne.io.Raw` and :class:`~mne.Epochs` data aren't loaded into memory(they're accessed from disk only when needed), but here we'll force loadinginto memory using the ``preload=True`` parameter so that we can see theresults of the rejection criteria being applied:<jupyter_code>epochs = mne.Epochs(raw, events, event_id=event_dict, tmin=-0.2, tmax=0.5,
                    reject=reject_criteria, preload=True)<jupyter_output><empty_output><jupyter_text>Since we are only interested in the EEG channels, we will select them. __Note__ the command works "in-place", i.e. it changes the object in memory.<jupyter_code>epochs.pick_types(meg=False, eeg=True)<jupyter_output><empty_output><jupyter_text>Next we'll pool across left/right stimulus presentations so we can compareauditory versus visual responses. To avoid biasing our signals to theleft or right, we'll use :meth:`~mne.Epochs.equalize_event_counts` first torandomly sample epochs from each condition to match the number of epochspresent in the condition with the fewest good epochs.<jupyter_code>conds_we_care_about = ['auditory/left', 'auditory/right',
                       'visual/left', 'visual/right']
epochs.equalize_event_counts(conds_we_care_about)  # this operates in-place
aud_epochs = epochs['auditory']
vis_epochs = epochs['visual']<jupyter_output><empty_output><jupyter_text>Like :class:`~mne.io.Raw` objects, :class:`~mne.Epochs` objects also have anumber of built-in plotting methods. One is :meth:`~mne.Epochs.plot_image`,which shows each epoch as one row of an image map, with color representingsignal magnitude; the average evoked response and the sensor location areshown below the image:  Evoked plots<jupyter_code>aud_left_evoked = aud_epochs['left'].average()
aud_right_evoked = aud_epochs['right'].average()
aud_left_evoked.plot_joint();
aud_right_evoked.plot_joint();<jupyter_output><empty_output><jupyter_text>Now we make a "difference" object be combining the _left_ and _right_ stimuli<jupyter_code>diff_evoked = mne.combine_evoked([aud_left_evoked, -aud_right_evoked], weights='equal')
diff_evoked.plot_joint()<jupyter_output><empty_output><jupyter_text>Is the peak changed by a lower low pass setting, e.g. 30 Hz?To test this we first need to find the peak for both and then compare it.We can use the "get_peak" method for that. Remember that help and default values can be found with "?", e.g.: ```aud_left_evoked.get_peak?```<jupyter_code>aud_left_evoked.get_peak(return_amplitude=True)
aud_right_evoked.get_peak(return_amplitude=True)<jupyter_output><empty_output><jupyter_text>Now to compare to the different low pass we can filter the epochs. __NOTE__ This can potentially lead to massive artifacts and should __NOT__ be done. But since it is just a quick excerise we will still do it.<jupyter_code>epochs_30 = epochs.copy()  # make a copy to not overwrite the origianls 
epochs_30.filter(0, 30)<jupyter_output><empty_output><jupyter_text>And just to show it, we can chain a lot of the commands together. It is better to be explicit with your own code.<jupyter_code>epochs_30['auditory/left'].average().get_peak(return_amplitude=True)
epochs_30['auditory/right'].average().get_peak(return_amplitude=True)<jupyter_output><empty_output><jupyter_text>The are excatly identical to the orignals, so it does not appear that the low pass changes the peak time. However the amplitute is not the same, __why__?  Global field power<jupyter_code>aud_left_data = aud_epochs['left'].get_data()
aud_right_data = aud_epochs['right'].get_data()<jupyter_output><empty_output><jupyter_text>We have defined GFP as the standard deviation of the evoked signal over time, so if we take the ```np.std``` of the evoked data we get the GFP.<jupyter_code>left_gfp = aud_left_evoked._data.std(axis=0)
right_gfp = aud_right_evoked._data.std(axis=0)<jupyter_output><empty_output><jupyter_text>We will need to the times of the samples to make a useful plot, we can get those from any of the object with ```.times```, e.g. ```aud_left_evoked.times```<jupyter_code>times = aud_left_evoked.times<jupyter_output><empty_output><jupyter_text>Now to plot it<jupyter_code>plt.plot(times, left_gfp, 'r', label='Left');  # the 'r' makes the line red
plt.plot(times, right_gfp, 'b', label='right');  # label is what goes in the legend
plt.grid();
plt.legend();<jupyter_output><empty_output><jupyter_text>Permutation test<jupyter_code>from mne.stats import spatio_temporal_cluster_test
from mne.channels import find_ch_adjacency
from mpl_toolkits.axes_grid1 import make_axes_locatable<jupyter_output><empty_output><jupyter_text>The permutation test functions in MNE-python does not work on epochs or evoked directly so we have to extract the data we want to use first. We will call the data tensor for __X__<jupyter_code>X = [aud_epochs[condition].get_data() for condition in ['left', 'right']]  # as 3D matrix
X = [np.transpose(x, (0, 2, 1)) for x in X]  # transpose for clustering<jupyter_output><empty_output><jupyter_text>In order to let the algorithm know which sensors are neighbors we first generate an adjacency matrix.<jupyter_code>adjacency, ch_names = find_ch_adjacency(epochs.info, ch_type='eeg')

print(type(adjacency))  # it's a sparse matrix!

plt.imshow(adjacency.toarray(), cmap='binary', origin='lower',
           interpolation='nearest')
plt.xlabel('{} EEG channels'.format(len(ch_names)))
plt.ylabel('{} EEG channels '.format(len(ch_names)))
plt.grid();
plt.title('Between-sensor adjacency')
# set cluster threshold
# set family-wise p-value
p_accept = 0.01

cluster_stats = spatio_temporal_cluster_test(X,
                                             n_permutations=5000,
                                             threshold=None,
                                             tail=0,
                                             n_jobs=1,
                                             buffer_size=None,
                                             adjacency=adjacency)
T_obs, clusters, p_values, h_zero = cluster_stats
good_cluster_inds = np.where(p_values < p_accept)[0]<jupyter_output><empty_output><jupyter_text>Here we are going to visualize the results. There is a lot going on in the code so do not worry about it.<jupyter_code># configure variables for visualization
colors = {"left": "crimson", "right": 'steelblue'}

# organize data for plotting
evokeds = {cond: epochs[cond].average() for cond in ['left', 'right']}

# loop over clusters
for i_clu, clu_idx in enumerate(good_cluster_inds):
    # unpack cluster information, get unique indices
    time_inds, space_inds = np.squeeze(clusters[clu_idx])
    ch_inds = np.unique(space_inds)
    time_inds = np.unique(time_inds)

    # get topography for F stat
    f_map = T_obs[time_inds, ...].mean(axis=0)

    # get signals at the sensors contributing to the cluster
    sig_times = epochs.times[time_inds]

    # create spatial mask
    mask = np.zeros((f_map.shape[0], 1), dtype=bool)
    mask[ch_inds, :] = True

    # initialize figure
    fig, ax_topo = plt.subplots(1, 1, figsize=(10, 3))

    # plot average test statistic and mark significant sensors
    f_evoked = mne.EvokedArray(f_map[:, np.newaxis], epochs.info, tmin=0)
    f_evoked.plot_topomap(times=0, mask=mask, axes=ax_topo, cmap='Reds',
                          vmin=np.min, vmax=np.max, show=False,
                          colorbar=False, mask_params=dict(markersize=10))
    image = ax_topo.images[0]

    # create additional axes (for ERF and colorbar)
    divider = make_axes_locatable(ax_topo)

    # add axes for colorbar
    ax_colorbar = divider.append_axes('right', size='5%', pad=0.05)
    plt.colorbar(image, cax=ax_colorbar)
    ax_topo.set_xlabel(
        'Averaged F-map ({:0.3f} - {:0.3f} s)'.format(*sig_times[[0, -1]]))

    # add new axis for time courses and plot time courses
    ax_signals = divider.append_axes('right', size='300%', pad=1.2)
    title = 'Cluster #{0}, {1} sensor'.format(i_clu + 1, len(ch_inds))
    if len(ch_inds) > 1:
        title += "s (mean)"
    mne.viz.plot_compare_evokeds(evokeds, title=title, picks=ch_inds, axes=ax_signals,
                         colors=colors, show=False,
                         split_legend=True, truncate_yaxis='auto')

    # plot temporal cluster extent
    ymin, ymax = ax_signals.get_ylim()
    ax_signals.fill_betweenx((ymin, ymax), sig_times[0], sig_times[-1],
                             color='orange', alpha=0.3)

    # clean up viz
    mne.viz.tight_layout(fig=fig)
    fig.subplots_adjust(bottom=.05)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10062340153452685, 'fraction_numerical': 0.0167039641943734, 'mean_word_length': 4.447540269917283, 'pattern_counts': {'":': 2, '<': 86, '<?xml version=': 0, '>': 86, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '349446', 'license': 'MIT', 'n_tokens_mistral': 4245, 'n_tokens_neox': 3954, 'n_words': 1209, 'path': 'practical_classes/Week_2/week_2_solutions.ipynb', 'repo_name': 'victor-m-p/adv_cognitive_neuroscience', 'size': '12512'}