bpp-phyl  2.2.0
RN95.cpp
Go to the documentation of this file.
1 //
2 // File: RN95.cpp
3 // Created by: Laurent Guéguen
4 // Created on: jeudi 24 février 2011, à 20h 42
5 //
6 
7 /*
8  Copyright or © or Copr. CNRS, (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 "RN95.h"
41 
42 #include <Bpp/Numeric/Matrix/MatrixTools.h>
43 
44 // From SeqLib:
45 #include <Bpp/Seq/Container/SequenceContainerTools.h>
46 
47 // From the STL:
48 #include <cmath>
49 
50 using namespace bpp;
51 using namespace std;
52 
53 /******************************************************************************/
54 
55 RN95::RN95(
56  const NucleicAlphabet* alphabet,
57  double alpha,
58  double beta,
59  double gamma,
60  double delta,
61  double epsilon,
62  double kappa,
63  double lambda,
64  double sigma) :
65  AbstractParameterAliasable("RN95."),
66  AbstractSubstitutionModel(alphabet, new CanonicalStateMap(alphabet, false), "RN95."),
67  alpha_(),
68  beta_(),
69  gamma_(),
70  delta_(),
71  epsilon_(),
72  kappa_(),
73  lambda_(),
74  sigma_(),
75  r_(),
76  c1_(),
77  c2_(),
78  c3_(),
79  c4_(),
80  c5_(),
81  c6_(),
82  c7_(),
83  c8_(),
84  c9_(),
85  p_(size_, size_),
86  exp1_(),
87  exp3_(),
88  exp6_(),
89  l_()
90 {
91  double f = gamma + lambda + delta + kappa;
92 
93  alpha_ = alpha / f;
94  beta_ = beta / f;
95  gamma_ = gamma / f;
96  delta_ = delta / f;
97  epsilon_ = epsilon / f;
98  kappa_ = kappa / f;
99  lambda_ = lambda / f;
100  sigma_ = sigma / f;
101 
102  double thetaR = delta_ + kappa_;
103  double thetaC = (gamma_ * thetaR + sigma_ * (1 - thetaR)) / (beta_ + sigma_ + thetaR) / (1 - thetaR);
104  double thetaG = (alpha_ * thetaR + kappa_ * (1 - thetaR)) / (alpha_ + epsilon_ + 1 - thetaR) / thetaR;
105  double kappaP = kappa_ / thetaR;
106  double gammaP = gamma_ / (1 - thetaR);
107  double alphaP = (alpha_ * (1 - thetaG) + (thetaG < kappaP ? thetaG : kappaP) * (1 - thetaR)) / (thetaG * (1 - thetaR));
108  double sigmaP = (sigma_ * (1 - thetaC) + (thetaC < gammaP ? thetaC : gammaP) * thetaR) / (thetaC * thetaR);
109  addParameter_(new Parameter("RN95.thetaR", thetaR, &Parameter::PROP_CONSTRAINT_EX));
110  addParameter_(new Parameter("RN95.thetaC", thetaC, &Parameter::PROP_CONSTRAINT_EX));
111  addParameter_(new Parameter("RN95.thetaG", thetaG, &Parameter::PROP_CONSTRAINT_EX));
112  addParameter_(new Parameter("RN95.gammaP", gammaP, &Parameter::PROP_CONSTRAINT_EX));
113  addParameter_(new Parameter("RN95.kappaP", kappaP, &Parameter::PROP_CONSTRAINT_EX));
114  addParameter_(new Parameter("RN95.alphaP", alphaP, new IntervalConstraint(1, 1, false), true));
115  addParameter_(new Parameter("RN95.sigmaP", sigmaP, new IntervalConstraint(1, 1, false), true));
116 
117  updateMatrices();
118 }
119 
120 /******************************************************************************/
121 void RN95::updateMatrices()
122 {
123  double alphaP = getParameterValue("alphaP");
124  double sigmaP = getParameterValue("sigmaP");
125  double thetaR = getParameterValue("thetaR");
126  double thetaC = getParameterValue("thetaC");
127  double thetaG = getParameterValue("thetaG");
128  double gammaP = getParameterValue("gammaP");
129  double kappaP = getParameterValue("kappaP");
130 
131  kappa_ = kappaP * thetaR;
132  gamma_ = gammaP * (1 - thetaR);
133  delta_ = thetaR - kappa_;
134  lambda_ = 1 - thetaR - gamma_;
135  alpha_ = (alphaP * (1 - thetaR) * thetaG - (thetaG < kappaP ? thetaG : kappaP) * (1 - thetaR)) / (1 - thetaG);
136  sigma_ = (sigmaP * thetaR * thetaC - (thetaC < gammaP ? thetaC : gammaP) * thetaR) / (1 - thetaC);
137  epsilon_ = (alpha_ * thetaR + kappa_ * (1 - thetaR)) / (thetaG * thetaR) - alpha_ - (1 - thetaR);
138  beta_ = (gamma_ * thetaR + sigma_ * (1 - thetaR)) / (thetaC * (1 - thetaR)) - sigma_ - thetaR;
139 
140  // stationnary frequencies
141 
142  freq_[0] = (1 - thetaG) * thetaR;
143  freq_[1] = thetaC * (1 - thetaR);
144  freq_[2] = thetaG * thetaR;
145  freq_[3] = (1 - thetaC) * (1 - thetaR);
146 
147  // Generator matrix:
148 
149  generator_(0, 1) = gamma_;
150  generator_(0, 2) = alpha_;
151  generator_(0, 3) = lambda_;
152 
153  generator_(0, 0) = -(gamma_ + alpha_ + lambda_);
154 
155  generator_(1, 0) = delta_;
156  generator_(1, 2) = kappa_;
157  generator_(1, 3) = beta_;
158 
159  generator_(1, 1) = -(delta_ + beta_ + kappa_);
160 
161  generator_(2, 0) = epsilon_;
162  generator_(2, 1) = gamma_;
163  generator_(2, 3) = lambda_;
164 
165  generator_(2, 2) = -(gamma_ + epsilon_ + lambda_);
166 
167  generator_(3, 0) = delta_;
168  generator_(3, 1) = sigma_;
169  generator_(3, 2) = kappa_;
170 
171  generator_(3, 3) = -(delta_ + sigma_ + kappa_);
172 
173  // Normalization
174 
175  double x = 0;
176  for (size_t i = 0; i < 4; i++)
177  {
178  x += generator_(i, i) * freq_[i];
179  }
180 
181  r_ = -1 / x;
182 
183  MatrixTools::scale(generator_, r_);
184  // variables for calculation purposes
185 
186  c1_ = 1;
187  c2_ = gamma_ + lambda_;
188  c3_ = alpha_ + gamma_ + epsilon_ + lambda_;
189  c4_ = kappa_ - alpha_;
190  c5_ = delta_ + kappa_;
191  c6_ = delta_ + kappa_ + sigma_ + beta_;
192  c7_ = gamma_ - sigma_;
193  c8_ = delta_ - epsilon_;
194  c9_ = lambda_ - beta_;
195 
196  // eigen vectors and values
197 
198  eigenValues_[0] = 0;
199  eigenValues_[1] = -c1_ * r_;
200  eigenValues_[2] = -c3_ * r_;
201  eigenValues_[3] = -c6_ * r_;
202 
203  rightEigenVectors_(0, 0) = 1.;
204  rightEigenVectors_(1, 0) = 1.;
205  rightEigenVectors_(2, 0) = 1.;
206  rightEigenVectors_(3, 0) = 1.;
207 
208  rightEigenVectors_(0, 1) = 1.;
209  rightEigenVectors_(1, 1) = -c5_ / c2_;
210  rightEigenVectors_(2, 1) = 1.;
211  rightEigenVectors_(3, 1) = -c5_ / c2_;
212 
213  rightEigenVectors_(0, 2) = (alpha_ * (c5_ - c3_) + kappa_ * c2_) / (delta_ * (c3_ - c2_) - epsilon_ * c5_);
214  rightEigenVectors_(1, 2) = 1.;
215  rightEigenVectors_(2, 2) = (-epsilon_ * (c5_ - c3_) - delta_ * c2_) / (delta_ * (c3_ - c2_) - epsilon_ * c5_);
216  rightEigenVectors_(3, 2) = 1.;
217 
218  rightEigenVectors_(0, 3) = 1.;
219  rightEigenVectors_(1, 3) = (beta_ * (c2_ - c6_) + lambda_ * c5_) / (gamma_ * (c6_ - c5_) - sigma_ * c2_);
220  rightEigenVectors_(2, 3) = 1;
221  rightEigenVectors_(3, 3) = (-sigma_ * (c2_ - c6_) - gamma_ * c5_) / (gamma_ * (c6_ - c5_) - sigma_ * c2_);
222 
223  // Need formula
224 
225  if (enableEigenDecomposition())
226  {
227  try
228  {
229  MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
230  isNonSingular_ = true;
231  isDiagonalizable_ = true;
232  for (size_t i = 0; i < size_ && isDiagonalizable_; i++)
233  {
234  if (abs(iEigenValues_[i]) > NumConstants::TINY())
235  isDiagonalizable_ = false;
236  }
237  }
238  catch (ZeroDivisionException& e)
239  {
240  ApplicationTools::displayMessage("Singularity during diagonalization of RN95. Taylor series used instead.");
241 
242  isNonSingular_ = false;
243  isDiagonalizable_ = false;
244  MatrixTools::Taylor(generator_, 30, vPowGen_);
245  }
246 
247  // and the exchangeability_
248  for (unsigned int i = 0; i < size_; i++)
249  for (unsigned int j = 0; j < size_; j++)
250  exchangeability_(i,j) = generator_(i,j) / freq_[j];
251  }
252 }
253 
254 /******************************************************************************/
255 double RN95::Pij_t(size_t i, size_t j, double d) const
256 {
257  l_ = rate_ * r_ * d;
258  exp1_ = exp(-c1_ * l_);
259  exp3_ = exp(-c3_ * l_);
260  exp6_ = exp(-c6_ * l_);
261 
262  switch (i)
263  {
264  {
265  // A
266  case 0: {
267  switch (j)
268  {
269  case 0: return freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
270  case 1: return freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
271  case 2: return freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
272  case 3: return freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
273  }
274  }
275  // C
276  case 1: {
277  switch (j)
278  {
279  case 0: return freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
280  case 1: return freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
281  case 2: return freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
282  case 3: return freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
283  }
284  }
285  // G
286  case 2: {
287  switch (j)
288  {
289  case 0: return freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
290  case 1: return freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
291  case 2: return freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
292  case 3: return freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
293  }
294  }
295  }
296  // T, U
297  case 3: {
298  switch (j)
299  {
300  case 0: return freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
301  case 1: return freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
302  case 2: return freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
303  case 3: return freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
304  }
305  }
306  }
307  return 0;
308 }
309 
310 /******************************************************************************/
311 double RN95::dPij_dt(size_t i, size_t j, double d) const
312 {
313  l_ = rate_ * r_ * d;
314  exp1_ = -c1_* rate_* r_* exp(-c1_ * l_);
315  exp3_ = -c3_* rate_* r_* exp(-c3_ * l_);
316  exp6_ = -c6_* rate_* r_* exp(-c6_ * l_);
317 
318  switch (i)
319  {
320  {
321  // A
322  case 0: {
323  switch (j)
324  {
325  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
326  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
327  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
328  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
329  }
330  }
331  // C
332  case 1: {
333  switch (j)
334  {
335  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
336  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
337  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
338  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
339  }
340  }
341  // G
342  case 2: {
343  switch (j)
344  {
345  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
346  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
347  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
348  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
349  }
350  }
351  }
352  // T, U
353  case 3: {
354  switch (j)
355  {
356  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
357  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
358  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
359  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
360  }
361  }
362  }
363  return 0;
364 }
365 
366 /******************************************************************************/
367 double RN95::d2Pij_dt2(size_t i, size_t j, double d) const
368 {
369  l_ = rate_ * r_ * d;
370  exp1_ = c1_ * rate_ * r_ * c1_ * rate_ * r_ * exp(-c1_ * l_);
371  exp3_ = c3_ * rate_ * r_ * c3_ * rate_ * r_ * exp(-c3_ * l_);
372  exp6_ = c6_ * rate_ * r_ * c6_ * rate_ * r_ * exp(-c6_ * l_);
373 
374  switch (i)
375  {
376  {
377  // A
378  case 0: {
379  switch (j)
380  {
381  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
382  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
383  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
384  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
385  }
386  }
387  // C
388  case 1: {
389  switch (j)
390  {
391  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
392  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
393  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
394  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
395  }
396  }
397  // G
398  case 2: {
399  switch (j)
400  {
401  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
402  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
403  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
404  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
405  }
406  }
407  }
408  // T, U
409  case 3: {
410  switch (j)
411  {
412  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
413  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
414  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
415  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
416  }
417  }
418  }
419  return 0;
420 }
421 
422 /******************************************************************************/
423 
424 const Matrix<double>& RN95::getPij_t(double d) const
425 {
426  l_ = rate_ * r_ * d;
427  exp1_ = exp(-c1_ * l_);
428  exp3_ = exp(-c3_ * l_);
429  exp6_ = exp(-c6_ * l_);
430 
431  // A
432  p_(0, 0) = freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
433  p_(0, 1) = freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
434  p_(0, 2) = freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
435  p_(0, 3) = freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
436  // C
437  p_(1, 0) = freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
438  p_(1, 1) = freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
439  p_(1, 2) = freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
440  p_(1, 3) = freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
441  // G
442  p_(2, 0) = freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
443  p_(2, 1) = freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
444  p_(2, 2) = freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
445  p_(2, 3) = freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
446  // T, U
447  p_(3, 0) = freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
448  p_(3, 1) = freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
449  p_(3, 2) = freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
450  p_(3, 3) = freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
451 
452  return p_;
453 }
454 
455 /******************************************************************************/
456 
457 const Matrix<double>& RN95::getdPij_dt(double d) const
458 {
459  l_ = rate_ * r_ * d;
460  exp1_ = -c1_* rate_* r_* exp(-c1_ * l_);
461  exp3_ = -c3_* rate_* r_* exp(-c3_ * l_);
462  exp6_ = -c6_* rate_* r_* exp(-c6_ * l_);
463 
464  // A
465  p_(0, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
466  p_(0, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
467  p_(0, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
468  p_(0, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
469  // C
470  p_(1, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
471  p_(1, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
472  p_(1, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
473  p_(1, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
474  // G
475  p_(2, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
476  p_(2, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
477  p_(2, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
478  p_(2, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
479  // T, U
480  p_(3, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
481  p_(3, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
482  p_(3, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
483  p_(3, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
484  return p_;
485 }
486 
487 /******************************************************************************/
488 
489 const Matrix<double>& RN95::getd2Pij_dt2(double d) const
490 {
491  l_ = rate_ * r_ * d;
492  exp1_ = c1_ * rate_ * r_ * c1_ * rate_ * r_ * exp(-c1_ * l_);
493  exp3_ = c3_ * rate_ * r_ * c3_ * rate_ * r_ * exp(-c3_ * l_);
494  exp6_ = c6_ * rate_ * r_ * c6_ * rate_ * r_ * exp(-c6_ * l_);
495 
496  // A
497  p_(0, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
498  p_(0, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
499  p_(0, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
500  p_(0, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
501  // C
502  p_(1, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
503  p_(1, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
504  p_(1, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
505  p_(1, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
506  // G
507  p_(2, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
508  p_(2, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
509  p_(2, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
510  p_(2, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
511  // T, U
512  p_(3, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
513  p_(3, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
514  p_(3, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
515  p_(3, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_;
516 
517  return p_;
518 }
519 
520 /******************************************************************************/
521 void RN95::setFreq(map<int, double>& freqs)
522 {
523  setParameterValue("thetaR", (freqs[0] + freqs[2]) / (freqs[0] + freqs[1] + freqs[2] + freqs[3]));
524  setParameterValue("thetaC", freqs[1] / (freqs[1] + freqs[3]));
525  setParameterValue("thetaG", freqs[2] / (freqs[0] + freqs[2]));
526 
527  updateMatrices();
528 }
529 
530 /******************************************************************************/
531 
bpp::RN95::exp6_
double exp6_
Definition: RN95.h:142
bpp::RN95::kappa_
double kappa_
Definition: RN95.h:135
bpp::AbstractSubstitutionModel::enableEigenDecomposition
bool enableEigenDecomposition()
Tell if eigenValues and Vectors must be computed.
Definition: AbstractSubstitutionModel.h:291
bpp::RN95::p_
RowMatrix< double > p_
Definition: RN95.h:141
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::RN95::c9_
double c9_
Definition: RN95.h:140
bpp::RN95::c7_
double c7_
Definition: RN95.h:140
bpp::AbstractSubstitutionModel::isNonSingular_
bool isNonSingular_
boolean value for non-singularity of rightEigenVectors_
Definition: AbstractSubstitutionModel.h:160
bpp::CanonicalStateMap
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:161
bpp::AbstractSubstitutionModel::rightEigenVectors_
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
Definition: AbstractSubstitutionModel.h:155
bpp::RN95::lambda_
double lambda_
Definition: RN95.h:135
bpp::AbstractSubstitutionModel::eigenValues_
Vdouble eigenValues_
The vector of eigen values.
Definition: AbstractSubstitutionModel.h:139
bpp::RN95::getdPij_dt
const Matrix< double > & getdPij_dt(double d) const
Definition: RN95.cpp:457
bpp::RN95::d2Pij_dt2
double d2Pij_dt2(size_t i, size_t j, double d) const
Definition: RN95.cpp:367
std
STL namespace.
bpp::AbstractSubstitutionModel::freq_
Vdouble freq_
The vector of equilibrium frequencies.
Definition: AbstractSubstitutionModel.h:113
bpp::RN95::beta_
double beta_
Definition: RN95.h:135
bpp::RN95::c1_
double c1_
Definition: RN95.h:140
bpp::AbstractSubstitutionModel
Partial implementation of the SubstitutionModel interface.
Definition: AbstractSubstitutionModel.h:78
bpp::RN95::c2_
double c2_
Definition: RN95.h:140
bpp::RN95::epsilon_
double epsilon_
Definition: RN95.h:135
bpp::RN95::delta_
double delta_
Definition: RN95.h:135
bpp::RN95::gamma_
double gamma_
Definition: RN95.h:135
bpp::RN95::c6_
double c6_
Definition: RN95.h:140
bpp::RN95::c5_
double c5_
Definition: RN95.h:140
bpp::RN95::Pij_t
double Pij_t(size_t i, size_t j, double d) const
Definition: RN95.cpp:255
bpp::RN95::getd2Pij_dt2
const Matrix< double > & getd2Pij_dt2(double d) const
Definition: RN95.cpp:489
bpp::RN95::alpha_
double alpha_
Definition: RN95.h:135
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::RN95::sigma_
double sigma_
Definition: RN95.h:135
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::RN95::RN95
RN95(const NucleicAlphabet *alphabet, double alpha=1, double beta=1, double gamma=1, double delta=1, double epsilon=1, double kappa=1, double lambda=1, double sigma=1)
Definition: RN95.cpp:55
bpp::AbstractSubstitutionModel::vPowGen_
std::vector< RowMatrix< double > > vPowGen_
vector of the powers of generator_ for Taylor development (if rightEigenVectors_ is singular)...
Definition: AbstractSubstitutionModel.h:172
bpp::RN95::updateMatrices
void updateMatrices()
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
Definition: RN95.cpp:121
bpp::RN95::c4_
double c4_
Definition: RN95.h:140
bpp::RN95::c8_
double c8_
Definition: RN95.h:140
bpp::AbstractSubstitutionModel::size_
size_t size_
The size of the generator, i.e. the number of states.
Definition: AbstractSubstitutionModel.h:96
bpp::RN95::exp1_
double exp1_
Definition: RN95.h:142
bpp::RN95::r_
double r_
Definition: RN95.h:136
bpp::RN95::getPij_t
const Matrix< double > & getPij_t(double d) const
Definition: RN95.cpp:424
bpp::RN95::c3_
double c3_
Definition: RN95.h:140
bpp::RN95::dPij_dt
double dPij_dt(size_t i, size_t j, double d) const
Definition: RN95.cpp:311
bpp::RN95::exp3_
double exp3_
Definition: RN95.h:142
bpp::AbstractSubstitutionModel::isDiagonalizable_
bool isDiagonalizable_
boolean value for diagonalizability in R of the generator_
Definition: AbstractSubstitutionModel.h:150
bpp::AbstractSubstitutionModel::iEigenValues_
Vdouble iEigenValues_
The vector of the imaginary part of the eigen values.
Definition: AbstractSubstitutionModel.h:144
bpp::RN95::l_
double l_
Definition: RN95.h:142
bpp::RN95::setFreq
void setFreq(std::map< int, double > &)
Set equilibrium frequencies.
Definition: RN95.cpp:521