Skip to content

Commit ded3a1a

Browse files
trongnhanuittrongnhanuit
committed
- Update refreshNonLowerLhsFromParent(); updateLowerLh(); updateLowerLhAvoidUsingUpperLRLh(); computeLhContribution()
1 parent 7a52459 commit ded3a1a

File tree

1 file changed

+137
-11
lines changed

1 file changed

+137
-11
lines changed

‎tree/tree.cpp‎

Lines changed: 137 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -5766,11 +5766,23 @@ void cmaple::Tree::refreshNonLowerLhsFromParent(Index& node_index,
57665766
const RealNumType threshold_prob = params->threshold_prob;
57675767
const Index parent_index = node.getNeighborIndex(TOP);
57685768
PhyloNode& parent_node = nodes[parent_index.getVectorIndex()];
5769-
const std::unique_ptr<SeqRegions>& parent_upper_lr_lh =
5770-
parent_node.getPartialLh(
5771-
parent_index
5772-
.getMiniIndex()); // node->neighbor->getPartialLhAtNode(aln,
5773-
// model, threshold_prob);
5769+
5770+
// 0. extract the mutations at the selected node
5771+
std::unique_ptr<SeqRegions>& selected_node_mutations =
5772+
node_mutations[node_index.getVectorIndex()];
5773+
// 1. create a new parent_upper_lr_lh that integrate the mutations, if any
5774+
std::unique_ptr<SeqRegions> mut_integrated_parent_upper_lr_lh =
5775+
(selected_node_mutations && selected_node_mutations->size())
5776+
? parent_node.getPartialLh(parent_index.getMiniIndex())
5777+
->integrateMutations<num_states>(selected_node_mutations, aln)
5778+
: nullptr;
5779+
// 2. create the pointer that points to the appropriate upper_lr_regions
5780+
const std::unique_ptr<SeqRegions>* parent_upper_lr_lh_ptr =
5781+
(selected_node_mutations && selected_node_mutations->size())
5782+
? &(mut_integrated_parent_upper_lr_lh)
5783+
: &(parent_node.getPartialLh(parent_index.getMiniIndex()));
5784+
// 3. create a reference from that pointer
5785+
auto& parent_upper_lr_lh = *parent_upper_lr_lh_ptr;
57745786

57755787
// update the total lh, total lh at the mid-branch point of the current node
57765788
if (node.getUpperLength() > 0) // node->length > 0)
@@ -7681,13 +7693,48 @@ template <const StateType num_states>
76817693
void cmaple::Tree::updateLowerLh(RealNumType& total_lh,
76827694
std::unique_ptr<SeqRegions>& new_lower_lh,
76837695
PhyloNode& node,
7684-
const std::unique_ptr<SeqRegions>& lower_lh_1,
7685-
const std::unique_ptr<SeqRegions>& lower_lh_2,
7696+
const std::unique_ptr<SeqRegions>& ori_lower_lh_1,
7697+
const std::unique_ptr<SeqRegions>& ori_lower_lh_2,
76867698
const Index neighbor_1_index,
76877699
PhyloNode& neighbor_1,
76887700
const Index neighbor_2_index,
76897701
PhyloNode& neighbor_2,
76907702
const PositionType& seq_length) {
7703+
7704+
// 0. extract the mutations at neighbor_1
7705+
std::unique_ptr<SeqRegions>& neighbor_1_mutations =
7706+
node_mutations[neighbor_1_index.getVectorIndex()];
7707+
// 1. create a new regions that de-integrate the mutations, if any
7708+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_1_regions =
7709+
(neighbor_1_mutations && neighbor_1_mutations->size())
7710+
? ori_lower_lh_1
7711+
->integrateMutations<num_states>(neighbor_1_mutations, aln, true)
7712+
: nullptr;
7713+
// 2. create the pointer that points to the appropriate regions
7714+
const std::unique_ptr<SeqRegions>* neighbor_1_regions_ptr =
7715+
(neighbor_1_mutations && neighbor_1_mutations->size())
7716+
? &(mut_integrated_neighbor_1_regions)
7717+
: &(ori_lower_lh_1);
7718+
// 3. create a reference from that pointer
7719+
auto& lower_lh_1 = *neighbor_1_regions_ptr;
7720+
7721+
// 0. extract the mutations at neighbor_2
7722+
std::unique_ptr<SeqRegions>& neighbor_2_mutations =
7723+
node_mutations[neighbor_2_index.getVectorIndex()];
7724+
// 1. create a new regions that de-integrate the mutations, if any
7725+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_2_regions =
7726+
(neighbor_2_mutations && neighbor_2_mutations->size())
7727+
? ori_lower_lh_2
7728+
->integrateMutations<num_states>(neighbor_2_mutations, aln, true)
7729+
: nullptr;
7730+
// 2. create the pointer that points to the appropriate regions
7731+
const std::unique_ptr<SeqRegions>* neighbor_2_regions_ptr =
7732+
(neighbor_2_mutations && neighbor_2_mutations->size())
7733+
? &(mut_integrated_neighbor_2_regions)
7734+
: &(ori_lower_lh_2);
7735+
// 3. create a reference from that pointer
7736+
auto& lower_lh_2 = *neighbor_2_regions_ptr;
7737+
76917738
lower_lh_1->mergeTwoLowers<num_states>(
76927739
new_lower_lh, neighbor_1.getUpperLength(), *lower_lh_2,
76937740
neighbor_2.getUpperLength(), aln, model, cumulative_rate,
@@ -7736,13 +7783,48 @@ void cmaple::Tree::updateLowerLhAvoidUsingUpperLRLh(
77367783
RealNumType& total_lh,
77377784
std::unique_ptr<SeqRegions>& new_lower_lh,
77387785
PhyloNode& node,
7739-
const std::unique_ptr<SeqRegions>& lower_lh_1,
7740-
const std::unique_ptr<SeqRegions>& lower_lh_2,
7786+
const std::unique_ptr<SeqRegions>& ori_lower_lh_1,
7787+
const std::unique_ptr<SeqRegions>& ori_lower_lh_2,
77417788
const Index neighbor_1_index,
77427789
PhyloNode& neighbor_1,
77437790
const Index neighbor_2_index,
77447791
PhyloNode& neighbor_2,
77457792
const PositionType& seq_length) {
7793+
7794+
// 0. extract the mutations at neighbor_1
7795+
std::unique_ptr<SeqRegions>& neighbor_1_mutations =
7796+
node_mutations[neighbor_1_index.getVectorIndex()];
7797+
// 1. create a new regions that de-integrate the mutations, if any
7798+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_1_regions =
7799+
(neighbor_1_mutations && neighbor_1_mutations->size())
7800+
? ori_lower_lh_1
7801+
->integrateMutations<num_states>(neighbor_1_mutations, aln, true)
7802+
: nullptr;
7803+
// 2. create the pointer that points to the appropriate regions
7804+
const std::unique_ptr<SeqRegions>* neighbor_1_regions_ptr =
7805+
(neighbor_1_mutations && neighbor_1_mutations->size())
7806+
? &(mut_integrated_neighbor_1_regions)
7807+
: &(ori_lower_lh_1);
7808+
// 3. create a reference from that pointer
7809+
auto& lower_lh_1 = *neighbor_1_regions_ptr;
7810+
7811+
// 0. extract the mutations at neighbor_2
7812+
std::unique_ptr<SeqRegions>& neighbor_2_mutations =
7813+
node_mutations[neighbor_2_index.getVectorIndex()];
7814+
// 1. create a new regions that de-integrate the mutations, if any
7815+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_2_regions =
7816+
(neighbor_2_mutations && neighbor_2_mutations->size())
7817+
? ori_lower_lh_2
7818+
->integrateMutations<num_states>(neighbor_2_mutations, aln, true)
7819+
: nullptr;
7820+
// 2. create the pointer that points to the appropriate regions
7821+
const std::unique_ptr<SeqRegions>* neighbor_2_regions_ptr =
7822+
(neighbor_2_mutations && neighbor_2_mutations->size())
7823+
? &(mut_integrated_neighbor_2_regions)
7824+
: &(ori_lower_lh_2);
7825+
// 3. create a reference from that pointer
7826+
auto& lower_lh_2 = *neighbor_2_regions_ptr;
7827+
77467828
lower_lh_1->mergeTwoLowers<num_states>(
77477829
new_lower_lh, neighbor_1.getUpperLength(), *lower_lh_2,
77487830
neighbor_2.getUpperLength(), aln, model, cumulative_rate,
@@ -7787,18 +7869,57 @@ void cmaple::Tree::computeLhContribution(
77877869
RealNumType& total_lh,
77887870
std::unique_ptr<SeqRegions>& new_lower_lh,
77897871
PhyloNode& node,
7790-
const std::unique_ptr<SeqRegions>& lower_lh_1,
7791-
const std::unique_ptr<SeqRegions>& lower_lh_2,
7872+
const std::unique_ptr<SeqRegions>& ori_lower_lh_1,
7873+
const std::unique_ptr<SeqRegions>& ori_lower_lh_2,
77927874
const Index neighbor_1_index,
77937875
PhyloNode& neighbor_1,
77947876
const Index neighbor_2_index,
77957877
PhyloNode& neighbor_2,
77967878
const PositionType& seq_length) {
7879+
7880+
// 0. extract the mutations at neighbor_1
7881+
std::unique_ptr<SeqRegions>& neighbor_1_mutations =
7882+
node_mutations[neighbor_1_index.getVectorIndex()];
7883+
// 1. create a new regions that de-integrate the mutations, if any
7884+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_1_regions =
7885+
(neighbor_1_mutations && neighbor_1_mutations->size())
7886+
? ori_lower_lh_1
7887+
->integrateMutations<num_states>(neighbor_1_mutations, aln, true)
7888+
: nullptr;
7889+
// 2. create the pointer that points to the appropriate regions
7890+
const std::unique_ptr<SeqRegions>* neighbor_1_regions_ptr =
7891+
(neighbor_1_mutations && neighbor_1_mutations->size())
7892+
? &(mut_integrated_neighbor_1_regions)
7893+
: &(ori_lower_lh_1);
7894+
// 3. create a reference from that pointer
7895+
auto& lower_lh_1 = *neighbor_1_regions_ptr;
7896+
7897+
// 0. extract the mutations at neighbor_2
7898+
std::unique_ptr<SeqRegions>& neighbor_2_mutations =
7899+
node_mutations[neighbor_2_index.getVectorIndex()];
7900+
// 1. create a new regions that de-integrate the mutations, if any
7901+
std::unique_ptr<SeqRegions> mut_integrated_neighbor_2_regions =
7902+
(neighbor_2_mutations && neighbor_2_mutations->size())
7903+
? ori_lower_lh_2
7904+
->integrateMutations<num_states>(neighbor_2_mutations, aln, true)
7905+
: nullptr;
7906+
// 2. create the pointer that points to the appropriate regions
7907+
const std::unique_ptr<SeqRegions>* neighbor_2_regions_ptr =
7908+
(neighbor_2_mutations && neighbor_2_mutations->size())
7909+
? &(mut_integrated_neighbor_2_regions)
7910+
: &(ori_lower_lh_2);
7911+
// 3. create a reference from that pointer
7912+
auto& lower_lh_2 = *neighbor_2_regions_ptr;
7913+
77977914
RealNumType lh_contribution = lower_lh_1->mergeTwoLowers<num_states>(
77987915
new_lower_lh, neighbor_1.getUpperLength(), *lower_lh_2,
77997916
neighbor_2.getUpperLength(), aln, model, cumulative_rate,
78007917
params->threshold_prob, true);
78017918
total_lh += lh_contribution;
7919+
7920+
/*std::cout << "lh_contribution " << neighbor_1_index.getVectorIndex()
7921+
<< " " << neighbor_2_index.getVectorIndex() <<": "
7922+
<< lh_contribution << std::endl;*/
78027923

78037924
// record the likelihood contribution at this node
78047925
// if likelihood contribution of this node has not yet existed -> add a new
@@ -11560,6 +11681,11 @@ template <const StateType num_states>
1156011681
auto cmaple::Tree::makeReferenceNode(PhyloNode& node, const NumSeqsType& node_vec_index,
1156111682
const int old_num_desc) -> void
1156211683
{
11684+
// Show debug info
11685+
if (cmaple::verbose_mode >= cmaple::VB_DEBUG) {
11686+
std::cout << "Make a local reference at node " << node_vec_index << std::endl;
11687+
}
11688+
1156311689
// dummy variables
1156411690
const PositionType seq_length = static_cast<PositionType>(aln->ref_seq.size());
1156511691
const RealNumType threshold_prob = params->threshold_prob;

0 commit comments

Comments
 (0)