Skip to content

Commit

Permalink
Fixed issues with multiple calibration constraints in PhyTime
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Jun 20, 2018
1 parent 2b76ab5 commit eddd319
Show file tree
Hide file tree
Showing 5 changed files with 127 additions and 51 deletions.
26 changes: 7 additions & 19 deletions src/date.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -444,14 +444,14 @@ int DATE_Check_Calibration_Constraints(t_tree *tree)

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

// Check that time constraints are satisfied. Note: it is a
// requirement to verify that the age of the root node is also
// within correct boundaries
int DATE_Check_Time_Constraints(t_tree *tree)
{
int i;

For(i,2*tree->n_otu-1)
for(int i=0;i<2*tree->n_otu-1;++i)
{
if(tree->a_nodes[i] != tree->n_root && tree->a_nodes[i]->tax == NO)
if(tree->a_nodes[i]->tax == NO)
{
if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] ||
tree->rates->nd_t[i] < tree->rates->t_prior_min[i])
Expand Down Expand Up @@ -605,7 +605,6 @@ phydbl *DATE_MCMC(t_tree *tree)
else
for(i=0;i<2*tree->n_otu-1;++i) PhyML_Fprintf(fp_stats,"nr%d\t",i);


for(i=0;i<2*tree->n_otu-1;++i) if(tree->a_nodes[i]->tax == NO) PhyML_Fprintf(fp_stats,"t%d\t",i);

PhyML_Fprintf(fp_stats,"accRT\t");
Expand Down Expand Up @@ -690,16 +689,6 @@ phydbl *DATE_MCMC(t_tree *tree)
else continue;

/* 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))
{
Expand Down Expand Up @@ -783,11 +772,10 @@ phydbl *DATE_MCMC(t_tree *tree)
for(i=0;i<tree->rates->n_cal;++i)
{
t_cal *cal = tree->rates->a_cal[i];
/* PhyML_Fprintf(fp_stats,"%d\t",cal->current_clade_idx); */
PhyML_Fprintf(fp_stats,"%s\t",cal->clade_list[cal->current_clade_idx]->id);
PhyML_Fprintf(fp_stats,"%d\t",cal->current_clade_idx);
/* PhyML_Fprintf(fp_stats,"%s\t",cal->clade_list[cal->current_clade_idx]->id); */
}


if(tree->rates->model == THORNE ||
tree->rates->model == LOGNORMAL ||
tree->rates->model == STRICTCLOCK)
Expand Down
19 changes: 16 additions & 3 deletions src/mcmc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4674,7 +4674,7 @@ void MCMC_Death_Rate(t_tree *tree)
tree->extra_tree->eval_glnL = YES;
tree->extra_tree->rates->c_lnL_rates = UNLIKELY;
tree->extra_tree->c_lnL = UNLIKELY;

TIMES_Randomize_Tree_With_Time_Constraints(tree->extra_tree->rates->a_cal[0],tree->extra_tree);
tree->extra_tree->rates->birth_rate = tree->rates->birth_rate;
tree->extra_tree->rates->death_rate = new_death_rate;
Expand All @@ -4684,8 +4684,21 @@ void MCMC_Death_Rate(t_tree *tree)
{
PhyML_Fprintf(stderr,"\n. glnL=%f",tree->extra_tree->rates->c_lnL_times);
PhyML_Fprintf(stderr,"\n. death=%G birth=%G [%G]",new_death_rate,tree->rates->birth_rate,tree->extra_tree->rates->birth_rate);
PhyML_Fprintf(stderr,"\n");
for(int i=0;i<tree->extra_tree->rates->n_cal;++i)
{
t_cal *cal = tree->extra_tree->rates->a_cal[i];
PhyML_Fprintf(stderr,"\n. calibration %s applies to clade %s\t",cal->id,cal->clade_list[cal->current_clade_idx]->id);
for(int j=0;j<cal->clade_list_size;++j)
{
t_clad *clade = cal->clade_list[j];
PhyML_Fprintf(stderr,"time:%G\t",tree->extra_tree->rates->nd_t[clade->target_nd->num]);
}
}

PhyML_Fprintf(stderr,"\n");
TIMES_Lk_Times(YES,tree->extra_tree);
Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
assert(FALSE);
}

