# ICA: EOG correction

%matplotlib qt

import mne
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from mne.preprocessing import ICA # for ICA


data_path = '/Users/kevinhsu/Documents/D/000_experiment/compounding_VR/2013_megdata'

sid = 2
# import raw CON

confile = data_path + '/de%.3dnn.con' % sid  #
#   /Users/kevinhsu/Documents/D/000_experiment/compounding_VR/2013_megdata/de001nn.con

raw = mne.io.kit.read_raw_kit(confile, stim = [192, 193, 194], slope = '+')
events = mne.find_events(raw, stim_channel='STI 014')

# preload raw MEG data
raw.load_data()

# duration of epochs
tmin = -0.1         # pre stimulis interval (in seconds) #
tmax = 0.7          # post stimulus interval #

# artifact rejection criteria
reject = dict(mag=1.5e-12)

# using strings (words) to label conditions
event_id = {'words':1, 'pseudo': 4}

# keep MEG channels
picks = mne.pick_types(raw.info, meg= True, stim = False)    # channels to use in epochs #

# baseline correction
baseline = (None, 0) # what to use as baseline comparison - here it's pre-stim interval #

# do the actual epoching processes
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj = True,
                    picks = picks, baseline=baseline, preload = True,
                    reject=None)

Extracting SQD Parameters from /Users/kevinhsu/Documents/D/000_experiment/compounding_VR/2013_megdata/de002nn.con...
Creating Raw.info structure...
Setting channel info structure...
Creating Info structure...
Current compensation grade : 0
Ready.
250 events found
Event IDs: [1 4]
Reading 0 ... 1334999 = 0.000 ... 1334.999 secs...
250 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 250 events and 801 original time points ...
0 bad epochs dropped

epochs

<Epochs | 250 events (all good), -0.1 - 0.7 sec, baseline [None, 0], ~240.2 MB, data loaded,
'pseudo': 125
'words': 125>

# ICA

kit_lay = mne.channels.layout.read_layout(data_path + '/KIT-TW-157.lout')

method = 'extended-infomax'
n_components = 15
max_pca_components = 15
random_state = 23

decim = 3

picks = mne.pick_types(epochs.info, meg=True, eeg=False, stim=False, exclude='bads')

ica = ICA(n_components=n_components, max_pca_components=max_pca_components, 
          random_state=random_state, method=method)

ica.fit(epochs, decim=decim, picks=picks)

ica.plot_components(layout = kit_lay)

Fitting ICA to data using 157 channels (please be patient, this may take a while)

<ipython-input-4-49160ed4f368>:15: DeprecationWarning: method='extended-infomax' is deprecated and will be removed in 0.19. If you want to use Extended Infomax, specify method='infomax' together with fit_params=dict(extended=True).
random_state=random_state, method=method)

Selection by number: 15 components
Computing Extended Infomax ICA
Fitting ICA took 9.0s.

[<Figure size 1500x1100 with 15 Axes>]

ic_ = 0
ica.plot_properties(epochs, ic_, topomap_args={'layout': kit_lay})
ica.plot_sources(epochs)

Using multitaper spectrum estimation with 7 DPSS windows


ica.plot_overlay(epochs.average(), exclude=[0, 2])

Transforming to ICA space (15 components)
Zeroing out 2 ICA components


epochs = ica.apply(epochs, exclude= [0, 2])
epochs.filter(0, 30, phase= 'zero-double').drop_bad(reject=reject)
epochs

Transforming to ICA space (15 components)
Zeroing out 2 ICA components
Setting up low-pass filter at 30 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-12 dB cutoff frequency: 33.75 Hz)
- Filter length: 441 samples (0.441 sec)

Rejecting epoch based on MAG : ['MEG 123', 'MEG 134', 'MEG 135', 'MEG 136']
Rejecting epoch based on MAG : ['MEG 153']
...
Rejecting epoch based on MAG : ['MEG 135']
Rejecting epoch based on MAG : ['MEG 135', 'MEG 136']
Rejecting epoch based on MAG : ['MEG 155']
111 bad epochs dropped

<Epochs | 139 events (all good), -0.1 - 0.7 sec, baseline [None, 0], ~133.7 MB, data loaded,
'pseudo': 70
'words': 69>

evoked1 = epochs['words'].average().apply_baseline((-0.1, 0.0))
evoked2 = epochs['pseudo'].average().apply_baseline((-0.1, 0.0))

evoked1.plot()

Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)


tmp = [evoked1, evoked2]
colors = 'green', 'red'
mne.viz.plot_evoked_topo(tmp, color = colors, title = 'LDT / lexicality', layout = kit_lay)


# saving an epoch object into a fif file
epochs.save(data_path + '/%.3dnn_ica_epo.fif' % sid, overwrite= True)


# saving multiple evoked objects into a fif file
mne.write_evokeds('%.3d-ica_ave.fif' % sid, [evoked1, evoked2])