Merge branch 'hotfix/2.0.7'

Author: oliche

.flake8

@@ -1,5 +1,5 @@
 [flake8]
-max-line-length = 99
+max-line-length = 130
 ignore = W504, W503, E266
 exclude =
     .git,

brainbox/plot.py

@@ -618,7 +618,7 @@ def peri_event_time_histogram(
 
 
 def driftmap(ts, feat, ax=None, plot_style='bincount',
-             t_bin=0.01, d_bin=20, weights=None, **kwargs):
+             t_bin=0.01, d_bin=20, weights=None, vmax=None, **kwargs):
     """
     Plots the values of a spike feature array (y-axis) over time (x-axis).
     Two arguments can be given for the plot_style of the drift map:
@@ -636,6 +636,7 @@ def driftmap(ts, feat, ax=None, plot_style='bincount',
     t_bin: time bin used when plot_style='bincount'
     d_bin: depth bin used when plot_style='bincount'
     plot_style: 'scatter', 'bincount'
+    **kwargs: matplotlib.imshow arguments
 
     Returns
     -------
@@ -674,7 +675,7 @@ def driftmap(ts, feat, ax=None, plot_style='bincount',
         R, times, depths = bincount2D(
             ts[iok], feat[iok], t_bin, d_bin, weights=weights)
         # plot raster map
-        ax.imshow(R, aspect='auto', cmap='binary', vmin=0, vmax=np.std(R) * 4,
+        ax.imshow(R, aspect='auto', cmap='binary', vmin=0, vmax=vmax or np.std(R) * 4,
                   extent=np.r_[times[[0, -1]], depths[[0, -1]]], origin='lower', **kwargs)
     ax.set_xlabel('time (secs)')
     ax.set_ylabel('depth (um)')

brainbox/task/passive.py

@@ -115,12 +115,17 @@ def get_rf_map_over_depth(rf_map_times, rf_map_pos, rf_stim_frames, spike_times,
             stim_on_times = rf_map_times[stim_frame[0]]
             stim_intervals = np.c_[stim_on_times - pre_stim, stim_on_times + post_stim]
 
-            idx_intervals = np.searchsorted(times, stim_intervals)
-
-            stim_trials = np.zeros((depths.shape[0], n_bins, idx_intervals.shape[0]))
-            for i, on in enumerate(idx_intervals):
-                stim_trials[:, :, i] = binned_array[:, on[0]:on[1]]
-            avg_stim_trials = np.mean(stim_trials, axis=2)
+            out_intervals = stim_intervals[:, 1] > times[-1]
+            idx_intervals = np.searchsorted(times, stim_intervals)[np.invert(out_intervals)]
+
+            # Case when no spikes during the passive period
+            if idx_intervals.shape[0] == 0:
+                avg_stim_trials = np.zeros((depths.shape[0], n_bins))
+            else:
+                stim_trials = np.zeros((depths.shape[0], n_bins, idx_intervals.shape[0]))
+                for i, on in enumerate(idx_intervals):
+                    stim_trials[:, :, i] = binned_array[:, on[0]:on[1]]
+                avg_stim_trials = np.mean(stim_trials, axis=2)
 
             _rf_map[:, x_pos, y_pos, :] = avg_stim_trials
 
@@ -197,8 +202,10 @@ def get_stim_aligned_activity(stim_events, spike_times, spike_depths, z_score_fl
 
         stim_intervals = np.c_[stim_times - pre_stim, stim_times + post_stim]
         base_intervals = np.c_[stim_times - base_stim, stim_times - pre_stim]
-        idx_stim = np.searchsorted(times, stim_intervals)
-        idx_base = np.searchsorted(times, base_intervals)
+        out_intervals = stim_intervals[:, 1] > times[-1]
+
+        idx_stim = np.searchsorted(times, stim_intervals)[np.invert(out_intervals)]
+        idx_base = np.searchsorted(times, base_intervals)[np.invert(out_intervals)]
 
         stim_trials = np.zeros((depths.shape[0], n_bins, idx_stim.shape[0]))
         noise_trials = np.zeros((depths.shape[0], n_bins_base, idx_stim.shape[0]))

ibllib/io/extractors/camera.py

@@ -91,6 +91,14 @@ def _extract(self, sync=None, chmap=None, video_path=None,
         """
         fpga_times = extract_camera_sync(sync=sync, chmap=chmap)
         count, (*_, gpio) = raw.load_embedded_frame_data(self.session_path, self.label)
