Missed-Events transition densities \({}^eG(t)\)¶

The callable object HJCFIT::MissedEventsG provides an interface to compute the likelihood \({}^eG(t)\) of open and shut events as a function of their lengths, for a fixed \(Q\)-matrix. It has the ability to compute both exact and approximate missed-events likelihood, returning one or the other depending on a given time resolution.

The asymptotic expression of the transition densities can be computed from the knowledge of the roots of a specific equations. On the one hand, root-finding can be a fairly difficult numerical operation. On the other, it would be more convenient if we can initialize HJCFIT::MissedEventsG directly from a \(Q\)-matrix object. As such, there are several means to initialize the functor:

from the knowledge of the roots and the determinant equations;
directly from a \(Q\)-matrix, using the default root-finding mechanism;
from a \(Q\)-matrix, using a custom root-finding mechanism (C++ only).

Initialization from a \(Q\)-matrix¶

Python:

from numpy import all, abs, arange
from HJCFIT.likelihood import QMatrix, DeterminantEq, MissedEventsG

# Define parameters.
qmatrix = QMatrix([ [-3050,        50,  3000,      0,    0], 
                    [2./3., -1502./3.,     0,    500,    0], 
                    [   15,         0, -2065,     50, 2000], 
                    [    0,     15000,  4000, -19000,    0], 
                    [    0,         0,    10,      0,  -10] ], 2)
tau = 1e-4
 
eG_automatic = MissedEventsG(qmatrix, tau)

C++11:

#include <iostream>
#include <exception>

#include <likelihood/missed_eventsG.h>
#include <likelihood/root_finder.h>
 
int main() {

  // Define parameters.
  HJCFIT::t_rmatrix matrix(5 ,5);
  matrix << -3050,        50,  3000,      0,    0, 
            2./3., -1502./3.,     0,    500,    0,  
               15,         0, -2065,     50, 2000,  
                0,     15000,  4000, -19000,    0,  
                0,         0,    10,      0,  -10;
  HJCFIT::QMatrix qmatrix(matrix, /*nopen=*/2);
  HJCFIT::t_real const tau(1e-4); // in seconds

  HJCFIT::MissedEventsG eG_automatic(qmatrix, tau);

  return 0;
}

A fair amount of work goes on behind the scene. First reasonable upper and lower bounds for the roots obtained (HJCFIT::find_lower_bound_for_roots(), and HJCFIT::find_upper_bound_for_roots()). Then intervals for each roots are computed (HJCFIT::find_root_intervals()). And finally, the roots themselves are obtained (HJCFIT::brentq()). All this work is done automatically in the case of this particular instantiation. A few extra parameters to control the root-finding process can be passed to the C++ and Python constructors.

Initialization from the roots and determinant equations¶

Python:

determinant_eq = DeterminantEq(qmatrix, tau)
af_roots = [( -3045.285776037674, 1), (-162.92946543451328, 1)]
fa_roots = [(-17090.192769236815, 1), (-2058.0812921673496, 1), (-0.24356535498785126, 1)]
eG_from_roots = MissedEventsG(determinant_eq, af_roots, determinant_eq.transpose(), fa_roots)

C++11:

  HJCFIT::DeterminantEq determinant_eq(qmatrix, tau);
  std::vector<HJCFIT::Root> af_roots{
    { /*root=*/ -3045.285776037674,   /*multiplicity=*/ 1}, 
    { /*root=*/ -162.92946543451328,  /*multiplicity=*/ 1}
  };
  std::vector<HJCFIT::Root> fa_roots{
    { /*root=*/ -17090.192769236815,      /*multiplicity=*/ 1},
    { /*root=*/  -2058.0812921673496,     /*multiplicity=*/ 1},
    { /*root=*/     -0.24356535498785126, /*multiplicity=*/ 1}
  };
  HJCFIT::MissedEventsG eG_from_roots( determinant_eq, af_roots, 
                                        determinant_eq.transpose(), fa_roots );

In this case, it is expected the roots are known, somehow, as well as their multiplicity. This format allows external root-finding methods to be interfaced with the packages.

Initialization from the \(Q\)-matrix and a root finding function¶

Given a root-finding function, it is possible to instantiate \({}^eG\). The root finding function should take a determinant equation as input, and return a vector of HJCFIT::Root as output. In the code below, we show how the prior initialization could be recreated.

  auto find_roots = [](HJCFIT::DeterminantEq const &_det) {
    return HJCFIT::find_roots(_det, 1e-12, 1e-12, 100, HJCFIT::quiet_nan, HJCFIT::quiet_nan);
  };
  HJCFIT::MissedEventsG eG_from_func(qmatrix, tau, find_roots);

This is mostly a convenience function, to make it slightly easier to interface with other root-finding methods in C++. This interface is not explicitly available in Python, although it can be created with ease.

Computing \({}^eG_{AF}(t)\) and \({}^eG_{FA}(t)\)¶

The computation of the likelihood matrices is illustrated by comparing the three initialization methods (two in python). It is clear that all three yield the same function.

python:

assert all(abs(eG_from_roots.af(tau) - eG_automatic.af(tau)) < 1e-8)
assert all(abs(eG_from_roots.fa(tau) - eG_automatic.fa(tau)) < 1e-8)

# Checks the three initialization are equivalent at different times
# The functions can be applied to arrays. 
times = arange(tau, 10*tau, 0.1*tau)
assert eG_from_roots.af(times).shape == (len(times), 2, 3)
assert eG_from_roots.fa(times).shape == (len(times), 3, 2)
assert all(abs(eG_from_roots.af(times) - eG_automatic.af(times)) < 1e-8)
assert all(abs(eG_from_roots.fa(times) - eG_automatic.fa(times)) < 1e-8)

c++11:

  for(HJCFIT::t_real t(tau);  t < 10*tau; t += tau * 0.1) {

    if(    ((eG_from_roots.af(t) - eG_from_func.af(t)).array().abs() > 1e-8).any()  
        or ((eG_from_roots.fa(t) - eG_from_func.fa(t)).array().abs() > 1e-8).any() )
      throw HJCFIT::errors::Runtime("root != func");

    if(    ((eG_from_roots.af(t) - eG_automatic.af(t)).array().abs() > 1e-8).any() 
        or ((eG_from_roots.fa(t) - eG_automatic.fa(t)).array().abs() > 1e-8).any() )
      throw HJCFIT::errors::Runtime("root != automatic");
  }

The python bindings accept any input that can be transformed to a numpy array of reals. If the input is a scalar, then the AF and FA blocks are returned. If the input is an array, then an array of similar shape is returned, where each component is a matrix.

The HJCFIT::MissedEventsG provides further functionality. For instance, the cutoff point between exact and asymptotic calculations can be set explicitly (it defaults to \(t < 3\tau\)). And the likelihood can be computed in Laplace space (see HJCFIT::MissedEventsG::laplace_af and HJCFIT::MissedEventsG::laplace_fa). We invite users to turn to the python and the c++ API for more details.