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