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