← Back to branch summary

~ubuntu-branches/ubuntu/precise/libbpp-phyl/precise

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Likelihood/RHomogeneousMixedTreeLikelihood.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Julien Dutheil
  • Date: 2011-06-09 11:00:00 UTC
  • Revision ID: james.westby@ubuntu.com-20110609110000-yvx78svv6w7xxgph
Tags: upstream-2.0.2
Import upstream version 2.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Bpp/Phyl/Likelihood/RHomogeneousMixedTreeLikelihood.cpp
 
1
//
 
2
// File: RHomogeneousMixedTreeLikelihood.h
 
3
// Created by: David Fournier, Laurent Gueguen
 
4
//
 
5
 
 
6
/*
 
7
   Copyright or © or Copr. CNRS, (November 16, 2004)
 
8
 
 
9
   This software is a computer program whose purpose is to provide classes
 
10
   for phylogenetic data analysis.
 
11
 
 
12
   This software is governed by the CeCILL  license under French law and
 
13
   abiding by the rules of distribution of free software.  You can  use,
 
14
   modify and/ or redistribute the software under the terms of the CeCILL
 
15
   license as circulated by CEA, CNRS and INRIA at the following URL
 
16
   "http://www.cecill.info".
 
17
 
 
18
   As a counterpart to the access to the source code and  rights to copy,
 
19
   modify and redistribute granted by the license, users are provided only
 
20
   with a limited warranty  and the software's author,  the holder of the
 
21
   economic rights,  and the successive licensors  have only  limited
 
22
   liability.
 
23
 
 
24
   In this respect, the user's attention is drawn to the risks associated
 
25
   with loading,  using,  modifying and/or developing or reproducing the
 
26
   software by the user in light of its specific status of free software,
 
27
   that may mean  that it is complicated to manipulate,  and  that  also
 
28
   therefore means  that it is reserved for developers  and  experienced
 
29
   professionals having in-depth computer knowledge. Users are therefore
 
30
   encouraged to load and test the software's suitability as regards their
 
31
   requirements in conditions enabling the security of their systems and/or
 
32
   data to be ensured and,  more generally, to use and operate it in the
 
33
   same conditions as regards security.
 
34
 
 
35
   The fact that you are presently reading this means that you have had
 
36
   knowledge of the CeCILL license and that you accept its terms.
 
37
 */
 
38
 
 
39
#include "RHomogeneousMixedTreeLikelihood.h"
 
40
 
 
41
 
 
42
// From the STL:
 
43
#include <iostream>
 
44
 
 
45
#include <math.h>
 
46
#include "../PatternTools.h"
 
47
 
 
48
#include <Bpp/Numeric/VectorTools.h>
 
49
#include <Bpp/App/ApplicationTools.h>
 
50
 
 
51
using namespace bpp;
 
52
using namespace std;
 
53
 
 
54
RHomogeneousMixedTreeLikelihood::RHomogeneousMixedTreeLikelihood(
 
55
  const Tree& tree,
 
56
  SubstitutionModel* model,
 
57
  DiscreteDistribution* rDist,
 
58
  bool checkRooted,
 
59
  bool verbose,
 
60
  bool usePatterns) throw (Exception) :
 
61
  RHomogeneousTreeLikelihood(tree, model, rDist, checkRooted, verbose, usePatterns),
 
62
  treeLikelihoodsContainer_(),
 
63
  probas_()
 
64
{
 
65
  MixedSubstitutionModel* mixedmodel;
 
66
  if ((mixedmodel = dynamic_cast<MixedSubstitutionModel*>(model_)) == 0)
 
67
    throw Exception("Bad model: RHomogeneousMixedTreeLikelihood needs a MixedSubstitutionModel.");
 
68
  unsigned int s = mixedmodel->getNumberOfModels();
 
69
  for (unsigned int i = 0; i < s; i++)
 
70
  {
 
71
    treeLikelihoodsContainer_.push_back(
 
72
      new RHomogeneousTreeLikelihood(tree, mixedmodel->getNModel(i), rDist, checkRooted, false, usePatterns));
 
73
    probas_.push_back(mixedmodel->getNProbability(i));
 
74
  }
 
75
}
 
76
 
 
77
RHomogeneousMixedTreeLikelihood::RHomogeneousMixedTreeLikelihood(
 
78
  const Tree& tree,
 
79
  const SiteContainer& data,
 
80
  SubstitutionModel* model,
 
81
  DiscreteDistribution* rDist,
 
82
  bool checkRooted,
 
83
  bool verbose,
 
84
  bool usePatterns) throw (Exception) :
 
85
  RHomogeneousTreeLikelihood(tree, model, rDist, checkRooted, verbose, usePatterns),
 
86
  treeLikelihoodsContainer_(),
 
87
  probas_()
 
