bpp-phyl3  3.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
gBGC.cpp
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: The Bio++ Development Group
2 //
3 // SPDX-License-Identifier: CECILL-2.1
4 
5 #include "gBGC.h"
6 
7 // From the STL:
8 #include <cmath>
9 
10 using namespace bpp;
11 
15 
16 /******************************************************************************/
17 
19  std::shared_ptr<const NucleicAlphabet> alph,
20  std::unique_ptr<NucleotideSubstitutionModelInterface> pm,
21  double B) :
23  AbstractNucleotideSubstitutionModel(alph, pm->getStateMap(), "gBGC."),
24  model_(std::move(pm)),
25  nestedPrefix_(model_->getNamespace()),
26  B_(B)
27 {
28  model_->setNamespace("gBGC." + nestedPrefix_);
29  model_->enableEigenDecomposition(0);
30  model_->computeFrequencies(false);
31 
32  addParameters_(model_->getParameters());
33  addParameter_(new Parameter("gBGC.B", B_, std::make_shared<IntervalConstraint>(-999, 10, true, true)));
34 
35  computeFrequencies(true);
37 }
38 
39 gBGC::gBGC(const gBGC& gbgc) :
42  model_(gbgc.model_->clone()),
43  nestedPrefix_(gbgc.nestedPrefix_),
44  B_(gbgc.B_)
45 {}
46 
47 gBGC& gBGC::operator=(const gBGC& gbgc)
48 {
51  model_.reset(gbgc.model_->clone());
53  B_ = gbgc.B_;
54  return *this;
55 }
56 
58 {
60  model_->matchParametersValues(parameters);
62 }
63 
65 {
66  B_ = getParameterValue("B");
67  unsigned int i, j;
68  // Generator:
69 
70  for (i = 0; i < 4; ++i)
71  {
72  for (j = 0; j < 4; ++j)
73  {
74  generator_(i, j) = model_->Qij(i, j);
75  }
76  }
77 
78  if (B_ != 0)
79  {
80  double bp = B_ / (1 - exp(-B_));
81  double bm = B_ / (exp(B_) - 1);
82 
83  generator_(0, 0) -= (generator_(0, 1) + generator_(0, 2)) * (bp - 1);
84  generator_(1, 1) -= (generator_(1, 0) + generator_(1, 3)) * (bm - 1);
85  generator_(2, 2) -= (generator_(2, 0) + generator_(2, 3)) * (bm - 1);
86  generator_(3, 3) -= (generator_(3, 1) + generator_(3, 2)) * (bp - 1);
87 
88  generator_(0, 1) *= bp;
89  generator_(0, 2) *= bp;
90  generator_(3, 1) *= bp;
91  generator_(3, 2) *= bp;
92  generator_(1, 0) *= bm;
93  generator_(2, 0) *= bm;
94  generator_(1, 3) *= bm;
95  generator_(2, 3) *= bm;
96  }
97 
99  {
100  // calcul spectral
101 
105 
106  rightEigenVectors_ = ev.getV();
107  try
108  {
110  isNonSingular_ = true;
111  isDiagonalizable_ = true;
112 
113  for (i = 0; i < 4 && isDiagonalizable_; i++)
114  {
116  isDiagonalizable_ = false;
117  }
118 
119  // frequence stationnaire
120 
121  if (isDiagonalizable_)
122  {
123  size_t nulleigen = 0;
124  double val;
125  isNonSingular_ = false;
126 
127  while (nulleigen < 4)
128  {
129  if (abs(eigenValues_[nulleigen]) < 0.000001 && abs(iEigenValues_[nulleigen]) < 0.000001)
130  {
131  val = rightEigenVectors_(0, nulleigen);
132  i = 1;
133  while (i < 4)
134  {
135  if (abs(rightEigenVectors_(i, nulleigen) - val) > NumConstants::SMALL())
136  break;
137  i++;
138  }
139 
140  if (i < 4)
141  nulleigen++;
142  else
143  {
144  isNonSingular_ = true;
145  break;
146  }
147  }
148  else
149  nulleigen++;
150  }
151 
152  if (isNonSingular_)
153  {
154  eigenValues_[nulleigen] = 0; // to avoid approximation errors on long long branches
155  iEigenValues_[nulleigen] = 0; // to avoid approximation errors on long long branches
156 
157  for (i = 0; i < 4; i++)
158  {
159  freq_[i] = leftEigenVectors_(nulleigen, i);
160  }
161 
162  val = 0;
163  for (i = 0; i < 4; i++)
164  {
165  val += freq_[i];
166  }
167 
168  for (i = 0; i < 4; i++)
169  {
170  freq_[i] /= val;
171  }
172  }
173  else
174  {
175  ApplicationTools::displayMessage("Unable to find eigenvector for eigenvalue 1 in gBGC. Taylor series used instead.");
176  isDiagonalizable_ = false;
177  }
178  }
179  }
180  catch (ZeroDivisionException& e)
181  {
182  ApplicationTools::displayMessage("Singularity during diagonalization of gBGC in gBGC. Taylor series used instead.");
183  isNonSingular_ = false;
184  isDiagonalizable_ = false;
185  }
186 
187  if (!isNonSingular_)
188  {
189  double min = generator_(0, 0);
190  for (i = 1; i < 4; i++)
191  {
192  if (min > generator_(i, i))
193  min = generator_(i, i);
194  }
195 
196  setScale(-1 / min);
197 
198  if (vPowGen_.size() == 0)
199  vPowGen_.resize(30);
200 
201  MatrixTools::getId(4, tmpMat_); // to compute the equilibrium frequency (Q+Id)^256
204 
205  for (i = 0; i < 4; ++i)
206  {
207  freq_[i] = vPowGen_[0](0, i);
208  }
209 
211  }
212 
213  // mise a l'echelle
214 
215  normalize();
216 
217 
218  if (!isNonSingular_)
220  }
221 }
222 
223 void gBGC::setNamespace(const std::string& prefix)
224 {
226  // We also need to update the namespace of the nested model:
227  model_->setNamespace(prefix + nestedPrefix_);
228 }
Specialisation abstract class for nucleotide substitution model.
void addParameters_(const ParameterList &parameters)
void addParameter_(Parameter *parameter)
AbstractParameterAliasable & operator=(const AbstractParameterAliasable &ap)
void setNamespace(const std::string &prefix)
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_
AbstractSubstitutionModel & operator=(const AbstractSubstitutionModel &model)
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.
void normalize()
normalize the generator
Vdouble iEigenValues_
The vector of the imaginary part of the eigen values.
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
bool enableEigenDecomposition()
Tell if eigenValues and Vectors must be computed.
void setScale(double scale)
Multiplies the current generator by the given scale.
RowMatrix< double > tmpMat_
For computational issues.
virtual void fireParameterChanged(const ParameterList &parameters) override
Tells the model that a parameter value has changed.
Vdouble freq_
The vector of equilibrium frequencies.
static void displayMessage(const std::string &text)
const RowMatrix< Real > & getV() const
const std::vector< Real > & getImagEigenValues() const
const std::vector< Real > & getRealEigenValues() const
static Scalar inv(const Matrix< Scalar > &A, Matrix< Scalar > &O)
static void Taylor(const Matrix &A, size_t p, std::vector< RowMatrix< Scalar >> &vO)
static void pow(const Matrix &A, size_t p, Matrix &O)
static void add(MatrixA &A, const MatrixB &B)
static void getId(size_t n, Matrix &O)
static double TINY()
static double SMALL()
gBGC model.
Definition: gBGC.h:48
void setNamespace(const std::string &) override
Definition: gBGC.cpp:223
void updateMatrices_() override
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
Definition: gBGC.cpp:64
std::string nestedPrefix_
Definition: gBGC.h:51
std::unique_ptr< NucleotideSubstitutionModelInterface > model_
Definition: gBGC.h:50
void fireParameterChanged(const ParameterList &) override
Tells the model that a parameter value has changed.
Definition: gBGC.cpp:57
gBGC & operator=(const gBGC &gbgc)
Definition: gBGC.cpp:47
gBGC(std::shared_ptr< const NucleicAlphabet >, std::unique_ptr< NucleotideSubstitutionModelInterface >, double B=0)
Build a new gBGC substitution model.
Definition: gBGC.cpp:18
double B_
the value of the bias.
Definition: gBGC.h:56
Defines the basic types of data flow nodes.
ExtendedFloat exp(const ExtendedFloat &ef)
ExtendedFloat abs(const ExtendedFloat &ef)