Coverage for ibllib/io/extractors/ephys_passive.py: 71%

1#!/usr/bin/env python 

2# -*- coding:utf-8 -*- 

3# @Author: Niccolò Bonacchi 

4# @Date: Monday, September 7th 2020, 11:51:17 am 

5import json (empty)

6import logging (empty)

7from pathlib import Path (empty)

8from typing import Tuple (empty)

9 

10import matplotlib.pyplot as plt (empty)

11import numpy as np (empty)

12import pandas as pd (empty)

13 

14import ibllib.io.raw_data_loaders as rawio (empty)

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

16from ibllib.io.extractors.base import BaseExtractor (empty)

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

18 plot_audio_times, 

19 plot_gabor_times, 

20 plot_passive_periods, 

21 plot_rfmapping, 

22 plot_stims_times, 

23 plot_sync_channels, 

24 plot_valve_times, 

25) 

26 

27log = logging.getLogger("ibllib") (empty)

28 

29# hardcoded var 

30FRAME_FS = 60 # Sampling freq of the ipad screen, in Hertz (empty)

31FS_FPGA = 30000 # Sampling freq of the neural recording system screen, in Hertz (empty)

32NVALVE = 40 # number of expected valve clicks (empty)

33NGABOR = 20 + 20 * 4 * 2 # number of expected Gabor patches (empty)

34NTONES = 40 (empty)

35NNOISES = 40 (empty)

36DEBUG_PLOTS = False (empty)

37 

38dataset_types = [ (empty)

39 "_spikeglx_sync.times", 

40 "_spikeglx_sync.channels", 

41 "_spikeglx_sync.polarities", 

42 "_iblrig_RFMapStim.raw", 

43 "_iblrig_stimPositionScreen.raw", 

44 "_iblrig_syncSquareUpdate.raw", 

45 "ephysData.raw.meta", 

46 "_iblrig_taskSettings.raw", 

47 "_iblrig_taskData.raw", 

48] 

49 

50min_dataset_types = [ (empty)

51 "_spikeglx_sync.times", 

52 "_spikeglx_sync.channels", 

53 "_spikeglx_sync.polarities", 

54 "_iblrig_RFMapStim.raw", 

55 "ephysData.raw.meta", 

56 "_iblrig_taskSettings.raw", 

57 "_iblrig_taskData.raw", 

58] 

59 

60 

61# load session fixtures 

62def _load_passive_session_fixtures(session_path: str, task_collection: str = 'raw_passive_data') -> dict: (empty)

63 """load_passive_session_fixtures Loads corresponding ephys session fixtures 

64 

65 :param session_path: the path to a session 

66 :type session_path: str 

67 :return: position contrast phase delays and stim id's 

68 :rtype: dict 

69 """ 

70 

71 # THIS CAN BE PREGENERATED SESSION NO 

72 settings = rawio.load_settings(session_path, task_collection=task_collection) 4 ctx1abcd

73 ses_nb = settings['PREGENERATED_SESSION_NUM'] 4 ctx1abcd

74 session_order = settings.get('SESSION_ORDER', None) 4 ctx1abcd

75 if session_order: # TODO test this out and make sure it okay 4 ctx1abcd

76 assert settings["SESSION_ORDER"][settings["SESSION_IDX"]] == ses_nb 3 ctx1acd

77 

78 path_fixtures = Path(ephys_fpga.__file__).parent.joinpath("ephys_sessions") 4 ctx1abcd

79 

80 fixture = { 4 ctx1abcd

81 "pcs": np.load(path_fixtures.joinpath(f"session_{ses_nb}_passive_pcs.npy")), 

82 "delays": np.load(path_fixtures.joinpath(f"session_{ses_nb}_passive_stimDelays.npy")), 

83 "ids": np.load(path_fixtures.joinpath(f"session_{ses_nb}_passive_stimIDs.npy")), 

84 } 

85 

86 return fixture 4 ctx1abcd

87 

88 

89def _load_task_version(session_path: str, task_collection: str = 'raw_passive_data') -> str: (empty)

90 """Find the IBL rig version used for the session 

91 

92 :param session_path: the path to a session 

93 :type session_path: str 

94 :return: ibl rig task protocol version 

95 :rtype: str 

96 

97 """ 

98 settings = rawio.load_settings(session_path, task_collection=task_collection) 4 ctx1abcd

99 ses_ver = settings["IBLRIG_VERSION"] 4 ctx1abcd

100 

101 return ses_ver 4 ctx1abcd

102 

103 

104def skip_task_replay(session_path: str, task_collection: str = 'raw_passive_data') -> bool: (empty)

105 """Find whether the task replay portion of the passive stimulus has been shown 

106 

107 :param session_path: the path to a session 

108 :type session_path: str 

109 :param task_collection: collection containing task data 

110 :type task_collection: str 

111 :return: whether or not the task replay has been run 

112 :rtype: bool 

113 """ 

114 

115 settings = rawio.load_settings(session_path, task_collection=task_collection) 4 ctx1abcd

