Coverage for ibllib/qc/task_metrics.py: 97%
410 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-07-02 18:55 +0100
« prev ^ index » next coverage.py v7.9.1, created at 2025-07-02 18:55 +0100
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 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 'stimOff_itiIn_delays': {'PASS': 0.99, 'WARNING': 0},
75 'positive_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
76 'negative_feedback_stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
77 'wheel_move_during_closed_loop': {'PASS': 0.99, 'WARNING': 0},
78 'response_stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0},
79 'detected_wheel_moves': {'PASS': 0.99, 'WARNING': 0},
80 'trial_length': {'PASS': 0.99, 'WARNING': 0},
81 'goCue_delays': {'PASS': 0.99, 'WARNING': 0},
82 'errorCue_delays': {'PASS': 0.99, 'WARNING': 0},
83 'stimOn_delays': {'PASS': 0.99, 'WARNING': 0},
84 'stimOff_delays': {'PASS': 0.99, 'WARNING': 0},
85 'stimFreeze_delays': {'PASS': 0.99, 'WARNING': 0},
86 'iti_delays': {'NOT_SET': 0},
87 'passed_trial_checks': {'NOT_SET': 0}
88}
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, qc.namespace).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, namespace='task', **kwargs):
210 """
211 Task QC for training, biased, and ephys choice world.
213 Parameters
214 ----------
215 session_path_or_eid : str, uuid.UUID, pathlib.Path
216 A session eid or path. If a path, it must be a valid ALFPath.
217 namespace : str
218 The namespace of the QC fields in the Alyx JSON field.
219 log : logging.Logger, optional
220 A logging.Logger instance, if None the 'ibllib.qc.base' logger is used.
221 one : one.api.OneAlyx, optional
222 An ONE instance for fetching and setting the QC on Alyx.
223 """
224 super().__init__(session_path_or_eid, **kwargs) 1ahcbd
226 # Data
227 self.extractor = None 1ahcbd
229 # Metrics and passed trials
230 self.metrics = None 1ahcbd
231 self.passed = None 1ahcbd
233 # Criteria (initialize as outcomes vary by class, task, and hardware)
234 self.namespace = namespace 1ahcbd
235 self.criteria = {k if k == 'default' else f'_{namespace}_{k}': v 1ahcbd
236 for k, v in BWM_CRITERIA.items()}
238 def compute(self, **kwargs):
239 """Compute and store the QC metrics.
241 Runs the QC on the session and stores a map of the metrics for each datapoint for each
242 test, and a map of which datapoints passed for each test.
244 Parameters
245 ----------
246 bpod_only : bool
247 If True no data is extracted from the FPGA for ephys sessions.
248 """
249 assert self.extractor is not None 1acbdfe
250 ns = kwargs.pop('namespace', self.namespace) 1acbdfe
251 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1acbdfe
252 if version.parse(ver) >= version.parse('8.0.0'): 1acbdfe
253 self.criteria[f'_{ns}_iti_delays'] = {'PASS': 0.99, 'WARNING': 0}
254 self.criteria[f'_{ns}_passed_trial_checks'] = {'PASS': 0.7, 'WARNING': 0}
255 else:
256 self.criteria[f'_{ns}_iti_delays'] = {'NOT_SET': 0} 1acbdfe
257 self.criteria[f'_{ns}_passed_trial_checks'] = {'NOT_SET': 0} 1acbdfe
259 self.log.info(f'Session {self.session_path}: Running QC on behavior data...') 1acbdfe
260 self.get_bpodqc_metrics_frame( 1acbdfe
261 self.extractor.data,
262 wheel_gain=self.extractor.settings['STIM_GAIN'], # The wheel gain
263 photodiode=self.extractor.frame_ttls,
264 audio=self.extractor.audio_ttls,
265 re_encoding=self.extractor.wheel_encoding or 'X1',
266 min_qt=self.extractor.settings.get('QUIESCENT_PERIOD') or 0.2,
267 audio_output=self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown')
268 )
270 def _get_checks(self):
271 """
272 Find all methods that begin with 'check_'.
274 Returns
275 -------
276 Dict[str, function]
277 A map of QC check function names and the corresponding functions that return `metric`
278 (any), `passed` (bool).
279 """
280 def is_metric(x): 1athcbdfe
281 return isfunction(x) and x.__name__.startswith('check_') 1athcbdfe
283 return dict(getmembers(sys.modules[__name__], is_metric)) 1athcbdfe
285 def get_bpodqc_metrics_frame(self, data, **kwargs):
286 """
287 Evaluate task QC metrics.
289 Evaluates all the QC metric functions in this module (those starting with 'check') and
290 returns the results. The optional kwargs listed below are passed to each QC metric function.
292 Parameters
293 ----------
294 data : dict
295 The extracted task data.
296 re_encoding : str {'X1', 'X2', 'X4'}
297 The encoding configuration of the rotary encoder.
298 enc_res : int
299 The rotary encoder resolution as number of fronts per revolution.
300 wheel_gain : float
301 The STIM_GAIN task parameter.
302 photodiode : dict
303 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel.
304 audio : dict
305 The fronts from Bpod's BNC2 input FPGA audio sync channel.
306 min_qt : float
307 The QUIESCENT_PERIOD task parameter.
309 Returns
310 -------
311 dict
312 Map of checks and their QC metric values (1 per trial).
313 dict
314 Map of checks and a float array of which samples passed.
315 """
316 # Find all methods that begin with 'check_'
317 checks = self._get_checks() 1acbdfe
318 prefix = f'_{kwargs.pop("namespace", self.namespace)}_' # Extended QC fields will start with this 1acbdfe
319 # Method 'check_foobar' stored with key '_task_foobar' in metrics map
320 qc_metrics_map = {prefix + k[6:]: fn(data, **kwargs) for k, fn in checks.items()} 1acbdfe
322 # Split metrics and passed frames
323 self.metrics = {} 1acbdfe
324 self.passed = {} 1acbdfe
325 for k in qc_metrics_map: 1acbdfe
326 self.metrics[k], self.passed[k] = qc_metrics_map[k] 1acbdfe
328 # Add a check for trial level pass: did a given trial pass all checks?
329 n_trials = data['intervals'].shape[0] 1acbdfe
330 # Trial-level checks return an array the length that equals the number of trials
331 trial_level_passed = [m for m in self.passed.values() if isinstance(m, Sized) and len(m) == n_trials] 1acbdfe
332 name = prefix + 'passed_trial_checks' 1acbdfe
333 self.metrics[name] = reduce(np.logical_and, trial_level_passed or (None, None)) 1acbdfe
334 self.passed[name] = self.metrics[name].astype(float) if trial_level_passed else None 1acbdfe
336 def run(self, update=False, **kwargs):
337 """
338 Compute the QC outcomes and return overall task QC outcome.
340 Parameters
341 ----------
342 update : bool
343 If True, updates the session QC fields on Alyx.
344 bpod_only : bool
345 If True no data is extracted from the FPGA for ephys sessions.
347 Returns
348 -------
349 str
350 Overall task QC outcome.
351 dict
352 A map of QC tests and the proportion of data points that passed them.
353 """
354 if self.metrics is None: 1acbdi
355 self.compute(**kwargs) 1acbd
356 outcome, results, _ = self.compute_session_status() 1acbdi
357 if update: 1acbdi
358 self.update_extended_qc(results) 1ci
359 self.update(outcome, kwargs.get('namespace', self.namespace)) 1ci
360 return outcome, results 1acbdi
362 def compute_session_status(self):
363 """
364 Compute the overall session QC for each key and aggregates in a single value.
366 Returns
367 -------
368 str
369 Overall session QC outcome.
370 dict
371 A map of QC tests and the proportion of data points that passed them.
372 dict
373 A map of QC tests and their outcomes.
374 """
375 if self.passed is None: 1agcbdei
376 raise AttributeError('passed is None; compute QC first') 1e
377 # Get mean passed of each check, or None if passed is None or all NaN
378 results = {k: None if v is None or np.isnan(v).all() else np.nanmean(v) 1agcbdei
379 for k, v in self.passed.items()}
380 session_outcome, outcomes = compute_session_status_from_dict(results, self.criteria) 1agcbdei
381 return session_outcome, results, outcomes 1agcbdei
383 @staticmethod
384 def compute_dataset_qc_status(outcomes, namespace='task'):
385 """Return map of dataset specific QC values.
387 Parameters
388 ----------
389 outcomes : dict
390 Map of checks and their individual outcomes.
391 namespace : str
392 The namespace of the QC fields in the Alyx JSON field.
394 Returns
395 -------
396 dict
397 Map of dataset names and their outcome.
398 """
399 trials_table_outcomes = { 1thc
400 'intervals': outcomes.get(f'_{namespace}_iti_delays', spec.QC.NOT_SET),
401 'goCue_times': outcomes.get(f'_{namespace}_goCue_delays', spec.QC.NOT_SET),
402 'response_times': spec.QC.NOT_SET, 'choice': spec.QC.NOT_SET,
403 'stimOn_times': outcomes.get(f'_{namespace}_stimOn_delays', spec.QC.NOT_SET),
404 'contrastLeft': spec.QC.NOT_SET, 'contrastRight': spec.QC.NOT_SET,
405 'feedbackType': spec.QC.NOT_SET, 'probabilityLeft': spec.QC.NOT_SET,
406 'feedback_times': outcomes.get(f'_{namespace}_errorCue_delays', spec.QC.NOT_SET),
407 'firstMovement_times': spec.QC.NOT_SET
408 }
409 reward_checks = (f'_{namespace}_reward_volumes', f'_{namespace}_reward_volume_set') 1thc
410 trials_table_outcomes['rewardVolume']: TaskQC.overall_outcome( 1thc
411 (outcomes.get(x, spec.QC.NOT_SET) for x in reward_checks)
412 )
413 dataset_outcomes = { 1thc
414 '_ibl_trials.stimOff_times': outcomes.get(f'_{namespace}_stimOff_delays', spec.QC.NOT_SET),
415 '_ibl_trials.table': trials_table_outcomes,
416 }
417 return dataset_outcomes 1thc
420class HabituationQC(TaskQC):
421 """Task QC for habituation choice world."""
423 def compute(self, **kwargs):
424 """Compute and store the QC metrics.
426 Runs the QC on the session and stores a map of the metrics for each datapoint for each
427 test, and a map of which datapoints passed for each test.
428 """
429 assert self.extractor is not None 1gb
430 self.log.info(f'Session {self.session_path}: Running QC on habituation data...') 1gb
432 # Initialize checks
433 prefix = f'_{kwargs.pop("namespace", self.namespace)}_' 1gb
434 data = self.extractor.data 1gb
435 audio_output = self.extractor.settings.get('device_sound', {}).get('OUTPUT', 'unknown') 1gb
436 metrics = {} 1gb
437 passed = {} 1gb
439 # Modify criteria based on version
440 ver = self.extractor.settings.get('IBLRIG_VERSION', '') or '0.0.0' 1gb
441 is_v8 = version.parse(ver) >= version.parse('8.0.0') 1gb
442 self.criteria['_task_iti_delays'] = {'PASS': 0.99, 'WARNING': 0} if is_v8 else {'NOT_SET': 0} 1gb
444 # Check all reward volumes == 3.0ul
445 check = prefix + 'reward_volumes' 1gb
446 metrics[check] = data['rewardVolume'] 1gb
447 passed[check] = metrics[check] == 3.0 1gb
449 # Check session durations are increasing in steps >= 12 minutes
450 check = prefix + 'habituation_time' 1gb
451 if not self.one or not self.session_path: 1gb
452 self.log.warning('unable to determine session trials without ONE')
453 metrics[check] = passed[check] = None
454 else:
455 subject, session_date = self.session_path.parts[-3:-1] 1gb
456 # compute from the date specified
457 date_minus_week = ( 1gb
458 datetime.strptime(session_date, '%Y-%m-%d') - timedelta(days=7)
459 ).strftime('%Y-%m-%d')
460 sessions = self.one.alyx.rest('sessions', 'list', subject=subject, 1gb
461 date_range=[date_minus_week, session_date],
462 task_protocol='habituation')
463 # Remove the current session if already registered
464 if sessions and sessions[0]['start_time'].startswith(session_date): 1gb
465 sessions = sessions[1:]
466 metric = ([0, data['intervals'][-1, 1] - data['intervals'][0, 0]] + 1gb
467 [(datetime.fromisoformat(x['end_time']) -
468 datetime.fromisoformat(x['start_time'])).total_seconds()
469 for x in [self.one.alyx.get(s['url']) for s in sessions]])
471 # The duration from raw trial data
472 # duration = map(float, self.extractor.raw_data[-1]['elapsed_time'].split(':'))
473 # duration = timedelta(**dict(zip(('hours', 'minutes', 'seconds'),
474 # duration))).total_seconds()
475 metrics[check] = np.array(metric) / 60 1gb
476 passed[check] = np.diff(metric) >= 12 1gb
478 # Check event orders: trial_start < stim on < stim center < feedback < stim off
479 check = prefix + 'trial_event_sequence' 1gb
480 nans = ( 1gb
481 np.isnan(data['intervals'][:, 0]) | # noqa
482 np.isnan(data['stimOn_times']) | # noqa
483 np.isnan(data['stimCenter_times']) |
484 np.isnan(data['valveOpen_times']) | # noqa
485 np.isnan(data['stimOff_times'])
486 )
487 a = np.less(data['intervals'][:, 0], data['stimOn_times'], where=~nans) 1gb
488 b = np.less(data['stimOn_times'], data['stimCenter_times'], where=~nans) 1gb
489 c = np.less(data['stimCenter_times'], data['valveOpen_times'], where=~nans) 1gb
490 d = np.less(data['valveOpen_times'], data['stimOff_times'], where=~nans) 1gb
492 metrics[check] = a & b & c & d & ~nans 1gb
493 passed[check] = metrics[check].astype(float) 1gb
495 # Check that the time difference between the visual stimulus center-command being
496 # triggered and the stimulus effectively appearing in the center is smaller than 150 ms.
497 check = prefix + 'stimCenter_delays' 1gb
498 metric = np.nan_to_num(data['stimCenter_times'] - data['stimCenterTrigger_times'], 1gb
499 nan=np.inf)
500 passed[check] = (metric <= 0.15) & (metric > 0) 1gb
501 metrics[check] = metric 1gb
503 # Phase check
504 check = prefix + 'phase' 1gb
505 metric = data['phase'] 1gb
506 passed[check] = (metric <= 2 * np.pi) & (metric >= 0) 1gb
507 metrics[check] = metric 1gb
509 # This is not very useful as a check because there are so few trials
510 check = prefix + 'phase_distribution' 1gb
511 metric, _ = np.histogram(data['phase']) 1gb
512 _, p = chisquare(metric) 1gb
513 passed[check] = p < 0.05 if len(data['phase']) >= 400 else None # skip if too few trials 1gb
514 metrics[check] = metric 1gb
516 # Check that the period of gray screen between stim off and the start of the next trial is
517 # 1s +/- 10%.
518 check = prefix + 'iti_delays' 1gb
519 iti = (np.roll(data['stimOn_times'], -1) - data['stimOff_times'])[:-1] 1gb
520 metric = np.r_[np.nan_to_num(iti, nan=np.inf), np.nan] - 1. 1gb
521 passed[check] = np.abs(metric) <= 0.1 1gb
522 passed[check][-1] = np.nan 1gb
523 metrics[check] = metric 1gb
525 # Check stim on go cue delay
526 # Unlike in ChoiceWorld, the go cue is in a race condition with the stimulus onset and
527 # therefore usually presents slightly before the stimulus. For now we will keep the QC
528 # thresholds the same as in ChoiceWorld which means this will always fail
529 check = prefix + 'stimOn_goCue_delays' 1gb
530 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
531 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1gb
532 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 1gb
533 passed[check] = np.abs(metric) < threshold 1gb
534 metrics[check] = metric 1gb
536 # Checks common to training QC
537 checks = [check_goCue_delays, check_stimOn_delays, check_stimOff_delays] 1gb
538 for fcn in checks: 1gb
539 check = prefix + fcn.__name__[6:] 1gb
540 metrics[check], passed[check] = fcn(data, audio_output=audio_output) 1gb
542 self.metrics, self.passed = (metrics, passed) 1gb
545# SINGLE METRICS
546# ---------------------------------------------------------------------------- #
548# === Delays between events checks ===
550def check_stimOn_goCue_delays(data, audio_output='harp', **_):
551 """
552 Check the go cue tone occurs less than 10ms before stimulus on.
554 Checks that the time difference between the onset of the visual stimulus
555 and the onset of the go cue tone is positive and less than 10ms.
557 Metric:
558 M = stimOn_times - goCue_times
560 Criteria:
561 0 < M < 0.010 s
563 Units:
564 seconds [s]
566 :param data: dict of trial data with keys ('goCue_times', 'stimOn_times', 'intervals')
567 :param audio_output: audio output device name.
569 Notes
570 -----
571 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
572 percentile of delays over 500 training sessions using the Xonar soundcard.
573 """
574 # Calculate the difference between stimOn and goCue times.
575 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
576 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1aBcbdfe
577 metric = np.nan_to_num(data['goCue_times'] - data['stimOn_times'], nan=np.inf) 1aBcbdfe
578 passed = (metric < threshold) & (metric > 0) 1aBcbdfe
579 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aBcbdfe
580 return metric, passed 1aBcbdfe
583def check_response_feedback_delays(data, audio_output='harp', **_):
584 """
585 Check the feedback delivered within 10ms of the response threshold.
587 Checks that the time difference between the response and the feedback onset
588 (error sound or valve) is positive and less than 10ms.
590 Metric:
591 M = feedback_time - response_time
593 Criterion:
594 0 < M < 0.010 s
596 Units:
597 seconds [s]
599 :param data: dict of trial data with keys ('feedback_times', 'response_times', 'intervals')
600 :param audio_output: audio output device name.
602 Notes
603 -----
604 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
605 percentile of delays over 500 training sessions using the Xonar soundcard.
606 """
607 threshold = 0.01 if audio_output.lower() == 'harp' else 0.053 1aCcbdfe
608 metric = np.nan_to_num(data['feedback_times'] - data['response_times'], nan=np.inf) 1aCcbdfe
609 passed = (metric < threshold) & (metric > 0) 1aCcbdfe
610 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aCcbdfe
611 return metric, passed 1aCcbdfe
614def check_response_stimFreeze_delays(data, **_):
615 """
616 Check the stimulus freezes within 100ms of the expected time.
618 Checks that the time difference between the visual stimulus freezing and the
619 response is positive and less than 100ms.
621 Metric:
622 M = (stimFreeze_times - response_times)
624 Criterion:
625 0 < M < 0.100 s
627 Units:
628 seconds [s]
630 :param data: dict of trial data with keys ('stimFreeze_times', 'response_times', 'intervals',
631 'choice')
632 """
633 # Calculate the difference between stimOn and goCue times.
634 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
635 metric = np.nan_to_num(data['stimFreeze_times'] - data['response_times'], nan=np.inf) 1aDcbdfe
636 # Test for valid values
637 passed = ((metric < 0.1) & (metric > 0)).astype(float) 1aDcbdfe
638 # Finally remove no_go trials (stimFreeze triggered differently in no_go trials)
639 # These values are ignored in calculation of proportion passed
640 passed[data['choice'] == 0] = np.nan 1aDcbdfe
641 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aDcbdfe
642 return metric, passed 1aDcbdfe
645def check_stimOff_itiIn_delays(data, **_):
646 """Check that the start of the trial interval is within 10ms of the visual stimulus turning off.
648 Metric:
649 M = itiIn_times - stimOff_times
651 Criterion:
652 0 < M < 0.010 s
654 Units:
655 seconds [s]
657 :param data: dict of trial data with keys ('stimOff_times', 'itiIn_times', 'intervals',
658 'choice')
659 """
660 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
661 metric = np.nan_to_num(data['itiIn_times'] - data['stimOff_times'], nan=np.inf) 1aEcbdfe
662 passed = ((metric < 0.01) & (metric >= 0)).astype(float) 1aEcbdfe
663 # Remove no_go trials (stimOff triggered differently in no_go trials)
664 # NaN values are ignored in calculation of proportion passed
665 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan 1aEcbdfe
666 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aEcbdfe
667 return metric, passed 1aEcbdfe
670def check_iti_delays(data, subtract_pauses=False, iti_delay_secs=ITI_DELAY_SECS,
671 feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_):
672 """
673 Check the open-loop grey screen period is approximately 1 second.
675 Check that the period of grey screen between stim off and the start of the next trial is
676 1s +/- 10%. If the trial was paused during this time, the check will account for that
678 Metric:
679 M = stimOff (n) - trialStart (n+1) - 1.
681 Criterion:
682 |M| < 0.1
684 Units:
685 seconds [s]
687 Parameters
688 ----------
689 data : dict
690 Trial data with keys ('stimOff_times', 'intervals', 'pause_duration').
691 subtract_pauses: bool
692 If True, account for experimenter-initiated pauses between trials; if False, trials where
693 the experimenter paused the task may fail this check.
695 Returns
696 -------
697 numpy.array
698 An array of metric values to threshold.
699 numpy.array
700 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
701 """
702 # Initialize array the length of completed trials
703 metric = np.full(data['intervals'].shape[0], np.nan) 1aucbdfe
704 passed = metric.copy() 1aucbdfe
705 pauses = (data['pause_duration'] if subtract_pauses else np.zeros_like(metric))[:-1] 1aucbdfe
706 # Get the difference between stim off and the start of the next trial
707 # Missing data are set to Inf, except for the last trial which is a NaN
708 metric[:-1] = np.nan_to_num( 1aucbdfe
709 data['intervals'][1:, 0] - data['stimOff_times'][:-1] - iti_delay_secs - pauses,
710 nan=np.inf
711 )
712 metric[data['choice'] == 0] = metric[data['choice'] == 0] - feedback_nogo_delay_secs 1aucbdfe
713 passed[:-1] = np.abs(metric[:-1]) < (iti_delay_secs / 10) # Last trial is not counted 1aucbdfe
714 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aucbdfe
715 return metric, passed 1aucbdfe
718def check_positive_feedback_stimOff_delays(data, **_):
719 """
720 Check stimulus offset occurs approximately 1 second after reward delivered.
722 Check that the time difference between the valve onset and the visual stimulus turning off
723 is 1 ± 0.150 seconds.
725 Metric:
726 M = stimOff_times - feedback_times - 1s
728 Criterion:
729 |M| < 0.150 s
731 Units:
732 seconds [s]
734 :param data: dict of trial data with keys ('stimOff_times', 'feedback_times', 'intervals',
735 'correct')
736 """
737 # If either are NaN, the result will be Inf to ensure that it crosses the failure threshold.
738 metric = np.nan_to_num(data['stimOff_times'] - data['feedback_times'] - 1, nan=np.inf) 1aFcbdfe
739 passed = (np.abs(metric) < 0.15).astype(float) 1aFcbdfe
740 # NaN values are ignored in calculation of proportion passed; ignore incorrect trials here
741 metric[~data['correct']] = passed[~data['correct']] = np.nan 1aFcbdfe
742 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aFcbdfe
743 return metric, passed 1aFcbdfe
746def check_negative_feedback_stimOff_delays(data, feedback_nogo_delay_secs=FEEDBACK_NOGO_DELAY_SECS, **_):
747 """
748 Check the stimulus offset occurs approximately 2 seconds after negative feedback delivery.
750 Check that the time difference between the error sound and the visual stimulus
751 turning off is 2 ± 0.150 seconds.
753 Metric:
754 M = stimOff_times - errorCue_times - 2s
756 Criterion:
757 |M| < 0.150 s
759 Units:
760 seconds [s]
762 :param data: dict of trial data with keys ('stimOff_times', 'errorCue_times', 'intervals')
763 """
764 metric = np.nan_to_num(data['stimOff_times'] - data['errorCue_times'] - 2, nan=np.inf) 1azcbdfe
765 # for the nogo trials, the feedback is the same as the stimOff
766 metric[data['choice'] == 0] = metric[data['choice'] == 0] + feedback_nogo_delay_secs 1azcbdfe
767 # Apply criteria
768 passed = (np.abs(metric) < 0.15).astype(float) 1azcbdfe
769 # Remove none negative feedback trials
770 metric[data['correct']] = passed[data['correct']] = np.nan 1azcbdfe
771 assert data['intervals'].shape[0] == len(metric) == len(passed) 1azcbdfe
772 return metric, passed 1azcbdfe
775# === Wheel movement during trial checks ===
777def check_wheel_move_before_feedback(data, **_):
778 """Check that the wheel does move within 100ms of the feedback onset (error sound or valve).
780 Metric:
781 M = (w_t - 0.05) - (w_t + 0.05), where t = feedback_times
783 Criterion:
784 M != 0
786 Units:
787 radians
789 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice',
790 'intervals', 'feedback_times')
791 """
792 # Get tuple of wheel times and positions within 100ms of feedback
793 traces = traces_by_trial( 1aocbdfe
794 data['wheel_timestamps'],
795 data['wheel_position'],
796 start=data['feedback_times'] - 0.05,
797 end=data['feedback_times'] + 0.05,
798 )
799 metric = np.zeros_like(data['feedback_times']) 1aocbdfe
800 # For each trial find the displacement
801 for i, trial in enumerate(traces): 1aocbdfe
802 pos = trial[1] 1aocbdfe
803 if pos.size > 1: 1aocbdfe
804 metric[i] = pos[-1] - pos[0] 1aocbdfe
806 # except no-go trials
807 metric[data['choice'] == 0] = np.nan # NaN = trial ignored for this check 1aocbdfe
808 nans = np.isnan(metric) 1aocbdfe
809 passed = np.zeros_like(metric) * np.nan 1aocbdfe
811 passed[~nans] = (metric[~nans] != 0).astype(float) 1aocbdfe
812 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aocbdfe
813 return metric, passed 1aocbdfe
816def _wheel_move_during_closed_loop(re_ts, re_pos, data, wheel_gain=None, tol=1, **_):
817 """
818 Check the wheel moves the correct amount to reach threshold.
820 Check that the wheel moves by approximately 35 degrees during the closed-loop period
821 on trials where a feedback (error sound or valve) is delivered.
823 Metric:
824 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
825 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
826 move 35 visual degrees
828 Criterion:
829 displacement < tol visual degree
831 Units:
832 degrees angle of wheel turn
834 :param re_ts: extracted wheel timestamps in seconds
835 :param re_pos: extracted wheel positions in radians
836 :param data: a dict with the keys (goCueTrigger_times, response_times, feedback_times,
837 position, choice, intervals)
838 :param wheel_gain: the 'STIM_GAIN' task setting
839 :param tol: the criterion in visual degrees
840 """
841 if wheel_gain is None: 1ajcbdfe
842 _log.warning('No wheel_gain input in function call, returning None')
843 return None, None
845 # Get tuple of wheel times and positions over each trial's closed-loop period
846 traces = traces_by_trial(re_ts, re_pos, 1ajcbdfe
847 start=data['goCueTrigger_times'],
848 end=data['response_times'])
850 metric = np.zeros_like(data['feedback_times']) 1ajcbdfe
851 # For each trial find the absolute displacement
852 for i, trial in enumerate(traces): 1ajcbdfe
853 t, pos = trial 1ajcbdfe
854 if pos.size != 0: 1ajcbdfe
855 # Find the position of the preceding sample and subtract it
856 idx = np.abs(re_ts - t[0]).argmin() - 1 1ajcbdfe
857 origin = re_pos[idx] 1ajcbdfe
858 metric[i] = np.abs(pos - origin).max() 1ajcbdfe
860 # Load wheel_gain and thresholds for each trial
861 wheel_gain = np.array([wheel_gain] * len(data['position'])) 1ajcbdfe
862 thresh = data['position'] 1ajcbdfe
863 # abs displacement, s, in mm required to move 35 visual degrees
864 s_mm = np.abs(thresh / wheel_gain) # don't care about direction 1ajcbdfe
865 criterion = cm_to_rad(s_mm * 1e-1) # convert abs displacement to radians (wheel pos is in rad) 1ajcbdfe
866 metric = metric - criterion # difference should be close to 0 1ajcbdfe
867 rad_per_deg = cm_to_rad(1 / wheel_gain * 1e-1) 1ajcbdfe
868 passed = (np.abs(metric) < rad_per_deg * tol).astype(float) # less than 1 visual degree off 1ajcbdfe
869 metric[data['choice'] == 0] = passed[data['choice'] == 0] = np.nan # except no-go trials 1ajcbdfe
870 assert data['intervals'].shape[0] == len(metric) == len(passed) 1ajcbdfe
871 return metric, passed 1ajcbdfe
874def check_wheel_move_during_closed_loop(data, wheel_gain=None, **_):
875 """
876 Check the wheel moves the correct amount to reach threshold.
878 Check that the wheel moves by approximately 35 degrees during the closed-loop period
879 on trials where a feedback (error sound or valve) is delivered.
881 Metric:
882 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
883 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
884 move 35 visual degrees
886 Criterion:
887 displacement < 3 visual degrees
889 Units:
890 degrees angle of wheel turn
892 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'choice',
893 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position')
894 :param wheel_gain: the 'STIM_GAIN' task setting
895 """
896 # Get the Bpod extracted wheel data
897 timestamps = data['wheel_timestamps'] 1ajcbdfe
898 position = data['wheel_position'] 1ajcbdfe
900 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=3) 1ajcbdfe
903def check_wheel_move_during_closed_loop_bpod(data, wheel_gain=None, **_):
904 """
905 Check the wheel moves the correct amount to reach threshold.
907 Check that the wheel moves by approximately 35 degrees during the closed-loop period
908 on trials where a feedback (error sound or valve) is delivered. This check uses the Bpod
909 wheel data (measured at a lower resolution) with a stricter tolerance (1 visual degree).
911 Metric:
912 M = abs(w_resp - w_t0) - threshold_displacement, where w_resp = position at response
913 time, w_t0 = position at go cue time, threshold_displacement = displacement required to
914 move 35 visual degrees.
916 Criterion:
917 displacement < 1 visual degree
919 Units:
920 degrees angle of wheel turn
922 :param data: dict of trial data with keys ('wheel_timestamps(_bpod)', 'wheel_position(_bpod)',
923 'choice', 'intervals', 'goCueTrigger_times', 'response_times', 'feedback_times', 'position')
924 :param wheel_gain: the 'STIM_GAIN' task setting
925 """
926 # Get the Bpod extracted wheel data
927 timestamps = data.get('wheel_timestamps_bpod', data['wheel_timestamps']) 1acbdfe
928 position = data.get('wheel_position_bpod', data['wheel_position']) 1acbdfe
930 return _wheel_move_during_closed_loop(timestamps, position, data, wheel_gain, tol=1) 1acbdfe
933def check_wheel_freeze_during_quiescence(data, **_):
934 """
935 Check the wheel is indeed still during the quiescent period.
937 Check that the wheel does not move more than 2 degrees in each direction during the
938 quiescence interval before the stimulus appears.
940 Metric:
941 M = |max(W) - min(W)| where W is wheel pos over quiescence interval;
942 interval = [stimOnTrigger_times - quiescent_duration, stimOnTrigger_times]
944 Criterion:
945 M < 2 degrees
947 Units:
948 degrees angle of wheel turn
950 :param data: dict of trial data with keys ('wheel_timestamps', 'wheel_position', 'quiescence',
951 'intervals', 'stimOnTrigger_times')
952 """
953 assert np.all(np.diff(data['wheel_timestamps']) >= 0) 1akcbdfe
954 assert data['quiescence'].size == data['stimOnTrigger_times'].size 1akcbdfe
955 # Get tuple of wheel times and positions over each trial's quiescence period
956 qevt_start_times = data['stimOnTrigger_times'] - data['quiescence'] 1akcbdfe
957 traces = traces_by_trial( 1akcbdfe
958 data['wheel_timestamps'],
959 data['wheel_position'],
960 start=qevt_start_times,
961 end=data['stimOnTrigger_times']
962 )
964 metric = np.zeros((len(data['quiescence']), 2)) # (n_trials, n_directions) 1akcbdfe
965 for i, trial in enumerate(traces): 1akcbdfe
966 t, pos = trial 1akcbdfe
967 # Get the last position before the period began
968 if pos.size > 0: 1akcbdfe
969 # Find the position of the preceding sample and subtract it
970 idx = np.abs(data['wheel_timestamps'] - t[0]).argmin() - 1 1akcbdfe
971 origin = data['wheel_position'][idx if idx != -1 else 0] 1akcbdfe
972 # Find the absolute min and max relative to the last sample
973 metric[i, :] = np.abs([np.min(pos - origin), np.max(pos - origin)]) 1akcbdfe
974 # Reduce to the largest displacement found in any direction
975 metric = np.max(metric, axis=1) 1akcbdfe
976 metric = 180 * metric / np.pi # convert to degrees from radians 1akcbdfe
977 criterion = 2 # Position shouldn't change more than 2 in either direction 1akcbdfe
978 passed = metric < criterion 1akcbdfe
979 assert data['intervals'].shape[0] == len(metric) == len(passed) 1akcbdfe
980 return metric, passed 1akcbdfe
983def check_detected_wheel_moves(data, min_qt=0, **_):
984 """Check that the detected first movement times are reasonable.
986 Metric:
987 M = firstMovement times
989 Criterion:
990 (goCue trigger time - min quiescent period) < M < response time
992 Units:
993 Seconds [s]
995 :param data: dict of trial data with keys ('firstMovement_times', 'goCueTrigger_times',
996 'response_times', 'choice', 'intervals')
997 :param min_qt: the minimum possible quiescent period (the QUIESCENT_PERIOD task parameter)
998 """
999 # Depending on task version this may be a single value or an array of quiescent periods
1000 min_qt = np.array(min_qt) 1ascbdfe
1001 if min_qt.size > data['intervals'].shape[0]: 1ascbdfe
1002 min_qt = min_qt[:data['intervals'].shape[0]] 1d
1004 metric = data['firstMovement_times'] 1ascbdfe
1005 qevt_start = data['goCueTrigger_times'] - np.array(min_qt) 1ascbdfe
1006 response = data['response_times'] 1ascbdfe
1007 # First movement time for each trial should be after the quiescent period and before feedback
1008 passed = np.array([a < m < b for m, a, b in zip(metric, qevt_start, response)], dtype=float) 1ascbdfe
1009 nogo = data['choice'] == 0 1ascbdfe
1010 passed[nogo] = np.nan # No go trial may have no movement times and that's fine 1ascbdfe
1011 return metric, passed 1ascbdfe
1014# === Sequence of events checks ===
1016def check_error_trial_event_sequence(data, **_):
1017 """
1018 Check trial events occur in correct order for negative feedback trials.
1020 Check that on incorrect / miss trials, there are exactly:
1021 2 audio events (go cue sound and error sound) and 2 Bpod events (trial start, ITI), occurring
1022 in the correct order
1024 Metric:
1025 M = Bpod (trial start) > audio (go cue) > audio (error) > Bpod (ITI) > Bpod (trial end)
1027 Criterion:
1028 M == True
1030 Units:
1031 -none-
1033 :param data: dict of trial data with keys ('errorCue_times', 'goCue_times', 'intervals',
1034 'itiIn_times', 'correct')
1035 """
1036 # An array the length of N trials where True means at least one event time was NaN (bad)
1037 nans = ( 1aqcbdfe
1038 np.isnan(data['intervals'][:, 0]) |
1039 np.isnan(data['goCue_times']) | # noqa
1040 np.isnan(data['errorCue_times']) | # noqa
1041 np.isnan(data['itiIn_times']) | # noqa
1042 np.isnan(data['intervals'][:, 1])
1043 )
1045 # For each trial check that the events happened in the correct order (ignore NaN values)
1046 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1aqcbdfe
1047 b = np.less(data['goCue_times'], data['errorCue_times'], where=~nans) # Go cue < error cue 1aqcbdfe
1048 c = np.less(data['errorCue_times'], data['itiIn_times'], where=~nans) # Error cue < ITI start 1aqcbdfe
1049 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1aqcbdfe
1051 # For each trial check all events were in order AND all event times were not NaN
1052 metric = a & b & c & d & ~nans 1aqcbdfe
1054 passed = metric.astype(float) 1aqcbdfe
1055 passed[data['correct']] = np.nan # Look only at incorrect trials 1aqcbdfe
1056 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aqcbdfe
1057 return metric, passed 1aqcbdfe
1060def check_correct_trial_event_sequence(data, **_):
1061 """
1062 Check trial events occur in correct order for positive feedback trials.
1064 Check that on correct trials, there are exactly:
1065 1 audio events and 3 Bpod events (valve open, trial start, ITI), occurring in the correct order
1067 Metric:
1068 M = Bpod (trial start) > audio (go cue) > Bpod (valve) > Bpod (ITI) > Bpod (trial end)
1070 Criterion:
1071 M == True
1073 Units:
1074 -none-
1076 :param data: dict of trial data with keys ('valveOpen_times', 'goCue_times', 'intervals',
1077 'itiIn_times', 'correct')
1078 """
1079 # An array the length of N trials where True means at least one event time was NaN (bad)
1080 nans = ( 1arcbdfe
1081 np.isnan(data['intervals'][:, 0]) |
1082 np.isnan(data['goCue_times']) | # noqa
1083 np.isnan(data['valveOpen_times']) |
1084 np.isnan(data['itiIn_times']) | # noqa
1085 np.isnan(data['intervals'][:, 1])
1086 )
1088 # For each trial check that the events happened in the correct order (ignore NaN values)
1089 a = np.less(data['intervals'][:, 0], data['goCue_times'], where=~nans) # Start time < go cue 1arcbdfe
1090 b = np.less(data['goCue_times'], data['valveOpen_times'], where=~nans) # Go cue < feedback 1arcbdfe
1091 c = np.less(data['valveOpen_times'], data['itiIn_times'], where=~nans) # Feedback < ITI start 1arcbdfe
1092 d = np.less(data['itiIn_times'], data['intervals'][:, 1], where=~nans) # ITI start < end time 1arcbdfe
1094 # For each trial True means all events were in order AND all event times were not NaN
1095 metric = a & b & c & d & ~nans 1arcbdfe
1097 passed = metric.astype(float) 1arcbdfe
1098 passed[~data['correct']] = np.nan # Look only at correct trials 1arcbdfe
1099 assert data['intervals'].shape[0] == len(metric) == len(passed) 1arcbdfe
1100 return metric, passed 1arcbdfe
1103def check_n_trial_events(data, **_):
1104 """Check that the number events per trial is correct.
1106 Within every trial interval there should be one of each trial event, except for
1107 goCueTrigger_times which should only be defined for incorrect trials
1109 Metric:
1110 M = all(start < event < end) for all event times except errorCueTrigger_times where
1111 start < error_trigger < end if not correct trial, else error_trigger == NaN
1113 Criterion:
1114 M == True
1116 Units:
1117 -none-, boolean
1119 :param data: dict of trial data with keys ('intervals', 'stimOnTrigger_times',
1120 'stimOffTrigger_times', 'stimOn_times', 'stimOff_times',
1121 'stimFreezeTrigger_times', 'errorCueTrigger_times', 'itiIn_times',
1122 'goCueTrigger_times', 'goCue_times', 'response_times', 'feedback_times')
1123 """
1124 intervals = data['intervals'] 1apcbdfe
1125 correct = data['correct'] 1apcbdfe
1126 err_trig = data['errorCueTrigger_times'] 1apcbdfe
1128 # Exclude these fields; valve and errorCue times are the same as feedback_times and we must
1129 # test errorCueTrigger_times separately
1130 # stimFreeze_times fails often due to TTL flicker
1131 exclude = ['camera_timestamps', 'errorCueTrigger_times', 'errorCue_times', 1apcbdfe
1132 'wheelMoves_peakVelocity_times', 'valveOpen_times', 'wheelMoves_peakAmplitude',
1133 'wheelMoves_intervals', 'wheel_timestamps', 'stimFreeze_times']
1134 events = [k for k in data.keys() if k.endswith('_times') and k not in exclude] 1apcbdfe
1135 metric = np.zeros(data['intervals'].shape[0], dtype=bool) 1apcbdfe
1137 # For each trial interval check that one of each trial event occurred. For incorrect trials,
1138 # check the error cue trigger occurred within the interval, otherwise check it is nan.
1139 for i, (start, end) in enumerate(intervals): 1apcbdfe
1140 metric[i] = (all([start < data[k][i] < end for k in events]) and 1apcbdfe
1141 (np.isnan(err_trig[i]) if correct[i] else start < err_trig[i] < end))
1142 passed = metric.astype(bool) 1apcbdfe
1143 assert intervals.shape[0] == len(metric) == len(passed) 1apcbdfe
1144 return metric, passed 1apcbdfe
1147def check_trial_length(data, **_):
1148 """
1149 Check open-loop duration positive and <= 1 minute.
1151 Check that the time difference between the onset of the go cue sound
1152 and the feedback (error sound or valve) is positive and smaller than 60.1 s.
1154 Metric:
1155 M = feedback_times - goCue_times
1157 Criteria:
1158 0 < M < 60.1 s
1160 Units:
1161 seconds [s]
1163 :param data: dict of trial data with keys ('feedback_times', 'goCue_times', 'intervals')
1164 """
1165 # NaN values are usually ignored so replace them with Inf so they fail the threshold
1166 metric = np.nan_to_num(data['feedback_times'] - data['goCue_times'], nan=np.inf) 1aHcbdfe
1167 passed = (metric < 60.1) & (metric > 0) 1aHcbdfe
1168 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aHcbdfe
1169 return metric, passed 1aHcbdfe
1172# === Trigger-response delay checks ===
1174def check_goCue_delays(data, audio_output='harp', **_):
1175 """
1176 Check the go cue tone occurs within 1ms of the intended time.
1178 Check that the time difference between the go cue sound being triggered and
1179 effectively played is smaller than 1ms.
1181 Metric:
1182 M = goCue_times - goCueTrigger_times
1184 Criterion:
1185 0 < M <= 0.0015 s
1187 Units:
1188 seconds [s]
1190 :param data: dict of trial data with keys ('goCue_times', 'goCueTrigger_times', 'intervals').
1191 :param audio_output: audio output device name.
1193 Notes
1194 -----
1195 - For non-harp sound card the permissible delay is 0.053s. This was chosen by taking the 99.5th
1196 percentile of delays over 500 training sessions using the Xonar soundcard.
1197 """
1198 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.053 1agGcbdfe
1199 metric = np.nan_to_num(data['goCue_times'] - data['goCueTrigger_times'], nan=np.inf) 1agGcbdfe
1200 passed = (metric <= threshold) & (metric > 0) 1agGcbdfe
1201 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agGcbdfe
1202 return metric, passed 1agGcbdfe
1205def check_errorCue_delays(data, audio_output='harp', **_):
1206 """
1207 Check the error tone occurs within 1.5ms of the intended time.
1209 Check that the time difference between the error sound being triggered and
1210 effectively played is smaller than 1ms.
1212 Metric:
1213 M = errorCue_times - errorCueTrigger_times
1215 Criterion:
1216 0 < M <= 0.0015 s
1218 Units:
1219 seconds [s]
1221 :param data: dict of trial data with keys ('errorCue_times', 'errorCueTrigger_times',
1222 'intervals', 'correct')
1223 :param audio_output: audio output device name.
1225 Notes
1226 -----
1227 - For non-harp sound card the permissible delay is 0.062s. This was chosen by taking the 99.5th
1228 percentile of delays over 500 training sessions using the Xonar soundcard.
1229 """
1230 threshold = 0.0015 if audio_output.lower() == 'harp' else 0.062 1aAcbdfe
1231 metric = np.nan_to_num(data['errorCue_times'] - data['errorCueTrigger_times'], nan=np.inf) 1aAcbdfe
1232 passed = ((metric <= threshold) & (metric > 0)).astype(float) 1aAcbdfe
1233 passed[data['correct']] = metric[data['correct']] = np.nan 1aAcbdfe
1234 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aAcbdfe
1235 return metric, passed 1aAcbdfe
1238def check_stimOn_delays(data, **_):
1239 """
1240 Check the visual stimulus onset occurs within 150ms of the intended time.
1242 Check that the time difference between the visual stimulus onset-command being triggered
1243 and the stimulus effectively appearing on the screen is smaller than 150 ms.
1245 Metric:
1246 M = stimOn_times - stimOnTrigger_times
1248 Criterion:
1249 0 < M < 0.15 s
1251 Units:
1252 seconds [s]
1254 :param data: dict of trial data with keys ('stimOn_times', 'stimOnTrigger_times',
1255 'intervals')
1256 """
1257 metric = np.nan_to_num(data['stimOn_times'] - data['stimOnTrigger_times'], nan=np.inf) 1agIcbdfe
1258 passed = (metric <= 0.15) & (metric > 0) 1agIcbdfe
1259 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agIcbdfe
1260 return metric, passed 1agIcbdfe
1263def check_stimOff_delays(data, **_):
1264 """
1265 Check stimulus offset occurs within 150ms of the intended time.
1267 Check that the time difference between the visual stimulus offset-command
1268 being triggered and the visual stimulus effectively turning off on the screen
1269 is smaller than 150 ms.
1271 Metric:
1272 M = stimOff_times - stimOffTrigger_times
1274 Criterion:
1275 0 < M < 0.15 s
1277 Units:
1278 seconds [s]
1280 :param data: dict of trial data with keys ('stimOff_times', 'stimOffTrigger_times',
1281 'intervals')
1282 """
1283 metric = np.nan_to_num(data['stimOff_times'] - data['stimOffTrigger_times'], nan=np.inf) 1agJcbdfe
1284 passed = (metric <= 0.15) & (metric > 0) 1agJcbdfe
1285 assert data['intervals'].shape[0] == len(metric) == len(passed) 1agJcbdfe
1286 return metric, passed 1agJcbdfe
1289def check_stimFreeze_delays(data, **_):
1290 """Check the stimulus freezes within 150ms of the intended time.
1292 Check that the time difference between the visual stimulus freeze-command
1293 being triggered and the visual stimulus effectively freezing on the screen
1294 is smaller than 150 ms.
1296 Metric:
1297 M = stimFreeze_times - stimFreezeTrigger_times
1299 Criterion:
1300 0 < M < 0.15 s
1302 Units:
1303 seconds [s]
1305 :param data: dict of trial data with keys ('stimFreeze_times', 'stimFreezeTrigger_times',
1306 'intervals')
1307 """
1308 metric = np.nan_to_num(data['stimFreeze_times'] - data['stimFreezeTrigger_times'], nan=np.inf) 1aKcbdfe
1309 passed = (metric <= 0.15) & (metric > 0) 1aKcbdfe
1310 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aKcbdfe
1311 return metric, passed 1aKcbdfe
1314# === Data integrity checks ===
1316def check_reward_volumes(data, **_):
1317 """Check that the reward volume is between 1.5 and 3 uL for correct trials, 0 for incorrect.
1319 Metric:
1320 M = reward volume
1322 Criterion:
1323 1.5 <= M <= 3 if correct else M == 0
1325 Units:
1326 uL
1328 :param data: dict of trial data with keys ('rewardVolume', 'correct', 'intervals')
1329 """
1330 metric = data['rewardVolume'] 1aycbdfe
1331 correct = data['correct'] 1aycbdfe
1332 passed = np.zeros_like(metric, dtype=bool) 1aycbdfe
1333 # Check correct trials within correct range
1334 passed[correct] = (1.5 <= metric[correct]) & (metric[correct] <= 3.) 1aycbdfe
1335 # Check incorrect trials are 0
1336 passed[~correct] = metric[~correct] == 0 1aycbdfe
1337 assert data['intervals'].shape[0] == len(metric) == len(passed) 1aycbdfe
1338 return metric, passed 1aycbdfe
1341def check_reward_volume_set(data, **_):
1342 """Check that there is only two reward volumes within a session, one of which is 0.
1344 Metric:
1345 M = set(rewardVolume)
1347 Criterion:
1348 (0 < len(M) <= 2) and 0 in M
1350 :param data: dict of trial data with keys ('rewardVolume')
1351 """
1352 metric = data['rewardVolume'] 1aLcbdfe
1353 passed = 0 < len(set(metric)) <= 2 and 0. in metric 1aLcbdfe
1354 return metric, passed 1aLcbdfe
1357def check_wheel_integrity(data, re_encoding='X1', enc_res=None, **_):
1358 """
1359 Check wheel position sampled at the expected resolution.
1361 Check that the difference between wheel position samples is close to the encoder resolution
1362 and that the wheel timestamps strictly increase.
1364 Metric:
1365 M = (absolute difference of the positions < 1.5 * encoder resolution)
1366 + 1 if (difference of timestamps <= 0) else 0
1368 Criterion:
1369 M ~= 0
1371 Units:
1372 arbitrary (radians, sometimes + 1)
1374 :param data: dict of wheel data with keys ('wheel_timestamps', 'wheel_position')
1375 :param re_encoding: the encoding of the wheel data, X1, X2 or X4
1376 :param enc_res: the rotary encoder resolution (default 1024 ticks per revolution)
1378 Notes
1379 -----
1380 - At high velocities some samples are missed due to the scanning frequency of the DAQ.
1381 This checks for more than 1 missing sample in a row (i.e. the difference between samples >= 2)
1383 """
1384 if isinstance(re_encoding, str): 1avcbdfe
1385 re_encoding = int(re_encoding[-1]) 1avcbdfe
1386 # The expected difference between samples in the extracted units
1387 resolution = 1 / (enc_res or ephys_fpga.WHEEL_TICKS 1avcbdfe
1388 ) * np.pi * 2 * ephys_fpga.WHEEL_RADIUS_CM / re_encoding
1389 # We expect the difference of neighbouring positions to be close to the resolution
1390 pos_check = np.abs(np.diff(data['wheel_position'])) 1avcbdfe
1391 # Timestamps should be strictly increasing
1392 ts_check = np.diff(data['wheel_timestamps']) <= 0. 1avcbdfe
1393 metric = pos_check + ts_check.astype(float) # all values should be close to zero 1avcbdfe
1394 passed = metric < 1.5 * resolution 1avcbdfe
1395 return metric, passed 1avcbdfe
1398# === Pre-stimulus checks ===
1399def check_stimulus_move_before_goCue(data, photodiode=None, **_):
1400 """
1401 Check there are no stimulus events before the go cue tone.
1403 Check that there are no visual stimulus change(s) between the start of the trial and the
1404 go cue sound onset, except for stim on.
1406 Metric:
1407 M = number of visual stimulus change events between trial start and goCue_times
1409 Criterion:
1410 M == 1
1412 Units:
1413 -none-, integer
1415 Parameters
1416 ----------
1417 data : dict
1418 Trial data with keys ('goCue_times', 'intervals', 'choice').
1419 photodiode : dict
1420 The fronts from Bpod's BNC1 input or FPGA frame2ttl channel.
1422 Returns
1423 -------
1424 numpy.array
1425 An array of metric values to threshold.
1426 numpy.array
1427 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
1429 Notes
1430 -----
1431 - There should be exactly 1 stimulus change before goCue; stimulus onset. Even if the stimulus
1432 contrast is 0, the sync square will still flip at stimulus onset, etc.
1433 - If there are no goCue times (all are NaN), the status should be NOT_SET.
1434 """
1435 if photodiode is None: 1acbdfe
1436 _log.warning('No photodiode TTL input in function call, returning None')
1437 return None
1438 photodiode_clean = ephys_fpga._clean_frame2ttl(photodiode) 1acbdfe
1439 s = photodiode_clean['times'] 1acbdfe
1440 s = s[~np.isnan(s)] # Remove NaNs 1acbdfe
1441 metric = np.array([]) 1acbdfe
1442 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1acbdfe
1443 metric = np.append(metric, np.count_nonzero(s[s > i] < c)) 1acbdfe
1445 passed = (metric == 1).astype(float) 1acbdfe
1446 passed[np.isnan(data['goCue_times'])] = np.nan 1acbdfe
1447 assert data['intervals'].shape[0] == len(metric) == len(passed) 1acbdfe
1448 return metric, passed 1acbdfe
1451def check_audio_pre_trial(data, audio=None, **_):
1452 """
1453 Check no audio stimuli before the go cue.
1455 Check there are no audio outputs between the start of the trial and the go cue sound onset - 20 ms.
1457 Metric:
1458 M = sum(start_times < audio TTL < (goCue_times - 20ms))
1460 Criterion:
1461 M == 0
1463 Units:
1464 -none-, integer
1466 Parameters
1467 ----------
1468 data : dict
1469 Trial data with keys ('goCue_times', 'intervals').
1470 audio : dict
1471 The fronts from Bpod's BNC2 input FPGA audio sync channel.
1473 Returns
1474 -------
1475 numpy.array
1476 An array of metric values to threshold.
1477 numpy.array
1478 An array of boolean values, 1 per trial, where True means trial passes QC threshold.
1479 """
1480 if audio is None: 1awcbdfe
1481 _log.warning('No BNC2 input in function call, retuning None')
1482 return None, None
1483 s = audio['times'][~np.isnan(audio['times'])] # Audio TTLs with NaNs removed 1awcbdfe
1484 metric = np.array([], dtype=np.int8) 1awcbdfe
1485 for i, c in zip(data['intervals'][:, 0], data['goCue_times']): 1awcbdfe
1486 metric = np.append(metric, sum(s[s > i] < (c - 0.02))) 1awcbdfe
1487 passed = metric == 0 1awcbdfe
1488 assert data['intervals'].shape[0] == len(metric) == len(passed) 1awcbdfe
1489 return metric, passed 1awcbdfe