Coverage for ibllib/ephys/ephysqc.py: 49%

287 statements  

« prev     ^ index     » next       coverage.py v7.5.4, created at 2024-07-08 17:16 +0100

1""" 

2Quality control of raw Neuropixel electrophysiology data. 

3""" 

4from pathlib import Path 

5import logging 

6 

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 

16 

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 

23 

24 

25_logger = logging.getLogger(__name__) 

26 

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 

34 

35 

36class EphysQC(base.QC): 

37 """ 

38 A class for computing Ephys QC metrics. 

39 

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 """ 

44 

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) 

48 

49 if self.use_alyx: 

50 super().__init__(probe_id, endpoint='insertions', **kwargs) 

51 self._outcome = 'NOT_SET' 

52 self.pid = probe_id 

53 

54 self.session_path = session_path 

55 keys = ('ap', 'ap_meta', 'lf', 'lf_meta') 

56 self.data = Bunch.fromkeys(keys) 

57 self.metrics = {} 

58 

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

76 

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() 

84 

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) 

103 

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 

123 

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. 

127 

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 

137 

138 # Load data 

139 self.load_data() 

140 self.out_path = kwargs.get('out_path', self.probe_path) 

141 

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 

159 

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) 

167 

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") 

171 

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)) 

176 

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)) 

204 

205 return qc_files 

206 

207 

208def rmsmap(sglx): 

209 """ 

210 Computes RMS map in time domain and spectra for each channel of Neuropixel probe 

211 

212 :param sglx: Open spikeglx reader 

213 :return: a dictionary with amplitudes in channeltime space, channelfrequency space, time 

214 and frequency scales 

215 """ 

216 rms_win_length_samples = 2 ** np.ceil(np.log2(sglx.fs * RMS_WIN_LENGTH_SECS)) 

217 # the window generator will generates window indices 

218 wingen = utils.WindowGenerator(ns=sglx.ns, nswin=rms_win_length_samples, overlap=0) 

219 # pre-allocate output dictionary of numpy arrays 

220 win = {'TRMS': np.zeros((wingen.nwin, sglx.nc)), 

221 'nsamples': np.zeros((wingen.nwin,)), 

222 'fscale': fourier.fscale(WELCH_WIN_LENGTH_SAMPLES, 1 / sglx.fs, one_sided=True), 

223 'tscale': wingen.tscale(fs=sglx.fs)} 

224 win['spectral_density'] = np.zeros((len(win['fscale']), sglx.nc)) 

225 # loop through the whole session 

226 with tqdm(total=wingen.nwin) as pbar: 

227 for first, last in wingen.firstlast: 

228 D = sglx.read_samples(first_sample=first, last_sample=last)[0].transpose() 

229 # remove low frequency noise below 1 Hz 

230 D = fourier.hp(D, 1 / sglx.fs, [0, 1]) 

231 iw = wingen.iw 

232 win['TRMS'][iw, :] = utils.rms(D) 

233 win['nsamples'][iw] = D.shape[1] 

234 # the last window may be smaller than what is needed for welch 

235 if last - first < WELCH_WIN_LENGTH_SAMPLES: 

236 continue 

237 # compute a smoothed spectrum using welch method 

238 _, w = signal.welch( 

239 D, fs=sglx.fs, window='hann', nperseg=WELCH_WIN_LENGTH_SAMPLES, 

240 detrend='constant', return_onesided=True, scaling='density', axis=-1 

241 ) 

242 win['spectral_density'] += w.T 

243 # print at least every 20 windows 

244 if (iw % min(20, max(int(np.floor(wingen.nwin / 75)), 1))) == 0: 

245 pbar.update(iw) 

246 sglx.close() 

247 return win 

248 

249 

250def extract_rmsmap(sglx, out_folder=None, overwrite=False): 

251 """ 

252 Wrapper for rmsmap that outputs _ibl_ephysRmsMap and _ibl_ephysSpectra ALF files 

253 

254 :param sglx: Open spikeglx Reader with data for which to compute rmsmap 

255 :param out_folder: folder in which to store output ALF files. Default uses the folder in which 

256 the `fbin` file lives. 

257 :param overwrite: do not re-extract if all ALF files already exist 

258 :param label: string or list of strings that will be appended to the filename before extension 

259 :return: None 

260 """ 

261 if out_folder is None: 

262 out_folder = sglx.file_bin.parent 

263 else: 

264 out_folder = Path(out_folder) 

