Bounding H2MM#
See also
This can also be viewed as a Jupyter Notebook
Downlaod H2MM_Bounds_Tutorial.ipynb
The data file can be downloaded here: sample_data_3det.txt
Let’s get our obligatory imports in order, this time we’ll start with the 3 detector data.
import os
import numpy as np
from matplotlib import pyplot as plt
import H2MM_C as hm
# load the data
color3 = list()
times3 = list()
i = 0
with open('sample_data_3det.txt','r') as f:
for line in f:
if i % 2 == 0:
times3.append(np.array([int(x) for x in line.split()],dtype='Q'))
else:
color3.append(np.array([int(x) for x in line.split()],dtype='L'))
i += 1
Bounding Parameter Values#
Sometimes we may want to restrain the range of possible values in an H2MM model (the h2mm_model).
For instance, we may want to keep transition rates (h2mm_model.trans) within values reasonable for the duration of the experiment, or because you want restrain the values in the emission probability matrix (h2mm_model.obs) based on how we know the system should behave.
Bounds are defined using h2mm_limits objects.
This object is passed to EM_H2MM_C() (or h2mm_model.optimize() in an object-oriented approach) through the bounds keyword argument.
When this is specified, it is also necessary to specify another keyword argument, bounds_func.
So let’s see an example
alt_period = 4000 # a fake alternation period
us_bounds = hm.h2mm_limits(max_trans = 1/(alt_period))
prior = np.array([1/4, 1/4, 1/4, 1/4])
trans = np.array([[1-3e-6, 1e-6, 1e-6, 1e-6],
[1e-6, 1-3e-6, 1e-6, 1e-6],
[1e-6, 1e-6, 1-3e-6, 1e-6],
[1e-6, 1e-6, 1e-6, 1-3e-6]])
obs = np.array([[0.4,0.4,0.2],
[0.3,0.1,0.6],
[0.2,0.4,0.4],
[0.1,0.1,0.8]])
imodel_4s3d = hm.h2mm_model(prior, trans, obs)
us_opt_model4 = hm.EM_H2MM_C(imodel_4s3d, color3, times3, bounds_func='revert', bounds=us_bounds)
us_opt_model4
So, what did we just to?
The h2mm_limits object us_bounds prevents any value (except on-diagonal values) of the transition probability matrix (h2mm_model.trans) from ever being larger (i.e. faster transition rate) than 1/4000
Bounds Process#
When you use a bounds method, each iteration goes through the following steps:
Calculate the loglikelihood and a new model
Check if the optimization has converged
Analyze the new model, and correct if necessary.
Check if any values are smaller or larger than the limits set in the
boundsargument.If values are out of bounds, apply a correction, based on the method defined by the argument passed to
bounds_func
Repeat optimization (return to step 1)
When creating a h2mm_limits object, all arguments are passed as keyword arguments.
They come in the form of min/max_[name] where [name] is prior, trans, or obs (the core arrays of h2mm_model objects).
These specify the minimum/maximum values in the respective array.
If no value is specified for a given min/max array, then the values of that array can be as small or as large as possible for an unbounded model.
In all cases, these values can be specified as either a float or an array.
- If specified as a float, then the value is universal for all values in the given array (this is most useful for the h2mm_model.trans matrix). This means less granularity, but then the h2mm_limits object can be used for any number of states/detectors
- If specified as an array, then the values in the array are matched with the cooresponding array in h2mm_model. This means greater granularity, but then you are locked into using the h2mm_limits object only for optimizations with a specific number of states/detectors.
As mentioned in the above outline, you also need to specify the bounds_func argument when using bounds.
There are 3 options for this:
'minmax': The shallowest correction, sets the out of bounds value to its minimum or maximum'revert': The preferred method, sets the out of bounds value to the value it was in the previous model'revert_old': A more extreme form of'revert'which goes to the model before the last in the optimization, and sets the out of bounds value to that “older” value.
Using factory_h2mm_model() with Bounds#
You will note that in the previous example, we specified the h2mm_model explicity, instead of using the factory_h2mm_model() function.
This is because it is possible that the factory_h2mm_model() could return a h2mm_model whose values are already out of bounds based for the h2mm_limits object.
This could create odd behaviour during optimization.
There is a way around this: you can call factory_h2mm_model() with the h2mm_limits object passed through the bounds keyword argument, and the function will automatically return a h2mm_model object that is within the bound provided.
Note
See the documentation of factory_h2mm_model() for a fill list of options for customizing the function’s output.
us_bounds = hm.h2mm_limits(max_trans = 1/4000)
# make factory_h2mm_model make a model within bounds
imodel = hm.factory_h2mm_model(3,3, bounds=us_bounds)
us_model = hm.EM_H2MM_C(imodel, color3, times3, bounds=us_bounds, bounds_func='revert')
Custom Bounds#
Finally, it is possible to supply a custom bounding function to bounds_func.
Note
This feature was designed to allow the user to handle things/circumstances that the writers of H2MM_C had not anticipated. Therefore this example is very simple, and does not show a useful method.
This function must is called at the end of the optimization loop, after a given model’s loglik has been calculated (and the standard H2MM next model for the next iteration produced).
This function takes the signature
bounds_func(new:h2mm_model, current:h2mm_model, old:h2mm_model, *bounds, **bounds_func)->h2mm_model|int|tuple[h2mm_model,int]
new, current and old are the h2mm_model s of the current iteration.
newis the model suggested/produced by the current iterationcurrentis the model whose loglik was just computedoldis the model computed in the previous iteration
These are always supplied each iteration. bounds and bounds_kwargs come from the
identically named keyword arguments in EM_H2MM_C(). bounds by default is None, which
is internally converted to a 0-size (empty) tuple, likewise bounds_kwargs is by default None
and is internally converted into an empty dict.
The return value can either or both specify
The “bounded” new model
If the optimization has converged
If only the new model is specified, convergence will be determined like all other optimizaztions,
by the difference in loglik of current and old.
However, if the bounds function returns a value specifying if the model has converged, then
EM_H2MM_C() will not separately check if the optimization has converged.
Note
max_iter and max_time are enforced separetely from bounds_func.
If specifying the converged state, this can be either a bool or 0, 1, 2.
As a bool, True indicates that the optimization has converged, and thus
can stop, the old model will be returned as the “optimal” model. False
will allow the optimization to proceed using the new model.
If specifying as 0 is equivalent to False, 1 to True, and 2 will return
the current model as the optimal model.
If both the new model and converged state are specified, this must be done by returning a 2-tuple
of (new, converged_state).
Warning
Make sure the model you return makes sense, otherwise the optimization will proceed unpredictably. Think of this as the “gloves off” approach, you might have a very powerful new method, or you might get something meaningless depending on how you code it. That’s your responsibility.
Below is a function that that re-implements the behavior of "minmax" but now the limits
normally specified with a h2mm_limits object supplied to bounds are replaced with kwargs:
def minmax_py(new, current, old, converged_min=1e-9,
min_prior=None, max_prior=None,
min_trans=None, max_trans=None,
min_obs=None, max_obs=None):
# bounding of trans matrix
if min_trans is not None or max_trans is not None:
trans = new.trans
idxs = np.arange(new.nstate)
if isinstance(min_trans, float):
trans[trans < min_trans*(~np.eye(new.nstate, dtype=np.bool_))] = min_trans
elif isinstance(min_trans, np.ndarray):
mask = trans < min_trans
trans[mask] = min_trans[mask]
if isinstance(max_trans, float):
trans[trans > max_trans*(~np.eye(new.nstate, dtype=np.bool_))] = max_trans
elif isinstance(max_trans, np.ndarray):
mask = trans > max_trans
trans[mask] = max_trans[mask]
for i in range(trans.shape[0]):
trans[i,i] = 1.0 - trans[i, idxs!=i].sum()
new.trans = trans
# bounding of obs matrix
if min_obs is not None or max_obs is not None:
obs = new.obs
if min_obs is not None:
minmask = obs < min_obs
obs[minmask] = min_obs[minmask]
else:
minmask = np.zeros(obs.shape, dtype=np.bool_)
if max_obs is not None:
maxmask = obs > max_obs
obs[maxmask] = max_obs[maxmask]
else:
maxmask = np.zeros(obs.shape, dtype=np.bool_)
obsmask = minmask | maxmask
for i in range(obs.shape[0]):
obs[i,~obsmask[i,:]] += (1-obs[i,:].sum()) / (~obsmask).sum()
new.obs = obs
if min_prior is not None or max_prior is not None:
prior = new.prior
if min_prior is not None:
minpmask = prior < min_prior
prior[minpmask] = min_prior[minpmask]
else:
minpmask = np.zeros(new.nstate, base=np.bool_)
if max_prior is not None:
maxpmask = prior > max_prior
prior[maxpmask] = max_prior[maxpmask]
else:
maxpmask = np.zeros(new.nstate, base=np.bool_)
pmask = minpmask | maxpmask
prior[~pmask] += (1-prior.sum()) / (~pmask).sum()
new.prior = prior
return new
Now let’s see it in action
prior = np.array([1/4, 1/4, 1/4, 1/4])
trans = np.array([[1-3e-6, 1e-6, 1e-6, 1e-6],
[1e-6, 1-3e-6, 1e-6, 1e-6],
[1e-6, 1e-6, 1-3e-6, 1e-6],
[1e-6, 1e-6, 1e-6, 1-3e-6]])
obs = np.array([[0.09, 0.01, 0.9],
[0.3, 0.1, 0.6],
[0.2, 0.4, 0.4],
[0.1, 0.1, 0.8]])
imodel4s3d = hm.h2mm_model(prior, trans, obs)
us_opt_model4 = imodel4s3d.optimize(color3, times3, bounds_func=minmax_py, bounds_kwargs=dict(max_trans=1e-4))
us_opt_model4
Now let’s re-implement how convergence of the optimization is handled
def limit_converged(new, current, old, conv_min):
if current.loglik < old.loglik:
return 1
if (current.loglik - old.loglik) < conv_min:
return 2
return 0
us_opt_model4 = imodel_4s3d.optimize(color3, times3, bounds_func=limit_converged, bounds=5e-8)
us_opt_model4
Finally, bellow is an example that re-implements the min-max procedure and checking for convergence
def minmax_conv_py(new, current, old, conv_min, **kwargs):
return minmax_py(new, current, old, **kwargs), limit_converged(new, current, old, conv_min)
us_opt_model4 = imodel_4s3d.optimize(color3, times3, bounds_func=minmax_conv_py, bounds=5e-8, bounds_kwargs=dict(max_trans=1e-4))
us_opt_model4