Build the log histogram on dcio (M3)

ekdist/ekplot.py

@@ -6,6 +6,10 @@
 from matplotlib import transforms as mtransforms
 from matplotlib import ticker
 
+from dcio.analysis.histogram import (
+    bins_per_decade, log_bin_histogram, staircase,
+)
+
 from ekdist import eklib
 from ekdist import exponentials
 
@@ -73,21 +77,31 @@ def burst_number_of_openings(nops):
 ###############################################################################
 # Dwell time histograms: x-log / y-sqrt
 def __histogram_bins_per_decade(X):
-    nbdec = 12
-    if (len(X) <= 300): nbdec = 5 
-    if (len(X) > 300) and (len(X) <= 1000): nbdec = 8
-    if (len(X) > 1000) and (len(X) <= 3000): nbdec = 10
-    return nbdec
+    return bins_per_decade(len(X))
 
 def __bin_width(X):
-    return math.exp(math.log(10.0) / float(__histogram_bins_per_decade(X)))
+    return 10.0 ** (1.0 / bins_per_decade(len(X)))
 
 def __exponential_scale_factor(X, pdf, tres):
+    """Normalise a fitted exponential mixture onto the histogram.
+
+    Not the same quantity as HJCFIT's ideal_pdf_scale_factor, which
+    renormalises an ideal pdf onto the resolved intervals from the Q matrix.
+    This one scales a mixture already fitted to these data, so it stays here.
+    """
     return (len(X) * math.log10(__bin_width(X)) * math.log(10) *
             (1 / np.sum(pdf.area * np.exp(-tres / pdf.tau))))
 
 def prepare_xlog_hist(X, tres):
-    """ Prepare x-log histogram.     
+    """ Prepare x-log histogram.
+
+    Binning is dcio.analysis.histogram; this is the staircase it produces.
+
+    The bin edges used to be built here, and the upper limit was computed as
+    exp(ceil(log(max(X)))) -- a power of e, where the whole scheme is bins per
+    decade. That limit can fall below the longest interval, and np.histogram
+    then drops the tail of the distribution silently: across 400 randomly drawn
+    exponential samples it happened in 119 of them.
 
     Parameters
     ----------
@@ -101,14 +115,9 @@ def prepare_xlog_hist(X, tres):
     xout, yout :  list of scalar
         x and y values to plot histogram.
     """
-    dx = __bin_width(X)
-    xend = math.exp(math.ceil(math.log(max(X)))) # round up maximum value
-    nbin = int(math.log(xend / tres) / math.log(dx)) # number of bins
-    my_bins = tres * np.array([dx**i for i in range(nbin+1)])
-    hist, bin_edges = np.histogram(X, bins=my_bins)
-    xout = [x for pair in zip(bin_edges, bin_edges) for x in pair]
-    yout = [0] + [y for pair in zip(hist, hist) for y in pair] + [0]
-    return xout, yout
+    counts, edges, _ = log_bin_histogram(X, tres)
+    xout, yout = staircase(edges, counts)
+    return list(xout), list(yout)
 
 def histogram_xlog_ysqrt_data(X, tres, pdf=None, tcrit=None, xlabel='Dwell times'):
     """ Plot dwell time histogram in log x and square root y. """

setup.py

@@ -12,6 +12,10 @@
       install_requires=[
           'numpy',
           'matplotlib',
+          # The shared record layer. Not on PyPI yet, so a source install
+          # needs it first:
+          #   pip install git+https://github.com/DCPROGS/dcio
+          'dcio',
           #'PyQt5',
       ],
       zip_safe=False)

tests/test_ekplot_histogram.py

@@ -0,0 +1,70 @@
+"""Log-binned histogram, now built on dcio.analysis.histogram.
+
+The binning arithmetic moved down to dcio (M3 of the record-layer plan). These
+tests pin the behaviour EKDIST exposes through prepare_xlog_hist, including the
+truncation bug that moving it fixed.
+"""
+
+import math
+
+import numpy as np
+import pytest
+
+from ekdist.ekplot import prepare_xlog_hist
+
+
+def _ek_old_edges(X, tres):
+    """The bin edges EKDIST built before the move, for comparison.
+
+    The upper limit is exp(ceil(ln(max))) -- a power of e, not of ten.
+    """
+    n = len(X)
+    nbdec = 5 if n <= 300 else 8 if n <= 1000 else 10 if n <= 3000 else 12
+    dx = math.exp(math.log(10.0) / float(nbdec))
+    xend = math.exp(math.ceil(math.log(max(X))))
+    nbin = int(math.log(xend / tres) / math.log(dx))
+    return tres * np.array([dx ** i for i in range(nbin + 1)])
+
+
+class TestPrepareXlogHist:
+
+    def test_returns_staircase_lists(self):
+        x = np.random.default_rng(0).exponential(0.005, 500) + 2.5e-5
+        xout, yout = prepare_xlog_hist(x, 2.5e-5)
+        assert isinstance(xout, list) and isinstance(yout, list)
+        assert len(xout) == len(yout)
+
+    def test_closes_to_zero_at_both_ends(self):
+        x = np.random.default_rng(1).exponential(0.005, 500) + 2.5e-5
+        _, yout = prepare_xlog_hist(x, 2.5e-5)
+        assert yout[0] == 0 and yout[-1] == 0
+
+    def test_starts_at_tres(self):
+        x = np.random.default_rng(2).exponential(0.005, 500) + 2.5e-5
+        xout, _ = prepare_xlog_hist(x, 2.5e-5)
+        assert xout[0] == pytest.approx(2.5e-5)
+
+    def test_counts_every_interval(self):
+        """The regression the move fixed.
+
+        The old upper limit could fall below max(X), and np.histogram then
+        dropped those intervals with no warning. Over 400 random exponential
+        samples that happened in 119 of them."""
+        rng = np.random.default_rng(0)
+        old_would_drop = 0
+        for _ in range(150):
+            n = int(rng.integers(50, 6000))
+            x = rng.exponential(0.005, n) + 2.5e-5
+            _, yout = prepare_xlog_hist(x, 2.5e-5)
+            # each bar contributes its count twice
+            assert sum(yout) // 2 == n
+            if x.max() > _ek_old_edges(x, 2.5e-5)[-1]:
+                old_would_drop += 1
+        assert old_would_drop > 0, "sample no longer exercises the old bug"
+
+    def test_bin_ratio_is_a_decade_root(self):
+        x = np.random.default_rng(3).exponential(0.005, 2000) + 2.5e-5
+        xout, _ = prepare_xlog_hist(x, 2.5e-5)
+        edges = np.array(xout[::2])
+        ratio = edges[1] / edges[0]
+        assert ratio == pytest.approx(10.0 ** (1.0 / 10))   # 2000 -> 10/decade