116 # Attempt to see if SKIP_EVENT_REPLAY is available, if not assume we do have task replay 

117 skip_replay = settings.get('SKIP_EVENT_REPLAY', False) 4 ctx1abcd

118 

119 return skip_replay 4 ctx1abcd

120 

121 

122def _load_passive_stim_meta() -> dict: (empty)

123 """load_passive_stim_meta Loads the passive protocol metadata 

124 

125 :return: metadata about passive protocol stimulus presentation 

126 :rtype: dict 

127 """ 

128 path_fixtures = Path(ephys_fpga.__file__).parent.joinpath("ephys_sessions") 6 ctx1haebcd

129 with open(path_fixtures.joinpath("passive_stim_meta.json"), "r") as f: 6 ctx1haebcd

130 meta = json.load(f) 6 ctx1haebcd

131 

132 return meta 6 ctx1haebcd

133 

134 

135# 1/3 Define start and end times of the 3 passive periods 

136def _get_spacer_times(spacer_template, jitter, ttl_signal, t_quiet, thresh=3.0): (empty)

137 """ 

138 Find timestamps of spacer signal. 

139 :param spacer_template: list of indices where ttl signal changes 

140 :type spacer_template: array-like 

141 :param jitter: jitter (in seconds) for matching ttl_signal with spacer_template 

142 :type jitter: float 

143 :param ttl_signal: 

144 :type ttl_signal: array-like 

145 :param t_quiet: seconds between spacer and next stim 

146 :type t_quiet: float 

147 :param thresh: threshold value for the fttl convolved signal (with spacer template) to pass over to detect a spacer 

148 :type thresh: float 

149 :return: times of spacer onset/offset 

150 :rtype: n_spacer x 2 np.ndarray; first col onset times, second col offset 

151 """ 

152 diff_spacer_template = np.diff(spacer_template) 5 ctx1aebcd

153 # add jitter; 

154 # remove extreme values 

155 spacer_model = jitter + diff_spacer_template[2:-2] 5 ctx1aebcd

156 # diff ttl signal to compare to spacer_model 

157 dttl = np.diff(ttl_signal) 5 ctx1aebcd

158 # remove diffs larger than max diff in model to clean up signal 

159 dttl[dttl > np.max(spacer_model)] = 0 5 ctx1aebcd

160 # convolve cleaned diff ttl signal w/ spacer model 

161 conv_dttl = np.correlate(dttl, spacer_model, mode="full") 5 ctx1aebcd

162 # find spacer location 

163 idxs_spacer_middle = np.where( 5 ctx1aebcd

164 (conv_dttl[1:-2] < thresh) & (conv_dttl[2:-1] > thresh) & (conv_dttl[3:] < thresh) 

165 )[0] 

166 # adjust indices for 

167 # - `np.where` call above 

168 # - length of spacer_model 

169 spacer_around = int((np.floor(len(spacer_model) / 2))) 5 ctx1aebcd

170 idxs_spacer_middle += 2 - spacer_around 5 ctx1aebcd

171 

172 # for each spacer make sure the times are monotonically non-decreasing before 

173 # and monotonically non-increasing afterwards 

174 is_valid = np.zeros((idxs_spacer_middle.size), dtype=bool) 5 ctx1aebcd

175 for i, t in enumerate(idxs_spacer_middle): 5 ctx1aebcd

176 before = all(np.diff(dttl[t - spacer_around:t]) >= 0) 5 ctx1aebcd

177 after = all(np.diff(dttl[t + 1:t + 1 + spacer_around]) <= 0) 5 ctx1aebcd

178 is_valid[i] = np.bitwise_and(before, after) 5 ctx1aebcd

179 

180 idxs_spacer_middle = idxs_spacer_middle[is_valid] 5 ctx1aebcd

181 

182 # pull out spacer times (middle) 

183 ts_spacer_middle = ttl_signal[idxs_spacer_middle] 5 ctx1aebcd

184 # put beginning/end of spacer times into an array 

185 spacer_length = np.max(spacer_template) 5 ctx1aebcd

186 spacer_times = np.zeros(shape=(ts_spacer_middle.shape[0], 2)) 5 ctx1aebcd

187 for i, t in enumerate(ts_spacer_middle): 5 ctx1aebcd

188 spacer_times[i, 0] = t - (spacer_length / 2) - t_quiet 5 ctx1aebcd

189 spacer_times[i, 1] = t + (spacer_length / 2) + t_quiet 5 ctx1aebcd

190 return spacer_times, conv_dttl 5 ctx1aebcd

191 

192 

193def _get_passive_spacers(session_path, sync_collection='raw_ephys_data', (empty)

194 sync=None, sync_map=None, tmin=None, tmax=None): 

195 """ 

196 load and get spacer information, do corr to find spacer timestamps 

197 returns t_passive_starts, t_starts, t_ends 

198 """ 

199 if sync is None or sync_map is None: 5 ctx1aebcd

