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