<jupyter_start><jupyter_text>Quantifying eye movements from EOG recordings Eye movements are a major source of artifacts during M/EEG experiments. This is especially relevant in the context of naturalistic viewing paradigms because unconstrained eye movements might differ between experimental conditions, especially when we compare viewing of different video content. In my case, I want to compare brain networks dynamics during eyes-open resting state, viewing of abstract movie and viewing of naturalistic movie (i.e. involving complex scenes, characters and social interactions).Because we wanted to limit eye-movements during the experiment without degrading too much the viewing conditions, we asked ourselves if we should present the video at full size or reduced. To decide between these options, we presented the videos twice to our pilot subject, once at full size and once at a size reduced by 50%. That way, I can use the EOG data to see if presenting the video at 50% size led to reductions in eye-movements quantity and amplitude. If not, we will keep presenting at full size because it's more comfortable for the subject. If yes, we might decide to present at 50% or run more analysis to see if the eye-movements are problematic for the testing of our research hypothesis.**_As a quick note, I realize that I don't have a good criterion that will help me decide if I should go for 50% or 100% videos. That's what makes these analysis "exploratory"._ **Whatever. **Here are some important factors to consider, that are specific to our experimental setup :**- Distance to screen : 100cm- Screen width : 43cm- Screen height : 33cm**Here are the different kinds of eye-movements we might want to quantify :**- Blinks : the eyelids closes for a split of a seconde, generating a huge stereotyped potential mostly visible on vertical EOG- Saccades : both eyes roll quickly in one direction on another in order to fixate a new location. Most saccades are not solely vertical nor horizontal, so they should be seen on both vEOG and hEOG)- Smooth poursuit eye movements (SPEM): these are the continuous movements performed to track a moving object. They are usually more difficult to identify because they produce less stereotyped activity. They are very important in the context of naturalistic viewing but we might be forced to drop their analysis if I don't find a good approach to identify them from EOG (which looks difficult if not impossible).  Importing the dataFirst we start by importing the data using MNE-python<jupyter_code>import mne
from mne.io import read_raw_ctf

FILEPATH = '/home/hyruuk/pCloudDrive/science/neuromod/inscapes_MEG/raw/20190529/'

filename_100 = 'P03_NEUROMOD-jboyle_20190529_03.ds'
filename_50 = 'P03_NEUROMOD-jboyle_20190529_05.ds'
fname_100 = FILEPATH + filename_100
fname_50 = FILEPATH + filename_50

# Open the files and load the data
raw_100 = read_raw_ctf(fname_100, preload=True).resample(200, npad='auto')
raw_50 = read_raw_ctf(fname_50, preload=True).resample(200, npad='auto')

# Here we pick only the relevant channels : EEG057 is vEOG, and EEG058 is hEOG
picks = mne.pick_channels(raw_100.info['ch_names'], ['EEG057', 'EEG058']) # For these analysis we only need the EOG chans
EOG_100 = raw_100.get_data(picks=picks)
EOG_50 = raw_50.get_data(picks=picks)
print(raw_100.info['sfreq'])
# Close the raw objects in order to free up some memory
raw_100.close()
raw_50.close()<jupyter_output>ds directory : /home/hyruuk/pCloudDrive/science/neuromod/inscapes_MEG/raw/20190529/P03_NEUROMOD-jboyle_20190529_03.ds
    res4 data read.
    hc data read.
    Separate EEG position data file read.
    Quaternion matching (desired vs. transformed):
      -2.12   72.19    0.00 mm <->   -2.12   72.19    0.00 mm (orig :  -52.71   50.46 -260.44 mm) diff =    0.000 mm
       2.12  -72.19    0.00 mm <->    2.12  -72.19    0.00 mm (orig :   52.79  -48.18 -258.67 mm) diff =    0.000 mm
      90.98    0.00    0.00 mm <->   90.98   -0.00    0.00 mm (orig :   64.10   65.71 -261.52 mm) diff =    0.000 mm
    Coordinate transformations established.
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
    Measurement info composed.
Finding samples for /home/hyruuk/pCloudDrive/science/neuromod/inscapes_MEG/raw/20190529/P03_NEUROMOD-jboyle_20190529_03.ds/P03_NEUROMOD-jboyle_20190529_03.meg4: 
    System clock channel is available, checking which samples are val[...]<jupyter_text>VisualisationNow lets compare visually the two signals, just to get an idea of what they are about.<jupyter_code>import matplotlib