i = 0;
Expand Down Expand Up @@ -4855,7 +4868,7 @@ void MCMC_Birth_Death_Updown(t_tree *tree)
PhyML_Fprintf(stderr,"\n. glnL=%f",tree->extra_tree->rates->c_lnL_times);
PhyML_Fprintf(stderr,"\n. death=%G birth=%G [%G]",new_death_rate,tree->rates->birth_rate,tree->extra_tree->rates->birth_rate);
TIMES_Lk_Times(YES,tree->extra_tree);
Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
assert(FALSE);
}


Expand Down
93 changes: 67 additions & 26 deletions src/times.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -697,7 +697,6 @@ phydbl TIMES_Lk_Times(int verbose, t_tree *tree)
DATE_Assign_Primary_Calibration(tree);
DATE_Update_T_Prior_MinMax(tree);
tree->rates->c_lnL_times = TIMES_Lk_Birth_Death(verbose,tree);
/* if(tree->rates->c_lnL_times > UNLIKELY) tree->rates->c_lnL_times = -1.; */
return(tree->rates->c_lnL_times);
}

Expand Down Expand Up @@ -1381,17 +1380,35 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)

if(b < d)
{
if(verbose) printf("\n. b: %G d: %G",b,d);
tree->rates->c_lnL_times = UNLIKELY;
if(verbose) PhyML_Printf("\n. b < d");
return UNLIKELY;
}

// Verify that calibration constraints are satisfied
for(i=0;i<2*tree->n_otu-1;++i)
if(tree->a_nodes[i]->tax == NO)
{
if(tree->rates->nd_t[i] < tree->rates->t_prior_min[i] ||
tree->rates->nd_t[i] > tree->rates->t_prior_max[i])
{
if(verbose) printf("\n. node %d @ time %f min: %f max: %f",i,tree->rates->nd_t[i],tree->rates->t_prior_min[i],tree->rates->t_prior_max[i]);
tree->rates->c_lnL_times = UNLIKELY;
if(verbose)
{
PhyML_Printf("\n. Time outside calibration range : %G [%G,%G]",tree->rates->nd_t[i],tree->rates->t_prior_min[i],tree->rates->t_prior_max[i]);
PhyML_Printf("\n. Clade incriminated: ");
if(tree->a_nodes[i] == tree->n_root)
{
List_Taxa_In_Clade(tree->n_root,tree->n_root->v[1],tree);
List_Taxa_In_Clade(tree->n_root,tree->n_root->v[2],tree);
}
else
{
assert(tree->a_nodes[i]->anc);
List_Taxa_In_Clade(tree->a_nodes[i]->anc,tree->a_nodes[i],tree);
}
PhyML_Printf("\n");
}
return UNLIKELY;
}
}
Expand All @@ -1417,9 +1434,18 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)
// Equation 6 in Yang and Rannala, 1997 with rho=1
logp_1t = 2.*logbmd - 2.*log(b-d*exp((d-b)*t)) + (d-b)*t;
lnL += logb + logp_1t - lognut1;


if(!(lnL > UNLIKELY))
{
PhyML_Printf("\n. logb: %G pt: %G nut1: %G lognut1: %G t: %G logp_1t: %G\n",logb,pt,nut1,lognut1,t,logp_1t);
tree->rates->c_lnL_times = UNLIKELY;
return UNLIKELY;
}
}

lnL += LnGamma(n-1);


