Merge branch 'release/2.7.1'

Author: juhuntenburg

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,26 +445,37 @@ 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.]')
     plt.xlabel('time [sec]')
     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]')

ibllib/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "2.7.0"
+__version__ = "2.7.1"
 import warnings
 
 from ibllib.misc import logger_config

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

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,14 +401,16 @@ 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')
         else:
             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')

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