Skip to content

Commit e9512fc

Browse files
trongnhanuittrongnhanuit
committed
de-integrate all local references before computing the absolute Lh at root
1 parent a5f261d commit e9512fc

File tree

2 files changed

+123
-38
lines changed

2 files changed

+123
-38
lines changed

‎tree/tree.cpp‎

Lines changed: 107 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -1432,11 +1432,14 @@ RealNumType cmaple::Tree::computeLhTemplate() {
14321432
refreshAllLhs<num_states>();
14331433

14341434
// initialize the total_lh by the likelihood from root
1435-
RealNumType total_lh =
1435+
/*RealNumType total_lh =
14361436
nodes[root_vector_index]
14371437
.getPartialLh(TOP)
14381438
->computeAbsoluteLhAtRoot<num_states>(node_mutations[root_vector_index], aln,
1439-
model, cumulative_base);
1439+
model, cumulative_base);*/
1440+
RealNumType total_lh = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
1441+
nodes[root_vector_index].getPartialLh(TOP),
1442+
cmaple::Index(root_vector_index, TOP));
14401443

14411444
// perform a DFS to add likelihood contributions from each internal nodes
14421445
total_lh += performDFS<&cmaple::Tree::computeLhContribution<num_states>>();
@@ -4637,16 +4640,20 @@ void cmaple::Tree::placeSubTreeAtNode(
46374640
const std::unique_ptr<SeqRegions>& lower_regions =
46384641
selected_node.getPartialLh(
46394642
TOP); // ->getPartialLhAtNode(aln, model, threshold_prob);
4640-
old_root_lh = lower_regions->computeAbsoluteLhAtRoot<num_states>(
4641-
node_mutations[root_vector_index], aln, model, cumulative_base);
4643+
/*old_root_lh = lower_regions->computeAbsoluteLhAtRoot<num_states>(
4644+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
4645+
old_root_lh = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
4646+
lower_regions, cmaple::Index(selected_node_vec, TOP));
46424647

46434648
// merge 2 lower vector into one
46444649
best_parent_lh = lower_regions->mergeTwoLowers<num_states>(
46454650
best_parent_regions, default_blength, *subtree_regions,
46464651
new_branch_length, aln, model, cumulative_rate, threshold_prob, true);
46474652

4648-
best_parent_lh += best_parent_regions->computeAbsoluteLhAtRoot<num_states>(
4649-
node_mutations[root_vector_index], aln, model, cumulative_base);
4653+
/*best_parent_lh += best_parent_regions->computeAbsoluteLhAtRoot<num_states>(
4654+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
4655+
best_parent_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
4656+
best_parent_regions, cmaple::Index(selected_node_vec, TOP));
46504657

46514658
// Try shorter branch lengths at root
46524659
best_root_blength = default_blength;
@@ -5294,9 +5301,11 @@ void cmaple::Tree::tryShorterBranchAtRoot(
52945301
new_root_lh = lower_regions->mergeTwoLowers<num_states>(
52955302
merged_root_sample_regions, new_blength, *sample, fixed_blength, aln,
52965303
model, cumulative_rate, params->threshold_prob, true);
5297-
new_root_lh +=
5304+
/*new_root_lh +=
52985305
merged_root_sample_regions->computeAbsoluteLhAtRoot<num_states>(
5299-
node_mutations[root_vector_index], aln, model, cumulative_base);
5306+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
5307+
new_root_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
5308+
merged_root_sample_regions, cmaple::Index(root_vector_index, TOP));
53005309

53015310
if (new_root_lh > best_parent_lh) {
53025311
best_parent_lh = new_root_lh;
@@ -5333,8 +5342,10 @@ bool cmaple::Tree::tryShorterNewBranchAtRoot(
53335342
new_root_lh = lower_regions->mergeTwoLowers<num_states>(
53345343
new_root_lower_regions, fixed_blength, *sample, new_blength, aln, model,
53355344
cumulative_rate, params->threshold_prob, true);
5336-
new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5337-
node_mutations[root_vector_index], aln, model, cumulative_base);
5345+
/*new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5346+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
5347+
new_root_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
5348+
new_root_lower_regions, cmaple::Index(root_vector_index, TOP));
53385349

53395350
if (new_root_lh > best_parent_lh) {
53405351
best_parent_lh = new_root_lh;
@@ -5373,8 +5384,10 @@ bool cmaple::Tree::tryLongerNewBranchAtRoot(
53735384
new_root_lh = lower_regions->mergeTwoLowers<num_states>(
53745385
new_root_lower_regions, fixed_blength, *sample, new_blength, aln, model,
53755386
cumulative_rate, params->threshold_prob, true);
5376-
new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5377-
node_mutations[root_vector_index], aln, model, cumulative_base);
5387+
/*new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5388+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
5389+
new_root_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
5390+
new_root_lower_regions, cmaple::Index(root_vector_index, TOP));
53785391

53795392
if (new_root_lh > best_parent_lh) {
53805393
best_parent_lh = new_root_lh;
@@ -5411,9 +5424,11 @@ void cmaple::Tree::estimateLengthNewBranchAtRoot(
54115424
new_root_lower_regions, fixed_blength, *sample, best_length, aln, model,
54125425
cumulative_rate, params->threshold_prob, true);
54135426

5414-
best_parent_lh +=
5427+
/*best_parent_lh +=
54155428
new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5416-
node_mutations[root_vector_index], aln, model, cumulative_base);
5429+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
5430+
best_parent_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
5431+
new_root_lower_regions, cmaple::Index(root_vector_index, TOP));
54175432

54185433
// replacePartialLH(best_parent_regions, new_root_lower_regions);
54195434
best_parent_regions = std::move(new_root_lower_regions);
@@ -5444,8 +5459,10 @@ void cmaple::Tree::estimateLengthNewBranchAtRoot(
54445459
// it happens when the min blength is too large
54455460
if (new_root_lower_regions != nullptr)
54465461
{
5447-
new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5448-
node_mutations[root_vector_index], aln, model, cumulative_base);
5462+
/*new_root_lh += new_root_lower_regions->computeAbsoluteLhAtRoot<num_states>(
5463+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
5464+
new_root_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
5465+
new_root_lower_regions, cmaple::Index(root_vector_index, TOP));
54495466

54505467
if (new_root_lh > best_parent_lh) {
54515468
best_length = -1;
@@ -8047,11 +8064,14 @@ void cmaple::Tree::calculate_aRLT(const bool allow_replacing_ML_tree) {
80478064
node_lhs.reserve(static_cast<std::vector<cmaple::NodeLh>::size_type>(nodes.size() * 0.5));
80488065

80498066
// compute the likelihood at root
8050-
RealNumType lh_at_root =
8067+
/*RealNumType lh_at_root =
80518068
nodes[root_vector_index]
80528069
.getPartialLh(TOP)
80538070
->computeAbsoluteLhAtRoot<num_states>(node_mutations[root_vector_index], aln,
8054-
model, cumulative_base);
8071+
model, cumulative_base);*/
8072+
RealNumType lh_at_root = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
8073+
nodes[root_vector_index].getPartialLh(TOP),
8074+
cmaple::Index(root_vector_index, TOP));
80558075

