bpp-phyl3  3.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
POMO.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 "POMO.h"
6 
7 using namespace bpp;
8 using namespace std;
9 
10 /******************************************************************************/
11 
13  std::shared_ptr<const AllelicAlphabet> allAlph,
14  std::unique_ptr<SubstitutionModelInterface> pmodel,
15  std::unique_ptr<FrequencySetInterface> pfitness) :
17  AbstractSubstitutionModel(allAlph, make_shared<CanonicalStateMap>(allAlph, false), "POMO."),
18  nbAlleles_(allAlph->getNbAlleles()),
19  pmodel_(std::move(pmodel)),
20  pfitness_(std::move(pfitness))
21 {
22  const auto& alph = allAlph->getStateAlphabet();
23 
24  if (alph->getAlphabetType() != pmodel_->alphabet().getAlphabetType())
25  throw AlphabetMismatchException("POMO mismatch state alphabet for model.", alph.get(), &pmodel_->alphabet());
26 
27  if (pfitness_ && alph->getAlphabetType() != pfitness_->alphabet().getAlphabetType())
28  throw AlphabetMismatchException("POMO mismatch state alphabet for fitness.", alph.get(), &pfitness_->alphabet());
29 
30  if (pfitness_)
31  pfitness_->setNamespace("POMO.fit_" + pfitness_->getNamespace());
32  pmodel_->setNamespace("POMO." + pmodel_->getNamespace());
33 
34  pmodel_->enableEigenDecomposition(pmodel_->computeFrequencies()); // enable eigen if needed for pmodel->freq_
35  if (pfitness_)
36  addParameters_(pfitness_->getParameters());
37  addParameters_(pmodel_->getParameters());
38 
39  computeFrequencies(false); // freq_ analytically defined
41 }
42 
44 {
45  auto nbStates = pmodel_->getNumberOfStates();
46  auto nbAlleles = allelicAlphabet().getNbAlleles();
47 
48  const auto& Q = pmodel_->generator();
49 
50  const Vdouble* fit = pfitness_ ? &pfitness_->getFrequencies() : 0;
51 
52  // Per couple of alleles
53 
54  // position of the block of alleles
55  size_t nbloc = nbStates;
56 
57  // for all couples starting with i
58  for (size_t i = 0; i < nbStates; i++)
59  {
60  double phi_i = fit ? (*fit)[i] : 1. / (double)nbStates;
61  // then for all couples ending with j
62  for (size_t j = i + 1; j < nbStates; j++)
63  {
64  double phi_j = fit ? (*fit)[j] : 1. / (double)nbStates;
65 
66  // mutations
67  generator_(i, nbloc) = nbAlleles * Q(i, j);
68  generator_(j, nbloc + nbAlleles - 2) = nbAlleles * Q(j, i);
69 
70  // drift + selection
71  if (std::abs(generator_(i, nbloc)) < NumConstants::TINY() &&
72  std::abs(generator_(j, nbloc + nbAlleles - 2)) < NumConstants::TINY()) // No mutation between states
73  {
74  for (size_t a = 1; a < nbAlleles; a++)
75  {
76  generator_(nbloc + a - 1, (a > 1 ? nbloc + a - 2 : i)) = 0;
77  generator_(nbloc + a - 1, (a < nbAlleles - 1 ? nbloc + a : j)) = 0;
78  }
79  }
80  else
81  for (size_t a = 1; a < nbAlleles; a++)
82  {
83  double rap = (double)(a * (nbAlleles - a )) / ((double)a * phi_i + (double)(nbAlleles - a) * phi_j);
84 
85  // drift towards i: a -> a+1
86  generator_(nbloc + a - 1, (a > 1 ? nbloc + a - 2 : i)) = phi_i * rap;
87 
88  // drift towards j: nbAlleles - a -> nbAlleles -a +1
89  generator_(nbloc + a - 1, (a < nbAlleles - 1 ? nbloc + a : j)) = phi_j * rap;
90  }
91  // setting for the next couple of alleles
92  nbloc += nbAlleles - 1;
93  }
94  }
95 
96  setDiagonal();
97 
98  // Stationnary distribution
99  // and variables used for rate of substitutions
100 
101  const auto& pFreq = pmodel_->getFrequencies();
102 
103  Vdouble pN(nbStates), pNm(nbStates), pNm2(nbStates), p(nbStates);
104  for (size_t i = 0; i < nbStates; i++)
105  {
106  pNm2[i] = fit ? std::pow((*fit)[i], nbAlleles - 2) : 1;
107  pNm[i] = fit ? pNm2[i] * (*fit)[i] : 1;
108  pN[i] = fit ? pNm[i] * (*fit)[i] : 1;
109  freq_[i] = pFreq[i] * pNm[i];
110  }
111 
112  size_t k = nbStates;
113  double sden = 0;
114  double snum = 0;
115 
116  for (size_t i = 0; i < nbStates - 1; ++i)
117  {
118  double phi_i = fit ? (*fit)[i] : 1. / (double)nbStates;
119  // then for all couples ending with j
120  for (size_t j = i + 1; j < nbStates; j++)
121  {
122  double phi_j = fit ? (*fit)[j] : 1. / (double)nbStates;
123 
124  auto mu = Q(i, j) * pFreq[i]; // alleles i & j
125 
126  double rfreq = pNm2[i];
127 
128  double rat = fit ? (*fit)[j] / (*fit)[i] : 1;
129 
130  for (size_t n = 1; n < nbAlleles; n++) // i_{N-n}j_{n}
131  {
132  freq_[k++] = mu * rfreq * ((int)(nbAlleles - n) * phi_i + (int)n * phi_j) * nbAlleles / ((int)n * (int)(nbAlleles - n));
133  rfreq *= rat;
134  }
135 
136  snum += 2 * mu * pNm2[i] * pN[j] * (pN[i] != pN[j] ? (phi_i - phi_j) / (pN[i] - pN[j]) : (1. / nbAlleles));
137  sden += mu * 2 * pNm[i] + ((phi_i + phi_j) *
138  ((phi_i != phi_j) ? ((pNm[i] - pNm[j]) / (phi_i - phi_j)) : (nbAlleles - 1)));
139  }
140  }
141 
142  // stationary freq
143  double x = VectorTools::sum(freq_);
144  freq_ /= x;
145 
146 
147  // Specific normalization in numbers of substitutions (from appendix D of Genetics. 2019 Aug; 212(4): 1321–1336.)
148  // s is the probability of substitution on a duration of 1 generation (ie the actual scale time of the model).
149  double s = snum / sden;
150 
151  // And everything for exponential
153 
154  setScale(1 / s);
155 }
156 
157 
159 {
160  if (pfitness_)
161  pfitness_->matchParametersValues(parameters);
162  pmodel_->matchParametersValues(parameters);
163 
164  updateMatrices_();
165 }
166 
167 
168 void POMO::setFreq(map<int, double>& frequencies)
169 {
170  // pfreqset_->setFrequenciesFromAlphabetStatesFrequencies(frequencies);
171  // matchParametersValues(pfreqset_->getParameters());
172 }
void addParameters_(const ParameterList &parameters)
RowMatrix< double > generator_
The generator matrix of the model.
void setDiagonal()
set the diagonal of the generator such that sum on each line equals 0.
virtual void updateMatrices_()
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
void setScale(double scale)
Multiplies the current generator by the given scale.
Vdouble freq_
The vector of equilibrium frequencies.
unsigned int getNbAlleles() const
This class implements a state map where all resolved states are modeled.
Definition: StateMap.h:168
static double TINY()
const AllelicAlphabet & allelicAlphabet() const
Definition: POMO.h:141
void fireParameterChanged(const ParameterList &parameters) override
Tells the model that a parameter value has changed.
Definition: POMO.cpp:158
std::unique_ptr< FrequencySetInterface > pfitness_
Definition: POMO.h:92
std::unique_ptr< SubstitutionModelInterface > pmodel_
Definition: POMO.h:91
void updateMatrices_() override
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
Definition: POMO.cpp:43
void setFreq(std::map< int, double > &frequencies) override
Set equilibrium frequencies.
Definition: POMO.cpp:168
POMO(std::shared_ptr< const AllelicAlphabet > allAlph, std::unique_ptr< SubstitutionModelInterface > pmodel, std::unique_ptr< FrequencySetInterface > pfitness)
Build a POMO instance.
Definition: POMO.cpp:12
static T sum(const std::vector< T > &v1)
Defines the basic types of data flow nodes.
std::vector< double > Vdouble
ExtendedFloat pow(const ExtendedFloat &ef, double exp)
ExtendedFloat abs(const ExtendedFloat &ef)