42 #include "../PatternTools.h" 44 #include <Bpp/Text/TextTools.h> 45 #include <Bpp/App/ApplicationTools.h> 59 DiscreteDistribution* rDist,
62 bool reparametrizeRoot)
68 if (!modelSet->isFullySetUpFor(tree))
69 throw Exception(
"RNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
77 const SiteContainer& data,
79 DiscreteDistribution* rDist,
82 bool reparametrizeRoot)
88 if (!modelSet->isFullySetUpFor(tree))
89 throw Exception(
"RNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
100 rateDistribution_->getNumberOfCategories(),
110 minusLogLik_(lik.minusLogLik_)
140 if (data_)
delete data_;
142 if (verbose_) ApplicationTools::displayTask(
"Initializing data structure");
143 likelihoodData_->initLikelihoods(*data_, *modelSet_->getModel(0));
145 if (verbose_) ApplicationTools::displayTaskDone();
147 nbSites_ = likelihoodData_->getNumberOfSites();
148 nbDistinctSites_ = likelihoodData_->getNumberOfDistinctSites();
149 nbStates_ = likelihoodData_->getNumberOfStates();
151 if (verbose_) ApplicationTools::displayResult(
"Number of distinct sites",
152 TextTools::toString(nbDistinctSites_));
153 initialized_ =
false;
161 for (
size_t i = 0; i <
nbSites_; i++)
174 for (
size_t i = 0; i <
nbSites_; i++)
178 sort(la.begin(), la.end());
179 for (
size_t i =
nbSites_; i > 0; i--)
255 throw (ParameterNotFoundException, ConstraintException)
257 setParametersValues(parameters);
266 if (params.getCommonParametersWith(
rateDistribution_->getIndependentParameters()).size() > 0)
274 for (
size_t i = 0; i < tmp.size(); i++)
277 ids = VectorTools::vectorUnion(ids, tmpv);
280 vector<const Node*> nodes;
281 for (
size_t i = 0; i < ids.size(); i++)
285 vector<const Node*> tmpv;
287 for (
size_t i = 0; i < tmp.size(); i++)
289 if (tmp[i] ==
"BrLenRoot" || tmp[i] ==
"RootPosition")
293 tmpv.push_back(
tree_->getRootNode()->getSon(0));
294 tmpv.push_back(
tree_->getRootNode()->getSon(1));
299 tmpv.push_back(
nodes_[TextTools::to < size_t > (tmp[i].substr(5))]);
301 nodes = VectorTools::vectorUnion(nodes, tmpv);
303 for (
size_t i = 0; i < nodes.size(); i++)
319 if (!
isInitialized())
throw Exception(
"RNonHomogeneousTreeLikelihood::getValue(). Instance is not initialized.");
329 size_t rateClass)
const 367 for (
size_t i = 0; i <
nbSites_; i++)
379 if (!hasParameter(variable))
380 throw ParameterNotFoundException(
"RNonHomogeneousTreeLikelihood::getFirstOrderDerivative().", variable);
381 if (getRateDistributionParameters().hasParameter(variable))
383 throw Exception(
"Derivatives respective to rate distribution parameter are not implemented.");
385 if (getSubstitutionModelParameters().hasParameter(variable))
387 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
391 return -getDLogLikelihood();
398 if (variable ==
"BrLenRoot")
400 const Node* father =
tree_->getRootNode();
405 size_t nbSites = _dLikelihoods_father->size();
406 for (
size_t i = 0; i < nbSites; i++)
408 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
411 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
414 (*_dLikelihoods_father_i_c)[s] = 1.;
420 for (
size_t l = 0; l < nbNodes; l++)
432 double pos = getParameterValue(
"RootPosition");
438 for (
size_t i = 0; i < nbSites; i++)
440 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
441 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
442 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
445 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
446 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
447 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
448 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
449 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
450 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
451 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
454 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
455 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
456 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
457 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
458 double dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
461 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
462 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
463 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
464 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
466 double dl = pos * dl1 * l2 + (1. - pos) * dl2 * l1;
467 (*_dLikelihoods_father_i_c)[x] *= dl;
482 VVVdouble* pxy__son = &
pxy_[son->
getId()];
483 for (
size_t i = 0; i < nbSites; i++)
485 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
486 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
489 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
490 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
491 VVdouble* pxy__son_c = &(*pxy__son)[c];
495 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
498 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
500 (*_dLikelihoods_father_i_c)[x] *= dl;
508 else if (variable ==
"RootPosition")
510 const Node* father =
tree_->getRootNode();
515 size_t nbSites = _dLikelihoods_father->size();
516 for (
size_t i = 0; i < nbSites; i++)
518 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
521 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
524 (*_dLikelihoods_father_i_c)[s] = 1.;
530 for (
size_t l = 0; l < nbNodes; l++)
542 double len = getParameterValue(
"BrLenRoot");
548 for (
size_t i = 0; i < nbSites; i++)
550 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
551 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
552 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
555 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
556 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
557 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
558 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
559 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
560 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
561 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
564 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
565 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
566 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
567 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
568 double dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
571 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
572 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
573 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
574 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
576 double dl = len * (dl1 * l2 - dl2 * l1);
577 (*_dLikelihoods_father_i_c)[x] *= dl;
592 VVVdouble* pxy__son = &
pxy_[son->
getId()];
593 for (
size_t i = 0; i < nbSites; i++)
595 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
596 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
599 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
600 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
601 VVdouble* pxy__son_c = &(*pxy__son)[c];
605 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
608 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
610 (*_dLikelihoods_father_i_c)[x] *= dl;
620 size_t brI = TextTools::to<size_t>(variable.substr(5));
627 size_t nbSites = _dLikelihoods_father->size();
628 for (
size_t i = 0; i < nbSites; i++)
630 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
633 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
636 (*_dLikelihoods_father_i_c)[s] = 1.;
642 for (
size_t l = 0; l < nbNodes; l++)
652 for (
size_t i = 0; i < nbSites; i++)
654 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
655 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
658 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
659 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
660 VVdouble* dpxy__son_c = &(*dpxy__son)[c];
664 Vdouble* dpxy__son_c_x = &(*dpxy__son_c)[x];
667 dl += (*dpxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
669 (*_dLikelihoods_father_i_c)[x] *= dl;
676 VVVdouble* pxy__son = &
pxy_[son->
getId()];
677 for (
size_t i = 0; i < nbSites; i++)
679 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
680 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
683 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
684 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
685 VVdouble* pxy__son_c = &(*pxy__son)[c];
689 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
692 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
694 (*_dLikelihoods_father_i_c)[x] *= dl;
712 if (father == 0)
return;
717 size_t nbSites = _dLikelihoods_father->size();
718 for (
size_t i = 0; i < nbSites; i++)
720 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
723 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
726 (*_dLikelihoods_father_i_c)[s] = 1.;
732 for (
size_t l = 0; l < nbNodes; l++)
736 VVVdouble* pxy__son = &
pxy_[son->
getId()];
742 for (
size_t i = 0; i < nbSites; i++)
744 VVdouble* _dLikelihoods_son_i = &(*_dLikelihoods_son)[(*_patternLinks_father_son)[i]];
745 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
748 Vdouble* _dLikelihoods_son_i_c = &(*_dLikelihoods_son_i)[c];
749 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
750 VVdouble* pxy__son_c = &(*pxy__son)[c];
754 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
757 dl += (*pxy__son_c_x)[y] * (*_dLikelihoods_son_i_c)[y];
759 (*_dLikelihoods_father_i_c)[x] *= dl;
767 for (
size_t i = 0; i < nbSites; i++)
769 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
770 VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
773 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
774 Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
775 VVdouble* pxy__son_c = &(*pxy__son)[c];
779 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
782 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
784 (*_dLikelihoods_father_i_c)[x] *= dl;
801 size_t rateClass)
const 839 for (
size_t i = 0; i <
nbSites_; i++)
851 if (!hasParameter(variable))
852 throw ParameterNotFoundException(
"RNonHomogeneousTreeLikelihood::getSecondOrderDerivative().", variable);
853 if (getRateDistributionParameters().hasParameter(variable))
855 throw Exception(
"Derivatives respective to rate distribution parameter are not implemented.");
857 if (getSubstitutionModelParameters().hasParameter(variable))
859 throw Exception(
"Derivatives respective to substitution model parameters are not implemented.");
863 return -getD2LogLikelihood();
870 if (variable ==
"BrLenRoot")
872 const Node* father =
tree_->getRootNode();
877 size_t nbSites = _d2Likelihoods_father->size();
878 for (
size_t i = 0; i < nbSites; i++)
880 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
883 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
886 (*_d2Likelihoods_father_i_c)[s] = 1.;
892 for (
size_t l = 0; l < nbNodes; l++)
904 double pos = getParameterValue(
"RootPosition");
912 for (
size_t i = 0; i < nbSites; i++)
914 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
915 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
916 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
919 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
920 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
921 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
922 VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
923 VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
924 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
925 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
926 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
927 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
930 Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
931 Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
932 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
933 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
934 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
935 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
936 double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
939 d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
940 d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
941 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
942 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
943 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
944 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
946 double d2l = pos * pos * d2l1 * l2 + (1. - pos) * (1. - pos) * d2l2 * l1 + 2 * pos * (1. - pos) * dl1 * dl2;
947 (*_d2Likelihoods_father_i_c)[x] *= d2l;
962 VVVdouble* pxy__son = &
pxy_[son->
getId()];
963 for (
size_t i = 0; i < nbSites; i++)
965 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
966 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
969 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
970 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
971 VVdouble* pxy__son_c = &(*pxy__son)[c];
975 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
978 d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
980 (*_d2Likelihoods_father_i_c)[x] *= d2l;
988 else if (variable ==
"RootPosition")
990 const Node* father =
tree_->getRootNode();
995 size_t nbSites = _d2Likelihoods_father->size();
996 for (
size_t i = 0; i < nbSites; i++)
998 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1001 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1004 (*_d2Likelihoods_father_i_c)[s] = 1.;
1010 for (
size_t l = 0; l < nbNodes; l++)
1022 double len = getParameterValue(
"BrLenRoot");
1030 for (
size_t i = 0; i < nbSites; i++)
1032 VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
1033 VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
1034 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1037 Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
1038 Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
1039 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1040 VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
1041 VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
1042 VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
1043 VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
1044 VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
1045 VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
1048 Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
1049 Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
1050 Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
1051 Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
1052 Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
1053 Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
1054 double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
1057 d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1058 d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1059 dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1060 dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1061 l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1062 l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1064 double d2l = len * len * (d2l1 * l2 + d2l2 * l1 - 2 * dl1 * dl2);
1065 (*_d2Likelihoods_father_i_c)[x] *= d2l;
1080 VVVdouble* pxy__son = &
pxy_[son->
getId()];
1081 for (
size_t i = 0; i < nbSites; i++)
1083 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1084 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1087 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1088 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1089 VVdouble* pxy__son_c = &(*pxy__son)[c];
1093 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1096 d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1098 (*_d2Likelihoods_father_i_c)[x] *= d2l;
1108 size_t brI = TextTools::to<size_t>(variable.substr(5));
1115 size_t nbSites = _d2Likelihoods_father->size();
1116 for (
size_t i = 0; i < nbSites; i++)
1118 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1121 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1124 (*_d2Likelihoods_father_i_c)[s] = 1.;
1130 for (
size_t l = 0; l < nbNodes; l++)
1140 for (
size_t i = 0; i < nbSites; i++)
1142 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1143 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1146 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1147 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1148 VVdouble* d2pxy__son_c = &(*d2pxy__son)[c];
1152 Vdouble* d2pxy__son_c_x = &(*d2pxy__son_c)[x];
1155 d2l += (*d2pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1157 (*_d2Likelihoods_father_i_c)[x] *= d2l;
1164 VVVdouble* pxy__son = &
pxy_[son->
getId()];
1165 for (
size_t i = 0; i < nbSites; i++)
1167 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1168 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1171 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1172 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1173 VVdouble* pxy__son_c = &(*pxy__son)[c];
1177 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1180 d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1182 (*_d2Likelihoods_father_i_c)[x] *= d2l;
1200 if (father == 0)
return;
1205 size_t nbSites = _d2Likelihoods_father->size();
1206 for (
size_t i = 0; i < nbSites; i++)
1208 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1211 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1214 (*_d2Likelihoods_father_i_c)[s] = 1.;
1220 for (
size_t l = 0; l < nbNodes; l++)
1224 VVVdouble* pxy__son = &
pxy_[son->
getId()];
1230 for (
size_t i = 0; i < nbSites; i++)
1232 VVdouble* _d2Likelihoods_son_i = &(*_d2Likelihoods_son)[(*_patternLinks_father_son)[i]];
1233 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1236 Vdouble* _d2Likelihoods_son_i_c = &(*_d2Likelihoods_son_i)[c];
1237 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1238 VVdouble* pxy__son_c = &(*pxy__son)[c];
1242 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1245 d2l += (*pxy__son_c_x)[y] * (*_d2Likelihoods_son_i_c)[y];
1247 (*_d2Likelihoods_father_i_c)[x] *= d2l;
1255 for (
size_t i = 0; i < nbSites; i++)
1257 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1258 VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1261 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1262 Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1263 VVdouble* pxy__son_c = &(*pxy__son)[c];
1267 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1270 dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1272 (*_d2Likelihoods_father_i_c)[x] *= dl;
1294 if (node->
isLeaf())
return;
1301 for (
size_t i = 0; i < nbSites; i++)
1304 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
1308 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
1312 (*_likelihoods_node_i_c)[x] = 1.;
1317 for (
size_t l = 0; l < nbNodes; l++)
1325 VVVdouble* pxy__son = &
pxy_[son->
getId()];
1329 for (
size_t i = 0; i < nbSites; i++)
1332 VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_node_son)[i]];
1333 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
1337 Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1338 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
1339 VVdouble* pxy__son_c = &(*pxy__son)[c];
1343 Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1344 double likelihood = 0;
1347 likelihood += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1349 (*_likelihoods_node_i_c)[x] *= likelihood;
1361 cout <<
"Likelihoods at node " << node->
getName() <<
": " << endl;
1363 cout <<
" ***" << endl;
virtual void computeDownSubtreeDLikelihood(const Node *)
This class implement the 'traditional' way of computing likelihood for a tree, allowing for non-homog...
Substitution models manager for non-homogeneous / non-reversible models of evolution.
Partial implementation for branch non-homogeneous models of the TreeLikelihood interface.
virtual void computeTreeD2Likelihood(const std::string &variable)
ParameterList brLenParameters_
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
virtual std::string getName() const
Get the name associated to this node, if there is one, otherwise throw a NodeException.
size_t getRootArrayPosition(size_t currentPosition) const
void setTree(const TreeTemplate< Node > *tree)
Set the tree associated to the data.
SubstitutionModelSet * modelSet_
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
double getLikelihoodForASite(size_t site) const
Get the likelihood for a site.
virtual double getDLogLikelihood() const
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
std::vector< double > getRootFrequencies() const
double getFirstOrderDerivative(const std::string &variable) const
void fireParameterChanged(const ParameterList ¶ms)
discrete Rate Across Sites, (simple) Recursive likelihood data structure.
virtual ~RNonHomogeneousTreeLikelihood()
virtual const Node * getFather() const
Get the father of this node is there is one.
Interface for phylogenetic tree objects.
RNonHomogeneousTreeLikelihood & operator=(const RNonHomogeneousTreeLikelihood &lik)
static void displayLikelihoodArray(const VVVdouble &likelihoodArray)
Print the likelihood array to terminal (debugging tool).
std::vector< double > rootFreqs_
virtual bool isLeaf() const
ParameterList getNodeParameters() const
Get the parameters corresponding attached to the nodes of the tree.
virtual double getD2LogLikelihood() const
AbstractNonHomogeneousTreeLikelihood & operator=(const AbstractNonHomogeneousTreeLikelihood &lik)
Assignation operator.
VVVdouble & getDLikelihoodArray(int nodeId)
virtual void computeTreeLikelihood()
double getLogLikelihoodForASite(size_t site) const
Get the logarithm of the likelihood for a site.
virtual int getId() const
Get the node's id.
virtual void computeDownSubtreeD2Likelihood(const Node *)
virtual double getD2LikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
double getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the likelihood for a site knowing its rate class.
virtual double getD2LikelihoodForASite(size_t site) const
void setData(const SiteContainer &sites)
Set the dataset for which the likelihood must be evaluated.
std::map< int, VVVdouble > d2pxy_
RNonHomogeneousTreeLikelihood(const Tree &tree, SubstitutionModelSet *modelSet, DiscreteDistribution *rDist, bool verbose=true, bool usePatterns=true, bool reparametrizeRoot=false)
Build a new NonHomogeneousTreeLikelihood object without data.
virtual void computeTreeDLikelihood(const std::string &variable)
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
virtual void computeSubtreeLikelihood(const Node *node)
Compute the likelihood for a subtree defined by the Tree::Node node.
TreeTemplate< Node > * tree_
virtual double getDLikelihoodForASite(size_t site) const
std::map< int, const Node * > idToNode_
An index linking nodes to their id, for faster access than the getNode() method.
DRASRTreeLikelihoodData * likelihoodData_
virtual double getDLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
virtual void displayLikelihood(const Node *node)
This method is mainly for debugging purpose.
The phylogenetic node class.
virtual double getDLogLikelihoodForASite(size_t site) const
double getSecondOrderDerivative(const std::string &variable) const
void init_(bool usePatterns)
Method called by constructors.
DRASRTreeLikelihoodData * clone() const
void setParameters(const ParameterList ¶meters)
Implements the Function interface.
double getLikelihood() const
Get the likelihood for the whole dataset.
double getLogLikelihood() const
Get the logarithm of the likelihood for the whole dataset.
std::map< int, VVVdouble > dpxy_
DiscreteDistribution * rateDistribution_
virtual size_t getNumberOfSons() const
virtual double getD2LogLikelihoodForASite(size_t site) const
const std::vector< size_t > & getArrayPositions(int parentId, int sonId) const
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.
double getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the logarithm of the likelihood for a site knowing its rate class.
VVVdouble & getLikelihoodArray(int nodeId)
bool isInitialized() const
std::vector< int > getNodesWithParameter(const std::string &name) const
VVVdouble & getD2LikelihoodArray(int nodeId)
std::map< int, VVVdouble > pxy_
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...