← 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/AbstractNonHomogeneousTreeLikelihood.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/AbstractNonHomogeneousTreeLikelihood.cpp
 
1
//
 
2
// File: AbstractNonHomogeneousTreeLikelihood.cpp
 
3
// Created by: Julien Dutheil
 
4
// Created on: Tue Oct 09 16:07 2007
 
5
// From file: AbstractHomogeneousTreeLikelihood.cpp
 
6
//
 
7
 
 
8
/*
 
9
Copyright or © or Copr. CNRS, (November 16, 2004)
 
10
 
 
11
This software is a computer program whose purpose is to provide classes
 
12
for phylogenetic data analysis.
 
13
 
 
14
This software is governed by the CeCILL  license under French law and
 
15
abiding by the rules of distribution of free software.  You can  use, 
 
16
modify and/ or redistribute the software under the terms of the CeCILL
 
17
license as circulated by CEA, CNRS and INRIA at the following URL
 
18
"http://www.cecill.info". 
 
19
 
 
20
As a counterpart to the access to the source code and  rights to copy,
 
21
modify and redistribute granted by the license, users are provided only
 
22
with a limited warranty  and the software's author,  the holder of the
 
23
economic rights,  and the successive licensors  have only  limited
 
24
liability. 
 
25
 
 
26
In this respect, the user's attention is drawn to the risks associated
 
27
with loading,  using,  modifying and/or developing or reproducing the
 
28
software by the user in light of its specific status of free software,
 
29
that may mean  that it is complicated to manipulate,  and  that  also
 
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 "AbstractNonHomogeneousTreeLikelihood.h"
 
40
#include "../PatternTools.h"
 
41
 
 
42
//From SeqLib:
 
43
#include <Bpp/Seq/SiteTools.h>
 
44
#include <Bpp/Seq/Container/SequenceContainerTools.h>
 
45
 
 
46
#include <Bpp/Text/TextTools.h>
 
47
#include <Bpp/App/ApplicationTools.h>
 
48
 
 
49
using namespace bpp;
 
50
 
 
51
// From the STL:
 
52
#include <iostream>
 
53
 
 
54
using namespace std;
 
55
 
 
56
/******************************************************************************/
 
57
 
 
58
AbstractNonHomogeneousTreeLikelihood::AbstractNonHomogeneousTreeLikelihood(
 
59
  const Tree& tree,
 
60
  SubstitutionModelSet* modelSet,
 
61
  DiscreteDistribution* rDist,
 
62
  bool verbose)
 
63
  throw (Exception) :
 
64
  AbstractDiscreteRatesAcrossSitesTreeLikelihood(rDist, verbose),
 
65
  modelSet_(0),
 
66
  brLenParameters_(),
 
67
  pxy_(),
 
68
  dpxy_(),
 
69
  d2pxy_(),
 
70
  rootFreqs_(),
 
71
  nodes_(),
 
72
  idToNode_(),
 
73
  nbSites_(),
 
74
  nbDistinctSites_(),
 
75
  nbClasses_(),
 
76
  nbStates_(),
 
77
  nbNodes_(),
 
78
  verbose_(),
 
79
  minimumBrLen_(),
 
80
  brLenConstraint_(0),
 
81
  root1_(),
 
82
  root2_()
 
83
{
 
84
  init_(tree, modelSet, rDist, verbose);
 
85
}
 
86
 
 
87
/******************************************************************************/
 
88
 
 
89
AbstractNonHomogeneousTreeLikelihood::AbstractNonHomogeneousTreeLikelihood(
 
90
    const AbstractNonHomogeneousTreeLikelihood& lik) :
 
91
  AbstractDiscreteRatesAcrossSitesTreeLikelihood(lik),
 
92
  modelSet_(lik.modelSet_),
 
93
  brLenParameters_(lik.brLenParameters_),
 
94
  pxy_(lik.pxy_),
 
95
  dpxy_(lik.dpxy_),
 
96
  d2pxy_(lik.d2pxy_),
 
97
  rootFreqs_(lik.rootFreqs_),
 
98
  nodes_(),
 
99
  idToNode_(),
 
100
        nbSites_(lik.nbSites_),
 
