Mixlen refactoring: Remove dependency on ModelFactoryMixlen and cppoptlib, add support for the -ft option

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 &params, 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 &params, 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 &params,
     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(&params);
     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 &params, 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) {

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

model/modelfactory.cpp

@@ -41,6 +41,7 @@
 #include "rateheterotachyinvar.h"
 //#include "ngs.h"
 #include <string>
+#include "tree/phylotreemixlen.h"
 #include "utils/timeutil.h"
 #include "nclextra/myreader.h"
 #include <sstream>
@@ -1041,7 +1042,6 @@ ModelFactory::ModelFactory(Params &params, 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();
 }
 

model/modelfactory.h

@@ -266,7 +266,11 @@ class ModelFactory : public unordered_map<int, double*>, 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