Fitting a mechanism

Log10Likelihood takes a Q matrix and returns a number. Between that and a person sits this module: a record in the form the likelihood wants, a search over the free rate constants of a mechanism, and a result worth printing.

It needs a mechanism, and HJCFIT does not define one — mechanisms, rate 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, A fit from a file 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:

from HJCFIT.likelihood.fitting import HJCFitter, Record

That is not a stylistic choice. Exposing it from the package __init__ made HJCFIT.likelihood import itself while partially initialised, and every CI job failed on it.

Nothing here imports scalcs. The fitter is duck-typed on a mechanism: anything offering theta(), theta_unsqueeze(), Rates, kA, set_eff() and Q will do, which in practice means scalcs.mechanism.Mechanism. That is what keeps pip install hjcfit wanting nothing but numpy and scipy.

A whole fit

import HJCFIT
from HJCFIT.likelihood.fitting import HJCFitter, Record
from scalcs.samples import samples

bursts = HJCFIT.read_idealized_bursts("CH82", tau=1e-4, tcrit=4e-3)
record = Record(conc=100e-9, tres=1e-4, tcrit=4e-3,
                groups=tuple(tuple(b) for b in bursts))

mechanism = samples.CH82()
mechanism.set_eff("c", 100e-9)

result = HJCFitter(mechanism, [record]).fit()
print(result)

log10(L) = 2289.1326   1313 evaluations in 1.1 s
  beta1            11.31
  beta2        1.328e+04
  alpha1            3668
  alpha2           424.4
  k(-1)             1605
  2k(-2)            4410
  2k(+1)       1.471e+04
  k(+2)        5.288e+08

Several concentrations at once

The likelihood takes one Q matrix per call, so a fit to several concentrations builds one instance per record and adds the logarithms. Pass them all and HJCFitter does that:

fitter = HJCFitter(mechanism, [rec_10uM, rec_30uM, rec_100uM])

The mechanism is set to each record’s own concentration before that record is evaluated, so the rate constants are shared and the concentrations are not.

The record

class HJCFIT.likelihood.fitting.Record(conc: float, groups: tuple, tres: float, tcrit: float = None, n_raw: int = None, n_apparent: int = None)[source]

One idealised experiment, in the form the likelihood wants.

Parameters:

  • conc – Agonist concentration [M].

  • groups – Alternating open/shut intervals, each group starting and ending with an opening. One group per burst or cluster; a single group when the whole record is treated as coming from one channel.

  • tres – Dead time already imposed [s].

  • tcrit – What the likelihood is given as tcritical. A number asks for CHS vectors (Colquhoun, Hawkes & Srodzinski 1996); None asks for equilibrium vectors (Colquhoun & Hawkes 1982). Note that this is the magnitude: elsewhere in the stack a negative tcrit is a flag selecting equilibrium vectors, and segmentation takes abs.

  • n_raw – Intervals before the dead time was imposed, if known.

  • n_apparent – Intervals after it. Defaults to the number actually held.

as_lists()[source]

The form Log10Likelihood wants: a list of lists.

check()[source]

Raise if the groups are not what the likelihood requires.

Called by HJCFitter on construction, because an even group does not fail loudly – it silently asks the likelihood for a product of matrices that does not alternate.

Raises:

ValueError – on an empty record or an even-length group.

property n_intervals

Intervals across every group.

property n_openings

Openings across every group; each group starts and ends on one.

HJCFIT.likelihood.fitting.trim_to_openings(intervals, amplitudes)[source]

Trim an alternating record so it starts and ends with an opening.

The likelihood requires an odd number of intervals per group. A record divided into bursts already satisfies that; a whole record fitted as a single group generally does not.

Parameters:

  • intervals – Durations, alternating open and shut.

  • amplitudes – Matching amplitudes; zero means shut.

Returns:

A view of intervals beginning and ending on an opening.

Every group must have an odd number of intervals — it must start and end with an opening — because the likelihood is a product of matrices alternating \(A \rightarrow F\) and \(F \rightarrow A\). Burst segmentation in dcio already gives odd groups; a whole record fitted as one group generally does not, which is what trim_to_openings() is for.

An even group does not fail loudly. Record.check() is called on construction of the fitter for that reason.

The fitter

class HJCFIT.likelihood.fitting.HJCFitter(mec, records, log_params=True, store_path=False, store_evaluations=False, solver=None)[source]

Maximise the HJC likelihood of one or more records over a mechanism.

Parameters:

  • mec – A mechanism already carrying its constraints and its initial guess. Duck-typed: it needs theta(), theta_unsqueeze(), Rates, kA and set_eff(). scalcs.mechanism.Mechanism provides them.

  • records – A sequence of Record, fitted simultaneously. The likelihood takes one Q matrix per call, so each record is evaluated at its own concentration and the logarithms added.

  • log_params (bool) – Search the logarithms of the rate constants rather than the rates. This is HJCFIT’s own default and three to four times faster (Colquhoun, Hatton & Hawkes 2003, p. 702). It also cannot produce a negative rate. The fits of that paper’s Figures 2-5 and 12-13 were made over the rates themselves, which is what the resetting below is for.

  • store_path (bool) – Keep the best vertex at each iteration. Cheap.

  • store_evaluations (bool) – Keep every point evaluated. Not cheap across many fits.

  • solver (dict) – Root-finding options for every record’s likelihood, keyed by the names in SOLVER_OPTIONS. Options left out keep Log10Likelihood’s defaults. To reproduce a published value, pass the settings it was computed with.

