From 2875cd90301ce877dde87b422477a3c2b5206417 Mon Sep 17 00:00:00 2001 From: Remis Date: Fri, 11 Sep 2026 21:14:33 +0100 Subject: [PATCH] A fit specification, with a command and a notebook over it Step 4 of the fitting layer. HJCFitter and FitResult made a fit expressible in Python; this makes it expressible in a file, so that the same description drives a notebook, a command, a batch runner and anything else -- including, one day, a desktop interface, which is a way of editing a fit specification. hjcfit template -o my-fit.toml hjcfit check my-fit.toml everything a fit does except the search hjcfit fit my-fit.toml -o out.json the search, and a result with provenance fitspec.py is the description and executes nothing: it imports only the standard library, and a test asserts that from a subprocess. runner.py is the half that executes -- records through dcio, mechanism through scalcs, search through HJCFitter. cli.py is the command; examples/fit_template.ipynb and examples/CH82.toml are the same fit as a notebook. TOML rather than the YAML the plan named. A specification is mostly rate constants, and PyYAML implements YAML 1.1, whose resolver wants a signed exponent: yaml.safe_load reads 1e8 as the string '1e8', and so does 1.0e8, and so does 1e+8. Only 1.0e+8 is a float. TOML takes all of them and is in the standard library from 3.11; tomli covers 3.10 and is in the extra. Two defects this surfaced, both silent, both guarded in build_mechanism: * mr naming a rate outside its cycle. SCALCS warns on stderr and carries on, which is worse than refusing: the rate leaves the free-parameter list while the cycle's constraint stays where it was, so it is neither fitted nor computed. On CH82 with mr = "k(-1)", eight free parameters became seven with k(-1) stuck at its guess. * a second rate of a cycle left flagged. A cycle determines one rate and set_mr does not clear the previous one -- update_mr carries the comment "TODO: check for consistency between cycle.mrconstr and rate.mr". CH82's sample already has 2k*(-2) under microscopic reversibility, so mr = "beta1" gave seven free parameters rather than eight with 2k*(-2) stuck at 0.411. Also reported rather than swallowed: rates that end on a limit. A rate outside its limits is reset before the likelihood sees it, as HJCFIT did, and nothing announced it -- so the fitted value read like an estimate when it is a statement that the model forbids where the likelihood wanted to go. SCALCS gives every rate default limits, so this happens with no limits set. The notebook draws the predicted distributions over the observed histograms and checks each curve against a quantity it did not come from, because a wrong curve there looks like a finding about the fit rather than a bug. The open density integrates to 1.00000 with a mean of 5.705 ms against 5.703 observed. The shut density has to be conditioned on [tres, tcrit]: a burst fit never saw a shut time above tcrit, only 31.6% of the predicted density lies below it, and the fitted mechanism's slowest shut component is 23 800 s in a record two minutes long. Conditioned, its mean is 0.1944 ms against 0.1938 observed. CI now runs the whole tests directory rather than two named files. tests/test_fitting.py -- the 36 tests written for the fitting layer in #194 -- was run by no job from the day it was added. A list of test files is one more piece of declared metadata that nothing exercises, which is the fifth time that pattern has cost something here. Verified: the fit through the specification reproduces the one in fitting.rst (2286.9746 at the guess, 2289.1326 in 1313 evaluations); load_records gives groups identical, element for element, to read_idealized_bursts; and a wheel built by scikit-build-core from this pyproject writes the hjcfit console script into entry_points.txt with the marked tomli requirement intact. 164 tests, all passing except the console-script one, which needs the package installed from this pyproject and so passes only in CI. Co-Authored-By: Claude Opus 5 --- .github/workflows/ci.yml | 22 +- documentation/source/api/python.rst | 1 + documentation/source/api/python/fitspec.rst | 234 ++++++ documentation/source/api/python/fitting.rst | 4 + examples/CH82.toml | 31 + examples/fit_template.ipynb | 316 ++++++++ likelihood/python/CMakeLists.txt | 2 +- likelihood/python/cli.py | 239 ++++++ likelihood/python/fitspec.py | 776 ++++++++++++++++++++ likelihood/python/fitting.py | 5 +- likelihood/python/runner.py | 542 ++++++++++++++ pyproject.toml | 10 + tests/test_fitspec.py | 298 ++++++++ tests/test_runner_cli.py | 432 +++++++++++ 14 files changed, 2902 insertions(+), 10 deletions(-) create mode 100644 documentation/source/api/python/fitspec.rst create mode 100644 examples/CH82.toml create mode 100644 examples/fit_template.ipynb create mode 100644 likelihood/python/cli.py create mode 100644 likelihood/python/fitspec.py create mode 100644 likelihood/python/runner.py create mode 100644 tests/test_fitspec.py create mode 100644 tests/test_runner_cli.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9a5ad31..dc14bc7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -222,14 +222,22 @@ jobs: # Run from outside the source tree, so the installed package is what is # imported rather than the checkout. # - # test_simplex_hjc.py rides along here rather than in the build matrix: it is - # pure Python and needs no dcio, but it does need HJCFIT *installed*, and - # this is the only job that installs it. The matrix runs ctest against - # the build tree instead. + # The whole directory, not a list of files. This job is the only one that + # installs the package, so it is the only place the pure-Python tests can + # run at all -- the matrix runs ctest against the build tree -- and the + # list that used to be here had fallen behind: tests/test_fitting.py, the + # 36 tests written for the fitting layer, was never run by any job from + # the day it was added. A named list of test files is one more piece of + # declared metadata that nothing exercises. + # + # examples/ goes too, beside tests/ as it is in the repository, because + # the notebook template and the specification it drives are tested. - name: Tests env: MPLBACKEND: Agg run: | - cp tests/test_dcio_integration.py tests/test_simplex_hjc.py /tmp/ - cd /tmp - pytest test_dcio_integration.py test_simplex_hjc.py -v + pip install nbclient nbformat ipykernel + mkdir -p /tmp/hjcfit-tests + cp -r tests examples /tmp/hjcfit-tests/ + cd /tmp/hjcfit-tests + pytest tests -v diff --git a/documentation/source/api/python.rst b/documentation/source/api/python.rst index e4a5431..c13cfb2 100644 --- a/documentation/source/api/python.rst +++ b/documentation/source/api/python.rst @@ -14,6 +14,7 @@ Python API Reference python/roots.rst python/optimization.rst python/fitting.rst + python/fitspec.rst Extras diff --git a/documentation/source/api/python/fitspec.rst b/documentation/source/api/python/fitspec.rst new file mode 100644 index 0000000..1f105b5 --- /dev/null +++ b/documentation/source/api/python/fitspec.rst @@ -0,0 +1,234 @@ +.. _python_fitspec_api: + +A fit from a file +----------------- + +:ref:`python_fitting_api` is the fitting layer as Python objects. This is the +same fit written down: a **fit specification**, which says which records, which +mechanism and how to search, and runs nothing. + +That separation is what lets the same description drive a notebook, a command, +a batch runner and -- if one is ever built -- a desktop interface, because a +user interface is a way of editing a fit specification. It is also what makes a +fit reproducible by somebody else: a file they can read, diff and keep beside +the result. + +:: + + hjcfit template -o my-fit.toml a specification to edit + hjcfit check my-fit.toml say what it would do; run nothing + hjcfit fit my-fit.toml -o out.json run it, and keep the result + +``examples/fit_template.ipynb`` is the same thing as a notebook, over +``examples/CH82.toml``. + +Everything here needs the ``[fitting]`` extra:: + + pip install 'hjcfit[fitting]' + +The file +"""""""" + +.. code-block:: toml + + title = "CH82 sample record at 100 nM" + + [[data]] + record = "CH82" # a sample record, or a path to an .scn file + conc = 1e-07 # M + tres = 0.0001 # dead time to impose, s + tcrit = 0.004 # critical shut time dividing the record, s + vectors = "chs" # or "equilibrium" + + [mechanism] + sample = "CH82" # a factory in scalcs.samples.samples + nfree = 8 # free parameters expected; asserted, not applied + + [mechanism.rates] # the initial guess, by rate name + beta1 = 15.0 + beta2 = 15000.0 + + [search] + method = "simplex" + log_params = true + maxfev = 20000 + +One ``[[data]]`` section per record. Several means several records contributing +to a single likelihood, so the rate constants are shared and the +concentrations are not. + +``hjcfit template`` writes this file with every option in it, commented. + +Why TOML and not YAML +""""""""""""""""""""" + +A specification is mostly rate constants, and rate constants get written +``1e8``. PyYAML implements YAML 1.1, whose resolver requires an exponent to +carry a sign, so ``yaml.safe_load`` reads ``1e8`` as the **string** ``'1e8'`` +-- and so does ``1.0e8``, and so does ``1e+8``. Only ``1.0e+8`` becomes a +float. A rate constant arriving silently as a string is the worst failure this +file could have. + +TOML has one number syntax, accepts every form above, and is in the standard +library from Python 3.11; ``tomli`` covers 3.10 and is in the ``[fitting]`` +extra. The ``.yaml`` file in ``scalcs/samples`` is a third thing again -- a +pickled Python object graph needing ``unsafe_load`` -- so YAML in this stack +already means something other than a document a person edits. + +Check it before you fit it +"""""""""""""""""""""""""" + +``hjcfit check`` does everything a fit does except the search: it finds the +records, builds the mechanism with every constraint applied, and evaluates the +likelihood once at the initial guess. + +.. code-block:: text + + $ hjcfit check examples/CH82.toml + CH82 sample record at 100 nM + CH82 at 100 nM: 100 us dead time, groups at tcrit 4 ms, chs vectors + mechanism: CH82 + guess: beta1 = 15, beta2 = 15000 + search: simplex, over logarithms, up to 20000 evaluations + + 100 nM: 4312 -> 1100 intervals at 100 us -> 572 groups, 836 openings, CHS vectors + + CH82: 5 states, 2 open, 10 rate constants + 8 free: beta1, beta2, alpha1, alpha2, k(-1), 2k(-2), 2k(+1), k(+2) + 2 fixed or constrained: k*(+2), 2k*(-2) + + log10L at the guess: 2286.9746 + nothing was fitted; run: hjcfit fit examples/CH82.toml + +A misspelled rate name, a critical time below the dead time, or a +free-parameter count nobody expected all surface here in a second, rather than +twenty minutes into a search or in a set of estimates that look plausible. + +Two fields rather than one signed number +"""""""""""""""""""""""""""""""""""""""" + +Elsewhere in the stack a negative critical time is a flag selecting +equilibrium vectors while its magnitude is still the time that divides the +record -- see :py:func:`HJCFIT.read_idealized_bursts`. Here they are +``tcrit`` and ``vectors``, because the critical time does two jobs and only +one of them is usually being changed: + +* **``tcrit`` with ``vectors = "chs"``** -- groups cut at a critical time + chosen to separate the activations of one channel, started and ended with + the CHS vectors of Colquhoun, Hawkes & Srodzinski (1996). The usual case. +* **``tcrit`` with ``vectors = "equilibrium"``** -- still groups, but + equilibrium vectors: clusters at a concentration high enough to desensitise. +* **no ``tcrit``** -- the whole record as a single group, which assumes one + channel throughout. ``vectors`` must then be ``"equilibrium"``, because CHS + vectors are defined between groups. + +Rates are named, never numbered +""""""""""""""""""""""""""""""" + +Every reference to a rate constant in a specification is by name. SCALCS +addresses rates by index into ``mec.Rates`` -- ``set_mr(True, 5, 0)`` -- and an +index is what users get wrong silently, because rate 5 of a mechanism is +whatever the sample happened to list fifth. +:py:func:`~HJCFIT.likelihood.runner.build_mechanism` is the one place a name +becomes an index, and it lists the names that exist when one does not. + +It also refuses two things SCALCS allows: + +* **``mr`` naming a rate outside the chosen cycle.** SCALCS warns on stderr + and carries on, which is worse than nothing: the rate leaves the + free-parameter list while the cycle's constraint stays where it was, so the + rate is neither fitted nor computed -- it is frozen at its initial guess and + the fit has one parameter fewer than the file says. +* **A second rate of a cycle left flagged.** A cycle determines exactly one + rate, and ``set_mr`` does not clear the previous one; its own ``update_mr`` + carries the comment ``TODO: check for consistency between cycle.mrconstr and + rate.mr``. On CH82, whose sample already has ``2k*(-2)`` under microscopic + reversibility, asking for ``mr = "beta1"`` gave seven free parameters rather + than eight with ``2k*(-2)`` stuck at whatever the sample carried. + ``build_mechanism`` releases the stale one. + +Rates that end on a limit +""""""""""""""""""""""""" + +A rate driven outside its limits is reset before the likelihood sees it, which +is what HJCFIT did, and nothing announces it. Such a rate is not an estimate: +it is a statement that the likelihood wanted to go somewhere the model forbids. +:py:func:`~HJCFIT.likelihood.runner.against_limits` finds them and they are +printed apart from the fitted values. + +.. code-block:: text + + log10(L) = 2288.7861 1154 evaluations in 1.3 s + beta1 5 + ... + + Rates that ended on a limit -- these are not estimates: + beta1 5 at its upper limit of 5 + +Note that SCALCS gives every rate default limits -- 10\ :sup:`-15` to +10\ :sup:`9` for a concentration-dependent rate and 10\ :sup:`-15` to +10\ :sup:`6` for the rest -- so this can happen in a specification that sets +no limits at all. + +The specification +""""""""""""""""" + +.. currentmodule:: HJCFIT.likelihood.fitspec + +.. autoclass:: FitSpec + :members: from_toml, from_dict, to_toml, write_toml, as_dict, validate + +.. autoclass:: DataSpec + :members: from_dict, validate + +.. autoclass:: MechanismSpec + :members: from_dict, validate + +.. autoclass:: SearchSpec + :members: from_dict, validate + +.. autoexception:: SpecError + +.. autofunction:: load_toml_bytes + +.. autodata:: TEMPLATE + :annotation: + +Carrying it out +""""""""""""""" + +.. currentmodule:: HJCFIT.likelihood.runner + +.. autofunction:: run + +.. autoclass:: Outcome + +.. autofunction:: load_records + +.. autofunction:: build_mechanism + +.. autofunction:: against_limits + +.. autofunction:: provenance + +.. autofunction:: result_as_dict + +.. autofunction:: write_result + +.. autoexception:: RunnerError + +The command +""""""""""" + +.. currentmodule:: HJCFIT.likelihood.cli + +.. autofunction:: main + +Also reachable as ``python -m HJCFIT.likelihood.cli`` when the console script +is not on the path, which is the usual state of affairs inside a conda +environment on Windows. + +Exit status is 0 on success, 1 on anything the command can explain -- one line +beginning with ``hjcfit:`` -- and 2 from ``fit`` when the search did not +converge. The estimates are still printed in that case; it is the status that +says so, because that is what a script reads. diff --git a/documentation/source/api/python/fitting.rst b/documentation/source/api/python/fitting.rst index d8baebf..5adbfd5 100644 --- a/documentation/source/api/python/fitting.rst +++ b/documentation/source/api/python/fitting.rst @@ -15,6 +15,10 @@ constants and constraints live in SCALCS. Install it with the extra:: pip install hjcfit[fitting] +If you would rather write a file than write Python, :ref:`python_fitspec_api` +is this module driven from a specification -- by the ``hjcfit`` command, or by +the notebook template in ``examples/``. + This module is imported on demand rather than exposed from HJCFIT.likelihood itself:: diff --git a/examples/CH82.toml b/examples/CH82.toml new file mode 100644 index 0000000..e61a270 --- /dev/null +++ b/examples/CH82.toml @@ -0,0 +1,31 @@ +# The CH82 sample record that ships with HJCFIT, fitted with the CH82 +# mechanism of Colquhoun & Hawkes (1982). Runs as it stands: +# +# hjcfit check CH82.toml +# hjcfit fit CH82.toml -o result.json +# +# `hjcfit template` writes the same thing with every option in it, commented. + +title = "CH82 sample record at 100 nM" + +[[data]] +record = "CH82" # a sample record; otherwise a path to an .scn file +conc = 1e-07 # M +tres = 0.0001 # dead time to impose, s +tcrit = 0.004 # critical shut time dividing the record into bursts, s +vectors = "chs" # CHS vectors, which are the right ones between bursts + +[mechanism] +sample = "CH82" # a factory in scalcs.samples.samples +nfree = 8 # free parameters expected; asserted before the search + +# The initial guess. Rates left out keep the sample's own value, and +# `hjcfit check` prints every name this mechanism has. +[mechanism.rates] +beta1 = 15.0 +beta2 = 15000.0 + +[search] +method = "simplex" # HJCFIT's own simplex, over the logarithms +log_params = true +maxfev = 20000 diff --git a/examples/fit_template.ipynb b/examples/fit_template.ipynb new file mode 100644 index 0000000..503fd70 --- /dev/null +++ b/examples/fit_template.ipynb @@ -0,0 +1,316 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e2610620", + "metadata": {}, + "source": [ + "# Fitting a mechanism from a specification\n", + "\n", + "This is a template. The only thing you have to change is the file name in the\n", + "next cell, and then the specification it points at.\n", + "\n", + "A **fit specification** is a TOML file that says which records, which mechanism\n", + "and how to search. It describes a fit and runs nothing, so the same file drives\n", + "this notebook, the `hjcfit` command and anything else that ever wants to run\n", + "the fit:\n", + "\n", + " hjcfit check my-fit.toml # everything except the search\n", + " hjcfit fit my-fit.toml # the search too\n", + "\n", + "`examples/CH82.toml` fits the CH82 sample record that ships with HJCFIT, so\n", + "this notebook runs before you have changed anything. Write a new one with\n", + "`hjcfit template -o my-fit.toml`, which is the same file with every option in\n", + "it, commented.\n", + "\n", + "Needs the `[fitting]` extra: `pip install 'hjcfit[fitting]'`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a89c12fb", + "metadata": {}, + "outputs": [], + "source": [ + "SPEC = \"CH82.toml\" # <- your specification\n", + "\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "\n", + "from scipy.integrate import trapezoid\n", + "\n", + "from HJCFIT.likelihood import IdealG, QMatrix\n", + "from HJCFIT.likelihood._methods import (dwell_time_histogram, ideal_pdf,\n", + " ideal_pdf_scale_factor,\n", + " missed_events_pdf)\n", + "from HJCFIT.likelihood.fitspec import FitSpec\n", + "from HJCFIT.likelihood.runner import (build_mechanism, load_records, run,\n", + " write_result)\n", + "\n", + "spec = FitSpec.from_toml(SPEC)\n", + "print(spec)" + ] + }, + { + "cell_type": "markdown", + "id": "c6742794", + "metadata": {}, + "source": [ + "## Check it before fitting it\n", + "\n", + "This is what `hjcfit check` does, and it is worth the habit: it finds the\n", + "records, builds the mechanism with every constraint applied, and evaluates the\n", + "likelihood once at the initial guess. A misspelled rate name or a critical time\n", + "below the dead time surfaces here, in a second, rather than twenty minutes into\n", + "a search.\n", + "\n", + "The free-parameter list is the other reason to look. Constraints are easy to\n", + "get wrong in a direction that still runs, and a count you did not expect is the\n", + "symptom. Put the number you expect in the specification as `nfree` and it\n", + "becomes a refusal to start rather than something to notice later." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4e34716", + "metadata": {}, + "outputs": [], + "source": [ + "records = load_records(spec)\n", + "mec = build_mechanism(spec)\n", + "\n", + "for record in records:\n", + " print(record)\n", + "\n", + "free = mec.get_free_parameter_names()\n", + "print()\n", + "print(f\"{mec.k} states, {mec.kA} open, {len(mec.Rates)} rate constants\")\n", + "print(f\"{len(free)} free: {', '.join(free)}\")\n", + "held = [r.name for r in mec.Rates if not r.is_free]\n", + "if held:\n", + " print(f\"{len(held)} fixed or constrained: {', '.join(held)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "5192a401", + "metadata": {}, + "source": [ + "## Fit\n", + "\n", + "`run` reads the records, builds the mechanism, searches, and reports. The\n", + "default search is `simplex_hjc`, HJCFIT's own -- the search that produced every\n", + "published result -- over the logarithms of the rate constants.\n", + "\n", + "Watch two things in the output. **Evaluations that could not be computed** are\n", + "reported rather than swallowed, because a mechanism the likelihood struggles\n", + "with should be visible. And **rates that ended on a limit** are printed\n", + "separately from the fitted values: a rate outside its limits is reset before\n", + "the likelihood sees it, so such a rate is not an estimate but a statement that\n", + "the likelihood wanted to go somewhere the model forbids. SCALCS gives every\n", + "rate default limits, so this can happen in a specification that sets none." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "511ad73d", + "metadata": {}, + "outputs": [], + "source": [ + "outcome = run(spec, verbose=True)\n", + "print()\n", + "print(outcome)" + ] + }, + { + "cell_type": "markdown", + "id": "109b457c", + "metadata": {}, + "source": [ + "## Does the fit describe the record?\n", + "\n", + "The estimates above say what the maximum is. They do not say whether the\n", + "mechanism describes what was recorded, and that is the question an\n", + "experimenter actually faces. The display for it is the one Sigworth and Sine\n", + "introduced: dwell times binned uniformly in log time, on a square-root\n", + "ordinate, with the density the fit predicts drawn over them.\n", + "\n", + "The curves come from the **fitted** mechanism, so they are what this fit\n", + "predicts and not what any truth predicts.\n", + "\n", + "**The shut-time curve has to be conditioned, and the open-time curve does\n", + "not.** This record was divided into bursts at a critical shut time of 4 ms, so\n", + "the likelihood never saw a shut time longer than that: every long shut interval\n", + "lies *between* bursts and was thrown away with them. The openings were not\n", + "touched by that, so every opening in the record is in the fit. Drawing the\n", + "apparent shut-time density over a histogram of within-burst shut times\n", + "therefore means comparing a density over all shut times with a sample of the\n", + "short ones, and the curve would sit far below the histogram for no reason\n", + "having anything to do with the fit. Dividing it by its own mass in\n", + "[t_res, t_crit] makes it the density of the shut times that were actually\n", + "fitted.\n", + "\n", + "How far this matters here: the fitted mechanism's slowest predicted shut time\n", + "constant is 23 800 s -- six and a half hours -- in a record about two minutes\n", + "long. A burst fit does not identify that component and was never asked to, and\n", + "only 31.6% of the predicted shut-time density lies below t_crit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9441e280", + "metadata": {}, + "outputs": [], + "source": [ + "record = outcome.records[0]\n", + "mec = outcome.mec\n", + "mec.set_eff('c', record.conc) # the Q matrix at this record's concentration\n", + "qmatrix = QMatrix(mec.Q, mec.kA)\n", + "tres, tcrit = record.tres, record.tcrit\n", + "\n", + "groups = [np.asarray(g, dtype=float) for g in record.groups]\n", + "opens = np.concatenate([g[0::2] for g in groups])\n", + "shuts = np.concatenate([g[1::2] for g in groups if len(g) > 1])\n", + "\n", + "apparent_open = missed_events_pdf(qmatrix, tres, shut=False)\n", + "apparent_shut = missed_events_pdf(qmatrix, tres, shut=True)\n", + "ideal_open = ideal_pdf(qmatrix, shut=False)\n", + "idealG = IdealG(qmatrix)\n", + "\n", + "# The mass of the apparent shut density below tcrit: the fraction of shut times\n", + "# a burst fit can see. Integrated on the same log grid the curve is drawn on,\n", + "# because the density spans eight decades.\n", + "grid = np.logspace(np.log10(tres), np.log10(tcrit), 4000)\n", + "shut_mass = float(trapezoid(apparent_shut(grid), grid))\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", + "dwell_time_histogram(\n", + " opens, tres, ax=axes[0], pdf=apparent_open, ideal=ideal_open,\n", + " ideal_scale=ideal_pdf_scale_factor(tres, qmatrix.aa, idealG.initial_vectors),\n", + " xlabel='Apparent open time (s)')\n", + "axes[0].set_title(f'{len(opens)} openings')\n", + "\n", + "dwell_time_histogram(\n", + " shuts, tres, ax=axes[1], tcrit=tcrit,\n", + " pdf=lambda t: apparent_shut(t) / shut_mass,\n", + " xlabel='Apparent shut time (s)')\n", + "axes[1].set_title(f'{len(shuts)} shut times within bursts '\n", + " f'({100 * shut_mass:.1f}% of the density)')\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "e7b3ddde", + "metadata": {}, + "source": [ + "## Check the curves before believing them\n", + "\n", + "A predicted curve over an observed histogram is the one plot where a wrong\n", + "curve looks like a finding about the fit rather than a bug, so it is worth one\n", + "cell to check each one against a quantity it did not come from. A density must\n", + "integrate to one over the range it is a density of, and its mean must match the\n", + "mean of the sample drawn under it.\n", + "\n", + "If either line below disagrees, the curve above is wrong and the fit is not\n", + "what is in question." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5d7082f", + "metadata": {}, + "outputs": [], + "source": [ + "wide = np.logspace(np.log10(tres), 8, 20000)\n", + "open_area = float(trapezoid(apparent_open(wide), wide))\n", + "open_mean = float(trapezoid(wide * apparent_open(wide), wide))\n", + "shut_mean = float(trapezoid(grid * apparent_shut(grid), grid)) / shut_mass\n", + "\n", + "print(f\"open: pdf integrates to {open_area:.5f}, mean {1e3 * open_mean:.4g} ms \"\n", + " f\"against {1e3 * opens.mean():.4g} ms observed\")\n", + "print(f\"shut: {100 * shut_mass:.2f}% of the density below tcrit, conditional \"\n", + " f\"mean {1e3 * shut_mean:.4g} ms against {1e3 * shuts.mean():.4g} ms observed\")" + ] + }, + { + "cell_type": "markdown", + "id": "f6c1995d", + "metadata": {}, + "source": [ + "## Keep the result\n", + "\n", + "`write_result` writes the estimates, the specification they came from, and the\n", + "versions of everything involved, as JSON. The specification goes into the file\n", + "on purpose: a result and the description of what produced it belong in one\n", + "place, and a number in a paper that can be traced to a file is worth more than\n", + "one that has to be argued about.\n", + "\n", + "The records are not in it. They are thousands of floats, and they are in the\n", + ".scn file the specification names." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01d4d45b", + "metadata": {}, + "outputs": [], + "source": [ + "write_result(\"result.json\", outcome)\n", + "\n", + "import json\n", + "print(json.dumps(outcome.provenance, indent=2))" + ] + }, + { + "cell_type": "markdown", + "id": "ad54dd40", + "metadata": {}, + "source": [ + "## Your own data\n", + "\n", + "Change the specification, not this notebook. `hjcfit template` writes one with\n", + "every option in it; the three that matter first are:\n", + "\n", + "* **`record`** -- the path to your `.scn` file, in place of a sample name.\n", + "* **`tres` and `tcrit`** -- the dead time to impose and the critical shut time\n", + " that divides the record. `tcrit` is what decides whether an interval starts a\n", + " new group, and getting it wrong changes the answer, not just the display.\n", + " Leave `tcrit` out to fit the whole record as one group, which assumes one\n", + " channel throughout; then `vectors` must be `\"equilibrium\"`.\n", + "* **`[mechanism]`** -- `sample` names a mechanism from `scalcs.samples.samples`;\n", + " `mec_file` reads a DCprogs `.mec` file instead. Give the initial guess under\n", + " `[mechanism.rates]` by rate name -- `hjcfit check` prints every name the\n", + " mechanism has, so a typo is a list rather than a puzzle.\n", + "\n", + "**Several concentrations at once** is more than one `[[data]]` section. Each\n", + "record is evaluated at its own concentration and the logarithms are added, so\n", + "the rate constants are shared and the concentrations are not. The displays\n", + "above then want one column per record.\n", + "\n", + "**Which search.** `method = \"simplex\"` is the default and is what to quote.\n", + "Reach for `\"scipy\"` when the starting point is poor: a regular simplex needs\n", + "every one of its vertices to be evaluable, and on CH82's eight parameters only\n", + "6 of 200 random starting points give a finite likelihood at all." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/likelihood/python/CMakeLists.txt b/likelihood/python/CMakeLists.txt index bc891eb..7486a61 100644 --- a/likelihood/python/CMakeLists.txt +++ b/likelihood/python/CMakeLists.txt @@ -70,7 +70,7 @@ if(pythonBindings) # Named rather than globbed, because CMake will not re-run on a new file # appearing and a glob would make the wheel depend on when cmake last ran. set(PYTHON_MODULES __init__.py random.py _methods.py optimization.py - simplex_hjc.py fitting.py) + simplex_hjc.py fitting.py fitspec.py runner.py cli.py) install(FILES ${PYTHON_MODULES} ${CMAKE_CURRENT_BINARY_DIR}/likelihood.py diff --git a/likelihood/python/cli.py b/likelihood/python/cli.py new file mode 100644 index 0000000..13865c6 --- /dev/null +++ b/likelihood/python/cli.py @@ -0,0 +1,239 @@ +######################## +# HJCFIT computes missed-events likelihood as described in +# Hawkes, Jalali and Colquhoun (1990, 1992) +# +# Copyright (C) 2013 University College London +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +######################### + +""" The ``hjcfit`` command: a fit from a file, with no Python written. + + :: + + hjcfit template -o ch82.toml a specification to edit + hjcfit check ch82.toml say what it would do; run nothing + hjcfit fit ch82.toml -o out.json run it, and keep the result + + Also reachable as ``python -m HJCFIT.likelihood.cli`` when the console + script is not on the path, which is the usual state of affairs inside a + conda environment on Windows. + + ``check`` is the one worth the habit. It reads the specification, finds the + records, builds the mechanism with every constraint applied, and evaluates + the likelihood once at the initial guess -- everything a fit does except + the search. A misspelled rate name, a critical time below the dead time, a + mechanism with a free-parameter count nobody expected: all of them surface + in under a second, rather than twenty minutes into a search or, worse, in + a set of estimates that look plausible. + + Every error this prints is one line beginning with ``hjcfit:`` and the exit + status is 1, because this is a command that belongs in scripts. +""" +__docformat__ = "restructuredtext en" +__all__ = ['main'] + +import argparse +import sys + +from .fitspec import TEMPLATE, FitSpec, SpecError + + +def _add_spec_argument(parser): + parser.add_argument('spec', metavar='SPEC.toml', + help='the fit specification') + + +def _parser(): + parser = argparse.ArgumentParser( + prog='hjcfit', + description='Maximum likelihood fitting of ion-channel mechanisms, ' + 'driven by a specification file.', + epilog='Start with: hjcfit template -o my-fit.toml') + parser.add_argument('--version', action='store_true', + help='print the versions of everything involved') + sub = parser.add_subparsers(dest='command', metavar='COMMAND') + + template = sub.add_parser( + 'template', help='write a specification to start from') + template.add_argument('-o', '--output', metavar='FILE', default=None, + help='write here instead of to standard output') + template.add_argument('-f', '--force', action='store_true', + help='overwrite FILE if it exists') + + check = sub.add_parser( + 'check', help='say what a specification would do, and do nothing') + _add_spec_argument(check) + + fit = sub.add_parser('fit', help='run a specification') + _add_spec_argument(fit) + fit.add_argument('-o', '--output', metavar='FILE', default=None, + help='write the result, with its provenance, as JSON') + fit.add_argument('-q', '--quiet', action='store_true', + help='print the result and nothing before it') + + show = sub.add_parser( + 'show', help='read a specification and write it back, normalised') + _add_spec_argument(show) + return parser + + +def _fail(message): + """ One line on stderr, after whatever has already been printed. + + The flush is not decoration. Python buffers stdout when it is not a + terminal, so without it a message written to the unbuffered stderr + arrives before the output it is about, and a command that printed + three useful lines and then failed reads as though it failed first. + """ + sys.stdout.flush() + sys.stderr.write('hjcfit: {0}\n'.format(message)) + + +def _load(path): + """ Read a specification, or fail with one line. """ + try: + return FitSpec.from_toml(path) + except FileNotFoundError: + raise SpecError('no such file: {0}'.format(path)) + except OSError as error: + raise SpecError('cannot read {0}: {1}'.format(path, error)) + + +def _versions(): + from .runner import provenance + prov = provenance() + lines = ['python {0} on {1}'.format(prov['python'], prov['platform'])] + for package, version in prov['versions'].items(): + lines.append('{0:12s} {1}'.format( + package, version if version else 'not installed')) + return "\n".join(lines) + + +def _template(args): + if args.output is None: + sys.stdout.write(TEMPLATE) + return 0 + import os + if os.path.exists(args.output) and not args.force: + raise SpecError('{0} exists; pass --force to overwrite it' + .format(args.output)) + with open(args.output, 'w', encoding='utf-8', newline='\n') as handle: + handle.write(TEMPLATE) + print('wrote {0}'.format(args.output)) + print('edit it, then: hjcfit check {0}'.format(args.output)) + return 0 + + +def _check(args): + """ Everything a fit does except the search. """ + from .runner import build_mechanism, load_records + from .fitting import HJCFitter + + spec = _load(args.spec) + print(spec) + print('') + + records = load_records(spec) + for record in records: + print(record) + + mec = build_mechanism(spec) + free = list(mec.get_free_parameter_names()) + print('') + print('{0}: {1} states, {2} open, {3} rate constants'.format( + getattr(mec, 'mtitle', '') or 'mechanism', mec.k, mec.kA, + len(mec.Rates))) + print('{0} free: {1}'.format(len(free), ", ".join(free))) + held = [r.name for r in mec.Rates if not r.is_free] + if held: + print('{0} fixed or constrained: {1}'.format( + len(held), ", ".join(held))) + + # Evaluating the likelihood once is the whole point of check: it is the + # only thing here that proves the records, the mechanism and the dead time + # are mutually possible. A guess at which the likelihood cannot be + # computed is a guess no search can start from. + fitter = HJCFitter(mec, records, log_params=spec.search.log_params) + print('') + print('log10L at the guess: {0:.4f}'.format(fitter.log10_likelihood())) + print('nothing was fitted; run: hjcfit fit {0}'.format(args.spec)) + return 0 + + +def _fit(args): + from .runner import run, write_result + + spec = _load(args.spec) + if not args.quiet: + print(spec) + print('') + outcome = run(spec, verbose=not args.quiet) + if not args.quiet: + print('') + print(outcome) + if args.output: + write_result(args.output, outcome) + print('') + print('wrote {0}'.format(args.output)) + # A fit whose optimiser did not converge has still produced a number, and + # the number is printed; the exit status is what a script reads, so it says + # so. Rates sitting on a limit are not a failure of the search -- they are + # a statement about the model -- and are reported without changing it. + return 0 if outcome.result.success else 2 + + +def _show(args): + sys.stdout.write(_load(args.spec).to_toml()) + return 0 + + +def main(argv=None): + """ The ``hjcfit`` command. + + :param argv: Arguments, without the program name. None takes + ``sys.argv``. + :returns: 0 on success, 1 on anything this module can explain, and 2 + from ``fit`` when the search did not converge. + """ + parser = _parser() + args = parser.parse_args(argv) + + if args.version: + print(_versions()) + return 0 + if args.command is None: + parser.print_help() + return 0 + + handlers = {'template': _template, 'check': _check, 'fit': _fit, + 'show': _show} + try: + return handlers[args.command](args) + except SpecError as error: + _fail(error) + return 1 + except Exception as error: + # RunnerError and anything dcio or scalcs raises. One line, named, and + # no traceback: a traceback out of a command is noise to the person who + # mistyped a rate name. --debug is deliberately not offered, because + # `python -c "from HJCFIT.likelihood.runner import run"` is the same + # thing with a real traceback. + from .runner import RunnerError + if not isinstance(error, (RunnerError, OSError, ValueError, + ArithmeticError, KeyError, TypeError)): + raise + _fail('{0}: {1}'.format(type(error).__name__, error)) + return 1 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/likelihood/python/fitspec.py b/likelihood/python/fitspec.py new file mode 100644 index 0000000..5569cae --- /dev/null +++ b/likelihood/python/fitspec.py @@ -0,0 +1,776 @@ +######################## +# HJCFIT computes missed-events likelihood as described in +# Hawkes, Jalali and Colquhoun (1990, 1992) +# +# Copyright (C) 2013 University College London +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +######################### + +""" What a fit is, written down. + + A fit needs three things said about it: which records, which mechanism and + how to search. This module is those three things as data, and **nothing + else**. It reads no records, builds no mechanism and runs no search, and it + imports nothing but the standard library. + + That separation is the point. The same description drives + :py:mod:`~HJCFIT.likelihood.runner`, the ``hjcfit`` command, the notebook + template in ``examples/`` and -- if one is ever built -- a desktop + interface, because a user interface is a way of editing a fit + specification. It is also what makes a fit reproducible by somebody else: + a file they can read, diff and keep beside the result. + + A specification is a file:: + + title = "CH82 sample record at 100 nM" + + [[data]] + record = "CH82" + conc = 100e-9 + tres = 100e-6 + tcrit = 4e-3 + + [mechanism] + sample = "CH82" + + [search] + method = "simplex" + + and :py:meth:`FitSpec.from_toml` reads it back. + + **Why TOML rather than YAML.** A specification is mostly rate constants, + and rate constants get written ``1e8``. PyYAML implements YAML 1.1, whose + resolver requires an exponent to carry a sign: ``yaml.safe_load`` reads + ``1e8`` as the *string* ``'1e8'``, and so does ``1.0e8``, and so does + ``1e+8``. Only ``1.0e+8`` becomes a float. A rate constant silently + arriving as a string is the worst failure this file could have. TOML has + one number syntax, accepts every form above, and is in the standard + library from Python 3.11 -- ``tomli`` covers 3.10 and is in the + ``[fitting]`` extra. The ``.yaml`` file in ``scalcs/samples`` is a third + thing again, a pickled Python object graph needing ``unsafe_load``, so + YAML in this stack already means something other than a document a person + edits. + + **Rates are named, never numbered.** Every reference to a rate constant + here is by name. SCALCS addresses rates by index into ``mec.Rates`` -- + ``set_mr(True, 5, 0)`` -- and an index is what users get wrong silently, + because rate 5 of a mechanism is whatever the sample happened to list + fifth. :py:mod:`~HJCFIT.likelihood.runner` resolves names to indices in + one place, and says which names exist when one does not. +""" +__docformat__ = "restructuredtext en" +__all__ = ['DataSpec', 'MechanismSpec', 'SearchSpec', 'FitSpec', + 'SpecError', 'TEMPLATE', 'load_toml_bytes'] + +from dataclasses import dataclass, field + + +class SpecError(ValueError): + """ A specification that cannot mean anything. + + Raised by the ``from_dict`` constructors and by ``validate``. Every + message names the field it is about. + """ + + +#: Accepted values of :py:attr:`DataSpec.vectors`. +VECTORS = ('chs', 'equilibrium') + +#: Accepted values of :py:attr:`SearchSpec.method`. +METHODS = ('simplex', 'scipy') + + +def _number(value, where): + if isinstance(value, bool) or not isinstance(value, (int, float)): + raise SpecError("{0}: expected a number, got {1!r}".format(where, value)) + return float(value) + + +def _named_numbers(mapping, where): + """ A ``{rate name: value}`` table, checked and copied. """ + if mapping is None: + return {} + if not isinstance(mapping, dict): + raise SpecError("{0}: expected a table of 'rate name' = value, got " + "{1!r}".format(where, mapping)) + return {str(k): _number(v, "{0}.{1}".format(where, k)) + for k, v in mapping.items()} + + +def _unexpected(given, known, where): + """ Refuse a key nobody reads. + + A misspelled key in a specification is worse than a missing one: the + fit runs, silently, without the setting the file says it has. Every + table here is closed for that reason. + """ + extra = [k for k in given if k not in known] + if extra: + raise SpecError("{0}: unknown key{1} {2}; known keys are {3}".format( + where, "" if len(extra) == 1 else "s", + ", ".join(repr(k) for k in sorted(extra)), + ", ".join(sorted(known)))) + + +@dataclass(frozen=True) +class DataSpec: + """ One record, and how it is presented to the likelihood. + + :param record: + An ``*.scn`` file, or the name of a sample record shipped with + HJCFIT (``"CH82"``, ``"CO"``, ``"CCO"``). + :param conc: Agonist concentration [M]. + :param tres: Dead time to impose [s]. + :param tcrit: + Critical shut time dividing the record into groups [s]. Omitted, the + whole record is fitted as a single group, which assumes one channel + throughout. + :param vectors: + ``"chs"`` for the CHS vectors of Colquhoun, Hawkes & Srodzinski + (1996), ``"equilibrium"`` for the equilibrium vectors of Colquhoun & + Hawkes (1982). Groups cut at a critical time chosen to separate + activations of one channel want CHS; clusters at a concentration + high enough to desensitise want equilibrium, and so does a whole + record. Defaults to CHS when *tcrit* is given and to equilibrium + when it is not. + + *tcrit* and *vectors* are two fields rather than one signed number. + Elsewhere in the stack a negative critical time is a flag selecting + equilibrium vectors while its magnitude is still the time that divides + the record -- :py:func:`HJCFIT.read_idealized_bursts` documents that, + and :py:class:`~HJCFIT.likelihood.fitting.Record` takes the magnitude + or None. Two fields say the same thing without a sign to remember, and + they make the combination that means nothing -- one group, but vectors + defined *between* groups -- something :py:meth:`validate` can refuse. + """ + + record: str + conc: float + tres: float + tcrit: float = None + vectors: str = 'chs' + + KEYS = ('record', 'conc', 'tres', 'tcrit', 'vectors') + + @classmethod + def from_dict(cls, d, where='data'): + """ Build one from a parsed table, checking as it goes. """ + _unexpected(d, cls.KEYS, where) + for required in ('record', 'conc', 'tres'): + if required not in d: + raise SpecError("{0}: {1} is required".format(where, required)) + tcrit = d.get('tcrit') + return cls( + record=str(d['record']), + conc=_number(d['conc'], where + '.conc'), + tres=_number(d['tres'], where + '.tres'), + tcrit=None if tcrit is None else _number(tcrit, where + '.tcrit'), + vectors=str(d.get('vectors', + 'chs' if tcrit is not None else 'equilibrium')), + ).validate(where) + + def validate(self, where='data'): + """ Raise :py:exc:`SpecError` unless this record can be fitted. """ + if self.conc < 0.0: + raise SpecError("{0}.conc: {1} is negative".format(where, self.conc)) + if self.tres <= 0.0: + raise SpecError("{0}.tres: must be positive, not {1}" + .format(where, self.tres)) + if self.vectors not in VECTORS: + raise SpecError("{0}.vectors: must be one of {1}, not {2!r}".format( + where, " or ".join(VECTORS), self.vectors)) + if self.tcrit is None: + if self.vectors != 'equilibrium': + raise SpecError( + "{0}: with no tcrit the whole record is one group, and " + "CHS vectors are defined between groups; either give " + "tcrit or set vectors to 'equilibrium'".format(where)) + else: + if self.tcrit <= 0.0: + raise SpecError( + "{0}.tcrit: must be positive. The sign is not a flag " + "here -- vectors is.".format(where)) + if self.tcrit <= self.tres: + raise SpecError( + "{0}: tcrit ({1:g} s) must exceed tres ({2:g} s); a " + "critical shut time below the dead time divides nothing" + .format(where, self.tcrit, self.tres)) + return self + + def as_dict(self): + d = {'record': self.record, 'conc': self.conc, 'tres': self.tres} + if self.tcrit is not None: + d['tcrit'] = self.tcrit + d['vectors'] = self.vectors + return d + + def __str__(self): + conc = ("{0:g} nM".format(self.conc * 1e9) if self.conc < 1e-6 + else "{0:g} uM".format(self.conc * 1e6)) + groups = ("whole record as one group" if self.tcrit is None else + "groups at tcrit {0:g} ms".format(self.tcrit * 1e3)) + return "{0} at {1}: {2:g} us dead time, {3}, {4} vectors".format( + self.record, conc, self.tres * 1e6, groups, self.vectors) + + +@dataclass(frozen=True) +class MechanismSpec: + """ Which mechanism, with which constraints and which initial guess. + + :param sample: + Name of a factory in :py:mod:`scalcs.samples.samples`, for example + ``"CH82"`` or ``"load_AChR_diamond_independent_binding"``. Exactly + one of this and *mec_file*. + :param mec_file: + A DCprogs ``*.mec`` file, read with + :py:func:`scalcs.scalcsio.mec_load`. + :param mec_number: + Which mechanism in that file, by the sequence number + ``scalcsio.mec_get_list`` reports. Needed only when the file holds + more than one. + :param rates: + Initial guess, by rate name. Names left out keep the value the + sample or the file carries, so a guess can be one line. + :param fixed: + Rate constants held at a stated value throughout the fit. + :param limits: + ``'name' = [lower, upper]``. A rate driven outside its limits is + **reset** to the limit before the likelihood sees it, which is what + HJCFIT did; see + :py:class:`~HJCFIT.likelihood.fitting.HJCFitter`. Nothing announces + that it happened, so a rate that ends against a limit has to be + noticed by comparing the result with what is set here -- + :py:mod:`~HJCFIT.likelihood.runner` does that comparison and says so. + :param mr: + Name of the rate constant computed from microscopic reversibility + rather than fitted. + :param mr_cycle: + Which cycle, when the mechanism has more than one. Cycles have no + names in SCALCS, so this one is an index. + :param ec50: + ``{rate = 'name', value = EC50 in M}``. The option of Colquhoun, + Hatton & Hawkes (2003) p. 702: an independently measured EC50 is + supplied and one rate constant is computed from it and all the + others at each iteration, which removes a free parameter. An EC50 + given wrongly biases the estimates rather than failing -- that + paper's Figures 9 and 10 are what a factor of two either way does -- + so what was used belongs in the specification, and in the result. + :param nfree: + Free parameters expected once every constraint above is applied. + Not an input: the runner asserts it. A change in how SCALCS handles + constraints then shows up as a refusal to start rather than as a + puzzling set of estimates a day later. + """ + + sample: str = None + mec_file: str = None + mec_number: int = None + rates: dict = field(default_factory=dict) + fixed: dict = field(default_factory=dict) + limits: dict = field(default_factory=dict) + mr: str = None + mr_cycle: int = 0 + ec50: dict = None + nfree: int = None + + KEYS = ('sample', 'mec_file', 'mec_number', 'rates', 'fixed', 'limits', + 'mr', 'mr_cycle', 'ec50', 'nfree') + + @classmethod + def from_dict(cls, d, where='mechanism'): + """ Build one from a parsed table, checking as it goes. """ + _unexpected(d, cls.KEYS, where) + + limits = {} + for name, pair in (d.get('limits') or {}).items(): + if not isinstance(pair, (list, tuple)) or len(pair) != 2: + raise SpecError("{0}.limits.{1}: expected [lower, upper], got " + "{2!r}".format(where, name, pair)) + lower = _number(pair[0], "{0}.limits.{1}[0]".format(where, name)) + upper = _number(pair[1], "{0}.limits.{1}[1]".format(where, name)) + if not lower < upper: + raise SpecError("{0}.limits.{1}: lower ({2:g}) must be below " + "upper ({3:g})".format(where, name, lower, + upper)) + limits[str(name)] = (lower, upper) + + ec50 = d.get('ec50') + if ec50 is not None: + if not isinstance(ec50, dict): + raise SpecError("{0}.ec50: expected a table with rate and " + "value, got {1!r}".format(where, ec50)) + _unexpected(ec50, ('rate', 'value'), where + '.ec50') + for required in ('rate', 'value'): + if required not in ec50: + raise SpecError("{0}.ec50: {1} is required" + .format(where, required)) + ec50 = {'rate': str(ec50['rate']), + 'value': _number(ec50['value'], where + '.ec50.value')} + + return cls( + sample=None if d.get('sample') is None else str(d['sample']), + mec_file=None if d.get('mec_file') is None else str(d['mec_file']), + mec_number=(None if d.get('mec_number') is None + else int(d['mec_number'])), + rates=_named_numbers(d.get('rates'), where + '.rates'), + fixed=_named_numbers(d.get('fixed'), where + '.fixed'), + limits=limits, + mr=None if d.get('mr') is None else str(d['mr']), + mr_cycle=int(d.get('mr_cycle', 0)), + ec50=ec50, + nfree=None if d.get('nfree') is None else int(d['nfree']), + ).validate(where) + + def validate(self, where='mechanism'): + """ Raise :py:exc:`SpecError` unless this mechanism can be built. """ + if (self.sample is None) == (self.mec_file is None): + raise SpecError("{0}: give exactly one of sample and mec_file" + .format(where)) + if self.mec_file is None and self.mec_number is not None: + raise SpecError("{0}.mec_number: means nothing without mec_file" + .format(where)) + both = sorted(set(self.rates) & set(self.fixed)) + if both: + raise SpecError( + "{0}: {1} {2} in both rates and fixed. A fixed rate takes its " + "value from fixed; naming it in rates as well says the same " + "thing twice and invites the two to disagree.".format( + where, ", ".join(both), + "is" if len(both) == 1 else "are")) + if self.ec50 is not None: + if self.ec50['value'] <= 0.0: + raise SpecError("{0}.ec50.value: must be positive, not {1}" + .format(where, self.ec50['value'])) + if self.ec50['rate'] in self.fixed: + raise SpecError( + "{0}: {1} is both fixed and computed from the EC50; it " + "cannot be both".format(where, self.ec50['rate'])) + if self.nfree is not None and self.nfree < 1: + raise SpecError("{0}.nfree: must be at least 1, not {1}" + .format(where, self.nfree)) + return self + + def as_dict(self): + d = {} + for name in ('sample', 'mec_file', 'mec_number', 'mr', 'nfree'): + value = getattr(self, name) + if value is not None: + d[name] = value + if self.mr is not None and self.mr_cycle: + d['mr_cycle'] = self.mr_cycle + for name in ('rates', 'fixed'): + if getattr(self, name): + d[name] = dict(getattr(self, name)) + if self.limits: + d['limits'] = {k: list(v) for k, v in self.limits.items()} + if self.ec50 is not None: + d['ec50'] = dict(self.ec50) + return d + + +@dataclass(frozen=True) +class SearchSpec: + """ How to search, and how hard. + + :param method: + ``"simplex"`` is + :py:func:`~HJCFIT.likelihood.optimization.simplex_hjc`, HJCFIT's own + and the search that produced every published result. ``"scipy"`` is + SciPy's Nelder-Mead with the restart loop below; reach for it when + the starting point is poor, because a regular simplex needs every + one of its vertices to be evaluable. + :param log_params: + Search the logarithms of the rate constants. HJCFIT's default, three + to four times faster, and it cannot produce a negative rate. It + belongs in the specification rather than in the driver because it + changes which maximum a search falls into: Colquhoun, Hatton & + Hawkes (2003) say of their Figure 2 fits that the rate constants + themselves were the free parameters, and the same set of fits run + over logarithms reaches the second solution at a quite different + rate. + :param restarts: + ``"scipy"`` only: extra runs from the previous solution, stopping + once one gains less than *fatol*. A single SciPy pass leaves fits + stranded between maxima. Measured over the twelve stranded fits of + one reproduction scenario, restarting gained 106.3 + log\\ :sub:`10` units in total, resolved five of them onto a real + maximum, and never lost any. + :param maxfev: Evaluation budget; both searches. + :param maxiter: ``"scipy"`` only. + :param xatol: ``"scipy"`` only. + :param fatol: ``"scipy"`` only. + :param options: + Passed straight to whichever search is chosen, for deliberate + departures from its own defaults. + """ + + method: str = 'simplex' + log_params: bool = True + restarts: int = 0 + maxfev: int = 20000 + maxiter: int = 20000 + xatol: float = 1e-4 + fatol: float = 1e-4 + options: dict = field(default_factory=dict) + + KEYS = ('method', 'log_params', 'restarts', 'maxfev', 'maxiter', + 'xatol', 'fatol', 'options') + + @classmethod + def from_dict(cls, d, where='search'): + """ Build one from a parsed table, checking as it goes. """ + _unexpected(d, cls.KEYS, where) + options = d.get('options') or {} + if not isinstance(options, dict): + raise SpecError("{0}.options: expected a table, got {1!r}" + .format(where, options)) + return cls( + method=str(d.get('method', 'simplex')), + log_params=bool(d.get('log_params', True)), + restarts=int(d.get('restarts', 0)), + maxfev=int(d.get('maxfev', 20000)), + maxiter=int(d.get('maxiter', 20000)), + xatol=_number(d.get('xatol', 1e-4), where + '.xatol'), + fatol=_number(d.get('fatol', 1e-4), where + '.fatol'), + options=dict(options), + ).validate(where) + + def validate(self, where='search'): + """ Raise :py:exc:`SpecError` unless this search can be run. """ + if self.method not in METHODS: + raise SpecError("{0}.method: must be one of {1}, not {2!r}".format( + where, " or ".join(METHODS), self.method)) + if self.restarts and self.method != 'scipy': + raise SpecError( + "{0}.restarts: only the scipy search takes restarts. " + "simplex_hjc has its own restart rule, capped by nresmax, " + "which goes in options.".format(where)) + if self.restarts < 0: + raise SpecError("{0}.restarts: cannot be negative, got {1}" + .format(where, self.restarts)) + for name in ('maxfev', 'maxiter'): + if getattr(self, name) < 1: + raise SpecError("{0}.{1}: must be at least 1, not {2}".format( + where, name, getattr(self, name))) + return self + + def as_dict(self): + d = {'method': self.method, 'log_params': self.log_params, + 'maxfev': self.maxfev} + if self.method == 'scipy': + d.update(restarts=self.restarts, maxiter=self.maxiter, + xatol=self.xatol, fatol=self.fatol) + if self.options: + d['options'] = dict(self.options) + return d + + +@dataclass(frozen=True) +class FitSpec: + """ A whole fit, as data. + + :param data: + One :py:class:`DataSpec` per record. More than one means several + records contributing to a single likelihood: the rate constants are + shared and the concentrations are not. + :param mechanism: A :py:class:`MechanismSpec`. + :param search: A :py:class:`SearchSpec`. + :param title: One line for the person reading the result. + + Read one with :py:meth:`from_toml`, write one with + :py:meth:`to_toml`, and run one with + :py:func:`HJCFIT.likelihood.runner.run`. + """ + + data: tuple + mechanism: MechanismSpec = field(default_factory=MechanismSpec) + search: SearchSpec = field(default_factory=SearchSpec) + title: str = '' + + KEYS = ('data', 'mechanism', 'search', 'title') + + # -- reading ----------------------------------------------------------- + + @classmethod + def from_dict(cls, d, where='spec'): + """ Build one from a parsed document. + + :raises SpecError: naming the field, on anything unusable. + """ + if not isinstance(d, dict): + raise SpecError("{0}: expected a table, got {1!r}".format(where, d)) + _unexpected(d, cls.KEYS, where) + records = d.get('data') + if not records: + raise SpecError( + "{0}: at least one [[data]] section is needed -- a fit " + "without a record is not a fit".format(where)) + if isinstance(records, dict): # a single [data] rather than [[data]] + records = [records] + if not isinstance(records, (list, tuple)): + raise SpecError("{0}.data: expected one or more [[data]] " + "sections, got {1!r}".format(where, records)) + return cls( + data=tuple(DataSpec.from_dict(r, "data[{0}]".format(i)) + for i, r in enumerate(records)), + mechanism=MechanismSpec.from_dict(d.get('mechanism') or {}), + search=SearchSpec.from_dict(d.get('search') or {}), + title=str(d.get('title', '')), + ).validate(where) + + @classmethod + def from_toml(cls, path): + """ Read a specification from a TOML file. + + :param path: The file. Read as bytes, which is what TOML wants: + the format is UTF-8 by definition, so letting the platform + choose an encoding could only get it wrong. + :raises SpecError: on a file TOML cannot parse, or a field this + module cannot use. Both carry the file name. + """ + with open(path, 'rb') as handle: + raw = handle.read() + try: + return cls.from_dict(load_toml_bytes(raw)) + except SpecError as error: + raise SpecError("{0}: {1}".format(path, error)) + + def validate(self, where='spec'): + """ Raise :py:exc:`SpecError` unless the whole thing hangs together. """ + for i, record in enumerate(self.data): + record.validate("data[{0}]".format(i)) + self.mechanism.validate() + self.search.validate() + # Two records at the same concentration are not an error -- two patches + # at one concentration is an ordinary experiment, and the reproduction + # of Colquhoun, Hatton & Hawkes (2003) fits three concentrations of + # which none repeats. What would be an error is a mechanism with no + # concentration dependence fitted to several concentrations, and that + # cannot be seen from here: it needs the mechanism built. + return self + + # -- writing ----------------------------------------------------------- + + def as_dict(self): + """ A plain dictionary, ready for :py:func:`to_toml` or for JSON. + + Defaults are left out, so a specification written back is the + shortest file that means the same thing. + """ + d = {} + if self.title: + d['title'] = self.title + d['data'] = [r.as_dict() for r in self.data] + d['mechanism'] = self.mechanism.as_dict() + d['search'] = self.search.as_dict() + return d + + def to_toml(self): + """ The specification as TOML text. + + Written by hand rather than with a library: ``tomllib`` reads and + does not write, and the alternative is a dependency to emit forty + lines of a format with four value types in it. + """ + lines = [] + if self.title: + lines.append('title = {0}'.format(_toml_value(self.title))) + lines.append('') + for record in self.data: + lines.append('[[data]]') + lines.extend(_toml_table(record.as_dict())) + lines.append('') + for name in ('mechanism', 'search'): + table = getattr(self, name).as_dict() + lines.append('[{0}]'.format(name)) + lines.extend(_toml_table(table, name)) + lines.append('') + return "\n".join(lines).rstrip() + "\n" + + def write_toml(self, path): + """ Write :py:meth:`to_toml` to *path*, as UTF-8. """ + with open(path, 'w', encoding='utf-8', newline='\n') as handle: + handle.write(self.to_toml()) + + # -- looking at it ----------------------------------------------------- + + def __str__(self): + lines = [self.title or 'a fit'] + for record in self.data: + lines.append(' ' + str(record)) + mech = self.mechanism + source = (mech.sample if mech.sample is not None else + '{0}{1}'.format(mech.mec_file, + '' if mech.mec_number is None + else ' #{0}'.format(mech.mec_number))) + lines.append(' mechanism: {0}'.format(source)) + for label, table in (('guess', mech.rates), ('fixed', mech.fixed)): + if table: + lines.append(' {0}: {1}'.format(label, ", ".join( + '{0} = {1:g}'.format(k, v) + for k, v in sorted(table.items())))) + if mech.mr is not None: + lines.append(' {0} from microscopic reversibility' + .format(mech.mr)) + if mech.ec50 is not None: + lines.append(' {0} from an EC50 of {1:g} uM'.format( + mech.ec50['rate'], mech.ec50['value'] * 1e6)) + if mech.limits: + lines.append(' limits on {0}'.format( + ", ".join(sorted(mech.limits)))) + lines.append(' search: {0}, over {1}, up to {2} evaluations'.format( + self.search.method, + 'logarithms' if self.search.log_params else 'rate constants', + self.search.maxfev)) + return "\n".join(lines) + + +# ---------------------------------------------------------------- TOML I/O + +def load_toml_bytes(raw): + """ Parse TOML from bytes, whichever parser this Python has. + + ``tomllib`` is in the standard library from 3.11; ``tomli`` is the + same parser by the same author and covers 3.10, and is in the + ``[fitting]`` extra for that reason. + + :raises SpecError: on a file TOML cannot parse, or when neither + parser is importable. + """ + try: + import tomllib as toml + except ImportError: # Python 3.10 + try: + import tomli as toml + except ImportError: + raise SpecError( + "reading a specification needs a TOML parser. Python 3.11 and " + "later have one in the standard library; on 3.10, install " + "tomli -- pip install 'hjcfit[fitting]' does that.") + try: + return toml.loads(raw.decode('utf-8')) + except UnicodeDecodeError as error: + raise SpecError("not UTF-8, which TOML must be: {0}".format(error)) + except toml.TOMLDecodeError as error: + raise SpecError("not valid TOML: {0}".format(error)) + + +def _toml_value(value): + """ One TOML scalar, or an array of them. """ + if isinstance(value, bool): + return 'true' if value else 'false' + if isinstance(value, int): + return repr(value) + if isinstance(value, float): + # repr round-trips, and TOML accepts every form it produces except the + # bare infinities and nan, which cannot be a rate constant or a time. + text = repr(value) + if text in ('inf', '-inf', 'nan'): + raise SpecError("{0} cannot go in a specification".format(text)) + return text + if isinstance(value, (list, tuple)): + return '[' + ", ".join(_toml_value(v) for v in value) + ']' + text = str(value) + return '"' + text.replace('\\', '\\\\').replace('"', '\\"') + '"' + + +def _toml_key(name): + """ A bare key where TOML allows one, a quoted key otherwise. + + Rate constants are called things like ``k(-1)`` and ``2k(+2)``, which + are not bare keys, so most of these come back quoted. + """ + if name and all(c.isalnum() or c in '-_' for c in name): + return name + return _toml_value(name) + + +def _toml_table(table, prefix=''): + """ One table as lines, with its sub-tables after its own keys. + + TOML reads every key after a ``[header]`` as belonging to it, so a + scalar written below a sub-table would silently join the sub-table. + Scalars therefore come first, which is not a matter of taste. + """ + lines = [ + '{0} = {1}'.format(_toml_key(k), _toml_value(v)) + for k, v in table.items() if not isinstance(v, dict) + ] + for key, value in table.items(): + if isinstance(value, dict): + if not value: + continue + header = '{0}.{1}'.format(prefix, key) if prefix else key + lines.append('') + lines.append('[{0}]'.format(header)) + lines.extend(_toml_table(value, header)) + return lines + + +#: A specification to start from, as written by ``hjcfit template``. It fits +#: the CH82 sample record that ships with HJCFIT, so it runs as it stands. +TEMPLATE = '''\ +# A fit specification. Everything a fit needs, and nothing that runs. +# +# hjcfit check this-file.toml say what it would do, and do nothing +# hjcfit fit this-file.toml run it +# +# As it stands this fits the CH82 sample record that ships with HJCFIT, so it +# works before you change anything. Then point `record` at your own .scn file. + +title = "CH82 sample record at 100 nM" + +# One [[data]] section per record. Several means several records contributing +# to one likelihood: the rate constants are shared, the concentrations are not. +[[data]] +record = "CH82" # an .scn file, or a sample: CH82, CO, CCO +conc = 1e-07 # M +tres = 0.0001 # dead time to impose, s +tcrit = 0.004 # critical shut time dividing the record into groups, s +vectors = "chs" # "chs" between groups of one channel's activations; + # "equilibrium" for clusters, or for a whole record + # (omit tcrit for that, and vectors must be equilibrium) + +[mechanism] +sample = "CH82" # a factory in scalcs.samples.samples +# mec_file = "mechs.mec" # or a DCprogs .mec file, instead of sample +# mec_number = 1 # which mechanism in it +# mr = "k(-1)" # a rate computed from microscopic reversibility +# nfree = 8 # free parameters expected; asserted, not applied + +# The initial guess, by rate name. Names left out keep the sample's own value, +# so this can be one line. `hjcfit check` prints every name the mechanism has. +[mechanism.rates] +beta1 = 15.0 +beta2 = 15000.0 + +# Held at this value throughout. +# [mechanism.fixed] +# "k(-1)" = 2000.0 + +# A rate driven outside its limits is reset to the limit before the likelihood +# sees it, and nothing announces it. `hjcfit fit` reports any rate that ends +# against one, because such a rate is not an estimate. +# [mechanism.limits] +# beta1 = [0.1, 100000.0] + +# One rate computed from an independently measured EC50 at every iteration, +# instead of being fitted (Colquhoun, Hatton & Hawkes 2003, p. 702). +# [mechanism.ec50] +# rate = "k(+2)" +# value = 3.3e-06 + +[search] +method = "simplex" # HJCFIT's own; "scipy" is Nelder-Mead with restarts +log_params = true # search the logarithms: faster, and cannot go negative +maxfev = 20000 +''' diff --git a/likelihood/python/fitting.py b/likelihood/python/fitting.py index 064054f..3a54a9d 100644 --- a/likelihood/python/fitting.py +++ b/likelihood/python/fitting.py @@ -147,9 +147,10 @@ def __str__(self): vectors = ("CHS vectors" if self.tcrit is not None else "equilibrium vectors") raw = "" if self.n_raw is None else "{0} -> ".format(self.n_raw) - return ("{0}: {1}{2} intervals at {3:g} us -> {4} groups, {5} " - "openings, {6}".format(conc, raw, self.n_intervals, + return ("{0}: {1}{2} intervals at {3:g} us -> {4} group{5}, {6} " + "openings, {7}".format(conc, raw, self.n_intervals, self.tres * 1e6, len(self.groups), + "" if len(self.groups) == 1 else "s", self.n_openings, vectors)) diff --git a/likelihood/python/runner.py b/likelihood/python/runner.py new file mode 100644 index 0000000..9c0a275 --- /dev/null +++ b/likelihood/python/runner.py @@ -0,0 +1,542 @@ +######################## +# HJCFIT computes missed-events likelihood as described in +# Hawkes, Jalali and Colquhoun (1990, 1992) +# +# Copyright (C) 2013 University College London +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +######################### + +""" Carrying out what a :py:class:`~HJCFIT.likelihood.fitspec.FitSpec` says. + + :py:mod:`~HJCFIT.likelihood.fitspec` describes a fit and executes nothing. + This module is the half that executes: it reads the records through + ``dcio``, builds the mechanism through ``scalcs``, runs the search through + :py:class:`~HJCFIT.likelihood.fitting.HJCFitter`, and writes a result with + enough beside it to be believed a year later. + + :: + + from HJCFIT.likelihood.fitspec import FitSpec + from HJCFIT.likelihood.runner import run + + outcome = run(FitSpec.from_toml("ch82.toml")) + print(outcome.result) + + Needs the ``[fitting]`` extra: ``pip install 'hjcfit[fitting]'``. + + Three things this module does that a hand-written script usually does not, + each because it has gone wrong: + + * **It names rates rather than numbering them.** SCALCS addresses rates by + index into ``mec.Rates``. :py:func:`build_mechanism` is the one place a + name becomes an index, and it lists the names that exist when one does + not. + * **It reports rates that end against a limit.** A rate outside its limits + is reset before the likelihood sees it and nothing says so, which means + a fit can end with a rate sitting exactly on a bound and look like an + estimate. Note that SCALCS gives every rate default limits -- + 1e-15 to 1e9 for a concentration-dependent rate and 1e-15 to 1e6 for the + rest -- so this can happen in a specification that sets no limits at all. + * **It asserts the free-parameter count** when the specification states + one, before the search starts rather than after it. +""" +__docformat__ = "restructuredtext en" +__all__ = ['Outcome', 'build_mechanism', 'load_records', 'run', + 'against_limits', 'provenance', 'result_as_dict', 'write_result', + 'RunnerError'] + +import json +import os +import platform +import socket +import sys +from dataclasses import dataclass, field +from datetime import datetime, timezone + +import numpy as np + +from .fitting import HJCFitter, Record, trim_to_openings + + +class RunnerError(RuntimeError): + """ A specification that cannot be carried out. + + Distinct from :py:exc:`~HJCFIT.likelihood.fitspec.SpecError`, which is + a specification that cannot be *read*. A name that no rate of the + chosen mechanism carries is well-formed and still impossible, and it + is this. + """ + + +def _need(module, what): + """ Import an extra, and say what to install when it is not there. """ + try: + __import__(module) + except ImportError as error: + raise RunnerError( + "{0} needs {1}, which is not a dependency of HJCFIT -- the " + "likelihood itself wants only numpy and scipy. Install the extra: " + "pip install 'hjcfit[fitting]'. (Original error: {2})" + .format(what, module.split('.')[0], error)) + return sys.modules[module] + + +# ------------------------------------------------------------------ records + +def _resolve_record(name): + """ A path, or the name of a sample record shipped with HJCFIT. + + The same lookup :py:func:`HJCFIT.read_idealized_bursts` does, so a + specification can say ``record = "CH82"`` and run anywhere the wheel + is installed. + """ + if os.path.exists(name): + return name + from os.path import abspath, dirname, join, splitext + # HJCFIT/likelihood/runner.py -> HJCFIT/data + sample = join(dirname(dirname(abspath(__file__))), 'data', + '{0}.scn'.format(splitext(name)[0])) + if os.path.exists(sample): + return sample + raise RunnerError( + "no file or sample record {0!r}. The samples that ship with HJCFIT " + "are CH82, CO and CCO; anything else has to be a path to an .scn " + "file.".format(name)) + + +def load_records(spec): + """ Read every record a specification names. + + :param spec: A :py:class:`~HJCFIT.likelihood.fitspec.FitSpec`. + :returns: A list of :py:class:`~HJCFIT.likelihood.fitting.Record`, in + the order the specification gives them. + :raises RunnerError: on a record that cannot be found or read. + + Two paths, because a record divided into groups and a record fitted + whole are genuinely different things: + + * with a critical time, the record's **periods** are segmented by + :py:func:`dcio.analysis.bursts_from_record`. Periods rather than + resolved intervals because imposing a dead time emits a fresh open + interval at every change of fitted amplitude, so a record idealised + with sub-conductance levels does not alternate open and shut, and + the likelihood is a product of matrices that must. + :py:func:`HJCFIT.read_idealized_bursts` is the other caller of that + pairing and carries the argument in full; a test here requires the + two to produce identical groups rather than similar ones. + * without one, the whole record is a single group, trimmed to start + and end on an opening. That assumes one channel for the length of + the record, which is why it is the exception. + + The critical time is used twice and means two different things, which + is why the specification separates them: its magnitude divides the + record, and the choice of vectors decides what the likelihood is told. + Groups cut with ``vectors = "equilibrium"`` -- clusters at a + desensitising concentration -- pass None as *tcritical* while still + being groups. + """ + scn = _need('dcio.formats.scn', 'reading a record') + analysis = _need('dcio.analysis', 'reading a record') + + records = [] + for i, data in enumerate(spec.data): + path = _resolve_record(data.record) + try: + raw = scn.read(path) + screc = analysis.from_scn(raw, tres=data.tres) + except Exception as error: + raise RunnerError("data[{0}]: cannot read {1}: {2}: {3}".format( + i, path, type(error).__name__, error)) + + if data.tcrit is None: + periods = screc.periods + groups = (tuple(trim_to_openings(periods.intervals, + periods.amplitudes)),) + else: + groups = tuple( + tuple(np.asarray(g, dtype=float)) + for g in analysis.bursts_from_record(screc, data.tcrit, + intervals_only=True)) + if not groups or not groups[0]: + raise RunnerError( + "data[{0}]: {1} yields no intervals at a dead time of " + "{2:g} us{3}".format(i, path, data.tres * 1e6, + "" if data.tcrit is None else + " and a critical time of {0:g} ms" + .format(data.tcrit * 1e3))) + + records.append(Record( + conc=data.conc, + groups=groups, + tres=data.tres, + # The magnitude divided the record above; here it selects the + # vectors, and None is what asks for equilibrium ones. + tcrit=data.tcrit if data.vectors == 'chs' else None, + n_raw=int(len(raw.intervals)), + ).check()) + return records + + +# --------------------------------------------------------------- mechanism + +def _indices(mec, name, where): + """ Every index in ``mec.Rates`` whose rate carries *name*. """ + found = [i for i, rate in enumerate(mec.Rates) if rate.name == name] + if not found: + raise RunnerError( + "{0}: this mechanism has no rate called {1!r}. It has: {2}" + .format(where, name, + ", ".join(r.name for r in mec.Rates))) + return found + + +def _one_index(mec, name, where): + """ The single index of *name*, for an interface that takes one. """ + found = _indices(mec, name, where) + if len(found) > 1: + raise RunnerError( + "{0}: {1} names {2} rates of this mechanism, and this setting " + "applies to one".format(where, name, len(found))) + return found[0] + + +def _build_bare(mech): + """ The mechanism named by the specification, before any constraint. """ + if mech.sample is not None: + samples = _need('scalcs.samples.samples', 'building a mechanism') + factory = getattr(samples, mech.sample, None) + if factory is None or not callable(factory): + available = sorted( + name for name, value in vars(samples).items() + if callable(value) and not name.startswith('_') + and name not in ('apply_rates',)) + raise RunnerError( + "mechanism.sample: scalcs.samples.samples has no {0!r}. It " + "has: {1}".format(mech.sample, ", ".join(available))) + built = factory() + # Some factories -- the independent-binding loaders -- return the + # mechanism and the constraints they applied. The mechanism is first. + return built[0] if isinstance(built, tuple) else built + + scalcsio = _need('scalcs.scalcsio', 'reading a .mec file') + if not os.path.exists(mech.mec_file): + raise RunnerError("mechanism.mec_file: no such file {0!r}" + .format(mech.mec_file)) + _version, meclist, _max = scalcsio.mec_get_list(mech.mec_file) + if not meclist: + raise RunnerError("mechanism.mec_file: {0} holds no mechanisms" + .format(mech.mec_file)) + if mech.mec_number is None: + if len(meclist) > 1: + listing = "\n".join( + " {0} {1} {2}".format(entry[1], entry[2], entry[3]) + for entry in meclist) + raise RunnerError( + "mechanism.mec_number: {0} holds {1} mechanisms, so which one " + "has to be said:\n{2}".format(mech.mec_file, len(meclist), + listing)) + chosen = meclist[0] + else: + matching = [e for e in meclist if e[1] == mech.mec_number] + if not matching: + raise RunnerError( + "mechanism.mec_number: {0} has no mechanism {1}; it has {2}" + .format(mech.mec_file, mech.mec_number, + ", ".join(str(e[1]) for e in meclist))) + chosen = matching[-1] # the last rate set saved for it + return scalcsio.mec_load(mech.mec_file, chosen[0]) + + +def build_mechanism(spec): + """ Build the mechanism a specification names, with its constraints. + + :param spec: A :py:class:`~HJCFIT.likelihood.fitspec.FitSpec`. + :returns: A ``scalcs.mechanism.Mechanism``, carrying the initial + guess, ready to be fitted. + :raises RunnerError: on a name no rate has, or a free-parameter count + that does not match ``mechanism.nfree``. + + The order matters and is the reproduction's: the guess, then the fixed + rates, then the limits, then microscopic reversibility, then the + constraints are updated, and only then the EC50 -- which is computed + *from* the other rates and so has to come last. + """ + mech = spec.mechanism + mec = _build_bare(mech) + + for name, value in mech.rates.items(): + for i in _indices(mec, name, 'mechanism.rates'): + mec.Rates[i].rateconstants = value + + for name, value in mech.fixed.items(): + for i in _indices(mec, name, 'mechanism.fixed'): + mec.Rates[i].rateconstants = value + mec.Rates[i].fixed = True + + for name, (lower, upper) in mech.limits.items(): + for i in _indices(mec, name, 'mechanism.limits'): + mec.Rates[i].limits = [[lower, upper]] + + if mech.mr is not None: + if not getattr(mec, 'Cycles', None): + raise RunnerError( + "mechanism.mr: this mechanism has no cycles, so no rate of it " + "is determined by microscopic reversibility") + if not 0 <= mech.mr_cycle < len(mec.Cycles): + raise RunnerError( + "mechanism.mr_cycle: {0} is not a cycle of this mechanism, " + "which has {1}".format(mech.mr_cycle, len(mec.Cycles))) + index = _one_index(mec, mech.mr, 'mechanism.mr') + rate = mec.Rates[index] + cycle = mec.Cycles[mech.mr_cycle] + # Checked here rather than left to SCALCS, which warns on stderr and + # carries on. What it does then is worse than nothing: it marks the + # rate constrained, so it leaves the free-parameter list, while the + # cycle's own constraint stays on whichever pair it was already on. + # The named rate is not computed from anything -- it is silently + # frozen at its initial guess, and the fit has one parameter fewer than + # the specification says. Measured on CH82 with mr = "k(-1)", which is + # not in its cycle: eight free parameters became seven, k(-1) stayed at + # the guess, and the microscopic-reversibility constraint remained on + # 2k*(-2). + states = list(getattr(cycle, 'states', [])) + outside = [s for s in (rate.State1.name, rate.State2.name) + if s not in states] + if outside: + raise RunnerError( + "mechanism.mr: {0} connects {1} and {2}, and {3} not in cycle " + "{4} ({5}). A rate outside the cycle cannot be determined by " + "microscopic reversibility round it.".format( + mech.mr, rate.State1.name, rate.State2.name, + "{0} is".format(outside[0]) if len(outside) == 1 + else "neither is", + mech.mr_cycle, ", ".join(states))) + + # A cycle determines exactly one rate, so any other rate of this cycle + # already flagged is stale and has to be released. SCALCS does not do + # this: set_mr moves the cycle's own constraint and leaves every + # previous rate's mr flag set, and its update_mr carries the comment + # "TODO: check for consistency between cycle.mrconstr and rate.mr". + # The consequence is the same silent freeze as above, in the supported + # path rather than a misuse of it. On CH82, whose sample already has + # 2k*(-2) under microscopic reversibility, asking for mr = "beta1" + # left seven free parameters rather than eight, with 2k*(-2) neither + # fitted nor computed from the cycle -- stuck at 0.411, the value the + # sample happened to carry. + for i, other in enumerate(mec.Rates): + if i == index or not getattr(other, 'mr', False): + continue + if (other.State1.name in states and other.State2.name in states): + other.mr = False + + mec.set_mr(True, index, mech.mr_cycle) + + mec.update_constrains() + + if mech.ec50 is not None: + mec.set_EC50_constraint( + _one_index(mec, mech.ec50['rate'], 'mechanism.ec50.rate'), + mech.ec50['value']) + + free = list(mec.get_free_parameter_names()) + if not free: + raise RunnerError( + "every rate of this mechanism is fixed or constrained, so there " + "is nothing to fit") + if mech.nfree is not None and len(free) != mech.nfree: + raise RunnerError( + "mechanism.nfree says {0} free parameters; the mechanism has {1}: " + "{2}".format(mech.nfree, len(free), ", ".join(free))) + return mec + + +# ------------------------------------------------------------------ the fit + +def against_limits(mec, tolerance=1e-6): + """ Free rates that ended on one of their limits. + + :param mec: A mechanism, after a fit. + :param tolerance: Relative closeness that counts as "on". + :returns: A list of ``(name, value, 'lower'|'upper', limit)``. + + A rate driven outside its limits is reset to the limit before the + likelihood sees it, which is what HJCFIT did and what + :py:class:`~HJCFIT.likelihood.fitting.HJCFitter` does. Nothing + announces it. A rate sitting on a bound is not an estimate of that + rate -- it is a statement that the likelihood wanted to go somewhere + the model forbids -- so it has to be reported separately from the + fitted values, and this is what does that. + """ + on = [] + for rate in mec.Rates: + if not rate.is_free or not rate.limits: + continue + value = float(np.ravel(rate.rateconstants)[0]) + lower, upper = rate.limits[0][0], rate.limits[0][1] + for which, limit in (('lower', lower), ('upper', upper)): + if limit and abs(value - limit) <= tolerance * abs(limit): + on.append((rate.name, value, which, float(limit))) + return on + + +@dataclass +class Outcome: + """ Everything one run produced. + + :param spec: The specification it was run from. + :param result: The :py:class:`~HJCFIT.likelihood.fitting.FitResult`. + :param records: The records as they were fitted. + :param mec: The mechanism, carrying the **fitted** rate constants -- + which is what predicted distributions have to be drawn from. + :param limited: + What :py:func:`against_limits` found. Empty is the good case. + :param provenance: What :py:func:`provenance` recorded. + """ + + spec: object + result: object + records: list + mec: object = field(repr=False, default=None) + limited: list = field(default_factory=list) + provenance: dict = field(default_factory=dict, repr=False) + + def __str__(self): + lines = [str(r) for r in self.records] + lines.append(str(self.result)) + if self.limited: + lines.append("") + lines.append("Rates that ended on a limit -- these are not " + "estimates:") + for name, value, which, limit in self.limited: + lines.append(" {0:9s} {1:12.4g} at its {2} limit of {3:g}" + .format(name, value, which, limit)) + return "\n".join(lines) + + +def run(spec, verbose=False, store_path=False, x0=None): + """ Read the records, build the mechanism, fit, and report. + + :param spec: A :py:class:`~HJCFIT.likelihood.fitspec.FitSpec`. + :param verbose: Print the records and the starting likelihood first. + :param store_path: Keep the best point of every iteration. + :param x0: Starting point, in the space being searched. Defaults to + the mechanism's own guess, which is what the specification says. + :returns: An :py:class:`Outcome`. + """ + records = load_records(spec) + mec = build_mechanism(spec) + search = spec.search + + fitter = HJCFitter(mec, records, log_params=search.log_params, + store_path=store_path) + if verbose: + # Not the records: Outcome.__str__ prints those, so that a run with + # verbose off still says what was fitted. + print("{0} free parameters: {1}".format( + len(mec.get_free_parameter_names()), + ", ".join(mec.get_free_parameter_names()))) + print("log10L at the guess: {0:.4f}".format(fitter.log10_likelihood())) + + result = fitter.fit(x0=x0, search=search.method, + restarts=search.restarts, maxfev=search.maxfev, + maxiter=search.maxiter, xatol=search.xatol, + fatol=search.fatol, options=search.options or None) + return Outcome(spec=spec, result=result, records=records, mec=fitter.mec, + limited=against_limits(fitter.mec), + provenance=provenance()) + + +# ---------------------------------------------------------------- the record + +def provenance(): + """ Where and with what this was computed. + + Package versions, interpreter, host and time. Cheap, and the + difference between a number in a paper that can be traced and one that + cannot: the reproduction of Colquhoun, Hatton & Hawkes (2003) attaches + one of these to every cached result, and it is the reason a + disagreement there can be chased to a version rather than argued about. + + Deliberately not a git description of the caller's working tree. + Shelling out to git costs the best part of a second, needs a + repository to be there, and says nothing at all about a fit run from + an installed wheel -- which is how the people this is for will run it. + """ + from importlib.metadata import PackageNotFoundError + from importlib.metadata import version as _version + + versions = {} + for package in ('hjcfit', 'dcio', 'scalcs', 'numpy', 'scipy', + 'matplotlib'): + try: + versions[package] = _version(package) + except PackageNotFoundError: + versions[package] = None + return { + 'when': datetime.now(timezone.utc).isoformat(timespec='seconds'), + 'versions': versions, + 'python': sys.version.split()[0], + 'platform': platform.platform(), + 'host': socket.gethostname(), + } + + +def result_as_dict(outcome): + """ An outcome as plain JSON-ready data. + + The specification is written back into it, so the file is a complete + account: what was asked for, what came out, and what it was computed + with. Nothing here holds the records -- they are thousands of floats + and they are in the .scn file the specification names. + """ + result = outcome.result + return { + 'hjcfit_result': 1, + 'spec': outcome.spec.as_dict(), + 'fit': { + 'log10_likelihood': result.log10_likelihood, + 'ln_likelihood': result.log10_likelihood * float(np.log(10.0)), + 'free_names': list(result.free_names), + 'free_values': [float(v) for v in result.free_values], + 'rates': {k: float(v) for k, v in result.rates.items()}, + 'nevals': int(result.nevals), + 'niter': int(result.niter), + 'nfailures': int(result.nfailures), + 'seconds': float(result.seconds), + 'success': bool(result.success), + 'message': str(result.message), + }, + 'records': [ + {'conc': r.conc, 'tres': r.tres, 'tcrit': r.tcrit, + 'groups': len(r.groups), 'intervals': r.n_intervals, + 'openings': r.n_openings, 'raw_intervals': r.n_raw} + for r in outcome.records + ], + 'at_limits': [ + {'rate': name, 'value': value, 'which': which, 'limit': limit} + for name, value, which, limit in outcome.limited + ], + 'provenance': outcome.provenance, + } + + +def write_result(path, outcome): + """ Write :py:func:`result_as_dict` to *path* as JSON. + + JSON rather than TOML: this is written and read by programs, and TOML + has no null, which every absent version and every equilibrium-vector + critical time needs. + """ + with open(path, 'w', encoding='utf-8', newline='\n') as handle: + json.dump(result_as_dict(outcome), handle, indent=2, sort_keys=False) + handle.write("\n") + return path diff --git a/pyproject.toml b/pyproject.toml index 361d2f4..f52453b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -82,6 +82,9 @@ fitting = [ # 1.1.0 does not have: the looser pin let the extra # resolve and the fit then fail on a missing attribute "matplotlib>=3.6", # the displays; nothing in the likelihood draws + # A fit specification is TOML, which is in the standard library from 3.11. + # tomli is the same parser by the same author, for 3.10. + "tomli>=1.1; python_version < '3.11'", ] # EKDIST -- exponential mixtures and t_crit estimation -- is deliberately not # an extra. It is not on PyPI, so declaring it would make @@ -89,6 +92,13 @@ fitting = [ # EKDIST's own dcio dependency was unresolvable in exactly that way until dcio # was published. Install it from source when a display needs it. +# `hjcfit fit spec.toml` -- the command over a fit specification. It needs the +# [fitting] extra to do anything, and says so rather than failing on an import; +# `pip install hjcfit` alone still installs the script, which is better than a +# script whose presence depends on how the package was installed. +[project.scripts] +hjcfit = "HJCFIT.likelihood.cli:main" + [project.urls] Homepage = "https://github.com/DCPROGS/HJCFIT" Documentation = "http://dcprogs.github.io/HJCFIT" diff --git a/tests/test_fitspec.py b/tests/test_fitspec.py new file mode 100644 index 0000000..ed97006 --- /dev/null +++ b/tests/test_fitspec.py @@ -0,0 +1,298 @@ +"""The fit specification: parsing, validation and round-tripping. + +Pure Python. Nothing here imports dcio, scalcs or even the likelihood, because +fitspec does not either -- a specification describes a fit and executes +nothing, and this file is the check that it stays that way. + +test_runner_cli.py is the other half: carrying one out. +""" + +import subprocess +import sys + +import pytest + +from HJCFIT.likelihood.fitspec import ( + TEMPLATE, DataSpec, FitSpec, MechanismSpec, SearchSpec, SpecError, + load_toml_bytes) + + +MINIMAL = { + 'data': [{'record': 'CH82', 'conc': 1e-7, 'tres': 1e-4, 'tcrit': 4e-3}], + 'mechanism': {'sample': 'CH82'}, +} + + +# -------------------------------------------------------------------------- +# the module's one structural promise +# -------------------------------------------------------------------------- + +def test_fitspec_imports_nothing_but_the_standard_library(): + """A specification must be readable without the [fitting] extra. + + The whole point of separating the description from the run is that + anything can read one -- a batch runner, a user interface, a script that + only wants to list what a directory of specifications contains. A stray + import of scalcs here would make reading a file cost a mechanism library. + """ + code = ( + "import sys;" + "import HJCFIT.likelihood.fitspec;" + "bad = [m for m in ('scalcs', 'dcio', 'matplotlib', 'scipy')" + " if m in sys.modules];" + "print(bad)" + ) + out = subprocess.run([sys.executable, '-c', code], check=True, + capture_output=True, text=True) + assert out.stdout.strip() == '[]', out.stdout + + +# -------------------------------------------------------------------------- +# DataSpec +# -------------------------------------------------------------------------- + +class TestDataSpec: + + def test_vectors_default_follows_tcrit(self): + """CHS vectors are defined between groups, so no groups means no CHS.""" + assert DataSpec.from_dict( + {'record': 'r', 'conc': 1e-7, 'tres': 1e-4, + 'tcrit': 4e-3}).vectors == 'chs' + assert DataSpec.from_dict( + {'record': 'r', 'conc': 1e-7, 'tres': 1e-4}).vectors == 'equilibrium' + + def test_chs_without_groups_is_refused(self): + with pytest.raises(SpecError, match='CHS vectors are defined between'): + DataSpec.from_dict({'record': 'r', 'conc': 1e-7, 'tres': 1e-4, + 'vectors': 'chs'}) + + def test_tcrit_below_tres_is_refused(self): + with pytest.raises(SpecError, match='divides nothing'): + DataSpec.from_dict({'record': 'r', 'conc': 1e-7, 'tres': 1e-4, + 'tcrit': 5e-5}) + + def test_a_negative_tcrit_is_not_a_flag_here(self): + """Elsewhere in the stack the sign selects the vectors. Not here.""" + with pytest.raises(SpecError, match='sign is not a flag'): + DataSpec.from_dict({'record': 'r', 'conc': 1e-7, 'tres': 1e-4, + 'tcrit': -4e-3, 'vectors': 'equilibrium'}) + + def test_zero_concentration_is_allowed(self): + """A mechanism with no agonist is an experiment, not a mistake.""" + assert DataSpec.from_dict( + {'record': 'r', 'conc': 0.0, 'tres': 1e-4}).conc == 0.0 + + def test_a_missing_required_field_names_itself(self): + with pytest.raises(SpecError, match='tres is required'): + DataSpec.from_dict({'record': 'r', 'conc': 1e-7}) + + @pytest.mark.parametrize('bad', ['vector', 'tcritical', 'resolution']) + def test_an_unknown_key_is_refused_with_the_known_ones(self, bad): + """A misspelled key is worse than a missing one: the fit would run, + silently, without the setting the file says it has.""" + d = {'record': 'r', 'conc': 1e-7, 'tres': 1e-4, bad: 1.0} + with pytest.raises(SpecError) as error: + DataSpec.from_dict(d) + assert bad in str(error.value) and 'known keys are' in str(error.value) + + def test_a_string_where_a_number_belongs_is_refused(self): + """The failure PyYAML would hand us silently. See the module docstring: + yaml.safe_load reads 1e8 as the string '1e8'.""" + with pytest.raises(SpecError, match='expected a number'): + DataSpec.from_dict({'record': 'r', 'conc': '1e-7', 'tres': 1e-4}) + + def test_a_boolean_is_not_a_number(self): + with pytest.raises(SpecError, match='expected a number'): + DataSpec.from_dict({'record': 'r', 'conc': True, 'tres': 1e-4}) + + def test_str_is_readable(self): + text = str(DataSpec.from_dict(MINIMAL['data'][0])) + assert '100 nM' in text and '100 us' in text and '4 ms' in text + + +# -------------------------------------------------------------------------- +# MechanismSpec +# -------------------------------------------------------------------------- + +class TestMechanismSpec: + + def test_exactly_one_source(self): + with pytest.raises(SpecError, match='exactly one of sample and'): + MechanismSpec.from_dict({}) + with pytest.raises(SpecError, match='exactly one of sample and'): + MechanismSpec.from_dict({'sample': 'CH82', 'mec_file': 'a.mec'}) + + def test_mec_number_without_a_file_means_nothing(self): + with pytest.raises(SpecError, match='means nothing without mec_file'): + MechanismSpec.from_dict({'sample': 'CH82', 'mec_number': 2}) + + def test_a_rate_cannot_be_both_a_guess_and_fixed(self): + with pytest.raises(SpecError, match='both rates and fixed'): + MechanismSpec.from_dict({'sample': 'CH82', + 'rates': {'beta1': 15.0}, + 'fixed': {'beta1': 15.0}}) + + def test_a_rate_cannot_be_both_fixed_and_from_the_ec50(self): + with pytest.raises(SpecError, match='cannot be both'): + MechanismSpec.from_dict({ + 'sample': 'CH82', 'fixed': {'k(+2)': 1e8}, + 'ec50': {'rate': 'k(+2)', 'value': 3.3e-6}}) + + def test_limits_must_be_a_pair_in_order(self): + with pytest.raises(SpecError, match=r'expected \[lower, upper\]'): + MechanismSpec.from_dict({'sample': 'CH82', + 'limits': {'beta1': [1.0]}}) + with pytest.raises(SpecError, match='must be below'): + MechanismSpec.from_dict({'sample': 'CH82', + 'limits': {'beta1': [10.0, 1.0]}}) + + def test_ec50_needs_both_fields(self): + with pytest.raises(SpecError, match='value is required'): + MechanismSpec.from_dict({'sample': 'CH82', + 'ec50': {'rate': 'k(+2)'}}) + + def test_rate_names_survive_a_round_trip(self): + """Rate constants are called things like 2k(+1), which TOML will not + take as a bare key. The writer has to quote them.""" + names = {'k(-1)': 1e3, '2k(+1)': 1e7, "k*(+2)": 5e8, 'beta1': 15.0} + spec = FitSpec.from_dict(dict(MINIMAL, mechanism={ + 'sample': 'CH82', 'rates': names})) + back = FitSpec.from_dict(load_toml_bytes(spec.to_toml().encode())) + assert back.mechanism.rates == names + + +# -------------------------------------------------------------------------- +# SearchSpec +# -------------------------------------------------------------------------- + +class TestSearchSpec: + + def test_defaults_are_hjcfits_own(self): + s = SearchSpec.from_dict({}) + assert s.method == 'simplex' and s.log_params is True + + def test_an_unknown_method_is_refused(self): + with pytest.raises(SpecError, match='must be one of simplex or scipy'): + SearchSpec.from_dict({'method': 'nelder-mead'}) + + def test_restarts_belong_to_the_scipy_search_alone(self): + """simplex_hjc has its own restart rule, capped by nresmax. Accepting + restarts here would silently do nothing.""" + with pytest.raises(SpecError, match='only the scipy search'): + SearchSpec.from_dict({'method': 'simplex', 'restarts': 8}) + assert SearchSpec.from_dict( + {'method': 'scipy', 'restarts': 8}).restarts == 8 + + def test_a_zero_budget_is_refused(self): + with pytest.raises(SpecError, match='at least 1'): + SearchSpec.from_dict({'maxfev': 0}) + + +# -------------------------------------------------------------------------- +# FitSpec +# -------------------------------------------------------------------------- + +class TestFitSpec: + + def test_a_fit_needs_a_record(self): + with pytest.raises(SpecError, match='at least one'): + FitSpec.from_dict({'mechanism': {'sample': 'CH82'}}) + + def test_a_single_data_table_is_taken_as_one_record(self): + """[data] rather than [[data]] is the easiest TOML mistake to make and + means exactly one thing, so it is read rather than refused.""" + spec = FitSpec.from_dict(dict(MINIMAL, data=MINIMAL['data'][0])) + assert len(spec.data) == 1 + + def test_several_records_keep_their_order(self): + spec = FitSpec.from_dict({ + 'data': [{'record': 'a', 'conc': 1e-8, 'tres': 1e-5}, + {'record': 'b', 'conc': 1e-7, 'tres': 1e-5}, + {'record': 'c', 'conc': 1e-6, 'tres': 1e-5}], + 'mechanism': {'sample': 'CH82'}}) + assert [d.record for d in spec.data] == ['a', 'b', 'c'] + + def test_the_index_of_a_bad_record_is_in_the_message(self): + with pytest.raises(SpecError, match=r'data\[1\]'): + FitSpec.from_dict({ + 'data': [{'record': 'a', 'conc': 1e-8, 'tres': 1e-5}, + {'record': 'b', 'conc': 1e-8, 'tres': -1.0}], + 'mechanism': {'sample': 'CH82'}}) + + def test_round_trip_through_toml(self): + spec = FitSpec.from_dict(dict( + MINIMAL, title='a fit', + mechanism={'sample': 'CH82', 'rates': {'beta1': 15.0}, + 'fixed': {'k(-1)': 2000.0}, + 'limits': {'beta1': [0.1, 1e5]}, + 'mr': 'k(-1)', 'nfree': 7}, + search={'method': 'scipy', 'restarts': 4, 'log_params': False})) + assert FitSpec.from_dict( + load_toml_bytes(spec.to_toml().encode())) == spec + + def test_to_toml_writes_scalars_before_sub_tables(self): + """TOML reads every key after a [header] as belonging to it, so a + scalar written below a sub-table would silently join it. This is a + correctness property of the writer, not a matter of taste.""" + text = FitSpec.from_dict(dict(MINIMAL, mechanism={ + 'sample': 'CH82', 'nfree': 8, + 'rates': {'beta1': 15.0}})).to_toml() + lines = [l.strip() for l in text.splitlines()] + assert lines.index('sample = "CH82"') < lines.index('[mechanism.rates]') + assert lines.index('nfree = 8') < lines.index('[mechanism.rates]') + + def test_from_toml_names_the_file_in_its_error(self, tmp_path): + path = tmp_path / 'broken.toml' + path.write_text('[[data]]\nrecord = "r"\nconc = 1e-7\n', + encoding='utf-8') + with pytest.raises(SpecError) as error: + FitSpec.from_toml(str(path)) + assert 'broken.toml' in str(error.value) + assert 'tres is required' in str(error.value) + + def test_unparseable_toml_says_so(self, tmp_path): + path = tmp_path / 'nonsense.toml' + path.write_text('this is not toml = = =\n', encoding='utf-8') + with pytest.raises(SpecError, match='not valid TOML'): + FitSpec.from_toml(str(path)) + + def test_str_mentions_every_constraint(self): + text = str(FitSpec.from_dict(dict(MINIMAL, mechanism={ + 'sample': 'CH82', 'fixed': {'k(-1)': 2000.0}, 'mr': 'beta1', + 'limits': {'beta2': [1.0, 1e5]}, + 'ec50': {'rate': 'k(+2)', 'value': 3.3e-6}}))) + assert 'k(-1) = 2000' in text + assert 'microscopic reversibility' in text + assert 'EC50 of 3.3 uM' in text + assert 'limits on beta2' in text + + +# -------------------------------------------------------------------------- +# the template, and the numbers TOML gets right +# -------------------------------------------------------------------------- + +class TestTemplate: + + def test_the_template_is_a_valid_specification(self): + """It is what `hjcfit template` hands a new user, so it has to parse + and to describe a fit that runs.""" + spec = FitSpec.from_dict(load_toml_bytes(TEMPLATE.encode())) + assert spec.data[0].record == 'CH82' + assert spec.mechanism.sample == 'CH82' + assert spec.search.method == 'simplex' + + def test_the_template_round_trips(self): + spec = FitSpec.from_dict(load_toml_bytes(TEMPLATE.encode())) + assert FitSpec.from_dict( + load_toml_bytes(spec.to_toml().encode())) == spec + + @pytest.mark.parametrize('text', ['1e8', '1.0e8', '1e+8', '1.0e+8', + '100e-9', '1E8']) + def test_every_way_of_writing_a_rate_constant_is_a_number(self, text): + """The reason this format is TOML. PyYAML implements YAML 1.1, whose + resolver wants a signed exponent: yaml.safe_load reads three of these + six as strings, and a rate constant arriving as a string is the worst + failure a specification could have. + """ + value = load_toml_bytes('x = {0}\n'.format(text).encode())['x'] + assert isinstance(value, float) diff --git a/tests/test_runner_cli.py b/tests/test_runner_cli.py new file mode 100644 index 0000000..05c7d25 --- /dev/null +++ b/tests/test_runner_cli.py @@ -0,0 +1,432 @@ +"""Carrying out a specification: the runner, the hjcfit command, the notebook. + +Needs the [fitting] extra -- dcio for the records, scalcs for the mechanism -- +so it runs in the dcio-integration job rather than the build matrix, and skips +rather than fails where those are absent. + +The fit here is the CH82 sample record that ships with HJCFIT, which is the +same fit documented in documentation/source/api/python/fitting.rst: log10 L +goes from 2286.97 at the guess to 2289.13. +""" + +import json +import os +import sys +from pathlib import Path + +import numpy as np +import pytest + +dcio = pytest.importorskip('dcio', reason='needs the [fitting] extra') +scalcs = pytest.importorskip('scalcs', reason='needs the [fitting] extra') + +import HJCFIT +from HJCFIT.likelihood import cli +from HJCFIT.likelihood.fitspec import FitSpec, load_toml_bytes +from HJCFIT.likelihood.runner import (RunnerError, build_mechanism, + load_records, provenance, + result_as_dict, run, write_result) + +HERE = Path(__file__).resolve().parent +EXAMPLES = HERE.parent / 'examples' + +#: The specification every test below starts from, as a dictionary so that +#: individual tests can vary one field of it. +CH82 = { + 'title': 'CH82 at 100 nM', + 'data': [{'record': 'CH82', 'conc': 1e-7, 'tres': 1e-4, 'tcrit': 4e-3}], + 'mechanism': {'sample': 'CH82', 'rates': {'beta1': 15.0, 'beta2': 15000.0}}, +} + + +def spec_for(**changes): + return FitSpec.from_dict(dict(CH82, **changes)) + + +def write(tmp_path, spec, name='spec.toml'): + path = tmp_path / name + spec.write_toml(str(path)) + return str(path) + + +# -------------------------------------------------------------------------- +# records +# -------------------------------------------------------------------------- + +class TestLoadRecords: + + def test_a_sample_name_resolves(self): + record = load_records(spec_for())[0] + assert record.n_openings == 836 + assert len(record.groups) == 572 + assert record.n_raw == 4312 + + def test_groups_are_identical_to_read_idealized_bursts(self): + """Two callers of the same pairing of dcio functions, so the check is + identity rather than similarity. + + read_idealized_bursts carries the argument for segmenting the record's + *periods* rather than its resolved intervals -- imposing a dead time + emits a fresh open interval at every change of fitted amplitude, so a + record with sub-conductance levels does not alternate open and shut. + The runner has to do the same thing, and this is what says it does. + """ + mine = load_records(spec_for())[0].groups + theirs = HJCFIT.read_idealized_bursts('CH82', tau=1e-4, tcrit=4e-3) + assert len(mine) == len(theirs) + for a, b in zip(mine, theirs): + assert np.array_equal(np.asarray(a, dtype=float), + np.asarray(b, dtype=float)) + + def test_no_tcrit_gives_one_group_starting_and_ending_open(self): + record = load_records(spec_for(data=[ + {'record': 'CH82', 'conc': 1e-7, 'tres': 1e-4}]))[0] + assert len(record.groups) == 1 + assert len(record.groups[0]) % 2 == 1 # the likelihood requires it + assert record.tcrit is None # equilibrium vectors + assert record.n_openings == 836 # every opening, as in bursts + + def test_equilibrium_vectors_on_groups_still_gives_groups(self): + """The cluster case: divided at a critical time, but started and ended + with equilibrium vectors rather than CHS ones. tcrit does two jobs and + only one of them is being turned off.""" + record = load_records(spec_for(data=[ + {'record': 'CH82', 'conc': 1e-7, 'tres': 1e-4, 'tcrit': 4e-3, + 'vectors': 'equilibrium'}]))[0] + assert len(record.groups) == 572 + assert record.tcrit is None + + def test_a_record_that_is_not_there_says_which_samples_are(self): + with pytest.raises(RunnerError, match='CH82, CO and CCO'): + load_records(spec_for(data=[ + {'record': 'not-a-file.scn', 'conc': 1e-7, 'tres': 1e-4}])) + + def test_several_records_come_back_in_order(self): + records = load_records(spec_for(data=[ + {'record': 'CH82', 'conc': 1e-8, 'tres': 1e-4, 'tcrit': 4e-3}, + {'record': 'CH82', 'conc': 1e-7, 'tres': 1e-4, 'tcrit': 4e-3}])) + assert [r.conc for r in records] == [1e-8, 1e-7] + + +# -------------------------------------------------------------------------- +# the mechanism +# -------------------------------------------------------------------------- + +class TestBuildMechanism: + + def test_the_guess_is_applied(self): + mec = build_mechanism(spec_for()) + rates = {r.name: float(np.ravel(r.rateconstants)[0]) for r in mec.Rates} + assert rates['beta1'] == 15.0 and rates['beta2'] == 15000.0 + + def test_eight_free_parameters_for_ch82(self): + assert build_mechanism(spec_for()).get_free_parameter_names() == [ + 'beta1', 'beta2', 'alpha1', 'alpha2', 'k(-1)', '2k(-2)', + '2k(+1)', 'k(+2)'] + + def test_a_fixed_rate_leaves_the_free_list(self): + mec = build_mechanism(spec_for(mechanism={ + 'sample': 'CH82', 'fixed': {'k(-1)': 2000.0}})) + assert 'k(-1)' not in mec.get_free_parameter_names() + assert len(mec.get_free_parameter_names()) == 7 + + def test_limits_reach_the_rate(self): + mec = build_mechanism(spec_for(mechanism={ + 'sample': 'CH82', 'limits': {'beta1': [0.5, 50.0]}})) + beta1 = next(r for r in mec.Rates if r.name == 'beta1') + assert list(beta1.limits[0]) == [0.5, 50.0] + + def test_an_unknown_rate_name_lists_the_names_there_are(self): + """The interface problem this module exists to fix: SCALCS addresses + rates by index into mec.Rates, and an index is what gets got wrong + silently.""" + with pytest.raises(RunnerError) as error: + build_mechanism(spec_for(mechanism={ + 'sample': 'CH82', 'rates': {'beta9': 1.0}})) + assert "no rate called 'beta9'" in str(error.value) + assert 'beta1' in str(error.value) and '2k(+1)' in str(error.value) + + def test_an_unknown_sample_lists_the_samples_there_are(self): + with pytest.raises(RunnerError) as error: + build_mechanism(spec_for(mechanism={'sample': 'NotAMechanism'})) + assert 'CH82' in str(error.value) and 'GlyR_flip' in str(error.value) + + def test_nfree_is_asserted_before_the_search(self): + with pytest.raises(RunnerError, match='nfree says 7'): + build_mechanism(spec_for(mechanism={'sample': 'CH82', 'nfree': 7})) + build_mechanism(spec_for(mechanism={'sample': 'CH82', 'nfree': 8})) + + def test_fixing_everything_is_refused(self): + every = {r.name: 1.0 for r in build_mechanism(spec_for()).Rates} + with pytest.raises(RunnerError, match='nothing to fit'): + build_mechanism(spec_for(mechanism={'sample': 'CH82', + 'fixed': every})) + + def test_mr_constrains_a_rate_that_is_in_the_cycle(self): + mec = build_mechanism(spec_for(mechanism={ + 'sample': 'CH82', 'mr': 'beta1'})) + assert 'beta1' not in mec.get_free_parameter_names() + assert len(mec.get_free_parameter_names()) == 8 # 2k*(-2) freed + + def test_mr_on_a_rate_outside_the_cycle_is_refused(self): + """SCALCS warns on stderr and carries on, and what it does then is + worse than nothing: the rate leaves the free-parameter list while the + cycle's constraint stays where it was, so the rate is silently frozen + at its initial guess and the fit has one parameter fewer than the + specification says. On CH82, mr = "k(-1)" turned eight free parameters + into seven with k(-1) stuck at the guess. + """ + with pytest.raises(RunnerError) as error: + build_mechanism(spec_for(mechanism={'sample': 'CH82', + 'mr': 'k(-1)'})) + assert 'not in cycle 0' in str(error.value) + assert 'A2R*' in str(error.value) # the states it could have used + + def test_mr_cycle_out_of_range_is_refused(self): + with pytest.raises(RunnerError, match='not a cycle of this mechanism'): + build_mechanism(spec_for(mechanism={ + 'sample': 'CH82', 'mr': 'beta1', 'mr_cycle': 3})) + + +# -------------------------------------------------------------------------- +# the fit +# -------------------------------------------------------------------------- + +class TestRun: + + @pytest.fixture(scope='class') + def outcome(self): + return run(spec_for()) + + def test_it_reaches_the_documented_maximum(self, outcome): + """The fit in fitting.rst: 2286.97 at the guess, 2289.13 at the end.""" + assert outcome.result.log10_likelihood > 2289.0 + assert outcome.result.success + assert outcome.result.nfailures == 0 + + def test_the_free_values_are_rates_not_logarithms(self, outcome): + """log_params searches the logarithms; a result reporting them would + be a units error that looks like an implausible mechanism.""" + values = dict(zip(outcome.result.free_names, outcome.result.free_values)) + assert 1e4 < values['k(+2)'] < 1e10 + assert 1.0 < values['beta1'] < 100.0 + + def test_the_mechanism_carries_the_fitted_rates(self, outcome): + """Predicted distributions are drawn from this, so it must not still + hold the guess.""" + rates = {r.name: float(np.ravel(r.rateconstants)[0]) + for r in outcome.mec.Rates} + assert rates['beta1'] != 15.0 + assert rates['beta1'] == pytest.approx( + dict(zip(outcome.result.free_names, + outcome.result.free_values))['beta1']) + + def test_nothing_ends_on_a_limit_here(self, outcome): + assert outcome.limited == [] + + def test_a_rate_driven_onto_a_limit_is_reported(self): + """A rate outside its limits is reset before the likelihood sees it, + and nothing announces it, so the fitted value reads like an estimate + when it is a statement that the model forbids where the likelihood + wanted to go.""" + outcome = run(spec_for( + mechanism={'sample': 'CH82', 'rates': {'beta1': 4.0}, + 'limits': {'beta1': [1.0, 5.0]}}, + search={'maxfev': 3000})) + assert [name for name, _, _, _ in outcome.limited] == ['beta1'] + name, value, which, limit = outcome.limited[0] + assert which == 'upper' and value == pytest.approx(5.0) + assert 'not estimates' in str(outcome) + + def test_the_scipy_search_reaches_the_same_maximum(self): + """Different paths, same maximum. Agreement on the height while + disagreeing on where it sits is evidence about the likelihood rather + than about either search.""" + outcome = run(spec_for(search={'method': 'scipy', 'restarts': 4})) + assert outcome.result.log10_likelihood > 2289.0 + + def test_str_says_what_was_fitted_as_well_as_what_came_out(self, outcome): + text = str(outcome) + assert '572 groups' in text and '836 openings' in text + assert 'log10(L)' in text and 'beta1' in text + + def test_ln_and_log10_differ_by_ln_10(self, outcome): + """Anything treating the likelihood as a statistical quantity needs + natural logarithms; using log10 inflates every standard deviation by + sqrt(ln 10) = 1.517, which looks like a bad fit rather than a units + error.""" + d = result_as_dict(outcome) + assert d['fit']['ln_likelihood'] == pytest.approx( + d['fit']['log10_likelihood'] * np.log(10.0)) + + +# -------------------------------------------------------------------------- +# the written result +# -------------------------------------------------------------------------- + +class TestWrittenResult: + + def test_it_carries_the_specification_that_produced_it(self, tmp_path): + outcome = run(spec_for(search={'maxfev': 200})) + path = tmp_path / 'result.json' + write_result(str(path), outcome) + d = json.loads(path.read_text(encoding='utf-8')) + assert FitSpec.from_dict(d['spec']) == outcome.spec + + def test_it_carries_the_versions_it_was_computed_with(self, tmp_path): + outcome = run(spec_for(search={'maxfev': 200})) + d = result_as_dict(outcome) + assert d['provenance']['versions']['hjcfit'] + assert d['provenance']['versions']['scalcs'] + assert d['provenance']['python'] == sys.version.split()[0] + + def test_it_does_not_carry_the_records(self, tmp_path): + """Thousands of floats that are already in the .scn file the + specification names.""" + outcome = run(spec_for(search={'maxfev': 200})) + path = tmp_path / 'result.json' + write_result(str(path), outcome) + assert path.stat().st_size < 8000 + d = json.loads(path.read_text(encoding='utf-8')) + assert d['records'][0]['groups'] == 572 + assert d['records'][0]['intervals'] == 1100 + + def test_provenance_has_no_hole_where_a_missing_package_would_be(self): + prov = provenance() + assert set(prov['versions']) >= {'hjcfit', 'dcio', 'scalcs', 'numpy'} + + +# -------------------------------------------------------------------------- +# the command +# -------------------------------------------------------------------------- + +class TestCommand: + + def test_template_then_check_then_fit(self, tmp_path, capsys, monkeypatch): + """The path a new user takes, end to end.""" + monkeypatch.chdir(tmp_path) + path = str(tmp_path / 'mine.toml') + + assert cli.main(['template', '-o', path]) == 0 + assert 'wrote' in capsys.readouterr().out + + assert cli.main(['check', path]) == 0 + out = capsys.readouterr().out + assert '572 groups' in out + assert '8 free:' in out + assert 'log10L at the guess: 2286.9746' in out + assert 'nothing was fitted' in out + + assert cli.main(['fit', path, '-o', 'out.json']) == 0 + out = capsys.readouterr().out + assert 'log10(L) = 2289.13' in out + d = json.loads((tmp_path / 'out.json').read_text(encoding='utf-8')) + assert d['fit']['log10_likelihood'] > 2289.0 + + def test_template_refuses_to_overwrite_without_force(self, tmp_path, + capsys): + path = str(tmp_path / 'mine.toml') + assert cli.main(['template', '-o', path]) == 0 + assert cli.main(['template', '-o', path]) == 1 + assert 'pass --force' in capsys.readouterr().err + assert cli.main(['template', '-o', path, '--force']) == 0 + + def test_show_normalises_a_specification(self, tmp_path, capsys): + path = write(tmp_path, spec_for()) + assert cli.main(['show', path]) == 0 + out = capsys.readouterr().out + assert FitSpec.from_dict(load_toml_bytes(out.encode())) == spec_for() + + def test_a_search_that_does_not_converge_exits_2(self, tmp_path, capsys): + """The number is still printed -- it is the exit status that says so, + because that is what a script reads.""" + path = write(tmp_path, spec_for(search={'maxfev': 50})) + assert cli.main(['fit', path, '-q']) == 2 + assert 'log10(L)' in capsys.readouterr().out + + def test_a_bad_specification_is_one_line_on_stderr(self, tmp_path, capsys): + path = tmp_path / 'bad.toml' + path.write_text('[[data]]\nrecord = "CH82"\nconc = 1e-7\n', + encoding='utf-8') + assert cli.main(['check', str(path)]) == 1 + err = capsys.readouterr().err + assert err.startswith('hjcfit: ') and err.count('\n') == 1 + assert 'tres is required' in err + + def test_a_missing_file_is_not_a_traceback(self, tmp_path, capsys): + assert cli.main(['check', str(tmp_path / 'nope.toml')]) == 1 + assert capsys.readouterr().err == ( + 'hjcfit: no such file: {0}\n'.format(tmp_path / 'nope.toml')) + + def test_an_impossible_rate_name_is_one_line_too(self, tmp_path, capsys): + path = write(tmp_path, spec_for( + mechanism={'sample': 'CH82', 'rates': {'beta9': 1.0}})) + assert cli.main(['check', path]) == 1 + err = capsys.readouterr().err + assert err.count('\n') == 1 and 'beta9' in err + + def test_version_lists_everything_involved(self, capsys): + assert cli.main(['--version']) == 0 + out = capsys.readouterr().out + for package in ('hjcfit', 'dcio', 'scalcs', 'numpy', 'scipy'): + assert package in out + + def test_no_arguments_prints_help_rather_than_failing(self, capsys): + assert cli.main([]) == 0 + assert 'hjcfit' in capsys.readouterr().out + + def test_the_console_script_entry_point_resolves(self): + """pyproject declares hjcfit = HJCFIT.likelihood.cli:main. A declared + entry point nothing imports is untested and will one day be wrong.""" + from importlib.metadata import entry_points + scripts = {e.name: e for e in entry_points(group='console_scripts')} + assert 'hjcfit' in scripts + assert scripts['hjcfit'].load() is cli.main + + +# -------------------------------------------------------------------------- +# the shipped example, and the notebook over it +# -------------------------------------------------------------------------- + +@pytest.mark.skipif(not EXAMPLES.is_dir(), + reason='examples/ is not beside the tests') +class TestExamples: + + def test_the_example_specification_is_valid_and_its_nfree_holds(self): + """It states nfree = 8, which is an assertion about how SCALCS handles + constraints. If that changes, this is where it shows.""" + spec = FitSpec.from_toml(str(EXAMPLES / 'CH82.toml')) + assert spec.mechanism.nfree == 8 + assert len(build_mechanism(spec).get_free_parameter_names()) == 8 + + def test_the_notebook_template_runs(self, tmp_path, monkeypatch): + """A notebook is the de facto interface to this package, and an + example that has stopped working is worse than none. This is the one + thing a desktop interface could not have: it is testable. + """ + nbformat = pytest.importorskip('nbformat') + nbclient = pytest.importorskip('nbclient') + pytest.importorskip('matplotlib') + + monkeypatch.setenv('MPLBACKEND', 'Agg') + # Copied, because the notebook writes result.json beside itself and a + # test must not leave anything in the source tree. + import shutil + for name in ('fit_template.ipynb', 'CH82.toml'): + shutil.copy(EXAMPLES / name, tmp_path / name) + + nb = nbformat.read(str(tmp_path / 'fit_template.ipynb'), as_version=4) + nbclient.NotebookClient( + nb, timeout=900, kernel_name='python3', + resources={'metadata': {'path': str(tmp_path)}}).execute() + + text = "\n".join( + output.get('text', '') + for cell in nb.cells for output in cell.get('outputs', []) + if output.output_type == 'stream') + assert 'log10(L) = 2289.13' in text + # The two checks the notebook makes on its own curves. A predicted + # curve over an observed histogram is the one plot where a wrong curve + # looks like a finding about the fit rather than a bug. + assert 'pdf integrates to 1.00000' in text + assert os.path.exists(tmp_path / 'result.json')