bpp-phyl3  3.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
RN95.cpp
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: The Bio++ Development Group
2 //
3 // SPDX-License-Identifier: CECILL-2.1
4 
6 
7 #include "RN95.h"
8 
9 // From the STL:
10 #include <cmath>
11 
12 using namespace bpp;
13 using namespace std;
14 
15 /******************************************************************************/
16 
18  shared_ptr<const NucleicAlphabet> alphabet,
19  double alpha,
20  double beta,
21  double gamma,
22  double delta,
23  double epsilon,
24  double kappa,
25  double lambda,
26  double sigma) :
28  AbstractNucleotideSubstitutionModel(alphabet, make_shared<CanonicalStateMap>(alphabet, false), "RN95."),
29  alpha_(),
30  beta_(),
31  gamma_(),
32  delta_(),
33  epsilon_(),
34  kappa_(),
35  lambda_(),
36  sigma_()
37 {
38  double f = gamma + lambda + delta + kappa;
39 
40  alpha_ = alpha / f;
41  beta_ = beta / f;
42  gamma_ = gamma / f;
43  delta_ = delta / f;
44  epsilon_ = epsilon / f;
45  kappa_ = kappa / f;
46  lambda_ = lambda / f;
47  sigma_ = sigma / f;
48 
49  double thetaR = delta_ + kappa_;
50  double kappaP = kappa_ / thetaR;
51  double gammaP = gamma_ / (1 - thetaR);
52 
53  addParameter_(new Parameter("RN95.thetaR", thetaR, Parameter::PROP_CONSTRAINT_EX));
54  addParameter_(new Parameter("RN95.gammaP", gammaP, Parameter::PROP_CONSTRAINT_EX));
55  addParameter_(new Parameter("RN95.kappaP", kappaP, Parameter::PROP_CONSTRAINT_EX));
56 
61 
62  computeFrequencies(false);
64 }
65 
66 /******************************************************************************/
67 
69 {
70  double thetaR = getParameterValue("thetaR");
71  double gammaP = getParameterValue("gammaP");
72  double kappaP = getParameterValue("kappaP");
73 
74  alpha_ = getParameterValue("alpha");
75  sigma_ = getParameterValue("sigma");
76  beta_ = getParameterValue("beta");
77  epsilon_ = getParameterValue("epsilon");
78 
79  // thetaR = delta_ + kappa_
80  kappa_ = kappaP * thetaR;
81  gamma_ = gammaP * (1 - thetaR);
82  delta_ = thetaR - kappa_;
83  lambda_ = 1 - thetaR - gamma_;
84 
85  // variables for calculation purposes
86  double thetaY = 1 - thetaR;
87 
88  auto c_1 = thetaR + sigma_ + beta_;
89  auto c_2 = 1 - c_1;
90  auto c_3 = thetaY + alpha_ + epsilon_;
91  auto c_4 = 1 - c_3;
92 
93  // stationnary frequencies
94 
95  freq_[0] = (delta_ * thetaY + epsilon_ * thetaR) / c_3;
96  freq_[1] = (sigma_ * thetaY + gamma_ * thetaR) / c_1;
97  freq_[2] = (kappa_ * thetaY + alpha_ * thetaR) / c_3;
98  freq_[3] = (beta_ * thetaY + lambda_ * thetaR) / c_1;
99 
100  // Generator matrix:
101 
102  generator_(0, 1) = gamma_;
103  generator_(0, 2) = alpha_;
104  generator_(0, 3) = lambda_;
105 
106  generator_(0, 0) = -(gamma_ + alpha_ + lambda_);
107 
108  generator_(1, 0) = delta_;
109  generator_(1, 2) = kappa_;
110  generator_(1, 3) = beta_;
111 
112  generator_(1, 1) = -(delta_ + beta_ + kappa_);
113 
114  generator_(2, 0) = epsilon_;
115  generator_(2, 1) = gamma_;
116  generator_(2, 3) = lambda_;
117 
118  generator_(2, 2) = -(gamma_ + epsilon_ + lambda_);
119 
120  generator_(3, 0) = delta_;
121  generator_(3, 1) = sigma_;
122  generator_(3, 2) = kappa_;
123 
124  generator_(3, 3) = -(delta_ + sigma_ + kappa_);
125 
126  // Normalization
127 
128  double x = 0;
129  for (size_t i = 0; i < 4; i++)
130  {
131  x += generator_(i, i) * freq_[i];
132  }
133 
134  auto r_ = isScalable() ? -1 / x : 1;
135 
136  setScale(r_);
137 
138 
139  // eigen vectors and values
140 
141  eigenValues_[0] = -r_;
142  eigenValues_[1] = -c_3 * r_;
143  eigenValues_[2] = -c_1 * r_;
144  eigenValues_[3] = 0;
145 
146  rightEigenVectors_(0, 0) = thetaY;
147  rightEigenVectors_(1, 0) = -thetaR;
148  rightEigenVectors_(2, 0) = thetaY;
149  rightEigenVectors_(3, 0) = -thetaR;
150 
151  rightEigenVectors_(0, 1) = (kappa_ - alpha_) * thetaY + alpha_ * c_4;
153  rightEigenVectors_(2, 1) = (epsilon_ - delta_) * thetaY - epsilon_ * c_4;
155 
157  rightEigenVectors_(1, 2) = (beta_ - lambda_) * thetaR - beta_ * c_2;
159  rightEigenVectors_(3, 2) = (gamma_ - sigma_) * thetaR + sigma_ * c_2;
160 
161  rightEigenVectors_(0, 3) = 1.;
162  rightEigenVectors_(1, 3) = 1.;
163  rightEigenVectors_(2, 3) = 1.;
164  rightEigenVectors_(3, 3) = 1.;
165 
166  // Need formula
167 
168  if (abs(c_2) < NumConstants::TINY() || abs(c_4) < NumConstants::TINY())
169  {
170  ApplicationTools::displayMessage("Singularity during diagonalization of RN95. Taylor series used instead.");
171 
172  isNonSingular_ = false;
173  isDiagonalizable_ = false;
175  }
176  else
177  {
178  isNonSingular_ = true;
179  isDiagonalizable_ = true;
180 
181  leftEigenVectors_(0, 0) = (delta_ - epsilon_) / c_4;
182  leftEigenVectors_(0, 1) = (sigma_ - gamma_) / c_2;
183  leftEigenVectors_(0, 2) = (kappa_ - alpha_) / c_4;
184  leftEigenVectors_(0, 3) = (beta_ - lambda_) / c_2;
185 
186  leftEigenVectors_(1, 0) = 1 / (c_3 * c_4);
187  leftEigenVectors_(1, 1) = 0;
188  leftEigenVectors_(1, 2) = -1 / (c_3 * c_4);
189  leftEigenVectors_(1, 3) = 0;
190 
191  leftEigenVectors_(2, 0) = 0;
192  leftEigenVectors_(2, 1) = -1 / (c_1 * c_2);
193  leftEigenVectors_(2, 2) = 0;
194  leftEigenVectors_(2, 3) = 1 / (c_1 * c_2);
195 
196  leftEigenVectors_(3, 0) = (epsilon_ + (delta_ - epsilon_) * thetaY) / c_3;
197  leftEigenVectors_(3, 1) = (gamma_ + (sigma_ - gamma_) * thetaY) / c_1;
198  leftEigenVectors_(3, 2) = (alpha_ + (kappa_ - alpha_) * thetaY) / c_3;
199  leftEigenVectors_(3, 3) = (lambda_ + (beta_ - lambda_) * thetaY) / c_1;
200  }
201 
202  // and the exchangeability_
203  for (unsigned int i = 0; i < size_; ++i)
204  {
205  for (unsigned int j = 0; j < size_; ++j)
206  {
207  exchangeability_(i, j) = generator_(i, j) / freq_[j];
208  }
209  }
210 }
211 
212 
213 /******************************************************************************/
214 void RN95::setFreq(map<int, double>& freqs)
215 {
216  auto thetaR = (freqs[0] + freqs[2]) / (freqs[0] + freqs[1] + freqs[2] + freqs[3]);
217  setParameterValue("thetaR", thetaR);
218  setParameterValue("gammaP", gamma_ / (1 - thetaR));
219  setParameterValue("kappaP", kappa_ / thetaR);
220 
221  updateMatrices_();
222 }
223 
224 /******************************************************************************/
Specialisation abstract class for nucleotide substitution model.
void addParameter_(Parameter *parameter)
void setParameterValue(const std::string &name, double value) override
double getParameterValue(const std::string &name) const override
RowMatrix< double > generator_
The generator matrix of the model.
bool isDiagonalizable_
boolean value for diagonalizability in R of the generator_
bool isNonSingular_
boolean value for non-singularity of rightEigenVectors_
std::vector< RowMatrix< double > > vPowGen_
vector of the powers of generator_ for Taylor development (if rightEigenVectors_ is singular).
Vdouble eigenValues_
The vector of eigen values.
RowMatrix< double > exchangeability_
The exchangeability matrix of the model, defined as . When the model is reversible,...
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
virtual bool isScalable() const
returns if model is scalable
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
void setScale(double scale)
Multiplies the current generator by the given scale.
size_t size_
The number of states.
Vdouble freq_
The vector of equilibrium frequencies.
static void displayMessage(const std::string &text)
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:168
static void Taylor(const Matrix &A, size_t p, std::vector< RowMatrix< Scalar >> &vO)
static double TINY()
static const std::shared_ptr< IntervalConstraint > PROP_CONSTRAINT_EX
static const std::shared_ptr< IntervalConstraint > R_PLUS_STAR
double beta_
Definition: RN95.h:131
double sigma_
Definition: RN95.h:131
void setFreq(std::map< int, double > &) override
Set equilibrium frequencies.
Definition: RN95.cpp:214
RN95(std::shared_ptr< const NucleicAlphabet > alphabet, double alpha=1, double beta=1, double gamma=0.25, double delta=0.25, double epsilon=1, double kappa=0.25, double lambda=0.25, double sigma=1)
Definition: RN95.cpp:17
void updateMatrices_() override
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
Definition: RN95.cpp:68
double alpha_
Definition: RN95.h:131
double epsilon_
Definition: RN95.h:131
double kappa_
Definition: RN95.h:131
double lambda_
Definition: RN95.h:131
double gamma_
Definition: RN95.h:131
double delta_
Definition: RN95.h:131
Defines the basic types of data flow nodes.
ExtendedFloat abs(const ExtendedFloat &ef)