bpp-phyl  2.4.0
AbstractHomogeneousTreeLikelihood.cpp
Go to the documentation of this file.
1 //
2 // File: AbstractHomogeneousTreeLikelihood.cpp
3 // Created by: Julien Dutheil
4 // Created on: Thr Dec 23 12:03 2004
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 #include "../PatternTools.h"
42 
43 #include <Bpp/Text/TextTools.h>
45 
46 // From bpp-seq:
47 #include <Bpp/Seq/SiteTools.h>
49 
50 using namespace bpp;
51 
52 // From the STL:
53 #include <iostream>
54 
55 using namespace std;
56 
57 /******************************************************************************/
58 
60  const Tree& tree,
61  TransitionModel* model,
62  DiscreteDistribution* rDist,
63  bool checkRooted,
64  bool verbose) :
66  model_(0),
67  brLenParameters_(),
68  pxy_(),
69  dpxy_(),
70  d2pxy_(),
71  rootFreqs_(),
72  nodes_(),
73  nbSites_(),
74  nbDistinctSites_(),
75  nbClasses_(),
76  nbStates_(),
77  nbNodes_(),
78  verbose_(),
79  minimumBrLen_(),
80  maximumBrLen_(),
81  brLenConstraint_()
82 {
83  init_(tree, model, rDist, checkRooted, verbose);
84 }
85 
86 /******************************************************************************/
87 
91  model_(lik.model_),
93  pxy_(lik.pxy_),
94  dpxy_(lik.dpxy_),
95  d2pxy_(lik.d2pxy_),
97  nodes_(),
98  nbSites_(lik.nbSites_),
100  nbClasses_(lik.nbClasses_),
101  nbStates_(lik.nbStates_),
102  nbNodes_(lik.nbNodes_),
103  verbose_(lik.verbose_),
107 {
108  nodes_ = tree_->getNodes();
109  nodes_.pop_back(); // Remove the root node (the last added!).
110 }
111 
112 /******************************************************************************/
113 
116 {
118  model_ = lik.model_;
120  pxy_ = lik.pxy_;
121  dpxy_ = lik.dpxy_;
122  d2pxy_ = lik.d2pxy_;
123  rootFreqs_ = lik.rootFreqs_;
124  nodes_ = tree_->getNodes();
125  nodes_.pop_back(); // Remove the root node (the last added!).
126  nbSites_ = lik.nbSites_;
128  nbClasses_ = lik.nbClasses_;
129  nbStates_ = lik.nbStates_;
130  nbNodes_ = lik.nbNodes_;
131  verbose_ = lik.verbose_;
134  if (brLenConstraint_.get())
135  brLenConstraint_.release();
136  brLenConstraint_.reset(lik.brLenConstraint_->clone());
137  return *this;
138 }
139 
140 /******************************************************************************/
141 
143  const Tree& tree,
144  TransitionModel* model,
145  DiscreteDistribution* rDist,
146  bool checkRooted,
147  bool verbose)
148 {
149  TreeTools::checkIds(tree, true);
150  tree_ = new TreeTemplate<Node>(tree);
151  if (checkRooted && tree_->isRooted())
152  {
153  if (verbose)
154  ApplicationTools::displayWarning("Tree has been unrooted.");
155  tree_->unroot();
156  }
157  nodes_ = tree_->getNodes();
158  nodes_.pop_back(); // Remove the root node (the last added!).
159  nbNodes_ = nodes_.size();
161  setModel(model);
162 
163  verbose_ = verbose;
164 
165  minimumBrLen_ = 0.000001;
166  maximumBrLen_ = 10000;
168 }
169 
170 /******************************************************************************/
171 
173 {
174  // Check:
175  if (data_)
176  {
178  throw Exception("AbstractHomogeneousTreeLikelihood::setSubstitutionModel(). Model alphabet do not match existing data.");
179  }
180 
181  model_ = model;
182 
183  if (data_)
184  {
185  if (model->getNumberOfStates() != model_->getNumberOfStates())
186  setData(*data_); // Have to reinitialize the whole data structure.
187  }
188 
189  nbStates_ = model->getNumberOfStates();
190 
191  // Allocate transition probabilities arrays:
192  for (unsigned int l = 0; l < nbNodes_; l++)
193  {
194  // For each son node,i
195  Node* son = nodes_[l];
196 
197  VVVdouble* pxy__son = &pxy_[son->getId()];
198  pxy__son->resize(nbClasses_);
199  for (unsigned int c = 0; c < nbClasses_; c++)
200  {
201  VVdouble* pxy__son_c = &(*pxy__son)[c];
202  pxy__son_c->resize(nbStates_);
203  for (unsigned int x = 0; x < nbStates_; x++)
204  {
205  (*pxy__son_c)[x].resize(nbStates_);
206  }
207  }
208 
209  VVVdouble* dpxy__son = &dpxy_[son->getId()];
210  dpxy__son->resize(nbClasses_);
211  for (unsigned int c = 0; c < nbClasses_; c++)
212  {
213  VVdouble* dpxy__son_c = &(*dpxy__son)[c];
214  dpxy__son_c->resize(nbStates_);
215  for (unsigned int x = 0; x < nbStates_; x++)
216  {
217  (*dpxy__son_c)[x].resize(nbStates_);
218  }
219  }
220 
221  VVVdouble* d2pxy__son = &d2pxy_[son->getId()];
222  d2pxy__son->resize(nbClasses_);
223  for (unsigned int c = 0; c < nbClasses_; c++)
224  {
225  VVdouble* d2pxy__son_c = &(*d2pxy__son)[c];
226  d2pxy__son_c->resize(nbStates_);
227  for (unsigned int x = 0; x < nbStates_; x++)
228  {
229  (*d2pxy__son_c)[x].resize(nbStates_);
230  }
231  }
232  }
233 }
234 
235 /******************************************************************************/
236 
238 {
239  if (initialized_)
240  throw Exception("AbstractHomogeneousTreeLikelihood::initialize(). Object is already initialized.");
241  if (!data_)
242  throw Exception("AbstractHomogeneousTreeLikelihood::initialize(). Data are no set.");
243  initParameters();
244  initialized_ = true;
246 }
247 
248 /******************************************************************************/
249 
251 {
252  if (!initialized_)
253  throw Exception("AbstractHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
255 }
256 
257 /******************************************************************************/
258 
260 {
261  if (!initialized_)
262  throw Exception("AbstractHomogeneousTreeLikelihood::getSubstitutionModelParameters(). Object is not initialized.");
264 }
265 
266 /******************************************************************************/
267 
269 {
270  // Reset parameters:
272 
273  // Branch lengths:
276 
277  // Substitution model:
279 
280  // Rate distribution:
282 }
283 
284 /******************************************************************************/
285 
287 {
288  if (!initialized_)
289  throw Exception("AbstractHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
290  // Apply branch lengths:
291  // brLenParameters_.matchParametersValues(parameters_); Not necessary!
292  for (unsigned int i = 0; i < nbNodes_; i++)
293  {
294  const Parameter* brLen = &getParameter(string("BrLen") + TextTools::toString(i));
295  if (brLen)
296  nodes_[i]->setDistanceToFather(brLen->getValue());
297  }
298  // Apply substitution model parameters:
301  // Apply rate distribution parameters:
303 }
304 
305 /******************************************************************************/
306 
308 {
310  for (unsigned int i = 0; i < nbNodes_; i++)
311  {
312  double d = minimumBrLen_;
313  if (!nodes_[i]->hasDistanceToFather())
314  {
315  if (verbose)
316  ApplicationTools::displayWarning("Missing branch length " + TextTools::toString(i) + ". Value is set to " + TextTools::toString(minimumBrLen_));
317  nodes_[i]->setDistanceToFather(minimumBrLen_);
318  }
319  else
320  {
321  d = nodes_[i]->getDistanceToFather();
322  if (d < minimumBrLen_)
323  {
324  if (verbose)
325  ApplicationTools::displayWarning("Branch length " + TextTools::toString(i) + " is too small: " + TextTools::toString(d) + ". Value is set to " + TextTools::toString(minimumBrLen_));
326  nodes_[i]->setDistanceToFather(minimumBrLen_);
327  d = minimumBrLen_;
328  }
329  if (d > maximumBrLen_)
330  {
331  if (verbose)
332  ApplicationTools::displayWarning("Branch length " + TextTools::toString(i) + " is too big: " + TextTools::toString(d) + ". Value is set to " + TextTools::toString(maximumBrLen_));
333  nodes_[i]->setDistanceToFather(maximumBrLen_);
334  d = maximumBrLen_;
335  }
336  }
337  brLenParameters_.addParameter(Parameter("BrLen" + TextTools::toString(i), d, brLenConstraint_->clone(), true)); // Attach constraint to avoid clonage problems!
338  }
339 }
340 
341 /*******************************************************************************/
342 
344 {
345  for (unsigned int l = 0; l < nbNodes_; l++)
346  {
347  // For each node,
348  Node* node = nodes_[l];
350  }
352 }
353 
354 /*******************************************************************************/
355 
357 {
358  double l = node->getDistanceToFather();
359 
360  // Computes all pxy and pyx once for all:
361  VVVdouble* pxy__node = &pxy_[node->getId()];
362  for (unsigned int c = 0; c < nbClasses_; c++)
363  {
364  VVdouble* pxy__node_c = &(*pxy__node)[c];
366 
367  for (unsigned int x = 0; x < nbStates_; x++)
368  {
369  Vdouble* pxy__node_c_x = &(*pxy__node_c)[x];
370  for (unsigned int y = 0; y < nbStates_; y++)
371  {
372  (*pxy__node_c_x)[y] = Q(x, y);
373  }
374  }
375  }
376 
378  {
379  // Computes all dpxy/dt once for all:
380  VVVdouble* dpxy__node = &dpxy_[node->getId()];
381  for (unsigned int c = 0; c < nbClasses_; c++)
382  {
383  VVdouble* dpxy__node_c = &(*dpxy__node)[c];
384  double rc = rateDistribution_->getCategory(c);
385  RowMatrix<double> dQ = model_->getdPij_dt(l * rc);
386  for (unsigned int x = 0; x < nbStates_; x++)
387  {
388  Vdouble* dpxy__node_c_x = &(*dpxy__node_c)[x];
389  for (unsigned int y = 0; y < nbStates_; y++)
390  {
391  (*dpxy__node_c_x)[y] = rc * dQ(x, y);
392  }
393  }
394  }
395  }
396 
398  {
399  // Computes all d2pxy/dt2 once for all:
400  VVVdouble* d2pxy__node = &d2pxy_[node->getId()];
401  for (unsigned int c = 0; c < nbClasses_; c++)
402  {
403  VVdouble* d2pxy__node_c = &(*d2pxy__node)[c];
404  double rc = rateDistribution_->getCategory(c);
405  RowMatrix<double> d2Q = model_->getd2Pij_dt2(l * rc);
406  for (unsigned int x = 0; x < nbStates_; x++)
407  {
408  Vdouble* d2pxy__node_c_x = &(*d2pxy__node_c)[x];
409  for (unsigned int y = 0; y < nbStates_; y++)
410  {
411  (*d2pxy__node_c_x)[y] = rc * rc * d2Q(x, y);
412  }
413  }
414  }
415  }
416 }
417 
418 /*******************************************************************************/
virtual bool matchParametersValues(const ParameterList &parameters)=0
virtual void setData(const SiteContainer &sites)=0
Set the dataset for which the likelihood must be evaluated.
static void displayWarning(const std::string &text)
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
virtual double getDistanceToFather() const
Get the distance to the father node is there is one, otherwise throw a NodeException.
Definition: Node.h:283
virtual size_t getNumberOfCategories() const =0
void init_(const Tree &tree, TransitionModel *model, DiscreteDistribution *rDist, bool checkRooted, bool verbose)
Method called by constructor.
Partial implementation of the DiscreteRatesAcrossSitesTreeLikelihood interface.
virtual void reset()
AbstractHomogeneousTreeLikelihood & operator=(const AbstractHomogeneousTreeLikelihood &lik)
Assignation operator.
virtual const Matrix< double > & getd2Pij_dt2(double t) const =0
virtual double getValue() const
virtual const Vdouble & getFrequencies() const =0
virtual void initParameters()
This builds the parameters list from all parametrizable objects, i.e. substitution model...
ParameterList getSubstitutionModelParameters() const
Get the parameters associated to substitution model(s).
virtual const ParameterList & getParameters() const =0
STL namespace.
The phylogenetic tree class.
Interface for phylogenetic tree objects.
Definition: Tree.h:148
virtual int getId() const
Get the node&#39;s id.
Definition: Node.h:203
virtual const Matrix< double > & getPij_t(double t) const =0
virtual void fireParameterChanged(const ParameterList &parameters)
virtual const Matrix< double > & getdPij_dt(double t) const =0
ParameterList getBranchLengthsParameters() const
Get the branch lengths parameters.
HomogeneousTreeLikelihood * clone() const =0
static bool checkIds(const Tree &tree, bool throwException)
Check if the ids are uniques.
Definition: TreeTools.cpp:778
virtual double getCategory(size_t categoryIndex) const =0
std::vector< double > Vdouble
virtual const ParameterList & getIndependentParameters() const =0
virtual const Alphabet * getAlphabet() const =0
AbstractDiscreteRatesAcrossSitesTreeLikelihood & operator=(const AbstractDiscreteRatesAcrossSitesTreeLikelihood &tl)
The phylogenetic node class.
Definition: Node.h:90
virtual ParameterList getCommonParametersWith(const ParameterList &params) const
AbstractHomogeneousTreeLikelihood(const Tree &tree, TransitionModel *model, DiscreteDistribution *rDist, bool checkRooted=true, bool verbose=true)
virtual void addParameter(const Parameter &param)
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
const ParameterList & getParameters() const
const Parameter & getParameter(const std::string &name) const
virtual const Alphabet * getAlphabet() const =0
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
std::string toString(T t)
virtual size_t getNumberOfStates() const =0
Get the number of states.
Partial implementation for homogeneous model of the TreeLikelihood interface.
std::vector< VVdouble > VVVdouble
std::vector< Vdouble > VVdouble
virtual void addParameters_(const ParameterList &parameters)
virtual std::string getAlphabetType() const =0