Add NIDQ and LF stream support to SpikeGLX reader

.github/workflows/test-symmetry.yml

@@ -26,15 +26,24 @@ jobs:
       # Check out NDR-matlab at the same branch name as this workflow run
       # (e.g. the PR's head branch or the branch being pushed to main).
       # Developers are expected to keep matching branch names across the
-      # two repos when making cross-language changes, the same way the
-      # symmetry branches work. Falls back to main via the default so
-      # `main` pushes also work.
+      # two repos when making cross-language changes. Falls back to main
+      # when the matching branch does not exist in NDR-matlab.
+      - name: Resolve NDR-matlab branch
+        id: matlab-ref
+        run: |
+          BRANCH="${{ github.head_ref || github.ref_name }}"
+          if git ls-remote --exit-code --heads https://github.com/VH-Lab/NDR-matlab.git "$BRANCH" >/dev/null 2>&1; then
+            echo "ref=$BRANCH" >> "$GITHUB_OUTPUT"
+          else
+            echo "ref=main" >> "$GITHUB_OUTPUT"
+          fi
+
       - name: Check out NDR-matlab (matching branch)
         uses: actions/checkout@v4
         with:
           repository: VH-Lab/NDR-matlab
           path: NDR-matlab
-          ref: ${{ github.head_ref || github.ref_name }}
+          ref: ${{ steps.matlab-ref.outputs.ref }}
 
       # -- Runtime setup ----------------------------------------------------
       - name: Start virtual display server

src/ndr/format/neuropixelsGLX/header.py

@@ -8,6 +8,8 @@
 from pathlib import Path
 from typing import Any
 
+import numpy as np
+
 from ndr.format.neuropixelsGLX.readmeta import readmeta
 
 
@@ -49,8 +51,13 @@ def header(metafilename: str | Path) -> dict[str, Any]:
         Dictionary with fields:
         - sample_rate : float
         - n_saved_chans : int
-        - n_neural_chans : int
-        - n_sync_chans : int
+        - n_neural_chans : int — for NIDQ this is MN + MA + XA
+        - n_sync_chans : int — number of digital word int16 columns
+        - n_digital_word_cols : int — same as n_sync_chans
+        - n_digital_lines : int — number of individual bit lines exposed
+        - digital_line_col : ndarray — 0-based DW column offset per line
+        - digital_line_bit : ndarray — 0-based bit position per line
+        - digital_line_label : list of str — label per line
         - saved_chan_list : list of int (1-based)
         - voltage_range : tuple of (float, float)
         - max_int : int
@@ -79,7 +86,13 @@ def header(metafilename: str | Path) -> dict[str, Any]:
     # Number of saved channels
     info["n_saved_chans"] = int(meta["nSavedChans"])
 
-    # Parse snsApLfSy or snsMnMaXaDw to determine neural vs sync channels
+    # Parse snsApLfSy or snsMnMaXaDw to determine neural vs sync channels.
+    # Also compute digital line mapping:
+    #   n_digital_word_cols : number of int16 columns holding digital data
+    #   n_digital_lines     : number of single-bit digital lines exposed
+    #   digital_line_col    : (n_digital_lines,) 0-based DW column offset
+    #   digital_line_bit    : (n_digital_lines,) 0-based bit position
+    #   digital_line_label  : list of str labels per line
     if "snsApLfSy" in meta:
         counts = [int(x) for x in meta["snsApLfSy"].split(",")]
         # counts[0] = AP chans, counts[1] = LF chans, counts[2] = SY chans
@@ -92,15 +105,76 @@ def header(metafilename: str | Path) -> dict[str, Any]:
         else:
             info["stream_type"] = "ap"
             info["n_neural_chans"] = counts[0]
+        # IMEC sync: each sync column provides 16 bits
+        info["n_digital_word_cols"] = info["n_sync_chans"]
+        n_lines = 16 * info["n_sync_chans"]
+        info["n_digital_lines"] = n_lines
+        cols = np.zeros(n_lines, dtype=int)
+        bits = np.zeros(n_lines, dtype=int)
+        labels: list[str] = []
+        idx = 0
+        for c in range(info["n_sync_chans"]):
+            for b in range(16):
+                cols[idx] = c
+                bits[idx] = b
+                labels.append(f"SY{c}.{b}")
+                idx += 1
+        info["digital_line_col"] = cols
+        info["digital_line_bit"] = bits
+        info["digital_line_label"] = labels
     elif "snsMnMaXaDw" in meta:
         info["stream_type"] = "nidq"
         counts = [int(x) for x in meta["snsMnMaXaDw"].split(",")]