88
{
 
89
  MixedSubstitutionModel* mixedmodel;
 
90
 
 
91
  if ((mixedmodel = dynamic_cast<MixedSubstitutionModel*>(model_)) == 0)
 
92
    throw Exception("Bad model: RHomogeneousMixedTreeLikelihood needs a MixedSubstitutionModel.");
 
93
 
 
94
  unsigned int s = mixedmodel->getNumberOfModels();
 
95
 
 
96
  for (unsigned int i = 0; i < s; i++)
 
97
  {
 
98
    treeLikelihoodsContainer_.push_back(
 
99
      new RHomogeneousTreeLikelihood(tree, mixedmodel->getNModel(i), rDist, checkRooted, false, usePatterns));
 
100
    probas_.push_back(mixedmodel->getNProbability(i));
 
101
  }
 
102
  setData(data);
 
103
}
 
104
 
 
105
RHomogeneousMixedTreeLikelihood& RHomogeneousMixedTreeLikelihood::operator=(const RHomogeneousMixedTreeLikelihood& lik)
 
106
{
 
107
  RHomogeneousTreeLikelihood::operator=(lik);
 
108
 
 
109
  treeLikelihoodsContainer_.clear();
 
110
  probas_.clear();
 
111
 
 
112
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
113
  {
 
114
    treeLikelihoodsContainer_.push_back(lik.treeLikelihoodsContainer_[i]->clone());
 
115
    probas_.push_back(lik.probas_[i]);
 
116
  }
 
117
 
 
118
  return *this;
 
119
}
 
120
 
 
121
 
 
122
RHomogeneousMixedTreeLikelihood::RHomogeneousMixedTreeLikelihood(const RHomogeneousMixedTreeLikelihood& lik) :
 
123
  RHomogeneousTreeLikelihood(lik),
 
124
  treeLikelihoodsContainer_(lik.treeLikelihoodsContainer_.size()),
 
125
  probas_(lik.probas_.size())
 
126
{
 
127
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
128
  {
 
129
    treeLikelihoodsContainer_[i] = lik.treeLikelihoodsContainer_[i]->clone();
 
130
    probas_.push_back(lik.probas_[i]);
 
131
  }
 
132
}
 
133
 
 
134
RHomogeneousMixedTreeLikelihood::~RHomogeneousMixedTreeLikelihood()
 
135
{
 
136
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
137
  {
 
138
    delete treeLikelihoodsContainer_[i];
 
139
  }
 
140
}
 
141
 
 
142
 
 
143
void RHomogeneousMixedTreeLikelihood::initialize() throw (Exception)
 
144
{
 
145
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
146
  {
 
147
    treeLikelihoodsContainer_[i]->initialize();
 
148
  }
 
149
  RHomogeneousTreeLikelihood::initialize();
 
150
}
 
151
 
 
152
void RHomogeneousMixedTreeLikelihood::setData(const SiteContainer& sites) throw (Exception)
 
153
{
 
154
  RHomogeneousTreeLikelihood::setData(sites);
 
155
 
 
156
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
157
  {
 
158
    treeLikelihoodsContainer_[i]->setData(sites);
 
159
  }
 
160
}
 
161
 
 
162
 
 
163
void RHomogeneousMixedTreeLikelihood::fireParameterChanged(const ParameterList& params)
 
164
{
 
165
  applyParameters();
 
166
 
 
167
  MixedSubstitutionModel* mixedmodel = dynamic_cast<MixedSubstitutionModel*>(model_);
 
168
 
 
169
  unsigned int s = mixedmodel->getNumberOfModels();
 
170
 
 
171
  const SubstitutionModel* pm;
 
172
  for (unsigned int i = 0; i < s; i++)
 
173
  {
 
174
    ParameterList pl;
 
175
    pm = mixedmodel->getNModel(i);
 
176
    pl.addParameters(pm->getParameters());
 
177
    pl.includeParameters(getParameters());
 
178
    treeLikelihoodsContainer_[i]->matchParametersValues(pl);
 
179
  }
 
180
  probas_ = mixedmodel->getProbabilities();
 
181
 
 
182
  minusLogLik_ = -getLogLikelihood();
 
183
}
 
184
 
 
185
void RHomogeneousMixedTreeLikelihood::computeTreeLikelihood()
 
186
{
 
187
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
188
  {
 
189
    treeLikelihoodsContainer_[i]->computeTreeLikelihood();
 
190
  }
 
191
}
 
192
 
 
193
/******************************************************************************
 
194
*                           Likelihoods                          *
 
195
******************************************************************************/
 
196
double RHomogeneousMixedTreeLikelihood::getLikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
197
{
 
198
  double res = 0;
 
199
 
 
200
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
201
  {
 
202
    res += treeLikelihoodsContainer_[i]->getLikelihoodForASiteForARateClass(site, rateClass) * probas_[i];
 
203
  }
 
204
 
 
205
  return res;
 
206
}
 
