From b3f8bfab6ce698eaebf77fb4292a934c803a381a Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 1 Sep 2026 13:27:16 +0100 Subject: [PATCH 1/2] Add burst segmentation to dcio.analysis Move burst splitting into the record layer, so the one implementation sits below SCALCS, EKDIST and HJCFIT rather than beside them. This is additive: nothing is removed from SCALCS, whose scsim functions keep working unchanged until they can be turned into re-exports. The implementation is SCALCS's, chosen over EKDIST's because it is the one with flag handling, the one with tests, and the one that reproduces the 2017 segmentation exactly. Verified identical to it over 500 random records in both flag modes, and on examples/scn/glyr_experimental.scn at three critical times. Adds a record-level entry point, bursts_from_record, which segments the record's periods rather than its resolved intervals. impose_resolution emits a fresh open interval at every change of fitted amplitude, so a record idealised with sub-conductance levels does not alternate open/shut: the shipped experimental example has 1770 consecutive open-open pairs in 3232 resolved intervals, and segmenting those gives 30 of 63 bursts that end on a shut. The missed-events likelihood needs alternating sequences. Burst count and lengths are unaffected; openings per burst are 731 rather than 2501, counting an opening once instead of once per sub-level. No HJCFIT sample record contains a sub-conductance level, so this fixes nothing currently broken there -- it removes a way for it to break later. Tests carry over the convention cases that pinned this behaviour in SCALCS, and add golden values for the shipped experimental record. Co-Authored-By: Claude Opus 5 --- README.md | 35 ++++++ dcio/analysis/__init__.py | 8 ++ dcio/analysis/bursts.py | 243 ++++++++++++++++++++++++++++++++++++++ tests/test_bursts.py | 212 +++++++++++++++++++++++++++++++++ 4 files changed, 498 insertions(+) create mode 100644 dcio/analysis/bursts.py create mode 100644 tests/test_bursts.py diff --git a/README.md b/README.md index e319cd3..d0d70e9 100644 --- a/README.md +++ b/README.md @@ -68,6 +68,41 @@ write("output.scn", intervals_s, amplitudes, flags, | `shut_mask` | `ndarray[bool]` | True for shut intervals | | `usable_mask` | `ndarray[bool]` | True when `flags & 8 == 0` | +## Record analysis + +Beyond file I/O, `dcio.analysis` carries the record layer shared by the rest of +the DCPROGS stack: applying a dead time, grouping the result into open and shut +periods, and cutting it into bursts. + +```python +from dcio.formats.scn import read +from dcio.analysis import from_scn, bursts_from_record + +record = from_scn(read("myrecording.scn"), tres=25e-6) # 25 us dead time +lengths, n_openings = bursts_from_record(record, tcrit=4e-3) + +print(len(lengths), "bursts,", n_openings.mean(), "openings each") +``` + +`bursts_from_record` segments the record's **periods**, not its resolved +intervals. `impose_resolution` emits a fresh open interval at every change of +fitted amplitude, so a record idealised with sub-conductance levels contains +runs of consecutive openings; segmenting those directly gives bursts that do +not alternate open/shut. Burst count and lengths are the same either way, the +number of openings per burst is not. + +Two conventions are worth stating, because implementations in this stack have +differed on them: + +- no gap longer than `tcrit` is required before the first burst -- the first + defined opening starts one; +- an unusable interval ends the burst before it, and its (meaningless) + duration is never compared with `tcrit`. + +Both follow SCAN's time-course fitting, which leaves the final interval of a +record with no defined length. Dropping the runs at each end instead loses two +bursts from every record. + ## Running tests ```bash diff --git a/dcio/analysis/__init__.py b/dcio/analysis/__init__.py index 9b719fc..a98e09a 100644 --- a/dcio/analysis/__init__.py +++ b/dcio/analysis/__init__.py @@ -1,5 +1,10 @@ """Single-channel record analysis.""" +from dcio.analysis.bursts import ( + bursts_from_record, + extract_burst_intervals, + extract_bursts, +) from dcio.analysis.record import ( Periods, SingleChannelRecord, @@ -14,4 +19,7 @@ "from_scn", "impose_resolution", "set_periods", + "extract_bursts", + "extract_burst_intervals", + "bursts_from_record", ] diff --git a/dcio/analysis/bursts.py b/dcio/analysis/bursts.py new file mode 100644 index 0000000..4d8b82b --- /dev/null +++ b/dcio/analysis/bursts.py @@ -0,0 +1,243 @@ +"""Burst segmentation of an idealised single-channel record. + +A *burst* is a run of openings and short shuttings delimited by shut intervals +longer than a critical time ``tcrit``. Segmentation sits directly on top of +:mod:`dcio.analysis.record`: it takes the resolved intervals a dead time +produced and cuts them into bursts. + +Typical workflow:: + + from dcio.formats.scn import read + from dcio.analysis.record import from_scn + from dcio.analysis.bursts import extract_bursts + + scr = from_scn(read("file.scn"), tres=25e-6) + lengths, n_openings = extract_bursts( + scr.resolved_intervals, scr.resolved_amplitudes, + tcrit=4e-3, flags=scr.resolved_flags, + ) + +or, equivalently, straight from the record:: + + from dcio.analysis.bursts import bursts_from_record + + lengths, n_openings = bursts_from_record(scr, tcrit=4e-3) + +Convention +---------- +The convention is the one EKDIST states in ``Bursts.slice_bursts`` and dcpyps +followed: + +1. no gap longer than ``tcrit`` is required before the first burst of a + record -- the first defined opening is a valid burst start; +2. an unusable interval is a valid end of burst. + +Both matter. Time-course fitting in SCAN leaves the last interval of a record +with no defined length, flagged unusable; it still ends the burst before it. A +leading interval may likewise be bad, and is discarded, but the first defined +opening after it starts a real burst. Dropping the runs at both ends -- as +earlier versions of this code did -- loses two bursts from every record, which +at 30 uM in the Burzomato 2004 set is a third of the data. +""" + +from __future__ import annotations + +import numpy as np + +from dcio.formats.scn import FLAG_UNUSABLE + +__all__ = [ + "extract_bursts", + "extract_burst_intervals", + "bursts_from_record", +] + + +def _burst_segments(intervals, amplitudes, tcrit, flags=None): + """Split a record into bursts and return them as (interval, amplitude) pairs. + + Shared by :func:`extract_bursts` and :func:`extract_burst_intervals`, which + differ only in what they report about one and the same segmentation. + + Parameters + ---------- + intervals, amplitudes : array_like + Alternating interval record (ideal or apparent). Durations in + seconds, amplitudes in pA with 0 meaning shut. + tcrit : float + Critical shut time separating within- from between-burst gaps [s]. + Its magnitude only -- see :func:`bursts_from_record` for the sign + convention that HJCFIT places on the same number. + flags : array_like of int, optional + Per-interval SCN property flags. An interval is unusable when + ``flags & FLAG_UNUSABLE`` is set. Without them no interval is treated + as unusable, which is right for a simulated record and wrong for an + experimental one. + + Returns + ------- + list of list of (float, float) + One list of ``(interval, amplitude)`` pairs per burst, each starting + and ending on an opening. + """ + intervals = np.asarray(intervals, float) + amplitudes = np.asarray(amplitudes, float) + if intervals.size == 0: + return [] + + if flags is None: + unusable = np.zeros(intervals.shape, dtype=bool) + else: + unusable = (np.asarray(flags, int) & FLAG_UNUSABLE) != 0 + + # A burst ends at a between-burst gap or at an interval of unknown length. + # An unusable interval has no measured duration, so it is never compared + # with tcrit. + # Strictly greater: tcrit is the time such that gaps *longer* than it are + # between-burst. EKDIST's Bursts.slice_bursts uses the same test, and a + # shut interval exactly equal to tcrit is within-burst by that reading. + separator = ((amplitudes == 0.0) & (intervals > tcrit) & ~unusable) | unusable + + # Trim to the first and last defined opening, so the record begins and + # ends on one. What lies outside is a shut interval or an unusable one, + # and neither belongs to a burst. + opening = (amplitudes != 0.0) & ~unusable + if not opening.any(): + return [] + lo = int(np.argmax(opening)) + hi = int(len(opening) - np.argmax(opening[::-1])) + + raw, seg = [], [] + for t, a, sep in zip(intervals[lo:hi], amplitudes[lo:hi], separator[lo:hi]): + if sep: + if seg: + raw.append(seg) + seg = [] + else: + seg.append((t, a)) + if seg: + raw.append(seg) + + bursts = [] + for seg in raw: + while seg and seg[0][1] == 0.0: # trim leading shut + seg = seg[1:] + while seg and seg[-1][1] == 0.0: # trim trailing shut + seg = seg[:-1] + if seg: + bursts.append(seg) + return bursts + + +def extract_bursts(intervals, amplitudes, tcrit, flags=None): + """Split a record into bursts at shut intervals longer than ``tcrit``. + + Each burst is trimmed to start and end on an opening. + + Parameters + ---------- + intervals, amplitudes : array_like + Alternating interval record (ideal or apparent). + tcrit : float + Critical shut time separating within- from between-burst gaps [s]. + flags : array_like of int, optional + Per-interval SCN property flags; see :func:`_burst_segments`. + + Returns + ------- + lengths : ndarray + Burst lengths [s] (first opening start to last opening end). + n_openings : ndarray of int + Number of (apparent) openings in each burst. + + See Also + -------- + extract_burst_intervals : the same bursts, as interval sequences. + """ + bursts = _burst_segments(intervals, amplitudes, tcrit, flags) + lengths = [sum(t for t, _ in seg) for seg in bursts] + nops = [sum(1 for _, a in seg if a != 0.0) for seg in bursts] + return np.array(lengths), np.array(nops, dtype=int) + + +def extract_burst_intervals(intervals, amplitudes, tcrit, flags=None): + """Split a record into bursts and return the intervals of each. + + The segmentation is that of :func:`extract_bursts`, which reduces each + burst to its length and its number of openings. Maximum-likelihood + fitting of missed-events mechanisms needs the interval sequences + themselves: the HJC likelihood is a product of matrices, one per interval, + so the order and the individual durations both matter. + + Parameters + ---------- + intervals, amplitudes : array_like + Alternating interval record (ideal or apparent). + tcrit : float + Critical shut time separating within- from between-burst gaps [s]. + flags : array_like of int, optional + Per-interval SCN property flags; see :func:`_burst_segments`. + + Returns + ------- + list of ndarray + One array of interval durations [s] per burst, alternating open and + shut and both starting and ending with an opening -- so every array + has odd length, which is what the missed-events likelihood requires. + + See Also + -------- + extract_bursts : the same bursts, as lengths and opening counts. + """ + return [np.array([t for t, _ in seg], dtype=float) + for seg in _burst_segments(intervals, amplitudes, tcrit, flags)] + + +def bursts_from_record(record, tcrit, intervals_only=False): + """Segment a :class:`~dcio.analysis.record.SingleChannelRecord` into bursts. + + Convenience over :func:`extract_bursts` that takes the intervals, + amplitudes and flags from the record itself, so a caller can get neither + the flags nor the input series wrong. + + Segmentation runs on the record's **periods**, not on its resolved + intervals, and that distinction is not cosmetic. + :func:`~dcio.analysis.record.impose_resolution` emits a fresh open + interval whenever the fitted amplitude changes, so a record idealised with + sub-conductance levels contains runs of consecutive open intervals: the + experimental example shipped with dcio has 1770 such adjacencies in 3232 + resolved intervals. Segmenting those directly gives bursts that do not + alternate open/shut, and counts each sub-level as a separate opening. + :func:`~dcio.analysis.record.set_periods` merges them, which is what makes + the burst an alternating sequence the missed-events likelihood can consume. + + Burst count and burst lengths are identical either way; the number of + openings per burst is not. + + Parameters + ---------- + record : SingleChannelRecord + A record with a dead time already applied. + tcrit : float + Critical shut time [s]. Only its magnitude is used. The sign carries + no meaning here, but it does downstream: HJCFIT reads a negative + ``tcrit`` as a flag selecting equilibrium vectors (Colquhoun & Hawkes + 1982) over CHS vectors (Colquhoun, Hawkes & Srodzinski 1996), and the + same number is passed to both. Taking the magnitude here means a + caller can hand the same value to either without thinking about it. + intervals_only : bool, default False + Return the per-burst interval sequences instead of lengths and + opening counts. + + Returns + ------- + lengths, n_openings : ndarray, ndarray + When *intervals_only* is False. + list of ndarray + When *intervals_only* is True. + """ + periods = record.periods + args = (periods.intervals, periods.amplitudes, abs(tcrit), periods.flags) + if intervals_only: + return extract_burst_intervals(*args) + return extract_bursts(*args) diff --git a/tests/test_bursts.py b/tests/test_bursts.py new file mode 100644 index 0000000..455dc12 --- /dev/null +++ b/tests/test_bursts.py @@ -0,0 +1,212 @@ +"""Burst segmentation. + +The convention tests are the ones that pinned this behaviour in SCALCS; they +move here with the code. Each names the reading of the burst definition it +fixes, because the three implementations in the DCPROGS stack disagreed on +exactly these points until they were written down. +""" + +import numpy as np +import pytest + +from dcio.analysis.bursts import ( + bursts_from_record, + extract_burst_intervals, + extract_bursts, +) +from dcio.analysis.record import from_scn +from dcio.formats.scn import FLAG_UNUSABLE, read + + +class TestExtractBursts: + + def test_handbuilt_three_bursts(self): + """This record opens and closes with a long shut, so all three runs + are bounded by separators and all three are complete bursts. + + These ends used to be dropped unconditionally. Against the Burzomato + 2004 records that lost two bursts from every file -- at 30 uM, two of + six.""" + # sep A.o A.s A.o sep B.o B.s B.o B.s B.o sep C.o sep + t = np.array([1, 0.1, 0.01, 0.1, 1, 0.2, 0.02, 0.3, 0.01, 0.1, 1, 0.5, 1.0]) + a = np.array([0, 5, 0, 5, 0, 5, 0, 5, 0, 5, 0, 5, 0.0]) + lengths, nops = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 3 + assert lengths[0] == pytest.approx(0.21) # 0.1+0.01+0.1 + assert lengths[1] == pytest.approx(0.63) # 0.2+0.02+0.3+0.01+0.1 + assert lengths[2] == pytest.approx(0.5) + assert list(nops) == [2, 3, 1] + + def test_keeps_ends_bounded_by_separators(self): + """Both runs are bounded by long shuts, so both are complete.""" + t = np.array([1.0, 0.1, 1.0, 0.2, 1.0]) + a = np.array([0.0, 5.0, 0.0, 5.0, 0.0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 2 + + def test_first_opening_starts_a_burst(self): + """No gap is required before the first burst: a record beginning on a + defined opening begins on a burst. EKDIST states this convention in + Bursts.slice_bursts and dcpyps followed it.""" + t = np.array([0.1, 1.0, 0.2, 1.0, 0.3]) + a = np.array([5.0, 0.0, 5.0, 0.0, 5.0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 3 + assert lengths == pytest.approx([0.1, 0.2, 0.3]) + + def test_no_separators_at_all(self): + """With nothing to cut on, the record is a single burst.""" + t = np.array([0.1, 0.01, 0.1]) + a = np.array([5.0, 0.0, 5.0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 1 + assert lengths[0] == pytest.approx(0.21) + + def test_shut_exactly_tcrit_is_within_burst(self): + """tcrit is the time such that gaps *longer* than it end a burst, so a + shut interval exactly equal to it does not. No real record exercises + this -- none of the four Burzomato 2004 files contains a shut time + equal to its tcrit -- so the convention is pinned here instead. EKDIST + makes the same choice.""" + t = np.array([0.1, 0.5, 0.2]) # middle shut is exactly tcrit + a = np.array([5.0, 0.0, 5.0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 1 # one burst, not two + assert lengths[0] == pytest.approx(0.8) + + lengths, _ = extract_bursts(t, a, tcrit=0.4999) + assert len(lengths) == 2 # just above, and it cuts + + def test_unusable_interval_ends_a_burst(self): + """Time-course fitting leaves the last interval with no defined + length, flagged unusable. It still ends the burst before it, and its + meaningless duration is never compared with tcrit.""" + # o s o unusable shut + t = np.array([0.1, 0.01, 0.1, 0.00005]) + a = np.array([5.0, 0.0, 5.0, 0.0]) + flags = np.array([0, 0, 0, FLAG_UNUSABLE]) + lengths, _ = extract_bursts(t, a, tcrit=0.5, flags=flags) + assert len(lengths) == 1 + assert lengths[0] == pytest.approx(0.21) + + def test_unusable_shut_cuts_mid_record(self): + """An unusable interval in the middle ends one burst and the next + begins after it, however short its nominal duration.""" + t = np.array([0.1, 0.00004, 0.2]) + a = np.array([5.0, 0.0, 5.0]) + flags = np.array([0, FLAG_UNUSABLE, 0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5, flags=flags) + assert len(lengths) == 2 + assert lengths == pytest.approx([0.1, 0.2]) + + def test_leading_bad_interval_discarded_not_the_burst(self): + """A bad leading interval is dropped, but the first defined opening + after it still starts a burst.""" + t = np.array([0.05, 0.1, 1.0, 0.2]) + a = np.array([5.0, 5.0, 0.0, 5.0]) + flags = np.array([FLAG_UNUSABLE, 0, 0, 0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5, flags=flags) + assert len(lengths) == 2 + assert lengths == pytest.approx([0.1, 0.2]) + + def test_empty_record(self): + lengths, nops = extract_bursts(np.array([]), np.array([]), tcrit=0.5) + assert len(lengths) == 0 + assert len(nops) == 0 + + def test_no_openings_at_all(self): + t = np.array([1.0, 2.0, 3.0]) + a = np.array([0.0, 0.0, 0.0]) + lengths, _ = extract_bursts(t, a, tcrit=0.5) + assert len(lengths) == 0 + + +class TestExtractBurstIntervals: + + def test_same_segmentation_as_extract_bursts(self): + t = np.array([1, 0.1, 0.01, 0.1, 1, 0.2, 0.02, 0.3, 1.0]) + a = np.array([0, 5, 0, 5, 0, 5, 0, 5, 0.0]) + lengths, nops = extract_bursts(t, a, tcrit=0.5) + seqs = extract_burst_intervals(t, a, tcrit=0.5) + assert len(seqs) == len(lengths) + for seq, length in zip(seqs, lengths): + assert seq.sum() == pytest.approx(length) + + def test_every_burst_has_odd_length(self): + """Each burst starts and ends on an opening, so its interval sequence + alternates open/shut/.../open. The missed-events likelihood requires + this -- an even-length sequence would end on a shut.""" + t = np.array([1, 0.1, 0.01, 0.1, 1, 0.2, 0.02, 0.3, 1.0]) + a = np.array([0, 5, 0, 5, 0, 5, 0, 5, 0.0]) + for seq in extract_burst_intervals(t, a, tcrit=0.5): + assert len(seq) % 2 == 1 + + +class TestBurstsFromRecord: + """The record-level entry point, on the shipped experimental example.""" + + @pytest.fixture() + def record(self, legacy_experimental_scn): + return from_scn(read(legacy_experimental_scn), tres=25e-6) + + def test_golden_burst_count(self, record): + """Pins the segmentation of the shipped experimental record. + + Any change to impose_resolution or to the burst convention moves these + numbers; that is the point of having them. Measured at tres = 25 us, + tcrit = 4 ms on examples/scn/glyr_experimental.scn.""" + assert record.n_resolved == 3232 + assert record.periods.n_open == 731 + lengths, nops = bursts_from_record(record, tcrit=4e-3) + assert len(lengths) == 63 + assert nops.sum() == 731 + assert lengths.mean() == pytest.approx(98.629e-3, abs=1e-6) + + def test_tcrit_sign_is_ignored(self, record): + """HJCFIT reads a negative tcrit as a flag selecting equilibrium + vectors over CHS vectors; the magnitude is the critical time. The same + number is passed to both, so segmentation must take the magnitude.""" + pos, _ = bursts_from_record(record, tcrit=4e-3) + neg, _ = bursts_from_record(record, tcrit=-4e-3) + np.testing.assert_allclose(pos, neg) + + def test_segments_periods_not_resolved_intervals(self, record): + """impose_resolution emits a fresh open interval at every change of + fitted amplitude, so a record with sub-conductance levels does not + alternate open/shut -- this one has 1770 open-open adjacencies in 3232 + resolved intervals. Segmenting periods instead gives the same bursts + and the same lengths, but counts an opening once rather than once per + sub-level, and yields the alternating sequences the missed-events + likelihood requires.""" + ra = record.resolved_amplitudes != 0.0 + assert int((ra[:-1] & ra[1:]).sum()) == 1770 # the record does have them + + by_period, nops_period = bursts_from_record(record, tcrit=4e-3) + by_interval, nops_interval = extract_bursts( + record.resolved_intervals, record.resolved_amplitudes, 4e-3, + flags=record.resolved_flags, + ) + assert len(by_period) == len(by_interval) # same bursts + np.testing.assert_allclose(by_period, by_interval) # same lengths + assert nops_period.sum() == 731 # openings, merged + assert nops_interval.sum() == 2501 # sub-levels counted + + def test_flags_are_supplied_automatically(self, record): + """The record-level call must pass the flags through; omitting them + treats an experimental record as if time-course fitting had left no + undefined intervals in it.""" + p = record.periods + auto, _ = bursts_from_record(record, tcrit=4e-3) + explicit, _ = extract_bursts(p.intervals, p.amplitudes, 4e-3, flags=p.flags) + np.testing.assert_allclose(auto, explicit) + + def test_intervals_only_agrees_on_lengths(self, record): + lengths, _ = bursts_from_record(record, tcrit=4e-3) + seqs = bursts_from_record(record, tcrit=4e-3, intervals_only=True) + assert len(seqs) == len(lengths) + np.testing.assert_allclose([s.sum() for s in seqs], lengths) + + def test_every_burst_starts_and_ends_open(self, record): + for seq in bursts_from_record(record, tcrit=4e-3, intervals_only=True): + assert len(seq) % 2 == 1 + assert len(seq) >= 1 From a0c8585ecf48a794cea6f2aa0735c54309eab208 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 1 Sep 2026 13:46:26 +0100 Subject: [PATCH 2/2] Check the burst API in the wheel job too Deferred from #2, which was based on a master that did not have it. Co-Authored-By: Claude Opus 5 --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 86a4034..9abcb4b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -140,6 +140,7 @@ jobs: from dcio.analysis import ( Periods, SingleChannelRecord, from_scn, impose_resolution, set_periods, + extract_bursts, extract_burst_intervals, bursts_from_record, ) print("installed wheel imports formats and analysis") PY