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 a5b6dca..c0b6412 100644 --- a/scalcs/scalcslib.py +++ b/scalcs/scalcslib.py @@ -701,9 +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)) - # 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 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 @@ -745,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)) - # 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 adds tres + mean = np.dot(phiA, np.dot(np.dot(DARS, QexpQF), uF))[0] + tres return mean @@ -1217,11 +1223,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, 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) +#: 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..357f0fd --- /dev/null +++ b/tests/test_apparent_means.py @@ -0,0 +1,188 @@ +"""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) + + +# --------------------------------------------------------------------------- +# 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) 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():