diff --git a/configure.ac b/configure.ac index a22de8fa..1bd52ca1 100755 --- a/configure.ac +++ b/configure.ac @@ -1,7 +1,7 @@ # -*- Autoconf -*- # Process this file with autoconf to produce a configure script. -AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20211231" | tr -d '\n'"]),[guindon@lirmm.fr]) +AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20220408" | tr -d '\n'"]),[guindon@lirmm.fr]) AM_INIT_AUTOMAKE([foreign]) AC_CONFIG_SRCDIR([src/main.c],[doc/phyml-manual.tex]) AC_CONFIG_HEADERS([config.h]) diff --git a/src/init.c b/src/init.c index e387b1a5..823eacd3 100644 --- a/src/init.c +++ b/src/init.c @@ -862,7 +862,7 @@ void RATES_Init_Rate_Struct(t_rate *rates, t_rate *existing_rates, int n_otu) rates->clock_r_prior_var = 1.0E-8; rates->max_rate = 1.E+1; - rates->min_rate = 1.E-5; + rates->min_rate = 1.E-1; /* rates->max_rate = 5.; */ /* rates->min_rate = 0.2; */ @@ -3499,11 +3499,8 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t) { for(int i=0;i<2*n_otu-1;++i) t->sigsq_scale[i] = 1.0; - /* !!!!!!!!!!!!!!!!!!!! */ - t->sigsq_scale_min = 1.E-4; - t->sigsq_scale_max = 1.E+4; - /* t->sigsq_scale_min = 0.1; */ - /* t->sigsq_scale_max = 10.; */ + t->sigsq_scale_min = 0.1; + t->sigsq_scale_max = 10.; t->rrw_norm_fact = 1.0; diff --git a/src/io.c b/src/io.c index 8194e95a..800abf7a 100644 --- a/src/io.c +++ b/src/io.c @@ -6657,6 +6657,12 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree) PhyML_Fprintf(fp_stats,"%s\t","lnAlgn"); PhyML_Fprintf(fp_stats,"%s\t","lnSpac"); PhyML_Fprintf(fp_stats,"%s\t","lnRate"); + /* PhyML_Fprintf(fp_stats,"%s\t","lnRateAutoCor"); */ + /* PhyML_Fprintf(fp_stats,"%s\t","lnRateClock"); */ + /* PhyML_Fprintf(fp_stats,"%s\t","lnRateRoot1"); */ + /* PhyML_Fprintf(fp_stats,"%s\t","lnRateRoot2"); */ + /* PhyML_Fprintf(fp_stats,"%s\t","lnRater1"); */ + /* PhyML_Fprintf(fp_stats,"%s\t","lnRater2"); */ PhyML_Fprintf(fp_stats,"%s\t","lnTime"); PhyML_Fprintf(fp_stats,"%s\t","substRate"); for(int i=0;immod->n_dim;++i) PhyML_Fprintf(fp_stats,"%s%s\t","sigSq",(i==0)?("Lon"):((i==1)?("Lat"):("xx"))); @@ -6675,6 +6681,8 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree) /* for(int i=0;imod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"rr(%d)\t",i+1); */ /* } */ PhyML_Fprintf(fp_stats,"%s\t","nu"); + PhyML_Fprintf(fp_stats,"%s\t","rrNormFact"); + PhyML_Fprintf(fp_stats,"%s\t","rrwNormFact"); /* PhyML_Fprintf(fp_stats,"%s\t","treeLen"); */ /* PhyML_Fprintf(fp_stats,"%s\t","speed"); */ @@ -6738,6 +6746,12 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree) PhyML_Fprintf(fp_stats,"%.2f\t",tree->c_lnL); PhyML_Fprintf(fp_stats,"%.2f\t",tree->mmod->c_lnL); PhyML_Fprintf(fp_stats,"%.2f\t",tree->rates->c_lnL); + /* PhyML_Fprintf(fp_stats,"%f\t",RATES_Autocor_Prior(tree)); */ + /* PhyML_Fprintf(fp_stats,"%f\t",RATES_Clock_R_Prior(tree)); */ + /* PhyML_Fprintf(fp_stats,"%f\t",RATES_Lk_Core(-1.,tree->rates->br_r[tree->n_root->v[1]->num],-1,-1,-1,-1,-1,tree->times->nd_t[tree->n_root->v[1]->num] - tree->times->nd_t[tree->n_root->num],tree)); */ + /* PhyML_Fprintf(fp_stats,"%f\t",RATES_Lk_Core(-1.,tree->rates->br_r[tree->n_root->v[2]->num],-1,-1,-1,-1,-1,tree->times->nd_t[tree->n_root->v[2]->num] - tree->times->nd_t[tree->n_root->num],tree)); */ + /* PhyML_Fprintf(fp_stats,"%f\t",tree->times->nd_t[tree->n_root->v[1]->num]); */ + /* PhyML_Fprintf(fp_stats,"%f\t",tree->times->nd_t[tree->n_root->v[2]->num]); */ PhyML_Fprintf(fp_stats,"%.2f\t",tree->times->c_lnL); PhyML_Fprintf(fp_stats,"%g\t",tree->rates->clock_r); for(int i=0;immod->n_dim;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->sigsq[i]); @@ -6760,6 +6774,8 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree) /* for(int i=0;imod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mod->ras->gamma_rr->v[i]); */ /* } */ PhyML_Fprintf(fp_stats,"%g\t",tree->rates->nu); + PhyML_Fprintf(fp_stats,"%g\t",tree->rates->norm_fact); + PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->rrw_norm_fact); /* PhyML_Fprintf(fp_stats,"%g\t",Tree_Length(tree)); */ /* PhyML_Fprintf(fp_stats,"%g\t",difftime(tree->mcmc->time_end,tree->mcmc->time_beg)/(phydbl)tree->mcmc->sample_interval); */ diff --git a/src/mcmc.c b/src/mcmc.c index 6c3e2f64..31098d9f 100644 --- a/src/mcmc.c +++ b/src/mcmc.c @@ -1313,7 +1313,7 @@ void MCMC_Root_Time(t_tree *tree, int print) phydbl r1_new,r2_new; phydbl t2,t1; t_node *v2,*v1; - phydbl tune,K; + phydbl tune,K,eps; int move_num,err,iter,n_tot_iter; t_node *root; short int failed; @@ -1323,12 +1323,13 @@ void MCMC_Root_Time(t_tree *tree, int print) if(tree->eval_alnL == YES) Lk(NULL,tree); root = tree->n_root; - + eps = -1.; failed = NO; move_num = tree->mcmc->num_move_root_time; - tune = MIN(10.,tree->mcmc->tune_move[tree->mcmc->num_move_root_time]); + /* tune = MIN(10.,tree->mcmc->tune_move[tree->mcmc->num_move_root_time]); */ + tune = 0.2; r1_cur = -1.; r2_cur = -1.; @@ -1386,26 +1387,27 @@ void MCMC_Root_Time(t_tree *tree, int print) t0_cur = tree->times->nd_t[root->num]; /* t_min = -TWO_TO_THE_LARGE; */ - t_min = t0_cur - 1.*(t_max - t0_cur); + t_min = t0_cur - (t_max - t0_cur); + assert(tree->young_disk); + eps = 0.05 * (tree->young_disk->time - t_max); K = (tree->young_disk->time - t0_cur) * tune; - t0_new = Rnorm_Trunc(t0_cur,K,t_min,t_max,&err); - /* PhyML_Printf("\n. K=%f t0_cur: %f t0_new: %f t_max:%f t_min:%f",K,t0_cur,t0_new,t_max,t_min); */ - hr -= Log_Dnorm_Trunc(t0_new,t0_cur,K,t_min,t_max,&err); - t_min = t0_new - 1.*(t_max - t0_new); + t0_new = Rnorm_Trunc(t0_cur,K,t_min-eps,t_max,&err); + hr -= Log_Dnorm_Trunc(t0_new,t0_cur,K,t_min-eps,t_max,&err); K = (tree->young_disk->time - t0_new) * tune; - hr += Log_Dnorm_Trunc(t0_cur,t0_new,K,t_min,t_max,&err); + t_min = t0_new - (t_max - t0_new); + hr += Log_Dnorm_Trunc(t0_cur,t0_new,K,t_min-eps,t_max,&err); - t_min = t0_cur - 1.*(t_max - t0_cur); - + t_min = t0_cur - (t_max - t0_cur); + if(isnan(t0_new)) { PhyML_Printf("\n. t_max=%f t0_cur=%f t0_new=%f",t_max,t0_cur,t0_new); assert(FALSE); } - if(t0_new > t_min && t0_new < t_max) + if(t0_new > t_min-eps && t0_new < t_max) { tree->times->nd_t[root->num] = t0_new; @@ -1426,7 +1428,7 @@ void MCMC_Root_Time(t_tree *tree, int print) if(failed == NO) { - if(tree->eval_alnL == YES) hr += log(r1_new/r1_cur) + log(r2_new/r2_cur); + if(tree->eval_rlnL == YES) hr += log(r1_new/r1_cur) + log(r2_new/r2_cur); tree->rates->br_r[tree->n_root->v[1]->num] = r1_new; tree->rates->br_r[tree->n_root->v[2]->num] = r2_new; @@ -1465,22 +1467,28 @@ void MCMC_Root_Time(t_tree *tree, int print) /* if(tree->aux_tree != NULL) */ - /* { */ - /* PhyML_Printf("\n. [%4d] Root time t: %10f->%10f t1:%f t2: %f r1: %10f->%10f r2: %10f->%10f R: %12f alnL:%10f->%10f[%d] tlnL: %10f->%10f[%d] glnL: %10f->%10f[%d] rlnL: %10f->%10f[%d] tune: %10f hr: %f ratio: %10f K: %10f failed: %d", */ - /* tree->mcmc->run, */ - /* t0_cur,t0_new, */ - /* t1,t2, */ - /* r1_cur,r1_new, */ - /* r2_cur,r2_new, */ - /* tree->rates->norm_fact, */ - /* cur_lnL_seq,new_lnL_seq,tree->eval_alnL, */ - /* cur_lnL_time,new_lnL_time,tree->eval_glnL, */ - /* cur_lnL_loc,new_lnL_loc,tree->eval_glnL, */ - /* cur_lnL_rate,new_lnL_rate,tree->eval_rlnL, */ - /* tune, */ - /* hr, */ - /* ratio,K,failed); */ - /* } */ + { + PhyML_Printf("\n."); + PhyML_Printf("\n. %c [%4d] \n. t: %1f->%1f t1:%f t2: %f \n. r1: %10f->%10f r2: %10f->%10f \n. R: %12f alnL:%10f->%10f[%d] \n. tlnL: %10f->%10f[%d] \n. glnL: %10f->%10f[%d] \n. rlnL: %10f->%10f[%d] \n. tune: %10f hr: %f [%f,%f] [%f,%f,%f] ratio: %10f K: %10f failed: %d", + (tree->aux_tree != NULL) ? '*' : 'X', + tree->mcmc->run, + t0_cur,t0_new, + t1,t2, + r1_cur,r1_new, + r2_cur,r2_new, + tree->rates->norm_fact, + cur_lnL_seq,new_lnL_seq,tree->eval_alnL, + cur_lnL_time,new_lnL_time,tree->eval_glnL, + cur_lnL_loc,new_lnL_loc,tree->eval_glnL, + cur_lnL_rate,new_lnL_rate,tree->eval_rlnL, + tune, + hr, + Log_Dnorm_Trunc(t0_new,t0_cur,(tree->young_disk->time - t0_cur) * tune,t0_cur - (t_max - t0_cur)-eps,t_max,&err), + Log_Dnorm_Trunc(t0_cur,t0_new,(tree->young_disk->time - t0_new) * tune,t0_new - (t_max - t0_new)-eps,t_max,&err), + t0_cur - (t_max - t0_cur)-eps,t0_new - (t_max - t0_new)-eps,t_max, + ratio,K,failed); + PhyML_Printf("\n."); + } assert(isnan(u) == NO && isinf(fabs(u)) == NO); @@ -6989,8 +6997,8 @@ void MCMC_PHYREX_Sigsq_Scale(t_tree *tree, int print) sigsq_scale_bkp = (phydbl *)mCalloc(2*tree->n_otu-2,sizeof(phydbl)); /* K = 10.0/tree->n_otu; // Small so as to facilitate convergence of exchange algo */ - K = 1.0; - n_changes = (int)(1. + 0.2*tree->n_otu); // Same here + K = 0.3; + n_changes = (int)(1. + 0.1*tree->n_otu); // Same here hr = 0.0; @@ -8463,7 +8471,7 @@ void MCMC_PHYREX_Node_Times_Pre(t_ldsk *a_ldsk, t_ldsk *d_ldsk, t_tree *tree, in phydbl cur_tlnL,new_tlnL; phydbl ratio,alpha,hr; int i; - int move_num; + int move_num,err; t_ldsk *ldsk_next; phydbl r0_cur,r1_cur,r2_cur; phydbl r0_new,r1_new,r2_new; @@ -8524,10 +8532,10 @@ void MCMC_PHYREX_Node_Times_Pre(t_ldsk *a_ldsk, t_ldsk *d_ldsk, t_tree *tree, in assert(t_max > t_min); - new_t = Uni()*(t_max - t_min) + t_min; - /* new_t = Rnorm_Trunc(cur_t,(t_max - t_min)/50.,t_min,t_max,&err); */ - /* hr -= Log_Dnorm_Trunc(new_t,cur_t,(t_max - t_min)/50.,t_min,t_max,&err); */ - /* hr += Log_Dnorm_Trunc(cur_t,new_t,(t_max - t_min)/50.,t_min,t_max,&err); */ + /* new_t = Uni()*(t_max - t_min) + t_min; */ + new_t = Rnorm_Trunc(cur_t,(t_max - t_min)/20.,t_min,t_max,&err); + hr -= Log_Dnorm_Trunc(new_t,cur_t,(t_max - t_min)/20.,t_min,t_max,&err); + hr += Log_Dnorm_Trunc(cur_t,new_t,(t_max - t_min)/20.,t_min,t_max,&err); /* new_t = Rnorm_Trunc(cur_t,(t_max - t_min)/100.,t_min,t_max,&err); */ /* hr -= Log_Dnorm_Trunc(new_t,cur_t,(t_max - t_min)/100.,t_min,t_max,&err); */ /* hr += Log_Dnorm_Trunc(cur_t,new_t,(t_max - t_min)/100.,t_min,t_max,&err); */ @@ -8552,7 +8560,7 @@ void MCMC_PHYREX_Node_Times_Pre(t_ldsk *a_ldsk, t_ldsk *d_ldsk, t_tree *tree, in if(failed == NO) { - if(tree->eval_alnL == YES && tree->eval_rlnL == YES) hr += log(r0_new/r0_cur) + log(r1_new/r1_cur) + log(r2_new/r2_cur); + if(tree->eval_rlnL == YES) hr += log(r0_new/r0_cur) + log(r1_new/r1_cur) + log(r2_new/r2_cur); tree->rates->br_r[v0->num] = r0_new; tree->rates->br_r[v1->num] = r1_new; @@ -8622,6 +8630,19 @@ void MCMC_PHYREX_Node_Times_Pre(t_ldsk *a_ldsk, t_ldsk *d_ldsk, t_tree *tree, in cur_tlnL,new_tlnL, alpha); + + if(tree->aux_tree == NULL && (v0 == tree->n_root->v[1] || v0 == tree->n_root->v[2])) + { + PhyML_Printf("\n! cur_t: %f new_t: %f t_min: %f t_max: %f alnL:%f->%f glnL:%f->%f rlnL:%f->%f tlnL:%f->%f hr: %f alpha:%f", + cur_t,new_t,t_min,t_max, + cur_alnL,new_alnL, + cur_glnL,new_glnL, + cur_rlnL,new_rlnL, + cur_tlnL,new_tlnL, + hr, + alpha); + } + assert(isnan(u) == NO && isinf(fabs(u)) == NO); if(u > alpha) /* Reject */ diff --git a/src/phyrex.c b/src/phyrex.c index e0496087..53dabd97 100644 --- a/src/phyrex.c +++ b/src/phyrex.c @@ -559,6 +559,7 @@ void PHYREX_XML(char *xml_filename) PHYREX_Set_Default_Migrep_Mod(mixt_tree->n_otu,aux_tree->mmod); aux_tree->mmod->model_id = mixt_tree->mmod->model_id; aux_tree->mmod->use_locations = mixt_tree->mmod->use_locations; + aux_tree->mmod->integrateAncestralLocations = mixt_tree->mmod->integrateAncestralLocations; Copy_Tree(mixt_tree,aux_tree); diff --git a/src/rates.c b/src/rates.c index d9bd3432..c05c4d2c 100644 --- a/src/rates.c +++ b/src/rates.c @@ -35,36 +35,9 @@ phydbl RATES_Lk(t_tree *tree) err = NO; - /* switch(tree->rates->model_id) */ - /* { */ - /* case LOGNORMAL : */ - /* { */ - /* tree->rates->c_lnL += Log_Dnorm(log(tree->rates->br_r[tree->n_root->num]),-tree->rates->nu*tree->rates->nu/2.,tree->rates->nu,&err); */ - /* tree->rates->c_lnL -= log(tree->rates->br_r[tree->n_root->num]); */ - /* break; */ - /* } */ - /* case GAMMA : */ - /* { */ - /* tree->rates->c_lnL += log(Dgamma(tree->rates->br_r[tree->n_root->num],1./tree->rates->nu,tree->rates->nu)); */ - /* break; */ - /* } */ - /* case STRICTCLOCK : */ - /* { */ - /* tree->rates->c_lnL += 0.0; */ - /* break; */ - /* } */ - /* default : */ - /* { */ - /* PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__); */ - /* Warn_And_Exit(""); */ - /* } */ - /* } */ - - // Priors on clock rate and rate autocorrelation. // Should be separated from likelihood function... - if(tree->rates->clock_r_has_prior == YES) tree->rates->c_lnL += RATES_Clock_R_Prior(tree); - + tree->rates->c_lnL += RATES_Clock_R_Prior(tree); tree->rates->c_lnL += RATES_Autocor_Prior(tree); if(isnan(tree->rates->c_lnL) || err == YES) assert(false); @@ -80,7 +53,7 @@ phydbl RATES_Clock_R_Prior(t_tree *tree) phydbl mean,sd,lnP; int err; - if(tree->rates->clock_r_has_prior == NO) return(UNLIKELY); + if(tree->rates->clock_r_has_prior == NO) return(0.0); err = NO; lnP = 0.0; @@ -120,8 +93,6 @@ phydbl RATES_Autocor_Prior(t_tree *tree) tree->rates->model_id == LOGNORMAL || tree->rates->model_id == GAMMA) lnP += log(lbda) - lbda * tree->rates->nu; - - /* PhyML_Printf("\n. lbda: %f nu: %f lnP: %f",lbda,tree->rates->nu,lnP); */ return(lnP); } diff --git a/src/rrw.c b/src/rrw.c index 08bda3e9..60a295bb 100644 --- a/src/rrw.c +++ b/src/rrw.c @@ -226,7 +226,7 @@ phydbl RRW_Forward_Lk_Path(t_ldsk *a, t_ldsk *d, t_tree *tree) phydbl RRW_Independent_Contrasts(t_tree *tree) { phydbl lnL; - + RRW_Update_Normalization_Factor(tree); RATES_Record_Times(tree); RRW_Rescale_Times(YES,tree); @@ -256,10 +256,14 @@ void RRW_Rescale_Times_Pre(t_node *a, t_node *d, phydbl cur_ta, int prod, t_tree td = tree->times->nd_t[d->num]; ta = tree->times->nd_t[a->num]; + assert(!(cur_ta > td)); + if(prod == YES) tree->times->nd_t[d->num] = ta + (td-cur_ta) * (tree->mmod->sigsq_scale[d->num] * tree->mmod->rrw_norm_fact); else tree->times->nd_t[d->num] = ta + (td-cur_ta) / (tree->mmod->sigsq_scale[d->num] * tree->mmod->rrw_norm_fact); + + assert(!(tree->times->nd_t[d->num] < ta)); if(d->tax == YES) return; else for(i=0;i<3;++i) if(d->v[i] != a && d->b[i] != tree->e_root) RRW_Rescale_Times_Pre(d,d->v[i],td,prod,tree); diff --git a/src/utilities.c b/src/utilities.c index 373813fa..bd226a08 100644 --- a/src/utilities.c +++ b/src/utilities.c @@ -36,9 +36,9 @@ phydbl String_To_Dbl(char *string) if(string == endptr || errno == ERANGE) { - PhyML_Printf("\n. Error in translating string '%s' to double.",string); - PhyML_Printf("\n. %d",errno == ERANGE); - PhyML_Printf("\n. buff = %f",buff); + PhyML_Fprintf(stderr,"\n. Error in translating string '%s' to double.",string); + PhyML_Fprintf(stderr,"\n. %d",errno == ERANGE); + PhyML_Fprintf(stderr,"\n. buff = %f",buff); Generic_Exit(__FILE__,__LINE__,__FUNCTION__); } return buff;