import matplotlib.pyplot as plt
plt.figure(figsize=(20,10))
plt.plot(EOG_100[0], linewidth=0.5)
plt.plot(EOG_50[0], linewidth=0.5)
plt.title('Vertical EOG_057')
plt.legend(('100%', '50%'))
#plt.savefig('Inscapes_vEOG_057')
plt.figure(figsize=(20,10))
plt.plot(EOG_100[1], linewidth=0.5)
plt.plot(EOG_50[1], linewidth=0.5)
plt.title('Horizontal EOG_058')
plt.legend(('100%', '50%'))
#plt.savefig('Inscapes_hEOG_058')<jupyter_output><empty_output><jupyter_text>Let's look at the vertical EOG closer now<jupyter_code>plt.figure(figsize=(20,10))
plt.plot(EOG_100[0,40000:60000], linewidth=0.5)
plt.plot(EOG_50[0, 40000:60000], linewidth=0.5)
plt.title('Vertical EOG_057')
plt.legend(('100%', '50%'))
#plt.savefig('Inscapes_vEOG_057')<jupyter_output><empty_output><jupyter_text>Uh oh, there might be more eyeblinks in the 50% than in the 100% condition...? We can see that :- Blinks are clearly visible on the vEOG channel, but not necessarly on the hEOG channel.- Overall the blue signal (100%) look more noisy and seem to display more events than the orange signal (50%).- Finally, huge spikes of activity tend to appear in all signals. This is bad news because I cannot then use the variability of each signal as a marker of the quantity of blinks/saccades, because the spikes might greatly affect this measure.   Using univariate convolutional sparse coding (CSC) to detect "atoms" (stereotyped patterns of activity) and their occurence through the signal.Since I had no idea on how to approach the problem, I went to see my friend [<NAME>](https://jasmainak.github.io/), who advised me **against** using ICA to compare the weight of eye-related components in differents conditions (and so, from different decompositions). The reason being that ICA is a non-convex model, which seems to mean that there is no unique solution for a decomposition, and starting the process with different random seeds might lead to different results (according to his words).He suggested however that I use a new tool he contributed developping, to see if it could be of any help for my questions. This tool is called [alphaCSC](https://alphacsc.github.io/) and it looks very promising for many neuroscientific applications. I have to admit that I don't understand much of it's working for now, but I'm on it.Here are the few things I got from now :- It's based on a machine learning technique called sparse dictionnary learning- It works by learning "atoms", or prototypical events that are reproduced across the time series- When the atoms are worked out, the technique allows to decompose the signal into a sum of "activation functions" (?) that describe the temporal occurence of the different atoms.So, if I'm correct, we could try to decompose our signal, identify the atoms relating to blinks and saccades, and then compare the activation functions of both recordings. If the activation function of recording A shows more events than that of recording B for the "blink atom", I could then say that there were more blinks during the recording A.**One thing that should be noted here is that Mainak advised me to concatenate the signals from the two conditions in order to learn atoms that would be shared between the recordings. If I don't do it that way, the activation functions I will want to compare would have different meaning, and thus couldn't be compared.**  So, here is my first try :<jupyter_code>import numpy as np
from numpy import hamming
from alphacsc import learn_d_z
from mne.filter import filter_data

reg = 6.0
n_times = 2500
n_times_atom = 350
n_trials = 100
n_atoms = 3
n_iter = 60
random_state = 10
sfreq = 1200

# Let's just take the code from alphaCSC's tutorials and use our data instead.
X = np.concatenate((EOG_100[0], EOG_50[0]))
#X = filter_data(X.astype(np.float64), sfreq, l_freq=1, h_freq=None, fir_design='firwin')
#print(X)
#0/0

# First we need to segment the signal. I'm not sure why it's necessary
overlap = 0
starts = np.arange(0, X.shape[0] - n_times, n_times - overlap)
stops = np.arange(n_times, X.shape[0], n_times - overlap)
X_new = []
for idx, (start, stop) in enumerate(zip(starts, stops)):
    if idx >= n_trials:
        break
    X_new.append(X[start:stop])
X_new = np.vstack(X_new)
del X

# Then we normalize the signal
X_new -= np.mean(X_new)
X_new /= np.std(X_new)

# CSC is subject to edge effects, so this hamming window should help us prevent it
X_new *= hamming(n_times)[None, :]

pobj, times, d_hat, z_hat, reg = learn_d_z(X_new, n_atoms, n_times_atom,
                                           reg=reg, n_iter=n_iter,
                                           random_state=random_state, n_jobs=1)<jupyter_output><empty_output><jupyter_text>Let's plot it<jupyter_code>plt.figure()
plt.plot(d_hat.T)
plt.savefig('atoms_OE_100+50_hEOG')
plt.show()


import mne
from mne.io import read_raw_ctf

FILEPATH = '/home/hyruuk/pCloudDrive/science/neuromod/inscapes_MEG/raw/20190529/'

filename_100 = 'P03_NEUROMOD-jboyle_20190529_03.ds'
filename_50 = 'P03_NEUROMOD-jboyle_20190529_05.ds'
fname_100 = FILEPATH + filename_100
fname_50 = FILEPATH + filename_50

# Open the files and load the data
raw_100 = read_raw_ctf(fname_100, preload=True)
raw_50 = read_raw_ctf(fname_50, preload=True)

# Here we pick only the relevant channels : EEG057 is vEOG, and EEG058 is hEOG
picks = mne.pick_channels(raw_100.info['ch_names'], ['EEG057', 'EEG058']) # For these analysis we only need the EOG chans
EOG_100 = raw_100.get_data(picks=picks)
EOG_50 = raw_50.get_data(picks=picks)
print(raw_100.info['sfreq'])
# Close the raw objects in order to free up some memory
raw_100.close()
raw_50.close()

