Coverage for brainbox/behavior/training.py: 64%

1"""Computing and testing IBL training status criteria. 

2 

3For an in-depth description of each training status, see `Appendix 2`_ of the IBL Protocol For Mice 

4Training. 

5 

6.. _Appendix 2: https://figshare.com/articles/preprint/A_standardized_and_reproducible_method_to_\ 

7measure_decision-making_in_mice_Appendix_2_IBL_protocol_for_mice_training/11634729 

8 

9Examples 

10-------- 

11Plot the psychometric curve for a given session. 

12 

13>>> trials = ONE().load_object(eid, 'trials') 

14>>> fix, ax = plot_psychometric(trials) 

15 

16Compute 'response times', defined as the duration of open-loop for each contrast. 

17 

18>>> reaction_time, contrasts, n_contrasts = compute_reaction_time(trials) 

19 

20Compute 'reaction times', defined as the time between go cue and first detected movement. 

21NB: These may be negative! 

22 

23>>> reaction_time, contrasts, n_contrasts = compute_reaction_time( 

24... trials, stim_on_type='goCue_times', stim_off_type='firstMovement_times') 

25 

26Compute 'response times', defined as the time between first detected movement and response. 

27 

28>>> reaction_time, contrasts, n_contrasts = compute_reaction_time( 

29... trials, stim_on_type='firstMovement_times', stim_off_type='response_times') 

30 

31Compute 'movement times', defined as the time between last detected movement and response threshold. 

32 

33>>> import brainbox.behavior.wheel as wh 

34>>> wheel_moves = ONE().load_object(eid, 'wheeMoves') 

35>>> trials['lastMovement_times'] = wh.get_movement_onset(wheel_moves.intervals, trial_data.response_times) 

36>>> reaction_time, contrasts, n_contrasts = compute_reaction_time( 

37... trials, stim_on_type='lastMovement_times', stim_off_type='response_times') 

38 

39""" 

40import logging (empty)

41import datetime (empty)

42import re (empty)

43from enum import IntFlag, auto, unique (empty)

44 

45import numpy as np (empty)

46import matplotlib (empty)

47import matplotlib.pyplot as plt (empty)

48import seaborn as sns (empty)

49import pandas as pd (empty)

50from scipy.stats import bootstrap (empty)

51from iblutil.util import Bunch (empty)

52from one.api import ONE (empty)

53from one.alf.io import AlfBunch (empty)

54from one.alf.exceptions import ALFObjectNotFound (empty)

55import psychofit as psy (empty)

56 

57_logger = logging.getLogger('ibllib') (empty)

58 

59TRIALS_KEYS = ['contrastLeft', (empty)

60 'contrastRight', 

61 'feedbackType', 

62 'probabilityLeft', 

63 'choice', 

64 'response_times', 

65 'stimOn_times'] 

66"""list of str: The required keys in the trials object for computing training status.""" (empty)

67 

68 

69@unique (empty)

70class TrainingStatus(IntFlag): (empty)

71 """Standard IBL training criteria. 

72 

73 Enumeration allows for comparisons between training status. 

74 

75 Examples 

76 -------- 

77 >>> status = 'ready4delay' 

78 ... assert TrainingStatus[status.upper()] is TrainingStatus.READY4DELAY 

79 ... assert TrainingStatus[status.upper()] not in TrainingStatus.FAILED, 'Subject failed training' 

80 ... assert TrainingStatus[status.upper()] >= TrainingStatus.TRAINED, 'Subject untrained' 

81 ... assert TrainingStatus[status.upper()] > TrainingStatus.IN_TRAINING, 'Subject untrained' 

82 ... assert TrainingStatus[status.upper()] in ~TrainingStatus.FAILED, 'Subject untrained' 

83 ... assert TrainingStatus[status.upper()] in TrainingStatus.TRAINED ^ TrainingStatus.READY 

84 

85 Get the next training status 

86 

87 >>> next(member for member in sorted(TrainingStatus) if member > TrainingStatus[status.upper()]) 

88 <TrainingStatus.READY4RECORDING: 128> 

89 

90 Notes 

91 ----- 

92 - ~TrainingStatus.TRAINED means any status but trained 1a or trained 1b. 

93 - A subject may acheive both TRAINED_1A and TRAINED_1B within a single session, therefore it 

94 is possible to have skipped the TRAINED_1A session status. 

95 """ 

96 UNTRAINABLE = auto() (empty)

97 UNBIASABLE = auto() (empty)

98 IN_TRAINING = auto() (empty)

99 TRAINED_1A = auto() (empty)

100 TRAINED_1B = auto() (empty)

101 READY4EPHYSRIG = auto() (empty)

102 READY4DELAY = auto() (empty)

103 READY4RECORDING = auto() (empty)

104 # Compound training statuses for convenience 

105 FAILED = UNTRAINABLE | UNBIASABLE (empty)

106 READY = READY4EPHYSRIG | READY4DELAY | READY4RECORDING (empty)

107 TRAINED = TRAINED_1A | TRAINED_1B (empty)

108 

109 

110def get_lab_training_status(lab, date=None, details=True, one=None): (empty)

