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