Merge branch 'release/2.17.2'

Author: juhuntenburg

brainbox/behavior/wheel.py

@@ -4,7 +4,7 @@
 import numpy as np
 from numpy import pi
 import scipy.interpolate as interpolate
-from scipy.signal import convolve, windows
+import scipy.signal
 from scipy.linalg import hankel
 import matplotlib.pyplot as plt
 from matplotlib.collections import LineCollection
@@ -90,6 +90,21 @@ def velocity(re_ts, re_pos):
         return ifcn(re_ts)
 
 
+def velocity_filtered(pos, fs, corner_frequency=20, order=8):
+    """
+    Compute wheel velocity from uniformly sampled wheel data
+
+    :param pos: vector of uniformly sampled wheel positions
+    :param fs: scalar, sampling frequency
+    :param corner_frequency: scalar, corner frequency of low-pass filter
+    :param order: scalar, order of Butterworth filter
+    """
+    sos = scipy.signal.butter(**{'N': order, 'Wn': corner_frequency / fs * 2, 'btype': 'lowpass'}, output='sos')
+    vel = np.insert(np.diff(scipy.signal.sosfiltfilt(sos, pos)), 0, 0) * fs
+    acc = np.insert(np.diff(vel), 0, 0) * fs
+    return vel, acc
+
+
 def velocity_smoothed(pos, freq, smooth_size=0.03):
     """
     Compute wheel velocity from uniformly sampled wheel data
@@ -114,12 +129,12 @@ def velocity_smoothed(pos, freq, smooth_size=0.03):
     std_samps = np.round(smooth_size * freq)  # Standard deviation relative to sampling frequency
     N = std_samps * 6  # Number of points in the Gaussian covering +/-3 standard deviations
     gauss_std = (N - 1) / 6
-    win = windows.gaussian(N, gauss_std)
+    win = scipy.signal.windows.gaussian(N, gauss_std)
     win = win / win.sum()  # Normalize amplitude
 
     # Convolve and multiply by sampling frequency to restore original units
-    vel = np.insert(convolve(np.diff(pos), win, mode='same'), 0, 0) * freq
-    acc = np.insert(convolve(np.diff(vel), win, mode='same'), 0, 0) * freq
+    vel = np.insert(scipy.signal.convolve(np.diff(pos), win, mode='same'), 0, 0) * freq
+    acc = np.insert(scipy.signal.convolve(np.diff(vel), win, mode='same'), 0, 0) * freq
 
     return vel, acc
 
@@ -272,8 +287,8 @@ def movements(t, pos, freq=1000, pos_thresh=8, t_thresh=.2, min_gap=.1, pos_thre
     peak_amps = np.fromiter(peaks, dtype=float, count=onsets.size)
     N = 10  # Number of points in the Gaussian
     STDEV = 1.8  # Equivalent to a width factor (alpha value) of 2.5
-    gauss = windows.gaussian(N, STDEV)  # A 10-point Gaussian window of a given s.d.
-    vel = convolve(np.diff(np.insert(pos, 0, 0)), gauss, mode='same')
+    gauss = scipy.signal.windows.gaussian(N, STDEV)  # A 10-point Gaussian window of a given s.d.
+    vel = scipy.signal.convolve(np.diff(np.insert(pos, 0, 0)), gauss, mode='same')
     # For each movement period, find the timestamp where the absolute velocity was greatest
     peaks = (t[m + np.abs(vel[m:n]).argmax()] for m, n in zip(onset_samps, offset_samps))
     peak_vel_times = np.fromiter(peaks, dtype=float, count=onsets.size)

brainbox/io/one.py

@@ -21,10 +21,11 @@
 from ibllib.pipes import histology
 from ibllib.pipes.ephys_alignment import EphysAlignment
 
+import brainbox.plot
 from brainbox.core import TimeSeries
 from brainbox.processing import sync
 from brainbox.metrics.single_units import quick_unit_metrics
-from brainbox.behavior.wheel import interpolate_position, velocity_smoothed
+from brainbox.behavior.wheel import interpolate_position, velocity_filtered
 from brainbox.behavior.dlc import likelihood_threshold, get_pupil_diameter, get_smooth_pupil_diameter
 
 _logger = logging.getLogger('ibllib')
@@ -1082,6 +1083,30 @@ def samples2times(self, values, direction='forward'):
             }
         return self._sync[direction](values)
 
+    @property
+    def pid2ref(self):
+        return f"{self.one.eid2ref(self.eid, as_dict=False)}_{self.pname}"
+
+    def raster(self, spikes, save_dir=None):
+        """
+        :param save_dir: optional if specified
+        :return:
+        """
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(figsize=(16, 9))
+        brainbox.plot.driftmap(spikes['times'], spikes['depths'], t_bin=0.007, d_bin=10, vmax=0.5, ax=ax)
+        title_str = f"{self.pid} \n" \
+                    f"{self.pid2ref} \n" \
+                    f"{spikes.clusters.size:_} spikes, {np.unique(spikes.clusters).size:_} clusters"
+        ax.title.set_text(title_str)
+        ax.set_ylim(0, 3800)
+        if save_dir is not None:
+            png_file = save_dir.joinpath(f"{self.pid}_{self.pid2ref}_raster.png") if Path(save_dir).is_dir() else Path(save_dir)
+            fig.savefig(png_file)
+            plt.close(fig)
+        else:
+            return fig, ax
+
 
 @dataclass
 class SessionLoader:
@@ -1253,29 +1278,30 @@ def load_trials(self):
         self.one.wildcards = True
         self.data_info.loc[self.data_info['name'] == 'trials', 'is_loaded'] = True
 
-    def load_wheel(self, sampling_rate=1000, smooth_size=0.03):
+    def load_wheel(self, fs=1000, corner_frequency=20, order=8):
         """
         Function to load wheel data (position, velocity, acceleration) into SessionLoader.wheel. The wheel position
         is first interpolated to a uniform sampling rate. Then velocity and acceleration are computed, during which
