bpp-phyl3  3.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
K80.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 "K80.h"
6 
7 // From the STL:
8 #include <cmath>
9 
10 using namespace bpp;
11 
13 
14 using namespace std;
15 
16 /******************************************************************************/
17 
18 K80::K80(shared_ptr<const NucleicAlphabet> alpha, double kappa) :
20  AbstractReversibleNucleotideSubstitutionModel(alpha, make_shared<CanonicalStateMap>(alpha, false), "K80."),
21  kappa_(kappa), r_(), l_(), k_(), exp1_(), exp2_(), p_(size_, size_)
22 {
23  addParameter_(new Parameter("K80.kappa", kappa, Parameter::R_PLUS_STAR));
25 }
26 
27 /******************************************************************************/
28 
30 {
31  kappa_ = getParameterValue("kappa");
32  k_ = (kappa_ + 1.) / 2.;
33  r_ = isScalable() ? 4. / (kappa_ + 2.) : 4;
34 
35  // Frequencies:
36  freq_[0] = freq_[1] = freq_[2] = freq_[3] = 1. / 4.;
37 
38  // Generator:
39  generator_(0, 0) = -2. - kappa_;
40  generator_(1, 1) = -2. - kappa_;
41  generator_(2, 2) = -2. - kappa_;
42  generator_(3, 3) = -2. - kappa_;
43 
44  generator_(0, 1) = 1.;
45  generator_(0, 3) = 1.;
46  generator_(1, 0) = 1.;
47  generator_(1, 2) = 1.;
48  generator_(2, 1) = 1.;
49  generator_(2, 3) = 1.;
50  generator_(3, 0) = 1.;
51  generator_(3, 2) = 1.;
52 
53  generator_(0, 2) = kappa_;
54  generator_(1, 3) = kappa_;
55  generator_(2, 0) = kappa_;
56  generator_(3, 1) = kappa_;
57 
58  // Normalization:
59  setScale(r_ / 4);
60 
61  // Exchangeability:
64 
65  // Eigen values:
66  eigenValues_[0] = 0;
67  eigenValues_[1] = -r_ * (1. + kappa_) / 2;
68  eigenValues_[2] = -r_ * (1. + kappa_) / 2;
69  eigenValues_[3] = -r_;
70 
71  // Eigen vectors:
72  leftEigenVectors_(0, 0) = 1. / 4.;
73  leftEigenVectors_(0, 1) = 1. / 4.;
74  leftEigenVectors_(0, 2) = 1. / 4.;
75  leftEigenVectors_(0, 3) = 1. / 4.;
76  leftEigenVectors_(1, 0) = 0.;
77  leftEigenVectors_(1, 1) = 1. / 2.;
78  leftEigenVectors_(1, 2) = 0.;
79  leftEigenVectors_(1, 3) = -1. / 2.;
80  leftEigenVectors_(2, 0) = 1. / 2.;
81  leftEigenVectors_(2, 1) = 0.;
82  leftEigenVectors_(2, 2) = -1. / 2.;
83  leftEigenVectors_(2, 3) = 0.;
84  leftEigenVectors_(3, 0) = 1. / 4.;
85  leftEigenVectors_(3, 1) = -1. / 4.;
86  leftEigenVectors_(3, 2) = 1. / 4.;
87  leftEigenVectors_(3, 3) = -1. / 4.;
88 
89  rightEigenVectors_(0, 0) = 1.;
90  rightEigenVectors_(0, 1) = 0.;
91  rightEigenVectors_(0, 2) = 1.;
92  rightEigenVectors_(0, 3) = 1.;
93  rightEigenVectors_(1, 0) = 1.;
94  rightEigenVectors_(1, 1) = 1.;
95  rightEigenVectors_(1, 2) = 0.;
96  rightEigenVectors_(1, 3) = -1.;
97  rightEigenVectors_(2, 0) = 1.;
98  rightEigenVectors_(2, 1) = 0.;
99  rightEigenVectors_(2, 2) = -1.;
100  rightEigenVectors_(2, 3) = 1.;
101  rightEigenVectors_(3, 0) = 1.;
102  rightEigenVectors_(3, 1) = -1.;
103  rightEigenVectors_(3, 2) = 0;
104  rightEigenVectors_(3, 3) = -1.;
105 }
106 
107 /******************************************************************************/
108 
109 double K80::Pij_t(size_t i, size_t j, double d) const
110 {
111  l_ = rate_ * r_ * d;
112  exp1_ = exp(-l_);
113  exp2_ = exp(-k_ * l_);
114 
115  switch (i)
116  {
117  case 0: // A
118  switch (j)
119  {
120  case 0: return 0.25 * (1. + exp1_) + 0.5 * exp2_; // A
121  case 1: return 0.25 * (1. - exp1_); // C
122  case 2: return 0.25 * (1. + exp1_) - 0.5 * exp2_; // G
123  case 3: return 0.25 * (1. - exp1_); // T, U
124  default: return 0;
125  }
126  case 1: // C
127  switch (j)
128  {
129  case 0: return 0.25 * (1. - exp1_); // A
130  case 1: return 0.25 * (1. + exp1_) + 0.5 * exp2_; // C
131  case 2: return 0.25 * (1. - exp1_); // G
132  case 3: return 0.25 * (1. + exp1_) - 0.5 * exp2_; // T, U
133  default: return 0;
134  }
135  case 2: // G
136  switch (j)
137  {
138  case 0: return 0.25 * (1. + exp1_) - 0.5 * exp2_; // A
139  case 1: return 0.25 * (1. - exp1_); // C
140  case 2: return 0.25 * (1. + exp1_) + 0.5 * exp2_; // G
141  case 3: return 0.25 * (1. - exp1_); // T, U
142  default: return 0;
143  }
144  case 3: // T, U
145  switch (j)
146  {
147  case 0: return 0.25 * (1. - exp1_); // A
148  case 1: return 0.25 * (1. + exp1_) - 0.5 * exp2_; // C
149  case 2: return 0.25 * (1. - exp1_); // G
150  case 3: return 0.25 * (1. + exp1_) + 0.5 * exp2_; // T, U
151  default: return 0;
152  }
153  default: return 0;
154  }
155 }
156 
157 /******************************************************************************/
158 
159 double K80::dPij_dt(size_t i, size_t j, double d) const
160 {
161  l_ = rate_ * r_ * d;
162  exp1_ = exp(-l_);
163  exp2_ = exp(-k_ * l_);
164 
165  switch (i)
166  {
167  case 0: // A
168  switch (j)
169  {
170  case 0: return rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // A
171  case 1: return rate_ * r_ / 4. * ( exp1_); // C
172  case 2: return rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // G
173  case 3: return rate_ * r_ / 4. * ( exp1_); // T, U
174  default: return 0;
175  }
176  case 1: // C
177  switch (j)
178  {
179  case 0: return rate_ * r_ / 4. * ( exp1_); // A
180  case 1: return rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // C
181  case 2: return rate_ * r_ / 4. * ( exp1_); // G
182  case 3: return rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // T, U
183  default: return 0;
184  }
185  case 2: // G
186  switch (j)
187  {
188  case 0: return rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // A
189  case 1: return rate_ * r_ / 4. * ( exp1_); // C
190  case 2: return rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // G
191  case 3: return rate_ * r_ / 4. * ( exp1_); // T, U
192  default: return 0;
193  }
194  case 3: // T, U
195  switch (j)
196  {
197  case 0: return rate_ * r_ / 4. * ( exp1_); // A
198  case 1: return rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // C
199  case 2: return rate_ * r_ / 4. * ( exp1_); // G
200  case 3: return rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // T, U
201  default: return 0;
202  }
203  default: return 0;
204  }
205 }
206 
207 /******************************************************************************/
208 
209 double K80::d2Pij_dt2(size_t i, size_t j, double d) const
210 {
211  double k_2 = k_ * k_;
212  double r_2 = rate_ * rate_ * r_ * r_;
213  l_ = rate_ * r_ * d;
214  exp1_ = exp(-l_);
215  exp2_ = exp(-k_ * l_);
216 
217  switch (i)
218  {
219  case 0: // A
220  switch (j)
221  {
222  case 0: return r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // A
223  case 1: return r_2 / 4. * (-exp1_); // C
224  case 2: return r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // G
225  case 3: return r_2 / 4. * (-exp1_); // T, U
226  default: return 0;
227  }
228  case 1: // C
229  switch (j)
230  {
231  case 0: return r_2 / 4. * (-exp1_); // A
232  case 1: return r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // C
233  case 2: return r_2 / 4. * (-exp1_); // G
234  case 3: return r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // T, U
235  default: return 0;
236  }
237  case 2: // G
238  switch (j)
239  {
240  case 0: return r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // A
241  case 1: return r_2 / 4. * (-exp1_); // C
242  case 2: return r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // G
243  case 3: return r_2 / 4. * (-exp1_); // T, U
244  default: return 0;
245  }
246  case 3: // T, U
247  switch (j)
248  {
249  case 0: return r_2 / 4. * (-exp1_); // A
250  case 1: return r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // C
251  case 2: return r_2 / 4. * (-exp1_); // G
252  case 3: return r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // T, U
253  default: return 0;
254  }
255  default: return 0;
256  }
257  return 0;
258 }
259 
260 /******************************************************************************/
261 
262 const Matrix<double>& K80::getPij_t(double d) const
263 {
264  l_ = rate_ * r_ * d;
265  exp1_ = exp(-l_);
266  exp2_ = exp(-k_ * l_);
267 
268  // A
269  p_(0, 0) = 0.25 * (1. + exp1_) + 0.5 * exp2_; // A
270  p_(0, 1) = 0.25 * (1. - exp1_); // C
271  p_(0, 2) = 0.25 * (1. + exp1_) - 0.5 * exp2_; // G
272  p_(0, 3) = 0.25 * (1. - exp1_); // T, U
273 
274  // C
275  p_(1, 0) = 0.25 * (1. - exp1_); // A
276  p_(1, 1) = 0.25 * (1. + exp1_) + 0.5 * exp2_; // C
277  p_(1, 2) = 0.25 * (1. - exp1_); // G
278  p_(1, 3) = 0.25 * (1. + exp1_) - 0.5 * exp2_; // T, U
279 
280  // G
281  p_(2, 0) = 0.25 * (1. + exp1_) - 0.5 * exp2_; // A
282  p_(2, 1) = 0.25 * (1. - exp1_); // C
283  p_(2, 2) = 0.25 * (1. + exp1_) + 0.5 * exp2_; // G
284  p_(2, 3) = 0.25 * (1. - exp1_); // T, U
285 
286  // T, U
287  p_(3, 0) = 0.25 * (1. - exp1_); // A
288  p_(3, 1) = 0.25 * (1. + exp1_) - 0.5 * exp2_; // C
289  p_(3, 2) = 0.25 * (1. - exp1_); // G
290  p_(3, 3) = 0.25 * (1. + exp1_) + 0.5 * exp2_; // T, U
291 
292  return p_;
293 }
294 
295 const Matrix<double>& K80::getdPij_dt(double d) const
296 {
297  l_ = rate_ * r_ * d;
298  exp1_ = exp(-l_);
299  exp2_ = exp(-k_ * l_);
300 
301  p_(0, 0) = rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // A
302  p_(0, 1) = rate_ * r_ / 4. * ( exp1_); // C
303  p_(0, 2) = rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // G
304  p_(0, 3) = rate_ * r_ / 4. * ( exp1_); // T, U
305 
306  // C
307  p_(1, 0) = rate_ * r_ / 4. * ( exp1_); // A
308  p_(1, 1) = rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // C
309  p_(1, 2) = rate_ * r_ / 4. * ( exp1_); // G
310  p_(1, 3) = rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // T, U
311 
312  // G
313  p_(2, 0) = rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // A
314  p_(2, 1) = rate_ * r_ / 4. * ( exp1_); // C
315  p_(2, 2) = rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // G
316  p_(2, 3) = rate_ * r_ / 4. * ( exp1_); // T, U
317 
318  // T, U
319  p_(3, 0) = rate_ * r_ / 4. * ( exp1_); // A
320  p_(3, 1) = rate_ * r_ / 4. * (-exp1_ + 2. * k_ * exp2_); // C
321  p_(3, 2) = rate_ * r_ / 4. * ( exp1_); // G
322  p_(3, 3) = rate_ * r_ / 4. * (-exp1_ - 2. * k_ * exp2_); // T, U
323 
324  return p_;
325 }
326 
327 const Matrix<double>& K80::getd2Pij_dt2(double d) const
328 {
329  double k_2 = k_ * k_;
330  double r_2 = rate_ * rate_ * r_ * r_;
331  l_ = rate_ * r_ * d;
332  exp1_ = exp(-l_);
333  exp2_ = exp(-k_ * l_);
334 
335  p_(0, 0) = r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // A
336  p_(0, 1) = r_2 / 4. * (-exp1_); // C
337  p_(0, 2) = r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // G
338  p_(0, 3) = r_2 / 4. * (-exp1_); // T, U
339 
340  // C
341  p_(1, 0) = r_2 / 4. * (-exp1_); // A
342  p_(1, 1) = r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // C
343  p_(1, 2) = r_2 / 4. * (-exp1_); // G
344  p_(1, 3) = r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // T, U
345 
346  // G
347  p_(2, 0) = r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // A
348  p_(2, 1) = r_2 / 4. * (-exp1_); // C
349  p_(2, 2) = r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // G
350  p_(2, 3) = r_2 / 4. * (-exp1_); // T, U
351 
352  // T, U
353  p_(3, 0) = r_2 / 4. * (-exp1_); // A
354  p_(3, 1) = r_2 / 4. * ( exp1_ - 2. * k_2 * exp2_); // C
355  p_(3, 2) = r_2 / 4. * (-exp1_); // G
356  p_(3, 3) = r_2 / 4. * ( exp1_ + 2. * k_2 * exp2_); // T, U
357 
358  return p_;
359 }
360 
361 /******************************************************************************/
void addParameter_(Parameter *parameter)
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
double k_
Definition: K80.h:119
RowMatrix< double > p_
Definition: K80.h:120
const Matrix< double > & getd2Pij_dt2(double d) const override
Definition: K80.cpp:327
double exp2_
Definition: K80.h:119
double kappa_
Definition: K80.h:118
K80(std::shared_ptr< const NucleicAlphabet > alpha, double kappa=1.)
Definition: K80.cpp:18
double d2Pij_dt2(size_t i, size_t j, double d) const override
Definition: K80.cpp:209
double dPij_dt(size_t i, size_t j, double d) const override
Definition: K80.cpp:159
double Pij_t(size_t i, size_t j, double d) const override
Definition: K80.cpp:109
double l_
Definition: K80.h:119
double exp1_
Definition: K80.h:119
double r_
Definition: K80.h:118
const Matrix< double > & getPij_t(double d) const override
Definition: K80.cpp:262
void updateMatrices_() override
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: K80.cpp:29
const Matrix< double > & getdPij_dt(double d) const override
Definition: K80.cpp:295
static void scale(Matrix &A, Scalar a, Scalar b=0)
static const std::shared_ptr< IntervalConstraint > R_PLUS_STAR
Defines the basic types of data flow nodes.
ExtendedFloat exp(const ExtendedFloat &ef)