101
  nbDistinctSites_(lik.nbDistinctSites_),
 
102
        nbClasses_(lik.nbClasses_),
 
103
        nbStates_(lik.nbStates_),
 
104
        nbNodes_(lik.nbNodes_),
 
105
  verbose_(lik.verbose_),
 
106
  minimumBrLen_(lik.minimumBrLen_),
 
107
  brLenConstraint_(dynamic_cast<Constraint*>(lik.brLenConstraint_->clone())),
 
108
  root1_(lik.root1_),
 
109
  root2_(lik.root2_)
 
110
{ 
 
111
  nodes_ = tree_->getNodes();
 
112
  nodes_.pop_back(); //Remove the root node (the last added!).  
 
113
  //Rebuild nodes index:
 
114
  for (unsigned int i = 0; i < nodes_.size(); i++)
 
115
  {
 
116
    const Node* node = nodes_[i];
 
117
    idToNode_[node->getId()] = node;
 
118
  }
 
119
}
 
120
 
 
121
/******************************************************************************/
 
122
 
 
123
AbstractNonHomogeneousTreeLikelihood& AbstractNonHomogeneousTreeLikelihood::operator=(
 
124
    const AbstractNonHomogeneousTreeLikelihood& lik)
 
125
{
 
126
  AbstractDiscreteRatesAcrossSitesTreeLikelihood::operator=(lik);
 
127
  modelSet_        = lik.modelSet_;
 
128
  brLenParameters_ = lik.brLenParameters_;
 
129
  pxy_             = lik.pxy_;
 
130
  dpxy_            = lik.dpxy_;
 
131
  d2pxy_           = lik.d2pxy_;
 
132
  rootFreqs_       = lik.rootFreqs_;
 
133
  nodes_           = tree_->getNodes();
 
134
  nodes_.pop_back(); //Remove the root node (the last added!).  
 
135
        nbSites_         = lik.nbSites_;
 
136
  nbDistinctSites_ = lik.nbDistinctSites_;
 
137
        nbClasses_       = lik.nbClasses_;
 
138
        nbStates_        = lik.nbStates_;
 
139
        nbNodes_         = lik.nbNodes_;
 
140
  verbose_         = lik.verbose_;
 
141
  minimumBrLen_    = lik.minimumBrLen_;
 
142
  brLenConstraint_.reset(dynamic_cast<Constraint*>(lik.brLenConstraint_->clone()));
 
143
  root1_           = lik.root1_;
 
144
  root2_           = lik.root2_;
 
145
  //Rebuild nodes index:
 
146
  for( unsigned int i = 0; i < nodes_.size(); i++)
 
147
  {
 
148
    const Node * node = nodes_[i];
 
149
    idToNode_[node->getId()] = node;
 
150
  }
 
151
  return *this;
 
152
}
 
153
 
 
154
/******************************************************************************/
 
155
 
 
156
void AbstractNonHomogeneousTreeLikelihood::init_(
 
157
    const Tree& tree,
 
158
                        SubstitutionModelSet* modelSet,
 
159
                        DiscreteDistribution* rDist,
 
160
                        bool verbose) throw (Exception)
 
161
{
 
162
  TreeTools::checkIds(tree, true);
 
163
  tree_ = new TreeTemplate<Node>(tree);
 
164
  root1_ = tree_->getRootNode()->getSon(0)->getId();
 
165
  root2_ = tree_->getRootNode()->getSon(1)->getId();
 
166
  nodes_ = tree_->getNodes();
 
167
  nodes_.pop_back(); //Remove the root node (the last added!).  
 
168
  nbNodes_ = nodes_.size();
 
169
  //Build nodes index:
 
170
  for (unsigned int i = 0; i < nodes_.size(); i++)
 
171
  {
 
172
    const Node * node = nodes_[i];
 
173
    idToNode_[node->getId()] = node;
 
174
  }
 
175
  nbClasses_ = rateDistribution_->getNumberOfCategories();
 
176
 
 
177
  verbose_ = verbose;
 
178
 
 
179
  minimumBrLen_ = 0.000001;
 
180
  brLenConstraint_.reset(new IncludingPositiveReal(minimumBrLen_));
 
181
  setSubstitutionModelSet(modelSet);
 
182
}
 