-        info["n_neural_chans"] = counts[0]
-        info["n_sync_chans"] = 0
+        info["n_mn_chans"] = counts[0]  # multiplexed neural
+        info["n_ma_chans"] = counts[1]  # multiplexed analog
+        info["n_xa_chans"] = counts[2]  # non-multiplexed analog
+        info["n_dw_chans"] = counts[3]  # digital word int16 columns
+        info["n_neural_chans"] = counts[0] + counts[1] + counts[2]
+        info["n_sync_chans"] = counts[3]
+        info["n_digital_word_cols"] = counts[3]
+
+        # Bytes saved per port — each byte = 8 active digital lines
+        n_bytes_p0 = int(meta.get("niXDBytes1", "0"))
+        n_bytes_p1 = int(meta.get("niXDBytes2", "0"))
+
+        if n_bytes_p0 == 0 and n_bytes_p1 == 0:
+            # Fall back: assume all 16 bits of every DW column are active
+            n_lines_p0 = 16 * info["n_dw_chans"]
+            n_lines_p1 = 0
+        else:
+            n_lines_p0 = 8 * n_bytes_p0
+            n_lines_p1 = 8 * n_bytes_p1
+
+        n_lines = n_lines_p0 + n_lines_p1
+        info["n_digital_lines"] = n_lines
+        cols = np.zeros(n_lines, dtype=int)
+        bits = np.zeros(n_lines, dtype=int)
+        labels = []
+        idx = 0
+        for k in range(n_lines_p0):
+            abs_bit = k
+            cols[idx] = abs_bit // 16
+            bits[idx] = abs_bit % 16
+            labels.append(f"XD{k}")
+            idx += 1
+        for k in range(n_lines_p1):
+            abs_bit = n_bytes_p0 * 8 + k
+            cols[idx] = abs_bit // 16
+            bits[idx] = abs_bit % 16
+            labels.append(f"XD1.{k}")
+            idx += 1
+        info["digital_line_col"] = cols
+        info["digital_line_bit"] = bits
+        info["digital_line_label"] = labels
     else:
         info["stream_type"] = "unknown"
         info["n_neural_chans"] = info["n_saved_chans"] - 1
         info["n_sync_chans"] = 1
+        info["n_digital_word_cols"] = 1
+        info["n_digital_lines"] = 16
+        info["digital_line_col"] = np.zeros(16, dtype=int)
+        info["digital_line_bit"] = np.arange(16, dtype=int)
+        info["digital_line_label"] = [f"bit{b}" for b in range(16)]
 
     # Parse saved channel subset
     if "snsSaveChanSubset" in meta:
@@ -115,15 +189,27 @@ def header(metafilename: str | Path) -> dict[str, Any]:
         vmax = float(meta["imAiRangeMax"])
         vmin = float(meta["imAiRangeMin"])
         info["voltage_range"] = (vmin, vmax)
+    elif "niAiRangeMax" in meta:
+        vmax = float(meta["niAiRangeMax"])
+        vmin = float(meta["niAiRangeMin"])
+        info["voltage_range"] = (vmin, vmax)
     else:
         info["voltage_range"] = (-0.6, 0.6)  # Neuropixels 1.0 default
 
     # Max integer value
     if "imMaxInt" in meta:
         info["max_int"] = int(meta["imMaxInt"])
+    elif "niMaxInt" in meta:
+        info["max_int"] = int(meta["niMaxInt"])
     else:
         info["max_int"] = 512  # Neuropixels 1.0 default
 
+    # NI-DAQ gains
+    if "niMNGain" in meta:
+        info["ni_mn_gain"] = float(meta["niMNGain"])
+    if "niMAGain" in meta:
+        info["ni_ma_gain"] = float(meta["niMAGain"])
+
     # Bits per sample
     info["bits_per_sample"] = 16
 

src/ndr/reader/neuropixelsGLX.py

@@ -18,21 +18,27 @@
 
 
 class ndr_reader_neuropixelsGLX(ndr_reader_base):
-    """Reader for Neuropixels SpikeGLX AP-band data.
+    """Reader for SpikeGLX data (AP, LF, and NIDQ streams).
 
-    This class reads action-potential band data from Neuropixels probes
-    acquired with the SpikeGLX software. Each instance handles one probe's
-    AP stream (one .ap.bin / .ap.meta file pair per epoch).
+    This class reads data from Neuropixels probes and NI-DAQ devices
+    acquired with the SpikeGLX software. Each instance handles one
+    stream (one .bin / .meta file pair per epoch).
 
