From f2c4f9fa4cee0b6273fc812f43cda32cc44e6121 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 1 Sep 2026 14:01:32 +0100 Subject: [PATCH] Add log-binned histogram arithmetic to dcio.analysis M3 of the record-layer plan: the binning maths moves down, the figures stay where they are. EKDIST draws for a person, HJCFIT draws inside a notebook and SCALCS feeds a Qt canvas; those are legitimately different figures over the same bins, and only the bins were ever duplicated. Comparing the two existing copies turned up a defect rather than a duplicate. The bin edges agree exactly -- same start, same ratio -- but the upper limit does not. EKDIST computes it as exp(ceil(log(max))), which rounds up to the next power of e, while the whole scheme is bins per decade. That limit can land below the longest interval, and numpy.histogram then drops the tail of the distribution silently: across 400 randomly drawn exponential samples it happened in 119 of them. HJCFIT's 10 ** ceil(log10(max)) is right, and is what this implements. Also splits out what was tangled together. The two "scale factor" helpers are not copies of each other: EKDIST's normalises a fitted exponential mixture, HJCFIT's renormalises an ideal pdf onto resolved intervals from the Q matrix survival. Different quantities, both staying put. Exports bins_per_decade, log_bin_edges, log_bin_histogram and staircase -- counts, edges and step coordinates, no matplotlib anywhere. 29 new tests, including the truncation regression over 200 random samples. 262 passing. Co-Authored-By: Claude Opus 5 --- .github/workflows/ci.yml | 1 + README.md | 23 +++++ dcio/analysis/__init__.py | 10 +++ dcio/analysis/histogram.py | 171 +++++++++++++++++++++++++++++++++++++ tests/test_histogram.py | 131 ++++++++++++++++++++++++++++ 5 files changed, 336 insertions(+) create mode 100644 dcio/analysis/histogram.py create mode 100644 tests/test_histogram.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9abcb4b..1e3e5da 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -141,6 +141,7 @@ jobs: Periods, SingleChannelRecord, from_scn, impose_resolution, set_periods, extract_bursts, extract_burst_intervals, bursts_from_record, + log_bin_edges, log_bin_histogram, staircase, ) print("installed wheel imports formats and analysis") PY diff --git a/README.md b/README.md index d0d70e9..ddb4e79 100644 --- a/README.md +++ b/README.md @@ -103,6 +103,29 @@ 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. +### Dwell-time histograms + +`dcio.analysis.histogram` carries the log-binning arithmetic -- bin counts, bin +edges, and staircase coordinates -- and nothing else. Drawing stays with the +caller, because EKDIST, HJCFIT and SCALCS legitimately want different figures +over the same bins. + +```python +import numpy as np +from dcio.analysis.histogram import log_bin_histogram, staircase + +counts, edges, nbdec = log_bin_histogram(record.periods.open_intervals, tres=25e-6) +x, y = staircase(edges, counts) +ax.semilogx(x, np.sqrt(y)) # Sigworth-Sine square-root ordinate +``` + +Bins start at the resolution and each is `10 ** (1/nbdec)` wider than the last, +with `nbdec` chosen from the sample size (5 / 8 / 10 / 12). The last edge is +rounded up to a whole decade, so no interval falls outside the bins -- writing +that round-up with a natural log instead, as earlier code in this stack did, +gives a power of e and lets `numpy.histogram` drop the tail of the distribution +without saying so. + ## Running tests ```bash diff --git a/dcio/analysis/__init__.py b/dcio/analysis/__init__.py index a98e09a..c8caf2f 100644 --- a/dcio/analysis/__init__.py +++ b/dcio/analysis/__init__.py @@ -1,5 +1,11 @@ """Single-channel record analysis.""" +from dcio.analysis.histogram import ( + bins_per_decade, + log_bin_edges, + log_bin_histogram, + staircase, +) from dcio.analysis.bursts import ( bursts_from_record, extract_burst_intervals, @@ -22,4 +28,8 @@ "extract_bursts", "extract_burst_intervals", "bursts_from_record", + "bins_per_decade", + "log_bin_edges", + "log_bin_histogram", + "staircase", ] diff --git a/dcio/analysis/histogram.py b/dcio/analysis/histogram.py new file mode 100644 index 0000000..52d2655 --- /dev/null +++ b/dcio/analysis/histogram.py @@ -0,0 +1,171 @@ +"""Log-binned dwell-time histograms: the arithmetic, not the figure. + +Dwell times span decades, so they are histogrammed on a logarithmic time axis: +each bin is a fixed ratio wider than the last, and the ordinate is plotted as +the square root of the count (Sigworth & Sine 1987), which makes an exponential +component a recognisable peak rather than a featureless decay. + +This module carries only the numbers -- bin counts, bin edges, and the step +coordinates a staircase plot needs. Drawing stays with the caller: EKDIST +draws for a person, HJCFIT draws inside a notebook, SCALCS feeds a Qt canvas, +and those are legitimately different figures over the same bins. + + from dcio.analysis.histogram import log_bin_histogram, staircase + + counts, edges, nbdec = log_bin_histogram(open_periods, tres=25e-6) + x, y = staircase(edges, counts) + ax.semilogx(x, np.sqrt(y)) +""" + +from __future__ import annotations + +import numpy as np + +__all__ = [ + "bins_per_decade", + "log_bin_edges", + "log_bin_histogram", + "staircase", +] + + +def bins_per_decade(n): + """Bins per decade for a sample of *n* intervals. + + The DCprogs convention: widen the bins for smaller samples so the counts + in them stay usable. + + Parameters + ---------- + n : int + Number of intervals in the sample. + + Returns + ------- + int + 5 at 300 intervals or fewer, 8 to 1000, 10 to 3000, 12 above. + """ + if n <= 300: + return 5 + if n <= 1000: + return 8 + if n <= 3000: + return 10 + return 12 + + +def log_bin_edges(intervals, tres, nbdec=None): + """Geometric bin edges for a dwell-time histogram. + + Bins start at the resolution and each is ``10 ** (1 / nbdec)`` times wider + than the last, so *nbdec* of them span a decade. The last edge is at or + above the longest interval, rounded up to a whole decade. + + That rounding is the part worth stating. Earlier implementations in this + stack wrote the decade round-up as ``exp(ceil(log(max)))``, which rounds up + to the next power of *e* rather than of ten. Because it can land below the + longest interval, ``numpy.histogram`` then drops the tail of the + distribution without saying so -- in randomly drawn exponential samples + that happens about three times in ten. + + Parameters + ---------- + intervals : array_like + Observed dwell times, in seconds. Must contain a positive value. + tres : float + Resolution, in seconds. The histogram starts here. + nbdec : int, optional + Bins per decade. Chosen from the sample size when omitted; see + :func:`bins_per_decade`. + + Returns + ------- + edges : ndarray + Bin edges, in seconds, ascending from *tres*. + nbdec : int + The number of bins per decade actually used. + + Raises + ------ + ValueError + If *tres* is not positive, or no interval is positive. + """ + intervals = np.asarray(intervals, dtype=float) + if tres <= 0.0: + raise ValueError(f"tres must be positive, got {tres!r}") + + positive = intervals[intervals > 0.0] + if positive.size == 0: + raise ValueError("no positive interval to histogram") + + if nbdec is None: + nbdec = bins_per_decade(len(intervals)) + + ratio = 10.0 ** (1.0 / nbdec) + tmax = 10.0 ** np.ceil(np.log10(positive.max())) + nbin = int(np.log(tmax / tres) / np.log(ratio)) + 1 + return tres * ratio ** np.arange(nbin + 1), nbdec + + +def log_bin_histogram(intervals, tres, nbdec=None): + """Bin dwell times logarithmically. + + Parameters + ---------- + intervals : array_like + Observed dwell times, in seconds. + tres : float + Resolution, in seconds. + nbdec : int, optional + Bins per decade; see :func:`bins_per_decade`. + + Returns + ------- + counts : ndarray of int + Number of intervals in each bin. + edges : ndarray + Bin edges, in seconds. + nbdec : int + Bins per decade used. + + Notes + ----- + Intervals shorter than *tres* fall below the first edge and are not + counted; by construction of the record they should not exist. No interval + falls above the last edge -- see :func:`log_bin_edges`. + """ + intervals = np.asarray(intervals, dtype=float) + edges, nbdec = log_bin_edges(intervals, tres, nbdec) + counts, _ = np.histogram(intervals, bins=edges) + return counts, edges, nbdec + + +def staircase(edges, counts): + """Step-plot coordinates for a binned histogram. + + Turns *n* counts and *n + 1* edges into the vertex sequence a line plot + needs to render as a staircase, closed to zero at both ends. + + Parameters + ---------- + edges : array_like + Bin edges, length ``n + 1``. + counts : array_like + Bin counts, length ``n``. + + Returns + ------- + x, y : ndarray + Vertices, each of length ``2 * (n + 1)``. Plot them directly; take + ``sqrt(y)`` for the conventional square-root ordinate. + """ + edges = np.asarray(edges, dtype=float) + counts = np.asarray(counts) + if len(edges) != len(counts) + 1: + raise ValueError( + f"expected {len(counts) + 1} edges for {len(counts)} counts, " + f"got {len(edges)}" + ) + x = np.repeat(edges, 2) + y = np.concatenate(([0], np.repeat(counts, 2), [0])) + return x, y diff --git a/tests/test_histogram.py b/tests/test_histogram.py new file mode 100644 index 0000000..b79d776 --- /dev/null +++ b/tests/test_histogram.py @@ -0,0 +1,131 @@ +"""Log-binned dwell-time histogram arithmetic.""" + +import math + +import numpy as np +import pytest + +from dcio.analysis.histogram import ( + bins_per_decade, + log_bin_edges, + log_bin_histogram, + staircase, +) + + +class TestBinsPerDecade: + + @pytest.mark.parametrize("n,expected", [ + (1, 5), (300, 5), (301, 8), (1000, 8), + (1001, 10), (3000, 10), (3001, 12), (100000, 12), + ]) + def test_boundaries(self, n, expected): + assert bins_per_decade(n) == expected + + +class TestLogBinEdges: + + def test_ratio_is_constant(self): + edges, nbdec = log_bin_edges(np.full(500, 0.01), tres=1e-5) + ratios = edges[1:] / edges[:-1] + np.testing.assert_allclose(ratios, ratios[0]) + assert ratios[0] == pytest.approx(10.0 ** (1.0 / nbdec)) + + def test_nbdec_bins_span_a_decade(self): + edges, nbdec = log_bin_edges(np.full(500, 0.01), tres=1e-5) + assert edges[nbdec] / edges[0] == pytest.approx(10.0) + + def test_starts_at_tres(self): + edges, _ = log_bin_edges(np.full(100, 0.01), tres=2.5e-5) + assert edges[0] == pytest.approx(2.5e-5) + + def test_nbdec_override_is_honoured(self): + edges, nbdec = log_bin_edges(np.full(5000, 0.01), tres=1e-5, nbdec=6) + assert nbdec == 6 + assert edges[6] / edges[0] == pytest.approx(10.0) + + @pytest.mark.parametrize("largest", [5.0, 0.05, 0.5, 30.0, 0.0403, 1.0]) + def test_last_edge_reaches_a_whole_decade(self, largest): + """The upper limit rounds up to a power of ten, not a power of e. + + Writing that round-up as exp(ceil(log(max))) -- as earlier code in this + stack did -- gives a power of e, which can fall below the longest + interval.""" + x = np.array([1e-4, largest]) + edges, _ = log_bin_edges(x, tres=2.5e-5) + decade = 10.0 ** math.ceil(math.log10(largest)) + assert edges[-1] >= decade + + def test_no_interval_falls_outside_the_bins(self): + """The regression that motivated moving this down. + + EKDIST's upper limit could land below max(X); numpy.histogram then + drops those intervals with no warning. Across randomly drawn + exponential samples it happened in about 3 of every 10.""" + rng = np.random.default_rng(0) + for _ in range(200): + n = int(rng.integers(50, 6000)) + x = rng.exponential(0.005, n) + 2.5e-5 + edges, _ = log_bin_edges(x, tres=2.5e-5) + assert x.max() <= edges[-1] + assert x.min() >= edges[0] + + def test_rejects_nonpositive_tres(self): + with pytest.raises(ValueError, match="tres must be positive"): + log_bin_edges(np.array([0.01]), tres=0.0) + + def test_rejects_empty_sample(self): + with pytest.raises(ValueError, match="no positive interval"): + log_bin_edges(np.array([]), tres=1e-5) + + def test_ignores_negative_sentinels(self): + """impose_resolution marks an unfinished final opening with -1.""" + x = np.array([0.001, 0.002, -1.0]) + edges, _ = log_bin_edges(x, tres=1e-5) + assert edges[-1] >= 0.002 + + +class TestLogBinHistogram: + + def test_counts_every_interval(self): + rng = np.random.default_rng(1) + x = rng.exponential(0.005, 2000) + 2.5e-5 + counts, edges, _ = log_bin_histogram(x, tres=2.5e-5) + assert counts.sum() == len(x) + + def test_edges_are_one_longer_than_counts(self): + x = np.full(400, 0.01) + counts, edges, _ = log_bin_histogram(x, tres=1e-5) + assert len(edges) == len(counts) + 1 + + def test_all_equal_intervals_land_in_one_bin(self): + x = np.full(400, 0.01) + counts, _, _ = log_bin_histogram(x, tres=1e-5) + assert (counts > 0).sum() == 1 + assert counts.max() == 400 + + +class TestStaircase: + + def test_vertex_count(self): + counts = np.array([1, 2, 3]) + edges = np.array([1.0, 2.0, 3.0, 4.0]) + x, y = staircase(edges, counts) + assert len(x) == len(y) == 2 * len(edges) + + def test_closes_to_zero_at_both_ends(self): + x, y = staircase(np.array([1.0, 2.0, 3.0]), np.array([4, 5])) + assert y[0] == 0 + assert y[-1] == 0 + + def test_each_bar_is_flat_across_its_bin(self): + counts = np.array([7, 3]) + edges = np.array([1.0, 2.0, 4.0]) + x, y = staircase(edges, counts) + # vertices: (1,0) (1,7) (2,7) (2,3) (4,3) (4,0) + assert list(y) == [0, 7, 7, 3, 3, 0] + assert list(x) == [1.0, 1.0, 2.0, 2.0, 4.0, 4.0] + + def test_length_mismatch_raises(self): + with pytest.raises(ValueError, match="expected 3 edges"): + staircase(np.array([1.0, 2.0]), np.array([1, 2]))