max_ll_util.cpp.html

max_ll_util.cpp.html	mathcode2html
Source file: max_ll_util.cpp
Converted: Sun May 10 2015 at 16:06:29
This documentation file will not reflect any later changes in the source file.

$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

/*

See also:

Up: Index

*/

#ifdef _WIN32
  #include <windows.h>
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <iostream>
#include <fstream>
using std::cout;
using std::cerr;
using std::endl;
#include <set>
using std::set;
#include <map>
using std::map;
#include <string.h>

#include "max_ll_util.h"
#include "ublas_matrixutil.h"

//std::ofstream logg;

#define MIN_ACTIVE_Q 0.0

/*

The Q matrix of rate constants (probability per second) is computed from intrinsic rate constants K0 and K1, and potentially Ligand- and Voltage-sensitive. \[Q_{a,b} = K0_{a,b} * Ligand_{a,b} * e^{K1_{a,b} * Voltage_{a,b}}\] \[Q_{a,a} = - \sum_i Q_{a,i}\]

*/
QUBFAST_API void BuildQ(int Ns, matrix<double>& Q, matrix<double>& K0, matrix<double>& K1, matrix<int>& Ligand, matrix<int>& Voltage, double *Constants)
{
  if ( ((int) Q.size1()) != Ns || ((int) Q.size2()) != Ns )
    Q.resize(Ns, Ns);
  for (int a=0; a<Ns; ++a) {
    double rowsum = 0.0;
    for (int b=0; b<Ns; ++b) {
      if (a != b) {
	double k = K1(a,b);
	if ( Voltage(a,b) )	  k *= Constants[ Voltage(a,b) ];
	k = exp(k);
	k *= K0(a,b);
	if ( Ligand(a,b) ) {
	  double lig = Constants[ Ligand(a,b) ];
	  k *= lig ? lig : MIN_ACTIVE_Q;
	}
	Q(a,b) = k;
	rowsum += k;
      }
    }
    Q(a,a) = - rowsum;
  }
}

QUBFAST_API void BuildQ_p(int Ns, matrix<double>& Q, matrix<double>& K0, matrix<double>& K1, matrix<double>& K2, matrix<int>& Ligand, matrix<int>& Voltage, matrix<int>& Pressure, double *Constants)
{
  if ( ((int) Q.size1()) != Ns || ((int) Q.size2()) != Ns )
    Q.resize(Ns, Ns);
  for (int a=0; a<Ns; ++a) {
    double rowsum = 0.0;
    for (int b=0; b<Ns; ++b) {
      if (a != b) {
	double k = K0(a,b);
	if ( Ligand(a,b) ) {
	  double lig = Constants[ Ligand(a,b) ];
	  k *= lig ? lig : MIN_ACTIVE_Q;
	}
	double kExp = 0.0;
	if ( Voltage(a,b) )
	  kExp += K1(a,b) * Constants[Voltage(a,b)];
	if ( Pressure(a,b) )
	  kExp += K2(a,b) * Constants[Pressure(a,b)];
	k *= exp(kExp);
	Q(a,b) = k;
	rowsum += k;
      }
    }
    Q(a,a) = - rowsum;
  }
}