111 """ 

112 Computes the training status of all alive and water restricted subjects in a specified lab. 

113 

114 The response are printed to std out. 

115 

116 Parameters 

117 ---------- 

118 lab : str 

119 Lab name (must match the name registered on Alyx). 

120 date : str 

121 The ISO date from which to compute training status. If not specified will compute from the 

122 latest date with available data. Format should be 'YYYY-MM-DD'. 

123 details : bool 

124 Whether to display all information about training status computation e.g. performance, 

125 number of trials, psychometric fit parameters. 

126 one : one.api.OneAlyx 

127 An instance of ONE. 

128 

129 """ 

130 one = one or ONE() 

131 subj_lab = one.alyx.rest('subjects', 'list', lab=lab, alive=True, water_restricted=True) 

132 subjects = [subj['nickname'] for subj in subj_lab] 

133 for subj in subjects: 

134 get_subject_training_status(subj, date=date, details=details, one=one) 

135 

136 

137def get_subject_training_status(subj, date=None, details=True, one=None): (empty)

138 """ 

139 Computes the training status of specified subject and prints results to std out. 

140 

141 Parameters 

142 ---------- 

143 subj : str 

144 Subject nickname (must match the name registered on Alyx). 

145 date : str 

146 The ISO date from which to compute training status. If not specified will compute from the 

147 latest date with available data. Format should be 'YYYY-MM-DD'. 

148 details : bool 

149 Whether to display all information about training status computation e.g. performance, 

150 number of trials, psychometric fit parameters. 

151 one : one.api.OneAlyx 

152 An instance of ONE. 

153 """ 

154 one = one or ONE() 

155 

156 trials, task_protocol, ephys_sess, n_delay = get_sessions(subj, date=date, one=one) 

157 if not trials: 

158 return 

159 sess_dates = list(trials.keys()) 

160 status, info, _ = get_training_status(trials, task_protocol, ephys_sess, n_delay) 

161 

162 if details: 

163 if np.any(info.get('psych')): 

164 display_status(subj, sess_dates, status, perf_easy=info.perf_easy, 

165 n_trials=info.n_trials, psych=info.psych, rt=info.rt) 

166 elif np.any(info.get('psych_20')): 

167 display_status(subj, sess_dates, status, perf_easy=info.perf_easy, 

168 n_trials=info.n_trials, psych_20=info.psych_20, psych_80=info.psych_80, 

169 rt=info.rt) 

170 else: 

171 display_status(subj, sess_dates, status) 

172 

173 

174def get_sessions(subj, date=None, one=None): (empty)

175 """ 

176 Download and load in training data for a specified subject. If a date is given it will load 

177 data from the three (or as many as are available) previous sessions up to the specified date. 

178 If not it will load data from the last three training sessions that have data available. 

179 

180 Parameters 

181 ---------- 

182 subj : str 

183 Subject nickname (must match the name registered on Alyx). 

184 date : str 

185 The ISO date from which to compute training status. If not specified will compute from the 

186 latest date with available data. Format should be 'YYYY-MM-DD'. 

187 one : one.api.OneAlyx 

188 An instance of ONE. 

189 

190 Returns 

191 ------- 

192 iblutil.util.Bunch 

193 Dictionary of trials objects where each key is the ISO session date string. 

194 list of str 

195 List of the task protocol used for each of the sessions. 

196 list of str 

197 List of ISO date strings where training was conducted on ephys rig. Empty list if all 

198 sessions on training rig. 

199 n_delay : int 

200 Number of sessions on ephys rig that had delay prior to starting session > 15min. 

201 Returns 0 if no sessions detected. 

202 """ 

203 one = one or ONE() 

204 

205 if date is None: 

206 # compute from yesterday 

207 specified_date = (datetime.date.today() - datetime.timedelta(days=1)) 

208 latest_sess = specified_date.strftime("%Y-%m-%d") 

209 latest_minus_week = (datetime.date.today() - 

210 datetime.timedelta(days=8)).strftime("%Y-%m-%d") 

211 else: 

212 # compute from the date specified 

213 specified_date = datetime.datetime.strptime(date, '%Y-%m-%d') 

214 latest_minus_week = (specified_date - datetime.timedelta(days=7)).strftime("%Y-%m-%d") 

215 latest_sess = date 

216 

217 sessions = one.alyx.rest('sessions', 'list', subject=subj, date_range=[latest_minus_week, 

218 latest_sess], dataset_types='trials.goCueTrigger_times') 

219 

220 # If not enough sessions in the last week, then just fetch them all 

221 if len(sessions) < 3: 

222 specified_date_plus = (specified_date + datetime.timedelta(days=1)).strftime("%Y-%m-%d") 

223 django_query = 'start_time__lte,' + specified_date_plus 

224 sessions = one.alyx.rest('sessions', 'list', subject=subj, 

225 dataset_types='trials.goCueTrigger_times', django=django_query) 

226 

227 # If still 0 sessions then return with warning 

228 if len(sessions) == 0: 

229 _logger.warning(f"No training sessions detected for {subj}") 

230 return [None] * 4 

231 

232 trials = Bunch() 

233 task_protocol = [] 

234 sess_dates = [] 

