Merge branch 'develop' into single_source_version

Author: nbonacchi

CITATION.cff

@@ -0,0 +1,11 @@
+cff-version: 0.0.0
+message: "If you use this software, please cite it as below."
+authors:
+  - family-names: International Brain Laboratory
+    given-names: The
+    orcid: 
+title: "ibllib"
+version: 
+doi: 
+date-released: 2021-12-09
+url: "https://github.com/int-brain-lab/ibllib"

brainbox/behavior/dlc.py

@@ -4,7 +4,7 @@
 import numpy as np
 from numpy import pi
 import scipy.interpolate as interpolate
-from scipy.signal import convolve, gaussian
+from scipy.signal import convolve, windows
 from scipy.linalg import hankel
 import matplotlib.pyplot as plt
 from matplotlib.collections import LineCollection
@@ -114,7 +114,7 @@ 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 = gaussian(N, gauss_std)
+    win = windows.gaussian(N, gauss_std)
     win = win / win.sum()  # Normalize amplitude
 
     # Convolve and multiply by sampling frequency to restore original units
@@ -274,7 +274,7 @@ 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 = gaussian(N, STDEV)  # A 10-point Gaussian window of a given s.d.
+    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')
     # 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))

brainbox/behavior/wheel.py

@@ -12,7 +12,7 @@
 
 from ibllib.io import spikeglx
 from ibllib.io.extractors.training_wheel import extract_wheel_moves, extract_first_movement_times
-from ibllib.ephys.neuropixel import SITES_COORDINATES, TIP_SIZE_UM
+from ibllib.ephys.neuropixel import SITES_COORDINATES, TIP_SIZE_UM, trace_header
 from ibllib.atlas import atlas
 from ibllib.atlas import AllenAtlas
 from ibllib.pipes import histology
