Extremely slow access to epochs after channel remapping.

[KSuljic]

I have two epoch lists. In one of the lists all channels have to be mirrored to the other side of the hemisphere (flipped). After that I want to concatenate both lists to one epoch file. So I'm using mne.rename_channels() to remap the channels of the one list. After that I'm using mne.equalize_channels() to get a uniform channel mapping of the two lists. After that I'm using mne.concatenate_epochs() to concat the lists. It works perfectly fine but if I try to access a condition after the concatenation in the resulting epoch instance, things become really slow. Before the access is really quick like instant. After it takes up to 20s. How could I resolve that I tried a lot already. Thank you very much!

[agramfort]

it's really hard to debug/profile without a script to replicate.

[KSuljic]

Got it, I will post it tomorrow. Thanks for the reply.

[KSuljic]

Unfortunately, I can't replicate the issue with sample data. The reason is probably that in my data I have around 400 Event IDs whereas here one has only 4. I'm flipping 200 out of these 400 IDs. Here I'm only flipping 2 out of 4. After concatenation of the epoch list with the 400 IDs of which 200 have flipped channels, the issue arises when I try to create evokeds out of them.

The reason I'm flipping the channels is that it is a tactile paradigm and I want to pool right stimuli with left stimuli by flipping them.

The procedure here is the same as in my code:

MNE version: 0.20.7

####
# Remap-Dict
OriginalChannels = ['Fpz',
 'Fz',
 'FCz',
 'Cz',
 'CPz',
 'Pz',
 'POz',
 'Oz',
 'Fp1',
 'F1',
 'FC1',
 'C1',
 'CP1',
 'P1',
 'F3',
 'FC3',
 'C3',
 'CP3',
 'P3',
 'PO3',
 'O1',
 'F5',
 'FC5',
 'C5',
 'CP5',
 'P5',
 'PO7',
 'F7',
 'FT7',
 'T7',
 'TP7',
 'P7',
 'Fp2',
 'F2',
 'FC2',
 'C2',
 'CP2',
 'P2',
 'F4',
 'FC4',
 'C4',
 'CP4',
 'P4',
 'PO4',
 'O2',
 'F6',
 'FC6',
 'C6',
 'CP6',
 'P6',
 'PO8',
 'F8',
 'FT8',
 'T8',
 'TP8',
 'P8',
 'FT10',
 'AF8',
 'AF7',
 'FT9']

RemappedChannels = ['Fpz',
 'Fz',
 'FCz',
 'Cz',
 'CPz',
 'Pz',
 'POz',
 'Oz',
 'Fp2',
 'F2',
 'FC2',
 'C2',
 'CP2',
 'P2',
 'F4',
 'FC4',
 'C4',
 'CP4',
 'P4',
 'PO4',
 'O2',
 'F6',
 'FC6',
 'C6',
 'CP6',
 'P6',
 'PO8',
 'F8',
 'FT8',
 'T8',
 'TP8',
 'P8',
 'Fp1',
 'F1',
 'FC1',
 'C1',
 'CP1',
 'P1',
 'F3',
 'FC3',
 'C3',
 'CP3',
 'P3',
 'PO3',
 'O1',
 'F5',
 'FC5',
 'C5',
 'CP5',
 'P5',
 'PO7',
 'F7',
 'FT7',
 'T7',
 'TP7',
 'P7',
 'FT9',
 'AF7',
 'AF8',
 'FT10']

remappingRtoL = dict(zip(OriginalChannels, RemappedChannels))

####
# Prepare Example Data
from mne.datasets import sample
data_path = sample.data_path()
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
raw = mne.io.read_raw_fif(fname, preload=True)

info = mne.create_info(ch_names=OriginalChannels, sfreq=500, ch_types='eeg')

events = mne.find_events(raw, stim_channel='STI 014')
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4, 'smiley': 5, 'buttonpress': 32}

raw = raw.pick('eeg')
raw.info = info

epochs = mne.Epochs(raw, events, event_id=event_dict, preload=True)



####
# Mirror Channels
Epochs = []

