← Back to branch summary

~ubuntu-branches/ubuntu/raring/libbpp-phyl/raring

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Likelihood/DRHomogeneousTreeLikelihood.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/DRHomogeneousTreeLikelihood.cpp
 
1
//
 
2
// File: DRHomogeneousTreeLikelihood.cpp
 
3
// Created by: Julien Dutheil
 
4
// Created on: Fri Oct 17 18:14:51 2003
 
5
//
 
6
 
 
7
/*
 
8
   Copyright or © or Copr. CNRS, (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
 
 
40
#include "DRHomogeneousTreeLikelihood.h"
 
41
#include "../PatternTools.h"
 
42
 
 
43
//From SeqLib:
 
44
#include <Bpp/Seq/SiteTools.h>
 
45
#include <Bpp/Seq/Container/AlignedSequenceContainer.h>
 
46
#include <Bpp/Seq/Container/SequenceContainerTools.h>
 
47
 
 
48
#include <Bpp/Text/TextTools.h>
 
49
#include <Bpp/App/ApplicationTools.h>
 
50
 
 
51
using namespace bpp;
 
52
 
 
53
// From the STL:
 
54
#include <iostream>
 
55
 
 
56
using namespace std;
 
57
 
 
58
/******************************************************************************/
 
59
 
 
60
DRHomogeneousTreeLikelihood::DRHomogeneousTreeLikelihood(
 
61
  const Tree& tree,
 
62
  SubstitutionModel* model,
 
63
  DiscreteDistribution* rDist,
 
64
  bool checkRooted,
 
65
  bool verbose)
 
66
throw (Exception) :
 
67
  AbstractHomogeneousTreeLikelihood(tree, model, rDist, checkRooted, verbose),
 
68
  likelihoodData_(0),
 
69
  minusLogLik_(-1.)
 
70
{
 
71
  init_();
 
72
}
 
73
 
 
74
/******************************************************************************/
 
75
 
 
76
DRHomogeneousTreeLikelihood::DRHomogeneousTreeLikelihood(
 
77
  const Tree& tree,
 
78
  const SiteContainer& data,
 
79
  SubstitutionModel* model,
 
80
  DiscreteDistribution* rDist,
 
81
  bool checkRooted,
 
82
  bool verbose)
 
83
throw (Exception) :
 
84
  AbstractHomogeneousTreeLikelihood(tree, model, rDist, checkRooted, verbose),
 
85
  likelihoodData_(0),
 
86
  minusLogLik_(-1.)
 
87
{
 
88
  init_();
 
89
  setData(data);
 
90
}
 
91
 
 
92
/******************************************************************************/
 
93
 
 
94
void DRHomogeneousTreeLikelihood::init_() throw (Exception)
 
95
{
 
96
  likelihoodData_ = new DRASDRTreeLikelihoodData(
 
97
    tree_,
 
98
    rateDistribution_->getNumberOfCategories());
 
99
}
 
100
 
 
101
/******************************************************************************/
 
102
 
 
103
DRHomogeneousTreeLikelihood::DRHomogeneousTreeLikelihood(const DRHomogeneousTreeLikelihood& lik) :
 
104
  AbstractHomogeneousTreeLikelihood(lik),
 
105
  likelihoodData_(0),
 
106
  minusLogLik_(-1.)
 
107
{
 
108
  likelihoodData_ = dynamic_cast<DRASDRTreeLikelihoodData*>(lik.likelihoodData_->clone());
 
109
  likelihoodData_->setTree(tree_);
 
110
  minusLogLik_ = lik.minusLogLik_;
 
111
}
 
112
 
 
113
/******************************************************************************/
 
114
 
 
115
DRHomogeneousTreeLikelihood& DRHomogeneousTreeLikelihood::operator=(const DRHomogeneousTreeLikelihood& lik)
 
116
{
 
117
  AbstractHomogeneousTreeLikelihood::operator=(lik);
 
118
  if (likelihoodData_) delete likelihoodData_;
 
119
  likelihoodData_ = dynamic_cast<DRASDRTreeLikelihoodData*>(lik.likelihoodData_->clone());
 
120
  likelihoodData_->setTree(tree_);
 
121
  minusLogLik_ = lik.minusLogLik_;
 
122
  return *this;
 
123
}
 
124
 
 
125
/******************************************************************************/
 
126
 
 
127
DRHomogeneousTreeLikelihood::~DRHomogeneousTreeLikelihood()
 
128
{
 
129
  delete likelihoodData_;
 
130
}
 
131
 
 
132
/******************************************************************************/
 
133
 
 
134
void DRHomogeneousTreeLikelihood::setData(const SiteContainer& sites) throw (Exception)
 
135
{
 
136
  if (data_) delete data_;
 
137
  data_ = PatternTools::getSequenceSubset(sites, *tree_->getRootNode());
 
138
  if (verbose_) ApplicationTools::displayTask("Initializing data structure");
 
139
  likelihoodData_->initLikelihoods(*data_, *model_);
 
140
  if (verbose_) ApplicationTools::displayTaskDone();
 
141
 
 
142
  nbSites_ = likelihoodData_->getNumberOfSites();
 
143
  nbDistinctSites_ = likelihoodData_->getNumberOfDistinctSites();
 
144
  nbStates_ = likelihoodData_->getNumberOfStates();
 
145
 
 
146
  if (verbose_) ApplicationTools::displayResult("Number of distinct sites",
 
147
                                                TextTools::toString(nbDistinctSites_));
 
148
  initialized_ = false;
 
149
}
 
150
 
 
151
/******************************************************************************/
 
152
 
 
153
double DRHomogeneousTreeLikelihood::getLikelihood() const
 