265 _logger.info(f"Computing RMS map for .{sglx.type} data in {out_folder}") 

266 alf_object_time = f'ephysTimeRms{sglx.type.upper()}' 

267 alf_object_freq = f'ephysSpectralDensity{sglx.type.upper()}' 

268 files_time = list(out_folder.glob(f"_iblqc_{alf_object_time}*")) 

269 files_freq = list(out_folder.glob(f"_iblqc_{alf_object_freq}*")) 

270 if (len(files_time) == 2 == len(files_freq)) and not overwrite: 

271 _logger.warning(f'RMS map already exists for .{sglx.type} data in {out_folder}, skipping. Use overwrite option.') 

272 return files_time + files_freq 

273 # crunch numbers 

274 rms = rmsmap(sglx) 

275 # output ALF files, single precision with the optional label as suffix before extension 

276 if not out_folder.exists(): 

277 out_folder.mkdir() 

278 tdict = {'rms': rms['TRMS'].astype(np.single), 'timestamps': rms['tscale'].astype(np.single)} 

279 fdict = {'power': rms['spectral_density'].astype(np.single), 

280 'freqs': rms['fscale'].astype(np.single)} 

281 out_time = alfio.save_object_npy( 

282 out_folder, object=alf_object_time, dico=tdict, namespace='iblqc') 

283 out_freq = alfio.save_object_npy( 

284 out_folder, object=alf_object_freq, dico=fdict, namespace='iblqc') 

285 return out_time + out_freq 

286 

287 

288def raw_qc_session(session_path, overwrite=False): 

289 """ 

290 Wrapper that exectutes QC from a session folder and outputs the results whithin the same folder 

291 as the original raw data. 

292 :param session_path: path of the session (Subject/yyyy-mm-dd/number 

293 :param overwrite: bool (False) Force means overwriting an existing QC file 

294 :return: None 

295 """ 

296 efiles = spikeglx.glob_ephys_files(session_path) 

297 qc_files = [] 

298 for efile in efiles: 

299 if efile.get('ap') and efile.ap.exists(): 

300 qc_files.extend(extract_rmsmap(efile.ap, out_folder=None, overwrite=overwrite)) 

301 if efile.get('lf') and efile.lf.exists(): 

302 qc_files.extend(extract_rmsmap(efile.lf, out_folder=None, overwrite=overwrite)) 

303 return qc_files 

304 

305 

306def validate_ttl_test(ses_path, display=False): 

307 """ 

308 For a mock session on the Ephys Choice world task, check the sync channels for all 

309 device properly connected and perform a synchronization if dual probes to check that 

310 all channels are recorded properly 

311 :param ses_path: session path 

312 :param display: show the probe synchronization plot if several probes 

313 :return: True if tests pass, errors otherwise 

314 """ 

315 

316 def _single_test(assertion, str_ok, str_ko): 1idefbc

317 if assertion: 1defbc

318 _logger.info(str_ok) 1defbc

319 return True 1defbc

320 else: 

321 _logger.error(str_ko) 

322 return False 

323 

324 EXPECTED_RATES_HZ = {'left_camera': 60, 'right_camera': 150, 'body_camera': 30} 1idefbc

325 SYNC_RATE_HZ = 1 1idefbc

326 MIN_TRIALS_NB = 6 1idefbc

327 

328 ok = True 1idefbc

329 ses_path = Path(ses_path) 1idefbc

330 if not ses_path.exists(): 1idefbc

331 return False 

332 

333 # get the synchronization fronts (from the raw binary if necessary) 

334 ephys_fpga.extract_sync(session_path=ses_path, overwrite=False) 1idefbc

335 rawsync, sync_map = ephys_fpga.get_main_probe_sync(ses_path) 1idefbc

336 last_time = rawsync['times'][-1] 1defbc

337 

338 # get upgoing fronts for each 

339 sync = Bunch({}) 1defbc

340 for k in sync_map: 1defbc

341 fronts = ephys_fpga.get_sync_fronts(rawsync, sync_map[k]) 1defbc

342 sync[k] = fronts['times'][fronts['polarities'] == 1] 1defbc

343 wheel = ephys_fpga.extract_wheel_sync(rawsync, chmap=sync_map) 1defbc

344 

345 frame_rates = {'right_camera': np.round(1 / np.median(np.diff(sync.right_camera))), 1defbc

346 'left_camera': np.round(1 / np.median(np.diff(sync.left_camera))), 

347 'body_camera': np.round(1 / np.median(np.diff(sync.body_camera)))} 

