Coverage for ibllib/io/extractors/mesoscope.py: 85%

1"""Mesoscope (timeline) data extraction.""" 

2import logging (empty)

3from itertools import chain (empty)

4 

5import numpy as np (empty)

6from scipy.signal import find_peaks (empty)

7import one.alf.io as alfio (empty)

8from one.alf.path import session_path_parts (empty)

9from iblutil.util import ensure_list (empty)

10import matplotlib.pyplot as plt (empty)

11from packaging import version (empty)

12 

13from ibllib.plots.misc import squares, vertical_lines (empty)

14from ibllib.io.raw_daq_loaders import ( (empty)

15 extract_sync_timeline, timeline_get_channel, correct_counter_discontinuities, load_timeline_sync_and_chmap) 

16import ibllib.io.extractors.base as extractors_base (empty)

17from ibllib.io.extractors.ephys_fpga import ( (empty)

18 FpgaTrials, FpgaTrialsHabituation, WHEEL_TICKS, WHEEL_RADIUS_CM, _assign_events_to_trial) 

19from ibllib.io.extractors.training_wheel import extract_wheel_moves (empty)

20from ibllib.io.extractors.camera import attribute_times (empty)

21from brainbox.behavior.wheel import velocity_filtered (empty)

22 

23_logger = logging.getLogger(__name__) (empty)

24 

25 

26def patch_imaging_meta(meta: dict) -> dict: (empty)

27 """ 

28 Patch imaging metadata for compatibility across versions. 

29 

30 A copy of the dict is NOT returned. 

31 

32 Parameters 

33 ---------- 

34 meta : dict 

35 A folder path that contains a rawImagingData.meta file. 

36 

37 Returns 

38 ------- 

39 dict 

40 The loaded metadata file, updated to the most recent version. 

41 """ 

42 # 2023-05-17 (unversioned) adds nFrames, channelSaved keys, MM and Deg keys 

43 ver = version.parse(meta.get('version') or '0.0.0') 7 ctx1jnoldbc

44 if ver <= version.parse('0.0.0'): 7 ctx1jnoldbc

45 if 'channelSaved' not in meta: 3 ctx1jno

46 meta['channelSaved'] = next((x['channelIdx'] for x in meta['FOV'] if 'channelIdx' in x), []) 3 ctx1jno

47 fields = ('topLeft', 'topRight', 'bottomLeft', 'bottomRight') 3 ctx1jno

48 for fov in meta.get('FOV', []): 3 ctx1jno

49 for unit in ('Deg', 'MM'): 3 ctx1jno

50 if unit not in fov: # topLeftDeg, etc. -> Deg[topLeft] 3 ctx1jno

51 fov[unit] = {f: fov.pop(f + unit, None) for f in fields} 3 ctx1jno

52 elif ver == version.parse('0.1.0'): 5 ctx1jldbc

53 for fov in meta.get('FOV', []): 1 ctx1j

54 if 'roiUuid' in fov: 1 ctx1j

55 fov['roiUUID'] = fov.pop('roiUuid') 1 ctx1j

56 # 2024-09-17 Modified the 2 unit vectors for the positive ML axis and the positive AP axis, 

57 # which then transform [X,Y] coordinates (in degrees) to [ML,AP] coordinates (in MM). 

58 if ver < version.Version('0.1.5') and 'imageOrientation' in meta: 7 ctx1jnoldbc

59 pos_ml, pos_ap = meta['imageOrientation']['positiveML'], meta['imageOrientation']['positiveAP'] 4 ctx1jlbc

60 center_ml, center_ap = meta['centerMM']['ML'], meta['centerMM']['AP'] 4 ctx1jlbc

61 res = meta['scanImageParams']['objectiveResolution'] 4 ctx1jlbc

62 # previously [[0, res/1000], [-res/1000, 0], [0, 0]] 

63 TF = np.linalg.pinv(np.c_[np.vstack([pos_ml, pos_ap, [0, 0]]), [1, 1, 1]]) @ \ 4 ctx1jlbc

64 (np.array([[res / 1000, 0], [0, res / 1000], [0, 0]]) + np.array([center_ml, center_ap])) 

65 TF = np.round(TF, 3) # handle floating-point error by rounding 4 ctx1jlbc

66 if not np.allclose(TF, meta['coordsTF']): 4 ctx1jlbc

67 meta['coordsTF'] = TF.tolist() 4 ctx1jlbc

68 centerDegXY = np.array([meta['centerDeg']['x'], meta['centerDeg']['y']]) 4 ctx1jlbc

69 for fov in meta.get('FOV', []): 4 ctx1jlbc

70 fov['MM'] = {k: (np.r_[np.array(v) - centerDegXY, 1] @ TF).tolist() for k, v in fov['Deg'].items()} 4 ctx1jlbc

71 

72 assert 'nFrames' in meta, '"nFrames" key missing from meta data; rawImagingData.meta.json likely an old version' 7 ctx1jnoldbc

