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

1"""Behaviour QC 

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

3 

4Examples 

5-------- 

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

7 

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

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

10 

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

12 

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

14>>> qc = TaskQC(eid) 

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

16 

17Inspecting individual test outcomes 

18 

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

20>>> qc = TaskQC(eid) 

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

22 

23Running bpod QC on ephys session 

24 

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

26>>> qc = TaskQC(eid) 

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

28>>> bpod_qc = qc.run() 

29 

30Running bpod QC only, from training rig PC 

31 

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

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

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

35>>> qc = TaskQC(session_path) 

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

37>>> qc.run() 

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

39 

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

41 

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

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

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

45>>> qc = TaskQC(session_path) 

46>>> qc.run() 

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

48""" 

49import logging (empty)

50import sys (empty)

51from datetime import datetime, timedelta (empty)

52from inspect import getmembers, isfunction (empty)

53from functools import reduce (empty)

54from collections.abc import Sized (empty)

55 

56import numpy as np (empty)

57from scipy.stats import chisquare (empty)

58 

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

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

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

62from one.alf.spec import is_session_path (empty)

63from . import base (empty)

64 

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

66 

67 

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

69 """A class for computing task QC metrics""" 

70 

71 criteria = dict() (empty)

72 criteria['default'] = {"PASS": 0.99, "WARNING": 0.90, "FAIL": 0} # Note: WARNING was 0.95 prior to Aug 2022 (empty)

