Determinant Equation \(W(s) = sI - H(s)\)

The function \(H(s)\) is an integral defined as:

\[H(s) = sI - \mathcal{Q}_{AA} - \mathcal{Q}_{AF}\ \int_0^\tau e^{-st}e^{\mathcal{Q}_{FF}t}\partial\,t\ \mathcal{Q}_{FA}\]

It is possible the function \(H\) as well as its determinant using the DeterminantEq objects. This is the object used when solving for the approximate missed-events likelihood. The determinant equation is initialized in one of two ways, either from a matrix or QMatrix.

python:
from dcprogs.likelihood import DeterminantEq, QMatrix

# Define parameters.
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] ]
qmatrix = QMatrix(matrix, 2)
tau = 1e-4

det0 = DeterminantEq(qmatrix, 1e-4);
det1 = DeterminantEq(matrix, 2, 1e-4);

print("{0!s}\n\n{1!s}".format(det0, det1))
c++11:
#include <iostream>

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

  // Define parameters.
  DCProgs::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;
  DCProgs::QMatrix qmatrix(matrix, /*nopen=*/2);

  // Create determinant using a QMatrix and a matrix+nopen.
  DCProgs::DeterminantEq det0(qmatrix, 1e-4);
  DCProgs::DeterminantEq det1(matrix, 2, 1e-4);

  std::cout << det0 << "\n\n" << det1 << "\n";
            << "  * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";

  return 0;
}

With an object in hand, it is possible to compute \(\mathop{det}W(s)\) for any \(s\). In the following we demonstrate that the two initialization methods are equivalent and that the determinant is zero at the roots of \(W(s)\), per definition.

python:

The python bindings accept both scalars and arrays as input.

x = [0, -1, -1e2]
assert all(abs(det0(x) - det1(x)) < 1e-8)

roots = array([-3045.285776037674, -162.92946543451328])
assert all(abs(det0(roots)) < 1e-6 * abs(roots))
c++11:
  if( std::abs(det0(0) - det1(0)) > 1e-6 
      or std::abs(det0(-1) - det1(-1)) > 1e-6
      or std::abs(det0(-1e2) - det1(-1e2)) > 1e-6)
    throw DCProgs::errors::Runtime("instanciations differ.");

  if( std::abs( det0(-3045.285776037674) ) > 1e-6 * 3e3
      or std::abs( det0(-162.92946543451328) ) > 1e-6 * 2e2 )
    throw DCProgs::errors::Runtime("Roots are not roots.");

There exists a convenience function to transform a determinant equation into its “transpose”, e.g. one where A states become F states and F states become A states:

python:
transpose = det0.transpose()
roots = array( [-17090.192769236815, -2058.0812921673496, -0.24356535498785126],
               dtype=internal_dtype )
assert all(abs(transpose(roots)) < 1e-6 * abs(roots))

Note

Here we choose to create an input which has same internal type as the dcprogs package. This may result in faster code since no conversion are required.
c++11:
  DCProgs::DeterminantEq transpose = det0.transpose();
  if( std::abs( transpose(-17090.192769236815) ) > 1e-6 * 2e5
      or std::abs( transpose(-2058.0812921673496) ) > 1e-6 * 2e3
      or std::abs( transpose(-0.24356535498785126) ) > 1e-6  )
    throw DCProgs::errors::Runtime("Roots are not roots.");

Finally, it is possible to compute \(H(s)\) directly, as well as \(\frac{\partial W(s)}{\partial s}\), as demonstrated below.

python:
print("  * H({0}):\n{1}\n".format(0, det0.H(0)))
print("  * H({0}):\n{1}\n".format(-1e2, det0.H(-1e2)))
print("  * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(0, det0.s_derivative(0)))
print("  * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(-1e2, det0.s_derivative(-1e2)))
c++11:
  std::cout << "  * H(0):\n" << DCProgs::numpy_io(det0.H(0)) << "\n\n"
            << "  * H(-1e2):\n" << DCProgs::numpy_io(det0.H(-1e2)) << "\n\n";

  std::cout << "  * d[sI-H(s)]/ds for s=0:\n" << DCProgs::numpy_io(det0.s_derivative(0)) << "\n\n"
            << "  * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";