Skip to content

Commit

Permalink
Cleaning generalized mixture code in PhyML
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Jul 18, 2024
1 parent ccde895 commit 1ece9fa
Show file tree
Hide file tree
Showing 14 changed files with 241 additions and 128 deletions.
10 changes: 3 additions & 7 deletions src/cl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,7 +278,7 @@ int Read_Command_Line(option *io, int argc, char **argv)
#endif
case 74:
{
io->mod->l_var_sigma = String_To_Dbl(optarg);
io->mod->l_var_sigma->v = String_To_Dbl(optarg);
break;
}
case 73:
Expand Down Expand Up @@ -396,10 +396,6 @@ 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 [email protected] 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;
Expand Down Expand Up @@ -1823,8 +1819,8 @@ int Read_Command_Line(option *io, int argc, char **argv)
PhyML_Printf("\n\n. Command line: ");
for(i=0;i<argc;i++) PhyML_Printf("%s ",argv[i]);
}

return 1;
return 1;
}

//////////////////////////////////////////////////////////////
Expand Down
1 change: 1 addition & 0 deletions src/free.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -591,6 +591,7 @@ void Free_Model_Basic(t_mod *mixt_mod)
Free_Scalar_Dbl(mixt_mod->lambda);
Free_Scalar_Dbl(mixt_mod->br_len_mult);
Free_Scalar_Dbl(mixt_mod->br_len_mult_unscaled);
Free_Scalar_Dbl(mixt_mod->l_var_sigma);

Free_Rmat_Weights(mixt_mod);
Free_Efrq_Weights(mixt_mod);
Expand Down
4 changes: 2 additions & 2 deletions src/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -683,7 +683,7 @@ void Set_Defaults_Model(t_mod *mod)

mod->kappa->v = 4.0;
mod->lambda->v = 1.0;
mod->l_var_sigma = 1.E-1;
mod->l_var_sigma->v = 1.E-1;
mod->e_frq_weight->v = 1.0;
mod->r_mat_weight->v = 1.0;

Expand Down Expand Up @@ -797,7 +797,7 @@ void Set_Defaults_Optimiz(t_opt *s_opt)
s_opt->br_len_in_spr = 10;
s_opt->opt_free_mixt_rates = YES;
s_opt->constrained_br_len = NO;
s_opt->opt_gamma_br_len = NO;
s_opt->opt_gamma_br_len = YES;
s_opt->first_opt_free_mixt_rates = YES;

s_opt->wim_n_rgrft = -1;
Expand Down
34 changes: 26 additions & 8 deletions src/io.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2376,9 +2376,9 @@ void Print_Fp_Out(FILE *fp_out, time_t t_beg, time_t t_end, t_tree *tree, option
PhyML_Fprintf(fp_out,"\n. Substitution model: \t\t%s",tree->mod->modelname->s);
}

PhyML_Fprintf(fp_out,"\n. Integrated length (IL) model: \t%s",(tree->io->mod->gamma_mgf_bl == YES)?"yes":"no");
PhyML_Fprintf(fp_out,"\n. Integrated length (IL) model: \t%s",(tree->mod->gamma_mgf_bl == YES)?"yes":"no");

if(tree->io->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp_out,"\n. Variance of IL model: \t\t%.5f",tree->mod->l_var_sigma);
if(tree->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp_out,"\n. Variance of IL model: \t\t%f",tree->mod->l_var_sigma->v);


