42 #include "Bpp/Numeric/Matrix/MatrixTools.h" 43 #include "Bpp/Numeric/NumTools.h" 55 nbStates_(model->getNumberOfStates()),
56 bMatrices_(reg->getNumberOfSubstitutionTypes()),
58 s_(reg->getNumberOfSubstitutionTypes()),
60 counts_(reg->getNumberOfSubstitutionTypes()),
65 throw Exception(
"UniformizationSubstitutionCount (constructor): alphabets do not match between register and model.");
72 for (
unsigned int i = 0; i <
nbStates_; ++i) {
79 throw Exception(
"UniformizationSubstitutionCount::UniformizationSubstitutionCount The maximum diagonal values of generator is above 10000. Abort, chose another mapping method");
87 size_t nbTypes =
register_->getNumberOfSubstitutionTypes();
97 for (
size_t i = 0; i <
register_->getNumberOfSubstitutionTypes(); ++i) {
108 if (i > 0 && k != j) {
122 double lam =
miu_ * length;
126 MatrixTools::scale(R, 1. /
miu_);
127 MatrixTools::add(R, I);
132 size_t nMax =
static_cast<size_t>(ceil(4 + 6 * sqrt(lam) + lam));
137 for (
size_t i = 1; i < nMax + 1; ++i)
140 for (
size_t i = 0; i <
register_->getNumberOfSubstitutionTypes(); ++i) {
141 s_[i].resize(nMax + 1);
144 for (
size_t l = 1; l < nMax + 1; ++l) {
145 MatrixTools::mult(R,
s_[i][l - 1],
s_[i][l]);
147 MatrixTools::add(
s_[i][l], tmp);
149 MatrixTools::fill(
counts_[i], 0);
150 for (
size_t l = 0; l < nMax + 1; ++l) {
153 double logF =
static_cast<double>(l + 1) * log(lam) - lam - log(
miu_) - NumTools::logFact(static_cast<double>(l + 1));
154 MatrixTools::scale(tmp, exp(logF));
155 MatrixTools::add(
counts_[i], tmp);
162 for (
size_t i = 0; i <
register_->getNumberOfSubstitutionTypes(); i++) {
170 counts_[i](j, k) *=
weights_->getIndex(supportedStates[j], supportedStates[k]);
181 throw Exception(
"UniformizationSubstitutionCount::getAllNumbersOfSubstitutions. Negative branch length: " + TextTools::toString(length) +
".");
187 return new RowMatrix<double>(
counts_[type - 1]);
195 throw Exception(
"UniformizationSubstitutionCount::getNumbersOfSubstitutions. Negative branch length: " + TextTools::toString(length) +
".");
201 return counts_[type - 1](initialState, finalState);
209 throw Exception(
"UniformizationSubstitutionCount::getNumbersOfSubstitutions. Negative branch length: " + TextTools::toString(length) +
".");
217 v[t] =
counts_[t](initialState, finalState);
228 throw Exception(
"UniformizationSubstitutionCount::setSubstitutionModel: alphabets do not match between register and model.");
240 for (
unsigned int i = 0; i <
nbStates_; ++i) {
247 throw Exception(
"UniformizationSubstitutionCount::setSubstitutionModel The maximum diagonal values of generator is above 10000. Abort, chose another mapping method.");
259 throw Exception(
"UniformizationSubstitutionCount::substitutionRegisterHasChanged: alphabets do not match between register and model.");
274 if (
weights_->getAlphabet()->getAlphabetType() !=
register_->getAlphabet()->getAlphabetType())
275 throw Exception(
"UniformizationSubstitutionCount::weightsHaveChanged. Incorrect alphabet type.");
Interface for all substitution models.
std::auto_ptr< SubstitutionRegister > register_
virtual const std::vector< int > & getAlphabetStates() const =0
Basic implementation of the the SubstitutionCount interface.
The SubstitutionRegister interface.
virtual const Alphabet * getAlphabet() const =0
const AlphabetIndex2 * weights_
virtual const Matrix< double > & getPij_t(double t) const =0
virtual size_t getNumberOfSubstitutionTypes() const
Short cut function, equivalent to getSubstitutionRegister().getNumberOfSubstitutionTypes().
virtual double Qij(size_t i, size_t j) const =0
virtual const Matrix< double > & getGenerator() const =0
Partial implementation of the WeightedSubstitutionCount interface.
virtual const Alphabet * getAlphabet() const =0