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