42 #include <Bpp/Numeric/VectorTools.h> 43 #include <Bpp/Numeric/Matrix/MatrixTools.h> 44 #include <Bpp/Numeric/Matrix/EigenValue.h> 53 AbstractParameterAliasable(model),
55 stateMap_ (model.stateMap_),
56 nbStates_ (model.nbStates_),
57 nbRates_ (model.nbRates_),
58 rates_ (model.rates_),
59 ratesExchangeability_(model.ratesExchangeability_),
60 ratesFreq_ (model.ratesFreq_),
61 ratesGenerator_ (model.ratesGenerator_),
62 generator_ (model.generator_),
63 exchangeability_ (model.exchangeability_),
64 leftEigenVectors_ (model.leftEigenVectors_),
65 rightEigenVectors_ (model.rightEigenVectors_),
66 eigenValues_ (model.eigenValues_),
67 iEigenValues_ (model.iEigenValues_),
68 eigenDecompose_ (model.eigenDecompose_),
70 dpijt_ (model.dpijt_),
71 d2pijt_ (model.d2pijt_),
73 normalizeRateChanges_(model.normalizeRateChanges_),
74 nestedPrefix_ (model.nestedPrefix_)
80 AbstractParametrizable::operator=(model);
112 RowMatrix<double> Tmp1, Tmp2;
117 MatrixTools::MatrixTools::getId< RowMatrix<double> >(
nbStates_, Tmp1);
122 MatrixTools::mult(
rates_, Tmp1, Tmp2);
134 double scale = -VectorTools::scalar<double, double>(Tmp,
freq_);
151 for (
unsigned int i = 0; i <
nbStates_; i++)
153 RowMatrix<double> tmp =
rates_;
154 MatrixTools::scale(tmp, modelEigenValues[i]);
156 EigenValue<double> ev(tmp);
157 vector<double> values = ev.getRealEigenValues();
158 RowMatrix<double> vectors = ev.getV();
159 for (
size_t j = 0; j <
nbRates_; j++)
164 for (
unsigned int ii = 0; ii <
nbRates_; ii++)
166 double vii = vectors(ii, j);
167 for (
unsigned int jj = 0; jj <
nbStates_; jj++)
205 if (i >= (nbStates_ * nbRates_))
206 throw IndexOutOfBoundsException(
"MarkovModulatedSubstitutionModel::getInitValue", i, 0, nbStates_ * nbRates_ - 1);
207 if (state < 0 || !model_->getAlphabet()->isIntInAlphabet(state))
208 throw BadIntException(state,
"MarkovModulatedSubstitutionModel::getInitValue. Character " + model_->getAlphabet()->intToChar(state) +
" is not allowed in model.");
209 vector<int> states = model_->getAlphabet()->getAlias(state);
210 for (
size_t j = 0; j < states.size(); j++)
212 if (getAlphabetStateAsInt(i) == states[j])
222 AbstractParameterAliasable::setNamespace(prefix);
virtual void updateMatrices()
virtual const Matrix< double > & getExchangeabilityMatrix() const =0
Vdouble freq_
The vector of equilibrium frequencies.
const Matrix< double > & getd2Pij_dt2(double t) const
MarkovModulatedStateMap stateMap_
Vdouble iEigenValues_
The vector of imaginary parts of the eigen values (zero in case of reversible pmodel).
RowMatrix< double > ratesGenerator_
void setNamespace(const std::string &prefix)
virtual const Vdouble & getFrequencies() const =0
MarkovModulatedSubstitutionModel & operator=(const MarkovModulatedSubstitutionModel &model)
virtual const Vdouble & getEigenValues() const =0
RowMatrix< double > generator_
The generator matrix of the model.
Interface for reversible substitution models.
Vdouble eigenValues_
The vector of real parts of eigen values.
RowMatrix< double > pijt_
These ones are for bookkeeping:
RowMatrix< double > rates_
virtual const Matrix< double > & getGenerator() const =0
RowMatrix< double > leftEigenVectors_
The matrix made of left eigen vectors (by row).
double getInitValue(size_t i, int state) const
RowMatrix< double > dpijt_
const Matrix< double > & getdPij_dt(double t) const
bool eigenDecompose_
Tell if the eigen decomposition should be performed.
std::string nestedPrefix_
Partial implementation of the Markov-modulated class of substitution models.
ReversibleSubstitutionModel * model_
MarkovModulatedSubstitutionModel(ReversibleSubstitutionModel *model, unsigned int nbRates, bool normalizeRateChanges, const std::string &prefix)
Build a new MarkovModulatedSubstitutionModel object.
RowMatrix< double > d2pijt_
ReversibleSubstitutionModel * clone() const =0
RowMatrix< double > ratesExchangeability_
virtual const Matrix< double > & getColumnRightEigenVectors() const =0
bool normalizeRateChanges_
virtual size_t getNumberOfStates() const =0
Get the number of states.
RowMatrix< double > rightEigenVectors_
The matrix made of right eigen vectors (by column).
const Matrix< double > & getPij_t(double t) const
RowMatrix< double > exchangeability_
The exchangeability matrix of the model.