13 std::shared_ptr<const AllelicAlphabet> allAlph,
14 std::unique_ptr<SubstitutionModelInterface> pmodel,
15 std::unique_ptr<FrequencySetInterface> pfitness) :
18 nbAlleles_(allAlph->getNbAlleles()),
19 pmodel_(std::move(pmodel)),
20 pfitness_(std::move(pfitness))
22 const auto& alph = allAlph->getStateAlphabet();
24 if (alph->getAlphabetType() !=
pmodel_->alphabet().getAlphabetType())
45 auto nbStates =
pmodel_->getNumberOfStates();
48 const auto& Q =
pmodel_->generator();
55 size_t nbloc = nbStates;
58 for (
size_t i = 0; i < nbStates; i++)
60 double phi_i = fit ? (*fit)[i] : 1. / (double)nbStates;
62 for (
size_t j = i + 1; j < nbStates; j++)
64 double phi_j = fit ? (*fit)[j] : 1. / (double)nbStates;
68 generator_(j, nbloc + nbAlleles - 2) = nbAlleles * Q(j, i);
74 for (
size_t a = 1; a < nbAlleles; a++)
76 generator_(nbloc + a - 1, (a > 1 ? nbloc + a - 2 : i)) = 0;
77 generator_(nbloc + a - 1, (a < nbAlleles - 1 ? nbloc + a : j)) = 0;
81 for (
size_t a = 1; a < nbAlleles; a++)
83 double rap = (double)(a * (nbAlleles - a )) / ((
double)a * phi_i + (double)(nbAlleles - a) * phi_j);
86 generator_(nbloc + a - 1, (a > 1 ? nbloc + a - 2 : i)) = phi_i * rap;
89 generator_(nbloc + a - 1, (a < nbAlleles - 1 ? nbloc + a : j)) = phi_j * rap;
92 nbloc += nbAlleles - 1;
101 const auto& pFreq =
pmodel_->getFrequencies();
103 Vdouble pN(nbStates), pNm(nbStates), pNm2(nbStates), p(nbStates);
104 for (
size_t i = 0; i < nbStates; i++)
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];
116 for (
size_t i = 0; i < nbStates - 1; ++i)
118 double phi_i = fit ? (*fit)[i] : 1. / (double)nbStates;
120 for (
size_t j = i + 1; j < nbStates; j++)
122 double phi_j = fit ? (*fit)[j] : 1. / (double)nbStates;
124 auto mu = Q(i, j) * pFreq[i];
126 double rfreq = pNm2[i];
128 double rat = fit ? (*fit)[j] / (*fit)[i] : 1;
130 for (
size_t n = 1; n < nbAlleles; n++)
132 freq_[k++] = mu * rfreq * ((int)(nbAlleles - n) * phi_i + (int)n * phi_j) * nbAlleles / ((int)n * (
int)(nbAlleles - n));
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)));
149 double s = snum / sden;
161 pfitness_->matchParametersValues(parameters);
162 pmodel_->matchParametersValues(parameters);
void addParameters_(const ParameterList ¶meters)
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.
bool computeFrequencies() const
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.
const AllelicAlphabet & allelicAlphabet() const
void fireParameterChanged(const ParameterList ¶meters) override
Tells the model that a parameter value has changed.
std::unique_ptr< FrequencySetInterface > pfitness_
std::unique_ptr< SubstitutionModelInterface > pmodel_
void updateMatrices_() override
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
void setFreq(std::map< int, double > &frequencies) override
Set equilibrium frequencies.
POMO(std::shared_ptr< const AllelicAlphabet > allAlph, std::unique_ptr< SubstitutionModelInterface > pmodel, std::unique_ptr< FrequencySetInterface > pfitness)
Build a POMO instance.
Defines the basic types of data flow nodes.
std::vector< double > Vdouble
ExtendedFloat pow(const ExtendedFloat &ef, double exp)
ExtendedFloat abs(const ExtendedFloat &ef)