200 sync, sync_map = ephys_fpga.get_sync_and_chn_map(session_path, sync_collection=sync_collection) 

201 meta = _load_passive_stim_meta() 5 ctx1aebcd

202 # t_end_ephys = passive.ephysCW_end(session_path=session_path) 

203 fttl = ephys_fpga.get_sync_fronts(sync, sync_map["frame2ttl"], tmin=tmin, tmax=tmax) 5 ctx1aebcd

204 fttl = ephys_fpga._clean_frame2ttl(fttl, display=False) 5 ctx1aebcd

205 spacer_template = ( 5 ctx1aebcd

206 np.array(meta["VISUAL_STIM_0"]["ttl_frame_nums"], dtype=np.float32) / FRAME_FS 

207 ) 

208 jitter = 3 / FRAME_FS # allow for 3 screen refresh as jitter 5 ctx1aebcd

209 t_quiet = meta["VISUAL_STIM_0"]["delay_around"] 5 ctx1aebcd

210 spacer_times, conv_dttl = _get_spacer_times( 5 ctx1aebcd

211 spacer_template=spacer_template, jitter=jitter, ttl_signal=fttl["times"], t_quiet=t_quiet 

212 ) 

213 

214 # Check correct number of spacers found 

215 n_exp_spacer = np.sum(np.array(meta['STIM_ORDER']) == 0) # Hardcoded 0 for spacer 5 ctx1aebcd

216 if n_exp_spacer != np.size(spacer_times) / 2: 5 ctx1aebcd

217 error_nspacer = True 

218 # sometimes the first spacer is truncated 

219 # assess whether the first spacer is undetected, and then launch another spacer detection on truncated fttl 

220 # with a lower threshold value 

221 # Note: take *3 for some margin 

222 if spacer_times[0][0] > (spacer_template[-1] + jitter) * 3 and (np.size(spacer_times) / 2) == n_exp_spacer - 1: 

223 # Truncate signals 

224 fttl_t = fttl["times"][np.where(fttl["times"] < spacer_times[0][0])] 

225 conv_dttl_t = conv_dttl[np.where(fttl["times"] < spacer_times[0][0])] 

226 ddttl = np.diff(np.diff(fttl_t)) 

227 # Find spacer location 

228 # NB: cannot re-use the same algo for spacer detection as conv peaks towards spacer end 

229 # 1. Find time point at which conv raises above a given threshold value 

230 thresh = 2.0 

231 idx_nearend_spacer = int(np.where((conv_dttl_t[1:-2] < thresh) & (conv_dttl_t[2:-1] > thresh))[0]) 

232 ddttl = ddttl[0:idx_nearend_spacer] 

233 # 2. Find time point before this, for which fttl diff increase/decrease (this is the middle of spacer) 

234 indx_middle = np.where((ddttl[0:-1] > 0) & (ddttl[1:] < 0))[0] 

235 if len(indx_middle) == 1: 

236 # 3. Add 1/2 spacer to middle idx to get the spacer end indx 

237 spacer_around = int((np.floor(len(spacer_template) / 2))) 

238 idx_end = int(indx_middle + spacer_around) + 1 

239 spacer_times = np.insert(spacer_times, 0, np.array([fttl["times"][0], fttl["times"][idx_end]]), axis=0) 

240 error_nspacer = False 

241 

242 if error_nspacer: 

243 raise ValueError( 

244 f'The number of expected spacer ({n_exp_spacer}) ' 

245 f'is different than the one found on the raw ' 

246 f'trace ({int(np.size(spacer_times) / 2)})' 

247 ) 

248 

249 if tmax is None: 5 ctx1aebcd

250 tmax = sync['times'][-1] 5 ctx1aebcd

251 

252 spacer_times = np.r_[spacer_times.flatten(), tmax] 5 ctx1aebcd

253 return spacer_times[0], spacer_times[1::2], spacer_times[2::2] 5 ctx1aebcd

254 

255 

256# 2/3 RFMapping stimuli 

257def _interpolate_rf_mapping_stimulus(idxs_up, idxs_dn, times, Xq, t_bin): (empty)

258 """ 

259 Interpolate stimulus presentation times to screen refresh rate to match `frames` 

260 :param ttl_01: 

261 :type ttl_01: array-like 

262 :param times: array of stimulus switch times 

263 :type times: array-like 

264 :param Xq: number of times (found in frames) 

265 :type frames: array-like 

266 :param t_bin: period of screen refresh rate 

267 :type t_bin: float 

268 :return: tuple of (stim_times, stim_frames) 

269 """ 

270 

271 beg_extrap_val = -10001 6 ctx1gaebcd

272 end_extrap_val = -10000 6 ctx1gaebcd

273 

274 X = np.sort(np.concatenate([idxs_up, idxs_dn])) 6 ctx1gaebcd

275 # make left and right extrapolations distinctive to easily find later 

276 Tq = np.interp(Xq, X, times, left=beg_extrap_val, right=end_extrap_val) 6 ctx1gaebcd

