diff --git a/brainbox/behavior/dlc.py b/brainbox/behavior/dlc.py index c5451d561..097c78a30 100644 --- a/brainbox/behavior/dlc.py +++ b/brainbox/behavior/dlc.py @@ -41,6 +41,8 @@ def insert_idx(array, values): idx[np.where(abs(values - array[idx - 1]) < abs(values - array[idx]))] -= 1 # If 0 index was reduced, revert idx[idx == -1] = 0 + if np.all(idx == 0): + raise ValueError('Something is wrong, all values to insert are outside of the array.') return idx @@ -121,7 +123,7 @@ def get_feature_event_times(dlc, dlc_t, features): def get_licks(dlc, dlc_t): """ - Compute lick times from the toungue dlc points + Compute lick times from the tongue dlc points :param dlc: dlc pqt table :param dlc_t: dlc times :return: @@ -216,6 +218,9 @@ def get_smooth_pupil_diameter(diameter_raw, camera, std_thresh=5, nan_thresh=1): else: raise NotImplementedError("camera has to be 'left' or 'right") + # Raise error if too many NaN time points, in this case it doesn't make sense to interpolate + if np.mean(np.isnan(diameter_raw)) > 0.9: + raise ValueError(f"Raw pupil diameter for {camera} is too often NaN, cannot smooth.") # run savitzy-golay filter on non-nan time points to denoise diameter_smoothed = smooth_interpolate_savgol(diameter_raw, window=window, order=3, interp_kind='linear') @@ -440,19 +445,26 @@ def plot_motion_energy_hist(camera_dict, trials_df): 'body': '#035382'} start_window, end_window = plt_window(trials_df['stimOn_times']) + missing_data = [] for cam in camera_dict.keys(): - try: - motion_energy = zscore(camera_dict[cam]['motion_energy'], nan_policy='omit') - start_idx = insert_idx(camera_dict[cam]['times'], start_window) - end_idx = np.array(start_idx + int(WINDOW_LEN * SAMPLING[cam]), dtype='int64') - me_all = [motion_energy[start_idx[i]:end_idx[i]] for i in range(len(start_idx))] - times = np.arange(len(me_all[0])) / SAMPLING[cam] + WINDOW_LAG - me_mean = np.mean(me_all, axis=0) - me_std = np.std(me_all, axis=0) / np.sqrt(len(me_all)) - plt.plot(times, me_mean, label=f'{cam} cam', color=colors[cam], linewidth=2) - plt.fill_between(times, me_mean + me_std, me_mean - me_std, color=colors[cam], alpha=0.2) - except AttributeError: - logger.warning(f"Cannot load motion energy AND times data for {cam} camera") + if (camera_dict[cam]['motion_energy'] is not None and len(camera_dict[cam]['motion_energy']) > 0 + and camera_dict[cam]['times'] is not None and len(camera_dict[cam]['times']) > 0): + try: + motion_energy = zscore(camera_dict[cam]['motion_energy'], nan_policy='omit') + start_idx = insert_idx(camera_dict[cam]['times'], start_window) + end_idx = np.array(start_idx + int(WINDOW_LEN * SAMPLING[cam]), dtype='int64') + me_all = [motion_energy[start_idx[i]:end_idx[i]] for i in range(len(start_idx))] + times = np.arange(len(me_all[0])) / SAMPLING[cam] + WINDOW_LAG + me_mean = np.mean(me_all, axis=0) + me_std = np.std(me_all, axis=0) / np.sqrt(len(me_all)) + plt.plot(times, me_mean, label=f'{cam} cam', color=colors[cam], linewidth=2) + plt.fill_between(times, me_mean + me_std, me_mean - me_std, color=colors[cam], alpha=0.2) + except AttributeError: + logger.warning(f"Cannot load motion energy and/or times data for {cam} camera") + missing_data.append(cam) + else: + logger.warning(f"Data missing or empty for motion energy and/or times data for {cam} camera") + missing_data.append(cam) plt.xticks([-0.5, 0, 0.5, 1, 1.5]) plt.ylabel('z-scored motion energy [a.u.]') @@ -460,6 +472,10 @@ def plot_motion_energy_hist(camera_dict, trials_df): plt.axvline(x=0, label='stimOn', linestyle='--', c='k') plt.legend(loc='lower right') plt.title('Motion Energy') + if len(missing_data) > 0: + ax = plt.gca() + ax.text(.95, .35, f"Data incomplete for\n{' and '.join(missing_data)} camera", color='r', fontsize=10, + horizontalalignment='right', verticalalignment='center', transform=ax.transAxes) return plt.gca() @@ -477,6 +493,8 @@ def plot_speed_hist(dlc_df, cam_times, trials_df, feature='paw_r', cam='left', l """ # Threshold the dlc traces dlc_df = likelihood_threshold(dlc_df) + # For pre-GPIO sessions, remove the first few timestamps to match the number of frames + cam_times = cam_times[-len(dlc_df):] # Get speeds speeds = get_speed(dlc_df, cam_times, camera=cam, feature=feature) # Windows aligned to align_to @@ -495,7 +513,7 @@ def plot_speed_hist(dlc_df, cam_times, trials_df, feature='paw_r', cam='left', l plt.plot(times, pd.DataFrame.from_dict(dict(zip(incorrect.index, incorrect.values))).mean(axis=1), c='gray', label='incorrect trial') plt.axvline(x=0, label='stimOn', linestyle='--', c='r') - plt.title(f'{feature.split("_")[0].capitalize()} speed') + plt.title(f'{feature.split("_")[0].capitalize()} speed ({cam} cam)') plt.xticks([-0.5, 0, 0.5, 1, 1.5]) plt.xlabel('time [sec]') plt.ylabel('speed [px/sec]') @@ -531,7 +549,7 @@ def plot_pupil_diameter_hist(pupil_diameter, cam_times, trials_df, cam='left'): plt.plot(times, pupil_mean, label=align_to.split("_")[0], color=color) plt.fill_between(times, pupil_mean + pupil_std, pupil_mean - pupil_std, color=color, alpha=0.5) plt.axvline(x=0, linestyle='--', c='k') - plt.title('Pupil diameter') + plt.title(f'Pupil diameter ({cam} cam)') plt.xlabel('time [sec]') plt.xticks([-0.5, 0, 0.5, 1, 1.5]) plt.ylabel('pupil diameter [px]') diff --git a/ibllib/__init__.py b/ibllib/__init__.py index 74a94b448..e9a602837 100644 --- a/ibllib/__init__.py +++ b/ibllib/__init__.py @@ -1,4 +1,4 @@ -__version__ = "2.7.0" +__version__ = "2.7.1" import warnings from ibllib.misc import logger_config diff --git a/ibllib/pipes/ephys_preprocessing.py b/ibllib/pipes/ephys_preprocessing.py index 3c1c3b661..6fe4d0762 100644 --- a/ibllib/pipes/ephys_preprocessing.py +++ b/ibllib/pipes/ephys_preprocessing.py @@ -877,6 +877,9 @@ def _run(self): class EphysPostDLC(tasks.Task): """ The post_dlc task takes dlc traces as input and computes useful quantities, as well as qc. + + For creating the full dlc_qc_plot, several other inputs are required that can be found in the docstring of + :py:func:ibllib.plots.figures.dlc_qc_plot """ io_charge = 90 level = 3 @@ -894,6 +897,16 @@ class EphysPostDLC(tasks.Task): } def _run(self, overwrite=False, run_qc=True, plot_qc=True): + """ + Run the EphysPostDLC task. Returns a list of file locations for the output files in signature. The created plot + (dlc_qc_plot.png) is not returned, but saved in session_path/snapshots and uploaded to Alyx as a note. + + :param overwrite: bool, whether to recompute existing output files (default is False). + Note that the dlc_qc_plot will be computed even if overwrite = False + :param run_qc: bool, whether to run the DLC QC (default is True) + :param plot_qc: book, whether to create the dlc_qc_plot (default is True) + + """ # Check if output files exist locally exist, output_files = self.assert_expected(self.signature['output_files'], silent=True) if exist and not overwrite: @@ -916,15 +929,13 @@ def _run(self, overwrite=False, run_qc=True, plot_qc=True): dlc = pd.read_parquet(dlc_file) dlc_thresh = likelihood_threshold(dlc, 0.9) # try to load respective camera times - try: - dlc_t = np.load(next(Path(self.session_path).joinpath('alf').glob(f'_ibl_{cam}Camera.times.*npy'))) - times = True - except StopIteration: - _logger.error(f'No camera.times found for {cam} camera. ' + dlc_t = np.load(next(Path(self.session_path).joinpath('alf').glob(f'_ibl_{cam}Camera.times.*npy'))) + times = True + if dlc_t.shape[0] == 0: + _logger.error(f'camera.times empty for {cam} camera. ' f'Computations using camera.times will be skipped') self.status = -1 times = False - # These features are only computed from left and right cam if cam in ('left', 'right'): features = pd.DataFrame() @@ -937,8 +948,13 @@ def _run(self, overwrite=False, run_qc=True, plot_qc=True): # Compute pupil diameter, raw and smoothed _logger.info(f"Computing raw pupil diameter for {cam} camera.") features['pupilDiameter_raw'] = get_pupil_diameter(dlc_thresh) - _logger.info(f"Computing smooth pupil diameter for {cam} camera.") - features['pupilDiameter_smooth'] = get_smooth_pupil_diameter(features['pupilDiameter_raw'], cam) + try: + _logger.info(f"Computing smooth pupil diameter for {cam} camera.") + features['pupilDiameter_smooth'] = get_smooth_pupil_diameter(features['pupilDiameter_raw'], + cam) + except BaseException: + _logger.error(f"Computing smooth pupil diameter for {cam} camera failed, saving all NaNs.") + features['pupilDiameter_smooth'] = np.nan # Safe to pqt features_file = Path(self.session_path).joinpath('alf', f'_ibl_{cam}Camera.features.pqt') features.to_parquet(features_file) @@ -977,6 +993,7 @@ def _run(self, overwrite=False, run_qc=True, plot_qc=True): fig_path.parent.mkdir(parents=True, exist_ok=True) fig = dlc_qc_plot(self.one.path2eid(self.session_path), one=self.one) fig.savefig(fig_path) + fig.clf() snp = ReportSnapshot(self.session_path, session_id, one=self.one) snp.outputs = [fig_path] snp.register_images(widths=['orig'], @@ -1060,5 +1077,6 @@ def __init__(self, session_path=None, **kwargs): tasks["EphysCellsQc"] = EphysCellsQc(self.session_path, parents=[tasks["SpikeSorting"]]) tasks["EphysDLC"] = EphysDLC(self.session_path, parents=[tasks["EphysVideoCompress"]]) # level 3 - tasks["EphysPostDLC"] = EphysPostDLC(self.session_path, parents=[tasks["EphysDLC"]]) + tasks["EphysPostDLC"] = EphysPostDLC(self.session_path, parents=[tasks["EphysDLC"], tasks["EphysTrials"], + tasks["EphysVideoSyncQc"]]) self.tasks = tasks diff --git a/ibllib/plots/figures.py b/ibllib/plots/figures.py index 2eeb64c13..645dc4c43 100644 --- a/ibllib/plots/figures.py +++ b/ibllib/plots/figures.py @@ -3,6 +3,7 @@ """ import logging from pathlib import Path +import traceback from string import ascii_uppercase import numpy as np @@ -341,6 +342,11 @@ def dlc_qc_plot(eid, one=None): else: logger.warning(f"Could not load _ibl_{cam}Camera.{feat} some DLC QC plots have to be skipped.") data[f'{cam}_{feat}'] = None + # Sometimes there is a file but the object is empty + if data[f'{cam}_{feat}'] is not None and len(data[f'{cam}_{feat}']) == 0: + logger.warning(f"Object loaded from _ibl_{cam}Camera.{feat} is empty, some plots have to be skipped.") + data[f'{cam}_{feat}'] = None + # Session data for alf_object in ['trials', 'wheel', 'licks']: try: @@ -355,7 +361,9 @@ def dlc_qc_plot(eid, one=None): data[f'{alf_object}'] = None # Simplify to what we actually need data['licks'] = data['licks'].times if data['licks'] else None - data['left_pupil'] = data['left_features'].pupilDiameter_smooth if data['left_features'] is not None else None + data['left_pupil'] = data['left_features'].pupilDiameter_smooth if ( + data['left_features'] is not None and not np.all(np.isnan(data['left_features'].pupilDiameter_smooth)) + ) else None data['wheel_time'] = data['wheel'].timestamps if data['wheel'] is not None else None data['wheel_position'] = data['wheel'].position if data['wheel'] is not None else None if data['trials']: @@ -393,7 +401,8 @@ def dlc_qc_plot(eid, one=None): ax = plt.subplot(2, 5, i + 1) ax.text(-0.1, 1.15, ascii_uppercase[i], transform=ax.transAxes, fontsize=16, fontweight='bold') # Check if any of the inputs is None - if any([v is None for v in panel[1].values()]): + if any([v is None for v in panel[1].values()]) or any([v.values() is None for v in panel[1].values() + if isinstance(v, dict)]): ax.text(.5, .5, f"Data incomplete\n{panel[0].__name__}", color='r', fontweight='bold', fontsize=12, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes) plt.axis('off') @@ -401,6 +410,7 @@ def dlc_qc_plot(eid, one=None): try: panel[0](**panel[1]) except BaseException: + logger.error(f'Error in {panel[0].__name__}\n' + traceback.format_exc()) ax.text(.5, .5, f'Error in \n{panel[0].__name__}', color='r', fontweight='bold', fontsize=12, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes) plt.axis('off') diff --git a/release_notes.md b/release_notes.md index 70c4c5626..bb5b2823c 100644 --- a/release_notes.md +++ b/release_notes.md @@ -1,4 +1,9 @@ ## Release Note 2.7 + +### Release Notes 2.7.1 2022-01-05 + +- Fixes and better logging for EphysPostDLC task + ### Release Notes 2.7.0 2021-12-20 - Remove atlas instantiation from import of histology module