From efd0397ac2aacde635b256bb14cdfd648f75b973 Mon Sep 17 00:00:00 2001 From: Stephane Guindon Date: Wed, 10 Jul 2024 08:33:40 +0200 Subject: [PATCH] Put scaling factors back on --- src/cl.c | 8 +-- src/init.c | 4 +- src/lk.c | 35 ++++++++--- src/main.c | 2 +- src/models.c | 162 +++++++++++++++--------------------------------- src/models.h | 4 +- src/optimiz.c | 6 +- src/utilities.c | 17 +++-- 8 files changed, 94 insertions(+), 144 deletions(-) diff --git a/src/cl.c b/src/cl.c index ec5246d1..86cc9778 100644 --- a/src/cl.c +++ b/src/cl.c @@ -396,10 +396,10 @@ int Read_Command_Line(option *io, int argc, char **argv) } case 66: { - PhyML_Printf("\n. Integrated edge length model needs additional work."); - PhyML_Printf("\n. Please send me an email at guindon@lirmm.fr if you are"); - PhyML_Printf("\n. really keen on using this model."); - assert(false); + /* PhyML_Printf("\n. Integrated edge length model needs additional work."); */ + /* PhyML_Printf("\n. Please send me an email at guindon@lirmm.fr if you are"); */ + /* PhyML_Printf("\n. really keen on using this model."); */ + /* assert(false); */ io->mod->gamma_mgf_bl = YES; io->mod->s_opt->opt_gamma_br_len = YES; break; diff --git a/src/init.c b/src/init.c index 90d97541..f420ce3f 100644 --- a/src/init.c +++ b/src/init.c @@ -148,7 +148,7 @@ void Init_Tree(t_tree *tree, int n_otu) tree->eval_rlnL = YES; tree->eval_glnL = YES; tree->eval_tlnL = YES; - tree->scaling_method = NO; + tree->scaling_method = SCALE_FAST; tree->perform_spr_right_away = YES; tree->tip_root = 0; tree->n_edges_traversed = 0; @@ -682,7 +682,7 @@ void Set_Defaults_Model(t_mod *mod) mod->kappa->v = 4.0; mod->lambda->v = 1.0; - mod->l_var_sigma = 1.E-2; + mod->l_var_sigma = 1.E-1; mod->e_frq_weight->v = 1.0; mod->r_mat_weight->v = 1.0; diff --git a/src/lk.c b/src/lk.c index 38b98ecc..b6ceb4ca 100644 --- a/src/lk.c +++ b/src/lk.c @@ -598,7 +598,7 @@ phydbl Lk(t_edge *b, t_tree *tree) if(tree->mixt_tree != NULL) len *= tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number]; if(len < tree->mod->l_min) len = tree->mod->l_min; else if(len > tree->mod->l_max) len = tree->mod->l_max; - for(state=0;statemod->eigen->e_val[state],len); + for(state=0;statemod->eigen->e_val[state]*len); } } @@ -655,7 +655,7 @@ phydbl Lk(t_edge *b, t_tree *tree) phydbl dLk(phydbl *l, t_edge *b, t_tree *tree) { unsigned int catg,state,site; - phydbl len,rr; + phydbl len,rr,var; phydbl lk,dlk,dlnlk,lnlk; phydbl ev,expevlen; @@ -691,7 +691,8 @@ phydbl dLk(phydbl *l, t_edge *b, t_tree *tree) rr *= tree->mod->br_len_mult->v; len = (*l) * rr; - + var = tree->mod->l_var_sigma * rr*rr; + if(isinf(len) || isnan(len)) { PhyML_Fprintf(stderr,"\n. len=%f rr=%f l=%f",len,rr,*l); @@ -709,8 +710,17 @@ phydbl dLk(phydbl *l, t_edge *b, t_tree *tree) ev = tree->mod->eigen->e_val[state]; expevlen = exp(ev*len); - expl[catg*2*ns + 2*state] = expevlen; - expl[catg*2*ns + 2*state + 1] = expevlen*ev*rr; + if(tree->io->mod->gamma_mgf_bl == YES) + { + expl[catg*2*ns + 2*state] = POW(1. - ev*var/len,-len*len/var); + expl[catg*2*ns + 2*state + 1] = expl[catg*2*ns + 2*state]; + expl[catg*2*ns + 2*state + 1] *= -(ev/(1.-ev*var/len) + 2.*len*LOG(1.-ev*var/len)/var); + } + else + { + expl[catg*2*ns + 2*state] = expevlen; + expl[catg*2*ns + 2*state + 1] = expevlen*ev*rr; + } } } @@ -2163,7 +2173,7 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree) int i; phydbl len; phydbl l_min, l_max; - phydbl shape, scale, mean, var; + phydbl mean, var; assert(b_fcus); assert(tree); @@ -2223,7 +2233,8 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree) else if(len > l_max) len = l_max; mean = len; - var = MAX(0.0,b_fcus->l_var->v) * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2); + /* var = MAX(0.0,b_fcus->l_var->v) * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2); */ + var = tree->mod->l_var_sigma * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2); if(tree->mixt_tree) var *= POW(tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number],2); /* var = 1.E-10; */ @@ -2247,9 +2258,8 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree) Warn_And_Exit(TODO_BEAGLE); #endif - shape = mean*mean/var; - scale = var/mean; - PMat_MGF_Gamma(b_fcus->Pij_rr+tree->mod->ns*tree->mod->ns*i,shape,scale,1.0,tree->mod); + /* PMat_MGF_Gamma(b_fcus->Pij_rr+tree->mod->ns*tree->mod->ns*i,shape,scale,1.0,tree->mod); */ + PMat_MGF_Gamma(mean,var,tree->mod,i*tree->mod->ns*tree->mod->ns,b_fcus->Pij_rr,b_fcus->tPij_rr); } } @@ -2719,6 +2729,11 @@ void Pull_Scaling_Factors(int site, t_edge *b, t_tree *tree) tree->fact_sum_scale[site] = sum_scale_left + sum_scale_rght; + break; + } + default : + { + assert(FALSE); break; } } diff --git a/src/main.c b/src/main.c index e1163883..94948cbe 100644 --- a/src/main.c +++ b/src/main.c @@ -252,7 +252,7 @@ int main(int argc, char **argv) Set_Update_Eigen(YES,tree->mod); Lk(NULL,tree); Set_Update_Eigen(NO,tree->mod); - + if(tree->mod->s_opt->opt_topo) { Global_Spr_Search(tree); diff --git a/src/models.c b/src/models.c index 4e7da241..b9659be5 100644 --- a/src/models.c +++ b/src/models.c @@ -245,7 +245,7 @@ phydbl Get_Lambda_F84(phydbl *pi, phydbl *kappa) /* 8000 = 20x20x20 times the operation + */ /********************************************************************/ -void PMat_Empirical(const phydbl l, const t_mod *mod, const int pos, phydbl *Pij, phydbl *tPij) +void PMat_Empirical(phydbl l, const t_mod *mod, const int pos, phydbl *Pij, phydbl *tPij) { const unsigned int ns = mod->ns; unsigned int i, j, k; @@ -253,14 +253,14 @@ void PMat_Empirical(const phydbl l, const t_mod *mod, const int pos, phydbl *Pij phydbl *expt; phydbl *uexpt; phydbl sum; - + assert(Pij); expt = mod->eigen->e_val_im; uexpt = mod->eigen->r_e_vect_im; U = mod->eigen->r_e_vect; V = mod->eigen->l_e_vect; - R = mod->eigen->e_val; /* exponential of the eigen value matrix */ + R = mod->eigen->e_val; for(k=0;kns; - expt = mod->eigen->e_val_im; - uexpt = mod->eigen->r_e_vect_im; - U = mod->eigen->r_e_vect; - V = mod->eigen->l_e_vect; - R = mod->eigen->e_val; /* exponential of the eigen value matrix */ - - if(mod->ras->n_catg == 1) shape = 1.E+4; - else shape = mod->ras->alpha->v; - - - for(i=0;ins*i+k] = .0; - - if(shape < 1.E-10) - { - PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__); - Warn_And_Exit(""); - } - - /* Formula 13.42, page 220 of Felsenstein's book ``Inferring Phylogenies'' */ - for(k=0;kns*i+j] += (uexpt[i*n+k] * V[k*n+j]); - } - if(Pij[pos+mod->ns*i+j] < SMALL_PIJ) Pij[pos+mod->ns*i+j] = SMALL_PIJ; - } -#ifdef DEBUG - phydbl sum; - sum = .0; - For (j,n) sum += Pij[pos+mod->ns*i+j]; - if((sum > 1.+.0001) || (sum < 1.-.0001)) - { - PhyML_Printf("\n"); - PhyML_Printf("\n. Q\n"); - for(i=0;ieigen->q[i*n+j]); PhyML_Printf("\n"); } - PhyML_Printf("\n. U\n"); - for(i=0;ins*i+j]); PhyML_Printf("\n"); } - PhyML_Printf("\n. sum = %f",sum); - PhyML_Fprintf(stderr,"\n. l=%f",l); - PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); - Warn_And_Exit(""); - } -#endif - } + /* PhyML_Printf("\n. Pmat len: %f",l); */ + /* for(i=0;ins; @@ -901,7 +834,7 @@ int Update_Eigen(t_mod *mod) if(mod->update_eigen == YES) { - //Update the Q-matrix first before computing the Eigen(because the Eigen is computed based on the Q-matrix) + //Update the Q-matrix first before computing the eigen decomposition (because eigen values/vectors is computed based on the Q-matrix) if(mod->io->datatype == NT) { if(mod->whichmodel == GTR) @@ -1032,55 +965,60 @@ void Switch_From_M4mod_To_Mod(t_mod *mod) ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// -void PMat_MGF_Gamma(phydbl *Pij, phydbl shape, phydbl scale, phydbl scaling_fact, t_mod *mod) +void PMat_MGF_Gamma(phydbl mu, phydbl sigsq, const t_mod *mod, const int pos, phydbl *Pij, phydbl *tPij) { - int dim; - int i,j,k; + const unsigned int ns = mod->ns; + unsigned int i,j,k; phydbl *uexpt,*imbd; + + Pij = Pij + pos; + if(tPij != NULL) tPij = tPij + pos; - dim = mod->eigen->size; uexpt = mod->eigen->r_e_vect_im; imbd = mod->eigen->e_val_im; - /* Get the eigenvalues of Q (not the exponentials) */ - for(i=0;ieigen->e_val[i]); + for(i=0;ieigen->e_val[i] * sigsq / mu, -mu*mu / sigsq); + for(i=0;ieigen->r_e_vect[i*ns+k] * imbd[k]; - /* Multiply them by the scaling factor */ - for(i=0;ieigen->r_e_vect[i*dim+k] * imbd[k]; - - for(i=0;ieigen->l_e_vect[k*dim+j]); + Pij[j] += (uexpt[i*ns+k] * mod->eigen->l_e_vect[k*ns+j]); } - if(Pij[dim*i+j] < SMALL_PIJ) Pij[dim*i+j] = SMALL_PIJ; + if(Pij[j] < SMALL_PIJ) Pij[ns*i+j] = SMALL_PIJ; } + Pij += ns; } - /* printf("\n. shape = %G scale = %G %f",shape,scale,Pij[1]); */ - /* printf("\n. Pij: %f",Pij[1]); */ + Pij -= ns*ns; + + if(tPij != NULL) + { + for(i=0;imod->s_opt->brent_it_max,tree->mod->s_opt->min_diff_lk_local,tree); Update_PMat_At_Given_Edge(b,tree); @@ -634,7 +634,6 @@ phydbl Br_Len_Opt(phydbl *l, t_edge *b, t_tree *tree) Set_Update_Eigen_Lr(NO,tree); Set_Use_Eigen_Lr(NO,tree); - /* lk_begin = Lk(b,tree); */ /* tree->n_tot_bl_opt += Generic_Brent_Lk(l, */ /* tree->mod->l_min, */ @@ -643,12 +642,13 @@ phydbl Br_Len_Opt(phydbl *l, t_edge *b, t_tree *tree) /* tree->mod->s_opt->brent_it_max, */ /* tree->mod->s_opt->quickdirty, */ /* Wrap_Lk_At_Given_Edge, */ - /* b,tree,NULL,NO); */ + /* b,tree,NULL,NO,NO); */ /* printf("\n. b->num: %4d l=%12G lnL: %12G",b->num,b->l->v,tree->c_lnL); */ /* lk_end = Lk(b,tree); /\* We can't assume that the log-lk value is up-to-date *\/ */ + lk_end = tree->c_lnL; if(lk_end < lk_begin - tree->mod->s_opt->min_diff_lk_local) diff --git a/src/utilities.c b/src/utilities.c index af2ea344..1543c4bd 100644 --- a/src/utilities.c +++ b/src/utilities.c @@ -10473,10 +10473,6 @@ phydbl *Dist_Btw_Tips(t_tree *tree) } -////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////// - - ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -10502,15 +10498,15 @@ void Best_Root_Position_IL_Model(t_tree *tree) best_edge = NULL; best_lnL = UNLIKELY; - For(i,2*tree->n_otu-3) + for(i=0;i<2*tree->n_otu-3;++i) { PhyML_Printf("\n. Positionning root node on edge %4d",tree->a_edges[i]->num); Add_Root(tree->a_edges[i],tree); tree->ignore_root = NO; Set_Both_Sides(YES,tree); Lk(NULL,tree); - /* Optimize_Br_Len_Serie(2,tree); */ + PhyML_Printf("\n HERE"); Update_Partial_Lk(tree,tree->n_root->b[1],tree->n_root); Br_Len_Opt(&(tree->n_root->b[1]->l->v),tree->n_root->b[1],tree); @@ -10549,22 +10545,23 @@ void Set_Br_Len_Var(t_edge *b, t_tree *tree) if(tree->rates == NO && tree->mod->gamma_mgf_bl == YES) { - phydbl len; if(b == NULL) { int i; - For(i,2*tree->n_otu-1) + for(i=0;i<2*tree->n_otu-1;++i) { len = MAX(0.0,tree->a_edges[i]->l->v); - tree->a_edges[i]->l_var->v = POW(len,2)*tree->mod->l_var_sigma; + /* tree->a_edges[i]->l_var->v = POW(len,2)*tree->mod->l_var_sigma; */ + tree->a_edges[i]->l_var->v = tree->mod->l_var_sigma; } } else { len = MAX(0.0,b->l->v); - b->l_var->v = POW(len,2)*tree->mod->l_var_sigma; + /* b->l_var->v = POW(len,2)*tree->mod->l_var_sigma; */ + b->l_var->v = tree->mod->l_var_sigma; } } }