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

1"""Behaviour QC. 

2 

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

4 

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

6instantiated directly. 

7 

8Examples 

9-------- 

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

11 

12>>> from ibllib.qc.task_metrics import TaskQC 

13>>> TaskQC('path/to/session').run(update=True) 

14 

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

16 

17>>> from ibllib.qc.task_metrics import TaskQC 

18>>> qc = TaskQC(eid) 

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

20 

21Inspecting individual test outcomes 

22 

23>>> from ibllib.qc.task_metrics import TaskQC 

24>>> qc = TaskQC(eid) 

25>>> outcome, results, outcomes = qc.compute().compute_session_status() 

26 

27Running bpod QC on ephys session 

28 

29>>> from ibllib.qc.task_metrics import TaskQC 

30>>> qc = TaskQC(eid) 

31>>> qc.load_data(bpod_only=True) # Extract without FPGA 

32>>> bpod_qc = qc.run() 

33 

34Running bpod QC only, from training rig PC 

35 

36>>> from ibllib.qc.task_metrics import TaskQC 

37>>> from ibllib.qc.qcplots import plot_results 

38>>> session_path = r'/home/nico/Downloads/FlatIron/mrsicflogellab/Subjects/SWC_023/2020-02-14/001' 

39>>> qc = TaskQC(session_path) 

40>>> qc.load_data(bpod_only=True, download_data=False) # Extract without FPGA 

41>>> qc.run() 

42>>> plot_results(qc, save_path=session_path) 

43 

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

45 

46>>> from ibllib.qc.task_metrics import TaskQC 

47>>> from ibllib.qc.qcplots import plot_results 

48>>> session_path = r'/home/nico/Downloads/FlatIron/mrsicflogellab/Subjects/SWC_023/2020-02-14/001' 

49>>> qc = TaskQC(session_path) 

50>>> qc.run() 

51>>> plot_results(qc, save_path=session_path) 

52""" 

53import logging (empty)

54import sys (empty)

55from packaging import version (empty)

56from pathlib import Path, PurePosixPath (empty)

57from datetime import datetime, timedelta (empty)

58from inspect import getmembers, isfunction (empty)

59from functools import reduce (empty)

60from collections.abc import Sized (empty)

61 

62import numpy as np (empty)

63from scipy.stats import chisquare (empty)

64 

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

66from ibllib.qc.task_extractors import TaskQCExtractor (empty)

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

68from one.alf import spec (empty)

69from . import base (empty)

70 

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

72 

73 

74BWM_CRITERIA = { (empty)

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

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

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

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

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

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

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

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

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

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

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

86 '_task_stimOff_delays': {'PASS': 0.99, 'WARNING': 0}, 

87 '_task_stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0}, 

88 '_task_iti_delays': {'NOT_SET': 0}, 

89 '_task_passed_trial_checks': {'NOT_SET': 0} 

90} 

91 

92 

93def compute_session_status_from_dict(results, criteria=None): (empty)

94 """ 

95 Compute overall task QC value from QC check results. 

96 

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

98 in a single value. 

99 

100 Parameters 

101 ---------- 

102 results : dict 

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

104 criteria : dict 

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

106 

107 Returns 

108 ------- 

109 one.alf.spec.QC 

110 Overall session QC outcome. 

111 dict 

112 A map of QC tests and their outcomes. 

113 """ 

114 if not criteria: 11 ctx1ynmohcadfeb

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

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

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

118 

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

120 session_outcome = base.QC.overall_outcome(outcomes.values()) 10 ctx1nmohcadfeb

121 return session_outcome, outcomes 10 ctx1nmohcadfeb

122 

123 

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

125 """ 

126 Update QC values for individual datasets. 

127 

128 Parameters 

129 ---------- 

130 qc : ibllib.qc.task_metrics.TaskQC 

131 A TaskQC object that has been run. 

132 registered_datasets : list of dict 

133 A list of Alyx dataset records. 

134 one : one.api.OneAlyx 

135 An online instance of ONE. 

136 override : bool 

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

138 

139 Returns 

140 ------- 

141 list of dict 

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

143 """ 

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

145 stem2id = {PurePosixPath(dset['name']).stem: dset.get('id') for dset in registered_datasets} 3 ctx1jcb

146 # Ensure dataset stems are unique 

147 assert len(stem2id) == len(registered_datasets), 'ambiguous dataset names' 3 ctx1jcb

148 

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

150 *_, outcomes = qc.compute_session_status() 3 ctx1jcb

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

