Ideal Likelihood \(\mathcal{G}(t)\)

A wrapper around QMatrix is provided which allows the calculation of the ideal likelihood:

\[\begin{split}\mathcal{G}(t) = \left( \begin{eqnarray} e^{-t\mathcal{Q}_{af}} & 0 \\ 0& e^{-t\mathcal{Q}_{fa}} \\ \end{eqnarray} \right)\end{split}\]

This object can be initialized directly from a QMatrix.

python:
from dcprogs.likelihood import QMatrix, IdealG, expm

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)

idealG = IdealG(qmatrix)
print(idealG)
c++11:
#include <iostream>
#include <exception>

#include <likelihood/idealG.h>
#include <likelihood/errors.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);

  DCProgs::IdealG idealG(qmatrix);

  std::cout << idealG << std::endl;
  return 0;
}

It provides the ideal likelihood as a function of time, as well as the laplace transforms:

python:
from numpy import dot, identity, abs, all

idealG_fa = dot(expm(2e-4 * qmatrix.ff), qmatrix.fa)
assert all( abs(idealG_fa - idealG.fa(2e-4)) < 1e-8 )

inversion = -0.5 * identity(2) - qmatrix.aa
assert all( abs( dot(inversion, idealG.laplace_af(-0.5)) - qmatrix.af ) < 1e-8 )

Note

This package can be compiled to use real number of 128bits. However, despite providing arrays for reals of that size, numpy does not include the corresponding linear algebra functions. A few useful functions, such as expm in this example, are provided to remediate to this situation.
c++11:
  DCProgs::t_rmatrix const idealG_fa = (2e-4*qmatrix.ff()).exp()*qmatrix.fa();
  if( ((idealG.fa(2e-4) - idealG_fa).array().abs() > 1e-8).any() )
    throw DCProgs::errors::Runtime("Not so ideal idealG");

  DCProgs::t_rmatrix const inversion = -0.5 * DCProgs::t_rmatrix::Identity(2, 2) - qmatrix.aa();
  if( ((inversion * idealG.laplace_af(-0.5) - qmatrix.af()).array().abs() > 1e-8).any() )
    throw DCProgs::errors::Runtime("Not so ideal idealG");