Coverage for ibllib/qc/task_metrics.py: 97%

410 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-07-02 18:55 +0100

1"""Behaviour QC. 

2 

3This module runs a list of quality control metrics on the behaviour data. 

4 

5.. warning:: 

6 The QC should be loaded using :meth:`ibllib.pipes.base_tasks.BehaviourTask.run_qc` and not 

7 instantiated directly. 

8 

9Examples 

10-------- 

11Running on a behaviour rig computer and updating QC fields in Alyx: 

12 

13>>> from ibllib.qc.task_qc_viewer.task_qc import show_session_task_qc 

14>>> qc = show_session_task_qc(session_path, bpod_only=True, local=True) # must close Viewer window 

15>>> qc = qc.run(update=True) 

16 

17Downloading the required data and inspecting the QC on a different computer: 

18 

19>>> from ibllib.pipes.dynamic_pipeline import get_trials_tasks 

20>>> from one.api import ONE 

21>>> task = get_trials_tasks(session_path, one=ONE())[0] # get first task run 

22>>> task.location = 'remote' 

23>>> task.setUp() # download required data 

24>>> qc = task.run_qc(update=False) 

25>>> outcome, results = qc.run() 

26 

27Inspecting individual test outcomes 

28 

29>>> outcome, results, outcomes = qc.compute_session_status() 

30 

31Running bpod QC on ephys session (when not on behaviour rig PC) 

32 

33>>> from ibllib.qc.task_qc_viewer.task_qc import get_bpod_trials_task, get_trials_tasks 

34>>> from one.api import ONE 

35>>> tasks = get_trials_tasks(session_path, one=ONE()) 

36>>> task = get_bpod_trials_task(tasks[0]) # Ensure Bpod only on behaviour rig 

37>>> task.location = 'remote' 

38>>> task.setUp() # download required data 

39>>> qc = task.run_qc(update=False) 

40>>> outcome, results = qc.run() 

41 

42Running ephys QC, from local server PC (after ephys + bpod data have been copied to a same folder) 

43 

44>>> from ibllib.pipes.dynamic_pipeline import get_trials_tasks 

45>>> task = get_trials_tasks(session_path, one=ONE())[0] # get first task run 

46>>> qc = task.run_qc(update=False) 

47>>> outcome, results = qc.run() 

48""" 

49import logging 

50import sys 

51from packaging import version 

52from pathlib import PurePosixPath 

53from datetime import datetime, timedelta 

54from inspect import getmembers, isfunction 

55from functools import reduce 

56from collections.abc import Sized 

57 

58import numpy as np 

59from scipy.stats import chisquare 

60 

61from brainbox.behavior.wheel import cm_to_rad, traces_by_trial 

62from ibllib.io.extractors import ephys_fpga 

63from one.alf import spec 

64from . import base 

65 

66_log = logging.getLogger(__name__) 

67 

68# todo the 2 followint parameters should be read from the task parameters for each session 

69ITI_DELAY_SECS = .5 

70FEEDBACK_NOGO_DELAY_SECS = 2 

71 

72BWM_CRITERIA = { 

73 'default': {'PASS': 0.99, 'WARNING': 0.90, 'FAIL': 0}, # Note: WARNING was 0.95 prior to Aug 2022 

74 'stimOff_itiIn_delays': {'PASS': 0.99, 'WARNING': 0}, 

75 'positive_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0}, 

76 'negative_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0}, 

77 'wheel_move_during_closed_loop': {'PASS': 0.99, 'WARNING': 0}, 

78 'response_stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0}, 

79 'detected_wheel_moves': {'PASS': 0.99, 'WARNING': 0}, 

80 'trial_length': {'PASS': 0.99, 'WARNING': 0}, 

81 'goCue_delays': {'PASS': 0.99, 'WARNING': 0}, 

82 'errorCue_delays': {'PASS': 0.99, 'WARNING': 0}, 

83 'stimOn_delays': {'PASS': 0.99, 'WARNING': 0}, 

84 'stimOff_delays': {'PASS': 0.99, 'WARNING': 0}, 

85 'stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0}, 

86 'iti_delays': {'NOT_SET': 0}, 

87 'passed_trial_checks': {'NOT_SET': 0} 

88} 

89 

90 

91def compute_session_status_from_dict(results, criteria=None): 

92 """ 

93 Compute overall task QC value from QC check results. 

94 

95 Given a dictionary of results, computes the overall session QC for each key and aggregates 

96 in a single value. 

97 

98 Parameters 

99 ---------- 

100 results : dict 

101 A dictionary of QC keys containing (usually scalar) values. 

102 criteria : dict 

103 A dictionary of qc keys containing map of PASS, WARNING, FAIL thresholds. 

104 

105 Returns 

106 ------- 

107 one.alf.spec.QC 

108 Overall session QC outcome. 

109 dict 

110 A map of QC tests and their outcomes. 

111 """ 

112 if not criteria: 1axmlngcbdei

113 criteria = {'default': BWM_CRITERIA['default']} 

114 outcomes = {k: TaskQC.thresholding(v, thresholds=criteria.get(k, criteria['default'])) 1axmlngcbdei

115 for k, v in results.items()} 

116 

117 # Criteria map is in order of severity so the max index is our overall QC outcome 

118 session_outcome = base.QC.overall_outcome(outcomes.values()) 1amlngcbdei

119 return session_outcome, outcomes 1amlngcbdei

120 

121 

122def update_dataset_qc(qc, registered_datasets, one, override=False): 

123 """ 

124 Update QC values for individual datasets. 

125 

126 Parameters 

127 ---------- 

128 qc : ibllib.qc.task_metrics.TaskQC 

129 A TaskQC object that has been run. 

130 registered_datasets : list of dict 

131 A list of Alyx dataset records. 

132 one : one.api.OneAlyx 

133 An online instance of ONE. 

134 override : bool 

135 If True the QC field is updated even if new value is better than previous. 

136 

137 Returns 

138 ------- 

139 list of dict 

140 The list of registered datasets but with the 'qc' fields updated. 

141 """ 

142 # Create map of dataset name, sans extension, to dataset id 

143 stem2id = {PurePosixPath(dset['name']).stem: dset.get('id') for dset in registered_datasets} 1hc

144 # Ensure dataset stems are unique 

145 assert len(stem2id) == len(registered_datasets), 'ambiguous dataset names' 1hc

146 

147 # dict of QC check to outcome (as enum value) 

148 *_, outcomes = qc.compute_session_status() 1hc

149 # work over map of dataset name (sans extension) to outcome (enum or dict of columns: enum) 

150 for name, outcome in qc.compute_dataset_qc_status(outcomes, qc.namespace).items(): 1hc

151 # if outcome is a dict, calculate aggregate outcome for each column 

152 if isinstance(outcome, dict): 1hc

153 extended_qc = outcome 1hc

154 outcome = qc.overall_outcome(outcome.values()) 1hc

155 else: 

156 extended_qc = {} 1hc

157 # check if dataset was registered to Alyx 

158 if not (did := stem2id.get(name)): 1hc

159 _log.debug('dataset %s not registered, skipping', name) 1h

160 continue 1h

161 # update the dataset QC value on Alyx 

162 if outcome > spec.QC.NOT_SET or override: 1hc

163 dset_qc = base.QC(did, one=one, log=_log, endpoint='datasets') 1hc

164 dset = next(x for x in registered_datasets if did == x.get('id')) 1hc

165 dset['qc'] = dset_qc.update(outcome, namespace='', override=override).name 1hc

166 if extended_qc: 1hc

167 dset_qc.update_extended_qc(extended_qc) 1hc

168 return registered_datasets 1hc

169 

170 

171class TaskQC(base.QC): 

172 """Task QC for training, biased, and ephys choice world.""" 

173 

174 criteria = BWM_CRITERIA 

175 

176 extractor = None 

177 """ibllib.qc.task_extractors.TaskQCExtractor: A task extractor object containing raw and extracted data.""" 

178 

179 @staticmethod 

180 def thresholding(qc_value, thresholds=None) -> spec.QC: 