cost(x)[source]

Negative summed \(\log_{10}\) likelihood, for a minimiser.

A Q matrix the likelihood cannot handle raises ArithmeticError, or else returns nan silently. Both become FAILURE_COST.

This is not what HJCFIT itself did. It kept the best parameters so far and, on a failure, replaced the current ones with those plus a bounded random perturbation (p. 702) – which its own simplex allows, because it owns its vertices. SciPy’s Nelder-Mead does not expose them, so a large constant penalty is used instead. The failure count is reported with every fit, so it is visible if this ever matters; the paper had two numerical failures in nearly 50 000 fits.

fit(x0=None, search='simplex', restarts=0, xatol=0.0001, fatol=0.0001, maxfev=20000, maxiter=20000, options=None)[source]

Run the search and return a FitResult.

Parameters:

  • x0 – Starting point, in the space being searched. Defaults to the mechanism’s current free parameters.

  • search – "simplex" is simplex_hjc(), HJCFIT’s own – the search that produced every published result, and the default here for that reason. "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 and a random Q matrix usually cannot offer that.

  • options – Passed to whichever search is chosen. The simplex’s own defaults are HJCFIT’s, so this is for deliberate departures.

  • restarts – "scipy" only: extra runs from the previous solution, stopping early when one gains less than fatol. The simplex has its own restart rule, capped at nresmax.

  • xatol – "scipy" only.

  • fatol – "scipy" only.

  • maxfev – Evaluation budget, both searches.

  • maxiter – "scipy" only.

ln_likelihood(x=None)[source]

The same thing in natural logarithms.

Anything treating the log likelihood as a statistical quantity – a Hessian, and so the covariance matrix, the standard deviations and the likelihood intervals – needs natural logarithms. Getting it wrong is not obvious in the result: it inflates every standard deviation by exactly \(\sqrt{\ln 10} = 1.517\), which looks like a badly behaved fit rather than a units error.

log10_likelihood(x=None)[source]

Summed \(\log_{10}\) likelihood of every record.

The likelihood takes one Q matrix per call, so each record is evaluated at its own concentration and the logarithms added.

Raises rather than returning a sentinel; cost() catches.

Parameters:

x – Parameters, in the space being searched. None uses the mechanism as it stands.

Raises:

ArithmeticError – if any record’s likelihood is not finite.

HJCFIT.likelihood.fitting.FAILURE_COST = 10000000000.0

Returned by the cost function when the likelihood cannot be computed. Large, finite and constant – see HJCFitter.cost() for why it is a penalty rather than the perturbation HJCFIT itself used.

HJCFIT.likelihood.fitting.SOLVER_OPTIONS = ('nmax', 'xtol', 'rtol', 'itermax', 'lower_bound', 'upper_bound')

The root-finding options Log10Likelihood accepts beyond the record itself, and so what solver may hold. 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 – gives a natural log-likelihood 0.38 higher than the defaults (nmax=3, 1e-10) at the same rates.

HJCFitter(..., solver=dict(nmax=2, xtol=1e-12)) passes them to every record’s likelihood. Options left out keep the defaults of Log10Likelihood.

The result

class HJCFIT.likelihood.fitting.FitResult(rates: dict, free_names: tuple, free_values: ndarray, log10_likelihood: float, nevals: int, niter: int, nfailures: int, seconds: float, success: bool, message: str, path: list = <factory>, evaluations: list = <factory>)[source]

What one fit produced.

Parameters:

  • rates – Every rate constant by name, constrained ones included.

  • free_names – Names of the free parameters, in order.

  • free_values – Their fitted values, as rates rather than logs.

  • log10_likelihood – \(\log_{10} L_{max}\).

  • nevals – Likelihood evaluations used.

  • niter – Iterations the optimiser reported.

  • nfailures – Evaluations that could not be computed and cost FAILURE_COST. Reported with every fit rather than swallowed, so that a mechanism the likelihood struggles with is visible.

  • seconds – Wall-clock time of the search.

  • success – The optimiser’s own verdict.

  • message – The optimiser’s own words.

  • path – Best point per iteration, when asked for.

  • evaluations – Every point evaluated, when asked for.

Three things about a fit that are easy to get wrong

Log space is the default, and not only for speed. Searching the logarithms of the rate constants is three to four times faster (Colquhoun, Hatton & Hawkes 2003, p. 702) and cannot produce a negative rate. Searching the rates themselves can: on that paper’s AChR mechanism it gave four fits in 250 with negative rate constants until the out-of-range reset was added.

A rate that ends against its limit is not a fit. It is a statement that the likelihood wanted to go somewhere the model forbids. HJCFitter resets rather than constrains, as HJCFIT did, so nothing tells you this happened — compare the result against the limits you set.

Standard errors need natural logarithms. Anything treating the log likelihood as a statistical quantity — a Hessian, and so the covariance matrix, the standard deviations and the likelihood intervals — needs HJCFitter.ln_likelihood(), not the base-10 one. Using log10 inflates every standard deviation by exactly \(\sqrt{\ln 10} = 1.517\), which looks like a badly behaved fit rather than a units error.