73 return meta 7 ctx1jnoldbc

74 

75 

76def plot_timeline(timeline, channels=None, raw=True): (empty)

77 """ 

78 Plot the timeline data. 

79 

80 Parameters 

81 ---------- 

82 timeline : one.alf.io.AlfBunch 

83 The timeline data object. 

84 channels : list of str 

85 An iterable of channel names to plot. 

86 raw : bool 

87 If true, plot the raw DAQ samples; if false, apply TTL thresholds and plot changes. 

88 

89 Returns 

90 ------- 

91 matplotlib.pyplot.Figure 

92 The figure containing timeline subplots. 

93 list of matplotlib.pyplot.Axes 

94 The axes for each timeline channel plotted. 

95 """ 

96 meta = {x.copy().pop('name'): x for x in timeline['meta']['inputs']} 1 ctx1k

97 channels = channels or meta.keys() 1 ctx1k

98 fig, axes = plt.subplots(len(channels), 1, sharex=True) 1 ctx1k

99 axes = ensure_list(axes) 1 ctx1k

100 if not raw: 1 ctx1k

101 chmap = {ch: meta[ch]['arrayColumn'] for ch in channels} 1 ctx1k

102 sync = extract_sync_timeline(timeline, chmap=chmap) 1 ctx1k

103 for i, (ax, ch) in enumerate(zip(axes, channels)): 1 ctx1k

104 if raw: 1 ctx1k

105 # axesScale controls vertical scaling of each trace (multiplicative) 

106 values = timeline['raw'][:, meta[ch]['arrayColumn'] - 1] * meta[ch]['axesScale'] 1 ctx1k

107 ax.plot(timeline['timestamps'], values) 1 ctx1k

108 elif np.any(idx := sync['channels'] == chmap[ch]): 1 ctx1k

109 squares(sync['times'][idx], sync['polarities'][idx], ax=ax) 1 ctx1k

110 ax.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) 1 ctx1k

111 ax.spines['bottom'].set_visible(False), ax.spines['left'].set_visible(True) 1 ctx1k

112 ax.set_ylabel(ch, rotation=45, fontsize=8) 1 ctx1k

113 # Add back x-axis ticks to the last plot 

114 axes[-1].tick_params(axis='x', which='both', bottom=True, labelbottom=True) 1 ctx1k

115 axes[-1].spines['bottom'].set_visible(True) 1 ctx1k

116 plt.get_current_fig_manager().window.showMaximized() # full screen 1 ctx1k

117 fig.tight_layout(h_pad=0) 1 ctx1k

118 return fig, axes 1 ctx1k

119 

120 

121class TimelineTrials(FpgaTrials): (empty)

122 """Similar extraction to the FPGA, however counter and position channels are treated differently.""" 

123 

124 timeline = None (empty)

125 """one.alf.io.AlfBunch: The timeline data object.""" (empty)

126 

127 sync_field = 'itiIn_times' (empty)

128 """str: The trial event to synchronize (must be present in extracted trials).""" (empty)

129 

130 def __init__(self, *args, sync_collection='raw_sync_data', **kwargs): (empty)

131 """An extractor for all ephys trial data, in Timeline time""" 

132 super().__init__(*args, **kwargs) 4 ctx1giae

133 self.timeline = alfio.load_object(self.session_path / sync_collection, 'DAQdata', namespace='timeline') 4 ctx1giae

134 

135 def load_sync(self, sync_collection='raw_sync_data', chmap=None, **_): (empty)

136 """Load the DAQ sync and channel map data. 

137 

138 Parameters 

139 ---------- 

140 sync_collection : str 

141 The session subdirectory where the sync data are located. 

142 chmap : dict 

143 A map of channel names and their corresponding indices. If None, the channel map is 

144 loaded using the :func:`ibllib.io.raw_daq_loaders.timeline_meta2chmap` method. 

145 

146 Returns 

147 ------- 

148 one.alf.io.AlfBunch 

149 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses 

150 and the corresponding channel numbers. 

151 dict 

152 A map of channel names and their corresponding indices. 

153 """ 

154 if not self.timeline: 2 ctx1ae

155 self.timeline = alfio.load_object(self.session_path / sync_collection, 'DAQdata', namespace='timeline') 

156 sync, chmap = load_timeline_sync_and_chmap( 2 ctx1ae

157 self.session_path / sync_collection, timeline=self.timeline, chmap=chmap) 

158 return sync, chmap 2 ctx1ae

159 

160 def _extract(self, sync=None, chmap=None, sync_collection='raw_sync_data', **kwargs) -> dict: (empty)

161 trials = super()._extract(sync, chmap, sync_collection=sync_collection, **kwargs) 2 ctx1ae

162 if kwargs.get('display', False): 2 ctx1ae

163 plot_timeline(self.timeline, channels=chmap.keys(), raw=True) 

164 return trials 2 ctx1ae

