Merge pull request #1924 from cboulay/nev-reset-per-packet-type

Blackrock: detect NEV clock resets per packet type

Author: h-mayorquin

neo/rawio/blackrockrawio.py

@@ -2038,30 +2038,42 @@ def _get_event_segment_ids(self, raw_event_data, masks, spec):
             # consistency check for monotone increasing time stamps
             # - Use logical comparator (instead of np.diff) to avoid unsigned dtype issues.
             # - Only consider handled/known event types.
-            mask_handled = np.any([mask for mask in masks.values()], axis=0)
-            jump_ids_handled = (
-                np.where(
-                    raw_event_data["timestamp"][mask_handled][1:] < raw_event_data["timestamp"][mask_handled][:-1]
-                )[0]
-                + 1
-            )
-            jump_ids = np.where(mask_handled)[0][jump_ids_handled]  # jump ids in full set of events (incl. unhandled)
-            overlap = np.isin(jump_ids, reset_ev_ids)
-            if not all(overlap):
+            # - Compare only packets of the same type: the file does not order packets of different
+            #   types (e.g. a comment is stamped when the user starts typing but written when they
+            #   finish, so it can be later than the spikes written after it).
+            # jumps maps the position of a packet earlier than the previous packet of its type to the
+            # position of that previous packet, in the full set of events (incl. unhandled).
+            jumps = {}
+            for mask in masks.values():
+                ids = np.flatnonzero(mask)
+                ts = raw_event_data["timestamp"][ids]
+                earlier = np.flatnonzero(ts[1:] < ts[:-1])
+                jumps.update(zip(ids[earlier + 1].tolist(), ids[earlier].tolist()))
+
+            # A reset sends the clock of every packet type back, so each type jumps at its first packet
+            # after it. Count a jump only if the packet it is compared with is in the current segment,
+            # so that one reset is one boundary.
+            reset_set = set(reset_ev_ids.tolist())
+            boundaries = []
+            segment_start = 0
+            for position in sorted(reset_set | jumps.keys()):
+                if position in reset_set or jumps[position] >= segment_start:
+                    boundaries.append(position)
+                    segment_start = position
+            additional_ids = np.array([b for b in boundaries if b not in reset_set])
+            if len(additional_ids):
                 # additional resets occurred without a reset event being stored
-                additional_ids = jump_ids[np.invert(overlap)]
                 warnings.warn(
                     f"Detected {len(additional_ids)} undocumented segments within "
                     f"nev data after timestamps {additional_ids}."
                     ""
                 )
-                reset_ev_ids = sorted(np.unique(np.concatenate((reset_ev_ids, jump_ids))))
 
             event_segment_ids = np.zeros(len(raw_event_data), dtype=int)
-            for reset_event_id in reset_ev_ids:
-                event_segment_ids[reset_event_id:] += 1
+            event_segment_ids[boundaries] = 1
+            event_segment_ids = np.cumsum(event_segment_ids)
 
-            self._nb_segment_nev = len(reset_ev_ids) + 1
+            self._nb_segment_nev = len(boundaries) + 1
             return event_segment_ids
 
         else:

neo/test/rawiotest/test_blackrockrawio.py

@@ -628,6 +628,29 @@ def test_nev_spikes_at_segment_edges_ptp(self):
                     returned += reader.get_spike_timestamps(block_index, seg_index, unit_index, None, None).size
         self.assertEqual(returned, 999)
 
+    def test_nev_segments_ignore_order_between_packet_types(self):
+        """
+        Only a packet earlier than the previous packet of its own type is a clock reset.
+
+        Central writes a comment a few ms after stamping it, so in this file two comments are
+        written before 96 spikes stamped 2.1 ms earlier. Comparing packets of different types
+        split the nev into two segments, and with the ns5 only the 384 spikes of the first one
+        were returned.
+        """
+        filename = self.get_local_path(
+            "blackrock/version_3_0/time_structure/comment_written_before_earlier_spikes/"
+            "comment_written_before_earlier_spikes"
+        )
+        for nsx_to_load in ([], [5]):
+            reader = BlackrockRawIO(filename=filename, nsx_to_load=nsx_to_load)
+            reader.parse_header()
+            self.assertEqual(reader.segment_count(0), 1)
+            returned = sum(
+                reader.get_spike_timestamps(0, 0, unit_index, None, None).size
+                for unit_index in range(reader.spike_channels_count())
+            )
+            self.assertEqual(returned, 960)
+
 
 if __name__ == "__main__":
     unittest.main()