80568076
// LT1 = tree_total_lh = likelihood at root + total likelihood contribution at
80578077
// all internal nodes
@@ -8181,9 +8201,12 @@ void cmaple::Tree::calSiteLhDiffRoot(
81818201
RealNumType best_parent_lh = parent_new_lower_lh->mergeTwoLowers<num_states>(
81828202
new_parent_new_lower_lh, parent_new_blength, *child_1_lower_regions,
81838203
child_1_blength, aln, model, cumulative_rate, threshold_prob, true);
8184-
best_parent_lh +=
8204+
/*best_parent_lh +=
81858205
new_parent_new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
8186-
node_mutations[root_vector_index], aln, model, cumulative_base);
8206+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
8207+
best_parent_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
8208+
new_parent_new_lower_lh,
8209+
cmaple::Index(root_vector_index, TOP));
81878210
// Try shorter branch lengths at root
81888211
tryShorterBranchAtRoot<num_states>(
81898212
child_1_lower_regions, parent_new_lower_lh, new_parent_new_lower_lh,
@@ -8893,9 +8916,12 @@ bool cmaple::Tree::calculateNNILhRoot(
88938916
RealNumType best_parent_lh = parent_new_lower_lh->mergeTwoLowers<num_states>(
88948917
new_parent_new_lower_lh, parent_new_blength, *child_1_lower_regions,
88958918
child_1_blength, aln, model, cumulative_rate, threshold_prob, true);
8896-
best_parent_lh +=
8919+
/*best_parent_lh +=
88978920
new_parent_new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
8898-
node_mutations[root_vector_index], aln, model, cumulative_base);
8921+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
8922+
best_parent_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
8923+
new_parent_new_lower_lh,
8924+
cmaple::Index(root_vector_index, TOP));
88998925
// Try shorter branch lengths at root
89008926
tryShorterBranchAtRoot<num_states>(
89018927
child_1_lower_regions, parent_new_lower_lh, new_parent_new_lower_lh,
@@ -8980,9 +9006,12 @@ bool cmaple::Tree::calculateNNILhRoot(
89809006
cumulative_rate, threshold_prob, true) -
89819007
node_lhs[parent.getNodelhIndex()].getLhContribution();
89829008
// 7.3 the absolute likelihood at root
8983-
lh_diff += new_parent_new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
9009+
/*lh_diff += new_parent_new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
89849010
node_mutations[root_vector_index], aln, model, cumulative_base)
8985-
- lh_at_root;
9011+
- lh_at_root;*/
9012+
lh_diff += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
9013+
new_parent_new_lower_lh, cmaple::Index(root_vector_index, TOP))
9014+
- lh_at_root;
89869015

