Optimization

HJCFIT computes a likelihood; maximising it is the caller’s business. Two things here help with that: the simplex the Fortran HJCFIT used, and the mapping from a reaction graph to the free parameters of a fit.

The simplex HJCFIT used

Every published HJCFIT result, Colquhoun, Hatton & Hawkes (2003) among them, was obtained with the simplex in SIMPHJC.FOR. It is not scipy.optimize.minimize() with method='Nelder-Mead', and the differences are not cosmetic — see simplex_hjc() for what they are.

HJCFIT.likelihood.optimization.simplex_hjc(fun, x0, args=(), logfit=True, stpfac=None, confac=None, resfac=None, nresmax=None, errfac=None, maxfev=None, local_search='corrected', callback=None)[source]

Minimise fun the way the Fortran HJCFIT did.

Parameters:

  • fun – Callable, fun(x, *args) -> float. In HJCFIT this is minus the log-likelihood.

  • x0 – Starting point. In log mode this is ``log(rates)``: the caller takes logarithms before the call and exponentiates after, which is what the Fortran does around its own call to SIMPHJC.

  • args – Extra arguments passed to fun.

  • logfit – Use the log-space step and convergence rules. This does not transform x0; it says what x0 already is, so that the step can be the same for every parameter and the convergence test can be absolute in log(rate). HJCFIT’s own prompt, “Use log(rate constant) in Simplex”, defaults to yes in every version of the source in DCPROGS/DCFORTRAN, so True is the faithful default.

  • stpfac – The step factor as the program prompts for it: 5.0 in log mode, 0.2 in linear. In log mode its logarithm becomes the step, so every parameter starts a factor of stpfac from the guess.

  • confac – Contraction factor, also used for shrinkage.

  • resfac – On a restart the step is reset to resfac * crtstp.

  • nresmax – Most restarts allowed. 0 disables them.

  • errfac – Converged when the simplex spans less than this in every parameter; in log mode it is absolute in log(rate), so 1e-3 is a tenth of a per cent of the rate.

  • maxfev – Evaluation budget. Exhausting it returns success=False.

  • local_search – "corrected" (the default), "fortran" or "off". The original’s local search contains a bug; "fortran" reproduces it, for reproducing published results, and nothing else should use it. See the module docstring.

  • callback – Called as callback(x) with the best vertex after each iteration, as scipy.optimize.minimize() calls it.

Returns:

a SimplexHJCResult

from HJCFIT.likelihood.optimization import simplex_hjc

result = simplex_hjc(lambda x: -likelihood(np.exp(x)),
                     np.log(theta))
rates = np.exp(result.x)
class HJCFIT.likelihood.optimization.SimplexHJCResult(x, fun, nfev, nit, nrestarts, iconv, success, message)[source]

The result of a simplex_hjc() run.

The field names are SciPy’s where SciPy has one, so that a script can swap scipy.optimize.minimize() for simplex_hjc() and keep reading res.x, res.fun, res.nfev, res.nit, res.success and res.message.

Parameters:

  • x – the point returned, in the space searched

  • fun – the function value there

  • nfev – likelihood evaluations used

  • nit – iterations

  • nrestarts – restarts taken, at most nresmax

  • iconv – which of the four endings was taken; see ICONV. 6 means the evaluation budget ran out

  • success – False only for iconv == 6

  • message – iconv in words

HJCFIT.likelihood.optimization.SIMPLEX_HJC_DEFAULTS = {'confac': 0.5, 'errfac': 0.001, 'maxfev': 20000, 'nresmax': 3, 'resfac': 10.0, 'stpfac_lin': 0.2, 'stpfac_log': 5.0}

dict() -> new empty dictionary dict(mapping) -> new dictionary initialized from a mapping object’s

(key, value) pairs

dict(iterable) -> new dictionary initialized as if via:

d = {} for k, v in iterable:

d[k] = v

dict(**kwargs) -> new dictionary initialized with the name=value pairs

in the keyword argument list. For example: dict(one=1, two=2)

Keeping rate constants in range

SIMPHJC.FOR has no notion of bounds, and neither does the port: the limits belonged to the program around the subroutine. If you search the rate constants themselves you have to supply that yourself, or the search will happily propose a negative rate constant.

HJCFIT.likelihood.optimization.reset_out_of_range(fun, lower=None, upper=None, logfit=True)[source]

Wraps an objective so a parameter leaving its range is reset, not searched.

simplex_hjc() is a port of SIMPHJC.FOR and, faithfully, has no notion of bounds: the subroutine never had any. The limits belonged to the program around it. Colquhoun, Hatton & Hawkes (2003) p. 702 describes both halves – an upper limit “to prevent physically unrealistic values”, and a floor because “if a value of a rate constant should go negative during the fitting process, it can be reset to a value near zero”.

