Feature: support event extraction on the ABF format

neo/rawio/axonrawio.py

@@ -289,22 +289,40 @@ def _parse_header(self):
         signal_buffers = np.array([("Signals", buffer_id)], dtype=_signal_buffer_dtype)
         signal_streams = np.array([("Signals", stream_id, buffer_id)], dtype=_signal_stream_dtype)
 
-        # only one events channel : tag
-        # In ABF timstamps are not attached too any particular segment
-        # so each segment acess all event
-        timestamps = []
-        labels = []
-        comments = []
-        for i, tag in enumerate(info["listTag"]):
-            timestamps.append(tag["lTagTime"])
-            labels.append(str(tag["nTagType"]))
-            comments.append(str(clean_string(tag["sComment"])))
-        self._raw_ev_timestamps = np.array(timestamps)
-        self._ev_labels = np.array(labels, dtype="U")
-        self._ev_comments = np.array(comments, dtype="U")
-        event_channels = [("Tag", "", "event")]
+        # Tags are the only events in an ABF file. Each tag type present in the file gets its own
+        # event channel, and the tag comment is the label. Tags belong to the whole recording, so
+        # each segment gets the tags that fall inside its sweep, and a tag between sweeps goes to
+        # no segment.
+        # lTagTime counts in units of fSynchTimeUnit microseconds, or in samples summed over all
+        # channels when fSynchTimeUnit is 0. Keep the seconds per count to rescale the timestamps.
+        if version < 2.0:
+            fSynchTimeUnit = info["fSynchTimeUnit"]
+        elif version >= 2.0:
+            fSynchTimeUnit = info["protocol"]["fSynchTimeUnit"]
+        if fSynchTimeUnit == 0:
+            self._tag_time_unit = 1.0 / (self._sampling_rate * nbchannel)
+        else:
+            self._tag_time_unit = fSynchTimeUnit * 1e-6
+
+        tags = info["listTag"]
+        tag_times = tags["lTagTime"] * self._tag_time_unit
+        tag_types = [tag_type for tag_type in TAG_TYPE_NAMES if np.any(tags["nTagType"] == tag_type)]
+        event_channels = [(TAG_TYPE_NAMES[tag_type], TAG_TYPE_NAMES[tag_type], "event") for tag_type in tag_types]
         event_channels = np.array(event_channels, dtype=_event_channel_dtype)
 
+        self._raw_ev_timestamps = {}
+        self._ev_labels = {}
+        for seg_index in range(nb_segment):
+            seg_t_start = self._t_starts[seg_index]
+            num_samples = self._buffer_descriptions[0][seg_index][buffer_id]["shape"][0]
+            seg_t_stop = seg_t_start + num_samples / self._sampling_rate
+            in_segment = (tag_times >= seg_t_start) & (tag_times < seg_t_stop)
+            for chan_index, tag_type in enumerate(tag_types):
+                seg_tags = tags[in_segment & (tags["nTagType"] == tag_type)]
+                labels = [clean_string(comment).decode("utf-8", errors="replace") for comment in seg_tags["sComment"]]
+                self._raw_ev_timestamps[seg_index, chan_index] = seg_tags["lTagTime"]
+                self._ev_labels[seg_index, chan_index] = np.array(labels, dtype="U")
+
         # No spikes
         spike_channels = []
         spike_channels = np.array(spike_channels, dtype=_spike_channel_dtype)
@@ -334,10 +352,6 @@ def _parse_header(self):
             nADCNum = np.array([adc_nums[c] for c in range(signal_channels.size)])
             signal_an["__array_annotations__"]["nADCNum"] = nADCNum
 
-            for c in range(event_channels.size):
-                ev_ann = seg_annotations["events"][c]
-                ev_ann["comments"] = self._ev_comments
-
     def _source_name(self):
         return self.filename
 
