From 29bc4ee1ff2cb344b0e6d765d06be3b633cebb42 Mon Sep 17 00:00:00 2001
From: Thomas Wong <thomaskf@gmail.com>
Date: Thu, 22 Feb 2024 14:39:44 +1100
Subject: [PATCH] Update the default of minimum branch length for MAST model to
 0.001

---
 tree/iqtreemix.cpp | 77 +++++++++++++++++++++++++++++++++++++++++++---
 tree/iqtreemix.h   | 15 +++++++--
 2 files changed, 86 insertions(+), 6 deletions(-)

diff --git a/tree/iqtreemix.cpp b/tree/iqtreemix.cpp
index 872528f1c..aced37993 100644
--- a/tree/iqtreemix.cpp
+++ b/tree/iqtreemix.cpp
@@ -8,7 +8,7 @@
 #include "iqtreemix.h"
 const double MIN_PROP = 0.001;
 const double MAX_PROP = 1000.0;
-const double MIN_LEN = 1e-4;
+// const double MIN_LEN = 1e-3;
 const double MAX_LEN = 1.0;
 const double LIKE_THRES = 0.1; // 10% more in team of likelihood value
 const double WEIGHT_EPSILON = 0.001;
@@ -1564,8 +1564,10 @@ void IQTreeMix::restoreCheckpoint() {
 void IQTreeMix::setMinBranchLen(Params& params) {
     size_t i;
     int num_prec;
+    
     if (params.min_branch_length <= 0.0) {
         params.min_branch_length = MAST_MIN_BRANCH_LEN;
+        /*
         if (size() > 0) {
             if (!at(0)->isSuperTree() && at(0)->getAlnNSite() >= 100000) {
                 params.min_branch_length = MAST_MIN_BRANCH_LEN; // 0.1 / (at(0)->getAlnNSite());
@@ -1585,8 +1587,9 @@ void IQTreeMix::setMinBranchLen(Params& params) {
             cout << "NOTE: minimal branch length is increased to " << params.min_branch_length << " because PoMo infers number of mutations and frequency shifts" << endl;
             cout.precision(3);
         }
+        */
     }
-    cout << setprecision(7) << "Minimum branch length is set to " << params.min_branch_length << endl;
+    // cout << setprecision(7) << "Minimum branch length is set to " << params.min_branch_length << endl;
 }
 
 /** set pointer of params variable */
@@ -1900,6 +1903,64 @@ void showDoubleArrayContent(string name, int dim, double* arr) {
     cout << endl;
 }
 
+/**
+    Estimate the average branch length for each tree
+    Based on this value to set the minimum value of branch length
+ */
+void IQTreeMix::estimateAvgBrLen() {
+    size_t i,j,k;
+    int freq;
+    int* parsimony_scores;
+    double avgLen;
+    vector<int> sumParScores;
+    
+    if (!parsi_computed) {
+        // compute parsimony scores for each tree along the patterns
+        // results are stored in the array patn_parsimony
+        computeParsimony();
+    }
+    
+    estAvgBrlen.clear();
+    sumParScores.clear();
+    for (i=0; i<ntree; i++)
+        sumParScores.push_back(0);
+    
+    k=0;
+    for (i=0; i<nptn; i++) {
+        freq = ptn_freq[i];
+        for (j=0; j<ntree; j++) {
+            sumParScores[j] += patn_parsimony[k] * freq;
+            k++;
+        }
+    }
+    
+    /*
+    cout << "sum of parsimony scores along all sites for each tree: ";
+    for (i=0; i<ntree; i++) {
+        if (i>0)
+            cout << ",";
+        cout << sumParScores[i];
+    }
+    cout << endl;
+    */
+    
+    for (i=0; i<ntree; i++) {
+        avgLen = (double) sumParScores[i] / getAlnNSite() / nbranch;
+        estAvgBrlen.push_back(avgLen);
+    }
+    
+    /*
+    cout << "estimated avg branch lengths: ";
+    for (i=0; i<ntree; i++) {
+        if (i>0)
+            cout <<",";
+        cout << estAvgBrlen[i];
+    }
+    cout << endl;
+    */
+}
+
+
 /**
     Initialize the tree weights using parsimony scores
     Idea:
@@ -2203,7 +2264,15 @@ string IQTreeMix::optimizeModelParameters(bool printInfo, double logl_epsilon) {
     
     // set all the branches of the same group to their weighted average for initialization of the branch lengths
     checkBranchGrp();
-    
+
+    // adjust the minimum branch length if params.min_branch_length > 0.1 * estimate_avg_branch_len
+    estimateAvgBrLen();
+    for (i = 0; i < ntree; i++) {
+        if (params->min_branch_length > 0.1 * estAvgBrlen[i])
+            params->min_branch_length = 0.1 * estAvgBrlen[i];
+    }
+    cout << setprecision(7) << "Minimum branch length is set to " << params->min_branch_length << endl;
+
     // show trees
     // cout << "Initial trees:" << endl;
     // showTree();
@@ -2939,7 +3008,7 @@ void IQTreeMix::setBounds(double *lower_bound, double *upper_bound, bool* bound_
         // optimization on branch length
         ndim = branch_group.size();
         for (i=0; i<ndim; i++) {
-            lower_bound[i+1] = MIN_LEN;
+            lower_bound[i+1] = params->min_branch_length;
             upper_bound[i+1] = MAX_LEN;
             bound_check[i+1] = false;
         }
diff --git a/tree/iqtreemix.h b/tree/iqtreemix.h
index 81dd6f74e..5e7d9fa2a 100644
--- a/tree/iqtreemix.h
+++ b/tree/iqtreemix.h
@@ -22,8 +22,8 @@
 //#include "phylokernelnew.h"
 #include <vectorclass/vectorclass.h>
 
-// by default, the minimum branch length for MAST model is 0.000001
-#define MAST_MIN_BRANCH_LEN 1e-6;
+// by default, the minimum branch length for MAST model is 0.001
+#define MAST_MIN_BRANCH_LEN 1e-3;
 
 // for checking the scaling for the likelihood values
 #define TINY_SCALE_DIFF 0.5
@@ -188,6 +188,12 @@ class IQTreeMix : public IQTree, public vector<IQTree*> {
      */
     virtual int testNumThreads();
     
+    /**
+        Estimate the average branch length for each tree
+        Based on this value to set the minimum value of branch length
+     */
+    void estimateAvgBrLen();
+
     /**
      Initialize the tree weights using parsimony scores
      Idea:
@@ -386,6 +392,11 @@ class IQTreeMix : public IQTree, public vector<IQTree*> {
      whether the site rates are linked, if exists
      */
     bool isLinkSiteRate;
+    
+    /**
+     the estimated average branch length for each tree
+     */
+    vector<double> estAvgBrlen;
 
 protected: