bpp-phyl  2.2.0
F84.cpp
Go to the documentation of this file.
1 //
2 // File: F84.cpp
3 // Created by: Julien Dutheil
4 // Created on: Tue May 23 11:13 2005
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 "F84.h"
41 #include "../FrequenciesSet/NucleotideFrequenciesSet.h"
42 
43 #include <Bpp/Numeric/Matrix/MatrixTools.h>
44 
45 // From SeqLib:
46 #include <Bpp/Seq/Container/SequenceContainerTools.h>
47 
48 // From the STL:
49 #include <cmath>
50 
51 using namespace bpp;
52 using namespace std;
53 
54 /******************************************************************************/
55 
57  const NucleicAlphabet* alpha,
58  double kappa,
59  double piA,
60  double piC,
61  double piG,
62  double piT) :
63  AbstractParameterAliasable("F84."),
64  AbstractReversibleSubstitutionModel(alpha, new CanonicalStateMap(alpha, false), "F84."),
65  kappa_(kappa), piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
66  r_(), k1_(), k2_(), theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
67  l_(), exp1_(), exp2_(), p_(size_, size_)
68 {
69  addParameter_(new Parameter("F84.kappa", kappa, &Parameter::R_PLUS));
70  addParameter_(new Parameter("F84.theta" , theta_ , &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
71  addParameter_(new Parameter("F84.theta1", theta1_, &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
72  addParameter_(new Parameter("F84.theta2", theta2_, &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
74 }
75 
76 /******************************************************************************/
77 
79 {
80  kappa_ = getParameterValue("kappa");
81  theta_ = getParameterValue("theta");
82  theta1_ = getParameterValue("theta1");
83  theta2_ = getParameterValue("theta2");
84  piA_ = theta1_ * (1. - theta_);
85  piC_ = (1. - theta2_) * theta_;
86  piG_ = theta2_ * theta_;
87  piT_ = (1. - theta1_) * (1. - theta_);
88  piR_ = piA_ + piG_;
89  piY_ = piT_ + piC_;
90  r_ = 1. / (1 - piA_ * piA_ - piC_ * piC_ - piG_ * piG_ - piT_*piT_ + 2. * kappa_ * (piC_ * piT_ / piY_ + piA_ * piG_ / piR_));
91  k1_ = 1;
92  k2_ = kappa_ + 1;
93 
94  // Frequences:
95  freq_[0] = piA_;
96  freq_[1] = piC_;
97  freq_[2] = piG_;
98  freq_[3] = piT_;
99 
100  //Generator:
101  generator_(0, 0) = -( piC_ + (1 + kappa_ / piR_) * piG_ + piT_);
102  generator_(1, 1) = -( piA_ + piG_ + (1 + kappa_ / piY_) * piT_);
103  generator_(2, 2) = -((1 + kappa_ / piR_) * piA_ + piC_ + piT_);
104  generator_(3, 3) = -( piA_ + (1 + kappa_ / piY_) * piC_ + piG_ );
105 
106  generator_(1, 0) = piA_;
107  generator_(3, 0) = piA_;
108  generator_(0, 1) = piC_;
109  generator_(2, 1) = piC_;
110  generator_(1, 2) = piG_;
111  generator_(3, 2) = piG_;
112  generator_(0, 3) = piT_;
113  generator_(2, 3) = piT_;
114 
115  generator_(2, 0) = (1 + kappa_ / piR_) * piA_;
116  generator_(3, 1) = (1 + kappa_ / piY_) * piC_;
117  generator_(0, 2) = (1 + kappa_ / piR_) * piG_;
118  generator_(1, 3) = (1 + kappa_ / piY_) * piT_;
119 
120  // Normalization:
121  MatrixTools::scale(generator_, r_);
122 
123  // Exchangeability:
124  exchangeability_(0,0) = generator_(0,0) / piA_;
125  exchangeability_(0,1) = generator_(0,1) / piC_;
126  exchangeability_(0,2) = generator_(0,2) / piG_;
127  exchangeability_(0,3) = generator_(0,3) / piT_;
128 
129  exchangeability_(1,0) = generator_(1,0) / piA_;
130  exchangeability_(1,1) = generator_(1,1) / piC_;
131  exchangeability_(1,2) = generator_(1,2) / piG_;
132  exchangeability_(1,3) = generator_(1,3) / piT_;
133 
134  exchangeability_(2,0) = generator_(2,0) / piA_;
135  exchangeability_(2,1) = generator_(2,1) / piC_;
136  exchangeability_(2,2) = generator_(2,2) / piG_;
137  exchangeability_(2,3) = generator_(2,3) / piT_;
138 
139  exchangeability_(3,0) = generator_(3,0) / piA_;
140  exchangeability_(3,1) = generator_(3,1) / piC_;
141  exchangeability_(3,2) = generator_(3,2) / piG_;
142  exchangeability_(3,3) = generator_(3,3) / piT_;
143 
144  // Eigen values:
145  eigenValues_[0] = 0;
146  eigenValues_[1] = -r_ * (1+kappa_);
147  eigenValues_[2] = -r_ * (1+kappa_);
148  eigenValues_[3] = -r_;
149 
150  // Eigen vectors:
151  leftEigenVectors_(0,0) = piA_;
152  leftEigenVectors_(0,1) = piC_;
153  leftEigenVectors_(0,2) = piG_;
154  leftEigenVectors_(0,3) = piT_;
155 
156  leftEigenVectors_(1,0) = 0.;
157  leftEigenVectors_(1,1) = piT_ / piY_;
158  leftEigenVectors_(1,2) = 0.;
159  leftEigenVectors_(1,3) = -piT_ / piY_;
160 
161  leftEigenVectors_(2,0) = piG_ / piR_;
162  leftEigenVectors_(2,1) = 0.;
163  leftEigenVectors_(2,2) = -piG_ / piR_;
164  leftEigenVectors_(2,3) = 0.;
165 
166  leftEigenVectors_(3,0) = piA_*piY_ / piR_;
167  leftEigenVectors_(3,1) = -piC_;
168  leftEigenVectors_(3,2) = piG_*piY_ / piR_;
169  leftEigenVectors_(3,3) = -piT_;
170 
171  rightEigenVectors_(0,0) = 1.;
172  rightEigenVectors_(0,1) = 0.;
173  rightEigenVectors_(0,2) = 1.;
174  rightEigenVectors_(0,3) = 1.;
175 
176  rightEigenVectors_(1,0) = 1.;
177  rightEigenVectors_(1,1) = 1.;
178  rightEigenVectors_(1,2) = 0.;;
179  rightEigenVectors_(1,3) = -piR_ / piY_;
180 
181  rightEigenVectors_(2,0) = 1.;
182  rightEigenVectors_(2,1) = 0.;
183  rightEigenVectors_(2,2) = -piA_ / piG_;
184  rightEigenVectors_(2,3) = 1.;
185 
186  rightEigenVectors_(3,0) = 1.;
187  rightEigenVectors_(3,1) = -piC_ / piT_;
188  rightEigenVectors_(3,2) = 0.;
189  rightEigenVectors_(3,3) = -piR_ / piY_;
190 }
191 
192 /******************************************************************************/
193 
194 double F84::Pij_t(size_t i, size_t j, double d) const
195 {
196  l_ = rate_ * r_ * d;
197  exp1_ = exp(-k1_*l_);
198  exp2_ = exp(-k2_*l_);
199 
200  switch(i) {
201  //A
202  case 0 : {
203  switch(j) {
204  case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp2_; //A
205  case 1 : return piC_ * (1. - exp1_); //C
206  case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp2_; //G
207  case 3 : return piT_ * (1. - exp1_); //T, U
208  }
209  }
210  //C
211  case 1 : {
212  switch(j) {
213  case 0 : return piA_ * (1. - exp1_); //A
214  case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp2_; //C
215  case 2 : return piG_ * (1. - exp1_); //G
216  case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp2_; //T, U
217  }
218  }
219  //G
220  case 2 : {
221  switch(j) {
222  case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp2_; //A
223  case 1 : return piC_ * (1. - exp1_); //C
224  case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp2_; //G
225  case 3 : return piT_ * (1. - exp1_); //T, U
226  }
227  }
228  //T, U
229  case 3 : {
230  switch(j) {
231  case 0 : return piA_ * (1. - exp1_); //A
232  case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp2_; //C
233  case 2 : return piG_ * (1. - exp1_); //G
234  case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp2_; //T, U
235  }
236  }
237  }
238  return 0;
239 }
240 
241 /******************************************************************************/
242 
243 double F84::dPij_dt(size_t i, size_t j, double d) const
244 {
245  l_ = rate_ * r_ * d;
246  exp1_ = exp(-k1_*l_);
247  exp2_ = exp(-k2_*l_);
248 
249  switch(i) {
250  //A
251  case 0 : {
252  switch(j) {
253  case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp2_); //A
254  case 1 : return rate_ * r_ * (piC_ * exp1_); //C
255  case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp2_); //G
256  case 3 : return rate_ * r_ * (piT_ * exp1_); //T, U
257  }
258  }
259  //C
260  case 1 : {
261  switch(j) {
262  case 0 : return rate_ * r_ * (piA_ * exp1_); //A
263  case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k2_ * exp2_); //C
264  case 2 : return rate_ * r_ * (piG_ * exp1_); //G
265  case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k2_ * exp2_); //T, U
266  }
267  }
268  //G
269  case 2 : {
270  switch(j) {
271  case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp2_); //A
272  case 1 : return rate_ * r_ * (piC_ * exp1_); //C
273  case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp2_); //G
274  case 3 : return rate_ * r_ * (piT_ * exp1_); //T, U
275  }
276  }
277  //T, U
278  case 3 : {
279  switch(j) {
280  case 0 : return rate_ * r_ * (piA_ * exp1_); //A
281  case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k2_ * exp2_); //C
282  case 2 : return rate_ * r_ * (piG_ * exp1_); //G
283  case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k2_ * exp2_); //T, U
284  }
285  }
286  }
287  return 0;
288 }
289 
290 /******************************************************************************/
291 
292 double F84::d2Pij_dt2(size_t i, size_t j, double d) const
293 {
294  double r_2 = rate_ * rate_ * r_ * r_;
295  l_ = rate_ * r_ * d;
296  double k2_2 = k2_ * k2_;
297  exp1_ = exp(-k1_*l_);
298  exp2_ = exp(-k2_*l_);
299 
300  switch(i) {
301  //A
302  case 0 : {
303  switch(j) {
304  case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp2_); //A
305  case 1 : return r_2 * (piC_ * - exp1_); //C
306  case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp2_); //G
307  case 3 : return r_2 * (piT_ * - exp1_); //T, U
308  }
309  }
310  //C
311  case 1 : {
312  switch(j) {
313  case 0 : return r_2 * (piA_ * - exp1_); //A
314  case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k2_2 * exp2_); //C
315  case 2 : return r_2 * (piG_ * - exp1_); //G
316  case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k2_2 * exp2_); //T, U
317  }
318  }
319  //G
320  case 2 : {
321  switch(j) {
322  case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp2_); //A
323  case 1 : return r_2 * (piC_ * - exp1_); //C
324  case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp2_); //G
325  case 3 : return r_2 * (piT_ * - exp1_); //T, U
326  }
327  }
328  //T, U
329  case 3 : {
330  switch(j) {
331  case 0 : return r_2 * (piA_ * - exp1_); //A
332  case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k2_2 * exp2_); //C
333  case 2 : return r_2 * (piG_ * - exp1_); //G
334  case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k2_2 * exp2_); //T, U
335  }
336  }
337  }
338  return 0;
339 }
340 
341 /******************************************************************************/
342 
343 const Matrix<double> & F84::getPij_t(double d) const
344 {
345  l_ = rate_ * r_ * d;
346  exp1_ = exp(-k1_*l_);
347  exp2_ = exp(-k2_*l_);
348 
349  //A
350  p_(0, 0) = piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp2_; //A
351  p_(0, 1) = piC_ * (1. - exp1_); //C
352  p_(0, 2) = piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp2_; //G
353  p_(0, 3) = piT_ * (1. - exp1_); //T, U
354 
355  //C
356  p_(1, 0) = piA_ * (1. - exp1_); //A
357  p_(1, 1) = piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp2_; //C
358  p_(1, 2) = piG_ * (1. - exp1_); //G
359  p_(1, 3) = piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp2_; //T, U
360 
361  //G
362  p_(2, 0) = piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp2_; //A
363  p_(2, 1) = piC_ * (1. - exp1_); //C
364  p_(2, 2) = piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp2_; //G
365  p_(2, 3) = piT_ * (1. - exp1_); //T, U
366 
367  //T, U
368  p_(3, 0) = piA_ * (1. - exp1_); //A
369  p_(3, 1) = piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp2_; //C
370  p_(3, 2) = piG_ * (1. - exp1_); //G
371  p_(3, 3) = piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp2_; //T, U
372 
373  return p_;
374 }
375 
376 const Matrix<double> & F84::getdPij_dt(double d) const
377 {
378  l_ = rate_ * r_ * d;
379  exp1_ = exp(-k1_*l_);
380  exp2_ = exp(-k2_*l_);
381 
382  //A
383  p_(0, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp2_); //A
384  p_(0, 1) = rate_ * r_ * (piC_ * exp1_); //C
385  p_(0, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp2_); //G
386  p_(0, 3) = rate_ * r_ * (piT_ * exp1_); //T, U
387 
388  //C
389  p_(1, 0) = rate_ * r_ * (piA_ * exp1_); //A
390  p_(1, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k2_ * exp2_); //C
391  p_(1, 2) = rate_ * r_ * (piG_ * exp1_); //G
392  p_(1, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k2_ * exp2_); //T, U
393 
394  //G
395  p_(2, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp2_); //A
396  p_(2, 1) = rate_ * r_ * (piC_ * exp1_); //C
397  p_(2, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp2_); //G
398  p_(2, 3) = rate_ * r_ * (piT_ * exp1_); //T, U
399 
400  //T, U
401  p_(3, 0) = rate_ * r_ * (piA_ * exp1_); //A
402  p_(3, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k2_ * exp2_); //C
403  p_(3, 2) = rate_ * r_ * (piG_ * exp1_); //G
404  p_(3, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k2_ * exp2_); //T, U
405 
406  return p_;
407 }
408 
409 const Matrix<double> & F84::getd2Pij_dt2(double d) const
410 {
411  double r_2 = rate_ * rate_ * r_ * r_;
412  l_ = rate_ * r_ * d;
413  double k2_2 = k2_ * k2_;
414  exp1_ = exp(-k1_*l_);
415  exp2_ = exp(-k2_*l_);
416 
417  //A
418  p_(0, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp2_); //A
419  p_(0, 1) = r_2 * (piC_ * - exp1_); //C
420  p_(0, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp2_); //G
421  p_(0, 3) = r_2 * (piT_ * - exp1_); //T, U
422 
423  //C
424  p_(1, 0) = r_2 * (piA_ * - exp1_); //A
425  p_(1, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k2_2 * exp2_); //C
426  p_(1, 2) = r_2 * (piG_ * - exp1_); //G
427  p_(1, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k2_2 * exp2_); //T, U
428 
429  //G
430  p_(2, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp2_); //A
431  p_(2, 1) = r_2 * (piC_ * - exp1_); //C
432  p_(2, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp2_); //G
433  p_(2, 3) = r_2 * (piT_ * - exp1_); //T, U
434 
435  //T, U
436  p_(3, 0) = r_2 * (piA_ * - exp1_); //A
437  p_(3, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k2_2 * exp2_); //C
438  p_(3, 2) = r_2 * (piG_ * - exp1_); //G
439  p_(3, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k2_2 * exp2_); //T, U
440 
441  return p_;
442 }
443 
444 /******************************************************************************/
445 
446 void F84::setFreq(map<int, double>& freqs)
447 {
448  piA_ = freqs[0];
449  piC_ = freqs[1];
450  piG_ = freqs[2];
451  piT_ = freqs[3];
452  vector<string> thetas(3);
453  thetas[0] = getNamespace()+ "theta";
454  thetas[1] = getNamespace() + "theta1";
455  thetas[2] = getNamespace() + "theta2";
456  ParameterList pl = getParameters().subList(thetas);
457  pl[0].setValue(piC_ + piG_);
458  pl[1].setValue(piA_ / (piA_ + piT_));
459  pl[2].setValue(piG_ / (piC_ + piG_));
460  setParametersValues(pl);
461 }
462 
463 /******************************************************************************/
464 
double k2_
Definition: F84.h:183
RowMatrix< double > exchangeability_
The exchangeability matrix of the model, defined as . When the model is reversible, this matrix is symetric.
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:161
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
Vdouble eigenValues_
The vector of eigen values.
STL namespace.
Vdouble freq_
The vector of equilibrium frequencies.
const Matrix< double > & getPij_t(double d) const
Definition: F84.cpp:343
double d2Pij_dt2(size_t i, size_t j, double d) const
Definition: F84.cpp:292
double kappa_
Definition: F84.h:183
const Matrix< double > & getd2Pij_dt2(double d) const
Definition: F84.cpp:409
double dPij_dt(size_t i, size_t j, double d) const
Definition: F84.cpp:243
Partial implementation of the ReversibleSubstitutionModel interface.
static IntervalConstraint FREQUENCE_CONSTRAINT_SMALL
double theta1_
Definition: F84.h:183
double theta_
Definition: F84.h:183
double piC_
Definition: F84.h:183
RowMatrix< double > p_
Definition: F84.h:185
double Pij_t(size_t i, size_t j, double d) const
Definition: F84.cpp:194
double r_
Definition: F84.h:183
double l_
Definition: F84.h:184
double k1_
Definition: F84.h:183
const Matrix< double > & getdPij_dt(double d) const
Definition: F84.cpp:376
double piA_
Definition: F84.h:183
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
RowMatrix< double > generator_
The generator matrix of the model.
double piY_
Definition: F84.h:183
double theta2_
Definition: F84.h:183
double rate_
The rate of the model (default: 1). The generator (and all its vectorial components) is independent o...
double exp2_
Definition: F84.h:184
void setFreq(std::map< int, double > &)
This method is redefined to actualize the corresponding parameters piA, piT, piG and piC too...
Definition: F84.cpp:446
double piR_
Definition: F84.h:183
void updateMatrices()
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
Definition: F84.cpp:78
double piG_
Definition: F84.h:183
double exp1_
Definition: F84.h:184
double piT_
Definition: F84.h:183
F84(const NucleicAlphabet *alpha, double kappa=1., double piA=0.25, double piC=0.25, double piG=0.25, double piT=0.25)
Definition: F84.cpp:56