@@ -356,29 +370,27 @@ def _get_analogsignal_buffer_description(self, block_index, seg_index, buffer_id
         return self._buffer_descriptions[block_index][seg_index][buffer_id]
 
     def _event_count(self, block_index, seg_index, event_channel_index):
-        return self._raw_ev_timestamps.size
+        return self._raw_ev_timestamps[seg_index, event_channel_index].size
 
     def _get_event_timestamps(self, block_index, seg_index, event_channel_index, t_start, t_stop):
-        # In ABF timestamps are not attached too any particular segment
-        # so each segment accesses all events
-        timestamp = self._raw_ev_timestamps
-        labels = self._ev_labels
+        timestamp = self._raw_ev_timestamps[seg_index, event_channel_index]
+        labels = self._ev_labels[seg_index, event_channel_index]
         durations = None
 
         if t_start is not None:
-            keep = timestamp >= int(t_start * self._sampling_rate)
+            keep = timestamp * self._tag_time_unit >= t_start
             timestamp = timestamp[keep]
             labels = labels[keep]
 
         if t_stop is not None:
-            keep = timestamp <= int(t_stop * self._sampling_rate)
+            keep = timestamp * self._tag_time_unit <= t_stop
             timestamp = timestamp[keep]
             labels = labels[keep]
 
         return timestamp, durations, labels
 
     def _rescale_event_timestamp(self, event_timestamps, dtype, event_channel_index):
-        event_times = event_timestamps.astype(dtype) / self._sampling_rate
+        event_times = event_timestamps.astype(dtype) * self._tag_time_unit
         return event_times
 
     def read_raw_protocol(self):
@@ -699,18 +711,9 @@ def _parse_abf_v1(f, header_description):
     header["lFileStartTime"] += header["nFileStartMillisecs"] * 0.001
 
     # tags
-    listTag = []
-    for i in range(header["lNumTagEntries"]):
-        f.seek(header["lTagSectionPtr"] + i * 64)
-        tag = {}
-        for key, fmt in TagInfoDescription:
-            val = f.read_f(fmt)
-            if len(val) == 1:
-                tag[key] = val[0]
-            else:
-                tag[key] = np.array(val)
-        listTag.append(tag)
-    header["listTag"] = listTag
+    f.seek(header["lTagSectionPtr"] * BLOCKSIZE)
+    tag_bytes = f.read(header["lNumTagEntries"] * TagInfoDescription.itemsize)
+    header["listTag"] = np.frombuffer(tag_bytes, dtype=TagInfoDescription)
 
     # protocol name formatting. Decode to str (like the channel names) so consumers get a plain
     # string rather than a bytes value, whose str() would bake the "b'...'" repr into the path.
@@ -1024,19 +1027,9 @@ def _parse_abf_v2(f, header_description):
     header["sProtocolPath"] = strings[header["uProtocolPathIndex"]].decode("utf-8", errors="replace")
 
     # tags
-    listTag = []
-    for i in range(sections["TagSection"]["llNumEntries"]):
-        f.seek(sections["TagSection"]["uBlockIndex"] * BLOCKSIZE + sections["TagSection"]["uBytes"] * i)
-        tag = {}
-        for key, fmt in TagInfoDescription:
-            val = f.read_f(fmt)
-            if len(val) == 1:
-                tag[key] = val[0]
-            else:
-                tag[key] = np.array(val)
-        listTag.append(tag)
-
-    header["listTag"] = listTag
+    f.seek(sections["TagSection"]["uBlockIndex"] * BLOCKSIZE)
+    tag_bytes = f.read(sections["TagSection"]["llNumEntries"] * TagInfoDescription.itemsize)
+    header["listTag"] = np.frombuffer(tag_bytes, dtype=TagInfoDescription)
 
     # DAC sections
     header["listDACInfo"] = []
@@ -1341,12 +1334,24 @@ def safe_decode_units(s):
     ("lADCUnitsIndex", "i"),
 ]
 