epochs_L = epochs['left']
Epochs.append(epochs_L)

epochs_R = epochs['right']
epochs_R.rename_channels(remappingRtoL) # Mirror Channels 
Epochs.append(epochs_R)

Epochs_eq = mne.channels.equalize_channels(Epochs, verbose=True)
        
#epochs with right channels flipped to left side
epochs_remapped = mne.concatenate_epochs(Epochs_eq * 150) # repeating epochs to increase effect



#unflipped epochs
epochs_OG = mne.concatenate_epochs([epochs] * 150) # repeating epochs to increase effect 


####
# Evoking
import datetime

#unflipped
a = datetime.datetime.now()
evoked_a = epochs_OG['auditory'].average()
evoked_v = epochs_OG['visual'].average()
b = datetime.datetime.now()

print(f'Time to evoke unflipped epochs: {b-a}')


#flipped
a = datetime.datetime.now()
evoked_remap_a = epochs_remapped['auditory'].average()
evoked_remap_v = epochs_remapped['visual'].average()
b = datetime.datetime.now()

print(f'Time to evoke flipped epochs: {b-a}')

[agramfort]

@KSuljic you can make fake events using sample data to replicate the issue

[KSuljic]

Thank you for the hint. I could replicate the issue to some extend but in my code the difference is stronger:

Time to evoke unflipped epochs: 0:00:00.056981
Time to evoke flipped epochs: 0:00:12.620404

It seems not much but in my analysis it makes out of 0.8s -> 3.5h.

Thank you for your aid.

####
#Prepare Example Data
from mne.datasets import sample
data_path = sample.data_path()
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
raw = mne.io.read_raw_fif(fname, preload=True)

info = mne.create_info(ch_names=OriginalChannels, sfreq=500, ch_types='eeg')
events = mne.find_events(raw, stim_channel='STI 014')

#Example Factors
moving = ['moving_no', 'moving_yes']
attention = ['attention_no', 'attention_yes']
phase = ['active', 'passive', 'stationary']
stim = ['single', 'empty']
side = ['side_left', 'side_right']
closeness = ['stim_close', 'stim_far']

factors = [side, phase, stim, closeness, moving, attention]

conditions = list(itertools.product(*factors))
conditions_agg = ['/'.join(c) for c in conditions]
condition_no = np.arange(1, len(conditions_agg))
condition_ids_dict = dict(zip(conditions_agg, condition_no))

#Creating Fake Events
import random

Events = np.empty((0,3))

for i in np.arange(20):
    x = random.uniform(0.9, 1.1)
    events_sim = events.copy()
    events_sim[:,0] = events[:,0] * x
        
    Events = np.concatenate((Events, events_sim))

    
import math
no_rep = math.ceil(len(Events)/len(condition_no))

event_ids_list = [condition_no] * no_rep
event_ids_list = np.concatenate((event_ids_list))

Events[:,2] = event_ids_list[:len(Events)]
Events = Events.astype('int')    


#Create epochs
raw = raw.pick('eeg')
raw.info = info

epochs = mne.Epochs(raw, Events, event_id=condition_ids_dict, event_repeated='drop', preload=True)



#####
Epochs_list = []

for k in list(epochs.event_id.keys()):
        epoch = epochs[k]
        Epochs_list.append(epoch)

Epochs_concat = mne.concatenate_epochs(Epochs_list)

#####
#Evoking
import datetime

#untouched epochs
a = datetime.datetime.now()
evoked_l = epochs['side_left'].average()
evoked_r = epochs['side_right'].average()
b = datetime.datetime.now()

print(f'Time to evoke untouched epochs: {b-a}')

#concatenated epochs
a = datetime.datetime.now()
evoked_remap_l = Epochs_concat['side_left'].average()
evoked_remap_r = Epochs_concat['side_right'].average()
b = datetime.datetime.now()

print(f'Time to evoke concatenated epochs: {b-a}')

[agramfort]

I confirm the weirdness:

In [6]: len(Epochs_remapped)
Out[6]: 6271

In [7]: len(epochs)
Out[7]: 6271