165 

166 def get_bpod_event_times(self, sync, chmap, bpod_event_ttls=None, display=False, **kwargs): (empty)

167 """ 

168 Extract Bpod times from sync. 

169 

170 Unlike the superclass method. This one doesn't reassign the first trial pulse. 

171 

172 Parameters 

173 ---------- 

174 sync : dict 

175 A dictionary with keys ('times', 'polarities', 'channels'), containing the sync pulses 

176 and the corresponding channel numbers. Must contain a 'bpod' key. 

177 chmap : dict 

178 A map of channel names and their corresponding indices. 

179 bpod_event_ttls : dict of tuple 

180 A map of event names to (min, max) TTL length. 

181 

182 Returns 

183 ------- 

184 dict 

185 A dictionary with keys {'times', 'polarities'} containing Bpod TTL fronts. 

186 dict 

187 A dictionary of events (from `bpod_event_ttls`) and their intervals as an Nx2 array. 

188 """ 

189 # Assign the Bpod BNC2 events based on TTL length. The defaults are below, however these 

190 # lengths are defined by the state machine of the task protocol and therefore vary. 

191 if bpod_event_ttls is None: 1 ctx1a

192 # The trial start TTLs are often too short for the low sampling rate of the DAQ and are 

193 # therefore not used in extraction 

194 bpod_event_ttls = {'valve_open': (2.33e-4, 0.4), 'trial_end': (0.4, np.inf)} 1 ctx1a

195 bpod, bpod_event_intervals = super().get_bpod_event_times( 1 ctx1a

196 sync=sync, chmap=chmap, bpod_event_ttls=bpod_event_ttls, display=display, **kwargs) 

197 

198 # TODO Here we can make use of the 'bpod_rising_edge' channel, if available 

199 return bpod, bpod_event_intervals 1 ctx1a

200 

201 def build_trials(self, sync=None, chmap=None, **kwargs): (empty)

202 """ 

203 Extract task related event times from the sync. 

204 

205 The two major differences are that the sampling rate is lower for imaging so the short Bpod 

206 trial start TTLs are often absent. For this reason, the sync happens using the ITI_in TTL. 

207 

208 Second, the valve used at the mesoscope has a way to record the raw voltage across the 

209 solenoid, giving a more accurate readout of the valve's activity. If the reward_valve 

210 channel is present on the DAQ, this is used to extract the valve open times. 

211 

212 Parameters 

213 ---------- 

214 sync : dict 

215 'polarities' of fronts detected on sync trace for all 16 chans and their 'times' 

216 chmap : dict 

217 Map of channel names and their corresponding index. Default to constant. 

218 

219 Returns 

220 ------- 

221 dict 

222 A map of trial event timestamps. 

223 """ 

224 # Get the events from the sync. 

225 # Store the cleaned frame2ttl, audio, and bpod pulses as this will be used for QC 

226 self.frame2ttl = self.get_stimulus_update_times(sync, chmap, **kwargs) 1 ctx1a

227 self.audio, audio_event_intervals = self.get_audio_event_times(sync, chmap, **kwargs) 1 ctx1a

228 if not set(audio_event_intervals.keys()) >= {'ready_tone', 'error_tone'}: 1 ctx1a

229 raise ValueError( 

230 'Expected at least "ready_tone" and "error_tone" audio events.' 

231 '`audio_event_ttls` kwarg may be incorrect.') 

232 

233 self.bpod, bpod_event_intervals = self.get_bpod_event_times(sync, chmap, **kwargs) 1 ctx1a

234 if not set(bpod_event_intervals.keys()) >= {'valve_open', 'trial_end'}: 1 ctx1a

235 raise ValueError( 

236 'Expected at least "trial_end" and "valve_open" audio events. ' 

237 '`bpod_event_ttls` kwarg may be incorrect.') 

238 

239 t_iti_in, t_trial_end = bpod_event_intervals['trial_end'].T 1 ctx1a

240 fpga_events = alfio.AlfBunch({ 1 ctx1a

241 'itiIn_times': t_iti_in, 

242 'intervals_1': t_trial_end, 

243 'goCue_times': audio_event_intervals['ready_tone'][:, 0], 

244 'errorTone_times': audio_event_intervals['error_tone'][:, 0] 

245 }) 

246 

247 # Sync the Bpod clock to the DAQ 

248 self.bpod2fpga, drift_ppm, ibpod, ifpga = self.sync_bpod_clock(self.bpod_trials, fpga_events, self.sync_field) 1 ctx1a

249 

250 out = dict() 1 ctx1a

251 out.update({k: self.bpod_trials[k][ibpod] for k in self.bpod_fields}) 1 ctx1a

252 out.update({k: self.bpod2fpga(self.bpod_trials[k][ibpod]) for k in self.bpod_rsync_fields}) 1 ctx1a

253 

254 start_times = out['intervals'][:, 0] 1 ctx1a