152 for name, outcome in qc.compute_dataset_qc_status(outcomes).items(): 3 ctx1jcb

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

154 if isinstance(outcome, dict): 3 ctx1jcb

155 extended_qc = outcome 3 ctx1jcb

156 outcome = qc.overall_outcome(outcome.values()) 3 ctx1jcb

157 else: 

158 extended_qc = {} 3 ctx1jcb

159 # check if dataset was registered to Alyx 

160 if not (did := stem2id.get(name)): 3 ctx1jcb

161 _log.debug('dataset %s not registered, skipping', name) 3 ctx1jcb

162 continue 3 ctx1jcb

163 # update the dataset QC value on Alyx 

164 if outcome > spec.QC.NOT_SET or override: 3 ctx1jcb

165 dset_qc = base.QC(did, one=one, log=_log, endpoint='datasets') 3 ctx1jcb

166 dset = next(x for x in registered_datasets if did == x.get('id')) 3 ctx1jcb

167 dset['qc'] = dset_qc.update(outcome, namespace='', override=override).name 3 ctx1jcb

168 if extended_qc: 3 ctx1jcb

169 dset_qc.update_extended_qc(extended_qc) 3 ctx1jcb

170 return registered_datasets 3 ctx1jcb

171 

172 

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

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

175 

176 criteria = BWM_CRITERIA (empty)

177 

178 extractor = None (empty)

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

180 

181 @staticmethod (empty)

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

183 """ 

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

185 

186 Parameters 

187 ---------- 

188 qc_value : float 

189 Proportion of passing qcs, between 0 and 1. 

190 thresholds : dict 

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

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

193 

194 Returns 

195 ------- 

196 one.alf.spec.QC 

197 The outcome. 

198 """ 

199 thresholds = {spec.QC.validate(k): v for k, v in thresholds.items() or {}} 11 ctx1ynmohcadfeb

200 MAX_BOUND, MIN_BOUND = (1, 0) 11 ctx1ynmohcadfeb

201 if qc_value is None or np.isnan(qc_value): 11 ctx1ynmohcadfeb

202 return spec.QC.NOT_SET 8 ctx1nmohafeb

203 elif (qc_value > MAX_BOUND) or (qc_value < MIN_BOUND): 11 ctx1ynmohcadfeb

204 raise ValueError('Values out of bound') 1 ctx1y

205 for crit in filter(None, sorted(spec.QC)): 10 ctx1nmohcadfeb

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

207 return crit 10 ctx1nmohcadfeb

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

209 return spec.QC.NOT_SET 8 ctx1mhcadfeb

210 

211 def __init__(self, session_path_or_eid, **kwargs): (empty)

212 """ 

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

214 

215 :param session_path_or_eid: A session eid or path 

216 :param log: A logging.Logger instance, if None the 'ibllib' logger is used 

217 :param one: An ONE instance for fetching and setting the QC on Alyx 

218 """ 

219 # When an eid is provided, we will download the required data by default (if necessary) 

220 self.download_data = not spec.is_session_path(Path(session_path_or_eid)) 6 ctx1ijcadb

221 super().__init__(session_path_or_eid, **kwargs) 6 ctx1ijcadb

222 

223 # Data 

224 self.extractor = None 6 ctx1ijcadb

225 

226 # Metrics and passed trials 

227 self.metrics = None 6 ctx1ijcadb

228 self.passed = None 6 ctx1ijcadb

229 

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

231 self.criteria = BWM_CRITERIA.copy() 6 ctx1ijcadb

232 

233 def load_data(self, bpod_only=False, download_data=True): (empty)

234 """Extract the data from raw data files. 

235 

236 Extracts all the required task data from the raw data files. 

237 

238 Parameters 

239 ---------- 

240 bpod_only : bool 

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

242 download_data : bool 

243 If True, any missing raw data is downloaded via ONE. By default data are not downloaded 

244 if a session path was provided to the constructor. 

245 """ 

246 self.extractor = TaskQCExtractor( (empty)

247 self.session_path, one=self.one, download_data=download_data, bpod_only=bpod_only) 

248 

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

250 """Compute and store the QC metrics. 

251 

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

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

254 

255 Parameters 

256 ---------- 

257 bpod_only : bool 

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

259 download_data : bool 

260 If True, any missing raw data is downloaded via ONE. By default data are not downloaded 

261 if a session path was provided to the constructor. 

262 """ 

