Merge pull request #38 from int-brain-lab/npultra_geom

uhd geometry

src/ibldsp/plots.py

@@ -0,0 +1,40 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal
+
+
+def show_channels_labels(raw, fs, channel_labels, xfeats):
+    """
+    Shows the features side by side a snippet of raw data
+    :param sr:
+    :return:
+    """
+    nc, ns = raw.shape
+    ns_plot = np.minimum(ns, 3000)
+    vaxis_uv = 75
+    sos_hp = scipy.signal.butter(**{"N": 3, "Wn": 300 / fs * 2, "btype": "highpass"}, output="sos")
+    butt = scipy.signal.sosfiltfilt(sos_hp, raw)
+    fig, ax = plt.subplots(1, 5, figsize=(18, 6), gridspec_kw={'width_ratios': [1, 1, 1, 8, .2]})
+    ax[0].plot(xfeats['xcor_hf'], np.arange(nc))
+    ax[0].plot(xfeats['xcor_hf'][iko := channel_labels == 1], np.arange(nc)[iko], 'r*')
+    ax[0].plot([- .5, -.5], [0, nc], 'r--')
+    ax[0].set(ylabel='channel #', xlabel='high coherence', ylim=[0, nc], title='a) dead channel')
+    ax[1].plot(xfeats['psd_hf'], np.arange(nc))
+    ax[1].plot(xfeats['psd_hf'][iko := channel_labels == 2], np.arange(nc)[iko], 'r*')
+    ax[1].plot([.02, .02], [0, nc], 'r--')
+
+    ax[1].set(yticklabels=[], xlabel='PSD', ylim=[0, nc], title='b) noisy channel')
+    ax[1].sharey(ax[0])
+    ax[2].plot(xfeats['xcor_lf'], np.arange(nc))
+    ax[2].plot(xfeats['xcor_lf'][iko := channel_labels == 3], np.arange(nc)[iko], 'r*')
+    ax[2].plot([-.75, -.75], [0, nc], 'r--')
+    ax[2].set(yticklabels=[], xlabel='low coherence', ylim=[0, nc], title='c) outside')
+    ax[2].sharey(ax[0])
+    im = ax[3].imshow(butt[:, :ns_plot] * 1e6, origin='lower', cmap='PuOr', aspect='auto',
+                      vmin=-vaxis_uv, vmax=vaxis_uv, extent=[0, ns_plot / fs * 1e3, 0, nc])
+    ax[3].set(yticklabels=[], title='d) Raw data', xlabel='time (ms)', ylim=[0, nc])
+    ax[3].grid(False)
+    ax[3].sharey(ax[0])
+    plt.colorbar(im, cax=ax[4], shrink=0.8).ax.set(ylabel='(uV)')
+    fig.tight_layout()
+    return fig, ax

src/ibldsp/voltage.py

@@ -15,6 +15,7 @@
 
 import ibldsp.fourier as fourier
 import ibldsp.utils as utils
+import ibldsp.plots
 
 
 def agc(x, wl=0.5, si=0.002, epsilon=1e-8, gpu=False):
@@ -639,7 +640,7 @@ def my_function(i_chunk, n_chunk):
         np.save(output_qc_path.joinpath("_iblqc_ephysSaturation.samples.npy"), saturation_data)
 
 
-def detect_bad_channels(raw, fs, similarity_threshold=(-0.5, 1), psd_hf_threshold=None):
+def detect_bad_channels(raw, fs, similarity_threshold=(-0.5, 1), psd_hf_threshold=None, display=False):
     """
     Bad channels detection for Neuropixel probes
     Labels channels
@@ -651,6 +652,7 @@ def detect_bad_channels(raw, fs, similarity_threshold=(-0.5, 1), psd_hf_threshol
     :param fs: sampling frequency
     :param similarity_threshold:
     :param psd_hf_threshold:
+    :param display: optinal (False) will show a plot of features alongside a raw data snippet
     :return: labels (numpy vector [nc]), xfeats: dictionary of features [nc]
     """
 
@@ -751,6 +753,8 @@ def nxcor(x, ref):
     ichannels[inoisy] = 2
     # from ibllib.plots.figures import ephys_bad_channels
     # ephys_bad_channels(x, 30000, ichannels, xfeats)
+    if display:
+        ibldsp.plots.show_channels_labels(raw, fs, ichannels, xfeats)
     return ichannels, xfeats
 
 