255 last_trial_end = out['intervals'][-1, 1] 1 ctx1a

256 

257 def assign_to_trial(events, take='last', starts=start_times, **kwargs): 1 ctx1a

258 """Assign DAQ events to trials. 

259 

260 Because we may not have trial start TTLs on the DAQ (because of the low sampling rate), 

261 there may be an extra last trial that's not in the Bpod intervals as the extractor 

262 ignores the last trial. This function trims the input array before assigning so that 

263 the last trial's events are correctly assigned. 

264 """ 

265 return _assign_events_to_trial(starts, events[events <= last_trial_end], take, **kwargs) 1 ctx1a

266 out['itiIn_times'] = assign_to_trial(fpga_events['itiIn_times'][ifpga]) 1 ctx1a

267 

268 # Extract valve open times from the DAQ 

269 valve_driver_ttls = bpod_event_intervals['valve_open'] 1 ctx1a

270 correct = self.bpod_trials['feedbackType'] == 1 1 ctx1a

271 # If there is a reward_valve channel, the voltage across the valve has been recorded and 

272 # should give a more accurate readout of the valve's activity. 

273 if any(ch['name'] == 'reward_valve' for ch in self.timeline['meta']['inputs']): 1 ctx1a

274 # TODO Let's look at the expected open length based on calibration and reward volume 

275 # import scipy.interpolate 

276 # # FIXME support v7 settings? 

277 # fcn_vol2time = scipy.interpolate.pchip( 

278 # self.bpod_extractor.settings['device_valve']['WATER_CALIBRATION_WEIGHT_PERDROP'], 

279 # self.bpod_extractor.settings['device_valve']['WATER_CALIBRATION_OPEN_TIMES'] 

280 # ) 

281 # reward_time = fcn_vol2time(self.bpod_extractor.settings.get('REWARD_AMOUNT_UL')) / 1e3 

282 

283 # Use the driver TTLs to find the valve open times that correspond to the valve opening 

284 valve_intervals, valve_open_times = self.get_valve_open_times(driver_ttls=valve_driver_ttls) 1 ctx1a

285 if valve_open_times.size != np.sum(correct): 1 ctx1a

286 _logger.warning( 

287 'Number of valve open times does not equal number of correct trials (%i != %i)', 

288 valve_open_times.size, np.sum(correct)) 

289 

290 out['valveOpen_times'] = assign_to_trial(valve_open_times) 1 ctx1a

291 else: 

292 # Use the valve controller TTLs recorded on the Bpod channel as the reward time 

293 out['valveOpen_times'] = assign_to_trial(valve_driver_ttls[:, 0]) 

294 

295 # Stimulus times extracted based on trigger times 

296 # When assigning events all start times must not be NaN so here we substitute freeze 

297 # trigger times on nogo trials for stim on trigger times, then replace with NaN again 

298 go_trials = np.where(out['choice'] != 0)[0] 1 ctx1a

299 lims = np.copy(out['stimOnTrigger_times']) 1 ctx1a

300 lims[go_trials] = out['stimFreezeTrigger_times'][go_trials] 1 ctx1a

301 out['stimFreeze_times'] = assign_to_trial( 1 ctx1a

302 self.frame2ttl['times'], 'last', 

303 starts=lims, t_trial_end=out['stimOffTrigger_times']) 

304 out['stimFreeze_times'][out['choice'] == 0] = np.nan 1 ctx1a

305 

306 # Here we do the same but use stim off trigger times 

307 lims = np.copy(out['stimOffTrigger_times']) 1 ctx1a

308 lims[go_trials] = out['stimFreezeTrigger_times'][go_trials] 1 ctx1a

309 out['stimOn_times'] = assign_to_trial( 1 ctx1a

310 self.frame2ttl['times'], 'first', 

311 starts=out['stimOnTrigger_times'], t_trial_end=lims) 

312 out['stimOff_times'] = assign_to_trial( 1 ctx1a

313 self.frame2ttl['times'], 'first', 

314 starts=out['stimOffTrigger_times'], t_trial_end=out['intervals'][:, 1] 

315 ) 

316 

317 # Audio times 

318 error_cue = fpga_events['errorTone_times'] 1 ctx1a

319 if error_cue.size != np.sum(~correct): 1 ctx1a

320 _logger.warning( 

321 'N detected error tones does not match number of incorrect trials (%i != %i)', 

322 error_cue.size, np.sum(~correct)) 

323 go_cue = fpga_events['goCue_times'] 1 ctx1a

324 out['goCue_times'] = assign_to_trial(go_cue, take='first') 1 ctx1a

325 out['errorCue_times'] = assign_to_trial(error_cue) 1 ctx1a

326 

327 if go_cue.size > start_times.size: 1 ctx1a

328 _logger.warning( 

329 'More go cue tones detected than trials! (%i vs %i)', go_cue.size, start_times.size) 

330 elif go_cue.size < start_times.size: 1 ctx1a