277 # uniform spacing outside boundaries of ttl signal 

278 # first values 

279 n_beg = len(np.where(Tq == beg_extrap_val)[0]) 6 ctx1gaebcd

280 if 0 < n_beg < Tq.shape[0]: 6 ctx1gaebcd

281 Tq[:n_beg] = times[0] - np.arange(n_beg, 0, -1) * t_bin 5 ctx1aebcd

282 # end values 

283 n_end = len(np.where(Tq == end_extrap_val)[0]) 6 ctx1gaebcd

284 if 0 < n_end < Tq.shape[0]: 6 ctx1gaebcd

285 Tq[-n_end:] = times[-1] + np.arange(1, n_end + 1) * t_bin 6 ctx1gaebcd

286 return Tq 6 ctx1gaebcd

287 

288 

289def _get_id_raisefall_from_analogttl(ttl_01): (empty)

290 """ 

291 Get index of raise/fall from analog continuous TTL signal (0-1 values) 

292 :param ttl_01: analog continuous TTL signal (0-1 values) 

293 :return: index up (0>1), index down (1>0), number of ttl transition 

294 """ 

295 # Check values are 0, 1, -1 

296 if not np.all(np.isin(np.unique(ttl_01), [-1, 0, 1])): 5 ctx1aebcd

297 raise ValueError("Values in input must be 0, 1, -1") 

298 else: 

299 # Find number of passage from [0 1] and [0 -1] 

300 d_ttl_01 = np.diff(ttl_01) 5 ctx1aebcd

301 id_up = np.where(np.logical_and(ttl_01 == 0, np.append(d_ttl_01, 0) == 1))[0] 5 ctx1aebcd

302 id_dw = np.where(np.logical_and(ttl_01 == 0, np.append(d_ttl_01, 0) == -1))[0] 5 ctx1aebcd

303 n_ttl_expected = 2 * (len(id_up) + len(id_dw)) # *2 for rise/fall of ttl pulse 5 ctx1aebcd

304 return id_up, id_dw, n_ttl_expected 5 ctx1aebcd

305 

306 

307def _reshape_RF(RF_file, meta_stim): (empty)

308 """ 

309 Reshape Receptive Field (RF) matrix. Take data associated to corner 

310 where frame2ttl placed to create TTL trace. 

311 :param RF_file: vector to be reshaped, containing RF info 

312 :param meta_stim: variable containing metadata information on RF 

313 :return: frames (reshaped RF), analog trace (0-1 values) 

314 """ 

315 frame_array = np.fromfile(RF_file, dtype="uint8") 5 ctx1aebcd

316 y_pix, x_pix, _ = meta_stim["stim_file_shape"] 5 ctx1aebcd

317 frames = np.transpose(np.reshape(frame_array, [y_pix, x_pix, -1], order="F"), [2, 1, 0]) 5 ctx1aebcd

318 ttl_trace = frames[:, 0, 0] 5 ctx1aebcd

319 # Convert values to 0,1,-1 for simplicity 

320 ttl_analogtrace_01 = np.zeros(np.size(ttl_trace)) 5 ctx1aebcd

321 ttl_analogtrace_01[np.where(ttl_trace == 0)] = -1 5 ctx1aebcd

322 ttl_analogtrace_01[np.where(ttl_trace == 255)] = 1 5 ctx1aebcd

323 return frames, ttl_analogtrace_01 5 ctx1aebcd

324 

325 

326# 3/3 Replay of task stimuli 

327def _extract_passiveGabor_df(fttl: dict, session_path: str, task_collection: str = 'raw_passive_data') -> pd.DataFrame: (empty)

328 # At this stage we want to define what pulses are and not quality control them. 

329 # Pulses are strictly alternating with intervals 

330 # find min max lengths for both (we don't know which are pulses and which are intervals yet) 

331 # trim edges of pulses 

332 diff0 = (np.min(np.diff(fttl["times"])[2:-2:2]), np.max(np.diff(fttl["times"])[2:-1:2])) 4 ctx1abcd

333 diff1 = (np.min(np.diff(fttl["times"])[3:-2:2]), np.max(np.diff(fttl["times"])[3:-1:2])) 4 ctx1abcd

334 # Highest max is of the intervals 

335 if max(diff0 + diff1) in diff0: 4 ctx1abcd

336 thresh = diff0[0] 

337 elif max(diff0 + diff1) in diff1: 4 ctx1abcd

338 thresh = diff1[0] 4 ctx1abcd

339 # Anything lower than the min length of intervals is a pulse 

340 idx_start_stims = np.where((np.diff(fttl["times"]) < thresh) & (np.diff(fttl["times"]) > 0.1))[0] 4 ctx1abcd

341 # Check if any pulse has been missed 

342 # i.e. expected length (without first pulse) and that it's alternating 

343 if len(idx_start_stims) < NGABOR - 1 and np.any(np.diff(idx_start_stims) > 2): 4 ctx1abcd

344 log.warning("Looks like one or more pulses were not detected, trying to extrapolate...") 

