Merge branch 'main' into extract_wf

release_notes.md

@@ -1,3 +1,9 @@
+# 0.10.0
+## 0.10.1 2024-03-19
+-  ensure compatibility with spikeglx 202309 metadata coordinates
+## 0.10.0 2024-03-14
+-  add support for online spikeglx reader
+
 # 0.9.0
 ## 0.9.2 2024-02-08
 -   `neurodsp` is now `ibldsp`. Drop-in replacement of the package name is all that is required to update. The `neurodsp` name will disappear on 01-Sep-2024; until then both names will work.

setup.py

@@ -3,12 +3,12 @@
 with open("README.md", "r", encoding="utf-8") as fh:
     long_description = fh.read()
 
-with open('requirements.txt') as f:
-    require = [x.strip() for x in f.readlines() if not x.startswith('git+')]
+with open("requirements.txt") as f:
+    require = [x.strip() for x in f.readlines() if not x.startswith("git+")]
 
 setuptools.setup(
     name="ibl-neuropixel",
-    version="0.9.2",
+    version="0.10.1",
     author="The International Brain Laboratory",
     description="Collection of tools for Neuropixel 1.0 and 2.0 probes data",
     long_description=long_description,
@@ -23,9 +23,9 @@
         "Operating System :: OS Independent",
     ],
     install_requires=require,
-    package_dir={'': 'src'},
-    packages=setuptools.find_packages(where="src", exclude=['tests']),
+    package_dir={"": "src"},
+    packages=setuptools.find_packages(where="src", exclude=["tests"]),
     include_package_data=True,
-    py_modules=['spikeglx', 'neuropixel'],
+    py_modules=["spikeglx", "neuropixel"],
     python_requires=">=3.8",
 )

src/ibldsp/cadzow.py

@@ -99,7 +99,17 @@ def denoise(WAV, x, y, r, imax=None, niter=1):
     return WAV_
 
 
-def cadzow_np1(wav, fs=30000, rank=5, niter=1, fmax=7500, h=None, ovx=int(16), nswx=int(32), npad=int(0)):
+def cadzow_np1(
+    wav,
+    fs=30000,
+    rank=5,
+    niter=1,
+    fmax=7500,
+    h=None,
+    ovx=int(16),
+    nswx=int(32),
+    npad=int(0),
+):
     """
     Apply Fxy rank-denoiser to a full recording of Neuropixel 1 probe geometry
     ntr - nswx has to be a multiple of (nswx - ovx)
@@ -124,11 +134,17 @@ def cadzow_np1(wav, fs=30000, rank=5, niter=1, fmax=7500, h=None, ovx=int(16), n
     imax = np.searchsorted(fscale, fmax)
     WAV = scipy.fft.rfft(wav[:, :])
     padgain = scipy.signal.windows.hann(npad * 2)[:npad]
-    WAV = np.r_[np.flipud(WAV[1:npad + 1, :]) * padgain[:, np.newaxis],
-                WAV,
-                np.flipud(WAV[-npad - 2: - 1, :]) * np.flipud(np.r_[padgain, 1])[:, np.newaxis]]  # apply padding
-    x = np.r_[np.flipud(h['x'][1:npad + 1]), h['x'], np.flipud(h['x'][-npad - 2: - 1])]
-    y = np.r_[np.flipud(h['y'][1:npad + 1]) - 120, h['y'], np.flipud(h['y'][-npad - 2: - 1]) + 120]
+    WAV = np.r_[
+        np.flipud(WAV[1: npad + 1, :]) * padgain[:, np.newaxis],
+        WAV,
+        np.flipud(WAV[-npad - 2: -1, :]) * np.flipud(np.r_[padgain, 1])[:, np.newaxis],
+    ]  # apply padding
+    x = np.r_[np.flipud(h["x"][1: npad + 1]), h["x"], np.flipud(h["x"][-npad - 2: -1])]
+    y = np.r_[
+        np.flipud(h["y"][1: npad + 1]) - 120,
+        h["y"],
+        np.flipud(h["y"][-npad - 2: -1]) + 120,
+    ]
     WAV_ = np.zeros_like(WAV)
     gain = np.zeros(ntr + npad * 2 + 1)
     hanning = scipy.signal.windows.hann(ovx * 2 - 1)[0:ovx]
@@ -144,9 +160,11 @@ def cadzow_np1(wav, fs=30000, rank=5, niter=1, fmax=7500, h=None, ovx=int(16), n
             gw = gain_window
         gain[firstx:lastx] += gw
         array = WAV[firstx:lastx, :]
-        array = denoise(array, x=x[firstx:lastx], y=y[firstx:lastx], r=rank, imax=imax, niter=niter)
+        array = denoise(
+            array, x=x[firstx:lastx], y=y[firstx:lastx], r=rank, imax=imax, niter=niter
+        )
         WAV_[firstx:lastx, :] += array * gw[:, np.newaxis]
 
-    WAV_ = WAV_[npad:-npad - 1]  # remove padding
+    WAV_ = WAV_[npad: -npad - 1]  # remove padding
     wav_ = scipy.fft.irfft(WAV_)
     return wav_

src/ibldsp/cuda_tools.py

@@ -4,7 +4,7 @@
 from iblutil.util import Bunch
 
 
-REGEX_PATTERN = r'const int\s+\S+\s+=\s+\S+.+'
+REGEX_PATTERN = r"const int\s+\S+\s+=\s+\S+.+"
 
 
 def get_cuda(fn, **constants):
@@ -15,10 +15,10 @@ def get_cuda(fn, **constants):
     :return: code: String
              constants: Bunch
     """
