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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
// $Id: computePijComponent.cpp 9253 2011-01-31 01:37:21Z rubi $
#include "definitions.h"
#include "treeIt.h"
#include "computePijComponent.h"
#include "logFile.h"
void computePijHomSpec::fillPij(const MDOUBLE dis, const stochasticProcess& sp, int derivationOrder, bool isReversible)
{
if (!(isReversible && sp.isReversible())) // if one is false
isReversible = false;
resize(sp.alphabetSize());
int i,j;
for (i=0; i<sp.alphabetSize(); i++) {
switch (derivationOrder) {
case 0:
_V[i][i] = sp.Pij_t(i,i,dis);
break;
case 1:
_V[i][i] = sp.dPij_dt(i,i,dis);
break;
case 2:
_V[i][i] = sp.d2Pij_dt2(i,i,dis);
break;
default:
errorMsg::reportError("error in function fillPij - derivationOrder must be 0, 1 or 2");
}
for (j=i+1; j<sp.alphabetSize(); j++) {
switch (derivationOrder) {
case 0:
_V[i][j] = sp.Pij_t(i,j,dis);
if ((_V[i][j] == 0 )&& (dis !=0)){
_V[i][j] = EPSILON;
}
break;
case 1:
_V[i][j] = sp.dPij_dt(i,j,dis);
break;
case 2:
_V[i][j] = sp.d2Pij_dt2(i,j,dis);
break;
default:
errorMsg::reportError("error in function fillPij - derivationOrder must be 0, 1 or 2");
}
if (sp.freq(j) == 0.0) {
if (isReversible) {
errorMsg::reportError("error in function fillPij");
}
}
// else {
if (isReversible){
_V[j][i] = _V[i][j]* sp.freq(i)/sp.freq(j);
}
else {
switch (derivationOrder) {
case 0:
_V[j][i] = sp.Pij_t(j,i,dis);
if ((_V[j][i] == 0 )&& (dis !=0))
_V[j][i] = EPSILON;
break;
case 1:
_V[j][i] = sp.dPij_dt(j,i,dis);
break;
case 2:
_V[j][i] = sp.d2Pij_dt2(j,i,dis);
break;
default:
errorMsg::reportError("error in function fillPij - derivationOrder must be 0, 1 or 2");
}
}
// }
}
}
}
void computePijHom::fillPij(const tree& et, const stochasticProcess& sp, int derivationOrder, bool isReversible) {
_V.resize(et.getNodesNum());
treeIterTopDownConst tIt(et);
tree::nodeP myNode = tIt.first();
{// skipping the root, but allocating place for the root pij even if they are not use
// to maintain that all arrays have the same size.
_V[myNode->id()].resize(sp.alphabetSize());
}
LOGDO(50,et.output(myLog::LogFile(),tree::ANCESTOR));
LOGDO(50,et.output(myLog::LogFile(),tree::PHYLIP));
for (; myNode != tIt.end(); myNode = tIt.next()) {
if (!(myNode->isRoot()))
_V[myNode->id()].fillPij(myNode->dis2father()*sp.getGlobalRate(),sp,derivationOrder,isReversible);
// else
// myLog::LogFile()<<"ROOT IS "<<myNode->name()<<endl;
}
}
void computePijGam::fillPij(const tree& et, const stochasticProcess& sp, int derivationOrder, bool isReversible) {
_V.resize(sp.categories());
for (int i=0; i < _V.size(); ++i) {
tree cp = et;
cp.multipleAllBranchesByFactor(sp.rates(i)/sp.getGlobalRate());// the global rate is taken care of in the hom pij.
_V[i].fillPij(cp,sp,derivationOrder,isReversible);
}
}
|
