42 #include <Bpp/Text/TextTools.h> 43 #include <Bpp/Numeric/VectorTools.h> 44 #include <Bpp/Numeric/Matrix/MatrixTools.h> 45 #include <Bpp/Numeric/Matrix/EigenValue.h> 46 #include <Bpp/Numeric/NumConstants.h> 49 #include <Bpp/Seq/Container/SequenceContainerTools.h> 57 AbstractParameterAliasable(prefix),
60 size_(alpha->getSize()),
62 generator_(size_, size_),
64 exchangeability_(size_, size_),
67 d2pijt_(size_, size_),
68 eigenDecompose_(true),
71 isDiagonalizable_(false),
72 rightEigenVectors_(size_, size_),
73 isNonSingular_(false),
74 leftEigenVectors_(size_, size_),
78 for (
size_t i = 0; i <
size_; i++)
92 if (!dynamic_cast<AbstractReversibleSubstitutionModel*>(
this)) {
93 for (
size_t i = 0; i <
size_; i++)
95 for (
size_t j = 0; j <
size_; j++)
120 catch (ZeroDivisionException& e)
122 ApplicationTools::displayMessage(
"Singularity during diagonalization. Taylor series used instead.");
148 std::vector<double> vdia(
size_);
149 std::vector<double> vup(
size_ - 1);
150 std::vector<double> vlo(
size_ - 1);
152 double l =
rate_ * t;
153 for (
size_t i = 0; i <
size_; i++)
160 vup[i] = vdia[i] * s;
163 vdia[i + 1] = vdia[i];
182 double v =
rate_ * t;
189 for (
size_t i = 1; i <
vPowGen_.size(); i++)
191 s *= v /
static_cast<double>(i);
215 std::vector<double> vdia(
size_);
216 std::vector<double> vup(
size_ - 1);
217 std::vector<double> vlo(
size_ - 1);
219 double l =
rate_ * t;
220 for (
size_t i = 0; i <
size_; i++)
230 vdia[i + 1] = vdia[i];
249 double v =
rate_ * t;
256 for (
size_t i = 1; i <
vPowGen_.size(); i++)
258 s *= v /
static_cast<double>(i);
284 std::vector<double> vdia(
size_);
285 std::vector<double> vup(
size_ - 1);
286 std::vector<double> vlo(
size_ - 1);
288 double l =
rate_ * t;
289 for (
size_t i = 0; i <
size_; i++)
296 vdia[i] = NumTools::sqr(
rate_)
299 vup[i] = NumTools::sqr(
rate_)
303 vdia[i + 1] = vdia[i];
322 double v =
rate_ * t;
329 for (
size_t i = 1; i <
vPowGen_.size(); i++)
331 s *= v /
static_cast<double>(i);
352 throw IndexOutOfBoundsException(
"AbstractSubstitutionModel::getInitValue", i, 0, size_ - 1);
353 if (state < 0 || !alphabet_->isIntInAlphabet(state))
354 throw BadIntException(state,
"AbstractSubstitutionModel::getInitValue. Character " + alphabet_->intToChar(state) +
" is not allowed in model.");
355 vector<int> states = alphabet_->getAlias(state);
356 for (
size_t j = 0; j < states.size(); j++)
358 if (getAlphabetStateAsInt(i) == states[j])
368 map<int, int> counts;
369 SequenceContainerTools::getCounts(data, counts);
371 map<int, double> freqs;
373 for (
int i = 0; i < static_cast<int>(
size_); i++)
375 t += (counts[i] + pseudoCount);
377 for (
int i = 0; i < static_cast<int>(
size_); i++)
379 freqs[i] = (
static_cast<double>(counts[i]) + pseudoCount) / t;
390 for (
size_t i = 0; i <
size_; ++i)
392 freq_[i] = freqs[
static_cast<int>(i)];
404 return -VectorTools::scalar<double, double>(v,
freq_);
425 throw Exception(
"Bad value for rate: " + TextTools::toString(rate));
427 if (hasParameter(
"rate"))
428 setParameterValue(
"rate",
rate_);
435 addParameter_(
new Parameter(getNamespace() +
"rate",
rate_, &Parameter::R_PLUS_STAR));
442 RowMatrix<double> Pi;
443 MatrixTools::diag(
freq_, Pi);
446 for (
size_t i = 0; i <
size_; i++)
449 for (
size_t j = 0; j <
size_; j++)
463 for (
size_t i = 0; i <
size_; i++)
virtual void updateMatrices()
Diagonalize the matrix, and fill the eigenValues_, iEigenValues_, leftEigenVectors_ and rightEigenVe...
bool enableEigenDecomposition()
Tell if eigenValues and Vectors must be computed.
RowMatrix< double > exchangeability_
The exchangeability matrix of the model, defined as . When the model is reversible, this matrix is symetric.
virtual double getRate() const
Get the rate.
bool isNonSingular_
boolean value for non-singularity of rightEigenVectors_
virtual const Matrix< double > & getPij_t(double t) const
RowMatrix< double > d2pijt_
RowMatrix< double > pijt_
These ones are for bookkeeping:
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
Vdouble eigenValues_
The vector of eigen values.
Vdouble freq_
The vector of equilibrium frequencies.
double getInitValue(size_t i, int state) const
virtual void updateMatrices()
Compute and diagonalize the matrix, and fill the eigenValues_, leftEigenVectors_ and rightEigenVecto...
double getScale() const
Get the scalar product of diagonal elements of the generator and the frequencies vector. If the generator is normalized, then scale=1. Otherwise each element must be multiplied by 1/scale.
void setFreqFromData(const SequenceContainer &data, double pseudoCount=0)
Set equilibrium frequencies equal to the frequencies estimated from the data.
void setScale(double scale)
Multiplies the current generator by the given scale.
virtual const Matrix< double > & getdPij_dt(double t) const
virtual const Matrix< double > & getd2Pij_dt2(double t) const
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row) if rightEigenVectors_ is non-singular.
RowMatrix< double > generator_
The generator matrix of the model.
double rate_
The rate of the model (default: 1). The generator (and all its vectorial components) is independent o...
std::vector< RowMatrix< double > > vPowGen_
vector of the powers of generator_ for Taylor development (if rightEigenVectors_ is singular)...
virtual void setRate(double rate)
Set the rate of the model (must be positive).
RowMatrix< double > dpijt_
size_t size_
The size of the generator, i.e. the number of states.
Map the states of a given alphabet which have a model state.
RowMatrix< double > tmpMat_
For computational issues.
bool isDiagonalizable_
boolean value for diagonalizability in R of the generator_
Vdouble iEigenValues_
The vector of the imaginary part of the eigen values.
AbstractSubstitutionModel(const Alphabet *alpha, StateMap *stateMap, const std::string &prefix)
virtual void setFreq(std::map< int, double > &)
Set equilibrium frequencies.