348 

349 # check the camera frame rates 

350 for lab in frame_rates: 1defbc

351 expect = EXPECTED_RATES_HZ[lab] 1defbc

352 ok &= _single_test(assertion=abs((1 - frame_rates[lab] / expect)) < 0.1, 1defbc

353 str_ok=f'PASS: {lab} frame rate: {frame_rates[lab]} = {expect} Hz', 

354 str_ko=f'FAILED: {lab} frame rate: {frame_rates[lab]} != {expect} Hz') 

355 

356 # check that the wheel has a minimum rate of activity on both channels 

357 re_test = abs(1 - sync.rotary_encoder_1.size / sync.rotary_encoder_0.size) < 0.1 1defbc

358 re_test &= len(wheel[1]) / last_time > 5 1defbc

359 ok &= _single_test(assertion=re_test, 1defbc

360 str_ok="PASS: Rotary encoder", str_ko="FAILED: Rotary encoder") 

361 # check that the frame 2 ttls has a minimum rate of activity 

362 ok &= _single_test(assertion=len(sync.frame2ttl) / last_time > 0.2, 1defbc

363 str_ok="PASS: Frame2TTL", str_ko="FAILED: Frame2TTL") 

364 # the audio has to have at least one event per trial 

365 ok &= _single_test(assertion=len(sync.bpod) > len(sync.audio) > MIN_TRIALS_NB, 1defbc

366 str_ok="PASS: audio", str_ko="FAILED: audio") 

367 # the bpod has to have at least twice the amount of min trial pulses 

368 ok &= _single_test(assertion=len(sync.bpod) > MIN_TRIALS_NB * 2, 1defbc

369 str_ok="PASS: Bpod", str_ko="FAILED: Bpod") 

370 try: 1defbc

371 # note: tried to depend as little as possible on the extraction code but for the valve... 

372 extractor = ephys_fpga.FpgaTrials(ses_path) 1defbc

373 _, bpod_intervals = extractor.get_bpod_event_times(rawsync, sync_map) 1defbc

374 t_valve_open = bpod_intervals['valve_open'][:, 0] 1defbc

375 res = t_valve_open.size > 1 1defbc

376 except AssertionError: 

377 res = False 

378 # check that the reward valve is actionned at least once 

379 ok &= _single_test(assertion=res, 1defbc

380 str_ok="PASS: Valve open", str_ko="FAILED: Valve open not detected") 

381 _logger.info('ALL CHECKS PASSED !') 1defbc

382 

383 # the imec sync is for 3B Probes only 

384 if sync.get('imec_sync') is not None: 1defbc

385 ok &= _single_test(assertion=np.all(1 - SYNC_RATE_HZ * np.diff(sync.imec_sync) < 0.1), 1bc

386 str_ok="PASS: imec sync", str_ko="FAILED: imec sync") 

387 

388 # second step is to test that we can make the sync. Assertions are whithin the synch code 

389 if sync.get('imec_sync') is not None: 1defbc

390 sync_result, _ = sync_probes.version3B(ses_path, display=display) 1bc

391 else: 

392 sync_result, _ = sync_probes.version3A(ses_path, display=display) 1def

393 

394 ok &= _single_test(assertion=sync_result, str_ok="PASS: synchronisation", 1defbc

395 str_ko="FAILED: probe synchronizations threshold exceeded") 

396 

397 if not ok: 1defbc

398 raise ValueError('FAILED TTL test') 

399 return ok 1defbc

400 

401 

402def spike_sorting_metrics_ks2(ks2_path=None, m=None, save=True, save_path=None): 

403 """ 

404 Given a path containing kilosort 2 output, compute quality metrics and optionally save them 

405 to a clusters_metric.csv file 

406 :param ks2_path: 

407 :param save 

408 :param save_path: If not given will save into the path given as ks2_path 

409 :return: 

410 """ 

411 

412 save_path = save_path or ks2_path 

413 

414 # ensure that either a ks2_path or a phylib `TemplateModel` object with unit info is given 

415 assert not (ks2_path is None and m is None), 'Must either specify a path to a ks2 output ' \ 

416 'directory, or a phylib `TemplateModel` object' 

417 # create phylib `TemplateModel` if not given 

418 m = phy_model_from_ks2_path(ks2_path) if None else m 

