bpp-phyl3  3.0.0
HKY85.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 "HKY85.h"
9 
10 using namespace bpp;
11 
12 // From the STL:
13 #include <cmath>
14 
15 using namespace std;
16 
17 /******************************************************************************/
18 
19 HKY85::HKY85(
20  shared_ptr<const NucleicAlphabet> alpha,
21  double kappa,
22  double piA,
23  double piC,
24  double piG,
25  double piT) :
26  AbstractParameterAliasable("HKY85."),
27  AbstractReversibleNucleotideSubstitutionModel(alpha, make_shared<CanonicalStateMap>(alpha, false), "HKY85."),
28  kappa_(kappa), k1_(), k2_(), r_(),
29  piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
30  theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
31  exp1_(), exp21_(), exp22_(), l_(), p_(size_, size_)
32 {
33  addParameter_(new Parameter("HKY85.kappa", kappa, Parameter::R_PLUS_STAR));
34  addParameter_(new Parameter("HKY85.theta", theta_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
35  addParameter_(new Parameter("HKY85.theta1", theta1_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
36  addParameter_(new Parameter("HKY85.theta2", theta2_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
37  updateMatrices_();
38 }
39 
40 /******************************************************************************/
41 
42 void HKY85::updateMatrices_()
43 {
44  kappa_ = getParameterValue("kappa");
45  theta_ = getParameterValue("theta");
46  theta1_ = getParameterValue("theta1");
47  theta2_ = getParameterValue("theta2");
48  piA_ = theta1_ * (1. - theta_);
49  piC_ = (1. - theta2_) * theta_;
50  piG_ = theta2_ * theta_;
51  piT_ = (1. - theta1_) * (1. - theta_);
52  piR_ = piA_ + piG_;
53  piY_ = piT_ + piC_;
54  k1_ = kappa_ * piY_ + piR_;
55  k2_ = kappa_ * piR_ + piY_;
56 
57  freq_[0] = piA_;
58  freq_[1] = piC_;
59  freq_[2] = piG_;
60  freq_[3] = piT_;
61 
62  generator_(0, 0) = -( piC_ + kappa_ * piG_ + piT_);
63  generator_(1, 1) = -( piA_ + piG_ + kappa_ * piT_);
64  generator_(2, 2) = -(kappa_ * piA_ + piC_ + piT_);
65  generator_(3, 3) = -( piA_ + kappa_ * piC_ + piG_ );
66 
67  generator_(1, 0) = piA_;
68  generator_(3, 0) = piA_;
69  generator_(0, 1) = piC_;
70  generator_(2, 1) = piC_;
71  generator_(1, 2) = piG_;
72  generator_(3, 2) = piG_;
73  generator_(0, 3) = piT_;
74  generator_(2, 3) = piT_;
75 
76  generator_(2, 0) = kappa_ * piA_;
77  generator_(3, 1) = kappa_ * piC_;
78  generator_(0, 2) = kappa_ * piG_;
79  generator_(1, 3) = kappa_ * piT_;
80 
81  // Normalization:
82  r_ = isScalable() ? 1. / (2. * (piA_ * piC_ + piC_ * piG_ + piA_ * piT_ + piG_ * piT_ + kappa_ * (piC_ * piT_ + piA_ * piG_))) : 1;
83 
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_ * (kappa_ * piY_ + piR_);
110  eigenValues_[2] = -r_ * (kappa_ * piR_ + piY_);
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 
157 double HKY85::Pij_t(size_t i, size_t j, double d) const
158 {
159  l_ = rate_ * r_ * d;
160  exp1_ = exp(-l_);
161  exp22_ = exp(-k2_ * l_);
162  exp21_ = exp(-k1_ * l_);
163 
164  switch (i)
165  {
166  case 0: // A
167  switch (j)
168  {
169  case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp22_; // A
170  case 1: return piC_ * (1. - exp1_); // C
171  case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp22_; // G
172  case 3: return piT_ * (1. - exp1_); // T, U
173  default: return 0;
174  }
175  case 1: // C
176  switch (j)
177  {
178  case 0: return piA_ * (1. - exp1_); // A
179  case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp21_; // C
180  case 2: return piG_ * (1. - exp1_); // G
181  case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp21_; // T, U
182  default: return 0;
183  }
184  case 2: // G
185  switch (j)
186  {
187  case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp22_; // A
188  case 1: return piC_ * (1. - exp1_); // C
189  case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp22_; // G
190  case 3: return piT_ * (1. - exp1_); // T, U
191  default: return 0;
192  }
193  case 3: // T, U
194  switch (j)
195  {
196  case 0: return piA_ * (1. - exp1_); // A
197  case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp21_; // C
198  case 2: return piG_ * (1. - exp1_); // G
199  case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp21_; // T, U
200  default: return 0;
201  }
202  default: return 0;
203  }
204 }
205 
206 /******************************************************************************/
207 
208 double HKY85::dPij_dt(size_t i, size_t j, double d) const
209 {
210  l_ = rate_ * r_ * d;
211  exp1_ = exp(-l_);
212  exp22_ = exp(-k2_ * l_);
213  exp21_ = exp(-k1_ * l_);
214 
215  switch (i)
216  {
217  case 0: // A
218  switch (j)
219  {
220  case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_ * exp22_); // A
221  case 1: return rate_ * r_ * (piC_ * exp1_); // C
222  case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_ * exp22_); // G
223  case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
224  default: return 0;
225  }
226  case 1: // C
227  switch (j)
228  {
229  case 0: return rate_ * r_ * (piA_ * exp1_); // A
230  case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_ * exp21_); // C
231  case 2: return rate_ * r_ * (piG_ * exp1_); // G
232  case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_ * exp21_); // T, U
233  default: return 0;
234  }
235  case 2: // G
236  switch (j)
237  {
238  case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_ * exp22_); // A
239  case 1: return rate_ * r_ * (piC_ * exp1_); // C
240  case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_ * exp22_); // G
241  case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
242  default: return 0;
243  }
244  case 3: // T, U
245  switch (j)
246  {
247  case 0: return rate_ * r_ * (piA_ * exp1_); // A
248  case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_ * exp21_); // C
249  case 2: return rate_ * r_ * (piG_ * exp1_); // G
250  case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_ * exp21_); // T, U
251  default: return 0;
252  }
253  default: return 0;
254  }
255 }
256 
257 /******************************************************************************/
258 
259 double HKY85::d2Pij_dt2(size_t i, size_t j, double d) const
260 {
261  double r_2 = rate_ * rate_ * r_ * r_;
262  l_ = rate_ * r_ * d;
263  double k1_2 = k1_ * k1_;
264  double k2_2 = k2_ * k2_;
265  exp1_ = exp(-l_);
266  exp22_ = exp(-k2_ * l_);
267  exp21_ = exp(-k1_ * l_);
268 
269  switch (i)
270  {
271  case 0: // A
272  switch (j)
273  {
274  case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_2 * exp22_); // A
275  case 1: return r_2 * (piC_ * -exp1_); // C
276  case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_2 * exp22_); // G
277  case 3: return r_2 * (piT_ * -exp1_); // T, U
278  default: return 0;
279  }
280  case 1: // C
281  switch (j)
282  {
283  case 0: return r_2 * (piA_ * -exp1_); // A
284  case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_2 * exp21_); // C
285  case 2: return r_2 * (piG_ * -exp1_); // G
286  case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_2 * exp21_); // T, U
287  default: return 0;
288  }
289  case 2: // G
290  switch (j)
291  {
292  case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_2 * exp22_); // A
293  case 1: return r_2 * (piC_ * -exp1_); // C
294  case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_2 * exp22_); // G
295  case 3: return r_2 * (piT_ * -exp1_); // T, U
296  default: return 0;
297  }
298  case 3: // T, U
299  switch (j)
300  {
301  case 0: return r_2 * (piA_ * -exp1_); // A
302  case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_2 * exp21_); // C
303  case 2: return r_2 * (piG_ * -exp1_); // G
304  case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_2 * exp21_); // T, U
305  default: return 0;
306  }
307  default: return 0;
308  }
309 }
310 
311 /******************************************************************************/
312 
313 const Matrix<double>& HKY85::getPij_t(double d) const
314 {
315  l_ = rate_ * r_ * d;
316  exp1_ = exp(-l_);
317  exp22_ = exp(-k2_ * l_);
318  exp21_ = exp(-k1_ * l_);
319 
320  // A
321  p_(0, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp22_; // A
322  p_(0, 1) = piC_ * (1. - exp1_); // C
323  p_(0, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp22_; // G
324  p_(0, 3) = piT_ * (1. - exp1_); // T, U
325 
326  // C
327  p_(1, 0) = piA_ * (1. - exp1_); // A
328  p_(1, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp21_; // C
329  p_(1, 2) = piG_ * (1. - exp1_); // G
330  p_(1, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp21_; // T, U
331 
332  // G
333  p_(2, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp22_; // A
334  p_(2, 1) = piC_ * (1. - exp1_); // C
335  p_(2, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp22_; // G
336  p_(2, 3) = piT_ * (1. - exp1_); // T, U
337 
338  // T, U
339  p_(3, 0) = piA_ * (1. - exp1_); // A
340  p_(3, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp21_; // C
341  p_(3, 2) = piG_ * (1. - exp1_); // G
342  p_(3, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp21_; // T, U
343 
344  return p_;
345 }
346 
347 const Matrix<double>& HKY85::getdPij_dt(double d) const
348 {
349  l_ = rate_ * r_ * d;
350  exp1_ = exp(-l_);
351  exp22_ = exp(-k2_ * l_);
352  exp21_ = exp(-k1_ * l_);
353 
354  // A
355  p_(0, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_ * exp22_); // A
356  p_(0, 1) = rate_ * r_ * (piC_ * exp1_); // C
357  p_(0, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_ * exp22_); // G
358  p_(0, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
359 
360  // C
361  p_(1, 0) = rate_ * r_ * (piA_ * exp1_); // A
362  p_(1, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_ * exp21_); // C
363  p_(1, 2) = rate_ * r_ * (piG_ * exp1_); // G
364  p_(1, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_ * exp21_); // T, U
365 
366  // G
367  p_(2, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_ * exp22_); // A
368  p_(2, 1) = rate_ * r_ * (piC_ * exp1_); // C
369  p_(2, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_ * exp22_); // G
370  p_(2, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
371 
372  // T, U
373  p_(3, 0) = rate_ * r_ * (piA_ * exp1_); // A
374  p_(3, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_ * exp21_); // C
375  p_(3, 2) = rate_ * r_ * (piG_ * exp1_); // G
376  p_(3, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_ * exp21_); // T, U
377 
378  return p_;
379 }
380 
381 const Matrix<double>& HKY85::getd2Pij_dt2(double d) const
382 {
383  double r_2 = rate_ * rate_ * r_ * r_;
384  l_ = rate_ * r_ * d;
385  double k1_2 = k1_ * k1_;
386  double k2_2 = k2_ * k2_;
387  exp1_ = exp(-l_);
388  exp22_ = exp(-k2_ * l_);
389  exp21_ = exp(-k1_ * l_);
390 
391  // A
392  p_(0, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_2 * exp22_); // A
393  p_(0, 1) = r_2 * (piC_ * -exp1_); // C
394  p_(0, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_2 * exp22_); // G
395  p_(0, 3) = r_2 * (piT_ * -exp1_); // T, U
396 
397  // C
398  p_(1, 0) = r_2 * (piA_ * -exp1_); // A
399  p_(1, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_2 * exp21_); // C
400  p_(1, 2) = r_2 * (piG_ * -exp1_); // G
401  p_(1, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_2 * exp21_); // T, U
402 
403  // G
404  p_(2, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_2 * exp22_); // A
405  p_(2, 1) = r_2 * (piC_ * -exp1_); // C
406  p_(2, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_2 * exp22_); // G
407  p_(2, 3) = r_2 * (piT_ * -exp1_); // T, U
408 
409  // T, U
410  p_(3, 0) = r_2 * (piA_ * -exp1_); // A
411  p_(3, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_2 * exp21_); // C
412  p_(3, 2) = r_2 * (piG_ * -exp1_); // G
413  p_(3, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_2 * exp21_); // T, U
414 
415  return p_;
416 }
417 
418 /******************************************************************************/
419 
420 void HKY85::setFreq(std::map<int, double>& freqs)
421 {
422  piA_ = freqs[0];
423  piC_ = freqs[1];
424  piG_ = freqs[2];
425  piT_ = freqs[3];
426  vector<string> thetas(3);
427  thetas[0] = getNamespace() + "theta";
428  thetas[1] = getNamespace() + "theta1";
429  thetas[2] = getNamespace() + "theta2";
430  ParameterList pl = getParameters().createSubList(thetas);
431  pl[0].setValue(piC_ + piG_);
432  pl[1].setValue(piA_ / (piA_ + piT_));
433  pl[2].setValue(piG_ / (piC_ + piG_));
434  setParametersValues(pl);
435 }
436 
437 /******************************************************************************/
bpp::AbstractParameterAliasable::addParameter_
void addParameter_(Parameter *parameter)
bpp::AbstractParameterAliasable::setParametersValues
void setParametersValues(const ParameterList &parameters) override
bpp::AbstractParameterAliasable::getNamespace
std::string getNamespace() const override
bpp::AbstractParameterAliasable::getParameters
const ParameterList & getParameters() const 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::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::HKY85::piY_
double piY_
Definition: HKY85.h:147
bpp::HKY85::piC_
double piC_
Definition: HKY85.h:147
bpp::HKY85::piG_
double piG_
Definition: HKY85.h:147
bpp::HKY85::theta1_
double theta1_
Definition: HKY85.h:147
bpp::HKY85::dPij_dt
double dPij_dt(size_t i, size_t j, double d) const
Definition: HKY85.cpp:208
bpp::HKY85::updateMatrices_
void updateMatrices_()
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: HKY85.cpp:42
bpp::HKY85::HKY85
HKY85(std::shared_ptr< const NucleicAlphabet > alpha, double kappa=1., double piA=0.25, double piC=0.25, double piG=0.25, double piT=0.25)
Definition: HKY85.cpp:19
bpp::HKY85::k1_
double k1_
Definition: HKY85.h:147
bpp::HKY85::Pij_t
double Pij_t(size_t i, size_t j, double d) const
Definition: HKY85.cpp:157
bpp::HKY85::getPij_t
const Matrix< double > & getPij_t(double d) const
Definition: HKY85.cpp:313
bpp::HKY85::k2_
double k2_
Definition: HKY85.h:147
bpp::HKY85::d2Pij_dt2
double d2Pij_dt2(size_t i, size_t j, double d) const
Definition: HKY85.cpp:259
bpp::HKY85::exp21_
double exp21_
Definition: HKY85.h:148
bpp::HKY85::kappa_
double kappa_
Definition: HKY85.h:147
bpp::HKY85::exp1_
double exp1_
Definition: HKY85.h:148
bpp::HKY85::l_
double l_
Definition: HKY85.h:148
bpp::HKY85::piA_
double piA_
Definition: HKY85.h:147
bpp::HKY85::theta2_
double theta2_
Definition: HKY85.h:147
bpp::HKY85::exp22_
double exp22_
Definition: HKY85.h:148
bpp::HKY85::p_
RowMatrix< double > p_
Definition: HKY85.h:149
bpp::HKY85::theta_
double theta_
Definition: HKY85.h:147
bpp::HKY85::setFreq
void setFreq(std::map< int, double > &freqs)
This method is redefined to actualize the corresponding parameters piA, piT, piG and piC too.
Definition: HKY85.cpp:420
bpp::HKY85::r_
double r_
Definition: HKY85.h:147
bpp::HKY85::piT_
double piT_
Definition: HKY85.h:147
bpp::HKY85::getdPij_dt
const Matrix< double > & getdPij_dt(double d) const
Definition: HKY85.cpp:347
bpp::HKY85::getd2Pij_dt2
const Matrix< double > & getd2Pij_dt2(double d) const
Definition: HKY85.cpp:381
bpp::HKY85::piR_
double piR_
Definition: HKY85.h:147
bpp::ParameterList::createSubList
virtual ParameterList createSubList(const std::vector< std::string > &names) const
bpp::Parameter::R_PLUS_STAR
static const std::shared_ptr< IntervalConstraint > R_PLUS_STAR
bpp
Defines the basic types of data flow nodes.
bpp::exp
ExtendedFloat exp(const ExtendedFloat &ef)
Definition: ExtendedFloat.h:474