@@ -783,7 +787,6 @@ def detect_bad_channels_cbin(bin_file, n_batches=10, batch_duration=0.3, display
     if display:
         raw = sr[sl, :nc].TO
         from ibllib.plots.figures import ephys_bad_channels
-
         ephys_bad_channels(raw, sr.fs, channel_flags, xfeats_med)
     return channel_flags
 

src/neuropixel.py

@@ -4,6 +4,7 @@
 from typing import Any
 import warnings
 import traceback
+import numbers
 
 import scipy.signal
 import numpy as np
@@ -60,7 +61,11 @@
 NC = 384
 SITES_COORDINATES: np.array
 # channel layouts for neuropixel probes as a function of the major version (1 or 2)
-CHANNEL_GRID = {1: dict(DX=16, X0=11, DY=20, Y0=20), 2: dict(DX=32, X0=27, DY=15, Y0=20)}
+CHANNEL_GRID = {
+    1: dict(DX=16, X0=11, DY=20, Y0=20),
+    2: dict(DX=32, X0=27, DY=15, Y0=20),
+    "NPultra": dict(DX=6, X0=0, DY=6, Y0=0)
+}
 
 
 def _deprecated_sites_coordinates() -> np.array:
@@ -97,7 +102,8 @@ def xy2rc(x, y, version=1):
     :param version: neuropixel major version 1 or 2
     :return: dictionary with keys x and y
     """
-    grid = CHANNEL_GRID[np.floor(version)]
+    version = np.floor(version) if isinstance(version, numbers.Number) else version
+    grid = CHANNEL_GRID[version]
     col = (x - grid['X0']) / grid['DX']
     row = (y - grid['Y0']) / grid['DY']
     return {"col": col, "row": row}
@@ -111,7 +117,8 @@ def rc2xy(row, col, version=1):
     :param version: neuropixel major version 1 or 2
     :return: dictionary with keys x and y
     """
-    grid = CHANNEL_GRID[np.floor(version)]
+    version = np.floor(version) if isinstance(version, numbers.Number) else version
+    grid = CHANNEL_GRID[version]
     x = col * grid['DX'] + grid['X0']
     y = row * grid['DY'] + grid['Y0']
     return {"x": x, "y": y}
@@ -131,6 +138,9 @@ def dense_layout(version=1, nshank=1):
 
     if version == 1:  # version 1 has a dense layout, checkerboard pattern
         ch.update({"col": np.tile(np.array([2, 0, 3, 1]), int(NC / 4))})
+    elif version == "NPultra":  # NPultra has 8 columns with square grid spacing
+        ch.update({"row": np.floor(np.arange(NC) / 8)})
+        ch.update({"col": np.tile(np.arange(8), int(NC / 8))})
     elif (
         np.floor(version) == 2 and nshank == 1
     ):  # single shank NP1 has 2 columns in a dense patter
@@ -189,7 +199,7 @@ def adc_shifts(version=1, nc=NC):
     :param version: neuropixel major version 1 or 2
     :param nc: number of channels
     """
-    if version == 1:
+    if version == 1 or version == "NPultra":
         adc_channels = 12
         n_cycles = 13
         # version 1 uses 32 ADC that sample 12 channels each

src/spikeglx.py

@@ -518,7 +518,7 @@ def _get_serial_number_from_meta(md):
 
 
 def _get_neuropixel_major_version_from_meta(md):
-    MAJOR_VERSION = {"3A": 1, "3B2": 1, "3B1": 1, "NP2.1": 2, "NP2.4": 2.4}
+    MAJOR_VERSION = {"3A": 1, "3B2": 1, "3B1": 1, "NP2.1": 2, "NP2.4": 2.4, "NPultra": "NPultra"}
     version = _get_neuropixel_version_from_meta(md)
     if version is not None:
         return MAJOR_VERSION[version]
@@ -563,6 +563,8 @@ def _get_neuropixel_version_from_meta(md):
     # Neuropixel 2.0 four shank
     elif prb_type == 24 or prb_type == 2013:
         return "NP2.4"
+    elif prb_type == 1100:
+        return "NPultra"
 
 
 def _get_sync_trace_indices_from_meta(md):

src/tests/unit/cpu/test_neuropixel.py

@@ -15,3 +15,73 @@ def test_adc_shifts():
     h24 = neuropixel.trace_header(version=2.4)
     np.testing.assert_equal(h24["sample_shift"], h21["sample_shift"])
     np.testing.assert_equal(np.unique(h21["sample_shift"] * 16), np.arange(16))
+    # test ADC shifts uhd
+    hUHD = neuropixel.trace_header(version="NPultra")
+    np.testing.assert_equal(hUHD["sample_shift"], h1["sample_shift"])
+
+
+def test_geom_np1():
+    gt = dict(
+        ind=np.arange(384),
+        shank=np.zeros(384),
+        row=np.repeat(np.arange(192), 2),
+        col=np.tile(np.array([2, 0, 3, 1]), 96),
+        x=np.tile(np.array([43, 11, 59, 27]), 96),
+        y=np.repeat(np.arange(0, 3840, 20), 2) + 20
+    )
+
+    h = neuropixel.trace_header(1)
+    for k, v in gt.items():
+        np.testing.assert_equal(v, h[k])
+
+
+def test_geom_np2_1shank():
+    gt = dict(
+        ind=np.arange(384),
+        shank=np.zeros(384),
+        row=np.repeat(np.arange(192), 2),
+        col=np.tile(np.array([0, 1]), 192),
+        x=np.tile(np.array([27, 59]), 192),
+        y=np.repeat(np.arange(0, 2880, 15), 2) + 20
+    )
+
+    h = neuropixel.trace_header(2, 1)
+    for k, v in gt.items():
+        np.testing.assert_equal(v, h[k])
+
+
+def test_geom_np2_4shank():
+    depth_blocks = np.vstack(
+        [
+            np.repeat(np.arange(24), 2),
+            np.repeat(np.arange(24, 48), 2)
+        ]
+    )
+    row_ind = np.concatenate([depth_blocks[i] for i in [0, 0, 1, 1, 0, 0, 1, 1]])
+    gt = dict(
+        ind=np.arange(384),
+        shank=np.repeat(np.array([0, 1, 0, 1, 2, 3, 2, 3]), 48),
+        row=row_ind,
+        col=np.tile(np.array([0, 1]), 192),
+        x=np.tile(np.array([27, 59]), 192),
+        y=row_ind * 15 + 20
+    )
+
+    h = neuropixel.trace_header(2, 4)
+    for k, v in gt.items():
+        np.testing.assert_equal(v, h[k])
+
+
+def test_geom_npultra():
+    gt = dict(
+        ind=np.arange(384),
+        shank=np.zeros(384),
+        row=np.repeat(np.arange(48), 8),
+        col=np.tile(np.arange(8), 48),
+        x=np.tile(np.arange(0, 48, 6), 48),
+        y=np.repeat(np.arange(0, 288, 6), 8)
+    )
+
+    h = neuropixel.trace_header("NPultra")
+    for k, v in gt.items():
+        np.testing.assert_equal(v, h[k])

src/tests/unit/cpu/test_spikeglx.py

@@ -404,6 +404,17 @@ def test_read_NP24(self):
             )
             self.assert_read_glx(bin_3b)
 
+    def test_read_NPultra(self):
+        with tempfile.TemporaryDirectory(prefix="glx_test") as tdir:
+            bin_3b = spikeglx._mock_spikeglx_file(
+                Path(tdir).joinpath("sampleNPultra_g0_t0.imec0.ap.bin"),
+                self.workdir / "sampleNPultra_g0_t0.imec0.ap.meta",
+                ns=32,
+                nc=385,
+                sync_depth=16,
+            )
+            self.assert_read_glx(bin_3b)
+
     def test_check_ephys_file(self):
         self.tdir = tempfile.TemporaryDirectory(prefix="glx_test")
         self.addCleanup(self.tdir.cleanup)
@@ -517,6 +528,12 @@ def testGetRevisionAndType(self):
                 self.assertTrue(len(md.keys()) >= 37)
 
                 if meta_data_file.name.split(".")[-2] in ["lf", "ap"]:
+                    # NPultra: non-numerical
+                    if "NPultra" in meta_data_file.name:
+                        self.assertEqual("NPultra", spikeglx._get_neuropixel_version_from_meta(md))
+                        self.assertEqual("NPultra", spikeglx._get_neuropixel_major_version_from_meta(md))
+                        continue
+
                     # for ap and lf look for version number
                     # test getting revision
                     revision = meta_data_file.name[6:8]