181 """ 

182 Compute the outcome of a single key by applying thresholding. 

183 

184 Parameters 

185 ---------- 

186 qc_value : float 

187 Proportion of passing qcs, between 0 and 1. 

188 thresholds : dict 

189 Dictionary with keys 'PASS', 'WARNING', 'FAIL', (or enum 

190 integers, c.f. one.alf.spec.QC). 

191 

192 Returns 

193 ------- 

194 one.alf.spec.QC 

195 The outcome. 

196 """ 

197 thresholds = {spec.QC.validate(k): v for k, v in thresholds.items() or {}} 1axmlngcbdei

198 MAX_BOUND, MIN_BOUND = (1, 0) 1axmlngcbdei

199 if qc_value is None or np.isnan(qc_value): 1axmlngcbdei

200 return spec.QC.NOT_SET 1amlngbei

201 elif (qc_value > MAX_BOUND) or (qc_value < MIN_BOUND): 1axmlngcbdei

202 raise ValueError('Values out of bound') 1x

203 for crit in filter(None, sorted(spec.QC)): 1amlngcbdei

204 if crit in thresholds.keys() and qc_value >= thresholds[crit]: 1amlngcbdei

205 return crit 1amlngcbdei

206 # if None of this applies, return 'NOT_SET' 

207 return spec.QC.NOT_SET 1algcbdei

208 

209 def __init__(self, session_path_or_eid, namespace='task', **kwargs): 

210 """ 

211 Task QC for training, biased, and ephys choice world. 

212 

213 Parameters 

214 ---------- 

215 session_path_or_eid : str, uuid.UUID, pathlib.Path 

216 A session eid or path. If a path, it must be a valid ALFPath. 

217 namespace : str 

218 The namespace of the QC fields in the Alyx JSON field. 

219 log : logging.Logger, optional 

220 A logging.Logger instance, if None the 'ibllib.qc.base' logger is used. 

221 one : one.api.OneAlyx, optional 

222 An ONE instance for fetching and setting the QC on Alyx. 

223 """ 

224 super().__init__(session_path_or_eid, **kwargs) 1ahcbd

225 

226 # Data 

227 self.extractor = None 1ahcbd

228 

229 # Metrics and passed trials 

230 self.metrics = None 1ahcbd

231 self.passed = None 1ahcbd

232 

233 # Criteria (initialize as outcomes vary by class, task, and hardware) 

234 self.namespace = namespace 1ahcbd

235 self.criteria = {k if k == 'default' else f'_{namespace}_{k}': v 1ahcbd

236 for k, v in BWM_CRITERIA.items()} 

237 

238 def compute(self, **kwargs): 

239 """Compute and store the QC metrics. 

240 

241 Runs the QC on the session and stores a map of the metrics for each datapoint for each 

242 test, and a map of which datapoints passed for each test. 

243 

244 Parameters 

245 ---------- 

246 bpod_only : bool 

247 If True no data is extracted from the FPGA for ephys sessions. 

248 """ 

249 assert self.extractor is not None 1acbdfe

250 ns = kwargs.pop('namespace', self.namespace) 1acbdfe

251 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1acbdfe

252 if version.parse(ver) >= version.parse('8.0.0'): 1acbdfe

253 self.criteria[f'_{ns}_iti_delays'] = {'PASS': 0.99, 'WARNING': 0} 

254 self.criteria[f'_{ns}_passed_trial_checks'] = {'PASS': 0.7, 'WARNING': 0} 

255 else: 

256 self.criteria[f'_{ns}_iti_delays'] = {'NOT_SET': 0} 1acbdfe

257 self.criteria[f'_{ns}_passed_trial_checks'] = {'NOT_SET': 0} 1acbdfe

258 

259 self.log.info(f'Session {self.session_path}: Running QC on behavior data...') 1acbdfe

260 self.get_bpodqc_metrics_frame( 1acbdfe

261 self.extractor.data, 

262 wheel_gain=self.extractor.settings['STIM_GAIN'], # The wheel gain 

263 photodiode=self.extractor.frame_ttls, 

264 audio=self.extractor.audio_ttls, 

265 re_encoding=self.extractor.wheel_encoding or 'X1', 

266 min_qt=self.extractor.settings.get('QUIESCENT_PERIOD') or 0.2, 

267 audio_output=self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown') 

268 ) 

269 

270 def _get_checks(self): 

271 """ 

272 Find all methods that begin with 'check_'. 

273 

274 Returns 

275 ------- 

276 Dict[str, function] 

277 A map of QC check function names and the corresponding functions that return `metric` 

278 (any), `passed` (bool). 

279 """ 

280 def is_metric(x): 1athcbdfe

281 return isfunction(x) and x.__name__.startswith('check_') 1athcbdfe

282 

283 return dict(getmembers(sys.modules[__name__], is_metric)) 1athcbdfe

284 

285 def get_bpodqc_metrics_frame(self, data, **kwargs): 

286 """ 

287 Evaluate task QC metrics. 

288 

289 Evaluates all the QC metric functions in this module (those starting with 'check') and 

290 returns the results. The optional kwargs listed below are passed to each QC metric function. 

291 

292 Parameters 

293 ---------- 

294 data : dict 

295 The extracted task data. 

296 re_encoding : str {'X1', 'X2', 'X4'} 

297 The encoding configuration of the rotary encoder. 

298 enc_res : int 

299 The rotary encoder resolution as number of fronts per revolution. 

300 wheel_gain : float 

301 The STIM_GAIN task parameter. 

302 photodiode : dict 

303 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel. 

304 audio : dict 

305 The fronts from Bpod's BNC2 input FPGA audio sync channel. 

306 min_qt : float 

307 The QUIESCENT_PERIOD task parameter. 

308 

309 Returns 

310 ------- 

311 dict 

312 Map of checks and their QC metric values (1 per trial). 

313 dict 

314 Map of checks and a float array of which samples passed. 

315 """ 

316 # Find all methods that begin with 'check_' 

317 checks = self._get_checks() 1acbdfe

318 prefix = f'_{kwargs.pop("namespace", self.namespace)}_' # Extended QC fields will start with this 1acbdfe

319 # Method 'check_foobar' stored with key '_task_foobar' in metrics map 

320 qc_metrics_map = {prefix + k[6:]: fn(data, **kwargs) for k, fn in checks.items()} 1acbdfe

321 

322 # Split metrics and passed frames 

323 self.metrics = {} 1acbdfe

324 self.passed = {} 1acbdfe

325 for k in qc_metrics_map: 1acbdfe

326 self.metrics[k], self.passed[k] = qc_metrics_map[k] 1acbdfe

327 

328 # Add a check for trial level pass: did a given trial pass all checks? 

329 n_trials = data['intervals'].shape[0] 1acbdfe

330 # Trial-level checks return an array the length that equals the number of trials 

331 trial_level_passed = [m for m in self.passed.values() if isinstance(m, Sized) and len(m) == n_trials] 1acbdfe

332 name = prefix + 'passed_trial_checks' 1acbdfe

333 self.metrics[name] = reduce(np.logical_and, trial_level_passed or (None, None)) 1acbdfe

334 self.passed[name] = self.metrics[name].astype(float) if trial_level_passed else None 1acbdfe

335 

336 def run(self, update=False, **kwargs): 

337 """ 

338 Compute the QC outcomes and return overall task QC outcome. 

339 

340 Parameters 

341 ---------- 

342 update : bool 

343 If True, updates the session QC fields on Alyx. 

344 bpod_only : bool 

345 If True no data is extracted from the FPGA for ephys sessions. 

346 

347 Returns 

348 ------- 

349 str 

350 Overall task QC outcome. 

351 dict 

352 A map of QC tests and the proportion of data points that passed them. 

353 """ 

354 if self.metrics is None: 1acbdi

355 self.compute(**kwargs) 1acbd

356 outcome, results, _ = self.compute_session_status() 1acbdi

357 if update: 1acbdi

358 self.update_extended_qc(results) 1ci

359 self.update(outcome, kwargs.get('namespace', self.namespace)) 1ci

360 return outcome, results 1acbdi

361 

362 def compute_session_status(self): 

363 """ 

364 Compute the overall session QC for each key and aggregates in a single value. 

365 

366 Returns 

367 ------- 

368 str 

369 Overall session QC outcome. 

370 dict 

371 A map of QC tests and the proportion of data points that passed them. 

372 dict 

373 A map of QC tests and their outcomes. 

374 """ 

375 if self.passed is None: 1agcbdei

