Occupancies

The start and end occupancies can be computed one of two way. They can be the equilibrium occupancies determined from the equation:

\[\begin{split}\phi_A = \phi_A {}^e\mathcal{G}_{AF} {}^e\mathcal{G}_{FA},\\ \phi_F = {}^e\mathcal{G}_{FA} {}^e\mathcal{G}_{AF} \phi_F,\\\end{split}\]

subject to the constraints,

\[\begin{split}\sum_i [\phi_A]_i = 1,\\ \sum_i [\phi_F]_i = 1.\end{split}\]

Where \({}^e\mathcal{G}_{AF}\) and \({}^e\mathcal{G}_{FA}\) are the laplacians of the missed-events likelihoods (or equivalently, the ideal likelihoods) for \(s=0\), and \([a]_i\) indicates the \(i^{th}\) component of vector \(a\).

Or they can be computed as CHS vectors, e.g. equation 5.11 and 5.8 from [CHS96].

The occupancies are accessed differently in c++ and in python.

python:

The equilibrium occupancies are accessed as properties of IdealG and MissedEventsG instances. The CHS vectors are functions of these same classes that take as arguments the critical time.

from dcprogs.likelihood import QMatrix, IdealG, 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 = MissedEventsG(qmatrix, tau)
idealG = IdealG(qmatrix)

tcritical = 5e-3

print("Equilibrium Occupancies\n"            \
      "=======================\n\n"          \
      "Ideal Likelihood\n"                   \
      "----------------\n\n"                 \
      "  * initial: {ideal_initial!r}\n"     \
      "  * final: {ideal_final!r}\n\n\n"     \
      "Missed-events Likelihood\n"           \
      "------------------------\n\n"         \
      "  * initial: {equi_initial!r}\n"      \
      "  * final: {equi_final!r}\n\n\n\n"    \
      "CHS Occupancies\n"                    \
      "===============\n\n"                  \
      "Missed-events Likelihood\n"           \
      "------------------------\n\n"         \
      "  * tcritical: {tcritical}\n"         \
      "  * initial: {chs_initial!r}\n"       \
      "  * final: {chs_final!r}"             \
      .format(
        ideal_initial = idealG.initial_occupancies,
        ideal_final   = idealG.final_occupancies,
        equi_initial  = eG.initial_occupancies,
        equi_final    = eG.final_occupancies,
        chs_initial   = eG.initial_CHS_occupancies(tcritical),
        chs_final     = eG.final_CHS_occupancies(tcritical),
        tcritical     = tcritical
      )
)
c++11:

Both equilibrium and CHS occupancies are accessed via function calls acting on the IdealG and MissedEventsG.

#include <iostream>

#include <likelihood/missed_eventsG.h>
#include <likelihood/idealG.h>
#include <likelihood/occupancies.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::t_real const tau(1e-4); // in seconds

  // Create missed-events G
  DCProgs::MissedEventsG eG(qmatrix, tau);
  // Create ideal G
  DCProgs::IdealG idealG(qmatrix);
  
  DCProgs::t_real const tcritical(5e-3);

  std::cout << "Equilibrium Occupancies\n"
            << "=======================\n\n"
            << "Ideal Likelihood\n"
            << "----------------\n\n"
            << "  * initial: " << DCProgs::occupancies(idealG) << "\n"
            << "  * final: "   << DCProgs::occupancies(idealG, false) << "\n\n\n"
            << "Missed-events Likelihood\n"
            << "------------------------\n\n"
            << "  * initial: " << DCProgs::occupancies(eG) << "\n"
            << "  * final: "   << DCProgs::occupancies(eG, false) << "\n\n\n\n"
            << "CHS Occupancies\n"
            << "===============\n\n"
            << "Missed-events Likelihood\n"
            << "------------------------\n\n"
            << "  * tcritical: " << tcritical << "\n"
            << "  * initial: " << DCProgs::CHS_occupancies(eG, tcritical) << "\n"
            << "  * final: "   << DCProgs::CHS_occupancies(eG, tcritical, false) << "\n";

  return 0;
}

In c++, the occupancies are kept outside of the classes because computing these values is outside the pure remit of the classes (which is to compute the likelihood). However, in python, practicality beats purity, and it makes practical sense to keep likelihood and occupancies together.

Note

Log10Likelihood() uses equilibrium occupancies depending on the value of its attribute tcritical:

  • if it is None, numpy.NaN, or negative, then the equilibrium occupancies are used
  • if it a strictly positive real number, then the CHS vectors are computed