41 #include "../PatternTools.h" 43 #include <Bpp/Text/TextTools.h> 44 #include <Bpp/App/ApplicationTools.h> 47 #include <Bpp/Seq/SiteTools.h> 48 #include <Bpp/Seq/Container/SequenceContainerTools.h> 62 DiscreteDistribution* rDist,
84 init_(tree, model, rDist, checkRooted, verbose);
93 brLenParameters_(lik.brLenParameters_),
97 rootFreqs_(lik.rootFreqs_),
99 nbSites_(lik.nbSites_),
100 nbDistinctSites_(lik.nbDistinctSites_),
101 nbClasses_(lik.nbClasses_),
102 nbStates_(lik.nbStates_),
103 nbNodes_(lik.nbNodes_),
104 verbose_(lik.verbose_),
105 minimumBrLen_(lik.minimumBrLen_),
106 maximumBrLen_(lik.maximumBrLen_),
107 brLenConstraint_(lik.brLenConstraint_->clone())
145 DiscreteDistribution* rDist,
147 bool verbose)
throw (Exception)
151 if (checkRooted && tree_->isRooted())
154 ApplicationTools::displayWarning(
"Tree has been unrooted.");
157 nodes_ = tree_->getNodes();
159 nbNodes_ = nodes_.size();
160 nbClasses_ = rateDistribution_->getNumberOfCategories();
161 setSubstitutionModel(model);
165 minimumBrLen_ = 0.000001;
166 maximumBrLen_ = 10000;
167 brLenConstraint_.reset(
new IntervalConstraint(minimumBrLen_, maximumBrLen_,
true,
true));
177 if (model->getAlphabet()->getAlphabetType() != data_->getAlphabet()->getAlphabetType())
178 throw Exception(
"AbstractHomogeneousTreeLikelihood::setSubstitutionModel(). Model alphabet do not match existing data.");
185 if (model->getNumberOfStates() != model_->getNumberOfStates())
189 nbStates_ = model->getNumberOfStates();
192 for (
unsigned int l = 0; l < nbNodes_; l++)
195 Node* son = nodes_[l];
197 VVVdouble* pxy__son = &pxy_[son->
getId()];
198 pxy__son->resize(nbClasses_);
199 for (
unsigned int c = 0; c < nbClasses_; c++)
201 VVdouble* pxy__son_c = &(*pxy__son)[c];
202 pxy__son_c->resize(nbStates_);
203 for (
unsigned int x = 0; x < nbStates_; x++)
205 (*pxy__son_c)[x].resize(nbStates_);
209 VVVdouble* dpxy__son = &dpxy_[son->
getId()];
210 dpxy__son->resize(nbClasses_);
211 for (
unsigned int c = 0; c < nbClasses_; c++)
213 VVdouble* dpxy__son_c = &(*dpxy__son)[c];
214 dpxy__son_c->resize(nbStates_);
215 for (
unsigned int x = 0; x < nbStates_; x++)
217 (*dpxy__son_c)[x].resize(nbStates_);
221 VVVdouble* d2pxy__son = &d2pxy_[son->
getId()];
222 d2pxy__son->resize(nbClasses_);
223 for (
unsigned int c = 0; c < nbClasses_; c++)
225 VVdouble* d2pxy__son_c = &(*d2pxy__son)[c];
226 d2pxy__son_c->resize(nbStates_);
227 for (
unsigned int x = 0; x < nbStates_; x++)
229 (*d2pxy__son_c)[x].resize(nbStates_);
240 throw Exception(
"AbstractHomogeneousTreeLikelihood::initialize(). Object is already initialized.");
242 throw Exception(
"AbstractHomogeneousTreeLikelihood::initialize(). Data are no set.");
245 fireParameterChanged(getParameters());
253 throw Exception(
"AbstractHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
262 throw Exception(
"AbstractHomogeneousTreeLikelihood::getSubstitutionModelParameters(). Object is not initialized.");
263 return model_->getParameters().getCommonParametersWith(getParameters());
278 addParameters_(
model_->getIndependentParameters());
289 throw Exception(
"AbstractHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
292 for (
unsigned int i = 0; i <
nbNodes_; i++)
294 const Parameter* brLen = &getParameter(
string(
"BrLen") + TextTools::toString(i));
296 nodes_[i]->setDistanceToFather(brLen->getValue());
299 model_->matchParametersValues(getParameters());
310 for (
unsigned int i = 0; i <
nbNodes_; i++)
313 if (!
nodes_[i]->hasDistanceToFather())
315 ApplicationTools::displayWarning(
"Missing branch length " + TextTools::toString(i) +
". Value is set to " + TextTools::toString(
minimumBrLen_));
320 d =
nodes_[i]->getDistanceToFather();
323 ApplicationTools::displayWarning(
"Branch length " + TextTools::toString(i) +
" is too small: " + TextTools::toString(d) +
". Value is set to " + TextTools::toString(
minimumBrLen_));
329 ApplicationTools::displayWarning(
"Branch length " + TextTools::toString(i) +
" is too big: " + TextTools::toString(d) +
". Value is set to " + TextTools::toString(
maximumBrLen_));
342 for (
unsigned int l = 0; l <
nbNodes_; l++)
358 VVVdouble* pxy__node = &
pxy_[node->
getId()];
361 VVdouble* pxy__node_c = &(*pxy__node)[c];
364 for (
unsigned int x = 0; x <
nbStates_; x++)
366 Vdouble* pxy__node_c_x = &(*pxy__node_c)[x];
367 for (
unsigned int y = 0; y <
nbStates_; y++)
369 (*pxy__node_c_x)[y] = Q(x, y);
377 VVVdouble* dpxy__node = &
dpxy_[node->
getId()];
380 VVdouble* dpxy__node_c = &(*dpxy__node)[c];
383 for (
unsigned int x = 0; x <
nbStates_; x++)
385 Vdouble* dpxy__node_c_x = &(*dpxy__node_c)[x];
386 for (
unsigned int y = 0; y <
nbStates_; y++)
388 (*dpxy__node_c_x)[y] = rc * dQ(x, y);
400 VVdouble* d2pxy__node_c = &(*d2pxy__node)[c];
403 for (
unsigned int x = 0; x <
nbStates_; x++)
405 Vdouble* d2pxy__node_c_x = &(*d2pxy__node_c)[x];
406 for (
unsigned int y = 0; y <
nbStates_; y++)
408 (*d2pxy__node_c_x)[y] = rc * rc * d2Q(x, y);
std::vector< double > rootFreqs_
Interface for all substitution models.
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
Partial implementation of the DiscreteRatesAcrossSitesTreeLikelihood interface.
std::map< int, VVVdouble > dpxy_
AbstractHomogeneousTreeLikelihood & operator=(const AbstractHomogeneousTreeLikelihood &lik)
Assignation operator.
virtual void initParameters()
This builds the parameters list from all parametrizable objects, i.e. substitution model...
The phylogenetic tree class.
Interface for phylogenetic tree objects.
AbstractHomogeneousTreeLikelihood(const Tree &tree, SubstitutionModel *model, DiscreteDistribution *rDist, bool checkRooted=true, bool verbose=true)
SubstitutionModel * model_
void init_(const Tree &tree, SubstitutionModel *model, DiscreteDistribution *rDist, bool checkRooted, bool verbose)
Method called by constructor.
virtual const Vdouble & getFrequencies() const =0
virtual const Matrix< double > & getd2Pij_dt2(double t) const =0
virtual double getDistanceToFather() const
Get the distance to the father node is there is one, otherwise throw a NodeException.
void initialize()
Init the likelihood object.
virtual const Matrix< double > & getPij_t(double t) const =0
virtual int getId() const
Get the node's id.
std::map< int, VVVdouble > d2pxy_
bool computeFirstOrderDerivatives_
virtual void initBranchLengthsParameters()
TreeTemplate< Node > * tree_
ParameterList getSubstitutionModelParameters() const
Get the parameters associated to substitution model(s).
AbstractDiscreteRatesAcrossSitesTreeLikelihood & operator=(const AbstractDiscreteRatesAcrossSitesTreeLikelihood &tl)
The phylogenetic node class.
void setSubstitutionModel(SubstitutionModel *model)
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
std::map< int, VVVdouble > pxy_
const SiteContainer * data_
virtual const Matrix< double > & getdPij_dt(double t) const =0
DiscreteDistribution * rateDistribution_
ParameterList getBranchLengthsParameters() const
Get the branch lengths parameters.
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
Partial implementation for homogeneous model of the TreeLikelihood interface.
bool computeSecondOrderDerivatives_
std::auto_ptr< Constraint > brLenConstraint_
ParameterList brLenParameters_