419 c, drift = spike_sorting_metrics(m.spike_times, m.spike_clusters, m.amplitudes, m.depths, 

420 cluster_ids=np.arange(m.clusters_channels.size)) 

421 # include the ks2 cluster contamination if `cluster_ContamPct` file exists 

422 file_contamination = ks2_path.joinpath('cluster_ContamPct.tsv') 

423 if file_contamination.exists(): 

424 contam = pd.read_csv(file_contamination, sep='\t') 

425 contam.rename(columns={'ContamPct': 'ks2_contamination_pct'}, inplace=True) 

426 c = c.set_index('cluster_id', drop=False).join(contam.set_index('cluster_id')) 

427 

428 # include the ks2 cluster labels if `cluster_KSLabel` file exists 

429 file_labels = ks2_path.joinpath('cluster_KSLabel.tsv') 

430 if file_labels.exists(): 

431 ks2_labels = pd.read_csv(file_labels, sep='\t') 

432 ks2_labels.rename(columns={'KSLabel': 'ks2_label'}, inplace=True) 

433 c = c.set_index('cluster_id', drop=False).join(ks2_labels.set_index('cluster_id')) 

434 

435 if save: 

436 Path(save_path).mkdir(exist_ok=True, parents=True) 

437 # the file name contains the label of the probe (directory name in this case) 

438 c.to_csv(Path(save_path).joinpath('cluster_metrics.csv')) 

439 

440 return c 

441 

442 

443def phy_model_from_ks2_path(ks2_path, bin_path, bin_file=None): 

444 if not bin_file: 

445 bin_file = next(bin_path.rglob('*.ap.*bin'), None) 

446 meta_file = next(bin_path.rglob('*.ap.meta'), None) 

447 if meta_file and meta_file.exists(): 

448 meta = spikeglx.read_meta_data(meta_file) 

449 fs = spikeglx._get_fs_from_meta(meta) 

450 nch = (spikeglx._get_nchannels_from_meta(meta) - 

451 len(spikeglx._get_sync_trace_indices_from_meta(meta))) 

452 else: 

453 fs = 30000 

454 nch = 384 

455 m = model.TemplateModel(dir_path=ks2_path, 

456 dat_path=bin_file, # this assumes the raw data is in the same folder 

457 sample_rate=fs, 

458 n_channels_dat=nch, 

459 n_closest_channels=NCH_WAVEFORMS) 

460 m.depths = m.get_depths() 

461 return m 

462 

463 

464# Make a bunch gathering all trial QC 

465def qc_fpga_task(fpga_trials, alf_trials): 

466 """ 

467 :fpga_task is the dictionary output of 

468 ibllib.io.extractors.ephys_fpga.extract_behaviour_sync 

469 : bpod_trials is the dictionary output of ibllib.io.extractors.ephys_trials.extract_all 

470 : alf_trials is the ALF _ibl_trials object after extraction (alfio.load_object) 

471 :return: qc_session, qc_trials, True means QC passes while False indicates a failure 

472 """ 

473 

474 GOCUE_STIMON_DELAY = 0.01 # -> 0.1 1g

475 FEEDBACK_STIMFREEZE_DELAY = 0.01 # -> 0.1 1g

476 VALVE_STIM_OFF_DELAY = 1 1g

477 VALVE_STIM_OFF_JITTER = 0.1 1g

478 ITI_IN_STIM_OFF_JITTER = 0.1 1g

479 ERROR_STIM_OFF_DELAY = 2 1g

480 ERROR_STIM_OFF_JITTER = 0.1 1g

481 RESPONSE_FEEDBACK_DELAY = 0.0005 1g

482 

483 def strictly_after(t0, t1, threshold): 1g

484 """ returns isafter, iswithinthreshold""" 

485 return (t1 - t0) > 0, np.abs((t1 - t0)) <= threshold 1g

486 

487 ntrials = fpga_trials['stimOn_times'].size 1g

488 qc_trials = Bunch({}) 1g

489 

490 """ 1g

491 First Check consistency of the dataset: whithin each trial, all events happen after trial 

492 start should not be NaNs and increasing. This is not a QC but an assertion. 

493 """ 

494 status = True 1g

495 for k in ['response_times', 'stimOn_times', 'response_times', 1g

496 'goCueTrigger_times', 'goCue_times', 'feedback_times']: 

497 if k.endswith('_bpod'): 1g

498 tstart = alf_trials['intervals_bpod'][:, 0] 

