From 1316e050a3d891d4a899edc636b092369e7a13f9 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 12:39:35 +0000 Subject: [PATCH 1/3] Let a mechanism be indexed by state name States are stored in the order Mechanism sorts them into, by state type: open first, then short-lived shut, then long-lived shut. That is rarely the order a paper tabulates them in. CH82 is stored as AR*, A2R*, AR, A2R, R, and the ten-state flip mechanism of Burzomato et al. (2004) differs from the published ordering by a swap of two states. Reading a Q matrix by position against a published table therefore invites an error that nothing catches, because the wrong entry is still a plausible number. Adds, all computed on demand and none of them caching anything: state_names names, in Q-matrix order state_index(name) row and column of a state rates_by_edge() {(from, to): Rate} for every transition rate(from, to) one transition, KeyError if there is no such step permutation(order) index array to reorder Q into a given order rate() raising rather than returning None is the point of it: asking for a step the mechanism does not have is a question about a different mechanism. When a name is simply unknown it says so, and lists the states there are. Reproducing a published paper previously meant writing a wrapper class to do this; that is what prompted it. Co-Authored-By: Claude Opus 5 --- scalcs/mechanism.py | 60 ++++++++++++++++++++++++++++++++++++++++ tests/test_mechanism.py | 61 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+) diff --git a/scalcs/mechanism.py b/scalcs/mechanism.py index 03adde4..0b16e53 100644 --- a/scalcs/mechanism.py +++ b/scalcs/mechanism.py @@ -780,6 +780,66 @@ def update_submat(self): self.QAI = self.Q[:self.kA, self.kA:self.k] self.QGG = self.Q[:self.kG, :self.kG] + # ----------------------------------------------------------------- + # Looking states and rates up by name + # + # States are stored in the order this class sorts them into, by state + # type: open states first, then short-lived shut, then long-lived shut. + # That is rarely the order a paper tabulates them in - CH82 is stored + # as AR*, A2R*, AR, A2R, R - so reading a Q matrix by position against + # a published table invites silent error. These let it be read by name. + # ----------------------------------------------------------------- + + @property + def state_names(self): + """Names of the states, in Q-matrix order.""" + return [state.name for state in self.States] + + def state_index(self, name): + """Row and column of `name` in the Q matrix.""" + for i, state in enumerate(self.States): + if state.name == name: + return i + raise KeyError('no state %r in this mechanism; it has %s' + % (name, ', '.join(self.state_names))) + + def rates_by_edge(self): + """{(from_name, to_name): Rate} for every transition. + + Built fresh on each call: a Rate's value can be changed in place, and + a stale cache of these would be worse than no cache at all. + """ + return {(rate.State1.name, rate.State2.name): rate + for rate in self.Rates} + + def rate(self, from_name, to_name): + """The Rate for one transition. + + Raises KeyError if the mechanism has no such transition, rather than + returning something for a step that does not exist. + """ + for rate in self.Rates: + if (rate.State1.name == from_name + and rate.State2.name == to_name): + return rate + # name an unknown state specifically; it is the likelier mistake + for name in (from_name, to_name): + self.state_index(name) + raise KeyError('%r -> %r is not a transition of this mechanism' + % (from_name, to_name)) + + def permutation(self, order): + """Index array that reorders the Q matrix into `order`. + + Use it as ``Q[np.ix_(perm, perm)]`` to put a Q matrix, or any matrix + indexed the same way, into the order a paper uses. + """ + if len(order) != self.k: + raise ValueError( + 'order lists %d states, but this mechanism has %d: %s' + % (len(order), self.k, ', '.join(self.state_names))) + return [self.state_index(name) for name in order] + def set_eff(self, eff, val): self.set_effdict({eff:val}) diff --git a/tests/test_mechanism.py b/tests/test_mechanism.py index ed17782..5a8bb08 100644 --- a/tests/test_mechanism.py +++ b/tests/test_mechanism.py @@ -410,3 +410,64 @@ def test_impose_limits_clips_out_of_range(self, cco): assert cco.Rates[0].unit_rate() <= cco.Rates[0].limits[0][1] # restore cco.Rates[0].rateconstants = original + + +# --------------------------------------------------------------------------- +# Looking states and rates up by name +# --------------------------------------------------------------------------- + +class TestStateLookup: + """Indexing a mechanism by state NAME rather than by position. + + The order SCALCS stores states in is its own - it sorts them by state type + - and is not the order a paper tabulates them in. CH82 is stored as + ['AR*', 'A2R*', 'AR', 'A2R', 'R']. Reading a Q matrix by position against a + published table is therefore a standing invitation to silent error; these + helpers exist so that it can be read by name instead. + """ + + def test_state_names_are_in_q_matrix_order(self, ch82): + assert ch82.state_names == ['AR*', 'A2R*', 'AR', 'A2R', 'R'] + + def test_state_index_agrees_with_the_q_matrix(self, ch82): + for i, name in enumerate(ch82.state_names): + assert ch82.state_index(name) == i + + def test_state_index_names_the_state_it_could_not_find(self, ch82): + with pytest.raises(KeyError, match='A3R'): + ch82.state_index('A3R') + + def test_rates_by_edge_has_one_entry_per_rate(self, ch82): + edges = ch82.rates_by_edge() + assert len(edges) == len(ch82.Rates) + assert ('R', 'AR') in edges and ('AR', 'R') in edges + + def test_rate_returns_the_transition_asked_for(self, ch82): + assert ch82.rate('R', 'AR').name.strip() == '2k(+1)' + assert ch82.rate('AR', 'R').name.strip() == 'k(-1)' + + def test_rate_refuses_a_transition_the_mechanism_does_not_have(self, ch82): + """R -> A2R is not a step of CH82: binding is one molecule at a time.""" + with pytest.raises(KeyError): + ch82.rate('R', 'A2R') + + def test_permutation_reindexes_q_into_the_given_order(self, ch82): + order = ['R', 'AR', 'A2R', 'AR*', 'A2R*'] + perm = ch82.permutation(order) + + Q = np.asarray(ch82.Q, dtype=float) + reordered = Q[np.ix_(perm, perm)] + + # element [0, 1] of the reordered matrix must now be R -> AR + assert reordered[0, 1] == pytest.approx( + Q[ch82.state_index('R'), ch82.state_index('AR')]) + # and the diagonal must still be the same set of exit rates + assert sorted(np.diag(reordered)) == pytest.approx(sorted(np.diag(Q))) + + def test_permutation_rejects_an_incomplete_order(self, ch82): + with pytest.raises(ValueError): + ch82.permutation(['R', 'AR']) + + def test_permutation_rejects_an_unknown_state(self, ch82): + with pytest.raises(KeyError, match='A3R'): + ch82.permutation(['R', 'AR', 'A2R', 'AR*', 'A3R']) From 4df79839f542a5d4b782f556a3bc2aa2be3846e8 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 12:39:35 +0000 Subject: [PATCH 2/3] Consolidate the transition probability and frequency matrices transition_probability had three implementations, two in scalcslib and one in scsim; transition_frequency had two. They had drifted: - scalcslib's copies carried empty docstrings; - scalcslib's transition_probability did Q.copy() without casting, so an integer Q matrix divided integer-wise and returned silently wrong numbers rather than raising; - scalcslib's transition_frequency computed pinf internally and so could only ever give equilibrium frequencies. Both now live in qmatlib, with the rest of the plain Q-matrix algebra, and are re-exported from scalcslib and scsim so nothing importing them breaks. transition_frequency now takes the occupancies as an optional argument. Called with one argument it still uses pinf, which is what every existing call site does; passing p gives the frequencies at an instant during a concentration jump, when the system is not at equilibrium. The tests pin the definitions, that neither function modifies its argument, that an integer Q is no longer truncated, and that at equilibrium detailed balance holds edge by edge. That last one is only approximate for CH82, whose MR-constrained rate is stored rounded to 0.66667, so its cycle closes to about 1e-9; a mechanism whose rates come from a fit closes to 1e-12. Co-Authored-By: Claude Opus 5 --- scalcs/qmatlib.py | 55 ++++++++++++++ scalcs/scalcslib.py | 29 +++----- scalcs/scsim.py | 18 ++--- tests/test_transition_matrices.py | 117 ++++++++++++++++++++++++++++++ 4 files changed, 186 insertions(+), 33 deletions(-) create mode 100644 tests/test_transition_matrices.py diff --git a/scalcs/qmatlib.py b/scalcs/qmatlib.py index 1e8abc4..8f7a326 100644 --- a/scalcs/qmatlib.py +++ b/scalcs/qmatlib.py @@ -215,6 +215,61 @@ def pinf(Q): pinf = np.append(pinf, 1 - np.sum(pinf)) return pinf +def transition_probability(Q): + r""" + Embedded-chain (jump) probability matrix of a Q matrix. + + .. math:: \pi_{ij} = \frac{q_{ij}}{-q_{ii}}, \quad \pi_{ii} = 0 + + Given that the channel leaves state *i*, :math:`\pi_{ij}` is the + probability that the jump is to *j*, so each row sums to 1. + + Parameters + ---------- + Q : array_like, shape (k, k) + + Returns + ------- + pi : ndarray, shape (k, k) + """ + Q = np.asarray(Q, dtype=float) + pi = Q / -np.diag(Q)[:, None] + np.fill_diagonal(pi, 0.0) + return pi + +def transition_frequency(Q, p=None): + r""" + Frequency of each transition: how many occur per second. + + .. math:: f_{ij} = p_i q_{ij}, \quad f_{ii} = 0 + + Parameters + ---------- + Q : array_like, shape (k, k) + p : array_like, shape (k,), optional + State occupancies. Defaults to the equilibrium occupancies, + :func:`pinf`. Pass the occupancies explicitly to get the frequencies + at some instant during a concentration jump, when the system is not + at equilibrium. + + Returns + ------- + f : ndarray, shape (k, k) + + Notes + ----- + At equilibrium a mechanism obeying microscopic reversibility satisfies + detailed balance, :math:`f_{ij} = f_{ji}` for every pair of states, so + the returned matrix is symmetric to within the precision in which the + constrained rate constants are stored. + """ + Q = np.asarray(Q, dtype=float) + if p is None: + p = pinf(Q) + f = np.asarray(p, dtype=float)[:, None] * Q + np.fill_diagonal(f, 0.0) + return f + def iGs(Q, kA, kB): r""" Calculate GBA and GAB matrices (Eq. 1.25, CH82). diff --git a/scalcs/scalcslib.py b/scalcs/scalcslib.py index c0b6412..f121654 100644 --- a/scalcs/scalcslib.py +++ b/scalcs/scalcslib.py @@ -1590,26 +1590,15 @@ def printout_distributions(mec, tres, eff='c'): return str -def transition_probability(Q): - """ - """ - k = Q.shape[0] - pi = Q.copy() - for i in range(k): - pi[i] = pi[i] / -Q[i,i] - pi[i,i] = 0 - return pi - -def transition_frequency(Q): - """ - """ - k = Q.shape[0] - pinf = qml.pinf(Q) - f = Q.copy().transpose() - for i in range(k): - f[i] = f[i] * pinf - f[i,i] = 0 - return f.transpose() +# Both of these now live in qmatlib, which is where the rest of the plain +# Q-matrix algebra is. They are re-exported here, undocumented duplicates +# having drifted apart: this one's transition_probability did not cast to +# float, so an integer Q divided integer-wise, and transition_frequency could +# only ever give equilibrium frequencies. qmatlib.transition_frequency takes +# the occupancies as an optional argument; called with one argument, as it +# always was here, it still uses pinf. +transition_probability = qml.transition_probability +transition_frequency = qml.transition_frequency def correlation_coefficient(cov, var1, var2): """ diff --git a/scalcs/scsim.py b/scalcs/scsim.py index b3c2927..171b281 100644 --- a/scalcs/scsim.py +++ b/scalcs/scsim.py @@ -36,24 +36,16 @@ from dcio.analysis import bursts as _bursts from dcio.formats import scn as _scn +from scalcs import qmatlib as _qml + #: SCN property flag marking an interval of unknown length. Re-exported from #: dcio, which is where it is defined. FLAG_UNUSABLE = _scn.FLAG_UNUSABLE -def transition_probability(Q): - r"""Embedded-chain (jump) probability matrix of a Q matrix. - - ``pi_ij = q_ij / -q_ii`` for ``i != j`` and ``pi_ii = 0``: given that the - channel leaves state ``i``, ``pi_ij`` is the probability the jump is to - ``j``. Each off-diagonal row therefore sums to 1. - """ - k = Q.shape[0] - pi = Q.copy().astype(float) - for i in range(k): - pi[i] = pi[i] / -Q[i, i] - pi[i, i] = 0.0 - return pi +# Defined in qmatlib with the rest of the Q-matrix algebra; re-exported here, +# where the simulation reaches for it, and where it was a second copy. +transition_probability = _qml.transition_probability def next_state(present, picum, tmean, kA, opamp): diff --git a/tests/test_transition_matrices.py b/tests/test_transition_matrices.py new file mode 100644 index 0000000..55d6319 --- /dev/null +++ b/tests/test_transition_matrices.py @@ -0,0 +1,117 @@ +"""Tests for the embedded-chain and transition-frequency matrices. + +Both quantities had grown three separate implementations - two in +``scalcslib`` and one in ``scsim`` - which is why they are consolidated in +``qmatlib`` and re-exported from the old homes. + +Reference mechanism: CH82 at 100 nM. +""" +import numpy as np +import pytest + +from scalcs import qmatlib as qml +from scalcs import scalcslib, scsim +from scalcs.samples.samples import CH82 + + +@pytest.fixture(scope="module") +def ch82(): + mec = CH82() + mec.set_eff('c', 100e-9) + return mec + + +@pytest.fixture(scope="module") +def Q(ch82): + return np.asarray(ch82.Q, dtype=float).copy() + + +class TestTransitionProbability: + r"""pi_ij = -q_ij / q_ii: given that the channel leaves state i, the + probability that the jump is to j.""" + + def test_rows_sum_to_one(self, Q): + pi = qml.transition_probability(Q) + assert pi.sum(axis=1) == pytest.approx(np.ones(Q.shape[0])) + + def test_diagonal_is_zero(self, Q): + pi = qml.transition_probability(Q) + assert np.diag(pi) == pytest.approx(np.zeros(Q.shape[0])) + + def test_matches_the_definition(self, Q): + pi = qml.transition_probability(Q) + for i in range(Q.shape[0]): + for j in range(Q.shape[0]): + if i != j: + assert pi[i, j] == pytest.approx(-Q[i, j] / Q[i, i]) + + def test_an_integer_q_matrix_is_not_truncated(self): + """A Q of integer dtype must not do integer division. + + Contrived, but it is the one way these functions could return + silently wrong numbers rather than raising. + """ + Qint = np.array([[-3, 3], [1, -1]]) + pi = qml.transition_probability(Qint) + assert pi[0, 1] == pytest.approx(1.0) + assert pi[1, 0] == pytest.approx(1.0) + + def test_does_not_modify_its_argument(self, Q): + before = Q.copy() + qml.transition_probability(Q) + assert Q == pytest.approx(before) + + +class TestTransitionFrequency: + """f_ij = p_i q_ij: transitions from i to j per second.""" + + def test_defaults_to_equilibrium_occupancies(self, Q): + assert qml.transition_frequency(Q) == pytest.approx( + qml.transition_frequency(Q, qml.pinf(Q))) + + def test_accepts_occupancies_that_are_not_at_equilibrium(self, Q): + """The reason for taking p as an argument: during a concentration + jump the occupancies are not the equilibrium ones.""" + p = np.zeros(Q.shape[0]) + p[-1] = 1.0 # everything in the last state + f = qml.transition_frequency(Q, p) + assert f[:-1, :] == pytest.approx(np.zeros_like(f[:-1, :])) + assert f[-1].sum() == pytest.approx(-Q[-1, -1]) + + def test_diagonal_is_zero(self, Q): + assert np.diag(qml.transition_frequency(Q)) == pytest.approx( + np.zeros(Q.shape[0])) + + def test_detailed_balance_holds_at_equilibrium(self, Q): + """f_ij = f_ji edge by edge, which is what microscopic reversibility + means. Not exact here: CH82's MR-constrained rate is stored rounded + to 0.66667, so its cycle closes only to about 1e-9. + """ + f = qml.transition_frequency(Q) + assert np.abs(f - f.T).max() / np.abs(f).max() < 1e-8 + + def test_does_not_modify_its_argument(self, Q): + before = Q.copy() + qml.transition_frequency(Q) + assert Q == pytest.approx(before) + + +class TestOldHomesStillWork: + """scalcslib and scsim keep their names, forwarding to qmatlib, so that + existing notebooks and scripts do not break.""" + + def test_scalcslib_forwards_both(self, Q): + assert scalcslib.transition_probability(Q) == pytest.approx( + qml.transition_probability(Q)) + assert scalcslib.transition_frequency(Q) == pytest.approx( + qml.transition_frequency(Q)) + + def test_scsim_forwards_the_probability_matrix(self, Q): + assert scsim.transition_probability(Q) == pytest.approx( + qml.transition_probability(Q)) + + def test_scalcslib_transition_frequency_still_takes_one_argument(self, Q): + """Its old signature computed pinf internally; that must still be + what a single-argument call gives.""" + assert scalcslib.transition_frequency(Q) == pytest.approx( + qml.transition_frequency(Q, qml.pinf(Q))) From 1f174c6ed811c0915f8b815be127b59296febcfe Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 12:51:18 +0000 Subject: [PATCH 3/3] Report whether a mechanism satisfies microscopic reversibility update_mr() imposes the constraint. Nothing reported whether it actually holds for the rate constants a mechanism is carrying, which is a different question and a useful one: a set that has been averaged, transcribed, or written by another program can fail it, and the size of the failure says how much to trust the set. cycle_ratios() [(states, forward/backward), ...], one per cycle reversibility_error() max |ratio - 1| over all cycles, 0.0 if none check_mr already returned the two products for one cycle, and these agree with it - there is a test pinning that so they cannot drift. What they add is the whole-mechanism summary, and failing loudly: check_mr looks each step up with a scan that contributes nothing when no rate matches, so a cycle naming a transition the mechanism does not have comes back as a perfectly closed cycle. These use rate(), which raises. Scale, for reading the number. The two ten-state datasets behind Colquhoun & Lape (2012): single fitted set 1.4e-07 (stored as 32-bit floats) averaged set 4.2e-02 (ratios 0.9576 and 0.9671) Averaging several fitted sets does not preserve cycle conditions, and that is what a violation of percents rather than parts in 10^7 looks like. CH82 gives 5e-06, because its constrained rate is stored rounded to 0.66667. Co-Authored-By: Claude Opus 5 --- scalcs/mechanism.py | 45 +++++++++++++++++++++++++++ tests/test_mechanism.py | 67 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 112 insertions(+) diff --git a/scalcs/mechanism.py b/scalcs/mechanism.py index 0b16e53..68e405f 100644 --- a/scalcs/mechanism.py +++ b/scalcs/mechanism.py @@ -840,6 +840,51 @@ def permutation(self, order): % (len(order), self.k, ', '.join(self.state_names))) return [self.state_index(name) for name in order] + # ----------------------------------------------------------------- + # Checking microscopic reversibility, as opposed to imposing it + # ----------------------------------------------------------------- + + def cycle_ratios(self): + """Forward/backward rate-constant product for each cycle. + + Returns a list of ``(states, ratio)``, one per cycle. Microscopic + reversibility requires every ratio to be 1. + + Unlike :meth:`check_mr`, a cycle that names a transition the + mechanism does not have raises rather than quietly contributing a + factor of 1, which would report a broken mechanism as reversible. + + The rate constants are taken at unit effector, so the concentration + cancels: a cycle crosses equally many association and dissociation + steps. + """ + out = [] + for cycle in self.Cycles: + states = list(cycle.states) + following = states[1:] + states[:1] + forward = backward = 1.0 + for a, b in zip(states, following): + forward *= self.rate(a, b).unit_rate() + backward *= self.rate(b, a).unit_rate() + out.append((tuple(states), float(forward / backward))) + return out + + def reversibility_error(self): + """How far the rate constants are from microscopic reversibility. + + ``max |forward/backward - 1|`` over all cycles; 0.0 for a mechanism + with no cycles, which cannot violate it. + + Expect roughly 1e-7 from a set of rate constants that was fitted + under the constraint and stored as 32-bit floats. Percents mean the + constraint was not applied - averaging several fitted sets, for + instance, does not preserve it. + """ + ratios = self.cycle_ratios() + if not ratios: + return 0.0 + return max(abs(ratio - 1.0) for _, ratio in ratios) + def set_eff(self, eff, val): self.set_effdict({eff:val}) diff --git a/tests/test_mechanism.py b/tests/test_mechanism.py index 5a8bb08..6e673da 100644 --- a/tests/test_mechanism.py +++ b/tests/test_mechanism.py @@ -471,3 +471,70 @@ def test_permutation_rejects_an_incomplete_order(self, ch82): def test_permutation_rejects_an_unknown_state(self, ch82): with pytest.raises(KeyError, match='A3R'): ch82.permutation(['R', 'AR', 'A2R', 'AR*', 'A3R']) + + +# --------------------------------------------------------------------------- +# Checking microscopic reversibility rather than imposing it +# --------------------------------------------------------------------------- + +class TestMicroscopicReversibility: + """Whether the rate constants a mechanism actually holds satisfy + microscopic reversibility. + + update_mr() imposes the constraint; these report whether it holds. That + is a different question, and a useful one: a set of rate constants that + has been averaged, or transcribed, or read from a file written by some + other program, can fail it, and the size of the failure says how much to + trust the set. + """ + + def test_ch82_is_reversible(self, ch82): + """Not exactly, though. CH82's MR-constrained rate is stored rounded + to 0.66667 rather than 2/3, which is a relative error of 5e-6, and + that is what comes back. The point is the order of magnitude: a set + that genuinely violates reversibility is out by percents. + """ + assert ch82.reversibility_error() < 1e-4 + + def test_a_mechanism_with_no_cycles_cannot_violate_it(self, cco): + """CCO is a chain. With no cycle there is nothing to close.""" + assert cco.Cycles == [] + assert cco.reversibility_error() == 0.0 + + def test_cycle_ratios_reports_one_ratio_per_cycle(self, ch82): + ratios = ch82.cycle_ratios() + assert len(ratios) == len(ch82.Cycles) + states, ratio = ratios[0] + assert set(states) == set(ch82.Cycles[0].states) + assert ratio == pytest.approx(1.0, abs=1e-4) + + def test_it_agrees_with_check_mr(self, ch82): + """check_mr is the existing, lower-level way to ask. The two must + not drift apart.""" + for cycle, (_, ratio) in zip(ch82.Cycles, ch82.cycle_ratios()): + fprod, bprod = ch82.check_mr(cycle) + assert ratio == pytest.approx(fprod / bprod) + + def test_a_broken_rate_constant_is_detected(self): + """Double one rate in the cycle and the cycle stops closing. + + Built fresh rather than from the fixture, which is module-scoped and + shared with every other test here. + """ + mec = CH82() + before = mec.reversibility_error() + rate = mec.rate('A2R*', 'AR*') + rate.rateconstants = rate.rateconstants * 2.0 + + assert before < 1e-4 + assert mec.reversibility_error() > 0.1 + + def test_a_cycle_naming_a_transition_that_does_not_exist_raises(self): + """A cycle must be made of real transitions. Silently treating a + missing step as a factor of 1 would report a broken mechanism as + perfectly reversible. + """ + mec = CH82() + mec.Cycles = [Cycle(['R', 'A2R*', 'AR', 'A2R'])] # R -> A2R* is not a step + with pytest.raises(KeyError): + mec.reversibility_error()