376 raise AttributeError('passed is None; compute QC first') 1e

377 # Get mean passed of each check, or None if passed is None or all NaN 

378 results = {k: None if v is None or np.isnan(v).all() else np.nanmean(v) 1agcbdei

379 for k, v in self.passed.items()} 

380 session_outcome, outcomes = compute_session_status_from_dict(results, self.criteria) 1agcbdei

381 return session_outcome, results, outcomes 1agcbdei

382 

383 @staticmethod 

384 def compute_dataset_qc_status(outcomes, namespace='task'): 

385 """Return map of dataset specific QC values. 

386 

387 Parameters 

388 ---------- 

389 outcomes : dict 

390 Map of checks and their individual outcomes. 

391 namespace : str 

392 The namespace of the QC fields in the Alyx JSON field. 

393 

394 Returns 

395 ------- 

396 dict 

397 Map of dataset names and their outcome. 

398 """ 

399 trials_table_outcomes = { 1thc

400 'intervals': outcomes.get(f'_{namespace}_iti_delays', spec.QC.NOT_SET), 

401 'goCue_times': outcomes.get(f'_{namespace}_goCue_delays', spec.QC.NOT_SET), 

402 'response_times': spec.QC.NOT_SET, 'choice': spec.QC.NOT_SET, 

403 'stimOn_times': outcomes.get(f'_{namespace}_stimOn_delays', spec.QC.NOT_SET), 

404 'contrastLeft': spec.QC.NOT_SET, 'contrastRight': spec.QC.NOT_SET, 

405 'feedbackType': spec.QC.NOT_SET, 'probabilityLeft': spec.QC.NOT_SET, 

406 'feedback_times': outcomes.get(f'_{namespace}_errorCue_delays', spec.QC.NOT_SET), 

407 'firstMovement_times': spec.QC.NOT_SET 

408 } 

409 reward_checks = (f'_{namespace}_reward_volumes', f'_{namespace}_reward_volume_set') 1thc

410 trials_table_outcomes['rewardVolume']: TaskQC.overall_outcome( 1thc

411 (outcomes.get(x, spec.QC.NOT_SET) for x in reward_checks) 

412 ) 

413 dataset_outcomes = { 1thc

414 '_ibl_trials.stimOff_times': outcomes.get(f'_{namespace}_stimOff_delays', spec.QC.NOT_SET), 

415 '_ibl_trials.table': trials_table_outcomes, 

416 } 

417 return dataset_outcomes 1thc

418 

419 

420class HabituationQC(TaskQC): 

421 """Task QC for habituation choice world.""" 

422 

423 def compute(self, **kwargs): 

424 """Compute and store the QC metrics. 

425 

426 Runs the QC on the session and stores a map of the metrics for each datapoint for each 

427 test, and a map of which datapoints passed for each test. 

428 """ 

429 assert self.extractor is not None 1gb

430 self.log.info(f'Session {self.session_path}: Running QC on habituation data...') 1gb

431 

432 # Initialize checks 

433 prefix = f'_{kwargs.pop("namespace", self.namespace)}_' 1gb

434 data = self.extractor.data 1gb

435 audio_output = self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown') 1gb

436 metrics = {} 1gb

437 passed = {} 1gb

438 

439 # Modify criteria based on version 

440 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1gb

441 is_v8 = version.parse(ver) >= version.parse('8.0.0') 1gb

442 self.criteria['_task_iti_delays'] = {'PASS': 0.99, 'WARNING': 0} if is_v8 else {'NOT_SET': 0} 1gb

443 

444 # Check all reward volumes == 3.0ul 

445 check = prefix + 'reward_volumes' 1gb

446 metrics[check] = data['rewardVolume'] 1gb

447 passed[check] = metrics[check] == 3.0 1gb

448 

449 # Check session durations are increasing in steps >= 12 minutes 

450 check = prefix + 'habituation_time' 1gb

451 if not self.one or not self.session_path: 1gb

452 self.log.warning('unable to determine session trials without ONE') 

453 metrics[check] = passed[check] = None 

454 else: 

455 subject, session_date = self.session_path.parts[-3:-1] 1gb

456 # compute from the date specified 

457 date_minus_week = ( 1gb

458 datetime.strptime(session_date, '%Y-%m-%d') - timedelta(days=7) 

459 ).strftime('%Y-%m-%d') 

460 sessions = self.one.alyx.rest('sessions', 'list', subject=subject, 1gb

461 date_range=[date_minus_week, session_date], 

462 task_protocol='habituation') 

463 # Remove the current session if already registered 

464 if sessions and sessions[0]['start_time'].startswith(session_date): 1gb

465 sessions = sessions[1:] 

466 metric = ([0, data['intervals'][-1, 1] - data['intervals'][0, 0]] + 1gb

467 [(datetime.fromisoformat(x['end_time']) - 

468 datetime.fromisoformat(x['start_time'])).total_seconds() 

469 for x in [self.one.alyx.get(s['url']) for s in sessions]]) 

470 

471 # The duration from raw trial data 

472 # duration = map(float, self.extractor.raw_data[-1]['elapsed_time'].split(':')) 

473 # duration = timedelta(**dict(zip(('hours', 'minutes', 'seconds'), 

474 # duration))).total_seconds() 

475 metrics[check] = np.array(metric) / 60 1gb

476 passed[check] = np.diff(metric) >= 12 1gb

477 

478 # Check event orders: trial_start < stim on < stim center < feedback < stim off 

479 check = prefix + 'trial_event_sequence' 1gb

480 nans = ( 1gb

481 np.isnan(data['intervals'][:, 0]) | # noqa 

482 np.isnan(data['stimOn_times']) | # noqa 

483 np.isnan(data['stimCenter_times']) | 

484 np.isnan(data['valveOpen_times']) | # noqa 

485 np.isnan(data['stimOff_times']) 

486 ) 

487 a = np.less(data['intervals'][:, 0], data['stimOn_times'], where=~nans) 1gb

488 b = np.less(data['stimOn_times'], data['stimCenter_times'], where=~nans) 1gb

489 c = np.less(data['stimCenter_times'], data['valveOpen_times'], where=~nans) 1gb

490 d = np.less(data['valveOpen_times'], data['stimOff_times'], where=~nans) 1gb

491 

492 metrics[check] = a & b & c & d & ~nans 1gb

493 passed[check] = metrics[check].astype(float) 1gb

494 

495 # Check that the time difference between the visual stimulus center-command being 

496 # triggered and the stimulus effectively appearing in the center is smaller than 150 ms. 

497 check = prefix + 'stimCenter_delays' 1gb

498 metric = np.nan_to_num(data['stimCenter_times'] - data['stimCenterTrigger_times'], 1gb

499 nan=np.inf) 

500 passed[check] = (metric <= 0.15) & (metric > 0) 1gb

501 metrics[check] = metric 1gb

502 

503 # Phase check 

504 check = prefix + 'phase' 1gb

505 metric = data['phase'] 1gb

506 passed[check] = (metric <= 2 * np.pi) & (metric >= 0) 1gb

507 metrics[check] = metric 1gb

508 

509 # This is not very useful as a check because there are so few trials 

510 check = prefix + 'phase_distribution' 1gb

511 metric, _ = np.histogram(data['phase']) 1gb

512 _, p = chisquare(metric) 1gb

513 passed[check] = p < 0.05 if len(data['phase']) >= 400 else None # skip if too few trials 1gb

514 metrics[check] = metric 1gb

515 

516 # Check that the period of gray screen between stim off and the start of the next trial is 

517 # 1s +/- 10%. 

518 check = prefix + 'iti_delays' 1gb

519 iti = (np.roll(data['stimOn_times'], -1) - data['stimOff_times'])[:-1] 1gb

520 metric = np.r_[np.nan_to_num(iti, nan=np.inf), np.nan] - 1. 1gb

521 passed[check] = np.abs(metric) <= 0.1 1gb

522 passed[check][-1] = np.nan 1gb

523 metrics[check] = metric 1gb

524 

525 # Check stim on go cue delay 

526 # Unlike in ChoiceWorld, the go cue is in a race condition with the stimulus onset and 

527 # therefore usually presents slightly before the stimulus. For now we will keep the QC 

528 # thresholds the same as in ChoiceWorld which means this will always fail 

