1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
#include "bblLSProportionalEB.h"
#include "numRec.h"
#include "logFile.h"
#include "errorMsg.h"
bblLSProportionalEB::bblLSProportionalEB(tree& et, const vector<sequenceContainer>& sc, multipleStochasticProcess* msp, const gammaDistribution* pProportionDist, Vdouble& treeLikelihoodVec, const bool optimizeSelectedBranches, int maxIter, MDOUBLE epsilon)
{
_treeLikelihoodVec = optimizeBranches(et,sc,msp,pProportionDist,treeLikelihoodVec,optimizeSelectedBranches,maxIter,epsilon);
}
Vdouble bblLSProportionalEB::optimizeBranches(tree& et, const vector<sequenceContainer>& sc, multipleStochasticProcess* msp, const gammaDistribution* pProportionDist, Vdouble& inTreeLikelihoodVec, const bool optimizeSelectedBranches, int maxIter, MDOUBLE epsilon)
{
Vdouble treeLikelihoodVec;
if (inTreeLikelihoodVec.empty()){
treeLikelihoodVec = likelihoodComputation::getTreeLikelihoodProportionalAllPosAlphTheSame(et,sc,msp,pProportionDist);
}
else{
treeLikelihoodVec = inTreeLikelihoodVec;
}
MDOUBLE treeLikelihood = sumVdouble(treeLikelihoodVec);
LOGnOUT(5,<<"ll before bblLSr4sp"<<" logL="<<treeLikelihood<<endl);
vector<tree::nodeP> nodesV;
et.getAllNodes(nodesV,et.getRoot());
MDOUBLE prevIterL = VERYSMALL;
for (int iter = 0; iter < maxIter; ++iter)
{
if (treeLikelihood < prevIterL + epsilon){
treeLikelihoodVec = likelihoodComputation::getTreeLikelihoodProportionalAllPosAlphTheSame(et,sc,msp,pProportionDist);
return treeLikelihoodVec; //likelihood converged
}
prevIterL = treeLikelihood;
MDOUBLE paramFound;
MDOUBLE oldBl;
MDOUBLE newL;
for (int i=0; i<nodesV.size(); i++)
{
if (nodesV[i]->isRoot()) continue;
if((optimizeSelectedBranches) && (nodesV[i]->getComment() != "1")) continue; //only selected branhes will be optimized
oldBl = nodesV[i]->dis2father();
newL = -brent(0.0,oldBl,MAX_BRANCH_LENGTH,evalR4SPBranch(nodesV[i],et,sc,msp,pProportionDist),epsilon,¶mFound);
LOGnOUT(4,<<"oldL="<<treeLikelihood<<" newL="<<newL<<" BL="<<nodesV[i]->dis2father()<<endl);
if (newL >= treeLikelihood)
{
treeLikelihood = newL;
nodesV[i]->setDisToFather(paramFound);
}
else //likelihood went down!
{
nodesV[i]->setDisToFather(oldBl); //return to previous BL
LOGnOUT(4,<<"Likelihood went down. oldL="<<treeLikelihood<<" newL="<<newL<<" Do not update tree"<<endl);
}
}
}
treeLikelihoodVec = likelihoodComputation::getTreeLikelihoodProportionalAllPosAlphTheSame(et,sc,msp,pProportionDist);
return treeLikelihoodVec;
}
|