/* t = FABS(tree->rates->nd_t[tree->n_root->num]); */
/* lnL += -bmd*t - log(b-d*pow(expmbmd,t)); */
Expand Down Expand Up @@ -1447,6 +1473,13 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)
// Equation 6 in Yang and Rannala, 1997 with rho=1
logp_1t = - b*t;
lnL += logb + logp_1t - lognut1;

if(!(lnL > UNLIKELY))
{
PhyML_Printf("\n. lognut1: %G t: %G logp_1t: %G",lognut1,t,logp_1t);
tree->rates->c_lnL_times = UNLIKELY;
return UNLIKELY;
}
}

lnL += LnGamma(n-1);
Expand All @@ -1466,7 +1499,8 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)
}
else if(b < bmin && d > dmin)
{
if(verbose) printf("\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
PhyML_Printf("\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
tree->rates->c_lnL_times = UNLIKELY;
return UNLIKELY;
}
else if(Are_Equal(bmd,0.0,bmin/10.) == YES) // Critical birth-death process
Expand All @@ -1484,6 +1518,13 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)

// Equation 7 in Yang and Rannala, 1997 with rho=1
lnL += log((1.+d)/pow(1.+d*t,2));

if(!(lnL > UNLIKELY))
{
PhyML_Printf("\n. logb: %G t: %G",logb,t);
tree->rates->c_lnL_times = UNLIKELY;
return UNLIKELY;
}
}

lnL += LnGamma(n-1);
Expand All @@ -1502,29 +1543,30 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)
}
else if(b < bmin && d < dmin) // Birth and death rates are below their limits
{
if(verbose) printf("\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
PhyML_Fprintf(stderr,"\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
tree->rates->c_lnL_times = UNLIKELY;
return -INFINITY;
}
else if(b > bmax && d > dmax)
{
if(verbose) printf("\n. b: %G bmax: %G d: %G dmax: %G",b,bmax,d,dmax);
PhyML_Fprintf(stderr,"\n. b: %G bmax: %G d: %G dmax: %G",b,bmax,d,dmax);
tree->rates->c_lnL_times = UNLIKELY;
return -INFINITY;
}
else
{
assert(FALSE);
}

if(isnan(lnL) || isinf(FABS(lnL)))
if(isnan(lnL) || isinf(fabs(lnL)) || !(lnL > UNLIKELY))
{
if(verbose) printf("\n. lnL: %f",lnL);
PhyML_Fprintf(stderr,"\n. lnL times: %f",lnL);
tree->rates->c_lnL_times = UNLIKELY;
return UNLIKELY;
}

tree->rates->c_lnL_times = lnL;


return(lnL);
}

Expand Down Expand Up @@ -1644,7 +1686,8 @@ void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tr
int i,j,nd_num,*cal_ordering,n_cal,swap,list_size,tmp,orig_is_mixt_tree,repeat,n_max_repeats,tip_num,*no_cal_tip_num,n_no_cal_tip_num,*permut;
t_cal *cal;
t_clad *clade;

int idx;

assert(mixt_tree->rates);

tips = (t_node **)mCalloc(mixt_tree->n_otu,sizeof(t_node *));
Expand All @@ -1656,6 +1699,7 @@ void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tr
n_max_repeats = 1000;
cal = NULL;
clade = NULL;
idx = -1;

// List node indices that are not in any calibration set
no_cal_tip_num = NULL;
Expand All @@ -1679,7 +1723,7 @@ void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tr
n_no_cal_tip_num++;
}
}


