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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#ifndef ___1_2_more_STATE_MODEL
#define ___1_2_more_STATE_MODEL
#include "definitions.h"
#include "replacementModel.h"
#include "fromQtoPt.h"
#include "errorMsg.h"
#include "matrixUtils.h"
class oneTwoMoreModel : public replacementModel {
public:
explicit oneTwoMoreModel(const MDOUBLE m1, const MDOUBLE m2,
const MDOUBLE m3, const MDOUBLE m4,const Vdouble &freq, bool useMarkovLimiting = true);
oneTwoMoreModel(const oneTwoMoreModel& other) {*this = other;}
virtual oneTwoMoreModel& operator=(const oneTwoMoreModel &other);
virtual oneTwoMoreModel* clone() const { return new oneTwoMoreModel(*this); }
virtual ~oneTwoMoreModel() {}
const int alphabetSize() const {return 3;} // two states and an intermediate (both states at once)
const MDOUBLE err_allow_for_pijt_function() const {return 1e-4;} // same as q2p definitions
const MDOUBLE Pij_t(const int i,const int j, const MDOUBLE d) const ;
const MDOUBLE dPij_dt(const int i,const int j, const MDOUBLE d) const{
if (d==0.0)
return _Q[i][j];
errorMsg::reportError("Error in oneTwoMoreModel, dPij_dt called");
return 0.0; // not supposed to be here
}
const MDOUBLE d2Pij_dt2(const int i,const int j, const MDOUBLE d) const{
errorMsg::reportError("Error in oneTwoMoreModel, d2Pij_dt2 called");
return 0.0; // not supposed to be here
}
const MDOUBLE freq(const int i) const {
if (i >= _freq.size())
errorMsg::reportError("Error in oneTwoMoreModel::freq, i > size of frequency vector");
return _freq[i];
}
const Vdouble getFreqs() const {return _freq;}
void setFreq(const Vdouble &freq);
void setMu1(const MDOUBLE val) ;
void setMu2(const MDOUBLE val) ;
void setMu3(const MDOUBLE val) ;
void setMu4(const MDOUBLE val) ;
const MDOUBLE getMu1() const {return _gain;}
const MDOUBLE getMu2() const {return _more;}
const MDOUBLE getMu3() const {return _less;}
const MDOUBLE getMu4() const {return _loss;}
void computeMarkovLimitingDistribution(); // compute P(infinity), which specifies the stationary distribution
MDOUBLE sumPijQij();
void norm(const MDOUBLE scale);
private:
virtual void updateQ();
void setEpsilonForZeroParams();
bool checkIsNullModel();
bool pijt_is_prob_value(MDOUBLE val) const;
bool areFreqsValid(Vdouble freq) const; // tests if frequencies are valid (>0, sum=1)
private:
MDOUBLE _gain; // _Q[0][1] not _Q[0][2]
MDOUBLE _more; // _Q[1][2]
MDOUBLE _less; // _Q[2][1] not _Q[2][0]
MDOUBLE _loss; // _Q[2][1]
VVdouble _Q;
Vdouble _freq;
bool _useMarkovLimiting; // should the markov limiting distribution be used to estimate the root frequencies
mutable bool _bQchanged; //indicates whether the Q matrix was changed after the last Pij_t call
mutable MDOUBLE _lastTcalculated;
mutable VVdouble _lastPtCalculated;
};
/*class gainLossModel : public replacementModel {
public:
explicit gainLossModel(const MDOUBLE m1, const MDOUBLE m2, const Vdouble freq);
virtual replacementModel* clone() const { return new gainLossModel(*this); }
gainLossModel(const gainLossModel& other): _q2pt(NULL) {*this = other;}
virtual gainLossModel& operator=(const gainLossModel &other);
virtual ~gainLossModel() {if (_q2pt) delete _q2pt;}
const int alphabetSize() const {return 3;} // two states and an intermediate (both states at once)
const MDOUBLE err_allow_for_pijt_function() const {return 1e-4;} // same as q2p definitions
const MDOUBLE Pij_t(const int i,const int j, const MDOUBLE d) const {
return _q2pt->Pij_t(i,j,d);
}
const MDOUBLE dPij_dt(const int i,const int j, const MDOUBLE d) const{
return _q2pt->dPij_dt(i,j,d);
}
const MDOUBLE d2Pij_dt2(const int i,const int j, const MDOUBLE d) const{
return _q2pt->d2Pij_dt2(i,j,d);
}
const MDOUBLE freq(const int i) const {
if (i >= _freq.size())
errorMsg::reportError("Error in gainLossModel::freq, i > size of frequency vector");
return _freq[i];
}
void setMu1(const MDOUBLE val, bool isReversible=true) { _gain = val; updateQ(isReversible);}
void setMu2(const MDOUBLE val,bool isReversible=true) { _more = val; updateQ(isReversible);}
const MDOUBLE getMu1() const {return _gain;}
const MDOUBLE getMu2() const {return _more;}
protected:
virtual void updateQ(bool isReversible=true);
virtual void normalizeQ();
protected:
Vdouble _freq;
MDOUBLE _gain;
MDOUBLE _more;
VVdouble _Q;
q2pt *_q2pt;
};
*/
/*
Q is a matrix of the following form:
0 1 01
0 1-m1 0 m1
1 0 1-m2 m2
01 (filled in assuming reversibility)
i.e. no direct change from state 0 to state 1 is allowed
*/
#endif // ___3STATE_MODEL
|
