Skip to content

Commit

Permalink
Fixed issues with mutation mapping. Works now in PhyTime
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Jun 19, 2018
1 parent 66b622f commit 376867f
Show file tree
Hide file tree
Showing 13 changed files with 573 additions and 326 deletions.
36 changes: 24 additions & 12 deletions src/date.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,8 @@ phydbl *DATE_MCMC(t_tree *tree)
fp_stats = tree->io->fp_out_stats;
fp_tree = tree->io->fp_out_tree;

if(tree->io->mutmap == YES) Make_MutMap(tree);

TIMES_Randomize_Tree_With_Time_Constraints(tree->rates->a_cal[0],tree);

mcmc = MCMC_Make_MCMC_Struct();
Expand All @@ -495,6 +497,7 @@ phydbl *DATE_MCMC(t_tree *tree)
MCMC_Randomize_Rate_Across_Sites(tree);
MCMC_Randomize_Rates(tree);


n_vars = 10;
adjust_len = tree->io->mcmc->chain_len_burnin;
mcmc->sample_interval = tree->io->mcmc->sample_interval;
Expand All @@ -519,7 +522,7 @@ phydbl *DATE_MCMC(t_tree *tree)
/* tree->bl_ndigits = 7; */
RATES_Update_Cur_Bl(tree);


PhyML_Printf("\n. AVX enabled: %s",
#if defined(__AVX__)
"yes"
Expand All @@ -541,6 +544,8 @@ phydbl *DATE_MCMC(t_tree *tree)
PhyML_Printf("\n. log(Pr(Seq|Tree)) = %f",tree->c_lnL);
PhyML_Printf("\n. log(Pr(Tree)) = %f",tree->rates->c_lnL_times);



tree->extra_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
tree->extra_tree->mod = tree->mod;
Copy_Tree(tree,tree->extra_tree);
Expand All @@ -558,8 +563,7 @@ phydbl *DATE_MCMC(t_tree *tree)
tree->extra_tree->mcmc = mcmc;
MCMC_Init_MCMC_Struct(NULL,NULL,mcmc);
MCMC_Complete_MCMC(mcmc,tree->extra_tree);



PhyML_Fprintf(fp_stats,"\n%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t",
"sample",
"lnL(posterior)",
Expand Down Expand Up @@ -647,9 +651,6 @@ phydbl *DATE_MCMC(t_tree *tree)
i = 0;
do { MCMC_Clock_R(tree); } while(i++ < 100);




do
{
if(tree->mcmc->run > adjust_len) for(i=0;i<tree->mcmc->n_moves;i++) tree->mcmc->adjust_tuning[i] = NO;
Expand All @@ -673,9 +674,9 @@ phydbl *DATE_MCMC(t_tree *tree)
else if(!strcmp(tree->mcmc->move_name[move],"times")) MCMC_Times_All(tree);
else if(!strcmp(tree->mcmc->move_name[move],"times_and_rates")) MCMC_Times_And_Rates_All(tree);

else if(!strcmp(tree->mcmc->move_name[move],"spr")) { MCMC_Prune_Regraft(tree); }
else if(!strcmp(tree->mcmc->move_name[move],"spr_local")) { MCMC_Prune_Regraft_Local(tree);}
else if(!strcmp(tree->mcmc->move_name[move],"spr_weighted")) { MCMC_Prune_Regraft_Weighted(tree);}
else if(!strcmp(tree->mcmc->move_name[move],"spr")) MCMC_Prune_Regraft(tree);
else if(!strcmp(tree->mcmc->move_name[move],"spr_local")) MCMC_Prune_Regraft_Local(tree);
else if(!strcmp(tree->mcmc->move_name[move],"spr_weighted")) MCMC_Prune_Regraft_Weighted(tree);

else if(!strcmp(tree->mcmc->move_name[move],"updown_t_cr")) MCMC_Updown_T_Cr(tree);
else if(!strcmp(tree->mcmc->move_name[move],"kappa")) MCMC_Kappa(tree);
Expand All @@ -688,8 +689,18 @@ phydbl *DATE_MCMC(t_tree *tree)
else if(!strcmp(tree->mcmc->move_name[move],"clade_change")) MCMC_Clade_Change(tree);
else continue;

/* PhyML_Fprintf(stderr,"\n. move: %s",tree->mcmc->move_name[move]); */

/* PhyML_Fprintf(stderr,"\n. move: %s",tree->mcmc->move_name[move]); */

/* PhyML_Printf("\n. %d move: %s",tree->mcmc->run,tree->mcmc->move_name[move]); */
/* PhyML_Printf("\n. root: %d v1: %d v2: %d b1->left: %d b1->rght: %d b2->left: %d b2->rght: %d", */
/* tree->n_root->num, */
/* tree->n_root->v[1]->num, */
/* tree->n_root->v[2]->num, */
/* tree->n_root->b[1]->left->num, */
/* tree->n_root->b[1]->rght->num, */
/* tree->n_root->b[2]->left->num, */
/* tree->n_root->b[2]->rght->num); */

if(!RATES_Check_Edge_Length_Consistency(tree))
{
PhyML_Fprintf(stderr,"\n. Issue detected by RATES_Check_Edge_Length_Consistency(tree).");
Expand Down Expand Up @@ -832,8 +843,9 @@ phydbl *DATE_MCMC(t_tree *tree)
fflush(NULL);
RATES_Update_Cur_Bl(tree);

if(tree->io->mutmap == YES) Sample_Ancestral_Seq(YES,NO,tree);

tree->mcmc->sample_num++;

}

tree->mcmc->run++;
Expand Down
14 changes: 5 additions & 9 deletions src/io.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3955,21 +3955,17 @@ void Read_Clade_Priors(char *file_name, t_tree *tree)
PhyML_Printf("\n. WARNING: could not find any clade in the tree referred to with the following taxon names:");
for(i=0;i<clade_size;i++) PhyML_Printf("\n. \"%s\"",clade_list[i]);
PhyML_Printf("\n. .................................................................");
/* sleep(3); */
}
else
{
tree->rates->t_has_prior[node_num] = YES;
tree->rates->t_prior_min[node_num] = MIN(prior_low,prior_up);
tree->rates->t_prior_max[node_num] = MAX(prior_low,prior_up);

tree->rates->t_prior_min[node_num] = -MAX(prior_low,prior_up);
tree->rates->t_prior_max[node_num] = -MIN(prior_low,prior_up);

/* Sergei to add stuff here re calibration object */


if(FABS(prior_low - prior_up) < 1.E-6 && tree->a_nodes[node_num]->tax == YES)
tree->rates->nd_t[node_num] = prior_low;

if(fabs(prior_low - prior_up) < 1.E-6 && tree->a_nodes[node_num]->tax == YES)
tree->rates->nd_t[node_num] = prior_low;

PhyML_Printf("\n");
PhyML_Printf("\n. [%3d]..................................................................",n_clade_priors);
PhyML_Printf("\n. Node %4d matches the clade referred to with the following taxon names:",node_num);
Expand Down
Loading

0 comments on commit 376867f

Please sign in to comment.