45 void TreeLikelihoodTools::getAncestralFrequencies(
48 std::map<
int, std::vector<double> >& frequencies,
49 bool alsoForLeaves)
throw (Exception)
51 int currentId = tl.getTree().getRootId();
52 vector<double> currentFreqs = tl.getRootFrequencies(tl.getSiteIndex(site));
53 getAncestralFrequencies_(tl, tl.getSiteIndex(site), currentId, currentFreqs, frequencies, alsoForLeaves);
56 void TreeLikelihoodTools::getAncestralFrequencies(
58 std::map<
int, std::vector<double> >& frequencies,
59 bool alsoForLeaves)
throw (Exception)
61 size_t n = tl.getLikelihoodData()->getNumberOfDistinctSites();
62 size_t ns = tl.getNumberOfStates();
64 map<int, vector<double> > siteFrequencies;
65 for (
size_t i = 0; i < n; ++i)
67 w = tl.getLikelihoodData()->getWeight(i);
70 for (
size_t i = 0; i < n; ++i)
72 w = tl.getLikelihoodData()->getWeight(i);
73 getAncestralFrequencies(tl, i, siteFrequencies, alsoForLeaves);
77 frequencies = siteFrequencies;
79 for (map<
int, vector<double> >::iterator it = frequencies.begin(); it != frequencies.end(); it++)
80 VectorTools::fill(it->second, 0.);
82 map<int, vector<double> >::iterator it = frequencies.begin();
83 map<int, vector<double> >::iterator itSite = siteFrequencies.begin();
84 for (
size_t j = 0; j < frequencies.size(); ++j)
86 for (
size_t k = 0; k < ns; ++k)
87 it->second[k] += itSite->second[k] * w / sumw;
94 void TreeLikelihoodTools::getAncestralFrequencies_(
98 const std::vector<double>& ancestralFrequencies,
99 std::map<
int,std::vector<double> >& frequencies,
100 bool alsoForLeaves)
throw (Exception)
102 if (!tl.getTree().isLeaf(parentId) || alsoForLeaves)
103 frequencies[parentId] = ancestralFrequencies;
104 vector<int> sonsId = tl.getTree().getSonsId(parentId);
105 for (
size_t i = 0; i < sonsId.size(); i++)
107 vector<double> sonFrequencies(tl.getNumberOfStates());
108 VVdouble pijt = tl.getTransitionProbabilities(sonsId[i], siteIndex);
109 for (
size_t j = 0; j < tl.getNumberOfStates(); j++)
112 for (
size_t k = 0; k < tl.getNumberOfStates(); k++)
114 x += pijt[k][j] * ancestralFrequencies[k];
116 sonFrequencies[j] = x;
118 getAncestralFrequencies_(tl, siteIndex, sonsId[i], sonFrequencies, frequencies, alsoForLeaves);
The TreeLikelihood interface.