529 check = prefix + 'stimOn_goCue_delays' 1gb

530 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold. 

531 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1gb

532 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 1gb

533 passed[check] = np.abs(metric) < threshold 1gb

534 metrics[check] = metric 1gb

535 

536 # Checks common to training QC 

537 checks = [check_goCue_delays, check_stimOn_delays, check_stimOff_delays] 1gb

538 for fcn in checks: 1gb

539 check = prefix + fcn.__name__[6:] 1gb

540 metrics[check], passed[check] = fcn(data, audio_output=audio_output) 1gb

541 

542 self.metrics, self.passed = (metrics, passed) 1gb

543 

544 

545# SINGLE METRICS 

546# ---------------------------------------------------------------------------- # 

547 

548# === Delays between events checks === 

549 

550def check_stimOn_goCue_delays(data, audio_output='harp', **_): 

551 """ 

552 Check the go cue tone occurs less than 10ms before stimulus on. 

553 

554 Checks that the time difference between the onset of the visual stimulus 

555 and the onset of the go cue tone is positive and less than 10ms. 

556 

557 Metric: 

558 M = stimOn_times - goCue_times 

559 

560 Criteria: 

561 0 < M < 0.010 s 

562 

563 Units: 

564 seconds [s] 

565 

566 :param data: dict of trial data with keys ('goCue_times', 'stimOn_times', 'intervals') 

567 :param audio_output: audio output device name. 

568 

569 Notes 

570 ----- 

571 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th 

572 percentile of delays over 500 training sessions using the Xonar soundcard. 

573 """ 

574 # Calculate the difference between stimOn and goCue times. 

575 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold. 

576 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1aBcbdfe

577 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 1aBcbdfe

578 passed = (metric < threshold) & (metric > 0) 1aBcbdfe

579 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aBcbdfe

580 return metric, passed 1aBcbdfe

581 

582 

583def check_response_feedback_delays(data, audio_output='harp', **_): 

584 """ 

585 Check the feedback delivered within 10ms of the response threshold. 

586 

587 Checks that the time difference between the response and the feedback onset 

588 (error sound or valve) is positive and less than 10ms. 

589 

590 Metric: 

591 M = feedback_time - response_time 

592 

593 Criterion: 

594 0 < M < 0.010 s 

595 

596 Units: 

597 seconds [s] 

598 

599 :param data: dict of trial data with keys ('feedback_times', 'response_times', 'intervals') 

600 :param audio_output: audio output device name. 

601 

602 Notes 

603 ----- 

604 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th 

605 percentile of delays over 500 training sessions using the Xonar soundcard. 

606 """ 

607 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1aCcbdfe

608 metric = np.nan_to_num(data['feedback_times'] - data['response_times'], nan=np.inf) 1aCcbdfe

609 passed = (metric < threshold) & (metric > 0) 1aCcbdfe

610 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aCcbdfe

611 return metric, passed 1aCcbdfe

612 

613 

614def check_response_stimFreeze_delays(data, **_): 

615 """ 

616 Check the stimulus freezes within 100ms of the expected time. 

617 

618 Checks that the time difference between the visual stimulus freezing and the 

619 response is positive and less than 100ms. 

620 

621 Metric: 

622 M = (stimFreeze_times - response_times) 

623 

624 Criterion: 

625 0 < M < 0.100 s 

626 

627 Units: 

628 seconds [s] 

629 

630 :param data: dict of trial data with keys ('stimFreeze_times', 'response_times', 'intervals', 

631 'choice') 

632 """ 

633 # Calculate the difference between stimOn and goCue times. 

634 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold. 

635 metric = np.nan_to_num(data['stimFreeze_times'] - data['response_times'], nan=np.inf) 1aDcbdfe

636 # Test for valid values 

637 passed = ((metric < 0.1) & (metric > 0)).astype(float) 1aDcbdfe

638 # Finally remove no_go trials (stimFreeze triggered differently in no_go trials) 

639 # These values are ignored in calculation of proportion passed 

640 passed[data['choice'] == 0] = np.nan 1aDcbdfe

641 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aDcbdfe

642 return metric, passed 1aDcbdfe

643 

644 

645def check_stimOff_itiIn_delays(data, **_): 

646 """Check that the start of the trial interval is within 10ms of the visual stimulus turning off. 

647 

648 Metric: 

649 M = itiIn_times - stimOff_times 

650 

651 Criterion: 

652 0 < M < 0.010 s 

653 

654 Units: 

655 seconds [s] 

656 

657 :param data: dict of trial data with keys ('stimOff_times', 'itiIn_times', 'intervals', 

658 'choice') 

659 """ 

660 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold. 

661 metric = np.nan_to_num(data['itiIn_times'] - data['stimOff_times'], nan=np.inf) 1aEcbdfe

662 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 1aEcbdfe

663 # Remove no_go trials (stimOff triggered differently in no_go trials) 

664 # NaN values are ignored in calculation of proportion passed 

665 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan 1aEcbdfe

666 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aEcbdfe

667 return metric, passed 1aEcbdfe

668 

669 

670def check_iti_delays(data, subtract_pauses=False, iti_delay_secs=ITI_DELAY_SECS, 

671 feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_): 

672 """ 

673 Check the open-loop grey screen period is approximately 1 second. 

674 

675 Check that the period of grey screen between stim off and the start of the next trial is 

676 1s +/- 10%. If the trial was paused during this time, the check will account for that 

677 

678 Metric: 

679 M = stimOff (n) - trialStart (n+1) - 1. 

680 

681 Criterion: 

682 |M| < 0.1 

683 

684 Units: 

685 seconds [s] 

686 

687 Parameters 

688 ---------- 

689 data : dict 

690 Trial data with keys ('stimOff_times', 'intervals', 'pause_duration'). 

691 subtract_pauses: bool 

692 If True, account for experimenter-initiated pauses between trials; if False, trials where 

693 the experimenter paused the task may fail this check. 

694 

695 Returns 

696 ------- 

697 numpy.array 

698 An array of metric values to threshold. 

699 numpy.array 

700 An array of boolean values, 1 per trial, where True means trial passes QC threshold. 

701 """ 

702 # Initialize array the length of completed trials 

703 metric = np.full(data['intervals'].shape[0], np.nan) 1aucbdfe

704 passed = metric.copy() 1aucbdfe

705 pauses = (data['pause_duration'] if subtract_pauses else np.zeros_like(metric))[:-1] 1aucbdfe

706 # Get the difference between stim off and the start of the next trial 

707 # Missing data are set to Inf, except for the last trial which is a NaN 

708 metric[:-1] = np.nan_to_num( 1aucbdfe

709 data['intervals'][1:, 0] - data['stimOff_times'][:-1] - iti_delay_secs - pauses, 

710 nan=np.inf 

711 ) 

712 metric[data['choice'] == 0] = metric[data['choice'] == 0] - feedback_nogo_delay_secs 1aucbdfe

713 passed[:-1] = np.abs(metric[:-1]) < (iti_delay_secs / 10) # Last trial is not counted 1aucbdfe

714 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aucbdfe

715 return metric, passed 1aucbdfe

716 

717 

718def check_positive_feedback_stimOff_delays(data, **_): 

719 """ 

720 Check stimulus offset occurs approximately 1 second after reward delivered. 

721 

722 Check that the time difference between the valve onset and the visual stimulus turning off 

723 is 1 ± 0.150 seconds. 

724 

725 Metric: 

726 M = stimOff_times - feedback_times - 1s 

727 

728 Criterion: 

729 |M| < 0.150 s 

730 

731 Units: 

732 seconds [s] 

733 

734 :param data: dict of trial data with keys ('stimOff_times', 'feedback_times', 'intervals', 

735 'correct') 

736 """ 

737 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold. 

738 metric = np.nan_to_num(data['stimOff_times'] - data['feedback_times'] - 1, nan=np.inf) 1aFcbdfe

739 passed = (np.abs(metric) < 0.15).astype(float) 1aFcbdfe

740 # NaN values are ignored in calculation of proportion passed; ignore incorrect trials here 

741 metric[~data['correct']] = passed[~data['correct']] = np.nan 1aFcbdfe

742 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aFcbdfe

743 return metric, passed 1aFcbdfe

744 

745 