331 """ 

332 If the error cues are all assigned and some go cues are missed it may be that some 

333 responses were so fast that the go cue and error tone merged, or the go cue TTL was too 

334 long. 

335 """ 

336 _logger.warning('%i go cue tones missed', start_times.size - go_cue.size) 

337 err_trig = self.bpod2fpga(self.bpod_trials['errorCueTrigger_times']) 

338 go_trig = self.bpod2fpga(self.bpod_trials['goCueTrigger_times']) 

339 assert not np.any(np.isnan(go_trig)) 

340 assert err_trig.size == go_trig.size # should be length of n trials with NaNs 

341 

342 # Find which trials are missing a go cue 

343 _go_cue = assign_to_trial(go_cue, take='first') 

344 error_cue = assign_to_trial(error_cue) 

345 missing = np.isnan(_go_cue) 

346 

347 # Get all the DAQ timestamps where audio channel was HIGH 

348 raw = timeline_get_channel(self.timeline, 'audio') 

349 raw = (raw - raw.min()) / (raw.max() - raw.min()) # min-max normalize 

350 ups = self.timeline.timestamps[raw > .5] # timestamps where input HIGH 

351 

352 # Get the timestamps of the first HIGH after the trigger times (allow up to 200ms after). 

353 # Indices of ups directly following a go trigger, or -1 if none found (or trigger NaN) 

354 idx = attribute_times(ups, go_trig, tol=0.2, take='after') 

355 # Trial indices that didn't have detected goCue and now has been assigned an `ups` index 

356 assigned = np.where(idx != -1 & missing)[0] # ignore unassigned 

357 _go_cue[assigned] = ups[idx[assigned]] 

358 

359 # Remove mis-assigned error tone times (i.e. those that have now been assigned to goCue) 

360 error_cue_without_trig, = np.where(~np.isnan(error_cue) & np.isnan(err_trig)) 

361 i_to_remove = np.intersect1d(assigned, error_cue_without_trig, assume_unique=True) 

362 error_cue[i_to_remove] = np.nan 

363 

364 # For those trials where go cue was merged with the error cue and therefore mis-assigned, 

365 # we must re-assign the error cue times as the first HIGH after the error trigger. 

366 idx = attribute_times(ups, err_trig, tol=0.2, take='after') 

367 assigned = np.where(idx != -1 & missing)[0] # ignore unassigned 

368 error_cue[assigned] = ups[idx[assigned]] 

369 out['goCue_times'] = _go_cue 

370 out['errorCue_times'] = error_cue 

371 

372 # Because we're not 

373 assert np.intersect1d(out['goCue_times'], out['errorCue_times']).size == 0, \ 1 ctx1a

374 'audio tones not assigned correctly; tones likely missed' 

375 

376 # Feedback times 

377 out['feedback_times'] = np.copy(out['valveOpen_times']) 1 ctx1a

378 ind_err = np.isnan(out['valveOpen_times']) 1 ctx1a

379 out['feedback_times'][ind_err] = out['errorCue_times'][ind_err] 1 ctx1a

380 

381 return out 1 ctx1a

382 

383 def extract_wheel_sync(self, ticks=WHEEL_TICKS, radius=WHEEL_RADIUS_CM, coding='x4', tmin=None, tmax=None): (empty)

384 """ 

385 Gets the wheel position from Timeline counter channel. 

386 

387 Parameters 

388 ---------- 

389 ticks : int 

390 Number of ticks corresponding to a full revolution (1024 for IBL rotary encoder). 

391 radius : float 

392 Radius of the wheel. Defaults to 1 for an output in radians. 

393 coding : str {'x1', 'x2', 'x4'} 

394 Rotary encoder encoding (IBL default is x4). 

395 tmin : float 

396 The minimum time from which to extract the sync pulses. 

397 tmax : float 

398 The maximum time up to which we extract the sync pulses. 

399 

400 Returns 

401 ------- 

402 np.array 

403 Wheel timestamps in seconds. 

404 np.array 

405 Wheel positions in radians. 

406 

407 See Also 

408 -------- 

409 ibllib.io.extractors.ephys_fpga.extract_wheel_sync 

410 """ 

411 if coding not in ('x1', 'x2', 'x4'): 2 ctx1ia

412 raise ValueError('Unsupported coding; must be one of x1, x2 or x4') 1 ctx1i

413 raw = correct_counter_discontinuities(timeline_get_channel(self.timeline, 'rotary_encoder')) 2 ctx1ia

414 

415 # Timeline evenly samples counter so we extract only change points 

416 d = np.diff(raw) 2 ctx1ia

417 ind, = np.where(~np.isclose(d, 0)) 2 ctx1ia

418 pos = raw[ind + 1] 2 ctx1ia

419 pos -= pos[0] # Start from zero 2 ctx1ia

420 pos = pos / ticks * np.pi * 2 * radius / int(coding[1]) # Convert to radians 2 ctx1ia

