bpp-phyl  2.4.0
AbstractDiscreteRatesAcrossSitesTreeLikelihood.cpp
Go to the documentation of this file.
1 //
2 // File: AbstractDiscreteRatesAcrossSitesTreeLikelihood.cpp
3 // Created by: Julien Dutheil
4 // Created on: Wue Jun 15 09:42 2005
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (November 16, 2004)
9 
10  This software is a computer program whose purpose is to provide classes
11  for phylogenetic data analysis.
12 
13  This software is governed by the CeCILL license under French law and
14  abiding by the rules of distribution of free software. You can use,
15  modify and/ or redistribute the software under the terms of the CeCILL
16  license as circulated by CEA, CNRS and INRIA at the following URL
17  "http://www.cecill.info".
18 
19  As a counterpart to the access to the source code and rights to copy,
20  modify and redistribute granted by the license, users are provided only
21  with a limited warranty and the software's author, the holder of the
22  economic rights, and the successive licensors have only limited
23  liability.
24 
25  In this respect, the user's attention is drawn to the risks associated
26  with loading, using, modifying and/or developing or reproducing the
27  software by the user in light of its specific status of free software,
28  that may mean that it is complicated to manipulate, and that also
29  therefore means that it is reserved for developers and experienced
30  professionals having in-depth computer knowledge. Users are therefore
31  encouraged to load and test the software's suitability as regards their
32  requirements in conditions enabling the security of their systems and/or
33  data to be ensured and, more generally, to use and operate it in the
34  same conditions as regards security.
35 
36  The fact that you are presently reading this means that you have had
37  knowledge of the CeCILL license and that you accept its terms.
38  */
39 
41 
43 
44 using namespace bpp;
45 
46 // From the STL:
47 #include <iostream>
48 
49 using namespace std;
50 
51 /******************************************************************************/
52 
54  DiscreteDistribution* rDist,
55  bool verbose) :
56  rateDistribution_(rDist)
57 {
59 }
60 
61 /******************************************************************************/
62 
64 {
65  if (!initialized_)
66  throw Exception("AbstractDiscreteRatesAcrossSitesTreeLikelihood::getRateDistributionParameters(). Object is not initialized.");
68 }
69 
70 /******************************************************************************/
71 
73 {
74  if (!initialized_)
75  throw Exception("AbstractDiscreteRatesAcrossSitesTreeLikelihood::getDerivableParameters(). Object is not initialized.");
77 }
78 
79 /******************************************************************************/
80 
82 {
83  if (!initialized_)
84  throw Exception("AbstractDiscreteRatesAcrossSitesTreeLikelihood::getNonDerivableParameters(). Object is not initialized.");
87  return tmp;
88 }
89 
90 /******************************************************************************/
91 
93 {
94  size_t nbSites = getNumberOfSites();
95  size_t nbClasses = getNumberOfClasses();
96  VVdouble l(nbSites);
97  for (size_t i = 0; i < nbSites; i++)
98  {
99  l[i].resize(nbClasses);
100  for (size_t j = 0; j < nbClasses; j++)
101  {
102  l[i][j] = getLikelihoodForASiteForARateClass(i, j);
103  }
104  }
105  return l;
106 }
107 
108 /******************************************************************************/
109 
111 {
112  size_t nbClasses = getNumberOfClasses();
113  double l = 0;
114  for (size_t i = 0; i < nbClasses; i++)
115  {
117  }
118  return l;
119 }
120 
121 /******************************************************************************/
122 
124 {
125  size_t nbClasses = getNumberOfClasses();
126  double l = 0;
127  for (size_t i = 0; i < nbClasses; i++)
128  {
130  }
131  // if(l <= 0.) cerr << "WARNING!!! Negative likelihood." << endl;
132  return log(l);
133 }
134 
135 /******************************************************************************/
136 
138 {
139  size_t nbSites = getNumberOfSites();
140  size_t nbClasses = getNumberOfClasses();
141  VVdouble l(nbSites);
142  for (size_t i = 0; i < nbSites; i++)
143  {
144  l[i] = Vdouble(nbClasses);
145  for (size_t j = 0; j < nbClasses; j++)
146  {
148  }
149  }
150  return l;
151 }
152 
153 /******************************************************************************/
154 
156 {
157  size_t nbSites = getNumberOfSites();
158  size_t nbClasses = getNumberOfClasses();
159  size_t nbStates = getNumberOfStates();
160  VVVdouble l(nbSites);
161  for (size_t i = 0; i < nbSites; i++)
162  {
163  l[i].resize(nbClasses);
164  for (size_t j = 0; j < nbClasses; j++)
165  {
166  l[i][j].resize(nbStates);
167  for (size_t x = 0; x < nbStates; x++)
168  {
169  l[i][j][x] = getLikelihoodForASiteForARateClassForAState(i, j, static_cast<int>(x));
170  }
171  }
172  }
173  return l;
174 }
175 
176 /******************************************************************************/
177 
179 {
180  size_t nbSites = getNumberOfSites();
181  size_t nbClasses = getNumberOfClasses();
182  size_t nbStates = getNumberOfStates();
183  VVVdouble l(nbSites);
184  for (size_t i = 0; i < nbSites; i++)
185  {
186  l[i].resize(nbClasses);
187  for (size_t j = 0; j < nbClasses; j++)
188  {
189  l[i][j].resize(nbStates);
190  for (size_t x = 0; x < nbStates; x++)
191  {
192  l[i][j][x] = getLogLikelihoodForASiteForARateClassForAState(i, j, static_cast<int>(x));
193  }
194  }
195  }
196  return l;
197 }
198 
199 /*******************************************************************************/
200 
202 {
203  size_t nbSites = getNumberOfSites();
204  size_t nbClasses = getNumberOfClasses();
207  for (size_t i = 0; i < nbSites; i++)
208  {
209  for (size_t j = 0; j < nbClasses; j++)
210  {
211  pb[i][j] = pb[i][j] * rateDistribution_->getProbability(j) / l[i];
212  }
213  }
214  return pb;
215 }
216 
217 /******************************************************************************/
218 
220 {
221  size_t nbSites = getNumberOfSites();
222  size_t nbClasses = getNumberOfClasses();
225  Vdouble rates(nbSites, 0.);
226  for (size_t i = 0; i < nbSites; i++)
227  {
228  for (size_t j = 0; j < nbClasses; j++)
229  {
230  rates[i] += (lr[i][j] / l[i]) * rateDistribution_->getProbability(j) * rateDistribution_->getCategory(j);
231  }
232  }
233  return rates;
234 }
235 
236 /******************************************************************************/
237 
239 {
240  size_t nbSites = getNumberOfSites();
242  vector<size_t> classes(nbSites);
243  for (size_t i = 0; i < nbSites; i++)
244  {
245  classes[i] = VectorTools::whichMax<double>(l[i]);
246  }
247  return classes;
248 }
249 
250 /******************************************************************************/
251 
253 {
254  size_t nbSites = getNumberOfSites();
256  Vdouble rates(nbSites);
257  for (size_t i = 0; i < nbSites; i++)
258  {
259  rates[i] = rateDistribution_->getCategory(VectorTools::whichMax<double>(l[i]));
260  }
261  return rates;
262 }
263 
264 /******************************************************************************/
265 
267  VVVdouble& likelihoodArray)
268 {
269  size_t nbSites = likelihoodArray.size();
270  size_t nbClasses = likelihoodArray[0].size();
271  size_t nbStates = likelihoodArray[0][0].size();
272  for (size_t i = 0; i < nbSites; i++)
273  {
274  for (size_t c = 0; c < nbClasses; c++)
275  {
276  for (size_t s = 0; s < nbStates; s++)
277  {
278  likelihoodArray[i][c][s] = 1.;
279  }
280  }
281  }
282 }
283 
284 /******************************************************************************/
285 
287  const VVVdouble& likelihoodArray)
288 {
289  size_t nbSites = likelihoodArray.size();
290  size_t nbClasses = likelihoodArray[0].size();
291  size_t nbStates = likelihoodArray[0][0].size();
292  for (size_t i = 0; i < nbSites; i++)
293  {
294  cout << "Site " << i << ":" << endl;
295  for (size_t c = 0; c < nbClasses; c++)
296  {
297  cout << "Rate class " << c;
298  for (size_t s = 0; s < nbStates; s++)
299  {
300  cout << "\t" << likelihoodArray[i][c][s];
301  }
302  cout << endl;
303  }
304  cout << endl;
305  }
306 }
307 
308 /******************************************************************************/
309 
311 {
312  VVVdouble p3 = getTransitionProbabilitiesPerRateClass(nodeId, siteIndex);
313  VVdouble p2;
315  p2.resize(getNumberOfStates());
316  for (size_t i = 0; i < p2.size(); i++)
317  {
318  p2[i].resize(getNumberOfStates());
319  for (size_t j = 0; j < p2.size(); j++)
320  {
321  for (size_t k = 0; k < getNumberOfClasses(); k++)
322  {
323  p2[i][j] += p3[k][i][j] * probs[k];
324  }
325  }
326  }
327  return p2;
328 }
329 
330 /******************************************************************************/
331 
virtual ParameterList getBranchLengthsParameters() const =0
Get the branch lengths parameters.
VVVdouble getLogLikelihoodForEachSiteForEachRateClassForEachState() const
Get the logarithm of the likelihood for each site and each rate class and each state.
VVVdouble getLikelihoodForEachSiteForEachRateClassForEachState() const
Get the likelihood for each site and each rate class and each state.
virtual double getProbability(size_t categoryIndex) const =0
virtual const ParameterList & getParameters() const =0
void enableDerivatives(bool yn)
Tell if derivatives must be computed.
STL namespace.
virtual size_t getNumberOfStates() const =0
VVdouble getLikelihoodForEachSiteForEachRateClass() const
Get the likelihood for each site and each rate class.
virtual double getLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const =0
Get the likelihood for a site knowing its rate class and its ancestral state.
std::vector< size_t > getRateClassWithMaxPostProbOfEachSite() const
Get the posterior rate class (the one with maximum posterior probability) for each site...
static void displayLikelihoodArray(const VVVdouble &likelihoodArray)
Print the likelihood array to terminal (debugging tool).
virtual double getLogLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const =0
Get the logarithm of the likelihood for a site knowing its rate class and its ancestral state...
double getLikelihoodForASiteForAState(size_t site, int state) const
Get the likelihood for a site and for a state.
virtual double getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const =0
Get the likelihood for a site knowing its rate class.
virtual void addParameters(const ParameterList &params)
virtual double getCategory(size_t categoryIndex) const =0
VVdouble getLogLikelihoodForEachSiteForEachRateClass() const
Get the logarithm of the likelihood for each site and each rate class.
VVdouble getTransitionProbabilities(int nodeId, size_t siteIndex) const
Retrieves all Pij(t) for a particular branch, defined by the upper node and site. ...
std::vector< double > Vdouble
Vdouble getRateWithMaxPostProbOfEachSite() const
Get the posterior rate (the one with maximum posterior probability) for each site.
ParameterList getNonDerivableParameters() const
All non derivable parameters.
static void resetLikelihoodArray(VVVdouble &likelihoodArray)
Set all conditional likelihoods to 1.
virtual double getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const =0
Get the logarithm of the likelihood for a site knowing its rate class.
virtual ParameterList getCommonParametersWith(const ParameterList &params) const
virtual VVVdouble getTransitionProbabilitiesPerRateClass(int nodeId, size_t siteIndex) const =0
Retrieves all Pij(t) for a particular branch, defined by the upper node.
VVdouble getPosteriorProbabilitiesOfEachRate() const
Get the posterior probability for each site of belonging to a particular rate class.
Vdouble getLikelihoodForEachSite() const
Get the likelihood for each site.
size_t getNumberOfSites() const
Get the number of sites in the dataset.
const ParameterList & getParameters() const
ParameterList getRateDistributionParameters() const
Get the parameters associated to the rate distirbution.
AbstractDiscreteRatesAcrossSitesTreeLikelihood(DiscreteDistribution *rDist, bool verbose=true)
virtual Vdouble getProbabilities() const =0
Vdouble getPosteriorRateOfEachSite() const
Get the posterior rate, i.e. averaged over all classes and weighted with posterior probabilities...
double getLogLikelihoodForASiteForAState(size_t site, int state) const
Get the logarithm of the likelihood for a site and for a state.
std::vector< VVdouble > VVVdouble
std::vector< Vdouble > VVdouble
virtual ParameterList getSubstitutionModelParameters() const =0
Get the parameters associated to substitution model(s).