Fixed bug in Map_Mutations

stephaneguindon · Jun 19, 2018 · 7c63f9c · 7c63f9c
1 parent 376867f
commit 7c63f9c
Show file tree

Hide file tree

Showing 5 changed files with 24 additions and 15 deletions.
diff --git a/src/cl.c b/src/cl.c
@@ -879,86 +879,103 @@ int Read_Command_Line(option *io, int argc, char **argv)
 	    else if(strcmp(optarg, "DAYHOFF") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = DAYHOFF;
 	      }
 	    else if(strcmp (optarg, "JTT") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = JTT;
 	      }
 	    else if(strcmp(optarg, "MTREV") == 0)
 	      {
 		io->datatype             = AA;
+		io->mod->ns              = 20;
 		io->mod->whichmodel      = MTREV;
 	      }
 	    else if(strcmp (optarg, "LG") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = LG;
 	      }
 	    else if(strcmp (optarg, "WAG") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = WAG;
 	      }
 	    else if(strcmp(optarg, "DCMUT") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = DCMUT;
 	      }
 	    else if(strcmp (optarg, "RTREV") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = RTREV;
 	      }
 	    else if(strcmp(optarg, "CPREV") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = CPREV;
 	      }
 	    else if(strcmp(optarg, "VT") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = VT;
 	      }
 	    else if(strcmp(optarg, "BLOSUM62") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = BLOSUM62;
 	      }
 	    else if(strcmp(optarg, "MTMAM") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = MTMAM;
 	      }
 	    else if (strcmp(optarg,"MTART") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = MTART;
 	      }
 	    else if (strcmp(optarg,"HIVW") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = HIVW;
 	      }
 	    else if(strcmp(optarg, "HIVB") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = HIVB;
 	      }
 	    else if(strcmp(optarg, "AB") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = AB;
 	      }
 	    else if (strcmp(optarg, "CUSTOM") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = CUSTOMAA;
 	      }
 	    else if(strcmp(optarg, "FLU") == 0)
 	      {
 		io->datatype              = AA;
+		io->mod->ns               = 20;
 		io->mod->whichmodel       = FLU;
 	      }
 	    else
@@ -1577,7 +1594,7 @@ int Read_Command_Line(option *io, int argc, char **argv)
       io->mod->s_opt->opt_rr = 1;
     }
 
-
+  
   return 1;
 }
 

diff --git a/src/lk.c b/src/lk.c
@@ -3195,35 +3195,29 @@ void Map_Mutations(t_node *a, t_node *d, int sa, int sd, t_edge *b, int site, in
   int *mut; // Array of mutations
   int thismut;
   int n_mut_branch;
-  phydbl cr,rr;
   int n_iter;
   int first_mut;
   phydbl T;
   int ns,tax_idx;
-
+  phydbl rr;
+
   ns = tree->mod->ns;
 
   all_probs = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
   mut = (int *)mCalloc(ns*ns,sizeof(int));
 
   // Site relative rate
   rr = tree->mod->ras->gamma_rr->v[rcat];
-
-#ifdef PHYTIME
-  cr = tree->rates->clock_r;
-#elif defined(PHYML)
-  cr = 0.0;
-#endif
-
+
   // Rate matrix
   Q = tree->mod->r_mat->qmat->v;
 
   // Length of the 'time' interval considered.
   T = 0.0;
 #ifdef PHYTIME
-  T = tree->rates->cur_l[d->num];
+  T = tree->rates->cur_l[d->num]*rr;
 #elif defined(PHYML)
-  T = b->l->v;
+  T = b->l->v*rr;
 #endif
 
 

diff --git a/src/main.c b/src/main.c
@@ -115,6 +115,7 @@ int main(int argc, char **argv)
       mod = io->mod;
       orig_random_input_tree = io->mod->s_opt->random_input_tree;
 
+
       if(io->data)
         {
           if(io->n_data_sets > 1) PhyML_Printf("\n. Data set [#%d]\n",num_data_set+1);

diff --git a/src/make.c b/src/make.c
@@ -1001,7 +1001,6 @@ t_rmat *Make_Rmat(int ns)
   r_mat->qmat->v                = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
   r_mat->qmat_buff->v           = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
 
-
   return(r_mat);
 }
 

diff --git a/src/times.c b/src/times.c
@@ -1954,8 +1954,6 @@ void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tr
         }
       while(tree);
     }
-
-  t_tree *tree = mixt_tree->next;
 
   Free(tips);
   Free(nd_list);
-Original file line number
+Diff line change
@@ Expand Up / @@ -1001,7 +1001,6 @@ t_rmat *Make_Rmat(int ns) @@
       r_mat->qmat->v                = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
       r_mat->qmat_buff->v           = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
       return(r_mat);
     }
@@ Expand Down @@