PhyML_Fprintf(fp_out,"\n. Number of taxa: \t\t\t%d",tree->n_otu);/*added FLT*/
Expand Down Expand Up @@ -2421,8 +2421,6 @@ void Print_Fp_Out(FILE *fp_out, time_t t_beg, time_t t_end, t_tree *tree, option

if(tree->mod->ras->invar) PhyML_Fprintf(fp_out,"\n. Proportion of invariant: \t\t%.3f",tree->mod->ras->pinvar->v);

if(tree->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp_out,"\n. Variance of branch lengths: \t\t%f",tree->mod->l_var_sigma);


/*was before Discrete gamma model ; moved here FLT*/
if((tree->mod->whichmodel == K80) ||
Expand Down Expand Up @@ -3130,9 +3128,9 @@ void Print_Banner(FILE *fp)
PhyML_Fprintf(fp," A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood \n");
PhyML_Fprintf(fp," Stephane Guindon & Olivier Gascuel \n");
PhyML_Fprintf(fp," \n");
PhyML_Fprintf(fp," https://github.com/stephaneguindon/phyml/tree/master \n");
PhyML_Fprintf(fp," https://github.com/stephaneguindon/phyml/ \n");
PhyML_Fprintf(fp," \n");
PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier \n");
PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier \n");
PhyML_Fprintf(fp," \n");
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
}
Expand All @@ -3146,7 +3144,7 @@ void Print_Banner_Small(FILE *fp)
PhyML_Fprintf(fp,"\n");
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
PhyML_Fprintf(fp," --- PhyML %s --- \n",VERSION);
PhyML_Fprintf(fp," https://github.com/stephaneguindon/phyml/tree/master \n");
PhyML_Fprintf(fp," https://github.com/stephaneguindon/phyml/ \n");
PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier \n");
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
}
Expand Down Expand Up @@ -4289,6 +4287,10 @@ void Print_Data_Structure(int final, FILE *fp, t_tree *mixt_tree)

PhyML_Fprintf(fp,"\n Equ. freq. weight:\t\t%20f",tree->mod->e_frq_weight->v / e_frq_weight_sum);

PhyML_Fprintf(fp,"\n Integrated length (IL) model:%20s",(tree->mod->gamma_mgf_bl == YES)?"yes":"no");

if(tree->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp,"\n Variance of IL model:\t%20.2g",tree->mod->l_var_sigma->v);

class++;

tree = tree->next;
Expand Down Expand Up @@ -5073,6 +5075,22 @@ void Make_Ratematrix_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
{
mod->s_opt->opt_rmat_weight = NO;
}


char *il_model = NULL;
il_model = XML_Get_Attribute_Value(instance,"integrated.lens");

if(il_model)
{
if(!strcmp(il_model,"yes") || !strcmp(il_model,"true"))
{
mod->gamma_mgf_bl = YES;
}
else
{
mod->gamma_mgf_bl = NO;
}
}
}

//////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -5214,7 +5232,7 @@ void Make_Topology_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
// Estimate tree topology
char *optimise = NULL;

optimise = XML_Get_Attribute_Value(instance,"optimise.tree");
optimise = XML_Get_Attribute_Value(instance,"optimise.topology");

if(optimise)
{
Expand Down
19 changes: 7 additions & 12 deletions src/lk.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -688,11 +688,12 @@ phydbl dLk(phydbl *l, t_edge *b, t_tree *tree)
for(catg=0;catg<ncatg;catg++)
{
rr = tree->mod->ras->gamma_rr->v[catg];
rr *= tree->mod->br_len_mult->v;
if(tree->mixt_tree) rr = tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
rr *= tree->mod->br_len_mult->v;

len = (*l) * rr;
/* var = tree->mod->l_var_sigma * rr*rr; */
var = (*l) * tree->mod->l_var_sigma * rr*rr;
var = (*l) * tree->mod->l_var_sigma->v * rr*rr;

if(isinf(len) || isnan(len))
{
Expand All @@ -711,11 +712,11 @@ phydbl dLk(phydbl *l, t_edge *b, t_tree *tree)
ev = tree->mod->eigen->e_val[state];
expevlen = exp(ev*len);

if(tree->io->mod->gamma_mgf_bl == YES)
if(tree->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);
expl[catg*2*ns + 2*state + 1] *= -(ev * rr/(1.-ev*var/len) + 2.*len*rr*LOG(1.-ev*var/len)/var);
}
else
{
Expand Down Expand Up @@ -2202,7 +2203,7 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree)

if(tree->mixt_tree != NULL) assert(tree->mod->ras->n_catg == 1);

if(tree->io->mod->gamma_mgf_bl == YES) Set_Br_Len_Var(b_fcus,tree);
if(tree->mod->gamma_mgf_bl == YES) Set_Br_Len_Var(b_fcus,tree);

l_min = tree->mod->l_min;
l_max = tree->mod->l_max;
Expand Down Expand Up @@ -2236,13 +2237,8 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree)
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 = tree->mod->l_var_sigma * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2); */
var = MAX(0.0,b_fcus->l->v) * tree->mod->l_var_sigma * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2);
var = MAX(0.0,b_fcus->l->v) * tree->mod->l_var_sigma->v * 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; */

