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 # The pregenerated session number has had many parameter names 

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

73 pars = map(settings.get, ['PRELOADED_SESSION_NUM', 'PREGENERATED_SESSION_NUM', 'SESSION_TEMPLATE_ID']) 4 ctx1abcd

74 ses_nb = next((k for k in pars if k is not None), None) 4 ctx1abcd

75 

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

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

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

79 

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

81 

82 fixture = { 4 ctx1abcd

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

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

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

86 } 

87 

88 return fixture 4 ctx1abcd

89 

90 

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

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

93 

94 :param session_path: the path to a session 

95 :type session_path: str 

96 :return: ibl rig task protocol version 

97 :rtype: str 

98 

99 """ 

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

101 ses_ver = settings["IBLRIG_VERSION"] 4 ctx1abcd

102 

103 return ses_ver 4 ctx1abcd

104 

105 

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

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

108 

109 :param session_path: the path to a session 

110 :type session_path: str 

111 :param task_collection: collection containing task data 

112 :type task_collection: str 

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

114 :rtype: bool 

115 """ 

116 

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

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

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

120 

121 return skip_replay 4 ctx1abcd

122 

123 

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

125 """load_passive_stim_meta Loads the passive protocol metadata 

126 

127 :return: metadata about passive protocol stimulus presentation 

128 :rtype: dict 

129 """ 

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

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

132 meta = json.load(f) 6 ctx1haebcd

133 

134 return meta 6 ctx1haebcd

135 

136 

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

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

139 """ 

140 Find timestamps of spacer signal. 

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

142 :type spacer_template: array-like 

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

144 :type jitter: float 

145 :param ttl_signal: 

146 :type ttl_signal: array-like 

147 :param t_quiet: seconds between spacer and next stim 

148 :type t_quiet: float 

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

150 :type thresh: float 

151 :return: times of spacer onset/offset 

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

153 """ 

154 diff_spacer_template = np.diff(spacer_template) 5 ctx1aebcd

155 # add jitter; 

156 # remove extreme values 

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

158 # diff ttl signal to compare to spacer_model 

159 dttl = np.diff(ttl_signal) 5 ctx1aebcd

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

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

162 # convolve cleaned diff ttl signal w/ spacer model 

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

164 # find spacer location 

165 idxs_spacer_middle = np.where( 5 ctx1aebcd

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

167 )[0] 

168 # adjust indices for 

169 # - `np.where` call above 

170 # - length of spacer_model 

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

172 idxs_spacer_middle += 2 - spacer_around 5 ctx1aebcd

173 

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

175 # and monotonically non-increasing afterwards 

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

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

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

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

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

181 

182 idxs_spacer_middle = idxs_spacer_middle[is_valid] 5 ctx1aebcd

183 

184 # pull out spacer times (middle) 

185 ts_spacer_middle = ttl_signal[idxs_spacer_middle] 5 ctx1aebcd

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

187 spacer_length = np.max(spacer_template) 5 ctx1aebcd

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

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

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

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

192 return spacer_times, conv_dttl 5 ctx1aebcd

193 

194 

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

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

197 """ 

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

199 returns t_passive_starts, t_starts, t_ends 

200 """ 

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

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

203 meta = _load_passive_stim_meta() 5 ctx1aebcd

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

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

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

207 spacer_template = ( 5 ctx1aebcd

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

209 ) 

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

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

212 spacer_times, conv_dttl = _get_spacer_times( 5 ctx1aebcd

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

214 ) 

215 

216 # Check correct number of spacers found 

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

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

219 error_nspacer = True 

220 # sometimes the first spacer is truncated 

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

222 # with a lower threshold value 

223 # Note: take *3 for some margin 

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

225 # Truncate signals 

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

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

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

229 # Find spacer location 

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

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

232 thresh = 2.0 

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

234 ddttl = ddttl[0:idx_nearend_spacer] 

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

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

237 if len(indx_middle) == 1: 

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

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

240 idx_end = int(indx_middle + spacer_around) + 1 

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

242 error_nspacer = False 

243 

244 if error_nspacer: 

245 raise ValueError( 

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

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

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

249 ) 

250 

251 if tmax is None: 5 ctx1aebcd

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

253 

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

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

256 

257 

258# 2/3 RFMapping stimuli 

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

260 """ 

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

262 :param ttl_01: 

263 :type ttl_01: array-like 

264 :param times: array of stimulus switch times 

265 :type times: array-like 

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

267 :type frames: array-like 

268 :param t_bin: period of screen refresh rate 

269 :type t_bin: float 

270 :return: tuple of (stim_times, stim_frames) 

271 """ 

272 

273 beg_extrap_val = -10001 6 ctx1gaebcd

274 end_extrap_val = -10000 6 ctx1gaebcd

275 

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

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

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

279 # uniform spacing outside boundaries of ttl signal 

280 # first values 

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

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

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

284 # end values 

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

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

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

288 return Tq 6 ctx1gaebcd

289 

290 

291def _get_id_raisefall_from_analogttl(ttl_01): (empty)

292 """ 

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

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

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

296 """ 

297 # Check values are 0, 1, -1 

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

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

300 else: 

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

302 d_ttl_01 = np.diff(ttl_01) 5 ctx1aebcd

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

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

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

306 return id_up, id_dw, n_ttl_expected 5 ctx1aebcd

307 

308 

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

310 """ 

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

312 where frame2ttl placed to create TTL trace. 

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

314 :param meta_stim: variable containing metadata information on RF 

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

316 """ 

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

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

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

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

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

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

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

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

325 return frames, ttl_analogtrace_01 5 ctx1aebcd

326 

327 

328# 3/3 Replay of task stimuli 

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

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

331 # Pulses are strictly alternating with intervals 

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

333 # trim edges of pulses 

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

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

336 # Highest max is of the intervals 

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

338 thresh = diff0[0] 

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

340 thresh = diff1[0] 4 ctx1abcd

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

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

343 # Check if any pulse has been missed 

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

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

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

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

348 insert_where = missing_where + 1 

349 missing_value = idx_start_stims[missing_where] + 2 

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

351 

352 idx_end_stims = idx_start_stims + 1 4 ctx1abcd

353 

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

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

356 # Check if we missed the first stim 

357 if len(start_times) < NGABOR: 4 ctx1abcd

358 first_stim_off_idx = idx_start_stims[0] - 1 4 ctx1abcd

359 # first_stim_on_idx = first_stim_off_idx - 1 

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

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

362 

363 # intervals dstype requires reshaping of start and end times 

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

365 

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

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

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

369 