235 if len(sessions) < 3: 

236 for n, _ in enumerate(sessions): 

237 try: 

238 trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials') 

239 except ALFObjectNotFound: 

240 trials_ = None 

241 

242 if trials_: 

243 task_protocol.append(re.search('tasks_(.*)Choice', 

244 sessions[n]['task_protocol']).group(1)) 

245 sess_dates.append(sessions[n]['start_time'][:10]) 

246 trials[sessions[n]['start_time'][:10]] = trials_ 

247 

248 else: 

249 n = 0 

250 while len(trials) < 3: 

251 print(sessions[n]['url'].split('/')[-1]) 

252 try: 

253 trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials') 

254 except ALFObjectNotFound: 

255 trials_ = None 

256 

257 if trials_: 

258 task_protocol.append(re.search('tasks_(.*)Choice', 

259 sessions[n]['task_protocol']).group(1)) 

260 sess_dates.append(sessions[n]['start_time'][:10]) 

261 trials[sessions[n]['start_time'][:10]] = trials_ 

262 

263 n += 1 

264 

265 if not np.any(np.array(task_protocol) == 'training'): 

266 ephys_sess = one.alyx.rest('sessions', 'list', subject=subj, 

267 date_range=[sess_dates[-1], sess_dates[0]], 

268 django='location__name__icontains,ephys') 

269 if len(ephys_sess) > 0: 

270 ephys_sess_dates = [sess['start_time'][:10] for sess in ephys_sess] 

271 

272 n_delay = len(one.alyx.rest('sessions', 'list', subject=subj, 

273 date_range=[sess_dates[-1], sess_dates[0]], 

274 django='json__SESSION_DELAY_START__gte,900')) 

275 else: 

276 ephys_sess_dates = [] 

277 n_delay = 0 

278 else: 

279 ephys_sess_dates = [] 

280 n_delay = 0 

281 

282 return trials, task_protocol, ephys_sess_dates, n_delay 

283 

284 

285def get_training_status(trials, task_protocol, ephys_sess_dates, n_delay): (empty)

286 """ 

287 Compute training status of a subject from consecutive training datasets. 

288 

289 For IBL, training status is calculated using trials from the last three consecutive sessions. 

290 

291 Parameters 

292 ---------- 

293 trials : dict of str 

294 Dictionary of trials objects where each key is the ISO session date string. 

295 task_protocol : list of str 

296 Task protocol used for each training session in `trials`, can be 'training', 'biased' or 

297 'ephys'. 

298 ephys_sess_dates : list of str 

299 List of ISO date strings where training was conducted on ephys rig. Empty list if all 

300 sessions on training rig. 

301 n_delay : int 

302 Number of sessions on ephys rig that had delay prior to starting session > 15min. 

303 Returns 0 if no sessions detected. 

304 

305 Returns 

306 ------- 

307 str 

308 Training status of the subject. 

309 iblutil.util.Bunch 

310 Bunch containing performance metrics that decide training status i.e. performance on easy 

311 trials, number of trials, psychometric fit parameters, reaction time. 

312 """ 

313 

314 info = Bunch() 9 ctx1edcgfhjba

315 trials_all = concatenate_trials(trials) 9 ctx1edcgfhjba

316 info.session_dates = list(trials.keys()) 9 ctx1edcgfhjba

317 info.protocols = [p for p in task_protocol] 9 ctx1edcgfhjba

318 

319 # Case when all sessions are trainingChoiceWorld 

320 if np.all(np.array(task_protocol) == 'training'): 9 ctx1edcgfhjba

321 signed_contrast = np.unique(get_signed_contrast(trials_all)) 4 ctx1efha

322 (info.perf_easy, info.n_trials, 4 ctx1efha

323 info.psych, info.rt) = compute_training_info(trials, trials_all) 

324 

325 pass_criteria, criteria = criterion_1b(info.psych, info.n_trials, info.perf_easy, info.rt, 4 ctx1efha

326 signed_contrast) 

327 if pass_criteria: 4 ctx1efha

328 failed_criteria = Bunch() 2 ctx1ha

329 failed_criteria['NBiased'] = {'val': info.protocols, 'pass': False} 2 ctx1ha

330 failed_criteria['Criteria'] = {'val': 'ready4ephysrig', 'pass': False} 2 ctx1ha

331 status = 'trained 1b' 2 ctx1ha

332 else: 

333 failed_criteria = criteria 3 ctx1efa

334 pass_criteria, criteria = criterion_1a(info.psych, info.n_trials, info.perf_easy, signed_contrast) 3 ctx1efa

335 if pass_criteria: 3 ctx1efa

336 status = 'trained 1a' 2 ctx1fa

337 else: 

338 failed_criteria = criteria 2 ctx1ea

339 status = 'in training' 2 ctx1ea

340 

341 return status, info, failed_criteria 4 ctx1efha

342 

343 # Case when there are < 3 biasedChoiceWorld sessions after reaching trained_1b criterion 

344 if ~np.all(np.array(task_protocol) == 'training') and \ 6 ctx1dcgjba

345 np.any(np.array(task_protocol) == 'training'): 

