41 #include "../PatternTools.h" 43 #include <Bpp/Text/TextTools.h> 44 #include <Bpp/App/ApplicationTools.h> 58 DiscreteDistribution* rDist,
74 const SiteContainer& data,
76 DiscreteDistribution* rDist,
95 rateDistribution_->getNumberOfCategories(),
105 minusLogLik_(lik.minusLogLik_)
135 if (data_)
delete data_;
137 if (verbose_) ApplicationTools::displayTask(
"Initializing data structure");
138 likelihoodData_->initLikelihoods(*data_, *model_);
139 if (verbose_) ApplicationTools::displayTaskDone();
141 nbSites_ = likelihoodData_->getNumberOfSites();
142 nbDistinctSites_ = likelihoodData_->getNumberOfDistinctSites();
143 nbStates_ = likelihoodData_->getNumberOfStates();
145 if (verbose_) ApplicationTools::displayResult(
"Number of distinct sites",
146 TextTools::toString(nbDistinctSites_));
147 initialized_ =
false;
155 for (
size_t i = 0; i <
nbSites_; i++)
168 for (
size_t i = 0; i <
nbSites_; i++)
172 sort(la.begin(), la.end());
173 for (
size_t i =
nbSites_; i > 0; i--)
251 throw (ParameterNotFoundException, ConstraintException)
253 setParametersValues(parameters);
262 if (
rateDistribution_->getParameters().getCommonParametersWith(params).size() > 0
263 ||
model_->getParameters().getCommonParametersWith(params).size() > 0)
268 else if (params.size() > 0)
271 for (
size_t i = 0; i < params.size(); i++)
273 string s = params[i].getName();
274 if (s.substr(0, 5) ==
"BrLen")
293 if (!
isInitialized())
throw Exception(
"RHomogeneousTreeLikelihood::getValue(). Instance is not initialized.");
303 size_t rateClass)
const 341 for (
size_t i = 0; i <
nbSites_; i++)
353 if (!hasParameter(variable))
354 throw ParameterNotFoundException(
"RHomogeneousTreeLikelihood::getFirstOrderDerivative().", variable);
355 if (getRateDistributionParameters().hasParameter(variable))
357 throw Exception(
"Derivatives respective to rate distribution parameter are not implemented.");
359 if (getSubstitutionModelParameters().hasParameter(variable))
361 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
365 return -getDLogLikelihood();
373 size_t brI = TextTools::to<size_t>(variable.substr(5));
380 size_t nbSites = _dLikelihoods_father->size();
381 for (
size_t i = 0; i < nbSites; i++)
383 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
386 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
389 (*_dLikelihoods_father_i_c)[s] = 1.;
395 for (
size_t l = 0; l < nbNodes; l++)
405 for (
size_t i = 0; i < nbSites; i++)
407 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
408 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
411 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
412 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
413 VVdouble* dpxy__son_c = &(*dpxy__son)[c];
417 Vdouble* dpxy__son_c_x = &(*dpxy__son_c)[x];
420 dl += (*dpxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
422 (*_dLikelihoods_father_i_c)[x] *= dl;
429 VVVdouble* pxy__son = &
pxy_[son->
getId()];
430 for (
size_t i = 0; i < nbSites; i++)
432 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
433 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
436 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
437 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
438 VVdouble* pxy__son_c = &(*pxy__son)[c];
442 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
445 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
447 (*_dLikelihoods_father_i_c)[x] *= dl;
465 if (father == NULL)
return;
470 size_t nbSites = _dLikelihoods_father->size();
471 for (
size_t i = 0; i < nbSites; i++)
473 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
476 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
479 (*_dLikelihoods_father_i_c)[s] = 1.;
485 for (
size_t l = 0; l < nbNodes; l++)
489 VVVdouble* pxy__son = &
pxy_[son->
getId()];
495 for (
size_t i = 0; i < nbSites; i++)
497 VVdouble* _dLikelihoods_son_i = &(*_dLikelihoods_son)[(*_patternLinks_father_son)[i]];
498 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
501 Vdouble* _dLikelihoods_son_i_c = &(*_dLikelihoods_son_i)[c];
502 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
503 VVdouble* pxy__son_c = &(*pxy__son)[c];
507 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
510 dl += (*pxy__son_c_x)[y] * (*_dLikelihoods_son_i_c)[y];
512 (*_dLikelihoods_father_i_c)[x] *= dl;
520 for (
size_t i = 0; i < nbSites; i++)
522 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
523 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
526 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
527 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
528 VVdouble* pxy__son_c = &(*pxy__son)[c];
532 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
535 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
537 (*_dLikelihoods_father_i_c)[x] *= dl;
554 size_t rateClass)
const 592 for (
size_t i = 0; i <
nbSites_; i++)
604 if (!hasParameter(variable))
605 throw ParameterNotFoundException(
"RHomogeneousTreeLikelihood::getSecondOrderDerivative().", variable);
606 if (getRateDistributionParameters().hasParameter(variable))
608 throw Exception(
"Derivatives respective to rate distribution parameter are not implemented.");
610 if (getSubstitutionModelParameters().hasParameter(variable))
612 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
616 return -getD2LogLikelihood();
624 size_t brI = TextTools::to<size_t>(variable.substr(5));
631 size_t nbSites = _d2Likelihoods_father->size();
632 for (
size_t i = 0; i < nbSites; i++)
634 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
637 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
640 (*_d2Likelihoods_father_i_c)[s] = 1.;
646 for (
size_t l = 0; l < nbNodes; l++)
656 for (
size_t i = 0; i < nbSites; i++)
658 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
659 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
662 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
663 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
664 VVdouble* d2pxy__son_c = &(*d2pxy__son)[c];
668 Vdouble* d2pxy__son_c_x = &(*d2pxy__son_c)[x];
671 d2l += (*d2pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
673 (*_d2Likelihoods_father_i_c)[x] *= d2l;
680 VVVdouble* pxy__son = &
pxy_[son->
getId()];
681 for (
size_t i = 0; i < nbSites; i++)
683 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
684 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
687 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
688 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
689 VVdouble* pxy__son_c = &(*pxy__son)[c];
693 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
696 d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
698 (*_d2Likelihoods_father_i_c)[x] *= d2l;
716 if (father == NULL)
return;
721 size_t nbSites = _d2Likelihoods_father->size();
722 for (
size_t i = 0; i < nbSites; i++)
724 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
727 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
730 (*_d2Likelihoods_father_i_c)[s] = 1.;
736 for (
size_t l = 0; l < nbNodes; l++)
740 VVVdouble* pxy__son = &
pxy_[son->
getId()];
746 for (
size_t i = 0; i < nbSites; i++)
748 VVdouble* _d2Likelihoods_son_i = &(*_d2Likelihoods_son)[(*_patternLinks_father_son)[i]];
749 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
752 Vdouble* _d2Likelihoods_son_i_c = &(*_d2Likelihoods_son_i)[c];
753 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
754 VVdouble* pxy__son_c = &(*pxy__son)[c];
758 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
761 d2l += (*pxy__son_c_x)[y] * (*_d2Likelihoods_son_i_c)[y];
763 (*_d2Likelihoods_father_i_c)[x] *= d2l;
771 for (
size_t i = 0; i < nbSites; i++)
773 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
774 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
777 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
778 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
779 VVdouble* pxy__son_c = &(*pxy__son)[c];
783 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
786 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
788 (*_d2Likelihoods_father_i_c)[x] *= dl;
810 if (node->
isLeaf())
return;
817 for (
size_t i = 0; i < nbSites; i++)
820 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
824 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
828 (*_likelihoods_node_i_c)[x] = 1.;
833 for (
size_t l = 0; l < nbNodes; l++)
841 VVVdouble* pxy__son = &
pxy_[son->
getId()];
845 for (
size_t i = 0; i < nbSites; i++)
848 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_node_son)[i]];
849 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
853 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
854 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
855 VVdouble* pxy__son_c = &(*pxy__son)[c];
859 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
860 double likelihood = 0;
863 likelihood += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
865 (*_likelihoods_node_i_c)[x] *= likelihood;
876 cout <<
"Likelihoods at node " << node->
getName() <<
": " << endl;
878 cout <<
" ***" << endl;
virtual void computeDownSubtreeDLikelihood(const Node *)
std::vector< double > rootFreqs_
Interface for all substitution models.
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
RHomogeneousTreeLikelihood(const Tree &tree, SubstitutionModel *model, DiscreteDistribution *rDist, bool checkRooted=true, bool verbose=true, bool usePatterns=true)
Build a new RHomogeneousTreeLikelihood object without data.
double getLogLikelihoodForASite(size_t site) const
Get the logarithm of the likelihood for a site.
virtual std::string getName() const
Get the name associated to this node, if there is one, otherwise throw a NodeException.
double getLogLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
Get the logarithm of the likelihood for a site knowing its rate class and its ancestral state...
size_t getRootArrayPosition(size_t currentPosition) const
std::map< int, VVVdouble > dpxy_
DRASRTreeLikelihoodData * likelihoodData_
AbstractHomogeneousTreeLikelihood & operator=(const AbstractHomogeneousTreeLikelihood &lik)
Assignation operator.
void setTree(const TreeTemplate< Node > *tree)
Set the tree associated to the data.
virtual double getD2LogLikelihoodForASite(size_t site) const
virtual void computeDownSubtreeD2Likelihood(const Node *)
virtual double getD2LikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
virtual void computeTreeD2Likelihood(const std::string &variable)
virtual const Node * getSon(size_t pos) const
discrete Rate Across Sites, (simple) Recursive likelihood data structure.
double getLikelihood() const
Get the likelihood for the whole dataset.
virtual double getD2LogLikelihood() const
virtual const Node * getFather() const
Get the father of this node is there is one.
Interface for phylogenetic tree objects.
virtual ~RHomogeneousTreeLikelihood()
static void displayLikelihoodArray(const VVVdouble &likelihoodArray)
Print the likelihood array to terminal (debugging tool).
This class implement the 'traditional' way of computing likelihood for a tree.
SubstitutionModel * model_
virtual const Vdouble & getFrequencies() const =0
virtual bool isLeaf() const
virtual double getDLogLikelihoodForASite(size_t site) const
RHomogeneousTreeLikelihood & operator=(const RHomogeneousTreeLikelihood &lik)
VVVdouble & getDLikelihoodArray(int nodeId)
void fireParameterChanged(const ParameterList ¶ms)
virtual int getId() const
Get the node's id.
double getLogLikelihood() const
Get the logarithm of the likelihood for the whole dataset.
std::map< int, VVVdouble > d2pxy_
TreeTemplate< Node > * tree_
virtual double getD2LikelihoodForASite(size_t site) const
The phylogenetic node class.
double getLikelihoodForASite(size_t site) const
Get the likelihood for a site.
void setParameters(const ParameterList ¶meters)
Implements the Function interface.
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
double getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the logarithm of the likelihood for a site knowing its rate class.
DRASRTreeLikelihoodData * clone() const
virtual void computeTreeDLikelihood(const std::string &variable)
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
virtual void computeSubtreeLikelihood(const Node *node)
Compute the likelihood for a subtree defined by the Tree::Node node.
std::map< int, VVVdouble > pxy_
virtual double getDLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
DiscreteDistribution * rateDistribution_
double getSecondOrderDerivative(const std::string &variable) const
virtual size_t getNumberOfSons() const
virtual double getDLikelihoodForASite(size_t site) const
double getFirstOrderDerivative(const std::string &variable) const
virtual void displayLikelihood(const Node *node)
This method is mainly for debugging purpose.
const std::vector< size_t > & getArrayPositions(int parentId, int sonId) const
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
double getLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
Get the likelihood for a site knowing its rate class and its ancestral state.
void computeTreeLikelihood()
VVVdouble & getLikelihoodArray(int nodeId)
Partial implementation for homogeneous model of the TreeLikelihood interface.
bool isInitialized() const
VVVdouble & getD2LikelihoodArray(int nodeId)
void init_(bool usePatterns)
Method called by constructors.
double getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the likelihood for a site knowing its rate class.
void setData(const SiteContainer &sites)
Set the dataset for which the likelihood must be evaluated.
virtual double getDLogLikelihood() const