From 5ea5b4ea11ddf48a9b0016afd035b5ed7685d6a2 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 1 Sep 2026 14:14:24 +0100 Subject: [PATCH] Build the log histogram on dcio.analysis.histogram Second half of M3: the binning arithmetic now comes from the shared record layer, and ekplot keeps the drawing. This is a bug fix as well as a deduplication. The bin edges were 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 land below the longest interval, and np.histogram then drops the tail of the distribution without saying so. Across 400 randomly drawn exponential samples the top edge fell below max(X) in 119 of them. __exponential_scale_factor stays. It looked like a third copy of HJCFIT's ideal_pdf_scale_factor and is a different quantity: this one normalises an exponential mixture already fitted to the data, HJCFIT's renormalises an ideal pdf onto the resolved intervals from the Q-matrix survival. prepare_xlog_hist keeps its signature and its list return type, so callers and notebooks are unaffected. dcio is declared in install_requires. It is not on PyPI yet, so a source install needs `pip install git+https://github.com/DCPROGS/dcio` first; the comment says so. Co-Authored-By: Claude Opus 5 --- ekdist/ekplot.py | 39 +++++++++++-------- setup.py | 4 ++ tests/test_ekplot_histogram.py | 70 ++++++++++++++++++++++++++++++++++ 3 files changed, 98 insertions(+), 15 deletions(-) create mode 100644 tests/test_ekplot_histogram.py diff --git a/ekdist/ekplot.py b/ekdist/ekplot.py index 3b716a7..8962b5e 100644 --- a/ekdist/ekplot.py +++ b/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. """ diff --git a/setup.py b/setup.py index d40d11b..cad02c0 100644 --- a/setup.py +++ b/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) \ No newline at end of file diff --git a/tests/test_ekplot_histogram.py b/tests/test_ekplot_histogram.py new file mode 100644 index 0000000..a201361 --- /dev/null +++ b/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