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/AlignedSequenceContainer.h> 49 #include <Bpp/Seq/Container/SequenceContainerTools.h> 63 DiscreteDistribution* rDist,
65 bool reparametrizeRoot)
71 if (!modelSet->isFullySetUpFor(tree))
72 throw Exception(
"DRNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
80 const SiteContainer& data,
82 DiscreteDistribution* rDist,
84 bool reparametrizeRoot)
90 if (!modelSet->isFullySetUpFor(tree))
91 throw Exception(
"DRNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
102 rateDistribution_->getNumberOfCategories());
110 minusLogLik_(lik.minusLogLik_)
144 ApplicationTools::displayTask(
"Initializing data structure");
145 likelihoodData_->initLikelihoods(*data_, *modelSet_->getModel(0));
148 ApplicationTools::displayTaskDone();
150 nbSites_ = likelihoodData_->getNumberOfSites();
151 nbDistinctSites_ = likelihoodData_->getNumberOfDistinctSites();
152 nbStates_ = likelihoodData_->getNumberOfStates();
155 ApplicationTools::displayResult(
"Number of distinct sites",
156 TextTools::toString(nbDistinctSites_));
157 initialized_ =
false;
169 l *= std::pow((*lik)[i], (
int)(*w)[i]);
184 la[i] = (*w)[i] * log((*lik)[i]);
186 sort(la.begin(), la.end());
238 throw (ParameterNotFoundException, ConstraintException)
240 setParametersValues(parameters);
249 if (params.getCommonParametersWith(
rateDistribution_->getIndependentParameters()).size() > 0)
257 for (
size_t i = 0; i < tmp.size(); i++)
260 ids = VectorTools::vectorUnion(ids, tmpv);
263 vector<const Node*> nodes;
264 for (
size_t i = 0; i < ids.size(); i++)
268 vector<const Node*> tmpv;
270 for (
size_t i = 0; i < tmp.size(); i++)
272 if (tmp[i] ==
"BrLenRoot" || tmp[i] ==
"RootPosition")
276 tmpv.push_back(
tree_->getRootNode()->getSon(0));
277 tmpv.push_back(
tree_->getRootNode()->getSon(1));
282 tmpv.push_back(
nodes_[TextTools::to < size_t > (tmp[i].substr(5))]);
284 nodes = VectorTools::vectorUnion(nodes, tmpv);
286 for (
size_t i = 0; i < nodes.size(); i++)
309 throw Exception(
"DRNonHomogeneousTreeLikelihood::getValue(). Instance is not initialized.");
321 VVVdouble* pxy__node = &
pxy_[node->
getId()];
322 VVVdouble* dpxy__node = &
dpxy_[node->
getId()];
327 double dLi, dLic, dLicx, numerator, denominator;
330 VVdouble* _likelihoods_father_node_i = &(*_likelihoods_father_node)[i];
331 VVdouble* larray_i = &larray[i];
335 Vdouble* _likelihoods_father_node_i_c = &(*_likelihoods_father_node_i)[c];
336 Vdouble* larray_i_c = &(*larray_i)[c];
337 VVdouble* pxy__node_c = &(*pxy__node)[c];
338 VVdouble* dpxy__node_c = &(*dpxy__node)[c];
344 Vdouble* pxy__node_c_x = &(*pxy__node_c)[x];
345 Vdouble* dpxy__node_c_x = &(*dpxy__node_c)[x];
349 numerator += (*dpxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
350 denominator += (*pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
352 dLicx = denominator == 0. ? 0. : (*larray_i_c)[x] * numerator / denominator;
357 (*_dLikelihoods_node)[i] = dLi / (*rootLikelihoodsSR)[i];
365 for (
size_t k = 0; k <
nbNodes_; k++)
376 if (!hasParameter(variable))
377 throw ParameterNotFoundException(
"DRNonHomogeneousTreeLikelihood::getFirstOrderDerivative().", variable);
378 if (getRateDistributionParameters().hasParameter(variable))
380 throw Exception(
"Derivatives respective to rate distribution parameters are not implemented.");
382 if (getSubstitutionModelParameters().hasParameter(variable))
384 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
390 const vector<unsigned int>* w = &likelihoodData_->getWeights();
391 Vdouble* _dLikelihoods_branch;
392 if (variable ==
"BrLenRoot")
394 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(root1_);
396 for (
size_t i = 0; i < nbDistinctSites_; i++)
398 d1 -= (*w)[i] * (*_dLikelihoods_branch)[i];
400 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(root2_);
402 for (
size_t i = 0; i < nbDistinctSites_; i++)
404 d2 -= (*w)[i] * (*_dLikelihoods_branch)[i];
406 double pos = getParameterValue(
"RootPosition");
407 return pos * d1 + (1. - pos) * d2;
409 else if (variable ==
"RootPosition")
411 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(root1_);
413 for (
size_t i = 0; i < nbDistinctSites_; i++)
415 d1 -= (*w)[i] * (*_dLikelihoods_branch)[i];
417 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(root2_);
419 for (
size_t i = 0; i < nbDistinctSites_; i++)
421 d2 -= (*w)[i] * (*_dLikelihoods_branch)[i];
423 double len = getParameterValue(
"BrLenRoot");
424 return len * (d1 - d2);
429 size_t brI = TextTools::to<size_t>(variable.substr(5));
430 const Node* branch = nodes_[brI];
431 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(branch->getId());
433 for (
size_t i = 0; i < nbDistinctSites_; i++)
435 d += (*w)[i] * (*_dLikelihoods_branch)[i];
449 VVVdouble* pxy__node = &
pxy_[node->
getId()];
455 double d2Li, d2Lic, d2Licx, numerator, denominator;
459 VVdouble* _likelihoods_father_node_i = &(*_likelihoods_father_node)[i];
460 VVdouble* larray_i = &larray[i];
464 Vdouble* _likelihoods_father_node_i_c = &(*_likelihoods_father_node_i)[c];
465 Vdouble* larray_i_c = &(*larray_i)[c];
466 VVdouble* pxy__node_c = &(*pxy__node)[c];
467 VVdouble* d2pxy__node_c = &(*d2pxy__node)[c];
473 Vdouble* pxy__node_c_x = &(*pxy__node_c)[x];
474 Vdouble* d2pxy__node_c_x = &(*d2pxy__node_c)[x];
478 numerator += (*d2pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
479 denominator += (*pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
481 d2Licx = denominator == 0. ? 0. : (*larray_i_c)[x] * numerator / denominator;
486 (*_d2Likelihoods_node)[i] = d2Li / (*rootLikelihoodsSR)[i];
494 for (
size_t k = 0; k <
nbNodes_; k++)
505 if (!hasParameter(variable))
506 throw ParameterNotFoundException(
"DRNonHomogeneousTreeLikelihood::getSecondOrderDerivative().", variable);
507 if (getRateDistributionParameters().hasParameter(variable))
509 throw Exception(
"Derivatives respective to rate distribution parameters are not implemented.");
511 if (getSubstitutionModelParameters().hasParameter(variable))
513 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
520 const vector<unsigned int>* w = &likelihoodData_->getWeights();
524 if (variable ==
"BrLenRoot")
526 const Node* father = tree_->getRootNode();
530 VVVdouble dLikelihoods_father(nbDistinctSites_);
531 VVVdouble d2Likelihoods_father(nbDistinctSites_);
532 for (
size_t i = 0; i < nbDistinctSites_; i++)
534 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
535 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
536 dLikelihoods_father_i->resize(nbClasses_);
537 d2Likelihoods_father_i->resize(nbClasses_);
538 for (
size_t c = 0; c < nbClasses_; c++)
540 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
541 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
542 dLikelihoods_father_i_c->resize(nbStates_);
543 d2Likelihoods_father_i_c->resize(nbStates_);
544 for (
size_t s = 0; s < nbStates_; s++)
546 (*dLikelihoods_father_i_c)[s] = 1.;
547 (*d2Likelihoods_father_i_c)[s] = 1.;
553 for (
size_t l = 0; l < nbNodes; l++)
557 if (son->
getId() == root1_)
559 VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(father->
getId(), root1_);
560 VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(father->
getId(), root2_);
561 double pos = getParameterValue(
"RootPosition");
563 VVVdouble* d2pxy_root1_ = &d2pxy_[root1_];
564 VVVdouble* d2pxy_root2_ = &d2pxy_[root2_];
565 VVVdouble* dpxy_root1_ = &dpxy_[root1_];
566 VVVdouble* dpxy_root2_ = &dpxy_[root2_];
567 VVVdouble* pxy_root1_ = &pxy_[root1_];
568 VVVdouble* pxy_root2_ = &pxy_[root2_];
569 for (
size_t i = 0; i < nbDistinctSites_; i++)
571 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[i];
572 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[i];
573 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
574 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
575 for (
size_t c = 0; c < nbClasses_; c++)
577 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
578 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
579 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
580 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
581 VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
582 VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
583 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
584 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
585 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
586 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
587 for (
size_t x = 0; x < nbStates_; x++)
589 Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
590 Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
591 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
592 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
593 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
594 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
595 double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
596 for (
size_t y = 0; y < nbStates_; y++)
598 d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
599 d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
600 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
601 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
602 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
603 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
605 double dl = pos * dl1 * l2 + (1. - pos) * dl2 * l1;
606 double d2l = pos * pos * d2l1 * l2 + (1. - pos) * (1. - pos) * d2l2 * l1 + 2 * pos * (1. - pos) * dl1 * dl2;
607 (*dLikelihoods_father_i_c)[x] *= dl;
608 (*d2Likelihoods_father_i_c)[x] *= d2l;
613 else if (son->
getId() == root2_)
620 VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(father->
getId(), son->
getId());
622 VVVdouble* pxy__son = &pxy_[son->
getId()];
623 for (
size_t i = 0; i < nbDistinctSites_; i++)
625 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[i];
626 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
627 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
628 for (
size_t c = 0; c < nbClasses_; c++)
630 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
631 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
632 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
633 VVdouble* pxy__son_c = &(*pxy__son)[c];
634 for (
size_t x = 0; x < nbStates_; x++)
637 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
638 for (
size_t y = 0; y < nbStates_; y++)
640 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
642 (*dLikelihoods_father_i_c)[x] *= dl;
643 (*d2Likelihoods_father_i_c)[x] *= dl;
649 Vdouble* rootLikelihoodsSR = &likelihoodData_->getRootRateSiteLikelihoodArray();
650 double d2l = 0, dlx, d2lx;
651 for (
size_t i = 0; i < nbDistinctSites_; i++)
653 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
654 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
656 for (
size_t c = 0; c < nbClasses_; c++)
658 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
659 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
660 for (
size_t x = 0; x < nbStates_; x++)
662 dlx += rateDistribution_->getProbability(c) * rootFreqs_[x] * (*dLikelihoods_father_i_c)[x];
663 d2lx += rateDistribution_->getProbability(c) * rootFreqs_[x] * (*d2Likelihoods_father_i_c)[x];
666 d2l += (*w)[i] * (d2lx / (*rootLikelihoodsSR)[i] - pow(dlx / (*rootLikelihoodsSR)[i], 2));
670 else if (variable ==
"RootPosition")
672 const Node* father = tree_->getRootNode();
676 VVVdouble dLikelihoods_father(nbDistinctSites_);
677 VVVdouble d2Likelihoods_father(nbDistinctSites_);
678 for (
size_t i = 0; i < nbDistinctSites_; i++)
680 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
681 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
682 dLikelihoods_father_i->resize(nbClasses_);
683 d2Likelihoods_father_i->resize(nbClasses_);
684 for (
size_t c = 0; c < nbClasses_; c++)
686 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
687 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
688 dLikelihoods_father_i_c->resize(nbStates_);
689 d2Likelihoods_father_i_c->resize(nbStates_);
690 for (
size_t s = 0; s < nbStates_; s++)
692 (*dLikelihoods_father_i_c)[s] = 1.;
693 (*d2Likelihoods_father_i_c)[s] = 1.;
699 for (
size_t l = 0; l < nbNodes; l++)
703 if (son->
getId() == root1_)
705 VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(father->
getId(), root1_);
706 VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(father->
getId(), root2_);
707 double len = getParameterValue(
"BrLenRoot");
709 VVVdouble* d2pxy_root1_ = &d2pxy_[root1_];
710 VVVdouble* d2pxy_root2_ = &d2pxy_[root2_];
711 VVVdouble* dpxy_root1_ = &dpxy_[root1_];
712 VVVdouble* dpxy_root2_ = &dpxy_[root2_];
713 VVVdouble* pxy_root1_ = &pxy_[root1_];
714 VVVdouble* pxy_root2_ = &pxy_[root2_];
715 for (
size_t i = 0; i < nbDistinctSites_; i++)
717 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[i];
718 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[i];
719 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
720 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
721 for (
size_t c = 0; c < nbClasses_; c++)
723 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
724 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
725 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
726 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
727 VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
728 VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
729 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
730 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
731 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
732 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
733 for (
size_t x = 0; x < nbStates_; x++)
735 Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
736 Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
737 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
738 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
739 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
740 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
741 double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
742 for (
size_t y = 0; y < nbStates_; y++)
744 d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
745 d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
746 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
747 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
748 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
749 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
751 double dl = len * (dl1 * l2 - dl2 * l1);
752 double d2l = len * len * (d2l1 * l2 + d2l2 * l1 - 2 * dl1 * dl2);
753 (*dLikelihoods_father_i_c)[x] *= dl;
754 (*d2Likelihoods_father_i_c)[x] *= d2l;
759 else if (son->
getId() == root2_)
766 VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(father->
getId(), son->
getId());
768 VVVdouble* pxy__son = &pxy_[son->
getId()];
769 for (
size_t i = 0; i < nbDistinctSites_; i++)
771 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[i];
772 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
773 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
774 for (
size_t c = 0; c < nbClasses_; c++)
776 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
777 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
778 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
779 VVdouble* pxy__son_c = &(*pxy__son)[c];
780 for (
size_t x = 0; x < nbStates_; x++)
783 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
784 for (
size_t y = 0; y < nbStates_; y++)
786 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
788 (*dLikelihoods_father_i_c)[x] *= dl;
789 (*d2Likelihoods_father_i_c)[x] *= dl;
795 Vdouble* rootLikelihoodsSR = &likelihoodData_->getRootRateSiteLikelihoodArray();
796 double d2l = 0, dlx, d2lx;
797 for (
size_t i = 0; i < nbDistinctSites_; i++)
799 VVdouble* dLikelihoods_father_i = &dLikelihoods_father[i];
800 VVdouble* d2Likelihoods_father_i = &d2Likelihoods_father[i];
802 for (
size_t c = 0; c < nbClasses_; c++)
804 Vdouble* dLikelihoods_father_i_c = &(*dLikelihoods_father_i)[c];
805 Vdouble* d2Likelihoods_father_i_c = &(*d2Likelihoods_father_i)[c];
806 for (
size_t x = 0; x < nbStates_; x++)
808 dlx += rateDistribution_->getProbability(c) * rootFreqs_[x] * (*dLikelihoods_father_i_c)[x];
809 d2lx += rateDistribution_->getProbability(c) * rootFreqs_[x] * (*d2Likelihoods_father_i_c)[x];
812 d2l += (*w)[i] * (d2lx / (*rootLikelihoodsSR)[i] - pow(dlx / (*rootLikelihoodsSR)[i], 2));
818 Vdouble* _dLikelihoods_branch;
819 Vdouble* _d2Likelihoods_branch;
821 size_t brI = TextTools::to<size_t>(variable.substr(5));
822 const Node* branch = nodes_[brI];
823 _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(branch->getId());
824 _d2Likelihoods_branch = &likelihoodData_->getD2LikelihoodArray(branch->getId());
826 for (
size_t i = 0; i < nbDistinctSites_; i++)
828 d2l += (*w)[i] * ((*_d2Likelihoods_branch)[i] - pow((*_dLikelihoods_branch)[i], 2));
873 for (
size_t l = 0; l < nbNodes; l++)
878 VVVdouble* _likelihoods_node_son = &(*_likelihoods_node)[son->
getId()];
886 Vdouble* _likelihoods_leaf_i = &(*_likelihoods_leaf)[i];
887 VVdouble* _likelihoods_node_son_i = &(*_likelihoods_node_son)[i];
891 Vdouble* _likelihoods_node_son_i_c = &(*_likelihoods_node_son_i)[c];
895 (*_likelihoods_node_son_i_c)[x] = (*_likelihoods_leaf_i)[x];
906 vector<const VVVdouble*> iLik(nbSons);
907 vector<const VVVdouble*> tProb(nbSons);
908 for (
size_t n = 0; n < nbSons; n++)
912 iLik[n] = &(*_likelihoods_son)[sonSon->
getId()];
928 for (
size_t n = 0; n < nbSons; n++)
939 VVVdouble* _likelihoods_node_father = &(*_likelihoods_node)[father->
getId()];
952 Vdouble* _likelihoods_leaf_i = &(*_likelihoods_leaf)[i];
953 VVdouble* _likelihoods_node_father_i = &(*_likelihoods_node_father)[i];
957 Vdouble* _likelihoods_node_father_i_c = &(*_likelihoods_node_father_i)[c];
961 (*_likelihoods_node_father_i_c)[x] = (*_likelihoods_leaf_i)[x];
968 vector<const Node*> nodes;
971 for (
size_t n = 0; n < nbFatherSons; n++)
975 nodes.push_back(son);
981 size_t nbSons = nodes.size();
983 vector<const VVVdouble*> iLik(nbSons);
984 vector<const VVVdouble*> tProb(nbSons);
985 for (
size_t n = 0; n < nbSons; n++)
987 const Node* fatherSon = nodes[n];
989 iLik[n] = &(*_likelihoods_father)[fatherSon->
getId()];
1008 VVdouble* _likelihoods_node_father_i = &(*_likelihoods_node_father)[i];
1011 Vdouble* _likelihoods_node_father_i_c = &(*_likelihoods_node_father_i)[c];
1014 (*_likelihoods_node_father_i_c)[x] *=
rootFreqs_[x];
1022 for (
size_t i = 0; i < nbNodeSons; i++)
1033 const Node* root =
tree_->getRootNode();
1041 VVdouble* rootLikelihoods_i = &(*rootLikelihoods)[i];
1042 Vdouble* leavesLikelihoods_root_i = &(*leavesLikelihoods_root)[i];
1045 Vdouble* rootLikelihoods_i_c = &(*rootLikelihoods_i)[c];
1048 (*rootLikelihoods_i_c)[x] = (*leavesLikelihoods_root_i)[x];
1060 vector<const VVVdouble*> iLik(nbNodes);
1061 vector<const VVVdouble*> tProb(nbNodes);
1062 for (
size_t n = 0; n < nbNodes; n++)
1066 iLik[n] = &(*likelihoods_root)[son->
getId()];
1076 VVdouble* rootLikelihoods_i = &(*rootLikelihoods)[i];
1077 Vdouble* rootLikelihoodsS_i = &(*rootLikelihoodsS)[i];
1078 (*rootLikelihoodsSR)[i] = 0;
1082 Vdouble* rootLikelihoods_i_c = &(*rootLikelihoods_i)[c];
1083 double* rootLikelihoodsS_i_c = &(*rootLikelihoodsS_i)[c];
1084 (*rootLikelihoodsS_i_c) = 0;
1088 (*rootLikelihoodsS_i_c) +=
rootFreqs_[x] * (*rootLikelihoods_i_c)[x];
1090 (*rootLikelihoodsSR)[i] += p[c] * (*rootLikelihoodsS_i_c);
1094 if ((*rootLikelihoodsSR)[i] < 0)
1095 (*rootLikelihoodsSR)[i] = 0.;
1104 int nodeId = node->
getId();
1114 VVdouble* likelihoodArray_i = &likelihoodArray[i];
1115 Vdouble* leavesLikelihoods_node_i = &(*leavesLikelihoods_node)[i];
1119 Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
1123 (*likelihoodArray_i_c)[x] = (*leavesLikelihoods_node_i)[x];
1134 VVdouble* likelihoodArray_i = &likelihoodArray[i];
1138 Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
1142 (*likelihoodArray_i_c)[x] = 1.;
1150 vector<const VVVdouble*> iLik(nbNodes);
1151 vector<const VVVdouble*> tProb(nbNodes);
1152 for (
size_t n = 0; n < nbNodes; n++)
1156 iLik[n] = &(*likelihoods_node)[son->
getId()];
1171 VVdouble* likelihoodArray_i = &likelihoodArray[i];
1174 Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
1187 const vector<const VVVdouble*>& iLik,
1188 const vector<const VVVdouble*>& tProb,
1191 size_t nbDistinctSites,
1199 for (
size_t n = 0; n < nbNodes; n++)
1201 const VVVdouble* pxy_n = tProb[n];
1202 const VVVdouble* iLik_n = iLik[n];
1204 for (
size_t i = 0; i < nbDistinctSites; i++)
1207 const VVdouble* iLik_n_i = &(*iLik_n)[i];
1208 VVdouble* oLik_i = &(oLik)[i];
1210 for (
size_t c = 0; c < nbClasses; c++)
1213 const Vdouble* iLik_n_i_c = &(*iLik_n_i)[c];
1214 Vdouble* oLik_i_c = &(*oLik_i)[c];
1215 const VVdouble* pxy_n_c = &(*pxy_n)[c];
1216 for (
size_t x = 0; x < nbStates; x++)
1219 const Vdouble* pxy_n_c_x = &(*pxy_n_c)[x];
1220 double likelihood = 0;
1221 for (
size_t y = 0; y < nbStates; y++)
1223 likelihood += (*pxy_n_c_x)[y] * (*iLik_n_i_c)[y];
1226 (*oLik_i_c)[x] *= likelihood;
1236 const vector<const VVVdouble*>& iLik,
1237 const vector<const VVVdouble*>& tProb,
1238 const VVVdouble* iLikR,
1239 const VVVdouble* tProbR,
1242 size_t nbDistinctSites,
1250 for (
size_t n = 0; n < nbNodes; n++)
1252 const VVVdouble* pxy_n = tProb[n];
1253 const VVVdouble* iLik_n = iLik[n];
1255 for (
size_t i = 0; i < nbDistinctSites; i++)
1258 const VVdouble* iLik_n_i = &(*iLik_n)[i];
1259 VVdouble* oLik_i = &(oLik)[i];
1261 for (
size_t c = 0; c < nbClasses; c++)
1264 const Vdouble* iLik_n_i_c = &(*iLik_n_i)[c];
1265 Vdouble* oLik_i_c = &(*oLik_i)[c];
1266 const VVdouble* pxy_n_c = &(*pxy_n)[c];
1267 for (
size_t x = 0; x < nbStates; x++)
1270 const Vdouble* pxy_n_c_x = &(*pxy_n_c)[x];
1271 double likelihood = 0;
1272 for (
size_t y = 0; y < nbStates; y++)
1276 likelihood += (*pxy_n_c_x)[y] * (*iLik_n_i_c)[y];
1280 (*oLik_i_c)[x] *= likelihood;
1287 for (
size_t i = 0; i < nbDistinctSites; i++)
1290 const VVdouble* iLikR_i = &(*iLikR)[i];
1291 VVdouble* oLik_i = &(oLik)[i];
1293 for (
size_t c = 0; c < nbClasses; c++)
1296 const Vdouble* iLikR_i_c = &(*iLikR_i)[c];
1297 Vdouble* oLik_i_c = &(*oLik_i)[c];
1298 const VVdouble* pxyR_c = &(*tProbR)[c];
1299 for (
size_t x = 0; x < nbStates; x++)
1301 double likelihood = 0;
1302 for (
size_t y = 0; y < nbStates; y++)
1305 const Vdouble* pxyR_c_y = &(*pxyR_c)[y];
1306 likelihood += (*pxyR_c_y)[x] * (*iLikR_i_c)[y];
1309 (*oLik_i_c)[x] *= likelihood;
1319 cout <<
"Likelihoods at node " << node->
getId() <<
": " << endl;
1323 cout <<
"Array for sub-node " << subNode->
getId() << endl;
1329 cout <<
"Array for father node " << father->
getId() << endl;
1332 cout <<
" ***" << endl;
void setData(const SiteContainer &sites)
Set the dataset for which the likelihood must be evaluated.
Substitution models manager for non-homogeneous / non-reversible models of evolution.
double getLogLikelihoodForASite(size_t site) const
Get the logarithm of the likelihood for a site.
Partial implementation for branch non-homogeneous models of the TreeLikelihood interface.
DRNonHomogeneousTreeLikelihood(const Tree &tree, SubstitutionModelSet *modelSet, DiscreteDistribution *rDist, bool verbose=true, bool reparametrizeRoot=false)
Build a new DRNonHomogeneousTreeLikelihood object without data.
ParameterList brLenParameters_
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
VVVdouble & getRootLikelihoodArray()
double getLogLikelihood() const
Get the logarithm of the likelihood for the whole dataset.
VVdouble & getRootSiteLikelihoodArray()
virtual void computeTreeD2Likelihoods()
SubstitutionModelSet * modelSet_
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
virtual void computeTreeDLikelihoodAtNode(const Node *node)
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
virtual const Node * getSon(size_t pos) const
void resetLikelihoodArrays(const Node *node)
std::vector< double > getRootFrequencies() const
Vdouble & getRootRateSiteLikelihoodArray()
double getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the likelihood for a site knowing its rate class.
const std::map< int, VVVdouble > & getLikelihoodArrays(int nodeId) const
double getFirstOrderDerivative(const std::string &variable) const
Vdouble & getDLikelihoodArray(int nodeId)
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.
const std::vector< unsigned int > & getWeights() const
virtual const Node * getFather() const
Get the father of this node is there is one.
Vdouble & getD2LikelihoodArray(int nodeId)
Interface for phylogenetic tree objects.
static void displayLikelihoodArray(const VVVdouble &likelihoodArray)
Print the likelihood array to terminal (debugging tool).
std::vector< double > rootFreqs_
virtual bool hasFather() const
Tell if this node has a father node.
virtual bool isLeaf() const
virtual void computeSubtreeLikelihoodPostfix(const Node *node)
ParameterList getNodeParameters() const
Get the parameters corresponding attached to the nodes of the tree.
DRASDRTreeLikelihoodData * likelihoodData_
AbstractNonHomogeneousTreeLikelihood & operator=(const AbstractNonHomogeneousTreeLikelihood &lik)
Assignation operator.
virtual ~DRNonHomogeneousTreeLikelihood()
virtual void computeRootLikelihood()
VVdouble & getLeafLikelihoods(int nodeId)
virtual int getId() const
Get the node's id.
void init_()
Method called by constructors.
double getValue() const
Function and NNISearchable interface.
bool computeFirstOrderDerivatives_
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...
void fireParameterChanged(const ParameterList ¶ms)
std::map< int, VVVdouble > d2pxy_
static void resetLikelihoodArray(VVVdouble &likelihoodArray)
Set all conditional likelihoods to 1.
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
TreeTemplate< Node > * tree_
VVVdouble & getLikelihoodArray(int parentId, int neighborId)
double getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the logarithm of the likelihood for a site knowing its rate class.
std::map< int, const Node * > idToNode_
An index linking nodes to their id, for faster access than the getNode() method.
The phylogenetic node class.
DRNonHomogeneousTreeLikelihood & operator=(const DRNonHomogeneousTreeLikelihood &lik)
virtual void displayLikelihood(const Node *node)
This method is mainly for debugging purpose.
static void computeLikelihoodFromArrays(const std::vector< const VVVdouble *> &iLik, const std::vector< const VVVdouble *> &tProb, VVVdouble &oLik, size_t nbNodes, size_t nbDistinctSites, size_t nbClasses, size_t nbStates, bool reset=true)
Compute conditional likelihoods.
void computeTreeLikelihood()
double getLikelihood() const
Get the likelihood for the whole dataset.
void setParameters(const ParameterList ¶meters)
Implements the Function interface.
virtual void computeSubtreeLikelihoodPrefix(const Node *node)
This class implements the likelihood computation for a tree using the double-recursive algorithm...
std::map< int, VVVdouble > dpxy_
DiscreteDistribution * rateDistribution_
virtual size_t getNumberOfSons() const
virtual void computeTreeDLikelihoods()
Likelihood data structure for rate across sites models, using a double-recursive algorithm.
DRASDRTreeLikelihoodData * clone() const
virtual void computeLikelihoodAtNode_(const Node *node, VVVdouble &likelihoodArray) const
virtual void computeTreeD2LikelihoodAtNode(const Node *node)
double getLikelihoodForASite(size_t site) const
Get the likelihood for a site.
bool isInitialized() const
std::vector< int > getNodesWithParameter(const std::string &name) const
bool computeSecondOrderDerivatives_
void setTree(const TreeTemplate< Node > *tree)
Set the tree associated to the data.
std::map< int, VVVdouble > pxy_
size_t getRootArrayPosition(const size_t site) const
double getSecondOrderDerivative(const std::string &variable) const