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