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 
10 using namespace bpp;
11 
12 // From the STL:
13 #include <cmath>
14 
15 using namespace std;
16 
17 /******************************************************************************/
18 
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_;
48  piG_ = 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 
152 double 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 
205 double 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 
258 double 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 
312 const 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 
344 const 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 
376 const 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 
411 void 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
F81(std::shared_ptr< const NucleicAlphabet > alpha, double piA=0.25, double piC=0.25, double piG=0.25, double piT=0.25)
Definition: F81.cpp:19
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)