bpp-phyl  2.2.0
K80.cpp
Go to the documentation of this file.
1 //
2 // File: K80.cpp
3 // Created by: Julien Dutheil
4 // Created on: Tue May 27 15:24:30 2003
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (November 16, 2004)
9 
10  This software is a computer program whose purpose is to provide classes
11  for phylogenetic data analysis.
12 
13  This software is governed by the CeCILL license under French law and
14  abiding by the rules of distribution of free software. You can use,
15  modify and/ or redistribute the software under the terms of the CeCILL
16  license as circulated by CEA, CNRS and INRIA at the following URL
17  "http://www.cecill.info".
18 
19  As a counterpart to the access to the source code and rights to copy,
20  modify and redistribute granted by the license, users are provided only
21  with a limited warranty and the software's author, the holder of the
22  economic rights, and the successive licensors have only limited
23  liability.
24 
25  In this respect, the user's attention is drawn to the risks associated
26  with loading, using, modifying and/or developing or reproducing the
27  software by the user in light of its specific status of free software,
28  that may mean that it is complicated to manipulate, and that also
29  therefore means that it is reserved for developers and experienced
30  professionals having in-depth computer knowledge. Users are therefore
31  encouraged to load and test the software's suitability as regards their
32  requirements in conditions enabling the security of their systems and/or
33  data to be ensured and, more generally, to use and operate it in the
34  same conditions as regards security.
35 
36  The fact that you are presently reading this means that you have had
37  knowledge of the CeCILL license and that you accept its terms.
38 */
39 
40 #include "K80.h"
41 
42 // From the STL:
43 #include <cmath>
44 
45 using namespace bpp;
46 
47 #include <Bpp/Numeric/Matrix/MatrixTools.h>
48 
49 using namespace std;
50 
51 /******************************************************************************/
52 
53 K80::K80(const NucleicAlphabet* alpha, double kappa) :
54 AbstractParameterAliasable("K80."),
55  AbstractReversibleSubstitutionModel(alpha, new CanonicalStateMap(alpha, false), "K80."),
56  kappa_(kappa), r_(), l_(), k_(), exp1_(), exp2_(), p_(size_, size_)
57 {
58  addParameter_(new Parameter("K80.kappa", kappa, &Parameter::R_PLUS_STAR));
59  updateMatrices();
60 }
61 
62 /******************************************************************************/
63 
64 void K80::updateMatrices()
65 {
66  kappa_ = getParameterValue("kappa");
67  k_ = (kappa_ + 1.) / 2.;
68  r_ = 4. / (kappa_ + 2.);
69 
70  // Frequences:
71  freq_[0] = freq_[1] = freq_[2] = freq_[3] = 1. / 4.;
72 
73  // Generator:
74  generator_(0, 0) = -2. - kappa_;
75  generator_(1, 1) = -2. - kappa_;
76  generator_(2, 2) = -2. - kappa_;
77  generator_(3, 3) = -2. - kappa_;
78 
79  generator_(0, 1) = 1.;
80  generator_(0, 3) = 1.;
81  generator_(1, 0) = 1.;
82  generator_(1, 2) = 1.;
83  generator_(2, 1) = 1.;
84  generator_(2, 3) = 1.;
85  generator_(3, 0) = 1.;
86  generator_(3, 2) = 1.;
87 
88  generator_(0, 2) = kappa_;
89  generator_(1, 3) = kappa_;
90  generator_(2, 0) = kappa_;
91  generator_(3, 1) = kappa_;
92 
93  // Normalization:
94  MatrixTools::scale(generator_, r_/4);
95 
96  // Exchangeability:
97  exchangeability_ = generator_;
98  MatrixTools::scale(exchangeability_, 4.);
99 
100  // Eigen values:
101  eigenValues_[0] = 0;
102  eigenValues_[1] = -r_ * (1. + kappa_)/2;
103  eigenValues_[2] = -r_ * (1. + kappa_)/2;
104  eigenValues_[3] = -r_;
105 
106  // Eigen vectors:
107  leftEigenVectors_(0,0) = 1. / 4.;
108  leftEigenVectors_(0,1) = 1. / 4.;
109  leftEigenVectors_(0,2) = 1. / 4.;
110  leftEigenVectors_(0,3) = 1. / 4.;
111  leftEigenVectors_(1,0) = 0.;
112  leftEigenVectors_(1,1) = 1. / 2.;
113  leftEigenVectors_(1,2) = 0.;
114  leftEigenVectors_(1,3) = -1. / 2.;
115  leftEigenVectors_(2,0) = 1. / 2.;
116  leftEigenVectors_(2,1) = 0.;
117  leftEigenVectors_(2,2) = -1. / 2.;
118  leftEigenVectors_(2,3) = 0.;
119  leftEigenVectors_(3,0) = 1. / 4.;
120  leftEigenVectors_(3,1) = -1. / 4.;
121  leftEigenVectors_(3,2) = 1. / 4.;
122  leftEigenVectors_(3,3) = -1. / 4.;
123 
124  rightEigenVectors_(0,0) = 1.;
125  rightEigenVectors_(0,1) = 0.;
126  rightEigenVectors_(0,2) = 1.;
127  rightEigenVectors_(0,3) = 1.;
128  rightEigenVectors_(1,0) = 1.;
129  rightEigenVectors_(1,1) = 1.;
130  rightEigenVectors_(1,2) = 0.;
131  rightEigenVectors_(1,3) = -1.;
132  rightEigenVectors_(2,0) = 1.;
133  rightEigenVectors_(2,1) = 0.;
134  rightEigenVectors_(2,2) = -1.;
135  rightEigenVectors_(2,3) = 1.;
136  rightEigenVectors_(3,0) = 1.;
137  rightEigenVectors_(3,1) = -1.;
138  rightEigenVectors_(3,2) = 0;
139  rightEigenVectors_(3,3) = -1.;
140 }
141 
142 /******************************************************************************/
143 
144 double K80::Pij_t(size_t i, size_t j, double d) const
145 {
146  l_ = rate_ * r_ * d;
147  exp1_ = exp(-l_);
148  exp2_ = exp(-k_ * l_);
149 
150  switch(i) {
151  //A
152  case 0 : {
153  switch(j) {
154  case 0 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //A
155  case 1 : return 0.25 * (1. - exp1_); //C
156  case 2 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //G
157  case 3 : return 0.25 * (1. - exp1_); //T, U
158  }
159  }
160  //C
161  case 1 : {
162  switch(j) {
163  case 0 : return 0.25 * (1. - exp1_); //A
164  case 1 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //C
165  case 2 : return 0.25 * (1. - exp1_); //G
166  case 3 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //T, U
167  }
168  }
169  //G
170  case 2 : {
171  switch(j) {
172  case 0 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //A
173  case 1 : return 0.25 * (1. - exp1_); //C
174  case 2 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //G
175  case 3 : return 0.25 * (1. - exp1_); //T, U
176  }
177  }
178  //T, U
179  case 3 : {
180  switch(j) {
181  case 0 : return 0.25 * (1. - exp1_); //A
182  case 1 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //C
183  case 2 : return 0.25 * (1. - exp1_); //G
184  case 3 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //T, U
185  }
186  }
187  }
188  return 0;
189 }
190 
191 /******************************************************************************/
192 
193 double K80::dPij_dt(size_t i, size_t j, double d) const
194 {
195  l_ = rate_ * r_ * d;
196  exp1_ = exp(-l_);
197  exp2_ = exp(-k_ * l_);
198 
199  switch(i) {
200  //A
201  case 0 : {
202  switch(j) {
203  case 0 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //A
204  case 1 : return rate_ * r_/4. * ( exp1_); //C
205  case 2 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //G
206  case 3 : return rate_ * r_/4. * ( exp1_); //T, U
207  }
208  }
209  //C
210  case 1 : {
211  switch(j) {
212  case 0 : return rate_ * r_/4. * ( exp1_); //A
213  case 1 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //C
214  case 2 : return rate_ * r_/4. * ( exp1_); //G
215  case 3 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //T, U
216  }
217  }
218  //G
219  case 2 : {
220  switch(j) {
221  case 0 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //A
222  case 1 : return rate_ * r_/4. * ( exp1_); //C
223  case 2 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //G
224  case 3 : return rate_ * r_/4. * ( exp1_); //T, U
225  }
226  }
227  //T, U
228  case 3 : {
229  switch(j) {
230  case 0 : return rate_ * r_/4. * ( exp1_); //A
231  case 1 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //C
232  case 2 : return rate_ * r_/4. * ( exp1_); //G
233  case 3 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //T, U
234  }
235  }
236  }
237  return 0;
238 }
239 
240 /******************************************************************************/
241 
242 double K80::d2Pij_dt2(size_t i, size_t j, double d) const
243 {
244  double k_2 = k_ * k_;
245  double r_2 = rate_ * rate_ * r_ * r_;
246  l_ = rate_ * r_ * d;
247  exp1_ = exp(-l_);
248  exp2_ = exp(-k_ * l_);
249 
250  switch(i) {
251  //A
252  case 0 : {
253  switch(j) {
254  case 0 : return r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //A
255  case 1 : return r_2/4. * (- exp1_); //C
256  case 2 : return r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //G
257  case 3 : return r_2/4. * (- exp1_); //T, U
258  }
259  }
260  //C
261  case 1 : {
262  switch(j) {
263  case 0 : return r_2/4. * (- exp1_); //A
264  case 1 : return r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //C
265  case 2 : return r_2/4. * (- exp1_); //G
266  case 3 : return r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //T, U
267  }
268  }
269  //G
270  case 2 : {
271  switch(j) {
272  case 0 : return r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //A
273  case 1 : return r_2/4. * (- exp1_); //C
274  case 2 : return r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //G
275  case 3 : return r_2/4. * (- exp1_); //T, U
276  }
277  }
278  //T, U
279  case 3 : {
280  switch(j) {
281  case 0 : return r_2/4. * (- exp1_); //A
282  case 1 : return r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //C
283  case 2 : return r_2/4. * (- exp1_); //G
284  case 3 : return r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //T, U
285  }
286  }
287  }
288  return 0;
289 }
290 
291 /******************************************************************************/
292 
293 const Matrix<double> & K80::getPij_t(double d) const
294 {
295  l_ = rate_ * r_ * d;
296  exp1_ = exp(-l_);
297  exp2_ = exp(-k_ * l_);
298 
299  //A
300  p_(0, 0) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //A
301  p_(0, 1) = 0.25 * (1. - exp1_); //C
302  p_(0, 2) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //G
303  p_(0, 3) = 0.25 * (1. - exp1_); //T, U
304 
305  //C
306  p_(1, 0) = 0.25 * (1. - exp1_); //A
307  p_(1, 1) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //C
308  p_(1, 2) = 0.25 * (1. - exp1_); //G
309  p_(1, 3) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //T, U
310 
311  //G
312  p_(2, 0) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //A
313  p_(2, 1) = 0.25 * (1. - exp1_); //C
314  p_(2, 2) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //G
315  p_(2, 3) = 0.25 * (1. - exp1_); //T, U
316 
317  //T, U
318  p_(3, 0) = 0.25 * (1. - exp1_); //A
319  p_(3, 1) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //C
320  p_(3, 2) = 0.25 * (1. - exp1_); //G
321  p_(3, 3) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //T, U
322 
323  return p_;
324 }
325 
326 const Matrix<double> & K80::getdPij_dt(double d) const
327 {
328  l_ = rate_ * r_ * d;
329  exp1_ = exp(-l_);
330  exp2_ = exp(-k_ * l_);
331 
332  p_(0, 0) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //A
333  p_(0, 1) = rate_ * r_/4. * ( exp1_); //C
334  p_(0, 2) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //G
335  p_(0, 3) = rate_ * r_/4. * ( exp1_); //T, U
336 
337  //C
338  p_(1, 0) = rate_ * r_/4. * ( exp1_); //A
339  p_(1, 1) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //C
340  p_(1, 2) = rate_ * r_/4. * ( exp1_); //G
341  p_(1, 3) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //T, U
342 
343  //G
344  p_(2, 0) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //A
345  p_(2, 1) = rate_ * r_/4. * ( exp1_); //C
346  p_(2, 2) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //G
347  p_(2, 3) = rate_ * r_/4. * ( exp1_); //T, U
348 
349  //T, U
350  p_(3, 0) = rate_ * r_/4. * ( exp1_); //A
351  p_(3, 1) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //C
352  p_(3, 2) = rate_ * r_/4. * ( exp1_); //G
353  p_(3, 3) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //T, U
354 
355  return p_;
356 }
357 
358 const Matrix<double> & K80::getd2Pij_dt2(double d) const
359 {
360  double k_2 = k_ * k_;
361  double r_2 = rate_ * rate_ * r_ * r_;
362  l_ = rate_ * r_ * d;
363  exp1_ = exp(-l_);
364  exp2_ = exp(-k_ * l_);
365 
366  p_(0, 0) = r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //A
367  p_(0, 1) = r_2/4. * (- exp1_); //C
368  p_(0, 2) = r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //G
369  p_(0, 3) = r_2/4. * (- exp1_); //T, U
370 
371  //C
372  p_(1, 0) = r_2/4. * (- exp1_); //A
373  p_(1, 1) = r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //C
374  p_(1, 2) = r_2/4. * (- exp1_); //G
375  p_(1, 3) = r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //T, U
376 
377  //G
378  p_(2, 0) = r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //A
379  p_(2, 1) = r_2/4. * (- exp1_); //C
380  p_(2, 2) = r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //G
381  p_(2, 3) = r_2/4. * (- exp1_); //T, U
382 
383  //T, U
384  p_(3, 0) = r_2/4. * (- exp1_); //A
385  p_(3, 1) = r_2/4. * ( exp1_ - 2. * k_2 * exp2_); //C
386  p_(3, 2) = r_2/4. * (- exp1_); //G
387  p_(3, 3) = r_2/4. * ( exp1_ + 2. * k_2 * exp2_); //T, U
388 
389  return p_;
390 }
391 
392 /******************************************************************************/
393 
bpp::K80::getdPij_dt
const Matrix< double > & getdPij_dt(double d) const
Definition: K80.cpp:326
bpp::K80::d2Pij_dt2
double d2Pij_dt2(size_t i, size_t j, double d) const
Definition: K80.cpp:242
bpp::K80::dPij_dt
double dPij_dt(size_t i, size_t j, double d) const
Definition: K80.cpp:193
bpp::K80::kappa_
double kappa_
Definition: K80.h:155
bpp::K80::k_
double k_
Definition: K80.h:156
bpp::AbstractSubstitutionModel::exchangeability_
RowMatrix< double > exchangeability_
The exchangeability matrix of the model, defined as . When the model is reversible, this matrix is symetric.
Definition: AbstractSubstitutionModel.h:122
bpp::CanonicalStateMap
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:161
bpp::K80::exp2_
double exp2_
Definition: K80.h:156
bpp::AbstractSubstitutionModel::rightEigenVectors_
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
Definition: AbstractSubstitutionModel.h:155
bpp::AbstractSubstitutionModel::eigenValues_
Vdouble eigenValues_
The vector of eigen values.
Definition: AbstractSubstitutionModel.h:139
bpp::K80::exp1_
double exp1_
Definition: K80.h:156
std
STL namespace.
bpp::AbstractSubstitutionModel::freq_
Vdouble freq_
The vector of equilibrium frequencies.
Definition: AbstractSubstitutionModel.h:113
bpp::K80::r_
double r_
Definition: K80.h:155
bpp::AbstractReversibleSubstitutionModel
Partial implementation of the ReversibleSubstitutionModel interface.
Definition: AbstractSubstitutionModel.h:372
bpp::K80::p_
RowMatrix< double > p_
Definition: K80.h:157
bpp::K80::K80
K80(const NucleicAlphabet *alpha, double kappa=1.)
Definition: K80.cpp:53
bpp::AbstractSubstitutionModel::leftEigenVectors_
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
Definition: AbstractSubstitutionModel.h:166
bpp::AbstractSubstitutionModel::generator_
RowMatrix< double > generator_
The generator matrix of the model.
Definition: AbstractSubstitutionModel.h:108
bpp::K80::Pij_t
double Pij_t(size_t i, size_t j, double d) const
Definition: K80.cpp:144
bpp::AbstractSubstitutionModel::rate_
double rate_
The rate of the model (default: 1). The generator (and all its vectorial components) is independent o...
Definition: AbstractSubstitutionModel.h:103
bpp::K80::l_
double l_
Definition: K80.h:156
bpp::K80::updateMatrices
void updateMatrices()
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: K80.cpp:64
bpp::K80::getPij_t
const Matrix< double > & getPij_t(double d) const
Definition: K80.cpp:293
bpp::K80::getd2Pij_dt2
const Matrix< double > & getd2Pij_dt2(double d) const
Definition: K80.cpp:358