346 status = 'trained 1b' 2 ctx1ja

347 (info.perf_easy, info.n_trials, 2 ctx1ja

348 info.psych, info.rt) = compute_training_info(trials, trials_all) 

349 

350 criteria = Bunch() 2 ctx1ja

351 criteria['NBiased'] = {'val': info.protocols, 'pass': False} 2 ctx1ja

352 criteria['Criteria'] = {'val': 'ready4ephysrig', 'pass': False} 2 ctx1ja

353 

354 return status, info, criteria 2 ctx1ja

355 

356 # Case when there is biasedChoiceWorld or ephysChoiceWorld in last three sessions 

357 if not np.any(np.array(task_protocol) == 'training'): 5 ctx1dcgba

358 

359 (info.perf_easy, info.n_trials, 5 ctx1dcgba

360 info.psych_20, info.psych_80, 

361 info.rt) = compute_bias_info(trials, trials_all) 

362 

363 n_ephys = len(ephys_sess_dates) 5 ctx1dcgba

364 info.n_ephys = n_ephys 5 ctx1dcgba

365 info.n_delay = n_delay 5 ctx1dcgba

366 

367 # Criterion recording 

368 pass_criteria, criteria = criteria_recording(n_ephys, n_delay, info.psych_20, info.psych_80, info.n_trials, 5 ctx1dcgba

369 info.perf_easy, info.rt) 

370 if pass_criteria: 5 ctx1dcgba

371 # Here the criteria doesn't actually fail but we have no other criteria to meet so we return this 

372 failed_criteria = criteria 2 ctx1ga

373 status = 'ready4recording' 2 ctx1ga

374 else: 

375 failed_criteria = criteria 4 ctx1dcba

376 assert all(date in trials for date in ephys_sess_dates) 4 ctx1dcba

377 perf_ephys_easy = np.array([compute_performance_easy(trials[k]) for k in 4 ctx1dcba

378 ephys_sess_dates]) 

379 n_ephys_trials = np.array([compute_n_trials(trials[k]) for k in ephys_sess_dates]) 4 ctx1dcba

380 

381 pass_criteria, criteria = criterion_delay(n_ephys_trials, perf_ephys_easy, n_ephys=n_ephys) 4 ctx1dcba

382 

383 if pass_criteria: 4 ctx1dcba

384 status = 'ready4delay' 2 ctx1da

385 else: 

386 failed_criteria = criteria 3 ctx1cba

387 pass_criteria, criteria = criterion_ephys(info.psych_20, info.psych_80, info.n_trials, 3 ctx1cba

388 info.perf_easy, info.rt) 

389 if pass_criteria: 3 ctx1cba

390 status = 'ready4ephysrig' 2 ctx1ca

391 else: 

392 failed_criteria = criteria 2 ctx1ba

393 status = 'trained 1b' 2 ctx1ba

394 

395 return status, info, failed_criteria 5 ctx1dcgba

396 

397 

398def display_status(subj, sess_dates, status, perf_easy=None, n_trials=None, psych=None, (empty)

399 psych_20=None, psych_80=None, rt=None): 

400 """ 

401 Display training status of subject to terminal. 

402 

403 Parameters 

404 ---------- 

405 subj : str 

406 Subject nickname (must match the name registered on Alyx). 

407 sess_dates : list of str 

408 ISO date strings of training sessions used to determine training status. 

409 status : str 

410 Training status of subject. 

411 perf_easy : numpy.array 

412 Proportion of correct high contrast trials for each training session. 

413 n_trials : numpy.array 

414 Total number of trials for each training session. 

415 psych : numpy.array 

416 Psychometric parameters fit to data from all training sessions - bias, threshold, lapse 

417 high, lapse low. 

418 psych_20 : numpy.array 

419 The fit psychometric parameters for the blocks where probability of a left stimulus is 0.2. 

420 psych_80 : numpy.array 

421 The fit psychometric parameters for the blocks where probability of a left stimulus is 0.8. 

422 rt : float 

423 The median response time for zero contrast trials across all training sessions. NaN 

424 indicates no zero contrast stimuli in training sessions. 

425 

426 """ 

427 

428 if perf_easy is None: 

429 print(f"\n{subj} : {status} \nSession dates=[{sess_dates[0]}, {sess_dates[1]}, " 

430 f"{sess_dates[2]}]") 

431 elif psych_20 is None: 

432 print(f"\n{subj} : {status} \nSession dates={[x for x in sess_dates]}, " 

433 f"Perf easy={[np.around(pe, 2) for pe in perf_easy]}, " 

434 f"N trials={[nt for nt in n_trials]} " 

435 f"\nPsych fit over last 3 sessions: " 

436 f"bias={np.around(psych[0], 2)}, thres={np.around(psych[1], 2)}, " 

437 f"lapse_low={np.around(psych[2], 2)}, lapse_high={np.around(psych[3], 2)} " 

438 f"\nMedian reaction time at 0 contrast over last 3 sessions = " 

439 f"{np.around(rt, 2)}") 

440 

441 else: 

