Feature/mechanism lookup

scalcs/mechanism.py

@@ -780,6 +780,111 @@ 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]
+
+    # -----------------------------------------------------------------
+    # 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})
 

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).

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):
     """

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):

tests/test_mechanism.py

@@ -410,3 +410,131 @@ 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'])
+
+
+# ---------------------------------------------------------------------------
+# 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()