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

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 (empty)

50import sys (empty)

51from packaging import version (empty)

52from pathlib import PurePosixPath (empty)

53from datetime import datetime, timedelta (empty)

54from inspect import getmembers, isfunction (empty)

55from functools import reduce (empty)

56from collections.abc import Sized (empty)

57 

58import numpy as np (empty)

59from scipy.stats import chisquare (empty)

60 

61from brainbox.behavior.wheel import cm_to_rad, traces_by_trial (empty)

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

63from one.alf import spec (empty)

64from . import base (empty)

65 

66_log = logging.getLogger(__name__) (empty)

67 

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

69ITI_DELAY_SECS = .5 (empty)

70FEEDBACK_NOGO_DELAY_SECS = 2 (empty)

71 

72BWM_CRITERIA = { (empty)

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): (empty)

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: 11 ctx1axmlngcbdei

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

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

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()) 10 ctx1amlngcbdei

119 return session_outcome, outcomes 10 ctx1amlngcbdei

120 

121 

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

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} 2 ctx1hc

144 # Ensure dataset stems are unique 

145 assert len(stem2id) == len(registered_datasets), 'ambiguous dataset names' 2 ctx1hc

146 

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

148 *_, outcomes = qc.compute_session_status() 2 ctx1hc

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(): 2 ctx1hc

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

152 if isinstance(outcome, dict): 2 ctx1hc

153 extended_qc = outcome 2 ctx1hc

154 outcome = qc.overall_outcome(outcome.values()) 2 ctx1hc

155 else: 

156 extended_qc = {} 2 ctx1hc

157 # check if dataset was registered to Alyx 

158 if not (did := stem2id.get(name)): 2 ctx1hc

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

160 continue 1 ctx1h

161 # update the dataset QC value on Alyx 

162 if outcome > spec.QC.NOT_SET or override: 2 ctx1hc

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

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

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

166 if extended_qc: 2 ctx1hc

167 dset_qc.update_extended_qc(extended_qc) 2 ctx1hc

168 return registered_datasets 2 ctx1hc

169 

170 

171class TaskQC(base.QC): (empty)

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

173 

174 criteria = BWM_CRITERIA (empty)

175 

176 extractor = None (empty)

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

178 

179 @staticmethod (empty)

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

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 {}} 11 ctx1axmlngcbdei

198 MAX_BOUND, MIN_BOUND = (1, 0) 11 ctx1axmlngcbdei

199 if qc_value is None or np.isnan(qc_value): 11 ctx1axmlngcbdei

200 return spec.QC.NOT_SET 8 ctx1amlngbei

201 elif (qc_value > MAX_BOUND) or (qc_value < MIN_BOUND): 11 ctx1axmlngcbdei

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

203 for crit in filter(None, sorted(spec.QC)): 10 ctx1amlngcbdei

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

205 return crit 10 ctx1amlngcbdei

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

207 return spec.QC.NOT_SET 8 ctx1algcbdei

208 

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

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) 5 ctx1ahcbd

225 

226 # Data 

227 self.extractor = None 5 ctx1ahcbd

228 

229 # Metrics and passed trials 

230 self.metrics = None 5 ctx1ahcbd

231 self.passed = None 5 ctx1ahcbd

232 

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

234 self.namespace = namespace 5 ctx1ahcbd

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

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

237 

238 def compute(self, **kwargs): (empty)

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 6 ctx1acbdfe

250 ns = kwargs.pop('namespace', self.namespace) 6 ctx1acbdfe

251 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 6 ctx1acbdfe

252 if version.parse(ver) >= version.parse('8.0.0'): 6 ctx1acbdfe

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} 6 ctx1acbdfe

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

258 

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

260 self.get_bpodqc_metrics_frame( 6 ctx1acbdfe

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): (empty)

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): 8 ctx1athcbdfe

281 return isfunction(x) and x.__name__.startswith('check_') 8 ctx1athcbdfe

282 

283 return dict(getmembers(sys.modules[__name__], is_metric)) 8 ctx1athcbdfe

284 

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

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() 6 ctx1acbdfe

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

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()} 6 ctx1acbdfe

321 

322 # Split metrics and passed frames 

323 self.metrics = {} 6 ctx1acbdfe

324 self.passed = {} 6 ctx1acbdfe

325 for k in qc_metrics_map: 6 ctx1acbdfe

326 self.metrics[k], self.passed[k] = qc_metrics_map[k] 6 ctx1acbdfe

327 

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

329 n_trials = data['intervals'].shape[0] 6 ctx1acbdfe

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] 6 ctx1acbdfe

332 name = prefix + 'passed_trial_checks' 6 ctx1acbdfe

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

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

335 

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

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: 5 ctx1acbdi

355 self.compute(**kwargs) 4 ctx1acbd

356 outcome, results, _ = self.compute_session_status() 5 ctx1acbdi

357 if update: 5 ctx1acbdi

358 self.update_extended_qc(results) 2 ctx1ci

359 self.update(outcome, kwargs.get('namespace', self.namespace)) 2 ctx1ci

360 return outcome, results 5 ctx1acbdi

361 

362 def compute_session_status(self): (empty)

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: 7 ctx1agcbdei

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

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) 7 ctx1agcbdei

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

380 session_outcome, outcomes = compute_session_status_from_dict(results, self.criteria) 7 ctx1agcbdei

381 return session_outcome, results, outcomes 7 ctx1agcbdei

382 