442 print(f"\n{subj} : {status} \nSession dates={[x for x in sess_dates]}, " 

443 f"Perf easy={[np.around(pe, 2) for pe in perf_easy]}, " 

444 f"N trials={[nt for nt in n_trials]} " 

445 f"\nPsych fit over last 3 sessions (20): " 

446 f"bias={np.around(psych_20[0], 2)}, thres={np.around(psych_20[1], 2)}, " 

447 f"lapse_low={np.around(psych_20[2], 2)}, lapse_high={np.around(psych_20[3], 2)} " 

448 f"\nPsych fit over last 3 sessions (80): bias={np.around(psych_80[0], 2)}, " 

449 f"thres={np.around(psych_80[1], 2)}, lapse_low={np.around(psych_80[2], 2)}, " 

450 f"lapse_high={np.around(psych_80[3], 2)} " 

451 f"\nMedian reaction time at 0 contrast over last 3 sessions = " 

452 f"{np.around(rt, 2)}") 

453 

454 

455def concatenate_trials(trials): (empty)

456 """ 

457 Concatenate trials from different training sessions. 

458 

459 Parameters 

460 ---------- 

461 trials : dict of str 

462 Dictionary of trials objects where each key is the ISO session date string. 

463 

464 Returns 

465 ------- 

466 one.alf.io.AlfBunch 

467 Trials object with data concatenated over three training sessions. 

468 """ 

469 trials_all = AlfBunch() 12 ctx1kemdcgfhjbna

470 for k in TRIALS_KEYS: 12 ctx1kemdcgfhjbna

471 trials_all[k] = np.concatenate(list(trials[kk][k] for kk in trials.keys())) 12 ctx1kemdcgfhjbna

472 

473 return trials_all 12 ctx1kemdcgfhjbna

474 

475 

476def compute_training_info(trials, trials_all): (empty)

477 """ 

478 Compute all relevant performance metrics for when subject is on trainingChoiceWorld. 

479 

480 Parameters 

481 ---------- 

482 trials : dict of str 

483 Dictionary of trials objects where each key is the ISO session date string. 

484 trials_all : one.alf.io.AlfBunch 

485 Trials object with data concatenated over three training sessions. 

486 

487 Returns 

488 ------- 

489 numpy.array 

490 Proportion of correct high contrast trials for each session. 

491 numpy.array 

492 Total number of trials for each training session. 

493 numpy.array 

494 Array of psychometric parameters fit to `all_trials` - bias, threshold, lapse high, 

495 lapse low. 

496 float 

497 The median response time for all zero-contrast trials across all sessions. Returns NaN if 

498 no trials zero-contrast trials). 

499 """ 

500 

501 signed_contrast = get_signed_contrast(trials_all) 5 ctx1efhja

502 perf_easy = np.array([compute_performance_easy(trials[k]) for k in trials.keys()]) 5 ctx1efhja

503 n_trials = np.array([compute_n_trials(trials[k]) for k in trials.keys()]) 5 ctx1efhja

504 psych = compute_psychometric(trials_all, signed_contrast=signed_contrast) 5 ctx1efhja

505 rt = compute_median_reaction_time(trials_all, contrast=0, signed_contrast=signed_contrast) 5 ctx1efhja

506 

507 return perf_easy, n_trials, psych, rt 5 ctx1efhja

508 

509 

510def compute_bias_info(trials, trials_all): (empty)

511 """ 

512 Compute all relevant performance metrics for when subject is on biasedChoiceWorld 

513 

514 :param trials: dict containing trials objects from three consecutive training sessions, 

515 keys are session dates 

516 :type trials: Bunch 

517 :param trials_all: trials object with data concatenated over three training sessions 

518 :type trials_all: Bunch 

519 :returns: 

520 - perf_easy - performance of easy trials for each session 

521 - n_trials - number of trials in each session 

522 - psych_20 - parameters for psychometric curve fit to trials in 20 block over all sessions 

523 - psych_80 - parameters for psychometric curve fit to trials in 80 block over all sessions 

524 - rt - median reaction time for zero contrast stimuli over all sessions 

525 """ 

526 

527 signed_contrast = get_signed_contrast(trials_all) 5 ctx1dcgba

528 perf_easy = np.array([compute_performance_easy(trials[k]) for k in trials.keys()]) 5 ctx1dcgba

529 n_trials = np.array([compute_n_trials(trials[k]) for k in trials.keys()]) 5 ctx1dcgba

530 psych_20 = compute_psychometric(trials_all, signed_contrast=signed_contrast, block=0.2) 5 ctx1dcgba

531 psych_80 = compute_psychometric(trials_all, signed_contrast=signed_contrast, block=0.8) 5 ctx1dcgba

532 rt = compute_median_reaction_time(trials_all, contrast=0, signed_contrast=signed_contrast) 5 ctx1dcgba

533 

534 return perf_easy, n_trials, psych_20, psych_80, rt 5 ctx1dcgba

535 

536 

537def get_signed_contrast(trials): (empty)

538 """ 

539 Compute signed contrast from trials object 

540 

541 :param trials: trials object that must contain contrastLeft and contrastRight keys 

542 :type trials: dict 

543 returns: array of signed contrasts in percent, where -ve values are on the left 

544 """ 

