Skip to content

Commit

Permalink
New release
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Apr 8, 2022
1 parent 7ec459e commit 1d08eee
Show file tree
Hide file tree
Showing 8 changed files with 89 additions and 79 deletions.
2 changes: 1 addition & 1 deletion configure.ac
View file
Original file line number Diff line number Diff line change
@@ -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'"]),[[email protected]])
AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20220408" | tr -d '\n'"]),[[email protected]])
AM_INIT_AUTOMAKE([foreign])
AC_CONFIG_SRCDIR([src/main.c],[doc/phyml-manual.tex])
AC_CONFIG_HEADERS([config.h])
Expand Down
9 changes: 3 additions & 6 deletions src/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; */

Expand Down Expand Up @@ -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;

Expand Down
16 changes: 16 additions & 0 deletions src/io.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;i<tree->mmod->n_dim;++i) PhyML_Fprintf(fp_stats,"%s%s\t","sigSq",(i==0)?("Lon"):((i==1)?("Lat"):("xx")));
Expand All @@ -6675,6 +6681,8 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
/* for(int i=0;i<tree->mod->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"); */

Expand Down Expand Up @@ -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;i<tree->mmod->n_dim;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->sigsq[i]);
Expand All @@ -6760,6 +6774,8 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
/* for(int i=0;i<tree->mod->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); */
Expand Down
95 changes: 58 additions & 37 deletions src/mcmc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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.;
Expand Down Expand Up @@ -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;

Expand All @@ -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;
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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); */
Expand All @@ -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;
Expand Down Expand Up @@ -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 */
Expand Down
1 change: 1 addition & 0 deletions src/phyrex.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
33 changes: 2 additions & 31 deletions src/rates.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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;
Expand Down Expand Up @@ -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);
}
Expand Down
6 changes: 5 additions & 1 deletion src/rrw.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
6 changes: 3 additions & 3 deletions src/utilities.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 1d08eee

Please sign in to comment.