-        Gaussian smoothing is applied.
+        a Butterworth low-pass filter is applied.
 
         Parameters
         ----------
-        sampling_rate: float
-            Rate at which to sample the wheel position, default is 1000 Hz
-        smooth_size: float
-            Size of Gaussian smoothing window in seconds, default is 0.03
+        fs: int, float
+            Sampling frequency for the wheel position, default is 1000 Hz
+        corner_frequency: int, float
+            Corner frequency of Butterworth low-pass filter, default is 20
+        order: int, float
+            Order of Butterworth low_pass filter, default is 8
         """
         wheel_raw = self.one.load_object(self.eid, 'wheel')
-        # TODO: Fix this instead of raising error?
         if wheel_raw['position'].shape[0] != wheel_raw['timestamps'].shape[0]:
             raise ValueError("Length mismatch between 'wheel.position' and 'wheel.timestamps")
         # resample the wheel position and compute velocity, acceleration
         self.wheel = pd.DataFrame(columns=['times', 'position', 'velocity', 'acceleration'])
         self.wheel['position'], self.wheel['times'] = interpolate_position(
-            wheel_raw['timestamps'], wheel_raw['position'], freq=sampling_rate)
-        self.wheel['velocity'], self.wheel['acceleration'] = velocity_smoothed(
-            self.wheel['position'], freq=sampling_rate, smooth_size=smooth_size)
+            wheel_raw['timestamps'], wheel_raw['position'], freq=fs)
+        self.wheel['velocity'], self.wheel['acceleration'] = velocity_filtered(
+            self.wheel['position'], fs=fs, corner_frequency=corner_frequency, order=order)
         self.wheel = self.wheel.apply(np.float32)
         self.data_info.loc[self.data_info['name'] == 'wheel', 'is_loaded'] = True
 

brainbox/tests/test_behavior.py

@@ -104,11 +104,10 @@ def test_direction_changes(self):
         on, off, *_ = self.test_data[0][1]
         vel, _ = wheel.velocity_smoothed(pos, 1000)
         times, indices = wheel.direction_changes(t, vel, np.c_[on, off])
-
+        # import matplotlib.pyplot as plt
+        # plt.plot(np.diff(pos) * 1000)
+        # plt.plot(vel)
         self.assertTrue(len(times) == len(indices) == 14, 'incorrect number of arrays returned')
-        # Check first arrays
-        np.testing.assert_allclose(times[0], [21.86593334, 22.12693334, 22.20193334, 22.66093334])
-        np.testing.assert_array_equal(indices[0], [21809, 22070, 22145, 22604])
 
 
 class TestTraining(unittest.TestCase):

ibllib/__init__.py

@@ -1,5 +1,5 @@
 """Library implementing the International Brain Laboratory data pipeline."""
-__version__ = "2.17.1"
+__version__ = "2.17.2"
 import warnings
 
 from iblutil.util import get_logger

ibllib/atlas/atlas.py

@@ -431,12 +431,18 @@ def plot_tilted_slice(self, xyz, axis, volume='image', cmap=None, ax=None, **kwa
         return ax
 
     @staticmethod
-    def _plot_slice(im, extent, ax=None, cmap=None, **kwargs):
+    def _plot_slice(im, extent, ax=None, cmap=None, volume=None, **kwargs):
         if not ax:
             ax = plt.gca()
             ax.axis('equal')
         if not cmap:
             cmap = plt.get_cmap('bone')
+
+        if volume == 'boundary':
+            imb = np.zeros((*im.shape[:2], 4), dtype=np.uint8)
+            imb[im == 1] = np.array([0, 0, 0, 255])
+            im = imb
+
         ax.imshow(im, extent=extent, cmap=cmap, **kwargs)
         return ax
 
@@ -534,9 +540,13 @@ def compute_boundaries(self, values):
         :param values:
         :return:
         """