263 if self.extractor is None: 7 ctx1cadgfeb

264 kwargs['download_data'] = kwargs.pop('download_data', self.download_data) 

265 self.load_data(**kwargs) 

266 

267 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 7 ctx1cadgfeb

268 if version.parse(ver) >= version.parse('8.0.0'): 7 ctx1cadgfeb

269 self.criteria['_task_iti_delays'] = {'PASS': 0.99, 'WARNING': 0} 

270 self.criteria['_task_passed_trial_checks'] = {'PASS': 0.7, 'WARNING': 0} 

271 else: 

272 self.criteria['_task_iti_delays'] = {'NOT_SET': 0} 7 ctx1cadgfeb

273 self.criteria['_task_passed_trial_checks'] = {'NOT_SET': 0} 7 ctx1cadgfeb

274 

275 self.log.info(f'Session {self.session_path}: Running QC on behavior data...') 7 ctx1cadgfeb

276 self.get_bpodqc_metrics_frame( 7 ctx1cadgfeb

277 self.extractor.data, 

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

279 photodiode=self.extractor.frame_ttls, 

280 audio=self.extractor.audio_ttls, 

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

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

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

284 ) 

285 

286 def _get_checks(self): (empty)

287 """ 

288 Find all methods that begin with 'check_'. 

289 

290 Returns 

291 ------- 

292 Dict[str, function] 

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

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

295 """ 

296 def is_metric(x): 9 ctx1ujcadgfeb

297 return isfunction(x) and x.__name__.startswith('check_') 9 ctx1ujcadgfeb

298 

299 return dict(getmembers(sys.modules[__name__], is_metric)) 9 ctx1ujcadgfeb

300 

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

302 """ 

303 Evaluate task QC metrics. 

304 

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

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

307 

308 Parameters 

309 ---------- 

310 data : dict 

311 The extracted task data. 

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

313 The encoding configuration of the rotary encoder. 

314 enc_res : int 

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

316 wheel_gain : float 

317 The STIM_GAIN task parameter. 

318 photodiode : dict 

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

320 audio : dict 

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

322 min_qt : float 

323 The QUIESCENT_PERIOD task parameter. 

324 

325 Returns 

326 ------- 

327 dict 

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

329 dict 

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

331 """ 

332 # Find all methods that begin with 'check_' 

333 checks = self._get_checks() 7 ctx1cadgfeb

334 prefix = '_task_' # Extended QC fields will start with this 7 ctx1cadgfeb

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

336 qc_metrics_map = {prefix + k[6:]: fn(data, **kwargs) for k, fn in checks.items()} 7 ctx1cadgfeb

337 

338 # Split metrics and passed frames 

339 self.metrics = {} 7 ctx1cadgfeb

340 self.passed = {} 7 ctx1cadgfeb

341 for k in qc_metrics_map: 7 ctx1cadgfeb

342 self.metrics[k], self.passed[k] = qc_metrics_map[k] 7 ctx1cadgfeb

343 

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

345 n_trials = data['intervals'].shape[0] 7 ctx1cadgfeb

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

347 trial_level_passed = [m for m in self.passed.values() if isinstance(m, Sized) and len(m) == n_trials] 7 ctx1cadgfeb

348 name = prefix + 'passed_trial_checks' 7 ctx1cadgfeb

349 self.metrics[name] = reduce(np.logical_and, trial_level_passed or (None, None)) 7 ctx1cadgfeb

350 self.passed[name] = self.metrics[name].astype(float) if trial_level_passed else None 7 ctx1cadgfeb

351 

352 def run(self, update=False, namespace='task', **kwargs): (empty)

353 """ 

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

355 

356 Parameters 

357 ---------- 

358 update : bool 

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

360 namespace : str 

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

362 bpod_only : bool 

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

364 download_data : bool 

365 If True, any missing raw data is downloaded via ONE. By default data are not downloaded 

366 if a session path was provided to the constructor. 

367 

368 Returns 

369 ------- 

370 str 

371 Overall task QC outcome. 

372 dict 

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

374 """ 

375 if self.metrics is None: 5 ctx1cadeb

376 self.compute(**kwargs) 5 ctx1cadeb

377 outcome, results, _ = self.compute_session_status() 5 ctx1cadeb

378 if update: 5 ctx1cadeb

379 self.update_extended_qc(results) 3 ctx1ceb

380 self.update(outcome, namespace) 3 ctx1ceb

381 return outcome, results 5 ctx1cadeb

