40 #include "../PatternTools.h" 43 #include <Bpp/Seq/SiteTools.h> 44 #include <Bpp/Seq/Container/SequenceContainerTools.h> 46 #include <Bpp/Text/TextTools.h> 47 #include <Bpp/App/ApplicationTools.h> 61 DiscreteDistribution* rDist,
63 bool reparametrizeRoot)
83 reparametrizeRoot_(reparametrizeRoot),
87 init_(tree, modelSet, rDist, verbose);
95 modelSet_(lik.modelSet_),
96 brLenParameters_(lik.brLenParameters_),
100 rootFreqs_(lik.rootFreqs_),
103 nbSites_(lik.nbSites_),
104 nbDistinctSites_(lik.nbDistinctSites_),
105 nbClasses_(lik.nbClasses_),
106 nbStates_(lik.nbStates_),
107 nbNodes_(lik.nbNodes_),
108 verbose_(lik.verbose_),
109 minimumBrLen_(lik.minimumBrLen_),
110 maximumBrLen_(lik.maximumBrLen_),
111 brLenConstraint_(dynamic_cast<Constraint*>(lik.brLenConstraint_->clone())),
112 reparametrizeRoot_(lik.reparametrizeRoot_),
119 for (
unsigned int i = 0; i <
nodes_.size(); i++)
154 for(
unsigned int i = 0; i <
nodes_.size(); i++)
167 DiscreteDistribution* rDist,
168 bool verbose)
throw (Exception)
173 root2_ = tree_->getRootNode()->getSon(1)->getId();
174 nodes_ = tree_->getNodes();
176 nbNodes_ = nodes_.size();
178 for (
unsigned int i = 0; i < nodes_.size(); i++)
180 const Node * node = nodes_[i];
181 idToNode_[node->
getId()] = node;
183 nbClasses_ = rateDistribution_->getNumberOfCategories();
187 minimumBrLen_ = 0.000001;
188 maximumBrLen_ = 10000;
189 brLenConstraint_.reset(
new IntervalConstraint(minimumBrLen_, maximumBrLen_,
true,
true));
190 setSubstitutionModelSet(modelSet);
200 if (modelSet->getAlphabet()->getAlphabetType() != data_->getAlphabet()->getAlphabetType())
201 throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::setSubstitutionModelSet(). Model alphabet do not match existing data.");
204 modelSet_ = modelSet;
208 if (modelSet->getNumberOfStates() != modelSet_->getNumberOfStates())
212 nbStates_ = modelSet->getNumberOfStates();
215 for (
unsigned int l = 0; l < nbNodes_; l++)
218 Node* son = nodes_[l];
220 VVVdouble* pxy__son = & pxy_[son->
getId()];
221 pxy__son->resize(nbClasses_);
222 for (
unsigned int c = 0; c < nbClasses_; c++)
224 VVdouble * pxy__son_c = & (* pxy__son)[c];
225 pxy__son_c->resize(nbStates_);
226 for(
unsigned int x = 0; x < nbStates_; x++)
228 (*pxy__son_c)[x].resize(nbStates_);
232 VVVdouble* dpxy__son = & dpxy_[son->
getId()];
233 dpxy__son->resize(nbClasses_);
234 for (
unsigned int c = 0; c < nbClasses_; c++)
236 VVdouble * dpxy__son_c = & (* dpxy__son)[c];
237 dpxy__son_c->resize(nbStates_);
238 for(
unsigned int x = 0; x < nbStates_; x++)
240 (* dpxy__son_c)[x].resize(nbStates_);
244 VVVdouble* d2pxy__son = & d2pxy_[son->
getId()];
245 d2pxy__son->resize(nbClasses_);
246 for (
unsigned int c = 0; c < nbClasses_; c++)
248 VVdouble * d2pxy__son_c = & (* d2pxy__son)[c];
249 d2pxy__son_c->resize(nbStates_);
250 for(
unsigned int x = 0; x < nbStates_; x++)
252 (* d2pxy__son_c)[x].resize(nbStates_);
261 computeAllTransitionProbabilities();
262 fireParameterChanged(getParameters());
270 if (
initialized_)
throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::initialize(). Object is already initialized.");
271 if (!
data_)
throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::initialize(). Data are no set.");
275 fireParameterChanged(getParameters());
282 if (!
initialized_)
throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
290 if(!
initialized_)
throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::getSubstitutionModelParameters(). Object is not initialized.");
291 return modelSet_->getParameters().getCommonParametersWith(getParameters());
306 addParameters_(
modelSet_->getIndependentParameters());
316 if (!
initialized_)
throw Exception(
"AbstractBranchNonHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
318 for (
unsigned int i = 0; i <
nbNodes_; i++)
320 int id =
nodes_[i]->getId();
323 const Parameter* rootBrLen = &getParameter(
"BrLenRoot");
324 const Parameter* rootPos = &getParameter(
"RootPosition");
325 nodes_[i]->setDistanceToFather(rootBrLen->getValue() * rootPos->getValue());
329 const Parameter* rootBrLen = &getParameter(
"BrLenRoot");
330 const Parameter* rootPos = &getParameter(
"RootPosition");
331 nodes_[i]->setDistanceToFather(rootBrLen->getValue() * (1. - rootPos->getValue()));
335 const Parameter* brLen = &getParameter(
string(
"BrLen") + TextTools::toString(i));
336 if (brLen)
nodes_[i]->setDistanceToFather(brLen->getValue());
341 modelSet_->matchParametersValues(getParameters());
351 double l1 = 0, l2 = 0;
352 for (
unsigned int i = 0; i <
nbNodes_; i++)
355 if (!
nodes_[i]->hasDistanceToFather())
357 ApplicationTools::displayWarning(
"Missing branch length " + TextTools::toString(i) +
". Value is set to " + TextTools::toString(
minimumBrLen_));
362 d =
nodes_[i]->getDistanceToFather();
365 ApplicationTools::displayWarning(
"Branch length " + TextTools::toString(i) +
" is too small: " + TextTools::toString(d) +
". Value is set to " + TextTools::toString(
minimumBrLen_));
371 ApplicationTools::displayWarning(
"Branch length " + TextTools::toString(i) +
" is too big: " + TextTools::toString(d) +
". Value is set to " + TextTools::toString(
maximumBrLen_));
387 brLenParameters_.addParameter(Parameter(
"RootPosition", l1 / (l1 + l2), &Parameter::PROP_CONSTRAINT_EX));
395 for(
unsigned int l = 0; l <
nbNodes_; l++)
412 VVVdouble * pxy__node = &
pxy_[node->
getId()];
415 VVdouble * pxy__node_c = & (* pxy__node)[c];
417 for(
unsigned int x = 0; x <
nbStates_; x++)
419 Vdouble * pxy__node_c_x = & (* pxy__node_c)[x];
420 for(
unsigned int y = 0; y <
nbStates_; y++)
422 (* pxy__node_c_x)[y] = Q(x, y);
430 VVVdouble * dpxy__node = &
dpxy_[node->
getId()];
434 VVdouble * dpxy__node_c = & (* dpxy__node)[c];
437 RowMatrix<double> dQ = model->
getdPij_dt(l * rc);
439 for(
unsigned int x = 0; x <
nbStates_; x++)
441 Vdouble * dpxy__node_c_x = & (* dpxy__node_c)[x];
442 for(
unsigned int y = 0; y <
nbStates_; y++)
443 (* dpxy__node_c_x)[y] = rc * dQ(x, y);
454 VVdouble * d2pxy__node_c = & (* d2pxy__node)[c];
457 for(
unsigned int x = 0; x <
nbStates_; x++)
459 Vdouble * d2pxy__node_c_x = & (* d2pxy__node_c)[x];
460 for(
unsigned int y = 0; y <
nbStates_; y++)
462 (* d2pxy__node_c_x)[y] = rc * rc * d2Q(x, y);
virtual N * getRootNode()
Substitution models manager for non-homogeneous / non-reversible models of evolution.
AbstractNonHomogeneousTreeLikelihood(const Tree &tree, SubstitutionModelSet *modelSet, DiscreteDistribution *rDist, bool verbose=true, bool reparametrizeRoot=true)
ParameterList getSubstitutionModelParameters() const
Get the parameters associated to substitution model(s).
Interface for all substitution models.
Partial implementation for branch non-homogeneous models of the TreeLikelihood interface.
const SubstitutionModel * getModelForNode(int nodeId) const
Get the model associated to a particular node id.
ParameterList brLenParameters_
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
Partial implementation of the DiscreteRatesAcrossSitesTreeLikelihood interface.
SubstitutionModelSet * modelSet_
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
std::vector< double > getRootFrequencies() const
The phylogenetic tree class.
Interface for phylogenetic tree objects.
std::vector< double > rootFreqs_
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.
AbstractNonHomogeneousTreeLikelihood & operator=(const AbstractNonHomogeneousTreeLikelihood &lik)
Assignation operator.
virtual const Matrix< double > & getPij_t(double t) const =0
virtual void initParameters()
This builds the parameters list from all parametrizable objects, i.e. substitution model...
virtual int getId() const
Get the node's id.
std::auto_ptr< Constraint > brLenConstraint_
bool computeFirstOrderDerivatives_
std::map< int, VVVdouble > d2pxy_
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
TreeTemplate< Node > * tree_
virtual void initBranchLengthsParameters()
std::map< int, const Node * > idToNode_
An index linking nodes to their id, for faster access than the getNode() method.
AbstractDiscreteRatesAcrossSitesTreeLikelihood & operator=(const AbstractDiscreteRatesAcrossSitesTreeLikelihood &tl)
The phylogenetic node class.
ParameterList getBranchLengthsParameters() const
Get the branch lengths parameters.
void setSubstitutionModelSet(SubstitutionModelSet *modelSet)
std::map< int, VVVdouble > dpxy_
const SiteContainer * data_
virtual const Matrix< double > & getdPij_dt(double t) const =0
DiscreteDistribution * rateDistribution_
void initialize()
Init the likelihood object.
void init_(const Tree &tree, SubstitutionModelSet *modelSet, DiscreteDistribution *rDist, bool verbose)
Method called by constructor.
bool computeSecondOrderDerivatives_
std::map< int, VVVdouble > pxy_