Coverage for ibllib/io/extractors/ephys_fpga.py: 91%
551 statements
« prev ^ index » next coverage.py v7.5.4, created at 2024-07-08 17:16 +0100
« prev ^ index » next coverage.py v7.5.4, created at 2024-07-08 17:16 +0100
1"""Data extraction from raw FPGA output.
3The behaviour extraction happens in the following stages:
5 1. The NI DAQ events are extracted into a map of event times and TTL polarities.
6 2. The Bpod trial events are extracted from the raw Bpod data, depending on the task protocol.
7 3. As protocols may be chained together within a given recording, the period of a given task
8 protocol is determined using the 'spacer' DAQ signal (see `get_protocol_period`).
9 4. Physical behaviour events such as stim on and reward time are separated out by TTL length or
10 sequence within the trial.
11 5. The Bpod clock is sync'd with the FPGA using one of the extracted trial events.
12 6. The Bpod software events are then converted to FPGA time.
14Examples
15--------
16For simple extraction, use the FPGATrials class:
18>>> extractor = FpgaTrials(session_path)
19>>> trials, _ = extractor.extract(update=False, save=False)
21Notes
22-----
23Sync extraction in this module only supports FPGA data acquired with an NI DAQ as part of a
24Neuropixels recording system, however a sync and channel map extracted from a different DAQ format
25can be passed to the FpgaTrials class.
27See Also
28--------
29For dynamic pipeline sessions it is best to call the extractor via the BehaviorTask class.
31TODO notes on subclassing various methods of FpgaTrials for custom hardware.
32"""
33import logging
34from itertools import cycle
35from pathlib import Path
36import uuid
37import re
38import warnings
39from functools import partial
41import matplotlib.pyplot as plt
42from matplotlib.colors import TABLEAU_COLORS
43import numpy as np
44from packaging import version
46import spikeglx
47import ibldsp.utils
48import one.alf.io as alfio
49from one.alf.files import filename_parts
50from iblutil.util import Bunch
51from iblutil.spacer import Spacer
53import ibllib.exceptions as err
54from ibllib.io import raw_data_loaders as raw, session_params
55from ibllib.io.extractors.bpod_trials import extract_all as bpod_extract_all
56import ibllib.io.extractors.base as extractors_base
57from ibllib.io.extractors.training_wheel import extract_wheel_moves
58from ibllib import plots
59from ibllib.io.extractors.default_channel_maps import DEFAULT_MAPS
61_logger = logging.getLogger(__name__)
63SYNC_BATCH_SIZE_SECS = 100
64"""int: Number of samples to read at once in bin file for sync."""
66WHEEL_RADIUS_CM = 1 # stay in radians
67"""float: The radius of the wheel used in the task. A value of 1 ensures units remain in radians."""
69WHEEL_TICKS = 1024
70"""int: The number of encoder pulses per channel for one complete rotation."""
72BPOD_FPGA_DRIFT_THRESHOLD_PPM = 150
73"""int: Throws an error if Bpod to FPGA clock drift is higher than this value."""
75CHMAPS = {'3A':
76 {'ap':
77 {'left_camera': 2,
78 'right_camera': 3,
79 'body_camera': 4,
80 'bpod': 7,
81 'frame2ttl': 12,
82 'rotary_encoder_0': 13,
83 'rotary_encoder_1': 14,
84 'audio': 15
85 }
86 },
87 '3B':
88 {'nidq':
89 {'left_camera': 0,
90 'right_camera': 1,
91 'body_camera': 2,
92 'imec_sync': 3,
93 'frame2ttl': 4,
94 'rotary_encoder_0': 5,
95 'rotary_encoder_1': 6,
96 'audio': 7,
97 'bpod': 16,
98 'laser': 17,
99 'laser_ttl': 18},
100 'ap':
101 {'imec_sync': 6}
102 },
103 }
104"""dict: The default channel indices corresponding to various devices for different recording systems."""
107def data_for_keys(keys, data):
108 """Check keys exist in 'data' dict and contain values other than None."""
109 return data is not None and all(k in data and data.get(k, None) is not None for k in keys) 1sdNUOPQDRSETz2Vcelmnkhb
112def get_ibl_sync_map(ef, version):
113 """
114 Gets default channel map for the version/binary file type combination
115 :param ef: ibllib.io.spikeglx.glob_ephys_file dictionary with field 'ap' or 'nidq'
116 :return: channel map dictionary
117 """
118 # Determine default channel map
119 if version == '3A': 13sJKFGdpqrNOPQDRSETzHI2Vacfelmnkjhbg
120 default_chmap = CHMAPS['3A']['ap'] 132aclmnb
121 elif version == '3B': 13sJKFGdpqrNOPQDRSETzHIVafekjhbg
122 if ef.get('nidq', None): 13sJKFGdpqrNOPQDRSETzHIVafekjhbg
123 default_chmap = CHMAPS['3B']['nidq'] 13sJKFGdpqrNOPQDRSETzHIVafekjhbg
124 elif ef.get('ap', None): 13pqrVkj
125 default_chmap = CHMAPS['3B']['ap'] 13pqrVkj
126 # Try to load channel map from file
127 chmap = spikeglx.get_sync_map(ef['path']) 13sJKFGdpqrNOPQDRSETzHI2Vacfelmnkjhbg
128 # If chmap provided but not with all keys, fill up with default values
129 if not chmap: 13sJKFGdpqrNOPQDRSETzHI2Vacfelmnkjhbg
130 return default_chmap 13JKFGpqrHI2Vafkjg
131 else:
132 if data_for_keys(default_chmap.keys(), chmap): 1sdNOPQDRSETz2Vcelmnkhb
133 return chmap 1sdz2Vcelmnkhb
134 else:
135 _logger.warning("Keys missing from provided channel map, " 1NOPQDRSET2lmb
136 "setting missing keys from default channel map")
137 return {**default_chmap, **chmap} 1NOPQDRSET2lmb
140def _sync_to_alf(raw_ephys_apfile, output_path=None, save=False, parts=''):
141 """
142 Extracts sync.times, sync.channels and sync.polarities from binary ephys dataset
144 :param raw_ephys_apfile: bin file containing ephys data or spike
145 :param output_path: output directory
146 :param save: bool write to disk only if True
147 :param parts: string or list of strings that will be appended to the filename before extension
148 :return:
149 """
150 # handles input argument: support ibllib.io.spikeglx.Reader, str and pathlib.Path
151 if isinstance(raw_ephys_apfile, spikeglx.Reader): 1wxypqrABCj
152 sr = raw_ephys_apfile 1wxypqrABCj
153 else:
154 raw_ephys_apfile = Path(raw_ephys_apfile)
155 sr = spikeglx.Reader(raw_ephys_apfile)
156 if not (opened := sr.is_open): 1wxypqrABCj
157 sr.open()
158 # if no output, need a temp folder to swap for big files
159 if not output_path: 1wxypqrABCj
160 output_path = raw_ephys_apfile.parent
161 file_ftcp = Path(output_path).joinpath(f'fronts_times_channel_polarity{uuid.uuid4()}.bin') 1wxypqrABCj
163 # loop over chunks of the raw ephys file
164 wg = ibldsp.utils.WindowGenerator(sr.ns, int(SYNC_BATCH_SIZE_SECS * sr.fs), overlap=1) 1wxypqrABCj
165 fid_ftcp = open(file_ftcp, 'wb') 1wxypqrABCj
166 for sl in wg.slice: 1wxypqrABCj
167 ss = sr.read_sync(sl) 1wxypqrABCj
168 ind, fronts = ibldsp.utils.fronts(ss, axis=0) 1wxypqrABCj
169 # a = sr.read_sync_analog(sl)
170 sav = np.c_[(ind[0, :] + sl.start) / sr.fs, ind[1, :], fronts.astype(np.double)] 1wxypqrABCj
171 sav.tofile(fid_ftcp) 1wxypqrABCj
172 # close temp file, read from it and delete
173 fid_ftcp.close() 1wxypqrABCj
174 tim_chan_pol = np.fromfile(str(file_ftcp)) 1wxypqrABCj
175 tim_chan_pol = tim_chan_pol.reshape((int(tim_chan_pol.size / 3), 3)) 1wxypqrABCj
176 file_ftcp.unlink() 1wxypqrABCj
177 sync = {'times': tim_chan_pol[:, 0], 1wxypqrABCj
178 'channels': tim_chan_pol[:, 1],
179 'polarities': tim_chan_pol[:, 2]}
180 # If opened Reader was passed into function, leave open
181 if not opened: 1wxypqrABCj
182 sr.close()
183 if save: 1wxypqrABCj
184 out_files = alfio.save_object_npy(output_path, sync, 'sync', 1wxypqrABCj
185 namespace='spikeglx', parts=parts)
186 return Bunch(sync), out_files 1wxypqrABCj
187 else:
188 return Bunch(sync)
191def _rotary_encoder_positions_from_fronts(ta, pa, tb, pb, ticks=WHEEL_TICKS, radius=WHEEL_RADIUS_CM, coding='x4'):
192 """
193 Extracts the rotary encoder absolute position as function of time from fronts detected
194 on the 2 channels. Outputs in units of radius parameters, by default radians
195 Coding options detailed here: http://www.ni.com/tutorial/7109/pt/
196 Here output is clockwise from subject perspective
198 :param ta: time of fronts on channel A
199 :param pa: polarity of fronts on channel A
200 :param tb: time of fronts on channel B
201 :param pb: polarity of fronts on channel B
202 :param ticks: number of ticks corresponding to a full revolution (1024 for IBL rotary encoder)
203 :param radius: radius of the wheel. Defaults to 1 for an output in radians
204 :param coding: x1, x2 or x4 coding (IBL default is x4)
205 :return: indices vector (ta) and position vector
206 """
207 if coding == 'x1': 1$Z%4Lsdacfelmnkjhbg
208 ia = np.searchsorted(tb, ta[pa == 1]) 1ZL
209 ia = ia[ia < ta.size] 1ZL
210 ia = ia[pa[ia] == 1] 1ZL
211 ib = np.searchsorted(ta, tb[pb == 1]) 1ZL
212 ib = ib[ib < tb.size] 1ZL
213 ib = ib[pb[ib] == 1] 1ZL
214 t = np.r_[ta[ia], tb[ib]] 1ZL
215 p = np.r_[ia * 0 + 1, ib * 0 - 1] 1ZL
216 ordre = np.argsort(t) 1ZL
217 t = t[ordre] 1ZL
218 p = p[ordre] 1ZL
219 p = np.cumsum(p) / ticks * np.pi * 2 * radius 1ZL
220 return t, p 1ZL
221 elif coding == 'x2': 1$%4Lsdacfelmnkjhbg
222 p = pb[np.searchsorted(tb, ta) - 1] * pa 1$%L
223 p = - np.cumsum(p) / ticks * np.pi * 2 * radius / 2 1$%L
224 return ta, p 1$%L
225 elif coding == 'x4': 14Lsdacfelmnkjhbg
226 p = np.r_[pb[np.searchsorted(tb, ta) - 1] * pa, -pa[np.searchsorted(ta, tb) - 1] * pb] 14Lsdacfelmnkjhbg
227 t = np.r_[ta, tb] 14Lsdacfelmnkjhbg
228 ordre = np.argsort(t) 14Lsdacfelmnkjhbg
229 t = t[ordre] 14Lsdacfelmnkjhbg
230 p = p[ordre] 14Lsdacfelmnkjhbg
231 p = - np.cumsum(p) / ticks * np.pi * 2 * radius / 4 14Lsdacfelmnkjhbg
232 return t, p 14Lsdacfelmnkjhbg
235def _assign_events_to_trial(t_trial_start, t_event, take='last', t_trial_end=None):
236 """
237 Assign events to a trial given trial start times and event times.
239 Trials without an event result in nan value in output time vector.
240 The output has a consistent size with t_trial_start and ready to output to alf.
242 Parameters
243 ----------
244 t_trial_start : numpy.array
245 An array of start times, used to bin edges for assigning values from `t_event`.
246 t_event : numpy.array
247 An array of event times to assign to trials.
248 take : str {'first', 'last'}, int
249 'first' takes first event > t_trial_start; 'last' takes last event < the next
250 t_trial_start; an int defines the index to take for events within trial bounds. The index
251 may be negative.
252 t_trial_end : numpy.array
253 Optional array of end times, used to bin edges for assigning values from `t_event`.
255 Returns
256 -------
257 numpy.array
258 An array the length of `t_trial_start` containing values from `t_event`. Unassigned values
259 are replaced with np.nan.
261 See Also
262 --------
263 FpgaTrials._assign_events - Assign trial events based on TTL length.
264 """
265 # make sure the events are sorted
266 try: 1Mdiacfehbg
267 assert np.all(np.diff(t_trial_start) >= 0) 1Mdiacfehbg
268 except AssertionError: 1M
269 raise ValueError('Trial starts vector not sorted') 1M
270 try: 1Mdiacfehbg
271 assert np.all(np.diff(t_event) >= 0) 1Mdiacfehbg
272 except AssertionError: 1M
273 raise ValueError('Events vector is not sorted') 1M
275 # remove events that happened before the first trial start
276 remove = t_event < t_trial_start[0] 1Mdiacfehbg
277 if t_trial_end is not None: 1Mdiacfehbg
278 if not np.all(np.diff(t_trial_end) >= 0): 1diacfehbg
279 raise ValueError('Trial end vector not sorted')
280 if not np.all(t_trial_end[:-1] < t_trial_start[1:]): 1diacfehbg
281 raise ValueError('Trial end times must not overlap with trial start times')
282 # remove events between end and next start, and after last end
283 remove |= t_event > t_trial_end[-1] 1diacfehbg
284 for e, s in zip(t_trial_end[:-1], t_trial_start[1:]): 1diacfehbg
285 remove |= np.logical_and(s > t_event, t_event >= e) 1diacfehbg
286 t_event = t_event[~remove] 1Mdiacfehbg
287 ind = np.searchsorted(t_trial_start, t_event) - 1 1Mdiacfehbg
288 t_event_nans = np.zeros_like(t_trial_start) * np.nan 1Mdiacfehbg
289 # select first or last element matching each trial start
290 if take == 'last': 1Mdiacfehbg
291 iall, iu = np.unique(np.flip(ind), return_index=True) 1Mdiacfehbg
292 t_event_nans[iall] = t_event[- (iu - ind.size + 1)] 1Mdiacfehbg
293 elif take == 'first': 1Mdiacfehbg
294 iall, iu = np.unique(ind, return_index=True) 1Mdiacfehbg
295 t_event_nans[iall] = t_event[iu] 1Mdiacfehbg
296 else: # if the index is arbitrary, needs to be numeric (could be negative if from the end)
297 iall = np.unique(ind) 1Ma
298 minsize = take + 1 if take >= 0 else - take 1Ma
299 # for each trial, take the take nth element if there are enough values in trial
300 for iu in iall: 1Ma
301 match = t_event[iu == ind] 1Ma
302 if len(match) >= minsize: 1Ma
303 t_event_nans[iu] = match[take] 1Ma
304 return t_event_nans 1Mdiacfehbg
307def get_sync_fronts(sync, channel_nb, tmin=None, tmax=None):
308 """
309 Return the sync front polarities and times for a given channel.
311 Parameters
312 ----------
313 sync : dict
314 'polarities' of fronts detected on sync trace for all 16 channels and their 'times'.
315 channel_nb : int
316 The integer corresponding to the desired sync channel.
317 tmin : float
318 The minimum time from which to extract the sync pulses.
319 tmax : float
320 The maximum time up to which we extract the sync pulses.
322 Returns
323 -------
324 Bunch
325 Channel times and polarities.
326 """
327 selection = sync['channels'] == channel_nb 1vutsJKFGdpqriNUWXYOPQDRSETzHIVacfelmnkjhbg
328 selection = np.logical_and(selection, sync['times'] <= tmax) if tmax else selection 1vutsJKFGdpqriNUWXYOPQDRSETzHIVacfelmnkjhbg
329 selection = np.logical_and(selection, sync['times'] >= tmin) if tmin else selection 1vutsJKFGdpqriNUWXYOPQDRSETzHIVacfelmnkjhbg
330 return Bunch({'times': sync['times'][selection], 1vutsJKFGdpqriNUWXYOPQDRSETzHIVacfelmnkjhbg
331 'polarities': sync['polarities'][selection]})
334def _clean_audio(audio, display=False):
335 """
336 one guy wired the 150 Hz camera output onto the soundcard. The effect is to get 150 Hz periodic
337 square pulses, 2ms up and 4.666 ms down. When this happens we remove all of the intermediate
338 pulses to repair the audio trace
339 Here is some helper code
340 dd = np.diff(audio['times'])
341 1 / np.median(dd[::2]) # 2ms up
342 1 / np.median(dd[1::2]) # 4.666 ms down
343 1 / (np.median(dd[::2]) + np.median(dd[1::2])) # both sum to 150 Hz
344 This only runs on sessions when the bug is detected and leaves others untouched
345 """
346 DISCARD_THRESHOLD = 0.01 1v1diacfehbg
347 average_150_hz = np.mean(1 / np.diff(audio['times'][audio['polarities'] == 1]) > 140) 1v1diacfehbg
348 naudio = audio['times'].size 1v1diacfehbg
349 if average_150_hz > 0.7 and naudio > 100: 1v1diacfehbg
350 _logger.warning('Soundcard signal on FPGA seems to have been mixed with 150Hz camera') 11
351 keep_ind = np.r_[np.diff(audio['times']) > DISCARD_THRESHOLD, False] 11
352 keep_ind = np.logical_and(keep_ind, audio['polarities'] == -1) 11
353 keep_ind = np.where(keep_ind)[0] 11
354 keep_ind = np.sort(np.r_[0, keep_ind, keep_ind + 1, naudio - 1]) 11
356 if display: # pragma: no cover 11
357 from ibllib.plots import squares
358 squares(audio['times'], audio['polarities'], ax=None, yrange=[-1, 1])
359 squares(audio['times'][keep_ind], audio['polarities'][keep_ind], yrange=[-1, 1])
360 audio = {'times': audio['times'][keep_ind], 11
361 'polarities': audio['polarities'][keep_ind]}
362 return audio 1v1diacfehbg
365def _clean_frame2ttl(frame2ttl, threshold=0.01, display=False):
366 """
367 Clean the frame2ttl events.
369 Frame 2ttl calibration can be unstable and the fronts may be flickering at an unrealistic
370 pace. This removes the consecutive frame2ttl pulses happening too fast, below a threshold
371 of F2TTL_THRESH.
373 Parameters
374 ----------
375 frame2ttl : dict
376 A dictionary of frame2TTL events, with keys {'times', 'polarities'}.
377 threshold : float
378 Consecutive pulses occurring with this many seconds ignored.
379 display : bool
380 If true, plots the input TTLs and the cleaned output.
382 Returns
383 -------
385 """
386 dt = np.diff(frame2ttl['times']) 16FGd8izHIacfeh9!#7bg
387 iko = np.where(np.logical_and(dt < threshold, frame2ttl['polarities'][:-1] == -1))[0] 16FGd8izHIacfeh9!#7bg
388 iko = np.unique(np.r_[iko, iko + 1]) 16FGd8izHIacfeh9!#7bg
389 frame2ttl_ = {'times': np.delete(frame2ttl['times'], iko), 16FGd8izHIacfeh9!#7bg
390 'polarities': np.delete(frame2ttl['polarities'], iko)}
391 if iko.size > (0.1 * frame2ttl['times'].size): 16FGd8izHIacfeh9!#7bg
392 _logger.warning(f'{iko.size} ({iko.size / frame2ttl["times"].size:.2%}) ' 16ic7
393 f'frame to TTL polarity switches below {threshold} secs')
394 if display: # pragma: no cover 16FGd8izHIacfeh9!#7bg
395 fig, (ax0, ax1) = plt.subplots(2, sharex=True)
396 plots.squares(frame2ttl['times'] * 1000, frame2ttl['polarities'], yrange=[0.1, 0.9], ax=ax0)
397 plots.squares(frame2ttl_['times'] * 1000, frame2ttl_['polarities'], yrange=[1.1, 1.9], ax=ax1)
398 import seaborn as sns
399 sns.displot(dt[dt < 0.05], binwidth=0.0005)
401 return frame2ttl_ 16FGd8izHIacfeh9!#7bg
404def extract_wheel_sync(sync, chmap=None, tmin=None, tmax=None):
405 """
406 Extract wheel positions and times from sync fronts dictionary for all 16 channels.
407 Output position is in radians, mathematical convention.
409 Parameters
410 ----------
411 sync : dict
412 'polarities' of fronts detected on sync trace for all 16 chans and their 'times'
413 chmap : dict
414 Map of channel names and their corresponding index. Default to constant.
415 tmin : float
416 The minimum time from which to extract the sync pulses.
417 tmax : float
418 The maximum time up to which we extract the sync pulses.
420 Returns
421 -------
422 numpy.array
423 Wheel timestamps in seconds.
424 numpy.array
425 Wheel positions in radians.
426 """
427 # Assume two separate edge count channels
428 assert chmap.keys() >= {'rotary_encoder_0', 'rotary_encoder_1'} 1sdacfelmnkjhbg
429 channela = get_sync_fronts(sync, chmap['rotary_encoder_0'], tmin=tmin, tmax=tmax) 1sdacfelmnkjhbg
430 channelb = get_sync_fronts(sync, chmap['rotary_encoder_1'], tmin=tmin, tmax=tmax) 1sdacfelmnkjhbg
431 re_ts, re_pos = _rotary_encoder_positions_from_fronts( 1sdacfelmnkjhbg
432 channela['times'], channela['polarities'], channelb['times'], channelb['polarities'],
433 ticks=WHEEL_TICKS, radius=WHEEL_RADIUS_CM, coding='x4')
434 return re_ts, re_pos 1sdacfelmnkjhbg
437def extract_sync(session_path, overwrite=False, ephys_files=None, namespace='spikeglx'):
438 """
439 Reads ephys binary file (s) and extract sync within the binary file folder
440 Assumes ephys data is within a `raw_ephys_data` folder
442 :param session_path: '/path/to/subject/yyyy-mm-dd/001'
443 :param overwrite: Bool on re-extraction, forces overwrite instead of loading existing files
444 :return: list of sync dictionaries
445 """
446 session_path = Path(session_path) 1wxypqrABC0lmnkj
447 if not ephys_files: 1wxypqrABC0lmnkj
448 ephys_files = spikeglx.glob_ephys_files(session_path) 1wxy0lmnkj
449 syncs = [] 1wxypqrABC0lmnkj
450 outputs = [] 1wxypqrABC0lmnkj
451 for efi in ephys_files: 1wxypqrABC0lmnkj
452 bin_file = efi.get('ap', efi.get('nidq', None)) 1wxypqrABClmnkj
453 if not bin_file: 1wxypqrABClmnkj
454 continue
455 alfname = dict(object='sync', namespace=namespace) 1wxypqrABClmnkj
456 if efi.label: 1wxypqrABClmnkj
457 alfname['extra'] = efi.label 1pqrlmnkj
458 file_exists = alfio.exists(bin_file.parent, **alfname) 1wxypqrABClmnkj
459 if not overwrite and file_exists: 1wxypqrABClmnkj
460 _logger.warning(f'Skipping raw sync: SGLX sync found for {efi.label}!') 1wxylmnk
461 sync = alfio.load_object(bin_file.parent, **alfname) 1wxylmnk
462 out_files, _ = alfio._ls(bin_file.parent, **alfname) 1wxylmnk
463 else:
464 sr = spikeglx.Reader(bin_file) 1wxypqrABCj
465 sync, out_files = _sync_to_alf(sr, bin_file.parent, save=True, parts=efi.label) 1wxypqrABCj
466 sr.close() 1wxypqrABCj
467 outputs.extend(out_files) 1wxypqrABClmnkj
468 syncs.extend([sync]) 1wxypqrABClmnkj
470 return syncs, outputs 1wxypqrABC0lmnkj
473def _get_all_probes_sync(session_path, bin_exists=True):
474 # round-up of all bin ephys files in the session, infer revision and get sync map
475 ephys_files = spikeglx.glob_ephys_files(session_path, bin_exists=bin_exists) 1sJKDEac0lmnkjb
476 version = spikeglx.get_neuropixel_version_from_files(ephys_files) 1sJKDEac0lmnkjb
477 # attach the sync information to each binary file found
478 for ef in ephys_files: 1sJKDEac0lmnkjb
479 ef['sync'] = alfio.load_object(ef.path, 'sync', namespace='spikeglx', short_keys=True) 1sJKDEaclmnkjb
480 ef['sync_map'] = get_ibl_sync_map(ef, version) 1sJKDEaclmnkjb
481 return ephys_files 1sJKDEac0lmnkjb
484def get_wheel_positions(sync, chmap, tmin=None, tmax=None):
485 """
486 Gets the wheel position from synchronisation pulses
488 Parameters
489 ----------
490 sync : dict
491 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses and
492 the corresponding channel numbers.
493 chmap : dict[str, int]
494 A map of channel names and their corresponding indices.
495 tmin : float
496 The minimum time from which to extract the sync pulses.
497 tmax : float
498 The maximum time up to which we extract the sync pulses.
500 Returns
501 -------
502 Bunch
503 A dictionary with keys ('timestamps', 'position'), containing the wheel event timestamps and
504 position in radians
505 Bunch
506 A dictionary of detected movement times with keys ('intervals', 'peakAmplitude', 'peakVelocity_times').
507 """
508 ts, pos = extract_wheel_sync(sync=sync, chmap=chmap, tmin=tmin, tmax=tmax) 1sdacfehbg
509 moves = Bunch(extract_wheel_moves(ts, pos)) 1sdacfehbg
510 wheel = Bunch({'timestamps': ts, 'position': pos}) 1sdacfehbg
511 return wheel, moves 1sdacfehbg
514def get_main_probe_sync(session_path, bin_exists=False):
515 """
516 From 3A or 3B multiprobe session, returns the main probe (3A) or nidq sync pulses
517 with the attached channel map (default chmap if none)
519 Parameters
520 ----------
521 session_path : str, pathlib.Path
522 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
523 bin_exists : bool
524 Whether there is a .bin file present.
526 Returns
527 -------
528 one.alf.io.AlfBunch
529 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses and
530 the corresponding channel numbers.
531 dict
532 A map of channel names and their corresponding indices.
533 """
534 ephys_files = _get_all_probes_sync(session_path, bin_exists=bin_exists) 1sJKDEac0lmnkjb
535 if not ephys_files: 1sJKDEac0lmnkjb
536 raise FileNotFoundError(f"No ephys files found in {session_path}") 10
537 version = spikeglx.get_neuropixel_version_from_files(ephys_files) 1sJKDEaclmnkjb
538 if version == '3A': 1sJKDEaclmnkjb
539 # the sync master is the probe with the most sync pulses
540 sync_box_ind = np.argmax([ef.sync.times.size for ef in ephys_files]) 1aclmnb
541 elif version == '3B': 1sJKDEkj
542 # the sync master is the nidq breakout box
543 sync_box_ind = np.argmax([1 if ef.get('nidq') else 0 for ef in ephys_files]) 1sJKDEkj
544 sync = ephys_files[sync_box_ind].sync 1sJKDEaclmnkjb
545 sync_chmap = ephys_files[sync_box_ind].sync_map 1sJKDEaclmnkjb
546 return sync, sync_chmap 1sJKDEaclmnkjb
549def get_protocol_period(session_path, protocol_number, bpod_sync):
550 """
552 Parameters
553 ----------
554 session_path : str, pathlib.Path
555 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
556 protocol_number : int
557 The order that the protocol was run in.
558 bpod_sync : dict
559 The sync times and polarities for Bpod BNC1.
561 Returns
562 -------
563 float
564 The time of the detected spacer for the protocol number.
565 float, None
566 The time of the next detected spacer or None if this is the last protocol run.
567 """
568 # The spacers are TTLs generated by Bpod at the start of each protocol
569 spacer_times = Spacer().find_spacers_from_fronts(bpod_sync) 1z
570 # Ensure that the number of detected spacers matched the number of expected tasks
571 if acquisition_description := session_params.read_params(session_path): 1z
572 n_tasks = len(acquisition_description.get('tasks', [])) 1z
573 assert n_tasks == len(spacer_times), f'expected {n_tasks} spacers, found {len(spacer_times)}' 1z
574 assert n_tasks > protocol_number >= 0, f'protocol number must be between 0 and {n_tasks}' 1z
575 else:
576 assert protocol_number < len(spacer_times)
577 start = spacer_times[int(protocol_number)] 1z
578 end = None if len(spacer_times) - 1 == protocol_number else spacer_times[int(protocol_number + 1)] 1z
579 return start, end 1z
582class FpgaTrials(extractors_base.BaseExtractor):
583 save_names = ('_ibl_trials.goCueTrigger_times.npy', None, None, None, None, None, None, None,
584 '_ibl_trials.stimOff_times.npy', None, None, None, '_ibl_trials.quiescencePeriod.npy',
585 '_ibl_trials.table.pqt', '_ibl_wheel.timestamps.npy',
586 '_ibl_wheel.position.npy', '_ibl_wheelMoves.intervals.npy',
587 '_ibl_wheelMoves.peakAmplitude.npy')
588 var_names = ('goCueTrigger_times', 'stimOnTrigger_times',
589 'stimOffTrigger_times', 'stimFreezeTrigger_times', 'errorCueTrigger_times',
590 'errorCue_times', 'itiIn_times', 'stimFreeze_times', 'stimOff_times',
591 'valveOpen_times', 'phase', 'position', 'quiescence', 'table',
592 'wheel_timestamps', 'wheel_position',
593 'wheelMoves_intervals', 'wheelMoves_peakAmplitude')
595 bpod_rsync_fields = ('intervals', 'response_times', 'goCueTrigger_times',
596 'stimOnTrigger_times', 'stimOffTrigger_times',
597 'stimFreezeTrigger_times', 'errorCueTrigger_times')
598 """tuple of str: Fields from Bpod extractor that we want to re-sync to FPGA."""
600 bpod_fields = ('feedbackType', 'choice', 'rewardVolume', 'contrastLeft', 'contrastRight',
601 'probabilityLeft', 'phase', 'position', 'quiescence')
602 """tuple of str: Fields from bpod extractor that we want to save."""
604 sync_field = 'intervals_0' # trial start events
605 """str: The trial event to synchronize (must be present in extracted trials)."""
607 bpod = None
608 """dict of numpy.array: The Bpod out TTLs recorded on the DAQ. Used in the QC viewer plot."""
610 def __init__(self, *args, bpod_trials=None, bpod_extractor=None, **kwargs):
611 """An extractor for ephysChoiceWorld trials data, in FPGA time.
613 This class may be subclassed to handle moderate variations in hardware and task protocol,
614 however there is flexible
615 """
616 super().__init__(*args, **kwargs) 15vutd'(iacfelmnkjhbg
617 self.bpod2fpga = None 15vutd'(iacfelmnkjhbg
618 self.bpod_trials = bpod_trials 15vutd'(iacfelmnkjhbg
619 self.frame2ttl = self.audio = self.bpod = self.settings = None 15vutd'(iacfelmnkjhbg
620 if bpod_extractor: 15vutd'(iacfelmnkjhbg
621 self.bpod_extractor = bpod_extractor 1diacfebg
622 self._update_var_names() 1diacfebg
624 def _update_var_names(self, bpod_fields=None, bpod_rsync_fields=None):
625 """
626 Updates this object's attributes based on the Bpod trials extractor.
628 Fields updated: bpod_fields, bpod_rsync_fields, save_names, and var_names.
630 Parameters
631 ----------
632 bpod_fields : tuple
633 A set of Bpod trials fields to keep.
634 bpod_rsync_fields : tuple
635 A set of Bpod trials fields to sync to the DAQ times.
636 """
637 if self.bpod_extractor: 1diacfebg
638 for var_name, save_name in zip(self.bpod_extractor.var_names, self.bpod_extractor.save_names): 1diacfebg
639 if var_name not in self.var_names: 1diacfebg
640 self.var_names += (var_name,) 1diacfebg
641 self.save_names += (save_name,) 1diacfebg
643 # self.var_names = self.bpod_extractor.var_names
644 # self.save_names = self.bpod_extractor.save_names
645 self.settings = self.bpod_extractor.settings # This is used by the TaskQC 1diacfebg
646 self.bpod_rsync_fields = bpod_rsync_fields 1diacfebg
647 if self.bpod_rsync_fields is None: 1diacfebg
648 self.bpod_rsync_fields = tuple(self._time_fields(self.bpod_extractor.var_names)) 1diacfebg
649 if 'table' in self.bpod_extractor.var_names: 1diacfebg
650 if not self.bpod_trials: 1diacfebg
651 self.bpod_trials = self.bpod_extractor.extract(save=False)
652 table_keys = alfio.AlfBunch.from_df(self.bpod_trials['table']).keys() 1diacfebg
653 self.bpod_rsync_fields += tuple(self._time_fields(table_keys)) 1diacfebg
654 elif bpod_rsync_fields:
655 self.bpod_rsync_fields = bpod_rsync_fields
656 excluded = (*self.bpod_rsync_fields, 'table') 1diacfebg
657 if bpod_fields: 1diacfebg
658 assert not set(self.bpod_fields).intersection(excluded), 'bpod_fields must not also be bpod_rsync_fields'
659 self.bpod_fields = bpod_fields
660 elif self.bpod_extractor: 1diacfebg
661 self.bpod_fields = tuple(x for x in self.bpod_extractor.var_names if x not in excluded) 1diacfebg
662 if 'table' in self.bpod_extractor.var_names: 1diacfebg
663 if not self.bpod_trials: 1diacfebg
664 self.bpod_trials = self.bpod_extractor.extract(save=False)
665 table_keys = alfio.AlfBunch.from_df(self.bpod_trials['table']).keys() 1diacfebg
666 self.bpod_fields += tuple([x for x in table_keys if x not in excluded]) 1diacfebg
668 @staticmethod
669 def _time_fields(trials_attr) -> set:
670 """
671 Iterates over Bpod trials attributes returning those that correspond to times for syncing.
673 Parameters
674 ----------
675 trials_attr : iterable of str
676 The Bpod field names.
678 Returns
679 -------
680 set
681 The field names that contain timestamps.
682 """
683 FIELDS = ('times', 'timestamps', 'intervals') 1)diacfebg
684 pattern = re.compile(fr'^[_\w]*({"|".join(FIELDS)})[_\w]*$') 1)diacfebg
685 return set(filter(pattern.match, trials_attr)) 1)diacfebg
687 def load_sync(self, sync_collection='raw_ephys_data', **kwargs):
688 """Load the DAQ sync and channel map data.
690 This method may be subclassed for novel DAQ systems. The sync must contain the following
691 keys: 'times' - an array timestamps in seconds; 'polarities' - an array of {-1, 1}
692 corresponding to TTL LOW and TTL HIGH, respectively; 'channels' - an array of ints
693 corresponding to channel number.
695 Parameters
696 ----------
697 sync_collection : str
698 The session subdirectory where the sync data are located.
699 kwargs
700 Optional arguments used by subclass methods.
702 Returns
703 -------
704 one.alf.io.AlfBunch
705 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
706 and the corresponding channel numbers.
707 dict
708 A map of channel names and their corresponding indices.
709 """
710 return get_sync_and_chn_map(self.session_path, sync_collection)
712 def _extract(self, sync=None, chmap=None, sync_collection='raw_ephys_data',
713 task_collection='raw_behavior_data', **kwargs) -> dict:
714 """Extracts ephys trials by combining Bpod and FPGA sync pulses.
716 It is essential that the `var_names`, `bpod_rsync_fields`, `bpod_fields`, and `sync_field`
717 attributes are all correct for the bpod protocol used.
719 Below are the steps involved:
720 0. Load sync and bpod trials, if required.
721 1. Determine protocol period and discard sync events outside the task.
722 2. Classify multiplexed TTL events based on length (see :meth:`FpgaTrials.build_trials`).
723 3. Sync the Bpod clock to the DAQ clock using one of the assigned trial events.
724 4. Assign classified TTL events to trial events based on order within the trial.
725 4. Convert Bpod software event times to DAQ clock.
726 5. Extract the wheel from the DAQ rotary encoder signal, if required.
728 Parameters
729 ----------
730 sync : dict
731 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
732 and the corresponding channel numbers. If None, the sync is loaded using the
733 `load_sync` method.
734 chmap : dict
735 A map of channel names and their corresponding indices. If None, the channel map is
736 loaded using the :meth:`FpgaTrials.load_sync` method.
737 sync_collection : str
738 The session subdirectory where the sync data are located. This is only used if the
739 sync or channel maps are not provided.
740 task_collection : str
741 The session subdirectory where the raw Bpod data are located. This is used for loading
742 the task settings and extracting the bpod trials, if not already done.
743 protocol_number : int
744 The protocol number if multiple protocols were run during the session. If provided, a
745 spacer signal must be present in order to determine the correct period.
746 kwargs
747 Optional arguments for subclass methods to use.
749 Returns
750 -------
751 dict
752 A dictionary of numpy arrays with `FpgaTrials.var_names` as keys.
753 """
754 if sync is None or chmap is None: 1diacfehbg
755 _sync, _chmap = self.load_sync(sync_collection) 1i
756 sync = sync or _sync 1i
757 chmap = chmap or _chmap 1i
759 if not self.bpod_trials: # extract the behaviour data from bpod 1diacfehbg
760 self.bpod_trials, *_ = bpod_extract_all(
761 session_path=self.session_path, task_collection=task_collection, save=False,
762 extractor_type=kwargs.get('extractor_type'))
764 # Explode trials table df
765 if 'table' in self.var_names: 1diacfehbg
766 trials_table = alfio.AlfBunch.from_df(self.bpod_trials.pop('table')) 1diacfehbg
767 table_columns = trials_table.keys() 1diacfehbg
768 self.bpod_trials.update(trials_table) 1diacfehbg
769 else:
770 if 'table' in self.bpod_trials: 1a
771 _logger.error(
772 '"table" found in Bpod trials but missing from `var_names` attribute and will'
773 'therefore not be extracted. This is likely in error.')
774 table_columns = None 1a
776 bpod = get_sync_fronts(sync, chmap['bpod']) 1diacfehbg
777 # Get the spacer times for this protocol
778 if any(arg in kwargs for arg in ('tmin', 'tmax')): 1diacfehbg
779 tmin, tmax = kwargs.get('tmin'), kwargs.get('tmax')
780 elif (protocol_number := kwargs.get('protocol_number')) is not None: # look for spacer 1diacfehbg
781 # The spacers are TTLs generated by Bpod at the start of each protocol
782 tmin, tmax = get_protocol_period(self.session_path, protocol_number, bpod)
783 else:
784 # Older sessions don't have protocol spacers so we sync the Bpod intervals here to
785 # find the approximate end time of the protocol (this will exclude the passive signals
786 # in ephysChoiceWorld that tend to ruin the final trial extraction).
787 _, trial_ints = self.get_bpod_event_times(sync, chmap, **kwargs) 1diacfehbg
788 t_trial_start = trial_ints.get('trial_start', np.array([[np.nan, np.nan]]))[:, 0] 1diacfehbg
789 bpod_start = self.bpod_trials['intervals'][:, 0] 1diacfehbg
790 if len(t_trial_start) > len(bpod_start) / 2: # if least half the trial start TTLs detected 1diacfehbg
791 _logger.warning('Attempting to get protocol period from aligning trial start TTLs') 1dacfehbg
792 fcn, *_ = ibldsp.utils.sync_timestamps(bpod_start, t_trial_start) 1dacfehbg
793 buffer = 2.5 # the number of seconds to include before/after task 1dacfehbg
794 start, end = fcn(self.bpod_trials['intervals'].flat[[0, -1]]) 1dacfehbg
795 tmin = min(sync['times'][0], start - buffer) 1dacfehbg
796 tmax = max(sync['times'][-1], end + buffer) 1dacfehbg
797 else: # This type of alignment fails for some sessions, e.g. mesoscope
798 tmin = tmax = None 1ia
800 # Remove unnecessary data from sync
801 selection = np.logical_and( 1diacfehbg
802 sync['times'] <= (tmax if tmax is not None else sync['times'][-1]),
803 sync['times'] >= (tmin if tmin is not None else sync['times'][0]),
804 )
805 sync = alfio.AlfBunch({k: v[selection] for k, v in sync.items()}) 1diacfehbg
806 _logger.debug('Protocol period from %.2fs to %.2fs (~%.0f min duration)', 1diacfehbg
807 *sync['times'][[0, -1]], np.diff(sync['times'][[0, -1]]) / 60)
809 # Get the trial events from the DAQ sync TTLs, sync clocks and build final trials datasets
810 out = self.build_trials(sync=sync, chmap=chmap, **kwargs) 1diacfehbg
812 # extract the wheel data
813 if any(x.startswith('wheel') for x in self.var_names): 1diacfehbg
814 wheel, moves = self.get_wheel_positions(sync=sync, chmap=chmap, tmin=tmin, tmax=tmax) 1diacfehbg
815 from ibllib.io.extractors.training_wheel import extract_first_movement_times 1diacfehbg
816 if not self.settings: 1diacfehbg
817 self.settings = raw.load_settings(session_path=self.session_path, task_collection=task_collection) 1h
818 min_qt = self.settings.get('QUIESCENT_PERIOD', None) 1diacfehbg
819 first_move_onsets, *_ = extract_first_movement_times(moves, out, min_qt=min_qt) 1diacfehbg
820 out.update({'firstMovement_times': first_move_onsets}) 1diacfehbg
821 out.update({f'wheel_{k}': v for k, v in wheel.items()}) 1diacfehbg
822 out.update({f'wheelMoves_{k}': v for k, v in moves.items()}) 1diacfehbg
824 # Re-create trials table
825 if table_columns: 1diacfehbg
826 trials_table = alfio.AlfBunch({x: out.pop(x) for x in table_columns}) 1diacfehbg
827 out['table'] = trials_table.to_df() 1diacfehbg
829 out = alfio.AlfBunch({k: out[k] for k in self.var_names if k in out}) # Reorder output 1diacfehbg
830 assert self.var_names == tuple(out.keys()) 1diacfehbg
831 return out 1diacfehbg
833 def _is_trials_object_attribute(self, var_name, variable_length_vars=None):
834 """
835 Check if variable name is expected to have the same length as trials.intervals.
837 Parameters
838 ----------
839 var_name : str
840 The variable name to check.
841 variable_length_vars : list
842 Set of variable names that are not expected to have the same length as trials.intervals.
843 This list may be passed by superclasses.
845 Returns
846 -------
847 bool
848 True if variable is a trials dataset.
850 Examples
851 --------
852 >>> assert self._is_trials_object_attribute('stimOnTrigger_times') is True
853 >>> assert self._is_trials_object_attribute('wheel_position') is False
854 """
855 save_name = self.save_names[self.var_names.index(var_name)] if var_name in self.var_names else None 15dacfehbg
856 if save_name: 15dacfehbg
857 return filename_parts(save_name)[1] == 'trials' 15dacfehbg
858 else:
859 return var_name not in (variable_length_vars or []) 15dacfehbg
861 def build_trials(self, sync, chmap, display=False, **kwargs):
862 """
863 Extract task related event times from the sync.
865 The trial start times are the shortest Bpod TTLs and occur at the start of the trial. The
866 first trial start TTL of the session is longer and must be handled differently. The trial
867 start TTL is used to assign the other trial events to each trial.
869 The trial end is the end of the so-called 'ITI' Bpod event TTL (classified as the longest
870 of the three Bpod event TTLs). Go cue audio TTLs are the shorter of the two expected audio
871 tones. The first of these after each trial start is taken to be the go cue time. Error
872 tones are longer audio TTLs and assigned as the last of such occurrence after each trial
873 start. The valve open Bpod TTLs are medium-length, the last of which is used for each trial.
874 The feedback times are times of either valve open or error tone as there should be only one
875 such event per trial.
877 The stimulus times are taken from the frame2ttl events (with improbably high frequency TTLs
878 removed): the first TTL after each trial start is assumed to be the stim onset time; the
879 second to last and last are taken as the stimulus freeze and offset times, respectively.
881 Parameters
882 ----------
883 sync : dict
884 'polarities' of fronts detected on sync trace for all 16 chans and their 'times'
885 chmap : dict
886 Map of channel names and their corresponding index. Default to constant.
887 display : bool, matplotlib.pyplot.Axes
888 Show the full session sync pulses display.
890 Returns
891 -------
892 dict
893 A map of trial event timestamps.
894 """
895 # Get the events from the sync.
896 # Store the cleaned frame2ttl, audio, and bpod pulses as this will be used for QC
897 self.frame2ttl = self.get_stimulus_update_times(sync, chmap, **kwargs) 1dacfehbg
898 self.audio, audio_event_intervals = self.get_audio_event_times(sync, chmap, **kwargs) 1dacfehbg
899 if not set(audio_event_intervals.keys()) >= {'ready_tone', 'error_tone'}: 1dacfehbg
900 raise ValueError(
901 'Expected at least "ready_tone" and "error_tone" audio events.'
902 '`audio_event_ttls` kwarg may be incorrect.')
903 self.bpod, bpod_event_intervals = self.get_bpod_event_times(sync, chmap, **kwargs) 1dacfehbg
904 if not set(bpod_event_intervals.keys()) >= {'trial_start', 'valve_open', 'trial_end'}: 1dacfehbg
905 raise ValueError(
906 'Expected at least "trial_start", "trial_end", and "valve_open" audio events. '
907 '`bpod_event_ttls` kwarg may be incorrect.')
909 t_iti_in, t_trial_end = bpod_event_intervals['trial_end'].T 1dacfehbg
910 fpga_events = alfio.AlfBunch({ 1dacfehbg
911 'goCue_times': audio_event_intervals['ready_tone'][:, 0],
912 'errorCue_times': audio_event_intervals['error_tone'][:, 0],
913 'valveOpen_times': bpod_event_intervals['valve_open'][:, 0],
914 'valveClose_times': bpod_event_intervals['valve_open'][:, 1],
915 'itiIn_times': t_iti_in,
916 'intervals_0': bpod_event_intervals['trial_start'][:, 0],
917 'intervals_1': t_trial_end
918 })
920 # Sync the Bpod clock to the DAQ.
921 # NB: The Bpod extractor typically drops the final, incomplete, trial. Hence there is
922 # usually at least one extra FPGA event. This shouldn't affect the sync. The final trial is
923 # dropped after assigning the FPGA events, using the `ifpga` index. Doing this after
924 # assigning the FPGA trial events ensures the last trial has the correct timestamps.
925 self.bpod2fpga, drift_ppm, ibpod, ifpga = self.sync_bpod_clock(self.bpod_trials, fpga_events, self.sync_field) 1dacfehbg
927 if np.any(np.diff(ibpod) != 1) and self.sync_field == 'intervals_0': 1dacfehbg
928 # One issue is that sometimes pulses may not have been detected, in this case
929 # add the events that have not been detected and re-extract the behaviour sync.
930 # This is only really relevant for the Bpod interval events as the other TTLs are
931 # from devices where a missing TTL likely means the Bpod event was truly absent.
932 _logger.warning('Missing Bpod TTLs; reassigning events using aligned Bpod start times')
933 bpod_start = self.bpod_trials['intervals'][:, 0]
934 missing_bpod = self.bpod2fpga(bpod_start[np.setxor1d(ibpod, np.arange(len(bpod_start)))])
935 t_trial_start = np.sort(np.r_[fpga_events['intervals_0'][:, 0], missing_bpod])
936 else:
937 t_trial_start = fpga_events['intervals_0'] 1dacfehbg
939 out = alfio.AlfBunch() 1dacfehbg
940 # Add the Bpod trial events, converting the timestamp fields to FPGA time.
941 # NB: The trial intervals are by default a Bpod rsync field.
942 out.update({k: self.bpod_trials[k][ibpod] for k in self.bpod_fields}) 1dacfehbg
943 for k in self.bpod_rsync_fields: 1dacfehbg
944 # Some personal projects may extract non-trials object datasets that may not have 1 event per trial
945 idx = ibpod if self._is_trials_object_attribute(k) else np.arange(len(self.bpod_trials[k]), dtype=int) 1dacfehbg
946 out[k] = self.bpod2fpga(self.bpod_trials[k][idx]) 1dacfehbg
948 f2ttl_t = self.frame2ttl['times'] 1dacfehbg
949 # Assign the FPGA events to individual trials
950 fpga_trials = { 1dacfehbg
951 'goCue_times': _assign_events_to_trial(t_trial_start, fpga_events['goCue_times'], take='first'),
952 'errorCue_times': _assign_events_to_trial(t_trial_start, fpga_events['errorCue_times']),
953 'valveOpen_times': _assign_events_to_trial(t_trial_start, fpga_events['valveOpen_times']),
954 'itiIn_times': _assign_events_to_trial(t_trial_start, fpga_events['itiIn_times']),
955 'stimOn_times': np.full_like(t_trial_start, np.nan),
956 'stimOff_times': np.full_like(t_trial_start, np.nan),
957 'stimFreeze_times': np.full_like(t_trial_start, np.nan)
958 }
960 # f2ttl times are unreliable owing to calibration and Bonsai sync square update issues.
961 # Take the first event after the FPGA aligned stimulus trigger time.
962 fpga_trials['stimOn_times'][ibpod] = _assign_events_to_trial( 1dacfehbg
963 out['stimOnTrigger_times'], f2ttl_t, take='first', t_trial_end=out['stimOffTrigger_times'])
964 fpga_trials['stimOff_times'][ibpod] = _assign_events_to_trial( 1dacfehbg
965 out['stimOffTrigger_times'], f2ttl_t, take='first', t_trial_end=out['intervals'][:, 1])
966 # For stim freeze we take the last event before the stim off trigger time.
967 # To avoid assigning early events (e.g. for sessions where there are few flips due to
968 # mis-calibration), we discount events before stim freeze trigger times (or stim on trigger
969 # times for versions below 6.2.5). We take the last event rather than the first after stim
970 # freeze trigger because often there are multiple flips after the trigger, presumably
971 # before the stim actually stops.
972 stim_freeze = np.copy(out['stimFreezeTrigger_times']) 1dacfehbg
973 go_trials = np.where(out['choice'] != 0)[0] 1dacfehbg
974 # NB: versions below 6.2.5 have no trigger times so use stim on trigger times
975 lims = np.copy(out['stimOnTrigger_times']) 1dacfehbg
976 if not np.isnan(stim_freeze).all(): 1dacfehbg
977 # Stim freeze times are NaN for nogo trials, but for all others use stim freeze trigger
978 # times. _assign_events_to_trial requires ascending timestamps so no NaNs allowed.
979 lims[go_trials] = stim_freeze[go_trials] 1dcehb
980 # take last event after freeze/stim on trigger, before stim off trigger
981 stim_freeze = _assign_events_to_trial(lims, f2ttl_t, take='last', t_trial_end=out['stimOffTrigger_times']) 1dacfehbg
982 fpga_trials['stimFreeze_times'][go_trials] = stim_freeze[go_trials] 1dacfehbg
984 # Feedback times are valve open on correct trials and error tone in on incorrect trials
985 fpga_trials['feedback_times'] = np.copy(fpga_trials['valveOpen_times']) 1dacfehbg
986 ind_err = np.isnan(fpga_trials['valveOpen_times']) 1dacfehbg
987 fpga_trials['feedback_times'][ind_err] = fpga_trials['errorCue_times'][ind_err] 1dacfehbg
989 out.update({k: fpga_trials[k][ifpga] for k in fpga_trials.keys()}) 1dacfehbg
991 if display: # pragma: no cover 1dacfehbg
992 width = 0.5
993 ymax = 5
994 if isinstance(display, bool):
995 plt.figure('Bpod FPGA Sync')
996 ax = plt.gca()
997 else:
998 ax = display
999 plots.squares(self.bpod['times'], self.bpod['polarities'] * 0.4 + 1, ax=ax, color='k')
1000 plots.squares(self.frame2ttl['times'], self.frame2ttl['polarities'] * 0.4 + 2, ax=ax, color='k')
1001 plots.squares(self.audio['times'], self.audio['polarities'] * 0.4 + 3, ax=ax, color='k')
1002 color_map = TABLEAU_COLORS.keys()
1003 for (event_name, event_times), c in zip(fpga_events.items(), cycle(color_map)):
1004 plots.vertical_lines(event_times, ymin=0, ymax=ymax, ax=ax, color=c, label=event_name, linewidth=width)
1005 # Plot the stimulus events along with the trigger times
1006 stim_events = filter(lambda t: 'stim' in t[0], fpga_trials.items())
1007 for (event_name, event_times), c in zip(stim_events, cycle(color_map)):
1008 plots.vertical_lines(
1009 event_times, ymin=0, ymax=ymax, ax=ax, color=c, label=event_name, linewidth=width, linestyle='--')
1010 nm = event_name.replace('_times', 'Trigger_times')
1011 plots.vertical_lines(
1012 out[nm], ymin=0, ymax=ymax, ax=ax, color=c, label=nm, linewidth=width, linestyle=':')
1013 ax.legend()
1014 ax.set_yticks([0, 1, 2, 3])
1015 ax.set_yticklabels(['', 'bpod', 'f2ttl', 'audio'])
1016 ax.set_ylim([0, 5])
1017 return out 1dacfehbg
1019 def get_wheel_positions(self, *args, **kwargs):
1020 """Extract wheel and wheelMoves objects.
1022 This method is called by the main extract method and may be overloaded by subclasses.
1023 """
1024 return get_wheel_positions(*args, **kwargs) 1dacfehbg
1026 def get_stimulus_update_times(self, sync, chmap, display=False, **_):
1027 """
1028 Extract stimulus update times from sync.
1030 Gets the stimulus times from the frame2ttl channel and cleans the signal.
1032 Parameters
1033 ----------
1034 sync : dict
1035 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
1036 and the corresponding channel numbers.
1037 chmap : dict
1038 A map of channel names and their corresponding indices. Must contain a 'frame2ttl' key.
1039 display : bool
1040 If true, plots the input TTLs and the cleaned output.
1042 Returns
1043 -------
1044 dict
1045 A dictionary with keys {'times', 'polarities'} containing stimulus TTL fronts.
1046 """
1047 frame2ttl = get_sync_fronts(sync, chmap['frame2ttl']) 1diacfehbg
1048 frame2ttl = _clean_frame2ttl(frame2ttl, display=display) 1diacfehbg
1049 return frame2ttl 1diacfehbg
1051 def get_audio_event_times(self, sync, chmap, audio_event_ttls=None, display=False, **_):
1052 """
1053 Extract audio times from sync.
1055 Gets the TTL times from the 'audio' channel, cleans the signal, and classifies each TTL
1056 event by length.
1058 Parameters
1059 ----------
1060 sync : dict
1061 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
1062 and the corresponding channel numbers.
1063 chmap : dict
1064 A map of channel names and their corresponding indices. Must contain an 'audio' key.
1065 audio_event_ttls : dict
1066 A map of event names to (min, max) TTL length.
1067 display : bool
1068 If true, plots the input TTLs and the cleaned output.
1070 Returns
1071 -------
1072 dict
1073 A dictionary with keys {'times', 'polarities'} containing audio TTL fronts.
1074 dict
1075 A dictionary of events (from `audio_event_ttls`) and their intervals as an Nx2 array.
1076 """
1077 audio = get_sync_fronts(sync, chmap['audio']) 1vdiacfehbg
1078 audio = _clean_audio(audio) 1vdiacfehbg
1080 if audio['times'].size == 0: 1vdiacfehbg
1081 _logger.error('No audio sync fronts found.')
1083 if audio_event_ttls is None: 1vdiacfehbg
1084 # For training/biased/ephys protocols, the ready tone should be below 110 ms. The error
1085 # tone should be between 400ms and 1200ms
1086 audio_event_ttls = {'ready_tone': (0, 0.11), 'error_tone': (0.4, 1.2)} 1vdiacfehbg
1087 audio_event_intervals = self._assign_events(audio['times'], audio['polarities'], audio_event_ttls, display=display) 1vdiacfehbg
1089 return audio, audio_event_intervals 1vdiacfehbg
1091 def get_bpod_event_times(self, sync, chmap, bpod_event_ttls=None, display=False, **kwargs):
1092 """
1093 Extract Bpod times from sync.
1095 Gets the Bpod TTL times from the sync 'bpod' channel and classifies each TTL event by
1096 length. NB: The first trial has an abnormal trial_start TTL that is usually mis-assigned.
1097 This method accounts for this.
1099 Parameters
1100 ----------
1101 sync : dict
1102 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
1103 and the corresponding channel numbers. Must contain a 'bpod' key.
1104 chmap : dict
1105 A map of channel names and their corresponding indices.
1106 bpod_event_ttls : dict of tuple
1107 A map of event names to (min, max) TTL length.
1109 Returns
1110 -------
1111 dict
1112 A dictionary with keys {'times', 'polarities'} containing Bpod TTL fronts.
1113 dict
1114 A dictionary of events (from `bpod_event_ttls`) and their intervals as an Nx2 array.
1115 """
1116 bpod = get_sync_fronts(sync, chmap['bpod']) 1utdiacfelmnkjhbg
1117 if bpod.times.size == 0: 1utdiacfelmnkjhbg
1118 raise err.SyncBpodFpgaException('No Bpod event found in FPGA. No behaviour extraction. '
1119 'Check channel maps.')
1120 # Assign the Bpod BNC2 events based on TTL length. The defaults are below, however these
1121 # lengths are defined by the state machine of the task protocol and therefore vary.
1122 if bpod_event_ttls is None: 1utdiacfelmnkjhbg
1123 # For training/biased/ephys protocols, the trial start TTL length is 0.1ms but this has
1124 # proven to drift on some Bpods and this is the highest possible value that
1125 # discriminates trial start from valve. Valve open events are between 50ms to 300 ms.
1126 # ITI events are above 400 ms.
1127 bpod_event_ttls = { 1utdacfelmnkjhbg
1128 'trial_start': (0, 2.33e-4), 'valve_open': (2.33e-4, 0.4), 'trial_end': (0.4, np.inf)}
1129 bpod_event_intervals = self._assign_events( 1utdiacfelmnkjhbg
1130 bpod['times'], bpod['polarities'], bpod_event_ttls, display=display)
1132 if 'trial_start' not in bpod_event_intervals or bpod_event_intervals['trial_start'].size == 0: 1utdiacfelmnkjhbg
1133 return bpod, bpod_event_intervals 1i
1135 # The first trial pulse is longer and often assigned to another event.
1136 # Here we move the earliest non-trial_start event to the trial_start array.
1137 t0 = bpod_event_intervals['trial_start'][0, 0] # expect 1st event to be trial_start 1utdacfelmnkjhbg
1138 pretrial = [(k, v[0, 0]) for k, v in bpod_event_intervals.items() if v.size and v[0, 0] < t0] 1utdacfelmnkjhbg
1139 if pretrial: 1utdacfelmnkjhbg
1140 (pretrial, _) = sorted(pretrial, key=lambda x: x[1])[0] # take the earliest event 1utdacfelmnkjhbg
1141 dt = np.diff(bpod_event_intervals[pretrial][0, :]) * 1e3 # record TTL length to log 1utdacfelmnkjhbg
1142 _logger.debug('Reassigning first %s to trial_start. TTL length = %.3g ms', pretrial, dt) 1utdacfelmnkjhbg
1143 bpod_event_intervals['trial_start'] = np.r_[ 1utdacfelmnkjhbg
1144 bpod_event_intervals[pretrial][0:1, :], bpod_event_intervals['trial_start']
1145 ]
1146 bpod_event_intervals[pretrial] = bpod_event_intervals[pretrial][1:, :] 1utdacfelmnkjhbg
1148 return bpod, bpod_event_intervals 1utdacfelmnkjhbg
1150 @staticmethod
1151 def _assign_events(ts, polarities, event_lengths, precedence='shortest', display=False):
1152 """
1153 Classify TTL events by length.
1155 Outputs the synchronisation events such as trial intervals, valve opening, and audio.
1157 Parameters
1158 ----------
1159 ts : numpy.array
1160 Numpy vector containing times of TTL fronts.
1161 polarities : numpy.array
1162 Numpy vector containing polarity of TTL fronts (1 rise, -1 fall).
1163 event_lengths : dict of tuple
1164 A map of TTL events and the range of permissible lengths, where l0 < ttl <= l1.
1165 precedence : str {'shortest', 'longest', 'dict order'}
1166 In the case of overlapping event TTL lengths, assign shortest/longest first or go by
1167 the `event_lengths` dict order.
1168 display : bool
1169 If true, plots the TTLs with coloured lines delineating the assigned events.
1171 Returns
1172 -------
1173 Dict[str, numpy.array]
1174 A dictionary of events and their intervals as an Nx2 array.
1176 See Also
1177 --------
1178 _assign_events_to_trial - classify TTLs by event order within a given trial period.
1179 """
1180 event_intervals = dict.fromkeys(event_lengths) 1vutdiacfelmnkjhbg
1181 assert 'unassigned' not in event_lengths.keys() 1vutdiacfelmnkjhbg
1183 if len(ts) == 0: 1vutdiacfelmnkjhbg
1184 return {k: np.array([[], []]).T for k in (*event_lengths.keys(), 'unassigned')}
1186 # make sure that there are no 2 consecutive fall or consecutive rise events
1187 assert np.all(np.abs(np.diff(polarities)) == 2) 1vutdiacfelmnkjhbg
1188 if polarities[0] == -1: 1vutdiacfelmnkjhbg
1189 ts = np.delete(ts, 0) 1vt
1190 if polarities[-1] == 1: # if the final TTL is left HIGH, insert a NaN 1vutdiacfelmnkjhbg
1191 ts = np.r_[ts, np.nan] 1t
1192 # take only even time differences: i.e. from rising to falling fronts
1193 dt = np.diff(ts)[::2] 1vutdiacfelmnkjhbg
1195 # Assign events from shortest TTL to largest
1196 assigned = np.zeros(ts.shape, dtype=bool) 1vutdiacfelmnkjhbg
1197 if precedence.lower() == 'shortest': 1vutdiacfelmnkjhbg
1198 event_items = sorted(event_lengths.items(), key=lambda x: np.diff(x[1])) 1vutdiacfelmnkjhbg
1199 elif precedence.lower() == 'longest':
1200 event_items = sorted(event_lengths.items(), key=lambda x: np.diff(x[1]), reverse=True)
1201 elif precedence.lower() == 'dict order':
1202 event_items = event_lengths.items()
1203 else:
1204 raise ValueError(f'Precedence must be one of "shortest", "longest", "dict order", got "{precedence}".')
1205 for event, (min_len, max_len) in event_items: 1vutdiacfelmnkjhbg
1206 _logger.debug('%s: %.4G < ttl <= %.4G', event, min_len, max_len) 1vutdiacfelmnkjhbg
1207 i_event = np.where(np.logical_and(dt > min_len, dt <= max_len))[0] * 2 1vutdiacfelmnkjhbg
1208 i_event = i_event[np.where(~assigned[i_event])[0]] # remove those already assigned 1vutdiacfelmnkjhbg
1209 event_intervals[event] = np.c_[ts[i_event], ts[i_event + 1]] 1vutdiacfelmnkjhbg
1210 assigned[np.r_[i_event, i_event + 1]] = True 1vutdiacfelmnkjhbg
1212 # Include the unassigned events for convenience and debugging
1213 event_intervals['unassigned'] = ts[~assigned].reshape(-1, 2) 1vutdiacfelmnkjhbg
1215 # Assert that event TTLs mutually exclusive
1216 all_assigned = np.concatenate(list(event_intervals.values())).flatten() 1vutdiacfelmnkjhbg
1217 assert all_assigned.size == np.unique(all_assigned).size, 'TTLs assigned to multiple events' 1vutdiacfelmnkjhbg
1219 # some debug plots when needed
1220 if display: # pragma: no cover 1vutdiacfelmnkjhbg
1221 plt.figure()
1222 plots.squares(ts, polarities, label='raw fronts')
1223 for event, intervals in event_intervals.items():
1224 plots.vertical_lines(intervals[:, 0], ymin=-0.2, ymax=1.1, linewidth=0.5, label=event)
1225 plt.legend()
1227 # Return map of event intervals in the same order as `event_lengths` dict
1228 return {k: event_intervals[k] for k in (*event_lengths, 'unassigned')} 1vutdiacfelmnkjhbg
1230 @staticmethod
1231 def sync_bpod_clock(bpod_trials, fpga_trials, sync_field):
1232 """
1233 Sync the Bpod clock to FPGA one using the provided trial event.
1235 It assumes that `sync_field` is in both `fpga_trials` and `bpod_trials`. Syncing on both
1236 intervals is not supported so to sync on trial start times, `sync_field` should be
1237 'intervals_0'.
1239 Parameters
1240 ----------
1241 bpod_trials : dict
1242 A dictionary of extracted Bpod trial events.
1243 fpga_trials : dict
1244 A dictionary of TTL events extracted from FPGA sync (see `extract_behaviour_sync`
1245 method).
1246 sync_field : str
1247 The trials key to use for syncing clocks. For intervals (i.e. Nx2 arrays) append the
1248 column index, e.g. 'intervals_0'.
1250 Returns
1251 -------
1252 function
1253 Interpolation function such that f(timestamps_bpod) = timestamps_fpga.
1254 float
1255 The clock drift in parts per million.
1256 numpy.array of int
1257 The indices of the Bpod trial events in the FPGA trial events array.
1258 numpy.array of int
1259 The indices of the FPGA trial events in the Bpod trial events array.
1261 Raises
1262 ------
1263 ValueError
1264 The key `sync_field` was not found in either the `bpod_trials` or `fpga_trials` dicts.
1265 """
1266 _logger.info(f'Attempting to align Bpod clock to DAQ using trial event "{sync_field}"') 1diacfehbg
1267 bpod_fpga_timestamps = [None, None] 1diacfehbg
1268 for i, trials in enumerate((bpod_trials, fpga_trials)): 1diacfehbg
1269 if sync_field not in trials: 1diacfehbg
1270 # handle syncing on intervals
1271 if not (m := re.match(r'(.*)_(\d)', sync_field)): 1dacfehbg
1272 # If missing from bpod trials, either the sync field is incorrect,
1273 # or the Bpod extractor is incorrect. If missing from the fpga events, check
1274 # the sync field and the `extract_behaviour_sync` method.
1275 raise ValueError(
1276 f'Sync field "{sync_field}" not in extracted {"fpga" if i else "bpod"} events')
1277 _sync_field, n = m.groups() 1dacfehbg
1278 bpod_fpga_timestamps[i] = trials[_sync_field][:, int(n)] 1dacfehbg
1279 else:
1280 bpod_fpga_timestamps[i] = trials[sync_field] 1diacfehbg
1282 # Sync the two timestamps
1283 fcn, drift, ibpod, ifpga = ibldsp.utils.sync_timestamps(*bpod_fpga_timestamps, return_indices=True) 1diacfehbg
1285 # If it's drifting too much throw warning or error
1286 _logger.info('N trials: %i bpod, %i FPGA, %i merged, sync %.5f ppm', 1diacfehbg
1287 *map(len, bpod_fpga_timestamps), len(ibpod), drift)
1288 if drift > 200 and bpod_fpga_timestamps[0].size != bpod_fpga_timestamps[1].size: 1diacfehbg
1289 raise err.SyncBpodFpgaException('sync cluster f*ck')
1290 elif drift > BPOD_FPGA_DRIFT_THRESHOLD_PPM: 1diacfehbg
1291 _logger.warning('BPOD/FPGA synchronization shows values greater than %.2f ppm',
1292 BPOD_FPGA_DRIFT_THRESHOLD_PPM)
1294 return fcn, drift, ibpod, ifpga 1diacfehbg
1297class FpgaTrialsHabituation(FpgaTrials):
1298 """Extract habituationChoiceWorld trial events from an NI DAQ."""
1300 save_names = ('_ibl_trials.stimCenter_times.npy', '_ibl_trials.feedbackType.npy', '_ibl_trials.rewardVolume.npy',
1301 '_ibl_trials.stimOff_times.npy', '_ibl_trials.contrastLeft.npy', '_ibl_trials.contrastRight.npy',
1302 '_ibl_trials.feedback_times.npy', '_ibl_trials.stimOn_times.npy', '_ibl_trials.stimOnTrigger_times.npy',
1303 '_ibl_trials.intervals.npy', '_ibl_trials.goCue_times.npy', '_ibl_trials.goCueTrigger_times.npy',
1304 None, None, None, None, None)
1305 """tuple of str: The filenames of each extracted dataset, or None if array should not be saved."""
1307 var_names = ('stimCenter_times', 'feedbackType', 'rewardVolume', 'stimOff_times', 'contrastLeft',
1308 'contrastRight', 'feedback_times', 'stimOn_times', 'stimOnTrigger_times', 'intervals',
1309 'goCue_times', 'goCueTrigger_times', 'itiIn_times', 'stimOffTrigger_times',
1310 'stimCenterTrigger_times', 'position', 'phase')
1311 """tuple of str: A list of names for the extracted variables. These become the returned output keys."""
1313 bpod_rsync_fields = ('intervals', 'stimOn_times', 'feedback_times', 'stimCenterTrigger_times',
1314 'goCue_times', 'itiIn_times', 'stimOffTrigger_times', 'stimOff_times',
1315 'stimCenter_times', 'stimOnTrigger_times', 'goCueTrigger_times')
1316 """tuple of str: Fields from Bpod extractor that we want to re-sync to FPGA."""
1318 bpod_fields = ('feedbackType', 'rewardVolume', 'contrastLeft', 'contrastRight', 'position', 'phase')
1319 """tuple of str: Fields from Bpod extractor that we want to save."""
1321 sync_field = 'feedback_times' # valve open events
1322 """str: The trial event to synchronize (must be present in extracted trials)."""
1324 def _extract(self, sync=None, chmap=None, sync_collection='raw_ephys_data',
1325 task_collection='raw_behavior_data', **kwargs) -> dict:
1326 """
1327 Extract habituationChoiceWorld trial events from an NI DAQ.
1329 It is essential that the `var_names`, `bpod_rsync_fields`, `bpod_fields`, and `sync_field`
1330 attributes are all correct for the bpod protocol used.
1332 Unlike FpgaTrials, this class assumes different Bpod TTL events and syncs the Bpod clock
1333 using the valve open times, instead of the trial start times.
1335 Parameters
1336 ----------
1337 sync : dict
1338 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
1339 and the corresponding channel numbers. If None, the sync is loaded using the
1340 `load_sync` method.
1341 dict
1342 A map of channel names and their corresponding indices. If None, the channel map is
1343 loaded using the `load_sync` method.
1344 sync_collection : str
1345 The session subdirectory where the sync data are located. This is only used if the
1346 sync or channel maps are not provided.
1347 task_collection : str
1348 The session subdirectory where the raw Bpod data are located. This is used for loading
1349 the task settings and extracting the bpod trials, if not already done.
1350 protocol_number : int
1351 The protocol number if multiple protocols were run during the session. If provided, a
1352 spacer signal must be present in order to determine the correct period.
1353 kwargs
1354 Optional arguments for class methods, e.g. 'display', 'bpod_event_ttls'.
1356 Returns
1357 -------
1358 dict
1359 A dictionary of numpy arrays with `FpgaTrialsHabituation.var_names` as keys.
1360 """
1361 # Version check: the ITI in TTL was added in a later version
1362 if not self.settings: 1a
1363 self.settings = raw.load_settings(session_path=self.session_path, task_collection=task_collection)
1364 iblrig_version = version.parse(self.settings.get('IBL_VERSION', '0.0.0')) 1a
1365 if version.parse('8.9.3') <= iblrig_version < version.parse('8.12.6'): 1a
1366 """A second 1s TTL was added in this version during the 'iti' state, however this is
1367 unrelated to the trial ITI and is unfortunately the same length as the trial start TTL."""
1368 raise NotImplementedError('Ambiguous TTLs in 8.9.3 >= version < 8.12.6')
1370 trials = super()._extract(sync=sync, chmap=chmap, sync_collection=sync_collection, 1a
1371 task_collection=task_collection, **kwargs)
1373 return trials 1a
1375 def get_bpod_event_times(self, sync, chmap, bpod_event_ttls=None, display=False, **kwargs):
1376 """
1377 Extract Bpod times from sync.
1379 Currently (at least v8.12 and below) there is no trial start or end TTL, only an ITI pulse.
1380 Also the first trial pulse is incorrectly assigned due to its abnormal length.
1382 Parameters
1383 ----------
1384 sync : dict
1385 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses
1386 and the corresponding channel numbers. Must contain a 'bpod' key.
1387 chmap : dict
1388 A map of channel names and their corresponding indices.
1389 bpod_event_ttls : dict of tuple
1390 A map of event names to (min, max) TTL length.
1392 Returns
1393 -------
1394 dict
1395 A dictionary with keys {'times', 'polarities'} containing Bpod TTL fronts.
1396 dict
1397 A dictionary of events (from `bpod_event_ttls`) and their intervals as an Nx2 array.
1398 """
1399 bpod = get_sync_fronts(sync, chmap['bpod']) 1a
1400 if bpod.times.size == 0: 1a
1401 raise err.SyncBpodFpgaException('No Bpod event found in FPGA. No behaviour extraction. '
1402 'Check channel maps.')
1403 # Assign the Bpod BNC2 events based on TTL length. The defaults are below, however these
1404 # lengths are defined by the state machine of the task protocol and therefore vary.
1405 if bpod_event_ttls is None: 1a
1406 # Currently (at least v8.12 and below) there is no trial start or end TTL, only an ITI pulse
1407 bpod_event_ttls = {'trial_iti': (1, 1.1), 'valve_open': (0, 0.4)} 1a
1408 bpod_event_intervals = self._assign_events( 1a
1409 bpod['times'], bpod['polarities'], bpod_event_ttls, display=display)
1411 # The first trial pulse is shorter and assigned to valve_open. Here we remove the first
1412 # valve event, prepend a 0 to the trial_start events, and drop the last trial if it was
1413 # incomplete in Bpod.
1414 bpod_event_intervals['trial_iti'] = np.r_[bpod_event_intervals['valve_open'][0:1, :], 1a
1415 bpod_event_intervals['trial_iti']]
1416 bpod_event_intervals['valve_open'] = bpod_event_intervals['valve_open'][1:, :] 1a
1418 return bpod, bpod_event_intervals 1a
1420 def build_trials(self, sync, chmap, display=False, **kwargs):
1421 """
1422 Extract task related event times from the sync.
1424 This is called by the superclass `_extract` method. The key difference here is that the
1425 `trial_start` LOW->HIGH is the trial end, and HIGH->LOW is trial start.
1427 Parameters
1428 ----------
1429 sync : dict
1430 'polarities' of fronts detected on sync trace for all 16 chans and their 'times'
1431 chmap : dict
1432 Map of channel names and their corresponding index. Default to constant.
1433 display : bool, matplotlib.pyplot.Axes
1434 Show the full session sync pulses display.
1436 Returns
1437 -------
1438 dict
1439 A map of trial event timestamps.
1440 """
1441 # Get the events from the sync.
1442 # Store the cleaned frame2ttl, audio, and bpod pulses as this will be used for QC
1443 self.frame2ttl = self.get_stimulus_update_times(sync, chmap, **kwargs) 1a
1444 self.audio, audio_event_intervals = self.get_audio_event_times(sync, chmap, **kwargs) 1a
1445 self.bpod, bpod_event_intervals = self.get_bpod_event_times(sync, chmap, **kwargs) 1a
1446 if not set(bpod_event_intervals.keys()) >= {'valve_open', 'trial_iti'}: 1a
1447 raise ValueError(
1448 'Expected at least "trial_iti" and "valve_open" Bpod events. `bpod_event_ttls` kwarg may be incorrect.')
1450 fpga_events = alfio.AlfBunch({ 1a
1451 'feedback_times': bpod_event_intervals['valve_open'][:, 0],
1452 'valveClose_times': bpod_event_intervals['valve_open'][:, 1],
1453 'intervals_0': bpod_event_intervals['trial_iti'][:, 1],
1454 'intervals_1': bpod_event_intervals['trial_iti'][:, 0],
1455 'goCue_times': audio_event_intervals['ready_tone'][:, 0]
1456 })
1458 # Sync the Bpod clock to the DAQ.
1459 self.bpod2fpga, drift_ppm, ibpod, ifpga = self.sync_bpod_clock(self.bpod_trials, fpga_events, self.sync_field) 1a
1461 out = alfio.AlfBunch() 1a
1462 # Add the Bpod trial events, converting the timestamp fields to FPGA time.
1463 # NB: The trial intervals are by default a Bpod rsync field.
1464 out.update({k: self.bpod_trials[k][ibpod] for k in self.bpod_fields}) 1a
1465 out.update({k: self.bpod2fpga(self.bpod_trials[k][ibpod]) for k in self.bpod_rsync_fields}) 1a
1467 # Assigning each event to a trial ensures exactly one event per trial (missing events are NaN)
1468 assign_to_trial = partial(_assign_events_to_trial, fpga_events['intervals_0']) 1a
1469 trials = alfio.AlfBunch({ 1a
1470 'goCue_times': assign_to_trial(fpga_events['goCue_times'], take='first'),
1471 'feedback_times': assign_to_trial(fpga_events['feedback_times']),
1472 'stimCenter_times': assign_to_trial(self.frame2ttl['times'], take=-2),
1473 'stimOn_times': assign_to_trial(self.frame2ttl['times'], take='first'),
1474 'stimOff_times': assign_to_trial(self.frame2ttl['times']),
1475 })
1476 out.update({k: trials[k][ifpga] for k in trials.keys()}) 1a
1478 # If stim on occurs before trial end, use stim on time. Likewise for trial end and stim off
1479 to_correct = ~np.isnan(out['stimOn_times']) & (out['stimOn_times'] < out['intervals'][:, 0]) 1a
1480 if np.any(to_correct): 1a
1481 _logger.warning('%i/%i stim on events occurring outside trial intervals', sum(to_correct), len(to_correct))
1482 out['intervals'][to_correct, 0] = out['stimOn_times'][to_correct]
1483 to_correct = ~np.isnan(out['stimOff_times']) & (out['stimOff_times'] > out['intervals'][:, 1]) 1a
1484 if np.any(to_correct): 1a
1485 _logger.debug( 1a
1486 '%i/%i stim off events occurring outside trial intervals; using stim off times as trial end',
1487 sum(to_correct), len(to_correct))
1488 out['intervals'][to_correct, 1] = out['stimOff_times'][to_correct] 1a
1490 if display: # pragma: no cover 1a
1491 width = 0.5
1492 ymax = 5
1493 if isinstance(display, bool):
1494 plt.figure('Bpod FPGA Sync')
1495 ax = plt.gca()
1496 else:
1497 ax = display
1498 plots.squares(self.bpod['times'], self.bpod['polarities'] * 0.4 + 1, ax=ax, color='k')
1499 plots.squares(self.frame2ttl['times'], self.frame2ttl['polarities'] * 0.4 + 2, ax=ax, color='k')
1500 plots.squares(self.audio['times'], self.audio['polarities'] * 0.4 + 3, ax=ax, color='k')
1501 color_map = TABLEAU_COLORS.keys()
1502 for (event_name, event_times), c in zip(trials.to_df().items(), cycle(color_map)):
1503 plots.vertical_lines(event_times, ymin=0, ymax=ymax, ax=ax, color=c, label=event_name, linewidth=width)
1504 ax.legend()
1505 ax.set_yticks([0, 1, 2, 3])
1506 ax.set_yticklabels(['', 'bpod', 'f2ttl', 'audio'])
1507 ax.set_ylim([0, 4])
1509 return out 1a
1512def extract_all(session_path, sync_collection='raw_ephys_data', save=True, save_path=None,
1513 task_collection='raw_behavior_data', protocol_number=None, **kwargs):
1514 """
1515 For the IBL ephys task, reads ephys binary file and extract:
1516 - sync
1517 - wheel
1518 - behaviour
1520 These `extract_all` functions should be deprecated as they make assumptions about hardware
1521 parameters. Additionally the FpgaTrials class now automatically loads DAQ sync files, extracts
1522 the Bpod trials, and returns a dict instead of a tuple. Therefore this function is entirely
1523 redundant. See the examples for the correct way to extract NI DAQ behaviour sessions.
1525 Parameters
1526 ----------
1527 session_path : str, pathlib.Path
1528 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
1529 sync_collection : str
1530 The session subdirectory where the sync data are located.
1531 save : bool
1532 If true, save the extracted files to save_path.
1533 task_collection : str
1534 The location of the behaviour protocol data.
1535 save_path : str, pathlib.Path
1536 The save location of the extracted files, defaults to the alf directory of the session path.
1537 protocol_number : int
1538 The order that the protocol was run in.
1539 **kwargs
1540 Optional extractor keyword arguments.
1542 Returns
1543 -------
1544 list
1545 The extracted data.
1546 list of pathlib.Path, None
1547 If save is True, a list of file paths to the extracted data.
1548 """
1549 warnings.warn( 1h
1550 'ibllib.io.extractors.ephys_fpga.extract_all will be removed in future versions; '
1551 'use FpgaTrials instead. For reliable extraction, use the dynamic pipeline behaviour tasks.',
1552 FutureWarning)
1553 return_extractor = kwargs.pop('return_extractor', False) 1h
1554 # Extract Bpod trials
1555 bpod_raw = raw.load_data(session_path, task_collection=task_collection) 1h
1556 assert bpod_raw is not None, 'No task trials data in raw_behavior_data - Exit' 1h
1557 bpod_trials, bpod_wheel, *_ = bpod_extract_all( 1h
1558 session_path=session_path, bpod_trials=bpod_raw, task_collection=task_collection,
1559 save=False, extractor_type=kwargs.get('extractor_type'))
1561 # Sync Bpod trials to FPGA
1562 sync, chmap = get_sync_and_chn_map(session_path, sync_collection) 1h
1563 # sync, chmap = get_main_probe_sync(session_path, bin_exists=bin_exists)
1564 trials = FpgaTrials(session_path, bpod_trials={**bpod_trials, **bpod_wheel}) # py3.9 -> | 1h
1565 outputs, files = trials.extract( 1h
1566 save=save, sync=sync, chmap=chmap, path_out=save_path,
1567 task_collection=task_collection, protocol_number=protocol_number, **kwargs)
1568 if not isinstance(outputs, dict): 1h
1569 outputs = {k: v for k, v in zip(trials.var_names, outputs)}
1570 if return_extractor: 1h
1571 return outputs, files, trials 1h
1572 else:
1573 return outputs, files
1576def get_sync_and_chn_map(session_path, sync_collection):
1577 """
1578 Return sync and channel map for session based on collection where main sync is stored.
1580 Parameters
1581 ----------
1582 session_path : str, pathlib.Path
1583 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
1584 sync_collection : str
1585 The session subdirectory where the sync data are located.
1587 Returns
1588 -------
1589 one.alf.io.AlfBunch
1590 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses and
1591 the corresponding channel numbers.
1592 dict
1593 A map of channel names and their corresponding indices.
1594 """
1595 if sync_collection == 'raw_ephys_data': 1FGdNUWXYOPQRSTzHIacfehbg
1596 # Check to see if we have nidq files, if we do just go with this otherwise go into other function that deals with
1597 # 3A probes
1598 nidq_meta = next(session_path.joinpath(sync_collection).glob('*nidq.meta'), None) 1FGdNOPQRSTzHIacfehbg
1599 if not nidq_meta: 1FGdNOPQRSTzHIacfehbg
1600 sync, chmap = get_main_probe_sync(session_path) 1acb
1601 else:
1602 sync = load_sync(session_path, sync_collection) 1FGdNOPQRSTzHIafehbg
1603 ef = Bunch() 1FGdNOPQRSTzHIafehbg
1604 ef['path'] = session_path.joinpath(sync_collection) 1FGdNOPQRSTzHIafehbg
1605 ef['nidq'] = nidq_meta 1FGdNOPQRSTzHIafehbg
1606 chmap = get_ibl_sync_map(ef, '3B') 1FGdNOPQRSTzHIafehbg
1608 else:
1609 sync = load_sync(session_path, sync_collection) 1UWXY
1610 chmap = load_channel_map(session_path, sync_collection) 1UWXY
1612 return sync, chmap 1FGdNUWXYOPQRSTzHIacfehbg
1615def load_channel_map(session_path, sync_collection):
1616 """
1617 Load syncing channel map for session path and collection
1619 Parameters
1620 ----------
1621 session_path : str, pathlib.Path
1622 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
1623 sync_collection : str
1624 The session subdirectory where the sync data are located.
1626 Returns
1627 -------
1628 dict
1629 A map of channel names and their corresponding indices.
1630 """
1632 device = sync_collection.split('_')[1] 1UWXY
1633 default_chmap = DEFAULT_MAPS[device]['nidq'] 1UWXY
1635 # Try to load channel map from file
1636 chmap = spikeglx.get_sync_map(session_path.joinpath(sync_collection)) 1UWXY
1637 # If chmap provided but not with all keys, fill up with default values
1638 if not chmap: 1UWXY
1639 return default_chmap 1WXY
1640 else:
1641 if data_for_keys(default_chmap.keys(), chmap): 1U
1642 return chmap 1U
1643 else:
1644 _logger.warning('Keys missing from provided channel map, '
1645 'setting missing keys from default channel map')
1646 return {**default_chmap, **chmap}
1649def load_sync(session_path, sync_collection):
1650 """
1651 Load sync files from session path and collection.
1653 Parameters
1654 ----------
1655 session_path : str, pathlib.Path
1656 The absolute session path, i.e. '/path/to/subject/yyyy-mm-dd/nnn'.
1657 sync_collection : str
1658 The session subdirectory where the sync data are located.
1660 Returns
1661 -------
1662 one.alf.io.AlfBunch
1663 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses and
1664 the corresponding channel numbers.
1665 """
1666 sync = alfio.load_object(session_path.joinpath(sync_collection), 'sync', namespace='spikeglx', short_keys=True) 1FGdNUWXYOPQRSTzHIafehbg
1668 return sync 1FGdNUWXYOPQRSTzHIafehbg