746def check_negative_feedback_stimOff_delays(data, feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_): 

747 """ 

748 Check the stimulus offset occurs approximately 2 seconds after negative feedback delivery. 

749 

750 Check that the time difference between the error sound and the visual stimulus 

751 turning off is 2 ± 0.150 seconds. 

752 

753 Metric: 

754 M = stimOff_times - errorCue_times - 2s 

755 

756 Criterion: 

757 |M| < 0.150 s 

758 

759 Units: 

760 seconds [s] 

761 

762 :param data: dict of trial data with keys ('stimOff_times', 'errorCue_times', 'intervals') 

763 """ 

764 metric = np.nan_to_num(data['stimOff_times'] - data['errorCue_times'] - 2, nan=np.inf) 1azcbdfe

765 # for the nogo trials, the feedback is the same as the stimOff 

766 metric[data['choice'] == 0] = metric[data['choice'] == 0] + feedback_nogo_delay_secs 1azcbdfe

767 # Apply criteria 

768 passed = (np.abs(metric) < 0.15).astype(float) 1azcbdfe

769 # Remove none negative feedback trials 

770 metric[data['correct']] = passed[data['correct']] = np.nan 1azcbdfe

771 assert data['intervals'].shape[0] == len(metric) == len(passed) 1azcbdfe

772 return metric, passed 1azcbdfe

773 

774 

775# === Wheel movement during trial checks === 

776 

777def check_wheel_move_before_feedback(data, **_): 

778 """Check that the wheel does move within 100ms of the feedback onset (error sound or valve). 

779 

780 Metric: 

781 M = (w_t - 0.05) - (w_t + 0.05), where t = feedback_times 

782 

783 Criterion: 

784 M != 0 

785 

786 Units: 

787 radians 

788 

789 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice', 

790 'intervals', 'feedback_times') 

791 """ 

792 # Get tuple of wheel times and positions within 100ms of feedback 

793 traces = traces_by_trial( 1aocbdfe

794 data['wheel_timestamps'], 

795 data['wheel_position'], 

796 start=data['feedback_times'] - 0.05, 

797 end=data['feedback_times'] + 0.05, 

798 ) 

799 metric = np.zeros_like(data['feedback_times']) 1aocbdfe

800 # For each trial find the displacement 

801 for i, trial in enumerate(traces): 1aocbdfe

802 pos = trial[1] 1aocbdfe

803 if pos.size > 1: 1aocbdfe

804 metric[i] = pos[-1] - pos[0] 1aocbdfe

805 

806 # except no-go trials 

807 metric[data['choice'] == 0] = np.nan # NaN = trial ignored for this check 1aocbdfe

808 nans = np.isnan(metric) 1aocbdfe

809 passed = np.zeros_like(metric) * np.nan 1aocbdfe

810 

811 passed[~nans] = (metric[~nans] != 0).astype(float) 1aocbdfe

812 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aocbdfe

813 return metric, passed 1aocbdfe

814 

815 

816def _wheel_move_during_closed_loop(re_ts, re_pos, data, wheel_gain=None, tol=1, **_): 

817 """ 

818 Check the wheel moves the correct amount to reach threshold. 

819 

820 Check that the wheel moves by approximately 35 degrees during the closed-loop period 

821 on trials where a feedback (error sound or valve) is delivered. 

822 

823 Metric: 

824 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

825 time, w_t0 = position at go cue time, threshold_displacement = displacement required to 

826 move 35 visual degrees 

827 

828 Criterion: 

829 displacement < tol visual degree 

830 

831 Units: 

832 degrees angle of wheel turn 

833 

834 :param re_ts: extracted wheel timestamps in seconds 

835 :param re_pos: extracted wheel positions in radians 

836 :param data: a dict with the keys (goCueTrigger_times, response_times, feedback_times, 

837 position, choice, intervals) 

838 :param wheel_gain: the 'STIM_GAIN' task setting 

839 :param tol: the criterion in visual degrees 

840 """ 

841 if wheel_gain is None: 1ajcbdfe

842 _log.warning('No wheel_gain input in function call, returning None') 

843 return None, None 

844 

845 # Get tuple of wheel times and positions over each trial's closed-loop period 

846 traces = traces_by_trial(re_ts, re_pos, 1ajcbdfe

847 start=data['goCueTrigger_times'], 

848 end=data['response_times']) 

849 

850 metric = np.zeros_like(data['feedback_times']) 1ajcbdfe

851 # For each trial find the absolute displacement 

852 for i, trial in enumerate(traces): 1ajcbdfe

853 t, pos = trial 1ajcbdfe

854 if pos.size != 0: 1ajcbdfe

855 # Find the position of the preceding sample and subtract it 

856 idx = np.abs(re_ts - t[0]).argmin() - 1 1ajcbdfe

857 origin = re_pos[idx] 1ajcbdfe

858 metric[i] = np.abs(pos - origin).max() 1ajcbdfe

859 

860 # Load wheel_gain and thresholds for each trial 

861 wheel_gain = np.array([wheel_gain] * len(data['position'])) 1ajcbdfe

862 thresh = data['position'] 1ajcbdfe

863 # abs displacement, s, in mm required to move 35 visual degrees 

864 s_mm = np.abs(thresh / wheel_gain) # don't care about direction 1ajcbdfe

865 criterion = cm_to_rad(s_mm * 1e-1) # convert abs displacement to radians (wheel pos is in rad) 1ajcbdfe

866 metric = metric - criterion # difference should be close to 0 1ajcbdfe

867 rad_per_deg = cm_to_rad(1 / wheel_gain * 1e-1) 1ajcbdfe

868 passed = (np.abs(metric) < rad_per_deg * tol).astype(float) # less than 1 visual degree off 1ajcbdfe

869 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan # except no-go trials 1ajcbdfe

870 assert data['intervals'].shape[0] == len(metric) == len(passed) 1ajcbdfe

871 return metric, passed 1ajcbdfe

872 

873 

874def check_wheel_move_during_closed_loop(data, wheel_gain=None, **_): 

875 """ 

876 Check the wheel moves the correct amount to reach threshold. 

877 

878 Check that the wheel moves by approximately 35 degrees during the closed-loop period 

879 on trials where a feedback (error sound or valve) is delivered. 

880 

881 Metric: 

882 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

883 time, w_t0 = position at go cue time, threshold_displacement = displacement required to 

884 move 35 visual degrees 

885 

886 Criterion: 

887 displacement < 3 visual degrees 

888 

889 Units: 

890 degrees angle of wheel turn 

891 

892 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice', 

893 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position') 

894 :param wheel_gain: the 'STIM_GAIN' task setting 

895 """ 

896 # Get the Bpod extracted wheel data 

897 timestamps = data['wheel_timestamps'] 1ajcbdfe

898 position = data['wheel_position'] 1ajcbdfe

899 

900 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=3) 1ajcbdfe

901 

902 

903def check_wheel_move_during_closed_loop_bpod(data, wheel_gain=None, **_): 

904 """ 

905 Check the wheel moves the correct amount to reach threshold. 

906 

907 Check that the wheel moves by approximately 35 degrees during the closed-loop period 

908 on trials where a feedback (error sound or valve) is delivered. This check uses the Bpod 

909 wheel data (measured at a lower resolution) with a stricter tolerance (1 visual degree). 

910 

911 Metric: 

912 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

913 time, w_t0 = position at go cue time, threshold_displacement = displacement required to 

914 move 35 visual degrees. 

915 

916 Criterion: 

917 displacement < 1 visual degree 

918 

919 Units: 

920 degrees angle of wheel turn 

921 

922 :param data: dict of trial data with keys ('wheel_timestamps(_bpod)', 'wheel_position(_bpod)', 

923 'choice', 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position') 

924 :param wheel_gain: the 'STIM_GAIN' task setting 

925 """ 

926 # Get the Bpod extracted wheel data 

927 timestamps = data.get('wheel_timestamps_bpod', data['wheel_timestamps']) 1acbdfe

928 position = data.get('wheel_position_bpod', data['wheel_position']) 1acbdfe

929 

930 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=1) 1acbdfe

931 

932 

933def check_wheel_freeze_during_quiescence(data, **_): 

