42 #include <Bpp/Numeric/Matrix/MatrixTools.h> 52 size_t s = Q.getNumberOfRows();
53 RowMatrix<double> QL(s, s);
54 for (
size_t i = 0; i < s; i++)
56 for (
size_t j = 0; j < s; j++)
58 QL(i, j) = ((i == j) ? 0. : Q(i, j));
62 MatrixTools::fill(
m_, 0.);
63 RowMatrix<double> M2(s, s);
64 RowMatrix<double> M3(s, s);
65 RowMatrix<double> M4(s, s);
66 RowMatrix<double> M5(s, s);
67 for (
size_t n = 1; n <
cutOff_; ++n)
69 MatrixTools::fill(M2, 0.);
70 for (
size_t p = 0; p < n; ++p)
72 MatrixTools::pow(Q, p, M3);
73 MatrixTools::mult(M3, QL, M4);
74 MatrixTools::pow(Q, n - p - 1, M3);
75 MatrixTools::mult(M4, M3, M5);
76 MatrixTools::add(M2, M5);
78 MatrixTools::scale(M2, pow(length, static_cast<double>(n)) / static_cast<double>(NumTools::fact(n)));
79 MatrixTools::add(
m_, M2);
84 for (
size_t i = 0; i < s; i++)
86 for (
size_t j = 0; j < s; j++)
98 return m_(initialState, finalState);
99 if (length < 0.000001)
100 return initialState == finalState ? 0. : 1.;
105 return m_(initialState, finalState);
113 return new RowMatrix<double>(
m_);
114 if (length < 0.000001)
117 for (
size_t i = 0; i < s; i++)
119 for (
size_t j = 0; j < s; j++)
121 m_(i, j) = i == j ? 0. : 1.;
133 return new RowMatrix<double>(
m_);
void setSubstitutionModel(const SubstitutionModel *model)
Set the substitution model associated with this count, if relevent.
Interface for all substitution models.
void computeCounts(double length) const
virtual const Alphabet * getAlphabet() const =0
virtual const Matrix< double > & getPij_t(double t) const =0
virtual const Matrix< double > & getGenerator() const =0
Matrix< double > * getAllNumbersOfSubstitutions(double length, size_t type=1) const
Get the numbers of susbstitutions on a branch, for each initial and final states, and given the branc...
const SubstitutionModel * model_
double getNumberOfSubstitutions(size_t initialState, size_t finalState, double length, size_t type=1) const
Get the number of susbstitutions on a branch, given the initial and final states, and the branch leng...