+        raw_ts = fpga_times[self.label]
+
+        if video_path is None:
+            filename = f'_iblrig_{self.label}Camera.raw.mp4'
+            video_path = self.session_path.joinpath('raw_video_data', filename)
+        # Permit the video path to be the length for development and debugging purposes
+        length = (video_path if isinstance(video_path, int) else get_video_length(video_path))
+        _logger.debug(f'Number of video frames = {length}')
 
         if gpio is not None and gpio['indices'].size > 1:
             _logger.info('Aligning to audio TTLs')
@@ -111,19 +119,10 @@ def _extract(self, sync=None, chmap=None, video_path=None,
                 right at the end of the video.  We therefore simply shorten the arrays to match
                 the length of the video.
                 """
-                if video_path is None:
-                    filename = f'_iblrig_{self.label}Camera.raw.mp4'
-                    video_path = self.session_path.joinpath('raw_video_data', filename)
-                # Permit the video path to be the length for development and debugging purposes
-                length = (video_path
-                          if isinstance(video_path, int)
-                          else get_video_length(video_path))
-                _logger.debug(f'Number of video frames = {length}')
                 if count.size > length:
                     count = count[:length]
                 else:
                     assert length == count.size, 'fewer counts than frames'
-                raw_ts = fpga_times[self.label]
                 return align_with_audio(raw_ts, audio, gpio, count,
                                         display=display,
                                         extrapolate_missing=extrapolate_missing)
@@ -132,7 +131,11 @@ def _extract(self, sync=None, chmap=None, video_path=None,
 
         # If you reach here extracting using audio TTLs was not possible
         _logger.warning('Alignment by wheel data not yet implemented')
-        return fpga_times[self.label]
+        if length < raw_ts.size:
+            df = raw_ts.size - length
+            _logger.info(f'Discarding first {df} pulses')
+            raw_ts = raw_ts[df:]
+        return raw_ts
 
 
 class CameraTimestampsBpod(BaseBpodTrialsExtractor):
@@ -219,6 +222,8 @@ def _extract(self, video_path=None, display=False, extrapolate_missing=True):
             raw_ts = np.r_[raw_ts, to_app]  # Append the missing times
         elif n_missing < 0:
             _logger.warning(f'{abs(n_missing)} fewer frames than Bpod timestamps')
+            _logger.info(f'Discarding first {abs(n_missing)} pulses')
+            raw_ts = raw_ts[abs(n_missing):]
 
         return raw_ts
 

ibllib/io/extractors/ephys_fpga.py

@@ -245,28 +245,34 @@ def _rotary_encoder_positions_from_fronts(ta, pa, tb, pb, ticks=WHEEL_TICKS, rad
         return t, p
 
 
-def _assign_events_audio(audio_t, audio_polarities, return_indices=False):
+def _assign_events_audio(audio_t, audio_polarities, return_indices=False, display=False):
     """
     From detected fronts on the audio sync traces, outputs the synchronisation events
     related to tone in
 
     :param audio_t: numpy vector containing times of fronts
     :param audio_fronts: numpy vector containing polarity of fronts (1 rise, -1 fall)
     :param return_indices (False): returns indices of tones
+    :param display (False): for debug mode, displays the raw fronts overlaid with detections
     :return: numpy arrays t_ready_tone_in, t_error_tone_in
     :return: numpy arrays ind_ready_tone_in, ind_error_tone_in if return_indices=True
     """
     # make sure that there are no 2 consecutive fall or consecutive rise events
     assert(np.all(np.abs(np.diff(audio_polarities)) == 2))
     # take only even time differences: ie. from rising to falling fronts
-    i0 = 0 if audio_polarities[0] == 1 else 1
-    dt = np.diff(audio_t)[i0::2]
+    dt = np.diff(audio_t)
     # detect ready tone by length below 110 ms
-    i_ready_tone_in = np.r_[np.where(dt <= 0.11)[0] * 2]
+    i_ready_tone_in = np.where(np.logical_and(dt <= 0.11, audio_polarities[:-1] == 1))[0]
     t_ready_tone_in = audio_t[i_ready_tone_in]
-    # error tones are events lasting from 400ms to 600ms
-    i_error_tone_in = np.where(np.logical_and(0.4 < dt, dt < 1.2))[0] * 2
+    # error tones are events lasting from 400ms to 1200ms
+    i_error_tone_in = np.where(np.logical_and(np.logical_and(0.4 < dt, dt < 1.2), audio_polarities[:-1] == 1))[0]
     t_error_tone_in = audio_t[i_error_tone_in]
+    if display:  # pragma: no cover
+        from ibllib.plots import squares, vertical_lines
+        squares(audio_t, audio_polarities, yrange=[-1, 1],)
+        vertical_lines(t_ready_tone_in, ymin=-.8, ymax=.8)
+        vertical_lines(t_error_tone_in, ymin=-.8, ymax=.8)
+
     if return_indices:
         return t_ready_tone_in, t_error_tone_in, i_ready_tone_in, i_error_tone_in
     else:
@@ -280,7 +286,6 @@ def _assign_events_to_trial(t_trial_start, t_event, take='last'):
     Trials without an event
     result in nan value in output time vector.
     The output has a consistent size with t_trial_start and ready to output to alf.
-
     :param t_trial_start: numpy vector of trial start times
     :param t_event: numpy vector of event times to assign to trials
     :param take: 'last' or 'first' (optional, default 'last'): index to take in case of duplicates
@@ -325,6 +330,37 @@ def get_sync_fronts(sync, channel_nb, tmin=None, tmax=None):
                   'polarities': sync['polarities'][selection]})
 
 
+def _clean_audio(audio, display=False):
+    """
+    one guy wired the 150 Hz camera output onto the soundcard. The effect is to get 150 Hz periodic
+    square pulses, 2ms up and 4.666 ms down. When this happens we remove all of the intermediate
+    pulses to repair the audio trace
+    Here is some helper code
+        dd = np.diff(audio['times'])
+        1 / np.median(dd[::2]) # 2ms up
+        1 / np.median(dd[1::2])  # 4.666 ms down
+        1 / (np.median(dd[::2]) + np.median(dd[1::2])) # both sum to 150 Hx
+    This only runs on sessions when the bug is detected and leaves others untouched
+    """
+    DISCARD_THRESHOLD = 0.01
+    average_150_hz = np.mean(1 / np.diff(audio['times'][audio['polarities'] == 1]) > 140)
+    naudio = audio['times'].size
+    if average_150_hz > 0.7 and naudio > 100:
+        _logger.warning("Soundcard signal on FPGA seems to have been mixed with 150Hz camera")
+        keep_ind = np.r_[np.diff(audio['times']) > DISCARD_THRESHOLD, False]
+        keep_ind = np.logical_and(keep_ind, audio['polarities'] == -1)
+        keep_ind = np.where(keep_ind)[0]
+        keep_ind = np.sort(np.r_[0, keep_ind, keep_ind + 1, naudio - 1])
+
+        if display:  # pragma: no cover
+            from ibllib.plots import squares
+            squares(audio['times'], audio['polarities'], ax=None, yrange=[-1, 1])
+            squares(audio['times'][keep_ind], audio['polarities'][keep_ind], yrange=[-1, 1])
+        audio = {'times': audio['times'][keep_ind],
+                 'polarities': audio['polarities'][keep_ind]}
+    return audio
+
+
 def _clean_frame2ttl(frame2ttl, display=False):
     """
     Frame 2ttl calibration can be unstable and the fronts may be flickering at an unrealistic
@@ -387,6 +423,7 @@ def extract_behaviour_sync(sync, chmap=None, display=False, bpod_trials=None):
     frame2ttl = get_sync_fronts(sync, chmap['frame2ttl'])
     frame2ttl = _clean_frame2ttl(frame2ttl)
     audio = get_sync_fronts(sync, chmap['audio'])
+    audio = _clean_audio(audio)
     # extract events from the fronts for each trace
     t_trial_start, t_valve_open, t_iti_in = _assign_events_bpod(bpod['times'], bpod['polarities'])
     # one issue is that sometimes bpod pulses may not have been detected, in this case

ibllib/io/extractors/ephys_passive.py

@@ -79,6 +79,20 @@ def _load_passive_session_fixtures(session_path: str) -> dict:
     return fixture
 
 
+def _load_task_protocol(session_path: str) -> str:
+    """Find the IBL rig version used for the session
+
+    :param session_path: the path to a session
+    :type session_path: str
+    :return: ibl rig task protocol version
+    :rtype: str
+    """
+    settings = rawio.load_settings(session_path)
+    ses_ver = settings["IBLRIG_VERSION_TAG"]
+
+    return ses_ver
+
+
 def _load_passive_stim_meta() -> dict:
     """load_passive_stim_meta Loads the passive protocol metadata
 
@@ -321,38 +335,72 @@ def _extract_passiveValve_intervals(bpod: dict) -> np.array:
 
     # check all values are within bpod tolerance of 100µs
     assert np.allclose(
-        valveOff_times - valveOn_times, valveOff_times[0] - valveOn_times[0], atol=0.0001
+        valveOff_times - valveOn_times, valveOff_times[1] - valveOn_times[1], atol=0.0001
     ), "Some valve outputs are longer or shorter than others"
 
     return np.array([(x, y) for x, y in zip(valveOn_times, valveOff_times)])
 
 
-def _extract_passiveAudio_intervals(audio: dict) -> Tuple[np.array, np.array]:
-    # Get Tone and Noise cue intervals
+def _extract_passiveAudio_intervals(audio: dict, rig_version: str) -> Tuple[np.array, np.array]:
+
+    # make an exception for task version = 6.2.5 where things are strange but data is recoverable
+    if rig_version == '6.2.5':
+        # Get all sound onsets and offsets
+        soundOn_times = audio["times"][audio["polarities"] > 0]
+        soundOff_times = audio["times"][audio["polarities"] < 0]
 
-    # Get all sound onsets and offsets
-    soundOn_times = audio["times"][audio["polarities"] > 0]
-    soundOff_times = audio["times"][audio["polarities"] < 0]
-    # Check they are the correct number
-    assert len(soundOn_times) == NTONES + NNOISES, "Wrong number of sound ONSETS"
-    assert len(soundOff_times) == NTONES + NNOISES, "Wrong number of sound OFFSETS"
+        # Have a couple that are wayyy too long!
+        time_threshold = 10
+        diff = soundOff_times - soundOn_times
+        stupid = np.where(diff > time_threshold)[0]
+        NREMOVE = len(stupid)
+        not_stupid = np.where(diff < time_threshold)[0]
 
-    diff = soundOff_times - soundOn_times
-    # Tone is ~100ms so check if diff < 0.3
-    toneOn_times = soundOn_times[diff <= 0.3]
-    toneOff_times = soundOff_times[diff <= 0.3]
-    # Noise is ~500ms so check if diff > 0.3
-    noiseOn_times = soundOn_times[diff > 0.3]
-    noiseOff_times = soundOff_times[diff > 0.3]
+        assert len(soundOn_times) == NTONES + NNOISES - NREMOVE, "Wrong number of sound ONSETS"
+        assert len(soundOff_times) == NTONES + NNOISES - NREMOVE, "Wrong number of sound OFFSETS"
 
-    assert len(toneOn_times) == NTONES
-    assert len(toneOff_times) == NTONES
-    assert len(noiseOn_times) == NNOISES
-    assert len(noiseOff_times) == NNOISES
+        soundOn_times = soundOn_times[not_stupid]
+        soundOff_times = soundOff_times[not_stupid]
 
-    # Fixed delays from soundcard ~500µs
-    np.allclose(toneOff_times - toneOn_times, 0.1, atol=0.0006)
-    np.allclose(noiseOff_times - noiseOn_times, 0.5, atol=0.0006)
+        diff = soundOff_times - soundOn_times
+        # Tone is ~100ms so check if diff < 0.3
+        toneOn_times = soundOn_times[diff <= 0.3]
+        toneOff_times = soundOff_times[diff <= 0.3]
+        # Noise is ~500ms so check if diff > 0.3
+        noiseOn_times = soundOn_times[diff > 0.3]
+        noiseOff_times = soundOff_times[diff > 0.3]
+
+        # append with nans
+        toneOn_times = np.r_[toneOn_times, np.full((NTONES - len(toneOn_times)), np.NAN)]
+        toneOff_times = np.r_[toneOff_times, np.full((NTONES - len(toneOff_times)), np.NAN)]
+        noiseOn_times = np.r_[noiseOn_times, np.full((NNOISES - len(noiseOn_times)), np.NAN)]
+        noiseOff_times = np.r_[noiseOff_times, np.full((NNOISES - len(noiseOff_times)), np.NAN)]
+
+    else:
+        # Get all sound onsets and offsets
+        soundOn_times = audio["times"][audio["polarities"] > 0]
+        soundOff_times = audio["times"][audio["polarities"] < 0]
+
+        # Check they are the correct number
+        assert len(soundOn_times) == NTONES + NNOISES, "Wrong number of sound ONSETS"
+        assert len(soundOff_times) == NTONES + NNOISES, "Wrong number of sound OFFSETS"
+
+        diff = soundOff_times - soundOn_times
+        # Tone is ~100ms so check if diff < 0.3
+        toneOn_times = soundOn_times[diff <= 0.3]
+        toneOff_times = soundOff_times[diff <= 0.3]
+        # Noise is ~500ms so check if diff > 0.3
+        noiseOn_times = soundOn_times[diff > 0.3]
+        noiseOff_times = soundOff_times[diff > 0.3]
+
+        assert len(toneOn_times) == NTONES
+        assert len(toneOff_times) == NTONES
+        assert len(noiseOn_times) == NNOISES
+        assert len(noiseOff_times) == NNOISES
+
+        # Fixed delays from soundcard ~500µs
+        np.allclose(toneOff_times - toneOn_times, 0.1, atol=0.0006)
+        np.allclose(noiseOff_times - noiseOn_times, 0.5, atol=0.0006)
 
     passiveTone_intervals = np.append(
         toneOn_times.reshape((len(toneOn_times), 1)),
@@ -444,8 +492,10 @@ def extract_task_replay(
     bpod = ephys_fpga.get_sync_fronts(sync, sync_map["bpod"], tmin=treplay[0])
     passiveValve_intervals = _extract_passiveValve_intervals(bpod)
 
+    task_version = _load_task_protocol(session_path)
     audio = ephys_fpga.get_sync_fronts(sync, sync_map["audio"], tmin=treplay[0])
-    passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio)
+    passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio,
+                                                                                    task_version)
 
     passiveStims_df = np.concatenate(
         [passiveValve_intervals, passiveTone_intervals, passiveNoise_intervals], axis=1
@@ -493,8 +543,10 @@ def extract_replay_debug(
     passiveValve_intervals = _extract_passiveValve_intervals(bpod)
     plot_valve_times(passiveValve_intervals, ax=ax)
 
+    task_version = _load_task_protocol(session_path)
     audio = ephys_fpga.get_sync_fronts(sync, sync_map["audio"], tmin=treplay[0])
-    passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio)
+    passiveTone_intervals, passiveNoise_intervals = _extract_passiveAudio_intervals(audio,
+                                                                                    task_version)
     plot_audio_times(passiveTone_intervals, passiveNoise_intervals, ax=ax)
 
     passiveStims_df = np.concatenate(

ibllib/io/extractors/video_motion.py

@@ -154,6 +154,7 @@ def align_motion(self, period=(-np.inf, np.inf), side='left', sd_thresh=10, disp
             # TODO Add function arg to make grayscale
             self.alignment.frames = \
                 vidio.get_video_frames_preload(camera_path, frame_numbers, mask=roi)
+            assert self.alignment.frames.size != 0
         except AssertionError:
             self.log.error('Failed to open video')
             return None, None, None

ibllib/pipes/ephys_alignment.py

@@ -35,6 +35,13 @@ def __init__(self, xyz_picks, chn_depths=None, track_prev=None,
         self.xyz_samples = histology.interpolate_along_track(self.xyz_track,
                                                              self.sampling_trk -
                                                              self.sampling_trk[0])
+        # ensure none of the track is outside the y or x lim of atlas
+        xlim = np.bitwise_and(self.xyz_samples[:, 0] > self.brain_atlas.bc.xlim[0],
+                              self.xyz_samples[:, 0] < self.brain_atlas.bc.xlim[1])
+        ylim = np.bitwise_and(self.xyz_samples[:, 1] < self.brain_atlas.bc.ylim[0],
+                              self.xyz_samples[:, 1] > self.brain_atlas.bc.ylim[1])
+        rem = np.bitwise_and(xlim, ylim)
+        self.xyz_samples = self.xyz_samples[rem]
 
         self.region, self.region_label, self.region_colour, self.region_id\
             = self.get_histology_regions(self.xyz_samples, self.sampling_trk, self.brain_atlas)