Skip to content

Commit a60b61b

Browse files
bqminhbqminh
authored
Merge pull request #41 from iqtree/piyumal
Urgent fix for MCMCTree control file generation and Hessian calculation with ASC flag.
2 parents d06ec38 + 33f9580 commit a60b61b

File tree

2 files changed

+28
-22
lines changed

2 files changed

+28
-22
lines changed

‎main/phyloanalysis.cpp‎

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -2830,14 +2830,14 @@ void printTrees(vector<string> trees, Params &params, string suffix) {
28302830
treesOut.close();
28312831
}
28322832

2833-
void processDervMCMCTree(double *gradient_vector, size_t branchNum, size_t nPtn, double *G_matrix,
2834-
const Map<RowVectorXd> &ptn_freq_diagonal, RowVectorXd &gradient_vector_eigen,
2833+
void processDervMCMCTree(double *gradient_vector, size_t branchNum, size_t mem_size, double *G_matrix,
2834+
RowVectorXd &ptn_freq_diagonal, RowVectorXd &gradient_vector_eigen,
28352835
MatrixXd &hessian, double *hessian_diagonal) {
28362836

28372837
Map<RowVectorXd> gradient_vector_eigen_mapped(gradient_vector, branchNum);
28382838
gradient_vector_eigen = gradient_vector_eigen_mapped;
28392839

2840-
Map<Matrix<double, Dynamic, Dynamic, RowMajor>> G_matrix_eigen(G_matrix, branchNum, nPtn);
2840+
Map<Matrix<double, Dynamic, Dynamic, RowMajor>> G_matrix_eigen(G_matrix, branchNum, mem_size);
28412841
MatrixXd G_matrix_eigen_t = G_matrix_eigen.transpose();
28422842
hessian = G_matrix_eigen * ptn_freq_diagonal.asDiagonal() * G_matrix_eigen_t;
28432843
hessian = (-1) * hessian;
@@ -2861,7 +2861,12 @@ void printMCMCFileFormat(PhyloTree *tree, MatrixXd &hessian, stringstream &tree_
28612861
size_t max_orig_nptn = get_safe_upper_limit(orig_nptn);
28622862
size_t nPtn = max_orig_nptn + tree->getModelFactory()->unobserved_ptns.size();
28632863
size_t nBranches = tree->branchNum;
2864-
Map<RowVectorXd> ptn_freq_diagonal(tree->ptn_freq, nPtn);
2864+
size_t mem_size = get_safe_upper_limit(tree->getAlnNPattern()) + max(get_safe_upper_limit(tree->getModel()->num_states),
2865+
get_safe_upper_limit(
2866+
tree->getModelFactory()->unobserved_ptns.size()));
2867+
2868+
RowVectorXd ptn_freq_diagonal = RowVectorXd::Zero(mem_size);
2869+
ptn_freq_diagonal.head(nPtn) = Map<RowVectorXd>(tree->ptn_freq, nPtn);
28652870

28662871
// Retrieve the branch lengths for MCMCTREE output file
28672872
DoubleVector branchLengths;
@@ -2908,10 +2913,7 @@ void printMCMCFileFormat(PhyloTree *tree, MatrixXd &hessian, stringstream &tree_
29082913
}
29092914
partBrMmap = partBrMmap2;
29102915
}
2911-
int numStates = tree->getModel()->num_states;
2912-
size_t mem_size = get_safe_upper_limit(tree->getAlnNSite()) + max(get_safe_upper_limit(numStates),
2913-
get_safe_upper_limit(
2914-
tree->getModelFactory()->unobserved_ptns.size()));
2916+
29152917
size_t branchNum = (superTree->params->partition_type != BRLEN_OPTIMIZE) ? superTree->branchNum
29162918
: tree->branchNum;
29172919
size_t g_matrix_size = branchNum * mem_size;
@@ -2927,7 +2929,7 @@ void printMCMCFileFormat(PhyloTree *tree, MatrixXd &hessian, stringstream &tree_
29272929
for (auto mapping: partBrMmap) {
29282930
int stree_branch_id = mapping.first;
29292931
int part_branch_id = mapping.second;
2930-
memmove(G_matrix_part+stree_branch_id*nPtn, G_matrix_sub_tree+part_branch_id*nPtn, sizeof(double)*nPtn);
2932+
memmove(G_matrix_part+stree_branch_id*mem_size, G_matrix_sub_tree+part_branch_id*mem_size, sizeof(double)*mem_size);
29312933
gradient_vector_part[stree_branch_id] = tree->gradient_vector[part_branch_id];
29322934
hessian_diagonal_part[stree_branch_id] = tree->hessian_diagonal[part_branch_id];
29332935
}
@@ -2943,25 +2945,22 @@ void printMCMCFileFormat(PhyloTree *tree, MatrixXd &hessian, stringstream &tree_
29432945
continue;
29442946
}
29452947
brVisitedMap[part_branch_id] = 1;
2946-
memmove(G_matrix_part+branchCounter*nPtn, G_matrix_sub_tree+part_branch_id*nPtn, sizeof(double)*nPtn);
2948+
memmove(G_matrix_part+branchCounter*mem_size, G_matrix_sub_tree+part_branch_id*mem_size, sizeof(double)*mem_size);
29472949
gradient_vector_part[branchCounter] = tree->gradient_vector[part_branch_id];
29482950
hessian_diagonal_part[branchCounter] = tree->hessian_diagonal[part_branch_id];
29492951
branchLengths2.push_back(branchLengths[part_branch_id]);
29502952
branchCounter++;
29512953
}
29522954
saveTreeMCMCTree(branchLengths2, branch_lengths_vector, tree, tree_stream);
29532955
}
2954-
processDervMCMCTree(gradient_vector_part, branchNum, nPtn, G_matrix_part, ptn_freq_diagonal,
2956+
processDervMCMCTree(gradient_vector_part, branchNum, mem_size, G_matrix_part, ptn_freq_diagonal,
29552957
gradient_vector_eigen, hessian, hessian_diagonal_part);
29562958
aligned_free(G_matrix_part);
29572959
aligned_free(gradient_vector_part);
29582960
aligned_free(hessian_diagonal_part);
29592961

29602962
} else {
29612963

2962-
int numStates = tree->getModel()->num_states;
2963-
size_t mem_size = get_safe_upper_limit(tree->getAlnNSite()) + max(get_safe_upper_limit(numStates),
2964-
get_safe_upper_limit(tree->getModelFactory()->unobserved_ptns.size()));
29652964
BranchVector singleAlnBranches;
29662965
tree->getBranches(singleAlnBranches);
29672966
vector<int> branch_ids;
@@ -2987,15 +2986,15 @@ void printMCMCFileFormat(PhyloTree *tree, MatrixXd &hessian, stringstream &tree_
29872986
int br_counter = 0;
29882987
DoubleVector branchLengthsNew;
29892988
for (int br_id: branch_ids) {
2990-
memmove(G_matrix_new+br_counter*nPtn, G_matrix_tree+br_id*nPtn, sizeof(double)*nPtn);
2989+
memmove(G_matrix_new+br_counter*mem_size, G_matrix_tree+br_id*mem_size, sizeof(double)*mem_size);
29912990
gradient_vector_new[br_counter] = tree->gradient_vector[br_id];
29922991
hessian_diagonal_new[br_counter] = tree->hessian_diagonal[br_id];
29932992
branchLengthsNew.push_back(branchLengths[br_id]);
29942993
br_counter++;
29952994
}
29962995

29972996

2998-
processDervMCMCTree(gradient_vector_new, nBranches, nPtn, G_matrix_new, ptn_freq_diagonal,
2997+
processDervMCMCTree(gradient_vector_new, nBranches, mem_size, G_matrix_new, ptn_freq_diagonal,
29992998
gradient_vector_eigen, hessian, hessian_diagonal_new);
30002999
saveTreeMCMCTree(branchLengthsNew, branch_lengths_vector, tree, tree_stream);
30013000

@@ -3045,7 +3044,7 @@ void printMCMCTreeCtlFile(IQTree *iqtree, ofstream &ctl, ofstream &dummyAlignmen
30453044
<< "RootAge = <1.0 * safe constraint on root age, used if no fossil for root." << endl << endl
30463045

30473046
<< "BDparas = " << Params::getInstance().mcmc_bds << " * birth rate, death rate, sampling priors for sampling times" << endl
3048-
<< "finetune = 1: 0.1 0.1 0.1 0.01 .5 * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr"
3047+
<< "finetune = 1: 0.1 0.1 0.1 0.01 .5 * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr" << endl
30493048
<< "print = 1 * 1: normal output; 2: verbose output" << endl << endl
30503049

30513050

@@ -3207,7 +3206,8 @@ void printHessian(IQTree *iqtree, int partition_type) {
32073206
}
32083207
}
32093208

3210-
printMCMCFileFormat(iqtree, hessian, tree_stream, branch_lengths_vector, gradient_vector_eigen, NULL, 0, iqtree->leftSingleRoot);
3209+
printMCMCFileFormat(iqtree, hessian, tree_stream, branch_lengths_vector, gradient_vector_eigen, NULL, 0,
3210+
iqtree->leftSingleRoot);
32113211

32123212
outfile << endl << iqtree->aln->getNSeq() << endl << endl;
32133213
outfile << tree_stream.str() << endl << endl;

‎utils/tools.cpp‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1498,7 +1498,7 @@ void parseArg(int argc, char *argv[], Params &params) {
14981498
params.date_outlier = -1.0;
14991499
params.dating_mf = false;
15001500
params.mcmc_clock = CORRELATED;
1501-
params.mcmc_bds = "1,1,0.5";
1501+
params.mcmc_bds = "1 1 0.5";
15021502
params.mcmc_iter = "20000, 100, 20000";
15031503

15041504
// added by TD
@@ -5670,10 +5670,14 @@ void parseArg(int argc, char *argv[], Params &params) {
56705670
params.mcmc_bds = argv[cnt];
56715671
StrVector mcmc_bds_vec;
56725672
convert_string_vec(params.mcmc_bds.c_str(), mcmc_bds_vec, ',');
5673-
if (mcmc_bds_vec.size()!=3 || !strcmp(mcmc_bds_vec[2].c_str(), ""))
5673+
if (mcmc_bds_vec.size() != 3 || mcmc_bds_vec[0].empty() ||
5674+
mcmc_bds_vec[1].empty() ||
5675+
mcmc_bds_vec[2].empty())
56745676
{
5675-
throw "three parameters should be set for birth-death model of MCMCtree (birth-rate, death-rate and sampling-fraction)";
5677+
throw
5678+
"three parameters should be set for birth-death model of MCMCtree (birth-rate, death-rate and sampling-fraction)";
56765679
}
5680+
params.mcmc_bds = mcmc_bds_vec[0] + " " + mcmc_bds_vec[1] + " " + mcmc_bds_vec[2];
56775681
continue;
56785682
}
56795683

@@ -5682,7 +5686,9 @@ void parseArg(int argc, char *argv[], Params &params) {
56825686
params.mcmc_iter = argv[cnt];
56835687
StrVector mcmc_iter_vec;
56845688
convert_string_vec(params.mcmc_iter.c_str(), mcmc_iter_vec, ',');
5685-
if (mcmc_iter_vec.size()!=3 || !strcmp(mcmc_iter_vec[2].c_str(), ""))
5689+
if (mcmc_iter_vec.size() != 3 || mcmc_iter_vec[0].empty() ||
5690+
mcmc_iter_vec[1].empty() ||
5691+
mcmc_iter_vec[2].empty())
56865692
{
56875693
throw "three parameters should be set for MCMCtree dating (Burin, samplefreq and nsamples)";
56885694
}

0 commit comments

Comments
 (0)