89879016
// if we found an NNI neighbor with higher lh => replace the ML tree
89889017
if (lh_diff > 0) {
@@ -9324,9 +9353,11 @@ bool cmaple::Tree::calculateNNILhNonRoot(
93249353
else {
93259354
// re-calculate likelihood at root
93269355
// std::cout << "lh at root (before): " << lh_at_root << std::endl;
9327-
lh_diff += new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
9356+
/*lh_diff += new_lower_lh->computeAbsoluteLhAtRoot<num_states>(
93289357
node_mutations[root_vector_index], aln, model, cumulative_base) -
9329-
lh_at_root;
9358+
lh_at_root;*/
9359+
lh_diff += computeAbsLhAtRootDeintegratedAllMuts<num_states>(new_lower_lh,
9360+
cmaple::Index(root_vector_index, TOP)) - lh_at_root;
93309361
// std::cout << "lh at root (after): " <<
93319362
// new_lower_lh->computeAbsoluteLhAtRoot(num_states, model) <<
93329363
// std::endl;
@@ -9505,8 +9536,10 @@ void cmaple::Tree::replaceMLTreebyNNIRoot(
95059536
child_1_best_blength, aln, model, cumulative_rate, threshold_prob,
95069537
true));
95079538
// update the absolute likelihood at root
9508-
lh_at_root = parent_new_lower->computeAbsoluteLhAtRoot<num_states>(
9509-
node_mutations[root_vector_index], aln, model, cumulative_base);
9539+
/*lh_at_root = parent_new_lower->computeAbsoluteLhAtRoot<num_states>(
9540+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
9541+
lh_at_root = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
9542+
parent_new_lower, cmaple::Index(root_vector_index, TOP));
95109543

95119544
// traverse downward to update the upper_left/right_region until the changes
95129545
// is insignificant
@@ -9740,9 +9773,11 @@ void cmaple::Tree::replaceMLTreebyNNINonRoot(
97409773
else {
97419774
// re-calculate likelihood at root
97429775
// std::cout << "lh at root (before): " << lh_at_root << std::endl;
9743-
lh_at_root =
9776+
/*lh_at_root =
97449777
node.getPartialLh(TOP)->computeAbsoluteLhAtRoot<num_states>(
9745-
node_mutations[root_vector_index], aln, model, cumulative_base);
9778+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
9779+
lh_at_root = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
9780+
node.getPartialLh(TOP), cmaple::Index(node_vec, TOP));
97469781
// std::cout << "lh at root (after): " << lh_at_root << std::endl;
97479782

97489783
// stop traversing further
@@ -11347,9 +11382,11 @@ NumSeqsType cmaple::Tree::seekBestRoot()
1134711382
RealNumType lh_contribution_at_root = MIN_NEGATIVE;
1134811383
if (lower_regions_merged)
1134911384
{
11350-
lh_contribution_at_root = lower_regions_merged
11385+
/*lh_contribution_at_root = lower_regions_merged
1135111386
->computeAbsoluteLhAtRoot<num_states>(node_mutations[root_vector_index],
11352-
aln, model, cumulative_base);
11387+
aln, model, cumulative_base);*/
11388+
lh_contribution_at_root = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
11389+
lower_regions_merged, candidate_node.getNeighborIndex(TOP));
1135311390
}
1135411391

1135511392
// compute the total score, taking into account the likelihood deduction and contribution
@@ -11436,9 +11473,11 @@ void cmaple::Tree::addStartingRootCandidate(
1143611473

1143711474
// compute the current likelihood contribution by merging likelihood
1143811475
// at root with the state frequencies
11439-
RealNumType lh_contribution_at_root =
11476+
/*RealNumType lh_contribution_at_root =
1144011477
node.getPartialLh(TOP)->computeAbsoluteLhAtRoot<num_states>(
11441-
node_mutations[root_vector_index], aln, model, cumulative_base);
11478+
node_mutations[root_vector_index], aln, model, cumulative_base);*/
11479+
RealNumType lh_contribution_at_root = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
11480+
node.getPartialLh(TOP), cmaple::Index(node_vec_id, TOP));
1144211481