345 missing_where = np.where(np.diff(idx_start_stims) > 2)[0] 

346 insert_where = missing_where + 1 

347 missing_value = idx_start_stims[missing_where] + 2 

348 idx_start_stims = np.insert(idx_start_stims, insert_where, missing_value) 

349 

350 idx_end_stims = idx_start_stims + 1 4 ctx1abcd

351 

352 start_times = fttl["times"][idx_start_stims] 4 ctx1abcd

353 end_times = fttl["times"][idx_end_stims] 4 ctx1abcd

354 # Check if we missed the first stim 

355 if len(start_times) < NGABOR: 4 ctx1abcd

356 first_stim_off_idx = idx_start_stims[0] - 1 4 ctx1abcd

357 # first_stim_on_idx = first_stim_off_idx - 1 

358 end_times = np.insert(end_times, 0, fttl["times"][first_stim_off_idx]) 4 ctx1abcd

359 start_times = np.insert(start_times, 0, end_times[0] - 0.3) 4 ctx1abcd

360 

361 # intervals dstype requires reshaping of start and end times 

362 passiveGabor_intervals = np.array([(x, y) for x, y in zip(start_times, end_times)]) 4 ctx1abcd

363 

364 # Check length of presentation of stim is within 150ms of expected 

365 if not np.allclose([y - x for x, y in passiveGabor_intervals], 0.3, atol=0.15): 4 ctx1abcd

366 log.warning("Some Gabor presentation lengths seem wrong.") 

367 

368 assert ( 

369 len(passiveGabor_intervals) == NGABOR 

370 ), f"Wrong number of Gabor stimuli detected: {len(passiveGabor_intervals)} / {NGABOR}" 

371 fixture = _load_passive_session_fixtures(session_path, task_collection) 4 ctx1abcd

372 passiveGabor_properties = fixture["pcs"] 4 ctx1abcd

373 passiveGabor_table = np.append(passiveGabor_intervals, passiveGabor_properties, axis=1) 4 ctx1abcd

374 columns = ["start", "stop", "position", "contrast", "phase"] 4 ctx1abcd

375 passiveGabor_df = pd.DataFrame(passiveGabor_table, columns=columns) 4 ctx1abcd

376 return passiveGabor_df 4 ctx1abcd

377 

378 

379def _extract_passiveValve_intervals(bpod: dict) -> np.array: (empty)

380 # passiveValve.intervals 

381 # Get valve intervals from bpod channel 

382 # bpod channel should only contain valve output for passiveCW protocol 

383 # All high fronts == valve open times and low fronts == valve close times 

384 valveOn_times = bpod["times"][bpod["polarities"] > 0] 4 ctx1abcd

385 valveOff_times = bpod["times"][bpod["polarities"] < 0] 4 ctx1abcd

386 

387 assert len(valveOn_times) == NVALVE, "Wrong number of valve ONSET times" 4 ctx1abcd

388 assert len(valveOff_times) == NVALVE, "Wrong number of valve OFFSET times" 4 ctx1abcd

389 assert len(bpod["times"]) == NVALVE * 2, "Wrong number of valve FRONTS detected" # (40 * 2) 4 ctx1abcd

390 

391 # check all values are within bpod tolerance of 100µs 

392 assert np.allclose( 4 ctx1abcd

393 valveOff_times - valveOn_times, valveOff_times[1] - valveOn_times[1], atol=0.0001 

394 ), "Some valve outputs are longer or shorter than others" 

395 

396 return np.array([(x, y) for x, y in zip(valveOn_times, valveOff_times)]) 4 ctx1abcd

397 

398 

399def _extract_passiveAudio_intervals(audio: dict, rig_version: str) -> Tuple[np.array, np.array]: (empty)

400 

401 # make an exception for task version = 6.2.5 where things are strange but data is recoverable 

402 if rig_version == '6.2.5': 4 ctx1abcd

403 # Get all sound onsets and offsets 

404 soundOn_times = audio["times"][audio["polarities"] > 0] 

405 soundOff_times = audio["times"][audio["polarities"] < 0] 

406 

407 # Have a couple that are wayyy too long! 

408 time_threshold = 10 

409 diff = soundOff_times - soundOn_times 

410 stupid = np.where(diff > time_threshold)[0] 

411 NREMOVE = len(stupid) 

412 not_stupid = np.where(diff < time_threshold)[0] 

413 

414 assert len(soundOn_times) == NTONES + NNOISES - NREMOVE, "Wrong number of sound ONSETS" 

415 assert len(soundOff_times) == NTONES + NNOISES - NREMOVE, "Wrong number of sound OFFSETS" 

416 

417 soundOn_times = soundOn_times[not_stupid] 

418 soundOff_times = soundOff_times[not_stupid] 

419 

420 diff = soundOff_times - soundOn_times 

421 # Tone is ~100ms so check if diff < 0.3 

422 toneOn_times = soundOn_times[diff <= 0.3] 

423 toneOff_times = soundOff_times[diff <= 0.3] 