-TagInfoDescription = [
-    ("lTagTime", "i"),
-    ("sComment", "56s"),
-    ("nTagType", "h"),
-    ("nVoiceTagNumber_or_AnnotationIndex", "h"),
-]
+TagInfoDescription = np.dtype(
+    [
+        ("lTagTime", "<i4"),
+        ("sComment", "S56"),
+        ("nTagType", "<i2"),
+        ("nVoiceTagNumber_or_AnnotationIndex", "<i2"),
+    ]
+)
+
+# Event channel name and id for each tag type (nTagType), after the ABF_*TAG constants in ABFHEADR.H
+TAG_TYPE_NAMES = {
+    0: "time_tag",
+    1: "comment_tag",
+    2: "external_tag",
+    3: "voice_tag",
+    4: "new_file_tag",
+    5: "annotation_tag",
+}
 
 DACInfoDescription = [
     ("nDACNum", "h"),

neo/test/iotest/test_axonio.py

@@ -15,11 +15,12 @@ class TestAxonIO(BaseTestIO, unittest.TestCase):
     entities_to_test = TestAxonRawIO.entities_to_test
     ioclass = AxonIO
 
-    def test_annotations(self):
+    def test_tag_comment_is_event_label(self):
         reader = AxonIO(filename=self.get_local_path("axon/File_axon_2.abf"))
         bl = reader.read_block()
         ev = bl.segments[0].events[0]
-        assert "comments" in ev.annotations
+        assert ev.name == "comment_tag"
+        assert ev.labels[0] == "Clampex start acquisition"
 
     def test_read_protocol(self):
         for f in self.entities_to_test:

neo/test/rawiotest/test_axonrawio.py

@@ -30,6 +30,65 @@ def test_read_raw_protocol(self):
 
         reader.read_raw_protocol()
 
+    def test_tag_event_in_synch_time_units(self):
+        # ABF1 locates the tag section in 512-byte blocks, and lTagTime counts in fSynchTimeUnit
+        # microseconds (12.5 here, with four channels). The tag at 1.8 s is inside the first of
+        # 16 sweeps of 2 s, so only segment 0 gets it.
+        reader = AxonRawIO(filename=self.get_local_path("axon/intracellular_data/tags/episodic_inside_sweep.abf"))
+        reader.parse_header()
+        timestamps, _, labels = reader.get_event_timestamps(seg_index=0, event_channel_index=0)
+        times = reader.rescale_event_timestamp(timestamps)
+        self.assertEqual(labels.tolist(), ["C9, L3,  RMP -66.8 mv"])
+        self.assertAlmostEqual(times[0], 1.8)
+        counts = [reader.event_count(seg_index=seg_index, event_channel_index=0) for seg_index in range(16)]
+        self.assertEqual(counts, [1] + [0] * 15)
+
+    def test_tag_event_in_samples_when_synch_time_unit_is_zero(self):
+        # With fSynchTimeUnit = 0, lTagTime counts samples summed over all channels (two channels
+        # at 10 kHz here). Expected values from pyABF.
+        filename = self.get_local_path("axon/intracellular_data/tags/time_in_samples_two_channels.abf")
+        reader = AxonRawIO(filename=filename)
+        reader.parse_header()
+        timestamps, _, labels = reader.get_event_timestamps(event_channel_index=0)
+        times = reader.rescale_event_timestamp(timestamps)
+        self.assertEqual(labels.tolist(), ["Digital Outputs => 00000001"])
+        self.assertAlmostEqual(times[0], 0.9216)
+
+    def test_one_event_channel_per_tag_type(self):
+        # One time, comment, voice and annotation tag, written by the vendor library.
+        filename = self.get_local_path("axon/intracellular_data/tags/annotation_and_voice_uncompressed.abf")
+        reader = AxonRawIO(filename=filename)
+        reader.parse_header()
+        expected_names = ["time_tag", "comment_tag", "voice_tag", "annotation_tag"]
+        self.assertEqual(reader.header["event_channels"]["name"].tolist(), expected_names)
+        self.assertEqual(reader.header["event_channels"]["id"].tolist(), expected_names)
+        labels = [reader.get_event_timestamps(event_channel_index=index)[2].tolist() for index in range(4)]
+        self.assertEqual(labels, [[""], ["comment tag"], ["<Voice Tag>"], ["annotation tag"]])
+
+    def test_tag_between_sweeps_goes_to_no_segment(self):
+        # Two comment tags, at 173.98 s and 294.86 s, fall between sweeps of 10 ms that start
+        # every 2 s.
+        reader = AxonRawIO(filename=self.get_local_path("axon/intracellular_data/tags/abf2.abf"))
+        reader.parse_header()
+        counts = [reader.event_count(seg_index=seg_index) for seg_index in range(reader.segment_count(0))]
+        self.assertEqual(sum(counts), 0)
+
+    def test_tag_event_time_window_in_seconds(self):
+        # t_start and t_stop are in seconds and are compared with the converted tag time (0.92 s
+        # here, from a count of 18432), not with the raw count.
+        filename = self.get_local_path("axon/intracellular_data/tags/time_in_samples_two_channels.abf")
+        reader = AxonRawIO(filename=filename)
+        reader.parse_header()
+        timestamps, _, _ = reader.get_event_timestamps(event_channel_index=0, t_stop=1.0)
+        self.assertEqual(timestamps.size, 1)
+        timestamps, _, _ = reader.get_event_timestamps(event_channel_index=0, t_start=1.0)
+        self.assertEqual(timestamps.size, 0)
+
+    def test_file_without_tags_has_no_event_channel(self):
+        reader = AxonRawIO(filename=self.get_local_path("axon/File_axon_1.abf"))
+        reader.parse_header()
+        self.assertEqual(reader.event_channels_count(), 0)
+
     def test_empty_channel_name_gets_fallback(self):
         # Some ABF files store a blank ADC channel name, which collapses to "" after space
         # stripping and leaves the channel unaddressable by name. A positional fallback (ch{id})