-    SpikeGLX saves Neuropixels data as flat interleaved int16 binary files
-    with companion .meta text files. The binary files have no header.
-    Channel count, sample rate, and gain information are read from the
-    .meta file.
+    SpikeGLX saves data as flat interleaved int16 binary files with
+    companion .meta text files. The binary files have no header.
+    Channel count, sample rate, and gain information are read from
+    the .meta file.
 
     Channel mapping:
-        - Neural channels are exposed as 'analog_in' (ai1..aiN)
-        - The sync word is exposed as 'digital_in' (di1)
-        - A single time channel 't1' is always present
+        - Analog channels are exposed as 'analog_in' (ai1..aiN).
+          For AP streams these are neural probe channels; for NIDQ
+          streams these are NI-DAQ analog inputs (MN + MA + XA).
+        - Digital lines are exposed as 'digital_in' (di1..diM),
+          where each di channel is a single bit of the packed digital
+          word(s). For NIDQ streams the count comes from
+          ``8 * (niXDBytes1 + niXDBytes2)``; for IMEC streams it is
+          ``16 * n_sync_chans``.
+        - A single time channel 't1' is always present.
 
     Data is returned as int16 to preserve native precision. Use
     :func:`ndr.format.neuropixelsGLX.samples2volts` for voltage conversion.
@@ -75,8 +81,10 @@ def getchannelsepoch(
     ) -> list[dict[str, Any]]:
         """List channels available for a given epoch.
 
-        Neural channels are 'analog_in' (ai1..aiN), the sync channel is
-        'digital_in' (di1), and a time channel 't1' is always present.
+        Analog channels are 'analog_in' (ai1..aiN), digital lines are
+        'digital_in' (di1..diM) with one entry per single-bit line in
+        the packed digital word(s), and a time channel 't1' is always
+        present.
         """
         metafile = self.filenamefromepochfiles(epochstreams)
         info = header(metafile)
