42 #include "Bpp/Numeric/Matrix/MatrixTools.h" 53 nbStates_(model->getNumberOfStates()),
54 jMat_(nbStates_, nbStates_),
55 v_(nbStates_, nbStates_),
56 vInv_(nbStates_, nbStates_),
57 lambda_(nbStates_, nbStates_),
58 bMatrice_(nbStates_, nbStates_),
59 insideProduct_(nbStates_, nbStates_),
60 rewards_(nbStates_, nbStates_),
65 throw Exception(
"DecompositionReward (constructor): alphabets do not match between alphabet index and model.");
66 if (!dynamic_cast<const ReversibleSubstitutionModel*>(model))
67 throw Exception(
"DecompositionReward::DecompositionReward. Only works with declared reversible models for now.");
112 vector<double> expLam = VectorTools::exp(lambda * t);
113 for (
unsigned int i = 0; i <
nbStates_; ++i) {
114 for (
unsigned int j = 0; j <
nbStates_; ++j) {
115 double dd = lambda[i] - lambda[j];
117 result(i, j) = t * expLam[i];
119 result(i, j) = (expLam[i] - expLam[j]) / dd;
132 MatrixTools::mult(
v_, tmp1, tmp2);
151 throw Exception(
"DecompositionReward::getAllRewards. Negative branch length: " + TextTools::toString(length) +
".");
157 return new RowMatrix<double>(
rewards_);
165 throw Exception(
"DecompositionReward::getRewards. Negative branch length: " + TextTools::toString(length) +
".");
171 return rewards_(initialState, finalState);
180 throw Exception(
"DecompositionReward::setSubstitutionModel. Only works with reversible models for now.");
184 throw Exception(
"DecompositionReward::setSubstitutionModel: alphabets do not match between alphabet index and model.");
204 throw Exception(
"DecompositionReward::AlphabetIndexHasChanged: alphabets do not match between alphbaet index and model.");
RowMatrix< double > bMatrice_
RowMatrix< double > jMat_
Interface for all substitution models.
void alphabetIndexHasChanged()
AlphabetIndex1 * alphIndex_
RowMatrix< double > rewards_
virtual const std::vector< int > & getAlphabetStates() const =0
virtual const Alphabet * getAlphabet() const =0
Basic implementation of the the Reward interface.
const AlphabetIndex1 * getAlphabetIndex() const
DecompositionReward(const SubstitutionModel *model, AlphabetIndex1 *alphIndex)
virtual const Vdouble & getEigenValues() const =0
virtual const Matrix< double > & getPij_t(double t) const =0
Interface for reversible substitution models.
double getReward(size_t initialState, size_t finalState, double length) const
Get the reward of susbstitutions on a branch, given the initial and final states, and the branch leng...
const ReversibleSubstitutionModel * model_
void computeRewards_(double length) const
RowMatrix< double > insideProduct_
void jFunction_(const std::vector< double > &lambda, double t, RowMatrix< double > &result) const
std::vector< double > lambda_
virtual const Matrix< double > & getRowLeftEigenVectors() const =0
virtual const Matrix< double > & getColumnRightEigenVectors() const =0
Matrix< double > * getAllRewards(double length) const
Get the rewards on a branch, for each initial and final states, and given the branch length...
void setSubstitutionModel(const SubstitutionModel *model)
Set the substitution model.
RowMatrix< double > vInv_