370 assert ( 

371 len(passiveGabor_intervals) == NGABOR 

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

373 fixture = _load_passive_session_fixtures(session_path, task_collection) 4 ctx1abcd

374 passiveGabor_properties = fixture["pcs"] 4 ctx1abcd

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

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

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

378 return passiveGabor_df 4 ctx1abcd

379 

380 

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

382 # passiveValve.intervals 

383 # Get valve intervals from bpod channel 

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

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

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

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

388 

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

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

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

392 

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

394 assert np.allclose( 4 ctx1abcd

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

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

397 

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

399 

400 

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

402 

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

404 if rig_version == '6.2.5': 4 ctx1abcd

405 # Get all sound onsets and offsets 

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

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

408 

409 # Have a couple that are wayyy too long! 

410 time_threshold = 10 

411 diff = soundOff_times - soundOn_times 

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

413 NREMOVE = len(stupid) 

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

415 

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

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

418 

419 soundOn_times = soundOn_times[not_stupid] 

420 soundOff_times = soundOff_times[not_stupid] 

421 

422 diff = soundOff_times - soundOn_times 

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

424 toneOn_times = soundOn_times[diff <= 0.3] 

425 toneOff_times = soundOff_times[diff <= 0.3] 

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

427 noiseOn_times = soundOn_times[diff > 0.3] 

428 noiseOff_times = soundOff_times[diff > 0.3] 

429 

430 # append with nans 

431 toneOn_times = np.r_[toneOn_times, np.full((NTONES - len(toneOn_times)), np.nan)] 

432 toneOff_times = np.r_[toneOff_times, np.full((NTONES - len(toneOff_times)), np.nan)] 

433 noiseOn_times = np.r_[noiseOn_times, np.full((NNOISES - len(noiseOn_times)), np.nan)] 

434 noiseOff_times = np.r_[noiseOff_times, np.full((NNOISES - len(noiseOff_times)), np.nan)] 

435 

436 else: 

437 # Get all sound onsets and offsets 

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

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

440 

441 # Check they are the correct number 

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

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

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

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

446 

447 diff = soundOff_times - soundOn_times 4 ctx1abcd

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

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

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

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

452 noiseOn_times = soundOn_times[diff > 0.3] 4 ctx1abcd

453 noiseOff_times = soundOff_times[diff > 0.3] 4 ctx1abcd

454 

455 assert len(toneOn_times) == NTONES 4 ctx1abcd

456 assert len(toneOff_times) == NTONES 4 ctx1abcd

457 assert len(noiseOn_times) == NNOISES 4 ctx1abcd

458 assert len(noiseOff_times) == NNOISES 4 ctx1abcd

459 

460 # Fixed delays from soundcard ~500µs 

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

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

463 

464 passiveTone_intervals = np.append( 4 ctx1abcd

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

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

467 axis=1, 

468 ) 

469 passiveNoise_intervals = np.append( 4 ctx1abcd

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

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

472 axis=1, 

473 ) 

474 return passiveTone_intervals, passiveNoise_intervals 4 ctx1abcd

475 

476 

477# ------------------------------------------------------------------ 

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

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

480 

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

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

483 

484 t_start_passive, t_starts, t_ends = _get_passive_spacers( 5 ctx1aebcd

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

486 ) 

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

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

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

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

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

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

493 passivePeriods_df = pd.DataFrame( 5 ctx1aebcd

494 [t_starts_col, t_ends_col], 

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

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

497 ) 

498 return passivePeriods_df # _ibl_passivePeriods.intervalsTable.csv 5 ctx1aebcd

499 

500 

501def extract_rfmapping( (empty)

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

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

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

505 meta = _load_passive_stim_meta() 5 ctx1aebcd

506 mkey = ( 5 ctx1aebcd

507 "VISUAL_STIM_" 

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

509 ) 

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

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

512 if trfm is None: 5 ctx1aebcd

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

514 trfm = passivePeriods_df.RFM.values 

515 

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

517 fttl = ephys_fpga._clean_frame2ttl(fttl) 5 ctx1aebcd

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

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

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

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

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

523 rf_times_off_idx = rf_times_on_idx + 1 5 ctx1aebcd

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

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

526 # Interpolate times for RF before outputting dataset 

527 passiveRFM_times = _interpolate_rf_mapping_stimulus( 5 ctx1aebcd

528 idxs_up=rf_id_up, 

529 idxs_dn=rf_id_dw, 

530 times=RF_times_1, 

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

532 t_bin=1 / FRAME_FS, 

533 ) 

534 

535 return passiveRFM_times # _ibl_passiveRFM.times.npy 5 ctx1aebcd

536 

537 

538def extract_task_replay( (empty)

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

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

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

542 

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

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

545 

546 if treplay is None: 4 ctx1abcd

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

548 treplay = passivePeriods_df.taskReplay.values 

549 

550 # TODO need to check this is okay 

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

552 fttl = ephys_fpga._clean_frame2ttl(fttl) 4 ctx1abcd

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

554 

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

556 passiveValve_intervals = _extract_passiveValve_intervals(bpod) 4 ctx1abcd

557 

558 task_version = _load_task_version(session_path, task_collection) 4 ctx1abcd

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

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

561 

562 passiveStims_df = np.concatenate( 4 ctx1abcd

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

564 ) 

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

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

567 return ( 4 ctx1abcd

568 passiveGabor_df, 

569 passiveStims_df, 

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

571 

572 

573def extract_replay_debug( (empty)

574 session_path: str, 

575 sync_collection: str = 'raw_ephys_data', 

576 task_collection: str = 'raw_passive_data', 

577 sync: dict = None, 

578 sync_map: dict = None, 

579 treplay: np.array = None, 

580 ax: plt.axes = None, 

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

582 # Load sessions sync channels, map 

583 if sync is None or sync_map is None: 

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

585 

586 if treplay is None: 

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

588 treplay = passivePeriods_df.taskReplay.values 

589 

590 if ax is None: 

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

592 

593 f = ax.figure 

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

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

596 

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

598 treplay = passivePeriods_df.taskReplay.values 

599 

600 plot_passive_periods(passivePeriods_df, ax=ax) 

601 

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

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

604 plot_gabor_times(passiveGabor_df, ax=ax) 

605 

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

607 passiveValve_intervals = _extract_passiveValve_intervals(bpod) 

608 plot_valve_times(passiveValve_intervals, ax=ax) 

609 

610 task_version = _load_task_version(session_path, task_collection) 

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

612 passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio, task_version) 

613 plot_audio_times(passiveTone_intervals, passiveNoise_intervals, ax=ax) 

614 

615 passiveStims_df = np.concatenate( 

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

617 ) 

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

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

620 

621 return ( 

622 passiveGabor_df, 

623 passiveStims_df, 

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

625 

626 

627# Main passiveCW extractor, calls all others 

628class PassiveChoiceWorld(BaseExtractor): (empty)

629 save_names = ( (empty)

630 "_ibl_passivePeriods.intervalsTable.csv", 

631 "_ibl_passiveRFM.times.npy", 

632 "_ibl_passiveGabor.table.csv", 

633 "_ibl_passiveStims.table.csv", 

634 ) 

635 var_names = ( (empty)

636 "passivePeriods_df", 

637 "passiveRFM_times", 

638 "passiveGabor_df", 

639 "passiveStims_df", 

640 ) 

641 

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

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

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

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

646 

647 # Get the start and end times of this protocol 

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

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

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

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

652 else: 

653 tmin = tmax = None 4 ctx1aecd

654 

655 # Passive periods 

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

657 sync_map=sync_map, tmin=tmin, tmax=tmax)