From 335b905ff1881b4502c77582ce72e9d9338cb984 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 16:27:03 +0000 Subject: [PATCH 1/3] Draw a mechanism as a reaction scheme scalcs.schemes draws a mechanism in the style of Colquhoun & Lape (2012), Figs. 3, 7 and 8: states on a grid, each reversible step a close pair of harpoon arrows with the barb of each on the outside of the pair, and a rate constant, transition probability or equilibrium transition frequency written against each direction. schemes.draw(mec, values='probability', conc=10e-6, numbers=True) The arrows are the rates of the mechanism. They are never inferred from where the states are placed: joining neighbours on a grid is the obvious shortcut, and it once drew the ten-state flip mechanism with thirteen connections instead of eleven, joining open states that do not interconvert. The placement comes from layout= if given, otherwise mec.layout - the drawing stored with a .mec file - otherwise a ready-made layout matched on the state names, which covers every mechanism in scalcs.samples. Each arrow sits to its own left and carries its barb on its own left. The two arrows of a pair run in opposite directions, so that puts both barbs on the outside at any angle, with no notion of a forward direction needed. Text and line widths shrink with the axes, relative to the scale of the published Figs. 7 and 8, so a scheme in a small panel keeps the proportions of the published figure instead of having its names overrun its arrows. They are not enlarged in a bigger panel, where neighbouring numbers crowd. Labels go where the published figures put them, beside each arrow near its tail. On a diagonal they are anchored by the corner nearest the arrow, so none of the text lies across the shaft. A layout that leaves a state out, names one that does not exist, puts two in one place, or runs a step's arrows over a third state - which would read as two steps through it - is an error naming the states involved. geometry() returns every position without drawing, so all but the rendering is tested without matplotlib; the tests check the edges drawn against the rates, the barbs against the pair, clearance round the state names, and the values against qmatlib. Every sample mechanism is drawn. Co-Authored-By: Claude Opus 5.5 --- README.md | 27 ++ scalcs/schemes.py | 595 ++++++++++++++++++++++++++++++++++++++++++ tests/test_schemes.py | 425 ++++++++++++++++++++++++++++++ 3 files changed, 1047 insertions(+) create mode 100644 scalcs/schemes.py create mode 100644 tests/test_schemes.py diff --git a/README.md b/README.md index c5c5e53..abacdae 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,7 @@ Modules | `popen` | Equilibrium open probability | | `mechanism` | Ion channel mechanism (state/rate) representation | | `scalcsio` | HJCFIT / DC_PyPs file I/O | +| `schemes` | Draw a mechanism as a reaction scheme (needs `scalcs[plot]`) | ### First-latency pdf (`firstlatency`) @@ -52,6 +53,28 @@ eigvals, g00, g10, g11 = fl.gamma_coefficients(tres, phi_shut, mec) f_exact = fl.exact_pdf(t, tres, roots, areas, eigvals, g00, g10, g11) ``` +### Reaction schemes (`schemes`) + +Draws a mechanism in the style of Colquhoun & Lape (2012), Figs. 3, 7 and 8: +states on a grid, each reversible step a close pair of harpoon arrows, with a +rate constant, transition probability or equilibrium transition frequency +written against each direction. + +```python +from scalcs import schemes +from scalcs.samples.samples import GlyR_flip + +schemes.draw(GlyR_flip(), values='probability', conc=10e-6, numbers=True) +``` + +The arrows drawn are the rates of the mechanism, never inferred from where the +states are placed. The placement comes from `layout=` if given, otherwise from +the drawing stored with a mechanism read from a .mec file, otherwise from a +ready-made layout covering every mechanism in `scalcs.samples`. + +Drawing needs matplotlib, which is optional: `pip install scalcs[plot]`. +`schemes.geometry()` computes every position without it. + References ---------- @@ -79,6 +102,10 @@ A Q-Matrix Cookbook. In: Single-channel recording. 2nd ed. (Eds: Sakmann B, Neher E) Plenum Press, New York, pp. 589-633. +CL2012: Colquhoun D, Lape R (2012) +Allosteric coupling in ligand-gated ion channels. +J Gen Physiol 140, 599-612. + CHS96: Colquhoun D, Hawkes AG, Srodzinski K (1996) Joint distributions of apparent open and shut times of single-ion channels and maximum likelihood fitting of mechanisms. diff --git a/scalcs/schemes.py b/scalcs/schemes.py new file mode 100644 index 0000000..3ec551c --- /dev/null +++ b/scalcs/schemes.py @@ -0,0 +1,595 @@ +"""Draw kinetic mechanisms as reaction schemes. + +The style is the one Colquhoun & Lape (2012) use for Figs. 3, 7 and 8: the +states set out on a grid, each reversible step drawn as a pair of harpoon +arrows lying close together, with a number written against each direction - +a rate constant, a transition probability or a transition frequency. + + from scalcs import schemes + from scalcs.samples.samples import CH82 + + schemes.draw(CH82(), values='probability', conc=100e-9, numbers=True) + +Which states are connected is taken from the rates of the mechanism and from +nothing else. Where the states go - the layout - is the only thing the caller +supplies, and usually not even that: + + 1. ``layout=`` given explicitly, as rows of state names (top row first, + None for a gap) or as ``{name: (column, row)}``; + 2. otherwise ``mec.layout``, the drawing stored with a mechanism read from + a .mec file; + 3. otherwise a ready-made layout from :data:`LAYOUTS`, matched on the names + of the states, which covers every mechanism in ``scalcs.samples``. + +Inferring the connections from the layout instead - joining neighbours on the +grid - is the obvious shortcut and the wrong one: it once drew the ten-state +flip mechanism with thirteen connections instead of eleven, joining open +states that do not interconvert. + +Everything except :func:`draw` is plain arithmetic, so the geometry can be +computed, inspected and tested without matplotlib, which is an optional +dependency: ``pip install scalcs[plot]``. +""" +import math +import re +from dataclasses import dataclass + +import numpy as np + +from scalcs import qmatlib as qml + + +# --------------------------------------------------------------------------- +# Style +# --------------------------------------------------------------------------- + +# Inches of figure per data unit at which the style's font sizes were tuned: +# Figs. 7 and 8 of Colquhoun & Lape (2012) set as two panels of a 7.6 x 9.4 +# inch figure. +REFERENCE_SCALE = 0.734 + +@dataclass +class Style: + """How a scheme is drawn. + + The defaults were measured from the published Figs. 3, 7 and 8 of + Colquhoun & Lape (2012). Every length is a fraction of `spacing`, the + distance between neighbouring states in data units, so changing + `spacing` rescales the whole drawing without retuning anything else. + """ + spacing: float = 1.4 # between neighbouring columns + row_spacing: float = None # between rows; None means the same as spacing. + # The published Fig. 3 uses 1.5 times spacing. + + pair_offset: float = 0.026 # half the gap between the two arrows of a pair + label_gap: float = 0.30 # clearance left round each state name + head_length: float = 0.20 + head_width: float = 0.046 + label_pad: float = 0.065 # from an arrow to the number written by it + line_width: float = 1.1 + + # Font sizes and the line width are those of the published figures, at + # REFERENCE_SCALE. Lengths are in data units but text is in points, so a + # scheme drawn into a smaller panel would have its names overrun the + # arrows; with scale_text they shrink with the axes instead. They are + # never enlarged. + state_fontsize: float = 19 + value_fontsize: float = 9 + number_fontsize: float = 9 + scale_text: bool = True + + open_colour: str = '#ed1c24' # the paper's red + shut_colour: str = 'black' + arrow_colour: str = 'black' + number_colour: str = '#1c9ad6' # the paper's blue state numbering + + def __post_init__(self): + if self.row_spacing is None: + self.row_spacing = self.spacing + + +# --------------------------------------------------------------------------- +# Ready-made layouts +# --------------------------------------------------------------------------- + +# Rows of state names, top row first, None for an empty cell. They follow the +# convention DC used in the drawings stored in his .mec files: binding from +# left to right along the resting row, a change of conformation downwards, +# open states at the bottom. +LAYOUTS = { + 'CO': [['C'], + ['O']], + 'RFO': [['R', 'F'], + [None, 'F*']], + 'RFFOO': [['R', 'F', 'FF'], + [None, 'F*', 'FF*']], + 'CCO': [['R', 'AR'], + [None, 'AR*']], + 'CCOB': [['R', 'AR', None], + [None, 'AR*', 'AR*B']], + 'CCOD': [['R', 'AR', None], + [None, 'AR*', 'AD']], + 'CH82': [['R', 'AR', 'A2R'], + [None, 'AR*', 'A2R*']], + 'CHME97': [['R', 'AR', 'A2R', 'A2D'], + [None, None, 'A2R*', None]], + 'GlyR_flip': [['R', 'AR', 'A2R', 'A3R'], + [None, 'AF', 'A2F', 'A3F'], + [None, 'AF*', 'A2F*', 'A3F*']], + 'AChR_diamond': [[None, 'AR*a', None, None], + [None, 'ARa', None, None], + ['R', None, 'A2R', 'A2R*'], + [None, 'ARb', None, None], + [None, 'AR*b', None, None]], + 'AChR_diamond_desens': [[None, 'AR*a', None, None, None], + [None, 'ARa', None, None, None], + ['R', None, 'A2R', 'A2R*', 'A2D'], + [None, 'ARb', None, None, None], + [None, 'AR*b', None, None, None]], + # Every state joined to every other, so the diagonals cross. Tried as a + # triangle with a state at its centre, which avoids the crossing but + # leaves three steps in a 53-degree corner, where their numbers collide. + # The square leaves 45 degrees between neighbouring steps. + 'fully_connected_cycle': [['C1', 'C2'], + ['O1', 'O2']], + 'six_cycles_mec': [['A', 'B', 'C', 'D'], + ['E', 'F', 'G', 'H'], + ['I', 'J', 'K', 'L']], +} + + +# --------------------------------------------------------------------------- +# Text +# --------------------------------------------------------------------------- + +def format_number(value): + """A number written the way the published figures write them. + + Plain decimals with trailing zeros removed, 1.0 as a bare 1; thousands + printed whole, which is what makes a frequency of 6471.73 per second come + out as 6472. Values of 1e5 and above - association rate constants, in + practice - get a power of ten, as do values too small for seven decimal + places. + """ + value = float(value) + if value == 0: + return '0' + if abs(value) >= 1e5: + return _power_of_ten(value) + if abs(value) >= 1000: + return '%.0f' % value + text = '%.3g' % value + if 'e' in text: # below 1e-4 + text = '%.7f' % value + if '.' in text: + text = text.rstrip('0').rstrip('.') + if text in ('0', '-0'): # below 5e-8 + return _power_of_ten(value) + return text + + +def _power_of_ten(value): + mantissa, exponent = ('%.3g' % value).split('e') + return r'$%s\times10^{%d}$' % (mantissa, int(exponent)) + + +def state_label(name): + """A state name as mathtext: digits after a letter become subscripts and + an asterisk a superscript, so 'A2F*' is written A2F* with a subscript 2 + and a raised star. Names with any other character are returned as they + are, rather than risk mathtext misreading them. + """ + if not re.fullmatch(r'[A-Za-z0-9*]+', name): + return name + body = re.sub(r'([A-Za-z])(\d+)', r'\1_{\2}', name) + body = body.replace('*', '^{*}') + return r'$\mathrm{%s}$' % body + + +def state_colours(mec, colours=None, style=None): + """{name: colour}: open states in the style's open colour, the rest in + its shut colour, then any given in `colours` - which is how the flipped + conformation of Figs. 7 and 8 gets its green.""" + style = style or Style() + out = {state.name: (style.open_colour if state.statetype == 'A' + else style.shut_colour) + for state in mec.States} + if colours: + _check_names(mec, colours, 'colours') + out.update(colours) + return out + + +def state_numbers(mec, numbers): + """{name: number} to write by each state. + + False writes none. True numbers the states in Q-matrix order, from 1, + which for a mechanism read from a .mec file is the numbering of the file. + A dict gives any other numbering - a paper's own, which need not be + SCALCS's: the published flip mechanism numbers two of its states the + other way round. + """ + if not numbers: + return {} + if numbers is True: + return {name: i + 1 for i, name in enumerate(mec.state_names)} + _check_names(mec, numbers, 'numbers') + return dict(numbers) + + +def _check_names(mec, given, what): + unknown = sorted(set(given) - set(mec.state_names)) + if unknown: + raise ValueError('%s names %s, which %s not a state of this mechanism ' + '(its states are %s)' + % (what, ', '.join(unknown), + 'is' if len(unknown) == 1 else 'are', + ', '.join(mec.state_names))) + + +# --------------------------------------------------------------------------- +# Geometry +# --------------------------------------------------------------------------- + +@dataclass +class Arrow: + """One direction of one step: source -> target. + + `barb` is the filled triangle at the head: the tip of the arrow, the + point of the barb, and where the barb rejoins the shaft. + """ + source: str + target: str + start: tuple + end: tuple + barb: tuple + label: str = None + label_xy: tuple = None + label_ha: str = 'center' + label_va: str = 'center' + + +@dataclass +class SchemeGeometry: + positions: dict # {name: (x, y)} in data units, y upwards + arrows: list # one Arrow per rate of the mechanism + + +def geometry(mec, layout=None, values=None, conc=None, style=None): + """Where every state, arrow and number goes. + + Parameters + ---------- + mec : Mechanism + layout : list of rows, or dict, optional + Rows of state names, top row first, None for an empty cell; or + ``{name: (column, row)}`` with row 0 at the top. If omitted, the + layout stored with the mechanism is used, then a ready-made one from + :data:`LAYOUTS`. + values : None, 'rates', 'probability', 'frequency', or a function + What to write against each arrow. 'rates' writes the rate constants; + 'probability' the transition probabilities pi_ij, Eq. 2 of + Colquhoun & Lape (2012); 'frequency' the transition frequencies at + equilibrium, f_ij = p_i q_ij, their Eq. 3. A function + ``f(source, target)`` may return any text, or None to leave that + direction unlabelled. + conc : float, optional + Agonist concentration (M) for 'probability' and 'frequency'. Needed + if any rate depends on concentration. The mechanism is returned to + the concentration it had. + style : Style, optional + + Returns + ------- + SchemeGeometry + """ + style = style or Style() + cells = _resolve_layout(mec, layout) + positions = {name: (column * style.spacing, -row * style.row_spacing) + for name, (column, row) in cells.items()} + _check_no_step_crosses_a_state(mec, positions, style) + + label = _labeller(mec, values, conc) + arrows = [_arrow(source, target, positions, style, + label(source, target) if label else None) + for source, target in mec.rates_by_edge()] + return SchemeGeometry(positions, arrows) + + +def _arrow(source, target, positions, style, label): + """Every arrow sits to its own left of the line between the two states + and carries its barb on its own left. The two arrows of a pair run in + opposite directions, so that puts each barb on the outside of the pair, + at whatever angle the step runs, with no need to decide which direction + is 'forward'. + """ + (xa, ya), (xb, yb) = positions[source], positions[target] + length = math.hypot(xb - xa, yb - ya) + ux, uy = (xb - xa) / length, (yb - ya) / length # along the arrow + nx, ny = -uy, ux # to its left + + gap = style.label_gap * style.spacing + offset = style.pair_offset * style.spacing + start = (xa + ux * gap + nx * offset, ya + uy * gap + ny * offset) + end = (xb - ux * gap + nx * offset, yb - uy * gap + ny * offset) + + head = style.head_length * style.spacing + width = style.head_width * style.spacing + back = (end[0] - ux * head, end[1] - uy * head) + barb = (end, + (back[0] + nx * width, back[1] + ny * width), + (end[0] - ux * head * 0.40, end[1] - uy * head * 0.40)) + + arrow = Arrow(source, target, start, end, barb, label) + if label is not None: + _place_label(arrow, ux, uy, nx, ny, style) + return arrow + + +# An arrow within about 25 degrees of horizontal or vertical has its label +# placed as Figs. 7 and 8 place them; anything steeper is a diagonal. +_NEARLY_ALIGNED = 0.9 + + +def _place_label(arrow, ux, uy, nx, ny, style): + """Beside the arrow, on the outside, near its tail - where Figs. 7 and 8 + put them. + + Horizontal: centred a fifth of the way along, above or below. Vertical: + beside the arrow, starting level with the tail and reading towards the + head. Diagonal: a fifth of the way along, anchored by the corner of the + text nearest the arrow, so that all of the text lies outside the pair - + centred, half of it would lie across the shaft. + """ + (sx, sy), (ex, ey) = arrow.start, arrow.end + pad = style.label_pad * style.spacing + if abs(ux) >= _NEARLY_ALIGNED: + along = 0.20 + arrow.label_ha = 'center' + arrow.label_va = 'bottom' if ny > 0 else 'top' + elif abs(uy) >= _NEARLY_ALIGNED: + along = 0.04 + arrow.label_ha = 'left' if nx > 0 else 'right' + arrow.label_va = 'top' if uy < 0 else 'bottom' + else: + along = 0.20 + arrow.label_ha = 'left' if nx > 0 else 'right' + arrow.label_va = 'bottom' if ny > 0 else 'top' + arrow.label_xy = (sx + along * (ex - sx) + nx * pad, + sy + along * (ey - sy) + ny * pad) + + +# --------------------------------------------------------------------------- +# Layout +# --------------------------------------------------------------------------- + +def _resolve_layout(mec, layout): + if layout is not None: + cells = _cells(layout) + elif getattr(mec, 'layout', None) is not None: + cells = _compact(mec.layout) + else: + cells = _ready_made(mec) + if cells is None: + raise ValueError( + 'no layout for this mechanism. Pass one as layout=[[...], ' + '[...]], rows of state names with the top row first and None ' + 'for an empty cell, or as layout={name: (column, row)}. Its ' + 'states are %s.' % ', '.join(mec.state_names)) + _check_layout(mec, cells) + return cells + + +def _cells(layout): + """{name: (column, row)} from either form of layout.""" + if isinstance(layout, dict): + return {name: tuple(cell) for name, cell in layout.items()} + cells = {} + for row, names in enumerate(layout): + for column, name in enumerate(names): + if name is None: + continue + if name in cells: + raise ValueError('%s appears twice in the layout' % name) + cells[name] = (column, row) + return cells + + +def _compact(board_layout): + """A layout read from a .mec drawing, on a grid of unit steps. + + The drawings put states on every other cell, with the connecting lines + in between. The common step is divided out of each axis, so a drawing + and the equivalent rows of names give the same positions. + """ + columns = [column for column, _ in board_layout.values()] + rows = [row for _, row in board_layout.values()] + left, top = min(columns), min(rows) + step_x = math.gcd(*[c - left for c in columns]) or 1 + step_y = math.gcd(*[r - top for r in rows]) or 1 + return {name: ((column - left) // step_x, (row - top) // step_y) + for name, (column, row) in board_layout.items()} + + +def _ready_made(mec): + """The entry of LAYOUTS with exactly this mechanism's state names.""" + names = set(mec.state_names) + for rows in LAYOUTS.values(): + cells = _cells(rows) + if set(cells) == names: + return cells + return None + + +def _check_layout(mec, cells): + names = set(mec.state_names) + missing = [n for n in mec.state_names if n not in cells] + if missing: + raise ValueError('the layout gives no position for %s' + % ', '.join(missing)) + _check_names(mec, cells, 'the layout') + taken = {} + for name, cell in cells.items(): + if cell in taken: + raise ValueError('%s and %s are in the same place, %s' + % (taken[cell], name, cell)) + taken[cell] = name + + +def _check_no_step_crosses_a_state(mec, positions, style): + """A step whose arrows pass over a third state would be read as two + steps through it - a transition the mechanism does not have, drawn.""" + clearance = style.label_gap * style.spacing + done = set() + for source, target in mec.rates_by_edge(): + pair = frozenset((source, target)) + if pair in done: + continue + done.add(pair) + a, b = np.array(positions[source]), np.array(positions[target]) + for name, point in positions.items(): + if name in pair: + continue + t = np.dot(np.array(point) - a, b - a) / np.dot(b - a, b - a) + if 0 < t < 1: + nearest = a + t * (b - a) + if np.hypot(*(np.array(point) - nearest)) < clearance: + raise ValueError( + 'the step between %s and %s would be drawn through ' + '%s; move one of them in the layout' + % (source, target, name)) + + +# --------------------------------------------------------------------------- +# Values written against the arrows +# --------------------------------------------------------------------------- + +VALUES = ('rates', 'probability', 'frequency') + + +def _labeller(mec, values, conc): + if values is None: + return None + if callable(values): + return values + if values not in VALUES: + raise ValueError("values must be one of %s, or a function " + "f(source, target) returning text; got %r" + % (', '.join(repr(v) for v in VALUES), values)) + if values == 'rates': + rates = mec.rates_by_edge() + return lambda a, b: format_number(rates[(a, b)].unit_rate()) + + Q = _q_at(mec, conc) + if values == 'probability': + matrix = qml.transition_probability(Q) + else: + matrix = qml.transition_frequency(Q) + return lambda a, b: format_number( + matrix[mec.state_index(a), mec.state_index(b)]) + + +def _q_at(mec, conc): + """The Q matrix at `conc`, leaving the mechanism as it was.""" + if conc is None: + if mec._effdict: + raise ValueError( + 'the rates of this mechanism depend on concentration, so ' + 'transition probabilities and frequencies need one: pass ' + 'conc= in M') + return np.array(mec.Q, dtype=float, copy=True) + saved = dict(mec._effdict) + try: + mec.set_eff('c', conc) + return np.array(mec.Q, dtype=float, copy=True) + finally: + mec.set_effdict(saved) + + +# --------------------------------------------------------------------------- +# Drawing +# --------------------------------------------------------------------------- + +def draw(mec, ax=None, layout=None, values=None, conc=None, colours=None, + numbers=False, style=None): + """Draw the mechanism; return the axes. + + Parameters are those of :func:`geometry`, plus: + + ax : matplotlib Axes, optional + Drawn into a new figure, sized to the scheme, if not given. Text and + lines are scaled to the size of the axes at the time of drawing; see + :class:`Style`. + colours : dict, optional + {name: colour} for any states not to be drawn in the default open or + shut colour. + numbers : bool or dict + State numbers to write in blue beside each state; see + :func:`state_numbers`. + + Requires matplotlib (``pip install scalcs[plot]``). + """ + import matplotlib.pyplot as plt + + style = style or Style() + scheme = geometry(mec, layout=layout, values=values, conc=conc, + style=style) + colour = state_colours(mec, colours, style) + number = state_numbers(mec, numbers) + + xs = [x for x, _ in scheme.positions.values()] + ys = [y for _, y in scheme.positions.values()] + margin = 0.6 * style.spacing + width = max(xs) - min(xs) + 2 * margin + height = max(ys) - min(ys) + 2 * margin + if ax is None: + # axes filling a figure sized so the scheme is at the reference scale + figure = plt.figure(figsize=(width * REFERENCE_SCALE, + height * REFERENCE_SCALE)) + ax = figure.add_axes([0, 0, 1, 1]) + + # Limits and aspect first, so the size of the axes - and from it the + # scale of the text - is known before anything is drawn. Equal aspect, + # or the harpoons shear and a pair stops looking parallel. + ax.set_aspect('equal') + ax.set_xlim(min(xs) - margin, max(xs) + margin) + ax.set_ylim(min(ys) - margin, max(ys) + margin) + ax.axis('off') + scale = _text_scale(ax, width) if style.scale_text else 1.0 + + for arrow in scheme.arrows: + (sx, sy), (ex, ey) = arrow.start, arrow.end + ax.plot([sx, ex], [sy, ey], color=style.arrow_colour, + lw=style.line_width * scale, solid_capstyle='butt', zorder=2) + ax.add_patch(plt.Polygon(arrow.barb, closed=True, + facecolor=style.arrow_colour, + edgecolor='none', zorder=3)) + if arrow.label is not None: + ax.text(*arrow.label_xy, arrow.label, ha=arrow.label_ha, + va=arrow.label_va, fontsize=style.value_fontsize * scale) + + for name, (x, y) in scheme.positions.items(): + ax.text(x, y, state_label(name), ha='center', va='center', + fontsize=style.state_fontsize * scale, color=colour[name]) + if name in number: + ax.text(x + 0.093 * style.spacing, y - 0.079 * style.spacing, + str(number[name]), ha='left', va='top', + fontsize=style.number_fontsize * scale, + color=style.number_colour) + return ax + + +def _text_scale(ax, width): + """Inches per data unit of `ax`, relative to REFERENCE_SCALE, and never + more than 1. + + Text shrinks with a panel smaller than the reference, so names do not + overrun the arrows; it is not enlarged in a bigger one, where it would + only make the numbers of neighbouring steps crowd each other. + + apply_aspect() resolves the box that equal aspect leaves the axes, so the + horizontal scale is also the vertical one. + """ + ax.apply_aspect() + inches = ax.get_position().width * ax.figure.get_figwidth() + return min(1.0, inches / width / REFERENCE_SCALE) diff --git a/tests/test_schemes.py b/tests/test_schemes.py new file mode 100644 index 0000000..ecf8a4a --- /dev/null +++ b/tests/test_schemes.py @@ -0,0 +1,425 @@ +"""Tests for scalcs.schemes, the kinetic-scheme drawer. + +Most of these test the geometry - where the states go and where every arrow +starts, ends and points - which is plain arithmetic and needs no matplotlib. +Only TestDrawing renders anything. + +The drawing style is the one Colquhoun & Lape (2012) use for Figs. 3, 7 and +8: closely paired harpoon arrows, the single barb of each on the outside of +its pair. +""" +import inspect +import math +import subprocess +import sys + +import numpy as np +import pytest + +from scalcs import qmatlib as qml +from scalcs import schemes +from scalcs.mechanism import Mechanism, Rate, State +from scalcs.samples import samples +from scalcs.samples.samples import CCO, CH82, GlyR_flip + + +def _sample_mechanisms(): + """Every mechanism in scalcs.samples that builds without arguments, so a + sample added later without a layout fails a test here.""" + out = [] + for name, factory in inspect.getmembers(samples, inspect.isfunction): + try: + mec = factory() + except TypeError: + continue + if isinstance(mec, Mechanism): + out.append((name, mec)) + return out + + +SAMPLES = _sample_mechanisms() + + +def _connections(mec): + return {frozenset(edge) for edge in mec.rates_by_edge()} + + +def _pairs(geometry): + by_edge = {} + for arrow in geometry.arrows: + by_edge.setdefault(frozenset((arrow.source, arrow.target)), []).append(arrow) + return [arrows for arrows in by_edge.values() if len(arrows) == 2] + + +def _unit(start, end): + dx, dy = end[0] - start[0], end[1] - start[1] + length = math.hypot(dx, dy) + return dx / length, dy / length + + +def _side(point, start, end): + """Signed distance of `point` to the left of the line start -> end.""" + ux, uy = _unit(start, end) + return (point[0] - start[0]) * -uy + (point[1] - start[1]) * ux + + +def _two_state(names=('X', 'Y')): + a, b = State('B', names[0], 0.0), State('A', names[1], 50e-12) + return Mechanism([Rate(10.0, a, b, name='on'), Rate(100.0, b, a, name='off')]) + + +# --------------------------------------------------------------------------- + +class TestTopologyComesFromTheMechanism: + """The connections drawn are the rates of the mechanism, never inferred + from where the states have been put.""" + + def test_one_arrow_per_rate(self): + mec = GlyR_flip() + assert len(schemes.geometry(mec).arrows) == len(mec.Rates) + + def test_the_flip_mechanism_has_eleven_connections_not_thirteen(self): + """The open states A F*, A2F*, A3F* sit side by side but do not + interconvert. A drawer that joined neighbours on a grid drew two + binding steps between them that the mechanism does not have.""" + drawn = {frozenset((a.source, a.target)) + for a in schemes.geometry(GlyR_flip()).arrows} + assert len(drawn) == 11 + assert frozenset(('AF*', 'A2F*')) not in drawn + assert frozenset(('A2F*', 'A3F*')) not in drawn + + @pytest.mark.parametrize('name, mec', SAMPLES, ids=[n for n, _ in SAMPLES]) + def test_arrows_join_exactly_the_states_the_rates_join(self, name, mec): + drawn = {(a.source, a.target) for a in schemes.geometry(mec).arrows} + assert drawn == set(mec.rates_by_edge()) + + +class TestPairs: + """The two arrows of a reversible step, at whatever angle it runs.""" + + @pytest.fixture(scope='class', params=['CH82', 'GlyR_flip', + 'AChR_diamond', + 'fully_connected_cycle']) + def geometry(self, request): + return schemes.geometry(getattr(samples, request.param)()) + + def test_the_two_arrows_run_in_opposite_directions(self, geometry): + for first, second in _pairs(geometry): + ux, uy = _unit(first.start, first.end) + vx, vy = _unit(second.start, second.end) + assert ux * vx + uy * vy == pytest.approx(-1.0) + + def test_they_are_separated_by_twice_the_pair_offset(self, geometry): + style = schemes.Style() + expected = 2 * style.pair_offset * style.spacing + for first, second in _pairs(geometry): + assert abs(_side(second.start, first.start, first.end)) == \ + pytest.approx(expected) + + def test_each_barb_is_on_the_outside_of_its_pair(self, geometry): + """The barb of each arrow points away from the other arrow.""" + for first, second in _pairs(geometry): + for arrow, other in ((first, second), (second, first)): + tip_of_barb = arrow.barb[1] + other_side = _side(other.start, arrow.start, arrow.end) + barb_side = _side(tip_of_barb, arrow.start, arrow.end) + assert barb_side * other_side < 0 + + def test_arrows_stop_short_of_the_state_labels(self, geometry): + style = schemes.Style() + clearance = style.label_gap * style.spacing + for arrow in geometry.arrows: + for point, name in ((arrow.start, arrow.source), + (arrow.end, arrow.target)): + centre = geometry.positions[name] + distance = math.hypot(point[0] - centre[0], point[1] - centre[1]) + assert distance >= clearance * 0.99 + + + def test_labels_on_a_diagonal_are_anchored_away_from_the_arrow(self): + """Centred above a 45-degree shaft, half of the text would lie on the + arrow. A diagonal label is anchored by the corner nearest its arrow, + so all of it lies outside the pair.""" + geometry = schemes.geometry(samples.AChR_diamond(), + values=lambda a, b: 'x') + diagonals = 0 + for arrow in geometry.arrows: + ux, uy = _unit(arrow.start, arrow.end) + if min(abs(ux), abs(uy)) < 0.5: + continue + diagonals += 1 + nx, ny = -uy, ux + assert arrow.label_ha == ('left' if nx > 0 else 'right') + assert arrow.label_va == ('bottom' if ny > 0 else 'top') + assert diagonals == 8 + + def test_labels_on_grid_steps_are_placed_as_in_the_paper(self): + """Horizontal: centred, above or below. Vertical: beside the arrow, + reading from its tail towards its head.""" + geometry = schemes.geometry(CH82(), values=lambda a, b: 'x') + by = {(a.source, a.target): a for a in geometry.arrows} + assert (by[('R', 'AR')].label_ha, by[('R', 'AR')].label_va) == ('center', 'bottom') + assert (by[('AR', 'R')].label_ha, by[('AR', 'R')].label_va) == ('center', 'top') + assert (by[('AR', 'AR*')].label_ha, by[('AR', 'AR*')].label_va) == ('left', 'top') + assert (by[('AR*', 'AR')].label_ha, by[('AR*', 'AR')].label_va) == ('right', 'bottom') + + +class TestLayout: + + def test_rows_of_names_are_given_top_row_first(self): + geometry = schemes.geometry(CCO(), layout=[['R', 'AR'], [None, 'AR*']]) + assert geometry.positions['AR'][0] == geometry.positions['AR*'][0] + assert geometry.positions['AR'][1] > geometry.positions['AR*'][1] + + def test_a_dict_of_column_and_row_works_the_same_way(self): + rows = schemes.geometry(CCO(), layout=[['R', 'AR'], [None, 'AR*']]) + cells = schemes.geometry(CCO(), layout={'R': (0, 0), 'AR': (1, 0), + 'AR*': (1, 1)}) + assert rows.positions == cells.positions + + def test_the_drawing_stored_with_a_mec_file_is_used(self): + """A .mec board puts states on every other cell. That spacing is + taken out, so a board layout and the equivalent rows agree.""" + mec = CCO() + mec.layout = {'R': (0, 0), 'AR': (2, 0), 'AR*': (2, 2)} + from_file = schemes.geometry(mec) + mec.layout = None + given = schemes.geometry(mec, layout=[['R', 'AR'], [None, 'AR*']]) + assert from_file.positions == given.positions + + def test_a_layout_given_explicitly_wins_over_the_stored_one(self): + mec = CCO() + mec.layout = {'R': (0, 0), 'AR': (2, 0), 'AR*': (2, 2)} + geometry = schemes.geometry(mec, layout=[['AR*', 'AR', 'R']]) + assert geometry.positions['AR*'][0] < geometry.positions['R'][0] + + @pytest.mark.parametrize('name, mec', SAMPLES, ids=[n for n, _ in SAMPLES]) + def test_every_sample_mechanism_has_a_ready_made_layout(self, name, mec): + assert mec.layout is None + geometry = schemes.geometry(mec) + assert set(geometry.positions) == set(mec.state_names) + + def test_a_state_left_off_the_layout_is_named(self): + with pytest.raises(ValueError, match='AR\\*'): + schemes.geometry(CCO(), layout=[['R', 'AR']]) + + def test_a_name_that_is_not_a_state_is_named(self): + with pytest.raises(ValueError, match='A2R'): + schemes.geometry(CCO(), layout=[['R', 'AR', 'A2R'], [None, 'AR*']]) + + def test_two_states_in_one_place_is_an_error(self): + with pytest.raises(ValueError, match='same place'): + schemes.geometry(CCO(), layout={'R': (0, 0), 'AR': (1, 0), + 'AR*': (1, 0)}) + + def test_an_arrow_drawn_through_another_state_is_an_error(self): + """R -> AR would pass straight over AR*, and read as two steps that + the mechanism does not have.""" + with pytest.raises(ValueError, match='AR\\*'): + schemes.geometry(CCO(), layout=[['R', 'AR*', 'AR']]) + + def test_with_no_layout_anywhere_the_error_says_how_to_give_one(self): + with pytest.raises(ValueError, match='layout='): + schemes.geometry(_two_state()) + + +class TestValues: + + def test_nothing_is_written_by_default(self): + assert all(a.label is None for a in schemes.geometry(CH82()).arrows) + + def test_transition_probabilities(self): + mec = CH82() + geometry = schemes.geometry(mec, values='probability', conc=100e-9) + mec.set_eff('c', 100e-9) + pi = qml.transition_probability(mec.Q) + for arrow in geometry.arrows: + i, j = mec.state_index(arrow.source), mec.state_index(arrow.target) + assert arrow.label == schemes.format_number(pi[i, j]) + + def test_transition_frequencies_at_equilibrium(self): + mec = CH82() + geometry = schemes.geometry(mec, values='frequency', conc=100e-9) + mec.set_eff('c', 100e-9) + f = qml.transition_frequency(mec.Q) + for arrow in geometry.arrows: + i, j = mec.state_index(arrow.source), mec.state_index(arrow.target) + assert arrow.label == schemes.format_number(f[i, j]) + + def test_rate_constants_need_no_concentration(self): + mec = CH82() + geometry = schemes.geometry(mec, values='rates') + for arrow in geometry.arrows: + expected = mec.rate(arrow.source, arrow.target).unit_rate() + assert arrow.label == schemes.format_number(expected) + + def test_probabilities_need_a_concentration(self): + """Without one the Q matrix is at whatever concentration was last + set - 1 M if none was - and the figure would be silently wrong.""" + with pytest.raises(ValueError, match='conc'): + schemes.geometry(CH82(), values='probability') + + def test_a_mechanism_without_ligand_needs_none(self): + mec = samples.RFO() + geometry = schemes.geometry(mec, values='probability') + assert all(a.label is not None for a in geometry.arrows) + + def test_the_mechanisms_own_concentration_is_left_alone(self): + mec = CH82() + mec.set_eff('c', 1e-6) + before = np.array(mec.Q, copy=True) + schemes.geometry(mec, values='frequency', conc=100e-9) + assert np.array(mec.Q) == pytest.approx(before) + + def test_a_function_may_label_some_directions_and_not_others(self): + geometry = schemes.geometry( + CCO(), values=lambda a, b: 'on' if (a, b) == ('R', 'AR') else None) + labels = {(a.source, a.target): a.label for a in geometry.arrows} + assert labels[('R', 'AR')] == 'on' + assert labels[('AR', 'R')] is None + + def test_an_unknown_kind_of_value_is_an_error(self): + with pytest.raises(ValueError, match='probability'): + schemes.geometry(CH82(), values='probabilities') + + +class TestFormatNumber: + """Numbers as the published figures print them.""" + + @pytest.mark.parametrize('value, text', [ + (0.0, '0'), + (1.0, '1'), + (0.418, '0.418'), + (0.0573, '0.0573'), + (6471.7296, '6472'), # f43: thousands print whole + (1850.6, '1851'), + ]) + def test_as_in_the_paper(self, value, text): + assert schemes.format_number(value) == text + + def test_large_rate_constants_use_a_power_of_ten(self): + assert schemes.format_number(1.66282e8) == r'$1.66\times10^{8}$' + + +class TestStateLabels: + + @pytest.mark.parametrize('name, label', [ + ('R', r'$\mathrm{R}$'), + ('A2R', r'$\mathrm{A_{2}R}$'), + ('A2F*', r'$\mathrm{A_{2}F^{*}}$'), + ('AR*a', r'$\mathrm{AR^{*}a}$'), + ]) + def test_digits_are_subscripts_and_stars_superscripts(self, name, label): + assert schemes.state_label(name) == label + + def test_a_name_mathtext_would_mangle_is_left_as_it_is(self): + assert schemes.state_label('A_2 R') == 'A_2 R' + + +class TestColoursAndNumbers: + + def test_open_states_are_red_and_shut_ones_black(self): + colours = schemes.state_colours(CH82()) + style = schemes.Style() + assert colours['AR*'] == colours['A2R*'] == style.open_colour + assert colours['R'] == colours['AR'] == style.shut_colour + + def test_any_state_can_be_given_its_own_colour(self): + colours = schemes.state_colours(GlyR_flip(), + {'AF': 'green', 'A2F': 'green'}) + assert colours['AF'] == colours['A2F'] == 'green' + assert colours['A3F'] == schemes.Style().shut_colour + + def test_numbers_true_numbers_states_in_q_matrix_order(self): + numbers = schemes.state_numbers(CH82(), True) + assert numbers == {'AR*': 1, 'A2R*': 2, 'AR': 3, 'A2R': 4, 'R': 5} + + def test_a_paper_s_own_numbering_can_be_given(self): + numbers = schemes.state_numbers(GlyR_flip(), {'A3F': 4, 'AF': 6}) + assert numbers == {'A3F': 4, 'AF': 6} + + def test_no_numbers_by_default(self): + assert schemes.state_numbers(CH82(), False) == {} + + +class TestImport: + + def test_importing_the_module_does_not_import_matplotlib(self): + """matplotlib is an optional extra: everything but draw() has to work + without it.""" + code = ('import sys; import scalcs.schemes; ' + 'print("matplotlib" in sys.modules)') + out = subprocess.run([sys.executable, '-c', code], capture_output=True, + text=True, check=True) + assert out.stdout.strip() == 'False' + + +class TestDrawing: + + @pytest.fixture(autouse=True) + def _agg(self): + pytest.importorskip('matplotlib') + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + yield + plt.close('all') + + @pytest.mark.parametrize('name, mec', SAMPLES, ids=[n for n, _ in SAMPLES]) + def test_every_sample_mechanism_draws(self, name, mec): + import matplotlib.pyplot as plt + ax = schemes.draw(mec) + barbs = [p for p in ax.patches if isinstance(p, plt.Polygon)] + assert len(barbs) == len(mec.Rates) + texts = {t.get_text() for t in ax.texts} + assert {schemes.state_label(n) for n in mec.state_names} <= texts + + def test_values_numbers_and_colours_reach_the_figure(self): + mec = CH82() + ax = schemes.draw(mec, values='rates', numbers=True, + colours={'R': 'blue'}) + texts = {t.get_text(): t for t in ax.texts} + assert schemes.format_number(mec.rate('R', 'AR').unit_rate()) in texts + assert '5' in texts + assert texts[schemes.state_label('R')].get_color() == 'blue' + + def _state_fontsize(self, ax, name='R'): + return next(t.get_fontsize() for t in ax.texts + if t.get_text() == schemes.state_label(name)) + + def test_a_new_figure_is_drawn_at_the_reference_scale(self): + """The style's font sizes are the ones measured from the published + figures, and a figure draw() makes for itself reproduces them.""" + ax = schemes.draw(CH82()) + assert self._state_fontsize(ax) == pytest.approx( + schemes.Style().state_fontsize, rel=0.02) + + def test_text_shrinks_with_the_axes_it_is_drawn_into(self): + """Lengths are in data units and font sizes in points, so a scheme + squeezed into a small panel would otherwise have its names overrun + the arrows. Text is scaled with the axes instead.""" + import matplotlib.pyplot as plt + # both panels smaller than the reference, so neither is capped + big = schemes.draw(CH82(), ax=plt.figure(figsize=(3, 2.25)).add_subplot()) + small = schemes.draw(CH82(), ax=plt.figure(figsize=(1.5, 1.125)).add_subplot()) + assert self._state_fontsize(big) < schemes.Style().state_fontsize + assert self._state_fontsize(small) / self._state_fontsize(big) == pytest.approx(0.5, rel=0.02) + + def test_text_is_never_enlarged_beyond_the_published_size(self): + """A small scheme in a big panel keeps the published sizes; scaled up, + the numbers of neighbouring steps crowd each other.""" + import matplotlib.pyplot as plt + ax = schemes.draw(CCO(), ax=plt.figure(figsize=(10, 8)).add_subplot()) + assert self._state_fontsize(ax, 'AR') == schemes.Style().state_fontsize + + def test_scaling_can_be_turned_off(self): + import matplotlib.pyplot as plt + ax = schemes.draw(CH82(), ax=plt.figure(figsize=(4, 3)).add_subplot(), + style=schemes.Style(scale_text=False)) + assert self._state_fontsize(ax) == schemes.Style().state_fontsize + + def test_it_draws_into_the_axes_it_is_given(self): + import matplotlib.pyplot as plt + fig, (left, right) = plt.subplots(1, 2) + assert schemes.draw(CCO(), ax=right) is right + assert not left.texts From 50e639fa9a6b2c51a9940ade108759771fff0f08 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 16:27:03 +0000 Subject: [PATCH 2/3] Declare matplotlib as the plot extra, and test with it in CI plot = ["matplotlib>=3.5"] for scalcs.schemes.draw and the one drawing function in sccurves; the rest of the package only returns arrays gui now lists matplotlib, which the GUI imported without declaring CI installs the plot extra on every leg. Without it the scheme-drawing tests importorskip and pass without drawing anything. The wheel check, which runs in a clean environment with no matplotlib, now also imports scalcs.schemes, so the lazy import cannot regress unnoticed. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/ci.yml | 10 +++++++--- pyproject.toml | 6 +++++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index de07e05..d988d9f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -55,9 +55,11 @@ jobs: run: | python -m pip install --upgrade pip if [ "${{ matrix.gui }}" = "true" ]; then - pip install -e '.[dev,gui]' matplotlib + pip install -e '.[dev,gui,plot]' else - pip install -e '.[dev]' + # plot as well, or the scheme-drawing tests importorskip and + # pass without drawing anything + pip install -e '.[dev,plot]' fi pip list @@ -114,6 +116,8 @@ jobs: if missing: raise SystemExit('missing from the installed wheel: ' + ', '.join(missing)) - from scalcs import scsim, scburst, qmatlib, mechanism, scalcsio + # schemes too: it must import without matplotlib, which this + # environment does not have + from scalcs import scsim, scburst, qmatlib, mechanism, scalcsio, schemes print('installed wheel imports and carries its data files') PY diff --git a/pyproject.toml b/pyproject.toml index 152a1c1..5e01ff8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,7 +32,11 @@ classifiers = [ ] [project.optional-dependencies] -gui = ["PyQt5>=5.12"] +# scalcs.schemes.draw and the one plotting function in sccurves; the rest +# of the package only returns arrays and does not need it +plot = ["matplotlib>=3.5"] +# the GUI draws with matplotlib too, which it imported without declaring +gui = ["PyQt5>=5.12", "matplotlib>=3.5"] dev = ["pytest>=7", "pytest-cov"] [project.urls] From fa6d5cecada587d1c8554b67ce97146b9cde0645 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Thu, 12 Nov 2026 16:34:35 +0000 Subject: [PATCH 3/3] Skip the Qt plot-window tests when matplotlib has no Qt binding The guard was importorskip('matplotlib.backends.backend_qtagg'). Without a Qt binding that module exists but raises a plain ImportError, not ModuleNotFoundError, and from pytest 9.1 importorskip no longer counts that as a reason to skip unless passed exc_type=ImportError. The import error then interrupts collection, and pytest exits with code 2. Latent on master: its Linux CI legs install no matplotlib, so the module was skipped one line earlier, by importorskip('matplotlib'). Installing the plot extra there exposed it, failing every Linux leg with pytest 9.2.0. It would equally have hit anyone running the tests with matplotlib but no Qt. Reproduced locally by blocking every Qt binding and making pytest's deprecation warning an error, as 9.1 does: collection error and exit 2 before, 760 passed and this module skipped after. Co-Authored-By: Claude Opus 5.5 --- tests/test_gui_plotwindow.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/test_gui_plotwindow.py b/tests/test_gui_plotwindow.py index 13c14eb..6094ce9 100644 --- a/tests/test_gui_plotwindow.py +++ b/tests/test_gui_plotwindow.py @@ -19,8 +19,12 @@ class already receives `axis` and was simply not passing it on. matplotlib = pytest.importorskip('matplotlib') matplotlib.use('Agg') +# Without a Qt binding the backend module exists but raises a plain ImportError +# rather than ModuleNotFoundError, which importorskip no longer treats as a +# skip from pytest 9.1 unless told to. Linux CI legs have matplotlib but no Qt. pytest.importorskip('matplotlib.backends.backend_qtagg', - reason='no Qt binding available for matplotlib') + reason='no Qt binding available for matplotlib', + exc_type=ImportError) from matplotlib import scale as mscale # noqa: E402