From c9d2a74624ffdb7ff80299dddf9b93aabb37f64e Mon Sep 17 00:00:00 2001 From: Remis Date: Wed, 2 Sep 2026 11:48:46 +0100 Subject: [PATCH 1/2] Measure the correlation between an opening and its neighbouring shut times ekdist.adjacent measures from a record what scalcs.sccurves and scalcs.scalcslib give in theory, so that one can be superimposed on the other. It is what Figs 6D and 6E of Colquhoun, Hatton & Hawkes (2003) J Physiol 547:699-728 display, and what a fit to a single-channel record is checked against once the rate constants are in hand. pairs() puts each opening beside a neighbouring shut time, conditional_mean_open_time() gives the mean apparent open time for openings whose adjacent shut time falls in each of a list of ranges, and dependency() gives the observed dependency of Magleby & Song (1992) on a log-log grid. Every function takes side='next', 'previous' or 'adjacent', because the theory distinguishes them: scalcs returns the mean open time given the previous gap and given the next gap as two separate curves, and its HJC_dependency is the joint density of an open time and the shut time that FOLLOWS it. Passing the wrong one raises nothing; it quietly compares a measurement with the wrong curve. The dependency's marginals are taken over every pair, not over the pairs that fall inside the grid. Using the row and column sums of the joint histogram renormalises to whatever rectangle the grid covers, and a shut-time grid never spans the distribution, because the gaps between activations run to seconds. In a record of 200000 intervals only 7316 of 49236 pairs were inside the grid, and the mistake turned a dependency of zero into values of two and three. There is a regression test for it. 15 tests, on constructions whose answer is known by hand or by symmetry rather than by comparison with another library: independent times must give zero dependency, correlated ones a positive diagonal, and the conditional means are checked against sums done by hand. Agreement with the Q-matrix theory is checked separately, where a record can be simulated from a mechanism: on 150000 intervals of the CHH 2003 AChR mechanism the observed dependency agrees over 21 cells with a mean difference of -0.004, against a theory spanning -0.90 to +2.58. Co-Authored-By: Claude Opus 5 --- ekdist/adjacent.py | 198 +++++++++++++++++++++++++++++++++++++++++ tests/test_adjacent.py | 181 +++++++++++++++++++++++++++++++++++++ 2 files changed, 379 insertions(+) create mode 100644 ekdist/adjacent.py create mode 100644 tests/test_adjacent.py diff --git a/ekdist/adjacent.py b/ekdist/adjacent.py new file mode 100644 index 0000000..f7d107d --- /dev/null +++ b/ekdist/adjacent.py @@ -0,0 +1,198 @@ +"""Statistics of open times and the shut times next to them. + +Correlations between an opening and its neighbouring shuttings are a property +of the mechanism, not of the fitting, and they are what the conditional +displays of a single-channel record are for. This module measures them **from a +record**; the theoretical counterparts live in :mod:`scalcs.sccurves` and +:mod:`scalcs.scalcslib`, and the point of having both is to superimpose one on +the other. + +Two measurements, matching Figs 6D and 6E of Colquhoun, Hatton & Hawkes (2003): + +:func:`conditional_mean_open_time` + the mean apparent open time for openings whose adjacent shut time falls in + each of a list of ranges -- the diamonds of Fig. 6D. +:func:`dependency` + the observed dependency of Magleby & Song (1992) on a log-log grid -- + Fig. 6E. + +**Which neighbour.** An opening has a shut time on each side, and the theory +distinguishes them: ``scalcs`` returns the mean open time given the *previous* +gap and given the *next* gap as two different curves, and its ``HJC_dependency`` +is the joint density of an open time and the shut time that **follows** it. So +every function here takes ``side``: + +``'next'`` + the shut time that follows the opening. Compare with ``HJC_dependency``, + and with the second return value of + ``HJC_adjacent_mean_open_to_shut_time_pdf``. +``'previous'`` + the shut time before the opening. Compare with the first return value. +``'adjacent'`` + both, pooled, so an opening with a shut time on each side is counted twice. + This is what a legend means by "adjacent". + +Getting this wrong does not raise anything: it quietly compares a measurement +with the wrong curve. +""" + +from __future__ import annotations + +import numpy as np + +SIDES = ("next", "previous", "adjacent") + + +def pairs(tints, ampls, side="next"): + """Pair each opening with a neighbouring shut time. + + Parameters + ---------- + tints, ampls : array_like + An alternating record: interval durations, and amplitudes that are zero + for a shut interval and non-zero for an open one. This is what + ``scalcs.scsim.impose_resolution`` returns. + side : {'next', 'previous', 'adjacent'} + + Returns + ------- + topen, tshut : ndarray + Equal length. Openings at either end of the record with no neighbour on + the required side are dropped. + """ + if side not in SIDES: + raise ValueError(f"side must be one of {SIDES}, not {side!r}") + tints = np.asarray(tints, float) + ampls = np.asarray(ampls, float) + if tints.shape != ampls.shape: + raise ValueError("tints and ampls must have the same length") + + is_open = ampls != 0.0 + topen, tshut = [], [] + if side in ("next", "adjacent"): + take = is_open[:-1] & ~is_open[1:] + topen.append(tints[:-1][take]) + tshut.append(tints[1:][take]) + if side in ("previous", "adjacent"): + take = is_open[1:] & ~is_open[:-1] + topen.append(tints[1:][take]) + tshut.append(tints[:-1][take]) + return np.concatenate(topen), np.concatenate(tshut) + + +def conditional_mean_open_time(topen, tshut, ranges): + """Mean apparent open time for each range of adjacent shut time. + + The diamonds and their bars in Fig. 6D. Ranges must be given, not chosen + here, because they have to be wide enough to hold enough observations and + that is a judgement about the data. + + Parameters + ---------- + topen, tshut : array_like + Paired, as returned by :func:`pairs`. + ranges : sequence of (float, float) + Half-open shut-time ranges ``[lo, hi)``; ``hi`` may be ``inf``. + + Returns + ------- + dict of ndarray + ``n`` (openings in the range), ``shut_mean`` (mean of the shut times + that defined it, which is where the paper plots the point), + ``open_mean``, ``open_sd``, ``open_sem``. Empty ranges give ``nan`` + rather than being dropped, so the result lines up with ``ranges``. + """ + topen = np.asarray(topen, float) + tshut = np.asarray(tshut, float) + n = np.zeros(len(ranges), int) + shut_mean = np.full(len(ranges), np.nan) + open_mean = np.full(len(ranges), np.nan) + open_sd = np.full(len(ranges), np.nan) + + for i, (lo, hi) in enumerate(ranges): + take = (tshut >= lo) & (tshut < hi) + n[i] = int(take.sum()) + if n[i] == 0: + continue + shut_mean[i] = tshut[take].mean() + open_mean[i] = topen[take].mean() + if n[i] > 1: + open_sd[i] = topen[take].std(ddof=1) + + with np.errstate(invalid="ignore", divide="ignore"): + open_sem = open_sd / np.sqrt(n) + return dict(n=n, shut_mean=shut_mean, open_mean=open_mean, + open_sd=open_sd, open_sem=open_sem) + + +def log_edges(values, tres=None, per_decade=8): + """Log-spaced bin edges spanning ``values``, for :func:`dependency`.""" + values = np.asarray(values, float) + lo = tres if tres else values[values > 0].min() + hi = values.max() + n = max(int(np.ceil(np.log10(hi / lo) * per_decade)), 2) + return np.logspace(np.log10(lo), np.log10(hi), n + 1) + + +def dependency(topen, tshut, open_edges, shut_edges, min_count=10): + """The observed dependency, on a grid of open and shut time. + + Dependency is defined (eqn (8) of the paper, after Magleby & Song 1992) as + + d(to, ts) = f(to, ts) / (fo(to) fs(ts)) - 1, + + zero where the two are independent, positive where openings of that length + occur next to shuttings of that length more often than chance. Binned, the + ratio of densities is the ratio of the joint probability of a cell to the + product of its marginals, so with ``N`` observations in total, + + d = N_ij N / (N_i N_j) - 1. + + Parameters + ---------- + topen, tshut : array_like + Paired, as returned by :func:`pairs`. + open_edges, shut_edges : array_like + Bin edges; use :func:`log_edges`, since dwell times span decades. + min_count : int + Cells holding fewer than this are returned as ``nan``. A dependency + estimated from three intervals is noise, and the published figures + leave such regions blank. + + Returns + ------- + dict + ``dependency`` (2-D, open by shut, ``nan`` where sparse), ``counts`` + (the joint histogram), ``open_centres`` and ``shut_centres`` + (geometric), ``n`` (pairs in all, which is what the marginals are taken + over) and ``inside`` (pairs that fell within the grid). + """ + topen = np.asarray(topen, float) + tshut = np.asarray(tshut, float) + open_edges = np.asarray(open_edges, float) + shut_edges = np.asarray(shut_edges, float) + + counts, _, _ = np.histogram2d(topen, tshut, bins=[open_edges, shut_edges]) + if counts.sum() == 0: + raise ValueError("no pairs fell inside the grid") + + # The marginals must come from EVERY pair, not only from those inside the + # grid, because the definition is a ratio of the joint density to the + # product of the unconditional marginals. Taking the row and column sums of + # the joint histogram instead silently renormalises to whatever rectangle + # the grid happens to cover, which is wrong whenever the grid does not span + # the whole distribution -- and a shut-time grid never does, since the long + # shut times between activations run to seconds. + total = len(topen) + row = np.histogram(topen, bins=open_edges)[0][:, None].astype(float) + col = np.histogram(tshut, bins=shut_edges)[0][None, :].astype(float) + + with np.errstate(invalid="ignore", divide="ignore"): + dep = counts * total / (row * col) - 1.0 + dep[counts < min_count] = np.nan + + centres = lambda e: np.sqrt(e[:-1] * e[1:]) # noqa: E731 (geometric) + return dict(dependency=dep, counts=counts, n=total, + inside=int(counts.sum()), + open_centres=centres(open_edges), + shut_centres=centres(shut_edges)) diff --git a/tests/test_adjacent.py b/tests/test_adjacent.py new file mode 100644 index 0000000..354463b --- /dev/null +++ b/tests/test_adjacent.py @@ -0,0 +1,181 @@ +"""Statistics of open times and the shut times next to them. + +Self-contained: the references here are constructions whose answer is known by +hand or by symmetry, not values taken from another library. The agreement with +the Q-matrix theory is checked separately, where a record can be simulated from +a mechanism. +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from ekdist import adjacent + + +# a five-interval record: open shut open shut open +TINTS = np.array([1.0, 10.0, 2.0, 20.0, 3.0]) +AMPLS = np.array([5.0, 0.0, 5.0, 0.0, 5.0]) + + +# ------------------------------------------------------------------- pairs + +def test_pairs_next(): + """Each opening with the shut time after it; the last has none.""" + topen, tshut = adjacent.pairs(TINTS, AMPLS, side="next") + npt.assert_array_equal(topen, [1.0, 2.0]) + npt.assert_array_equal(tshut, [10.0, 20.0]) + + +def test_pairs_previous(): + """Each opening with the shut time before it; the first has none.""" + topen, tshut = adjacent.pairs(TINTS, AMPLS, side="previous") + npt.assert_array_equal(topen, [2.0, 3.0]) + npt.assert_array_equal(tshut, [10.0, 20.0]) + + +def test_pairs_adjacent_counts_the_middle_opening_twice(): + topen, tshut = adjacent.pairs(TINTS, AMPLS, side="adjacent") + assert len(topen) == 4 + assert sorted(topen.tolist()) == [1.0, 2.0, 2.0, 3.0] + + +def test_pairs_record_starting_shut(): + tints = np.array([10.0, 1.0, 20.0]) + ampls = np.array([0.0, 5.0, 0.0]) + npt.assert_array_equal(adjacent.pairs(tints, ampls, "next")[1], [20.0]) + npt.assert_array_equal(adjacent.pairs(tints, ampls, "previous")[1], [10.0]) + + +def test_pairs_rejects_nonsense(): + with pytest.raises(ValueError, match="side must be"): + adjacent.pairs(TINTS, AMPLS, side="sideways") + with pytest.raises(ValueError, match="same length"): + adjacent.pairs(TINTS, AMPLS[:-1]) + + +# ------------------------------------------- conditional mean open time + +def test_conditional_mean_open_time(): + topen = np.array([1.0, 3.0, 10.0, 20.0, 100.0]) + tshut = np.array([0.5, 0.7, 5.0, 6.0, 500.0]) + out = adjacent.conditional_mean_open_time( + topen, tshut, [(0.0, 1.0), (1.0, 10.0), (10.0, 100.0)]) + npt.assert_array_equal(out["n"], [2, 2, 0]) + npt.assert_allclose(out["open_mean"][:2], [2.0, 15.0]) + npt.assert_allclose(out["shut_mean"][:2], [0.6, 5.5]) + npt.assert_allclose(out["open_sd"][:2], + [np.std([1, 3], ddof=1), np.std([10, 20], ddof=1)]) + npt.assert_allclose(out["open_sem"][:2], out["open_sd"][:2] / np.sqrt(2)) + + +def test_empty_range_is_nan_and_keeps_its_place(): + """The result must line up with the ranges, or a plot misaligns.""" + out = adjacent.conditional_mean_open_time( + np.array([1.0]), np.array([1.0]), [(0.0, 0.5), (0.5, 2.0)]) + assert out["n"][0] == 0 + assert np.isnan(out["open_mean"][0]) and np.isnan(out["shut_mean"][0]) + npt.assert_allclose(out["open_mean"][1], 1.0) + + +def test_ranges_are_half_open(): + topen = np.array([1.0, 2.0]) + tshut = np.array([1.0, 2.0]) + out = adjacent.conditional_mean_open_time(topen, tshut, [(1.0, 2.0)]) + npt.assert_array_equal(out["n"], [1]) + npt.assert_allclose(out["open_mean"], [1.0]) + + +def test_infinite_upper_limit(): + out = adjacent.conditional_mean_open_time( + np.array([1.0, 2.0]), np.array([1.0, 1e9]), [(10.0, np.inf)]) + npt.assert_array_equal(out["n"], [1]) + npt.assert_allclose(out["open_mean"], [2.0]) + + +# ------------------------------------------------------------- dependency + +def lognormal_pair(n, rho=0.0, seed=0): + """``n`` pairs of positive times with a given correlation of their logs.""" + rng = np.random.default_rng(seed) + z1 = rng.normal(size=n) + z2 = rho * z1 + np.sqrt(1 - rho ** 2) * rng.normal(size=n) + return np.exp(z1 * 1.2 - 6.0), np.exp(z2 * 1.5 - 5.0) + + +def test_independent_times_give_zero_dependency(): + """The defining property: independence means dependency zero.""" + topen, tshut = lognormal_pair(400000, rho=0.0, seed=1) + edges_o = np.logspace(np.log10(np.quantile(topen, 0.05)), + np.log10(np.quantile(topen, 0.95)), 7) + edges_s = np.logspace(np.log10(np.quantile(tshut, 0.05)), + np.log10(np.quantile(tshut, 0.95)), 7) + out = adjacent.dependency(topen, tshut, edges_o, edges_s, min_count=50) + dep = out["dependency"][~np.isnan(out["dependency"])] + assert len(dep) > 20 + assert abs(dep.mean()) < 0.02, dep.mean() + assert np.abs(dep).max() < 0.15, np.abs(dep).max() + + +def test_marginals_are_taken_over_every_pair_not_only_the_grid(): + """Regression: a grid covering part of the distribution must not bias it. + + Taking the marginals as the row and column sums of the joint histogram + renormalises to whatever rectangle the grid covers. For independent times + that turned a dependency of zero into values of two or three; a shut-time + grid never spans the whole distribution, because the gaps between + activations run to seconds. + """ + topen, tshut = lognormal_pair(400000, rho=0.0, seed=2) + # deliberately narrow: most pairs fall outside + edges_o = np.logspace(np.log10(np.quantile(topen, 0.30)), + np.log10(np.quantile(topen, 0.60)), 5) + edges_s = np.logspace(np.log10(np.quantile(tshut, 0.30)), + np.log10(np.quantile(tshut, 0.60)), 5) + out = adjacent.dependency(topen, tshut, edges_o, edges_s, min_count=50) + assert out["inside"] < 0.2 * out["n"], "the grid should be narrow" + dep = out["dependency"][~np.isnan(out["dependency"])] + assert np.abs(dep).max() < 0.15, np.abs(dep).max() + + +def test_correlated_times_give_positive_dependency_on_the_diagonal(): + """Both short or both long together, one of each much less than chance.""" + topen, tshut = lognormal_pair(400000, rho=0.6, seed=3) + edges_o = np.logspace(np.log10(np.quantile(topen, 0.05)), + np.log10(np.quantile(topen, 0.95)), 7) + edges_s = np.logspace(np.log10(np.quantile(tshut, 0.05)), + np.log10(np.quantile(tshut, 0.95)), 7) + dep = adjacent.dependency(topen, tshut, edges_o, edges_s, + min_count=20)["dependency"] + assert dep[0, 0] > 0.5 and dep[-1, -1] > 0.5 + # the opposite corners are strongly negative, or so nearly empty that they + # come back nan -- which says the same thing + for corner in (dep[0, -1], dep[-1, 0]): + assert np.isnan(corner) or corner < -0.5, corner + # and the effect is monotone across the grid + assert np.nanmean(np.diag(dep)) > np.nanmean(np.diag(np.fliplr(dep))) + + +def test_sparse_cells_are_nan(): + topen, tshut = lognormal_pair(2000, seed=4) + edges_o = np.logspace(np.log10(topen.min()), np.log10(topen.max()), 30) + edges_s = np.logspace(np.log10(tshut.min()), np.log10(tshut.max()), 30) + out = adjacent.dependency(topen, tshut, edges_o, edges_s, min_count=10) + dep, counts = out["dependency"], out["counts"] + assert np.isnan(dep[counts < 10]).all() + assert not np.isnan(dep[counts >= 10]).any() + + +def test_empty_grid_raises(): + with pytest.raises(ValueError, match="no pairs"): + adjacent.dependency(np.array([1.0]), np.array([1.0]), + np.array([10.0, 20.0]), np.array([10.0, 20.0])) + + +def test_log_edges(): + values = np.array([1e-5, 1e-3, 1e-1]) + edges = adjacent.log_edges(values, tres=1e-5, per_decade=2) + assert edges[0] == pytest.approx(1e-5) + assert edges[-1] == pytest.approx(1e-1) + assert len(edges) == 9 + assert np.all(np.diff(edges) > 0) From 92321ee072988b04f50be1e4adea5ff4ad88a615 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 8 Sep 2026 15:40:24 +0100 Subject: [PATCH 2/2] clean-install: enumerate the modules instead of listing them The job is titled "Import every module with nothing else installed" and exists because setup.py has to name everything the package imports -- scipy was missing from install_requires for a while and only a clean environment could show it. But it imported a hand-written list, and the list had already fallen behind: ekdist.adjacent is not in it. So the job that catches an undeclared dependency would have missed one in any module added since the list was written. pkgutil.iter_modules finds whatever is actually installed; it reports 7 where the list had 6. Co-Authored-By: Claude Opus 5 --- .github/workflows/ci.yml | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8bd9e57..b878562 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -96,13 +96,30 @@ jobs: /tmp/clean/bin/pip install . /tmp/clean/bin/pip list + # Enumerated rather than listed by hand. The hand-written list said "every + # module" and had already fallen behind -- ekdist.adjacent was not in it -- + # so a new module with an undeclared dependency would have gone straight + # through the job meant to catch exactly that. pkgutil finds whatever is + # actually installed. - name: Import every module with nothing else installed env: MPLBACKEND: Agg run: | cd /tmp /tmp/clean/bin/python - <<'PY' - from ekdist import ekplot, ekrecord, ekscn, eklib, exponentials, errors - from ekdist.ekplot import prepare_xlog_hist - print("all ekdist modules import from a clean install") + import importlib + import pkgutil + + import ekdist + + names = sorted(m.name for m in pkgutil.iter_modules(ekdist.__path__)) + if not names: + raise SystemExit("no ekdist submodules found; the install is wrong") + for name in names: + importlib.import_module(f"ekdist.{name}") + print(f" ok ekdist.{name}") + + # one attribute, so that a module which imports but is empty shows up + from ekdist.ekplot import prepare_xlog_hist # noqa: F401 + print(f"all {len(names)} ekdist modules import from a clean install") PY