73 criteria['_task_stimOff_itiIn_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

74 criteria['_task_positive_feedback_stimOff_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

75 criteria['_task_negative_feedback_stimOff_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

76 criteria['_task_wheel_move_during_closed_loop'] = {"PASS": 0.99, "WARNING": 0} (empty)

77 criteria['_task_response_stimFreeze_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

78 criteria['_task_detected_wheel_moves'] = {"PASS": 0.99, "WARNING": 0} (empty)

79 criteria['_task_trial_length'] = {"PASS": 0.99, "WARNING": 0} (empty)

80 criteria['_task_goCue_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

81 criteria['_task_errorCue_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

82 criteria['_task_stimOn_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

83 criteria['_task_stimOff_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

84 criteria['_task_stimFreeze_delays'] = {"PASS": 0.99, "WARNING": 0} (empty)

85 criteria['_task_iti_delays'] = {"NOT_SET": 0} (empty)

86 criteria['_task_passed_trial_checks'] = {"NOT_SET": 0} (empty)

87 

88 @staticmethod (empty)

89 def _thresholding(qc_value, thresholds=None): (empty)

90 """ 

91 Computes the outcome of a single key by applying thresholding. 

92 :param qc_value: proportion of passing qcs, between 0 and 1 

93 :param thresholds: dictionary with keys 'PASS', 'WARNING', 'FAIL' 

94 (cf. TaskQC.criteria attribute) 

95 :return: int where -1: NOT_SET, 0: PASS, 1: WARNING, 2: FAIL 

96 """ 

97 MAX_BOUND, MIN_BOUND = (1, 0) 14 ctx1vonpkjfcabehgd

98 if not thresholds: 14 ctx1vonpkjfcabehgd

99 thresholds = TaskQC.criteria['default'].copy() 

100 if qc_value is None or np.isnan(qc_value): 14 ctx1vonpkjfcabehgd

101 return int(-1) 3 ctx1onp

102 elif (qc_value > MAX_BOUND) or (qc_value < MIN_BOUND): 14 ctx1vonpkjfcabehgd

103 raise ValueError("Values out of bound") 1 ctx1v

104 if 'PASS' in thresholds.keys() and qc_value >= thresholds['PASS']: 13 ctx1onpkjfcabehgd

105 return 0 13 ctx1onpkjfcabehgd

106 if 'WARNING' in thresholds.keys() and qc_value >= thresholds['WARNING']: 13 ctx1onpkjfcabehgd

107 return 1 11 ctx1opjfcabehgd

108 if 'FAIL' in thresholds and qc_value >= thresholds['FAIL']: 12 ctx1onkjfcabehgd

109 return 2 10 ctx1okjfcabhgd

110 if 'NOT_SET' in thresholds and qc_value >= thresholds['NOT_SET']: 9 ctx1nfcabehgd

111 return -1 9 ctx1nfcabehgd

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

113 return -1 

114 

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

116 """ 

117 :param session_path_or_eid: A session eid or path 

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

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

120 """ 

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

122 self.download_data = not is_session_path(session_path_or_eid) 8 ctx1ljfcabed

123 super().__init__(session_path_or_eid, **kwargs) 8 ctx1ljfcabed

124 

125 # Data 

126 self.extractor = None 8 ctx1ljfcabed

127 

128 # Metrics and passed trials 

129 self.metrics = None 8 ctx1ljfcabed

130 self.passed = None 8 ctx1ljfcabed

131 

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

133 """Extract the data from raw data files 

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

135 

136 :param bpod_only: if True no data is extracted from the FPGA for ephys sessions 

137 :param download_data: if True, any missing raw data is downloaded via ONE. 

138 """ 

139 self.extractor = TaskQCExtractor( 3 ctx1lab

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

141 

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

143 """Compute and store the QC metrics 

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

145 test, and a map of which datapoints passed for each test 

146 :param bpod_only: if True no data is extracted from the FPGA for ephys sessions 

147 :param download_data: if True, any missing raw data is downloaded via ONE. By default 

148 data are not downloaded if a session path was provided to the constructor. 

149 :return: 

150 """ 

151 if self.extractor is None: 9 ctx1fcabeihgd

152 kwargs['download_data'] = kwargs.pop('download_data', self.download_data) 2 ctx1ab

153 self.load_data(**kwargs) 2 ctx1ab

154 self.log.info(f"Session {self.session_path}: Running QC on behavior data...") 9 ctx1fcabeihgd

155 self.metrics, self.passed = get_bpodqc_metrics_frame( 9 ctx1fcabeihgd

156 self.extractor.data, 

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

158 photodiode=self.extractor.frame_ttls, 

159 audio=self.extractor.audio_ttls, 

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

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

162 ) 

163 return 9 ctx1fcabeihgd

164 

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

166 """ 

167 :param update: if True, updates the session QC fields on Alyx 

168 :param bpod_only: if True no data is extracted from the FPGA for ephys sessions 

169 :param download_data: if True, any missing raw data is downloaded via ONE. By default 

170 data are not downloaded if a session path was provided to the constructor. 

171 :return: QC outcome (str), a dict for extended QC 

172 """ 

173 if self.metrics is None: 8 ctx1jfcabegd

174 self.compute(**kwargs) 8 ctx1jfcabegd

175 outcome, results, _ = self.compute_session_status() 8 ctx1jfcabegd

176 if update: 8 ctx1jfcabegd

177 self.update_extended_qc(results) 4 ctx1fcgd

178 self.update(outcome, namespace) 4 ctx1fcgd

179 return outcome, results 8 ctx1jfcabegd

180 

181 @staticmethod (empty)

182 def compute_session_status_from_dict(results): (empty)

183 """ 

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

185 in a single value 

186 :param results: a dictionary of qc keys containing (usually scalar) values 

187 :return: Overall session QC outcome as a string 

188 :return: A dict of QC tests and their outcomes 

189 """ 

190 indices = np.zeros(len(results), dtype=int) 14 ctx1vonpkjfcabehgd

191 for i, k in enumerate(results): 14 ctx1vonpkjfcabehgd

192 if k in TaskQC.criteria.keys(): 14 ctx1vonpkjfcabehgd

193 indices[i] = TaskQC._thresholding(results[k], thresholds=TaskQC.criteria[k]) 12 ctx1npkjfcabehgd

194 else: 

195 indices[i] = TaskQC._thresholding(results[k], thresholds=TaskQC.criteria['default']) 14 ctx1vonpkjfcabehgd

196 

197 def key_map(x): 13 ctx1onpkjfcabehgd

198 return 'NOT_SET' if x < 0 else list(TaskQC.criteria['default'].keys())[x] 13 ctx1onpkjfcabehgd

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

200 session_outcome = key_map(max(indices)) 13 ctx1onpkjfcabehgd

201 outcomes = dict(zip(results.keys(), map(key_map, indices))) 13 ctx1onpkjfcabehgd

202 return session_outcome, outcomes 13 ctx1onpkjfcabehgd

203 

204 def compute_session_status(self): (empty)

205 """ 

206 Computes the overall session QC for each key and aggregates in a single value 

207 :return: Overall session QC outcome as a string 

208 :return: A dict of QC tests and the proportion of data points that passed them 

209 :return: A dict of QC tests and their outcomes 

210 """ 

211 if self.passed is None: 10 ctx1kjfcabehgd

212 raise AttributeError('passed is None; compute QC first') 1 ctx1h

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

214 results = {k: None if v is None or np.isnan(v).all() else np.nanmean(v) 10 ctx1kjfcabehgd

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

216 session_outcome, outcomes = self.compute_session_status_from_dict(results) 10 ctx1kjfcabehgd

217 return session_outcome, results, outcomes 10 ctx1kjfcabehgd

218 

219 

220class HabituationQC(TaskQC): (empty)

221 

222 def compute(self, download_data=None): (empty)

223 """Compute and store the QC metrics 

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

225 test, and a map of which datapoints passed for each test 

226 :return: 

227 """ 

228 if self.extractor is None: 2 ctx1kj

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

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

231 self.load_data(download_data=ensure_data) 

232 self.log.info(f"Session {self.session_path}: Running QC on habituation data...") 2 ctx1kj

233 

234 # Initialize checks 

235 prefix = '_task_' 2 ctx1kj

236 data = self.extractor.data 2 ctx1kj

237 metrics = {} 2 ctx1kj

238 passed = {} 2 ctx1kj

239 

240 # Check all reward volumes == 3.0ul 

241 check = prefix + 'reward_volumes' 2 ctx1kj

242 metrics[check] = data['rewardVolume'] 2 ctx1kj

243 passed[check] = metrics[check] == 3.0 2 ctx1kj

244 

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

246 check = prefix + 'habituation_time' 2 ctx1kj

247 if not self.one or not self.session_path: 2 ctx1kj

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

249 metrics[check] = passed[check] = None 

250 else: 

251 subject, session_date = self.session_path.parts[-3:-1] 2 ctx1kj

252 # compute from the date specified 

253 date_minus_week = ( 2 ctx1kj

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

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

256 sessions = self.one.alyx.rest('sessions', 'list', subject=subject, 2 ctx1kj

257 date_range=[date_minus_week, session_date], 

258 task_protocol='habituation') 

259 # Remove the current session if already registered 

260 if sessions and sessions[0]['start_time'].startswith(session_date): 2 ctx1kj

261 sessions = sessions[1:] 

262 metric = ([0, data['intervals'][-1, 1] - data['intervals'][0, 0]] + 2 ctx1kj

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

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

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

266 

267 # The duration from raw trial data 

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

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

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

271 metrics[check] = np.array(metric) 2 ctx1kj

272 passed[check] = np.diff(metric) >= 12 2 ctx1kj

273 

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

275 check = prefix + 'trial_event_sequence' 2 ctx1kj

276 nans = ( 2 ctx1kj

277 np.isnan(data["intervals"][:, 0]) | # noqa 

278 np.isnan(data["stimOn_times"]) | # noqa 

279 np.isnan(data["stimCenter_times"]) | 

280 np.isnan(data["valveOpen_times"]) | # noqa 

281 np.isnan(data["stimOff_times"]) 

282 ) 

283 a = np.less(data["intervals"][:, 0], data["stimOn_times"], where=~nans) 2 ctx1kj

284 b = np.less(data["stimOn_times"], data["stimCenter_times"], where=~nans) 2 ctx1kj

285 c = np.less(data["stimCenter_times"], data["valveOpen_times"], where=~nans) 2 ctx1kj

286 d = np.less(data["valveOpen_times"], data["stimOff_times"], where=~nans) 2 ctx1kj

287 

288 metrics[check] = a & b & c & d & ~nans 2 ctx1kj

289 passed[check] = metrics[check].astype(float) 2 ctx1kj

290 

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

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

293 check = prefix + 'stimCenter_delays' 2 ctx1kj

294 metric = np.nan_to_num(data["stimCenter_times"] - data["stimCenterTrigger_times"], 2 ctx1kj

295 nan=np.inf) 

296 passed[check] = (metric <= 0.15) & (metric > 0) 2 ctx1kj

297 metrics[check] = metric 2 ctx1kj

298 

299 # Phase check 

300 check = prefix + 'phase' 2 ctx1kj

301 metric = data['phase'] 2 ctx1kj

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

303 metrics[check] = metric 2 ctx1kj

304 

305 check = prefix + 'phase_distribution' 2 ctx1kj

306 metric, _ = np.histogram(data['phase']) 2 ctx1kj

307 _, p = chisquare(metric) 2 ctx1kj

308 passed[check] = p < 0.05 2 ctx1kj

309 metrics[check] = metric 2 ctx1kj

310 

311 # Checks common to training QC 

312 checks = [check_goCue_delays, check_stimOn_goCue_delays, 2 ctx1kj

313 check_stimOn_delays, check_stimOff_delays] 

314 for fcn in checks: 2 ctx1kj

315 check = prefix + fcn.__name__[6:] 2 ctx1kj

316 metrics[check], passed[check] = fcn(data) 2 ctx1kj

317 

318 self.metrics, self.passed = (metrics, passed) 2 ctx1kj

319 

320 

321def get_bpodqc_metrics_frame(data, **kwargs): (empty)

322 """ 

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

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

325 :param data: dict of extracted task data 

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

327 :param enc_res: the rotary encoder resolution 

328 :param wheel_gain: the STIM_GAIN task parameter 

329 :param photodiode: the fronts from Bpod's BNC1 input or FPGA frame2ttl channel 

330 :param audio: the fronts from Bpod's BNC2 input FPGA audio sync channel 

331 :param min_qt: the QUIESCENT_PERIOD task parameter 

332 :return metrics: dict of checks and their QC metrics 

333 :return passed: dict of checks and a float array of which samples passed 

334 """ 

335 def is_metric(x): 9 ctx1fcabeihgd

336 return isfunction(x) and x.__name__.startswith('check_') 9 ctx1fcabeihgd

337 # Find all methods that begin with 'check_' 

338 checks = getmembers(sys.modules[__name__], is_metric) 9 ctx1fcabeihgd

339 prefix = '_task_' # Extended QC fields will start with this 9 ctx1fcabeihgd

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

341 qc_metrics_map = {prefix + k[6:]: fn(data, **kwargs) for k, fn in checks} 9 ctx1fcabeihgd

342 

343 # Split metrics and passed frames 

344 metrics = {} 9 ctx1fcabeihgd

345 passed = {} 9 ctx1fcabeihgd

346 for k in qc_metrics_map: 9 ctx1fcabeihgd

347 metrics[k], passed[k] = qc_metrics_map[k] 9 ctx1fcabeihgd

348 

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

350 n_trials = data['intervals'].shape[0] 9 ctx1fcabeihgd

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

352 trial_level_passed = [m for m in passed.values() 9 ctx1fcabeihgd

353 if isinstance(m, Sized) and len(m) == n_trials] 

354 name = prefix + 'passed_trial_checks' 9 ctx1fcabeihgd

355 metrics[name] = reduce(np.logical_and, trial_level_passed or (None, None)) 9 ctx1fcabeihgd

356 passed[name] = metrics[name].astype(float) if trial_level_passed else None 9 ctx1fcabeihgd

357 

358 return metrics, passed 9 ctx1fcabeihgd

359 

360 

361# SINGLE METRICS 

362# ---------------------------------------------------------------------------- # 

363 

364# === Delays between events checks === 

365 

366def check_stimOn_goCue_delays(data, **_): (empty)

367 """ Checks that the time difference between the onset of the visual stimulus 

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

369 

370 Metric: M = stimOn_times - goCue_times 

371 Criteria: 0 < M < 0.010 s 

372 Units: seconds [s] 

373 

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

375 """ 

376 # Calculate the difference between stimOn and goCue times. 

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

378 metric = np.nan_to_num(data["goCue_times"] - data["stimOn_times"], nan=np.inf) 12 ctx1kGjfcabeihgd

379 passed = (metric < 0.01) & (metric > 0) 12 ctx1kGjfcabeihgd

380 assert data["intervals"].shape[0] == len(metric) == len(passed) 12 ctx1kGjfcabeihgd

381 return metric, passed 12 ctx1kGjfcabeihgd

382 

383 

384def check_response_feedback_delays(data, **_): (empty)

385 """ Checks that the time difference between the response and the feedback onset 

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

387 

388 Metric: M = feedback_time - response_time 

389 Criterion: 0 < M < 0.010 s 

390 Units: seconds [s] 

391 

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

393 """ 

394 metric = np.nan_to_num(data["feedback_times"] - data["response_times"], nan=np.inf) 10 ctx1Hfcabeihgd

395 passed = (metric < 0.01) & (metric > 0) 10 ctx1Hfcabeihgd

396 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Hfcabeihgd

397 return metric, passed 10 ctx1Hfcabeihgd

398 

399 

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

401 """ Checks that the time difference between the visual stimulus freezing and the 

402 response is positive and less than 100ms. 

403 

404 Metric: M = (stimFreeze_times - response_times) 

405 Criterion: 0 < M < 0.100 s 

406 Units: seconds [s] 

407 

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

409 'choice') 

410 """ 

411 # Calculate the difference between stimOn and goCue times. 

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

413 metric = np.nan_to_num(data["stimFreeze_times"] - data["response_times"], nan=np.inf) 10 ctx1Bfcabeihgd

414 # Test for valid values 

415 passed = ((metric < 0.1) & (metric > 0)).astype(float) 10 ctx1Bfcabeihgd

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

417 # These values are ignored in calculation of proportion passed 

418 passed[data["choice"] == 0] = np.nan 10 ctx1Bfcabeihgd

419 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Bfcabeihgd

420 return metric, passed 10 ctx1Bfcabeihgd

421 

422 

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

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

425 

426 Metric: M = itiIn_times - stimOff_times 

427 Criterion: 0 < M < 0.010 s 

428 Units: seconds [s] 

429 

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

431 'choice') 

432 """ 

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

434 metric = np.nan_to_num(data["itiIn_times"] - data["stimOff_times"], nan=np.inf) 10 ctx1Cfcabeihgd

435 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 10 ctx1Cfcabeihgd

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

437 # NaN values are ignored in calculation of proportion passed 

438 metric[data["choice"] == 0] = passed[data["choice"] == 0] = np.nan 10 ctx1Cfcabeihgd

439 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Cfcabeihgd

440 return metric, passed 10 ctx1Cfcabeihgd

441 

442 

443def check_iti_delays(data, **_): (empty)

444 """ Check that the period of gray screen between stim off and the start of the next trial is 

445 0.5s +/- 200%. 

446 

447 Metric: M = stimOff (n) - trialStart (n+1) - 0.5 

448 Criterion: |M| < 1 

449 Units: seconds [s] 

450 

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

452 """ 

453 # Initialize array the length of completed trials 

454 metric = np.full(data["intervals"].shape[0], np.nan) 10 ctx1Afcabeihgd

455 passed = metric.copy() 10 ctx1Afcabeihgd

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

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

458 metric[:-1] = \ 10 ctx1Afcabeihgd

459 np.nan_to_num(data["intervals"][1:, 0] - data["stimOff_times"][:-1] - 0.5, nan=np.inf) 

460 passed[:-1] = np.abs(metric[:-1]) < .5 # Last trial is not counted 10 ctx1Afcabeihgd

461 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Afcabeihgd

462 return metric, passed 10 ctx1Afcabeihgd

463 

464 

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

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

467 is 1 ± 0.150 seconds. 

468 

469 Metric: M = stimOff_times - feedback_times - 1s 

470 Criterion: |M| < 0.150 s 

471 Units: seconds [s] 

472 

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

474 'correct') 

475 """ 

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

477 metric = np.nan_to_num(data["stimOff_times"] - data["feedback_times"] - 1, nan=np.inf) 10 ctx1Dfcabeihgd

478 passed = (np.abs(metric) < 0.15).astype(float) 10 ctx1Dfcabeihgd

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

480 metric[~data["correct"]] = passed[~data["correct"]] = np.nan 10 ctx1Dfcabeihgd

481 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Dfcabeihgd

482 return metric, passed 10 ctx1Dfcabeihgd

483 

484 

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

486 """ Check that the time difference between the error sound and the visual stimulus 

487 turning off is 2 ± 0.150 seconds. 

488 

489 Metric: M = stimOff_times - errorCue_times - 2s 

490 Criterion: |M| < 0.150 s 

491 Units: seconds [s] 

492 

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

494 """ 

495 metric = np.nan_to_num(data["stimOff_times"] - data["errorCue_times"] - 2, nan=np.inf) 10 ctx1Efcabeihgd

496 # Apply criteria 

497 passed = (np.abs(metric) < 0.15).astype(float) 10 ctx1Efcabeihgd

498 # Remove none negative feedback trials 

499 metric[data["correct"]] = passed[data["correct"]] = np.nan 10 ctx1Efcabeihgd

500 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1Efcabeihgd

501 return metric, passed 10 ctx1Efcabeihgd

502 

503 

504# === Wheel movement during trial checks === 

505 

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

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

508 

509 Metric: M = (w_t - 0.05) - (w_t + 0.05), where t = feedback_times 

510 Criterion: M != 0 

511 Units: radians 

512 

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

514 'intervals', 'feedback_times') 

515 """ 

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

517 traces = traces_by_trial( 10 ctx1rfcabeihgd

518 data["wheel_timestamps"], 

519 data["wheel_position"], 

520 start=data["feedback_times"] - 0.05, 

521 end=data["feedback_times"] + 0.05, 

522 ) 

523 metric = np.zeros_like(data["feedback_times"]) 10 ctx1rfcabeihgd

524 # For each trial find the displacement 

525 for i, trial in enumerate(traces): 10 ctx1rfcabeihgd

526 pos = trial[1] 10 ctx1rfcabeihgd

527 if pos.size > 1: 10 ctx1rfcabeihgd

528 metric[i] = pos[-1] - pos[0] 10 ctx1rfcabeihgd

529 

530 # except no-go trials 

531 metric[data["choice"] == 0] = np.nan # NaN = trial ignored for this check 10 ctx1rfcabeihgd

532 nans = np.isnan(metric) 10 ctx1rfcabeihgd

533 passed = np.zeros_like(metric) * np.nan 10 ctx1rfcabeihgd

534 

535 passed[~nans] = (metric[~nans] != 0).astype(float) 10 ctx1rfcabeihgd

536 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1rfcabeihgd

537 return metric, passed 10 ctx1rfcabeihgd

538 

539 

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

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

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

543 

544 Metric: M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

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

546 move 35 visual degrees 

547 Criterion: displacement < tol visual degree 

548 Units: degrees angle of wheel turn 

549 

550 :param re_ts: extracted wheel timestamps in seconds 

551 :param re_pos: extracted wheel positions in radians 

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

553 position, choice, intervals) 

554 :param wheel_gain: the 'STIM_GAIN' task setting 

555 :param tol: the criterion in visual degrees 

556 """ 

557 if wheel_gain is None: 10 ctx1mfcabeihgd

558 _log.warning("No wheel_gain input in function call, returning None") 

559 return None, None 

560 

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

562 traces = traces_by_trial(re_ts, re_pos, 10 ctx1mfcabeihgd

563 start=data["goCueTrigger_times"], 

564 end=data["response_times"]) 

565 

566 metric = np.zeros_like(data["feedback_times"]) 10 ctx1mfcabeihgd

567 # For each trial find the absolute displacement 

568 for i, trial in enumerate(traces): 10 ctx1mfcabeihgd

569 t, pos = trial 10 ctx1mfcabeihgd

570 if pos.size != 0: 10 ctx1mfcabeihgd

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

572 idx = np.abs(re_ts - t[0]).argmin() - 1 10 ctx1mfcabeihgd

573 origin = re_pos[idx] 10 ctx1mfcabeihgd

574 metric[i] = np.abs(pos - origin).max() 10 ctx1mfcabeihgd

575 

576 # Load wheel_gain and thresholds for each trial 

577 wheel_gain = np.array([wheel_gain] * len(data["position"])) 10 ctx1mfcabeihgd

578 thresh = data["position"] 10 ctx1mfcabeihgd

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

580 s_mm = np.abs(thresh / wheel_gain) # don't care about direction 10 ctx1mfcabeihgd

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

582 metric = metric - criterion # difference should be close to 0 10 ctx1mfcabeihgd

583 rad_per_deg = cm_to_rad(1 / wheel_gain * 1e-1) 10 ctx1mfcabeihgd

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

585 metric[data["choice"] == 0] = passed[data["choice"] == 0] = np.nan # except no-go trials 10 ctx1mfcabeihgd

586 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1mfcabeihgd

587 return metric, passed 10 ctx1mfcabeihgd

588 

589 

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

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

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

593 

594 Metric: M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

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

596 move 35 visual degrees 

597 Criterion: displacement < 3 visual degrees 

598 Units: degrees angle of wheel turn 

599 

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

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

602 :param wheel_gain: the 'STIM_GAIN' task setting 

603 """ 

604 # Get the Bpod extracted wheel data 

605 timestamps = data['wheel_timestamps'] 10 ctx1mfcabeihgd

606 position = data['wheel_position'] 10 ctx1mfcabeihgd

607 

608 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=3) 10 ctx1mfcabeihgd

609 

610 

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

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

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

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

615 

616 Metric: M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response 

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

618 move 35 visual degrees 

619 Criterion: displacement < 1 visual degree 

620 Units: degrees angle of wheel turn 

621 

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

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

624 :param wheel_gain: the 'STIM_GAIN' task setting 

625 """ 

626 # Get the Bpod extracted wheel data 

627 timestamps = data.get('wheel_timestamps_bpod', data['wheel_timestamps']) 9 ctx1fcabeihgd

628 position = data.get('wheel_position_bpod', data['wheel_position']) 9 ctx1fcabeihgd

629 

630 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=1) 9 ctx1fcabeihgd

631 

632 

633def check_wheel_freeze_during_quiescence(data, **_): (empty)

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

635 quiescence interval before the stimulus appears. 

636 

637 Metric: M = |max(W) - min(W)| where W is wheel pos over quiescence interval 

638 interval = [stimOnTrigger_times - quiescent_duration, stimOnTrigger_times] 

639 Criterion: M < 2 degrees 

640 Units: degrees angle of wheel turn 

641 

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

643 'intervals', 'stimOnTrigger_times') 

644 """ 

645 assert np.all(np.diff(data["wheel_timestamps"]) >= 0) 10 ctx1qfcabeihgd

646 assert data["quiescence"].size == data["stimOnTrigger_times"].size 10 ctx1qfcabeihgd

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

648 qevt_start_times = data["stimOnTrigger_times"] - data["quiescence"] 10 ctx1qfcabeihgd

649 traces = traces_by_trial( 10 ctx1qfcabeihgd

650 data["wheel_timestamps"], 

651 data["wheel_position"], 

652 start=qevt_start_times, 

653 end=data["stimOnTrigger_times"] 

654 ) 

655 

656 metric = np.zeros((len(data["quiescence"]), 2)) # (n_trials, n_directions) 10 ctx1qfcabeihgd

657 for i, trial in enumerate(traces): 10 ctx1qfcabeihgd

658 t, pos = trial 10 ctx1qfcabeihgd

659 # Get the last position before the period began 

660 if pos.size > 0: 10 ctx1qfcabeihgd

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

662 idx = np.abs(data["wheel_timestamps"] - t[0]).argmin() - 1 6 ctx1qcabed

663 origin = data["wheel_position"][idx if idx != -1 else 0] 6 ctx1qcabed

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

665 metric[i, :] = np.abs([np.min(pos - origin), np.max(pos - origin)]) 6 ctx1qcabed

666 # Reduce to the largest displacement found in any direction 

667 metric = np.max(metric, axis=1) 10 ctx1qfcabeihgd

668 metric = 180 * metric / np.pi # convert to degrees from radians 10 ctx1qfcabeihgd

669 criterion = 2 # Position shouldn't change more than 2 in either direction 10 ctx1qfcabeihgd

670 passed = metric < criterion 10 ctx1qfcabeihgd

671 assert data["intervals"].shape[0] == len(metric) == len(passed) 10 ctx1qfcabeihgd

672 return metric, passed 10 ctx1qfcabeihgd

673 

674 

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

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

677 

678 Metric: M = firstMovement times 

679 Criterion: (goCue trigger time - min quiescent period) < M < response time 

680 Units: Seconds [s] 

681 

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

683 'response_times', 'choice', 'intervals') 

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

685 """ 

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

687 min_qt = np.array(min_qt) 10 ctx1wfcabeihgd

688 if min_qt.size > data["intervals"].shape[0]: 10 ctx1wfcabeihgd

689 min_qt = min_qt[:data["intervals"].shape[0]] 1 ctx1e

690 

691 metric = data['firstMovement_times'] 10 ctx1wfcabeihgd

692 qevt_start = data['goCueTrigger_times'] - np.array(min_qt) 10 ctx1wfcabeihgd

693 response = data['response_times'] 10 ctx1wfcabeihgd

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

695 passed = np.array([a < m < b for m, a, b in zip(metric, qevt_start, response)], dtype=float) 10 ctx1wfcabeihgd

696 nogo = data['choice'] == 0 10 ctx1wfcabeihgd

697 passed[nogo] = np.nan # No go trial may have no movement times and that's fine 10 ctx1wfcabeihgd

698 return metric, passed 10 ctx1wfcabeihgd

699 

700 

701# === Sequence of events checks === 

702 

703def check_error_trial_event_sequence(data, **_): (empty)

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

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

706 in the correct order 

707 

708 Metric: M = Bpod (trial start) > audio (go cue) > audio (error) > Bpod (ITI) > Bpod (trial end) 

709 Criterion: M == True 

710 Units: -none- 

711 

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

713 'itiIn_times', 'correct') 

714 """ 

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

716 nans = ( 10 ctx1tfcabeihgd

717 np.isnan(data["intervals"][:, 0]) | 

718 np.isnan(data["goCue_times"]) | # noqa 

719 np.isnan(data["errorCue_times"]) | # noqa 

720 np.isnan(data["itiIn_times"]) | # noqa 

721 np.isnan(data["intervals"][:, 1]) 

722 ) 

723 

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

725 a = np.less(data["intervals"][:, 0], data["goCue_times"], where=~nans) # Start time < go cue 10 ctx1tfcabeihgd