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 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