if(var > tree->mod->l_var_max) var = tree->mod->l_var_max;
if(var < tree->mod->l_var_min) var = tree->mod->l_var_min;
}

//Update the transition prob. matrix
Expand All @@ -2260,7 +2256,6 @@ void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree)
Warn_And_Exit(TODO_BEAGLE);
#endif

/* 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);
}
}
Expand Down
8 changes: 6 additions & 2 deletions src/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,7 +252,10 @@ int main(int argc, char **argv)
Set_Update_Eigen(YES,tree->mod);
Lk(NULL,tree);
Set_Update_Eigen(NO,tree->mod);



PhyML_Printf("\n. Init log-likelihood: %f",tree->c_lnL);

if(tree->mod->s_opt->opt_topo)
{
Global_Spr_Search(tree);
Expand Down Expand Up @@ -280,7 +283,8 @@ int main(int argc, char **argv)
Check_Br_Lens(tree);
Br_Len_Involving_Invar(tree);
Rescale_Br_Len_Multiplier_Tree(tree);



#elif defined EVOLVE
Evolve(tree->data,tree->mod,tree);
Exit("\n. Exiting 'evolve'\n");
Expand Down
3 changes: 3 additions & 0 deletions src/make.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -934,6 +934,9 @@ t_mod *Make_Model_Basic(void)
mod->aa_rate_mat_file = Make_String(T_MAX_FILE);
Init_String(mod->aa_rate_mat_file);

mod->l_var_sigma = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
Init_Scalar_Dbl(mod->l_var_sigma);

return mod;
}

Expand Down
30 changes: 21 additions & 9 deletions src/mixt.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2844,7 +2844,7 @@ phydbl MIXT_dLk(phydbl *l, t_edge *mixt_b, t_tree *mixt_tree)
phydbl log_site_lk,inv_site_lk;
int num_prec_issue;
phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas;
phydbl len;
phydbl len,var;
phydbl *expl,*dot_prod;
phydbl rr;
phydbl ev,expevlen;
Expand All @@ -2860,10 +2860,8 @@ phydbl MIXT_dLk(phydbl *l, t_edge *mixt_b, t_tree *mixt_tree)
cpy_mixt_b = mixt_b;
len = -1.;

if(mixt_tree->update_eigen_lr == YES)
{
MIXT_Update_Eigen_Lr(mixt_b,mixt_tree);
}
if(mixt_tree->update_eigen_lr == YES) MIXT_Update_Eigen_Lr(mixt_b,mixt_tree);


/*! Make sure that l is one of the lengths of mixt_b */
b = mixt_b;
Expand Down Expand Up @@ -2951,17 +2949,31 @@ phydbl MIXT_dLk(phydbl *l, t_edge *mixt_b, t_tree *mixt_tree)
{
len = 0.0;
}


if(length_found == YES)
var = (*l) * tree->mod->l_var_sigma->v * rr*rr;
else
var = b->l->v * tree->mod->l_var_sigma->v * rr*rr;


expl = mixt_tree->expl;

for(state=0;state<tree->mod->ns;++state)
{
ev = tree->mod->eigen->e_val[state];
expevlen = exp(ev*len);

expl[class*2*ns + 2*state] = expevlen;
expl[class*2*ns + 2*state + 1] = expevlen*ev*rr;

if(tree->mod->gamma_mgf_bl == YES)
{
expl[class*2*ns + 2*state] = POW(1. - ev*var/len,-len*len/var);
expl[class*2*ns + 2*state + 1] = expl[class*2*ns + 2*state];
expl[class*2*ns + 2*state + 1] *= -(ev * rr/(1.-ev*var/len) + 2.*len*rr*LOG(1.-ev*var/len)/var);
}
else
{
expl[class*2*ns + 2*state] = expevlen;
expl[class*2*ns + 2*state + 1] = expevlen*ev*rr;
}
}

class++;
Expand Down
Loading

0 comments on commit 1ece9fa

Please sign in to comment.