diff --git a/src/date.c b/src/date.c index 3ddc441e..c826cae0 100755 --- a/src/date.c +++ b/src/date.c @@ -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]) @@ -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"); @@ -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)) { @@ -783,11 +772,10 @@ phydbl *DATE_MCMC(t_tree *tree) for(i=0;irates->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) diff --git a/src/mcmc.c b/src/mcmc.c index 7ac2a0c3..64df5bad 100755 --- a/src/mcmc.c +++ b/src/mcmc.c @@ -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; @@ -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;iextra_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;jclade_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; @@ -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); } diff --git a/src/times.c b/src/times.c index c1d43fd1..9fac5d0b 100755 --- a/src/times.c +++ b/src/times.c @@ -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); } @@ -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; } } @@ -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)); */ @@ -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); @@ -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 @@ -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); @@ -1502,12 +1543,14 @@ 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 @@ -1515,16 +1558,15 @@ phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree) 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); } @@ -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 *)); @@ -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; @@ -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;repeatn_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;irates->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;irates->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)) { diff --git a/src/utilities.c b/src/utilities.c index f822266f..672fddb5 100755 --- a/src/utilities.c +++ b/src/utilities.c @@ -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); @@ -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); + } + } + } +} + +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ diff --git a/src/utilities.h b/src/utilities.h index 02072b2a..198ba045 100755 --- a/src/utilities.h +++ b/src/utilities.h @@ -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 @@ -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"