bpp-phyl  2.1.0
 All Classes Namespaces Files Functions Variables Friends Pages
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, "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  try
226  {
227  MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
228  isNonSingular_ = true;
229  isDiagonalizable_ = true;
230  for (size_t i = 0; i < size_ && isDiagonalizable_; i++)
231  {
232  if (abs(iEigenValues_[i]) > NumConstants::TINY())
233  isDiagonalizable_ = false;
234  }
235  }
236  catch (ZeroDivisionException& e)
237  {
238  ApplicationTools::displayMessage("Singularity during diagonalization. Taylor series used instead.");
239 
240  isNonSingular_ = false;
241  isDiagonalizable_ = false;
242  MatrixTools::Taylor(generator_, 30, vPowGen_);
243  }
244 
245  // and the exchangeability_
246  for (unsigned int i = 0; i < size_; i++)
247  for (unsigned int j = 0; j < size_; j++)
248  exchangeability_(i,j) = generator_(i,j) / freq_[j];
249 
250 }
251 
252 /******************************************************************************/
253 double RN95::Pij_t(int i, int j, double d) const
254 {
255  l_ = rate_ * r_ * d;
256  exp1_ = exp(-c1_ * l_);
257  exp3_ = exp(-c3_ * l_);
258  exp6_ = exp(-c6_ * l_);
259 
260  switch (i)
261  {
262  {
263  // A
264  case 0: {
265  switch (j)
266  {
267  case 0: return freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
268  case 1: return freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
269  case 2: return freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
270  case 3: return freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
271  }
272  }
273  // C
274  case 1: {
275  switch (j)
276  {
277  case 0: return freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
278  case 1: return freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
279  case 2: return freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
280  case 3: return freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
281  }
282  }
283  // G
284  case 2: {
285  switch (j)
286  {
287  case 0: return freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
288  case 1: return freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
289  case 2: return freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
290  case 3: return freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
291  }
292  }
293  }
294  // T, U
295  case 3: {
296  switch (j)
297  {
298  case 0: return freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
299  case 1: return freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
300  case 2: return freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
301  case 3: return freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
302  }
303  }
304  }
305  return 0;
306 }
307 
308 /******************************************************************************/
309 double RN95::dPij_dt(int i, int j, double d) const
310 {
311  l_ = rate_ * r_ * d;
312  exp1_ = -c1_* rate_* r_* exp(-c1_ * l_);
313  exp3_ = -c3_* rate_* r_* exp(-c3_ * l_);
314  exp6_ = -c6_* rate_* r_* exp(-c6_ * l_);
315 
316  switch (i)
317  {
318  {
319  // A
320  case 0: {
321  switch (j)
322  {
323  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
324  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
325  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
326  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
327  }
328  }
329  // C
330  case 1: {
331  switch (j)
332  {
333  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
334  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
335  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
336  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
337  }
338  }
339  // G
340  case 2: {
341  switch (j)
342  {
343  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
344  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
345  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
346  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
347  }
348  }
349  }
350  // T, U
351  case 3: {
352  switch (j)
353  {
354  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
355  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
356  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
357  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
358  }
359  }
360  }
361  return 0;
362 }
363 
364 /******************************************************************************/
365 double RN95::d2Pij_dt2(int i, int j, double d) const
366 {
367  l_ = rate_ * r_ * d;
368  exp1_ = c1_ * rate_ * r_ * c1_ * rate_ * r_ * exp(-c1_ * l_);
369  exp3_ = c3_ * rate_ * r_ * c3_ * rate_ * r_ * exp(-c3_ * l_);
370  exp6_ = c6_ * rate_ * r_ * c6_ * rate_ * r_ * exp(-c6_ * l_);
371 
372  switch (i)
373  {
374  {
375  // A
376  case 0: {
377  switch (j)
378  {
379  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
380  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
381  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
382  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
383  }
384  }
385  // C
386  case 1: {
387  switch (j)
388  {
389  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
390  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
391  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
392  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
393  }
394  }
395  // G
396  case 2: {
397  switch (j)
398  {
399  case 0: return -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
400  case 1: return c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
401  case 2: return -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
402  case 3: return c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
403  }
404  }
405  }
406  // T, U
407  case 3: {
408  switch (j)
409  {
410  case 0: return c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
411  case 1: return -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
412  case 2: return c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
413  case 3: return -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
414  }
415  }
416  }
417  return 0;
418 }
419 
420 /******************************************************************************/
421 
422 const Matrix<double>& RN95::getPij_t(double d) const
423 {
424  l_ = rate_ * r_ * d;
425  exp1_ = exp(-c1_ * l_);
426  exp3_ = exp(-c3_ * l_);
427  exp6_ = exp(-c6_ * l_);
428 
429  // A
430  p_(0, 0) = freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
431  p_(0, 1) = freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
432  p_(0, 2) = freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
433  p_(0, 3) = freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
434  // C
435  p_(1, 0) = freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
436  p_(1, 1) = freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
437  p_(1, 2) = freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
438  p_(1, 3) = freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
439  // G
440  p_(2, 0) = freq_[0] - c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
441  p_(2, 1) = freq_[1] + c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
442  p_(2, 2) = freq_[2] - c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
443  p_(2, 3) = freq_[3] + c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
444  // T, U
445  p_(3, 0) = freq_[0] + c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
446  p_(3, 1) = freq_[1] - c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
447  p_(3, 2) = freq_[2] + c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
448  p_(3, 3) = freq_[3] - c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
449 
450  return p_;
451 }
452 
453 /******************************************************************************/
454 
455 const Matrix<double>& RN95::getdPij_dt(double d) const
456 {
457  l_ = rate_ * r_ * d;
458  exp1_ = -c1_* rate_* r_* exp(-c1_ * l_);
459  exp3_ = -c3_* rate_* r_* exp(-c3_ * l_);
460  exp6_ = -c6_* rate_* r_* exp(-c6_ * l_);
461 
462  // A
463  p_(0, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
464  p_(0, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
465  p_(0, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
466  p_(0, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
467  // C
468  p_(1, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
469  p_(1, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
470  p_(1, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
471  p_(1, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
472  // G
473  p_(2, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
474  p_(2, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
475  p_(2, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
476  p_(2, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
477  // T, U
478  p_(3, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
479  p_(3, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
480  p_(3, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
481  p_(3, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // T
482  return p_;
483 }
484 
485 /******************************************************************************/
486 
487 const Matrix<double>& RN95::getd2Pij_dt2(double d) const
488 {
489  l_ = rate_ * r_ * d;
490  exp1_ = c1_ * rate_ * r_ * c1_ * rate_ * r_ * exp(-c1_ * l_);
491  exp3_ = c3_ * rate_ * r_ * c3_ * rate_ * r_ * exp(-c3_ * l_);
492  exp6_ = c6_ * rate_ * r_ * c6_ * rate_ * r_ * exp(-c6_ * l_);
493 
494  // A
495  p_(0, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // A
496  p_(0, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
497  p_(0, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (alpha_ * (c3_ - c1_) - c2_ * c4_) / (c3_ * (c3_ - c1_)) * exp3_; // G
498  p_(0, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
499  // C
500  p_(1, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
501  p_(1, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // C
502  p_(1, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
503  p_(1, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (beta_ * (c6_ - c1_) - c5_ * c9_) / (c6_ * (c6_ - c1_)) * exp6_; // T
504  // G
505  p_(2, 0) = -c2_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // A
506  p_(2, 1) = c2_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // C
507  p_(2, 2) = -c2_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (c2_ * c8_ - epsilon_ * (c3_ - c1_)) / (c3_ * (c3_ - c1_)) * exp3_; // G
508  p_(2, 3) = c2_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (lambda_ * sigma_ - gamma_ * beta_) / (c6_ * (c6_ - c1_)) * exp6_; // T, U
509  // T, U
510  p_(3, 0) = c5_ * c8_ / (c1_ * (c3_ - c1_)) * exp1_ + (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // A
511  p_(3, 1) = -c5_ * c7_ / (c1_ * (c6_ - c1_)) * exp1_ + (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_; // C
512  p_(3, 2) = c5_ * c4_ / (c1_ * (c3_ - c1_)) * exp1_ - (epsilon_ * kappa_ - delta_ * alpha_) / (c3_ * (c3_ - c1_)) * exp3_; // G
513  p_(3, 3) = -c5_ * c9_ / (c1_ * (c6_ - c1_)) * exp1_ - (c5_ * c7_ - sigma_ * (c6_ - c1_)) / (c6_ * (c6_ - c1_)) * exp6_;
514 
515  return p_;
516 }
517 
518 /******************************************************************************/
519 void RN95::setFreq(map<int, double>& freqs)
520 {
521  setParameterValue("thetaR", (freqs[0] + freqs[2]) / (freqs[0] + freqs[1] + freqs[2] + freqs[3]));
522  setParameterValue("thetaC", freqs[1] / (freqs[1] + freqs[3]));
523  setParameterValue("thetaG", freqs[2] / (freqs[0] + freqs[2]));
524 
525  updateMatrices();
526 }
527 
528 /******************************************************************************/
529