@@ -4224,10 +4224,13 @@ ModelPairSet PartitionFinder::getBetterPairsmAIC() {
42244224 ModelPairSet cur_better_pairs;
42254225 set<int > part_ids;
42264226
4227+ double cpu_time = getCPUTime ();
4228+ double real_time = getRealTime ();
4229+
42274230 // evaluate all pairs of partitions for merging
42284231 int index = 0 ;
42294232 for (auto it = sorted_pairs.begin (); it != sorted_pairs.end (); it++, index++) {
4230- // early stopping: if no better mAIC pairs until all pairs with cAIC < inf_score + 10 checked
4233+ // early stopping: if no better mAIC pairs until going through number pairs reach to number of current partitions
42314234 if (cur_better_pairs.size () == 0 && index > gene_sets.size ()) break ;
42324235
42334236 // check for compatibility
@@ -4270,6 +4273,13 @@ ModelPairSet PartitionFinder::getBetterPairsmAIC() {
42704273
42714274
42724275 cout << cur_better_pairs.size () << " compatible better partition pairs found based on mAIC" << endl;
4276+
4277+ cpu_time = getCPUTime () - cpu_time;
4278+ real_time = getRealTime () - real_time;
4279+ cout << " CPU time for computing mAIC scores: " << cpu_time << " seconds (" << convert_time (cpu_time) << " )" << endl;
4280+ cout << " Wall-clock time for computing mAIC scores: " << real_time << " seconds (" << convert_time (real_time) << " )" << endl;
4281+ cout << endl;
4282+
42734283 sorted_pairs.clear ();
42744284 return cur_better_pairs;
42754285}
@@ -5030,8 +5040,8 @@ void PartitionFinder::test_PartitionModel() {
50305040 better_pairs = getBetterPairsmAIC ();
50315041 bool is_maic_pairs_empty = better_pairs.empty ();
50325042
5033- if (switch_to_caic) {
5034- if (is_maic_pairs_empty ) {
5043+ if (is_maic_pairs_empty) {
5044+ if (switch_to_caic ) {
50355045 cout << " No better pairs based on mAIC after one round of cAIC merging and finish merging with previous optimal mAIC" << endl;
50365046 // load back back up
50375047 lhsum = lhsum_bu;
@@ -5044,12 +5054,6 @@ void PartitionFinder::test_PartitionModel() {
50445054 greedy_model_trees = greedy_model_trees_bu;
50455055 break ;
50465056 } else {
5047- cout << " Better pairs based on mAIC are found after one round of cAIC merging" << endl;
5048- switch_to_caic = false ;
5049- compatible_pairs = better_pairs;
5050- }
5051- } else {
5052- if (is_maic_pairs_empty) {
50535057 if (is_pairs_empty || gene_sets.size () == 2 ) {
50545058 cout << " No better pairs based on both mAIC and cAIC" << endl;
50555059 break ;
@@ -5066,10 +5070,15 @@ void PartitionFinder::test_PartitionModel() {
50665070 lenvec_bu = lenvec;
50675071 model_names_bu = model_names;
50685072 greedy_model_trees_bu = greedy_model_trees;
5069- } else {
5070- compatible_pairs = better_pairs;
50715073 }
5074+ } else {
5075+ if (switch_to_caic) {
5076+ cout << " Better pairs based on mAIC are found after one round of cAIC merging" << endl;
5077+ switch_to_caic = false ;
5078+ }
5079+ compatible_pairs = better_pairs;
50725080 }
5081+
50735082 }
50745083 // 2017-12-21: simultaneously merging better pairs
50755084 for (auto it_pair = compatible_pairs.begin (); it_pair != compatible_pairs.end (); it_pair++) {
@@ -5125,23 +5134,31 @@ void PartitionFinder::test_PartitionModel() {
51255134 // save and output mAIC after merging all pairs
51265135 if (params->marginal_lh_aic ) {
51275136 if (!switch_to_caic) {
5128- inf_score_maic = getmAICforMergeScheme (gene_sets, model_names, dfsum, true );
5129- cout << " Current partition model mAIC score: " << inf_score_maic
5130- << " (Marginal LnL: " << lh_marginal << " df:" << dfsum << " )" << endl;
5137+ double cur_score_maic = getmAICforMergeScheme (gene_sets, model_names, dfsum, true );
5138+ if (cur_score_maic < inf_score_maic) {
5139+ cout << " Current partition model mAIC score: " << cur_score_maic
5140+ << " (Marginal LnL: " << lh_marginal << " df:" << dfsum << " )" << endl;
5141+ } else {
5142+ cout << " [Rare case] Current partition model mAIC score: " << cur_score_maic
5143+ << " (Marginal LnL: " << lh_marginal << " df:" << dfsum << " ) is worse than "
5144+ << inf_score_maic << endl;
5145+ }
5146+ inf_score_maic = cur_score_maic;
51315147 } else {
51325148 double cur_score_maic = getmAICforMergeScheme (gene_sets, model_names, dfsum, true );
51335149 if (cur_score_maic < inf_score_maic) {
51345150 cout << " Current partition model mAIC score: " << cur_score_maic
51355151 << " (Marginal LnL: " << lh_marginal << " df:" << dfsum << " ) is better than "
51365152 << inf_score_maic << " , back to merging with mAIC since next round" << endl;
5137- cur_score_maic = inf_score_maic ;
5153+ inf_score_maic = cur_score_maic ;
51385154 switch_to_caic = false ;
51395155 } else {
51405156 cout << " Current partition model mAIC score: " << cur_score_maic
51415157 << " (Marginal LnL: " << lh_marginal << " df:" << dfsum << " ) is worse than "
51425158 << inf_score_maic << endl;
51435159 }
51445160 }
5161+ cout << endl;
51455162 }
51465163
51475164
0 commit comments