Coverage for ibllib/ephys/ephysqc.py: 48%
292 statements
« prev ^ index » next coverage.py v7.7.0, created at 2025-03-17 15:25 +0000
1"""
2Quality control of raw Neuropixel electrophysiology data.
3"""
4from pathlib import Path
5import logging
7import numpy as np
8import pandas as pd
9from scipy import signal, stats
10import one.alf.io as alfio
11from iblutil.util import Bunch
12import spikeglx
13import neuropixel
14from ibldsp import fourier, utils, voltage
15from tqdm import tqdm
17from brainbox.io.spikeglx import Streamer
18from brainbox.metrics.single_units import spike_sorting_metrics
19from ibllib.ephys import sync_probes, spikes
20from ibllib.qc import base
21from ibllib.io.extractors import ephys_fpga
22from phylib.io import model
25_logger = logging.getLogger(__name__)
27RMS_WIN_LENGTH_SECS = 3
28WELCH_WIN_LENGTH_SAMPLES = 1024
29NCH_WAVEFORMS = 32 # number of channels to be saved in templates.waveforms and channels.waveforms
30BATCHES_SPACING = 300
31TMIN = 40
32SAMPLE_LENGTH = 1
33SPIKE_THRESHOLD_UV = -50 # negative, the threshold used for spike detection on pre-processed raw data
36class EphysQC(base.QC):
37 """
38 A class for computing Ephys QC metrics.
40 :param probe_id: An existing and registered probe insertion ID.
41 :param one: An ONE instance pointing to the database the probe_id is registered with. Optional, will instantiate
42 default database if not given.
43 """
45 def __init__(self, probe_id, session_path=None, **kwargs):
46 self.use_alyx = kwargs.pop('use_alyx', True)
47 self.stream = kwargs.pop('stream', True)
49 if self.use_alyx:
50 super().__init__(probe_id, endpoint='insertions', **kwargs)
51 self._outcome = 'NOT_SET'
52 self.pid = probe_id
54 self.session_path = session_path
55 keys = ('ap', 'ap_meta', 'lf', 'lf_meta')
56 self.data = Bunch.fromkeys(keys)
57 self.metrics = {}
59 def _ensure_required_data(self):
60 """
61 Ensures the datasets required for QC are available locally or remotely.
62 """
63 assert self.one is not None, 'ONE instance is required to ensure required data' 1h
64 eid, pname = self.one.pid2eid(self.pid) 1h
65 if self.session_path is None: 1h
66 self.session_path = self.one.eid2path(eid) 1h
67 self.probe_path = Path(self.session_path).joinpath('raw_ephys_data', pname) 1h
68 # Check if there is at least one meta file available
69 meta_files = list(self.probe_path.rglob('*.meta')) 1h
70 assert len(meta_files) != 0, f'No meta files in {self.probe_path}' 1h
71 # Check if there is no more than one meta file per type
72 ap_meta = [meta for meta in meta_files if 'ap.meta' in meta.name] 1h
73 assert not len(ap_meta) > 1, f'More than one ap.meta file in {self.probe_path}. Remove redundant files to run QC' 1h
74 lf_meta = [meta for meta in meta_files if 'lf.meta' in meta.name] 1h
75 assert not len(lf_meta) > 1, f'More than one lf.meta file in {self.probe_path}. Remove redundant files to run QC' 1h
77 def load_data(self, ensure=True) -> None:
78 """
79 Load any locally available data.
80 """
81 # First sanity check
82 if self.use_alyx:
83 self._ensure_required_data()
85 _logger.info('Gathering data for QC')
86 # Load metadata and, if locally present, bin file
87 for dstype in ['ap', 'lf']:
88 # We already checked that there is not more than one meta file per type
89 meta_file = next(self.probe_path.rglob(f'*{dstype}.meta'), None)
90 if meta_file is None:
91 _logger.warning(f'No {dstype}.meta file in {self.probe_path}, skipping QC for {dstype} data.')
92 else:
93 self.data[f'{dstype}_meta'] = spikeglx.read_meta_data(meta_file)
94 bin_file = next(meta_file.parent.glob(f'*{dstype}.*bin'), None)
95 if not bin_file:
96 # we only stream the AP file, we won't stream the full LF file...
97 if dstype == 'ap':
98 self.data[f'{dstype}'] = Streamer(pid=self.pid, one=self.one, remove_cached=True)
99 else:
100 self.data[f'{dstype}'] = None
101 else:
102 self.data[f'{dstype}'] = spikeglx.Reader(bin_file, open=True)
104 @staticmethod
105 def _compute_metrics_array(raw, fs, h):
106 """
107 From a numpy array, computes rms on raw data, destripes, computes rms on destriped data
108 and performs a simple spike detection
109 :param raw: voltage numpy.array(ntraces, nsamples)
110 :param fs: sampling frequency (Hz)
111 :param h: dictionary containing sensor coordinates, see neuropixel.trace_header
112 :return: 3 numpy vectors nchannels length
113 """
114 destripe = voltage.destripe(raw, fs=fs, h=h)
115 rms_raw = utils.rms(raw)
116 rms_pre_proc = utils.rms(destripe)
117 detections = spikes.detection(data=destripe.T, fs=fs, h=h, detect_threshold=SPIKE_THRESHOLD_UV * 1e-6)
118 spike_rate = np.bincount(detections.trace, minlength=raw.shape[0]).astype(np.float32)
119 channel_labels, _ = voltage.detect_bad_channels(raw, fs=fs)
120 _, psd = signal.welch(destripe, fs=fs, window='hann', nperseg=WELCH_WIN_LENGTH_SAMPLES,
121 detrend='constant', return_onesided=True, scaling='density', axis=-1)
122 return rms_raw, rms_pre_proc, spike_rate, channel_labels, psd
124 def run(self, update: bool = False, overwrite: bool = True, stream: bool = None, **kwargs) -> (str, dict):
125 """
126 Run QC on samples of the .ap file, and on the entire file for .lf data if it is present.
128 :param update: bool, whether to update the qc json fields for this probe. Default is False.
129 :param overwrite: bool, whether to overwrite locally existing outputs of this function. Default is False.
130 :param stream: bool, whether to stream the samples of the .ap data if not locally available. Defaults to value
131 set in class init (True if none set).
132 :return: A list of QC output files. In case of a complete run that is one file for .ap and three files for .lf.
133 """
134 # If stream is explicitly given in run, overwrite value from init
135 if stream is not None:
136 self.stream = stream
138 # Load data
139 self.load_data()
140 self.out_path = kwargs.get('out_path', self.probe_path)
142 qc_files = []
143 # If ap meta file present, calculate median RMS per channel before and after destriping
144 # NB: ideally this should go a a separate function once we have a spikeglx.Streamer that behaves like the Reader
145 if self.data.ap_meta:
146 files = {'rms': self.out_path.joinpath("_iblqc_ephysChannels.apRMS.npy"),
147 'spike_rate': self.out_path.joinpath("_iblqc_ephysChannels.rawSpikeRates.npy"),
148 'channel_labels': self.out_path.joinpath("_iblqc_ephysChannels.labels.npy"),
149 'ap_freqs': self.out_path.joinpath("_iblqc_ephysSpectralDensityAP.freqs.npy"),
150 'ap_power': self.out_path.joinpath("_iblqc_ephysSpectralDensityAP.power.npy"),
151 }
152 if all([files[k].exists() for k in files]) and not overwrite:
153 _logger.warning(f'RMS map already exists for .ap data in {self.probe_path}, skipping. '
154 f'Use overwrite option.')
155 results = {k: np.load(files[k]) for k in files}
156 else:
157 sr = self.data['ap']
158 nc = sr.nc - sr.nsync
160 # verify that the channel layout is correct according to IBL layout
161 th = sr.geometry
162 if sr.meta.get('NP2.4_shank', None) is not None:
163 h = neuropixel.trace_header(sr.major_version, nshank=4)
164 h = neuropixel.split_trace_header(h, shank=int(sr.meta.get('NP2.4_shank')))
165 else:
166 h = neuropixel.trace_header(sr.major_version, nshank=np.unique(th['shank']).size)
168 if not (np.all(h['x'] == th['x']) and np.all(h['y'] == th['y'])):
169 _logger.critical("Channel geometry seems incorrect")
170 # raise ValueError("Wrong Neuropixel channel mapping used - ABORT")
172 t0s = np.arange(TMIN, sr.rl - SAMPLE_LENGTH, BATCHES_SPACING)
173 all_rms = np.zeros((2, nc, t0s.shape[0]))
174 all_srs, channel_ok = (np.zeros((nc, t0s.shape[0])) for _ in range(2))
175 psds = np.zeros((nc, fourier.fscale(WELCH_WIN_LENGTH_SAMPLES, 1, one_sided=True).size))
177 _logger.info(f'Computing RMS samples for .ap data {self.probe_path}')
178 for i, t0 in enumerate(t0s):
179 sl = slice(int(t0 * sr.fs), int((t0 + SAMPLE_LENGTH) * sr.fs))
180 raw = sr[sl, :-sr.nsync].T
181 all_rms[0, :, i], all_rms[1, :, i], all_srs[:, i], channel_ok[:, i], psd =\
182 self._compute_metrics_array(raw, sr.fs, h)
183 psds += psd
184 # Calculate the median RMS across all samples per channel
185 results = {'rms': np.median(all_rms, axis=-1),
186 'spike_rate': np.median(all_srs, axis=-1),
187 'channel_labels': stats.mode(channel_ok, axis=1)[0],
188 'ap_freqs': fourier.fscale(WELCH_WIN_LENGTH_SAMPLES, 1 / sr.fs, one_sided=True),
189 'ap_power': psds.T / len(t0s), # shape: (nfreqs, nchannels)
190 }
191 for k in files:
192 np.save(files[k], results[k])
193 qc_files.extend([files[k] for k in files])
194 for p in [10, 90]:
195 self.metrics[f'apRms_p{p}_raw'] = np.format_float_scientific(
196 np.percentile(results['rms'][0, :], p), precision=2)
197 self.metrics[f'apRms_p{p}_proc'] = np.format_float_scientific(
198 np.percentile(results['rms'][1, :], p), precision=2)
199 if update:
200 self.update_extended_qc(self.metrics)
201 # If lf meta and bin file present, run the old qc on LF data
202 if self.data.lf_meta and self.data.lf:
203 qc_files.extend(extract_rmsmap(self.data.lf, out_folder=self.out_path, overwrite=overwrite))
205 return qc_files
208def rmsmap(sglx, spectra=True, nmod=1):
209 """
210 Computes RMS map in time domain and spectra for each channel of Neuropixel probe
212 :param sglx: Open spikeglx reader
213 :param spectra: Whether to compute the spectra
214 :param nmod: take every nmod windows, in cases where we don't want to compute over the whole signal
215 :return: a dictionary with amplitudes in channeltime space, channelfrequency space, time
216 and frequency scales
217 """
218 rms_win_length_samples = 2 ** np.ceil(np.log2(sglx.fs * RMS_WIN_LENGTH_SECS))
219 # the window generator will generates window indices
220 wingen = utils.WindowGenerator(ns=sglx.ns, nswin=rms_win_length_samples, overlap=0)
221 nwin = np.ceil(wingen.nwin / nmod).astype(int)
222 # pre-allocate output dictionary of numpy arrays
223 win = {'TRMS': np.zeros((nwin, sglx.nc)),
224 'nsamples': np.zeros((nwin,)),
225 'fscale': fourier.fscale(WELCH_WIN_LENGTH_SAMPLES, 1 / sglx.fs, one_sided=True),
226 'tscale': wingen.tscale(fs=sglx.fs)[::nmod]}
227 win['spectral_density'] = np.zeros((len(win['fscale']), sglx.nc))
228 # loop through the whole session
229 with tqdm(total=wingen.nwin) as pbar:
230 for iwindow, (first, last) in enumerate(wingen.firstlast):
231 if np.mod(iwindow, nmod) != 0:
232 continue
234 D = sglx.read_samples(first_sample=first, last_sample=last)[0].transpose()
235 # remove low frequency noise below 1 Hz
236 D = fourier.hp(D, 1 / sglx.fs, [0, 1])
237 iw = np.floor(wingen.iw / nmod).astype(int)
238 win['TRMS'][iw, :] = utils.rms(D)
239 win['nsamples'][iw] = D.shape[1]
240 if spectra:
241 # the last window may be smaller than what is needed for welch
242 if last - first < WELCH_WIN_LENGTH_SAMPLES:
243 continue
244 # compute a smoothed spectrum using welch method
245 _, w = signal.welch(
246 D, fs=sglx.fs, window='hann', nperseg=WELCH_WIN_LENGTH_SAMPLES,
247 detrend='constant', return_onesided=True, scaling='density', axis=-1
248 )
249 win['spectral_density'] += w.T
250 # print at least every 20 windows
251 if (iw % min(20, max(int(np.floor(wingen.nwin / 75)), 1))) == 0:
252 pbar.update(iw)
253 sglx.close()
254 return win
257def extract_rmsmap(sglx, out_folder=None, overwrite=False, spectra=True, nmod=1):
258 """
259 Wrapper for rmsmap that outputs _ibl_ephysRmsMap and _ibl_ephysSpectra ALF files
261 :param sglx: Open spikeglx Reader with data for which to compute rmsmap
262 :param out_folder: folder in which to store output ALF files. Default uses the folder in which
263 the `fbin` file lives.
264 :param overwrite: do not re-extract if all ALF files already exist
265 :param spectra: Whether to compute the spectral density across the signal
266 :param nmod: take every nmod windows, in cases where we don't want to compute over the whole signal
267 :return: None
268 """
269 if out_folder is None:
270 out_folder = sglx.file_bin.parent
271 else:
272 out_folder = Path(out_folder)
273 _logger.info(f"Computing RMS map for .{sglx.type} data in {out_folder}")
274 alf_object_time = f'ephysTimeRms{sglx.type.upper()}'
275 alf_object_freq = f'ephysSpectralDensity{sglx.type.upper()}'
276 files_time = list(out_folder.glob(f"_iblqc_{alf_object_time}*"))
277 files_freq = list(out_folder.glob(f"_iblqc_{alf_object_freq}*"))
278 if (len(files_time) == 2 == len(files_freq)) and not overwrite:
279 _logger.warning(f'RMS map already exists for .{sglx.type} data in {out_folder}, skipping. Use overwrite option.')
280 return files_time + files_freq
281 # crunch numbers
282 rms = rmsmap(sglx, spectra=spectra, nmod=nmod)
283 # output ALF files, single precision with the optional label as suffix before extension
284 if not out_folder.exists():
285 out_folder.mkdir()
286 tdict = {'rms': rms['TRMS'].astype(np.single), 'timestamps': rms['tscale'].astype(np.single)}
287 out_time = alfio.save_object_npy(
288 out_folder, object=alf_object_time, dico=tdict, namespace='iblqc')
289 if spectra:
290 fdict = {'power': rms['spectral_density'].astype(np.single),
291 'freqs': rms['fscale'].astype(np.single)}
292 out_freq = alfio.save_object_npy(
293 out_folder, object=alf_object_freq, dico=fdict, namespace='iblqc')
294 return out_time + out_freq if spectra else out_time
297def raw_qc_session(session_path, overwrite=False):
298 """
299 Wrapper that exectutes QC from a session folder and outputs the results whithin the same folder
300 as the original raw data.
301 :param session_path: path of the session (Subject/yyyy-mm-dd/number
302 :param overwrite: bool (False) Force means overwriting an existing QC file
303 :return: None
304 """
305 efiles = spikeglx.glob_ephys_files(session_path)
306 qc_files = []
307 for efile in efiles:
308 if efile.get('ap') and efile.ap.exists():
309 qc_files.extend(extract_rmsmap(efile.ap, out_folder=None, overwrite=overwrite))
310 if efile.get('lf') and efile.lf.exists():
311 qc_files.extend(extract_rmsmap(efile.lf, out_folder=None, overwrite=overwrite))
312 return qc_files
315def validate_ttl_test(ses_path, display=False):
316 """
317 For a mock session on the Ephys Choice world task, check the sync channels for all
318 device properly connected and perform a synchronization if dual probes to check that
319 all channels are recorded properly
320 :param ses_path: session path
321 :param display: show the probe synchronization plot if several probes
322 :return: True if tests pass, errors otherwise
323 """
325 def _single_test(assertion, str_ok, str_ko): 1idefbc
326 if assertion: 1defbc
327 _logger.info(str_ok) 1defbc
328 return True 1defbc
329 else:
330 _logger.error(str_ko)
331 return False
333 EXPECTED_RATES_HZ = {'left_camera': 60, 'right_camera': 150, 'body_camera': 30} 1idefbc
334 SYNC_RATE_HZ = 1 1idefbc
335 MIN_TRIALS_NB = 6 1idefbc
337 ok = True 1idefbc
338 ses_path = Path(ses_path) 1idefbc
339 if not ses_path.exists(): 1idefbc
340 return False
342 # get the synchronization fronts (from the raw binary if necessary)
343 ephys_fpga.extract_sync(session_path=ses_path, overwrite=False) 1idefbc
344 rawsync, sync_map = ephys_fpga.get_main_probe_sync(ses_path) 1idefbc
345 last_time = rawsync['times'][-1] 1defbc
347 # get upgoing fronts for each
348 sync = Bunch({}) 1defbc
349 for k in sync_map: 1defbc
350 fronts = ephys_fpga.get_sync_fronts(rawsync, sync_map[k]) 1defbc
351 sync[k] = fronts['times'][fronts['polarities'] == 1] 1defbc
352 wheel = ephys_fpga.extract_wheel_sync(rawsync, chmap=sync_map) 1defbc
354 frame_rates = {'right_camera': np.round(1 / np.median(np.diff(sync.right_camera))), 1defbc
355 'left_camera': np.round(1 / np.median(np.diff(sync.left_camera))),
356 'body_camera': np.round(1 / np.median(np.diff(sync.body_camera)))}
358 # check the camera frame rates
359 for lab in frame_rates: 1defbc
360 expect = EXPECTED_RATES_HZ[lab] 1defbc
361 ok &= _single_test(assertion=abs((1 - frame_rates[lab] / expect)) < 0.1, 1defbc
362 str_ok=f'PASS: {lab} frame rate: {frame_rates[lab]} = {expect} Hz',
363 str_ko=f'FAILED: {lab} frame rate: {frame_rates[lab]} != {expect} Hz')
365 # check that the wheel has a minimum rate of activity on both channels
366 re_test = abs(1 - sync.rotary_encoder_1.size / sync.rotary_encoder_0.size) < 0.1 1defbc
367 re_test &= len(wheel[1]) / last_time > 5 1defbc
368 ok &= _single_test(assertion=re_test, 1defbc
369 str_ok="PASS: Rotary encoder", str_ko="FAILED: Rotary encoder")
370 # check that the frame 2 ttls has a minimum rate of activity
371 ok &= _single_test(assertion=len(sync.frame2ttl) / last_time > 0.2, 1defbc
372 str_ok="PASS: Frame2TTL", str_ko="FAILED: Frame2TTL")
373 # the audio has to have at least one event per trial
374 ok &= _single_test(assertion=len(sync.bpod) > len(sync.audio) > MIN_TRIALS_NB, 1defbc
375 str_ok="PASS: audio", str_ko="FAILED: audio")
376 # the bpod has to have at least twice the amount of min trial pulses
377 ok &= _single_test(assertion=len(sync.bpod) > MIN_TRIALS_NB * 2, 1defbc
378 str_ok="PASS: Bpod", str_ko="FAILED: Bpod")
379 try: 1defbc
380 # note: tried to depend as little as possible on the extraction code but for the valve...
381 extractor = ephys_fpga.FpgaTrials(ses_path) 1defbc
382 _, bpod_intervals = extractor.get_bpod_event_times(rawsync, sync_map) 1defbc
383 t_valve_open = bpod_intervals['valve_open'][:, 0] 1defbc
384 res = t_valve_open.size > 1 1defbc
385 except AssertionError:
386 res = False
387 # check that the reward valve is actionned at least once
388 ok &= _single_test(assertion=res, 1defbc
389 str_ok="PASS: Valve open", str_ko="FAILED: Valve open not detected")
390 _logger.info('ALL CHECKS PASSED !') 1defbc
392 # the imec sync is for 3B Probes only
393 if sync.get('imec_sync') is not None: 1defbc
394 ok &= _single_test(assertion=np.all(1 - SYNC_RATE_HZ * np.diff(sync.imec_sync) < 0.1), 1bc
395 str_ok="PASS: imec sync", str_ko="FAILED: imec sync")
397 # second step is to test that we can make the sync. Assertions are whithin the synch code
398 if sync.get('imec_sync') is not None: 1defbc
399 sync_result, _ = sync_probes.version3B(ses_path, display=display) 1bc
400 else:
401 sync_result, _ = sync_probes.version3A(ses_path, display=display) 1def
403 ok &= _single_test(assertion=sync_result, str_ok="PASS: synchronisation", 1defbc
404 str_ko="FAILED: probe synchronizations threshold exceeded")
406 if not ok: 1defbc
407 raise ValueError('FAILED TTL test')
408 return ok 1defbc
411def spike_sorting_metrics_ks2(ks2_path=None, m=None, save=True, save_path=None):
412 """
413 Given a path containing kilosort 2 output, compute quality metrics and optionally save them
414 to a clusters_metric.csv file
415 :param ks2_path:
416 :param save
417 :param save_path: If not given will save into the path given as ks2_path
418 :return:
419 """
421 save_path = save_path or ks2_path
423 # ensure that either a ks2_path or a phylib `TemplateModel` object with unit info is given
424 assert not (ks2_path is None and m is None), 'Must either specify a path to a ks2 output ' \
425 'directory, or a phylib `TemplateModel` object'
426 # create phylib `TemplateModel` if not given
427 m = phy_model_from_ks2_path(ks2_path) if None else m
428 c, drift = spike_sorting_metrics(m.spike_times, m.spike_clusters, m.amplitudes, m.depths,
429 cluster_ids=np.arange(m.clusters_channels.size))
430 # include the ks2 cluster contamination if `cluster_ContamPct` file exists
431 file_contamination = ks2_path.joinpath('cluster_ContamPct.tsv')
432 if file_contamination.exists():
433 contam = pd.read_csv(file_contamination, sep='\t')
434 contam.rename(columns={'ContamPct': 'ks2_contamination_pct'}, inplace=True)
435 c = c.set_index('cluster_id', drop=False).join(contam.set_index('cluster_id'))
437 # include the ks2 cluster labels if `cluster_KSLabel` file exists
438 file_labels = ks2_path.joinpath('cluster_KSLabel.tsv')
439 if file_labels.exists():
440 ks2_labels = pd.read_csv(file_labels, sep='\t')
441 ks2_labels.rename(columns={'KSLabel': 'ks2_label'}, inplace=True)
442 c = c.set_index('cluster_id', drop=False).join(ks2_labels.set_index('cluster_id'))
444 if save:
445 Path(save_path).mkdir(exist_ok=True, parents=True)
446 # the file name contains the label of the probe (directory name in this case)
447 c.to_csv(Path(save_path).joinpath('cluster_metrics.csv'))
449 return c
452def phy_model_from_ks2_path(ks2_path, bin_path, bin_file=None):
453 if not bin_file:
454 bin_file = next(bin_path.rglob('*.ap.*bin'), None)
455 meta_file = next(bin_path.rglob('*.ap.meta'), None)
456 if meta_file and meta_file.exists():
457 meta = spikeglx.read_meta_data(meta_file)
458 fs = spikeglx._get_fs_from_meta(meta)
459 nch = (spikeglx._get_nchannels_from_meta(meta) -
460 len(spikeglx._get_sync_trace_indices_from_meta(meta)))
461 else:
462 fs = 30000
463 nch = 384
464 m = model.TemplateModel(dir_path=ks2_path,
465 dat_path=bin_file, # this assumes the raw data is in the same folder
466 sample_rate=fs,
467 n_channels_dat=nch,
468 n_closest_channels=NCH_WAVEFORMS)
469 m.depths = m.get_depths()
470 return m
473# Make a bunch gathering all trial QC
474def qc_fpga_task(fpga_trials, alf_trials):
475 """
476 :fpga_task is the dictionary output of
477 ibllib.io.extractors.ephys_fpga.extract_behaviour_sync
478 : bpod_trials is the dictionary output of ibllib.io.extractors.ephys_trials.extract_all
479 : alf_trials is the ALF _ibl_trials object after extraction (alfio.load_object)
480 :return: qc_session, qc_trials, True means QC passes while False indicates a failure
481 """
483 GOCUE_STIMON_DELAY = 0.01 # -> 0.1 1g
484 FEEDBACK_STIMFREEZE_DELAY = 0.01 # -> 0.1 1g
485 VALVE_STIM_OFF_DELAY = 1 1g
486 VALVE_STIM_OFF_JITTER = 0.1 1g
487 ITI_IN_STIM_OFF_JITTER = 0.1 1g
488 ERROR_STIM_OFF_DELAY = 2 1g
489 ERROR_STIM_OFF_JITTER = 0.1 1g
490 RESPONSE_FEEDBACK_DELAY = 0.0005 1g
492 def strictly_after(t0, t1, threshold): 1g
493 """ returns isafter, iswithinthreshold"""
494 return (t1 - t0) > 0, np.abs((t1 - t0)) <= threshold 1g
496 ntrials = fpga_trials['stimOn_times'].size 1g
497 qc_trials = Bunch({}) 1g
499 """ 1g
500 First Check consistency of the dataset: whithin each trial, all events happen after trial
501 start should not be NaNs and increasing. This is not a QC but an assertion.
502 """
503 status = True 1g
504 for k in ['response_times', 'stimOn_times', 'response_times', 1g
505 'goCueTrigger_times', 'goCue_times', 'feedback_times']:
506 if k.endswith('_bpod'): 1g
507 tstart = alf_trials['intervals_bpod'][:, 0]
508 else:
509 tstart = alf_trials['intervals'][:, 0] 1g
510 selection = ~np.isnan(alf_trials[k]) 1g
511 status &= np.all(alf_trials[k][selection] - tstart[selection] > 0) 1g
512 status &= np.all(np.diff(alf_trials[k][selection]) > 0) 1g
513 assert status 1g
515 """ 1g
516 This part of the function uses only fpga_trials information
517 """
518 # check number of feedbacks: should always be one
519 qc_trials['n_feedback'] = (np.uint32(~np.isnan(fpga_trials['valveOpen_times'])) + 1g
520 np.uint32(~np.isnan(fpga_trials['errorCue_times'])))
522 # check for non-Nans
523 qc_trials['stimOn_times_nan'] = ~np.isnan(fpga_trials['stimOn_times']) 1g
524 qc_trials['goCue_times_nan'] = ~np.isnan(fpga_trials['goCue_times']) 1g
526 # stimOn before goCue
527 qc_trials['stimOn_times_before_goCue_times'], qc_trials['stimOn_times_goCue_times_delay'] =\ 1g
528 strictly_after(fpga_trials['stimOn_times'], fpga_trials['goCue_times'], GOCUE_STIMON_DELAY)
530 # stimFreeze before feedback
531 qc_trials['stim_freeze_before_feedback'], qc_trials['stim_freeze_feedback_delay'] = \ 1g
532 strictly_after(fpga_trials['stimFreeze_times'], fpga_trials['feedback_times'],
533 FEEDBACK_STIMFREEZE_DELAY)
535 # stimOff 1 sec after valve, with 0.1 as acceptable jitter
536 qc_trials['stimOff_delay_valve'] = np.less( 1g
537 np.abs(
538 fpga_trials['stimOff_times'] - fpga_trials['valveOpen_times'] - VALVE_STIM_OFF_DELAY
539 ),
540 VALVE_STIM_OFF_JITTER, out=np.ones(ntrials, dtype=bool),
541 where=~np.isnan(fpga_trials['valveOpen_times']))
543 # iti_in whithin 0.01 sec of stimOff
544 qc_trials['iti_in_delay_stim_off'] = \ 1g
545 np.abs(fpga_trials['stimOff_times'] - fpga_trials['itiIn_times']) < ITI_IN_STIM_OFF_JITTER
547 # stimOff 2 secs after errorCue_times with jitter
548 # noise off happens 2 secs after stimm, with 0.1 as acceptable jitter
549 qc_trials['stimOff_delay_noise'] = np.less( 1g
550 np.abs(
551 fpga_trials['stimOff_times'] - fpga_trials['errorCue_times'] - ERROR_STIM_OFF_DELAY
552 ),
553 ERROR_STIM_OFF_JITTER, out=np.ones(ntrials, dtype=bool),
554 where=~np.isnan(fpga_trials['errorCue_times']))
556 """ 1g
557 This part uses only alf_trials information
558 """
559 # TEST Response times (from session start) should be increasing continuously
560 # Note: RT are not durations but time stamps from session start
561 # 1. check for non-Nans
562 qc_trials['response_times_nan'] = ~np.isnan(alf_trials['response_times']) 1g
563 # 2. check for positive increase
564 qc_trials['response_times_increase'] = \ 1g
565 np.diff(np.append([0], alf_trials['response_times'])) > 0
566 # TEST Response times (from goCue) should be positive
567 qc_trials['response_times_goCue_times_diff'] = \ 1g
568 alf_trials['response_times'] - alf_trials['goCue_times'] > 0
569 # TEST 1. Response_times should be before feedback
570 qc_trials['response_before_feedback'] = \ 1g
571 alf_trials['feedback_times'] - alf_trials['response_times'] > 0
572 # 2. Delay between wheel reaches threshold (response time) and
573 # feedback is 100us, acceptable jitter 500 us
574 qc_trials['response_feedback_delay'] = \ 1g
575 alf_trials['feedback_times'] - alf_trials['response_times'] < RESPONSE_FEEDBACK_DELAY
577 # Test output at session level
578 qc_session = {k: np.all(qc_trials[k]) for k in qc_trials} 1g
580 return qc_session, qc_trials 1g