183
 
 
184
/******************************************************************************/
 
185
 
 
186
void AbstractNonHomogeneousTreeLikelihood::setSubstitutionModelSet(SubstitutionModelSet* modelSet) throw (Exception)
 
187
{
 
188
  //Check:
 
189
  if (data_)
 
190
  {
 
191
    if (modelSet->getAlphabet()->getAlphabetType() != data_->getAlphabet()->getAlphabetType())
 
192
      throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::setSubstitutionModelSet(). Model alphabet do not match existing data.");
 
193
  }
 
194
 
 
195
  modelSet_ = modelSet;
 
196
  
 
197
  if (data_)
 
198
  {
 
199
    if (modelSet->getNumberOfStates() != modelSet_->getNumberOfStates())
 
200
      setData(*data_); //Have to reinitialize the whole data structure.
 
201
  }
 
202
  
 
203
  nbStates_ = modelSet->getNumberOfStates();
 
204
 
 
205
  //Allocate transition probabilities arrays:
 
206
  for (unsigned int l = 0; l < nbNodes_; l++)
 
207
  {
 
208
    //For each son node,i
 
209
    Node* son = nodes_[l];
 
210
 
 
211
    VVVdouble* pxy__son = & pxy_[son->getId()];
 
212
    pxy__son->resize(nbClasses_);
 
213
    for (unsigned int c = 0; c < nbClasses_; c++)
 
214
    {
 
215
      VVdouble * pxy__son_c = & (* pxy__son)[c];
 
216
      pxy__son_c->resize(nbStates_);
 
217
      for(unsigned int x = 0; x < nbStates_; x++)
 
218
      {
 
219
        (*pxy__son_c)[x].resize(nbStates_);
 
220
      }
 
221
    }
 
222
  
 
223
    VVVdouble* dpxy__son = & dpxy_[son->getId()];
 
224
    dpxy__son->resize(nbClasses_);
 
225
    for (unsigned int c = 0; c < nbClasses_; c++)
 
226
    {
 
227
      VVdouble * dpxy__son_c = & (* dpxy__son)[c];
 
228
      dpxy__son_c->resize(nbStates_);
 
229
      for(unsigned int x = 0; x < nbStates_; x++)
 
230
      {
 
231
        (* dpxy__son_c)[x].resize(nbStates_);
 
232
      }
 
233
    }
 
234
      
 
235
    VVVdouble* d2pxy__son = & d2pxy_[son->getId()];
 
236
    d2pxy__son->resize(nbClasses_);
 
237
    for (unsigned int c = 0; c < nbClasses_; c++)
 
238
    {
 
239
      VVdouble * d2pxy__son_c = & (* d2pxy__son)[c];
 
240
      d2pxy__son_c->resize(nbStates_);
 
241
      for(unsigned int x = 0; x < nbStates_; x++)
 
242
      {
 
243
        (* d2pxy__son_c)[x].resize(nbStates_);
 
244
      }
 
245
    }
 
246
  }
 
247
 
 
248
  //We have to reset parameters. If the instance is not initialized, this will be done by the initialize method.
 
249
  if (initialized_) 
 
250
  {
 
251
    initParameters();
 
252
    computeAllTransitionProbabilities();
 
253
    fireParameterChanged(getParameters());
 
254
  }
 
255
}
 
256
 
 
257
/******************************************************************************/
 
258
 
 
259
void AbstractNonHomogeneousTreeLikelihood::initialize() throw (Exception)
 