421 

422 # Get timestamps of changes and trim based on protocol spacers 

423 ts = self.timeline['timestamps'][ind + 1] 2 ctx1ia

424 tmin = ts.min() if tmin is None else tmin 2 ctx1ia

425 tmax = ts.max() if tmax is None else tmax 2 ctx1ia

426 mask = np.logical_and(ts >= tmin, ts <= tmax) 2 ctx1ia

427 return ts[mask], pos[mask] 2 ctx1ia

428 

429 def get_wheel_positions(self, ticks=WHEEL_TICKS, radius=WHEEL_RADIUS_CM, coding='x4', (empty)

430 tmin=None, tmax=None, display=False, **kwargs): 

431 """ 

432 Gets the wheel position and detected movements from Timeline counter channel. 

433 

434 Called by the super class extractor (FPGATrials._extract). 

435 

436 Parameters 

437 ---------- 

438 ticks : int 

439 Number of ticks corresponding to a full revolution (1024 for IBL rotary encoder). 

440 radius : float 

441 Radius of the wheel. Defaults to 1 for an output in radians. 

442 coding : str {'x1', 'x2', 'x4'} 

443 Rotary encoder encoding (IBL default is x4). 

444 tmin : float 

445 The minimum time from which to extract the sync pulses. 

446 tmax : float 

447 The maximum time up to which we extract the sync pulses. 

448 display : bool 

449 If true, plot the wheel positions from bpod and the DAQ. 

450 

451 Returns 

452 ------- 

453 dict 

454 wheel object with keys ('timestamps', 'position'). 

455 dict 

456 wheelMoves object with keys ('intervals' 'peakAmplitude'). 

457 """ 

458 wheel = self.extract_wheel_sync(ticks=ticks, radius=radius, coding=coding, tmin=tmin, tmax=tmax) 2 ctx1ia

459 wheel = dict(zip(('timestamps', 'position'), wheel)) 2 ctx1ia

460 moves = extract_wheel_moves(wheel['timestamps'], wheel['position']) 2 ctx1ia

461 

462 if display: 2 ctx1ia

463 assert self.bpod_trials, 'no bpod trials to compare' 1 ctx1i

464 fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True) 1 ctx1i

465 bpod_ts = self.bpod_trials['wheel_timestamps'] 1 ctx1i

466 bpod_pos = self.bpod_trials['wheel_position'] 1 ctx1i

467 ax0.plot(self.bpod2fpga(bpod_ts), bpod_pos) 1 ctx1i

468 ax0.set_ylabel('Bpod wheel position / rad') 1 ctx1i

469 ax1.plot(wheel['timestamps'], wheel['position']) 1 ctx1i

470 ax1.set_ylabel('DAQ wheel position / rad'), ax1.set_xlabel('Time / s') 1 ctx1i

471 return wheel, moves 2 ctx1ia

472 

473 def get_valve_open_times(self, display=False, threshold=100, driver_ttls=None): (empty)

474 """ 

475 Get the valve open times from the raw timeline voltage trace. 

476 

477 Parameters 

478 ---------- 

479 display : bool 

480 Plot detected times on the raw voltage trace. 

481 threshold : float 

482 The threshold of voltage change to apply. The default was set by eye; units should be 

483 Volts per sample but doesn't appear to be. 

484 driver_ttls : numpy.array 

485 An optional array of driver TTLs to use for assigning with the valve times. 

486 

487 Returns 

488 ------- 

489 numpy.array 

490 The detected valve open intervals. 

491 numpy.array 

492 If driver_ttls is not None, returns an array of open times that occurred directly after 

493 the driver TTLs. 

494 """ 

495 WARN_THRESH = 10e-3 # open time threshold below which to log warning 3 ctx1gae

496 tl = self.timeline 3 ctx1gae

497 info = next(x for x in tl['meta']['inputs'] if x['name'] == 'reward_valve') 3 ctx1gae

498 values = tl['raw'][:, info['arrayColumn'] - 1] # Timeline indices start from 1 3 ctx1gae

499 

500 # The voltage changes over ~1ms and can therefore occur over two DAQ samples at 2kHz 

501 # making simple thresholding an issue. For this reason we convolve the signal with a 

502 # window and detect the peaks and troughs. 

503 if (Fs := tl['meta']['daqSampleRate']) != 2000: # e.g. 2kHz 3 ctx1gae

504 _logger.warning('Reward valve detection not tested with a DAQ sample rate of %i', Fs) 3 ctx1gae

505 dt = 1e-3 # change in voltage takes ~1ms when changing valve open state 3 ctx1gae

506 N = dt / (1 / Fs) # this means voltage change occurs over N samples 3 ctx1gae

507 vel, _ = velocity_filtered(values, int(Fs / N)) # filtered voltage change over time 3 ctx1gae

508 ups, _ = find_peaks(vel, height=threshold) # valve closes (-5V -> 0V) 3 ctx1gae