499 else: 

500 tstart = alf_trials['intervals'][:, 0] 1g

501 selection = ~np.isnan(alf_trials[k]) 1g

502 status &= np.all(alf_trials[k][selection] - tstart[selection] > 0) 1g

503 status &= np.all(np.diff(alf_trials[k][selection]) > 0) 1g

504 assert status 1g

505 

506 """ 1g

507 This part of the function uses only fpga_trials information 

508 """ 

509 # check number of feedbacks: should always be one 

510 qc_trials['n_feedback'] = (np.uint32(~np.isnan(fpga_trials['valveOpen_times'])) + 1g

511 np.uint32(~np.isnan(fpga_trials['errorCue_times']))) 

512 

513 # check for non-Nans 

514 qc_trials['stimOn_times_nan'] = ~np.isnan(fpga_trials['stimOn_times']) 1g

515 qc_trials['goCue_times_nan'] = ~np.isnan(fpga_trials['goCue_times']) 1g

516 

517 # stimOn before goCue 

518 qc_trials['stimOn_times_before_goCue_times'], qc_trials['stimOn_times_goCue_times_delay'] =\ 1g

519 strictly_after(fpga_trials['stimOn_times'], fpga_trials['goCue_times'], GOCUE_STIMON_DELAY) 

520 

521 # stimFreeze before feedback 

522 qc_trials['stim_freeze_before_feedback'], qc_trials['stim_freeze_feedback_delay'] = \ 1g

523 strictly_after(fpga_trials['stimFreeze_times'], fpga_trials['feedback_times'], 

524 FEEDBACK_STIMFREEZE_DELAY) 

525 

526 # stimOff 1 sec after valve, with 0.1 as acceptable jitter 

527 qc_trials['stimOff_delay_valve'] = np.less( 1g

528 np.abs( 

529 fpga_trials['stimOff_times'] - fpga_trials['valveOpen_times'] - VALVE_STIM_OFF_DELAY 

530 ), 

531 VALVE_STIM_OFF_JITTER, out=np.ones(ntrials, dtype=bool), 

532 where=~np.isnan(fpga_trials['valveOpen_times'])) 

533 

534 # iti_in whithin 0.01 sec of stimOff 

535 qc_trials['iti_in_delay_stim_off'] = \ 1g

536 np.abs(fpga_trials['stimOff_times'] - fpga_trials['itiIn_times']) < ITI_IN_STIM_OFF_JITTER 

537 

538 # stimOff 2 secs after errorCue_times with jitter 

539 # noise off happens 2 secs after stimm, with 0.1 as acceptable jitter 

540 qc_trials['stimOff_delay_noise'] = np.less( 1g

541 np.abs( 

542 fpga_trials['stimOff_times'] - fpga_trials['errorCue_times'] - ERROR_STIM_OFF_DELAY 

543 ), 

544 ERROR_STIM_OFF_JITTER, out=np.ones(ntrials, dtype=bool), 

545 where=~np.isnan(fpga_trials['errorCue_times'])) 

546 

547 """ 1g

548 This part uses only alf_trials information 

549 """ 

550 # TEST Response times (from session start) should be increasing continuously 

551 # Note: RT are not durations but time stamps from session start 

552 # 1. check for non-Nans 

553 qc_trials['response_times_nan'] = ~np.isnan(alf_trials['response_times']) 1g

554 # 2. check for positive increase 

555 qc_trials['response_times_increase'] = \ 1g

556 np.diff(np.append([0], alf_trials['response_times'])) > 0 

557 # TEST Response times (from goCue) should be positive 

558 qc_trials['response_times_goCue_times_diff'] = \ 1g

559 alf_trials['response_times'] - alf_trials['goCue_times'] > 0 

560 # TEST 1. Response_times should be before feedback 

561 qc_trials['response_before_feedback'] = \ 1g

562 alf_trials['feedback_times'] - alf_trials['response_times'] > 0 

563 # 2. Delay between wheel reaches threshold (response time) and 

564 # feedback is 100us, acceptable jitter 500 us 

565 qc_trials['response_feedback_delay'] = \ 1g

566 alf_trials['feedback_times'] - alf_trials['response_times'] < RESPONSE_FEEDBACK_DELAY 

567 

568 # Test output at session level 

569 qc_session = {k: np.all(qc_trials[k]) for k in qc_trials} 1g

570 

571 return qc_session, qc_trials 1g