-    path = Path(__file__).parent / (fn + '.cu')
+    path = Path(__file__).parent / (fn + ".cu")
     assert path.exists
     code = path.read_text()
-    code = code.replace('__global__ void', 'extern "C" __global__ void')
+    code = code.replace("__global__ void", 'extern "C" __global__ void')
     if not constants:
         return code, Bunch(extract_constants_from_cuda(code))
     return change_cuda_constants(code, constants)
@@ -33,9 +33,9 @@ def extract_constants_from_cuda(code):
     r = re.compile(REGEX_PATTERN)
     m = r.search(code)
     if m:
-        constants = m.group(0).replace('const int', '').replace(';', '').split(',')
+        constants = m.group(0).replace("const int", "").replace(";", "").split(",")
         for const in constants:
-            a, b = const.strip().split('=')
+            a, b = const.strip().split("=")
             yield a.strip(), int(b.strip())
 
 
@@ -49,14 +49,16 @@ def change_cuda_constants(code, constants):
     """
     r = re.compile(REGEX_PATTERN)
     m = r.match(code)
-    assert m, 'No constants found in CUDA code'
+    assert m, "No constants found in CUDA code"
     pattern_length = m.span(0)[1] - 1
-    default_constants_string = m.group(0).replace('const int', '').replace(';', '').split(',')
+    default_constants_string = (
+        m.group(0).replace("const int", "").replace(";", "").split(",")
+    )
     code_constants = {}
 
     # Find default constants in CUDA code
     for default_constants_string in default_constants_string:
-        name, value = default_constants_string.split('=')
+        name, value = default_constants_string.split("=")
         code_constants[name.strip()] = int(value.strip())
 
     # Replace default constants with the new user constants
@@ -65,9 +67,9 @@ def change_cuda_constants(code, constants):
 
     new_strings = []
     for name, value in code_constants.items():
-        new_strings.append(f'{name} = {value}')
-    new_constants_string = ', '.join(new_strings)
+        new_strings.append(f"{name} = {value}")
+    new_constants_string = ", ".join(new_strings)
 
-    new_code = f'const int  {new_constants_string}{code[pattern_length:]}'
+    new_code = f"const int  {new_constants_string}{code[pattern_length:]}"
 
     return new_code, Bunch(code_constants)

src/ibldsp/destripe_gpu.py

@@ -7,8 +7,16 @@
 from .voltage import _get_destripe_parameters, interpolate_bad_channels, kfilt
 
 
-def destripe_array(data, fs=30000, fshigh=300., taper_size=64, sample_shifts=None, channel_labels=None,
-                   channel_xcoords=None, channel_ycoords=None):
+def destripe_array(
+    data,
+    fs=30000,
+    fshigh=300.0,
+    taper_size=64,
+    sample_shifts=None,
+    channel_labels=None,
+    channel_xcoords=None,
+    channel_ycoords=None,
+):
     """
     Applies de-striping to a cupy array
     :param data: float32 cupy array, shape (n_channels, n_times)
