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 
5 #include <Bpp/Numeric/Matrix/MatrixTools.h>
6 
7 #include "../FrequencySet/NucleotideFrequencySet.h"
8 #include "F84.h"
9 
10 // From the STL:
11 #include <cmath>
12 
13 using namespace bpp;
14 using namespace std;
15 
16 /******************************************************************************/
17 
18 F84::F84(
19  shared_ptr<const NucleicAlphabet> alpha,
20  double kappa,
21  double piA,
22  double piC,
23  double piG,
24  double piT) :
25  AbstractParameterAliasable("F84."),
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));
32  addParameter_(new Parameter("F84.theta", theta_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
33  addParameter_(new Parameter("F84.theta1", theta1_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
34  addParameter_(new Parameter("F84.theta2", theta2_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
35  updateMatrices_();
36 }
37 
38 /******************************************************************************/
39 
40 void F84::updateMatrices_()
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_;
48  piG_ = 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 
157 double 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 
207 double 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 
257 double 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 
310 const 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 
343 const 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 
376 const 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 
413 void 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";
423  ParameterList pl = getParameters().createSubList(thetas);
424  pl[0].setValue(piC_ + piG_);
425  pl[1].setValue(piA_ / (piA_ + piT_));
426  pl[2].setValue(piG_ / (piC_ + piG_));
427  setParametersValues(pl);
428 }
429 
430 /******************************************************************************/
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::F84::p_
RowMatrix< double > p_
Definition: F84.h:147
bpp::F84::kappa_
double kappa_
Definition: F84.h:145
bpp::F84::getdPij_dt
const Matrix< double > & getdPij_dt(double d) const
Definition: F84.cpp:343
bpp::F84::k1_
double k1_
Definition: F84.h:145
bpp::F84::getd2Pij_dt2
const Matrix< double > & getd2Pij_dt2(double d) const
Definition: F84.cpp:376
bpp::F84::exp2_
double exp2_
Definition: F84.h:146
bpp::F84::piC_
double piC_
Definition: F84.h:145
bpp::F84::exp1_
double exp1_
Definition: F84.h:146
bpp::F84::getPij_t
const Matrix< double > & getPij_t(double d) const
Definition: F84.cpp:310
bpp::F84::updateMatrices_
void updateMatrices_()
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: F84.cpp:40
bpp::F84::theta_
double theta_
Definition: F84.h:145
bpp::F84::piG_
double piG_
Definition: F84.h:145
bpp::F84::Pij_t
double Pij_t(size_t i, size_t j, double d) const
Definition: F84.cpp:157
bpp::F84::piY_
double piY_
Definition: F84.h:145
bpp::F84::setFreq
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
bpp::F84::dPij_dt
double dPij_dt(size_t i, size_t j, double d) const
Definition: F84.cpp:207
bpp::F84::theta2_
double theta2_
Definition: F84.h:145
bpp::F84::piT_
double piT_
Definition: F84.h:145
bpp::F84::k2_
double k2_
Definition: F84.h:145
bpp::F84::piA_
double piA_
Definition: F84.h:145
bpp::F84::r_
double r_
Definition: F84.h:145
bpp::F84::F84
F84(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: F84.cpp:18
bpp::F84::l_
double l_
Definition: F84.h:146
bpp::F84::theta1_
double theta1_
Definition: F84.h:145
bpp::F84::piR_
double piR_
Definition: F84.h:145
bpp::F84::d2Pij_dt2
double d2Pij_dt2(size_t i, size_t j, double d) const
Definition: F84.cpp:257
bpp::FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL
static std::shared_ptr< IntervalConstraint > FREQUENCE_CONSTRAINT_SMALL
Definition: FrequencySet.h:90
bpp::ParameterList::createSubList
virtual ParameterList createSubList(const std::vector< std::string > &names) const
bpp::Parameter::R_PLUS
static const std::shared_ptr< IntervalConstraint > R_PLUS
bpp
Defines the basic types of data flow nodes.
bpp::exp
ExtendedFloat exp(const ExtendedFloat &ef)
Definition: ExtendedFloat.h:474