The \(Q\)-matrix¶

The \(Q\)-matrix represents the most basic data structure of this library. It holds the rates between different states in the mechanism, for open and shut states. In practice, it can be represented by a matrix of rates and the number of open states. As such, the initialization takes both of these input:

from numpy import sum, abs, all
from dcprogs.likelihood import QMatrix

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);
 
print(qmatrix)

#include <iostream>
#include <typeinfo>

#include <exception>
#include <likelihood/qmatrix.h>
 

int main() {

  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);
 
  std::cout << qmatrix << std::endl;
 
  return 0;
}

In both cases, the open-open transitions must be located in the top-left corner of the input matrix.

QMatrix is a structure instance, without private members. Hence both the matrix and number of open states can be accessed directly. In the following, we show access to nopen and to the rates. Among other things, the code below verifies that the rates of the CH82 model do sum to zero over rows (this condition is not enforced by the QMatrix object).

assert qmatrix.nopen == 2
assert qmatrix.nshut == 3
assert all(abs(qmatrix.matrix - matrix) < 1e-12)

  
# Compute sum over rows, row by row.
for i, row in enumerate(qmatrix.matrix): print("sum(row[{0}]) = {1}".format(i, sum(row)))

# Compute sum over rows, but let numpy do it.
print("Sum over rows: {0}".format(sum(qmatrix.matrix, axis=1)))
  
# Check that sum over rows returns zero.
assert all(abs(sum(qmatrix.matrix, axis=1)) < 1e-12)

  if(qmatrix.nopen != 2) return 1;
  if(qmatrix.nshut() != 3) return 1;
  if( ((qmatrix.matrix - matrix).array().abs() > 1e-12).any() ) return 1;
  
  // Compute sum over rows, row by row.
  for(DCProgs::t_int i(0); i < qmatrix.matrix.rows(); ++i) 
    std::cout << "sum(row[" << i << "]): " << qmatrix.matrix.row(i).sum() << std::endl;

  // Compute sum over rows, but let eigen do it.
  std::cout << "Sum over rows: " << qmatrix.matrix.rowwise().sum().transpose() << std::endl;
  
  // Check that sum over rows returns zero.
  if( (qmatrix.matrix.rowwise().sum().array().abs() > 1e-10).any() ) throw std::exception();

Finally, the blocks of the \(Q\) can be accessed efficiently via functions (or properties in python). These are wrappers around the rate matrix. These wrappers are numerically efficient. They can and should be used within numerical algorithms.

# Gets AA block.
print("AA block:\n---------\n{0}".format(qmatrix.aa))
                
# aa is simply wrapper around the matrix block 
qmatrix.aa[0, 0] = 42
assert abs(qmatrix[0, 0] - 42) < 1e-12

  // Gets AA block.
  std::cout << "AA block:\n---------\n" << qmatrix.aa() << std::endl;
                
  // aa returns a wrapper around the matrix block. Hence memory is the same.
  if(&(qmatrix.aa()(0,0)) != &(qmatrix(0, 0))) throw std::exception();