382 

383 def compute_session_status(self): (empty)

384 """ 

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

386 

387 Returns 

388 ------- 

389 str 

390 Overall session QC outcome. 

391 dict 

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

393 dict 

394 A map of QC tests and their outcomes. 

395 """ 

396 if self.passed is None: 7 ctx1hcadfeb

397 raise AttributeError('passed is None; compute QC first') 1 ctx1f

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

399 results = {k: None if v is None or np.isnan(v).all() else np.nanmean(v) 7 ctx1hcadfeb

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

401 session_outcome, outcomes = compute_session_status_from_dict(results, self.criteria) 7 ctx1hcadfeb

402 return session_outcome, results, outcomes 7 ctx1hcadfeb

403 

404 @staticmethod (empty)

405 def compute_dataset_qc_status(outcomes): (empty)

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

407 

408 Parameters 

409 ---------- 

410 outcomes : dict 

411 Map of checks and their individual outcomes. 

412 

413 Returns 

414 ------- 

415 dict 

416 Map of dataset names and their outcome. 

417 """ 

418 trials_table_outcomes = { 4 ctx1ujcb

419 'intervals': outcomes.get('_task_iti_delays', spec.QC.NOT_SET), 

420 'goCue_times': outcomes.get('_task_goCue_delays', spec.QC.NOT_SET), 

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

422 'stimOn_times': outcomes.get('_task_stimOn_delays', spec.QC.NOT_SET), 

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

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

425 'feedback_times': outcomes.get('_task_errorCue_delays', spec.QC.NOT_SET), 

426 'firstMovement_times': spec.QC.NOT_SET 

427 } 

428 reward_checks = ('_task_reward_volumes', '_task_reward_volume_set') 4 ctx1ujcb

429 trials_table_outcomes['rewardVolume']: TaskQC.overall_outcome( 4 ctx1ujcb

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

431 ) 

432 dataset_outcomes = { 4 ctx1ujcb

433 '_ibl_trials.stimOff_times': outcomes.get('_task_stimOff_delays', spec.QC.NOT_SET), 

434 '_ibl_trials.table': trials_table_outcomes, 

435 } 

436 return dataset_outcomes 4 ctx1ujcb

437 

438 

439class HabituationQC(TaskQC): (empty)

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

441 

442 def compute(self, download_data=None, **kwargs): (empty)

443 """Compute and store the QC metrics. 

444 

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

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

447 """ 

448 if self.extractor is None: 2 ctx1ha

449 # If download_data is None, decide based on whether eid or session path was provided 

450 ensure_data = self.download_data if download_data is None else download_data 

451 self.load_data(download_data=ensure_data, **kwargs) 

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

453 

454 # Initialize checks 

455 prefix = '_task_' 2 ctx1ha

456 data = self.extractor.data 2 ctx1ha

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

458 metrics = {} 2 ctx1ha

459 passed = {} 2 ctx1ha

460 

461 # Modify criteria based on version 

462 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 2 ctx1ha

463 is_v8 = version.parse(ver) >= version.parse('8.0.0') 2 ctx1ha

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

465 

466 # Check all reward volumes == 3.0ul 

467 check = prefix + 'reward_volumes' 2 ctx1ha

468 metrics[check] = data['rewardVolume'] 2 ctx1ha

469 passed[check] = metrics[check] == 3.0 2 ctx1ha

470 

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

472 check = prefix + 'habituation_time' 2 ctx1ha

473 if not self.one or not self.session_path: 2 ctx1ha

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

475 metrics[check] = passed[check] = None 

476 else: 

477 subject, session_date = self.session_path.parts[-3:-1] 2 ctx1ha

478 # compute from the date specified 

479 date_minus_week = ( 2 ctx1ha

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

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

482 sessions = self.one.alyx.rest('sessions', 'list', subject=subject, 2 ctx1ha

483 date_range=[date_minus_week, session_date], 

484 task_protocol='habituation') 

485 # Remove the current session if already registered 

486 if sessions and sessions[0]['start_time'].startswith(session_date): 2 ctx1ha

487 sessions = sessions[1:] 

488 metric = ([0, data['intervals'][-1, 1] - data['intervals'][0, 0]] + 2 ctx1ha

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

490 datetime.fromisoformat(x['start_time'])).total_seconds() / 60 

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

492 

493 # The duration from raw trial data 

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

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

496 # duration))).total_seconds() / 60 

497 metrics[check] = np.array(metric) 2 ctx1ha

498 passed[check] = np.diff(metric) >= 12 2 ctx1ha

499 

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

501 check = prefix + 'trial_event_sequence' 2 ctx1ha

502 nans = ( 2 ctx1ha

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

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

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

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

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

508 ) 

509 a = np.less(data['intervals'][:, 0], data['stimOn_times'], where=~nans) 2 ctx1ha

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

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

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

513 

514 metrics[check] = a & b & c & d & ~nans 2 ctx1ha

515 passed[check] = metrics[check].astype(float) 2 ctx1ha

516 

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

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

519 check = prefix + 'stimCenter_delays' 2 ctx1ha

520 metric = np.nan_to_num(data['stimCenter_times'] - data['stimCenterTrigger_times'], 2 ctx1ha

521 nan=np.inf) 

522 passed[check] = (metric <= 0.15) & (metric > 0) 2 ctx1ha

523 metrics[check] = metric 2 ctx1ha

524 

525 # Phase check 

526 check = prefix + 'phase' 2 ctx1ha

527 metric = data['phase'] 2 ctx1ha

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

529 metrics[check] = metric 2 ctx1ha

530 

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

532 check = prefix + 'phase_distribution' 2 ctx1ha

533 metric, _ = np.histogram(data['phase']) 2 ctx1ha

534 _, p = chisquare(metric) 2 ctx1ha

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

536 metrics[check] = metric 2 ctx1ha

537 

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

539 # 1s +/- 10%. 

540 check = prefix + 'iti_delays' 2 ctx1ha

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

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

543 passed[check] = np.abs(metric) <= 0.1 2 ctx1ha

544 passed[check][-1] = np.NaN 2 ctx1ha

545 metrics[check] = metric 2 ctx1ha

546 

547 # Checks common to training QC 

548 checks = [check_goCue_delays, check_stimOn_goCue_delays, 2 ctx1ha

549 check_stimOn_delays, check_stimOff_delays] 

550 for fcn in checks: 2 ctx1ha

551 check = prefix + fcn.__name__[6:] 2 ctx1ha

552 metrics[check], passed[check] = fcn(data, audio_output=audio_output) 2 ctx1ha

553 

554 self.metrics, self.passed = (metrics, passed) 2 ctx1ha

555 

556 

557# SINGLE METRICS 

558# ---------------------------------------------------------------------------- # 

559 

560# === Delays between events checks === 

561 

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

563 """ 

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

565 

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

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

568 

569 Metric: 

570 M = stimOn_times - goCue_times 

571 

572 Criteria: 

573 0 < M < 0.010 s 

574 

575 Units: 

576 seconds [s] 

577 

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

579 :param audio_output: audio output device name. 

580 

581 Notes 

582 ----- 

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

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

585 """ 

586 # Calculate the difference between stimOn and goCue times. 

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

588 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 9 ctx1hBcadgfeb

589 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 9 ctx1hBcadgfeb

590 passed = (metric < threshold) & (metric > 0) 9 ctx1hBcadgfeb

591 assert data['intervals'].shape[0] == len(metric) == len(passed) 9 ctx1hBcadgfeb

592 return metric, passed 9 ctx1hBcadgfeb

593 

594 

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

596 """ 

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

598 

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

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

601 

602 Metric: 

603 M = feedback_time - response_time 

604 

605 Criterion: 

606 0 < M < 0.010 s 

607 

608 Units: 

609 seconds [s] 

610 

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

612 :param audio_output: audio output device name. 

613 

614 Notes 

615 ----- 

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

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

618 """ 

619 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 8 ctx1Ccadgfeb

620 metric = np.nan_to_num(data['feedback_times'] - data['response_times'], nan=np.inf) 8 ctx1Ccadgfeb

621 passed = (metric < threshold) & (metric > 0) 8 ctx1Ccadgfeb

622 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1Ccadgfeb

623 return metric, passed 8 ctx1Ccadgfeb

624 

625 

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

627 """ 

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

629 

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

631 response is positive and less than 100ms. 

632 

633 Metric: 

634 M = (stimFreeze_times - response_times) 

635 

636 Criterion: 

637 0 < M < 0.100 s 

638 

639 Units: 

640 seconds [s] 

641 

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

643 'choice') 

