Coverage for ibllib/qc/task_metrics.py: 97%
404 statements
« prev ^ index » next coverage.py v7.7.0, created at 2025-03-17 15:25 +0000
1"""Behaviour QC.
3This module runs a list of quality control metrics on the behaviour data.
5.. warning::
6 The QC should be loaded using :meth:`ibllib.pipes.base_tasks.BehaviourTask.run_qc` and not
7 instantiated directly.
9Examples
10--------
11Running on a behaviour rig computer and updating QC fields in Alyx:
13>>> from ibllib.qc.task_qc_viewer.task_qc import show_session_task_qc
14>>> qc = show_session_task_qc(session_path, bpod_only=True, local=True) # must close Viewer window
15>>> qc = qc.run(update=True)
17Downloading the required data and inspecting the QC on a different computer:
19>>> from ibllib.pipes.dynamic_pipeline import get_trials_tasks
20>>> from one.api import ONE
21>>> task = get_trials_tasks(session_path, one=ONE())[0] # get first task run
22>>> task.location = 'remote'
23>>> task.setUp() # download required data
24>>> qc = task.run_qc(update=False)
25>>> outcome, results = qc.run()
27Inspecting individual test outcomes
29>>> outcome, results, outcomes = qc.compute_session_status()
31Running bpod QC on ephys session (when not on behaviour rig PC)
33>>> from ibllib.qc.task_qc_viewer.task_qc import get_bpod_trials_task, get_trials_tasks
34>>> from one.api import ONE
35>>> tasks = get_trials_tasks(session_path, one=ONE())
36>>> task = get_bpod_trials_task(tasks[0]) # Ensure Bpod only on behaviour rig
37>>> task.location = 'remote'
38>>> task.setUp() # download required data
39>>> qc = task.run_qc(update=False)
40>>> outcome, results = qc.run()
42Running ephys QC, from local server PC (after ephys + bpod data have been copied to a same folder)
44>>> from ibllib.pipes.dynamic_pipeline import get_trials_tasks
45>>> task = get_trials_tasks(session_path, one=ONE())[0] # get first task run
46>>> qc = task.run_qc(update=False)
47>>> outcome, results = qc.run()
48"""
49import logging
50import sys
51from packaging import version
52from pathlib import Path, PurePosixPath
53from datetime import datetime, timedelta
54from inspect import getmembers, isfunction
55from functools import reduce
56from collections.abc import Sized
58import numpy as np
59from scipy.stats import chisquare
61from brainbox.behavior.wheel import cm_to_rad, traces_by_trial
62from ibllib.io.extractors import ephys_fpga
63from one.alf import spec
64from . import base
66_log = logging.getLogger(__name__)
68# todo the 2 followint parameters should be read from the task parameters for each session
69ITI_DELAY_SECS = .5
70FEEDBACK_NOGO_DELAY_SECS = 2
72BWM_CRITERIA = {
73 'default': {'PASS': 0.99, 'WARNING': 0.90, 'FAIL': 0}, # Note: WARNING was 0.95 prior to Aug 2022
74 '_task_stimOff_itiIn_delays': {'PASS': 0.99, 'WARNING': 0},
75 '_task_positive_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
76 '_task_negative_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
77 '_task_wheel_move_during_closed_loop': {'PASS': 0.99, 'WARNING': 0},
78 '_task_response_stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0},
79 '_task_detected_wheel_moves': {'PASS': 0.99, 'WARNING': 0},
80 '_task_trial_length': {'PASS': 0.99, 'WARNING': 0},
81 '_task_goCue_delays': {'PASS': 0.99, 'WARNING': 0},
82 '_task_errorCue_delays': {'PASS': 0.99, 'WARNING': 0},
83 '_task_stimOn_delays': {'PASS': 0.99, 'WARNING': 0},
84 '_task_stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
85 '_task_stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0},
86 '_task_iti_delays': {'NOT_SET': 0},
87 '_task_passed_trial_checks': {'NOT_SET': 0}
88}
91def compute_session_status_from_dict(results, criteria=None):
92 """
93 Compute overall task QC value from QC check results.
95 Given a dictionary of results, computes the overall session QC for each key and aggregates
96 in a single value.
98 Parameters
99 ----------
100 results : dict
101 A dictionary of QC keys containing (usually scalar) values.
102 criteria : dict
103 A dictionary of qc keys containing map of PASS, WARNING, FAIL thresholds.
105 Returns
106 -------
107 one.alf.spec.QC
108 Overall session QC outcome.
109 dict
110 A map of QC tests and their outcomes.
111 """
112 if not criteria: 1axmlngcbdei
113 criteria = {'default': BWM_CRITERIA['default']}
114 outcomes = {k: TaskQC.thresholding(v, thresholds=criteria.get(k, criteria['default'])) 1axmlngcbdei
115 for k, v in results.items()}
117 # Criteria map is in order of severity so the max index is our overall QC outcome
118 session_outcome = base.QC.overall_outcome(outcomes.values()) 1amlngcbdei
119 return session_outcome, outcomes 1amlngcbdei
122def update_dataset_qc(qc, registered_datasets, one, override=False):
123 """
124 Update QC values for individual datasets.
126 Parameters
127 ----------
128 qc : ibllib.qc.task_metrics.TaskQC
129 A TaskQC object that has been run.
130 registered_datasets : list of dict
131 A list of Alyx dataset records.
132 one : one.api.OneAlyx
133 An online instance of ONE.
134 override : bool
135 If True the QC field is updated even if new value is better than previous.
137 Returns
138 -------
139 list of dict
140 The list of registered datasets but with the 'qc' fields updated.
141 """
142 # Create map of dataset name, sans extension, to dataset id
143 stem2id = {PurePosixPath(dset['name']).stem: dset.get('id') for dset in registered_datasets} 1hc
144 # Ensure dataset stems are unique
145 assert len(stem2id) == len(registered_datasets), 'ambiguous dataset names' 1hc
147 # dict of QC check to outcome (as enum value)
148 *_, outcomes = qc.compute_session_status() 1hc
149 # work over map of dataset name (sans extension) to outcome (enum or dict of columns: enum)
150 for name, outcome in qc.compute_dataset_qc_status(outcomes).items(): 1hc
151 # if outcome is a dict, calculate aggregate outcome for each column
152 if isinstance(outcome, dict): 1hc
153 extended_qc = outcome 1hc
154 outcome = qc.overall_outcome(outcome.values()) 1hc
155 else:
156 extended_qc = {} 1hc
157 # check if dataset was registered to Alyx
158 if not (did := stem2id.get(name)): 1hc
159 _log.debug('dataset %s not registered, skipping', name) 1h
160 continue 1h
161 # update the dataset QC value on Alyx
162 if outcome > spec.QC.NOT_SET or override: 1hc
163 dset_qc = base.QC(did, one=one, log=_log, endpoint='datasets') 1hc
164 dset = next(x for x in registered_datasets if did == x.get('id')) 1hc
165 dset['qc'] = dset_qc.update(outcome, namespace='', override=override).name 1hc
166 if extended_qc: 1hc
167 dset_qc.update_extended_qc(extended_qc) 1hc
168 return registered_datasets 1hc
171class TaskQC(base.QC):
172 """Task QC for training, biased, and ephys choice world."""
174 criteria = BWM_CRITERIA
176 extractor = None
177 """ibllib.qc.task_extractors.TaskQCExtractor: A task extractor object containing raw and extracted data."""
179 @staticmethod
180 def thresholding(qc_value, thresholds=None) -> spec.QC:
181 """
182 Compute the outcome of a single key by applying thresholding.
184 Parameters
185 ----------
186 qc_value : float
187 Proportion of passing qcs, between 0 and 1.
188 thresholds : dict
189 Dictionary with keys 'PASS', 'WARNING', 'FAIL', (or enum
190 integers, c.f. one.alf.spec.QC).
192 Returns
193 -------
194 one.alf.spec.QC
195 The outcome.
196 """
197 thresholds = {spec.QC.validate(k): v for k, v in thresholds.items() or {}} 1axmlngcbdei
198 MAX_BOUND, MIN_BOUND = (1, 0) 1axmlngcbdei
199 if qc_value is None or np.isnan(qc_value): 1axmlngcbdei
200 return spec.QC.NOT_SET 1amlngbei
201 elif (qc_value > MAX_BOUND) or (qc_value < MIN_BOUND): 1axmlngcbdei
202 raise ValueError('Values out of bound') 1x
203 for crit in filter(None, sorted(spec.QC)): 1amlngcbdei
204 if crit in thresholds.keys() and qc_value >= thresholds[crit]: 1amlngcbdei
205 return crit 1amlngcbdei
206 # if None of this applies, return 'NOT_SET'
207 return spec.QC.NOT_SET 1algcbdei
209 def __init__(self, session_path_or_eid, **kwargs):
210 """
211 Task QC for training, biased, and ephys choice world.
213 :param session_path_or_eid: A session eid or path
214 :param log: A logging.Logger instance, if None the 'ibllib' logger is used
215 :param one: An ONE instance for fetching and setting the QC on Alyx
216 """
217 # When an eid is provided, we will download the required data by default (if necessary)
218 self.download_data = not spec.is_session_path(Path(session_path_or_eid)) 1ahcbd
219 super().__init__(session_path_or_eid, **kwargs) 1ahcbd
221 # Data
222 self.extractor = None 1ahcbd
224 # Metrics and passed trials
225 self.metrics = None 1ahcbd
226 self.passed = None 1ahcbd
228 # Criteria (initialize as outcomes vary by class, task, and hardware)
229 self.criteria = BWM_CRITERIA.copy() 1ahcbd
231 def compute(self, **kwargs):
232 """Compute and store the QC metrics.
234 Runs the QC on the session and stores a map of the metrics for each datapoint for each
235 test, and a map of which datapoints passed for each test.
237 Parameters
238 ----------
239 bpod_only : bool
240 If True no data is extracted from the FPGA for ephys sessions.
241 """
242 assert self.extractor is not None 1acbdfe
244 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1acbdfe
245 if version.parse(ver) >= version.parse('8.0.0'): 1acbdfe
246 self.criteria['_task_iti_delays'] = {'PASS': 0.99, 'WARNING': 0}
247 self.criteria['_task_passed_trial_checks'] = {'PASS': 0.7, 'WARNING': 0}
248 else:
249 self.criteria['_task_iti_delays'] = {'NOT_SET': 0} 1acbdfe
250 self.criteria['_task_passed_trial_checks'] = {'NOT_SET': 0} 1acbdfe
252 self.log.info(f'Session {self.session_path}: Running QC on behavior data...') 1acbdfe
253 self.get_bpodqc_metrics_frame( 1acbdfe
254 self.extractor.data,
255 wheel_gain=self.extractor.settings['STIM_GAIN'], # The wheel gain
256 photodiode=self.extractor.frame_ttls,
257 audio=self.extractor.audio_ttls,
258 re_encoding=self.extractor.wheel_encoding or 'X1',
259 min_qt=self.extractor.settings.get('QUIESCENT_PERIOD') or 0.2,
260 audio_output=self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown')
261 )
263 def _get_checks(self):
264 """
265 Find all methods that begin with 'check_'.
267 Returns
268 -------
269 Dict[str, function]
270 A map of QC check function names and the corresponding functions that return `metric`
271 (any), `passed` (bool).
272 """
273 def is_metric(x): 1athcbdfe
274 return isfunction(x) and x.__name__.startswith('check_') 1athcbdfe
276 return dict(getmembers(sys.modules[__name__], is_metric)) 1athcbdfe
278 def get_bpodqc_metrics_frame(self, data, **kwargs):
279 """
280 Evaluate task QC metrics.
282 Evaluates all the QC metric functions in this module (those starting with 'check') and
283 returns the results. The optional kwargs listed below are passed to each QC metric function.
285 Parameters
286 ----------
287 data : dict
288 The extracted task data.
289 re_encoding : str {'X1', 'X2', 'X4'}
290 The encoding configuration of the rotary encoder.
291 enc_res : int
292 The rotary encoder resolution as number of fronts per revolution.
293 wheel_gain : float
294 The STIM_GAIN task parameter.
295 photodiode : dict
296 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel.
297 audio : dict
298 The fronts from Bpod's BNC2 input FPGA audio sync channel.
299 min_qt : float
300 The QUIESCENT_PERIOD task parameter.
302 Returns
303 -------
304 dict
305 Map of checks and their QC metric values (1 per trial).
306 dict
307 Map of checks and a float array of which samples passed.
308 """
309 # Find all methods that begin with 'check_'
310 checks = self._get_checks() 1acbdfe
311 prefix = '_task_' # Extended QC fields will start with this 1acbdfe
312 # Method 'check_foobar' stored with key '_task_foobar' in metrics map
313 qc_metrics_map = {prefix + k[6:]: fn(data, **kwargs) for k, fn in checks.items()} 1acbdfe
315 # Split metrics and passed frames
316 self.metrics = {} 1acbdfe
317 self.passed = {} 1acbdfe
318 for k in qc_metrics_map: 1acbdfe
319 self.metrics[k], self.passed[k] = qc_metrics_map[k] 1acbdfe
321 # Add a check for trial level pass: did a given trial pass all checks?
322 n_trials = data['intervals'].shape[0] 1acbdfe
323 # Trial-level checks return an array the length that equals the number of trials
324 trial_level_passed = [m for m in self.passed.values() if isinstance(m, Sized) and len(m) == n_trials] 1acbdfe
325 name = prefix + 'passed_trial_checks' 1acbdfe
326 self.metrics[name] = reduce(np.logical_and, trial_level_passed or (None, None)) 1acbdfe
327 self.passed[name] = self.metrics[name].astype(float) if trial_level_passed else None 1acbdfe
329 def run(self, update=False, namespace='task', **kwargs):
330 """
331 Compute the QC outcomes and return overall task QC outcome.
333 Parameters
334 ----------
335 update : bool
336 If True, updates the session QC fields on Alyx.
337 namespace : str
338 The namespace of the QC fields in the Alyx JSON field.
339 bpod_only : bool
340 If True no data is extracted from the FPGA for ephys sessions.
342 Returns
343 -------
344 str
345 Overall task QC outcome.
346 dict
347 A map of QC tests and the proportion of data points that passed them.
348 """
349 if self.metrics is None: 1acbdi
350 self.compute(**kwargs) 1acbd
351 outcome, results, _ = self.compute_session_status() 1acbdi
352 if update: 1acbdi
353 self.update_extended_qc(results) 1ci
354 self.update(outcome, namespace) 1ci
355 return outcome, results 1acbdi
357 def compute_session_status(self):
358 """
359 Compute the overall session QC for each key and aggregates in a single value.
361 Returns
362 -------
363 str
364 Overall session QC outcome.
365 dict
366 A map of QC tests and the proportion of data points that passed them.
367 dict
368 A map of QC tests and their outcomes.
369 """
370 if self.passed is None: 1agcbdei
371 raise AttributeError('passed is None; compute QC first') 1e
372 # Get mean passed of each check, or None if passed is None or all NaN
373 results = {k: None if v is None or np.isnan(v).all() else np.nanmean(v) 1agcbdei
374 for k, v in self.passed.items()}
375 session_outcome, outcomes = compute_session_status_from_dict(results, self.criteria) 1agcbdei
376 return session_outcome, results, outcomes 1agcbdei
378 @staticmethod
379 def compute_dataset_qc_status(outcomes):
380 """Return map of dataset specific QC values.
382 Parameters
383 ----------
384 outcomes : dict
385 Map of checks and their individual outcomes.
387 Returns
388 -------
389 dict
390 Map of dataset names and their outcome.
391 """
392 trials_table_outcomes = { 1thc
393 'intervals': outcomes.get('_task_iti_delays', spec.QC.NOT_SET),
394 'goCue_times': outcomes.get('_task_goCue_delays', spec.QC.NOT_SET),
395 'response_times': spec.QC.NOT_SET, 'choice': spec.QC.NOT_SET,
396 'stimOn_times': outcomes.get('_task_stimOn_delays', spec.QC.NOT_SET),
397 'contrastLeft': spec.QC.NOT_SET, 'contrastRight': spec.QC.NOT_SET,
398 'feedbackType': spec.QC.NOT_SET, 'probabilityLeft': spec.QC.NOT_SET,
399 'feedback_times': outcomes.get('_task_errorCue_delays', spec.QC.NOT_SET),
400 'firstMovement_times': spec.QC.NOT_SET
401 }
402 reward_checks = ('_task_reward_volumes', '_task_reward_volume_set') 1thc
403 trials_table_outcomes['rewardVolume']: TaskQC.overall_outcome( 1thc
404 (outcomes.get(x, spec.QC.NOT_SET) for x in reward_checks)
405 )
406 dataset_outcomes = { 1thc
407 '_ibl_trials.stimOff_times': outcomes.get('_task_stimOff_delays', spec.QC.NOT_SET),
408 '_ibl_trials.table': trials_table_outcomes,
409 }
410 return dataset_outcomes 1thc
413class HabituationQC(TaskQC):
414 """Task QC for habituation choice world."""
416 def compute(self, **kwargs):
417 """Compute and store the QC metrics.
419 Runs the QC on the session and stores a map of the metrics for each datapoint for each
420 test, and a map of which datapoints passed for each test.
421 """
422 assert self.extractor is not None 1gb
423 self.log.info(f'Session {self.session_path}: Running QC on habituation data...') 1gb
425 # Initialize checks
426 prefix = '_task_' 1gb
427 data = self.extractor.data 1gb
428 audio_output = self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown') 1gb
429 metrics = {} 1gb
430 passed = {} 1gb
432 # Modify criteria based on version
433 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1gb
434 is_v8 = version.parse(ver) >= version.parse('8.0.0') 1gb
435 self.criteria['_task_iti_delays'] = {'PASS': 0.99, 'WARNING': 0} if is_v8 else {'NOT_SET': 0} 1gb
437 # Check all reward volumes == 3.0ul
438 check = prefix + 'reward_volumes' 1gb
439 metrics[check] = data['rewardVolume'] 1gb
440 passed[check] = metrics[check] == 3.0 1gb
442 # Check session durations are increasing in steps >= 12 minutes
443 check = prefix + 'habituation_time' 1gb
444 if not self.one or not self.session_path: 1gb
445 self.log.warning('unable to determine session trials without ONE')
446 metrics[check] = passed[check] = None
447 else:
448 subject, session_date = self.session_path.parts[-3:-1] 1gb
449 # compute from the date specified
450 date_minus_week = ( 1gb
451 datetime.strptime(session_date, '%Y-%m-%d') - timedelta(days=7)
452 ).strftime('%Y-%m-%d')
453 sessions = self.one.alyx.rest('sessions', 'list', subject=subject, 1gb
454 date_range=[date_minus_week, session_date],
455 task_protocol='habituation')
456 # Remove the current session if already registered
457 if sessions and sessions[0]['start_time'].startswith(session_date): 1gb
458 sessions = sessions[1:]
459 metric = ([0, data['intervals'][-1, 1] - data['intervals'][0, 0]] + 1gb
460 [(datetime.fromisoformat(x['end_time']) -
461 datetime.fromisoformat(x['start_time'])).total_seconds() / 60
462 for x in [self.one.alyx.get(s['url']) for s in sessions]])
464 # The duration from raw trial data
465 # duration = map(float, self.extractor.raw_data[-1]['elapsed_time'].split(':'))
466 # duration = timedelta(**dict(zip(('hours', 'minutes', 'seconds'),
467 # duration))).total_seconds() / 60
468 metrics[check] = np.array(metric) 1gb
469 passed[check] = np.diff(metric) >= 12 1gb
471 # Check event orders: trial_start < stim on < stim center < feedback < stim off
472 check = prefix + 'trial_event_sequence' 1gb
473 nans = ( 1gb
474 np.isnan(data['intervals'][:, 0]) | # noqa
475 np.isnan(data['stimOn_times']) | # noqa
476 np.isnan(data['stimCenter_times']) |
477 np.isnan(data['valveOpen_times']) | # noqa
478 np.isnan(data['stimOff_times'])
479 )
480 a = np.less(data['intervals'][:, 0], data['stimOn_times'], where=~nans) 1gb
481 b = np.less(data['stimOn_times'], data['stimCenter_times'], where=~nans) 1gb
482 c = np.less(data['stimCenter_times'], data['valveOpen_times'], where=~nans) 1gb
483 d = np.less(data['valveOpen_times'], data['stimOff_times'], where=~nans) 1gb
485 metrics[check] = a & b & c & d & ~nans 1gb
486 passed[check] = metrics[check].astype(float) 1gb
488 # Check that the time difference between the visual stimulus center-command being
489 # triggered and the stimulus effectively appearing in the center is smaller than 150 ms.
490 check = prefix + 'stimCenter_delays' 1gb
491 metric = np.nan_to_num(data['stimCenter_times'] - data['stimCenterTrigger_times'], 1gb
492 nan=np.inf)
493 passed[check] = (metric <= 0.15) & (metric > 0) 1gb
494 metrics[check] = metric 1gb
496 # Phase check
497 check = prefix + 'phase' 1gb
498 metric = data['phase'] 1gb
499 passed[check] = (metric <= 2 * np.pi) & (metric >= 0) 1gb
500 metrics[check] = metric 1gb
502 # This is not very useful as a check because there are so few trials
503 check = prefix + 'phase_distribution' 1gb
504 metric, _ = np.histogram(data['phase']) 1gb
505 _, p = chisquare(metric) 1gb
506 passed[check] = p < 0.05 if len(data['phase']) >= 400 else None # skip if too few trials 1gb
507 metrics[check] = metric 1gb
509 # Check that the period of gray screen between stim off and the start of the next trial is
510 # 1s +/- 10%.
511 check = prefix + 'iti_delays' 1gb
512 iti = (np.roll(data['stimOn_times'], -1) - data['stimOff_times'])[:-1] 1gb
513 metric = np.r_[np.nan_to_num(iti, nan=np.inf), np.nan] - 1. 1gb
514 passed[check] = np.abs(metric) <= 0.1 1gb
515 passed[check][-1] = np.nan 1gb
516 metrics[check] = metric 1gb
518 # Checks common to training QC
519 checks = [check_goCue_delays, check_stimOn_goCue_delays, 1gb
520 check_stimOn_delays, check_stimOff_delays]
521 for fcn in checks: 1gb
522 check = prefix + fcn.__name__[6:] 1gb
523 metrics[check], passed[check] = fcn(data, audio_output=audio_output) 1gb
525 self.metrics, self.passed = (metrics, passed) 1gb
528# SINGLE METRICS
529# ---------------------------------------------------------------------------- #
531# === Delays between events checks ===
533def check_stimOn_goCue_delays(data, audio_output='harp', **_):
534 """
535 Check the go cue tone occurs less than 10ms before stimulus on.
537 Checks that the time difference between the onset of the visual stimulus
538 and the onset of the go cue tone is positive and less than 10ms.
540 Metric:
541 M = stimOn_times - goCue_times
543 Criteria:
544 0 < M < 0.010 s
546 Units:
547 seconds [s]
549 :param data: dict of trial data with keys ('goCue_times', 'stimOn_times', 'intervals')
550 :param audio_output: audio output device name.
552 Notes
553 -----
554 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
555 percentile of delays over 500 training sessions using the Xonar soundcard.
556 """
557 # Calculate the difference between stimOn and goCue times.
558 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
559 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1agBcbdfe
560 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 1agBcbdfe
561 passed = (metric < threshold) & (metric > 0) 1agBcbdfe
562 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agBcbdfe
563 return metric, passed 1agBcbdfe
566def check_response_feedback_delays(data, audio_output='harp', **_):
567 """
568 Check the feedback delivered within 10ms of the response threshold.
570 Checks that the time difference between the response and the feedback onset
571 (error sound or valve) is positive and less than 10ms.
573 Metric:
574 M = feedback_time - response_time
576 Criterion:
577 0 < M < 0.010 s
579 Units:
580 seconds [s]
582 :param data: dict of trial data with keys ('feedback_times', 'response_times', 'intervals')
583 :param audio_output: audio output device name.
585 Notes
586 -----
587 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
588 percentile of delays over 500 training sessions using the Xonar soundcard.
589 """
590 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1aCcbdfe
591 metric = np.nan_to_num(data['feedback_times'] - data['response_times'], nan=np.inf) 1aCcbdfe
592 passed = (metric < threshold) & (metric > 0) 1aCcbdfe
593 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aCcbdfe
594 return metric, passed 1aCcbdfe
597def check_response_stimFreeze_delays(data, **_):
598 """
599 Check the stimulus freezes within 100ms of the expected time.
601 Checks that the time difference between the visual stimulus freezing and the
602 response is positive and less than 100ms.
604 Metric:
605 M = (stimFreeze_times - response_times)
607 Criterion:
608 0 < M < 0.100 s
610 Units:
611 seconds [s]
613 :param data: dict of trial data with keys ('stimFreeze_times', 'response_times', 'intervals',
614 'choice')
615 """
616 # Calculate the difference between stimOn and goCue times.
617 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
618 metric = np.nan_to_num(data['stimFreeze_times'] - data['response_times'], nan=np.inf) 1aDcbdfe
619 # Test for valid values
620 passed = ((metric < 0.1) & (metric > 0)).astype(float) 1aDcbdfe
621 # Finally remove no_go trials (stimFreeze triggered differently in no_go trials)
622 # These values are ignored in calculation of proportion passed
623 passed[data['choice'] == 0] = np.nan 1aDcbdfe
624 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aDcbdfe
625 return metric, passed 1aDcbdfe
628def check_stimOff_itiIn_delays(data, **_):
629 """Check that the start of the trial interval is within 10ms of the visual stimulus turning off.
631 Metric:
632 M = itiIn_times - stimOff_times
634 Criterion:
635 0 < M < 0.010 s
637 Units:
638 seconds [s]
640 :param data: dict of trial data with keys ('stimOff_times', 'itiIn_times', 'intervals',
641 'choice')
642 """
643 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
644 metric = np.nan_to_num(data['itiIn_times'] - data['stimOff_times'], nan=np.inf) 1aEcbdfe
645 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 1aEcbdfe
646 # Remove no_go trials (stimOff triggered differently in no_go trials)
647 # NaN values are ignored in calculation of proportion passed
648 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan 1aEcbdfe
649 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aEcbdfe
650 return metric, passed 1aEcbdfe
653def check_iti_delays(data, subtract_pauses=False, iti_delay_secs=ITI_DELAY_SECS,
654 feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_):
655 """
656 Check the open-loop grey screen period is approximately 1 second.
658 Check that the period of grey screen between stim off and the start of the next trial is
659 1s +/- 10%. If the trial was paused during this time, the check will account for that
661 Metric:
662 M = stimOff (n) - trialStart (n+1) - 1.
664 Criterion:
665 |M| < 0.1
667 Units:
668 seconds [s]
670 Parameters
671 ----------
672 data : dict
673 Trial data with keys ('stimOff_times', 'intervals', 'pause_duration').
674 subtract_pauses: bool
675 If True, account for experimenter-initiated pauses between trials; if False, trials where
676 the experimenter paused the task may fail this check.
678 Returns
679 -------
680 numpy.array
681 An array of metric values to threshold.
682 numpy.array
683 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
684 """
685 # Initialize array the length of completed trials
686 metric = np.full(data['intervals'].shape[0], np.nan) 1aucbdfe
687 passed = metric.copy() 1aucbdfe
688 pauses = (data['pause_duration'] if subtract_pauses else np.zeros_like(metric))[:-1] 1aucbdfe
689 # Get the difference between stim off and the start of the next trial
690 # Missing data are set to Inf, except for the last trial which is a NaN
691 metric[:-1] = np.nan_to_num( 1aucbdfe
692 data['intervals'][1:, 0] - data['stimOff_times'][:-1] - iti_delay_secs - pauses,
693 nan=np.inf
694 )
695 metric[data['choice'] == 0] = metric[data['choice'] == 0] - feedback_nogo_delay_secs 1aucbdfe
696 passed[:-1] = np.abs(metric[:-1]) < (iti_delay_secs / 10) # Last trial is not counted 1aucbdfe
697 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aucbdfe
698 return metric, passed 1aucbdfe
701def check_positive_feedback_stimOff_delays(data, **_):
702 """
703 Check stimulus offset occurs approximately 1 second after reward delivered.
705 Check that the time difference between the valve onset and the visual stimulus turning off
706 is 1 ± 0.150 seconds.
708 Metric:
709 M = stimOff_times - feedback_times - 1s
711 Criterion:
712 |M| < 0.150 s
714 Units:
715 seconds [s]
717 :param data: dict of trial data with keys ('stimOff_times', 'feedback_times', 'intervals',
718 'correct')
719 """
720 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
721 metric = np.nan_to_num(data['stimOff_times'] - data['feedback_times'] - 1, nan=np.inf) 1aFcbdfe
722 passed = (np.abs(metric) < 0.15).astype(float) 1aFcbdfe
723 # NaN values are ignored in calculation of proportion passed; ignore incorrect trials here
724 metric[~data['correct']] = passed[~data['correct']] = np.nan 1aFcbdfe
725 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aFcbdfe
726 return metric, passed 1aFcbdfe
729def check_negative_feedback_stimOff_delays(data, feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_):
730 """
731 Check the stimulus offset occurs approximately 2 seconds after negative feedback delivery.
733 Check that the time difference between the error sound and the visual stimulus
734 turning off is 2 ± 0.150 seconds.
736 Metric:
737 M = stimOff_times - errorCue_times - 2s
739 Criterion:
740 |M| < 0.150 s
742 Units:
743 seconds [s]
745 :param data: dict of trial data with keys ('stimOff_times', 'errorCue_times', 'intervals')
746 """
747 metric = np.nan_to_num(data['stimOff_times'] - data['errorCue_times'] - 2, nan=np.inf) 1azcbdfe
748 # for the nogo trials, the feedback is the same as the stimOff
749 metric[data['choice'] == 0] = metric[data['choice'] == 0] + feedback_nogo_delay_secs 1azcbdfe
750 # Apply criteria
751 passed = (np.abs(metric) < 0.15).astype(float) 1azcbdfe
752 # Remove none negative feedback trials
753 metric[data['correct']] = passed[data['correct']] = np.nan 1azcbdfe
754 assert data['intervals'].shape[0] == len(metric) == len(passed) 1azcbdfe
755 return metric, passed 1azcbdfe
758# === Wheel movement during trial checks ===
760def check_wheel_move_before_feedback(data, **_):
761 """Check that the wheel does move within 100ms of the feedback onset (error sound or valve).
763 Metric:
764 M = (w_t - 0.05) - (w_t + 0.05), where t = feedback_times
766 Criterion:
767 M != 0
769 Units:
770 radians
772 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice',
773 'intervals', 'feedback_times')
774 """
775 # Get tuple of wheel times and positions within 100ms of feedback
776 traces = traces_by_trial( 1aocbdfe
777 data['wheel_timestamps'],
778 data['wheel_position'],
779 start=data['feedback_times'] - 0.05,
780 end=data['feedback_times'] + 0.05,
781 )
782 metric = np.zeros_like(data['feedback_times']) 1aocbdfe
783 # For each trial find the displacement
784 for i, trial in enumerate(traces): 1aocbdfe
785 pos = trial[1] 1aocbdfe
786 if pos.size > 1: 1aocbdfe
787 metric[i] = pos[-1] - pos[0] 1aocbdfe
789 # except no-go trials
790 metric[data['choice'] == 0] = np.nan # NaN = trial ignored for this check 1aocbdfe
791 nans = np.isnan(metric) 1aocbdfe
792 passed = np.zeros_like(metric) * np.nan 1aocbdfe
794 passed[~nans] = (metric[~nans] != 0).astype(float) 1aocbdfe
795 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aocbdfe
796 return metric, passed 1aocbdfe
799def _wheel_move_during_closed_loop(re_ts, re_pos, data, wheel_gain=None, tol=1, **_):
800 """
801 Check the wheel moves the correct amount to reach threshold.
803 Check that the wheel moves by approximately 35 degrees during the closed-loop period
804 on trials where a feedback (error sound or valve) is delivered.
806 Metric:
807 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
808 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
809 move 35 visual degrees
811 Criterion:
812 displacement < tol visual degree
814 Units:
815 degrees angle of wheel turn
817 :param re_ts: extracted wheel timestamps in seconds
818 :param re_pos: extracted wheel positions in radians
819 :param data: a dict with the keys (goCueTrigger_times, response_times, feedback_times,
820 position, choice, intervals)
821 :param wheel_gain: the 'STIM_GAIN' task setting
822 :param tol: the criterion in visual degrees
823 """
824 if wheel_gain is None: 1ajcbdfe
825 _log.warning('No wheel_gain input in function call, returning None')
826 return None, None
828 # Get tuple of wheel times and positions over each trial's closed-loop period
829 traces = traces_by_trial(re_ts, re_pos, 1ajcbdfe
830 start=data['goCueTrigger_times'],
831 end=data['response_times'])
833 metric = np.zeros_like(data['feedback_times']) 1ajcbdfe
834 # For each trial find the absolute displacement
835 for i, trial in enumerate(traces): 1ajcbdfe
836 t, pos = trial 1ajcbdfe
837 if pos.size != 0: 1ajcbdfe
838 # Find the position of the preceding sample and subtract it
839 idx = np.abs(re_ts - t[0]).argmin() - 1 1ajcbdfe
840 origin = re_pos[idx] 1ajcbdfe
841 metric[i] = np.abs(pos - origin).max() 1ajcbdfe
843 # Load wheel_gain and thresholds for each trial
844 wheel_gain = np.array([wheel_gain] * len(data['position'])) 1ajcbdfe
845 thresh = data['position'] 1ajcbdfe
846 # abs displacement, s, in mm required to move 35 visual degrees
847 s_mm = np.abs(thresh / wheel_gain) # don't care about direction 1ajcbdfe
848 criterion = cm_to_rad(s_mm * 1e-1) # convert abs displacement to radians (wheel pos is in rad) 1ajcbdfe
849 metric = metric - criterion # difference should be close to 0 1ajcbdfe
850 rad_per_deg = cm_to_rad(1 / wheel_gain * 1e-1) 1ajcbdfe
851 passed = (np.abs(metric) < rad_per_deg * tol).astype(float) # less than 1 visual degree off 1ajcbdfe
852 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan # except no-go trials 1ajcbdfe
853 assert data['intervals'].shape[0] == len(metric) == len(passed) 1ajcbdfe
854 return metric, passed 1ajcbdfe
857def check_wheel_move_during_closed_loop(data, wheel_gain=None, **_):
858 """
859 Check the wheel moves the correct amount to reach threshold.
861 Check that the wheel moves by approximately 35 degrees during the closed-loop period
862 on trials where a feedback (error sound or valve) is delivered.
864 Metric:
865 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
866 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
867 move 35 visual degrees
869 Criterion:
870 displacement < 3 visual degrees
872 Units:
873 degrees angle of wheel turn
875 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice',
876 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position')
877 :param wheel_gain: the 'STIM_GAIN' task setting
878 """
879 # Get the Bpod extracted wheel data
880 timestamps = data['wheel_timestamps'] 1ajcbdfe
881 position = data['wheel_position'] 1ajcbdfe
883 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=3) 1ajcbdfe
886def check_wheel_move_during_closed_loop_bpod(data, wheel_gain=None, **_):
887 """
888 Check the wheel moves the correct amount to reach threshold.
890 Check that the wheel moves by approximately 35 degrees during the closed-loop period
891 on trials where a feedback (error sound or valve) is delivered. This check uses the Bpod
892 wheel data (measured at a lower resolution) with a stricter tolerance (1 visual degree).
894 Metric:
895 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
896 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
897 move 35 visual degrees.
899 Criterion:
900 displacement < 1 visual degree
902 Units:
903 degrees angle of wheel turn
905 :param data: dict of trial data with keys ('wheel_timestamps(_bpod)', 'wheel_position(_bpod)',
906 'choice', 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position')
907 :param wheel_gain: the 'STIM_GAIN' task setting
908 """
909 # Get the Bpod extracted wheel data
910 timestamps = data.get('wheel_timestamps_bpod', data['wheel_timestamps']) 1acbdfe
911 position = data.get('wheel_position_bpod', data['wheel_position']) 1acbdfe
913 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=1) 1acbdfe
916def check_wheel_freeze_during_quiescence(data, **_):
917 """
918 Check the wheel is indeed still during the quiescent period.
920 Check that the wheel does not move more than 2 degrees in each direction during the
921 quiescence interval before the stimulus appears.
923 Metric:
924 M = |max(W) - min(W)| where W is wheel pos over quiescence interval;
925 interval = [stimOnTrigger_times - quiescent_duration, stimOnTrigger_times]
927 Criterion:
928 M < 2 degrees
930 Units:
931 degrees angle of wheel turn
933 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'quiescence',
934 'intervals', 'stimOnTrigger_times')
935 """
936 assert np.all(np.diff(data['wheel_timestamps']) >= 0) 1akcbdfe
937 assert data['quiescence'].size == data['stimOnTrigger_times'].size 1akcbdfe
938 # Get tuple of wheel times and positions over each trial's quiescence period
939 qevt_start_times = data['stimOnTrigger_times'] - data['quiescence'] 1akcbdfe
940 traces = traces_by_trial( 1akcbdfe
941 data['wheel_timestamps'],
942 data['wheel_position'],
943 start=qevt_start_times,
944 end=data['stimOnTrigger_times']
945 )
947 metric = np.zeros((len(data['quiescence']), 2)) # (n_trials, n_directions) 1akcbdfe
948 for i, trial in enumerate(traces): 1akcbdfe
949 t, pos = trial 1akcbdfe
950 # Get the last position before the period began
951 if pos.size > 0: 1akcbdfe
952 # Find the position of the preceding sample and subtract it
953 idx = np.abs(data['wheel_timestamps'] - t[0]).argmin() - 1 1akcbdfe
954 origin = data['wheel_position'][idx if idx != -1 else 0] 1akcbdfe
955 # Find the absolute min and max relative to the last sample
956 metric[i, :] = np.abs([np.min(pos - origin), np.max(pos - origin)]) 1akcbdfe
957 # Reduce to the largest displacement found in any direction
958 metric = np.max(metric, axis=1) 1akcbdfe
959 metric = 180 * metric / np.pi # convert to degrees from radians 1akcbdfe
960 criterion = 2 # Position shouldn't change more than 2 in either direction 1akcbdfe
961 passed = metric < criterion 1akcbdfe
962 assert data['intervals'].shape[0] == len(metric) == len(passed) 1akcbdfe
963 return metric, passed 1akcbdfe
966def check_detected_wheel_moves(data, min_qt=0, **_):
967 """Check that the detected first movement times are reasonable.
969 Metric:
970 M = firstMovement times
972 Criterion:
973 (goCue trigger time - min quiescent period) < M < response time
975 Units:
976 Seconds [s]
978 :param data: dict of trial data with keys ('firstMovement_times', 'goCueTrigger_times',
979 'response_times', 'choice', 'intervals')
980 :param min_qt: the minimum possible quiescent period (the QUIESCENT_PERIOD task parameter)
981 """
982 # Depending on task version this may be a single value or an array of quiescent periods
983 min_qt = np.array(min_qt) 1ascbdfe
984 if min_qt.size > data['intervals'].shape[0]: 1ascbdfe
985 min_qt = min_qt[:data['intervals'].shape[0]] 1d
987 metric = data['firstMovement_times'] 1ascbdfe
988 qevt_start = data['goCueTrigger_times'] - np.array(min_qt) 1ascbdfe
989 response = data['response_times'] 1ascbdfe
990 # First movement time for each trial should be after the quiescent period and before feedback
991 passed = np.array([a < m < b for m, a, b in zip(metric, qevt_start, response)], dtype=float) 1ascbdfe
992 nogo = data['choice'] == 0 1ascbdfe
993 passed[nogo] = np.nan # No go trial may have no movement times and that's fine 1ascbdfe
994 return metric, passed 1ascbdfe
997# === Sequence of events checks ===
999def check_error_trial_event_sequence(data, **_):
1000 """
1001 Check trial events occur in correct order for negative feedback trials.
1003 Check that on incorrect / miss trials, there are exactly:
1004 2 audio events (go cue sound and error sound) and 2 Bpod events (trial start, ITI), occurring
1005 in the correct order
1007 Metric:
1008 M = Bpod (trial start) > audio (go cue) > audio (error) > Bpod (ITI) > Bpod (trial end)
1010 Criterion:
1011 M == True
1013 Units:
1014 -none-
1016 :param data: dict of trial data with keys ('errorCue_times', 'goCue_times', 'intervals',
1017 'itiIn_times', 'correct')
1018 """
1019 # An array the length of N trials where True means at least one event time was NaN (bad)
1020 nans = ( 1aqcbdfe
1021 np.isnan(data['intervals'][:, 0]) |
1022 np.isnan(data['goCue_times']) | # noqa
1023 np.isnan(data['errorCue_times']) | # noqa
1024 np.isnan(data['itiIn_times']) | # noqa
1025 np.isnan(data['intervals'][:, 1])
1026 )
1028 # For each trial check that the events happened in the correct order (ignore NaN values)
1029 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1aqcbdfe
1030 b = np.less(data['goCue_times'], data['errorCue_times'], where=~nans) # Go cue < error cue 1aqcbdfe
1031 c = np.less(data['errorCue_times'], data['itiIn_times'], where=~nans) # Error cue < ITI start 1aqcbdfe
1032 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1aqcbdfe
1034 # For each trial check all events were in order AND all event times were not NaN
1035 metric = a & b & c & d & ~nans 1aqcbdfe
1037 passed = metric.astype(float) 1aqcbdfe
1038 passed[data['correct']] = np.nan # Look only at incorrect trials 1aqcbdfe
1039 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aqcbdfe
1040 return metric, passed 1aqcbdfe
1043def check_correct_trial_event_sequence(data, **_):
1044 """
1045 Check trial events occur in correct order for positive feedback trials.
1047 Check that on correct trials, there are exactly:
1048 1 audio events and 3 Bpod events (valve open, trial start, ITI), occurring in the correct order
1050 Metric:
1051 M = Bpod (trial start) > audio (go cue) > Bpod (valve) > Bpod (ITI) > Bpod (trial end)
1053 Criterion:
1054 M == True
1056 Units:
1057 -none-
1059 :param data: dict of trial data with keys ('valveOpen_times', 'goCue_times', 'intervals',
1060 'itiIn_times', 'correct')
1061 """
1062 # An array the length of N trials where True means at least one event time was NaN (bad)
1063 nans = ( 1arcbdfe
1064 np.isnan(data['intervals'][:, 0]) |
1065 np.isnan(data['goCue_times']) | # noqa
1066 np.isnan(data['valveOpen_times']) |
1067 np.isnan(data['itiIn_times']) | # noqa
1068 np.isnan(data['intervals'][:, 1])
1069 )
1071 # For each trial check that the events happened in the correct order (ignore NaN values)
1072 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1arcbdfe
1073 b = np.less(data['goCue_times'], data['valveOpen_times'], where=~nans) # Go cue < feedback 1arcbdfe
1074 c = np.less(data['valveOpen_times'], data['itiIn_times'], where=~nans) # Feedback < ITI start 1arcbdfe
1075 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1arcbdfe
1077 # For each trial True means all events were in order AND all event times were not NaN
1078 metric = a & b & c & d & ~nans 1arcbdfe
1080 passed = metric.astype(float) 1arcbdfe
1081 passed[~data['correct']] = np.nan # Look only at correct trials 1arcbdfe
1082 assert data['intervals'].shape[0] == len(metric) == len(passed) 1arcbdfe
1083 return metric, passed 1arcbdfe
1086def check_n_trial_events(data, **_):
1087 """Check that the number events per trial is correct.
1089 Within every trial interval there should be one of each trial event, except for
1090 goCueTrigger_times which should only be defined for incorrect trials
1092 Metric:
1093 M = all(start < event < end) for all event times except errorCueTrigger_times where
1094 start < error_trigger < end if not correct trial, else error_trigger == NaN
1096 Criterion:
1097 M == True
1099 Units:
1100 -none-, boolean
1102 :param data: dict of trial data with keys ('intervals', 'stimOnTrigger_times',
1103 'stimOffTrigger_times', 'stimOn_times', 'stimOff_times',
1104 'stimFreezeTrigger_times', 'errorCueTrigger_times', 'itiIn_times',
1105 'goCueTrigger_times', 'goCue_times', 'response_times', 'feedback_times')
1106 """
1107 intervals = data['intervals'] 1apcbdfe
1108 correct = data['correct'] 1apcbdfe
1109 err_trig = data['errorCueTrigger_times'] 1apcbdfe
1111 # Exclude these fields; valve and errorCue times are the same as feedback_times and we must
1112 # test errorCueTrigger_times separately
1113 # stimFreeze_times fails often due to TTL flicker
1114 exclude = ['camera_timestamps', 'errorCueTrigger_times', 'errorCue_times', 1apcbdfe
1115 'wheelMoves_peakVelocity_times', 'valveOpen_times', 'wheelMoves_peakAmplitude',
1116 'wheelMoves_intervals', 'wheel_timestamps', 'stimFreeze_times']
1117 events = [k for k in data.keys() if k.endswith('_times') and k not in exclude] 1apcbdfe
1118 metric = np.zeros(data['intervals'].shape[0], dtype=bool) 1apcbdfe
1120 # For each trial interval check that one of each trial event occurred. For incorrect trials,
1121 # check the error cue trigger occurred within the interval, otherwise check it is nan.
1122 for i, (start, end) in enumerate(intervals): 1apcbdfe
1123 metric[i] = (all([start < data[k][i] < end for k in events]) and 1apcbdfe
1124 (np.isnan(err_trig[i]) if correct[i] else start < err_trig[i] < end))
1125 passed = metric.astype(bool) 1apcbdfe
1126 assert intervals.shape[0] == len(metric) == len(passed) 1apcbdfe
1127 return metric, passed 1apcbdfe
1130def check_trial_length(data, **_):
1131 """
1132 Check open-loop duration positive and <= 1 minute.
1134 Check that the time difference between the onset of the go cue sound
1135 and the feedback (error sound or valve) is positive and smaller than 60.1 s.
1137 Metric:
1138 M = feedback_times - goCue_times
1140 Criteria:
1141 0 < M < 60.1 s
1143 Units:
1144 seconds [s]
1146 :param data: dict of trial data with keys ('feedback_times', 'goCue_times', 'intervals')
1147 """
1148 # NaN values are usually ignored so replace them with Inf so they fail the threshold
1149 metric = np.nan_to_num(data['feedback_times'] - data['goCue_times'], nan=np.inf) 1aHcbdfe
1150 passed = (metric < 60.1) & (metric > 0) 1aHcbdfe
1151 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aHcbdfe
1152 return metric, passed 1aHcbdfe
1155# === Trigger-response delay checks ===
1157def check_goCue_delays(data, audio_output='harp', **_):
1158 """
1159 Check the go cue tone occurs within 1ms of the intended time.
1161 Check that the time difference between the go cue sound being triggered and
1162 effectively played is smaller than 1ms.
1164 Metric:
1165 M = goCue_times - goCueTrigger_times
1167 Criterion:
1168 0 < M <= 0.0015 s
1170 Units:
1171 seconds [s]
1173 :param data: dict of trial data with keys ('goCue_times', 'goCueTrigger_times', 'intervals').
1174 :param audio_output: audio output device name.
1176 Notes
1177 -----
1178 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
1179 percentile of delays over 500 training sessions using the Xonar soundcard.
1180 """
1181 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.053 1agGcbdfe
1182 metric = np.nan_to_num(data['goCue_times'] - data['goCueTrigger_times'], nan=np.inf) 1agGcbdfe
1183 passed = (metric <= threshold) & (metric > 0) 1agGcbdfe
1184 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agGcbdfe
1185 return metric, passed 1agGcbdfe
1188def check_errorCue_delays(data, audio_output='harp', **_):
1189 """
1190 Check the error tone occurs within 1.5ms of the intended time.
1192 Check that the time difference between the error sound being triggered and
1193 effectively played is smaller than 1ms.
1195 Metric:
1196 M = errorCue_times - errorCueTrigger_times
1198 Criterion:
1199 0 < M <= 0.0015 s
1201 Units:
1202 seconds [s]
1204 :param data: dict of trial data with keys ('errorCue_times', 'errorCueTrigger_times',
1205 'intervals', 'correct')
1206 :param audio_output: audio output device name.
1208 Notes
1209 -----
1210 - For non-harp sound card the permissible delay is 0.062s. This was chosen by taking the 99.5th
1211 percentile of delays over 500 training sessions using the Xonar soundcard.
1212 """
1213 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.062 1aAcbdfe
1214 metric = np.nan_to_num(data['errorCue_times'] - data['errorCueTrigger_times'], nan=np.inf) 1aAcbdfe
1215 passed = ((metric <= threshold) & (metric > 0)).astype(float) 1aAcbdfe
1216 passed[data['correct']] = metric[data['correct']] = np.nan 1aAcbdfe
1217 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aAcbdfe
1218 return metric, passed 1aAcbdfe
1221def check_stimOn_delays(data, **_):
1222 """
1223 Check the visual stimulus onset occurs within 150ms of the intended time.
1225 Check that the time difference between the visual stimulus onset-command being triggered
1226 and the stimulus effectively appearing on the screen is smaller than 150 ms.
1228 Metric:
1229 M = stimOn_times - stimOnTrigger_times
1231 Criterion:
1232 0 < M < 0.15 s
1234 Units:
1235 seconds [s]
1237 :param data: dict of trial data with keys ('stimOn_times', 'stimOnTrigger_times',
1238 'intervals')
1239 """
1240 metric = np.nan_to_num(data['stimOn_times'] - data['stimOnTrigger_times'], nan=np.inf) 1agIcbdfe
1241 passed = (metric <= 0.15) & (metric > 0) 1agIcbdfe
1242 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agIcbdfe
1243 return metric, passed 1agIcbdfe
1246def check_stimOff_delays(data, **_):
1247 """
1248 Check stimulus offset occurs within 150ms of the intended time.
1250 Check that the time difference between the visual stimulus offset-command
1251 being triggered and the visual stimulus effectively turning off on the screen
1252 is smaller than 150 ms.
1254 Metric:
1255 M = stimOff_times - stimOffTrigger_times
1257 Criterion:
1258 0 < M < 0.15 s
1260 Units:
1261 seconds [s]
1263 :param data: dict of trial data with keys ('stimOff_times', 'stimOffTrigger_times',
1264 'intervals')
1265 """
1266 metric = np.nan_to_num(data['stimOff_times'] - data['stimOffTrigger_times'], nan=np.inf) 1agJcbdfe
1267 passed = (metric <= 0.15) & (metric > 0) 1agJcbdfe
1268 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agJcbdfe
1269 return metric, passed 1agJcbdfe
1272def check_stimFreeze_delays(data, **_):
1273 """Check the stimulus freezes within 150ms of the intended time.
1275 Check that the time difference between the visual stimulus freeze-command
1276 being triggered and the visual stimulus effectively freezing on the screen
1277 is smaller than 150 ms.
1279 Metric:
1280 M = stimFreeze_times - stimFreezeTrigger_times
1282 Criterion:
1283 0 < M < 0.15 s
1285 Units:
1286 seconds [s]
1288 :param data: dict of trial data with keys ('stimFreeze_times', 'stimFreezeTrigger_times',
1289 'intervals')
1290 """
1291 metric = np.nan_to_num(data['stimFreeze_times'] - data['stimFreezeTrigger_times'], nan=np.inf) 1aKcbdfe
1292 passed = (metric <= 0.15) & (metric > 0) 1aKcbdfe
1293 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aKcbdfe
1294 return metric, passed 1aKcbdfe
1297# === Data integrity checks ===
1299def check_reward_volumes(data, **_):
1300 """Check that the reward volume is between 1.5 and 3 uL for correct trials, 0 for incorrect.
1302 Metric:
1303 M = reward volume
1305 Criterion:
1306 1.5 <= M <= 3 if correct else M == 0
1308 Units:
1309 uL
1311 :param data: dict of trial data with keys ('rewardVolume', 'correct', 'intervals')
1312 """
1313 metric = data['rewardVolume'] 1aycbdfe
1314 correct = data['correct'] 1aycbdfe
1315 passed = np.zeros_like(metric, dtype=bool) 1aycbdfe
1316 # Check correct trials within correct range
1317 passed[correct] = (1.5 <= metric[correct]) & (metric[correct] <= 3.) 1aycbdfe
1318 # Check incorrect trials are 0
1319 passed[~correct] = metric[~correct] == 0 1aycbdfe
1320 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aycbdfe
1321 return metric, passed 1aycbdfe
1324def check_reward_volume_set(data, **_):
1325 """Check that there is only two reward volumes within a session, one of which is 0.
1327 Metric:
1328 M = set(rewardVolume)
1330 Criterion:
1331 (0 < len(M) <= 2) and 0 in M
1333 :param data: dict of trial data with keys ('rewardVolume')
1334 """
1335 metric = data['rewardVolume'] 1aLcbdfe
1336 passed = 0 < len(set(metric)) <= 2 and 0. in metric 1aLcbdfe
1337 return metric, passed 1aLcbdfe
1340def check_wheel_integrity(data, re_encoding='X1', enc_res=None, **_):
1341 """
1342 Check wheel position sampled at the expected resolution.
1344 Check that the difference between wheel position samples is close to the encoder resolution
1345 and that the wheel timestamps strictly increase.
1347 Metric:
1348 M = (absolute difference of the positions < 1.5 * encoder resolution)
1349 + 1 if (difference of timestamps <= 0) else 0
1351 Criterion:
1352 M ~= 0
1354 Units:
1355 arbitrary (radians, sometimes + 1)
1357 :param data: dict of wheel data with keys ('wheel_timestamps', 'wheel_position')
1358 :param re_encoding: the encoding of the wheel data, X1, X2 or X4
1359 :param enc_res: the rotary encoder resolution (default 1024 ticks per revolution)
1361 Notes
1362 -----
1363 - At high velocities some samples are missed due to the scanning frequency of the DAQ.
1364 This checks for more than 1 missing sample in a row (i.e. the difference between samples >= 2)
1366 """
1367 if isinstance(re_encoding, str): 1avcbdfe
1368 re_encoding = int(re_encoding[-1]) 1avcbdfe
1369 # The expected difference between samples in the extracted units
1370 resolution = 1 / (enc_res or ephys_fpga.WHEEL_TICKS 1avcbdfe
1371 ) * np.pi * 2 * ephys_fpga.WHEEL_RADIUS_CM / re_encoding
1372 # We expect the difference of neighbouring positions to be close to the resolution
1373 pos_check = np.abs(np.diff(data['wheel_position'])) 1avcbdfe
1374 # Timestamps should be strictly increasing
1375 ts_check = np.diff(data['wheel_timestamps']) <= 0. 1avcbdfe
1376 metric = pos_check + ts_check.astype(float) # all values should be close to zero 1avcbdfe
1377 passed = metric < 1.5 * resolution 1avcbdfe
1378 return metric, passed 1avcbdfe
1381# === Pre-stimulus checks ===
1382def check_stimulus_move_before_goCue(data, photodiode=None, **_):
1383 """
1384 Check there are no stimulus events before the go cue tone.
1386 Check that there are no visual stimulus change(s) between the start of the trial and the
1387 go cue sound onset, except for stim on.
1389 Metric:
1390 M = number of visual stimulus change events between trial start and goCue_times
1392 Criterion:
1393 M == 1
1395 Units:
1396 -none-, integer
1398 Parameters
1399 ----------
1400 data : dict
1401 Trial data with keys ('goCue_times', 'intervals', 'choice').
1402 photodiode : dict
1403 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel.
1405 Returns
1406 -------
1407 numpy.array
1408 An array of metric values to threshold.
1409 numpy.array
1410 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
1412 Notes
1413 -----
1414 - There should be exactly 1 stimulus change before goCue; stimulus onset. Even if the stimulus
1415 contrast is 0, the sync square will still flip at stimulus onset, etc.
1416 - If there are no goCue times (all are NaN), the status should be NOT_SET.
1417 """
1418 if photodiode is None: 1acbdfe
1419 _log.warning('No photodiode TTL input in function call, returning None')
1420 return None
1421 photodiode_clean = ephys_fpga._clean_frame2ttl(photodiode) 1acbdfe
1422 s = photodiode_clean['times'] 1acbdfe
1423 s = s[~np.isnan(s)] # Remove NaNs 1acbdfe
1424 metric = np.array([]) 1acbdfe
1425 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1acbdfe
1426 metric = np.append(metric, np.count_nonzero(s[s > i] < c)) 1acbdfe
1428 passed = (metric == 1).astype(float) 1acbdfe
1429 passed[np.isnan(data['goCue_times'])] = np.nan 1acbdfe
1430 assert data['intervals'].shape[0] == len(metric) == len(passed) 1acbdfe
1431 return metric, passed 1acbdfe
1434def check_audio_pre_trial(data, audio=None, **_):
1435 """
1436 Check no audio stimuli before the go cue.
1438 Check there are no audio outputs between the start of the trial and the go cue sound onset - 20 ms.
1440 Metric:
1441 M = sum(start_times < audio TTL < (goCue_times - 20ms))
1443 Criterion:
1444 M == 0
1446 Units:
1447 -none-, integer
1449 Parameters
1450 ----------
1451 data : dict
1452 Trial data with keys ('goCue_times', 'intervals').
1453 audio : dict
1454 The fronts from Bpod's BNC2 input FPGA audio sync channel.
1456 Returns
1457 -------
1458 numpy.array
1459 An array of metric values to threshold.
1460 numpy.array
1461 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
1462 """
1463 if audio is None: 1awcbdfe
1464 _log.warning('No BNC2 input in function call, retuning None')
1465 return None, None
1466 s = audio['times'][~np.isnan(audio['times'])] # Audio TTLs with NaNs removed 1awcbdfe
1467 metric = np.array([], dtype=np.int8) 1awcbdfe
1468 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1awcbdfe
1469 metric = np.append(metric, sum(s[s > i] < (c - 0.02))) 1awcbdfe
1470 passed = metric == 0 1awcbdfe
1471 assert data['intervals'].shape[0] == len(metric) == len(passed) 1awcbdfe
1472 return metric, passed 1awcbdfe