509 downs, _ = find_peaks(-1 * vel, height=threshold) # valve opens (0V -> -5V) 3 ctx1gae

510 

511 # Convert these times into intervals 

512 ixs = np.argsort(np.r_[downs, ups]) # sort indices 3 ctx1gae

513 times = tl['timestamps'][np.r_[downs, ups]][ixs] # ordered valve event times 3 ctx1gae

514 polarities = np.r_[np.zeros_like(downs) - 1, np.ones_like(ups)][ixs] # polarity sorted 3 ctx1gae

515 missing = np.where(np.diff(polarities) == 0)[0] # if some changes were missed insert NaN 3 ctx1gae

516 times = np.insert(times, missing + int(polarities[0] == -1), np.nan) 3 ctx1gae

517 if polarities[-1] == -1: # ensure ends with a valve close 3 ctx1gae

518 times = np.r_[times, np.nan] 2 ctx1ga

519 if polarities[0] == 1: # ensure starts with a valve open 3 ctx1gae

520 # It seems it can start out at -5V (open), then when the reward happens it closes and 

521 # immediately opens. In this case we insert discard the first open time. 

522 times = np.r_[np.nan, times] 

523 intervals = times.reshape(-1, 2) 3 ctx1gae

524 

525 # Log warning of improbably short intervals 

526 short = np.sum(np.diff(intervals) < WARN_THRESH) 3 ctx1gae

527 if short > 0: 3 ctx1gae

528 _logger.warning('%i valve open intervals shorter than %i ms', short, WARN_THRESH) 

529 

530 # The closing of the valve is noisy. Keep only the falls that occur immediately after a Bpod TTL 

531 if driver_ttls is not None: 3 ctx1gae

532 # Returns an array of open_times indices, one for each driver TTL 

533 ind = attribute_times(intervals[:, 0], driver_ttls[:, 0], tol=.1, take='after') 2 ctx1ae

534 open_times = intervals[ind[ind >= 0], 0] 2 ctx1ae

535 # TODO Log any > 40ms? Difficult to report missing valve times because of calibration 

536 

537 if display: 3 ctx1gae

538 fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True) 1 ctx1g

539 ax0.plot(tl['timestamps'], timeline_get_channel(tl, 'bpod'), color='grey', linestyle='-') 1 ctx1g

540 if driver_ttls is not None: 1 ctx1g

541 x = np.empty_like(driver_ttls.flatten()) 

542 x[0::2] = driver_ttls[:, 0] 

543 x[1::2] = driver_ttls[:, 1] 

544 y = np.ones_like(x) 

545 y[1::2] -= 2 

546 squares(x, y, ax=ax0, yrange=[0, 5]) 

547 # vertical_lines(driver_ttls, ymax=5, ax=ax0, linestyle='--', color='b') 

548 ax0.plot(open_times, np.ones_like(open_times) * 4.5, 'g*') 

549 ax1.plot(tl['timestamps'], values, 'k-o') 1 ctx1g

550 ax1.set_ylabel('Voltage / V'), ax1.set_xlabel('Time / s') 1 ctx1g

551 

552 ax2 = ax1.twinx() 1 ctx1g

553 ax2.set_ylabel('dV', color='grey') 1 ctx1g

554 ax2.plot(tl['timestamps'], vel, linestyle='-', color='grey') 1 ctx1g

555 ax2.plot(intervals[:, 1], np.ones(len(intervals)) * threshold, 'r*', label='close') 1 ctx1g

556 ax2.plot(intervals[:, 0], np.ones(len(intervals)) * threshold, 'g*', label='open') 1 ctx1g

557 return intervals if driver_ttls is None else (intervals, open_times) 3 ctx1gae

558 

559 def _assign_events_audio(self, audio_times, audio_polarities, display=False): (empty)

560 """ 

561 This is identical to ephys_fpga._assign_events_audio, except for the ready tone threshold. 

562 

563 Parameters 

564 ---------- 

565 audio_times : numpy.array 

566 An array of audio TTL front times. 

567 audio_polarities : numpy.array 

568 An array of audio TTL front polarities (1 for rises, -1 for falls). 

569 display : bool 

570 If true, display audio pulses and the assigned onsets. 

571 

572 Returns 

573 ------- 

574 numpy.array 

575 The times of the go cue onsets. 

576 numpy.array 

577 The times of the error tone onsets. 

578 """ 

579 # make sure that there are no 2 consecutive fall or consecutive rise events 

580 assert np.all(np.abs(np.diff(audio_polarities)) == 2) 

581 # take only even time differences: i.e. from rising to falling fronts 

582 dt = np.diff(audio_times) 

583 onsets = audio_polarities[:-1] == 1 

584 

585 # error tones are events lasting from 400ms to 1200ms 

586 i_error_tone_in = np.where(np.logical_and(0.4 < dt, dt < 1.2) & onsets)[0] 