545 # Replace NaNs with zeros, stack and take the difference 

546 contrast = np.nan_to_num(np.c_[trials['contrastLeft'], trials['contrastRight']]) 11 ctx1kemdcgfhjba

547 return np.diff(contrast).flatten() * 100 11 ctx1kemdcgfhjba

548 

549 

550def compute_performance_easy(trials): (empty)

551 """ 

552 Compute performance on easy trials (stimulus >= 50 %) from trials object 

553 

554 :param trials: trials object that must contain contrastLeft, contrastRight and feedbackType 

555 keys 

556 :type trials: dict 

557 returns: float containing performance on easy contrast trials 

558 """ 

559 signed_contrast = get_signed_contrast(trials) 10 ctx1kedcgfhjba

560 easy_trials = np.where(np.abs(signed_contrast) >= 50)[0] 10 ctx1kedcgfhjba

561 return np.sum(trials['feedbackType'][easy_trials] == 1) / easy_trials.shape[0] 10 ctx1kedcgfhjba

562 

563 

564def compute_performance(trials, signed_contrast=None, block=None, prob_right=False): (empty)

565 """ 

566 Compute performance on all trials at each contrast level from trials object 

567 

568 :param trials: trials object that must contain contrastLeft, contrastRight and feedbackType 

569 keys 

570 :type trials: dict 

571 returns: float containing performance on easy contrast trials 

572 """ 

573 if signed_contrast is None: 10 ctx1kedcgfhjba

574 signed_contrast = get_signed_contrast(trials) 2 ctx1ba

575 

576 if block is None: 10 ctx1kedcgfhjba

577 block_idx = np.full(trials.probabilityLeft.shape, True, dtype=bool) 7 ctx1kefhjba

578 else: 

579 block_idx = trials.probabilityLeft == block 5 ctx1dcgba

580 

581 if not np.any(block_idx): 10 ctx1kedcgfhjba

582 return np.nan * np.zeros(3) 1 ctx1b

583 

584 contrasts, n_contrasts = np.unique(signed_contrast[block_idx], return_counts=True) 10 ctx1kedcgfhjba

585 

586 if not prob_right: 10 ctx1kedcgfhjba

587 correct = trials.feedbackType == 1 

588 performance = np.vectorize(lambda x: np.mean(correct[(x == signed_contrast) & block_idx]))(contrasts) 

589 else: 

590 rightward = trials.choice == -1 10 ctx1kedcgfhjba

591 # Calculate the proportion rightward for each contrast type 

592 performance = np.vectorize(lambda x: np.mean(rightward[(x == signed_contrast) & block_idx]))(contrasts) 10 ctx1kedcgfhjba

593 

594 return performance, contrasts, n_contrasts 10 ctx1kedcgfhjba

595 

596 

597def compute_n_trials(trials): (empty)

598 """ 

599 Compute number of trials in trials object 

600 

601 :param trials: trials object 

602 :type trials: dict 

603 returns: int containing number of trials in session 

604 """ 

605 return trials['choice'].shape[0] 10 ctx1kedcgfhjba

606 

607 

608def compute_psychometric(trials, signed_contrast=None, block=None, plotting=False, compute_ci=False, alpha=.032): (empty)

609 """ 

610 Compute psychometric fit parameters for trials object. 

611 

612 Parameters 

613 ---------- 

614 trials : one.alf.io.AlfBunch 

615 An ALF trials object containing the keys {'probabilityLeft', 'contrastLeft', 

616 'contrastRight', 'feedbackType', 'choice', 'response_times', 'stimOn_times'}. 

617 signed_contrast : numpy.array 

618 An array of signed contrasts in percent the length of trials, where left contrasts are -ve. 

619 If None, these are computed from the trials object. 

620 block : float 

621 The block type to compute. If None, all trials are included, otherwise only trials where 

622 probabilityLeft matches this value are included. For biasedChoiceWorld, the 

623 probabilityLeft set is {0.5, 0.2, 0.8}. 

624 plotting : bool 

625 Which set of psychofit model parameters to use (see notes). 

626 compute_ci : bool 

627 If true, computes and returns the confidence intervals for response at each contrast. 

628 alpha : float, default=0.032 

629 Significance level for confidence interval. Must be in (0, 1). If `compute_ci` is false, 

630 this value is ignored. 

631 

632 Returns 

633 ------- 

634 numpy.array 

635 Array of psychometric fit parameters - bias, threshold, lapse high, lapse low. 

636 (tuple of numpy.array) 

637 If `compute_ci` is true, a tuple of 

638 

639 See Also 

640 -------- 

641 statsmodels.stats.proportion.proportion_confint - The function used to compute confidence 

642 interval. 

643 psychofit.mle_fit_psycho - The function used to fit the psychometric parameters. 

644 

645 Notes 

646 ----- 

647 The psychofit starting parameters and model constraints used for the fit when computing the 

648 training status (e.g. trained_1a, etc.) are sub-optimal and can produce a poor fit. To keep 

649 the precise criteria the same for all subjects, these parameters have not changed. To produce a 

650 better fit for plotting purposes, or to calculate the training status in a manner inconsistent 

651 with the IBL training pipeline, use plotting=True. 

652 """ 