for(repeat=0;repeat<n_max_repeats;++repeat)
{
Expand Down Expand Up @@ -1880,21 +1924,18 @@ void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tr
DATE_Assign_Primary_Calibration(mixt_tree);
DATE_Update_T_Prior_MinMax(mixt_tree);

/* { */
/* Print_Node(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree); */
/* Print_Node(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree); */
/* fflush(NULL); */

/* int i; */
/* for(i=0;i<mixt_tree->rates->n_cal;i++) */
/* { */
/* PhyML_Printf("\n. Node number to which calibration [%d] applies to is [%d]",i,Find_Clade(mixt_tree->rates->a_cal[i]->target_tax, */
/* mixt_tree->rates->a_cal[i]->n_target_tax, */
/* mixt_tree)); */
/* PhyML_Printf("\n. Lower bound set to: %15f time units.",mixt_tree->rates->a_cal[i]->lower); */
/* PhyML_Printf("\n. Upper bound set to: %15f time units.",mixt_tree->rates->a_cal[i]->upper); */
/* } */
/* } */

/* for(int i=0;i<mixt_tree->rates->n_cal;i++) */
/* { */
/* cal = mixt_tree->rates->a_cal[i]; */
/* clade = cal->clade_list[cal->current_clade_idx]; */
/* idx = Find_Clade(clade->tax_list,clade->n_tax,mixt_tree); */
/* PhyML_Printf("\n. Node number to which calibration [%d] applies to is [%d]",i,idx); */
/* PhyML_Printf("\n. Lower bound set to: %15f time units.",mixt_tree->rates->a_cal[i]->lower); */
/* PhyML_Printf("\n. Upper bound set to: %15f time units.",mixt_tree->rates->a_cal[i]->upper); */
/* PhyML_Printf("\n. t_prior_min: %G t_prior_max: %G",mixt_tree->rates->t_prior_min[idx],mixt_tree->rates->t_prior_max[idx]); */
/* PhyML_Printf("\n. Time set to %G",mixt_tree->rates->nd_t[idx]); */
/* } */

if(!DATE_Check_Calibration_Constraints(mixt_tree))
{
Expand Down
37 changes: 35 additions & 2 deletions src/utilities.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12068,10 +12068,12 @@ char *Mutation_Id(int mut_idx, t_tree *tree)

s = (char *)mCalloc(20,sizeof(char));

strcpy(s,"from ");
/* strcpy(s,"from "); */
strcpy(s," ");
c = Reciproc_Assign_State((int)(mut_idx/ns),tree->mod->io->datatype);
sprintf(s+strlen(s),"%c",c);
strcat(s," to ");
/* strcat(s," to "); */
strcat(s," ");
c = Reciproc_Assign_State((int)(mut_idx%ns),tree->mod->io->datatype);
sprintf(s+strlen(s),"%c",c);

Expand Down Expand Up @@ -12103,3 +12105,34 @@ void Random_Tax_Idx(t_node *a, t_node *d, int *idx, t_tree *tree)

/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/

void List_Taxa_In_Clade(t_node *a, t_node *d, t_tree *tree)
{
if(d->tax == YES)
{
PhyML_Printf("\n- [%50s]",d->name);
}
else
{
for(int i=0;i<3;++i)
{
if(d->v[i] != a && d->b[i] != tree->e_root)
{
List_Taxa_In_Clade(d,d->v[i],tree);
}
}
}
}

/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/
3 changes: 2 additions & 1 deletion src/utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,7 @@ extern int TIME;
#define GOLDEN_C (1.0-GOLDEN_R)
#define N_MAX_INSERT 20
#define N_MAX_OTU 4000
#define UNLIKELY -1.e10
#define UNLIKELY -1.e20
#define NJ_SEUIL 0.1
#define ROUND_MAX 100
#define DIST_MAX 2.00
Expand Down Expand Up @@ -2293,6 +2293,7 @@ void Remove_Duplicates(calign *data, t_mod *mod, option *io);
short int Are_Sequences_Identical(align *seq1, align *seq2);
char *Mutation_Id(int mut_idx, t_tree *tree);
void Random_Tax_Idx(t_node *a, t_node *d, int *idx, t_tree *tree);
void List_Taxa_In_Clade(t_node *a, t_node *d, t_tree *tree);

#include "xml.h"
#include "free.h"
Expand Down

0 comments on commit eddd319

Please sign in to comment.