From 1c2fbbb2b2401b62ecc9c4077d4379623dfafa38 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 8 Sep 2026 10:34:49 +0100 Subject: [PATCH] Check the simulator against Q-matrix theory for multi-state mechanisms test_scsim.py checks the two-state CO mechanism, where the means are 1/alpha and 1/beta. That is the easy case: with one open and one shut state the dwell times are single exponentials, and a simulator can get them right while getting the multi-state case wrong. These check CH82 and AChR_diamond, where the answer is a Q-matrix expression rather than a reciprocal rate: the ideal means phi_A (-Q_AA)^-1 u and phi_F (-Q_FF)^-1 u, P(open) against the equilibrium occupancies, and -- after a dead time is imposed -- the apparent means against scalcslib.exact_mean_open_shut_time. The last earns its keep. Nothing else ties the simulator to the missed-events code: test_apparent_means.py checks that code against closed forms and against itself by another route, but never against data. An error of two dead times in the apparent mean open time -- exactly what was there until the dARSdS sign was corrected -- is caught here even if every analytic cross-check had been adjusted to agree with it. Verified by mutation: putting the error back makes both mechanisms fail. How much power the check has was measured rather than asserted, because it depends on the dead time as a fraction of the mean open time. At 400000 intervals and tres = 25 us: mechanism one tres two tres actual disagreement CH82 3.4 SE 6.7 SE 1.3 SE AChR_diamond 5.6 SE 11.3 SE 1.0 SE So a two-dead-time error is caught on both, and a one-dead-time error on the AChR mechanism but not on CH82, whose open times are a hundred times the dead time. The test that guards this asserts the two-dead-time case everywhere and the one-dead-time case only where the measurement supports it; claiming otherwise would make a passing suite mean less than it appears to. An earlier draft did claim it, and failed honestly at 100000 intervals, which is how the numbers above came to be measured. Tolerances are multiples of the measured standard error rather than fixed percentages, so the tests neither pass by being loose on a sharp quantity nor fail by being tight on a noisy one. Seeds are fixed. Marked slow: 616 passed with them, 598 in 32 s without. Co-Authored-By: Claude Opus 5 --- tests/test_scsim_vs_theory.py | 192 ++++++++++++++++++++++++++++++++++ 1 file changed, 192 insertions(+) create mode 100644 tests/test_scsim_vs_theory.py diff --git a/tests/test_scsim_vs_theory.py b/tests/test_scsim_vs_theory.py new file mode 100644 index 0000000..9e63b9b --- /dev/null +++ b/tests/test_scsim_vs_theory.py @@ -0,0 +1,192 @@ +"""The simulator against Q-matrix theory, for mechanisms with several states. + +``test_scsim.py`` already checks the two-state CO mechanism, where the means +are ``1/alpha`` and ``1/beta``. That is the easy case: with one open and one +shut state the dwell times are single exponentials, and a simulator can get +them right while getting the multi-state case wrong. + +These check the cases where the answer is a Q-matrix expression rather than a +reciprocal rate: + +* the ideal means, ``phi_A (-Q_AA)^-1 u`` and ``phi_F (-Q_FF)^-1 u``, for + mechanisms with three or four states of each kind; +* ``P(open)`` against the equilibrium occupancies; +* and, after a dead time is imposed, the **apparent** means against + :func:`scalcslib.exact_mean_open_shut_time` -- the HJC missed-events theory. + +The last is the one that earns its keep. Nothing else ties the simulator to the +missed-events code: ``test_apparent_means.py`` checks that code against closed +forms and against itself by another route, but never against data. An error of +two dead times in the apparent mean open time -- exactly the error that was +there until the ``dARSdS`` sign was corrected -- would be caught here even if +every analytic cross-check had been adjusted to agree with it. + +How much power that check has was measured rather than assumed, because it +depends on the dead time as a fraction of the mean open time and so differs by +mechanism. At 400 000 intervals and ``tres = 25 us``: + +=============== =========== =========== ================== +mechanism one tres two tres actual disagreement +=============== =========== =========== ================== +CH82 3.4 SE 6.7 SE 1.3 SE +AChR_diamond 5.6 SE 11.3 SE 1.0 SE +=============== =========== =========== ================== + +So at the tolerance used here a **two**-dead-time error is caught on both, and +a one-dead-time error on the AChR mechanism but not on CH82, whose open times +are a hundred times the dead time. +:func:`test_the_apparent_open_time_check_has_power` asserts the two-dead-time +case, so a passing suite means the check could actually fail. + +Marked slow: those standard errors need of the order of 10^5 intervals. + +The tolerance is a multiple of the *measured* standard error rather than a +fixed percentage, so the tests neither pass by being loose on a sharp quantity +nor fail by being tight on a noisy one. Seeds are fixed, so a failure is a real +disagreement rather than a bad draw. +""" + +import numpy as np +import pytest + +from scalcs import qmatlib as qml +from scalcs import scalcslib as scl +from scalcs import scsim +from scalcs.samples import samples + +#: How many standard errors of disagreement to allow. With fixed seeds these +#: are deterministic; 4 leaves room for the estimator's own skew without +#: admitting an error of the size worth catching. +NSIGMA = 4.0 + +TRES = 25e-6 +NINT = 400000 + + +def mechanisms(): + ch82 = samples.CH82() + ch82.set_eff('c', 100e-9) + achr = samples.AChR_diamond() + achr.set_eff('c', 100e-9) + return [("CH82", ch82), ("AChR_diamond", achr)] + + +def ideal_means(mec): + """phi_A (-Q_AA)^-1 u and phi_F (-Q_FF)^-1 u, and P(open).""" + kA, kF = mec.kA, mec.k - mec.kA + QAA, QFF = mec.Q[:kA, :kA], mec.Q[kA:, kA:] + phiA = np.asarray(qml.phiA(mec)).reshape(1, kA) + phiF = np.asarray(qml.phiF(mec)).reshape(1, kF) + mopen = float(np.squeeze(phiA @ np.linalg.inv(-QAA) @ np.ones((kA, 1)))) + mshut = float(np.squeeze(phiF @ np.linalg.inv(-QFF) @ np.ones((kF, 1)))) + popen = float(np.asarray(qml.pinf(mec.Q)).ravel()[:kA].sum()) + return mopen, mshut, popen + + +def measured(tints, ampls): + """Mean open, mean shut and P(open) of a record, with standard errors.""" + op, sh = tints[ampls > 0], tints[ampls == 0] + return dict( + mopen=op.mean(), se_open=op.std(ddof=1) / np.sqrt(op.size), + mshut=sh.mean(), se_shut=sh.std(ddof=1) / np.sqrt(sh.size), + popen=op.sum() / tints.sum(), n_open=op.size) + + +@pytest.fixture(scope="module", params=mechanisms(), ids=lambda m: m[0]) +def record(request): + """One long simulated record per mechanism, reused by every test here.""" + name, mec = request.param + tints, ampls, _ = scsim.simulate_intervals(mec, nintmax=NINT, seed=20030547) + at, aa = scsim.impose_resolution(tints, ampls, TRES) + return name, mec, tints, ampls, at, aa + + +@pytest.mark.slow +def test_ideal_mean_open_time(record): + """phi_A (-Q_AA)^-1 u, not 1/alpha: several open states.""" + _, mec, ti, ai, _, _ = record + want, _, _ = ideal_means(mec) + got = measured(ti, ai) + assert abs(got["mopen"] - want) < NSIGMA * got["se_open"], ( + f"simulated {got['mopen']:.6g} s, theory {want:.6g} s, " + f"{abs(got['mopen'] - want) / got['se_open']:.1f} SE") + + +@pytest.mark.slow +def test_ideal_mean_shut_time(record): + _, mec, ti, ai, _, _ = record + _, want, _ = ideal_means(mec) + got = measured(ti, ai) + assert abs(got["mshut"] - want) < NSIGMA * got["se_shut"], ( + f"simulated {got['mshut']:.6g} s, theory {want:.6g} s") + + +@pytest.mark.slow +def test_popen(record): + """P(open) against the equilibrium occupancies of the open states.""" + _, mec, ti, ai, _, _ = record + _, _, want = ideal_means(mec) + got = measured(ti, ai)["popen"] + assert abs(got - want) / want < 0.02, f"simulated {got:.6g}, theory {want:.6g}" + + +@pytest.mark.slow +def test_apparent_mean_open_time_matches_the_missed_events_theory(record): + """The simulator against exact_mean_open_shut_time, after a dead time. + + This is the check with power over the missed-events code -- the only one + anywhere that tests it against simulated data rather than against another + expression. See the module docstring for how much power, which depends on + the mechanism. + """ + _, mec, _, _, at, aa = record + want, _ = scl.exact_mean_open_shut_time(mec, TRES) + got = measured(at, aa) + off = abs(got["mopen"] - want) + assert off < NSIGMA * got["se_open"], ( + f"simulated {got['mopen']*1e3:.6g} ms, theory {want*1e3:.6g} ms, " + f"{off / got['se_open']:.1f} SE; one dead time is " + f"{TRES / got['se_open']:.1f} SE") + + +@pytest.mark.slow +def test_apparent_mean_shut_time_matches_the_missed_events_theory(record): + """The same for shut times. + + Much less sensitive -- the apparent mean shut time is large and its + distribution broad, so a dead time is a small fraction of a standard error + -- but a gross error in the shut-time branch would still show. + """ + _, mec, _, _, at, aa = record + _, want = scl.exact_mean_open_shut_time(mec, TRES) + got = measured(at, aa) + assert abs(got["mshut"] - want) < NSIGMA * got["se_shut"], ( + f"simulated {got['mshut']*1e3:.6g} ms, theory {want*1e3:.6g} ms") + + +@pytest.mark.slow +def test_the_apparent_open_time_check_has_power(record): + """The check above is only worth having if it would fail when it should. + + A theory value wrong by **two** dead times -- the error that was actually + present -- must lie outside the tolerance on every mechanism here. One dead + time is asserted only where the measurement supports it: on CH82 it is 3.4 + SE, inside the 4 SE tolerance, because its open times are a hundred times + the dead time. Claiming otherwise would make a passing suite mean less than + it appears to. + """ + name, mec, _, _, at, aa = record + want, _ = scl.exact_mean_open_shut_time(mec, TRES) + got = measured(at, aa) + tol = NSIGMA * got["se_open"] + + for wrong in (want + 2 * TRES, want - 2 * TRES): + assert abs(got["mopen"] - wrong) > tol, ( + f"{name}: an answer wrong by two dead times would pass. The record " + f"is too short or the tolerance too loose to be worth anything.") + + sensitivity = TRES / got["se_open"] + if name == "AChR_diamond": + assert sensitivity > NSIGMA, ( + f"one dead time is only {sensitivity:.1f} SE here; it used to be " + f"{5.6:.1f}, so the record or the mechanism has changed")