@@ -86,13 +94,13 @@ def getchannelsepoch(
         # Time channel
         channels.append({"name": "t1", "type": "time", "time_channel": 1})
 
-        # Neural channels (analog_in)
+        # Analog channels (analog_in)
         for i in range(1, info["n_neural_chans"] + 1):
             channels.append({"name": f"ai{i}", "type": "analog_in", "time_channel": 1})
 
-        # Sync channel (digital_in)
-        if info["n_sync_chans"] > 0:
-            channels.append({"name": "di1", "type": "digital_in", "time_channel": 1})
+        # Digital lines (digital_in) — one per bit of the packed digital word(s)
+        for i in range(1, info["n_digital_lines"] + 1):
+            channels.append({"name": f"di{i}", "type": "digital_in", "time_channel": 1})
 
         return channels
 
@@ -145,9 +153,12 @@ def readchannels_epochsamples(
 
         Reads data between sample s0 and s1 (inclusive, 1-based).
 
-        For 'analog_in': returns int16 neural data.
+        For 'analog_in': returns int16 data (neural for AP, NI-DAQ analog
+        inputs for NIDQ).
         For 'time': returns float64 time stamps in seconds.
-        For 'digital_in': returns int16 sync word values.
+        For 'digital_in': returns int16 single-bit values (0 or 1)
+        extracted from the packed digital word(s). ``channel`` gives
+        the 1-based digital line(s).
         """
         metafile = self.filenamefromepochfiles(epochstreams)
         info = header(metafile)
@@ -178,18 +189,46 @@ def readchannels_epochsamples(
             return data
 
         elif ct in ("digital_in", "di"):
-            # Sync channel is the last channel in the file
-            sync_chan = np.array([info["n_saved_chans"]], dtype=np.int64)
-            data, _, _, _ = binarymatrix_read(
-                binfile,
-                info["n_saved_chans"],
-                sync_chan,
-                float(s0),
-                float(s1),
-                dataType="int16",
-                byteOrder="ieee-le",
-                headerSkip=0,
-            )
+            if isinstance(channel, int):
+                channel = [channel]
+            line_idx = np.array(channel, dtype=int)
+
+            if np.any(line_idx < 1) or np.any(line_idx > info["n_digital_lines"]):
+                raise ValueError(
+                    f"Digital line out of range; valid lines are 1..{info['n_digital_lines']}."
+                )
+
+            # Digital word columns are at the end of each sample row
+            first_dw_col = info["n_saved_chans"] - info["n_digital_word_cols"] + 1
+
+            # Look up the (column, bit) position for each requested line
+            # line_idx is 1-based; digital_line_col/bit arrays are 0-indexed
+            col_offsets = info["digital_line_col"][line_idx - 1]
+            bit_pos = info["digital_line_bit"][line_idx - 1]
+
+            n_samples = int(s1) - int(s0) + 1
+            data = np.zeros((n_samples, len(channel)), dtype=np.int16)
+
+            unique_cols = np.unique(col_offsets)
+            for uc in unique_cols:
+                file_col = np.array([first_dw_col + uc], dtype=np.int64)
+                raw, _, _, _ = binarymatrix_read(
+                    binfile,
+                    info["n_saved_chans"],
+                    file_col,
+                    float(s0),
+                    float(s1),
+                    dataType="int16",
+                    byteOrder="ieee-le",
+                    headerSkip=0,
+                )
+                # raw is (n_samples, 1) int16 — treat as uint16 for bit extraction
+                raw_uint = raw[:, 0].view(np.uint16)
+                mask = col_offsets == uc
+                indices = np.where(mask)[0]
+                for idx in indices:
+                    data[:, idx] = ((raw_uint >> bit_pos[idx]) & 1).astype(np.int16)
+
             return data
 
         else:
@@ -204,9 +243,9 @@ def readevents_epochsamples_native(
         t0: float,
         t1: float,
     ) -> tuple[np.ndarray, np.ndarray]:
-        """Read events or markers. Neuropixels has no native event channels.
+        """Read events or markers. SpikeGLX has no native event channels.
 
-        Returns empty arrays since Neuropixels data is purely
+        Returns empty arrays since SpikeGLX data is purely
         regularly-sampled.
         """
         return np.array([]), np.array([])
@@ -220,8 +259,8 @@ def samplerate(
     ) -> np.ndarray | float:
         """Get the sample rate for specified channels.
 
-        For Neuropixels AP-band, all channels share the same sample rate
-        (typically 30 kHz).
+        All channels in a single SpikeGLX binary file share one sample
+        rate: typically 30 kHz for AP, 2.5 kHz for LF, or ~25 kHz for NIDQ.
         """
         metafile = self.filenamefromepochfiles(epochstreams)
         info = header(metafile)
@@ -230,18 +269,39 @@ def samplerate(
         return info["sample_rate"]
 
     def filenamefromepochfiles(self, filename_array: list[str]) -> str:
-        """Identify the .ap.meta file from epoch file list.
+        """Identify the companion .meta file from the epoch file list.
 
-        Searches the list for a file matching *.ap.meta. Errors if zero
-        or more than one match is found.
-        """
-        matches = [f for f in filename_array if f.lower().endswith(".ap.meta")]
+        First searches for a ``.bin`` file and derives the ``.meta`` path
+        from it (replacing the last 3 characters). This allows multiple
+        ``.meta`` files (e.g. both AP and NIDQ) to coexist in the same
+        epoch file list — the ``.bin`` file disambiguates which stream
+        this reader instance handles.
 
-        if len(matches) == 0:
-            raise FileNotFoundError("No .ap.meta file found in the epoch file list.")
-        if len(matches) > 1:
-            raise ValueError("Multiple .ap.meta files found. Each epoch should have exactly one.")
-        return matches[0]
+        If no ``.bin`` file is present, falls back to finding a single
+        ``.meta`` file. Errors if zero or more than one ``.meta`` match
+        is found without a ``.bin`` to disambiguate.
+        """
+        # Primary path: find the .bin file and derive .meta from it
+        binfile = ""
+        for f in filename_array:
+            if f.lower().endswith(".bin"):
+                binfile = f
+                break
+
+        if binfile:
+            metafile = binfile[:-3] + "meta"
+            # Verify the .meta exists in the file list or on disk
+            if metafile not in filename_array and not Path(metafile).is_file():
+                raise FileNotFoundError(f"No companion .meta file found for {binfile}.")
+            return metafile
+
+        # Fallback: find a single .meta file
+        meta_matches = [f for f in filename_array if f.lower().endswith(".meta")]
+        if len(meta_matches) == 0:
+            raise FileNotFoundError("No .meta file found in the epoch file list.")
+        if len(meta_matches) > 1:
+            raise ValueError("Multiple .meta files found and no .bin file to disambiguate.")
+        return meta_matches[0]
 
     def daqchannels2internalchannels(
         self,