In [8]: %timeit epochs['side_left'].average()
754 ms ± 67.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [9]: %timeit Epochs_remapped['side_left'].average()
3.05 s ± 174 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [12]: len(epochs['side_left'])
Out[12]: 3193

In [13]: len(Epochs_remapped['side_left'])
Out[13]: 3193

@larsoner any hint?

[KSuljic]

The issue is also reproducible without mne.rename_channels(), only by appending the epochs to a list and concatinating with mne.concatenate_epochs().

[drammock]

The issue is also reproducible without mne.rename_channels(), only by appending the epochs to a list and concatinating with mne.concatenate_epochs().

In that case, can you edit your post above to simplify the code? The more minimal your example is, the easier it is for us to track down the problem.

[KSuljic]

Done. It seems the issue has something to do with appending the epochs and then using mne.concatenate_epochs(). Is there any way to mirror/flip the channels in place (meaning in the original epoch instance)? Then I wouldn't need concatenate_epochs(). I tried rename_channels() on the original epoch instance and reorder_channels(). After that I tried equalize_channels() in both cases. Both times it didn't work. The channels stay the same:

for k in list(epochs.event_id.keys()): 
        if 'side_right' in k:
            epochs[k].reorder_channels(RemappedChannels)

or

for k in list(epochs.event_id.keys()): 
        if 'side_right' in k:
            epochs[k].rename_channels(remappingRtoL)  

EDIT: I just realized that both functions don't operate in place.

[larsoner]

Here is an example that can be copy-pasted:

from copy import deepcopy
import datetime
import mne
data_path = mne.datasets.sample.data_path()
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
raw = mne.io.read_raw_fif(fname).pick_types(meg=False, eeg=True).load_data()
events = mne.make_fixed_length_events(raw, duration=0.5)
epochs = mne.Epochs(raw, events, preload=True)
epochs_list = [epochs[ii] for ii in range(len(epochs))]
epochs_concat = mne.concatenate_epochs(epochs_list)
assert len(epochs_concat) == len(epochs)

#untouched epochs
a = datetime.datetime.now()
evoked_l = epochs['1']
b = datetime.datetime.now()
print(f'Time untouched epochs:    {b-a}')

#concatenated epochs
a = datetime.datetime.now()
evoked_remap_l = epochs_concat['1']
b = datetime.datetime.now()
print(f'Time concatenated epochs: {b-a}')

a = datetime.datetime.now()
deepcopy(epochs_concat.drop_log)
b = datetime.datetime.now()
print(f'Time deepcopy drop_log:   {b-a}')

The output is:

Time untouched epochs:    0:00:00.046417
Time concatenated epochs: 0:00:00.637247
Time deepcopy drop_log:   0:00:00.632559

i.e., the problem is deepcopy on a list of list of str that has lots of elements.

One fix would be to make it a tuple of tuple of str, which is immutable. I actually like this change because epochs.drop_log is not meant to be modified directly by users anyway.

[larsoner]

(and it will speed up all epochs.copy operations, not just the one associated with epochs[...] -> epochs.__getitem__ -> epochs._getitem(..., copy=True) which is the cause of the slowdown here)

[agramfort]

+1 for using tuples then

…

[KSuljic]

Thank you very much for looking into it.
Is there any direct fix I could apply now?

[KSuljic]

@agramfort @larsoner Can I fix this by myself anyhow (with tuples?) or do I have to wait for an update? Unfortunately, due to this issue my analysis came to a stop. Thank you!

[larsoner]

You can try doing epochs.drop_log = [[]] * len(epochs) before subselection to see if it helps. If it doesn't, then you can do epochs.drop_log = []. Neither are actually correct, but they might allow you to continue your analysis until we fix the slowness.

[KSuljic]

Unfortunately, both didn't work. If I try to subselect it states:
IndexError: list assignment index out of range

[KSuljic]

I'll try to install the dev version in another conda env. It seems you fixed it there already.

[KSuljic]

In the dev version it works perfectly (roughly only double the time as you mentioned). Thank you very much!