Optimize FunDi rho from initial 0.5 and new options to set initial rho and branch length

Author: bqminh

tree/phylotree.cpp

@@ -2822,7 +2822,7 @@ double PhyloTree::computeFundiLikelihood() {
         taxa_set.insert(*it);
     }
 
-    cout << "rho = " << params->alisim_fundi_proportion << endl;
+    cout << "Fundi proportion rho: " << params->alisim_fundi_proportion << endl;
 
     findNodeNames(taxa_set, central_branch, root, nullptr);
     if (!central_branch.first) {
@@ -2910,6 +2910,16 @@ double PhyloTree::computeFundiLikelihood() {
         // by BFGS algorithm
         current_it = (PhyloNeighbor*)central_branch.second;
         current_it_back = (PhyloNeighbor*)central_branch.second->node->findNeighbor(central_branch.first);
+
+        params->alisim_fundi_proportion = params->fundi_init_proportion;
+        if (params->fundi_init_branch_length > 0.0) {
+            current_it->length = params->fundi_init_branch_length;
+            current_it_back->length = params->fundi_init_branch_length;
+        }
+
+        cout << "Init rho = " << params->alisim_fundi_proportion
+            << " and init branch length = " << current_it->length << endl;
+        
         variables[1] = params->alisim_fundi_proportion;
         variables[2] = current_it->length;
         lower_bound[1] = 0.0;

utils/tools.cpp

@@ -1571,7 +1571,9 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.alisim_length_ratio = 2;
     params.birth_rate = 0.8;
     params.death_rate = 0.2;
-    params.alisim_fundi_proportion = 0;
+    params.alisim_fundi_proportion = 0.0;
+    params.fundi_init_proportion = 0.5;
+    params.fundi_init_branch_length = 0.0;
     params.alisim_distribution_definitions = NULL;
     params.alisim_skip_checking_memory = false;
     params.alisim_write_internal_sequences = false;
@@ -2921,13 +2923,34 @@ void parseArg(int argc, char *argv[], Params &params) {
                     params.alisim_fundi_proportion = 0.0;
                 } else {
                     params.alisim_fundi_proportion = convert_double(fundi_input.c_str());
-                    if (params.alisim_fundi_proportion > 1 || params.alisim_fundi_proportion < 0)
+                    if (params.alisim_fundi_proportion > 1 || params.alisim_fundi_proportion <= 0)
                         throw "Proportion in FunDi model must be positive and not greater than 1";
                 }
 
                 continue;
             }
-			if (strcmp(argv[cnt], "-ngs_gap") == 0) {
+
+            if (strcmp(argv[cnt], "--fundi-init-rho") == 0) {
+                cnt++;
+                if (cnt >= argc)
+                    throw "Use --fundi-init-rho <proportion>";
+                params.fundi_init_proportion = convert_double(argv[cnt]);
+                if (params.fundi_init_proportion >= 1 || params.fundi_init_proportion <= 0)
+                    throw "Initial proportion in FunDi model must be positive and smaller than 1";
+                continue;
+            }
+
+            if (strcmp(argv[cnt], "--fundi-init-branch") == 0) {
+                cnt++;
+                if (cnt >= argc)
+                    throw "Use --fundi-init-branch <branch_legth>";
+                params.fundi_init_branch_length = convert_double(argv[cnt]);
+                if (params.fundi_init_branch_length >= params.max_branch_length || params.fundi_init_branch_length <= 0)
+                    throw "Initial branch length in FunDi model must be positive and smaller than 10";
+                continue;
+            }
+
+            if (strcmp(argv[cnt], "-ngs_gap") == 0) {
 				params.ngs_ignore_gaps = false;
 				continue;
 			}
@@ -4777,7 +4800,7 @@ void parseArg(int argc, char *argv[], Params &params) {
                 params.modelfinder_eps = convert_double(argv[cnt]);
                 if (params.modelfinder_eps <= 0.0)
                     throw "ModelFinder epsilon must be positive";
-                if (params.modelEps > 1.0)
+                if (params.modelfinder_eps > 1.0)
                     throw "ModelFinder epsilon must not be larger than 1.0";
                 continue;
             }

utils/tools.h

@@ -2629,7 +2629,17 @@ class Params {
     *  fundi model - proportion
     */
     double alisim_fundi_proportion;
-    
+
+    /**
+    *  fundi model - initial proportion for optimisation
+    */
+    double fundi_init_proportion;
+
+    /**
+    *  fundi model - initial branch length for optimisation
+    */
+    double fundi_init_branch_length;
+
     /**
     *  distribution_definition_file
     */