424 # Noise is ~500ms so check if diff > 0.3 

425 noiseOn_times = soundOn_times[diff > 0.3] 

426 noiseOff_times = soundOff_times[diff > 0.3] 

427 

428 # append with nans 

429 toneOn_times = np.r_[toneOn_times, np.full((NTONES - len(toneOn_times)), np.NAN)] 

430 toneOff_times = np.r_[toneOff_times, np.full((NTONES - len(toneOff_times)), np.NAN)] 

431 noiseOn_times = np.r_[noiseOn_times, np.full((NNOISES - len(noiseOn_times)), np.NAN)] 

432 noiseOff_times = np.r_[noiseOff_times, np.full((NNOISES - len(noiseOff_times)), np.NAN)] 

433 

434 else: 

435 # Get all sound onsets and offsets 

436 soundOn_times = audio["times"][audio["polarities"] > 0] 4 ctx1abcd

437 soundOff_times = audio["times"][audio["polarities"] < 0] 4 ctx1abcd

438 

439 # Check they are the correct number 

440 assert len(soundOn_times) == NTONES + NNOISES, f"Wrong number of sound ONSETS, " \ 4 ctx1abcd

441 f"{len(soundOn_times)}/{NTONES + NNOISES}" 

442 assert len(soundOff_times) == NTONES + NNOISES, f"Wrong number of sound OFFSETS, " \ 4 ctx1abcd

443 f"{len(soundOn_times)}/{NTONES + NNOISES}" 

444 

445 diff = soundOff_times - soundOn_times 4 ctx1abcd

446 # Tone is ~100ms so check if diff < 0.3 

447 toneOn_times = soundOn_times[diff <= 0.3] 4 ctx1abcd

448 toneOff_times = soundOff_times[diff <= 0.3] 4 ctx1abcd

449 # Noise is ~500ms so check if diff > 0.3 

450 noiseOn_times = soundOn_times[diff > 0.3] 4 ctx1abcd

451 noiseOff_times = soundOff_times[diff > 0.3] 4 ctx1abcd

452 

453 assert len(toneOn_times) == NTONES 4 ctx1abcd

454 assert len(toneOff_times) == NTONES 4 ctx1abcd

455 assert len(noiseOn_times) == NNOISES 4 ctx1abcd

456 assert len(noiseOff_times) == NNOISES 4 ctx1abcd

457 

458 # Fixed delays from soundcard ~500µs 

459 np.allclose(toneOff_times - toneOn_times, 0.1, atol=0.0006) 4 ctx1abcd

460 np.allclose(noiseOff_times - noiseOn_times, 0.5, atol=0.0006) 4 ctx1abcd

461 

462 passiveTone_intervals = np.append( 4 ctx1abcd

463 toneOn_times.reshape((len(toneOn_times), 1)), 

464 toneOff_times.reshape((len(toneOff_times), 1)), 

465 axis=1, 

466 ) 

467 passiveNoise_intervals = np.append( 4 ctx1abcd

468 noiseOn_times.reshape((len(noiseOn_times), 1)), 

469 noiseOff_times.reshape((len(noiseOff_times), 1)), 

470 axis=1, 

471 ) 

472 return passiveTone_intervals, passiveNoise_intervals 4 ctx1abcd

473 

474 

475# ------------------------------------------------------------------ 

476def extract_passive_periods(session_path: str, sync_collection: str = 'raw_ephys_data', sync: dict = None, (empty)

477 sync_map: dict = None, tmin=None, tmax=None) -> pd.DataFrame: 

478 

479 if sync is None or sync_map is None: 5 ctx1aebcd

480 sync, sync_map = ephys_fpga.get_sync_and_chn_map(session_path, sync_collection) 

481 

482 t_start_passive, t_starts, t_ends = _get_passive_spacers( 5 ctx1aebcd

483 session_path, sync_collection, sync=sync, sync_map=sync_map, tmin=tmin, tmax=tmax 

484 ) 

485 t_starts_col = np.insert(t_starts, 0, t_start_passive) 5 ctx1aebcd

486 t_ends_col = np.insert(t_ends, 0, t_ends[-1]) 5 ctx1aebcd

487 # tpassive_protocol = [t_start_passive, t_ends[-1]] 

488 # tspontaneous = [t_starts[0], t_ends[0]] 

489 # trfm = [t_starts[1], t_ends[1]] 

490 # treplay = [t_starts[2], t_ends[2]] 

491 passivePeriods_df = pd.DataFrame( 5 ctx1aebcd

492 [t_starts_col, t_ends_col], 

493 index=["start", "stop"], 

494 columns=["passiveProtocol", "spontaneousActivity", "RFM", "taskReplay"], 

495 ) 

496 return passivePeriods_df # _ibl_passivePeriods.intervalsTable.csv 5 ctx1aebcd

497 

498 

499def extract_rfmapping( (empty)

500 session_path: str, sync_collection: str = 'raw_ephys_data', task_collection: str = 'raw_passive_data', 

501 sync: dict = None, sync_map: dict = None, trfm: np.array = None 

502) -> Tuple[np.array, np.array]: 