207
 
 
208
double RHomogeneousMixedTreeLikelihood::getLogLikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
209
{
 
210
  double x = getLikelihoodForASiteForARateClass(site, rateClass);
 
211
  if (x < 0)
 
212
    x = 0;
 
213
  return log(x);
 
214
}
 
215
 
 
216
double RHomogeneousMixedTreeLikelihood::getLikelihoodForASiteForARateClassForAState(unsigned int site, unsigned int rateClass, int state) const
 
217
{
 
218
  double res = 0;
 
219
 
 
220
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
221
  {
 
222
    res += treeLikelihoodsContainer_[i]->getLikelihoodForASiteForARateClassForAState(site, rateClass, state) * probas_[i];
 
223
  }
 
224
 
 
225
  return res;
 
226
}
 
227
 
 
228
double RHomogeneousMixedTreeLikelihood::getLogLikelihoodForASiteForARateClassForAState(unsigned int site, unsigned int rateClass, int state) const
 
229
{
 
230
  double x = getLikelihoodForASiteForARateClassForAState(site, rateClass, state);
 
231
  if (x < 0)
 
232
    x = 0;
 
233
  return log(x);
 
234
}
 
235
 
 
236
 
 
237
/******************************************************************************
 
238
*                           First Order Derivatives                          *
 
239
******************************************************************************/
 
240
double RHomogeneousMixedTreeLikelihood::getDLikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
241
{
 
242
  double res = 0;
 
243
 
 
244
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
245
  {
 
246
    res += treeLikelihoodsContainer_[i]->getDLikelihoodForASiteForARateClass(site, rateClass) * probas_[i];
 
247
  }
 
248
 
 
249
  return res;
 
250
}
 
251
 
 
252
void RHomogeneousMixedTreeLikelihood::computeTreeDLikelihood(const string& variable)
 
253
{
 
254
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
255
  {
 
256
    treeLikelihoodsContainer_[i]->computeTreeDLikelihood(variable);
 
257
  }
 
258
}
 
259
 
 
260
/******************************************************************************
 
261
*                           Second Order Derivatives                          *
 
262
******************************************************************************/
 
263
double RHomogeneousMixedTreeLikelihood::getD2LikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
264
{
 
265
  double res = 0;
 
266
 
 
267
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
268
  {
 
269
    res += treeLikelihoodsContainer_[i]->getD2LikelihoodForASiteForARateClass(site, rateClass) * probas_[i];
 
270
  }
 
271
 
 
272
  return res;
 
273
}
 
274
 
 
275
 
 
276
void RHomogeneousMixedTreeLikelihood::computeTreeD2Likelihood(const string& variable)
 
277
{
 
278
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
279
  {
 
280
    treeLikelihoodsContainer_[i]->computeTreeD2Likelihood(variable);
 
281
  }
 
282
}
 
283
 
 
284
void RHomogeneousMixedTreeLikelihood::computeSubtreeLikelihood(const Node* node)
 
285
{
 
286
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
287
  {
 
288
    treeLikelihoodsContainer_[i]->computeSubtreeLikelihood(node);
 
289
  }
 
290
}
 
291
 
 
292
void RHomogeneousMixedTreeLikelihood::computeDownSubtreeDLikelihood(const Node* node)
 
293
{
 
294
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
295
  {
 
296
    treeLikelihoodsContainer_[i]->computeDownSubtreeDLikelihood(node);
 
297
  }
 
298
}
 
299
 
 
300
void RHomogeneousMixedTreeLikelihood::computeDownSubtreeD2Likelihood(const Node* node)
 
301
{
 
302
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
303
  {
 
304
    treeLikelihoodsContainer_[i]->computeDownSubtreeD2Likelihood(node);
 
305
  }
 
306
}
 
307
 
 
308
 
 
309
void RHomogeneousMixedTreeLikelihood::computeAllTransitionProbabilities()
 
310
{
 
311
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
312
  {
 
313
    treeLikelihoodsContainer_[i]->computeAllTransitionProbabilities();
 
314
  }
 
315
}
 
316
 
 
317
 
 
318
void RHomogeneousMixedTreeLikelihood::computeTransitionProbabilitiesForNode(const Node* node)
 
319
{
 
320
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
321
  {
 
322
    treeLikelihoodsContainer_[i]->computeTransitionProbabilitiesForNode(node);
 
323
  }
 
324
}
 
325
 
 
326
 
 
327
void RHomogeneousMixedTreeLikelihood::displayLikelihood(const Node* node)
 
328
{
 
329
  for (unsigned int i = 0; i < treeLikelihoodsContainer_.size(); i++)
 
330
  {
 
331
    treeLikelihoodsContainer_[i]->displayLikelihood(node);
 
332
  }
 
333
}