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