From 26501e945a608a046b3aca861b6f84e0d1a92a35 Mon Sep 17 00:00:00 2001 From: Remis Date: Wed, 2 Sep 2026 11:35:50 +0100 Subject: [PATCH 1/2] Correct the apparent mean open and shut times exact_mean_open_shut_time was too large by exactly twice the dead time, and HJC_adjacent_mean_open_to_shut_time_pdf by once. The algebra in neither was wrong; the offset assumed for the dARSdS term was. That term evaluates to the mean apparent duration plus tres, not minus it, so the mean is recovered by subtracting tres rather than adding it. Checked against references that owe nothing to the implementation: for a two-state channel it now reproduces eqns (121) and (122) of Colquhoun & Hawkes (1995b) exactly, and for CH82 and the AChR mechanism it agrees with the mean of the asymptotic HJC pdf, tres + sum(area * tau), which reaches the same quantity by a different route through the library. At tres = 0 the apparent means become the ideal ones, as they must. The root cause is the convention of qmatlib.dARSdS, which has not been re-derived against CHS96 eqn (3.6). The correction is applied at each call site instead, where what is wanted is unambiguous. Anyone calling dARSdS directly still inherits the same trap. Adds HJC_adjacent_mean_open_to_shut_range_mean, the mean apparent open time for openings whose adjacent shut time falls in a given range -- what is superimposed on the observed conditional means in Fig. 6D of Colquhoun, Hatton & Hawkes (2003), and not the same thing as the continuous relationship, which cannot be compared with data directly because the data must be binned. It needs the apparent shut time density as a weight, so HJC_shut_time_pdf is exposed on its own rather than remaining a denominator inside the conditional mean. The finite-dead-time reference values in test_real_dtype.py move with the correction. They also move for a second, smaller reason recorded there: popen.EC50 was made to converge on the concentration rather than stopping as soon as Popen was within 0.001 of half maximal. Nothing about the real-dtype cast those tests guard has changed. Co-Authored-By: Claude Opus 5 --- scalcs/scalcslib.py | 133 +++++++++++++++++++++++++++++++++-- tests/test_apparent_means.py | 112 +++++++++++++++++++++++++++++ tests/test_real_dtype.py | 27 ++++--- 3 files changed, 256 insertions(+), 16 deletions(-) create mode 100644 tests/test_apparent_means.py diff --git a/scalcs/scalcslib.py b/scalcs/scalcslib.py index a5b6dca..aded413 100644 --- a/scalcs/scalcslib.py +++ b/scalcs/scalcslib.py @@ -701,9 +701,17 @@ def exact_mean_open_shut_time(mec, tres): DARS = qml.dARSdS(tres, mec.QAA, mec.QFF, GAF, GFA, expQFF, mec.kA, mec.kF) DFRS = qml.dARSdS(tres, mec.QFF, mec.QAA, GFA, GAF, expQAA, mec.kF, mec.kA) uF, uA = np.ones((mec.kF, 1)), np.ones((mec.kA, 1)) - # meanOpenTime = tres + phiA * DARS * QexpQF * uF - meanA = tres + np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] - meanF = tres + np.dot(phiF, np.dot(np.dot(DFRS, QexpQA), uA))[0] + # The DARS term below evaluates to the mean apparent duration PLUS tres, not + # minus it, so the mean is recovered by subtracting tres rather than adding it. + # Corrected 2026-09-02, having been wrong by 2*tres. Checked two ways: for a + # two-state channel it now reproduces eqns (121) and (122) of Colquhoun & + # Hawkes (1995b) exactly, and for CH82 and the AChR mechanism it agrees with + # the mean of the asymptotic HJC pdf, tres + sum(area * tau), to six figures. + # The root cause is the convention of qmatlib.dARSdS, which has not been + # re-derived against CHS96 eqn (3.6); the correction is applied at each call + # site instead, where what is wanted is unambiguous. + meanA = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] - tres + meanF = np.dot(phiF, np.dot(np.dot(DFRS, QexpQA), uA))[0] - tres return meanA, meanF @@ -745,8 +753,8 @@ def exact_mean_time(tres, QAA, QFF, QAF, kA, kF, GAF, GFA): DARS = qml.dARSdS(tres, QAA, QFF, GAF, GFA, expQFF, kA, kF) uF = np.ones((kF, 1)) - # meanOpenTime = tres + phiA * DARS * QexpQF * uF - mean = tres + np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] + # see exact_mean_open_shut_time for why this subtracts tres + mean = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] - tres return mean @@ -1217,11 +1225,122 @@ def HJC_adjacent_mean_open_to_shut_time_pdf(sht, tres, Q, QAA, QAF, QFF, QFA): denom = np.dot(np.dot(phiF, eGFAt), uA)[0] nom1 = np.dot(np.dot(phiF, eGFAt), col1)[0] nom2 = np.dot(np.dot(row1, eGFAt), uA)[0] - mp.append(nom1 / denom) - mn.append(nom2 / denom) + # the same tres offset as in exact_mean_open_shut_time: these are + # conditional means of the apparent open time, and the DARS term + # carries an extra tres + mp.append(nom1 / denom - tres) + mn.append(nom2 / denom - tres) return np.array(mp), np.array(mn) +#: numpy renamed trapz to trapezoid in 2.0; keep working on both +_trapezoid = getattr(np, 'trapezoid', None) or np.trapz + + +def HJC_adjacent_mean_open_to_shut_range_mean(ranges, tres, Q, QAA, QAF, + QFF, QFA, points=64, + side='next'): + """ + Mean apparent open time for openings whose adjacent apparent shut time + falls in each of a series of ranges, allowing for missed events. + + This is what is superimposed on the observed conditional means of a + record -- the filled circles of Fig. 6D in Colquhoun, Hatton & Hawkes + (2003) J Physiol 547:699-728. It is not the same thing as the continuous + relationship of + :func:`HJC_adjacent_mean_open_to_shut_time_pdf`, which cannot be compared + with data directly because the data have to be binned into ranges wide + enough to hold enough observations. + + The average over a range is + + E[open | u1 < shut < u2] = Int m(t) f(t) dt / Int f(t) dt, + + where ``f(t)`` is the apparent shut time density and ``m(t)`` the + conditional mean at ``t``. Both integrands come from + :func:`HJC_adjacent_mean_open_to_shut_time_pdf`, whose ratio is ``m(t)``, + so the numerator and denominator are integrated separately by Simpson's + rule over a log-spaced grid. + + Parameters + ---------- + ranges : sequence of (float, float) + Shut time ranges [s]. An infinite upper limit is replaced by a + concentration of the density: 100 times the slowest shut time + constant. + tres : float + Time resolution (dead time). + Q, QAA, QAF, QFF, QFA : array_like + The Q matrix and its submatrices. + points : int + Grid points per range for the quadrature. + side : {'next', 'previous'} + Whether the shut time follows the opening ('next', the usual case) or + precedes it. + + Returns + ------- + means : ndarray + One mean open time [s] per range; ``nan`` for a range holding no + probability. + """ + if side not in ('next', 'previous'): + raise ValueError("side must be 'next' or 'previous'") + + Froots = asymptotic_roots(tres, QFF, QAA, QFA, QAF, + QFF.shape[0], QAA.shape[0]) + slowest = -1.0 / Froots.max() + + means = np.zeros(len(ranges)) + for i, (u1, u2) in enumerate(ranges): + hi = u2 if np.isfinite(u2) else max(u1 * 10.0, slowest * 100.0) + lo = max(u1, tres) + if not hi > lo: + means[i] = np.nan + continue + t = np.logspace(np.log10(lo), np.log10(hi), points) + mp, mn = HJC_adjacent_mean_open_to_shut_time_pdf(t, tres, Q, QAA, + QAF, QFF, QFA) + m = mn if side == 'next' else mp + # the denominator of m(t) is the apparent shut time density itself, + # so weighting by it and dividing recovers the range average + density = HJC_shut_time_pdf(t, tres, Q, QAA, QAF, QFF, QFA) + num = _trapezoid(m * density, t) + den = _trapezoid(density, t) + means[i] = num / den if den > 0 else np.nan + return means + + +def HJC_shut_time_pdf(t, tres, Q, QAA, QAF, QFF, QFA): + """ + Apparent (missed events) shut time probability density at times ``t``. + + The same quantity that appears as the denominator of + :func:`HJC_adjacent_mean_open_to_shut_time_pdf`, exposed on its own + because the range averages need it as a weight. + """ + kA, kF = QAA.shape[0], QFF.shape[0] + uA = np.ones((kA))[:, np.newaxis] + expQFF = qml.expQt(QFF, tres) + expQAA = qml.expQt(QAA, tres) + GAF, GFA = qml.iGs(Q, kA, kF) + eGAF = qml.eGs(GAF, GFA, kA, kF, expQFF) + eGFA = qml.eGs(GFA, GAF, kF, kA, expQAA) + phiF = qml.phiHJC(eGFA, eGAF, kF) + eigs, A = qml.eigs(-Q) + Feigvals, FZ00, FZ10, FZ11 = qml.Zxx(Q, eigs, A, kA, QAA, QFA, QAF, + expQAA, False) + Froots = asymptotic_roots(tres, QFF, QAA, QFA, QAF, kF, kA) + FR = qml.AR(Froots, tres, QFF, QAA, QFA, QAF, kF, kA) + + out = np.zeros(len(np.atleast_1d(t))) + for i, ti in enumerate(np.atleast_1d(t)): + eGFAt = qml.eGAF(ti, tres, Feigvals, FZ00, FZ10, FZ11, Froots, + FR, QFA, expQAA) + out[i] = np.dot(np.dot(phiF, eGFAt), uA)[0] + return out + + def adjacent_open_to_shut_range_pdf_components(u1, u2, QAA, QAF, QFF, QFA, phiA): """ Calculate time constants and areas for an ideal (no missed events) diff --git a/tests/test_apparent_means.py b/tests/test_apparent_means.py new file mode 100644 index 0000000..55b47fc --- /dev/null +++ b/tests/test_apparent_means.py @@ -0,0 +1,112 @@ +"""Mean apparent open and shut times, and conditional means. + +These pin the correction made on 2026-09-02, when +``exact_mean_open_shut_time`` was found to be too large by exactly twice the +dead time and ``HJC_adjacent_mean_open_to_shut_time_pdf`` by once. The error was +in neither function's algebra but in the offset assumed for the ``dARSdS`` +term, so both are checked here against references that owe nothing to the +implementation: + +* for a two-state channel, eqns (121) and (122) of Colquhoun & Hawkes (1995b), + which give the apparent means in closed form; +* for larger mechanisms, the mean of the asymptotic HJC pdf, ``tres + sum(area + * tau)``, computed from the roots and areas -- a different route through the + library entirely; +* and the limit ``tres -> 0``, where the apparent means must become the ideal + ones. +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from scalcs import qmatlib as qml +from scalcs import scalcslib as scl +from scalcs.samples import samples + +#: The C-O sample: mean open 1/50 s, mean shut 1/20 s. +CO_MEAN_OPEN, CO_MEAN_SHUT = 1 / 50.0, 1 / 20.0 + + +def blue_book_means(mo, ms, tres): + """Apparent mean open and shut times of a two-state channel. + + Eqns (121) and (122) of Colquhoun & Hawkes (1995b), p. 456, with the same + dead time imposed on openings and shuttings. + """ + emo = tres + (mo + ms) * np.exp(tres / ms) - (tres + ms) + ems = tres + (mo + ms) * np.exp(tres / mo) - (tres + mo) + return emo, ems + + +def asymptotic_mean_open(mec, tres): + """tres + sum(area * tau) for the asymptotic HJC open time pdf.""" + roots = scl.asymptotic_roots(tres, mec.QAA, mec.QFF, mec.QAF, mec.QFA, + mec.kA, mec.kF) + GAF, GFA = qml.iGs(mec.Q, mec.kA, mec.kF) + areas = scl.asymptotic_areas(tres, roots, mec.QAA, mec.QFF, mec.QAF, + mec.QFA, mec.kA, mec.kF, GAF, GFA) + return tres + float(np.sum(areas * (-1.0 / roots))) + + +# ------------------------------------------- two states, against the algebra + +@pytest.mark.parametrize('tres', [1e-3, 2e-3, 5e-3, 1e-2, 2e-2]) +def test_two_state_apparent_means(tres): + mec = samples.CO() + expected = blue_book_means(CO_MEAN_OPEN, CO_MEAN_SHUT, tres) + npt.assert_allclose(scl.exact_mean_open_shut_time(mec, tres), expected, + rtol=1e-9) + + +def test_two_state_conditional_mean_is_flat(): + """A two-state channel has no correlations, so the mean open time cannot + depend on the adjacent shut time, and must equal the apparent mean.""" + mec = samples.CO() + tres = 5e-3 + expected, _ = blue_book_means(CO_MEAN_OPEN, CO_MEAN_SHUT, tres) + shut = np.array([5e-3, 1e-2, 3e-2, 1e-1, 5e-1]) + mp, mn = scl.HJC_adjacent_mean_open_to_shut_time_pdf( + shut, tres, mec.Q, mec.QAA, mec.QAF, mec.QFF, mec.QFA) + npt.assert_allclose(mp, expected, rtol=1e-9) + npt.assert_allclose(mn, expected, rtol=1e-9) + + +# -------------------------------------- more states, against a second route + +@pytest.mark.parametrize('name, conc, tres', [ + ('CH82', 100e-9, 30e-6), + ('AChR_diamond', 30e-9, 25e-6), + ('AChR_diamond', 10e-6, 25e-6), +]) +def test_mean_open_matches_the_asymptotic_pdf(name, conc, tres): + mec = getattr(samples, name)() + mec.set_eff('c', conc) + mean_open, _ = scl.exact_mean_open_shut_time(mec, tres) + npt.assert_allclose(mean_open, asymptotic_mean_open(mec, tres), rtol=1e-5) + + +def test_conditional_mean_agrees_at_long_shut_times(): + """Far out, the conditional mean must lie within the range the + unconditional apparent mean and the ideal mean bracket.""" + mec = samples.AChR_diamond() + mec.set_eff('c', 30e-9) + tres = 25e-6 + mp, mn = scl.HJC_adjacent_mean_open_to_shut_time_pdf( + np.array([1.0]), tres, mec.Q, mec.QAA, mec.QAF, mec.QFF, mec.QFA) + assert 0.0 < mn[0] < scl.exact_mean_open_shut_time(mec, tres)[0] + assert 0.0 < mp[0] < scl.exact_mean_open_shut_time(mec, tres)[0] + + +# ------------------------------------------------------------- the tres = 0 limit + +@pytest.mark.parametrize('name, conc', [('CO', None), + ('CH82', 100e-9), + ('AChR_diamond', 30e-9)]) +def test_zero_dead_time_gives_the_ideal_means(name, conc): + mec = getattr(samples, name)() + if conc is not None: + mec.set_eff('c', conc) + ideal_open = float(np.sum(qml.phiA(mec) @ (-np.linalg.inv(mec.QAA)))) + mean_open, _ = scl.exact_mean_open_shut_time(mec, 1e-12) + npt.assert_allclose(mean_open, ideal_open, rtol=1e-6) diff --git a/tests/test_real_dtype.py b/tests/test_real_dtype.py index 64e12a0..42111b1 100644 --- a/tests/test_real_dtype.py +++ b/tests/test_real_dtype.py @@ -93,20 +93,29 @@ def test_maxPopen_with_dead_time_is_real_and_silent(): if issubclass(w.category, np.exceptions.ComplexWarning)] assert not np.iscomplexobj(maxP), 'maxPopen returned %r' % (maxP,) assert not complex_warnings, '%d ComplexWarning(s)' % len(complex_warnings) - assert float(maxP) == pytest.approx(0.953453875851444, rel=1e-9) - - -# The two EC50 reference values below changed on 2026-09-02, when popen.EC50 -# was made to converge on the concentration instead of stopping as soon as -# Popen was within 0.001 of half maximal. The old values, 2.3560804576936233e-06 -# and 2.4032020668474956e-06, sit at responses of 0.4994 and 0.4999 rather than -# 0.5. Nothing about the real-dtype cast these tests guard has changed. + assert float(maxP) == pytest.approx(0.970082292957507, rel=1e-9) + + +# Every finite-dead-time reference value in this file changed on 2026-09-02, +# for two reasons, neither of them to do with the real-dtype cast these tests +# guard. +# +# First, popen.EC50 was made to converge on the concentration instead of +# stopping as soon as Popen was within 0.001 of half maximal; the old values, +# 2.3560804576936233e-06 and 2.4032020668474956e-06, sat at responses of 0.4994 +# and 0.4999 rather than 0.5. +# +# Second, and larger, scalcslib.exact_mean_open_shut_time was found to be too +# big by twice the dead time, which fed straight into Popen at finite +# resolution. maxPopen at a 30 us dead time was 0.9534, BELOW its value of +# 0.9677 with no dead time -- the wrong way round, since missing brief +# shuttings lengthens apparent openings. It is now 0.9701. def test_EC50_with_dead_time_is_real(): m = samples.CH82() ec50 = popenlib.EC50(m, TRES) assert not np.iscomplexobj(ec50) - assert float(ec50) == pytest.approx(2.3588952204118414e-06, rel=1e-9) + assert float(ec50) == pytest.approx(2.4042844021469254e-06, rel=1e-9) def test_dead_time_results_are_unchanged_by_the_real_cast(): From 2396306adf924a7af258230628ef40106f60d02c Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Fri, 4 Sep 2026 14:06:45 +0100 Subject: [PATCH 2/2] Re-derive dARSdS: one sign, and the callers can add tres again The 2*tres offset patched at the call sites in the previous commit had a single cause, now found by doing the derivation the docstring only sketched. HJC give AR*(s) = [sI - QAA - QAF M(s) QFA]^-1 with M(s) = (I - exp(-(sI - QFF) tres)) (sI - QFF)^-1. Factorising as CHS96 eqn (3.6) does, AR*(s) = inv(VA(s)) inv(U(s)) with U(s) = sI - QAA and VA(s) = I - inv(U) QAF M(s) QFA, and differentiating the product of two inverses gives exactly the DARS formula already in the docstring. So the outer structure was right and only dVAds was in question. Differentiating M with N(s) = sI - QFF gives M'(0) = -tres expQFF inv(QFF) - SFF inv(QFF)^2, and hence dVAds = -inv(QAA) GAF SFF GFA - tres GAF expQFF GFA - GAF SFF inv(QFF) GFA Two of those three terms were coded correctly. The dead-time term was coded POSITIVE. That made dVAds too large by 2 tres GAF expQFF GFA, and every apparent mean duration built on it too long by 2 tres. With the sign fixed, the callers add tres rather than subtracting it, which is what the theory says: the density is written in u = t - tres and integrates to 1, since phi inv(VA) GAF expQFF u = 1. The numbers do not change -- all four independent checks in test_apparent_means.py still hold, and the 590 tests pass -- but the correction now lives where the error was instead of being cancelled three times downstream. Anyone calling dARSdS directly gets a correct answer for the first time. Three tests added. The one that matters differentiates AR*(s) numerically straight from its definition, using nothing else in the library, and compares with dARSdS at three dead times; a flipped sign fails it. The others pin the normalisation phi inv(VA) GAF expQFF u = 1 and the resulting mean = DARS term + tres. Co-Authored-By: Claude Opus 5 --- scalcs/qmatlib.py | 55 +++++++++++++++++++++----- scalcs/scalcslib.py | 32 +++++++-------- tests/test_apparent_means.py | 76 ++++++++++++++++++++++++++++++++++++ 3 files changed, 137 insertions(+), 26 deletions(-) diff --git a/scalcs/qmatlib.py b/scalcs/qmatlib.py index 7a000f7..56df848 100644 --- a/scalcs/qmatlib.py +++ b/scalcs/qmatlib.py @@ -524,14 +524,45 @@ def dARSdS(tres, QAA, QFF, GAF, GFA, expQFF, kA, kF): For same evaluation for shut states exhange A by F and F by A in function call. - SFF = I - exp(QFF * tres) - First evaluate [dVA(s) / ds] * s = 0. - dVAds = -inv(QAA) * GAF * SFF * GFA - GAF * SFF * inv(QFF) * GFA + - + tres * GAF * expQFF * GFA + The derivation, since a sign in it was wrong for a long time. - Then: DARS = inv(VA) * QAA^(-2) - inv(VA) * dVAds * inv(VA) * inv(QAA) = - = inv(VA) * [inv(QAA) - dVAds * inv(VA)] * inv(QAA) - where VA = I - GAF * SFF * GFA + HJC give the transform as + + AR*(s) = [sI - QAA - QAF M(s) QFA]^-1, + M(s) = (I - exp(-(sI - QFF) tres)) (sI - QFF)^-1. + + Factorise as CHS96 eqn (3.6) does, pulling out (sI - QAA): + + AR*(s) = VA(s)^-1 U(s)^-1, U(s) = sI - QAA, + VA(s) = I - U(s)^-1 QAF M(s) QFA. + + At s = 0, U = -QAA, U^-1 QAF = -GAF and M(0) QFA = SFF GFA, so + VA(0) = I - GAF SFF GFA, the VA below. Differentiating the product of + two inverses, with U' = I, + + DARS = [-d AR*/ds]_0 = inv(VA) QAA^-2 - inv(VA) dVAds inv(VA) inv(QAA) + = inv(VA) [inv(QAA) - dVAds inv(VA)] inv(QAA). + + For dVAds, with N(s) = sI - QFF and M = (I - exp(-N tres)) N^-1, + + M'(s) = tres exp(-N tres) N^-1 - (I - exp(-N tres)) N^-2, + M'(0) = -tres expQFF inv(QFF) - SFF inv(QFF)^2, + + and, using QAF = -QAA GAF and inv(QFF) QFA = -GFA, + + dVAds = -inv(QAA) GAF SFF GFA (from U^-2 QAF M(0) QFA) + - tres GAF expQFF GFA (from -U^-1 QAF M'(0) QFA) + - GAF SFF inv(QFF) GFA. + + **The middle term is negative.** It was coded positive, which made dVAds + too large by 2 tres GAF expQFF GFA and every apparent mean duration built + on it too long by 2 tres. Checked against eqns (121) and (122) of + Colquhoun & Hawkes (1995b) and against the mean of the asymptotic HJC pdf; + see tests/test_apparent_means.py. + + A caller wanting a mean apparent duration **adds** tres to + phi AR-term, because the density is written in u = t - tres and + phi inv(VA) GAF expQFF u = 1. Parameters ---------- @@ -565,8 +596,14 @@ def dARSdS(tres, QAA, QFF, GAF, GFA, expQFF, kA, kF): I = np.eye(kF) SFF = I - expQFF - #Q1 = tres * GAF * exp(QFF*tres) * GFA - Q1 = tres * np.dot(GAF, np.dot(expQFF, GFA)) + # Q1 = -tres * GAF * exp(QFF*tres) * GFA + # + # The sign of this term was wrong until the derivation in the docstring + # was carried out. M'(0) = -tres*expQFF*QFF^-1 - SFF*QFF^-2, so the + # contribution -GAF M'(0) QFA is -tres*GAF*expQFF*GFA and is negative. + # With it positive, dVAds was too large by 2*tres*GAF*expQFF*GFA, and + # every caller inherited an apparent mean duration 2*tres too long. + Q1 = -tres * np.dot(GAF, np.dot(expQFF, GFA)) #Q2 = GAF * SFF * inv(QFF) * GFA Q2 = np.dot(GAF, np.dot(SFF, np.dot(invQFF, GFA))) #Q3 = -inv(QAA) * GAF * SFF * GFA diff --git a/scalcs/scalcslib.py b/scalcs/scalcslib.py index aded413..c0b6412 100644 --- a/scalcs/scalcslib.py +++ b/scalcs/scalcslib.py @@ -701,17 +701,15 @@ def exact_mean_open_shut_time(mec, tres): DARS = qml.dARSdS(tres, mec.QAA, mec.QFF, GAF, GFA, expQFF, mec.kA, mec.kF) DFRS = qml.dARSdS(tres, mec.QFF, mec.QAA, GFA, GAF, expQAA, mec.kF, mec.kA) uF, uA = np.ones((mec.kF, 1)), np.ones((mec.kA, 1)) - # The DARS term below evaluates to the mean apparent duration PLUS tres, not - # minus it, so the mean is recovered by subtracting tres rather than adding it. - # Corrected 2026-09-02, having been wrong by 2*tres. Checked two ways: for a - # two-state channel it now reproduces eqns (121) and (122) of Colquhoun & - # Hawkes (1995b) exactly, and for CH82 and the AChR mechanism it agrees with - # the mean of the asymptotic HJC pdf, tres + sum(area * tau), to six figures. - # The root cause is the convention of qmatlib.dARSdS, which has not been - # re-derived against CHS96 eqn (3.6); the correction is applied at each call - # site instead, where what is wanted is unambiguous. - meanA = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] - tres - meanF = np.dot(phiF, np.dot(np.dot(DFRS, QexpQA), uA))[0] - tres + # The density is written in u = t - tres and integrates to 1, so a mean + # apparent duration is the DARS term PLUS tres. See the derivation in + # qmatlib.dARSdS: its dead-time term had the wrong sign until 2026-09-03, + # and while it did, this line had to subtract tres to cancel an error of + # 2*tres. Checked against eqns (121) and (122) of Colquhoun & Hawkes + # (1995b) and against the mean of the asymptotic HJC pdf; see + # tests/test_apparent_means.py. + meanA = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] + tres + meanF = np.dot(phiF, np.dot(np.dot(DFRS, QexpQA), uA))[0] + tres return meanA, meanF @@ -753,8 +751,8 @@ def exact_mean_time(tres, QAA, QFF, QAF, kA, kF, GAF, GFA): DARS = qml.dARSdS(tres, QAA, QFF, GAF, GFA, expQFF, kA, kF) uF = np.ones((kF, 1)) - # see exact_mean_open_shut_time for why this subtracts tres - mean = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] - tres + # see exact_mean_open_shut_time for why this adds tres + mean = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] + tres return mean @@ -1226,10 +1224,10 @@ def HJC_adjacent_mean_open_to_shut_time_pdf(sht, tres, Q, QAA, QAF, QFF, QFA): nom1 = np.dot(np.dot(phiF, eGFAt), col1)[0] nom2 = np.dot(np.dot(row1, eGFAt), uA)[0] # the same tres offset as in exact_mean_open_shut_time: these are - # conditional means of the apparent open time, and the DARS term - # carries an extra tres - mp.append(nom1 / denom - tres) - mn.append(nom2 / denom - tres) + # conditional means of the apparent open time, so the DARS term is + # short by tres + mp.append(nom1 / denom + tres) + mn.append(nom2 / denom + tres) return np.array(mp), np.array(mn) diff --git a/tests/test_apparent_means.py b/tests/test_apparent_means.py index 55b47fc..357f0fd 100644 --- a/tests/test_apparent_means.py +++ b/tests/test_apparent_means.py @@ -110,3 +110,79 @@ def test_zero_dead_time_gives_the_ideal_means(name, conc): ideal_open = float(np.sum(qml.phiA(mec) @ (-np.linalg.inv(mec.QAA)))) mean_open, _ = scl.exact_mean_open_shut_time(mec, 1e-12) npt.assert_allclose(mean_open, ideal_open, rtol=1e-6) + + +# --------------------------------------------------------------------------- +# dARSdS itself, added 2026-09-03 when the root cause was finally re-derived. +# +# The offset was traced to one sign in `qmatlib.dARSdS`: the dead-time term of +# dVAds is -tres*GAF*expQFF*GFA and had been coded positive, which made the +# function too large by 2*tres and left every caller to cancel it. Fixing it +# means callers now ADD tres, as the theory says, and it also means anyone +# calling dARSdS directly gets a correct answer -- which the tests above, +# working only through the callers, could never have shown. +# +# The check here owes nothing to the rest of the library: it differentiates +# the Laplace transform numerically, straight from its definition. + +def _AR_star(s, tres, QAA, QFF, QAF, QFA): + """AR*(s) from its definition, for numerical differentiation. + + AR*(s) = [sI - QAA - QAF M(s) QFA]^-1, + M(s) = (I - exp(-(sI - QFF) tres)) (sI - QFF)^-1 + """ + from scipy.linalg import expm + kA, kF = QAA.shape[0], QFF.shape[0] + N = s * np.eye(kF) - QFF + M = (np.eye(kF) - expm(-N * tres)) @ np.linalg.inv(N) + return np.linalg.inv(s * np.eye(kA) - QAA - QAF @ M @ QFA) + + +@pytest.mark.parametrize("tres", [10e-6, 25e-6, 100e-6]) +def test_dARSdS_is_minus_the_derivative_of_the_transform(tres): + """dARSdS must equal [-d AR*(s)/ds] at s = 0, by finite differences.""" + mec = samples.CH82() + mec.set_eff('c', 100e-9) + expQFF = qml.expQt(mec.QFF, tres) + GAF, GFA = qml.iGs(mec.Q, mec.kA, mec.kF) + + got = qml.dARSdS(tres, mec.QAA, mec.QFF, GAF, GFA, expQFF, + mec.kA, mec.kF) + + # central difference on s, in units set by the slowest rate + h = 1e-4 * abs(np.diag(mec.QAA)).min() + plus = _AR_star(h, tres, mec.QAA, mec.QFF, mec.QAF, mec.QFA) + minus = _AR_star(-h, tres, mec.QAA, mec.QFF, mec.QAF, mec.QFA) + expected = -(plus - minus) / (2 * h) + + npt.assert_allclose(got, expected, rtol=1e-5, atol=0) + + +def test_mean_is_the_dARSdS_term_plus_tres(): + """The documented relationship between dARSdS and the apparent mean. + + Callers add tres because the density is written in u = t - tres and + integrates to 1. Check both halves: the norm, and the mean. + """ + tres = 25e-6 + mec = samples.AChR_diamond() + mec.set_eff('c', 100e-9) + + expQFF = qml.expQt(mec.QFF, tres) + expQAA = qml.expQt(mec.QAA, tres) + GAF, GFA = qml.iGs(mec.Q, mec.kA, mec.kF) + eGAF = qml.eGs(GAF, GFA, mec.kA, mec.kF, expQFF) + eGFA = qml.eGs(GFA, GAF, mec.kF, mec.kA, expQAA) + phiA = qml.phiHJC(eGAF, eGFA, mec.kA) + uF = np.ones((mec.kF, 1)) + + SFF = np.eye(mec.kF) - expQFF + VA = np.eye(mec.kA) - GAF @ SFF @ GFA + norm = float((phiA @ np.linalg.inv(VA) @ GAF @ expQFF @ uF)[0]) + npt.assert_allclose(norm, 1.0, rtol=1e-12) + + DARS = qml.dARSdS(tres, mec.QAA, mec.QFF, GAF, GFA, expQFF, + mec.kA, mec.kF) + term = float((phiA @ (DARS @ (mec.QAF @ expQFF)) @ uF)[0]) + meanA, _ = scl.exact_mean_open_shut_time(mec, tres) + npt.assert_allclose(term + tres, meanA, rtol=1e-12)