587 t_error_tone_in = audio_times[i_error_tone_in] 

588 

589 # detect ready tone by length below 300 ms 

590 i_ready_tone_in = np.where(np.logical_and(dt <= 0.3, onsets))[0] 

591 t_ready_tone_in = audio_times[i_ready_tone_in] 

592 if display: # pragma: no cover 

593 fig, ax = plt.subplots(nrows=2, sharex=True) 

594 ax[0].plot(self.timeline.timestamps, timeline_get_channel(self.timeline, 'audio'), 'k-o') 

595 ax[0].set_ylabel('Voltage / V') 

596 squares(audio_times, audio_polarities, yrange=[-1, 1], ax=ax[1]) 

597 vertical_lines(t_ready_tone_in, ymin=-.8, ymax=.8, ax=ax[1], label='go cue') 

598 vertical_lines(t_error_tone_in, ymin=-.8, ymax=.8, ax=ax[1], label='error tone') 

599 ax[1].set_xlabel('Time / s') 

600 ax[1].legend() 

601 

602 return t_ready_tone_in, t_error_tone_in 

603 

604 

605class TimelineTrialsHabituation(FpgaTrialsHabituation, TimelineTrials): (empty)

606 """Habituation trials extraction from Timeline DAQ data.""" 

607 

608 sync_field = 'intervals_0' (empty)

609 

610 def build_trials(self, sync=None, chmap=None, **kwargs): (empty)

611 """ 

612 Extract task related event times from the sync. 

613 

614 The valve used at the mesoscope has a way to record the raw voltage across the solenoid, 

615 giving a more accurate readout of the valve's activity. If the reward_valve channel is 

616 present on the DAQ, this is used to extract the valve open times. 

617 

618 Parameters 

619 ---------- 

620 sync : dict 

621 'polarities' of fronts detected on sync trace for all 16 chans and their 'times' 

622 chmap : dict 

623 Map of channel names and their corresponding index. Default to constant. 

624 

625 Returns 

626 ------- 

627 dict 

628 A map of trial event timestamps. 

629 """ 

630 out = super().build_trials(sync, chmap, **kwargs) 1 ctx1e

631 

632 start_times = out['intervals'][:, 0] 1 ctx1e

633 last_trial_end = out['intervals'][-1, 1] 1 ctx1e

634 

635 # Extract valve open times from the DAQ 

636 _, bpod_event_intervals = self.get_bpod_event_times(sync, chmap, **kwargs) 1 ctx1e

637 bpod_feedback_times = self.bpod2fpga(self.bpod_trials['feedback_times']) 1 ctx1e

638 valve_driver_ttls = bpod_event_intervals['valve_open'] 1 ctx1e

639 # If there is a reward_valve channel, the voltage across the valve has been recorded and 

640 # should give a more accurate readout of the valve's activity. 

641 if any(ch['name'] == 'reward_valve' for ch in self.timeline['meta']['inputs']): 1 ctx1e

642 # Use the driver TTLs to find the valve open times that correspond to the valve opening 

643 valve_intervals, valve_open_times = self.get_valve_open_times(driver_ttls=valve_driver_ttls) 1 ctx1e

644 if valve_open_times.size != start_times.size: 1 ctx1e

645 _logger.warning( 1 ctx1e

646 'Number of valve open times does not equal number of correct trials (%i != %i)', 

647 valve_open_times.size, start_times.size) 

648 else: 

649 # Use the valve controller TTLs recorded on the Bpod channel as the reward time 

650 valve_open_times = valve_driver_ttls[:, 0] 

651 # there may be an extra last trial that's not in the Bpod intervals as the extractor ignores the last trial 

652 valve_open_times = valve_open_times[valve_open_times <= last_trial_end] 1 ctx1e

653 out['valveOpen_times'] = _assign_events_to_trial( 1 ctx1e

654 bpod_feedback_times, valve_open_times, take='first', t_trial_end=out['intervals'][:, 1]) 

655 

656 # Feedback times 

657 out['feedback_times'] = np.copy(out['valveOpen_times']) 1 ctx1e

658 

659 return out 1 ctx1e

660 

661 

662class MesoscopeSyncTimeline(extractors_base.BaseExtractor): (empty)

663 """Extraction of mesoscope imaging times.""" 

664 

665 var_names = ('mpci_times', 'mpciStack_timeshift') (empty)

666 save_names = ('mpci.times.npy', 'mpciStack.timeshift.npy') (empty)

667 

668 """one.alf.io.AlfBunch: The raw imaging meta data and frame times""" (empty)

669 rawImagingData = None (empty)

670 

671 def __init__(self, session_path, n_FOVs): (empty)

672 """ 

673 Extract the mesoscope frame times from DAQ data acquired through Timeline. 

674 

675 Parameters 

676 ---------- 

677 session_path : str, pathlib.Path 

678 The session path to extract times from. 

679 n_FOVs : int 

680 The number of fields of view acquired.