503 meta = _load_passive_stim_meta() 5 ctx1aebcd

504 mkey = ( 5 ctx1aebcd

505 "VISUAL_STIM_" 

506 + {v: k for k, v in meta["VISUAL_STIMULI"].items()}["receptive_field_mapping"] 

507 ) 

508 if sync is None or sync_map is None: 5 ctx1aebcd

509 sync, sync_map = ephys_fpga.get_sync_and_chn_map(session_path, sync_collection) 

510 if trfm is None: 5 ctx1aebcd

511 passivePeriods_df = extract_passive_periods(session_path, sync_collection, sync=sync, sync_map=sync_map) 

512 trfm = passivePeriods_df.RFM.values 

513 

514 fttl = ephys_fpga.get_sync_fronts(sync, sync_map["frame2ttl"], tmin=trfm[0], tmax=trfm[1]) 5 ctx1aebcd

515 fttl = ephys_fpga._clean_frame2ttl(fttl) 5 ctx1aebcd

516 RF_file = Path().joinpath(session_path, task_collection, "_iblrig_RFMapStim.raw.bin") 5 ctx1aebcd

517 passiveRFM_frames, RF_ttl_trace = _reshape_RF(RF_file=RF_file, meta_stim=meta[mkey]) 5 ctx1aebcd

518 rf_id_up, rf_id_dw, RF_n_ttl_expected = _get_id_raisefall_from_analogttl(RF_ttl_trace) 5 ctx1aebcd

519 meta[mkey]["ttl_num"] = RF_n_ttl_expected 5 ctx1aebcd

520 rf_times_on_idx = np.where(np.diff(fttl["times"]) < 1)[0] 5 ctx1aebcd

521 rf_times_off_idx = rf_times_on_idx + 1 5 ctx1aebcd

522 RF_times = fttl["times"][np.sort(np.concatenate([rf_times_on_idx, rf_times_off_idx]))] 5 ctx1aebcd

523 RF_times_1 = RF_times[0::2] 5 ctx1aebcd

524 # Interpolate times for RF before outputting dataset 

525 passiveRFM_times = _interpolate_rf_mapping_stimulus( 5 ctx1aebcd

526 idxs_up=rf_id_up, 

527 idxs_dn=rf_id_dw, 

528 times=RF_times_1, 

529 Xq=np.arange(passiveRFM_frames.shape[0]), 

530 t_bin=1 / FRAME_FS, 

531 ) 

532 

533 return passiveRFM_times # _ibl_passiveRFM.times.npy 5 ctx1aebcd

534 

535 

536def extract_task_replay( (empty)

537 session_path: str, sync_collection: str = 'raw_ephys_data', task_collection: str = 'raw_passive_data', 

538 sync: dict = None, sync_map: dict = None, treplay: np.array = None 

539) -> Tuple[pd.DataFrame, pd.DataFrame]: 

540 

541 if sync is None or sync_map is None: 4 ctx1abcd

542 sync, sync_map = ephys_fpga.get_sync_and_chn_map(session_path, sync_collection) 

543 

544 if treplay is None: 4 ctx1abcd

545 passivePeriods_df = extract_passive_periods(session_path, sync_collection, sync=sync, sync_map=sync_map) 

546 treplay = passivePeriods_df.taskReplay.values 

547 

548 # TODO need to check this is okay 

549 fttl = ephys_fpga.get_sync_fronts(sync, sync_map["frame2ttl"], tmin=treplay[0], tmax=treplay[1]) 4 ctx1abcd

550 fttl = ephys_fpga._clean_frame2ttl(fttl) 4 ctx1abcd

551 passiveGabor_df = _extract_passiveGabor_df(fttl, session_path, task_collection=task_collection) 4 ctx1abcd

552 

553 bpod = ephys_fpga.get_sync_fronts(sync, sync_map["bpod"], tmin=treplay[0], tmax=treplay[1]) 4 ctx1abcd

554 passiveValve_intervals = _extract_passiveValve_intervals(bpod) 4 ctx1abcd

555 

556 task_version = _load_task_version(session_path, task_collection) 4 ctx1abcd

557 audio = ephys_fpga.get_sync_fronts(sync, sync_map["audio"], tmin=treplay[0], tmax=treplay[1]) 4 ctx1abcd

558 passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio, task_version) 4 ctx1abcd

559 

560 passiveStims_df = np.concatenate( 4 ctx1abcd

561 [passiveValve_intervals, passiveTone_intervals, passiveNoise_intervals], axis=1 

562 ) 

563 columns = ["valveOn", "valveOff", "toneOn", "toneOff", "noiseOn", "noiseOff"] 4 ctx1abcd

564 passiveStims_df = pd.DataFrame(passiveStims_df, columns=columns) 4 ctx1abcd