644 """ 

645 # Calculate the difference between stimOn and goCue times. 

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

647 metric = np.nan_to_num(data['stimFreeze_times'] - data['response_times'], nan=np.inf) 8 ctx1Dcadgfeb

648 # Test for valid values 

649 passed = ((metric < 0.1) & (metric > 0)).astype(float) 8 ctx1Dcadgfeb

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

651 # These values are ignored in calculation of proportion passed 

652 passed[data['choice'] == 0] = np.nan 8 ctx1Dcadgfeb

653 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1Dcadgfeb

654 return metric, passed 8 ctx1Dcadgfeb

655 

656 

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

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

659 

660 Metric: 

661 M = itiIn_times - stimOff_times 

662 

663 Criterion: 

664 0 < M < 0.010 s 

665 

666 Units: 

667 seconds [s] 

668 

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

670 'choice') 

671 """ 

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

673 metric = np.nan_to_num(data['itiIn_times'] - data['stimOff_times'], nan=np.inf) 8 ctx1Ecadgfeb

674 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 8 ctx1Ecadgfeb

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

676 # NaN values are ignored in calculation of proportion passed 

677 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan 8 ctx1Ecadgfeb

678 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1Ecadgfeb

679 return metric, passed 8 ctx1Ecadgfeb

680 