653 

654 if signed_contrast is None: 11 ctx1kemdcgfhjba

655 signed_contrast = get_signed_contrast(trials) 4 ctx1kmba

656 

657 if block is None: 11 ctx1kemdcgfhjba

658 block_idx = np.full(trials.probabilityLeft.shape, True, dtype=bool) 6 ctx1kefhja

659 else: 

660 block_idx = trials.probabilityLeft == block 6 ctx1mdcgba

661 

662 if not np.any(block_idx): 11 ctx1kemdcgfhjba

663 return np.nan * np.zeros(4) 3 ctx1mba

664 

665 prob_choose_right, contrasts, n_contrasts = compute_performance( 10 ctx1kedcgfhjba

666 trials, signed_contrast=signed_contrast, block=block, prob_right=True) 

667 

668 if plotting: 10 ctx1kedcgfhjba

669 # These starting parameters and constraints tend to produce a better fit, and are therefore 

670 # used for plotting. 

671 psych, _ = psy.mle_fit_psycho( 1 ctx1b

672 np.vstack([contrasts, n_contrasts, prob_choose_right]), 

673 P_model='erf_psycho_2gammas', 

674 parstart=np.array([0., 40., 0.1, 0.1]), 

675 parmin=np.array([-50., 10., 0., 0.]), 

676 parmax=np.array([50., 50., 0.2, 0.2]), 

677 nfits=10) 

678 else: 

679 # These starting parameters and constraints are not ideal but are still used for computing 

680 # the training status for consistency. 

681 psych, _ = psy.mle_fit_psycho( 10 ctx1kedcgfhjba

682 np.vstack([contrasts, n_contrasts, prob_choose_right]), 

683 P_model='erf_psycho_2gammas', 

684 parstart=np.array([np.mean(contrasts), 20., 0.05, 0.05]), 

685 parmin=np.array([np.min(contrasts), 0., 0., 0.]), 

686 parmax=np.array([np.max(contrasts), 100., 1, 1])) 

687 

688 if compute_ci: 10 ctx1kedcgfhjba

689 import statsmodels.stats.proportion as smp # noqa 

690 # choice == -1 means contrast on right hand side 

691 n_right = np.vectorize(lambda x: np.sum(trials['choice'][(x == signed_contrast) & block_idx] == -1))(contrasts) 

692 ci = smp.proportion_confint(n_right, n_contrasts, alpha=alpha, method='normal') - prob_choose_right 

693 

694 return psych, ci 

695 else: 

696 return psych 10 ctx1kedcgfhjba

697 

698 

699def compute_median_reaction_time(trials, stim_on_type='stimOn_times', contrast=None, signed_contrast=None): (empty)

700 """ 

701 Compute median response time on zero contrast trials from trials object 

702 

703 Parameters 

704 ---------- 

705 trials : one.alf.io.AlfBunch 

706 An ALF trials object containing the keys {'probabilityLeft', 'contrastLeft', 

707 'contrastRight', 'feedbackType', 'choice', 'response_times', 'stimOn_times'}. 

708 stim_on_type : str, default='stimOn_times' 

709 The trials key to use when calculating the response times. The difference between this and 

710 'feedback_times' is used (see notes). 

711 contrast : float 

712 If None, the median response time is calculated for all trials, regardless of contrast, 

713 otherwise only trials where the matching signed percent contrast was presented are used. 

714 signed_contrast : numpy.array 

715 An array of signed contrasts in percent the length of trials, where left contrasts are -ve. 

716 If None, these are computed from the trials object. 

717 

718 Returns 

719 ------- 

720 float 

721 The median response time for trials with `contrast` (returns NaN if no trials matching 

722 `contrast` in trials object). 

723 

724 Notes 

725 ----- 

726 - The `stim_on_type` is 'stimOn_times' by default, however for IBL rig data, the photodiode is 

727 sometimes not calibrated properly which can lead to inaccurate (or absent, i.e. NaN) stim on 

728 times. Therefore, it is sometimes more accurate to use the 'stimOnTrigger_times' (the time of 

729 the stimulus onset command), if available, or the 'goCue_times' (the time of the soundcard 

730 output TTL when the audio go cue is played) or the 'goCueTrigger_times' (the time of the 

731 audio go cue command). 

732 - The response/reaction time here is defined as the time between stim on and feedback, i.e. the 

733 entire open-loop trial duration. 

734 """ 

735 if signed_contrast is None: 10 ctx1kedcgfhjba

736 signed_contrast = get_signed_contrast(trials) 3 ctx1kba

737 

738 if contrast is None: 10 ctx1kedcgfhjba

739 contrast_idx = np.full(trials.probabilityLeft.shape, True, dtype=bool) 2 ctx1ba

740 else: 

741 contrast_idx = signed_contrast == contrast 10 ctx1kedcgfhjba

742 

743 if np.any(contrast_idx): 10 ctx1kedcgfhjba

744 reaction_time = np.nanmedian((trials.response_times - trials[stim_on_type]) 10 ctx1kedcgfhjba

745 [contrast_idx]) 

746 else: 

747 reaction_time = np.nan 1 ctx1a

