A fit from a file

Fitting a mechanism 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
hjcfit sample my-fit.toml -o post  sample the posterior of the rates

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

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.

$ 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 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. 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. against_limits() finds them and they are printed apart from the fitted values.

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^-15 to 10⁹ for a concentration-dependent rate and 10^-15 to 10⁶ for the rest – so this can happen in a specification that sets no limits at all.

The specification

class HJCFIT.likelihood.fitspec.FitSpec(data: tuple, mechanism: MechanismSpec = <factory>, search: SearchSpec = <factory>, title: str = '', likelihood: LikelihoodSpec = <factory>, mcmc: MCMCSpec = <factory>)

A whole fit, as data.

Parameters:

• data – One DataSpec per record. More than one means several records contributing to a single likelihood: the rate constants are shared and the concentrations are not.

• mechanism – A MechanismSpec.

• search – A SearchSpec.

• likelihood – A LikelihoodSpec: root-finding settings, for fitting and sampling alike.

• mcmc – A MCMCSpec, read by hjcfit sample and ignored by hjcfit fit.

• title – One line for the person reading the result.

Read one with from_toml(), write one with to_toml(), and run one with HJCFIT.likelihood.runner.run().

as_dict()

A plain dictionary, ready for to_toml() or for JSON.

Defaults are left out, so a specification written back is the shortest file that means the same thing.

classmethod from_dict(d, where='spec')

Build one from a parsed document.

Raises:

SpecError – naming the field, on anything unusable.

classmethod from_toml(path)

Read a specification from a TOML file.

Parameters:

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.

to_toml()

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.

validate(where='spec')

Raise SpecError unless the whole thing hangs together.

write_toml(path)

Write to_toml() to path, as UTF-8.

class HJCFIT.likelihood.fitspec.DataSpec(record: str, conc: float, tres: float, tcrit: float = None, vectors: str = 'chs')

One record, and how it is presented to the likelihood.

Parameters:

• record – An *.scn file, or the name of a sample record shipped with HJCFIT ("CH82", "CO", "CCO").

• conc – Agonist concentration [M].

• tres – Dead time to impose [s].

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

• 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 – HJCFIT.read_idealized_bursts() documents that, and 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 validate() can refuse.

classmethod from_dict(d, where='data')

Build one from a parsed table, checking as it goes.

validate(where='data')

Raise SpecError unless this record can be fitted.

class HJCFIT.likelihood.fitspec.MechanismSpec(sample: str = None, mec_file: str = None, mec_number: int = None, rates: dict = <factory>, fixed: dict = <factory>, limits: dict = <factory>, mr: str = None, mr_cycle: int = 0, ec50: dict = None, nfree: int = None)

Which mechanism, with which constraints and which initial guess.

Parameters:

• sample – Name of a factory in scalcs.samples.samples, for example "CH82" or "load_AChR_diamond_independent_binding". Exactly one of this and mec_file.

• mec_file – A DCprogs *.mec file, read with scalcs.scalcsio.mec_load().

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

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

• fixed – Rate constants held at a stated value throughout the fit.

• 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 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 – runner does that comparison and says so.

• mr – Name of the rate constant computed from microscopic reversibility rather than fitted.

• mr_cycle – Which cycle, when the mechanism has more than one. Cycles have no names in SCALCS, so this one is an index.

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

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

classmethod from_dict(d, where='mechanism')

Build one from a parsed table, checking as it goes.

validate(where='mechanism')

Raise SpecError unless this mechanism can be built.

class HJCFIT.likelihood.fitspec.SearchSpec(method: str = 'simplex', log_params: bool = True, restarts: int = 0, maxfev: int = 20000, maxiter: int = 20000, xatol: float = 0.0001, fatol: float = 0.0001, options: dict = <factory>)

How to search, and how hard.

Parameters:

• method – "simplex" is 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.

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

