1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
// $Id: likeDistProp.h 962 2006-11-07 15:13:34Z privmane $
#ifndef ___LIKE_DIST_PROP
#define ___LIKE_DIST_PROP
#include "definitions.h"
#include "countTableComponent.h"
#include "stochasticProcess.h"
#include <cmath>
class likeDistProp {
private:
const int _alphabetSize;
const vector<stochasticProcess>& _s1;
const MDOUBLE _toll;
public:
const MDOUBLE giveDistance( const vector<countTableComponentGam>& ctc,
MDOUBLE& resL) const;
explicit likeDistProp(const int alphabetSize,
const vector<stochasticProcess>& s1,
const MDOUBLE toll) : _alphabetSize(alphabetSize), _s1(s1) ,_toll(toll){
}
};
class C_evallikeDistProp_d{ // derivative.
public:
C_evallikeDistProp_d(const vector<countTableComponentGam>& ctc,
const vector<stochasticProcess>& inS1) : _ctc(ctc), _sp(inS1) {};
private:
const vector<countTableComponentGam>& _ctc;
const vector<stochasticProcess>& _sp;
public:
MDOUBLE operator() (MDOUBLE dist) {
MDOUBLE sumDL=0.0;
const MDOUBLE epsilonPIJ = 1e-10;
for (int gene=0; gene < _ctc.size(); ++ gene) {
for (int alph1=0; alph1 < _ctc[gene].alphabetSize(); ++alph1){
for (int alph2=0; alph2 < _ctc[gene].alphabetSize(); ++alph2){
for (int rateCategor = 0; rateCategor<_sp[gene].categories(); ++rateCategor) {
MDOUBLE rate = _sp[gene].rates(rateCategor);
MDOUBLE pij= _sp[gene].Pij_t(alph1,alph2,dist*rate);
MDOUBLE dpij = _sp[gene].dPij_dt(alph1,alph2,dist*rate);
if (pij<epsilonPIJ) {
pij = epsilonPIJ;
dpij = epsilonPIJ;
}
sumDL+= _ctc[gene].getCounts(alph1,alph2,rateCategor)*dpij*_sp[gene].ratesProb(rateCategor)
*rate/pij;
}
}
}
}
return -sumDL;
}
};
class C_evallikeDistProp{
private:
const vector<countTableComponentGam>& _ctc;
const vector<stochasticProcess>& _sp;
public:
C_evallikeDistProp(const vector<countTableComponentGam>& ctc,
const vector<stochasticProcess>& inS1):_ctc(ctc), _sp(inS1) {};
MDOUBLE operator() (MDOUBLE dist) {
const MDOUBLE epsilonPIJ = 1e-10;
MDOUBLE sumL=0.0;
for (int gene=0; gene < _ctc.size(); ++ gene) {
for (int alph1=0; alph1 < _ctc[gene].alphabetSize(); ++alph1){
for (int alph2=0; alph2 < _ctc[gene].alphabetSize(); ++alph2){
for (int rateCategor = 0; rateCategor<_sp[gene].categories(); ++rateCategor) {
MDOUBLE rate = _sp[gene].rates(rateCategor);
MDOUBLE pij= _sp[gene].Pij_t(alph1,alph2,dist*rate);
if (pij<0) {
pij = epsilonPIJ;
}
sumL += _ctc[gene].getCounts(alph1,alph2,rateCategor)*(log(pij)-log(_sp[gene].freq(alph2)))*_sp[gene].ratesProb(rateCategor);
}
}
}
}
return -sumL;
}
};
#endif
|