@@ -39,17 +47,21 @@ def destripe_array(data, fs=30000, fshigh=300., taper_size=64, sample_shifts=Non
 
     # align channels if the time shifts are provided
     if sample_shifts is not None:
-        sample_shifts = cp.array(sample_shifts, dtype='float32')
+        sample_shifts = cp.array(sample_shifts, dtype="float32")
         data = channel_shift(data, sample_shifts)
 
     # apply spatial filter
     # TODO: Currently using default settings, should allow user to change this in function arguments
     kfilt_kwargs = _get_destripe_parameters(fs, None, None, True)[1]
 
     if channel_labels is not None:
-        data = interpolate_bad_channels(data, channel_labels, channel_xcoords, channel_ycoords, gpu=True)
+        data = interpolate_bad_channels(
+            data, channel_labels, channel_xcoords, channel_ycoords, gpu=True
+        )
         inside_brain = cp.where(channel_labels != 3)[0]
-        data[inside_brain, :] = kfilt(data[inside_brain, :], gpu=True, **kfilt_kwargs)  # apply the k-filter / CAR
+        data[inside_brain, :] = kfilt(
+            data[inside_brain, :], gpu=True, **kfilt_kwargs
+        )  # apply the k-filter / CAR
     else:
         data = kfilt(data, gpu=True, **kfilt_kwargs)  # apply the k-filter / CAR
 
@@ -65,6 +77,6 @@ def get_sos(fs, fshigh, fslow=None):
     :return: sos, second-order sections
     """
     if fslow and fslow < fs / 2:
-        return butter(3, (2 * fshigh / fs, 2 * fslow / fs), 'bandpass', output='sos')
+        return butter(3, (2 * fshigh / fs, 2 * fslow / fs), "bandpass", output="sos")
     else:
-        return butter(3, 2 * fshigh / fs, 'high', output='sos')
+        return butter(3, 2 * fshigh / fs, "high", output="sos")

src/ibldsp/filter_gpu.py

@@ -22,7 +22,7 @@ def sosfiltfilt_gpu(sos, x, axis=-1):
 
     n_sections, m = sos.shape
     if m != 6:
-        raise ValueError('sos array must be shape (n_sections, 6)')
+        raise ValueError("sos array must be shape (n_sections, 6)")
 
     ntaps = 2 * n_sections + 1
     ntaps -= min((sos[:, 2] == 0).sum(), (sos[:, 5] == 0).sum())
@@ -64,20 +64,20 @@ def sosfilt_gpu(sos, x, axis, zi):
 
     n_sections, m = sos.shape
     if m != 6:
-        raise ValueError('sos array must be shape (n_sections, 6)')
+        raise ValueError("sos array must be shape (n_sections, 6)")
 
     x_zi_shape = list(x.shape)
     x_zi_shape[axis] = 2
     x_zi_shape = tuple([n_sections] + x_zi_shape)
 
     if zi is not None:
-        assert zi.shape == x_zi_shape, f'zi has shape {zi.shape}, expected {x_zi_shape}'
-        zi = cp.array(zi, dtype='float32')
+        assert zi.shape == x_zi_shape, f"zi has shape {zi.shape}, expected {x_zi_shape}"
+        zi = cp.array(zi, dtype="float32")
     else:
-        zi = cp.zeros(x_zi_shape, dtype='float32')
+        zi = cp.zeros(x_zi_shape, dtype="float32")
 
-    sos = cp.array(sos, dtype='float32')
-    assert x.dtype == 'float32', f'Expected float32 data, got {x.dtype}'
+    sos = cp.array(sos, dtype="float32")
+    assert x.dtype == "float32", f"Expected float32 data, got {x.dtype}"
 
     axis = axis % x.ndim
     x = cp.ascontiguousarray(cp.moveaxis(x, axis, -1))
@@ -91,15 +91,15 @@ def sosfilt_gpu(sos, x, axis, zi):
 
 
 def _cuda_sosfilt(sos, x, zi):
-
     n_signals, n_samples = x.shape
     n_sections = sos.shape[0]
 
     n_blocks, n_threads = sosfilt_kernel_params(n_signals)
 
-    code, consts = get_cuda('sosfilt', n_signals=n_signals, n_samples=n_samples,
-                            n_sections=n_sections)
-    kernel = cp.RawKernel(code, 'sosfilt')
+    code, consts = get_cuda(
+        "sosfilt", n_signals=n_signals, n_samples=n_samples, n_sections=n_sections
+    )
+    kernel = cp.RawKernel(code, "sosfilt")
     kernel((n_blocks,), (n_threads,), (sos, x, zi))
 
 
@@ -152,17 +152,20 @@ def odd_ext(x, n, axis=-1):
     if n < 1:
         return x
     if n > x.shape[axis] - 1:
-        raise ValueError(("The extension length n (%d) is too big. " +
-                          "It must not exceed x.shape[axis]-1, which is %d.")
-                         % (n, x.shape[axis] - 1))
+        raise ValueError(
+            (
+                "The extension length n (%d) is too big. "
+                + "It must not exceed x.shape[axis]-1, which is %d."
+            )
+            % (n, x.shape[axis] - 1)
+        )
     left_end = axis_slice(x, start=0, stop=1, axis=axis)
     left_ext = axis_slice(x, start=n, stop=0, step=-1, axis=axis)
     right_end = axis_slice(x, start=-1, axis=axis)
     right_ext = axis_slice(x, start=-2, stop=-(n + 2), step=-1, axis=axis)
-    ext = cp.concatenate((2 * left_end - left_ext,
-                          x,
-                          2 * right_end - right_ext),
-                         axis=axis)
+    ext = cp.concatenate(
+        (2 * left_end - left_ext, x, 2 * right_end - right_ext), axis=axis
+    )
     return ext