• 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₁₀ units in total, resolved five of them onto a real maximum, and never lost any.

• maxfev – Evaluation budget; both searches.

• maxiter – "scipy" only.

• xatol – "scipy" only.

• fatol – "scipy" only.

• options – Passed straight to whichever search is chosen, for deliberate departures from its own defaults.

class HJCFIT.likelihood.fitspec.LikelihoodSpec(nmax: int = None, xtol: float = None, rtol: float = None, itermax: int = None, lower_bound: float = None, upper_bound: float = None)

Root-finding settings for every record’s likelihood.

The names are those of SOLVER_OPTIONS; anything left out keeps Log10Likelihood’s default. They are not cosmetic. On the three AChR records of Epstein et al. (2016), nmax = 2 with tolerances of 1e-12 – that paper’s settings – move the natural log-likelihood by 0.38 against the defaults at the same rates, so a published value is reproduced only with the settings it was computed with.

Applies to hjcfit fit and hjcfit sample alike.

solver()

The settings given, as the solver argument of the fitter.

class HJCFIT.likelihood.fitspec.MCMCSpec(start: str = 'fit', sampler: str = 'adaptive', n: int = 20000, burnin: int = 5000, pilot_n: int = 3000, covariance: str = 'hessian', mixture: str = 'sum', log_space: bool = None, prior: str = 'uniform', chains: int = 1, seed: int = 1)

How to sample the posterior, for hjcfit sample.

The method is that of Epstein, Calderhead, Girolami & Sivilotti (2016); see Bayesian inference.

Parameters:

• start – "fit" (the default) runs the [search] first and starts every chain at the maximum-likelihood estimate, with the Hessian there as the first proposal covariance. "guess" is the paper’s own scheme: a Metropolis-within-Gibbs pilot from the initial guess finds the mode instead. With the paper’s flat prior the two points coincide; the fit is usually the cheaper way there.

• sampler – "adaptive", adaptive Metropolis, or "mwg", Metropolis-within-Gibbs.

• n – Iterations per chain, burn-in included.

• burnin – Iterations during which step sizes are tuned; not kept.

• pilot_n – start = "guess" only: iterations of the pilot, half of them burn-in.

• covariance – "hessian": shape the adaptive sampler’s first proposals with the inverse Hessian at the start, at the optimal scale 2.38/sqrt(k). Where the Hessian cannot be used – a point that is not a maximum, a direction no record determines – the run says so and falls back to "identity", a small isotropic step, which is the paper’s.

• mixture – "sum" (the paper’s) or "choice"; see adaptive_sample().

• log_space – Walk the logarithms of the rates. Unset, each sampler keeps its own default: logarithms for "mwg", the rates themselves for "adaptive".

• prior – "uniform" (the paper’s) or "loguniform", in both cases between each free rate’s limits.

• chains – Independent chains, run in parallel processes. Chains after the first start from draws of the Gaussian approximation at the start point, when there is one, so that their agreement means something.

• seed – Seeds chain i with seed + i.

exception HJCFIT.likelihood.fitspec.SpecError

A specification that cannot mean anything.

Raised by the from_dict constructors and by validate. Every message names the field it is about.

HJCFIT.likelihood.fitspec.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.

HJCFIT.likelihood.fitspec.TEMPLATE

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.

Carrying it out

HJCFIT.likelihood.runner.run(spec, verbose=False, store_path=False, x0=None)

Read the records, build the mechanism, fit, and report.

Parameters:

• spec – A FitSpec.

• verbose – Print the records and the starting likelihood first.

• store_path – Keep the best point of every iteration.

• x0 – Starting point, in the space being searched. Defaults to the mechanism’s own guess, which is what the specification says.

Returns:

An Outcome.

class HJCFIT.likelihood.runner.Outcome(spec: object, result: object, records: list, mec: object = None, limited: list = <factory>, provenance: dict = <factory>)

Everything one run produced.

Parameters:

• spec – The specification it was run from.

• result – The FitResult.

• records – The records as they were fitted.

• mec – The mechanism, carrying the fitted rate constants – which is what predicted distributions have to be drawn from.

• limited – What against_limits() found. Empty is the good case.

• provenance – What provenance() recorded.

HJCFIT.likelihood.runner.load_records(spec)

Read every record a specification names.

Parameters:

spec – A FitSpec.

Returns:

A list of 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 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. 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.

HJCFIT.likelihood.runner.build_mechanism(spec)

Build the mechanism a specification names, with its constraints.

Parameters:

spec – A 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.

HJCFIT.likelihood.runner.against_limits(mec, tolerance=1e-06)

Free rates that ended on one of their limits.

Parameters:

• mec – A mechanism, after a fit.

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

HJCFIT.likelihood.runner.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.

HJCFIT.likelihood.runner.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.

HJCFIT.likelihood.runner.write_result(path, outcome)

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

HJCFIT.likelihood.runner.sample(spec, processes=None, verbose=False)

Sample the posterior a specification describes ([mcmc]).

Parameters:

• spec – A FitSpec.

• processes – Worker processes for the chains. None uses one per chain, up to the number of CPUs; 1 runs them one after another in this process.

• verbose – Print the stages as they happen.

Returns:

A SampleOutcome.

The C++ likelihood already spreads each evaluation over the cores with OpenMP, but past two or three threads it gains little. Measured on the AChR records of Epstein et al. (2016): 57 ms on one thread, 36 ms on two, 29 ms on four. So parallel chains are run as separate processes, each limited to its share of the logical CPUs.

What that buys depends on the machine. On a 4-core, 8-thread laptop, four CH82 chains ran 1.3-1.6 times faster in four processes than one after another in one. Fewer threads per process (one each) was slower there, because hyperthreading helps. More physical cores, or a likelihood that costs more per evaluation than CH82’s 1-2 ms, gain more.

Worker processes are spawned. As with any spawned process on Windows, a script that calls this with processes > 1 must do so under if __name__ == '__main__':. The hjcfit command and notebooks need nothing.

class HJCFIT.likelihood.runner.SampleOutcome(spec: object, records: list, names: tuple, start: ndarray, start_from: str, chains: list, summary: list, fit: object = None, pilot: object = None, approximation: object = None, notes: list = <factory>, seconds: float = 0.0, provenance: dict = <factory>)

Everything one sampling run produced.

Parameters:

• spec – The specification.

• records – The records as they were sampled against.

• start – Where the first chain started, as rates.

• start_from – 'fit' or 'pilot'.

• fit – The Outcome of the fit, for start = "fit".

• pilot – The pilot Chain, for start = "guess".

• approximation – The GaussianApproximation at the start, when the Hessian could be used.

• chains – One Chain each.

• summary – Per free rate: mean, sd, quantiles, ESS summed over chains and split-R-hat across them.

• notes – Anything the run had to decide that the reader should know, such as falling back from the Hessian.

• seconds – Wall-clock time of the whole run.

HJCFIT.likelihood.runner.samples_as_dict(outcome)

A sampling outcome as JSON-ready data: the specification, the start, the summary and the provenance. The chains themselves go in .npz files beside it (write_samples()).

HJCFIT.likelihood.runner.write_samples(prefix, outcome)

Write <prefix>.json (samples_as_dict()) and one <prefix>_chain<i>.npz per chain (save()).

Returns:

The paths written.

exception HJCFIT.likelihood.runner.RunnerError

A specification that cannot be carried out.

Distinct from 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.

The command

HJCFIT.likelihood.cli.main(argv=None)

The hjcfit command.

Parameters:

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 or from sample when the chains disagree (split-R-hat above 1.1).

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, or from sample when the chains disagree (split-R-hat above 1.1). The estimates are still printed in that case; it is the status that says so, because that is what a script reads.