This is that resetting, and it matters. Searching the rate constants themselves without it produced four fits in 250 with negative rate constants on the AChR mechanism of that paper. In log space it cannot happen, which is one reason the log search is the default there and here.

The reset is applied to the parameters handed to fun, not to the simplex’s own vertices, which is where the original put it: the search may still step outside, and simply learns that it gains nothing by doing so. Clamping the vertices instead would change the geometry of the simplex rather than the function it sees.

Parameters:

  • fun – The objective. Called with the reset parameters.

  • lower – Lower limits, in rate space, broadcast against the parameters. None leaves rates unbounded below – which, searching rates directly, means nothing stops a negative one.

  • upper – Upper limits, in rate space. Hjcfit1-09122003.for prompts for one and the paper’s fits used 1e6 for every rate.

  • logfit – True when the search is over the logarithms of the rate constants, as simplex_hjc() is by default. The limits are always given as rates; this says how to get from the search coordinates to them.

Returns:

A callable with the same signature as fun.

from HJCFIT.likelihood.optimization import (
    simplex_hjc, reset_out_of_range)

bounded = reset_out_of_range(cost, lower=1e-12, upper=1e6)
result = simplex_hjc(bounded, log(guess))

Note what this does not do. It does not report that a limit was reached, and a fit that ends against one is not a fit – it is a statement that the likelihood wanted to go somewhere the model forbids. Check the result against the limits you passed.

A worked shape, for one record and a mechanism whose free rates you can get at as a vector:

from numpy import log, exp
from HJCFIT.likelihood import Log10Likelihood
from HJCFIT.likelihood.optimization import simplex_hjc, reset_out_of_range

likelihood = Log10Likelihood(bursts, nopen=nopen, tau=tau, tcritical=tcrit)

def cost(x):
    """Negative log10 likelihood at log rate constants ``x``."""
    mechanism.set_free_rates(exp(x))
    return -likelihood(mechanism.Q)

result = simplex_hjc(reset_out_of_range(cost, lower=1e-12, upper=1e6),
                     log(guess))
print(result.fun, result.nfev, result.message)
rates = exp(result.x)

Three things about that are worth stating rather than leaving to be rediscovered.

The search is over log(rate) because that is HJCFIT’s own default, it is three to four times faster (Colquhoun, Hatton & Hawkes 2003, p. 702), and a logarithm cannot go negative. The fits of that paper’s Figures 2–5 and 12–13 were made over the rates themselves, which is why the resetting exists at all.

result.fun is the value of whatever you minimised. If that is the negative log10 likelihood, as above, then the log10 likelihood is -result.fun. Mixing the two conventions between the objective and the running best is an easy mistake and does not announce itself.

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

What it is not for

simplex_hjc is a local refiner started from a guess a person believes in, searching the logarithms of the rate constants. That is how HJCFIT used it, and it is where it works.

It is not a global search from a random Q matrix. Tried that way on CH82 with eight free parameters it does not converge at all, and the reason is not the simplex: of 200 random reduced coordinate vectors, 6 give a finite likelihood. The starting simplex has nine vertices, so almost every one of them lands somewhere the likelihood cannot be evaluated and there is nothing to walk on. A derivative-free simplex needs a valid starting simplex, not merely a valid starting point.

For a random multi-start over a mechanism you have no guess for, a method that handles constraints directly — COBYLA or SLSQP through scipy.optimize.minimize() — will get you into a plausible region. simplex_hjc is what to finish with, and what to quote.

Reducing a likelihood to its free parameters

HJCFIT.likelihood.optimization.reduce_likelihood(likelihood, graph_matrix)[source]

Maps likelihood to a set of variable components.

The goal is a callable that takes on input a numpy array with only variable components. It hides from the input the components that are fixed or can be obtained as a result of an expression.

param likelihood:

This should be a callable object that takes a QMatrix or numpy matrix on input.

param graph_matrix:

Defines the reaction graph. This is a list of list (matrix) where each element is either “V” (variable), 0 (no direct reaction between these two states), a number (Fixed reaction rate), or a string with a valid python expression (that eval understands). In the latter case, q will be replaced with the value of the qmatrix, i is set to the current row index, and j to the current column index. The open-states should be in the top-left corner.

Returns:

A callable from which the fixed components have been abstracted. It takes on input a numpy vector with as many components as there are variable components in graph_matrix.

For convenience, the callable has a to_reduced_coords method which takes a numpy matrix and returns a vector with only the variable components. It also sports a to_full_coords coords that maps back to the whole space.