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/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/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] 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_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 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