45 #include "../Likelihood/DRTreeLikelihoodTools.h" 46 #include "../Likelihood/MarginalAncestralStateReconstruction.h" 48 #include <Bpp/Text/TextTools.h> 49 #include <Bpp/App/ApplicationTools.h> 50 #include <Bpp/Numeric/Matrix/MatrixTools.h> 51 #include <Bpp/Numeric/DataTable.h> 52 #include <Bpp/Seq/AlphabetIndex/UserAlphabetIndex1.h> 65 const vector<int>& nodeIds,
67 bool verbose)
throw (Exception)
70 if (!drtl.isInitialized())
71 throw Exception(
"SubstitutionMappingTools::computeSubstitutionVectors(). Likelihood object is not initialized.");
76 const SiteContainer* sequences = drtl.getData();
77 const DiscreteDistribution* rDist = drtl.getRateDistribution();
79 size_t nbSites = sequences->getNumberOfSites();
80 size_t nbDistinctSites = drtl.getLikelihoodData()->getNumberOfDistinctSites();
81 size_t nbStates = sequences->getAlphabet()->getSize();
82 size_t nbClasses = rDist->getNumberOfCategories();
83 size_t nbTypes = substitutionCount.getNumberOfSubstitutionTypes();
84 vector<const Node*> nodes = tree.
getNodes();
85 const vector<size_t>* rootPatternLinks
86 = &drtl.getLikelihoodData()->getRootArrayPositions();
88 size_t nbNodes = nodes.size();
95 drtl.computeLikelihoodAtNode(tree.getRootId(), lik);
96 Vdouble Lr(nbDistinctSites, 0);
97 Vdouble rcProbs = rDist->getProbabilities();
98 Vdouble rcRates = rDist->getCategories();
99 for (
size_t i = 0; i < nbDistinctSites; i++)
101 VVdouble* lik_i = &lik[i];
102 for (
size_t c = 0; c < nbClasses; c++)
104 Vdouble* lik_i_c = &(*lik_i)[c];
105 double rc = rDist->getProbability(c);
106 for (
size_t s = 0; s < nbStates; s++)
108 Lr[i] += (*lik_i_c)[s] * rc;
115 ApplicationTools::displayTask(
"Compute joint node-pairs likelihood",
true);
117 for (
size_t l = 0; l < nbNodes; ++l)
120 const Node* currentNode = nodes[l];
121 if (nodeIds.size() > 0 && !VectorTools::contains(nodeIds, currentNode->
getId()))
129 ApplicationTools::displayGauge(l, nbNodes - 1);
130 VVdouble substitutionsForCurrentNode(nbDistinctSites);
131 for (
size_t i = 0; i < nbDistinctSites; ++i)
133 substitutionsForCurrentNode[i].resize(nbTypes);
137 VVVdouble likelihoodsFatherConstantPart(nbDistinctSites);
138 for (
size_t i = 0; i < nbDistinctSites; i++)
140 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
141 likelihoodsFatherConstantPart_i->resize(nbClasses);
142 for (
size_t c = 0; c < nbClasses; c++)
144 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
145 likelihoodsFatherConstantPart_i_c->resize(nbStates);
146 double rc = rDist->getProbability(c);
147 for (
size_t s = 0; s < nbStates; s++)
151 (*likelihoodsFatherConstantPart_i_c)[s] = rc;
158 for (
size_t n = 0; n < nbSons; n++)
161 if (currentSon->
getId() != currentNode->
getId())
163 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
166 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentSon->
getId()));
174 while (sit->hasNext())
176 size_t i = sit->next();
180 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentSon->
getId(), i);
183 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
184 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
185 for (
size_t c = 0; c < nbClasses; c++)
187 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
188 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
189 VVdouble* pxy_c = &pxy[c];
190 for (
size_t x = 0; x < nbStates; x++)
192 Vdouble* pxy_c_x = &(*pxy_c)[x];
193 double likelihood = 0.;
194 for (
size_t y = 0; y < nbStates; y++)
196 likelihood += (*pxy_c_x)[y] * (*likelihoodsFather_son_i_c)[y];
198 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
208 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
210 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(father->
getId()));
218 while (sit->hasNext())
220 size_t i = sit->next();
224 pxy = drtl.getTransitionProbabilitiesPerRateClass(father->
getId(), i);
227 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
228 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
229 for (
size_t c = 0; c < nbClasses; c++)
231 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
232 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
233 VVdouble* pxy_c = &pxy[c];
234 for (
size_t x = 0; x < nbStates; x++)
236 double likelihood = 0.;
237 for (
size_t y = 0; y < nbStates; y++)
239 Vdouble* pxy_c_x = &(*pxy_c)[y];
240 likelihood += (*pxy_c_x)[x] * (*likelihoodsFather_son_i_c)[y];
242 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
251 for (
size_t i = 0; i < nbDistinctSites; i++)
253 vector<double> freqs = drtl.getRootFrequencies(i);
254 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
255 for (
size_t c = 0; c < nbClasses; c++)
257 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
258 for (
size_t x = 0; x < nbStates; x++)
260 (*likelihoodsFatherConstantPart_i_c)[x] *= freqs[x];
272 const VVVdouble* likelihoodsFather_node = &(drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentNode->
getId()));
273 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentNode->
getId()));
279 substitutionCount.setSubstitutionModel(bmd->
getModel());
281 VVVVdouble nxy(nbClasses);
282 for (
size_t c = 0; c < nbClasses; ++c)
284 VVVdouble* nxy_c = &nxy[c];
285 double rc = rcRates[c];
286 nxy_c->resize(nbTypes);
287 for (
size_t t = 0; t < nbTypes; ++t)
289 VVdouble* nxy_c_t = &(*nxy_c)[t];
290 Matrix<double>* nijt = substitutionCount.getAllNumbersOfSubstitutions(d * rc, t + 1);
291 nxy_c_t->resize(nbStates);
292 for (
size_t x = 0; x < nbStates; ++x)
294 Vdouble* nxy_c_t_x = &(*nxy_c_t)[x];
295 nxy_c_t_x->resize(nbStates);
296 for (
size_t y = 0; y < nbStates; ++y)
298 (*nxy_c_t_x)[y] = (*nijt)(x, y);
308 while (sit->hasNext())
310 size_t i = sit->next();
314 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentNode->
getId(), i);
317 const VVdouble* likelihoodsFather_node_i = &(*likelihoodsFather_node)[i];
318 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
319 for (
size_t c = 0; c < nbClasses; ++c)
321 const Vdouble* likelihoodsFather_node_i_c = &(*likelihoodsFather_node_i)[c];
322 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
323 const VVdouble* pxy_c = &pxy[c];
324 VVVdouble* nxy_c = &nxy[c];
325 for (
size_t x = 0; x < nbStates; ++x)
327 double* likelihoodsFatherConstantPart_i_c_x = &(*likelihoodsFatherConstantPart_i_c)[x];
328 const Vdouble* pxy_c_x = &(*pxy_c)[x];
329 for (
size_t y = 0; y < nbStates; ++y)
331 double likelihood_cxy = (*likelihoodsFatherConstantPart_i_c_x)
333 * (*likelihoodsFather_node_i_c)[y];
335 for (
size_t t = 0; t < nbTypes; ++t)
338 substitutionsForCurrentNode[i][t] += likelihood_cxy * (*nxy_c)[t][x][y];
353 for (
size_t i = 0; i < nbSites; ++i)
355 for (
size_t t = 0; t < nbTypes; ++t)
357 (*substitutions)(l, i, t) = substitutionsForCurrentNode[(*rootPatternLinks)[i]][t] / Lr[(*rootPatternLinks)[i]];
363 if (ApplicationTools::message)
364 *ApplicationTools::message <<
" ";
365 ApplicationTools::displayTaskDone();
368 return substitutions;
376 const vector<int>& nodeIds,
378 bool verbose)
throw (Exception)
381 if (!drtl.isInitialized())
382 throw Exception(
"SubstitutionMappingTools::computeSubstitutionVectors(). Likelihood object is not initialized.");
387 const SiteContainer* sequences = drtl.getData();
388 const DiscreteDistribution* rDist = drtl.getRateDistribution();
390 size_t nbSites = sequences->getNumberOfSites();
391 size_t nbDistinctSites = drtl.getLikelihoodData()->getNumberOfDistinctSites();
392 size_t nbStates = sequences->getAlphabet()->getSize();
393 size_t nbClasses = rDist->getNumberOfCategories();
394 size_t nbTypes = substitutionCount.getNumberOfSubstitutionTypes();
395 vector<const Node*> nodes = tree.
getNodes();
396 const vector<size_t>* rootPatternLinks
397 = &drtl.getLikelihoodData()->getRootArrayPositions();
399 size_t nbNodes = nodes.size();
406 drtl.computeLikelihoodAtNode(tree.getRootId(), lik);
407 Vdouble Lr(nbDistinctSites, 0);
408 Vdouble rcProbs = rDist->getProbabilities();
409 Vdouble rcRates = rDist->getCategories();
410 for (
size_t i = 0; i < nbDistinctSites; i++)
412 VVdouble* lik_i = &lik[i];
413 for (
size_t c = 0; c < nbClasses; c++)
415 Vdouble* lik_i_c = &(*lik_i)[c];
416 double rc = rDist->getProbability(c);
417 for (
size_t s = 0; s < nbStates; s++)
419 Lr[i] += (*lik_i_c)[s] * rc;
426 ApplicationTools::displayTask(
"Compute joint node-pairs likelihood",
true);
428 for (
size_t l = 0; l < nbNodes; ++l)
431 const Node* currentNode = nodes[l];
432 if (nodeIds.size() > 0 && !VectorTools::contains(nodeIds, currentNode->
getId()))
440 ApplicationTools::displayGauge(l, nbNodes - 1);
441 VVdouble substitutionsForCurrentNode(nbDistinctSites);
442 for (
size_t i = 0; i < nbDistinctSites; ++i)
444 substitutionsForCurrentNode[i].resize(nbTypes);
448 VVVdouble likelihoodsFatherConstantPart(nbDistinctSites);
449 for (
size_t i = 0; i < nbDistinctSites; i++)
451 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
452 likelihoodsFatherConstantPart_i->resize(nbClasses);
453 for (
size_t c = 0; c < nbClasses; c++)
455 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
456 likelihoodsFatherConstantPart_i_c->resize(nbStates);
457 double rc = rDist->getProbability(c);
458 for (
size_t s = 0; s < nbStates; s++)
462 (*likelihoodsFatherConstantPart_i_c)[s] = rc;
469 for (
size_t n = 0; n < nbSons; n++)
472 if (currentSon->
getId() != currentNode->
getId())
474 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
477 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentSon->
getId()));
485 while (sit->hasNext())
487 size_t i = sit->next();
491 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentSon->
getId(), i);
494 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
495 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
496 for (
size_t c = 0; c < nbClasses; c++)
498 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
499 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
500 VVdouble* pxy_c = &pxy[c];
501 for (
size_t x = 0; x < nbStates; x++)
503 Vdouble* pxy_c_x = &(*pxy_c)[x];
504 double likelihood = 0.;
505 for (
size_t y = 0; y < nbStates; y++)
507 likelihood += (*pxy_c_x)[y] * (*likelihoodsFather_son_i_c)[y];
509 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
519 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
521 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(father->
getId()));
529 while (sit->hasNext())
531 size_t i = sit->next();
535 pxy = drtl.getTransitionProbabilitiesPerRateClass(father->
getId(), i);
538 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
539 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
540 for (
size_t c = 0; c < nbClasses; c++)
542 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
543 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
544 VVdouble* pxy_c = &pxy[c];
545 for (
size_t x = 0; x < nbStates; x++)
547 double likelihood = 0.;
548 for (
size_t y = 0; y < nbStates; y++)
550 Vdouble* pxy_c_x = &(*pxy_c)[y];
551 likelihood += (*pxy_c_x)[x] * (*likelihoodsFather_son_i_c)[y];
553 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
562 for (
size_t i = 0; i < nbDistinctSites; i++)
564 vector<double> freqs = drtl.getRootFrequencies(i);
565 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
566 for (
size_t c = 0; c < nbClasses; c++)
568 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
569 for (
size_t x = 0; x < nbStates; x++)
571 (*likelihoodsFatherConstantPart_i_c)[x] *= freqs[x];
583 const VVVdouble* likelihoodsFather_node = &(drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentNode->
getId()));
584 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentNode->
getId()));
590 substitutionCount.setSubstitutionModel(modelSet.getModelForNode(currentNode->
getId()));
593 VVVVdouble nxy(nbClasses);
594 for (
size_t c = 0; c < nbClasses; ++c)
596 VVVdouble* nxy_c = &nxy[c];
597 double rc = rcRates[c];
598 nxy_c->resize(nbTypes);
599 for (
size_t t = 0; t < nbTypes; ++t)
601 VVdouble* nxy_c_t = &(*nxy_c)[t];
602 Matrix<double>* nijt = substitutionCount.getAllNumbersOfSubstitutions(d * rc, t + 1);
603 nxy_c_t->resize(nbStates);
604 for (
size_t x = 0; x < nbStates; ++x)
606 Vdouble* nxy_c_t_x = &(*nxy_c_t)[x];
607 nxy_c_t_x->resize(nbStates);
608 for (
size_t y = 0; y < nbStates; ++y)
610 (*nxy_c_t_x)[y] = (*nijt)(x, y);
620 while (sit->hasNext())
622 size_t i = sit->next();
626 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentNode->
getId(), i);
629 const VVdouble* likelihoodsFather_node_i = &(*likelihoodsFather_node)[i];
630 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
631 for (
size_t c = 0; c < nbClasses; ++c)
633 const Vdouble* likelihoodsFather_node_i_c = &(*likelihoodsFather_node_i)[c];
634 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
635 const VVdouble* pxy_c = &pxy[c];
636 VVVdouble* nxy_c = &nxy[c];
637 for (
size_t x = 0; x < nbStates; ++x)
639 double* likelihoodsFatherConstantPart_i_c_x = &(*likelihoodsFatherConstantPart_i_c)[x];
640 const Vdouble* pxy_c_x = &(*pxy_c)[x];
641 for (
size_t y = 0; y < nbStates; ++y)
643 double likelihood_cxy = (*likelihoodsFatherConstantPart_i_c_x)
645 * (*likelihoodsFather_node_i_c)[y];
647 for (
size_t t = 0; t < nbTypes; ++t)
650 substitutionsForCurrentNode[i][t] += likelihood_cxy * (*nxy_c)[t][x][y];
665 for (
size_t i = 0; i < nbSites; ++i)
667 for (
size_t t = 0; t < nbTypes; ++t)
669 (*substitutions)(l, i, t) = substitutionsForCurrentNode[(*rootPatternLinks)[i]][t] / Lr[(*rootPatternLinks)[i]];
675 if (ApplicationTools::message)
676 *ApplicationTools::message <<
" ";
677 ApplicationTools::displayTaskDone();
680 return substitutions;
688 bool verbose)
throw (Exception)
691 if (!drtl.isInitialized())
692 throw Exception(
"SubstitutionMappingTools::computeSubstitutionVectorsNoAveraging(). Likelihood object is not initialized.");
696 const SiteContainer* sequences = drtl.getData();
697 const DiscreteDistribution* rDist = drtl.getRateDistribution();
699 size_t nbSites = sequences->getNumberOfSites();
700 size_t nbDistinctSites = drtl.getLikelihoodData()->getNumberOfDistinctSites();
701 size_t nbStates = sequences->getAlphabet()->getSize();
702 size_t nbClasses = rDist->getNumberOfCategories();
703 size_t nbTypes = substitutionCount.getNumberOfSubstitutionTypes();
704 vector<const Node*> nodes = tree.
getNodes();
705 const vector<size_t>* rootPatternLinks
706 = &drtl.getLikelihoodData()->getRootArrayPositions();
708 size_t nbNodes = nodes.size();
713 Vdouble rcRates = rDist->getCategories();
717 ApplicationTools::displayTask(
"Compute joint node-pairs likelihood",
true);
719 for (
size_t l = 0; l < nbNodes; ++l)
722 const Node* currentNode = nodes[l];
729 ApplicationTools::displayGauge(l, nbNodes - 1);
730 VVdouble substitutionsForCurrentNode(nbDistinctSites);
731 for (
size_t i = 0; i < nbDistinctSites; ++i)
733 substitutionsForCurrentNode[i].resize(nbTypes);
737 VVVdouble likelihoodsFatherConstantPart(nbDistinctSites);
738 for (
size_t i = 0; i < nbDistinctSites; ++i)
740 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
741 likelihoodsFatherConstantPart_i->resize(nbClasses);
742 for (
size_t c = 0; c < nbClasses; ++c)
744 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
745 likelihoodsFatherConstantPart_i_c->resize(nbStates);
746 double rc = rDist->getProbability(c);
747 for (
size_t s = 0; s < nbStates; ++s)
751 (*likelihoodsFatherConstantPart_i_c)[s] = rc;
758 for (
size_t n = 0; n < nbSons; ++n)
761 if (currentSon->
getId() != currentNode->
getId())
763 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
766 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentSon->
getId()));
774 while (sit->hasNext())
776 size_t i = sit->next();
780 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentSon->
getId(), i);
783 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
784 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
785 for (
size_t c = 0; c < nbClasses; ++c)
787 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
788 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
789 VVdouble* pxy_c = &pxy[c];
790 for (
size_t x = 0; x < nbStates; ++x)
792 Vdouble* pxy_c_x = &(*pxy_c)[x];
793 double likelihood = 0.;
794 for (
size_t y = 0; y < nbStates; ++y)
796 likelihood += (*pxy_c_x)[y] * (*likelihoodsFather_son_i_c)[y];
798 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
808 const VVVdouble* likelihoodsFather_son = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentSon->
getId());
810 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(father->
getId()));
818 while (sit->hasNext())
820 size_t i = sit->next();
824 pxy = drtl.getTransitionProbabilitiesPerRateClass(father->
getId(), i);
827 const VVdouble* likelihoodsFather_son_i = &(*likelihoodsFather_son)[i];
828 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
829 for (
size_t c = 0; c < nbClasses; ++c)
831 const Vdouble* likelihoodsFather_son_i_c = &(*likelihoodsFather_son_i)[c];
832 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
833 VVdouble* pxy_c = &pxy[c];
834 for (
size_t x = 0; x < nbStates; ++x)
836 double likelihood = 0.;
837 for (
size_t y = 0; y < nbStates; ++y)
839 Vdouble* pxy_c_x = &(*pxy_c)[y];
840 likelihood += (*pxy_c_x)[x] * (*likelihoodsFather_son_i_c)[y];
842 (*likelihoodsFatherConstantPart_i_c)[x] *= likelihood;
851 for (
size_t i = 0; i < nbDistinctSites; ++i)
853 vector<double> freqs = drtl.getRootFrequencies(i);
854 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
855 for (
size_t c = 0; c < nbClasses; ++c)
857 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
858 for (
size_t x = 0; x < nbStates; ++x)
860 (*likelihoodsFatherConstantPart_i_c)[x] *= freqs[x];
871 const VVVdouble* likelihoodsFather_node = &drtl.getLikelihoodData()->getLikelihoodArray(father->
getId(), currentNode->
getId());
872 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentNode->
getId()));
878 substitutionCount.setSubstitutionModel(bmd->
getModel());
880 VVVVdouble nxy(nbClasses);
881 for (
size_t c = 0; c < nbClasses; ++c)
883 double rc = rcRates[c];
884 VVVdouble* nxy_c = &nxy[c];
885 nxy_c->resize(nbTypes);
886 for (
size_t t = 0; t < nbTypes; ++t)
888 VVdouble* nxy_c_t = &(*nxy_c)[t];
889 nxy_c_t->resize(nbStates);
890 Matrix<double>* nijt = substitutionCount.getAllNumbersOfSubstitutions(d * rc, t + 1);
891 for (
size_t x = 0; x < nbStates; ++x)
893 Vdouble* nxy_c_t_x = &(*nxy_c_t)[x];
894 nxy_c_t_x->resize(nbStates);
895 for (
size_t y = 0; y < nbStates; ++y)
897 (*nxy_c_t_x)[y] = (*nijt)(x, y);
907 while (sit->hasNext())
909 size_t i = sit->next();
913 pxy = drtl.getTransitionProbabilitiesPerRateClass(currentNode->
getId(), i);
916 const VVdouble* likelihoodsFather_node_i = &(*likelihoodsFather_node)[i];
917 VVdouble* likelihoodsFatherConstantPart_i = &likelihoodsFatherConstantPart[i];
918 RowMatrix<double> pairProbabilities(nbStates, nbStates);
919 MatrixTools::fill(pairProbabilities, 0.);
920 VVVdouble subsCounts(nbStates);
921 for (
size_t j = 0; j < nbStates; ++j)
923 subsCounts[j].resize(nbStates);
924 for (
size_t k = 0; k < nbStates; ++k)
926 subsCounts[j][k].resize(nbTypes);
929 for (
size_t c = 0; c < nbClasses; ++c)
931 const Vdouble* likelihoodsFather_node_i_c = &(*likelihoodsFather_node_i)[c];
932 Vdouble* likelihoodsFatherConstantPart_i_c = &(*likelihoodsFatherConstantPart_i)[c];
933 const VVdouble* pxy_c = &pxy[c];
934 VVVdouble* nxy_c = &nxy[c];
935 for (
size_t x = 0; x < nbStates; ++x)
937 double* likelihoodsFatherConstantPart_i_c_x = &(*likelihoodsFatherConstantPart_i_c)[x];
938 const Vdouble* pxy_c_x = &(*pxy_c)[x];
939 for (
size_t y = 0; y < nbStates; ++y)
941 double likelihood_cxy = (*likelihoodsFatherConstantPart_i_c_x)
943 * (*likelihoodsFather_node_i_c)[y];
944 pairProbabilities(x, y) += likelihood_cxy;
945 for (
size_t t = 0; t < nbTypes; ++t)
947 subsCounts[x][y][t] += likelihood_cxy * (*nxy_c)[t][x][y];
955 vector<size_t> xy = MatrixTools::whichMax(pairProbabilities);
956 for (
size_t t = 0; t < nbTypes; ++t)
958 substitutionsForCurrentNode[i][t] += subsCounts[xy[0]][xy[1]][t] / pairProbabilities(xy[0], xy[1]);
963 for (
size_t i = 0; i < nbSites; i++)
965 for (
size_t t = 0; t < nbTypes; t++)
967 (*substitutions)(l, i, t) = substitutionsForCurrentNode[(*rootPatternLinks)[i]][t];
973 if (ApplicationTools::message)
974 *ApplicationTools::message <<
" ";
975 ApplicationTools::displayTaskDone();
977 return substitutions;
985 bool verbose)
throw (Exception)
988 if (!drtl.isInitialized())
989 throw Exception(
"SubstitutionMappingTools::computeSubstitutionVectorsNoAveragingMarginal(). Likelihood object is not initialized.");
994 const SiteContainer* sequences = drtl.getData();
995 const DiscreteDistribution* rDist = drtl.getRateDistribution();
996 const Alphabet* alpha = sequences->getAlphabet();
998 size_t nbSites = sequences->getNumberOfSites();
999 size_t nbDistinctSites = drtl.getLikelihoodData()->getNumberOfDistinctSites();
1000 size_t nbStates = alpha->getSize();
1001 size_t nbTypes = substitutionCount.getNumberOfSubstitutionTypes();
1002 vector<const Node*> nodes = tree.
getNodes();
1003 const vector<size_t>* rootPatternLinks
1004 = &drtl.getLikelihoodData()->getRootArrayPositions();
1006 size_t nbNodes = nodes.size();
1013 Vdouble rcRates = rDist->getCategories();
1017 ApplicationTools::displayTask(
"Compute marginal ancestral states");
1021 ApplicationTools::displayTaskDone();
1025 ApplicationTools::displayTask(
"Compute substitution vectors",
true);
1027 for (
size_t l = 0; l < nbNodes; l++)
1029 const Node* currentNode = nodes[l];
1035 vector<size_t> nodeStates = ancestors[currentNode->
getId()];
1036 vector<size_t> fatherStates = ancestors[father->
getId()];
1040 ApplicationTools::displayGauge(l, nbNodes - 1);
1041 VVdouble substitutionsForCurrentNode(nbDistinctSites);
1042 for (
size_t i = 0; i < nbDistinctSites; ++i)
1044 substitutionsForCurrentNode[i].resize(nbTypes);
1052 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentNode->
getId()));
1056 substitutionCount.setSubstitutionModel(bmd->
getModel());
1058 VVVdouble nxyt(nbTypes);
1059 for (
size_t t = 0; t < nbTypes; ++t)
1061 nxyt[t].resize(nbStates);
1062 Matrix<double>* nxy = substitutionCount.getAllNumbersOfSubstitutions(d, t + 1);
1063 for (
size_t x = 0; x < nbStates; ++x)
1065 nxyt[t][x].resize(nbStates);
1066 for (
size_t y = 0; y < nbStates; ++y)
1068 nxyt[t][x][y] = (*nxy)(x, y);
1075 while (sit->hasNext())
1077 size_t i = sit->next();
1078 size_t fatherState = fatherStates[i];
1079 size_t nodeState = nodeStates[i];
1080 if (fatherState >= nbStates || nodeState >= nbStates)
1081 for (
size_t t = 0; t < nbTypes; ++t)
1083 substitutionsForCurrentNode[i][t] = 0;
1086 for (
size_t t = 0; t < nbTypes; ++t)
1088 substitutionsForCurrentNode[i][t] = nxyt[t][fatherState][nodeState];
1094 for (
size_t i = 0; i < nbSites; i++)
1096 for (
size_t t = 0; t < nbTypes; t++)
1098 (*substitutions)(l, i, t) = substitutionsForCurrentNode[(*rootPatternLinks)[i]][t];
1104 if (ApplicationTools::message)
1105 *ApplicationTools::message <<
" ";
1106 ApplicationTools::displayTaskDone();
1108 return substitutions;
1116 bool verbose)
throw (Exception)
1119 if (!drtl.isInitialized())
1120 throw Exception(
"SubstitutionMappingTools::computeSubstitutionVectorsMarginal(). Likelihood object is not initialized.");
1125 const SiteContainer* sequences = drtl.getData();
1126 const DiscreteDistribution* rDist = drtl.getRateDistribution();
1128 size_t nbSites = sequences->getNumberOfSites();
1129 size_t nbDistinctSites = drtl.getLikelihoodData()->getNumberOfDistinctSites();
1130 size_t nbStates = sequences->getAlphabet()->getSize();
1131 size_t nbClasses = rDist->getNumberOfCategories();
1132 size_t nbTypes = substitutionCount.getNumberOfSubstitutionTypes();
1133 vector<const Node*> nodes = tree.
getNodes();
1134 const vector<size_t>* rootPatternLinks
1135 = &drtl.getLikelihoodData()->getRootArrayPositions();
1137 size_t nbNodes = nodes.size();
1144 Vdouble rcProbs = rDist->getProbabilities();
1145 Vdouble rcRates = rDist->getCategories();
1149 ApplicationTools::displayTask(
"Compute marginal node-pairs likelihoods",
true);
1151 for (
size_t l = 0; l < nbNodes; l++)
1153 const Node* currentNode = nodes[l];
1161 ApplicationTools::displayGauge(l, nbNodes - 1);
1162 VVdouble substitutionsForCurrentNode(nbDistinctSites);
1163 for (
size_t i = 0; i < nbDistinctSites; ++i)
1165 substitutionsForCurrentNode[i].resize(nbTypes);
1174 auto_ptr<TreeLikelihood::ConstBranchModelIterator> mit(drtl.getNewBranchModelIterator(currentNode->
getId()));
1175 while (mit->hasNext())
1178 substitutionCount.setSubstitutionModel(bmd->
getModel());
1180 VVVVdouble nxy(nbClasses);
1181 for (
size_t c = 0; c < nbClasses; ++c)
1183 VVVdouble* nxy_c = &nxy[c];
1184 double rc = rcRates[c];
1185 nxy_c->resize(nbTypes);
1186 for (
size_t t = 0; t < nbTypes; ++t)
1188 VVdouble* nxy_c_t = &(*nxy_c)[t];
1189 Matrix<double>* nijt = substitutionCount.getAllNumbersOfSubstitutions(d * rc, t + 1);
1190 nxy_c_t->resize(nbStates);
1191 for (
size_t x = 0; x < nbStates; ++x)
1193 Vdouble* nxy_c_t_x = &(*nxy_c_t)[x];
1194 nxy_c_t_x->resize(nbStates);
1195 for (
size_t y = 0; y < nbStates; ++y)
1197 (*nxy_c_t_x)[y] = (*nijt)(x, y);
1206 while (sit->hasNext())
1208 size_t i = sit->next();
1209 VVdouble* probsNode_i = &probsNode[i];
1210 VVdouble* probsFather_i = &probsFather[i];
1211 for (
size_t c = 0; c < nbClasses; ++c)
1213 Vdouble* probsNode_i_c = &(*probsNode_i)[c];
1214 Vdouble* probsFather_i_c = &(*probsFather_i)[c];
1215 VVVdouble* nxy_c = &nxy[c];
1216 for (
size_t x = 0; x < nbStates; ++x)
1218 for (
size_t y = 0; y < nbStates; ++y)
1220 double prob_cxy = (*probsFather_i_c)[x] * (*probsNode_i_c)[y];
1222 for (
size_t t = 0; t < nbTypes; ++t)
1224 substitutionsForCurrentNode[i][t] += prob_cxy * (*nxy_c)[t][x][y];
1239 for (
size_t i = 0; i < nbSites; ++i)
1241 for (
size_t t = 0; t < nbTypes; ++t)
1243 (*substitutions)(l, i, t) = substitutionsForCurrentNode[(*rootPatternLinks)[i]][t];
1249 if (ApplicationTools::message)
1250 *ApplicationTools::message <<
" ";
1251 ApplicationTools::displayTaskDone();
1253 return substitutions;
1260 const SiteContainer& sites,
1266 throw IOException(
"SubstitutionMappingTools::writeToFile. Can't write to stream.");
1269 for (
size_t i = 0; i < substitutions.getNumberOfSites(); i++)
1271 out <<
"\tSite" << sites.getSite(i).getPosition();
1275 for (
size_t j = 0; j < substitutions.getNumberOfBranches(); j++)
1277 out << substitutions.getNode(j)->getId() <<
"\t" << substitutions.getNode(j)->getDistanceToFather();
1278 for (
size_t i = 0; i < substitutions.getNumberOfSites(); i++)
1280 out <<
"\t" << substitutions(j, i, type);
1293 DataTable* data = DataTable::read(in,
"\t",
true, -1);
1294 vector<string> ids = data->getColumn(0);
1295 data->deleteColumn(0);
1296 data->deleteColumn(0);
1298 size_t nbSites = data->getNumberOfColumns();
1299 substitutions.setNumberOfSites(nbSites);
1300 size_t nbBranches = data->getNumberOfRows();
1301 for (
size_t i = 0; i < nbBranches; i++)
1303 int id = TextTools::toInt(ids[i]);
1304 size_t br = substitutions.getNodeIndex(
id);
1305 for (
size_t j = 0; j < nbSites; j++)
1307 substitutions(br, j, type) = TextTools::toDouble((*data)(i, j));
1311 for (
size_t i = 0; i < nbSites; i++)
1313 string siteTxt = data->getColumnName(i);
1315 if (siteTxt.substr(0, 4) ==
"Site")
1316 site = TextTools::to<int>(siteTxt.substr(4));
1318 site = TextTools::to<int>(siteTxt);
1319 substitutions.setSitePosition(i, site);
1324 catch (Exception& e)
1326 throw IOException(
string(
"Bad input file. ") + e.what());
1336 Vdouble v(nbBranches);
1337 for (
size_t l = 0; l < nbBranches; ++l)
1340 for (
size_t t = 0; t < nbTypes; ++t)
1342 v[l] += smap(l, siteIndex, t);
1352 double sumSquare = 0;
1358 sum += smap(l, siteIndex, t);
1360 sumSquare += sum * sum;
1362 return sqrt(sumSquare);
1371 Vdouble v(nbTypes, 0);
1372 for (
size_t i = 0; i < nbSites; ++i)
1374 for (
size_t t = 0; t < nbTypes; ++t)
1376 v[t] += smap(branchIndex, i, t);
1388 Vdouble v(nbTypes, 0);
1389 for (
size_t i = 0; i < nbBranches; ++i)
1391 for (
size_t t = 0; t < nbTypes; ++t)
1393 v[t] += smap(i, siteIndex, t);
1402 const vector<int>& ids,
1414 vector< vector<double> > counts(ids.size());
1418 for (
size_t k = 0; k < ids.size(); ++k)
1420 vector<double> countsf(nbTypes, 0);
1421 vector<double> tmp(nbTypes, 0);
1422 size_t nbIgnored = 0;
1424 for (
size_t i = 0; !error && i < nbSites; ++i)
1427 for (
size_t t = 0; t < nbTypes; ++t)
1429 tmp[t] = (*mapping)(mapping->
getNodeIndex(ids[k]), i, t);
1430 error = isnan(tmp[t]);
1455 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(ids[k]) +
", counts could not be computed.");
1456 for (
size_t t = 0; t < nbTypes; ++t)
1466 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(ids[k]) +
", " + TextTools::toString(nbIgnored) +
" sites (" + TextTools::toString(ceil(static_cast<double>(nbIgnored * 100) / static_cast<double>(nbSites))) +
"%) have been ignored because they are presumably saturated.");
1470 counts[k].resize(countsf.size());
1471 for (
size_t j = 0; j < countsf.size(); ++j)
1473 counts[k][j] = countsf[j];
1484 const vector<int>& ids,
1496 vector< vector<double> > counts(ids.size());
1500 for (
size_t k = 0; k < ids.size(); ++k)
1502 vector<double> countsf(nbTypes, 0);
1503 vector<double> tmp(nbTypes, 0);
1504 size_t nbIgnored = 0;
1506 for (
size_t i = 0; !error && i < nbSites; ++i)
1509 for (
size_t t = 0; t < nbTypes; ++t)
1511 tmp[t] = (*mapping)(mapping->
getNodeIndex(ids[k]), i, t);
1512 error = isnan(tmp[t]);
1537 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(ids[k]) +
", counts could not be computed.");
1538 for (
size_t t = 0; t < nbTypes; ++t)
1548 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(ids[k]) +
", " + TextTools::toString(nbIgnored) +
" sites (" + TextTools::toString(ceil(static_cast<double>(nbIgnored * 100) / static_cast<double>(nbSites))) +
"%) have been ignored because they are presumably saturated.");
1552 counts[k].resize(countsf.size());
1553 for (
size_t j = 0; j < countsf.size(); ++j)
1555 counts[k][j] = countsf[j];
1566 const vector<int>& ids,
1572 size_t nbStates = nullModel->
getAlphabet()->getSize();
1577 vector<UserAlphabetIndex1 > usai(nbTypes, UserAlphabetIndex1(nullModel->
getAlphabet()));
1579 for (
size_t i = 0; i < nbStates; i++)
1581 for (
size_t j = 0; j < nbStates; j++)
1585 size_t nbt = reg.
getType(i, j);
1587 usai[nbt - 1].setIndex(supportedStates[i], usai[nbt - 1].getIndex(supportedStates[i]) + nullModel->
Qij(i, j));
1593 vector< vector<double> > rewards(ids.size());
1595 for (
size_t k = 0; k < ids.size(); ++k)
1597 rewards[k].resize(nbTypes);
1600 for (
size_t nbt = 0; nbt < nbTypes; nbt++)
1606 for (
size_t k = 0; k < ids.size(); ++k)
1609 for (
size_t i = 0; i < nbSites; ++i)
1611 double tmp = (*mapping)(k, i);
1615 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(ids[k]) +
", reward for type " + reg.
getTypeName(nbt + 1) +
" could not be computed.");
1621 rewards[k][nbt] = s;
1633 const vector<int>& ids,
1639 size_t nbStates = nullModelSet->
getAlphabet()->getSize();
1645 vector< vector<double> > rewards(ids.size());
1647 for (
size_t k = 0; k < ids.size(); ++k)
1649 rewards[k].resize(nbTypes);
1652 vector<UserAlphabetIndex1 > usai(nbTypes, UserAlphabetIndex1(nullModelSet->
getAlphabet()));
1654 for (
size_t nbm = 0; nbm < nbModels; nbm++)
1656 vector<int> mids = VectorTools::vectorIntersection(ids, nullModelSet->
getNodesWithModel(nbm));
1663 for (
size_t nbt = 0; nbt < nbTypes; nbt++)
1664 for (
size_t i = 0; i < nbStates; i++)
1665 usai[nbt].setIndex(supportedStates[i], 0);
1667 for (
size_t i = 0; i < nbStates; i++)
1669 for (
size_t j = 0; j < nbStates; j++)
1673 size_t nbt = reg.
getType(i, j);
1675 usai[nbt - 1].setIndex(supportedStates[i], usai[nbt - 1].getIndex(supportedStates[i]) + modn->
Qij(i, j));
1680 for (
size_t nbt = 0; nbt < nbTypes; nbt++)
1686 for (
size_t k = 0; k < mids.size(); k++)
1689 for (
size_t i = 0; i < nbSites; ++i)
1691 double tmp = (*mapping)(mapping->
getNodeIndex(mids[k]), i);
1695 ApplicationTools::displayWarning(
"On branch " + TextTools::toString(mids[k]) +
", reward for type " + reg.
getTypeName(nbt + 1) +
" could not be computed.");
1703 rewards[VectorTools::which(ids, mids[k])][nbt] = s;
1718 const vector<int>& ids,
1724 vector< vector<double> > counts;
1725 vector< vector<double> > factors;
1727 counts = getCountsPerBranch(drtl, ids, model, reg, -1, verbose);
1728 factors = getNormalizationsPerBranch(drtl, ids, nullModel, reg, verbose);
1730 size_t nbTypes = counts[0].size();
1732 for (
size_t k = 0; k < ids.size(); ++k)
1734 for (
size_t t = 0; t < nbTypes; ++t)
1736 if (factors[k][t] != 0)
1737 counts[k][t] /= factors[k][t];
1745 for (
size_t k = 0; k < ids.size(); ++k)
1747 double l = tree.
getNode(ids[k])->getDistanceToFather();
1748 for (
size_t t = 0; t < nbTypes; ++t)
1761 const vector<int>& ids,
1767 vector< vector<double> > counts;
1768 vector< vector<double> > factors;
1770 counts = getCountsPerBranch(drtl, ids, modelSet->
getModel(0), reg, -1, verbose);
1771 factors = getNormalizationsPerBranch(drtl, ids, nullModelSet, reg, verbose);
1773 size_t nbTypes = counts[0].size();
1775 for (
size_t k = 0; k < ids.size(); ++k)
1777 for (
size_t t = 0; t < nbTypes; ++t)
1779 if (factors[k][t] != 0)
1780 counts[k][t] /= factors[k][t];
1788 for (
size_t k = 0; k < ids.size(); ++k)
1790 double l = tree.
getNode(ids[k])->getDistanceToFather();
1791 for (
size_t t = 0; t < nbTypes; ++t)
1804 const vector<int>& ids,
1810 vector< vector<double> > counts = getCountsPerBranch(drtl, ids, model, reg, threshold);
1819 catch (Exception& ex)
1821 throw Exception(
"The stationarity option can only be used with a category substitution register.");
1824 size_t nbTypes = counts[0].size();
1826 for (
size_t k = 0; k < ids.size(); ++k)
1831 for (
size_t i = 0; i < freqs.size(); ++i)
1834 freqsTypes[c - 1] += freqs[i];
1838 double s = VectorTools::sum(counts[k]);
1839 for (
size_t t = 0; t < nbTypes; ++t)
1844 double s2 = VectorTools::sum(counts[k]);
1846 counts[k] = (counts[k] / s2) * s;
1858 const vector<int>& ids,
1866 file.open(filename.c_str());
1869 size_t nbBr = ids.size();
1871 vector<size_t> sdi(nbBr);
1872 for (
size_t i = 0; i < nbBr; ++i)
1878 for (
size_t i = 0; i < nbBr; ++i)
1884 for (
size_t k = 0; k < nbSites; ++k)
1888 for (
size_t i = 0; i < nbBr; ++i)
1890 file <<
"\t" << countsf[sdi[i]];
1900 const string& filenamePrefix,
1902 const vector<int>& ids,
1912 size_t nbBr = ids.size();
1916 string path = filenamePrefix + TextTools::toString(i + 1) + string(
".count");
1917 ApplicationTools::displayResult(
string(
"Output counts of type ") + TextTools::toString(i + 1) +
string(
" to file"), path);
1918 file.open(path.c_str());
1921 for (
size_t k = 0; k < nbBr; ++k)
1927 for (
size_t j = 0; j < nbSites; ++j)
1930 for (
size_t k = 0; k < nbBr; ++k)
1932 file <<
"\t" << (*smap)(smap->
getNodeIndex(ids[k]), j, i);
A pair of SubstitutionModel / SiteIterator.
Likelihood ancestral states reconstruction: marginal method.
Substitution models manager for non-homogeneous / non-reversible models of evolution.
Interface for all substitution models.
virtual const SubstitutionModel * getModel() const =0
const SubstitutionModel * getModel(size_t i) const
Get one model from the set knowing its index.
virtual const std::vector< int > & getAlphabetStates() const =0
std::map< int, std::vector< size_t > > getAllAncestralStates() const
Get all ancestral states for all nodes.
virtual SiteIterator * getNewSiteIterator() const =0
virtual SubstitutionRegister * clone() const =0
Data storage class for probabilistic substitution mappings.
General interface for storing mapping data.
The SubstitutionRegister interface.
virtual const Node * getSon(size_t pos) const
virtual const Alphabet * getAlphabet() const =0
size_t getNumberOfSubstitutionTypes() const
Analytical reward using the eigen decomposition method.
The phylogenetic tree class.
virtual size_t getCategoryFrom(size_t type) const
size_t getNumberOfCategories() const
virtual const Node * getFather() const
Get the father of this node is there is one.
virtual size_t getNumberOfBranches() const =0
virtual bool hasFather() const
Tell if this node has a father node.
Gather states into defined categories, and count the changes between categories.
virtual double getDistanceToFather() const
Get the distance to the father node is there is one, otherwise throw a NodeException.
const Alphabet * getAlphabet() const
virtual size_t getNumberOfSites() const =0
Get the number of sites in the dataset.
virtual int getId() const
Get the node's id.
virtual double Qij(size_t i, size_t j) const =0
virtual const Tree & getTree() const =0
Get the tree (topology and branch lengths).
virtual ConstBranchModelDescription * next()=0
The SubstitutionsCount interface.
virtual size_t getCategory(size_t state) const
The phylogenetic node class.
virtual size_t getNumberOfSubstitutionTypes() const =0
virtual size_t getType(size_t fromState, size_t toState) const =0
Get the substitution type far a given pair of model states.
virtual std::string getTypeName(size_t type) const =0
Get the name of a given substitution type.
const std::vector< int > & getNodesWithModel(size_t i) const
Get a list of nodes id for which the given model is associated.
virtual size_t getNodeIndex(int nodeId) const
virtual size_t getNumberOfSons() const
Interface for double-recursive (DR) implementation of the likelihood computation. ...
size_t getNumberOfModels() const
virtual size_t getNumberOfSites() const =0
size_t getNumberOfSites() const
virtual bool hasNext() const =0
virtual size_t getNumberOfSubstitutionTypes() const =0
virtual N * getNode(int id, bool checkId=false)
virtual std::vector< const N * > getNodes() const