From 2a420d9f533fb7a75f4ca7cee8e74179f2f2adc0 Mon Sep 17 00:00:00 2001 From: Stephane Guindon Date: Tue, 24 Sep 2024 13:50:35 +0200 Subject: [PATCH] Hard-coded exp priors on IOU parameters --- src/init.c | 16 +-- src/iou.c | 29 +++-- src/mcmc.c | 297 ++++++++++++++++++++++++++++++++++++++---------- src/mcmc.h | 1 + src/phyrex.c | 11 +- src/rw.c | 2 +- src/utilities.h | 1 + src/velocity.c | 71 +++++++----- 8 files changed, 314 insertions(+), 114 deletions(-) diff --git a/src/init.c b/src/init.c index f53981fc..7bd398d7 100644 --- a/src/init.c +++ b/src/init.c @@ -3579,7 +3579,7 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t) t->model_id = -1; t->use_locations = -1; t->sampling_scheme = -1; - t->safe_phyrex = NO; + t->safe_phyrex = YES; t->dist_type = HAVERSINE; t->lim_up->lonlat[0] = 100.; @@ -3604,9 +3604,6 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t) t->prior_param_rad = 1.; t->update_rad = NO; - t->min_sigsq = 0.0; - t->max_sigsq = 1.E+3; - t->prior_param_sigsq = 10.0; t->min_veloc = -1.E+5; t->max_veloc = +1.E+5; @@ -3615,16 +3612,19 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t) t->min_omega = 0.0; t->max_omega = 1.E+3; - t->ou_theta = 1.0; - t->min_ou_theta = 1.0E-4; - t->max_ou_theta = 10.0; + t->ou_theta = 1.0E-1; + t->min_ou_theta = 0.0; + t->max_ou_theta = 1000.0; for(int i=0;in_dim;++i) t->ou_mu[i] = 0.0; t->min_ou_mu = -10.; t->max_ou_mu = 10.; + t->min_sigsq = 0.0; + t->max_sigsq = 1.E+3; + t->prior_param_sigsq = 10.0; assert(t->n_dim > 0); - for(int i=0;in_dim;++i) t->sigsq[i] = t->min_sigsq + (t->max_sigsq-t->min_sigsq)/20.; + for(int i=0;in_dim;++i) t->sigsq[i] = 1.E-1; t->nu = 1.0E-0; diff --git a/src/iou.c b/src/iou.c index 09b6e2c1..0266495c 100644 --- a/src/iou.c +++ b/src/iou.c @@ -65,7 +65,8 @@ phydbl IOU_Velocity_Variance_Along_Edge(t_node *d, short int dim, t_tree *tree) { PhyML_Printf("\n. ta: %f t: %f sg2: %f sinh: %f",ta,t,sg2,sinh(ta*t)); } - +// PhyML_Printf("\n. t: %f ta: %f sg2: %f var: %f",t,ta,sg2,var); + assert(isnan(var) == NO); assert(isinf(var) == NO); @@ -115,22 +116,24 @@ phydbl IOU_Location_Mean_Along_Edge(t_node *d, short int dim, t_tree *tree) at both extremities */ phydbl IOU_Location_Variance_Along_Edge(t_node *d, short int dim, t_tree *tree) { - phydbl var,t,ta,sg; + phydbl var,t,ta,sg2; assert(d != tree->n_root); t = fabs(tree->times->nd_t[d->num] - tree->times->nd_t[d->anc->num]); ta = tree->mmod->ou_theta; - sg = tree->mmod->sigsq[dim]; + sg2 = tree->mmod->sigsq[dim] * tree->mmod->sigsq_scale[d->num] * tree->mmod->sigsq_scale_norm_fact; if(t < SMALL) return(0.0); - /* var = sg/pow(ta,3)*(ta*t-2.*(cosh(ta*t)-1.)/sinh(ta*t)); */ - var = sg/pow(ta,3)*(ta*t-2.*(1./tanh(ta*t) - 1./sinh(ta*t))); +// var = sg/pow(ta,3)*(ta*t-2.*(cosh(ta*t)-1.)/sinh(ta*t)); + var = sg2/pow(ta,3)*(ta*t-2.*(1./tanh(ta*t) - 1./sinh(ta*t))); + + if(var > sg2*pow(t,3)/12) var = sg2*pow(t,3)/12; if(var < 0.0) var = 0.0; - + assert(isnan(var) == NO); return(var); @@ -156,9 +159,9 @@ phydbl IOU_Prior_Theta(t_tree *tree) { phydbl lbda,lnP; - lbda = 1.0; - lnP = log(lbda)-lbda*tree->mmod->ou_theta; - + lnP = 0.0; + lbda = 1.0E+4; + lnP += Log_Dexp(tree->mmod->ou_theta,lbda); return(lnP); } @@ -167,7 +170,13 @@ phydbl IOU_Prior_Theta(t_tree *tree) phydbl IOU_Prior_Mu(t_tree *tree) { - return(0.0); + phydbl lbda,lnP; + + lnP = 0.0; + lbda = 1.0E+3; + for (int i = 0; i < tree->mmod->n_dim; i++) + lnP += Log_Dexp(tree->mmod->ou_mu[i], lbda); + return (lnP); } ////////////////////////////////////////////////////////////// diff --git a/src/mcmc.c b/src/mcmc.c index ac66d050..6236ad31 100644 --- a/src/mcmc.c +++ b/src/mcmc.c @@ -3485,7 +3485,7 @@ void MCMC_Randomize_Sigsq(t_tree *tree) } case FLAT_PRIOR: { - tree->mmod->sigsq[i] = 50.; + tree->mmod->sigsq[i] = 1.E-1; break; } default : assert(false); @@ -3518,10 +3518,12 @@ void MCMC_Randomize_Veloc(t_tree *tree) for(j=0;jmmod->n_dim;++j) { - mean = tree->mmod->rw_root_mean[VELOCITY]; - var = tree->mmod->rw_root_var[VELOCITY]; + // mean = tree->mmod->rw_root_mean[VELOCITY]; + // var = tree->mmod->rw_root_var[VELOCITY]; - tree->n_root->ldsk->veloc->deriv[j] = Rnorm(mean,sqrt(var)); + // tree->n_root->ldsk->veloc->deriv[j] = Rnorm(mean,sqrt(var)); + + tree->n_root->ldsk->veloc->deriv[j] = 0.0; // Using a large initial root velocity -> PhyREX cannot sample root velocity properly assert(tree->n_root->ldsk->veloc->deriv[j] > tree->mmod->min_veloc); assert(tree->n_root->ldsk->veloc->deriv[j] < tree->mmod->max_veloc); @@ -6561,7 +6563,8 @@ void MCMC_Complete_MCMC(t_mcmc *mcmc, t_tree *tree) mcmc->num_move_rates_shrink = mcmc->n_moves; mcmc->n_moves += 1; mcmc->num_move_obs_var = mcmc->n_moves; mcmc->n_moves += 1; mcmc->num_move_phyrex_tip_loc = mcmc->n_moves; mcmc->n_moves += 1; - + mcmc->num_move_phyrex_iou_theta_sigsq = mcmc->n_moves; mcmc->n_moves += 1; + mcmc->run_move = (int *)mCalloc(mcmc->n_moves,sizeof(int)); mcmc->acc_move = (int *)mCalloc(mcmc->n_moves,sizeof(int)); mcmc->prev_run_move = (int *)mCalloc(mcmc->n_moves,sizeof(int)); @@ -6659,9 +6662,9 @@ void MCMC_Complete_MCMC(t_mcmc *mcmc, t_tree *tree) strcpy(mcmc->move_name[mcmc->num_move_phyrex_iou_theta],"phyrex_iou_theta"); strcpy(mcmc->move_name[mcmc->num_move_phyrex_iou_mu],"phyrex_iou_mu"); strcpy(mcmc->move_name[mcmc->num_move_rates_shrink],"rates_shrink"); - strcpy(mcmc->move_name[mcmc->num_move_obs_var],"observational_var"); - strcpy(mcmc->move_name[mcmc->num_move_phyrex_tip_loc],"tip_location"); - + strcpy(mcmc->move_name[mcmc->num_move_obs_var], "observational_var"); + strcpy(mcmc->move_name[mcmc->num_move_phyrex_tip_loc], "tip_location"); + strcpy(mcmc->move_name[mcmc->num_move_phyrex_iou_theta_sigsq], "phyrex_iou_theta_sigsq"); for(i=0;in_moves;i++) mcmc->move_type[i] = -1; mcmc->move_type[mcmc->num_move_nd_r] = MCMC_MOVE_SCALE_THORNE; @@ -6799,12 +6802,13 @@ void MCMC_Complete_MCMC(t_mcmc *mcmc, t_tree *tree) mcmc->move_prob[mcmc->num_move_phyrex_node_times] = 4.0; mcmc->move_prob[mcmc->num_move_phyrex_node_veloc] = 4.0; mcmc->move_prob[mcmc->num_move_phyrex_correlated_node_veloc] = 0.0; - mcmc->move_prob[mcmc->num_move_phyrex_all_veloc] = 1.0; + mcmc->move_prob[mcmc->num_move_phyrex_all_veloc] = 2.0; mcmc->move_prob[mcmc->num_move_phyrex_shuffle_node_times] = 1.0; mcmc->move_prob[mcmc->num_move_phyrex_iwn_omega] = 1.0; mcmc->move_prob[mcmc->num_move_phyrex_iou_theta] = 1.0; mcmc->move_prob[mcmc->num_move_phyrex_iou_mu] = 1.0; mcmc->move_prob[mcmc->num_move_phyrex_tip_loc] = 0.0; + mcmc->move_prob[mcmc->num_move_phyrex_iou_theta_sigsq] = 1.0; # else @@ -6849,6 +6853,7 @@ void MCMC_Complete_MCMC(t_mcmc *mcmc, t_tree *tree) mcmc->move_prob[mcmc->num_move_rates_shrink] = 0.0; mcmc->move_prob[mcmc->num_move_obs_var] = 0.0; mcmc->move_prob[mcmc->num_move_phyrex_tip_loc] = 0.0; + mcmc->move_prob[mcmc->num_move_phyrex_iou_theta_sigsq] = 0.0; #endif @@ -7022,6 +7027,7 @@ void MCMC_Run_Core(t_tree *tree) if(!strcmp(tree->mcmc->move_name[move],"clock")) MCMC_Clock_R(tree); if(!strcmp(tree->mcmc->move_name[move],"phyrex_iwn_omega")) MCMC_PHYREX_IWN_Update_Omega(tree); if(!strcmp(tree->mcmc->move_name[move],"phyrex_iou_theta")) MCMC_PHYREX_IOU_Update_Theta(tree); + if(!strcmp(tree->mcmc->move_name[move],"phyrex_iou_theta_sigsq")) MCMC_PHYREX_IOU_Update_Theta_Sigsq(tree); if(!strcmp(tree->mcmc->move_name[move],"phyrex_iou_mu")) MCMC_PHYREX_IOU_Update_Mu(tree); if(!strcmp(tree->mcmc->move_name[move],"updown_t_br")) MCMC_Updown_T_Br(tree); if(!strcmp(tree->mcmc->move_name[move],"observational_var")) MCMC_Obs_Var(tree); @@ -9324,7 +9330,7 @@ void MCMC_PHYREX_Node_Times(t_tree *tree, int print) Lk(NULL,tree); } - if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && tree->eval_glnL == YES) + if(tree->eval_glnL == YES) { tree->contmod->both_sides[LOCATION] = YES; tree->contmod->both_sides[VELOCITY] = YES; @@ -9354,13 +9360,14 @@ void MCMC_PHYREX_Node_Times(t_tree *tree, int print) if(tree->eval_alnL == YES) { Set_Both_Sides(NO,tree); - Lk(NULL,tree); /* Required in order to rescale all branch rates */ + Lk(NULL,tree); /* Required in order to rescale all branch rates */ } - if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && tree->eval_glnL == YES) + if(tree->eval_glnL == YES) { tree->contmod->both_sides[LOCATION] = NO; tree->contmod->both_sides[VELOCITY] = NO; + LOCATION_Lk(NULL,tree); } #ifdef PHYREX @@ -9468,7 +9475,7 @@ void MCMC_PHYREX_Node_Times_Pre(t_ldsk *a_ldsk, t_ldsk *d_ldsk, t_tree *tree, in new_glnL = (tree->mmod->use_locations == YES) ? LOCATION_Lk(d_ldsk->nd,tree) : 0.0; // new_glnL = (tree->mmod->use_locations == YES) ? LOCATION_Lk(NULL,tree) : 0.0; /* !!!!!!!!!!!!!!!!!! */ } - else + else // NOTE: should use node lk calculation for RRW here as well.... { new_glnL = (tree->mmod->use_locations == YES) ? LOCATION_Lk(NULL,tree) : 0.0; } @@ -13944,8 +13951,9 @@ void MCMC_PHYREX_Node_Velocity(t_tree *tree) MCMC_PHYREX_Node_Velocity_Pre(NULL,tree->n_root,tree); tree->contmod->both_sides[LOCATION] = NO; tree->contmod->both_sides[VELOCITY] = NO; + LOCATION_Lk(NULL,tree); -// /* Velocities */ + /* Velocities */ // int *permut; // int i, j; @@ -13987,9 +13995,9 @@ void MCMC_PHYREX_Node_Velocity(t_tree *tree) // /* hr -= Log_Dnorm(new,cur,0.1,&err); */ // /* hr += Log_Dnorm(cur,new,0.1,&err); */ -// new = Rnorm(cur, 1.E-0); -// hr -= Log_Dnorm(new, cur, 1.E-0, &err); -// hr += Log_Dnorm(cur, new, 1.E-0, &err); +// new = Rnorm(cur, 1.E-2); +// hr -= Log_Dnorm(new, cur, 1.E-2, &err); +// hr += Log_Dnorm(cur, new, 1.E-2, &err); // /* if(n->ldsk->prev != NULL) */ // /* { */ @@ -14129,12 +14137,15 @@ void MCMC_PHYREX_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree) if(dad != tree->n_root) { var_dad = VELOC_Velocity_Variance_Along_Edge(dad,i,tree); - mu_dad = dad->anc->ldsk->veloc->deriv[i]; + // mu_dad = dad->anc->ldsk->veloc->deriv[i]; + mu_dad = VELOC_Velocity_Mean_Along_Edge(dad,i,tree); } else { + // var_dad = tree->mmod->rw_root_var[VELOCITY]; + // mu_dad = tree->mmod->rw_root_mean[VELOCITY]; var_dad = tree->mmod->rw_root_var[VELOCITY]; - mu_dad = tree->mmod->rw_root_mean[VELOCITY]; + mu_dad = dad->ldsk->veloc->deriv[i]; } if(dad->tax == YES) @@ -14147,8 +14158,10 @@ void MCMC_PHYREX_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree) var_son = VELOC_Velocity_Variance_Along_Edge(son,i,tree); var_bro = VELOC_Velocity_Variance_Along_Edge(bro,i,tree); - mu_son = son->ldsk->veloc->deriv[i]; - mu_bro = bro->ldsk->veloc->deriv[i]; + // mu_son = son->ldsk->veloc->deriv[i]; + // mu_bro = bro->ldsk->veloc->deriv[i]; + mu_son = VELOC_Velocity_Mean_Along_Edge(son,i,tree); + mu_bro = VELOC_Velocity_Mean_Along_Edge(bro,i,tree); if(var_dad > SMALL && var_son > SMALL && var_bro > SMALL) { @@ -14182,18 +14195,16 @@ void MCMC_PHYREX_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree) if (isnan(var)) assert(false); - new_vel = Rnorm(mu, 1.5 * sqrt(var)); - hr -= Log_Dnorm(new_vel, mu, 1.5 * sqrt(var), &err); - hr += Log_Dnorm(cur_vel, mu, 1.5 * sqrt(var), &err); - - // if(dad == tree->n_root) - // { - // PhyML_Printf("\n. mu: %g var: %g v1: %g v2: %g new_vel: %g cur_vel: %g hr: %g", - // mu,var, - // var_son,var_bro, - // new_vel,cur_vel, - // hr); - // } + if (dad == tree->n_root) + { + new_vel = Rnorm(cur_vel,1.E-2); + } + else + { + new_vel = Rnorm(mu, 1. * sqrt(var)); + hr -= Log_Dnorm(new_vel, mu, 1. * sqrt(var), &err); + hr += Log_Dnorm(cur_vel, mu, 1. * sqrt(var), &err); + } d->ldsk->veloc->deriv[i] = new_vel; } @@ -14213,24 +14224,6 @@ void MCMC_PHYREX_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree) new_glnL = LOCATION_Lk(d,tree); // new_glnL = LOCATION_Lk(NULL,tree); /* !!!!!!!!!!!!!!!!! */ -// PhyML_Printf("\n. a: %3d d: %3d (%d%d%d) hr: %13f cur_glnL: %13f new_glnL: %13f vel: %12f %12f lnLv(%13f %13f) lnLl(%13f %13f) VEL down:(%13f %13f) up:(%13f %13f) dt: %f", -// a ? a->num : -1, d->num, -// d == tree->n_root ? 1 : 0, -// d == tree->n_root->v[1] ? 1 : 0, -// d == tree->n_root->v[2] ? 1 : 0, -// hr, -// cur_glnL, new_glnL, -// d->ldsk->veloc->deriv[0], -// d->ldsk->veloc->deriv[1], -// tree->contmod->lnL[VELOCITY * tree->mmod->n_dim + 0], -// tree->contmod->lnL[VELOCITY * tree->mmod->n_dim + 1], -// tree->contmod->lnL[LOCATION * tree->mmod->n_dim + 0], -// tree->contmod->lnL[LOCATION * tree->mmod->n_dim + 1], -// tree->contmod->lnL_down[Contmod_Start(VELOCITY, 0, tree) + d->num], -// tree->contmod->lnL_down[Contmod_Start(VELOCITY, 1, tree) + d->num], -// tree->contmod->lnL_up[Contmod_Start(VELOCITY, 0, tree) + d->num], -// tree->contmod->lnL_up[Contmod_Start(VELOCITY, 1, tree) + d->num], -// a ? (d->ldsk->disk->time - a->ldsk->disk->time) : -1.); ratio = 0.0; ratio += (new_glnL - cur_glnL); @@ -14239,6 +14232,26 @@ void MCMC_PHYREX_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree) ratio = exp(ratio); alpha = MIN(1.,ratio); +// if (dad == tree->n_root) +// PhyML_Printf("\n. a: %3d d: %3d (%d%d%d) hr: %13f alpha: %13f cur_glnL: %13f new_glnL: %13f vel: %12f %12f lnLv(%13f %13f) lnLl(%13f %13f) VEL down:(%13f %13f) up:(%13f %13f) dt: %f", +// a ? a->num : -1, d->num, +// d == tree->n_root ? 1 : 0, +// d == tree->n_root->v[1] ? 1 : 0, +// d == tree->n_root->v[2] ? 1 : 0, +// hr, +// alpha, +// cur_glnL, new_glnL, +// d->ldsk->veloc->deriv[0], +// d->ldsk->veloc->deriv[1], +// tree->contmod->lnL[VELOCITY * tree->mmod->n_dim + 0], +// tree->contmod->lnL[VELOCITY * tree->mmod->n_dim + 1], +// tree->contmod->lnL[LOCATION * tree->mmod->n_dim + 0], +// tree->contmod->lnL[LOCATION * tree->mmod->n_dim + 1], +// tree->contmod->lnL_down[Contmod_Start(VELOCITY, 0, tree) + d->num], +// tree->contmod->lnL_down[Contmod_Start(VELOCITY, 1, tree) + d->num], +// tree->contmod->lnL_up[Contmod_Start(VELOCITY, 0, tree) + d->num], +// tree->contmod->lnL_up[Contmod_Start(VELOCITY, 1, tree) + d->num], +// a ? (d->ldsk->disk->time - a->ldsk->disk->time) : -1.); u = Uni(); @@ -14491,6 +14504,7 @@ void MCMC_PHYREX_All_Velocities(t_tree *tree) phydbl cur_glnP, new_glnP; t_node *n; phydbl cur_sigsq, new_sigsq; + phydbl cur_iou_theta, new_iou_theta; /* Scale all velocities */ @@ -14514,6 +14528,9 @@ void MCMC_PHYREX_All_Velocities(t_tree *tree) cur_sigsq = tree->mmod->sigsq[j]; new_sigsq = tree->mmod->sigsq[j]; + cur_iou_theta = tree->mmod->ou_theta; + new_iou_theta = tree->mmod->ou_theta; + n = tree->n_root; u = Uni(); @@ -14521,7 +14538,12 @@ void MCMC_PHYREX_All_Velocities(t_tree *tree) n_nodes = 0; Scale_Subtree_Veloc(n, mult, &n_nodes, j, tree); - /* Add_Subtree_Veloc(n,new-cur,&n_nodes,j,tree); */ + + if (IOU_Is_Iou(tree->mmod) == YES) + { + new_iou_theta = cur_iou_theta * mult; + tree->mmod->ou_theta = new_iou_theta; + } new_sigsq = cur_sigsq * mult; tree->mmod->sigsq[j] = new_sigsq; @@ -14534,6 +14556,7 @@ void MCMC_PHYREX_All_Velocities(t_tree *tree) ratio += (new_glnP - cur_glnP); ratio += n_nodes * log(mult); ratio += log(mult); + if(IOU_Is_Iou(tree->mmod) == YES) ratio += log(mult); ratio = exp(ratio); alpha = MIN(1., ratio); @@ -14546,6 +14569,7 @@ void MCMC_PHYREX_All_Velocities(t_tree *tree) { PHYREX_Restore_All_Veloc(tree); tree->mmod->sigsq[j] = cur_sigsq; + tree->mmod->ou_theta = cur_iou_theta; Reset_Lk(tree); Reset_Prior(tree); } @@ -14588,16 +14612,86 @@ void MCMC_PHYREX_IWN_Update_Omega(t_tree *tree) #ifdef PHYREX void MCMC_PHYREX_IOU_Update_Theta(t_tree *tree) { - if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES) + if (VELOC_Is_Integrated_Velocity(tree->mmod) == YES && IOU_Is_Iou(tree->mmod) == YES) { - MCMC_Single_Param_Generic(&(tree->mmod->ou_theta), - tree->mmod->min_ou_theta, - tree->mmod->max_ou_theta, - tree->mcmc->num_move_phyrex_iou_theta, - &(tree->mmod->c_lnP),&(tree->mmod->c_lnL), - LOCATION_Wrap_Prior,LOCATION_Wrap_Lk, - tree->mcmc->move_type[tree->mcmc->num_move_phyrex_iou_theta], - NO,NO,NULL,tree,NULL); + // MCMC_Single_Param_Generic(&(tree->mmod->ou_theta), + // tree->mmod->min_ou_theta, + // tree->mmod->max_ou_theta, + // tree->mcmc->num_move_phyrex_iou_theta, + // &(tree->mmod->c_lnP),&(tree->mmod->c_lnL), + // LOCATION_Wrap_Prior,LOCATION_Wrap_Lk, + // tree->mcmc->move_type[tree->mcmc->num_move_phyrex_iou_theta], + // NO,NO,NULL,tree,NULL); + + int j; + phydbl u, alpha, ratio, mult; + phydbl cur_glnL, new_glnL; + phydbl cur_glnP, new_glnP; + phydbl cur_theta, new_theta; + + // LOCATION_Lk(NULL, tree); + // LOCATION_Prior(tree); + + for (j = 0; j < tree->mmod->n_dim; ++j) + { + tree->mcmc->run_move[tree->mcmc->num_move_phyrex_iou_theta]++; + + Set_Lk(tree); + Set_Prior(tree); + + cur_glnL = tree->mmod->c_lnL; + new_glnL = UNLIKELY; + + cur_glnP = tree->mmod->c_lnP; + new_glnP = UNLIKELY; + + cur_theta = tree->mmod->ou_theta; + new_theta = tree->mmod->ou_theta; + + u = Uni(); + mult = exp((u - 0.5)); + + new_theta = cur_theta * mult; + tree->mmod->ou_theta = new_theta; + + new_glnL = LOCATION_Lk(NULL, tree); + new_glnP = LOCATION_Prior(tree); + + ratio = 0.0; + ratio += (new_glnL - cur_glnL); + ratio += (new_glnP - cur_glnP); + ratio += 1 * log(mult); + + ratio = exp(ratio); + alpha = MIN(1., ratio); + + // PhyML_Printf("\n. lnL: %13f->%13f scale: %13f alpha: %13f theta: %13g->%13g [%13g %13g %13g %13g]", + // cur_glnL, new_glnL, + // mult, alpha, + // cur_theta, new_theta, + // tree->contmod->lnL[VELOCITY*tree->mmod->n_dim+0], + // tree->contmod->lnL[VELOCITY*tree->mmod->n_dim+1], + // tree->contmod->lnL[LOCATION*tree->mmod->n_dim+0], + // tree->contmod->lnL[LOCATION*tree->mmod->n_dim+1]); + + u = Uni(); + + if (u > alpha) /* Reject */ + { + tree->mmod->ou_theta = cur_theta; + Reset_Lk(tree); + Reset_Prior(tree); + } + else + { + tree->mcmc->acc_move[tree->mcmc->num_move_phyrex_iou_theta]++; + } + + tree->mcmc->run++; + + PHYREX_Print_MCMC_Tree(tree); + PHYREX_Print_MCMC_Summary(tree); + } } } @@ -14609,7 +14703,7 @@ void MCMC_PHYREX_IOU_Update_Theta(t_tree *tree) #ifdef PHYREX void MCMC_PHYREX_IOU_Update_Mu(t_tree *tree) { - if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES) + if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && IOU_Is_Iou(tree->mmod) == YES) { for(int i=0; immod->n_dim; ++i) MCMC_Single_Param_Generic(&(tree->mmod->ou_mu[i]), @@ -14628,6 +14722,85 @@ void MCMC_PHYREX_IOU_Update_Mu(t_tree *tree) /*//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////*/ +#ifdef PHYREX +void MCMC_PHYREX_IOU_Update_Theta_Sigsq(t_tree *tree) +{ + if (VELOC_Is_Integrated_Velocity(tree->mmod) == YES && IOU_Is_Iou(tree->mmod) == YES) + { + int j; + phydbl u, alpha, ratio, mult; + phydbl cur_glnL, new_glnL; + phydbl cur_glnP, new_glnP; + phydbl cur_sigsq, new_sigsq; + phydbl cur_theta, new_theta; + + for (j = 0; j < tree->mmod->n_dim; ++j) + { + tree->mcmc->run_move[tree->mcmc->num_move_phyrex_iou_theta_sigsq]++; + + Set_Lk(tree); + Set_Prior(tree); + + cur_glnL = tree->mmod->c_lnL; + new_glnL = UNLIKELY; + + cur_glnP = tree->mmod->c_lnP; + new_glnP = UNLIKELY; + + cur_sigsq = tree->mmod->sigsq[j]; + new_sigsq = tree->mmod->sigsq[j]; + + cur_theta = tree->mmod->ou_theta; + new_theta = tree->mmod->ou_theta; + + u = Uni(); + mult = exp((u - 0.5)); + + new_sigsq = cur_sigsq * mult; + tree->mmod->sigsq[j] = new_sigsq; + + new_theta = cur_theta * mult; + tree->mmod->ou_theta = new_theta; + + new_glnL = LOCATION_Lk(NULL, tree); + new_glnP = LOCATION_Prior(tree); + + ratio = 0.0; + ratio += (new_glnL - cur_glnL); + ratio += (new_glnP - cur_glnP); + ratio += 2 * log(mult); + + ratio = exp(ratio); + alpha = MIN(1., ratio); + + // PhyML_Printf("\n. lnL: %f->%f scale: %f alpha: %f",cur_glnL,new_glnL,mult,alpha); + + u = Uni(); + + if (u > alpha) /* Reject */ + { + tree->mmod->sigsq[j] = cur_sigsq; + tree->mmod->ou_theta = cur_theta; + Reset_Lk(tree); + Reset_Prior(tree); + } + else + { + tree->mcmc->acc_move[tree->mcmc->num_move_phyrex_iou_theta_sigsq]++; + } + + tree->mcmc->run++; + PHYREX_Print_MCMC_Stats(tree); + PHYREX_Print_MCMC_Tree(tree); + PHYREX_Print_MCMC_Summary(tree); + } + } +} +#endif + +/*//////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////*/ + #ifdef PHYREX void MCMC_Obs_Var(t_tree *tree) { diff --git a/src/mcmc.h b/src/mcmc.h index c2592f65..0f97d8cf 100644 --- a/src/mcmc.h +++ b/src/mcmc.h @@ -206,5 +206,6 @@ void MCMC_Crossvalidate_Locations(t_tree *tree); void MCMC_Randomize_Sigsq(t_tree *tree); void MCMC_PHYREX_Correlated_Node_Velocity(t_tree *tree); void MCMC_PHYREX_Correlated_Node_Velocity_Pre(t_node *a, t_node *d, t_tree *tree); +void MCMC_PHYREX_IOU_Update_Theta_Sigsq(t_tree *tree); #endif diff --git a/src/phyrex.c b/src/phyrex.c index 5b2dcba0..b8024eb3 100644 --- a/src/phyrex.c +++ b/src/phyrex.c @@ -896,7 +896,13 @@ void PHYREX_XML(char *xml_filename) MCMC_Randomize_Locations(mixt_tree); MCMC_Randomize_Veloc(mixt_tree); MCMC_Randomize_Contmod(mixt_tree); - + +// mixt_tree->mmod->ou_theta = 0.1; +// mixt_tree->mmod->ou_mu[0] = 0.1; +// mixt_tree->mmod->ou_mu[1] = 0.1; +// mixt_tree->mmod->sigsq[0] = 1.E-3; +// mixt_tree->mmod->sigsq[1] = 1.E-3; + Set_Ignore_Root(YES,mixt_tree); Set_Bl_From_Rt(YES,mixt_tree); @@ -921,8 +927,7 @@ void PHYREX_XML(char *xml_filename) Update_Ancestors(mixt_tree->aux_tree[i]->n_root,mixt_tree->aux_tree[i]->n_root->v[1],mixt_tree->aux_tree[i]->n_root->b[1],mixt_tree->aux_tree[i]); } } - - + assert(PHYREX_Check_Struct(mixt_tree,YES)); TIMES_Lk(mixt_tree); RATES_Lk(mixt_tree); diff --git a/src/rw.c b/src/rw.c index 245d47e1..d1428afc 100644 --- a/src/rw.c +++ b/src/rw.c @@ -237,7 +237,7 @@ void RW_Integrated_Lk_Down(phydbl son_a, phydbl son_b, phydbl son_mu_down, phydb logr = son_logrem + bro_logrem; logr -= log(fabs(bro_a*son_a)); logr += Log_Dnorm((son_mu_down-son_b)/son_a,(bro_mu_down-bro_b)/bro_a,sqrt((son_var_down+son_var)/pow(son_a,2)+(bro_var_down+bro_var)/pow(bro_a,2)),&err); - + if((son_var + son_var_down > SMALL) && (bro_var + bro_var_down > SMALL)) // Standard case { v = pow(son_a,2)/(son_var_down + son_var) + pow(bro_a,2)/(bro_var_down + bro_var); diff --git a/src/utilities.h b/src/utilities.h index 742cc589..1d899dc7 100644 --- a/src/utilities.h +++ b/src/utilities.h @@ -1872,6 +1872,7 @@ typedef struct __Tmcmc { int num_move_phyrex_shuffle_node_times; int num_move_phyrex_iwn_omega; int num_move_phyrex_iou_theta; + int num_move_phyrex_iou_theta_sigsq; int num_move_phyrex_iou_mu; int num_move_obs_var; int num_move_phyrex_tip_loc; diff --git a/src/velocity.c b/src/velocity.c index e7b3eb51..efcee5f6 100644 --- a/src/velocity.c +++ b/src/velocity.c @@ -461,15 +461,13 @@ phydbl VELOC_Integrated_Lk_Location(t_node *z, t_tree *tree) start = Contmod_Start(LOCATION,i,tree); - tree->contmod->lnL[LOCATION*tree->mmod->n_dim+i] = 0.0; tree->contmod->lnL[LOCATION*tree->mmod->n_dim+i] += tree->contmod->logrem_down[start + tree->n_root->num]; tree->contmod->lnL[LOCATION*tree->mmod->n_dim+i] += Log_Dnorm(tree->contmod->mu_down[start + tree->n_root->num], root_mean, sqrt(root_var+tree->contmod->var_down[start + tree->n_root->num]), &err); - - lnL += tree->contmod->lnL[LOCATION*tree->mmod->n_dim+i]; + lnL += tree->contmod->lnL[LOCATION * tree->mmod->n_dim + i]; } } else @@ -591,7 +589,6 @@ void VELOC_Update_Lk_Location_Up(t_node *a, t_node *d, t_tree *tree) phydbl dad_mu_up, dad_var_up, dad_logrem_up; phydbl son_a, son_b, son_var; phydbl bro_a, bro_b, bro_mu_down, bro_var_down, bro_var, bro_logrem_down; - phydbl dt_son, dt_bro; dad = a; son = d; @@ -607,8 +604,6 @@ void VELOC_Update_Lk_Location_Up(t_node *a, t_node *d, t_tree *tree) assert(bro->anc == dad); assert(bro != son); - dt_son = tree->times->nd_t[son->num] - tree->times->nd_t[dad->num]; - dt_bro = tree->times->nd_t[bro->num] - tree->times->nd_t[dad->num]; for(i=0;immod->n_dim;++i) { @@ -669,7 +664,7 @@ void VELOC_Update_Lk_Location_Up(t_node *a, t_node *d, t_tree *tree) assert(isnan(tree->contmod->logrem_up[start+d->num]) == NO); assert(isnan(tree->contmod->mu_up[start+d->num]) == NO); - assert(isnan(tree->contmod->var_up[start+d->num]) == NO && !(tree->contmod->var_down[start+d->num] < 0.0)); + assert(isnan(tree->contmod->var_up[start+d->num]) == NO && !(tree->contmod->var_up[start+d->num] < 0.0)); } @@ -688,7 +683,6 @@ void VELOC_Update_Lk_Location_Down(t_node *a, t_node *d, t_tree *tree) phydbl son_logrem_down,bro_logrem_down; phydbl son_a, bro_a; phydbl son_b, bro_b; - phydbl dt_son, dt_bro; if(d->tax) return; @@ -703,8 +697,6 @@ void VELOC_Update_Lk_Location_Down(t_node *a, t_node *d, t_tree *tree) } } - dt_son = tree->times->nd_t[son->num] - tree->times->nd_t[dad->num]; - dt_bro = tree->times->nd_t[bro->num] - tree->times->nd_t[dad->num]; for(i=0;immod->n_dim;++i) { @@ -741,15 +733,8 @@ void VELOC_Update_Lk_Location_Down(t_node *a, t_node *d, t_tree *tree) } else if(IWN_Is_Iwn(tree->mmod) == YES) { - IWN_Integrated_Location_Down(dt_son,dt_bro, - son_a,son_b,son_mu_down,son_var_down,son_var, - bro_a,bro_b,bro_mu_down,bro_var_down,bro_var, - son_logrem_down,bro_logrem_down, - d->ldsk->veloc->deriv[i],son->ldsk->veloc->deriv[i],bro->ldsk->veloc->deriv[i], - tree->mmod->omega, - tree->contmod->mu_down+start+dad->num, - tree->contmod->var_down+start+dad->num, - tree->contmod->logrem_down+start+dad->num); + PhyML_Printf("\n. Not implemented yet"); + assert(false); } else if(IOU_Is_Iou(tree->mmod) == YES) { @@ -949,18 +934,32 @@ void VELOC_Sample_Node_Locations_Joint(t_tree *tree) for(i=0;immod->n_dim;++i) { - start = Contmod_Start(LOCATION,i,tree); - - var = 1./tree->contmod->var_down[start+tree->n_root->num] + 1./root_var; - var = 1./var; - - root_mean = LOCATION_Mean_Lonlat(i,tree); + start = Contmod_Start(LOCATION, i, tree); - mean = (tree->contmod->mu_down[start+tree->n_root->num]/tree->contmod->var_down[start+tree->n_root->num] + - root_mean / root_var)*var; + root_mean = LOCATION_Mean_Lonlat(i, tree); + + if (tree->contmod->var_down[start + tree->n_root->num] > SMALL) + { + var = 1. / tree->contmod->var_down[start + tree->n_root->num] + 1. / root_var; + var = 1. / var; + mean = (tree->contmod->mu_down[start + tree->n_root->num] / tree->contmod->var_down[start + tree->n_root->num] + + root_mean / root_var) * + var; + } + else + { + var = 0.0; + mean = tree->contmod->mu_down[start + tree->n_root->num]; + } tree->n_root->ldsk->coord->lonlat[i] = Rnorm(mean,sqrt(var)); + assert(isnan(mean) == NO); + assert(isnan(var) == NO); + assert(isinf(mean) == NO); + assert(isinf(var) == NO); + + assert(isnan(tree->contmod->var_down[start + tree->n_root->num]) == NO); } @@ -980,17 +979,20 @@ void VELOC_Sample_Node_Locations_Joint_Post(t_node *a, t_node *d, t_tree *tree) { int start, i; phydbl son_a, son_b, son_var, son_var_down, son_mu, son_mu_down, var, mean; + t_node *son, *dad; for (int i = 0; i < 3; ++i) if (d->v[i] != a && !(a == tree->n_root && d->b[i] == tree->e_root)) VELOC_Sample_Node_Locations_Joint_Post(d, d->v[i], tree); + son = d; + dad = a; for (i = 0; i < tree->mmod->n_dim; ++i) { start = Contmod_Start(LOCATION, i, tree); son_a = 1.0; - son_b = 0.0; + son_b = VELOC_Location_Mean_Along_Edge(son, i, tree) - son_a * dad->ldsk->coord->lonlat[i]; son_var = VELOC_Location_Variance_Along_Edge(d, i, tree); son_var_down = tree->contmod->var_down[start + d->num]; @@ -1022,6 +1024,11 @@ void VELOC_Sample_Node_Locations_Joint_Post(t_node *a, t_node *d, t_tree *tree) } assert(!(var < 0.0)); + assert(isnan(mean) == NO); + assert(isnan(var) == NO); + assert(isinf(mean) == NO); + assert(isinf(var) == NO); + d->ldsk->coord->lonlat[i] = Rnorm(mean, sqrt(var)); } } @@ -1057,8 +1064,12 @@ phydbl VELOC_Mean_Speed(t_tree *tree) phydbl mean_speed; mean_speed = 0.0; - for(i=0;i<2*tree->n_otu-1;++i) mean_speed += VELOC_Veloc_To_Speed(tree->a_nodes[i]->ldsk->veloc,tree); - return(mean_speed / (phydbl)(2*tree->n_otu-1)); + // for(i=0;i<2*tree->n_otu-1;++i) mean_speed += VELOC_Veloc_To_Speed(tree->a_nodes[i]->ldsk->veloc,tree); + // return(mean_speed / (phydbl)(2*tree->n_otu-1)); + // for(i=0;i<2*tree->n_otu-1;++i) if(tree->a_nodes[i] != tree->n_root) mean_speed += VELOC_Veloc_To_Speed(tree->a_nodes[i]->ldsk->veloc,tree); + // return(mean_speed / (phydbl)(2*tree->n_otu-2)); + for(i=0;in_otu;++i) mean_speed += VELOC_Veloc_To_Speed(tree->a_nodes[i]->ldsk->veloc,tree); + return(mean_speed / (phydbl)(tree->n_otu)); } return(-1.); }