934 """ 

935 Check the wheel is indeed still during the quiescent period. 

936 

937 Check that the wheel does not move more than 2 degrees in each direction during the 

938 quiescence interval before the stimulus appears. 

939 

940 Metric: 

941 M = |max(W) - min(W)| where W is wheel pos over quiescence interval; 

942 interval = [stimOnTrigger_times - quiescent_duration, stimOnTrigger_times] 

943 

944 Criterion: 

945 M < 2 degrees 

946 

947 Units: 

948 degrees angle of wheel turn 

949 

950 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'quiescence', 

951 'intervals', 'stimOnTrigger_times') 

952 """ 

953 assert np.all(np.diff(data['wheel_timestamps']) >= 0) 1akcbdfe

954 assert data['quiescence'].size == data['stimOnTrigger_times'].size 1akcbdfe

955 # Get tuple of wheel times and positions over each trial's quiescence period 

956 qevt_start_times = data['stimOnTrigger_times'] - data['quiescence'] 1akcbdfe

957 traces = traces_by_trial( 1akcbdfe

958 data['wheel_timestamps'], 

959 data['wheel_position'], 

960 start=qevt_start_times, 

961 end=data['stimOnTrigger_times'] 

962 ) 

963 

964 metric = np.zeros((len(data['quiescence']), 2)) # (n_trials, n_directions) 1akcbdfe

965 for i, trial in enumerate(traces): 1akcbdfe

966 t, pos = trial 1akcbdfe

967 # Get the last position before the period began 

968 if pos.size > 0: 1akcbdfe

969 # Find the position of the preceding sample and subtract it 

970 idx = np.abs(data['wheel_timestamps'] - t[0]).argmin() - 1 1akcbdfe

971 origin = data['wheel_position'][idx if idx != -1 else 0] 1akcbdfe

972 # Find the absolute min and max relative to the last sample 

973 metric[i, :] = np.abs([np.min(pos - origin), np.max(pos - origin)]) 1akcbdfe

974 # Reduce to the largest displacement found in any direction 

975 metric = np.max(metric, axis=1) 1akcbdfe

976 metric = 180 * metric / np.pi # convert to degrees from radians 1akcbdfe

977 criterion = 2 # Position shouldn't change more than 2 in either direction 1akcbdfe

978 passed = metric < criterion 1akcbdfe

979 assert data['intervals'].shape[0] == len(metric) == len(passed) 1akcbdfe

980 return metric, passed 1akcbdfe

981 

982 

983def check_detected_wheel_moves(data, min_qt=0, **_): 

984 """Check that the detected first movement times are reasonable. 

985 

986 Metric: 

987 M = firstMovement times 

988 

989 Criterion: 

990 (goCue trigger time - min quiescent period) < M < response time 

991 

992 Units: 

993 Seconds [s] 

994 

995 :param data: dict of trial data with keys ('firstMovement_times', 'goCueTrigger_times', 

996 'response_times', 'choice', 'intervals') 

997 :param min_qt: the minimum possible quiescent period (the QUIESCENT_PERIOD task parameter) 

998 """ 

999 # Depending on task version this may be a single value or an array of quiescent periods 

1000 min_qt = np.array(min_qt) 1ascbdfe

1001 if min_qt.size > data['intervals'].shape[0]: 1ascbdfe

1002 min_qt = min_qt[:data['intervals'].shape[0]] 1d

1003 

1004 metric = data['firstMovement_times'] 1ascbdfe

1005 qevt_start = data['goCueTrigger_times'] - np.array(min_qt) 1ascbdfe

1006 response = data['response_times'] 1ascbdfe

1007 # First movement time for each trial should be after the quiescent period and before feedback 

1008 passed = np.array([a < m < b for m, a, b in zip(metric, qevt_start, response)], dtype=float) 1ascbdfe

1009 nogo = data['choice'] == 0 1ascbdfe

1010 passed[nogo] = np.nan # No go trial may have no movement times and that's fine 1ascbdfe

1011 return metric, passed 1ascbdfe

1012 

1013 

1014# === Sequence of events checks === 

1015 

1016def check_error_trial_event_sequence(data, **_): 

1017 """ 

1018 Check trial events occur in correct order for negative feedback trials. 

1019 

1020 Check that on incorrect / miss trials, there are exactly: 

1021 2 audio events (go cue sound and error sound) and 2 Bpod events (trial start, ITI), occurring 

1022 in the correct order 

1023 

1024 Metric: 

1025 M = Bpod (trial start) > audio (go cue) > audio (error) > Bpod (ITI) > Bpod (trial end) 

1026 

1027 Criterion: 

1028 M == True 

1029 

1030 Units: 

1031 -none- 

1032 

1033 :param data: dict of trial data with keys ('errorCue_times', 'goCue_times', 'intervals', 

1034 'itiIn_times', 'correct') 

1035 """ 

1036 # An array the length of N trials where True means at least one event time was NaN (bad) 

1037 nans = ( 1aqcbdfe

1038 np.isnan(data['intervals'][:, 0]) | 

1039 np.isnan(data['goCue_times']) | # noqa 

1040 np.isnan(data['errorCue_times']) | # noqa 

1041 np.isnan(data['itiIn_times']) | # noqa 

1042 np.isnan(data['intervals'][:, 1]) 

1043 ) 

1044 

1045 # For each trial check that the events happened in the correct order (ignore NaN values) 

1046 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1aqcbdfe

1047 b = np.less(data['goCue_times'], data['errorCue_times'], where=~nans) # Go cue < error cue 1aqcbdfe

1048 c = np.less(data['errorCue_times'], data['itiIn_times'], where=~nans) # Error cue < ITI start 1aqcbdfe

1049 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1aqcbdfe

1050 

1051 # For each trial check all events were in order AND all event times were not NaN 

1052 metric = a & b & c & d & ~nans 1aqcbdfe

1053 

1054 passed = metric.astype(float) 1aqcbdfe

1055 passed[data['correct']] = np.nan # Look only at incorrect trials 1aqcbdfe

1056 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aqcbdfe

1057 return metric, passed 1aqcbdfe

1058 

1059 

1060def check_correct_trial_event_sequence(data, **_): 

1061 """ 

1062 Check trial events occur in correct order for positive feedback trials. 

1063 

1064 Check that on correct trials, there are exactly: 

1065 1 audio events and 3 Bpod events (valve open, trial start, ITI), occurring in the correct order 

1066 

1067 Metric: 

1068 M = Bpod (trial start) > audio (go cue) > Bpod (valve) > Bpod (ITI) > Bpod (trial end) 

1069 

1070 Criterion: 

1071 M == True 

1072 

1073 Units: 

1074 -none- 

1075 

1076 :param data: dict of trial data with keys ('valveOpen_times', 'goCue_times', 'intervals', 

1077 'itiIn_times', 'correct') 

1078 """ 

1079 # An array the length of N trials where True means at least one event time was NaN (bad) 

1080 nans = ( 1arcbdfe

1081 np.isnan(data['intervals'][:, 0]) | 

1082 np.isnan(data['goCue_times']) | # noqa 

1083 np.isnan(data['valveOpen_times']) | 

1084 np.isnan(data['itiIn_times']) | # noqa 

1085 np.isnan(data['intervals'][:, 1]) 

1086 ) 

1087 

1088 # For each trial check that the events happened in the correct order (ignore NaN values) 

1089 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1arcbdfe

1090 b = np.less(data['goCue_times'], data['valveOpen_times'], where=~nans) # Go cue < feedback 1arcbdfe

1091 c = np.less(data['valveOpen_times'], data['itiIn_times'], where=~nans) # Feedback < ITI start 1arcbdfe

1092 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1arcbdfe

1093 

1094 # For each trial True means all events were in order AND all event times were not NaN 

1095 metric = a & b & c & d & ~nans 1arcbdfe

1096 

1097 passed = metric.astype(float) 1arcbdfe

1098 passed[~data['correct']] = np.nan # Look only at correct trials 1arcbdfe

1099 assert data['intervals'].shape[0] == len(metric) == len(passed) 1arcbdfe

1100 return metric, passed 1arcbdfe

1101 

1102 

1103def check_n_trial_events(data, **_): 

