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
11using namespace bpp;
12
13// From the STL:
14#include <iostream>
15
16using namespace std;
17
18/******************************************************************************/
19
20RNonHomogeneousTreeLikelihood::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
38RNonHomogeneousTreeLikelihood::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
58void RNonHomogeneousTreeLikelihood::init_(bool usePatterns)
59{
60 likelihoodData_ = make_shared<DRASRTreeLikelihoodData>(
61 tree_,
62 rateDistribution_->getNumberOfCategories(),
63 usePatterns);
64}
65
66/******************************************************************************/
67
68RNonHomogeneousTreeLikelihood::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
80RNonHomogeneousTreeLikelihood& 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
92RNonHomogeneousTreeLikelihood::~RNonHomogeneousTreeLikelihood() {}
93
94/******************************************************************************/
95
96void 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
118double 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
130double 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
148double 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
160double 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
175double 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
188double 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
201double 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
208double 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
215void RNonHomogeneousTreeLikelihood::setParameters(const ParameterList& parameters)
216{
217 setParametersValues(parameters);
218}
219
220/******************************************************************************/
221
222void 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
276double 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 ******************************************************************************/
286double 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
301double 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
314double 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
322double 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
335double 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
354void 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
665void 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 ******************************************************************************/
757double 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
772double 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
785double RNonHomogeneousTreeLikelihood::getD2LogLikelihoodForASite(size_t site) const
786{
787 return getD2LikelihoodForASite(site) / getLikelihoodForASite(site)
788 - pow( getDLikelihoodForASite(site) / getLikelihoodForASite(site), 2);
789}
790
791/******************************************************************************/
792
793double 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
806double 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
825void 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
1152void 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
1243void RNonHomogeneousTreeLikelihood::computeTreeLikelihood()
1244{
1245 computeSubtreeLikelihood(tree_->getRootNode());
1246}
1247
1248/******************************************************************************/
1249
1250void 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
1319void 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::getFather
virtual const Node * getFather() const
Get the father of this node is there is one.
Definition: Node.h:306
bpp::Node::getSon
virtual const Node * getSon(size_t pos) const
Definition: Node.h:362
bpp::Node::isLeaf
virtual bool isLeaf() const
Definition: Node.h:679
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