154
{
 
155
  double l = 1.;
 
156
  Vdouble* lik = &likelihoodData_->getRootRateSiteLikelihoodArray();
 
157
  const vector<unsigned int> * w = &likelihoodData_->getWeights();
 
158
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
159
  {
 
160
    l *= std::pow((*lik)[i], (int)(*w)[i]);
 
161
  }
 
162
  return l;
 
163
}
 
164
 
 
165
/******************************************************************************/
 
166
 
 
167
double DRHomogeneousTreeLikelihood::getLogLikelihood() const
 
168
{
 
169
  double ll = 0;
 
170
  Vdouble* lik = &likelihoodData_->getRootRateSiteLikelihoodArray();
 
171
  const vector<unsigned int> * w = &likelihoodData_->getWeights();
 
172
  vector<double> la(nbDistinctSites_);
 
173
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
174
  {
 
175
    la[i] = (*w)[i] * log((*lik)[i]);
 
176
  }
 
177
  sort(la.begin(), la.end());
 
178
  for (unsigned int i = nbDistinctSites_; i > 0; i--)
 
179
  {
 
180
    ll += la[i - 1];
 
181
  }
 
182
  return ll;
 
183
}
 
184
 
 
185
/******************************************************************************/
 
186
 
 
187
double DRHomogeneousTreeLikelihood::getLikelihoodForASite(unsigned int site) const
 
188
{
 
189
  return likelihoodData_->getRootRateSiteLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)];
 
190
}
 
191
 
 
192
/******************************************************************************/
 
193
 
 
194
double DRHomogeneousTreeLikelihood::getLogLikelihoodForASite(unsigned int site) const
 
195
{
 
196
  return log(likelihoodData_->getRootRateSiteLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)]);
 
197
}
 
198
 
 
199
/******************************************************************************/
 
200
double DRHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
201
{
 
202
  return likelihoodData_->getRootSiteLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)][rateClass];
 
203
}
 
204
 
 
205
/******************************************************************************/
 
206
 
 
207
double DRHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClass(unsigned int site, unsigned int rateClass) const
 
208
{
 
209
  return log(likelihoodData_->getRootSiteLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)][rateClass]);
 
210
}
 
211
 
 
212
/******************************************************************************/
 
213
 
 
214
double DRHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClassForAState(unsigned int site, unsigned int rateClass, int state) const
 
215
{
 
216
  return likelihoodData_->getRootLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)][rateClass][state];
 
217
}
 
218
 
 
219
/******************************************************************************/
 
220
 
 
221
double DRHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClassForAState(unsigned int site, unsigned int rateClass, int state) const
 
222
{
 
223
  return log(likelihoodData_->getRootLikelihoodArray()[likelihoodData_->getRootArrayPosition(site)][rateClass][state]);
 
224
}
 
225
 
 
226
/******************************************************************************/
 
227
 
 
228
void DRHomogeneousTreeLikelihood::setParameters(const ParameterList& parameters)
 
229
throw (ParameterNotFoundException, ConstraintException)
 
230
{
 
231
  setParametersValues(parameters);
 
232
}
 
233
 
 
234
/******************************************************************************/
 
235
 
 
236
void DRHomogeneousTreeLikelihood::fireParameterChanged(const ParameterList& params)
 
237
{
 
238
  applyParameters();
 
239
 
 
240
  if (rateDistribution_->getParameters().getCommonParametersWith(params).size() > 0
 
241
      || model_->getParameters().getCommonParametersWith(params).size() > 0)
 
242
  {
 
243
    //Rate parameter changed, need to recompute all probs:
 
244
    computeAllTransitionProbabilities();
 
245
  }
 
246
  else if (params.size() > 0)
 
247
  {
 
248
    //We may save some computations:
 
249
    for (unsigned int i = 0; i < params.size(); i++)
 
250
    {
 
251
      string s = params[i].getName();
 
252
      if (s.substr(0, 5) == "BrLen")
 
253
      {
 
254
        //Branch length parameter:
 
255
        computeTransitionProbabilitiesForNode(nodes_[TextTools::to < unsigned int > (s.substr(5))]);
 
256
      }
 
257
    }
 
258
  }
 
259
 
 
260
  computeTreeLikelihood();
 
261
  if (computeFirstOrderDerivatives_)
 
262
  {
 
263
    computeTreeDLikelihoods();
 
264
  }
 
265
  if (computeSecondOrderDerivatives_)
 
266
  {
 
267
    computeTreeD2Likelihoods();
 
268
  }
 
269
 
 
270
  minusLogLik_ = -getLogLikelihood();
 
271
}
 
272
 
 
273
/******************************************************************************/
 
274
 
 
275
double DRHomogeneousTreeLikelihood::getValue() const
 
276
throw (Exception)
 
277
{
 
278
  if (!isInitialized()) throw Exception("DRHomogeneousTreeLikelihood::getValue(). Instance is not initialized.");
 
279
  return minusLogLik_;
 
280
}
 
281
 
 
282
/******************************************************************************
 
283
*                           First Order Derivatives                          *
 
284
******************************************************************************/
 
285
 
 
286
void DRHomogeneousTreeLikelihood::computeTreeDLikelihoodAtNode(const Node* node)
 