260
{
 
261
  if (initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::initialize(). Object is already initialized.");
 
262
  if (!data_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::initialize(). Data are no set.");
 
263
  initParameters();
 
264
  initialized_ = true;
 
265
  computeAllTransitionProbabilities();
 
266
  fireParameterChanged(getParameters());
 
267
}
 
268
 
 
269
/******************************************************************************/
 
270
 
 
271
ParameterList AbstractNonHomogeneousTreeLikelihood::getBranchLengthsParameters() const
 
272
{
 
273
  if(!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
 
274
  return brLenParameters_.getCommonParametersWith(getParameters());
 
275
}
 
276
 
 
277
/******************************************************************************/
 
278
 
 
279
ParameterList AbstractNonHomogeneousTreeLikelihood::getSubstitutionModelParameters() const
 
280
{
 
281
  if(!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::getSubstitutionModelParameters(). Object is not initialized.");
 
282
  return modelSet_->getParameters().getCommonParametersWith(getParameters());
 
283
}
 
284
 
 
285
/******************************************************************************/
 
286
 
 
287
void AbstractNonHomogeneousTreeLikelihood::initParameters()
 
288
{
 
289
  // Reset parameters:
 
290
  resetParameters_();
 
291
  
 
292
  // Branch lengths:
 
293
  initBranchLengthsParameters();
 
294
  addParameters_(brLenParameters_);
 
295
  
 
296
  // Substitution model:
 
297
  addParameters_(modelSet_->getIndependentParameters());
 
298
  
 
299
  // Rate distribution:
 
300
  addParameters_(rateDistribution_->getIndependentParameters());
 
301
}
 
302
 
 
303
/******************************************************************************/
 
304
 
 
305
void AbstractNonHomogeneousTreeLikelihood::applyParameters() throw (Exception)
 
306
{
 
307
  if(!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
 
308
  //Apply branch lengths:
 
309
  for(unsigned int i = 0; i < nbNodes_; i++)
 
310
  {
 
311
    int id = nodes_[i]->getId();
 
312
    if(id == root1_)
 
313
    {
 
314
      const Parameter * rootBrLen = &getParameter("BrLenRoot");
 
315
      const Parameter * rootPos = &getParameter("RootPosition");
 
316
      nodes_[i]->setDistanceToFather(rootBrLen->getValue() * rootPos->getValue());
 
317
    }
 
318
    else if(id == root2_)
 
319
    {
 
320
      const Parameter * rootBrLen = &getParameter("BrLenRoot");
 
321
      const Parameter * rootPos = &getParameter("RootPosition");
 
322
      nodes_[i]->setDistanceToFather(rootBrLen->getValue() * (1. - rootPos->getValue()));
 
323
    }
 
324
    else
 
325
    {
 
326
      const Parameter * brLen = &getParameter(string("BrLen") + TextTools::toString(i));
 
327
      if(brLen != NULL) nodes_[i]->setDistanceToFather(brLen->getValue());
 
328
    }
 
329
  }
 
330
  //Apply substitution model parameters:
 
331
  modelSet_->matchParametersValues(getParameters());
 
332
  //Apply rate distribution parameters:
 
333
  rateDistribution_->matchParametersValues(getParameters());
 
334
}
 
335
 
 
336
/******************************************************************************/
 
337
 
 
338
void AbstractNonHomogeneousTreeLikelihood::initBranchLengthsParameters()
 
339
{
 
340
  brLenParameters_.reset();
 
341
  double l1 = 0, l2 = 0;
 
342
  for (unsigned int i = 0; i < nbNodes_; i++)
 
343
  {
 
344
    double d = minimumBrLen_;
 
345
    if (!nodes_[i]->hasDistanceToFather())
 
346
    {
 
347
      ApplicationTools::displayWarning("Missing branch length " + TextTools::toString(i) + ". Value is set to " + TextTools::toString(minimumBrLen_));
 
348
      nodes_[i]->setDistanceToFather(minimumBrLen_);
 
349
    }
 
350
    else
 
351
    {
 
352
      d = nodes_[i]->getDistanceToFather();
 
353
      if (d < minimumBrLen_)
 
354
      {
 
355
        ApplicationTools::displayWarning("Branch length " + TextTools::toString(i) + " is too small: " + TextTools::toString(d) + ". Value is set to " + TextTools::toString(minimumBrLen_));
 
356
        nodes_[i]->setDistanceToFather(minimumBrLen_);
 
357
        d = minimumBrLen_;
 
358
      }
 
359
    }
 
360
    if (nodes_[i]->getId() == root1_)
 
361
      l1 = d;
 
362
    else if(nodes_[i]->getId() == root2_)
 
363
      l2 = d;
 
364
    else
 
365
    {
 
366
      brLenParameters_.addParameter(Parameter("BrLen" + TextTools::toString(i), d, brLenConstraint_->clone(), true)); //Attach constraint to avoid clonage problems!
 
367
    }
 
368
  }
 
369
  brLenParameters_.addParameter(Parameter("BrLenRoot", l1 + l2, brLenConstraint_->clone(), true));
 
370
  brLenParameters_.addParameter(Parameter("RootPosition", l1 / (l1 + l2), &Parameter::PROP_CONSTRAINT_EX));
 
371
}
 
372
 
 
373
/*******************************************************************************/
 
374
 
 
375
void AbstractNonHomogeneousTreeLikelihood::computeAllTransitionProbabilities()
 
376
{
 
377
  for(unsigned int l = 0; l < nbNodes_; l++)
 
378
  {
 
379
    //For each node,
 
380
    Node * node = nodes_[l];
 
381
    computeTransitionProbabilitiesForNode(node);
 
382
  }
 
383
  rootFreqs_ = modelSet_->getRootFrequencies();
 
384
}
 
385
 
 
386
/*******************************************************************************/
 
387
 
 
388
void AbstractNonHomogeneousTreeLikelihood::computeTransitionProbabilitiesForNode(const Node* node)
 
389
{
 
390
  const SubstitutionModel* model = modelSet_->getModelForNode(node->getId());
 
391
  double l = node->getDistanceToFather(); 
 
392
 
 
393
  //Computes all pxy and pyx once for all:
 
394
  VVVdouble * pxy__node = & pxy_[node->getId()];
 
395
  for(unsigned int c = 0; c < nbClasses_; c++)
 
396
  {
 
397
    VVdouble * pxy__node_c = & (* pxy__node)[c];
 
398
    RowMatrix<double> Q = model->getPij_t(l * rateDistribution_->getCategory(c));
 
399
    for(unsigned int x = 0; x < nbStates_; x++)
 
400
    {
 
401
      Vdouble * pxy__node_c_x = & (* pxy__node_c)[x];
 
402
      for(unsigned int y = 0; y < nbStates_; y++)
 
403
      {
 
404
        (* pxy__node_c_x)[y] = Q(x, y);
 
405
      }
 
406
    }
 
407
  }
 
408
  
 
409
  if(computeFirstOrderDerivatives_)
 
410
  {
 
411
    //Computes all dpxy/dt once for all:
 
412
    VVVdouble * dpxy__node = & dpxy_[node->getId()];
 
413
 
 
414
    for(unsigned int c = 0; c < nbClasses_; c++)
 
415
    {
 
416
      VVdouble * dpxy__node_c = & (* dpxy__node)[c];
 
417
      double rc = rateDistribution_->getCategory(c);
 
418
 
 
419
      RowMatrix<double> dQ = model->getdPij_dt(l * rc);  
 
420
 
 
421
      for(unsigned int x = 0; x < nbStates_; x++)
 
422
      {
 
423
        Vdouble * dpxy__node_c_x = & (* dpxy__node_c)[x];
 
424
        for(unsigned int y = 0; y < nbStates_; y++)
 
425
          (* dpxy__node_c_x)[y] = rc * dQ(x, y); 
 
426
      }
 
427
    }
 
428
  }
 
429
      
 
430
  if(computeSecondOrderDerivatives_)
 
431
  {
 
432
    //Computes all d2pxy/dt2 once for all:
 
433
    VVVdouble * d2pxy__node = & d2pxy_[node->getId()];
 
434
    for(unsigned int c = 0; c < nbClasses_; c++)
 
435
    {
 
436
      VVdouble * d2pxy__node_c = & (* d2pxy__node)[c];
 
437
      double rc =  rateDistribution_->getCategory(c);
 
438
      RowMatrix<double> d2Q = model->getd2Pij_dt2(l * rc);
 
439
      for(unsigned int x = 0; x < nbStates_; x++)
 
440
      {
 
441
        Vdouble * d2pxy__node_c_x = & (* d2pxy__node_c)[x];
 
442
        for(unsigned int y = 0; y < nbStates_; y++)
 
443
        {
 
444
          (* d2pxy__node_c_x)[y] = rc * rc * d2Q(x, y);
 
445
        }
 
446
      }
 
447
    }
 
448
  }
 
449
}
 
450
 
 
451
/*******************************************************************************/
 
452