383 @staticmethod (empty)

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

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 = { 3 ctx1thc

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') 3 ctx1thc

410 trials_table_outcomes['rewardVolume']: TaskQC.overall_outcome( 3 ctx1thc

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

412 ) 

413 dataset_outcomes = { 3 ctx1thc

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 3 ctx1thc

418 

419 

420class HabituationQC(TaskQC): (empty)

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

422 

423 def compute(self, **kwargs): (empty)

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 2 ctx1gb

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

431 

432 # Initialize checks 

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

434 data = self.extractor.data 2 ctx1gb

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

436 metrics = {} 2 ctx1gb

437 passed = {} 2 ctx1gb

438 

439 # Modify criteria based on version 

440 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 2 ctx1gb

441 is_v8 = version.parse(ver) >= version.parse('8.0.0') 2 ctx1gb

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

443 

444 # Check all reward volumes == 3.0ul 

445 check = prefix + 'reward_volumes' 2 ctx1gb

446 metrics[check] = data['rewardVolume'] 2 ctx1gb

447 passed[check] = metrics[check] == 3.0 2 ctx1gb

448 

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

450 check = prefix + 'habituation_time' 2 ctx1gb

451 if not self.one or not self.session_path: 2 ctx1gb

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] 2 ctx1gb

456 # compute from the date specified 

457 date_minus_week = ( 2 ctx1gb

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, 2 ctx1gb

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): 2 ctx1gb

465 sessions = sessions[1:] 

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

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 2 ctx1gb

476 passed[check] = np.diff(metric) >= 12 2 ctx1gb

477 

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

479 check = prefix + 'trial_event_sequence' 2 ctx1gb

480 nans = ( 2 ctx1gb

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) 2 ctx1gb

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

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

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

491 

492 metrics[check] = a & b & c & d & ~nans 2 ctx1gb

493 passed[check] = metrics[check].astype(float) 2 ctx1gb

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' 2 ctx1gb

498 metric = np.nan_to_num(data['stimCenter_times'] - data['stimCenterTrigger_times'], 2 ctx1gb

499 nan=np.inf) 

500 passed[check] = (metric <= 0.15) & (metric > 0) 2 ctx1gb

501 metrics[check] = metric 2 ctx1gb

502 

503 # Phase check 

504 check = prefix + 'phase' 2 ctx1gb

505 metric = data['phase'] 2 ctx1gb

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

507 metrics[check] = metric 2 ctx1gb

508 

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

510 check = prefix + 'phase_distribution' 2 ctx1gb

511 metric, _ = np.histogram(data['phase']) 2 ctx1gb

512 _, p = chisquare(metric) 2 ctx1gb

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

514 metrics[check] = metric 2 ctx1gb

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' 2 ctx1gb

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

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

521 passed[check] = np.abs(metric) <= 0.1 2 ctx1gb

522 passed[check][-1] = np.nan 2 ctx1gb

523 metrics[check] = metric 2 ctx1gb

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' 2 ctx1gb

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 2 ctx1gb

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

533 passed[check] = np.abs(metric) < threshold 2 ctx1gb

534 metrics[check] = metric 2 ctx1gb

535 

536 # Checks common to training QC 

537 checks = [check_goCue_delays, check_stimOn_delays, check_stimOff_delays] 2 ctx1gb

538 for fcn in checks: 2 ctx1gb

539 check = prefix + fcn.__name__[6:] 2 ctx1gb

540 metrics[check], passed[check] = fcn(data, audio_output=audio_output) 2 ctx1gb

541 

542 self.metrics, self.passed = (metrics, passed) 2 ctx1gb

543 

544 

545# SINGLE METRICS 

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

547 

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

549 

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

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 7 ctx1aBcbdfe

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

578 passed = (metric < threshold) & (metric > 0) 7 ctx1aBcbdfe

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

580 return metric, passed 7 ctx1aBcbdfe

581 

582 

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

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 7 ctx1aCcbdfe

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

609 passed = (metric < threshold) & (metric > 0) 7 ctx1aCcbdfe

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

611 return metric, passed 7 ctx1aCcbdfe

612 

613 

614def check_response_stimFreeze_delays(data, **_): (empty)

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) 7 ctx1aDcbdfe

636 # Test for valid values 

637 passed = ((metric < 0.1) & (metric > 0)).astype(float) 7 ctx1aDcbdfe

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 7 ctx1aDcbdfe

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

642 return metric, passed 7 ctx1aDcbdfe

643 

644 

645def check_stimOff_itiIn_delays(data, **_): (empty)

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) 7 ctx1aEcbdfe

662 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 7 ctx1aEcbdfe

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 7 ctx1aEcbdfe

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

667 return metric, passed 7 ctx1aEcbdfe

668 

669 

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

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) 7 ctx1aucbdfe

704 passed = metric.copy() 7 ctx1aucbdfe

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

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( 7 ctx1aucbdfe

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 7 ctx1aucbdfe

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

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

715 return metric, passed 7 ctx1aucbdfe

716 

717 

718def check_positive_feedback_stimOff_delays(data, **_): (empty)

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) 7 ctx1aFcbdfe

739 passed = (np.abs(metric) < 0.15).astype(float) 7 ctx1aFcbdfe

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

741 metric[~data['correct']] = passed[~data['correct']] = np.nan 7 ctx1aFcbdfe

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

743 return metric, passed 7 ctx1aFcbdfe

744 

745 

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

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) 7 ctx1azcbdfe

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 7 ctx1azcbdfe

767 # Apply criteria 

768 passed = (np.abs(metric) < 0.15).astype(float) 7 ctx1azcbdfe

769 # Remove none negative feedback trials 

770 metric[data['correct']] = passed[data['correct']] = np.nan 7 ctx1azcbdfe

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

772 return metric, passed 7 ctx1azcbdfe

773 

774 

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

776 

777def check_wheel_move_before_feedback(data, **_): (empty)