Add log-binned histogram arithmetic (M3)

.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

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

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",
 ]

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

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