diff --git a/main/phyloanalysis.cpp b/main/phyloanalysis.cpp index c2dac05a5..a70cfe61c 100644 --- a/main/phyloanalysis.cpp +++ b/main/phyloanalysis.cpp @@ -57,7 +57,7 @@ #include "whtest/whtest_wrapper.h" #include "model/partitionmodel.h" #include "model/modelmixture.h" -#include "model/modelfactorymixlen.h" +#include "model/modelfactory.h" //#include "guidedbootstrap.h" #include "model/modelset.h" #include "utils/timeutil.h" @@ -3441,9 +3441,31 @@ void startTreeReconstruction(Params ¶ms, IQTree* &iqtree, ModelCheckpoint &m string best_subst_name, best_rate_name; map > nest_network; runModelFinder(params, *iqtree, model_info, best_subst_name, best_rate_name, nest_network); - runMixtureFinder(params, iqtree, model_info); - + + // ModelFinder or MixtureFinder could selected a heterotachy rate model, + // so we need to check and possibly switch the tree type + if (posRateHeterotachy(iqtree->aln->model_name) != string::npos && !iqtree->isMixlen()) { + ASSERT(!iqtree->isSuperTree() && !iqtree->isTreeMix()); + // create a new instance + IQTree* iqtree_new = new PhyloTreeMixlen(iqtree->aln); + iqtree_new->setCheckpoint(iqtree->getCheckpoint()); + if (!iqtree->constraintTree.empty()) { + iqtree_new->constraintTree.readConstraint(iqtree->constraintTree); + } + iqtree_new->removed_seqs = iqtree->removed_seqs; + iqtree_new->twin_seqs = iqtree->twin_seqs; + if (params.start_tree == STT_PLL_PARSIMONY || params.start_tree == STT_RANDOM_TREE || params.pll) { + /* Initialized all data structure for PLL*/ + iqtree_new->initializePLL(params); + } + iqtree_new->setParams(¶ms); + iqtree_new->copyPhyloTree(iqtree, false); + // replace iqtree object + delete iqtree; + iqtree = iqtree_new; + } + // if users want to perform tree dating (with mcmc) // and if ModelFinder was run, the traversal starting node was incidently deleted (after copyTree and restoreCheckpoint) // we have to delete tree nodes to force IQ-TREE to re-read the tree from the treefile @@ -3494,11 +3516,6 @@ void optimizeConTree(Params ¶ms, IQTree *tree) { tree->getCheckpoint()->put("contree", contree); } -// check whether it is tree-mixture model -bool isTreeMixture(Params& params) { - return (params.model_name.find("+T") != string::npos); -} - void runTreeReconstruction(Params ¶ms, IQTree* &iqtree) { iqtree->aln->checkAbsentStates("alignment"); @@ -3519,40 +3536,12 @@ void runTreeReconstruction(Params ¶ms, IQTree* &iqtree) { if ((params.start_tree == STT_BIONJ || params.iqp || params.leastSquareBranch) && !iqtree->root) { computeInitialDist(params, *iqtree); } - + /******************** Pass the parameter object params to IQTree *******************/ iqtree->setParams(¶ms); - - ModelsBlock *models_block = readModelsDefinition(params); initializeParams(params, *iqtree); - if (posRateHeterotachy(iqtree->aln->model_name) != string::npos && !iqtree->isMixlen()) { - // rate mixture - - if (iqtree->isTreeMix()) { - outError("Tree-mixture model does not work with rate mixture"); - } - - // create a new instance - IQTree* iqtree_new = new PhyloTreeMixlen(iqtree->aln, 0); - iqtree_new->setCheckpoint(iqtree->getCheckpoint()); - if (!iqtree->constraintTree.empty()) - iqtree_new->constraintTree.readConstraint(iqtree->constraintTree); - iqtree_new->removed_seqs = iqtree->removed_seqs; - iqtree_new->twin_seqs = iqtree->twin_seqs; - if (params.start_tree == STT_PLL_PARSIMONY || params.start_tree == STT_RANDOM_TREE || params.pll) { - /* Initialized all data structure for PLL*/ - iqtree_new->initializePLL(params); - } - iqtree_new->setParams(¶ms); - iqtree_new->copyPhyloTree(iqtree, false); - - // replace iqtree object - delete iqtree; - iqtree = iqtree_new; - } - if (!params.compute_ml_tree_only) { iqtree->setRootNode(params.root); } @@ -3567,17 +3556,13 @@ void runTreeReconstruction(Params ¶ms, IQTree* &iqtree) { iqtree->initSettings(params); /*********************** INITIAL MODEL OPTIMIZATION *****************/ - - if (!iqtree->getModelFactory()) { + ModelsBlock *models_block = readModelsDefinition(params); iqtree->initializeModel(params, iqtree->aln->model_name, models_block); + delete models_block; } - if (iqtree->getRate()->isHeterotachy() && !iqtree->isMixlen()) { - ASSERT(0 && "Heterotachy tree not properly created"); - } -// iqtree.restoreCheckpoint(); - delete models_block; +// iqtree.restoreCheckpoint(); // UpperBounds analysis. Here, to analyse the initial tree without any tree search or optimization /* @@ -4181,32 +4166,26 @@ void runMultipleTreeReconstruction(Params ¶ms, Alignment *alignment, IQTree init_random(params.ran_seed); IQTree *iqtree; - if (alignment->isSuperAlignment()){ - if(params.partition_type != BRLEN_OPTIMIZE){ - iqtree = new PhyloSuperTreePlen((SuperAlignment*) alignment, (PhyloSuperTree*) tree); + if (alignment->isSuperAlignment()) { + if (params.partition_type != BRLEN_OPTIMIZE) { + iqtree = new PhyloSuperTreePlen((SuperAlignment*)alignment, (PhyloSuperTree*)tree); } else { - iqtree = new PhyloSuperTree((SuperAlignment*) alignment, (PhyloSuperTree*) tree); + iqtree = new PhyloSuperTree((SuperAlignment*)alignment, (PhyloSuperTree*)tree); } } else { - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy(alignment->model_name); - - if (params.num_mixlen > 1) { - iqtree = new PhyloTreeMixlen(alignment, params.num_mixlen); - } else if (pos != string::npos) { - iqtree = new PhyloTreeMixlen(alignment, 0); - } else + if (posRateHeterotachy(alignment->model_name) != string::npos) { + iqtree = new PhyloTreeMixlen(alignment); + } else { iqtree = new IQTree(alignment); + } } - if (!tree->constraintTree.empty()) { iqtree->constraintTree.readConstraint(tree->constraintTree); } - // set checkpoint iqtree->setCheckpoint(tree->getCheckpoint()); iqtree->num_precision = tree->num_precision; - + // run reconstruction for this tree runTreeReconstruction(params, iqtree); // read in the output tree file stringstream ss; @@ -4264,8 +4243,7 @@ void runMultipleTreeReconstruction(Params ¶ms, Alignment *alignment, IQTree cout << "Run " << best_run+1 << " gave best log-likelihood: " << runLnL[best_run] << endl; - // initialize tree and model strucgture - ModelsBlock *models_block = readModelsDefinition(params); + // initialize tree and model structure tree->setParams(¶ms); // restore thread count for tree search after ModelFinder's per-partition reduction if (params.num_threads_max > params.num_threads) { @@ -4276,12 +4254,10 @@ void runMultipleTreeReconstruction(Params ¶ms, Alignment *alignment, IQTree } tree->setNumThreads(params.num_threads); if (!tree->getModelFactory()) { + ModelsBlock *models_block = readModelsDefinition(params); tree->initializeModel(params, tree->aln->model_name, models_block); + delete models_block; } - if (tree->getRate()->isHeterotachy() && !tree->isMixlen()) { - ASSERT(0 && "Heterotachy tree not properly created"); - } - delete models_block; // restore the tree and model from the best run tree->getCheckpoint()->startStruct("run" + convertIntToString(best_run+1)); @@ -4594,52 +4570,50 @@ void runStandardBootstrap(Params ¶ms, Alignment *alignment, IQTree *tree) { finish_random(); randstream = saved_randstream; - if (params.print_tree_lh && MPIHelper::getInstance().isMaster()) { - double prob; - bootstrap_alignment->multinomialProb(*alignment, prob); - ofstream boot_lh; - if (sample == 0) - boot_lh.open(bootlh_name.c_str()); - else - boot_lh.open(bootlh_name.c_str(), ios_base::out | ios_base::app); - boot_lh << "0\t" << prob << endl; - boot_lh.close(); + if (MPIHelper::getInstance().isMaster()) { + if (params.print_tree_lh) { + double prob; + bootstrap_alignment->multinomialProb(*alignment, prob); + ofstream boot_lh; + if (sample == 0) { + boot_lh.open(bootlh_name.c_str()); + } else { + boot_lh.open(bootlh_name.c_str(), ios_base::out | ios_base::app); + } + boot_lh << sample << "\t" << prob << endl; + boot_lh.close(); + } + if (params.print_bootaln) { + bootstrap_alignment->printAlignment(params.aln_output_format, bootaln_name.c_str(), true); + } + if (params.print_boot_site_freq) { + string sample_name = (string)params.out_prefix + "." + convertIntToString(sample); + printSiteStateFreq((sample_name + ".bootsitefreq").c_str(), bootstrap_alignment); + bootstrap_alignment->printAlignment(params.aln_output_format, (sample_name + ".bootaln").c_str()); + } } + IQTree *boot_tree; - if (alignment->isSuperAlignment()){ - if(params.partition_type != BRLEN_OPTIMIZE){ - boot_tree = new PhyloSuperTreePlen((SuperAlignment*) bootstrap_alignment, (PhyloSuperTree*) tree); + if (alignment->isSuperAlignment()) { + if (params.partition_type != BRLEN_OPTIMIZE) { + boot_tree = new PhyloSuperTreePlen((SuperAlignment*)bootstrap_alignment, (PhyloSuperTree*)tree); } else { - boot_tree = new PhyloSuperTree((SuperAlignment*) bootstrap_alignment, (PhyloSuperTree*) tree); + boot_tree = new PhyloSuperTree((SuperAlignment*)bootstrap_alignment, (PhyloSuperTree*)tree); } } else { - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy(alignment->model_name); - - if (params.num_mixlen > 1) { - boot_tree = new PhyloTreeMixlen(bootstrap_alignment, params.num_mixlen); - } else if (pos != string::npos) { - boot_tree = new PhyloTreeMixlen(bootstrap_alignment, 0); - } else + if (posRateHeterotachy(alignment->model_name) != string::npos) { + boot_tree = new PhyloTreeMixlen(bootstrap_alignment); + } else { boot_tree = new IQTree(bootstrap_alignment); + } } - if (params.print_bootaln && MPIHelper::getInstance().isMaster()) { - bootstrap_alignment->printAlignment(params.aln_output_format, bootaln_name.c_str(), true); - } - - if (params.print_boot_site_freq && MPIHelper::getInstance().isMaster()) { - printSiteStateFreq((((string)params.out_prefix)+"."+convertIntToString(sample)+".bootsitefreq").c_str(), bootstrap_alignment); - bootstrap_alignment->printAlignment(params.aln_output_format, (((string)params.out_prefix)+"."+convertIntToString(sample)+".bootaln").c_str()); - } - if (!tree->constraintTree.empty()) { boot_tree->constraintTree.readConstraint(tree->constraintTree); } - // set checkpoint boot_tree->setCheckpoint(tree->getCheckpoint()); boot_tree->num_precision = tree->num_precision; - + // run reconstruction for this tree runTreeReconstruction(params, boot_tree); // read in the output tree file stringstream ss; @@ -4832,24 +4806,26 @@ void convertAlignment(Params ¶ms, IQTree *iqtree) { 2016-08-04: compute a site frequency model for profile mixture model */ void computeSiteFrequencyModel(Params ¶ms, Alignment *alignment) { - - cout << endl << "===> COMPUTING SITE FREQUENCY MODEL BASED ON TREE FILE " << params.tree_freq_file << endl; ASSERT(params.tree_freq_file); - PhyloTree *tree = new PhyloTree(alignment); + cout << endl << "===> COMPUTING SITE FREQUENCY MODEL BASED ON TREE FILE " << params.tree_freq_file << endl; + double saved_min_branch_length = params.min_branch_length; + params.min_branch_length = 1e-6; + IQTree *tree; + if (posRateHeterotachy(alignment->model_name) != string::npos) { + tree = new PhyloTreeMixlen(alignment); + } else { + tree = new IQTree(alignment); + } tree->setParams(¶ms); + tree->setLikelihoodKernel(params.SSE); + tree->setNumThreads(params.num_threads); bool myrooted = params.is_rooted; tree->readTree(params.tree_freq_file, myrooted); - tree->setAlignment(alignment); tree->setRootNode(params.root); - + tree->setAlignment(alignment); ModelsBlock *models_block = readModelsDefinition(params); tree->setModelFactory(new ModelFactory(params, alignment->model_name, tree, models_block)); delete models_block; - tree->setModel(tree->getModelFactory()->model); - tree->setRate(tree->getModelFactory()->site_rate); - tree->setLikelihoodKernel(params.SSE); - tree->setNumThreads(params.num_threads); - if (!tree->getModel()->isMixture()) outError("No mixture model was specified!"); uint64_t mem_size = tree->getMemoryRequired(); @@ -4880,10 +4856,9 @@ void computeSiteFrequencyModel(Params ¶ms, Alignment *alignment) { } printSiteStateFreq(((string)params.out_prefix+".sitefreq").c_str(), tree, ptn_state_freq); params.print_site_state_freq = WSF_NONE; - delete [] ptn_state_freq; delete tree; - + params.min_branch_length = saved_min_branch_length; cout << endl << "===> CONTINUE ANALYSIS USING THE INFERRED SITE FREQUENCY MODEL" << endl; } @@ -4894,47 +4869,41 @@ void computeSiteFrequencyModel(Params ¶ms, Alignment *alignment) { IQTree *newIQTree(Params ¶ms, Alignment *alignment) { IQTree *tree; - - // check whether it is a tree mixture model - bool isTreeMix = isTreeMixture(params); - + // check whether it is a multitree mixture model + bool isTreeMix = (params.model_name.find("+T") != string::npos); if (alignment->isSuperAlignment()) { if (params.partition_type == TOPO_UNLINKED) { - if (isTreeMix) - outError("Error! The topology-unlinked option is not supported for multitree mixture model"); + if (isTreeMix) { + outError("Topology-unlinked partition model is not supported for multitree mixture model"); + } tree = new PhyloSuperTreeUnlinked((SuperAlignment*)alignment); - } else if(params.partition_type != BRLEN_OPTIMIZE){ - if (isTreeMix) - outError("Error! The proportional edge option is not supported for multitree mixture model"); - // initialize supertree - Proportional Edges case + } else if (params.partition_type != BRLEN_OPTIMIZE) { + if (isTreeMix) { + outError("Edge-linked partition models are not supported for multitree mixture model"); + } tree = new PhyloSuperTreePlen((SuperAlignment*)alignment, params.partition_type); } else { if (isTreeMix) { tree = new PhyloSuperHmm((SuperAlignment*)alignment, params); } else { - // initialize supertree stuff if user specifies partition file with -sp option tree = new PhyloSuperTree((SuperAlignment*)alignment); } } - // this alignment will actually be of type SuperAlignment - // alignment = tree->aln; - if (((PhyloSuperTree*)tree)->rescale_codon_brlen) + if (((PhyloSuperTree*)tree)->rescale_codon_brlen) { cout << "NOTE: Mixed codon and other data, branch lengths of codon partitions are rescaled by 3!" << endl; - + } } else { - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy(alignment->model_name); - - if (params.num_mixlen > 1) { - tree = new PhyloTreeMixlen(alignment, params.num_mixlen); - } else if (pos != string::npos) { - tree = new PhyloTreeMixlen(alignment, 0); + if (posRateHeterotachy(alignment->model_name) != string::npos) { + if (isTreeMix) { + outError("Heterotachy rate model is not supported for multitree mixture model"); + } + tree = new PhyloTreeMixlen(alignment); } else if (isTreeMix) { tree = new IQTreeMixHmm(params, alignment); - } else + } else { tree = new IQTree(alignment); + } } - return tree; } @@ -5633,80 +5602,65 @@ bool runCMaple(Params ¶ms) are unlinked between partitions */ void runUnlinkedPhyloAnalysis(Params ¶ms, Checkpoint *checkpoint) { - SuperAlignment *super_aln; - ASSERT(params.partition_file); - /****************** read in alignment **********************/ - // Partition model analysis - super_aln = new SuperAlignmentUnlinked(params); + SuperAlignment *super_aln = new SuperAlignmentUnlinked(params); PhyloSuperTree *super_tree = new PhyloSuperTree(super_aln); - /**** do separate tree reconstruction for each partition ***/ - MTreeSet part_trees; - if (params.user_file) { // reading user tree file for all partitions bool is_rooted = false; part_trees.readTrees(params.user_file, is_rooted, 0, super_aln->partitions.size()); - if (is_rooted) + if (is_rooted) { outError("Rooted trees not allowed: ", params.user_file); - if (part_trees.size() != super_aln->partitions.size()) + } + if (part_trees.size() != super_aln->partitions.size()) { outError("User tree file does not have the same number of trees as partitions"); + } params.user_file = nullptr; } - ModelCheckpoint *model_info = new ModelCheckpoint; - int part = 0; - for (auto alnit = super_aln->partitions.begin(); alnit != super_aln->partitions.end(); alnit++, part++) { - - checkpoint->startStruct((*alnit)->name); - - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy((*alnit)->model_name); + for (size_t part = 0; part < super_aln->partitions.size(); ++part) { + Alignment *part_aln = super_aln->partitions[part]; + checkpoint->startStruct(part_aln->name); IQTree *tree; - - if (params.num_mixlen > 1) { - tree = new PhyloTreeMixlen((*alnit), params.num_mixlen); - } else if (pos != string::npos) { - tree = new PhyloTreeMixlen((*alnit), 0); - } else - tree = new IQTree((*alnit)); - + if (posRateHeterotachy(part_aln->model_name) != string::npos) { + tree = new PhyloTreeMixlen(part_aln); + } else { + tree = new IQTree(part_aln); + } tree->setCheckpoint(checkpoint); if (checkpoint->getBool("finished")) { tree->restoreCheckpoint(); } else { - if (!part_trees.empty()) + if (!part_trees.empty()) { tree->copyTree(part_trees[part]); - + } startTreeReconstruction(params, tree, *model_info); // call main tree reconstruction - if (params.num_runs == 1) + if (params.num_runs == 1) { runTreeReconstruction(params, tree); - else + } else { runMultipleTreeReconstruction(params, tree->aln, tree); + } checkpoint->putBool("finished", true); checkpoint->dump(); } - super_tree->at(part)->copyTree(tree); - delete tree; checkpoint->endStruct(); } - IQTree *iqtree = super_tree; super_tree->setCheckpoint(checkpoint); startTreeReconstruction(params, iqtree, *model_info); runTreeReconstruction(params, iqtree); - if (MPIHelper::getInstance().isMaster()) + if (MPIHelper::getInstance().isMaster()) { reportPhyloAnalysis(params, *iqtree, *model_info); - + } + delete model_info; delete super_tree; delete super_aln; - delete model_info; } void assignBranchSupportNew(Params ¶ms) { diff --git a/main/phylotesting.cpp b/main/phylotesting.cpp index 34cfde377..481658cd1 100644 --- a/main/phylotesting.cpp +++ b/main/phylotesting.cpp @@ -36,7 +36,7 @@ #include "model/modelliemarkov.h" #include "model/modelpomo.h" #include "utils/timeutil.h" -#include "model/modelfactorymixlen.h" +#include "model/modelfactory.h" #include "tree/phylosupertreeplen.h" #include "tree/phylosupertreeunlinked.h" @@ -757,26 +757,25 @@ string computeFastMLTree(Params ¶ms, Alignment *aln, if (aln->isSuperAlignment()) { SuperAlignment *saln = (SuperAlignment*)aln; - if (params.partition_type == TOPO_UNLINKED) + if (params.partition_type == TOPO_UNLINKED) { iqtree = new PhyloSuperTreeUnlinked(saln); - else if (params.partition_type == BRLEN_OPTIMIZE) + } else if (params.partition_type == BRLEN_OPTIMIZE) { iqtree = new PhyloSuperTree(saln); - else + } else { iqtree = new PhyloSuperTreePlen(saln, brlen_type); + } for (int part = 0; part != subst_names.size(); part++) { saved_model_names.push_back(saln->partitions[part]->model_name); saln->partitions[part]->model_name = subst_names[part] + rate_names[part]; } } else if (posRateHeterotachy(rate_names[0]) != string::npos) { - iqtree = new PhyloTreeMixlen(aln, 0); + iqtree = new PhyloTreeMixlen(aln); } else { iqtree = new IQTree(aln); } - if (params.constraint_tree_file) { iqtree->constraintTree.readConstraint(params.constraint_tree_file, aln->getSeqNames()); } - if ((params.start_tree == STT_PLL_PARSIMONY || params.start_tree == STT_RANDOM_TREE || params.pll) && !iqtree->isInitializedPLL()) { /* Initialized all data structure for PLL*/ iqtree->initializePLL(params); @@ -819,9 +818,9 @@ string computeFastMLTree(Params ¶ms, Alignment *aln, // disable thorough I+G optimization params.opt_gammai = false; initTree = iqtree->optimizeModelParameters(false, params.modelEps*50.0); - if (iqtree->isMixlen()) - initTree = ((ModelFactoryMixlen*)iqtree->getModelFactory())->sortClassesByTreeLength(); - + if (iqtree->isMixlen()) { + initTree = iqtree->getModelFactory()->sortClassesByTreeLength(); + } // do quick NNI search if (params.start_tree != STT_USER_TREE) { cout << "Perform nearest neighbor interchange..." << endl; @@ -1936,27 +1935,30 @@ string CandidateModel::evaluate(Params ¶ms, IQTree *iqtree = nullptr; if (in_aln->isSuperAlignment()) { SuperAlignment *saln = (SuperAlignment*)in_aln; - if (params.partition_type == BRLEN_OPTIMIZE) + if (params.partition_type == BRLEN_OPTIMIZE) { iqtree = new PhyloSuperTree(saln); - else + } else { iqtree = new PhyloSuperTreePlen(saln, brlen_type); + } StrVector subst_names; StrVector rate_names; convert_string_vec(subst_name.c_str(), subst_names); convert_string_vec(rate_name.c_str(), rate_names); ASSERT(subst_names.size() == rate_names.size()); - for (int part = 0; part != subst_names.size(); part++) + for (int part = 0; part != subst_names.size(); part++) { saln->partitions[part]->model_name = subst_names[part]+rate_names[part]; - } else if (posRateHeterotachy(getName()) != string::npos) - iqtree = new PhyloTreeMixlen(in_aln, 0); - else + } + } else if (posRateHeterotachy(getName()) != string::npos) { + iqtree = new PhyloTreeMixlen(in_aln); + } else { iqtree = new IQTree(in_aln); + } iqtree->setParams(¶ms); iqtree->setLikelihoodKernel(params.SSE); iqtree->optimize_by_newton = params.optimize_by_newton; iqtree->setNumThreads(num_threads); - iqtree->setCheckpoint(&in_model_info); + #ifdef _OPENMP #pragma omp critical #endif @@ -7435,18 +7437,15 @@ void runMixtureFinder(Params ¶ms, IQTree* &iqtree, ModelCheckpoint &model_in // create a new IQTree object for this mixture model - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy(aln->model_name); - if (params.num_mixlen > 1) { - new_iqtree = new PhyloTreeMixlen(aln, params.num_mixlen); - } else if (pos != string::npos) { - new_iqtree = new PhyloTreeMixlen(aln, 0); + if (posRateHeterotachy(aln->model_name) != string::npos) { + new_iqtree = new PhyloTreeMixlen(aln); } else { new_iqtree = new IQTree(aln); } new_iqtree->setCheckpoint(iqtree->getCheckpoint()); - if (!iqtree->constraintTree.empty()) + if (!iqtree->constraintTree.empty()) { new_iqtree->constraintTree.readConstraint(iqtree->constraintTree); + } new_iqtree->removed_seqs = iqtree->removed_seqs; new_iqtree->twin_seqs = iqtree->twin_seqs; if (params.start_tree == STT_PLL_PARSIMONY || params.start_tree == STT_RANDOM_TREE || params.pll) { diff --git a/model/CMakeLists.txt b/model/CMakeLists.txt index d1284e50b..74eef08ea 100644 --- a/model/CMakeLists.txt +++ b/model/CMakeLists.txt @@ -29,7 +29,6 @@ modelunrest.cpp modelunrest.h modelliemarkov.cpp modelliemarkov.h modelpomo.cpp modelpomo.h modelpomomixture.cpp modelpomomixture.h -modelfactorymixlen.cpp modelfactorymixlen.h modelhmm.cpp modelhmm.h modelhmmgm.cpp modelhmmgm.h modelhmmtm.cpp modelhmmtm.h diff --git a/model/modelfactory.cpp b/model/modelfactory.cpp index a45f1a842..689bffc9f 100644 --- a/model/modelfactory.cpp +++ b/model/modelfactory.cpp @@ -41,6 +41,7 @@ #include "rateheterotachyinvar.h" //#include "ngs.h" #include +#include "tree/phylotreemixlen.h" #include "utils/timeutil.h" #include "nclextra/myreader.h" #include @@ -1041,7 +1042,6 @@ ModelFactory::ModelFactory(Params ¶ms, string &model_name, PhyloTree *tree, } catch (const char* str) { outError(str); } - } void ModelFactory::setCheckpoint(Checkpoint *checkpoint) { @@ -1234,11 +1234,11 @@ void ModelFactory::initFromClassMinusOne(double init_weight) { } int ModelFactory::getNParameters(int brlen_type) { - int df = model->getNDim() + model->getNDimFreq() + site_rate->getNDim() + - site_rate->getTree()->getNBranchParameters(brlen_type); - + int df = model->getNDim() + model->getNDimFreq() + + site_rate->getNDim() + site_rate->getTree()->getNBranchParameters(brlen_type); return df; } + /* double ModelFactory::initGTRGammaIParameters(RateHeterogeneity *rate, ModelSubst *model, double initAlpha, double initPInvar, double *initRates, double *initStateFreqs) { @@ -1554,13 +1554,13 @@ double ModelFactory::optimizeParameters(int fixed_len, bool write_info, double logl_epsilon, double gradient_epsilon) { ASSERT(model); ASSERT(site_rate); - -// double defaultEpsilon = logl_epsilon; - - double begin_time = getRealTime(); - double cur_lh; PhyloTree *tree = site_rate->getTree(); ASSERT(tree); + if (tree->isMixlen()) { + ((PhyloTreeMixlen*)tree)->initializeMixlen(logl_epsilon, write_info); + } + double begin_time = getRealTime(); + double cur_lh; stopStoringTransMatrix(); @@ -1830,9 +1830,84 @@ ModelFactory::~ModelFactory() clear(); } +string ModelFactory::sortClassesByTreeLength() { + ASSERT(model); + ASSERT(site_rate); + PhyloTree *tree = site_rate->getTree(); + ASSERT(tree); + if (!tree->isMixlen()) { + return tree->getTreeString(); + } + DoubleVector brlen; + tree->saveBranchLengths(brlen); + size_t nmixlen = tree->getNMixlen(); + ASSERT(brlen.size() == tree->branchNum * nmixlen); + // compute tree lengths + int index[nmixlen]; + double treelen[nmixlen]; + memset(treelen, 0, nmixlen*sizeof(double)); + for (size_t i = 0; i < nmixlen; ++i) { + index[i] = i; + } + for (size_t b = 0; b < brlen.size(); ++b) { + treelen[b%nmixlen] += brlen[b]; + } + // sort tree lengths and reorder mixlen classes + quicksort(treelen, 0, nmixlen-1, index); + bool sorted = true; + for (size_t i = 0; i < nmixlen; ++i) { + sorted &= (index[i] == i); + } + if (!sorted) { + double score = tree->curScore; + cout << "Reordering rate classes by tree lengths" << endl; + // reorder mixlen class branch lengths + DoubleVector sorted_brlen(brlen.size(), 0.0); + for (size_t b = 0; b < tree->branchNum; ++b) { + for (size_t i = 0; i < nmixlen; ++i) { + sorted_brlen[b*nmixlen + i] = brlen[b*nmixlen + index[i]]; + } + } + tree->restoreBranchLengths(sorted_brlen); + // reoder mixlen class weights + ASSERT(site_rate->getNRate() == nmixlen); + double prop[nmixlen]; + for (size_t i = 0; i < nmixlen; ++i) { + prop[i] = site_rate->getProp(index[i]); + } + for (size_t i = 0; i < nmixlen; ++i) { + site_rate->setProp(i, prop[i]); + } + // reorder fused mixture models + if (fused_mix_rate) { + ASSERT(model->getNMixtures() == nmixlen); + ModelSubst *models[nmixlen]; + for (size_t i = 0; i < nmixlen; ++i) { + ASSERT(model->getMixtureWeight(i) == 1.0); + models[i] = model->getMixtureClass(index[i]); + } + for (size_t i = 0; i < nmixlen; ++i) { + model->setMixtureClass(i, models[i]); + } + int num_states = model->num_states; + for (size_t i = 0; i < nmixlen; ++i) { + ((ModelMarkov*)model->getMixtureClass(i))->setEigenvalues(&model->getEigenvalues()[i*num_states]); + ((ModelMarkov*)model->getMixtureClass(i))->setEigenvectors(&model->getEigenvectors()[i*num_states*num_states]); + ((ModelMarkov*)model->getMixtureClass(i))->setInverseEigenvectors(&model->getInverseEigenvectors()[i*num_states*num_states]); + ((ModelMarkov*)model->getMixtureClass(i))->setInverseEigenvectorsTransposed(&model->getInverseEigenvectorsTransposed()[i*num_states*num_states]); + } + model->decomposeRateMatrix(); + site_rate->writeInfo(cout); + } + tree->clearAllPartialLH(); + ASSERT(fabs(score - tree->computeLikelihood()) < 0.1); + } + return tree->getTreeString(); +} + /************* FOLLOWING SERVE FOR JOINT OPTIMIZATION OF MODEL AND RATE PARAMETERS *******/ -int ModelFactory::getNDim() -{ + +int ModelFactory::getNDim() { return model->getNDim() + site_rate->getNDim(); } diff --git a/model/modelfactory.h b/model/modelfactory.h index 8b29cb062..dc9fcf205 100644 --- a/model/modelfactory.h +++ b/model/modelfactory.h @@ -266,7 +266,11 @@ class ModelFactory : public unordered_map, public Optimization, pu */ double optimizeParametersOnly(int num_steps, double gradient_epsilon, double cur_logl); - /************* FOLLOWING FUNCTIONS SERVE FOR JOINT OPTIMIZATION OF MODEL AND RATE PARAMETERS *******/ + /** + * Sort mixlen classes in ascending order of tree lengths + * @return Tree string with sorted branch lengths + */ + string sortClassesByTreeLength(); /** * TRUE to optimize all parameters simultaneously, default: FALSE diff --git a/model/modelfactorymixlen.cpp b/model/modelfactorymixlen.cpp deleted file mode 100644 index 79998d794..000000000 --- a/model/modelfactorymixlen.cpp +++ /dev/null @@ -1,158 +0,0 @@ -/* - * modelfactorymixlen.cpp - * - * Created on: Sep 2, 2015 - * Author: minh - */ - -#include "tree/phylotreemixlen.h" -#include "utils/timeutil.h" -#include "utils/tools.h" -#include "model/modelfactorymixlen.h" -#include "model/modelmarkov.h" -#include "model/modelmixture.h" -#include "rateheterotachy.h" - -ModelFactoryMixlen::ModelFactoryMixlen(Params ¶ms, string &model_name, PhyloTree *tree, ModelsBlock *models_block) : - ModelFactory(params, model_name, tree, models_block) { - if (!tree->isMixlen()) { - cerr << "Please add '-mixlen " << site_rate->getNRate() << "' option into the command line" << endl; - outError("Sorry for the inconvience, please rerun IQ-TREE with option above"); - } - if (tree->getMixlen() != site_rate->getNRate()) { - ((PhyloTreeMixlen*)tree)->setMixlen(site_rate->getNRate()); -// outError("#heterotachy classes and #mixture branch lengths do not match"); - } - if (fused_mix_rate) { - // fix the rate of heterotachy -// RateHeterotachy *hrate = (RateHeterotachy*)site_rate; -// ModelMixture *mmodel = (ModelMixture*)model; - - /* - if (site_rate->getFixParams() == 1) { - // swap the weights between model and site_rate - for (int i = 0; i < site_rate->getNRate(); i++) { - model->setMixtureWeight(i, site_rate->getProp(i)); - site_rate->setProp(i, 1.0); - } - } else { - // fix the weight of heterotachy model to 1 - site_rate->setFixParams(2); - double fix_prop = (1.0-site_rate->getPInvar()); - for (int i = 0; i < site_rate->getNRate(); i++) - site_rate->setProp(i, fix_prop); - } - */ - } -} - -double ModelFactoryMixlen::optimizeParameters(int fixed_len, bool write_info, double logl_epsilon, double gradient_epsilon) { - - PhyloTreeMixlen *tree = (PhyloTreeMixlen*)site_rate->getTree(); - ASSERT(tree); - - tree->initializeMixlen(logl_epsilon, write_info); - - double score = ModelFactory::optimizeParameters(fixed_len, write_info, logl_epsilon, gradient_epsilon); - - return score; -} - -string ModelFactoryMixlen::sortClassesByTreeLength() { - - PhyloTreeMixlen *tree = (PhyloTreeMixlen*)site_rate->getTree(); - - // now sort the classes by tree lengths - DoubleVector brlen; - tree->saveBranchLengths(brlen); - ASSERT(brlen.size() == tree->branchNum * tree->mixlen); - - // compute tree lengths - double treelen[tree->mixlen]; - int index[tree->mixlen]; - memset(treelen, 0, sizeof(double)*tree->mixlen); - int i, j; - for (i = 0; i < tree->mixlen; i++) - index[i] = i; - for (i = 0, j = 0; i < brlen.size(); i++, j++) { - if (j == tree->mixlen) j = 0; - treelen[j] += brlen[i]; - } - - // sort tree lengths and reorder branch lengths - quicksort(treelen, 0, tree->mixlen-1, index); - bool sorted = true; - for (j = 0; j < tree->mixlen; j++) - if (index[j] != j) { sorted = false; break; }; - if (!sorted) { - double score = tree->curScore; - cout << "Reordering classes by tree lengths" << endl; - DoubleVector sorted_brlen; - sorted_brlen.resize(brlen.size()); - for (i = 0; i < tree->branchNum; i++) - for (j = 0; j < tree->mixlen; j++) - sorted_brlen[i*tree->mixlen + j] = brlen[i*tree->mixlen + index[j]]; - tree->restoreBranchLengths(sorted_brlen); - - ASSERT(tree->mixlen == site_rate->getNRate()); - // reoder class weights - double prop[site_rate->getNRate()]; - for (j = 0; j < site_rate->getNRate(); j++) - prop[j] = site_rate->getProp(index[j]); - for (j = 0; j < site_rate->getNRate(); j++) - site_rate->setProp(j, prop[j]); - - // reorder mixture models - if (fused_mix_rate) { - ASSERT(model->getNMixtures() == site_rate->getNRate()); -// ModelMixture *mixmodel = (ModelMixture*)model; - int nmix = model->getNMixtures(); - ModelSubst *models[nmix]; - for (j = 0; j < nmix; j++) - models[j] = model->getMixtureClass(index[j]); - for (j = 0; j < nmix; j++) - model->setMixtureClass(j, models[j]); - - for (j = 0; j < site_rate->getNRate(); j++) - prop[j] = model->getMixtureWeight(index[j]); - for (j = 0; j < site_rate->getNRate(); j++) - model->setMixtureWeight(j, prop[j]); - - // assigning memory for individual models - int m = 0; - int num_states = model->num_states; - for (m = 0; m < nmix; m++) { - ((ModelMarkov*)model->getMixtureClass(m))->setEigenvalues(&model->getEigenvalues()[m*num_states]); - ((ModelMarkov*)model->getMixtureClass(m))->setEigenvectors(&model->getEigenvectors()[m*num_states*num_states]); - ((ModelMarkov*)model->getMixtureClass(m))->setInverseEigenvectors(&model->getInverseEigenvectors()[m*num_states*num_states]); - ((ModelMarkov*)model->getMixtureClass(m))->setInverseEigenvectorsTransposed(&model->getInverseEigenvectorsTransposed()[m*num_states*num_states]); - } - model->decomposeRateMatrix(); - site_rate->writeInfo(cout); - } - tree->clearAllPartialLH(); - ASSERT(fabs(score - tree->computeLikelihood()) < 0.1); - } - - // update relative_rate - /* - double sum = 0.0; - for (j = 0; j < tree->mixlen; j++) - sum += treelen[j]; - sum = tree->mixlen/sum; - for (j = 0; j < tree->mixlen; j++) - tree->relative_rate->setRate(j, treelen[j]*sum); - */ - return tree->getTreeString(); -} - -int ModelFactoryMixlen::getNParameters(int brlen_type) { - int df = ModelFactory::getNParameters(brlen_type); - if (brlen_type == BRLEN_OPTIMIZE) { - df += site_rate->phylo_tree->branchNum * (site_rate->phylo_tree->getMixlen() - 1); - } - else if (brlen_type == BRLEN_SCALE) { - df += (site_rate->phylo_tree->getMixlen() - 1); - } - return df; -} diff --git a/model/modelfactorymixlen.h b/model/modelfactorymixlen.h deleted file mode 100644 index 71603bdbc..000000000 --- a/model/modelfactorymixlen.h +++ /dev/null @@ -1,46 +0,0 @@ -/* - * modelfactorymixlen.h - * - * Created on: Sep 2, 2015 - * Author: minh - */ - - -#include "modelfactory.h" - -class ModelFactoryMixlen : public ModelFactory { - -public: - - /** - constructor - create substitution model with possible rate heterogeneity. Create proper class objects - for two variables: model and site_rate. It takes the following field of params into account: - model_name, num_rate_cats, freq_type, store_trans_matrix - @param params program parameters - @param tree associated phylogenetic tree - */ - ModelFactoryMixlen(Params ¶ms, string &model_name, PhyloTree *tree, ModelsBlock *models_block); - - /** - optimize model parameters and tree branch lengths - @param fixed_len TRUE to fix branch lengths, default is false - @return the best likelihood - */ - virtual double optimizeParameters(int fixed_len = BRLEN_OPTIMIZE, bool write_info = true, - double logl_epsilon = 0.1, double gradient_epsilon = 0.0001); - - /** - sort classes in ascending order of tree lengths - @return tree string with sorted branch lengths - */ - string sortClassesByTreeLength(); - - - /** - * @param brlen_type either BRLEN_OPTIMIZE, BRLEN_FIX or BRLEN_SCALE - * @return #parameters of the model + # branches - */ - virtual int getNParameters(int brlen_type); - -}; diff --git a/model/rateheterotachy.cpp b/model/rateheterotachy.cpp index 38e8994f9..00c43e0cd 100644 --- a/model/rateheterotachy.cpp +++ b/model/rateheterotachy.cpp @@ -10,6 +10,7 @@ #include "rateheterotachy.h" RateHeterotachy::RateHeterotachy(int ncat, string params, PhyloTree *tree) : RateHeterogeneity() { + ASSERT(tree->isMixlen()); phylo_tree = tree; prop = nullptr; fix_params = 0; diff --git a/model/rateheterotachy.h b/model/rateheterotachy.h index 77b87245b..024eacdcc 100644 --- a/model/rateheterotachy.h +++ b/model/rateheterotachy.h @@ -18,9 +18,6 @@ class PhyloTree; rate-heterotachy model, allowing for mixed branch lengths */ class RateHeterotachy: virtual public RateHeterogeneity { - - friend class ModelFactoryMixlen; - public: /** constructor diff --git a/simulator/alisimulator.cpp b/simulator/alisimulator.cpp index ce44c7c1d..14d55b59e 100644 --- a/simulator/alisimulator.cpp +++ b/simulator/alisimulator.cpp @@ -75,10 +75,9 @@ AliSimulator::AliSimulator(Params *input_params, IQTree *iq_tree, int expected_n /** * initialize an IQTree instance from input file */ -void AliSimulator::initializeIQTreeFromTreeFile() -{ - // handle the case with partition models +void AliSimulator::initializeIQTreeFromTreeFile() { if (params->partition_file) { + /* partitioned model case */ // initilize partition alignments Alignment *aln; if (params->partition_type == TOPO_UNLINKED) @@ -130,12 +129,16 @@ void AliSimulator::initializeIQTreeFromTreeFile() cout<<" Loading partition trees one by one. Each tree should be specified in a single line in the input tree file."<at(i); + // load the tree for the current partition + IQTree *¤t_tree = (IQTree*&)((PhyloSuperTree*)tree)->at(i); + if (posRateHeterotachy(current_tree->aln->model_name) != string::npos) { + // heterotachy rate model -> reset the tree as PhyloTreeMixlen + IQTree *new_tree = new PhyloTreeMixlen(current_tree->aln); + delete current_tree; + current_tree = new_tree; + } bool is_rooted = false; current_tree->readTree(params->user_file, is_rooted, tree_line_index); - // update the alignment for the current partition int expected_num_states_current_tree = current_tree->aln->getNSite(); ASSERT(expected_num_states_current_tree == 0); @@ -151,29 +154,8 @@ void AliSimulator::initializeIQTreeFromTreeFile() pat.frequency = expected_num_states_current_tree; pat.flag = PAT_INVARIANT; current_tree->aln->addPattern(pat); - // initialize the model for the current partition initializeModel(current_tree, current_tree->aln->model_name); - - // if a Heterotachy model is used -> re-read the PhyloTreeMixlen from file - if (current_tree->getRate()->isHeterotachy()) - { - // initialize a new PhyloTreeMixlen - IQTree* new_tree = new PhyloTreeMixlen(current_tree->aln, current_tree->getRate()->getNRate()); - - // delete the old tree - delete current_tree; - - // set the new PhyloTreeMixlen to the new tree - current_tree = new_tree; - - // re-load the tree/branch-lengths from the file - current_tree->IQTree::readTree(params->user_file, is_rooted, tree_line_index); - - // re-initialize the model - initializeModel(current_tree, current_tree->aln->model_name); - } - // set partition rate if (params->partition_type == BRLEN_SCALE) { @@ -265,41 +247,22 @@ void AliSimulator::initializeIQTreeFromTreeFile() ((PhyloSuperTree*) tree)->part_info[i].part_rate *= sum; } } - } - // other cases without partition models - else - { + } else { + /* non-partitioned model case */ // initialize tree - tree = new IQTree(); + if (posRateHeterotachy(params->model_name) != string::npos) { + tree = new PhyloTreeMixlen(); + } else { + tree = new IQTree(); + } bool is_rooted = false; tree->readTree(params->user_file, is_rooted); tree->setParams(params); - // initialize alignment tree->aln = new Alignment(); initializeAlignment(tree, params->model_name); - - // inittialize model + // initialize model initializeModel(tree, params->model_name); - - // if a Heterotachy model is used -> re-read the PhyloTreeMixlen from file - if (tree->getRate()->isHeterotachy()) - { - // initialize a new PhyloTreeMixlen - IQTree* new_tree = new PhyloTreeMixlen(tree->aln, tree->getRate()->getNRate()); - - // delete the old tree - delete tree; - - // set the new PhyloTreeMixlen to the new tree - tree = new_tree; - - // re-load the tree/branch-lengths from the file - tree->IQTree::readTree(params->user_file, is_rooted); - - // re-initialize the model - initializeModel(tree, params->model_name); - } } } @@ -2488,7 +2451,13 @@ void AliSimulator::branchSpecificEvolutionMasterThread(int sequence_length, doub string model_full_name = (*it)->attributes["model"]; // convert separator from "/" to "," std::replace(model_full_name.begin(), model_full_name.end(), '/', ','); - IQTree *tmp_tree = new IQTree(); + // initialize tree + IQTree *tmp_tree; + if (tree->getRate()->isHeterotachy()) { + tmp_tree = new PhyloTreeMixlen(); + } else { + tmp_tree = new IQTree(); + } tmp_tree->copyPhyloTree(tree, true); initializeModel(tmp_tree, model_full_name); diff --git a/tree/iqtree.cpp b/tree/iqtree.cpp index 05064355f..47de30d76 100644 --- a/tree/iqtree.cpp +++ b/tree/iqtree.cpp @@ -24,13 +24,11 @@ #include "phylosupertree.h" #include "phylosupertreeplen.h" #include "model/partitionmodelplen.h" -#include "model/modelfactorymixlen.h" +#include "model/modelfactory.h" #include "mexttree.h" #include "utils/timeutil.h" #include "model/modelmarkov.h" #include "model/rategamma.h" -//#include "phylotreemixlen.h" -//#include "model/modelfactorymixlen.h" #include #include "utils/tools.h" #include "utils/MPIHelper.h" @@ -2374,7 +2372,8 @@ string IQTree::ensureModelParametersAreSet(double initEpsilon) { initTree = getTreeString(); cout << "CHECKPOINT: Model parameters restored, LogL: " << getCurScore() << endl; } else { - // for mixtureFinder, verify whether the likelihood is the same as the best likelihood obtained in mixtureFinder + // for mixtureFinder, verify that the likelihood is the same as + // the best likelihood obtained in mixtureFinder double mixFinderLogL; if (getCheckpoint()->get("MixFinderLogL", mixFinderLogL)) { double allowableDiff = 0.01; @@ -2382,10 +2381,9 @@ string IQTree::ensureModelParametersAreSet(double initEpsilon) { ASSERT(fabs(currLogL - mixFinderLogL) < allowableDiff); getCheckpoint()->eraseKeyPrefix("MixFinderLogL"); } - initTree = optimizeModelParameters(true, initEpsilon); if (isMixlen()) { - initTree = ((ModelFactoryMixlen*)getModelFactory())->sortClassesByTreeLength(); + initTree = getModelFactory()->sortClassesByTreeLength(); } saveCheckpoint(); getModelFactory()->saveCheckpoint(); @@ -2965,20 +2963,15 @@ void IQTree::refineBootTrees() { // create bootstrap tree IQTree *boot_tree; - if (aln->isSuperAlignment()){ - if(params->partition_type != BRLEN_OPTIMIZE){ - boot_tree = new PhyloSuperTreePlen((SuperAlignment*) bootstrap_alignment, (PhyloSuperTree*) this); + if (aln->isSuperAlignment()) { + if (params->partition_type != BRLEN_OPTIMIZE) { + boot_tree = new PhyloSuperTreePlen((SuperAlignment*)bootstrap_alignment, (PhyloSuperTree*)this); } else { - boot_tree = new PhyloSuperTree((SuperAlignment*) bootstrap_alignment, (PhyloSuperTree*) this); + boot_tree = new PhyloSuperTree((SuperAlignment*)bootstrap_alignment, (PhyloSuperTree*)this); } } else { - // allocate heterotachy tree if neccessary - int pos = posRateHeterotachy(aln->model_name); - - if (params->num_mixlen > 1) { - boot_tree = new PhyloTreeMixlen(bootstrap_alignment, params->num_mixlen); - } else if (pos != string::npos) { - boot_tree = new PhyloTreeMixlen(bootstrap_alignment, 0); + if (posRateHeterotachy(aln->model_name) != string::npos) { + boot_tree = new PhyloTreeMixlen(bootstrap_alignment); } else { boot_tree = new IQTree(bootstrap_alignment); } @@ -4983,10 +4976,10 @@ int PhyloTree::testNumThreads() { cout << "BEST NUMBER OF THREADS: " << bestProc+1 << endl << endl; setNumThreads(bestProc+1); - - // clear the relative treelength arrays if it is GHOST model + + // clear the relative_treelen vector if it is GHOST model if (isMixlen()) { - ((PhyloTreeMixlen*)this)->clear_relative_treelen(); + ((PhyloTreeMixlen*)this)->clearRelativeTreelen(); } return bestProc+1; diff --git a/tree/phylokernelnew.h b/tree/phylokernelnew.h index 6f2704c60..8969beb24 100644 --- a/tree/phylokernelnew.h +++ b/tree/phylokernelnew.h @@ -2082,7 +2082,7 @@ void PhyloTree::computeLikelihoodBufferGenericSIMD(PhyloNeighbor *dad_branch, Ph // reserve 3*block for computeLikelihoodDerv double *buffer_partial_lh_ptr = buffer_partial_lh + 3*get_safe_upper_limit(block); if (isMixlen()) { - size_t nmix = getMixlen(); + size_t nmix = getNMixlen(); buffer_partial_lh_ptr += nmix*(nmix+1)*VectorClass::size() + (nmix+3)*nmix*VectorClass::size()*num_packets; } @@ -2367,7 +2367,7 @@ void PhyloTree::computeLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch, Phyl VectorClass *all_dfvec = nullptr; VectorClass *all_ddfvec = nullptr; - size_t nmixlen = getMixlen(), nmixlen2 = nmixlen*nmixlen; + size_t nmixlen = getNMixlen(), nmixlen2 = nmixlen*nmixlen; if (isMixlen()) { ASSERT(nmixlen == ncat); all_dfvec = (VectorClass*)buffer_partial_lh_ptr; diff --git a/tree/phylosupertree.cpp b/tree/phylosupertree.cpp index a4b496ea4..098b0dce7 100644 --- a/tree/phylosupertree.cpp +++ b/tree/phylosupertree.cpp @@ -299,36 +299,31 @@ void PhyloSuperTree::readTreeString(const string &tree_string) { } - -/** - * save branch lengths into a vector - */ -void PhyloSuperTree::saveBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { - ASSERT(getMixlen() == 1); // supertree and treemixlen not allowed together - int totalBranchNum = branchNum * getMixlen(); - iterator it; - for (it = begin(); it != end(); it++) { - totalBranchNum += (*it)->branchNum * (*it)->getMixlen(); - } - lenvec.resize(startid + totalBranchNum); - - PhyloTree::saveBranchLengths(lenvec, startid); - startid += branchNum * getMixlen(); - for (iterator it = begin(); it != end(); it++) { - (*it)->saveBranchLengths(lenvec, startid); - startid += (*it)->branchNum * (*it)->getMixlen(); - } +void PhyloSuperTree::saveBranchLengths(DoubleVector &lenvec, size_t startid, PhyloNode *, PhyloNode *) { + if (!isSuperTreeUnlinked()) { + PhyloTree::saveBranchLengths(lenvec, startid); + startid += branchNum * getNMixlen(); + } + size_t min_size = startid; + for (iterator it = begin(); it != end(); ++it) { + min_size += (*it)->branchNum * (*it)->getNMixlen(); + } + lenvec.resize(max(lenvec.size(), min_size), 0.0); + for (iterator it = begin(); it != end(); ++it) { + (*it)->saveBranchLengths(lenvec, startid); + startid += (*it)->branchNum * (*it)->getNMixlen(); + } } -/** - * restore branch lengths from a vector previously called with saveBranchLengths - */ -void PhyloSuperTree::restoreBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { - PhyloTree::restoreBranchLengths(lenvec, startid); - startid += branchNum * getMixlen(); - for (iterator it = begin(); it != end(); it++) { - (*it)->restoreBranchLengths(lenvec, startid); - startid += (*it)->branchNum * (*it)->getMixlen(); - } + +void PhyloSuperTree::restoreBranchLengths(DoubleVector &lenvec, size_t startid, PhyloNode *, PhyloNode *) { + if (!isSuperTreeUnlinked()) { + PhyloTree::restoreBranchLengths(lenvec, startid); + startid += branchNum * getNMixlen(); + } + for (iterator it = begin(); it != end(); ++it) { + (*it)->restoreBranchLengths(lenvec, startid); + startid += (*it)->branchNum * (*it)->getNMixlen(); + } } int PhyloSuperTree::collapseInternalBranches(Node *node, Node *dad, double threshold) { diff --git a/tree/phylosupertree.h b/tree/phylosupertree.h index 168d76892..593a38a49 100644 --- a/tree/phylosupertree.h +++ b/tree/phylosupertree.h @@ -159,13 +159,17 @@ class PhyloSuperTree : public IQTree, public vector virtual void readTreeString(const string &tree_string) override; /** - * save branch lengths into a vector + * Save branch lengths into a vector */ - virtual void saveBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; + virtual void saveBranchLengths(DoubleVector &lenvec, size_t startid = 0, + PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; + /** - * restore branch lengths from a vector previously called with saveBranchLengths + * Restore branch lengths from a vector that + * was previously filled by saveBranchLengths() */ - virtual void restoreBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; + virtual void restoreBranchLengths(DoubleVector &lenvec, size_t startid = 0, + PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; /** Collapse all internal branches with length <= threshold diff --git a/tree/phylosupertreeunlinked.cpp b/tree/phylosupertreeunlinked.cpp index 768e6dff4..683a90101 100644 --- a/tree/phylosupertreeunlinked.cpp +++ b/tree/phylosupertreeunlinked.cpp @@ -123,32 +123,6 @@ void PhyloSuperTreeUnlinked::restoreCheckpoint() { } } -/** - * save branch lengths into a vector - */ -void PhyloSuperTreeUnlinked::saveBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { - int totalBranchNum = 0; - iterator it; - for (it = begin(); it != end(); it++) { - totalBranchNum += (*it)->branchNum * (*it)->getMixlen(); - } - lenvec.resize(startid + totalBranchNum); - - for (iterator it = begin(); it != end(); it++) { - (*it)->saveBranchLengths(lenvec, startid); - startid += (*it)->branchNum * (*it)->getMixlen(); - } -} -/** - * restore branch lengths from a vector previously called with saveBranchLengths - */ -void PhyloSuperTreeUnlinked::restoreBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { - for (iterator it = begin(); it != end(); it++) { - (*it)->restoreBranchLengths(lenvec, startid); - startid += (*it)->branchNum * (*it)->getMixlen(); - } -} - void PhyloSuperTreeUnlinked::setRootNode(const char *my_root, bool multi_taxa) { // DOES NOTHING } diff --git a/tree/phylosupertreeunlinked.h b/tree/phylosupertreeunlinked.h index f92136a48..07e1e1b67 100644 --- a/tree/phylosupertreeunlinked.h +++ b/tree/phylosupertreeunlinked.h @@ -94,16 +94,7 @@ class PhyloSuperTreeUnlinked : public PhyloSuperTree { restore object from the checkpoint */ virtual void restoreCheckpoint() override; - - /** - * save branch lengths into a vector - */ - virtual void saveBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; - /** - * restore branch lengths from a vector previously called with saveBranchLengths - */ - virtual void restoreBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr) override; - + /** set the root by name @param my_root root node name @param multi_taxa TRUE if my_root is a comma-separated list of nodes diff --git a/tree/phylotree.cpp b/tree/phylotree.cpp index f5b4ed4ff..8374d508f 100644 --- a/tree/phylotree.cpp +++ b/tree/phylotree.cpp @@ -376,14 +376,11 @@ void PhyloTree::copyPhyloTree(PhyloTree *tree, bool borrowSummary) { } } -void PhyloTree::copyPhyloTreeMixlen(PhyloTree *tree, int mix, bool borrowSummary) { - if (tree->isMixlen()) { - ((PhyloTreeMixlen*)tree)->cur_mixture = mix; - } +void PhyloTree::copyPhyloTreeMixlen(PhyloTree *tree, int c, bool borrowSummary) { + ASSERT(tree->getCurMixture() == -1); + tree->setCurMixture(c); copyPhyloTree(tree, borrowSummary); - if (tree->isMixlen()) { - ((PhyloTreeMixlen*)tree)->cur_mixture = -1; - } + tree->setCurMixture(-1); } #define FAST_NAME_CHECK 1 @@ -730,31 +727,34 @@ string PhyloTree::getModelNameParams(bool show_fixed_params) { return name; } -void PhyloTree::saveBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { +void PhyloTree::saveBranchLengths(DoubleVector &lenvec, size_t startid, PhyloNode *node, PhyloNode *dad) { if (!node) { - node = (PhyloNode*) root; - ASSERT(branchNum == nodeNum-1); - if (lenvec.empty()) lenvec.resize((branchNum*getMixlen()) + startid); + node = (PhyloNode*)root; + ASSERT(branchNum == nodeNum - 1); + size_t min_size = startid + (branchNum * getNMixlen()); + lenvec.resize(max(lenvec.size(), min_size), 0.0); } - FOR_NEIGHBOR_IT(node, dad, it){ - (*it)->getLength(lenvec, (*it)->id*getMixlen() + startid); - PhyloTree::saveBranchLengths(lenvec, startid, (PhyloNode*) (*it)->node, node); + FOR_NEIGHBOR_IT(node, dad, it) { + size_t this_startid = startid + ((*it)->id * getNMixlen()); + (*it)->getLength(lenvec, this_startid); + PhyloTree::saveBranchLengths(lenvec, startid, (PhyloNode*)(*it)->node, node); } } -void PhyloTree::restoreBranchLengths(DoubleVector &lenvec, int startid, PhyloNode *node, PhyloNode *dad) { +void PhyloTree::restoreBranchLengths(DoubleVector &lenvec, size_t startid, PhyloNode *node, PhyloNode *dad) { if (!node) { - node = (PhyloNode*) root; - ASSERT(!lenvec.empty()); + node = (PhyloNode*)root; + size_t min_size = startid + (branchNum * getNMixlen()); + ASSERT(lenvec.size() >= min_size); } - FOR_NEIGHBOR_IT(node, dad, it){ - (*it)->setLength(lenvec, (*it)->id*getMixlen() + startid, getMixlen()); - (*it)->node->findNeighbor(node)->setLength(lenvec, (*it)->id*getMixlen() + startid, getMixlen()); - PhyloTree::restoreBranchLengths(lenvec, startid, (PhyloNode*) (*it)->node, node); + FOR_NEIGHBOR_IT(node, dad, it) { + size_t this_startid = startid + ((*it)->id * getNMixlen()); + (*it)->setLength(lenvec, this_startid, getNMixlen()); + (*it)->node->findNeighbor(node)->setLength(lenvec, this_startid, getNMixlen()); + PhyloTree::restoreBranchLengths(lenvec, startid, (PhyloNode*)(*it)->node, node); } } - /**************************************************************************** Parsimony function ****************************************************************************/ @@ -893,7 +893,7 @@ size_t PhyloTree::getBufferPartialLhSize() { buffer_size += get_safe_upper_limit(3*block*model->num_states); if (isMixlen()) { - size_t nmix = max(getMixlen(), getRate()->getNRate()); + size_t nmix = max(getNMixlen(), getRate()->getNRate()); buffer_size += nmix*(nmix+1)*VECTOR_SIZE + (nmix+3)*nmix*VECTOR_SIZE*num_packets; } return buffer_size; @@ -2596,29 +2596,31 @@ void PhyloTree::optimizePatternRates(DoubleVector &pattern_rates) { } int PhyloTree::getNBranchParameters(int brlen_type) { - if (params->fixed_branch_length || brlen_type == BRLEN_FIX) + if (params->fixed_branch_length) { return 0; - + } int df = 0; - - if (brlen_type == BRLEN_OPTIMIZE) { - df = branchNum - (int)rooted; - // If model is Lie-Markov, and is in fact time reversible, one of the - // degrees of freedom is illusary. (Of the two edges coming from the - // root, only sum of their lenghts affects likelihood.) - // So correct for this. Without this correction, K2P and RY2.2b - // would not be synonymous, for example. - -// string className(typeid(*model).name()); -// if (className.find("ModelLieMarkov")!=string::npos && model->isReversible()) -// df--; - - // BQM 2017-04-28, alternatively, check if there is a virtual_root and model is reversible - if (rooted && model && model->isReversible()) - df--; - - } else if (brlen_type == BRLEN_SCALE) + if (brlen_type == BRLEN_SCALE) { + // one scaler for all branches and mixlen categories df = 1; + } else if (brlen_type == BRLEN_OPTIMIZE) { + df = branchNum; + if (rooted) { + // a rooted tree always has an additional root leaf, + // its branch length does not affect the likelihood + df--; + if (model && model->isReversible()) { + // if a reversible model is used with a rooted tree, then + // of the two edges coming from the root, only the sum + // of their lenghts affects the likelihood + // Example: for Lie-Markov models, without this correction + // K2P and RY2.2b would not be synonymous + df--; + } + } + // each mixlen category has its own set of branch lengths + df *= getNMixlen(); + } return df; } @@ -4325,7 +4327,7 @@ NNIMove PhyloTree::getBestNNIForBran(PhyloNode *node1, PhyloNode *node2, NNIMove ((PhyloNeighbor*)*saved_it[id])->size = 0; } - int nni5_num_eval = max(params->nni5_num_eval, getMixlen()); + int nni5_num_eval = max(params->nni5_num_eval, getNMixlen()); for (int step = 0; step < nni5_num_eval; step++) { diff --git a/tree/phylotree.h b/tree/phylotree.h index 2b752e4b1..e096445a7 100644 --- a/tree/phylotree.h +++ b/tree/phylotree.h @@ -341,21 +341,18 @@ Phylogenetic Tree class @author BUI Quang Minh, Steffen Klaere, Arndt von Haeseler */ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { - - friend class PhyloSuperTree; - friend class PhyloSuperTreePlen; - friend class RateGamma; - friend class RateGammaInvar; - friend class RateKategory; + friend class PhyloSuperTree; + friend class PhyloSuperTreePlen; + friend class RateGamma; + friend class RateGammaInvar; + friend class RateKategory; friend class ModelMixture; friend class RateFree; friend class RateHeterotachy; friend class PhyloTreeMixlen; - friend class ModelFactoryMixlen; friend class MemSlotVector; friend class ModelFactory; friend class IQTreeMix; - public: /** default constructor ( everything is initialized to nullptr) @@ -439,22 +436,21 @@ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { void copyConstraintTree(MTree *tree, IntVector &taxon_order, int *rand_stream); /** - copy the phylogenetic tree structure into this tree, designed specifically for PhyloTree. - So there is some distinction with copyTree. - @param tree the tree to copy - @param borrowSummary true if the alignment summary of the other tree is to be "borrowed" + * Copy the phylogenetic tree structure into this tree. + * Designed specifically for copying PhyloTree objects + * @param tree The tree to copy + * @param borrowSummary TRUE to use the alignment summary of tree */ void copyPhyloTree(PhyloTree *tree, bool borrowSummary); /** - copy the phylogenetic tree structure into this tree, designed specifically for PhyloTree. - So there is some distinction with copyTree. - @param tree the tree to copy - @param mix mixture ID of branch lengths - @param borrowSummary true if the alignment summary of the other tree is to be "borrowed" + * Copy the phylogenetic tree structure into this tree. + * Designed specifically for copying PhyloTree objects + * @param tree The tree to copy + * @param c The mixlen category to copy + * @param borrowSummary TRUE to use the alignment summary of tree */ - virtual void copyPhyloTreeMixlen(PhyloTree *tree, int mix, bool borrowSummary); - + virtual void copyPhyloTreeMixlen(PhyloTree *tree, int c, bool borrowSummary); /** Set the alignment, important to compute parsimony or likelihood score @@ -540,11 +536,6 @@ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { return false; } - /** - @return true if this is a tree with mixture branch lengths, default: false - */ - virtual bool isMixlen() { return false; } - /** @return true if this is a mixture of trees, default: false */ @@ -556,9 +547,15 @@ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { virtual bool isHMM() { return false; } /** - @return number of mixture branch lengths, default: 1 - */ - virtual int getMixlen() { return 1; } + * Functions for the mixlen stuff + */ + virtual bool isMixlen() const { return false; } + + virtual int getNMixlen() const { return 1; } + + virtual int getCurMixture() const { return -1; } + + virtual void setCurMixture(int c) {} /** allocate a new node. Override this if you have an inherited Node class. @@ -591,13 +588,16 @@ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { } /** - * save branch lengths into a vector + * Save branch lengths into a vector */ - virtual void saveBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr); + virtual void saveBranchLengths(DoubleVector &lenvec, size_t startid = 0, + PhyloNode *node = nullptr, PhyloNode *dad = nullptr); /** - * restore branch lengths from a vector previously called with saveBranchLengths + * Restore branch lengths from a vector that + * was previously filled by saveBranchLengths() */ - virtual void restoreBranchLengths(DoubleVector &lenvec, int startid = 0, PhyloNode *node = nullptr, PhyloNode *dad = nullptr); + virtual void restoreBranchLengths(DoubleVector &lenvec, size_t startid = 0, + PhyloNode *node = nullptr, PhyloNode *dad = nullptr); /**************************************************************************** Dot product @@ -1311,9 +1311,6 @@ class PhyloTree : public MTree, public Optimization, public CheckpointFactory { template void computeLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch, PhyloNode *dad, double *df, double *ddf); - /** For Mixlen stuffs */ - virtual int getCurMixture() { return 0; } - template void computeLikelihoodDervMixlenSIMD(PhyloNeighbor *dad_branch, PhyloNode *dad, double &df, double &ddf); diff --git a/tree/phylotreemixlen.cpp b/tree/phylotreemixlen.cpp index cc2e6f45b..e5316a3d7 100644 --- a/tree/phylotreemixlen.cpp +++ b/tree/phylotreemixlen.cpp @@ -5,121 +5,70 @@ // Created by Minh Bui on 24/08/15. // // - #include "phylotreemixlen.h" #include "phylonodemixlen.h" -#include "model/modelfactorymixlen.h" -#include "model/modelmixture.h" #include "model/ratefree.h" #include "utils/MPIHelper.h" -#ifdef USE_CPPOPTLIB -#include "cppoptlib/solver/newtondescentsolver.h" -#include "cppoptlib/solver/lbfgsbsolver.h" -#endif - -PhyloTreeMixlen::PhyloTreeMixlen() : IQTree() -#ifdef USE_CPPOPTLIB -, cppoptlib::BoundedProblem() -#endif -{ - mixlen = 1; - cur_mixture = -1; -// relative_treelen = nullptr; - initializing_mixlen = false; -} +PhyloTreeMixlen::PhyloTreeMixlen() +: IQTree(), mixlen(1), cur_mixture(-1), initializing_mixlen(false) {} -PhyloTreeMixlen::PhyloTreeMixlen(Alignment *aln, int mixlen) : IQTree(aln) -#ifdef USE_CPPOPTLIB -, cppoptlib::BoundedProblem(mixlen) -#endif -{ -// cout << "Initializing heterotachy mixture branch lengths" << endl; - cur_mixture = -1; -// relative_treelen = nullptr; - initializing_mixlen = false; - setMixlen(mixlen); -} +PhyloTreeMixlen::PhyloTreeMixlen(Alignment *aln) +: IQTree(aln), mixlen(1), cur_mixture(-1), initializing_mixlen(false) {} -PhyloTreeMixlen::~PhyloTreeMixlen() { -// if (relative_treelen) -// aligned_free(relative_treelen); +void PhyloTreeMixlen::setCurMixture(int c) { + ASSERT(c >= -1 && c < (int)mixlen); + cur_mixture = c; } void PhyloTreeMixlen::startCheckpoint() { - if (mixlen > 0) - checkpoint->startStruct("PhyloTreeMixlen" + convertIntToString(getMixlen())); - else + if (mixlen > 1) { + checkpoint->startStruct("PhyloTreeMixlen" + convertIntToString(mixlen)); + } else { PhyloTree::startCheckpoint(); + } } void PhyloTreeMixlen::saveCheckpoint() { - if (mixlen > 0) { + if (mixlen > 1) { startCheckpoint(); - if (this->relative_treelen.size() > 0) { - ASSERT(mixlen == this->relative_treelen.size()); - double relative_treelen[mixlen]; - for (int i = 0; i < mixlen; i++) - relative_treelen[i] = this->relative_treelen[i]; - CKP_ARRAY_SAVE(mixlen, relative_treelen); + if (relative_treelen.size() > 0) { + ASSERT(mixlen == relative_treelen.size()); + CKP_ARRAY_SAVE(mixlen, &relative_treelen[0]); } endCheckpoint(); } IQTree::saveCheckpoint(); } -/** - restore object from the checkpoint -*/ void PhyloTreeMixlen::restoreCheckpoint() { - if (mixlen > 0) { + if (mixlen > 1) { startCheckpoint(); - double relative_treelen[mixlen]; - if (CKP_ARRAY_RESTORE(mixlen, relative_treelen)) { - this->relative_treelen.resize(mixlen); - for (int i = 0; i < mixlen; i++) - this->relative_treelen[i] = relative_treelen[i]; + DoubleVector restored_relative_treelen(mixlen, 0.0); + if (CKP_ARRAY_RESTORE(mixlen, &restored_relative_treelen[0])) { + relative_treelen = restored_relative_treelen; } endCheckpoint(); } IQTree::restoreCheckpoint(); - if (!root) { - // if not success, try to restore from PhyloTree - int orig_mixlen = mixlen; - mixlen = 0; - PhyloTree::restoreCheckpoint(); - mixlen = orig_mixlen; - } } Node* PhyloTreeMixlen::newNode(int node_id, const char* node_name) { - return (Node*) (new PhyloNodeMixlen(node_id, node_name)); + return (Node*)(new PhyloNodeMixlen(node_id, node_name)); } Node* PhyloTreeMixlen::newNode(int node_id, int node_name) { - return (Node*) (new PhyloNodeMixlen(node_id, node_name)); -} - -void PhyloTreeMixlen::setMixlen(int mixlen) { - this->mixlen = mixlen; + return (Node*)(new PhyloNodeMixlen(node_id, node_name)); } void PhyloTreeMixlen::readTreeString(const string &tree_string) { IQTree::readTreeString(tree_string); - treeLengths(relative_treelen); - if (mixlen > 0 && relative_treelen[0] == 0.0) - relative_treelen.clear(); -} - -void PhyloTreeMixlen::initializeModel(Params ¶ms, string model_name, ModelsBlock *models_block) { - try { - if (!getModelFactory()) { - setModelFactory(new ModelFactoryMixlen(params, model_name, this, models_block)); + if (mixlen > 1) { + treeLengths(relative_treelen); + if (relative_treelen[0] == 0.0) { + relative_treelen.clear(); } - } catch (string & str) { - outError(str); } - IQTree::initializeModel(params, model_name, models_block); } void PhyloTreeMixlen::treeLengths(DoubleVector &lenvec, Node *node, Node *dad) { @@ -133,6 +82,15 @@ void PhyloTreeMixlen::treeLengths(DoubleVector &lenvec, Node *node, Node *dad) { } } +void PhyloTreeMixlen::setModelFactory(ModelFactory *model_fac) { + IQTree::setModelFactory(model_fac); + if (site_rate) { + ASSERT(site_rate->isHeterotachy()); + mixlen = site_rate->getNRate(); + // clear, as we now have new weights and mixlen + relative_treelen.clear(); + } +} void PhyloTreeMixlen::initializeMixBranches(PhyloNode *node, PhyloNode *dad) { if (!node) { @@ -297,6 +255,7 @@ void PhyloTreeMixlen::initializeMixlen(double tolerance, bool write_info) { } if (((PhyloNeighborMixlen*)root->neighbors[0])->lengths.size() != mixlen) { + ASSERT(relative_treelen.size() == mixlen); // assign branch length from rate model DoubleVector saved_treelen = relative_treelen; DoubleVector lenvec; @@ -349,31 +308,6 @@ void PhyloTreeMixlen::optimizeOneBranch(PhyloNode *node1, PhyloNode *node2, bool theta_computed = false; -#ifdef USE_CPPOPTLIB - if (params->optimize_alg_mixlen.find("cppopt") != string::npos) { - //----- using cppoptlib ------// - - TVector lower_bound(mixlen), upper_bound(mixlen), variables(mixlen); - - // variables.resize(mixlen); - for (i = 0; i < mixlen; i++) { - lower_bound[i] = params->min_branch_length; - variables[i] = current_it->getLength(i); - upper_bound[i] = params->max_branch_length; - } - - setBoxConstraint(lower_bound, upper_bound); - - cppoptlib::NewtonDescentSolver solver; - // cppoptlib::LbfgsbSolver solver; - solver.minimize(*this, variables); - for (i = 0; i < mixlen; i++) { - current_it->setLength(i, variables[i]); - current_it_back->setLength(i, variables[i]); - } - } else -#endif - if (params->optimize_alg_mixlen.find("newton") != string::npos) { //----- Newton-Raphson -----// @@ -687,50 +621,6 @@ void PhyloTreeMixlen::printResultTree(string suffix) { cout << "Best tree printed to " << tree_file_name << endl; } - -/*************** Using cppoptlib for branch length optimization ***********/ - -#ifdef USE_CPPOPTLIB -double PhyloTreeMixlen::value(const TVector &x) { - double xx[mixlen+1]; - for (int i = 0; i < mixlen; i++) - xx[i+1] = x(i); - return targetFunk(xx); -} - -void PhyloTreeMixlen::gradient(const TVector &x, TVector &grad) { - int i; - double xx[mixlen]; - for (i = 0; i < mixlen; i++) - xx[i] = x(i); - double df[mixlen+1], ddf[mixlen*mixlen]; - computeFuncDervMulti(xx, df, ddf); - for (i = 0; i < mixlen; i++) - grad(i) = df[i]; -} - -void PhyloTreeMixlen::hessian(const TVector &x, THessian &hessian) { - int i, j; - double xx[mixlen]; - for (i = 0; i < mixlen; i++) - xx[i] = x(i); - int mixlen2 = mixlen*mixlen; - double df[mixlen+1], ddf[mixlen2]; - computeFuncDervMulti(xx, df, ddf); - - for (i = 0; i < mixlen; i++) - for (j = 0; j < mixlen; j++) - hessian(i, j) = ddf[i*mixlen+j]; -} -#endif - -/** - * clear the array "relative_treelen" - */ -void PhyloTreeMixlen::clear_relative_treelen() { - relative_treelen.clear(); -} - // defining log-likelihood derivative function for EM algorithm void PhyloTreeMixlen::computeFuncDerv(double value, double &df, double &ddf) { if (initializing_mixlen) { diff --git a/tree/phylotreemixlen.h b/tree/phylotreemixlen.h index 01e97dd5b..81b037e16 100644 --- a/tree/phylotreemixlen.h +++ b/tree/phylotreemixlen.h @@ -5,42 +5,30 @@ // Created by Minh Bui on 24/08/15. // // +#ifndef PHYLOTREEMIXLEN_H_ +#define PHYLOTREEMIXLEN_H_ -#ifndef __iqtree__phylotreemixlen__ -#define __iqtree__phylotreemixlen__ - -#include -#ifdef USE_CPPOPTLIB -#include "cppoptlib/meta.h" -#include "cppoptlib/boundedproblem.h" -#endif #include "iqtree.h" +/** + * Phylogenetic tree with a mixture of branch lengths + * Started within a joint project with Stephen Crotty + */ +class PhyloTreeMixlen : public IQTree { +public: + PhyloTreeMixlen(); + PhyloTreeMixlen(Alignment *aln); -/** - Phylogenetic tree with mixture of branch lengths - Started within joint project with Stephen Crotty -*/ -#ifdef USE_CPPOPTLIB -class PhyloTreeMixlen : public IQTree, public cppoptlib::BoundedProblem -#else -class PhyloTreeMixlen : public IQTree -#endif -{ - - friend class ModelFactoryMixlen; + virtual bool isMixlen() const override { return !initializing_mixlen; } -public: + virtual int getNMixlen() const override { return (initializing_mixlen) ? 1 : mixlen; } - /** - default constructor - */ - PhyloTreeMixlen(); + virtual int getCurMixture() const override { return cur_mixture; } - PhyloTreeMixlen(Alignment *aln, int mixlen); + virtual void setCurMixture(int c) override; - virtual ~PhyloTreeMixlen() override; + void clearRelativeTreelen() { relative_treelen.clear(); } /** start structure for checkpointing @@ -80,28 +68,6 @@ class PhyloTreeMixlen : public IQTree */ virtual void readTreeString(const string &tree_string) override; - virtual void initializeModel(Params ¶ms, string model_name, ModelsBlock *models_block) override; - - /** - @return true if this is a tree with mixture branch lengths, default: false - */ - virtual bool isMixlen() override { return !initializing_mixlen; } - - /** - @return number of mixture branch lengths, default: 1 - */ - virtual int getMixlen() override { - if (initializing_mixlen) - return 1; - else - return mixlen; - } - - /** - set number of mixture branch lengths - */ - void setMixlen(int mixlen); - /** @param[out] lenvec tree lengths for each class in mixlen model @param node the starting node, nullptr to start from the root @@ -109,20 +75,6 @@ class PhyloTreeMixlen : public IQTree */ virtual void treeLengths(DoubleVector &lenvec, Node *node = nullptr, Node *dad = nullptr) override; - /** - * assign branch length as mean over all branch lengths of categories - */ - void assignMeanMixBranches(Node *node = nullptr, Node *dad = nullptr); - - - /** - parse the string containing branch length(s) - by default, this will parse just one length - @param lenstr string containing branch length(s) - @param[out] branch_len output branch length(s) - */ -// virtual void parseBranchLength(string &lenstr, DoubleVector &branch_len); - /** * internal function called by printTree to print branch length * @param out output stream @@ -137,9 +89,10 @@ class PhyloTreeMixlen : public IQTree virtual void printResultTree(string suffix = "") override; /** - initialize mixture branch lengths - */ - void initializeMixBranches(PhyloNode *node = nullptr, PhyloNode *dad = nullptr); + * Set the model factory + * @param model_fac Model factory + */ + virtual void setModelFactory(ModelFactory *model_fac) override; /** initialize parameters if necessary */ void initializeMixlen(double tolerance, bool write_info); @@ -211,9 +164,6 @@ class PhyloTreeMixlen : public IQTree */ virtual double derivativeFunk(double x[], double dfx[]) override; - /** For Mixlen stuffs */ - virtual int getCurMixture() override { return cur_mixture; } - /** * Optimize current tree using NNI * @@ -222,58 +172,29 @@ class PhyloTreeMixlen : public IQTree */ virtual pair optimizeNNI(bool speedNNI = true) override; - /** number of mixture categories */ - int mixlen; - - /** current category, for optimizing branch length */ - int cur_mixture; - - -/*************** Using cppoptlib for branch length optimization ***********/ - -#ifdef USE_CPPOPTLIB - -// using typename BoundedProblem::TVector; -// using typename BoundedProblem::THessian; - - /** - * @brief returns objective value in x - * @details [long description] - * - * @param x [description] - * @return [description] - */ - double value(const TVector &x); - - /** - * @brief returns gradient in x as reference parameter - * @details should be overwritten by symbolic gradient - * - * @param grad [description] - */ - void gradient(const TVector &x, TVector &grad); - +protected: /** - * @brief This computes the hessian - * @details should be overwritten by symbolic hessian, if solver relies on hessian - */ - void hessian(const TVector &x, THessian &hessian); - -#endif + * Initialize mixlen branch lengths from branch lengths and relative_treelen + */ + void initializeMixBranches(PhyloNode *node = nullptr, PhyloNode *dad = nullptr); /** - * clear the array "relative_treelen" + * Set branch lengths as respective mean mixlen branch lengths */ - void clear_relative_treelen(); + void assignMeanMixBranches(Node *node = nullptr, Node *dad = nullptr); protected: - + /** number of mixture categories */ + int mixlen; + + /** current category, for optimizing a branch length */ + int cur_mixture; + /** relative rate, used to initialize branch lengths */ DoubleVector relative_treelen; - /** true if during initialization phase */ + /** true if within initializeMixlen() */ bool initializing_mixlen; - }; -#endif /* defined(__iqtree__phylotreemixlen__) */ +#endif /* PHYLOTREEMIXLEN_H_ */ diff --git a/utils/tools.cpp b/utils/tools.cpp index 4bb0bfbcb..a1ef1134e 100644 --- a/utils/tools.cpp +++ b/utils/tools.cpp @@ -4980,16 +4980,6 @@ void parseArg(int argc, char *argv[], Params ¶ms) { continue; } - if (strcmp(argv[cnt], "-mixlen") == 0) { - cnt++; - if (cnt >= argc) - throw "Use -mixlen "; - params.num_mixlen = convert_int(argv[cnt]); - if (params.num_mixlen < 1) - throw("-mixlen must be >= 1"); - continue; - } - if (strcmp(argv[cnt], "--link-alpha") == 0) { params.link_alpha = true; continue; @@ -7437,7 +7427,6 @@ void Params::setDefault() { lmap_num_quartets = -1; lmap_cluster_file = nullptr; print_lmap_quartet_lh = false; - num_mixlen = 1; link_alpha = false; link_model = false; model_joint = ""; diff --git a/utils/tools.h b/utils/tools.h index 58473002e..acbad34e2 100644 --- a/utils/tools.h +++ b/utils/tools.h @@ -1684,8 +1684,6 @@ class Params { /** TRUE to optimize mixture model nucleotide/amino acide frequency */ bool optimize_mixmodel_freq; - /** number of mixture branch lengths, default 1 */ - int num_mixlen; /** TRUE to always optimize rate matrix even if user parameters are specified in e.g. GTR{1,2,3,4,5} */ bool optimize_rate_matrix;