bpp-phyl3 3.0.0
F81.cpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: The Bio++ Development Group
2//
3// SPDX-License-Identifier: CECILL-2.1
4
6
7#include "../FrequencySet/NucleotideFrequencySet.h"
8#include "F81.h"
9
10using namespace bpp;
11
12// From the STL:
13#include <cmath>
14
15using namespace std;
16
17/******************************************************************************/
18
19F81::F81(
20 shared_ptr<const NucleicAlphabet> alpha,
21 double piA,
22 double piC,
23 double piG,
24 double piT) :
26 AbstractReversibleNucleotideSubstitutionModel(alpha, make_shared<CanonicalStateMap>(alpha, false), "F81."),
27 r_(),
28 piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
29 theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
30 exp1_(), l_(), p_(size_, size_)
31{
35 computeFrequencies(false);
37}
38
39/******************************************************************************/
40
42{
43 theta_ = getParameterValue("theta");
44 theta1_ = getParameterValue("theta1");
45 theta2_ = getParameterValue("theta2");
46 piA_ = theta1_ * (1. - theta_);
47 piC_ = (1. - theta2_) * theta_;
49 piT_ = (1. - theta1_) * (1. - theta_);
50 piR_ = piA_ + piG_;
51 piY_ = piT_ + piC_;
52
53 freq_[0] = piA_;
54 freq_[1] = piC_;
55 freq_[2] = piG_;
56 freq_[3] = piT_;
57
58 generator_(0, 0) = piA_ - 1;
59 generator_(1, 1) = piC_ - 1;
60 generator_(2, 2) = piG_ - 1;
61 generator_(3, 3) = piT_ - 1;
62
63 generator_(1, 0) = piA_;
64 generator_(3, 0) = piA_;
65 generator_(0, 1) = piC_;
66 generator_(2, 1) = piC_;
67 generator_(1, 2) = piG_;
68 generator_(3, 2) = piG_;
69 generator_(0, 3) = piT_;
70 generator_(2, 3) = piT_;
71
72 generator_(2, 0) = piA_;
73 generator_(3, 1) = piC_;
74 generator_(0, 2) = piG_;
75 generator_(1, 3) = piT_;
76
77 // Normalization:
78 r_ = 1. / (2. * (piA_ * piC_ + piC_ * piG_ + piA_ * piT_ + piG_ * piT_ + (piC_ * piT_ + piA_ * piG_)));
80
81 // Exchangeability:
82 exchangeability_(0, 0) = generator_(0, 0) / piA_;
83 exchangeability_(0, 1) = generator_(0, 1) / piC_;
84 exchangeability_(0, 2) = generator_(0, 2) / piG_;
85 exchangeability_(0, 3) = generator_(0, 3) / piT_;
86
87 exchangeability_(1, 0) = generator_(1, 0) / piA_;
88 exchangeability_(1, 1) = generator_(1, 1) / piC_;
89 exchangeability_(1, 2) = generator_(1, 2) / piG_;
90 exchangeability_(1, 3) = generator_(1, 3) / piT_;
91
92 exchangeability_(2, 0) = generator_(2, 0) / piA_;
93 exchangeability_(2, 1) = generator_(2, 1) / piC_;
94 exchangeability_(2, 2) = generator_(2, 2) / piG_;
95 exchangeability_(2, 3) = generator_(2, 3) / piT_;
96
97 exchangeability_(3, 0) = generator_(3, 0) / piA_;
98 exchangeability_(3, 1) = generator_(3, 1) / piC_;
99 exchangeability_(3, 2) = generator_(3, 2) / piG_;
100 exchangeability_(3, 3) = generator_(3, 3) / piT_;
101
102 // Eigen values:
103 eigenValues_[0] = 0;
104 eigenValues_[1] = -r_;
105 eigenValues_[2] = -r_;
106 eigenValues_[3] = -r_;
107
108 // Eigen vectors:
109 leftEigenVectors_(0, 0) = piA_;
110 leftEigenVectors_(0, 1) = piC_;
111 leftEigenVectors_(0, 2) = piG_;
112 leftEigenVectors_(0, 3) = piT_;
113
114 leftEigenVectors_(1, 0) = 0.;
115 leftEigenVectors_(1, 1) = piT_ / piY_;
116 leftEigenVectors_(1, 2) = 0.;
117 leftEigenVectors_(1, 3) = -piT_ / piY_;
118
119 leftEigenVectors_(2, 0) = piG_ / piR_;
120 leftEigenVectors_(2, 1) = 0.;
121 leftEigenVectors_(2, 2) = -piG_ / piR_;
122 leftEigenVectors_(2, 3) = 0.;
123
124 leftEigenVectors_(3, 0) = piA_ * piY_ / piR_;
125 leftEigenVectors_(3, 1) = -piC_;
126 leftEigenVectors_(3, 2) = piG_ * piY_ / piR_;
127 leftEigenVectors_(3, 3) = -piT_;
128
129 rightEigenVectors_(0, 0) = 1.;
130 rightEigenVectors_(0, 1) = 0.;
131 rightEigenVectors_(0, 2) = 1.;
132 rightEigenVectors_(0, 3) = 1.;
133
134 rightEigenVectors_(1, 0) = 1.;
135 rightEigenVectors_(1, 1) = 1.;
136 rightEigenVectors_(1, 2) = 0.;
137 rightEigenVectors_(1, 3) = -piR_ / piY_;
138
139 rightEigenVectors_(2, 0) = 1.;
140 rightEigenVectors_(2, 1) = 0.;
141 rightEigenVectors_(2, 2) = -piA_ / piG_;
142 rightEigenVectors_(2, 3) = 1.;
143
144 rightEigenVectors_(3, 0) = 1.;
145 rightEigenVectors_(3, 1) = -piC_ / piT_;
146 rightEigenVectors_(3, 2) = 0.;
147 rightEigenVectors_(3, 3) = -piR_ / piY_;
148}
149
150/******************************************************************************/
151
152double F81::Pij_t(size_t i, size_t j, double d) const
153{
154 l_ = rate_ * r_ * d;
155 exp1_ = exp(-l_);
156
157 switch (i)
158 {
159 // A
160 case 0: {
161 switch (j)
162 {
163 case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp1_; // A
164 case 1: return piC_ * (1. - exp1_); // C
165 case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp1_; // G
166 case 3: return piT_ * (1. - exp1_); // T, U
167 }
168 }
169 // C
170 case 1: {
171 switch (j)
172 {
173 case 0: return piA_ * (1. - exp1_); // A
174 case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp1_; // C
175 case 2: return piG_ * (1. - exp1_); // G
176 case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp1_; // T, U
177 }
178 }
179 // G
180 case 2: {
181 switch (j)
182 {
183 case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp1_; // A
184 case 1: return piC_ * (1. - exp1_); // C
185 case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp1_; // G
186 case 3: return piT_ * (1. - exp1_); // T, U
187 }
188 }
189 // T, U
190 case 3: {
191 switch (j)
192 {
193 case 0: return piA_ * (1. - exp1_); // A
194 case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp1_; // C
195 case 2: return piG_ * (1. - exp1_); // G
196 case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp1_; // T, U
197 }
198 }
199 }
200 return 0;
201}
202
203/******************************************************************************/
204
205double F81::dPij_dt(size_t i, size_t j, double d) const
206{
207 l_ = rate_ * r_ * d;
208 exp1_ = exp(-l_);
209
210 switch (i)
211 {
212 // A
213 case 0: {
214 switch (j)
215 {
216 case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * exp1_); // A
217 case 1: return rate_ * r_ * (piC_ * exp1_); // C
218 case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * exp1_); // G
219 case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
220 }
221 }
222 // C
223 case 1: {
224 switch (j)
225 {
226 case 0: return rate_ * r_ * (piA_ * exp1_); // A
227 case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * exp1_); // C
228 case 2: return rate_ * r_ * (piG_ * exp1_); // G
229 case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * exp1_); // T, U
230 }
231 }
232 // G
233 case 2: {
234 switch (j)
235 {
236 case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * exp1_); // A
237 case 1: return rate_ * r_ * (piC_ * exp1_); // C
238 case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * exp1_); // G
239 case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
240 }
241 }
242 // T, U
243 case 3: {
244 switch (j)
245 {
246 case 0: return rate_ * r_ * (piA_ * exp1_); // A
247 case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * exp1_); // C
248 case 2: return rate_ * r_ * (piG_ * exp1_); // G
249 case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * exp1_); // T, U
250 }
251 }
252 }
253 return 0;
254}
255
256/******************************************************************************/
257
258double F81::d2Pij_dt2(size_t i, size_t j, double d) const
259{
260 double r_2 = rate_ * rate_ * r_ * r_;
261 l_ = rate_ * r_ * d;
262 exp1_ = exp(-l_);
263
264 switch (i)
265 {
266 // A
267 case 0: {
268 switch (j)
269 {
270 case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * exp1_); // A
271 case 1: return r_2 * (piC_ * -exp1_); // C
272 case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * exp1_); // G
273 case 3: return r_2 * (piT_ * -exp1_); // T, U
274 }
275 }
276 // C
277 case 1: {
278 switch (j)
279 {
280 case 0: return r_2 * (piA_ * -exp1_); // A
281 case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * exp1_); // C
282 case 2: return r_2 * (piG_ * -exp1_); // G
283 case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * exp1_); // T, U
284 }
285 }
286 // G
287 case 2: {
288 switch (j)
289 {
290 case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * exp1_); // A
291 case 1: return r_2 * (piC_ * -exp1_); // C
292 case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * exp1_); // G
293 case 3: return r_2 * (piT_ * -exp1_); // T, U
294 }
295 }
296 // T, U
297 case 3: {
298 switch (j)
299 {
300 case 0: return r_2 * (piA_ * -exp1_); // A
301 case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * exp1_); // C
302 case 2: return r_2 * (piG_ * -exp1_); // G
303 case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * exp1_); // T, U
304 }
305 }
306 }
307 return 0;
308}
309
310/******************************************************************************/
311
312const Matrix<double>& F81::getPij_t(double d) const
313{
314 l_ = rate_ * r_ * d;
315 exp1_ = exp(-l_);
316
317 // A
318 p_(0, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp1_; // A
319 p_(0, 1) = piC_ * (1. - exp1_); // C
320 p_(0, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp1_; // G
321 p_(0, 3) = piT_ * (1. - exp1_); // T, U
322
323 // C
324 p_(1, 0) = piA_ * (1. - exp1_); // A
325 p_(1, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp1_; // C
326 p_(1, 2) = piG_ * (1. - exp1_); // G
327 p_(1, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp1_; // T, U
328
329 // G
330 p_(2, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp1_; // A
331 p_(2, 1) = piC_ * (1. - exp1_); // C
332 p_(2, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp1_; // G
333 p_(2, 3) = piT_ * (1. - exp1_); // T, U
334
335 // T, U
336 p_(3, 0) = piA_ * (1. - exp1_); // A
337 p_(3, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp1_; // C
338 p_(3, 2) = piG_ * (1. - exp1_); // G
339 p_(3, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp1_; // T, U
340
341 return p_;
342}
343
344const Matrix<double>& F81::getdPij_dt(double d) const
345{
346 l_ = rate_ * r_ * d;
347 exp1_ = exp(-l_);
348
349 // A
350 p_(0, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * exp1_); // A
351 p_(0, 1) = rate_ * r_ * (piC_ * exp1_); // C
352 p_(0, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * exp1_); // G
353 p_(0, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
354
355 // C
356 p_(1, 0) = rate_ * r_ * (piA_ * exp1_); // A
357 p_(1, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * exp1_); // C
358 p_(1, 2) = rate_ * r_ * (piG_ * exp1_); // G
359 p_(1, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * exp1_); // T, U
360
361 // G
362 p_(2, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * exp1_); // A
363 p_(2, 1) = rate_ * r_ * (piC_ * exp1_); // C
364 p_(2, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * exp1_); // G
365 p_(2, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
366
367 // T, U
368 p_(3, 0) = rate_ * r_ * (piA_ * exp1_); // A
369 p_(3, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * exp1_); // C
370 p_(3, 2) = rate_ * r_ * (piG_ * exp1_); // G
371 p_(3, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * exp1_); // T, U
372
373 return p_;
374}
375
376const Matrix<double>& F81::getd2Pij_dt2(double d) const
377{
378 double r_2 = rate_ * rate_ * r_ * r_;
379 l_ = rate_ * r_ * d;
380 exp1_ = exp(-l_);
381
382 // A
383 p_(0, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * exp1_); // A
384 p_(0, 1) = r_2 * (piC_ * -exp1_); // C
385 p_(0, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * exp1_); // G
386 p_(0, 3) = r_2 * (piT_ * -exp1_); // T, U
387
388 // C
389 p_(1, 0) = r_2 * (piA_ * -exp1_); // A
390 p_(1, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * exp1_); // C
391 p_(1, 2) = r_2 * (piG_ * -exp1_); // G
392 p_(1, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * exp1_); // T, U
393
394 // G
395 p_(2, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * exp1_); // A
396 p_(2, 1) = r_2 * (piC_ * -exp1_); // C
397 p_(2, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * exp1_); // G
398 p_(2, 3) = r_2 * (piT_ * -exp1_); // T, U
399
400 // T, U
401 p_(3, 0) = r_2 * (piA_ * -exp1_); // A
402 p_(3, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * exp1_); // C
403 p_(3, 2) = r_2 * (piG_ * -exp1_); // G
404 p_(3, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * exp1_); // T, U
405
406 return p_;
407}
408
409/******************************************************************************/
410
411void F81::setFreq(std::map<int, double>& freqs)
412{
413 piA_ = freqs[0];
414 piC_ = freqs[1];
415 piG_ = freqs[2];
416 piT_ = freqs[3];
417 vector<string> thetas(3);
418 thetas[0] = getNamespace() + "theta";
419 thetas[1] = getNamespace() + "theta1";
420 thetas[2] = getNamespace() + "theta2";
422 pl[0].setValue(piC_ + piG_);
423 pl[1].setValue(piA_ / (piA_ + piT_));
424 pl[2].setValue(piG_ / (piC_ + piG_));
426}
427
428/******************************************************************************/
void addParameter_(Parameter *parameter)
void setParametersValues(const ParameterList &parameters) override
std::string getNamespace() const override
const ParameterList & getParameters() const override
double getParameterValue(const std::string &name) const override
Specialisation abstract class for reversible nucleotide substitution model.
RowMatrix< double > generator_
The generator matrix of the model.
Vdouble eigenValues_
The vector of eigen values.
RowMatrix< double > exchangeability_
The exchangeability matrix of the model, defined as . When the model is reversible,...
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
Vdouble freq_
The vector of equilibrium frequencies.
double rate_
The rate of the model (default: 1). The generator (and all its vectorial components) is independent o...
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:168
void updateMatrices_() override
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: F81.cpp:41
void setFreq(std::map< int, double > &freqs) override
This method is redefined to actualize the corresponding parameters piA, piT, piG and piC too.
Definition: F81.cpp:411
double theta1_
Definition: F81.h:147
double piT_
Definition: F81.h:147
double exp1_
Definition: F81.h:148
const Matrix< double > & getdPij_dt(double d) const override
Definition: F81.cpp:344
const Matrix< double > & getPij_t(double d) const override
Definition: F81.cpp:312
double dPij_dt(size_t i, size_t j, double d) const override
Definition: F81.cpp:205
double r_
Definition: F81.h:147
double piG_
Definition: F81.h:147
const Matrix< double > & getd2Pij_dt2(double d) const override
Definition: F81.cpp:376
double d2Pij_dt2(size_t i, size_t j, double d) const override
Definition: F81.cpp:258
double Pij_t(size_t i, size_t j, double d) const override
Definition: F81.cpp:152
double l_
Definition: F81.h:148
double piC_
Definition: F81.h:147
double theta_
Definition: F81.h:147
double piA_
Definition: F81.h:147
double theta2_
Definition: F81.h:147
RowMatrix< double > p_
Definition: F81.h:149
double piR_
Definition: F81.h:147
double piY_
Definition: F81.h:147
static std::shared_ptr< IntervalConstraint > FREQUENCE_CONSTRAINT_SMALL
Definition: FrequencySet.h:90
static void scale(Matrix &A, Scalar a, Scalar b=0)
virtual ParameterList createSubList(const std::vector< std::string > &names) const
Defines the basic types of data flow nodes.
ExtendedFloat exp(const ExtendedFloat &ef)