1104 """Check that the number events per trial is correct. 

1105 

1106 Within every trial interval there should be one of each trial event, except for 

1107 goCueTrigger_times which should only be defined for incorrect trials 

1108 

1109 Metric: 

1110 M = all(start < event < end) for all event times except errorCueTrigger_times where 

1111 start < error_trigger < end if not correct trial, else error_trigger == NaN 

1112 

1113 Criterion: 

1114 M == True 

1115 

1116 Units: 

1117 -none-, boolean 

1118 

1119 :param data: dict of trial data with keys ('intervals', 'stimOnTrigger_times', 

1120 'stimOffTrigger_times', 'stimOn_times', 'stimOff_times', 

1121 'stimFreezeTrigger_times', 'errorCueTrigger_times', 'itiIn_times', 

1122 'goCueTrigger_times', 'goCue_times', 'response_times', 'feedback_times') 

1123 """ 

1124 intervals = data['intervals'] 1apcbdfe

1125 correct = data['correct'] 1apcbdfe

1126 err_trig = data['errorCueTrigger_times'] 1apcbdfe

1127 

1128 # Exclude these fields; valve and errorCue times are the same as feedback_times and we must 

1129 # test errorCueTrigger_times separately 

1130 # stimFreeze_times fails often due to TTL flicker 

1131 exclude = ['camera_timestamps', 'errorCueTrigger_times', 'errorCue_times', 1apcbdfe

1132 'wheelMoves_peakVelocity_times', 'valveOpen_times', 'wheelMoves_peakAmplitude', 

1133 'wheelMoves_intervals', 'wheel_timestamps', 'stimFreeze_times'] 

1134 events = [k for k in data.keys() if k.endswith('_times') and k not in exclude] 1apcbdfe

1135 metric = np.zeros(data['intervals'].shape[0], dtype=bool) 1apcbdfe

1136 

1137 # For each trial interval check that one of each trial event occurred. For incorrect trials, 

1138 # check the error cue trigger occurred within the interval, otherwise check it is nan. 

1139 for i, (start, end) in enumerate(intervals): 1apcbdfe

1140 metric[i] = (all([start < data[k][i] < end for k in events]) and 1apcbdfe

1141 (np.isnan(err_trig[i]) if correct[i] else start < err_trig[i] < end)) 

1142 passed = metric.astype(bool) 1apcbdfe

1143 assert intervals.shape[0] == len(metric) == len(passed) 1apcbdfe

1144 return metric, passed 1apcbdfe

1145 

1146 

1147def check_trial_length(data, **_): 

1148 """ 

1149 Check open-loop duration positive and <= 1 minute. 

1150 

1151 Check that the time difference between the onset of the go cue sound 

1152 and the feedback (error sound or valve) is positive and smaller than 60.1 s. 

1153 

1154 Metric: 

1155 M = feedback_times - goCue_times 

1156 

1157 Criteria: 

1158 0 < M < 60.1 s 

1159 

1160 Units: 

1161 seconds [s] 

1162 

1163 :param data: dict of trial data with keys ('feedback_times', 'goCue_times', 'intervals') 

1164 """ 

1165 # NaN values are usually ignored so replace them with Inf so they fail the threshold 

1166 metric = np.nan_to_num(data['feedback_times'] - data['goCue_times'], nan=np.inf) 1aHcbdfe

1167 passed = (metric < 60.1) & (metric > 0) 1aHcbdfe

1168 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aHcbdfe

1169 return metric, passed 1aHcbdfe

1170 

1171 

1172# === Trigger-response delay checks === 

1173 

1174def check_goCue_delays(data, audio_output='harp', **_): 

1175 """ 

1176 Check the go cue tone occurs within 1ms of the intended time. 

1177 

1178 Check that the time difference between the go cue sound being triggered and 

1179 effectively played is smaller than 1ms. 

1180 

1181 Metric: 

1182 M = goCue_times - goCueTrigger_times 

1183 

1184 Criterion: 

1185 0 < M <= 0.0015 s 

1186 

1187 Units: 

1188 seconds [s] 

1189 

1190 :param data: dict of trial data with keys ('goCue_times', 'goCueTrigger_times', 'intervals'). 

1191 :param audio_output: audio output device name. 

1192 

1193 Notes 

1194 ----- 

1195 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th 

1196 percentile of delays over 500 training sessions using the Xonar soundcard. 

1197 """ 

1198 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.053 1agGcbdfe

1199 metric = np.nan_to_num(data['goCue_times'] - data['goCueTrigger_times'], nan=np.inf) 1agGcbdfe

1200 passed = (metric <= threshold) & (metric > 0) 1agGcbdfe

1201 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agGcbdfe

1202 return metric, passed 1agGcbdfe

1203 

1204 

1205def check_errorCue_delays(data, audio_output='harp', **_): 

1206 """ 

1207 Check the error tone occurs within 1.5ms of the intended time. 

1208 

1209 Check that the time difference between the error sound being triggered and 

1210 effectively played is smaller than 1ms. 

1211 

1212 Metric: 

1213 M = errorCue_times - errorCueTrigger_times 

1214 

1215 Criterion: 

1216 0 < M <= 0.0015 s 

1217 

1218 Units: 

1219 seconds [s] 

1220 

1221 :param data: dict of trial data with keys ('errorCue_times', 'errorCueTrigger_times', 

1222 'intervals', 'correct') 

1223 :param audio_output: audio output device name. 

1224 

1225 Notes 

1226 ----- 

1227 - For non-harp sound card the permissible delay is 0.062s. This was chosen by taking the 99.5th 

1228 percentile of delays over 500 training sessions using the Xonar soundcard. 

1229 """ 

1230 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.062 1aAcbdfe

1231 metric = np.nan_to_num(data['errorCue_times'] - data['errorCueTrigger_times'], nan=np.inf) 1aAcbdfe

1232 passed = ((metric <= threshold) & (metric > 0)).astype(float) 1aAcbdfe

1233 passed[data['correct']] = metric[data['correct']] = np.nan 1aAcbdfe

1234 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aAcbdfe

1235 return metric, passed 1aAcbdfe

1236 

1237 

1238def check_stimOn_delays(data, **_): 

1239 """ 

1240 Check the visual stimulus onset occurs within 150ms of the intended time. 

1241 

1242 Check that the time difference between the visual stimulus onset-command being triggered 

1243 and the stimulus effectively appearing on the screen is smaller than 150 ms. 

1244 

1245 Metric: 

1246 M = stimOn_times - stimOnTrigger_times 

1247 

1248 Criterion: 

1249 0 < M < 0.15 s 

1250 

1251 Units: 

1252 seconds [s] 

1253 

1254 :param data: dict of trial data with keys ('stimOn_times', 'stimOnTrigger_times', 

1255 'intervals') 

1256 """ 

1257 metric = np.nan_to_num(data['stimOn_times'] - data['stimOnTrigger_times'], nan=np.inf) 1agIcbdfe

1258 passed = (metric <= 0.15) & (metric > 0) 1agIcbdfe

1259 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agIcbdfe

1260 return metric, passed 1agIcbdfe

1261 

1262 

1263def check_stimOff_delays(data, **_): 

1264 """ 

1265 Check stimulus offset occurs within 150ms of the intended time. 

1266 

1267 Check that the time difference between the visual stimulus offset-command 

1268 being triggered and the visual stimulus effectively turning off on the screen 

1269 is smaller than 150 ms. 

1270 

1271 Metric: 

1272 M = stimOff_times - stimOffTrigger_times 

1273 

1274 Criterion: 

1275 0 < M < 0.15 s 

1276 

1277 Units: 

1278 seconds [s] 

1279 

1280 :param data: dict of trial data with keys ('stimOff_times', 'stimOffTrigger_times', 

1281 'intervals') 

1282 """ 

1283 metric = np.nan_to_num(data['stimOff_times'] - data['stimOffTrigger_times'], nan=np.inf) 1agJcbdfe

1284 passed = (metric <= 0.15) & (metric > 0) 1agJcbdfe

1285 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agJcbdfe

1286 return metric, passed 1agJcbdfe

1287 

1288 

1289def check_stimFreeze_delays(data, **_): 