565 return ( 4 ctx1abcd

566 passiveGabor_df, 

567 passiveStims_df, 

568 ) # _ibl_passiveGabor.table.csv, _ibl_passiveStims.times_table.csv 

569 

570 

571def extract_replay_debug( (empty)

572 session_path: str, 

573 sync_collection: str = 'raw_ephys_data', 

574 task_collection: str = 'raw_passive_data', 

575 sync: dict = None, 

576 sync_map: dict = None, 

577 treplay: np.array = None, 

578 ax: plt.axes = None, 

579) -> Tuple[pd.DataFrame, pd.DataFrame]: 

580 # Load sessions sync channels, map 

581 if sync is None or sync_map is None: 

582 sync, sync_map = ephys_fpga.get_sync_and_chn_map(session_path, sync_collection) 

583 

584 if treplay is None: 

585 passivePeriods_df = extract_passive_periods(session_path, sync_collection=sync_collection, sync=sync, sync_map=sync_map) 

586 treplay = passivePeriods_df.taskReplay.values 

587 

588 if ax is None: 

589 f, ax = plt.subplots(1, 1) 

590 

591 f = ax.figure 

592 f.suptitle("/".join(str(session_path).split("/")[-5:])) 

593 plot_sync_channels(sync=sync, sync_map=sync_map, ax=ax) 

594 

595 passivePeriods_df = extract_passive_periods(session_path, sync_collection=sync_collection, sync=sync, sync_map=sync_map) 

596 treplay = passivePeriods_df.taskReplay.values 

597 

598 plot_passive_periods(passivePeriods_df, ax=ax) 

599 

600 fttl = ephys_fpga.get_sync_fronts(sync, sync_map["frame2ttl"], tmin=treplay[0]) 

601 passiveGabor_df = _extract_passiveGabor_df(fttl, session_path, task_collection=task_collection) 

602 plot_gabor_times(passiveGabor_df, ax=ax) 

603 

604 bpod = ephys_fpga.get_sync_fronts(sync, sync_map["bpod"], tmin=treplay[0]) 

605 passiveValve_intervals = _extract_passiveValve_intervals(bpod) 

606 plot_valve_times(passiveValve_intervals, ax=ax) 

607 

608 task_version = _load_task_version(session_path, task_collection) 

609 audio = ephys_fpga.get_sync_fronts(sync, sync_map["audio"], tmin=treplay[0]) 

610 passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio, task_version) 

611 plot_audio_times(passiveTone_intervals, passiveNoise_intervals, ax=ax) 

612 

613 passiveStims_df = np.concatenate( 

614 [passiveValve_intervals, passiveTone_intervals, passiveNoise_intervals], axis=1 

615 ) 

616 columns = ["valveOn", "valveOff", "toneOn", "toneOff", "noiseOn", "noiseOff"] 

617 passiveStims_df = pd.DataFrame(passiveStims_df, columns=columns) 

618 

619 return ( 

620 passiveGabor_df, 

621 passiveStims_df, 

622 ) # _ibl_passiveGabor.table.csv, _ibl_passiveStims.table.csv 

623 

624 

625# Main passiveCW extractor, calls all others 

626class PassiveChoiceWorld(BaseExtractor): (empty)

627 save_names = ( (empty)

628 "_ibl_passivePeriods.intervalsTable.csv", 

629 "_ibl_passiveRFM.times.npy", 

630 "_ibl_passiveGabor.table.csv", 

631 "_ibl_passiveStims.table.csv", 

632 ) 

633 var_names = ( (empty)

634 "passivePeriods_df", 

635 "passiveRFM_times", 

636 "passiveGabor_df", 

637 "passiveStims_df", 

638 ) 

639 

640 def _extract(self, sync_collection: str = 'raw_ephys_data', task_collection: str = 'raw_passive_data', sync: dict = None, (empty)

641 sync_map: dict = None, plot: bool = False, **kwargs) -> tuple: 

642 if sync is None or sync_map is None: 5 ctx1aebcd

643 sync, sync_map = ephys_fpga.get_sync_and_chn_map(self.session_path, sync_collection) 5 ctx1aebcd

644 

645 # Get the start and end times of this protocol 

646 if (protocol_number := kwargs.get('protocol_number')) is not None: # look for spacer 5 ctx1aebcd

647 # The spacers are TTLs generated by Bpod at the start of each protocol 

648 bpod = ephys_fpga.get_sync_fronts(sync, sync_map['bpod']) 1 ctx1b

649 tmin, tmax = ephys_fpga.get_protocol_period(self.session_path, protocol_number, bpod) 1 ctx1b

650 else: 

651 tmin = tmax = None 4 ctx1aecd

652 

653 # Passive periods 

654 passivePeriods_df = extract_passive_periods(self.session_path, sync_collection=sync_collection, sync=sync, 5 ctx1aebcd

655 sync_map=sync_map, tmin=tmin, tmax=tmax) 

656 trfm = passivePeriods_df.RFM.values 5 ctx1aebcd

657 treplay = passivePeriods_df.taskReplay.values 5 ctx1aebcd