bpp-phyl3  3.0.0
TN93.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/EigenValue.h>
6 #include <Bpp/Numeric/Matrix/MatrixTools.h>
7 
8 #include "../FrequencySet/NucleotideFrequencySet.h"
9 #include "TN93.h"
10 
11 // From the STL:
12 #include <cmath>
13 
14 using namespace std;
15 
16 using namespace bpp;
17 
18 /******************************************************************************/
19 
20 TN93::TN93(
21  shared_ptr<const NucleicAlphabet> alpha,
22  double kappa1,
23  double kappa2,
24  double piA,
25  double piC,
26  double piG,
27  double piT) :
28  AbstractParameterAliasable("TN93."),
29  AbstractReversibleNucleotideSubstitutionModel(alpha, make_shared<CanonicalStateMap>(alpha, false), "TN93."),
30  kappa1_(kappa1), kappa2_(kappa2),
31  piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
32  r_(), k1_(), k2_(),
33  theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
34  exp1_(), exp21_(), exp22_(), l_(), p_(size_, size_)
35 {
36  addParameter_(new Parameter("TN93.kappa1", kappa1, Parameter::R_PLUS_STAR));
37  addParameter_(new Parameter("TN93.kappa2", kappa2, Parameter::R_PLUS_STAR));
38  addParameter_(new Parameter("TN93.theta", theta_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
39  addParameter_(new Parameter("TN93.theta1", theta1_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
40  addParameter_(new Parameter("TN93.theta2", theta2_, FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL));
41  p_.resize(size_, size_);
42  updateMatrices_();
43 }
44 
45 /******************************************************************************/
46 
47 void TN93::updateMatrices_()
48 {
49  kappa1_ = getParameterValue("kappa1");
50  kappa2_ = getParameterValue("kappa2");
51  theta_ = getParameterValue("theta" );
52  theta1_ = getParameterValue("theta1");
53  theta2_ = getParameterValue("theta2");
54  piA_ = theta1_ * (1. - theta_);
55  piC_ = (1. - theta2_) * theta_;
56  piG_ = theta2_ * theta_;
57  piT_ = (1. - theta1_) * (1. - theta_);
58  piR_ = piA_ + piG_;
59  piY_ = piT_ + piC_;
60  r_ = isScalable() ? 1. / (2. * (piA_ * piC_ + piC_ * piG_ + piA_ * piT_ + piG_ * piT_ + kappa2_ * piC_ * piT_ + kappa1_ * piA_ * piG_)) : 1;
61 
62  k1_ = kappa2_ * piY_ + piR_;
63  k2_ = kappa1_ * piR_ + piY_;
64 
65  freq_[0] = piA_;
66  freq_[1] = piC_;
67  freq_[2] = piG_;
68  freq_[3] = piT_;
69 
70  generator_(0, 0) = -( piC_ + kappa1_ * piG_ + piT_);
71  generator_(1, 1) = -( piA_ + piG_ + kappa2_ * piT_);
72  generator_(2, 2) = -(kappa1_ * piA_ + piC_ + piT_);
73  generator_(3, 3) = -( piA_ + kappa2_ * piC_ + piG_ );
74 
75  generator_(1, 0) = piA_;
76  generator_(3, 0) = piA_;
77  generator_(0, 1) = piC_;
78  generator_(2, 1) = piC_;
79  generator_(1, 2) = piG_;
80  generator_(3, 2) = piG_;
81  generator_(0, 3) = piT_;
82  generator_(2, 3) = piT_;
83 
84  generator_(2, 0) = kappa1_ * piA_;
85  generator_(3, 1) = kappa2_ * piC_;
86  generator_(0, 2) = kappa1_ * piG_;
87  generator_(1, 3) = kappa2_ * piT_;
88 
89  // Normalization:
90  setScale(r_);
91 
92  // Exchangeability:
93  exchangeability_(0, 0) = generator_(0, 0) / piA_;
94  exchangeability_(0, 1) = generator_(0, 1) / piC_;
95  exchangeability_(0, 2) = generator_(0, 2) / piG_;
96  exchangeability_(0, 3) = generator_(0, 3) / piT_;
97 
98  exchangeability_(1, 0) = generator_(1, 0) / piA_;
99  exchangeability_(1, 1) = generator_(1, 1) / piC_;
100  exchangeability_(1, 2) = generator_(1, 2) / piG_;
101  exchangeability_(1, 3) = generator_(1, 3) / piT_;
102 
103  exchangeability_(2, 0) = generator_(2, 0) / piA_;
104  exchangeability_(2, 1) = generator_(2, 1) / piC_;
105  exchangeability_(2, 2) = generator_(2, 2) / piG_;
106  exchangeability_(2, 3) = generator_(2, 3) / piT_;
107 
108  exchangeability_(3, 0) = generator_(3, 0) / piA_;
109  exchangeability_(3, 1) = generator_(3, 1) / piC_;
110  exchangeability_(3, 2) = generator_(3, 2) / piG_;
111  exchangeability_(3, 3) = generator_(3, 3) / piT_;
112 
113  // We are not sure that the values are computed in this order :p
114  // Eigen values:
115  // eigenValues_[0] = 0;
116  // eigenValues_[1] = -r * (kappa2 * piY + piR);
117  // eigenValues_[2] = -r * (kappa1 * piR + piY);
118  // eigenValues_[3] = -r;
119 
120  // Eigen vectors and values:
121  EigenValue<double> ev(generator_);
122  rightEigenVectors_ = ev.getV();
123  MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
124  eigenValues_ = ev.getRealEigenValues();
125 }
126 
127 /******************************************************************************/
128 
129 double TN93::Pij_t(size_t i, size_t j, double d) const
130 {
131  l_ = rate_ * r_ * d;
132  exp1_ = exp(-l_);
133  exp22_ = exp(-k2_ * l_);
134  exp21_ = exp(-k1_ * l_);
135 
136  switch (i)
137  {
138  case 0: // A
139  switch (j)
140  {
141  case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp22_; // A
142  case 1: return piC_ * (1. - exp1_); // C
143  case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp22_; // G
144  case 3: return piT_ * (1. - exp1_); // T, U
145  default: return 0;
146  }
147  case 1: // C
148  switch (j)
149  {
150  case 0: return piA_ * (1. - exp1_); // A
151  case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp21_; // C
152  case 2: return piG_ * (1. - exp1_); // G
153  case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp21_; // T, U
154  default: return 0;
155  }
156  case 2: // G
157  switch (j)
158  {
159  case 0: return piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp22_; // A
160  case 1: return piC_ * (1. - exp1_); // C
161  case 2: return piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp22_; // G
162  case 3: return piT_ * (1. - exp1_); // T, U
163  default: return 0;
164  }
165  case 3: // T, U
166  switch (j)
167  {
168  case 0: return piA_ * (1. - exp1_); // A
169  case 1: return piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp21_; // C
170  case 2: return piG_ * (1. - exp1_); // G
171  case 3: return piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp21_; // T, U
172  default: return 0;
173  }
174  default: return 0;
175  }
176 }
177 
178 /******************************************************************************/
179 
180 double TN93::dPij_dt(size_t i, size_t j, double d) const
181 {
182  l_ = rate_ * r_ * d;
183  exp1_ = exp(-l_);
184  exp22_ = exp(-k2_ * l_);
185  exp21_ = exp(-k1_ * l_);
186 
187  switch (i)
188  {
189  case 0: // A
190  switch (j)
191  {
192  case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_ * exp22_); // A
193  case 1: return rate_ * r_ * (piC_ * exp1_); // C
194  case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_ * exp22_); // G
195  case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
196  default: return 0;
197  }
198  case 1: // C
199  switch (j)
200  {
201  case 0: return rate_ * r_ * (piA_ * exp1_); // A
202  case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_ * exp21_); // C
203  case 2: return rate_ * r_ * (piG_ * exp1_); // G
204  case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_ * exp21_); // T, U
205  default: return 0;
206  }
207  case 2: // G
208  switch (j)
209  {
210  case 0: return rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_ * exp22_); // A
211  case 1: return rate_ * r_ * (piC_ * exp1_); // C
212  case 2: return rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_ * exp22_); // G
213  case 3: return rate_ * r_ * (piT_ * exp1_); // T, U
214  default: return 0;
215  }
216  case 3: // T, U
217  switch (j)
218  {
219  case 0: return rate_ * r_ * (piA_ * exp1_); // A
220  case 1: return rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_ * exp21_); // C
221  case 2: return rate_ * r_ * (piG_ * exp1_); // G
222  case 3: return rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_ * exp21_); // T, U
223  default: return 0;
224  }
225  default: return 0;
226  }
227 }
228 
229 /******************************************************************************/
230 
231 double TN93::d2Pij_dt2(size_t i, size_t j, double d) const
232 {
233  double r_2 = rate_ * rate_ * r_ * r_;
234  l_ = rate_ * r_ * d;
235  double k1_2 = k1_ * k1_;
236  double k2_2 = k2_ * k2_;
237  exp1_ = exp(-l_);
238  exp22_ = exp(-k2_ * l_);
239  exp21_ = exp(-k1_ * l_);
240 
241  switch (i)
242  {
243  case 0: // A
244  switch (j)
245  {
246  case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_2 * exp22_); // A
247  case 1: return r_2 * (piC_ * -exp1_); // C
248  case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_2 * exp22_); // G
249  case 3: return r_2 * (piT_ * -exp1_); // T, U
250  default: return 0;
251  }
252  case 1: // C
253  switch (j)
254  {
255  case 0: return r_2 * (piA_ * -exp1_); // A
256  case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_2 * exp21_); // C
257  case 2: return r_2 * (piG_ * -exp1_); // G
258  case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_2 * exp21_); // T, U
259  default: return 0;
260  }
261  case 2: // G
262  switch (j)
263  {
264  case 0: return r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_2 * exp22_); // A
265  case 1: return r_2 * (piC_ * -exp1_); // C
266  case 2: return r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_2 * exp22_); // G
267  case 3: return r_2 * (piT_ * -exp1_); // T, U
268  default: return 0;
269  }
270  case 3: // T, U
271  switch (j)
272  {
273  case 0: return r_2 * (piA_ * -exp1_); // A
274  case 1: return r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_2 * exp21_); // C
275  case 2: return r_2 * (piG_ * -exp1_); // G
276  case 3: return r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_2 * exp21_); // T, U
277  default: return 0;
278  }
279  default: return 0;
280  }
281 }
282 
283 /******************************************************************************/
284 
285 const Matrix<double>& TN93::getPij_t(double d) const
286 {
287  l_ = rate_ * r_ * d;
288  exp1_ = exp(-l_);
289  exp22_ = exp(-k2_ * l_);
290  exp21_ = exp(-k1_ * l_);
291 
292  // A
293  p_(0, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) + (piG_ / piR_) * exp22_; // A
294  p_(0, 1) = piC_ * (1. - exp1_); // C
295  p_(0, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) - (piG_ / piR_) * exp22_; // G
296  p_(0, 3) = piT_ * (1. - exp1_); // T, U
297 
298  // C
299  p_(1, 0) = piA_ * (1. - exp1_); // A
300  p_(1, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) + (piT_ / piY_) * exp21_; // C
301  p_(1, 2) = piG_ * (1. - exp1_); // G
302  p_(1, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) - (piT_ / piY_) * exp21_; // T, U
303 
304  // G
305  p_(2, 0) = piA_ * (1. + (piY_ / piR_) * exp1_) - (piA_ / piR_) * exp22_; // A
306  p_(2, 1) = piC_ * (1. - exp1_); // C
307  p_(2, 2) = piG_ * (1. + (piY_ / piR_) * exp1_) + (piA_ / piR_) * exp22_; // G
308  p_(2, 3) = piT_ * (1. - exp1_); // T, U
309 
310  // T, U
311  p_(3, 0) = piA_ * (1. - exp1_); // A
312  p_(3, 1) = piC_ * (1. + (piR_ / piY_) * exp1_) - (piC_ / piY_) * exp21_; // C
313  p_(3, 2) = piG_ * (1. - exp1_); // G
314  p_(3, 3) = piT_ * (1. + (piR_ / piY_) * exp1_) + (piC_ / piY_) * exp21_; // T, U
315 
316  return p_;
317 }
318 
319 const Matrix<double>& TN93::getdPij_dt(double d) const
320 {
321  l_ = rate_ * r_ * d;
322  exp1_ = exp(-l_);
323  exp22_ = exp(-k2_ * l_);
324  exp21_ = exp(-k1_ * l_);
325 
326  // A
327  p_(0, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_ * exp22_); // A
328  p_(0, 1) = rate_ * r_ * (piC_ * exp1_); // C
329  p_(0, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_ * exp22_); // G
330  p_(0, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
331 
332  // C
333  p_(1, 0) = rate_ * r_ * (piA_ * exp1_); // A
334  p_(1, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_ * exp21_); // C
335  p_(1, 2) = rate_ * r_ * (piG_ * exp1_); // G
336  p_(1, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_ * exp21_); // T, U
337 
338  // G
339  p_(2, 0) = rate_ * r_ * (piA_ * -(piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_ * exp22_); // A
340  p_(2, 1) = rate_ * r_ * (piC_ * exp1_); // C
341  p_(2, 2) = rate_ * r_ * (piG_ * -(piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_ * exp22_); // G
342  p_(2, 3) = rate_ * r_ * (piT_ * exp1_); // T, U
343 
344  // T, U
345  p_(3, 0) = rate_ * r_ * (piA_ * exp1_); // A
346  p_(3, 1) = rate_ * r_ * (piC_ * -(piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_ * exp21_); // C
347  p_(3, 2) = rate_ * r_ * (piG_ * exp1_); // G
348  p_(3, 3) = rate_ * r_ * (piT_ * -(piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_ * exp21_); // T, U
349 
350  return p_;
351 }
352 
353 const Matrix<double>& TN93::getd2Pij_dt2(double d) const
354 {
355  double r_2 = rate_ * rate_ * r_ * r_;
356  l_ = rate_ * r_ * d;
357  double k1_2 = k1_ * k1_;
358  double k2_2 = k2_ * k2_;
359  exp1_ = exp(-l_);
360  exp22_ = exp(-k2_ * l_);
361  exp21_ = exp(-k1_ * l_);
362 
363  // A
364  p_(0, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ + (piG_ / piR_) * k2_2 * exp22_); // A
365  p_(0, 1) = r_2 * (piC_ * -exp1_); // C
366  p_(0, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ - (piG_ / piR_) * k2_2 * exp22_); // G
367  p_(0, 3) = r_2 * (piT_ * -exp1_); // T, U
368 
369  // C
370  p_(1, 0) = r_2 * (piA_ * -exp1_); // A
371  p_(1, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ + (piT_ / piY_) * k1_2 * exp21_); // C
372  p_(1, 2) = r_2 * (piG_ * -exp1_); // G
373  p_(1, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ - (piT_ / piY_) * k1_2 * exp21_); // T, U
374 
375  // G
376  p_(2, 0) = r_2 * (piA_ * (piY_ / piR_) * exp1_ - (piA_ / piR_) * k2_2 * exp22_); // A
377  p_(2, 1) = r_2 * (piC_ * -exp1_); // C
378  p_(2, 2) = r_2 * (piG_ * (piY_ / piR_) * exp1_ + (piA_ / piR_) * k2_2 * exp22_); // G
379  p_(2, 3) = r_2 * (piT_ * -exp1_); // T, U
380 
381  // T, U
382  p_(3, 0) = r_2 * (piA_ * -exp1_); // A
383  p_(3, 1) = r_2 * (piC_ * (piR_ / piY_) * exp1_ - (piC_ / piY_) * k1_2 * exp21_); // C
384  p_(3, 2) = r_2 * (piG_ * -exp1_); // G
385  p_(3, 3) = r_2 * (piT_ * (piR_ / piY_) * exp1_ + (piC_ / piY_) * k1_2 * exp21_); // T, U
386 
387  return p_;
388 }
389 
390 /******************************************************************************/
391 
392 void TN93::setFreq(std::map<int, double>& freqs)
393 {
394  piA_ = freqs[0];
395  piC_ = freqs[1];
396  piG_ = freqs[2];
397  piT_ = freqs[3];
398  vector<string> thetas(3);
399  thetas[0] = getNamespace() + "theta";
400  thetas[1] = getNamespace() + "theta1";
401  thetas[2] = getNamespace() + "theta2";
402  ParameterList pl = getParameters().createSubList(thetas);
403  pl[0].setValue(piC_ + piG_);
404  pl[1].setValue(piA_ / (piA_ + piT_));
405  pl[2].setValue(piG_ / (piC_ + piG_));
406  setParametersValues(pl);
407 }
408 
409 /******************************************************************************/
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::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::EigenValue::getV
const RowMatrix< Real > & getV() const
bpp::EigenValue::getRealEigenValues
const std::vector< Real > & getRealEigenValues() const
bpp::FrequencySetInterface::FREQUENCE_CONSTRAINT_SMALL
static std::shared_ptr< IntervalConstraint > FREQUENCE_CONSTRAINT_SMALL
Definition: FrequencySet.h:90
bpp::MatrixTools::inv
static Scalar inv(const Matrix< Scalar > &A, Matrix< Scalar > &O)
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
RowMatrix< double >::resize
void resize(size_t nRows, size_t nCols)
bpp::TN93::Pij_t
double Pij_t(size_t i, size_t j, double d) const override
Definition: TN93.cpp:129
bpp::TN93::r_
double r_
Definition: TN93.h:97
bpp::TN93::k2_
double k2_
Definition: TN93.h:97
bpp::TN93::getdPij_dt
const Matrix< double > & getdPij_dt(double d) const override
Definition: TN93.cpp:319
bpp::TN93::l_
double l_
Definition: TN93.h:98
bpp::TN93::k1_
double k1_
Definition: TN93.h:97
bpp::TN93::kappa2_
double kappa2_
Definition: TN93.h:97
bpp::TN93::setFreq
void setFreq(std::map< int, double > &freqs) override
This method is over-defined to actualize the corresponding parameters piA, piT, piG and piC too.
Definition: TN93.cpp:392
bpp::TN93::piC_
double piC_
Definition: TN93.h:97
bpp::TN93::getd2Pij_dt2
const Matrix< double > & getd2Pij_dt2(double d) const override
Definition: TN93.cpp:353
bpp::TN93::exp22_
double exp22_
Definition: TN93.h:98
bpp::TN93::p_
RowMatrix< double > p_
Definition: TN93.h:99
bpp::TN93::piR_
double piR_
Definition: TN93.h:97
bpp::TN93::piY_
double piY_
Definition: TN93.h:97
bpp::TN93::piG_
double piG_
Definition: TN93.h:97
bpp::TN93::theta_
double theta_
Definition: TN93.h:97
bpp::TN93::updateMatrices_
void updateMatrices_() override
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: TN93.cpp:47
bpp::TN93::piT_
double piT_
Definition: TN93.h:97
bpp::TN93::exp21_
double exp21_
Definition: TN93.h:98
bpp::TN93::kappa1_
double kappa1_
Definition: TN93.h:97
bpp::TN93::piA_
double piA_
Definition: TN93.h:97
bpp::TN93::theta1_
double theta1_
Definition: TN93.h:97
bpp::TN93::dPij_dt
double dPij_dt(size_t i, size_t j, double d) const override
Definition: TN93.cpp:180
bpp::TN93::d2Pij_dt2
double d2Pij_dt2(size_t i, size_t j, double d) const override
Definition: TN93.cpp:231
bpp::TN93::theta2_
double theta2_
Definition: TN93.h:97
bpp::TN93::getPij_t
const Matrix< double > & getPij_t(double d) const override
Definition: TN93.cpp:285
bpp::TN93::exp1_
double exp1_
Definition: TN93.h:98
bpp
Defines the basic types of data flow nodes.
bpp::exp
ExtendedFloat exp(const ExtendedFloat &ef)
Definition: ExtendedFloat.h:474