287
{
 
288
  const Node* father = node->getFather();
 
289
  VVVdouble* _likelihoods_father_node = &likelihoodData_->getLikelihoodArray(father->getId(), node->getId());
 
290
  Vdouble* _dLikelihoods_node = &likelihoodData_->getDLikelihoodArray(node->getId());
 
291
  VVVdouble*  pxy__node = &pxy_[node->getId()];
 
292
  VVVdouble* dpxy__node = &dpxy_[node->getId()];
 
293
  VVVdouble larray;
 
294
  computeLikelihoodAtNode_(father, larray);
 
295
  Vdouble* rootLikelihoodsSR = &likelihoodData_->getRootRateSiteLikelihoodArray();
 
296
 
 
297
  double dLi, dLic, dLicx, numerator, denominator;
 
298
 
 
299
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
300
  {
 
301
    VVdouble* _likelihoods_father_node_i = &(*_likelihoods_father_node)[i];
 
302
    VVdouble* larray_i = &larray[i];
 
303
    dLi = 0;
 
304
    for (unsigned int c = 0; c < nbClasses_; c++)
 
305
    {
 
306
      Vdouble* _likelihoods_father_node_i_c = &(*_likelihoods_father_node_i)[c];
 
307
      Vdouble* larray_i_c = &(*larray_i)[c];
 
308
      VVdouble*  pxy__node_c = &(*pxy__node)[c];
 
309
      VVdouble* dpxy__node_c = &(*dpxy__node)[c];
 
310
      dLic = 0;
 
311
      for (unsigned int x = 0; x < nbStates_; x++)
 
312
      {
 
313
        numerator = 0;
 
314
        denominator = 0;
 
315
        Vdouble*  pxy__node_c_x = &(*pxy__node_c)[x];
 
316
        Vdouble* dpxy__node_c_x = &(*dpxy__node_c)[x];
 
317
        dLicx = 0;
 
318
        for (unsigned int y = 0; y < nbStates_; y++)
 
319
        {
 
320
          numerator   += (*dpxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
 
321
          denominator += (*pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
 
322
        }
 
323
        dLicx = denominator == 0. ? 0. : (*larray_i_c)[x] * numerator / denominator;
 
324
        //cout << i << "\t" << c << "\t" << x << "\t" << (*larray_i_c)[x] << "\t" << numerator << "\t" << denominator << endl;
 
325
        dLic += dLicx;
 
326
      }
 
327
      dLi += rateDistribution_->getProbability(c) * dLic;
 
328
    }
 
329
    (*_dLikelihoods_node)[i] = dLi / (*rootLikelihoodsSR)[i];
 
330
    //cout << dLi << "\t" << (*rootLikelihoodsSR)[i] << endl;
 
331
  }
 
332
}
 
333
 
 
334
/******************************************************************************/
 
335
 
 
336
void DRHomogeneousTreeLikelihood::computeTreeDLikelihoods()
 
337
{
 
338
  for (unsigned int k = 0; k < nbNodes_; k++)
 
339
  {
 
340
    computeTreeDLikelihoodAtNode(nodes_[k]);
 
341
  }
 
342
}
 
343
 
 
344
/******************************************************************************/
 
345
 
 
346
double DRHomogeneousTreeLikelihood::getFirstOrderDerivative(const std::string& variable) const
 
347
throw (Exception)
 
348
{
 
349
  if (!hasParameter(variable))
 
350
    throw ParameterNotFoundException("DRHomogeneousTreeLikelihood::getFirstOrderDerivative().", variable);
 
351
  if (getRateDistributionParameters().hasParameter(variable))
 
352
  {
 
353
    throw Exception("Derivatives respective to rate distribution parameters are not implemented.");
 
354
  }
 
355
  if (getSubstitutionModelParameters().hasParameter(variable))
 
356
  {
 
357
    throw Exception("Derivatives respective to substitution model parameters are not implemented.");
 
358
  }
 
359
 
 
360
  //
 
361
  // Computation for branch lengths:
 
362
  //
 
363
 
 
364
  // Get the node with the branch whose length must be derivated:
 
365
  unsigned int brI = TextTools::to<unsigned int>(variable.substr(5));
 
366
  const Node* branch = nodes_[brI];
 
367
  Vdouble* _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(branch->getId());
 
368
  double d = 0;
 
369
  const vector<unsigned int> * w = &likelihoodData_->getWeights();
 
370
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
371
  {
 
372
    d += (*w)[i] * (*_dLikelihoods_branch)[i];
 
373
  }
 
374
  return -d;
 
375
}
 
376
 
 
377
/******************************************************************************
 
378
*                           Second Order Derivatives                         *
 
379
******************************************************************************/
 
380
 
 
381
void DRHomogeneousTreeLikelihood::computeTreeD2LikelihoodAtNode(const Node* node)
 
382
{
 
383
  const Node* father = node->getFather();
 
384
  VVVdouble* _likelihoods_father_node = &likelihoodData_->getLikelihoodArray(father->getId(), node->getId());
 
385
  Vdouble* _d2Likelihoods_node = &likelihoodData_->getD2LikelihoodArray(node->getId());
 
386
  VVVdouble*   pxy__node = &pxy_[node->getId()];
 
387
  VVVdouble* d2pxy__node = &d2pxy_[node->getId()];
 
388
  VVVdouble larray;
 
389
  computeLikelihoodAtNode_(father, larray);
 
390
  Vdouble* rootLikelihoodsSR = &likelihoodData_->getRootRateSiteLikelihoodArray();
 
391
 
 
392
  double d2Li, d2Lic, d2Licx, numerator, denominator;
 
393
 
 
394
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
395
  {
 
396
    VVdouble* _likelihoods_father_node_i = &(*_likelihoods_father_node)[i];
 
397
    VVdouble* larray_i = &larray[i];
 
398
    d2Li = 0;
 
399
    for (unsigned int c = 0; c < nbClasses_; c++)
 
400
    {
 
401
      Vdouble* _likelihoods_father_node_i_c = &(*_likelihoods_father_node_i)[c];
 
402
      Vdouble* larray_i_c = &(*larray_i)[c];
 
403
      VVdouble*   pxy__node_c = &(*pxy__node)[c];
 
404
      VVdouble* d2pxy__node_c = &(*d2pxy__node)[c];
 
405
      d2Lic = 0;
 
406
      for (unsigned int x = 0; x < nbStates_; x++)
 
407
      {
 
408
        numerator = 0;
 
409
        denominator = 0;
 
410
        Vdouble*   pxy__node_c_x = &(*pxy__node_c)[x];
 
411
        Vdouble* d2pxy__node_c_x = &(*d2pxy__node_c)[x];
 
412
        d2Licx = 0;
 
413
        for (unsigned int y = 0; y < nbStates_; y++)
 
414
        {
 
415
          numerator   += (*d2pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
 
416
          denominator += (*pxy__node_c_x)[y] * (*_likelihoods_father_node_i_c)[y];
 
417
        }
 
418
        d2Licx = denominator == 0. ? 0. : (*larray_i_c)[x] * numerator / denominator;
 
419
        d2Lic += d2Licx;
 
420
      }
 
421
      d2Li += rateDistribution_->getProbability(c) * d2Lic;
 
422
    }
 
423
    (*_d2Likelihoods_node)[i] = d2Li / (*rootLikelihoodsSR)[i];
 
424
  }
 
425
}
 
426
 
 
427
/******************************************************************************/
 
428
 
 
429
void DRHomogeneousTreeLikelihood::computeTreeD2Likelihoods()
 
430
{
 
431
  for (unsigned int k = 0; k < nbNodes_; k++)
 
432
  {
 
433
    computeTreeD2LikelihoodAtNode(nodes_[k]);
 
434
  }
 
435
}
 
436
 
 
437
/******************************************************************************/
 
438
 
 
439
double DRHomogeneousTreeLikelihood::getSecondOrderDerivative(const std::string& variable) const
 
440
throw (Exception)
 
441
{
 
442
  if (!hasParameter(variable))
 
443
    throw ParameterNotFoundException("DRHomogeneousTreeLikelihood::getSecondOrderDerivative().", variable);
 
444
  if (getRateDistributionParameters().hasParameter(variable))
 
445
  {
 
446
    throw Exception("Derivatives respective to rate distribution parameters are not implemented.");
 
447
  }
 
448
  if (getSubstitutionModelParameters().hasParameter(variable))
 
449
  {
 
450
    throw Exception("Derivatives respective to substitution model parameters are not implemented.");
 
451
  }
 
452
 
 
453
  //
 
454
  // Computation for branch lengths:
 
455
  //
 
456
 
 
457
  // Get the node with the branch whose length must be derivated:
 
458
  unsigned int brI = TextTools::to<unsigned int>(variable.substr(5));
 
459
  const Node* branch = nodes_[brI];
 
460
  Vdouble* _dLikelihoods_branch = &likelihoodData_->getDLikelihoodArray(branch->getId());
 
461
  Vdouble* _d2Likelihoods_branch = &likelihoodData_->getD2LikelihoodArray(branch->getId());
 
462
  double d2 = 0;
 
463
  const vector<unsigned int> * w = &likelihoodData_->getWeights();
 
464
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
465
  {
 
466
    d2 += (*w)[i] * ((*_d2Likelihoods_branch)[i] - pow((*_dLikelihoods_branch)[i], 2));
 
467
  }
 
468
  return -d2;
 
469
}
 
470
 
 
471
/******************************************************************************/
 
472
 
 
473
void DRHomogeneousTreeLikelihood::resetLikelihoodArrays(const Node* node)
 
474
{
 
475
  for (unsigned int n = 0; n < node->getNumberOfSons(); n++)
 
476
  {
 
477
    const Node* subNode = node->getSon(n);
 
478
    resetLikelihoodArray(likelihoodData_->getLikelihoodArray(node->getId(), subNode->getId()));
 
479
  }
 
480
  if (node->hasFather())
 
481
  {
 
482
    const Node* father = node->getFather();
 
483
    resetLikelihoodArray(likelihoodData_->getLikelihoodArray(node->getId(), father->getId()));
 
484
  }
 
485
}
 
486
 
 
487
/******************************************************************************/
 
488
 
 
489
void DRHomogeneousTreeLikelihood::computeTreeLikelihood()
 
490
{
 
491
  computeSubtreeLikelihoodPostfix(tree_->getRootNode());
 
492
  computeSubtreeLikelihoodPrefix(tree_->getRootNode());
 
493
  computeRootLikelihood();
 
494
}
 
495
 
 
496
/******************************************************************************/
 
497
 
 
498
void DRHomogeneousTreeLikelihood::computeSubtreeLikelihoodPostfix(const Node* node)
 
499
{
 
500
//  if(node->isLeaf()) return;
 
501
//cout << node->getId() << "\t" << (node->hasName()?node->getName():"") << endl;
 
502
  if (node->getNumberOfSons() == 0) return;
 
503
 
 
504
  // Set all likelihood arrays to 1 for a start:
 
505
  resetLikelihoodArrays(node);
 
506
 
 
507
  map<int, VVVdouble> * _likelihoods_node = &likelihoodData_->getLikelihoodArrays(node->getId());
 
508
  unsigned int nbNodes = node->getNumberOfSons();
 
509
  for (unsigned int l = 0; l < nbNodes; l++)
 
510
  {
 
511
    //For each son node...
 
512
 
 
513
    const Node* son = node->getSon(l);
 
514
    VVVdouble* _likelihoods_node_son = &(*_likelihoods_node)[son->getId()];
 
515
 
 
516
    if (son->isLeaf())
 
517
    {
 
518
      VVdouble* _likelihoods_leaf = &likelihoodData_->getLeafLikelihoods(son->getId());
 
519
      for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
520
      {
 
521
        //For each site in the sequence,
 
522
        Vdouble* _likelihoods_leaf_i = &(*_likelihoods_leaf)[i];
 
523
        VVdouble* _likelihoods_node_son_i = &(*_likelihoods_node_son)[i];
 
524
        for (unsigned int c = 0; c < nbClasses_; c++)
 
525
        {
 
526
          //For each rate classe,
 
527
          Vdouble* _likelihoods_node_son_i_c = &(*_likelihoods_node_son_i)[c];
 
528
          for (unsigned int x = 0; x < nbStates_; x++)
 
529
          {
 
530
            //For each initial state,
 
531
            (*_likelihoods_node_son_i_c)[x] = (*_likelihoods_leaf_i)[x];
 
532
          }
 
533
        }
 
534
      }
 
535
    }
 
536
    else
 
537
    {
 
538
      computeSubtreeLikelihoodPostfix(son); //Recursive method:
 
539
      unsigned int nbSons = son->getNumberOfSons();
 
540
      map<int, VVVdouble> * _likelihoods_son = &likelihoodData_->getLikelihoodArrays(son->getId());
 
541
 
 
542
      vector<const VVVdouble*> iLik(nbSons);
 
543
      vector<const VVVdouble*> tProb(nbSons);
 
544
      for (unsigned int n = 0; n < nbSons; n++)
 
545
      {
 
546
        const Node* sonSon = son->getSon(n);
 
547
        tProb[n] = &pxy_[sonSon->getId()];
 
548
        iLik[n] = &(*_likelihoods_son)[sonSon->getId()];
 
549
      }
 
550
      computeLikelihoodFromArrays(iLik, tProb, *_likelihoods_node_son, nbSons, nbDistinctSites_, nbClasses_, nbStates_, false);
 
551
    }
 
552
  }
 
553
}
 
554
 
 
555
/******************************************************************************/
 
556
 
 
557
void DRHomogeneousTreeLikelihood::computeSubtreeLikelihoodPrefix(const Node* node)
 
558
{
 
559
  if (!node->hasFather())
 
560
  {
 
561
    // 'node' is the root of the tree.
 
562
    // Just call the method on each son node:
 
563
    unsigned int nbSons = node->getNumberOfSons();
 
564
    for (unsigned int n = 0; n < nbSons; n++)
 
565
    {
 
566
      computeSubtreeLikelihoodPrefix(node->getSon(n));
 
567
    }
 
568
    return;
 
569
  }
 
570
  else
 
571
  {
 
572
    const Node* father = node->getFather();
 
573
    map<int, VVVdouble> * _likelihoods_node = &likelihoodData_->getLikelihoodArrays(node->getId());
 
574
    map<int, VVVdouble> * _likelihoods_father = &likelihoodData_->getLikelihoodArrays(father->getId());
 
575
    VVVdouble* _likelihoods_node_father = &(*_likelihoods_node)[father->getId()];
 
576
    if (node->isLeaf())
 
577
    {
 
578
      resetLikelihoodArray(*_likelihoods_node_father);
 
579
    }
 
580
 
 
581
    if (father->isLeaf())
 
582
    {
 
583
      // If the tree is rooted by a leaf
 
584
      VVdouble* _likelihoods_leaf = &likelihoodData_->getLeafLikelihoods(father->getId());
 
585
      for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
586
      {
 
587
        //For each site in the sequence,
 
588
        Vdouble* _likelihoods_leaf_i = &(*_likelihoods_leaf)[i];
 
589
        VVdouble* _likelihoods_node_father_i = &(*_likelihoods_node_father)[i];
 
590
        for (unsigned int c = 0; c < nbClasses_; c++)
 
591
        {
 
592
          //For each rate classe,
 
593
          Vdouble* _likelihoods_node_father_i_c = &(*_likelihoods_node_father_i)[c];
 
594
          for (unsigned int x = 0; x < nbStates_; x++)
 
595
          {
 
596
            //For each initial state,
 
597
            (*_likelihoods_node_father_i_c)[x] = (*_likelihoods_leaf_i)[x];
 
598
          }
 
599
        }
 
600
      }
 
601
    }
 
602
    else
 
603
    {
 
604
      vector<const Node*> nodes;
 
605
      // Add brothers:
 
606
      unsigned int nbFatherSons = father->getNumberOfSons();
 
607
      for (unsigned int n = 0; n < nbFatherSons; n++)
 
608
      {
 
609
        const Node* son = father->getSon(n);
 
610
        if (son->getId() != node->getId()) nodes.push_back(son); //This is a real brother, not current node!
 
611
      }
 
612
      // Now the real stuff... We've got to compute the likelihoods for the
 
613
      // subtree defined by node 'father'.
 
614
      // This is the same as postfix method, but with different subnodes.
 
615
 
 
616
      unsigned int nbSons = nodes.size(); // In case of a bifurcating tree, this is equal to 1, excepted for the root.
 
617
 
 
618
      vector<const VVVdouble*> iLik(nbSons);
 
619
      vector<const VVVdouble*> tProb(nbSons);
 
620
      for (unsigned int n = 0; n < nbSons; n++)
 
621
      {
 
622
        const Node* fatherSon = nodes[n];
 
623
        tProb[n] = &pxy_[fatherSon->getId()];
 
624
        iLik[n] = &(*_likelihoods_father)[fatherSon->getId()];
 
625
      }
 
626
 
 
627
      if (father->hasFather())
 
628
      {
 
629
        const Node* fatherFather = father->getFather();
 
630
        computeLikelihoodFromArrays(iLik, tProb, &(*_likelihoods_father)[fatherFather->getId()], &pxy_[father->getId()], *_likelihoods_node_father, nbSons, nbDistinctSites_, nbClasses_, nbStates_, false);
 
631
      }
 
632
      else
 
633
      {
 
634
        computeLikelihoodFromArrays(iLik, tProb, *_likelihoods_node_father, nbSons, nbDistinctSites_, nbClasses_, nbStates_, false);
 
635
      }
 
636
    }
 
637
 
 
638
    if (!father->hasFather())
 
639
    {
 
640
      //We have to account for the root frequencies:
 
641
      for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
642
      {
 
643
        VVdouble* _likelihoods_node_father_i = &(*_likelihoods_node_father)[i];
 
644
        for (unsigned int c = 0; c < nbClasses_; c++)
 
645
        {
 
646
          Vdouble* _likelihoods_node_father_i_c = &(*_likelihoods_node_father_i)[c];
 
647
          for (unsigned int x = 0; x < nbStates_; x++)
 
648
          {
 
649
            (*_likelihoods_node_father_i_c)[x] *= rootFreqs_[x];
 
650
          }
 
651
        }
 
652
      }
 
653
    }
 
654
 
 
655
    // Call the method on each son node:
 
656
    unsigned int nbNodeSons = node->getNumberOfSons();
 
657
    for (unsigned int i = 0; i < nbNodeSons; i++)
 
658
    {
 
659
      computeSubtreeLikelihoodPrefix(node->getSon(i)); //Recursive method.
 
660
    }
 
661
  }
 
662
}
 
663
 
 
664
/******************************************************************************/
 
665
 
 
666
void DRHomogeneousTreeLikelihood::computeRootLikelihood()
 
667
{
 
668
  const Node* root = tree_->getRootNode();
 
669
  VVVdouble* rootLikelihoods = &likelihoodData_->getRootLikelihoodArray();
 
670
  // Set all likelihoods to 1 for a start:
 
671
  if (root->isLeaf())
 
672
  {
 
673
    VVdouble* leavesLikelihoods_root = &likelihoodData_->getLeafLikelihoods(root->getId());
 
674
    for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
675
    {
 
676
      VVdouble* rootLikelihoods_i = &(*rootLikelihoods)[i];
 
677
      Vdouble* leavesLikelihoods_root_i = &(*leavesLikelihoods_root)[i];
 
678
      for (unsigned int c = 0; c < nbClasses_; c++)
 
679
      {
 
680
        Vdouble* rootLikelihoods_i_c = &(*rootLikelihoods_i)[c];
 
681
        for (unsigned int x = 0; x < nbStates_; x++)
 
682
        {
 
683
          (*rootLikelihoods_i_c)[x] = (*leavesLikelihoods_root_i)[x];
 
684
        }
 
685
      }
 
686
    }
 
687
  }
 
688
  else
 
689
  {
 
690
    resetLikelihoodArray(*rootLikelihoods);
 
691
  }
 
692
 
 
693
  map<int, VVVdouble> * likelihoods_root = &likelihoodData_->getLikelihoodArrays(root->getId());
 
694
  unsigned int nbNodes = root->getNumberOfSons();
 
695
  vector<const VVVdouble*> iLik(nbNodes);
 
696
  vector<const VVVdouble*> tProb(nbNodes);
 
697
  for (unsigned int n = 0; n < nbNodes; n++)
 
698
  {
 
699
    const Node* son = root->getSon(n);
 
700
    tProb[n] = &pxy_[son->getId()];
 
701
    iLik[n] = &(*likelihoods_root)[son->getId()];
 
702
  }
 
703
  computeLikelihoodFromArrays(iLik, tProb, *rootLikelihoods, nbNodes, nbDistinctSites_, nbClasses_, nbStates_, false);
 
704
 
 
705
  Vdouble p = rateDistribution_->getProbabilities();
 
706
  VVdouble* rootLikelihoodsS  = &likelihoodData_->getRootSiteLikelihoodArray();
 
707
  Vdouble* rootLikelihoodsSR = &likelihoodData_->getRootRateSiteLikelihoodArray();
 
708
  for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
709
  {
 
710
    //For each site in the sequence,
 
711
    VVdouble* rootLikelihoods_i = &(*rootLikelihoods)[i];
 
712
    Vdouble* rootLikelihoodsS_i = &(*rootLikelihoodsS)[i];
 
713
    (*rootLikelihoodsSR)[i] = 0;
 
714
    for (unsigned int c = 0; c < nbClasses_; c++)
 
715
    {
 
716
      //For each rate classe,
 
717
      Vdouble* rootLikelihoods_i_c = &(*rootLikelihoods_i)[c];
 
718
      double* rootLikelihoodsS_i_c = &(*rootLikelihoodsS_i)[c];
 
719
      (*rootLikelihoodsS_i_c) = 0;
 
720
      for (unsigned int x = 0; x < nbStates_; x++)
 
721
      {
 
722
        //For each initial state,
 
723
        (*rootLikelihoodsS_i_c) += rootFreqs_[x] * (*rootLikelihoods_i_c)[x];
 
724
      }
 
725
      (*rootLikelihoodsSR)[i] += p[c] * (*rootLikelihoodsS_i_c);
 
726
    }
 
727
 
 
728
    //Final checking (for numerical errors):
 
729
    if ((*rootLikelihoodsSR)[i] < 0) (*rootLikelihoodsSR)[i] = 0.;
 
730
  }
 
731
}
 
732
 
 
733
/******************************************************************************/
 
734
 
 
735
void DRHomogeneousTreeLikelihood::computeLikelihoodAtNode_(const Node* node, VVVdouble& likelihoodArray) const
 
736
{
 
737
  // const Node * node = tree_->getNode(nodeId);
 
738
  int nodeId = node->getId();
 
739
  likelihoodArray.resize(nbDistinctSites_);
 
740
  map<int, VVVdouble> * likelihoods_node = &likelihoodData_->getLikelihoodArrays(nodeId);
 
741
 
 
742
  //Initialize likelihood array:
 
743
  if (node->isLeaf())
 
744
  {
 
745
    VVdouble* leavesLikelihoods_node = &likelihoodData_->getLeafLikelihoods(nodeId);
 
746
    for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
747
    {
 
748
      VVdouble* likelihoodArray_i = &likelihoodArray[i];
 
749
      Vdouble* leavesLikelihoods_node_i = &(*leavesLikelihoods_node)[i];
 
750
      likelihoodArray_i->resize(nbClasses_);
 
751
      for (unsigned int c = 0; c < nbClasses_; c++)
 
752
      {
 
753
        Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
 
754
        likelihoodArray_i_c->resize(nbStates_);
 
755
        for (unsigned int x = 0; x < nbStates_; x++)
 
756
        {
 
757
          (*likelihoodArray_i_c)[x] = (*leavesLikelihoods_node_i)[x];
 
758
        }
 
759
      }
 
760
    }
 
761
  }
 
762
  else
 
763
  {
 
764
    // Otherwise:
 
765
    // Set all likelihoods to 1 for a start:
 
766
    for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
767
    {
 
768
      VVdouble* likelihoodArray_i = &likelihoodArray[i];
 
769
      likelihoodArray_i->resize(nbClasses_);
 
770
      for (unsigned int c = 0; c < nbClasses_; c++)
 
771
      {
 
772
        Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
 
773
        likelihoodArray_i_c->resize(nbStates_);
 
774
        for (unsigned int x = 0; x < nbStates_; x++)
 
775
        {
 
776
          (*likelihoodArray_i_c)[x] = 1.;
 
777
        }
 
778
      }
 
779
    }
 
780
  }
 
781
 
 
782
  unsigned int nbNodes = node->getNumberOfSons();
 
783
 
 
784
  vector<const VVVdouble*> iLik(nbNodes);
 
785
  vector<const VVVdouble*> tProb(nbNodes);
 
786
  for (unsigned int n = 0; n < nbNodes; n++)
 
787
  {
 
788
    const Node* son = node->getSon(n);
 
789
    tProb[n] = &pxy_[son->getId()];
 
790
    iLik[n] = &(*likelihoods_node)[son->getId()];
 
791
  }
 
792
 
 
793
  if (node->hasFather())
 
794
  {
 
795
    const Node* father = node->getFather();
 
796
    computeLikelihoodFromArrays(iLik, tProb, &(*likelihoods_node)[father->getId()], &pxy_[nodeId], likelihoodArray, nbNodes, nbDistinctSites_, nbClasses_, nbStates_, false);
 
797
  }
 
798
  else
 
799
  {
 
800
    computeLikelihoodFromArrays(iLik, tProb, likelihoodArray, nbNodes, nbDistinctSites_, nbClasses_, nbStates_, false);
 
801
 
 
802
    //We have to account for the equilibrium frequencies:
 
803
    for (unsigned int i = 0; i < nbDistinctSites_; i++)
 
804
    {
 
805
      VVdouble* likelihoodArray_i = &likelihoodArray[i];
 
806
      for (unsigned int c = 0; c < nbClasses_; c++)
 
807
      {
 
808
        Vdouble* likelihoodArray_i_c = &(*likelihoodArray_i)[c];
 
809
        for (unsigned int x = 0; x < nbStates_; x++)
 
810
        {
 
811
          (*likelihoodArray_i_c)[x] *= rootFreqs_[x];
 
812
        }
 
813
      }
 
814
    }
 
815
  }
 
816
}
 
817
 
 
818
/******************************************************************************/
 
819
 
 
820
void DRHomogeneousTreeLikelihood::computeLikelihoodFromArrays(
 
821
  const vector<const VVVdouble*>& iLik,
 
822
  const vector<const VVVdouble*>& tProb,
 
823
  VVVdouble& oLik,
 
824
  unsigned int nbNodes,
 
825
  unsigned int nbDistinctSites,
 
826
  unsigned int nbClasses,
 
827
  unsigned int nbStates,
 
828
  bool reset)
 
829
{
 
830
  if (reset) resetLikelihoodArray(oLik);
 
831
 
 
832
  for (unsigned int n = 0; n < nbNodes; n++)
 
833
  {
 
834
    const VVVdouble* pxy_n = tProb[n];
 
835
    const VVVdouble* iLik_n = iLik[n];
 
836
 
 
837
    for (unsigned int i = 0; i < nbDistinctSites; i++)
 
838
    {
 
839
      //For each site in the sequence,
 
840
      const VVdouble* iLik_n_i = &(*iLik_n)[i];
 
841
      VVdouble* oLik_i = &(oLik)[i];
 
842
 
 
843
      for (unsigned int c = 0; c < nbClasses; c++)
 
844
      {
 
845
        //For each rate classe,
 
846
        const Vdouble* iLik_n_i_c = &(*iLik_n_i)[c];
 
847
        Vdouble* oLik_i_c = &(*oLik_i)[c];
 
848
        const VVdouble* pxy_n_c = &(*pxy_n)[c];
 
849
        for (unsigned int x = 0; x < nbStates; x++)
 
850
        {
 
851
          //For each initial state,
 
852
          const Vdouble* pxy_n_c_x = &(*pxy_n_c)[x];
 
853
          double likelihood = 0;
 
854
          for (unsigned int y = 0; y < nbStates; y++)
 
855
          {
 
856
            likelihood += (*pxy_n_c_x)[y] * (*iLik_n_i_c)[y];
 
857
          }
 
858
          // We store this conditionnal likelihood into the corresponding array:
 
859
          (*oLik_i_c)[x] *= likelihood;
 
860
        }
 
861
      }
 
862
    }
 
863
  }
 
864
}
 
865
 
 
866
/******************************************************************************/
 
867
 
 
868
void DRHomogeneousTreeLikelihood::computeLikelihoodFromArrays(
 
869
  const vector<const VVVdouble*>& iLik,
 
870
  const vector<const VVVdouble*>& tProb,
 
871
  const VVVdouble* iLikR,
 
872
  const VVVdouble* tProbR,
 
873
  VVVdouble& oLik,
 
874
  unsigned int nbNodes,
 
875
  unsigned int nbDistinctSites,
 
876
  unsigned int nbClasses,
 
877
  unsigned int nbStates,
 
878
  bool reset)
 
879
{
 
880
  if (reset) resetLikelihoodArray(oLik);
 
881
 
 
882
  for (unsigned int n = 0; n < nbNodes; n++)
 
883
  {
 
884
    const VVVdouble* pxy_n = tProb[n];
 
885
    const VVVdouble* iLik_n = iLik[n];
 
886
 
 
887
    for (unsigned int i = 0; i < nbDistinctSites; i++)
 
888
    {
 
889
      //For each site in the sequence,
 
890
      const VVdouble* iLik_n_i = &(*iLik_n)[i];
 
891
      VVdouble* oLik_i = &(oLik)[i];
 
892
 
 
893
      for (unsigned int c = 0; c < nbClasses; c++)
 
894
      {
 
895
        //For each rate classe,
 
896
        const Vdouble* iLik_n_i_c = &(*iLik_n_i)[c];
 
897
        Vdouble* oLik_i_c = &(*oLik_i)[c];
 
898
        const VVdouble* pxy_n_c = &(*pxy_n)[c];
 
899
        for (unsigned int x = 0; x < nbStates; x++)
 
900
        {
 
901
          //For each initial state,
 
902
          const Vdouble* pxy_n_c_x = &(*pxy_n_c)[x];
 
903
          double likelihood = 0;
 
904
          for (unsigned int y = 0; y < nbStates; y++)
 
905
          {
 
906
            //cout << "1:" << (* pxy_n_c_x)[y]  << endl;
 
907
            //cout << "2:" << (* iLik_n_i_c)[y] << endl;
 
908
            likelihood += (*pxy_n_c_x)[y] * (*iLik_n_i_c)[y];
 
909
            //cout << i << "\t" << c << "\t" << x << "\t" << y << "\t" <<  (* pxy__son_c_x)[y] << "\t" << (* likelihoods_root_son_i_c)[y] << endl;
 
910
          }
 
911
          // We store this conditionnal likelihood into the corresponding array:
 
912
          (*oLik_i_c)[x] *= likelihood;
 
913
        }
 
914
      }
 
915
    }
 
916
  }
 
917
 
 
918
  // Now deal with the subtree containing the root:
 
919
  for (unsigned int i = 0; i < nbDistinctSites; i++)
 
920
  {
 
921
    //For each site in the sequence,
 
922
    const VVdouble* iLikR_i = &(*iLikR)[i];
 
923
    VVdouble* oLik_i = &(oLik)[i];
 
924
 
 
925
    for (unsigned int c = 0; c < nbClasses; c++)
 
926
    {
 
927
      //For each rate classe,
 
928
      const Vdouble* iLikR_i_c = &(*iLikR_i)[c];
 
929
      Vdouble* oLik_i_c = &(*oLik_i)[c];
 
930
      const VVdouble* pxyR_c = &(*tProbR)[c];
 
931
      for (unsigned int x = 0; x < nbStates; x++)
 
932
      {
 
933
        double likelihood = 0;
 
934
        for (unsigned int y = 0; y < nbStates; y++)
 
935
        {
 
936
          //For each final state,
 
937
          const Vdouble* pxyR_c_y = &(*pxyR_c)[y];
 
938
          likelihood += (*pxyR_c_y)[x] * (*iLikR_i_c)[y];
 
939
        }
 
940
        // We store this conditionnal likelihood into the corresponding array:
 
941
        (*oLik_i_c)[x] *= likelihood;
 
942
      }
 
943
    }
 
944
  }
 
945
}
 
946
 
 
947
/******************************************************************************/
 
948
 
 
949
void DRHomogeneousTreeLikelihood::displayLikelihood(const Node* node)
 
950
{
 
951
  cout << "Likelihoods at node " << node->getId() << ": " << endl;
 
952
  for (unsigned int n = 0; n < node->getNumberOfSons(); n++)
 
953
  {
 
954
    const Node* subNode = node->getSon(n);
 
955
    cout << "Array for sub-node " << subNode->getId() << endl;
 
956
    displayLikelihoodArray(likelihoodData_->getLikelihoodArray(node->getId(), subNode->getId()));
 
957
  }
 
958
  if (node->hasFather())
 
959
  {
 
960
    const Node* father = node->getFather();
 
961
    cout << "Array for father node " << father->getId() << endl;
 
962
    displayLikelihoodArray(likelihoodData_->getLikelihoodArray(node->getId(), father->getId()));
 
963
  }
 
964
  cout << "                                         ***" << endl;
 
965
}
 
966
 
 
967
/*******************************************************************************/
 
968