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
14using namespace std;
15
16using namespace bpp;
17
18/******************************************************************************/
19
20TN93::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
47void 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
129double 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
180double 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
231double 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
285const 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
319const 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
353const 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
392void 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