bpp-phyl3  3.0.0
RNonHomogeneousTreeLikelihood.cpp
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: The Bio++ Development Group
2 //
3 // SPDX-License-Identifier: CECILL-2.1
4 
5 #include <Bpp/App/ApplicationTools.h>
6 #include <Bpp/Text/TextTools.h>
7 
8 #include "../../PatternTools.h"
9 #include "RNonHomogeneousTreeLikelihood.h"
10 
11 using namespace bpp;
12 
13 // From the STL:
14 #include <iostream>
15 
16 using namespace std;
17 
18 /******************************************************************************/
19 
20 RNonHomogeneousTreeLikelihood::RNonHomogeneousTreeLikelihood(
21  const Tree& tree,
22  std::shared_ptr<SubstitutionModelSet> modelSet,
23  std::shared_ptr<DiscreteDistributionInterface> rDist,
24  bool verbose,
25  bool usePatterns,
26  bool reparametrizeRoot) :
27  AbstractNonHomogeneousTreeLikelihood(tree, modelSet, rDist, verbose, reparametrizeRoot),
28  likelihoodData_(),
29  minusLogLik_(-1.)
30 {
31  if (!modelSet->isFullySetUpFor(tree))
32  throw Exception("RNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
33  init_(usePatterns);
34 }
35 
36 /******************************************************************************/
37 
38 RNonHomogeneousTreeLikelihood::RNonHomogeneousTreeLikelihood(
39  const Tree& tree,
40  const AlignmentDataInterface& data,
41  std::shared_ptr<SubstitutionModelSet> modelSet,
42  std::shared_ptr<DiscreteDistributionInterface> rDist,
43  bool verbose,
44  bool usePatterns,
45  bool reparametrizeRoot) :
46  AbstractNonHomogeneousTreeLikelihood(tree, modelSet, rDist, verbose, reparametrizeRoot),
47  likelihoodData_(),
48  minusLogLik_(-1.)
49 {
50  if (!modelSet->isFullySetUpFor(tree))
51  throw Exception("RNonHomogeneousTreeLikelihood(constructor). Model set is not fully specified.");
52  init_(usePatterns);
53  setData(data);
54 }
55 
56 /******************************************************************************/
57 
58 void RNonHomogeneousTreeLikelihood::init_(bool usePatterns)
59 {
60  likelihoodData_ = make_shared<DRASRTreeLikelihoodData>(
61  tree_,
62  rateDistribution_->getNumberOfCategories(),
63  usePatterns);
64 }
65 
66 /******************************************************************************/
67 
68 RNonHomogeneousTreeLikelihood::RNonHomogeneousTreeLikelihood(
69  const RNonHomogeneousTreeLikelihood& lik) :
70  AbstractNonHomogeneousTreeLikelihood(lik),
71  likelihoodData_(),
72  minusLogLik_(lik.minusLogLik_)
73 {
74  likelihoodData_ = shared_ptr<DRASRTreeLikelihoodData>(lik.likelihoodData_->clone());
75  likelihoodData_->setTree(tree_);
76 }
77 
78 /******************************************************************************/
79 
80 RNonHomogeneousTreeLikelihood& RNonHomogeneousTreeLikelihood::operator=(
81  const RNonHomogeneousTreeLikelihood& lik)
82 {
83  AbstractNonHomogeneousTreeLikelihood::operator=(lik);
84  likelihoodData_ = shared_ptr<DRASRTreeLikelihoodData>(lik.likelihoodData_->clone());
85  likelihoodData_->setTree(tree_);
86  minusLogLik_ = lik.minusLogLik_;
87  return *this;
88 }
89 
90 /******************************************************************************/
91 
92 RNonHomogeneousTreeLikelihood::~RNonHomogeneousTreeLikelihood() {}
93 
94 /******************************************************************************/
95 
96 void RNonHomogeneousTreeLikelihood::setData(const AlignmentDataInterface& sites)
97 {
98  data_ = PatternTools::getSequenceSubset(sites, *tree_->getRootNode());
99  if (verbose_)
100  ApplicationTools::displayTask("Initializing data structure");
101  likelihoodData_->initLikelihoods(*data_, *modelSet_->getModel(0)); // We assume here that all models have the same number of states, and that they have the same 'init' method,
102  // Which is a reasonable assumption as long as they share the same alphabet.
103  if (verbose_)
104  ApplicationTools::displayTaskDone();
105 
106  nbSites_ = likelihoodData_->getNumberOfSites();
107  nbDistinctSites_ = likelihoodData_->getNumberOfDistinctSites();
108  nbStates_ = likelihoodData_->getNumberOfStates();
109 
110  if (verbose_)
111  ApplicationTools::displayResult("Number of distinct sites",
112  TextTools::toString(nbDistinctSites_));
113  initialized_ = false;
114 }
115 
116 /******************************************************************************/
117 
118 double RNonHomogeneousTreeLikelihood::getLikelihood() const
119 {
120  double l = 1.;
121  for (size_t i = 0; i < nbSites_; i++)
122  {
123  l *= getLikelihoodForASite(i);
124  }
125  return l;
126 }
127 
128 /******************************************************************************/
129 
130 double RNonHomogeneousTreeLikelihood::getLogLikelihood() const
131 {
132  double ll = 0;
133  vector<double> la(nbSites_);
134  for (size_t i = 0; i < nbSites_; i++)
135  {
136  la[i] = getLogLikelihoodForASite(i);
137  }
138  sort(la.begin(), la.end());
139  for (size_t i = nbSites_; i > 0; i--)
140  {
141  ll += la[i - 1];
142  }
143  return ll;
144 }
145 
146 /******************************************************************************/
147 
148 double RNonHomogeneousTreeLikelihood::getLikelihoodForASite(size_t site) const
149 {
150  double l = 0;
151  for (size_t i = 0; i < nbClasses_; i++)
152  {
153  l += getLikelihoodForASiteForARateClass(site, i) * rateDistribution_->getProbability(i);
154  }
155  return l;
156 }
157 
158 /******************************************************************************/
159 
160 double RNonHomogeneousTreeLikelihood::getLogLikelihoodForASite(size_t site) const
161 {
162  double l = 0;
163  for (size_t i = 0; i < nbClasses_; i++)
164  {
165  l += getLikelihoodForASiteForARateClass(site, i) * rateDistribution_->getProbability(i);
166  }
167  // if(l <= 0.) cerr << "WARNING!!! Negative likelihood." << endl;
168  if (l < 0)
169  l = 0; // May happen because of numerical errors.
170  return log(l);
171 }
172 
173 /******************************************************************************/
174 
175 double RNonHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
176 {
177  double l = 0;
178  Vdouble* la = &likelihoodData_->getLikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass];
179  for (size_t i = 0; i < nbStates_; i++)
180  {
181  l += (*la)[i] * rootFreqs_[i];
182  }
183  return l;
184 }
185 
186 /******************************************************************************/
187 
188 double RNonHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
189 {
190  double l = 0;
191  Vdouble* la = &likelihoodData_->getLikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass];
192  for (size_t i = 0; i < nbStates_; i++)
193  {
194  l += (*la)[i] * rootFreqs_[i];
195  }
196  return log(l);
197 }
198 
199 /******************************************************************************/
200 
201 double RNonHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
202 {
203  return likelihoodData_->getLikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass][static_cast<size_t>(state)];
204 }
205 
206 /******************************************************************************/
207 
208 double RNonHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
209 {
210  return log(likelihoodData_->getLikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass][static_cast<size_t>(state)]);
211 }
212 
213 /******************************************************************************/
214 
215 void RNonHomogeneousTreeLikelihood::setParameters(const ParameterList& parameters)
216 {
217  setParametersValues(parameters);
218 }
219 
220 /******************************************************************************/
221 
222 void RNonHomogeneousTreeLikelihood::fireParameterChanged(const ParameterList& params)
223 {
224  applyParameters();
225 
226  if (params.getCommonParametersWith(rateDistribution_->getIndependentParameters()).size() > 0)
227  {
228  computeAllTransitionProbabilities();
229  }
230  else
231  {
232  vector<int> ids;
233  vector<string> tmp = params.getCommonParametersWith(modelSet_->getNodeParameters()).getParameterNames();
234  for (size_t i = 0; i < tmp.size(); i++)
235  {
236  vector<int> tmpv = modelSet_->getNodesWithParameter(tmp[i]);
237  ids = VectorTools::vectorUnion(ids, tmpv);
238  }
239  tmp = params.getCommonParametersWith(brLenParameters_).getParameterNames();
240  vector<const Node*> nodes;
241  for (size_t i = 0; i < ids.size(); i++)
242  {
243  nodes.push_back(idToNode_[ids[i]]);
244  }
245  vector<const Node*> tmpv;
246  bool test = false;
247  for (size_t i = 0; i < tmp.size(); i++)
248  {
249  if (tmp[i] == "BrLenRoot" || tmp[i] == "RootPosition")
250  {
251  if (!test)
252  {
253  tmpv.push_back(tree_->getRootNode()->getSon(0));
254  tmpv.push_back(tree_->getRootNode()->getSon(1));
255  test = true; // Add only once.
256  }
257  }
258  else
259  tmpv.push_back(nodes_[TextTools::to< size_t >(tmp[i].substr(5))]);
260  }
261  nodes = VectorTools::vectorUnion(nodes, tmpv);
262 
263  for (size_t i = 0; i < nodes.size(); i++)
264  {
265  computeTransitionProbabilitiesForNode(nodes[i]);
266  }
267  rootFreqs_ = modelSet_->getRootFrequencies();
268  }
269  computeTreeLikelihood();
270 
271  minusLogLik_ = -getLogLikelihood();
272 }
273 
274 /******************************************************************************/
275 
276 double RNonHomogeneousTreeLikelihood::getValue() const
277 {
278  if (!isInitialized())
279  throw Exception("RNonHomogeneousTreeLikelihood::getValue(). Instance is not initialized.");
280  return minusLogLik_;
281 }
282 
283 /******************************************************************************
284 * First Order Derivatives *
285 ******************************************************************************/
286 double RNonHomogeneousTreeLikelihood::getDLikelihoodForASiteForARateClass(
287  size_t site,
288  size_t rateClass) const
289 {
290  double dl = 0;
291  Vdouble* dla = &likelihoodData_->getDLikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass];
292  for (size_t i = 0; i < nbStates_; i++)
293  {
294  dl += (*dla)[i] * rootFreqs_[i];
295  }
296  return dl;
297 }
298 
299 /******************************************************************************/
300 
301 double RNonHomogeneousTreeLikelihood::getDLikelihoodForASite(size_t site) const
302 {
303  // Derivative of the sum is the sum of derivatives:
304  double dl = 0;
305  for (size_t i = 0; i < nbClasses_; i++)
306  {
307  dl += getDLikelihoodForASiteForARateClass(site, i) * rateDistribution_->getProbability(i);
308  }
309  return dl;
310 }
311 
312 /******************************************************************************/
313 
314 double RNonHomogeneousTreeLikelihood::getDLogLikelihoodForASite(size_t site) const
315 {
316  // d(f(g(x)))/dx = dg(x)/dx . df(g(x))/dg :
317  return getDLikelihoodForASite(site) / getLikelihoodForASite(site);
318 }
319 
320 /******************************************************************************/
321 
322 double RNonHomogeneousTreeLikelihood::getDLogLikelihood() const
323 {
324  // Derivative of the sum is the sum of derivatives:
325  double dl = 0;
326  for (size_t i = 0; i < nbSites_; i++)
327  {
328  dl += getDLogLikelihoodForASite(i);
329  }
330  return dl;
331 }
332 
333 /******************************************************************************/
334 
335 double RNonHomogeneousTreeLikelihood::getFirstOrderDerivative(const string& variable) const
336 {
337  if (!hasParameter(variable))
338  throw ParameterNotFoundException("RNonHomogeneousTreeLikelihood::getFirstOrderDerivative().", variable);
339  if (getRateDistributionParameters().hasParameter(variable))
340  {
341  throw Exception("Derivatives respective to rate distribution parameter are not implemented.");
342  }
343  if (getSubstitutionModelParameters().hasParameter(variable))
344  {
345  throw Exception("Derivatives respective to substitution model parameters are not implemented.");
346  }
347 
348  const_cast<RNonHomogeneousTreeLikelihood*>(this)->computeTreeDLikelihood(variable);
349  return -getDLogLikelihood();
350 }
351 
352 /******************************************************************************/
353 
354 void RNonHomogeneousTreeLikelihood::computeTreeDLikelihood(const string& variable)
355 {
356  if (variable == "BrLenRoot")
357  {
358  const Node* father = tree_->getRootNode();
359 
360  // Compute dLikelihoods array for the father node.
361  // Fist initialize to 1:
362  VVVdouble* _dLikelihoods_father = &likelihoodData_->getDLikelihoodArray(father->getId());
363  size_t nbSites = _dLikelihoods_father->size();
364  for (size_t i = 0; i < nbSites; i++)
365  {
366  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
367  for (size_t c = 0; c < nbClasses_; c++)
368  {
369  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
370  for (size_t s = 0; s < nbStates_; s++)
371  {
372  (*_dLikelihoods_father_i_c)[s] = 1.;
373  }
374  }
375  }
376 
377  size_t nbNodes = father->getNumberOfSons();
378  for (size_t l = 0; l < nbNodes; l++)
379  {
380  const Node* son = father->getSon(l);
381 
382  if (son->getId() == root1_)
383  {
384  const Node* root1 = father->getSon(0);
385  const Node* root2 = father->getSon(1);
386  vector<size_t>* _patternLinks_fatherroot1_ = &likelihoodData_->getArrayPositions(father->getId(), root1->getId());
387  vector<size_t>* _patternLinks_fatherroot2_ = &likelihoodData_->getArrayPositions(father->getId(), root2->getId());
388  VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(root1->getId());
389  VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(root2->getId());
390  double pos = getParameterValue("RootPosition");
391 
392  VVVdouble* dpxy_root1_ = &dpxy_[root1_];
393  VVVdouble* dpxy_root2_ = &dpxy_[root2_];
394  VVVdouble* pxy_root1_ = &pxy_[root1_];
395  VVVdouble* pxy_root2_ = &pxy_[root2_];
396  for (size_t i = 0; i < nbSites; i++)
397  {
398  VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
399  VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
400  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
401  for (size_t c = 0; c < nbClasses_; c++)
402  {
403  Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
404  Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
405  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
406  VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
407  VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
408  VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
409  VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
410  for (size_t x = 0; x < nbStates_; x++)
411  {
412  Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
413  Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
414  Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
415  Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
416  double dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
417  for (size_t y = 0; y < nbStates_; y++)
418  {
419  dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
420  dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
421  l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
422  l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
423  }
424  double dl = pos * dl1 * l2 + (1. - pos) * dl2 * l1;
425  (*_dLikelihoods_father_i_c)[x] *= dl;
426  }
427  }
428  }
429  }
430  else if (son->getId() == root2_)
431  {
432  // Do nothing, this was accounted when dealing with root1_
433  }
434  else
435  {
436  // Account for a putative multifurcation:
437  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
438  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
439 
440  VVVdouble* pxy__son = &pxy_[son->getId()];
441  for (size_t i = 0; i < nbSites; i++)
442  {
443  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
444  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
445  for (size_t c = 0; c < nbClasses_; c++)
446  {
447  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
448  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
449  VVdouble* pxy__son_c = &(*pxy__son)[c];
450  for (size_t x = 0; x < nbStates_; x++)
451  {
452  double dl = 0;
453  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
454  for (size_t y = 0; y < nbStates_; y++)
455  {
456  dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
457  }
458  (*_dLikelihoods_father_i_c)[x] *= dl;
459  }
460  }
461  }
462  }
463  }
464  return;
465  }
466  else if (variable == "RootPosition")
467  {
468  const Node* father = tree_->getRootNode();
469 
470  // Compute dLikelihoods array for the father node.
471  // Fist initialize to 1:
472  VVVdouble* _dLikelihoods_father = &likelihoodData_->getDLikelihoodArray(father->getId());
473  size_t nbSites = _dLikelihoods_father->size();
474  for (size_t i = 0; i < nbSites; i++)
475  {
476  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
477  for (size_t c = 0; c < nbClasses_; c++)
478  {
479  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
480  for (size_t s = 0; s < nbStates_; s++)
481  {
482  (*_dLikelihoods_father_i_c)[s] = 1.;
483  }
484  }
485  }
486 
487  size_t nbNodes = father->getNumberOfSons();
488  for (size_t l = 0; l < nbNodes; l++)
489  {
490  const Node* son = father->getSon(l);
491 
492  if (son->getId() == root1_)
493  {
494  const Node* root1 = father->getSon(0);
495  const Node* root2 = father->getSon(1);
496  vector<size_t>* _patternLinks_fatherroot1_ = &likelihoodData_->getArrayPositions(father->getId(), root1->getId());
497  vector<size_t>* _patternLinks_fatherroot2_ = &likelihoodData_->getArrayPositions(father->getId(), root2->getId());
498  VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(root1->getId());
499  VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(root2->getId());
500  double len = getParameterValue("BrLenRoot");
501 
502  VVVdouble* dpxy_root1_ = &dpxy_[root1_];
503  VVVdouble* dpxy_root2_ = &dpxy_[root2_];
504  VVVdouble* pxy_root1_ = &pxy_[root1_];
505  VVVdouble* pxy_root2_ = &pxy_[root2_];
506  for (size_t i = 0; i < nbSites; i++)
507  {
508  VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
509  VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
510  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
511  for (size_t c = 0; c < nbClasses_; c++)
512  {
513  Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
514  Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
515  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
516  VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
517  VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
518  VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
519  VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
520  for (size_t x = 0; x < nbStates_; x++)
521  {
522  Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
523  Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
524  Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
525  Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
526  double dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
527  for (size_t y = 0; y < nbStates_; y++)
528  {
529  dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
530  dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
531  l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
532  l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
533  }
534  double dl = len * (dl1 * l2 - dl2 * l1);
535  (*_dLikelihoods_father_i_c)[x] *= dl;
536  }
537  }
538  }
539  }
540  else if (son->getId() == root2_)
541  {
542  // Do nothing, this was accounted when dealing with root1_
543  }
544  else
545  {
546  // Account for a putative multifurcation:
547  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
548  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
549 
550  VVVdouble* pxy__son = &pxy_[son->getId()];
551  for (size_t i = 0; i < nbSites; i++)
552  {
553  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
554  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
555  for (size_t c = 0; c < nbClasses_; c++)
556  {
557  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
558  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
559  VVdouble* pxy__son_c = &(*pxy__son)[c];
560  for (size_t x = 0; x < nbStates_; x++)
561  {
562  double dl = 0;
563  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
564  for (size_t y = 0; y < nbStates_; y++)
565  {
566  dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
567  }
568  (*_dLikelihoods_father_i_c)[x] *= dl;
569  }
570  }
571  }
572  }
573  }
574  return;
575  }
576 
577  // Get the node with the branch whose length must be derivated:
578  size_t brI = TextTools::to<size_t>(variable.substr(5));
579  const Node* branch = nodes_[brI];
580  const Node* father = branch->getFather();
581  VVVdouble* _dLikelihoods_father = &likelihoodData_->getDLikelihoodArray(father->getId());
582 
583  // Compute dLikelihoods array for the father node.
584  // Fist initialize to 1:
585  size_t nbSites = _dLikelihoods_father->size();
586  for (size_t i = 0; i < nbSites; i++)
587  {
588  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
589  for (size_t c = 0; c < nbClasses_; c++)
590  {
591  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
592  for (size_t s = 0; s < nbStates_; s++)
593  {
594  (*_dLikelihoods_father_i_c)[s] = 1.;
595  }
596  }
597  }
598 
599  size_t nbNodes = father->getNumberOfSons();
600  for (size_t l = 0; l < nbNodes; l++)
601  {
602  const Node* son = father->getSon(l);
603 
604  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
605  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
606 
607  if (son == branch)
608  {
609  VVVdouble* dpxy__son = &dpxy_[son->getId()];
610  for (size_t i = 0; i < nbSites; i++)
611  {
612  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
613  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
614  for (size_t c = 0; c < nbClasses_; c++)
615  {
616  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
617  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
618  VVdouble* dpxy__son_c = &(*dpxy__son)[c];
619  for (size_t x = 0; x < nbStates_; x++)
620  {
621  double dl = 0;
622  Vdouble* dpxy__son_c_x = &(*dpxy__son_c)[x];
623  for (size_t y = 0; y < nbStates_; y++)
624  {
625  dl += (*dpxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
626  }
627  (*_dLikelihoods_father_i_c)[x] *= dl;
628  }
629  }
630  }
631  }
632  else
633  {
634  VVVdouble* pxy__son = &pxy_[son->getId()];
635  for (size_t i = 0; i < nbSites; i++)
636  {
637  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
638  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
639  for (size_t c = 0; c < nbClasses_; c++)
640  {
641  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
642  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
643  VVdouble* pxy__son_c = &(*pxy__son)[c];
644  for (size_t x = 0; x < nbStates_; x++)
645  {
646  double dl = 0;
647  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
648  for (size_t y = 0; y < nbStates_; y++)
649  {
650  dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
651  }
652  (*_dLikelihoods_father_i_c)[x] *= dl;
653  }
654  }
655  }
656  }
657  }
658 
659  // Now we go down the tree toward the root node:
660  computeDownSubtreeDLikelihood(father);
661 }
662 
663 /******************************************************************************/
664 
665 void RNonHomogeneousTreeLikelihood::computeDownSubtreeDLikelihood(const Node* node)
666 {
667  const Node* father = node->getFather();
668  // We assume that the _dLikelihoods array has been filled for the current node 'node'.
669  // We will evaluate the array for the father node.
670  if (father == 0)
671  return; // We reached the root!
672 
673  // Compute dLikelihoods array for the father node.
674  // Fist initialize to 1:
675  VVVdouble* _dLikelihoods_father = &likelihoodData_->getDLikelihoodArray(father->getId());
676  size_t nbSites = _dLikelihoods_father->size();
677  for (size_t i = 0; i < nbSites; i++)
678  {
679  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
680  for (size_t c = 0; c < nbClasses_; c++)
681  {
682  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
683  for (size_t s = 0; s < nbStates_; s++)
684  {
685  (*_dLikelihoods_father_i_c)[s] = 1.;
686  }
687  }
688  }
689 
690  size_t nbNodes = father->getNumberOfSons();
691  for (size_t l = 0; l < nbNodes; l++)
692  {
693  const Node* son = father->getSon(l);
694 
695  VVVdouble* pxy__son = &pxy_[son->getId()];
696  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
697 
698  if (son == node)
699  {
700  VVVdouble* _dLikelihoods_son = &likelihoodData_->getDLikelihoodArray(son->getId());
701  for (size_t i = 0; i < nbSites; i++)
702  {
703  VVdouble* _dLikelihoods_son_i = &(*_dLikelihoods_son)[(*_patternLinks_father_son)[i]];
704  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
705  for (size_t c = 0; c < nbClasses_; c++)
706  {
707  Vdouble* _dLikelihoods_son_i_c = &(*_dLikelihoods_son_i)[c];
708  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
709  VVdouble* pxy__son_c = &(*pxy__son)[c];
710  for (size_t x = 0; x < nbStates_; x++)
711  {
712  double dl = 0;
713  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
714  for (size_t y = 0; y < nbStates_; y++)
715  {
716  dl += (*pxy__son_c_x)[y] * (*_dLikelihoods_son_i_c)[y];
717  }
718  (*_dLikelihoods_father_i_c)[x] *= dl;
719  }
720  }
721  }
722  }
723  else
724  {
725  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
726  for (size_t i = 0; i < nbSites; i++)
727  {
728  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
729  VVdouble* _dLikelihoods_father_i = &(*_dLikelihoods_father)[i];
730  for (size_t c = 0; c < nbClasses_; c++)
731  {
732  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
733  Vdouble* _dLikelihoods_father_i_c = &(*_dLikelihoods_father_i)[c];
734  VVdouble* pxy__son_c = &(*pxy__son)[c];
735  for (size_t x = 0; x < nbStates_; x++)
736  {
737  double dl = 0;
738  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
739  for (size_t y = 0; y < nbStates_; y++)
740  {
741  dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
742  }
743  (*_dLikelihoods_father_i_c)[x] *= dl;
744  }
745  }
746  }
747  }
748  }
749 
750  // Next step: move toward grand father...
751  computeDownSubtreeDLikelihood(father);
752 }
753 
754 /******************************************************************************
755 * Second Order Derivatives *
756 ******************************************************************************/
757 double RNonHomogeneousTreeLikelihood::getD2LikelihoodForASiteForARateClass(
758  size_t site,
759  size_t rateClass) const
760 {
761  double d2l = 0;
762  Vdouble* d2la = &likelihoodData_->getD2LikelihoodArray(tree_->getRootNode()->getId())[likelihoodData_->getRootArrayPosition(site)][rateClass];
763  for (size_t i = 0; i < nbStates_; i++)
764  {
765  d2l += (*d2la)[i] * rootFreqs_[i];
766  }
767  return d2l;
768 }
769 
770 /******************************************************************************/
771 
772 double RNonHomogeneousTreeLikelihood::getD2LikelihoodForASite(size_t site) const
773 {
774  // Derivative of the sum is the sum of derivatives:
775  double d2l = 0;
776  for (size_t i = 0; i < nbClasses_; i++)
777  {
778  d2l += getD2LikelihoodForASiteForARateClass(site, i) * rateDistribution_->getProbability(i);
779  }
780  return d2l;
781 }
782 
783 /******************************************************************************/
784 
785 double RNonHomogeneousTreeLikelihood::getD2LogLikelihoodForASite(size_t site) const
786 {
787  return getD2LikelihoodForASite(site) / getLikelihoodForASite(site)
788  - pow( getDLikelihoodForASite(site) / getLikelihoodForASite(site), 2);
789 }
790 
791 /******************************************************************************/
792 
793 double RNonHomogeneousTreeLikelihood::getD2LogLikelihood() const
794 {
795  // Derivative of the sum is the sum of derivatives:
796  double dl = 0;
797  for (size_t i = 0; i < nbSites_; i++)
798  {
799  dl += getD2LogLikelihoodForASite(i);
800  }
801  return dl;
802 }
803 
804 /******************************************************************************/
805 
806 double RNonHomogeneousTreeLikelihood::getSecondOrderDerivative(const string& variable) const
807 {
808  if (!hasParameter(variable))
809  throw ParameterNotFoundException("RNonHomogeneousTreeLikelihood::getSecondOrderDerivative().", variable);
810  if (getRateDistributionParameters().hasParameter(variable))
811  {
812  throw Exception("Derivatives respective to rate distribution parameter are not implemented.");
813  }
814  if (getSubstitutionModelParameters().hasParameter(variable))
815  {
816  throw Exception("Derivatives respective to substitution model parameters are not implemented.");
817  }
818 
819  const_cast<RNonHomogeneousTreeLikelihood*>(this)->computeTreeD2Likelihood(variable);
820  return -getD2LogLikelihood();
821 }
822 
823 /******************************************************************************/
824 
825 void RNonHomogeneousTreeLikelihood::computeTreeD2Likelihood(const string& variable)
826 {
827  if (variable == "BrLenRoot")
828  {
829  const Node* father = tree_->getRootNode();
830 
831  // Compute dLikelihoods array for the father node.
832  // Fist initialize to 1:
833  VVVdouble* _d2Likelihoods_father = &likelihoodData_->getD2LikelihoodArray(father->getId());
834  size_t nbSites = _d2Likelihoods_father->size();
835  for (size_t i = 0; i < nbSites; i++)
836  {
837  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
838  for (size_t c = 0; c < nbClasses_; c++)
839  {
840  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
841  for (size_t s = 0; s < nbStates_; s++)
842  {
843  (*_d2Likelihoods_father_i_c)[s] = 1.;
844  }
845  }
846  }
847 
848  size_t nbNodes = father->getNumberOfSons();
849  for (size_t l = 0; l < nbNodes; l++)
850  {
851  const Node* son = father->getSon(l);
852 
853  if (son->getId() == root1_)
854  {
855  const Node* root1 = father->getSon(0);
856  const Node* root2 = father->getSon(1);
857  vector<size_t>* _patternLinks_fatherroot1_ = &likelihoodData_->getArrayPositions(father->getId(), root1->getId());
858  vector<size_t>* _patternLinks_fatherroot2_ = &likelihoodData_->getArrayPositions(father->getId(), root2->getId());
859  VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(root1->getId());
860  VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(root2->getId());
861  double pos = getParameterValue("RootPosition");
862 
863  VVVdouble* d2pxy_root1_ = &d2pxy_[root1_];
864  VVVdouble* d2pxy_root2_ = &d2pxy_[root2_];
865  VVVdouble* dpxy_root1_ = &dpxy_[root1_];
866  VVVdouble* dpxy_root2_ = &dpxy_[root2_];
867  VVVdouble* pxy_root1_ = &pxy_[root1_];
868  VVVdouble* pxy_root2_ = &pxy_[root2_];
869  for (size_t i = 0; i < nbSites; i++)
870  {
871  VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
872  VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
873  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
874  for (size_t c = 0; c < nbClasses_; c++)
875  {
876  Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
877  Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
878  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
879  VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
880  VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
881  VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
882  VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
883  VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
884  VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
885  for (size_t x = 0; x < nbStates_; x++)
886  {
887  Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
888  Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
889  Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
890  Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
891  Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
892  Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
893  double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
894  for (size_t y = 0; y < nbStates_; y++)
895  {
896  d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
897  d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
898  dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
899  dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
900  l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
901  l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
902  }
903  double d2l = pos * pos * d2l1 * l2 + (1. - pos) * (1. - pos) * d2l2 * l1 + 2 * pos * (1. - pos) * dl1 * dl2;
904  (*_d2Likelihoods_father_i_c)[x] *= d2l;
905  }
906  }
907  }
908  }
909  else if (son->getId() == root2_)
910  {
911  // Do nothing, this was accounted when dealing with root1_
912  }
913  else
914  {
915  // Account for a putative multifurcation:
916  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
917  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
918 
919  VVVdouble* pxy__son = &pxy_[son->getId()];
920  for (size_t i = 0; i < nbSites; i++)
921  {
922  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
923  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
924  for (size_t c = 0; c < nbClasses_; c++)
925  {
926  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
927  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
928  VVdouble* pxy__son_c = &(*pxy__son)[c];
929  for (size_t x = 0; x < nbStates_; x++)
930  {
931  double d2l = 0;
932  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
933  for (size_t y = 0; y < nbStates_; y++)
934  {
935  d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
936  }
937  (*_d2Likelihoods_father_i_c)[x] *= d2l;
938  }
939  }
940  }
941  }
942  }
943  return;
944  }
945  else if (variable == "RootPosition")
946  {
947  const Node* father = tree_->getRootNode();
948 
949  // Compute dLikelihoods array for the father node.
950  // Fist initialize to 1:
951  VVVdouble* _d2Likelihoods_father = &likelihoodData_->getD2LikelihoodArray(father->getId());
952  size_t nbSites = _d2Likelihoods_father->size();
953  for (size_t i = 0; i < nbSites; i++)
954  {
955  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
956  for (size_t c = 0; c < nbClasses_; c++)
957  {
958  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
959  for (size_t s = 0; s < nbStates_; s++)
960  {
961  (*_d2Likelihoods_father_i_c)[s] = 1.;
962  }
963  }
964  }
965 
966  size_t nbNodes = father->getNumberOfSons();
967  for (size_t l = 0; l < nbNodes; l++)
968  {
969  const Node* son = father->getSon(l);
970 
971  if (son->getId() == root1_)
972  {
973  const Node* root1 = father->getSon(0);
974  const Node* root2 = father->getSon(1);
975  vector<size_t>* _patternLinks_fatherroot1_ = &likelihoodData_->getArrayPositions(father->getId(), root1->getId());
976  vector<size_t>* _patternLinks_fatherroot2_ = &likelihoodData_->getArrayPositions(father->getId(), root2->getId());
977  VVVdouble* _likelihoodsroot1_ = &likelihoodData_->getLikelihoodArray(root1->getId());
978  VVVdouble* _likelihoodsroot2_ = &likelihoodData_->getLikelihoodArray(root2->getId());
979  double len = getParameterValue("BrLenRoot");
980 
981  VVVdouble* d2pxy_root1_ = &d2pxy_[root1_];
982  VVVdouble* d2pxy_root2_ = &d2pxy_[root2_];
983  VVVdouble* dpxy_root1_ = &dpxy_[root1_];
984  VVVdouble* dpxy_root2_ = &dpxy_[root2_];
985  VVVdouble* pxy_root1_ = &pxy_[root1_];
986  VVVdouble* pxy_root2_ = &pxy_[root2_];
987  for (size_t i = 0; i < nbSites; i++)
988  {
989  VVdouble* _likelihoodsroot1__i = &(*_likelihoodsroot1_)[(*_patternLinks_fatherroot1_)[i]];
990  VVdouble* _likelihoodsroot2__i = &(*_likelihoodsroot2_)[(*_patternLinks_fatherroot2_)[i]];
991  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
992  for (size_t c = 0; c < nbClasses_; c++)
993  {
994  Vdouble* _likelihoodsroot1__i_c = &(*_likelihoodsroot1__i)[c];
995  Vdouble* _likelihoodsroot2__i_c = &(*_likelihoodsroot2__i)[c];
996  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
997  VVdouble* d2pxy_root1__c = &(*d2pxy_root1_)[c];
998  VVdouble* d2pxy_root2__c = &(*d2pxy_root2_)[c];
999  VVdouble* dpxy_root1__c = &(*dpxy_root1_)[c];
1000  VVdouble* dpxy_root2__c = &(*dpxy_root2_)[c];
1001  VVdouble* pxy_root1__c = &(*pxy_root1_)[c];
1002  VVdouble* pxy_root2__c = &(*pxy_root2_)[c];
1003  for (size_t x = 0; x < nbStates_; x++)
1004  {
1005  Vdouble* d2pxy_root1__c_x = &(*d2pxy_root1__c)[x];
1006  Vdouble* d2pxy_root2__c_x = &(*d2pxy_root2__c)[x];
1007  Vdouble* dpxy_root1__c_x = &(*dpxy_root1__c)[x];
1008  Vdouble* dpxy_root2__c_x = &(*dpxy_root2__c)[x];
1009  Vdouble* pxy_root1__c_x = &(*pxy_root1__c)[x];
1010  Vdouble* pxy_root2__c_x = &(*pxy_root2__c)[x];
1011  double d2l1 = 0, d2l2 = 0, dl1 = 0, dl2 = 0, l1 = 0, l2 = 0;
1012  for (size_t y = 0; y < nbStates_; y++)
1013  {
1014  d2l1 += (*d2pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1015  d2l2 += (*d2pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1016  dl1 += (*dpxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1017  dl2 += (*dpxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1018  l1 += (*pxy_root1__c_x)[y] * (*_likelihoodsroot1__i_c)[y];
1019  l2 += (*pxy_root2__c_x)[y] * (*_likelihoodsroot2__i_c)[y];
1020  }
1021  double d2l = len * len * (d2l1 * l2 + d2l2 * l1 - 2 * dl1 * dl2);
1022  (*_d2Likelihoods_father_i_c)[x] *= d2l;
1023  }
1024  }
1025  }
1026  }
1027  else if (son->getId() == root2_)
1028  {
1029  // Do nothing, this was accounted when dealing with root1_
1030  }
1031  else
1032  {
1033  // Account for a putative multifurcation:
1034  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
1035  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
1036 
1037  VVVdouble* pxy__son = &pxy_[son->getId()];
1038  for (size_t i = 0; i < nbSites; i++)
1039  {
1040  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1041  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1042  for (size_t c = 0; c < nbClasses_; c++)
1043  {
1044  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1045  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1046  VVdouble* pxy__son_c = &(*pxy__son)[c];
1047  for (size_t x = 0; x < nbStates_; x++)
1048  {
1049  double d2l = 0;
1050  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1051  for (size_t y = 0; y < nbStates_; y++)
1052  {
1053  d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1054  }
1055  (*_d2Likelihoods_father_i_c)[x] *= d2l;
1056  }
1057  }
1058  }
1059  }
1060  }
1061  return;
1062  }
1063 
1064  // Get the node with the branch whose length must be derivated:
1065  size_t brI = TextTools::to<size_t>(variable.substr(5));
1066  const Node* branch = nodes_[brI];
1067  const Node* father = branch->getFather();
1068 
1069  // Compute dLikelihoods array for the father node.
1070  // Fist initialize to 1:
1071  VVVdouble* _d2Likelihoods_father = &likelihoodData_->getD2LikelihoodArray(father->getId());
1072  size_t nbSites = _d2Likelihoods_father->size();
1073  for (size_t i = 0; i < nbSites; i++)
1074  {
1075  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1076  for (size_t c = 0; c < nbClasses_; c++)
1077  {
1078  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1079  for (size_t s = 0; s < nbStates_; s++)
1080  {
1081  (*_d2Likelihoods_father_i_c)[s] = 1.;
1082  }
1083  }
1084  }
1085 
1086  size_t nbNodes = father->getNumberOfSons();
1087  for (size_t l = 0; l < nbNodes; l++)
1088  {
1089  const Node* son = father->getSon(l);
1090 
1091  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
1092  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
1093 
1094  if (son == branch)
1095  {
1096  VVVdouble* d2pxy__son = &d2pxy_[son->getId()];
1097  for (size_t i = 0; i < nbSites; i++)
1098  {
1099  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1100  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1101  for (size_t c = 0; c < nbClasses_; c++)
1102  {
1103  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1104  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1105  VVdouble* d2pxy__son_c = &(*d2pxy__son)[c];
1106  for (size_t x = 0; x < nbStates_; x++)
1107  {
1108  double d2l = 0;
1109  Vdouble* d2pxy__son_c_x = &(*d2pxy__son_c)[x];
1110  for (size_t y = 0; y < nbStates_; y++)
1111  {
1112  d2l += (*d2pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1113  }
1114  (*_d2Likelihoods_father_i_c)[x] *= d2l;
1115  }
1116  }
1117  }
1118  }
1119  else
1120  {
1121  VVVdouble* pxy__son = &pxy_[son->getId()];
1122  for (size_t i = 0; i < nbSites; i++)
1123  {
1124  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1125  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1126  for (size_t c = 0; c < nbClasses_; c++)
1127  {
1128  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1129  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1130  VVdouble* pxy__son_c = &(*pxy__son)[c];
1131  for (size_t x = 0; x < nbStates_; x++)
1132  {
1133  double d2l = 0;
1134  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1135  for (size_t y = 0; y < nbStates_; y++)
1136  {
1137  d2l += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1138  }
1139  (*_d2Likelihoods_father_i_c)[x] *= d2l;
1140  }
1141  }
1142  }
1143  }
1144  }
1145 
1146  // Now we go down the tree toward the root node:
1147  computeDownSubtreeD2Likelihood(father);
1148 }
1149 
1150 /******************************************************************************/
1151 
1152 void RNonHomogeneousTreeLikelihood::computeDownSubtreeD2Likelihood(const Node* node)
1153 {
1154  const Node* father = node->getFather();
1155  // We assume that the _dLikelihoods array has been filled for the current node 'node'.
1156  // We will evaluate the array for the father node.
1157  if (father == 0)
1158  return; // We reached the root!
1159 
1160  // Compute dLikelihoods array for the father node.
1161  // Fist initialize to 1:
1162  VVVdouble* _d2Likelihoods_father = &likelihoodData_->getD2LikelihoodArray(father->getId());
1163  size_t nbSites = _d2Likelihoods_father->size();
1164  for (size_t i = 0; i < nbSites; i++)
1165  {
1166  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1167  for (size_t c = 0; c < nbClasses_; c++)
1168  {
1169  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1170  for (size_t s = 0; s < nbStates_; s++)
1171  {
1172  (*_d2Likelihoods_father_i_c)[s] = 1.;
1173  }
1174  }
1175  }
1176 
1177  size_t nbNodes = father->getNumberOfSons();
1178  for (size_t l = 0; l < nbNodes; l++)
1179  {
1180  const Node* son = father->getSon(l);
1181 
1182  VVVdouble* pxy__son = &pxy_[son->getId()];
1183  vector<size_t>* _patternLinks_father_son = &likelihoodData_->getArrayPositions(father->getId(), son->getId());
1184 
1185  if (son == node)
1186  {
1187  VVVdouble* _d2Likelihoods_son = &likelihoodData_->getD2LikelihoodArray(son->getId());
1188  for (size_t i = 0; i < nbSites; i++)
1189  {
1190  VVdouble* _d2Likelihoods_son_i = &(*_d2Likelihoods_son)[(*_patternLinks_father_son)[i]];
1191  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1192  for (size_t c = 0; c < nbClasses_; c++)
1193  {
1194  Vdouble* _d2Likelihoods_son_i_c = &(*_d2Likelihoods_son_i)[c];
1195  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1196  VVdouble* pxy__son_c = &(*pxy__son)[c];
1197  for (size_t x = 0; x < nbStates_; x++)
1198  {
1199  double d2l = 0;
1200  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1201  for (size_t y = 0; y < nbStates_; y++)
1202  {
1203  d2l += (*pxy__son_c_x)[y] * (*_d2Likelihoods_son_i_c)[y];
1204  }
1205  (*_d2Likelihoods_father_i_c)[x] *= d2l;
1206  }
1207  }
1208  }
1209  }
1210  else
1211  {
1212  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
1213  for (size_t i = 0; i < nbSites; i++)
1214  {
1215  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_father_son)[i]];
1216  VVdouble* _d2Likelihoods_father_i = &(*_d2Likelihoods_father)[i];
1217  for (size_t c = 0; c < nbClasses_; c++)
1218  {
1219  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1220  Vdouble* _d2Likelihoods_father_i_c = &(*_d2Likelihoods_father_i)[c];
1221  VVdouble* pxy__son_c = &(*pxy__son)[c];
1222  for (size_t x = 0; x < nbStates_; x++)
1223  {
1224  double dl = 0;
1225  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1226  for (size_t y = 0; y < nbStates_; y++)
1227  {
1228  dl += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1229  }
1230  (*_d2Likelihoods_father_i_c)[x] *= dl;
1231  }
1232  }
1233  }
1234  }
1235  }
1236 
1237  // Next step: move toward grand father...
1238  computeDownSubtreeD2Likelihood(father);
1239 }
1240 
1241 /******************************************************************************/
1242 
1243 void RNonHomogeneousTreeLikelihood::computeTreeLikelihood()
1244 {
1245  computeSubtreeLikelihood(tree_->getRootNode());
1246 }
1247 
1248 /******************************************************************************/
1249 
1250 void RNonHomogeneousTreeLikelihood::computeSubtreeLikelihood(const Node* node)
1251 {
1252  if (node->isLeaf())
1253  return;
1254 
1255  size_t nbSites = likelihoodData_->getLikelihoodArray(node->getId()).size();
1256  size_t nbNodes = node->getNumberOfSons();
1257 
1258  // Must reset the likelihood array first (i.e. set all of them to 1):
1259  VVVdouble* _likelihoods_node = &likelihoodData_->getLikelihoodArray(node->getId());
1260  for (size_t i = 0; i < nbSites; i++)
1261  {
1262  // For each site in the sequence,
1263  VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
1264  for (size_t c = 0; c < nbClasses_; c++)
1265  {
1266  // For each rate class,
1267  Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
1268  for (size_t x = 0; x < nbStates_; x++)
1269  {
1270  // For each initial state,
1271  (*_likelihoods_node_i_c)[x] = 1.;
1272  }
1273  }
1274  }
1275 
1276  for (size_t l = 0; l < nbNodes; l++)
1277  {
1278  // For each son node,
1279 
1280  const Node* son = node->getSon(l);
1281 
1282  computeSubtreeLikelihood(son); // Recursive method:
1283 
1284  VVVdouble* pxy__son = &pxy_[son->getId()];
1285  vector<size_t>* _patternLinks_node_son = &likelihoodData_->getArrayPositions(node->getId(), son->getId());
1286  VVVdouble* _likelihoods_son = &likelihoodData_->getLikelihoodArray(son->getId());
1287 
1288  for (size_t i = 0; i < nbSites; i++)
1289  {
1290  // For each site in the sequence,
1291  VVdouble* _likelihoods_son_i = &(*_likelihoods_son)[(*_patternLinks_node_son)[i]];
1292  VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
1293  for (size_t c = 0; c < nbClasses_; c++)
1294  {
1295  // For each rate class,
1296  Vdouble* _likelihoods_son_i_c = &(*_likelihoods_son_i)[c];
1297  Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
1298  VVdouble* pxy__son_c = &(*pxy__son)[c];
1299 
1300  for (size_t x = 0; x < nbStates_; x++)
1301  {
1302  // For each initial state,
1303  Vdouble* pxy__son_c_x = &(*pxy__son_c)[x];
1304  double likelihood = 0;
1305  for (size_t y = 0; y < nbStates_; y++)
1306  {
1307  likelihood += (*pxy__son_c_x)[y] * (*_likelihoods_son_i_c)[y];
1308  }
1309  (*_likelihoods_node_i_c)[x] *= likelihood;
1310  }
1311  }
1312  }
1313  }
1314 }
1315 
1316 
1317 /******************************************************************************/
1318 
1319 void RNonHomogeneousTreeLikelihood::displayLikelihood(const Node* node)
1320 {
1321  cout << "Likelihoods at node " << node->getName() << ": " << endl;
1322  displayLikelihoodArray(likelihoodData_->getLikelihoodArray(node->getId()));
1323  cout << " ***" << endl;
1324 }
bpp::AbstractDiscreteRatesAcrossSitesTreeLikelihood::displayLikelihoodArray
static void displayLikelihoodArray(const VVVdouble &likelihoodArray)
Print the likelihood array to terminal (debugging tool).
Definition: AbstractDiscreteRatesAcrossSitesTreeLikelihood.cpp:251
bpp::AbstractNonHomogeneousTreeLikelihood
Partial implementation for branch non-homogeneous models of the TreeLikelihood interface.
Definition: AbstractNonHomogeneousTreeLikelihood.h:25
bpp::AbstractNonHomogeneousTreeLikelihood::operator=
AbstractNonHomogeneousTreeLikelihood & operator=(const AbstractNonHomogeneousTreeLikelihood &lik)
Assignation operator.
Definition: AbstractNonHomogeneousTreeLikelihood.cpp:93
bpp::AbstractNonHomogeneousTreeLikelihood::computeAllTransitionProbabilities
virtual void computeAllTransitionProbabilities()
Fill the pxy_, dpxy_ and d2pxy_ arrays for all nodes.
Definition: AbstractNonHomogeneousTreeLikelihood.cpp:370
bpp::AbstractNonHomogeneousTreeLikelihood::getSubstitutionModelParameters
ParameterList getSubstitutionModelParameters() const override
Get the parameters associated to substitution model(s).
Definition: AbstractNonHomogeneousTreeLikelihood.cpp:258
bpp::AbstractNonHomogeneousTreeLikelihood::computeTransitionProbabilitiesForNode
virtual void computeTransitionProbabilitiesForNode(const Node *node)
Fill the pxy_, dpxy_ and d2pxy_ arrays for one node.
Definition: AbstractNonHomogeneousTreeLikelihood.cpp:383
bpp::AbstractNonHomogeneousTreeLikelihood::idToNode_
std::map< int, const Node * > idToNode_
An index linking nodes to their id, for faster access than the getNode() method.
Definition: AbstractNonHomogeneousTreeLikelihood.h:79
bpp::AbstractNonHomogeneousTreeLikelihood::nodes_
std::vector< Node * > nodes_
Pointer toward all nodes in the tree.
Definition: AbstractNonHomogeneousTreeLikelihood.h:74
bpp::AbstractNonHomogeneousTreeLikelihood::applyParameters
virtual void applyParameters()
All parameters are stored in a parameter list. This function apply these parameters to the substituti...
Definition: AbstractNonHomogeneousTreeLikelihood.cpp:285
bpp::AbstractNonHomogeneousTreeLikelihood::getRateDistributionParameters
ParameterList getRateDistributionParameters() const override
Get the parameters associated to the rate distribution.
Definition: AbstractNonHomogeneousTreeLikelihood.h:153
bpp::AbstractParametrizable::setParametersValues
void setParametersValues(const ParameterList &parameters) override
bpp::AbstractParametrizable::hasParameter
bool hasParameter(const std::string &name) const override
bpp::AbstractParametrizable::getParameterValue
double getParameterValue(const std::string &name) const override
bpp::AbstractTreeLikelihood::data
const AlignmentDataInterface & data() const
Get the dataset for which the likelihood must be evaluated.
Definition: AbstractTreeLikelihood.h:262
bpp::AbstractTreeLikelihood::tree
const Tree & tree() const
Get the tree (topology and branch lengths).
Definition: AbstractTreeLikelihood.h:269
bpp::AbstractTreeLikelihood::data_
std::unique_ptr< const AlignmentDataInterface > data_
Definition: AbstractTreeLikelihood.h:208
bpp::AbstractTreeLikelihood::tree_
std::shared_ptr< TreeTemplate< Node > > tree_
Definition: AbstractTreeLikelihood.h:209
bpp::ApplicationTools::displayTask
static void displayTask(const std::string &text, bool eof=false)
bpp::ApplicationTools::displayTaskDone
static void displayTaskDone()
bpp::ApplicationTools::displayResult
static void displayResult(const std::string &text, const T &result)
bpp::Node
The phylogenetic node class.
Definition: Node.h:59
bpp::Node::getName
virtual std::string getName() const
Get the name associated to this node, if there is one, otherwise throw a NodeException.
Definition: Node.h:203
bpp::Node::getId
virtual int getId() const
Get the node's id.
Definition: Node.h:170
bpp::Node::getSon
virtual const Node * getSon(size_t pos) const
Definition: Node.h:362
bpp::Node::getFather
virtual const Node * getFather() const
Get the father of this node is there is one.
Definition: Node.h:306
bpp::Node::isLeaf
virtual bool isLeaf() const
Definition: Node.h:667
bpp::Node::getNumberOfSons
virtual size_t getNumberOfSons() const
Definition: Node.h:355
bpp::ParameterList::size
size_t size() const
bpp::ParameterList::getCommonParametersWith
virtual ParameterList getCommonParametersWith(const ParameterList &params) const
bpp::ParameterList::getParameterNames
virtual std::vector< std::string > getParameterNames() const
bpp::PatternTools::getSequenceSubset
static std::unique_ptr< AlignmentDataInterface > getSequenceSubset(const AlignmentDataInterface &sequenceSet, const std::shared_ptr< N > node, const AssociationTreeGraphImplObserver< N, E, I > &tree)
Extract the sequences corresponding to a given subtree.
Definition: PatternTools.h:44
bpp::RNonHomogeneousTreeLikelihood
This class implement the 'traditional' way of computing likelihood for a tree, allowing for non-homog...
Definition: RNonHomogeneousTreeLikelihood.h:51
bpp::RNonHomogeneousTreeLikelihood::getFirstOrderDerivative
double getFirstOrderDerivative(const std::string &variable) const
Definition: RNonHomogeneousTreeLikelihood.cpp:335
bpp::RNonHomogeneousTreeLikelihood::init_
void init_(bool usePatterns)
Method called by constructors.
Definition: RNonHomogeneousTreeLikelihood.cpp:58
bpp::RNonHomogeneousTreeLikelihood::getD2LikelihoodForASite
virtual double getD2LikelihoodForASite(size_t site) const
Definition: RNonHomogeneousTreeLikelihood.cpp:772
bpp::RNonHomogeneousTreeLikelihood::setParameters
void setParameters(const ParameterList &parameters)
Implements the Function interface.
Definition: RNonHomogeneousTreeLikelihood.cpp:215
bpp::RNonHomogeneousTreeLikelihood::getLikelihood
double getLikelihood() const
Get the likelihood for the whole dataset.
Definition: RNonHomogeneousTreeLikelihood.cpp:118
bpp::RNonHomogeneousTreeLikelihood::computeTreeDLikelihood
virtual void computeTreeDLikelihood(const std::string &variable)
Definition: RNonHomogeneousTreeLikelihood.cpp:354
bpp::RNonHomogeneousTreeLikelihood::operator=
RNonHomogeneousTreeLikelihood & operator=(const RNonHomogeneousTreeLikelihood &lik)
Definition: RNonHomogeneousTreeLikelihood.cpp:80
bpp::RNonHomogeneousTreeLikelihood::getLikelihoodForASite
double getLikelihoodForASite(size_t site) const
Get the likelihood for a site.
Definition: RNonHomogeneousTreeLikelihood.cpp:148
bpp::RNonHomogeneousTreeLikelihood::RNonHomogeneousTreeLikelihood
RNonHomogeneousTreeLikelihood(const Tree &tree, std::shared_ptr< SubstitutionModelSet > modelSet, std::shared_ptr< DiscreteDistributionInterface > rDist, bool verbose=true, bool usePatterns=true, bool reparametrizeRoot=false)
Build a new NonHomogeneousTreeLikelihood object without data.
Definition: RNonHomogeneousTreeLikelihood.cpp:20
bpp::RNonHomogeneousTreeLikelihood::getLogLikelihood
double getLogLikelihood() const
Get the logarithm of the likelihood for the whole dataset.
Definition: RNonHomogeneousTreeLikelihood.cpp:130
bpp::RNonHomogeneousTreeLikelihood::fireParameterChanged
void fireParameterChanged(const ParameterList &params)
Definition: RNonHomogeneousTreeLikelihood.cpp:222
bpp::RNonHomogeneousTreeLikelihood::getDLikelihoodForASite
virtual double getDLikelihoodForASite(size_t site) const
Definition: RNonHomogeneousTreeLikelihood.cpp:301
bpp::RNonHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClassForAState
double getLogLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
Get the logarithm of the likelihood for a site knowing its rate class and its ancestral state.
Definition: RNonHomogeneousTreeLikelihood.cpp:208
bpp::RNonHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClassForAState
double getLikelihoodForASiteForARateClassForAState(size_t site, size_t rateClass, int state) const
Get the likelihood for a site knowing its rate class and its ancestral state.
Definition: RNonHomogeneousTreeLikelihood.cpp:201
bpp::RNonHomogeneousTreeLikelihood::getD2LogLikelihoodForASite
virtual double getD2LogLikelihoodForASite(size_t site) const
Definition: RNonHomogeneousTreeLikelihood.cpp:785
bpp::RNonHomogeneousTreeLikelihood::getLikelihoodForASiteForARateClass
double getLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the likelihood for a site knowing its rate class.
Definition: RNonHomogeneousTreeLikelihood.cpp:175
bpp::RNonHomogeneousTreeLikelihood::setData
void setData(const AlignmentDataInterface &sites)
Set the dataset for which the likelihood must be evaluated.
Definition: RNonHomogeneousTreeLikelihood.cpp:96
bpp::RNonHomogeneousTreeLikelihood::getDLikelihoodForASiteForARateClass
virtual double getDLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Definition: RNonHomogeneousTreeLikelihood.cpp:286
bpp::RNonHomogeneousTreeLikelihood::likelihoodData_
std::shared_ptr< DRASRTreeLikelihoodData > likelihoodData_
Definition: RNonHomogeneousTreeLikelihood.h:53
bpp::RNonHomogeneousTreeLikelihood::getDLogLikelihoodForASite
virtual double getDLogLikelihoodForASite(size_t site) const
Definition: RNonHomogeneousTreeLikelihood.cpp:314
bpp::RNonHomogeneousTreeLikelihood::displayLikelihood
virtual void displayLikelihood(const Node *node)
This method is mainly for debugging purpose.
Definition: RNonHomogeneousTreeLikelihood.cpp:1319
bpp::RNonHomogeneousTreeLikelihood::computeDownSubtreeDLikelihood
virtual void computeDownSubtreeDLikelihood(const Node *)
Definition: RNonHomogeneousTreeLikelihood.cpp:665
bpp::RNonHomogeneousTreeLikelihood::getD2LikelihoodForASiteForARateClass
virtual double getD2LikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Definition: RNonHomogeneousTreeLikelihood.cpp:757
bpp::RNonHomogeneousTreeLikelihood::computeTreeD2Likelihood
virtual void computeTreeD2Likelihood(const std::string &variable)
Definition: RNonHomogeneousTreeLikelihood.cpp:825
bpp::RNonHomogeneousTreeLikelihood::getLogLikelihoodForASite
double getLogLikelihoodForASite(size_t site) const
Get the logarithm of the likelihood for a site.
Definition: RNonHomogeneousTreeLikelihood.cpp:160
bpp::RNonHomogeneousTreeLikelihood::getLogLikelihoodForASiteForARateClass
double getLogLikelihoodForASiteForARateClass(size_t site, size_t rateClass) const
Get the logarithm of the likelihood for a site knowing its rate class.
Definition: RNonHomogeneousTreeLikelihood.cpp:188
bpp::RNonHomogeneousTreeLikelihood::getSecondOrderDerivative
double getSecondOrderDerivative(const std::string &variable) const
Definition: RNonHomogeneousTreeLikelihood.cpp:806
bpp::RNonHomogeneousTreeLikelihood::computeSubtreeLikelihood
virtual void computeSubtreeLikelihood(const Node *node)
Compute the likelihood for a subtree defined by the Tree::Node node.
Definition: RNonHomogeneousTreeLikelihood.cpp:1250
bpp::Tree
Interface for phylogenetic tree objects.
Definition: Tree.h:115
bpp::VectorTools::vectorUnion
static std::vector< T > vectorUnion(const std::vector< T > &vec1, const std::vector< T > &vec2)
bpp::TextTools::toString
std::string toString(T t)
bpp
Defines the basic types of data flow nodes.
bpp::log
double log(const ExtendedFloat &ef)
Definition: ExtendedFloat.h:443
bpp::Vdouble
std::vector< double > Vdouble
bpp::pow
ExtendedFloat pow(const ExtendedFloat &ef, double exp)
Definition: ExtendedFloat.h:469
bpp::VVVdouble
std::vector< VVdouble > VVVdouble
bpp::VVdouble
std::vector< Vdouble > VVdouble