1290 """Check the stimulus freezes within 150ms of the intended time. 

1291 

1292 Check that the time difference between the visual stimulus freeze-command 

1293 being triggered and the visual stimulus effectively freezing on the screen 

1294 is smaller than 150 ms. 

1295 

1296 Metric: 

1297 M = stimFreeze_times - stimFreezeTrigger_times 

1298 

1299 Criterion: 

1300 0 < M < 0.15 s 

1301 

1302 Units: 

1303 seconds [s] 

1304 

1305 :param data: dict of trial data with keys ('stimFreeze_times', 'stimFreezeTrigger_times', 

1306 'intervals') 

1307 """ 

1308 metric = np.nan_to_num(data['stimFreeze_times'] - data['stimFreezeTrigger_times'], nan=np.inf) 1aKcbdfe

1309 passed = (metric <= 0.15) & (metric > 0) 1aKcbdfe

1310 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aKcbdfe

1311 return metric, passed 1aKcbdfe

1312 

1313 

1314# === Data integrity checks === 

1315 

1316def check_reward_volumes(data, **_): 

1317 """Check that the reward volume is between 1.5 and 3 uL for correct trials, 0 for incorrect. 

1318 

1319 Metric: 

1320 M = reward volume 

1321 

1322 Criterion: 

1323 1.5 <= M <= 3 if correct else M == 0 

1324 

1325 Units: 

1326 uL 

1327 

1328 :param data: dict of trial data with keys ('rewardVolume', 'correct', 'intervals') 

1329 """ 

1330 metric = data['rewardVolume'] 1aycbdfe

1331 correct = data['correct'] 1aycbdfe

1332 passed = np.zeros_like(metric, dtype=bool) 1aycbdfe

1333 # Check correct trials within correct range 

1334 passed[correct] = (1.5 <= metric[correct]) & (metric[correct] <= 3.) 1aycbdfe

1335 # Check incorrect trials are 0 

1336 passed[~correct] = metric[~correct] == 0 1aycbdfe

1337 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aycbdfe

1338 return metric, passed 1aycbdfe

1339 

1340 

1341def check_reward_volume_set(data, **_): 

1342 """Check that there is only two reward volumes within a session, one of which is 0. 

1343 

1344 Metric: 

1345 M = set(rewardVolume) 

1346 

1347 Criterion: 

1348 (0 < len(M) <= 2) and 0 in M 

1349 

1350 :param data: dict of trial data with keys ('rewardVolume') 

1351 """ 

1352 metric = data['rewardVolume'] 1aLcbdfe

1353 passed = 0 < len(set(metric)) <= 2 and 0. in metric 1aLcbdfe

1354 return metric, passed 1aLcbdfe

1355 

1356 

1357def check_wheel_integrity(data, re_encoding='X1', enc_res=None, **_): 

1358 """ 

1359 Check wheel position sampled at the expected resolution. 

1360 

1361 Check that the difference between wheel position samples is close to the encoder resolution 

1362 and that the wheel timestamps strictly increase. 

1363 

1364 Metric: 

1365 M = (absolute difference of the positions < 1.5 * encoder resolution) 

1366 + 1 if (difference of timestamps <= 0) else 0 

1367 

1368 Criterion: 

1369 M ~= 0 

1370 

1371 Units: 

1372 arbitrary (radians, sometimes + 1) 

1373 

1374 :param data: dict of wheel data with keys ('wheel_timestamps', 'wheel_position') 

1375 :param re_encoding: the encoding of the wheel data, X1, X2 or X4 

1376 :param enc_res: the rotary encoder resolution (default 1024 ticks per revolution) 

1377 

1378 Notes 

1379 ----- 

1380 - At high velocities some samples are missed due to the scanning frequency of the DAQ. 

1381 This checks for more than 1 missing sample in a row (i.e. the difference between samples >= 2) 

1382 

1383 """ 

1384 if isinstance(re_encoding, str): 1avcbdfe

1385 re_encoding = int(re_encoding[-1]) 1avcbdfe

1386 # The expected difference between samples in the extracted units 

1387 resolution = 1 / (enc_res or ephys_fpga.WHEEL_TICKS 1avcbdfe

1388 ) * np.pi * 2 * ephys_fpga.WHEEL_RADIUS_CM / re_encoding 

1389 # We expect the difference of neighbouring positions to be close to the resolution 

1390 pos_check = np.abs(np.diff(data['wheel_position'])) 1avcbdfe

1391 # Timestamps should be strictly increasing 

1392 ts_check = np.diff(data['wheel_timestamps']) <= 0. 1avcbdfe

1393 metric = pos_check + ts_check.astype(float) # all values should be close to zero 1avcbdfe

1394 passed = metric < 1.5 * resolution 1avcbdfe

1395 return metric, passed 1avcbdfe

1396 

1397 

1398# === Pre-stimulus checks === 

1399def check_stimulus_move_before_goCue(data, photodiode=None, **_): 

1400 """ 

1401 Check there are no stimulus events before the go cue tone. 

1402 

1403 Check that there are no visual stimulus change(s) between the start of the trial and the 

1404 go cue sound onset, except for stim on. 

1405 

1406 Metric: 

1407 M = number of visual stimulus change events between trial start and goCue_times 

1408 

1409 Criterion: 

1410 M == 1 

1411 

1412 Units: 

1413 -none-, integer 

1414 

1415 Parameters 

1416 ---------- 

1417 data : dict 

1418 Trial data with keys ('goCue_times', 'intervals', 'choice'). 

1419 photodiode : dict 

1420 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel. 

1421 

1422 Returns 

1423 ------- 

1424 numpy.array 

1425 An array of metric values to threshold. 

1426 numpy.array 

1427 An array of boolean values, 1 per trial, where True means trial passes QC threshold. 

1428 

1429 Notes 

1430 ----- 

1431 - There should be exactly 1 stimulus change before goCue; stimulus onset. Even if the stimulus 

1432 contrast is 0, the sync square will still flip at stimulus onset, etc. 

1433 - If there are no goCue times (all are NaN), the status should be NOT_SET. 

1434 """ 

1435 if photodiode is None: 1acbdfe

1436 _log.warning('No photodiode TTL input in function call, returning None') 

1437 return None 

1438 photodiode_clean = ephys_fpga._clean_frame2ttl(photodiode) 1acbdfe

1439 s = photodiode_clean['times'] 1acbdfe

1440 s = s[~np.isnan(s)] # Remove NaNs 1acbdfe

1441 metric = np.array([]) 1acbdfe

1442 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1acbdfe

1443 metric = np.append(metric, np.count_nonzero(s[s > i] < c)) 1acbdfe

1444 

1445 passed = (metric == 1).astype(float) 1acbdfe

1446 passed[np.isnan(data['goCue_times'])] = np.nan 1acbdfe

1447 assert data['intervals'].shape[0] == len(metric) == len(passed) 1acbdfe

1448 return metric, passed 1acbdfe

1449 

1450 

1451def check_audio_pre_trial(data, audio=None, **_): 

1452 """ 

1453 Check no audio stimuli before the go cue. 

1454 

1455 Check there are no audio outputs between the start of the trial and the go cue sound onset - 20 ms. 

1456 

1457 Metric: 

1458 M = sum(start_times < audio TTL < (goCue_times - 20ms)) 

1459 

1460 Criterion: 

1461 M == 0 

1462 

1463 Units: 

1464 -none-, integer 

1465 

1466 Parameters 

1467 ---------- 

1468 data : dict 

1469 Trial data with keys ('goCue_times', 'intervals'). 

1470 audio : dict 

1471 The fronts from Bpod's BNC2 input FPGA audio sync channel. 

1472 

1473 Returns 

1474 ------- 

1475 numpy.array 

1476 An array of metric values to threshold. 

1477 numpy.array 

1478 An array of boolean values, 1 per trial, where True means trial passes QC threshold. 

1479 """ 

1480 if audio is None: 1awcbdfe

1481 _log.warning('No BNC2 input in function call, retuning None') 

1482 return None, None 

1483 s = audio['times'][~np.isnan(audio['times'])] # Audio TTLs with NaNs removed 1awcbdfe

1484 metric = np.array([], dtype=np.int8) 1awcbdfe

1485 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1awcbdfe

1486 metric = np.append(metric, sum(s[s > i] < (c - 0.02))) 1awcbdfe

1487 passed = metric == 0 1awcbdfe

1488 assert data['intervals'].shape[0] == len(metric) == len(passed) 1awcbdfe

1489 return metric, passed 1awcbdfe