### This saves the EOG data
from scipy.io import savemat, loadmat

savemat('EOG_100.mat', {'EOG_100':EOG_100, 'orientation':['vertical', 'horizontal']})
savemat('EOG_50.mat', {'EOG_50':EOG_50, 'orientation':['vertical', 'horizontal']})
EOG_concat = np.concatenate((EOG_100[0], EOG_50[0]))
a = loadmat('EOG_100.mat')
import numpy as np
from scipy.io import loadmat

import matplotlib.pyplot as plt

from mne.filter import filter_data
from alphacsc import GreedyCDL

n_atoms = 2
reg = 0.5
n_iter = 60
n_times_atom = 400
random_state = 42

#X = loadmat('EOG_50.mat')['EOG_50']

#X = X[0, 50000:100000][None, :]  # crop signal to demo

X = np.concatenate((EOG_100[0], EOG_50[0])).reshape(1, EOG_100[0].shape[0] + EOG_50[0].shape[0])
print(X.shape)

# filter data
sfreq = 200.
X = filter_data(X, sfreq, l_freq=1.0, h_freq=40.)

# z-score to have sane reg parameter
X -= np.mean(X)
X /= np.std(X)

cdl = GreedyCDL(
    # Shape of the dictionary
    n_atoms=n_atoms,
    n_times_atom=n_times_atom,
    # Request a rank1 dictionary with unit norm temporal and spatial maps
    rank1=False,
    uv_constraint='joint',
    # apply a temporal window reparametrization
    window=True,
    # at the end, refit the activations with fixed support and no reg to unbias
    unbiased_z_hat=True,
    # Initialize the dictionary with random chunk from the data
    D_init='chunk',
    # rescale the regularization parameter to be a percentage of lambda_max
    lmbd_max="scaled",
    reg=reg,
    # Number of iteration for the alternate minimization and cvg threshold
    n_iter=100,
    eps=1e-6,
    # solver for the z-step
    solver_z="lgcd",
    solver_z_kwargs={'tol': 1e-3,
                     'max_iter': 100000},
    # solver for the d-step
    solver_d='alternate_adaptive',
    solver_d_kwargs={'max_iter': 300},
    # sort atoms by explained variances
    sort_atoms=True,
    # Technical parameters
    verbose=1,
    random_state=0,
    n_jobs=1)

cdl.fit(X[None, ...])
plt.plot(cdl.D_hat_[0][0])
plt.figure()
plt.plot(cdl.z_hat_[0][0], '-o')
plt.plot(X.T, alpha=0.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Detecting EOG events with MNE-PythonThe MNE-Python toolbox conveniently provides a function that can identify EOG events, including blinks and maybe some of the biggest saccades. From what I understood, it works by filtering the EOG channels and then detecting peaks in the signal.Let's see if that function detects more events for the 100% than for the 50% condition :<jupyter_code>from mne.preprocessing import find_eog_events

raw_100 = read_raw_ctf(fname_100, preload=True, verbose=False)
raw_50 = read_raw_ctf(fname_50, preload=True, verbose=False)
eog_events_100 = find_eog_events(raw_100, ch_name='EEG057', verbose=False)
eog_events_50 = find_eog_events(raw_50, ch_name='EEG057', verbose=False)
print('Number of events in 100% condition : ' + str(len(eog_events_100)))
print('Number of events in 50% condition : ' + str(len(eog_events_50)))<jupyter_output>Number of events in 100% condition : 172
Number of events in 50% condition : 80<jupyter_text>Ok so, looks like there are clearly more events in the 100% condition... But what is the impact of these events on the MEG data ?To answer that question, we can segment the data around the events, and average them to produce ERFs. If the ERFs are strong, we can assume that the signal will be very contamined by occurence of a blink.<jupyter_code>from hytools.meg_utils import get_ch_pos, array_topoplot

baseline=(None, 0.0)
tmin=-0.2
tmax=0.8

event_id = 998
epochs_100 = mne.Epochs(raw_100, events=eog_events_100[:80], event_id=event_id, tmin=tmin,
                    tmax=tmax, baseline=baseline, reject=None, picks=None, preload=True)
epochs_50 = mne.Epochs(raw_50, events=eog_events_50, event_id=event_id, tmin=tmin,
                    tmax=tmax, baseline=baseline, reject=None, picks=None, preload=True)

evoked_100 = epochs_100.average()
evoked_50 = epochs_50.average()

evoked_100.plot();
evoked_50.plot();<jupyter_output>80 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 80 events and 1201 original time points ...
0 bad epochs dropped
80 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 80 events and 1201 original time points ...
0 bad epochs dropped
{'fraction_non_alphanumeric': 0.08692081308603107, 'fraction_numerical': 0.05057833139415509, 'mean_word_length': 4.587762906309751, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '463485', 'license': 'Unlicense', 'n_tokens_mistral': 5242, 'n_tokens_neox': 4647, 'n_words': 1741, 'path': 'eog_analysis_alphaCSC.ipynb', 'repo_name': 'hyruuk/eog_analysis', 'size': '14615'}