1144311482
// compute the likelihood contribution
1144411483
// by merging two lower regions from the two children
@@ -11794,3 +11833,37 @@ auto cmaple::Tree::makeReferenceNode(PhyloNode& node, const NumSeqsType& node_ve
1179411833
}
1179511834
}
1179611835
}
11836+
11837+
template <const StateType num_states>
11838+
auto cmaple::Tree::computeAbsLhAtRootDeintegratedAllMuts(
11839+
const std::unique_ptr<SeqRegions>& regions,
11840+
cmaple::Index node_index) -> RealNumType
11841+
{
11842+
// clone the original regions
11843+
std::unique_ptr<SeqRegions> mut_deintegrated_regions =
11844+
cmaple::make_unique<SeqRegions>(regions);
11845+
11846+
// traverse upward to de-integrate all mutations from local references
11847+
while (node_index.getMiniIndex() != UNDEFINED)
11848+
{
11849+
const NumSeqsType node_vec_index = node_index.getVectorIndex();
11850+
PhyloNode& node = nodes[node_vec_index];
11851+
11852+
// extract the mutations at the this node
11853+
std::unique_ptr<SeqRegions>& this_node_mutations =
11854+
node_mutations[node_vec_index];
11855+
// de-integrate mutations, if any
11856+
if (this_node_mutations && this_node_mutations->size())
11857+
{
11858+
mut_deintegrated_regions = mut_deintegrated_regions
11859+
->integrateMutations<num_states>(this_node_mutations, aln, true);
11860+
}
11861+
11862+
// move upward
11863+
node_index = node.getNeighborIndex(TOP);
11864+
}
11865+
11866+
// compute the absolute lh value after de-integrating all local references
11867+
return mut_deintegrated_regions->computeAbsoluteLhAtRoot<num_states>(
11868+
aln, model, cumulative_base);
11869+
}

‎tree/tree.h‎

Lines changed: 16 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2124,6 +2124,14 @@ bool isDiffFromOrigPlacement(
21242124
template <const StateType num_states>
21252125
auto makeReferenceNode(PhyloNode& node, const NumSeqsType& node_vec_index,
21262126
const int old_num_desc) -> void;
2127+
2128+
/**
2129+
Compute the absolute likelihood at root, de-integrating all local references
2130+
*/
2131+
template <const StateType num_states>
2132+
auto computeAbsLhAtRootDeintegratedAllMuts(
2133+
const std::unique_ptr<SeqRegions>& regions,
2134+
cmaple::Index node_index) -> RealNumType;
21272135
};
21282136

21292137
/*!
@@ -2668,17 +2676,21 @@ void cmaple::Tree::placeNewSampleAtNode(const Index selected_node_index,
26682676
selected_node.getPartialLh(TOP);
26692677

26702678
// compute the absolute lh
2671-
old_root_lh = lower_regions->computeAbsoluteLhAtRoot<num_states>(
2672-
selected_node_mutations, aln, model, cumulative_base);
2679+
/* old_root_lh = lower_regions->computeAbsoluteLhAtRoot<num_states>(
2680+
selected_node_mutations, aln, model, cumulative_base);*/
2681+
old_root_lh = computeAbsLhAtRootDeintegratedAllMuts<num_states>(
2682+
lower_regions, selected_node_index);
26732683

26742684
// merge 2 lower vector into one
26752685
RealNumType new_root_lh = lower_regions->mergeTwoLowers<num_states>(
26762686
best_parent_regions, default_blength, *sample, default_blength, aln,
26772687
model, cumulative_rate, threshold_prob, true);
26782688

26792689
// compute the absolute lh
2680-
new_root_lh += best_parent_regions->computeAbsoluteLhAtRoot<num_states>(
2681-
selected_node_mutations, aln, model, cumulative_base);
2690+
/*new_root_lh += best_parent_regions->computeAbsoluteLhAtRoot<num_states>(
2691+
selected_node_mutations, aln, model, cumulative_base);*/
2692+
new_root_lh += computeAbsLhAtRootDeintegratedAllMuts<num_states>(
2693+
best_parent_regions, selected_node_index);
26822694

26832695
best_parent_lh = new_root_lh;
26842696

0 commit comments

Comments
 (0)