Merge pull request #1862 from h-mayorquin/add_intan_test

Add a real regression test for one-file-per-signal Intan digital streams

Author: alejoe91

neo/rawio/intanrawio.py

@@ -571,13 +571,15 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
 
         """
         dtype = np.uint16  # We fix this to match the memmap dtype
+        # Slice to the requested window once up front so the bitwise unpacking below only runs over
+        # the requested samples rather than the whole recording for every channel.
+        raw_digital_data = raw_digital_data[i_start:i_stop]
         output = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
 
         for channel_index, channel_id in enumerate(channel_ids):
             native_order = self.native_channel_order[channel_id]
             mask = 1 << native_order
-            demultiplex_data = np.bitwise_and(raw_digital_data, mask) > 0
-            output[:, channel_index] = demultiplex_data[i_start:i_stop].flatten()
+            output[:, channel_index] = (np.bitwise_and(raw_digital_data, mask) > 0).flatten()
 
         return output
 

neo/test/rawiotest/test_intanrawio.py

@@ -18,6 +18,7 @@ class TestIntanRawIO(
         "intan/rhs_fpc_multistim_240514_082243/rhs_fpc_multistim_240514_082243.rhs",  # Format header-attached newer version
         "intan/intan_fpc_test_231117_052630/info.rhd",  # Format one-file-per-channel
         "intan/intan_fps_test_231117_052500/info.rhd",  # Format one file per signal
+        "intan/intan_fps_multiple_digital_channels/info.rhd",  # one-file-per-signal with multiple packed digital channels (issue #1853)
         "intan/intan_fpc_rhs_test_240329_091637/info.rhs",  # Format one-file-per-channel
         "intan/intan_fps_rhs_test_240329_091536/info.rhs",  # Format one-file-per-signal
         "intan/rhd_fpc_multistim_240514_082044/info.rhd",  # Multiple digital channels one-file-per-channel rhd
@@ -218,25 +219,58 @@ def test_correct_decoding_of_stimulus_current(self):
 
         assert np.isclose(duration_of_positive_pulse, expected_duration)
 
+    def test_reading_one_file_per_signal_multiple_digital_channels(self):
+        "Regression test for https://github.com/NeuralEnsemble/python-neo/issues/1853"
+        # One-file-per-signal recording with 16 digital-input and 16 digital-output channels.
+        # All channels of a digital stream are packed as bit positions of a single 16-bit word
+        # per sample, so the per-file sample count must divide by one word, not by the channel
+        # count. Earlier fixtures had a single digital channel per stream, so the division by one
+        # happened to be correct and never exercised this path.
+        file_path = Path(self.get_local_path("intan/intan_fps_multiple_digital_channels/info.rhd"))
+        intan_reader = IntanRawIO(filename=file_path)
+        intan_reader.parse_header()
 
-class TestIntanDigitalDemultiplexShape(unittest.TestCase):
-    """Regression coverage for https://github.com/NeuralEnsemble/python-neo/issues/1853."""
-
-    def test_demultiplex_handles_packed_digital_buffer(self):
-        # Digital streams pack all channels into one uint16 word per timestamp.
-        n_samples = 8
-        packed = np.array([0, 1, 16, 17, 0, 16, 1, 17], dtype=np.uint16).reshape(n_samples, 1)
-        expected_bit0 = np.array([0, 1, 0, 1, 0, 0, 1, 1], dtype=np.uint16)
-        expected_bit4 = np.array([0, 0, 1, 1, 0, 1, 0, 1], dtype=np.uint16)
-
-        reader = IntanRawIO.__new__(IntanRawIO)
-        reader.native_channel_order = {"DIN-00": 0, "DIN-04": 4}
+        signal_streams = intan_reader.header["signal_streams"]
+        signal_channels = intan_reader.header["signal_channels"]
+        stream_names = signal_streams["name"].tolist()
+        stream_ids = signal_streams["id"].tolist()
 
-        chunk = reader._demultiplex_digital_data(packed, ["DIN-00", "DIN-04"], 0, n_samples)
+        amplifier_stream_index = stream_names.index("RHD2000 amplifier channel")
+        expected_num_samples = intan_reader.get_signal_size(
+            block_index=0, seg_index=0, stream_index=amplifier_stream_index
+        )
 
-        self.assertEqual(chunk.shape, (n_samples, 2))
-        np.testing.assert_array_equal(chunk[:, 0], expected_bit0)
-        np.testing.assert_array_equal(chunk[:, 1], expected_bit4)
+        folder_path = file_path.parent
+        digital_stream_to_raw_file = {
+            "USB board digital input channel": folder_path / "digitalin.dat",
+            "USB board digital output channel": folder_path / "digitalout.dat",
+        }
+        for stream_name, raw_file_path in digital_stream_to_raw_file.items():
+            stream_index = stream_names.index(stream_name)
+            stream_id = stream_ids[stream_index]
+            channel_ids = [channel["name"] for channel in signal_channels if channel["stream_id"] == stream_id]
+            assert len(channel_ids) == 16
+
+            # The packed word holds all channels, so the stream reports the full sample count.
+            # Before the fix it was divided by the channel count and get_analogsignal_chunk crashed.
+            num_samples = intan_reader.get_signal_size(block_index=0, seg_index=0, stream_index=stream_index)
+            assert num_samples == expected_num_samples
+
+            chunk = intan_reader.get_analogsignal_chunk(stream_index=stream_index, channel_ids=channel_ids)
+            assert chunk.shape == (expected_num_samples, len(channel_ids))
+
+            # Each channel is one bit of the packed word; demultiplex the raw file and compare.
+            # Intan RHD Application Note: Data File Formats (digitalin.dat / digitalout.dat):
+            # "All 16 digital inputs are encoded bit-by-bit in each 16-bit word. For example, if
+            # digital inputs 0, 4, and 5 are high and the rest low, the uint16 value for this sample
+            # time will be 2^0 + 2^4 + 2^5 = 1 + 16 + 32 = 49." The note isolates a channel with the
+            # MATLAB recipe "digital_input_ch = (bitand(digital_word, 2^ch) > 0)", which is what the
+            # bitwise_and below reproduces (2^ch -> 1 << bit_position).
+            packed_words = np.fromfile(raw_file_path, dtype=np.uint16)
+            for channel_index, channel_id in enumerate(channel_ids):
+                bit_position = int(channel_id.split("-")[-1])
+                expected = (np.bitwise_and(packed_words, np.uint16(1 << bit_position)) > 0).astype(np.uint16)
+                np.testing.assert_array_equal(chunk[:, channel_index], expected)
 
 
 if __name__ == "__main__":