@@ -221,7 +221,7 @@ def _load_channels_locations_from_disk(eid, collection=None, one=None, revision=
     return channels
 
 
-def channel_locations_interpolation(channels_aligned, channels, brain_regions=None):
+def channel_locations_interpolation(channels_aligned, channels=None, brain_regions=None):
     """
     oftentimes the channel map for different spike sorters may be different so interpolate the alignment onto
     if there is no spike sorting in the base folder, the alignment doesn't have the localCoordinates field
@@ -238,16 +238,18 @@ def channel_locations_interpolation(channels_aligned, channels, brain_regions=No
       'x', 'y', 'z', 'acronym', 'atlas_id', 'axial_um', 'lateral_um'
     :return: Bunch or dictionary of channels with brain coordinates keys
     """
+    NEUROPIXEL_VERSION = 1
+    h = trace_header(version=NEUROPIXEL_VERSION)
+    if channels is None:
+        channels = {'localCoordinates': np.c_[h['x'], h['y']]}
     nch = channels['localCoordinates'].shape[0]
     if set(['x', 'y', 'z']).issubset(set(channels_aligned.keys())):
         channels_aligned = _channels_bunch2alf(channels_aligned)
     if 'localCoordinates' in channels_aligned.keys():
         aligned_depths = channels_aligned['localCoordinates'][:, 1]
     else:  # this is a edge case for a few spike sorting sessions
         assert channels_aligned['mlapdv'].shape[0] == 384
-        NEUROPIXEL_VERSION = 1
-        from ibllib.ephys.neuropixel import trace_header
-        aligned_depths = trace_header(version=NEUROPIXEL_VERSION)['y']
+        aligned_depths = h['y']
     depth_aligned, ind_aligned = np.unique(aligned_depths, return_index=True)
     depths, ind, iinv = np.unique(channels['localCoordinates'][:, 1], return_index=True, return_inverse=True)
     channels['mlapdv'] = np.zeros((nch, 3))

brainbox/io/one.py

@@ -781,7 +781,7 @@ def xyz2ccf(self, xyz, ccf_order='mlapdv'):
         """
         Converts coordinates to the CCF coordinates, which is assumed to be the cube indices
         times the spacing.
-        :param xyz: mlapdv coordinates in um, origin Bregma
+        :param xyz: mlapdv coordinates in meters, origin Bregma
         :param ccf_order: order that you want values returned 'mlapdv' (ibl) or 'apdvml'
         (Allen mcc vertices)
         :return: coordinates in CCF space um, origin is the front left top corner of the data

ibllib/atlas/atlas.py

@@ -109,6 +109,37 @@ def fk(x, si=.002, dx=1, vbounds=None, btype='highpass', ntr_pad=0, ntr_tap=None
     return xf / gain
 
 
+def car(x, collection=None, lagc=300, butter_kwargs=None):
+    """
+    Applies common average referencing with optional automatic gain control
+    :param x: the input array to be filtered. dimension, the filtering is considering
+    axis=0: spatial dimension, axis=1 temporal dimension. (ntraces, ns)
+    :param collection:
+    :param lagc: window size for time domain automatic gain control (no agc otherwise)
+    :param butter_kwargs: filtering parameters: defaults: {'N': 3, 'Wn': 0.1, 'btype': 'highpass'}
+    :return:
+    """
+    if butter_kwargs is None:
+        butter_kwargs = {'N': 3, 'Wn': 0.1, 'btype': 'highpass'}
+    if collection is not None:
+        xout = np.zeros_like(x)
+        for c in np.unique(collection):
+            sel = collection == c
+            xout[sel, :] = kfilt(x=x[sel, :], ntr_pad=0, ntr_tap=None, collection=None,
+                                 butter_kwargs=butter_kwargs)
+        return xout
+
+    # apply agc and keep the gain in handy
+    if not lagc:
+        xf = np.copy(x)
+        gain = 1
+    else:
+        xf, gain = agc(x, wl=lagc, si=1.0)
+    # apply CAR and then un-apply the gain
+    xf = xf - np.median(xf, axis=0)
+    return xf / gain
+
+
 def kfilt(x, collection=None, ntr_pad=0, ntr_tap=None, lagc=300, butter_kwargs=None):
     """
     Applies a butterworth filter on the 0-axis with tapering / padding
@@ -209,13 +240,19 @@ def destripe(x, fs, neuropixel_version=1, butter_kwargs=None, k_kwargs=None, cha
        True: deduces the bad channels from the data provided
     :param butter_kwargs: (optional, None) butterworth params, see the code for the defaults dict
     :param k_kwargs: (optional, None) K-filter params, see the code for the defaults dict
+        can also be set to 'car', in which case the median accross channels will be subtracted
     :return: x, filtered array
     """
     if butter_kwargs is None:
         butter_kwargs = {'N': 3, 'Wn': 300 / fs * 2, 'btype': 'highpass'}
     if k_kwargs is None:
         k_kwargs = {'ntr_pad': 60, 'ntr_tap': 0, 'lagc': 3000,
                     'butter_kwargs': {'N': 3, 'Wn': 0.01, 'btype': 'highpass'}}
+        spatial_fcn = lambda dat: kfilt(dat, **k_kwargs)  # noqa
+    elif isinstance(k_kwargs, dict):
+        spatial_fcn = lambda dat: kfilt(dat, **k_kwargs)  # noqa
+    else:
+        spatial_fcn = lambda dat: car(dat, lagc=int(0.1 * fs))  # noqa
     h = neuropixel.trace_header(version=neuropixel_version)
     if channel_labels is True:
         channel_labels, _ = detect_bad_channels(x, fs)
@@ -231,9 +268,9 @@ def destripe(x, fs, neuropixel_version=1, butter_kwargs=None, k_kwargs=None, cha
     if channel_labels is not None:
         x = interpolate_bad_channels(x, channel_labels, h)
         inside_brain = np.where(channel_labels != 3)[0]
-        x[inside_brain, :] = kfilt(x[inside_brain, :], **k_kwargs)  # apply the k-filter
+        x[inside_brain, :] = spatial_fcn(x[inside_brain, :])  # apply the k-filter
     else:
-        x = kfilt(x, **k_kwargs)
+        x = spatial_fcn(x)
     return x
 
 
@@ -245,7 +282,7 @@ def decompress_destripe_cbin(sr_file, output_file=None, h=None, wrot=None, appen
     Production version with optimized FFTs - requires pyfftw
     :param sr: seismic reader object (spikeglx.Reader)
     :param output_file: (optional, defaults to .bin extension of the compressed bin file)
-    :param h: (optional)
+    :param h: (optional) neuropixel trace header. Dictionary with key 'sample_shift'
     :param wrot: (optional) whitening matrix [nc x nc] or amplitude scalar to apply to the output
     :param append: (optional, False) for chronic recordings, append to end of file
     :param nc_out: (optional, True) saves non selected channels (synchronisation trace) in output
@@ -273,7 +310,7 @@ def decompress_destripe_cbin(sr_file, output_file=None, h=None, wrot=None, appen
     k_kwargs = {'ntr_pad': 60, 'ntr_tap': 0, 'lagc': 3000,
                 'butter_kwargs': {'N': 3, 'Wn': 0.01, 'btype': 'highpass'}}
     h = neuropixel.trace_header(version=1) if h is None else h
-    ncv = h['x'].size  # number of channels
+    ncv = h['sample_shift'].size  # number of channels
     output_file = sr.file_bin.with_suffix('.bin') if output_file is None else output_file
     assert output_file != sr.file_bin
     taper = np.r_[0, scipy.signal.windows.cosine((SAMPLES_TAPER - 1) * 2), 0]
@@ -504,7 +541,7 @@ def nxcor(x, ref):
     xcor = channels_similarity(raw)
     fscale, psd = scipy.signal.welch(raw * 1e6, fs=fs)  # units; uV ** 2 / Hz
 
-    sos_hp = scipy.signal.butter(**{'N': 3, 'Wn': 1000 / fs / 2, 'btype': 'highpass'}, output='sos')
+    sos_hp = scipy.signal.butter(**{'N': 3, 'Wn': 300 / fs * 2, 'btype': 'highpass'}, output='sos')
     hf = scipy.signal.sosfiltfilt(sos_hp, raw)
     xcorf = channels_similarity(hf)
 
@@ -513,7 +550,7 @@ def nxcor(x, ref):
         'rms_raw': rms(raw),  # very similar to the rms avfter butterworth filter
         'xcor_hf': detrend(xcor, 11),
         'xcor_lf': xcorf - detrend(xcorf, 11) - 1,
-        'psd_hf': np.mean(psd[:, fscale > 12000], axis=-1),
+        'psd_hf': np.mean(psd[:, fscale > (fs / 2 * 0.8)], axis=-1),  # 80% nyquists
     })
 
     # make recommendation

ibllib/dsp/voltage.py

@@ -553,7 +553,7 @@ def groom_pin_state(gpio, audio, ts, tolerance=2., display=False, take='first',
         downs = ts[high2low] - ts[high2low][0]
         offsets = audio_times[1::2] - audio_times[1]
         assigned = attribute_times(offsets, downs, tol=tolerance, take=take)
-        unassigned = np.setdiff1d(np.arange(onsets.size), assigned[assigned > -1])
+        unassigned = np.setdiff1d(np.arange(offsets.size), assigned[assigned > -1])
         if unassigned.size > 0:
             _logger.debug(f'{unassigned.size} audio TTL falls were not detected by the camera')
         # Check that all pin state downticks could be attributed to an offset TTL

ibllib/io/extractors/camera.py

@@ -478,6 +478,11 @@ def _geometry_from_meta(meta_data):
     """
     cm = _map_channels_from_meta(meta_data)
     major_version = _get_neuropixel_major_version_from_meta(meta_data)
+    if cm is None:
+        _logger.warning("Meta data doesn't have geometry (snsShankMap field), returning defaults")
+        th = neuropixel.trace_header(version=major_version)
+        th['flag'] = th['x'] * 0 + 1.
+        return th
     th = cm.copy()
     if major_version == 1:
         # the spike sorting channel maps have a flipped version of the channel map