681 

682def check_iti_delays(data, subtract_pauses=False, **_): (empty)

683 """ 

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

685 

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

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

688 

689 Metric: 

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

691 

692 Criterion: 

693 |M| < 0.1 

694 

695 Units: 

696 seconds [s] 

697 

698 Parameters 

699 ---------- 

700 data : dict 

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

702 subtract_pauses: bool 

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

704 the experimenter paused the task may fail this check. 

705 

706 Returns 

707 ------- 

708 numpy.array 

709 An array of metric values to threshold. 

710 numpy.array 

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

712 """ 

713 # Initialize array the length of completed trials 

714 ITI = 1. 8 ctx1vcadgfeb

715 metric = np.full(data['intervals'].shape[0], np.nan) 8 ctx1vcadgfeb

716 passed = metric.copy() 8 ctx1vcadgfeb

717 pauses = (data['pause_duration'] if subtract_pauses else np.zeros_like(metric))[:-1] 8 ctx1vcadgfeb

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

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

720 metric[:-1] = \ 8 ctx1vcadgfeb

721 np.nan_to_num(data['intervals'][1:, 0] - data['stimOff_times'][:-1] - ITI - pauses, nan=np.inf) 

722 passed[:-1] = np.abs(metric[:-1]) < (ITI / 10) # Last trial is not counted 8 ctx1vcadgfeb

723 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1vcadgfeb

724 return metric, passed 8 ctx1vcadgfeb

725 

726 

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

728 """ 

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

730 

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

732 is 1 ± 0.150 seconds. 

733 

734 Metric: 

735 M = stimOff_times - feedback_times - 1s 

736 

737 Criterion: 

738 |M| < 0.150 s 

739 

740 Units: 

741 seconds [s] 

742 

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

744 'correct') 

745 """ 

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

747 metric = np.nan_to_num(data['stimOff_times'] - data['feedback_times'] - 1, nan=np.inf) 8 ctx1Fcadgfeb

748 passed = (np.abs(metric) < 0.15).astype(float) 8 ctx1Fcadgfeb

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

750 metric[~data['correct']] = passed[~data['correct']] = np.nan 8 ctx1Fcadgfeb

751 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1Fcadgfeb

752 return metric, passed 8 ctx1Fcadgfeb

753 

754 

755def check_negative_feedback_stimOff_delays(data, **_): (empty)

756 """ 

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

758 

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

760 turning off is 2 ± 0.150 seconds. 

761 

762 Metric: 

763 M = stimOff_times - errorCue_times - 2s 

764 

765 Criterion: 

766 |M| < 0.150 s 

767 

768 Units: 

769 seconds [s] 

770 

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

772 """ 

773 metric = np.nan_to_num(data['stimOff_times'] - data['errorCue_times'] - 2, nan=np.inf) 8 ctx1Gcadgfeb

774 # Apply criteria 

775 passed = (np.abs(metric) < 0.15).astype(float) 8 ctx1Gcadgfeb

776 # Remove none negative feedback trials 

777 metric[data['correct']] = passed[data['correct']] = np.nan 8 ctx1Gcadgfeb

778 assert data['intervals'].shape[0] == len(metric) == len(passed) 8 ctx1Gcadgfeb

779 return metric, passed 8 ctx1Gcadgfeb

780 

781 

782# === Wheel movement during trial checks === 

783 

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

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

786 

787 Metric: 

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