/*

"Qe" aka ${}^eQ$ is the "apparent" rate constant matrix, given that events with duration <= tdead are not recorded. MIL computes the probability of staying in class a for duration t using the submatrix ${}^eQ_{aa}$ \[A(a, t) = e^{{}^eQ_{aa} (t-tdead)}\]

*/
QUBFAST_API void QtoQe(matrix<double>& Q, int *clazz, double tdead, matrix<double>& eQ)
{
  /*

eQ is partitioned into its Na x Nb submatrices, for classes a and b (likewise partition Q) each partition is computed separately. Notation: going from class a to class b, z is "any class but a", and c is "any class but a or b". Iz is the identity matrix with dimensions of Qzz.

$a = b$: ${}^eQ_{aa} = Q_{aa} - Q_{az} * (I_z - e^{Q_{zz} * tdead}) * Q_{za}$

(accounts for going to another class for a short, then returning undetected)

$a \ne b$: ${}^eQ_{ab} = (Q_{ab} - (Q_{ac} * (I_c - e^{Q_{cc} * tdead}) * Q^{-1}_{cc} * Q_{cb}) * e^{Q_{bb} * tdead}$

(accounts for going to class c for a short, on the way from a to b, then staying in b long enough to be detected)

  */

  int n = (int) Q.size1();
  if ( ((int) eQ.size1()) != n || ((int) eQ.size2()) != n )
    eQ.resize(n, n);

  if ( tdead == 0.0 ) {
    for (int i=0; i<n; ++i)
      for (int j=0; j<n; ++j)
	eQ(i,j) = Q(i,j);
    return;
  }

  set<int> classet; // no duplicates, sorted increasing, for iteration
  for (int i=0; i<n; ++i)
    classet.insert(clazz[i]);
  for (set<int>::iterator cls_i = classet.begin(); cls_i != classet.end(); ++cls_i) {
    std::vector<int> states, outs;
    for (int i=0; i<n; ++i)
      if ( *cls_i == clazz[i] )
	states.push_back(i);
      else
	outs.push_back(i);
    int Na = (int) states.size();
    int Nz = (int) outs.size();
    int *sts = & (states[0]);
    int *ots = & (outs[0]);
    matrix<double> Qaa(Na, Na);        matrix_to_partition(Q, Na, Na, sts, sts, Qaa);
    matrix<double> Qaz(Na, Nz);        matrix_to_partition(Q, Na, Nz, sts, ots, Qaz);
    matrix<double> Qza(Nz, Na);        matrix_to_partition(Q, Nz, Na, ots, sts, Qza);
    matrix<double> Qzz(Nz, Nz);        matrix_to_partition(Q, Nz, Nz, ots, ots, Qzz);
    matrix<double> Iz(Nz, Nz);         fill_identity(Iz);
    matrix<double> Az(Nz, Nz);         expm(Qzz, tdead, Az);
    matrix<double> Qzi;                invm(Qzz, Qzi);
    matrix<double> missed_returning;   missed_returning.resize(Na, Na);
    noalias(missed_returning) = prod(matrix<double>(prod(matrix<double>(prod(Qaz, matrix<double>(Iz - Az))), Qzi)), Qza);
    matrix<double> eQaa(Na, Na);
    noalias(eQaa) = Qaa - missed_returning;
    matrix_from_partition(eQ, Na, Na, sts, sts, eQaa);
    
    for (set<int>::iterator cls_j = classet.begin(); cls_j != classet.end(); ++cls_j) {
      if ( *cls_j != *cls_i ) {
	std::vector<int> bees, cees;
	for (int i=0; i<n; ++i)
	  if ( *cls_j == clazz[i] )
	    bees.push_back(i);
	  else if ( *cls_i != clazz[i] )
	    cees.push_back(i);
	int Nb = (int) bees.size();
	int Nc = (int) cees.size();
	int *bs = Nb ? &(bees[0]) : NULL;
	int *cs = Nc ? &(cees[0]) : NULL;
	matrix<double> Qab(Na, Nb);    matrix_to_partition(Q, Na, Nb, sts, bs, Qab);
	matrix<double> Qbb(Nb, Nb);    matrix_to_partition(Q, Nb, Nb, bs, bs, Qbb);
	matrix<double> escaped_indirect(Na, Nb);
	if ( Nc ) {
	  matrix<double> Qac(Na, Nc);  matrix_to_partition(Q, Na, Nc, sts, cs, Qac);
	  matrix<double> Qcc(Nc, Nc);  matrix_to_partition(Q, Nc, Nc, cs, cs, Qcc);
	  matrix<double> Qcb(Nc, Nb);  matrix_to_partition(Q, Nc, Nb, cs, bs, Qcb);
	  matrix<double> Ic(Nc, Nc);   fill_identity(Ic);
	  matrix<double> Ac(Nc, Nc);   expm(Qcc, tdead, Ac);
	  matrix<double> Qci(Nc, Nc);  invm(Qcc, Qci);
	  noalias(escaped_indirect)  = prod( matrix<double>(prod( matrix<double>(prod(Qac, Ic - Ac)), Qci )), Qcb);
	}
	else {
	  // no third class? can't escape indirectly
	  fill_const(escaped_indirect, 0.0);
	}
	matrix<double> Ab(Nb, Nb);     expm(Qbb, tdead, Ab);
	matrix<double> eQab(Na, Nb);
	noalias(eQab) = prod(Qab - escaped_indirect, Ab);
	matrix_from_partition(eQ, Na, Nb, sts, bs, eQab);
      }
    }

  }
}

/*

As discussed in (Milescu 2005), it can be preferable to start an experiment at equilibrium. Entry probabilities (P0) are then no longer constant, but a function of rate constants. We use the "direct" method given in Neher and Sakmann:

Defining $S = [Q | 1]$ and $u$ a row vector of ones, \[P_{eq} = u \cdot (S \cdot S^T)^{-1}\]

*/
QUBFAST_API void QtoPe(matrix<double>& Q, matrix<double>& Pe)
{
  int N = (int) Q.size1();
  if ( ((int) Pe.size1()) != 1 || ((int) Pe.size2()) != N )
    Pe.resize(1, N);
  matrix<double> S(N, N+1);
  matrix<double> U(1, N);
  for (int i=0; i<N; ++i) {
    for (int j=0; j<=N; ++j)
      S(i,j) = (j < N) ? Q(i,j) : 0.0;
    S(i,N) = 1.0;
    U(0,i) = 1.0;
  }
  matrix<double> SST(N, N);     noalias(SST) = prod(S, trans(S));
  matrix<double> SSTI(N, N);    invm(SST, SSTI);
  noalias(Pe) = prod(U, SSTI);
}



QUBFAST_API double alpha_logscale(vector<double>& alpha)
{
  double tot = 0.0;
  for (vector<double>::iterator ii = alpha.begin(); ii != alpha.end(); ++ii)
    tot += *ii;
  double scale = tot ? (1.0 / tot) : 1e10;
  for (vector<double>::iterator ii = alpha.begin(); ii != alpha.end(); ++ii)
    *ii *= scale;
  // std::cerr << scale << std::endl;
  return - log(scale);
}