748 

749 return reaction_time 10 ctx1kedcgfhjba

750 

751 

752def compute_reaction_time(trials, stim_on_type='stimOn_times', stim_off_type='response_times', signed_contrast=None, block=None, (empty)

753 compute_ci=False, alpha=0.32): 

754 """ 

755 Compute median response time for all contrasts. 

756 

757 Parameters 

758 ---------- 

759 trials : one.alf.io.AlfBunch 

760 An ALF trials object containing the keys {'probabilityLeft', 'contrastLeft', 

761 'contrastRight', 'feedbackType', 'choice', 'response_times', 'stimOn_times'}. 

762 stim_on_type : str, default='stimOn_times' 

763 The trials key to use when calculating the response times. The difference between this and 

764 `stim_off_type` is used (see notes). 

765 stim_off_type : str, default='response_times' 

766 The trials key to use when calculating the response times. The difference between this and 

767 `stim_on_type` is used (see notes). 

768 signed_contrast : numpy.array 

769 An array of signed contrasts in percent the length of trials, where left contrasts are -ve. 

770 If None, these are computed from the trials object. 

771 block : float 

772 The block type to compute. If None, all trials are included, otherwise only trials where 

773 probabilityLeft matches this value are included. For biasedChoiceWorld, the 

774 probabilityLeft set is {0.5, 0.2, 0.8}. 

775 compute_ci : bool 

776 If true, computes and returns the confidence intervals for response time at each contrast. 

777 alpha : float, default=0.32 

778 Significance level for confidence interval. Must be in (0, 1). If `compute_ci` is false, 

779 this value is ignored. 

780 

781 Returns 

782 ------- 

783 numpy.array 

784 The median response times for each unique signed contrast. 

785 numpy.array 

786 The set of unique signed contrasts. 

787 numpy.array 

788 The number of trials for each unique signed contrast. 

789 (numpy.array) 

790 If `compute_ci` is true, an array of confidence intervals is return in the shape (n_trials, 

791 2). 

792 

793 Notes 

794 ----- 

795 - The response/reaction time by default is the time between stim on and response, i.e. the 

796 entire open-loop trial duration. One could use 'stimOn_times' and 'firstMovement_times' to 

797 get the true reaction time, or 'firstMovement_times' and 'response_times' to get the true 

798 response times, or calculate the last movement onset times and calculate the true movement 

799 times. See module examples for how to calculate this. 

800 

801 See Also 

802 -------- 

803 scipy.stats.bootstrap - the function used to compute the confidence interval. 

804 """ 

805 

806 if signed_contrast is None: 

807 signed_contrast = get_signed_contrast(trials) 

808 

809 if block is None: 

810 block_idx = np.full(trials.probabilityLeft.shape, True, dtype=bool) 

811 else: 

812 block_idx = trials.probabilityLeft == block 

813 

814 contrasts, n_contrasts = np.unique(signed_contrast[block_idx], return_counts=True) 

815 reaction_time = np.vectorize( 

816 lambda x: np.nanmedian((trials[stim_off_type] - trials[stim_on_type])[(x == signed_contrast) & block_idx]), 

817 otypes=[float] 

818 )(contrasts) 

819 

820 if compute_ci: 

821 ci = np.full((contrasts.size, 2), np.nan) 

822 for i, x in enumerate(contrasts): 

823 data = (trials[stim_off_type] - trials[stim_on_type])[(x == signed_contrast) & block_idx] 

824 bt = bootstrap((data,), np.nanmedian, confidence_level=1 - alpha) 

825 ci[i, 0] = bt.confidence_interval.low 

826 ci[i, 1] = bt.confidence_interval.high 

827 

828 return reaction_time, contrasts, n_contrasts, ci 

829 else: 

830 return reaction_time, contrasts, n_contrasts, 

831 

832 

833def criterion_1a(psych, n_trials, perf_easy, signed_contrast): (empty)

834 """ 

835 Returns bool indicating whether criteria for status 'trained_1a' are met. 

836 

837 Criteria 

838 -------- 

839 - Bias is less than 16 

840 - Threshold is less than 19 

841 - Lapse rate on both sides is less than 0.2 

842 - The total number of trials is greater than 200 for each session 

843 - Performance on easy contrasts > 80% for all sessions 

844 - Zero contrast trials must be present 

845 

846 Parameters 

847 ---------- 

848 psych : numpy.array 

849 The fit psychometric parameters three consecutive sessions. Parameters are bias, threshold, 

850 lapse high, lapse low. 

851 n_trials : numpy.array of int 

852 The number for trials for each session. 

853 perf_easy : numpy.array of float 

854 The proportion of correct high contrast trials for each session. 

855 signed_contrast: numpy.array 

856 Unique list of contrasts displayed 

857 

858 Returns 

859 ------- 

860 bool 

861 True if the criteria are met for 'trained_1a'. 

862 Bunch 

863 Bunch containing breakdown of the passing/ failing critieria 

864 

865 Notes 

866 ----- 

867 The parameter thresholds chosen here were originally determined by averaging the parameter fits 

868 for a number of sessions determined to be of 'good' performance by an experimenter.