From e638d117e84ecd9ad458a70898dcb5945600df98 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Sun, 8 Nov 2026 22:24:37 -0600 Subject: [PATCH] Reconstruct ABF epoch shapes, the last epoch level and ABF1 protocols Draw ramps and pulse, triangle, cosine and biphasic trains, start each sweep at the last epoch level when the protocol asks for it, read the protocol of ABF1 files from version 1.6 on, and warn when a file uses user lists, alternating outputs, stimulus-file waveforms or conditioning trains, which are not reconstructed yet. --- neo/io/axonio.py | 5 +- neo/rawio/axonrawio.py | 275 ++++++++++++++++++++++----- neo/test/rawiotest/test_axonrawio.py | 79 ++++++++ 3 files changed, 309 insertions(+), 50 deletions(-) diff --git a/neo/io/axonio.py b/neo/io/axonio.py index 9e556e91d..7e33739af 100644 --- a/neo/io/axonio.py +++ b/neo/io/axonio.py @@ -45,9 +45,8 @@ def __init__(self, filename): def read_protocol(self): """ - Read the protocol waveform of the file, if present; - function works with ABF2 only. Protocols can be reconstructed - from the ABF1 header. + Read the protocol waveform of the file, if present. + See ``AxonRawIO.read_raw_protocol`` for what is reconstructed. Returns ------- diff --git a/neo/rawio/axonrawio.py b/neo/rawio/axonrawio.py index b0603e197..b21f5272b 100644 --- a/neo/rawio/axonrawio.py +++ b/neo/rawio/axonrawio.py @@ -49,6 +49,7 @@ import os import struct import datetime +import warnings from io import open, BufferedReader import numpy as np @@ -383,7 +384,7 @@ def _rescale_event_timestamp(self, event_timestamps, dtype, event_channel_index) def read_raw_protocol(self): """ - Read stimulus waveforms (DAC output protocols) from ABF2 files. + Read stimulus waveforms (DAC output protocols) from ABF files. Purpose ------- @@ -404,10 +405,11 @@ def read_raw_protocol(self): Limitations ----------- - - ABF2 files only: ABF1 files don't store protocol sections, though - protocols can theoretically be reconstructed from ABF1 headers + - ABF1 files need version 1.6 or later when a DAC plays a waveform - Waveforms are RECONSTRUCTED from epoch definitions, not recorded data - - Some complex protocols may not reconstruct perfectly (see TODO in code) + - User lists, alternating outputs, stimulus-file waveforms and conditioning + trains are not reconstructed yet: a warning names them and the waveform + follows the epoch table only Returns ------- @@ -443,60 +445,230 @@ def read_raw_protocol(self): Notes ----- - Waveforms are constructed from: - - Holding levels (fDACHoldingLevel) - - Epoch definitions (EpochPerDACSection) describing step changes + Waveforms are constructed the way the vendor library (``ABFH_GetWaveform`` + in the ABF File Support Pack) does, from: + - Holding levels (fDACHoldingLevel), for the first and last 1/64 of the + sweep + - Epoch definitions (EpochPerDACSection) in their shape: step, ramp, or + pulse, triangle, cosine or biphasic train - Episode number (waveforms can change across episodes via increment parameters) + - The "use last epoch level" option (nInterEpisodeLevel), where a sweep + starts at the level the previous one ended + A DAC with no waveform, and every DAC outside episodic stimulation mode, + stays at its holding level. Author: JS Nowacki """ - info = self._axon_info - - if info["fFileVersionNumber"] < 2.0: - raise IOError("Protocol section is only present in ABF2 files.") - - nADC = info["sections"]["ADCSection"]["llNumEntries"] # Number of ADC channels - nDAC = info["sections"]["DACSection"]["llNumEntries"] # Number of DAC channels - nSam = int(info["protocol"]["lNumSamplesPerEpisode"] / nADC) # Number of samples per episode - nEpi = info["lActualEpisodes"] # Actual number of episodes + protocol = _get_protocol(self._axon_info) + unsupported = _unsupported_protocol_features(protocol) + if unsupported: + warnings.warn( + f"read_raw_protocol does not reconstruct {', '.join(unsupported)} yet, so the waveforms of " + f"{self.filename} follow the epoch table only and can differ from what was played." + ) - # Make a list of segments with analog signals with just holding levels - # List of segments relates to number of episodes, as for recorded data + n_channels = protocol["n_channels"] sigs_by_segments = [] - for epiNum in range(nEpi): - # One analog signal for each DAC in segment (episode) + for episode in range(1, protocol["n_episodes"] + 1): signals = [] - for DACNum in range(nDAC): - sig = np.ones(nSam) * info["listDACInfo"][DACNum]["fDACHoldingLevel"] - # If there are epoch infos for this DAC - if DACNum in info["dictEpochInfoPerDAC"]: - # Save last sample index - i_last = int(nSam * 15625 / 10**6) - # TODO guess for first holding - # Go over EpochInfoPerDAC and change the analog signal - # according to the epochs - epochInfo = info["dictEpochInfoPerDAC"][DACNum] - for epochNum, epoch in epochInfo.items(): - i_begin = i_last - i_end = i_last + epoch["lEpochInitDuration"] + epoch["lEpochDurationInc"] * epiNum - dif = i_end - i_begin - sig[i_begin:i_end] = np.ones(dif) * ( - epoch["fEpochInitLevel"] + epoch["fEpochLevelInc"] * epiNum - ) - i_last += epoch["lEpochInitDuration"] + epoch["lEpochDurationInc"] * epiNum - signals.append(sig) + for dac in range(protocol["n_dacs"]): + if protocol["operation_mode"] == 5 and protocol["waveform_enable"][dac]: + waveform = _epoch_table_waveform(protocol, dac, episode) + else: + waveform = np.full(protocol["samples_per_episode"], protocol["holding_level"][dac]) + # the waveform is built on the interleaved samples of all channels, keep one per sampling sequence + signals.append(waveform[::n_channels]) sigs_by_segments.append(signals) - sig_names = [] - sig_units = [] - for DACNum in range(nDAC): - name = info["listDACInfo"][DACNum]["DACChNames"].decode("utf-8") - units = safe_decode_units(info["listDACInfo"][DACNum]["DACChUnits"]) - sig_names.append(name) - sig_units.append(units) + return sigs_by_segments, protocol["dac_names"], protocol["dac_units"] + + +ABF_EPOCH_COUNT = 10 +EPOCH_KEYS = [ + "nEpochType", + "fEpochInitLevel", + "fEpochLevelInc", + "lEpochInitDuration", + "lEpochDurationInc", + "lEpochPulsePeriod", + "lEpochPulseWidth", +] - return sigs_by_segments, sig_names, sig_units + +def _get_protocol(info): + """ + Collect the fields that define the command waveforms, in one layout for ABF1 and ABF2. + + Per-DAC epoch fields are arrays of shape (n_dacs, ABF_EPOCH_COUNT). Sample counts are in + samples of the interleaved data (all channels), as in the header. + """ + version = info["fFileVersionNumber"] + if version < 2.0: + # ABF1 has four DACs, and the first two can play a waveform + n_dacs = 4 + epochs = {} + for key in EPOCH_KEYS: + epochs[key] = np.zeros((n_dacs, ABF_EPOCH_COUNT), dtype=np.asarray(info[key]).dtype) + epochs[key][:2] = np.reshape(info[key], (2, ABF_EPOCH_COUNT)) + + def per_dac(key): + values = np.zeros(n_dacs, dtype=np.asarray(info[key]).dtype) + values[: len(info[key])] = info[key] + return values + + protocol = { + "operation_mode": info["nOperationMode"], + "n_channels": info["nADCNumChannels"], + "samples_per_episode": info["lNumSamplesPerEpisode"], + "episodes_per_run": info["lEpisodesPerRun"], + "dac_names": [name.decode("utf-8", errors="replace").strip() for name in info["sDACChannelName"]], + "dac_units": [safe_decode_units(units) for units in info["sDACChannelUnits"]], + "holding_level": per_dac("fDACHoldingLevel"), + "waveform_enable": per_dac("nWaveformEnable"), + "waveform_source": per_dac("nWaveformSource"), + "inter_episode_level": per_dac("nInterEpisodeLevel"), + "condit_enable": per_dac("nConditEnable"), + "user_lists": bool(np.any(info["nULEnable"])), + "alternate_dac_outputs": bool(info["nAlternateDACOutputState"]), + } + if version < 1.6: + # Headers older than 1.6 are 2 KB long and keep a single waveform in other fields, so + # the fields after byte 2048 hold unrelated bytes. No file we have uses that waveform. + if info["_nWaveformSource"]: + raise NeoReadWriteError("Protocol reconstruction needs ABF 1.6 or later when a DAC plays a waveform") + protocol["waveform_enable"] = np.zeros(n_dacs, dtype=int) + else: + dac_infos = info["listDACInfo"] + n_dacs = len(dac_infos) + epochs = {key: np.zeros((n_dacs, ABF_EPOCH_COUNT)) for key in EPOCH_KEYS} + for dac, epochs_of_dac in info["dictEpochInfoPerDAC"].items(): + for epoch_num, epoch in epochs_of_dac.items(): + for key in EPOCH_KEYS: + epochs[key][dac, epoch_num] = epoch[key] + + def per_dac(key): + return np.array([dac_info[key] for dac_info in dac_infos]) + + protocol_section = info["protocol"] + protocol = { + "operation_mode": protocol_section["nOperationMode"], + "n_channels": info["sections"]["ADCSection"]["llNumEntries"], + "samples_per_episode": protocol_section["lNumSamplesPerEpisode"], + "episodes_per_run": protocol_section["lEpisodesPerRun"], + "dac_names": [dac_info["DACChNames"].decode("utf-8") for dac_info in dac_infos], + "dac_units": [safe_decode_units(dac_info["DACChUnits"]) for dac_info in dac_infos], + "holding_level": per_dac("fDACHoldingLevel"), + "waveform_enable": per_dac("nWaveformEnable"), + "waveform_source": per_dac("nWaveformSource"), + "inter_episode_level": per_dac("nInterEpisodeLevel"), + "condit_enable": per_dac("nConditEnable"), + # the vendor reader plays every user list in the section, whatever its nULEnable says + "user_lists": info["sections"]["UserListSection"]["llNumEntries"] > 0, + "alternate_dac_outputs": bool(protocol_section["nAlternateDACOutputState"]), + } + protocol.update(epochs) + protocol["n_dacs"] = n_dacs + protocol["n_episodes"] = info["lActualEpisodes"] + return protocol + + +def _unsupported_protocol_features(protocol): + """The protocol features that change a waveform but are not reconstructed yet.""" + playing = np.flatnonzero(protocol["waveform_enable"]) + if protocol["operation_mode"] != 5 or playing.size == 0: + return [] + features = [] + if protocol["user_lists"]: + features.append("user lists") + if protocol["alternate_dac_outputs"]: + features.append("alternating outputs") + if np.any(protocol["waveform_source"][playing] == 2): + features.append("stimulus-file waveforms") + if np.any(protocol["condit_enable"][playing]): + features.append("conditioning trains") + return features + + +def _holding_length(protocol): + """Samples (interleaved) of the first holding: 1/64 of the sweep, rounded down to whole + sampling sequences and at least one.""" + n_channels = protocol["n_channels"] + length = protocol["samples_per_episode"] // 64 + length -= length % n_channels + return max(length, n_channels) + + +def _epoch_level(protocol, dac, episode, epoch): + return protocol["fEpochInitLevel"][dac, epoch] + (episode - 1) * protocol["fEpochLevelInc"][dac, epoch] + + +def _holding_level(protocol, dac, episode): + """The level before the first epoch: with the "use last epoch level" option, the level of the + last epoch of the previous sweep.""" + if not protocol["inter_episode_level"][dac] or episode < 2: + return protocol["holding_level"][dac] + enabled = np.flatnonzero(protocol["nEpochType"][dac]) + if enabled.size == 0: + return protocol["holding_level"][dac] + return _epoch_level(protocol, dac, episode - 1, enabled[-1]) + + +def _epoch_shape(epoch_type, n, start_level, level, period, width): + """One epoch's samples, from the previous level (start_level) to this epoch's level.""" + position = np.arange(n) + if epoch_type == 1: # step + return np.full(n, level, dtype="float64") + if epoch_type == 2: # ramp + return start_level + (level - start_level) * (position + 1) / n + phase = position % period + if epoch_type == 3: # rectangular pulse train + return np.where(phase < width, level, start_level) + if epoch_type == 4: # triangle train: up over the width, down over the rest of the period + rising = start_level + (level - start_level) * phase / width + falling = level - (level - start_level) * (phase - width) / (period - width) + return np.where(phase < width, rising, falling) + if epoch_type == 5: # cosine train + return start_level + (level - start_level) * (1 - np.cos(2 * np.pi * phase / period)) / 2 + if epoch_type == 7: # biphasic train: level for half the width, its mirror for the other half + return np.where(phase < width // 2, level, np.where(phase < width, 2 * start_level - level, start_level)) + raise NeoReadWriteError(f"Epoch type {epoch_type} is not supported") + + +def _epoch_table_waveform(protocol, dac, episode): + """Port of GenerateWaveform (ABFHWAVE.CPP): the interleaved waveform of one DAC in one sweep.""" + n_samples = protocol["samples_per_episode"] + n_channels = protocol["n_channels"] + holding_level = _holding_level(protocol, dac, episode) + waveform = np.full(n_samples, holding_level, dtype="float64") + # sweeps added after the protocol ran (off-line) stay at holding + if episode > protocol["episodes_per_run"]: + return waveform + + position = _holding_length(protocol) + start_level = holding_level + for epoch in range(ABF_EPOCH_COUNT): + epoch_type = protocol["nEpochType"][dac, epoch] + if epoch_type == 0: + continue + duration = ( + protocol["lEpochInitDuration"][dac, epoch] + (episode - 1) * protocol["lEpochDurationInc"][dac, epoch] + ) + n = int(duration) * n_channels + if n <= 0: + continue + if position + n > n_samples: + raise NeoReadWriteError(f"The epochs of DAC {dac} in sweep {episode} are longer than the sweep") + level = _epoch_level(protocol, dac, episode, epoch) + period = int(protocol["lEpochPulsePeriod"][dac, epoch]) * n_channels + width = int(protocol["lEpochPulseWidth"][dac, epoch]) * n_channels + waveform[position : position + n] = _epoch_shape(epoch_type, n, start_level, level, period, width) + position += n + start_level = level + + end_level = start_level if protocol["inter_episode_level"][dac] else protocol["holding_level"][dac] + waveform[position:] = end_level + return waveform def parse_axon_soup(filename): @@ -1173,7 +1345,11 @@ def safe_decode_units(s): ("fInstrumentOffset", 986, "16f"), ("fSignalGain", 1050, "16f"), ("fSignalOffset", 1114, "16f"), + ("sDACChannelName", 1306, "10s" * 4), + ("sDACChannelUnits", 1346, "8s" * 4), + ("fDACHoldingLevel", 1394, "4f"), ("nDigitalEnable", 1436, "h"), + ("_nWaveformSource", 1438, "h"), ("nActiveDACChannel", 1440, "h"), ("nDigitalHolding", 1584, "h"), ("nDigitalInterEpisode", 1586, "h"), @@ -1182,6 +1358,8 @@ def safe_decode_units(s): ("lDACFileNumEpisodes", 2056, "2i"), ("fDACCalibrationFactor", 2074, "4f"), ("fDACCalibrationOffset", 2090, "4f"), + ("lEpochPulsePeriod", 2136, "20i"), + ("lEpochPulseWidth", 2216, "20i"), ("nWaveformEnable", 2296, "2h"), ("nWaveformSource", 2300, "2h"), ("nInterEpisodeLevel", 2304, "2h"), @@ -1190,9 +1368,12 @@ def safe_decode_units(s): ("fEpochLevelInc", 2428, "20f"), ("lEpochInitDuration", 2508, "20i"), ("lEpochDurationInc", 2588, "20i"), + ("nConditEnable", 3260, "2h"), + ("nULEnable", 3360, "4h"), ("nTelegraphEnable", 4512, "16h"), ("fTelegraphAdditGain", 4576, "16f"), ("sProtocolPath", 4898, "384s"), + ("nAlternateDACOutputState", 5876, "h"), ] headerDescriptionV2 = [ diff --git a/neo/test/rawiotest/test_axonrawio.py b/neo/test/rawiotest/test_axonrawio.py index d742dae02..d3cb776a4 100644 --- a/neo/test/rawiotest/test_axonrawio.py +++ b/neo/test/rawiotest/test_axonrawio.py @@ -30,6 +30,85 @@ def test_read_raw_protocol(self): reader.read_raw_protocol() + def test_read_raw_protocol_epoch_shapes(self): + # One file per epoch shape, on DAC 1. The first holding is 1/64 of the sweep (312 of 20000 + # samples), then a step to -20 mV, a step to -10 mV and the shape from -10 mV to 25 mV over + # samples 5312 to 15312, with trains of period 1000 and width 200. Expected values from the + # vendor library (ABF_GetWaveform, ABF File Support Pack 2.03). + start = 5312 + expected = { + "step": {0: 25.0, 9999: 25.0}, + "ramp": {0: -9.99825, 5000: 7.50175, 9999: 24.99825}, + "pulse_train": {0: 25.0, 199: 25.0, 200: -10.0, 1000: 25.0}, + "triangle_train": {100: 7.5, 200: 25.0, 600: 7.5}, + "cosine_train": {0: -10.0, 250: 7.5, 500: 25.0}, + "biphasic_train": {50: 25.0, 150: -45.0, 250: -10.0}, + } + for shape, values in expected.items(): + with self.subTest(shape=shape): + filename = self.get_local_path(f"axon/intracellular_data/read_raw_protocol/{shape}.abf") + reader = AxonRawIO(filename=filename) + reader.parse_header() + sigs_by_segments, _, _ = reader.read_raw_protocol() + waveform = sigs_by_segments[0][1] + self.assertEqual(waveform[311], -10.0) + self.assertEqual(waveform[312], -20.0) + self.assertEqual(waveform[start - 1], -10.0) + for offset, value in values.items(): + self.assertAlmostEqual(waveform[start + offset], value, places=4) + self.assertEqual(waveform[start + 10000], -10.0) + + def test_read_raw_protocol_increments(self): + # A level that grows by 10 pA each sweep (DAC 0), and a ramp on DAC 1 whose duration grows + # by 20 samples each sweep, so the step after it starts later. Expected values from the + # vendor library. + reader = AxonRawIO( + filename=self.get_local_path("axon/intracellular_data/read_raw_protocol/level_increment.abf") + ) + reader.parse_header() + sigs_by_segments, _, _ = reader.read_raw_protocol() + self.assertEqual([sigs[0][2937] for sigs in sigs_by_segments], [-50.0, -40.0, -30.0]) + + filename = self.get_local_path("axon/intracellular_data/read_raw_protocol/duration_increment.abf") + reader = AxonRawIO(filename=filename) + reader.parse_header() + sigs_by_segments, _, _ = reader.read_raw_protocol() + for sweep, sigs in enumerate(sigs_by_segments): + step_start = 4170 + 20 * sweep + self.assertLess(sigs[1][step_start - 2], 1.0 + sweep) + self.assertEqual(sigs[1][step_start], 1.0 + sweep) + + def test_read_raw_protocol_last_epoch_level(self): + # A ramp to 50 mV that grows by 20 mV each sweep, with the level between sweeps set to the + # last epoch's level, so each sweep starts where the previous one ended. + reader = AxonRawIO( + filename=self.get_local_path("axon/intracellular_data/read_raw_protocol/last_epoch_level_ramp.abf") + ) + reader.parse_header() + sigs_by_segments, _, _ = reader.read_raw_protocol() + self.assertEqual([sigs[0][0] for sigs in sigs_by_segments], [0.0, 50.0, 70.0, 90.0]) + self.assertAlmostEqual(sigs_by_segments[0][0][-1], 50.0, places=4) + + def test_read_raw_protocol_abf1(self): + # ABF 1.83, four channels, a step on DAC 1 from 0 to -400 pA growing by 50 pA each sweep. + # Expected values from the vendor library. + reader = AxonRawIO(filename=self.get_local_path("axon/intracellular_data/tags/episodic_inside_sweep.abf")) + reader.parse_header() + sigs_by_segments, names, units = reader.read_raw_protocol() + self.assertEqual((names[1], units[1]), ("Cmd 1", "pA")) + waveform = sigs_by_segments[2][1] + self.assertEqual(waveform.size, 20000) + self.assertEqual(waveform[1311], 0.0) + self.assertEqual(waveform[1312], -300.0) + self.assertEqual(waveform[7311], -300.0) + self.assertEqual(waveform[7312], 0.0) + + def test_read_raw_protocol_warns_on_unsupported_features(self): + reader = AxonRawIO(filename=self.get_local_path("axon/intracellular_data/read_raw_protocol/user_list.abf")) + reader.parse_header() + with self.assertWarnsRegex(UserWarning, "user lists"): + reader.read_raw_protocol() + 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})