-        boundary = np.diff(values, axis=0, append=0)
-        boundary = boundary + np.diff(values, axis=1, append=0)
+        boundary = np.abs(np.diff(values, axis=0, prepend=0))
+        boundary = boundary + np.abs(np.diff(values, axis=1, prepend=0))
+        boundary = boundary + np.abs(np.diff(values, axis=1, append=0))
+        boundary = boundary + np.abs(np.diff(values, axis=0, append=0))
+
         boundary[boundary != 0] = 1
+
         return boundary
 
     def plot_slices(self, xyz, *args, **kwargs):
@@ -580,7 +590,7 @@ def plot_cslice(self, ap_coordinate, volume='image', mapping=None, region_values
         """
 
         cslice = self.slice(ap_coordinate, axis=1, volume=volume, mapping=mapping, region_values=region_values)
-        return self._plot_slice(np.moveaxis(cslice, 0, 1), extent=self.extent(axis=1), **kwargs)
+        return self._plot_slice(np.moveaxis(cslice, 0, 1), extent=self.extent(axis=1), volume=volume, **kwargs)
 
     def plot_hslice(self, dv_coordinate, volume='image', mapping=None, region_values=None, **kwargs):
         """
@@ -604,7 +614,7 @@ def plot_hslice(self, dv_coordinate, volume='image', mapping=None, region_values
         """
 
         hslice = self.slice(dv_coordinate, axis=2, volume=volume, mapping=mapping, region_values=region_values)
-        return self._plot_slice(hslice, extent=self.extent(axis=2), **kwargs)
+        return self._plot_slice(hslice, extent=self.extent(axis=2), volume=volume, **kwargs)
 
     def plot_sslice(self, ml_coordinate, volume='image', mapping=None, region_values=None, **kwargs):
         """
@@ -628,7 +638,7 @@ def plot_sslice(self, ml_coordinate, volume='image', mapping=None, region_values
         """
 
         sslice = self.slice(ml_coordinate, axis=0, volume=volume, mapping=mapping, region_values=region_values)
-        return self._plot_slice(np.swapaxes(sslice, 0, 1), extent=self.extent(axis=0), **kwargs)
+        return self._plot_slice(np.swapaxes(sslice, 0, 1), extent=self.extent(axis=0), volume=volume, **kwargs)
 
     def plot_top(self, volume='annotation', mapping=None, region_values=None, ax=None, **kwargs):
         """
@@ -668,7 +678,7 @@ def plot_top(self, volume='annotation', mapping=None, region_values=None, ax=Non
         elif volume == 'boundary':
             im = self.compute_boundaries(regions)
 
-        return self._plot_slice(im, self.extent(axis=2), ax=ax, **kwargs)
+        return self._plot_slice(im, self.extent(axis=2), ax=ax, volume=volume, **kwargs)
 
 
 @dataclass
@@ -1164,7 +1174,7 @@ def plot_flatmap(self, depth=0, volume='annotation', mapping='Allen', region_val
         if not ax:
             ax = plt.gca()
 
-        return self._plot_slice(im, self.extent_flmap(), ax=ax, **kwargs)
+        return self._plot_slice(im, self.extent_flmap(), ax=ax, volume=volume, **kwargs)
 
     def extent_flmap(self):
         extent = np.r_[0, self.flatmap.shape[1], 0, self.flatmap.shape[0]]

ibllib/atlas/plots.py

@@ -344,12 +344,6 @@ def _plot_slice(coord, slice, region_values, vol_type, background='boundary', ma
     else:
         fig, ax = plt.subplots()
 
-    if background == 'boundary':
-        cmap_bound = matplotlib.cm.get_cmap("bone_r").copy()
-        cmap_bound.set_under([1, 1, 1], 0)
-    else:
-        cmap_bound = None
-
     if slice == 'coronal':
         if background == 'image':
             ba.plot_cslice(coord, volume='image', mapping=map, ax=ax)
@@ -358,7 +352,7 @@ def _plot_slice(coord, slice, region_values, vol_type, background='boundary', ma
         else:
             ba.plot_cslice(coord, volume=vol_type, region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
                            vmax=clevels[1], ax=ax)
-            ba.plot_cslice(coord, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+            ba.plot_cslice(coord, volume='boundary', mapping=map, ax=ax)
 
     elif slice == 'sagittal':
         if background == 'image':
@@ -368,7 +362,7 @@ def _plot_slice(coord, slice, region_values, vol_type, background='boundary', ma
         else:
             ba.plot_sslice(coord, volume=vol_type, region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
                            vmax=clevels[1], ax=ax)
-            ba.plot_sslice(coord, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+            ba.plot_sslice(coord, volume='boundary', mapping=map, ax=ax)
 
     elif slice == 'horizontal':
         if background == 'image':
@@ -378,7 +372,7 @@ def _plot_slice(coord, slice, region_values, vol_type, background='boundary', ma
         else:
             ba.plot_hslice(coord, volume=vol_type, region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
                            vmax=clevels[1], ax=ax)
-            ba.plot_hslice(coord, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+            ba.plot_hslice(coord, volume='boundary', mapping=map, ax=ax)
 
     elif slice == 'top':
         if background == 'image':
@@ -388,8 +382,7 @@ def _plot_slice(coord, slice, region_values, vol_type, background='boundary', ma
         else:
             ba.plot_top(volume=vol_type, region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
                         vmax=clevels[1], ax=ax)
-            ba.plot_top(volume='boundary', mapping=map, ax=ax,
-                        cmap=cmap_bound, vmin=0.01, vmax=0.8)
+            ba.plot_top(volume='boundary', mapping=map, ax=ax)
 
     if show_cbar:
         norm = matplotlib.colors.Normalize(vmin=clevels[0], vmax=clevels[1], clip=False)

ibllib/ephys/ephysqc.py

@@ -209,7 +209,7 @@ def rmsmap(sglx):
            'tscale': wingen.tscale(fs=sglx.fs)}
     win['spectral_density'] = np.zeros((len(win['fscale']), sglx.nc))
     # loop through the whole session
-    with tqdm(total=wingen.firstlast) as pbar:
+    with tqdm(total=wingen.nwin) as pbar:
         for first, last in wingen.firstlast:
             D = sglx.read_samples(first_sample=first, last_sample=last)[0].transpose()
             # remove low frequency noise below 1 Hz

ibllib/io/extractors/camera.py

@@ -160,7 +160,17 @@ def _extract(self, sync=None, chmap=None, video_path=None,
         video_frames = get_video_length(self.session_path.joinpath('raw_video_data', f'_iblrig_{self.label}Camera.raw.mp4'))
         raw_ts = fpga_times[self.label]
 
-        assert video_frames == raw_ts.size, 'dimension mismatch between video frames and TTL pulses'
+        # For left camera sometimes we have one extra pulse than video frame
+        if (raw_ts.size - video_frames) == 1:
+            _logger.warning(f'One extra sync pulse detected for {self.label} camera')
+            raw_ts = raw_ts[:-1]
+        elif (raw_ts.size - video_frames) == -1:
+            _logger.warning(f'One extra video frame detected for {self.label} camera')
+            med_time = np.median(np.diff(raw_ts))
+            raw_ts = np.r_[raw_ts, np.array([raw_ts[-1] + med_time])]
+
+        assert video_frames == raw_ts.size, f'dimension mismatch between video frames and TTL pulses for {self.label} camera' \
+                                            f'by {np.abs(video_frames - raw_ts.size)} frames'
 
         return raw_ts
 

ibllib/io/extractors/video_motion.py

@@ -136,6 +136,34 @@ def _set_eid_or_path(self, session_path_or_eid):
             raise ValueError("'session' must be a valid session path or uuid")
 
     def align_motion(self, period=(-np.inf, np.inf), side='left', sd_thresh=10, display=False):
+        """
+        Align video to the wheel using cross-correlation of the video motion signal and the rotary
+        encoder.
+
+        Parameters
+        ----------
+        period : (float, float)
+            The time period over which to do the alignment.
+        side : {'left', 'right'}
+            With which camera to perform the alignment.
+        sd_thresh : float
+            For plotting where the motion energy goes above this standard deviation threshold.
+        display : bool
+            When true, displays the aligned wheel motion energy along with the rotary encoder
+            signal.
+
+        Returns
+        -------
+        int
+            Frame offset, i.e. by how many frames the video was shifted to match the rotary encoder
+            signal.  Negative values mean the video was shifted backwards with respect to the wheel
+            timestamps.
+        float
+            The peak cross-correlation.
+        numpy.ndarray
+            The motion energy used in the cross-correlation, i.e. the frame difference for the
+            period given.
+        """
         # Get data samples within period
         wheel = self.data['wheel']
         self.alignment.label = side

ibllib/io/extractors/widefield.py

@@ -121,7 +121,7 @@ def _save(self, data=None, path_out=None):
     def preprocess(self, fs=30, functional_channel=0, nbaseline_frames=30, k=200, nchannels=2):
 
         # MOTION CORRECTION
-        wfield_cli._motion(str(self.data_path), nchannels=nchannels)
+        wfield_cli._motion(str(self.data_path), nchannels=nchannels, plot_ext='.png')
         # COMPUTE AVERAGE FOR BASELINE
         wfield_cli._baseline(str(self.data_path), nbaseline_frames, nchannels=nchannels)
         # DATA REDUCTION
@@ -131,23 +131,23 @@ def preprocess(self, fs=30, functional_channel=0, nbaseline_frames=30, k=200, nc
         dat = wfield_cli.load_stack(str(self.data_path), nchannels=nchannels)
         if dat.shape[1] == 2:
             del dat
-            wfield_cli._hemocorrect(str(self.data_path), fs=fs, functional_channel=functional_channel)
+            wfield_cli._hemocorrect(str(self.data_path), fs=fs, functional_channel=functional_channel, plot_ext='.png')
 
     def remove_files(self, file_prefix='motion'):
         motion_files = self.data_path.glob(f'{file_prefix}*')
         for file in motion_files:
             _logger.info(f'Removing {file}')
             file.unlink()
 
-    def sync_timestamps(self, bin_exists=False, save=False, save_paths=None, **kwargs):
+    def sync_timestamps(self, bin_exists=False, save=False, save_paths=None, sync_collection='raw_sync_data', **kwargs):
 
         if save and save_paths:
             assert len(save_paths) == 3, 'Must provide save_path as list with 3 paths'
             for save_path in save_paths:
                 Path(save_path).parent.mkdir(parents=True, exist_ok=True)
 
         # Load in fpga sync
-        fpga_sync, chmap = get_sync_and_chn_map(self.session_path, 'raw_widefield_data')
+        fpga_sync, chmap = get_sync_and_chn_map(self.session_path, sync_collection)
         fpga_led = get_sync_fronts(fpga_sync, chmap['frame_trigger'])
         fpga_led_up = fpga_led['times'][fpga_led['polarities'] == 1]  # only consider up pulse times