Skip to content

Commit

Permalink
Change the default minimum branch length for MAST model to 1e-4
Browse files Browse the repository at this point in the history
  • Loading branch information
thomaskf committed Feb 24, 2024
1 parent 2d2f473 commit 03ed4fd
Show file tree
Hide file tree
Showing 2 changed files with 3 additions and 75 deletions.
68 changes: 1 addition & 67 deletions tree/iqtreemix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1589,7 +1589,7 @@ void IQTreeMix::setMinBranchLen(Params& params) {
}
*/
}
// 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 */
Expand Down Expand Up @@ -1903,64 +1903,6 @@ 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:
Expand Down Expand Up @@ -2265,14 +2207,6 @@ 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();
Expand Down
10 changes: 2 additions & 8 deletions tree/iqtreemix.h
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@
//#include "phylokernelnew.h"
#include <vectorclass/vectorclass.h>

// by default, the minimum branch length for MAST model is 0.001
#define MAST_MIN_BRANCH_LEN 1e-3;
// by default, the minimum branch length for MAST model is 0.0001
#define MAST_MIN_BRANCH_LEN 1e-4;

// for checking the scaling for the likelihood values
#define TINY_SCALE_DIFF 0.5
Expand Down Expand Up @@ -188,12 +188,6 @@ 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:
Expand Down

0 comments on commit 03ed4fd

Please sign in to comment.