Skip to content

Commit

Permalink
Put scaling factors back on
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Jul 10, 2024
1 parent fb891e2 commit efd0397
Show file tree
Hide file tree
Showing 8 changed files with 94 additions and 144 deletions.
8 changes: 4 additions & 4 deletions src/cl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 [email protected] 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 [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
4 changes: 2 additions & 2 deletions src/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down
35 changes: 25 additions & 10 deletions src/lk.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;state<ns;++state) expl[catg*ns+state] = (phydbl)pow(tree->mod->eigen->e_val[state],len);
for(state=0;state<ns;++state) expl[catg*ns+state] = exp(tree->mod->eigen->e_val[state]*len);
}
}

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

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

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

Expand Down Expand Up @@ -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;
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
162 changes: 50 additions & 112 deletions src/models.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,22 +245,22 @@ 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;
const phydbl *U,*V,*R;
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;k<ns;++k) expt[k] = exp(R[k]*l);

Expand All @@ -281,6 +281,7 @@ void PMat_Empirical(const phydbl l, const t_mod *mod, const int pos, phydbl *Pij
}
if(Pij[j] < SMALL_PIJ) Pij[j] = SMALL_PIJ;
}
/* Below is required ? */
sum = 0.0;
for(j=0;j<ns;++j) sum += Pij[j];
for(j=0;j<ns;++j) Pij[j] /= sum;
Expand All @@ -300,90 +301,22 @@ void PMat_Empirical(const phydbl l, const t_mod *mod, const int pos, phydbl *Pij
}
}
}
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void PMat_Gamma(phydbl l, t_mod *mod, int pos, phydbl *Pij)
{
int n;
int i, j, k;
phydbl *U,*V,*R;
phydbl *expt;
phydbl *uexpt;
phydbl shape;


n = mod->ns;
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;i<n;i++) for(k=0;k<n;k++) Pij[pos+mod->ns*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;k<n;k++) expt[k] = POW(shape/(shape-log(R[k])*l),shape);

/* multiply Vr*expt*Vi into Pij */
for(i=0;i<n;i++) for(k=0;k<n;k++) uexpt[i*n+k] = U[i*n+k] * expt[k];

For (i,n)
{
For (j,n)
{
for(k=0;k<n;k++)
{
Pij[pos+mod->ns*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;i<n;i++) { for(j=0;j<n;j++) PhyML_Printf("%7.3f ",mod->eigen->q[i*n+j]); PhyML_Printf("\n"); }
PhyML_Printf("\n. U\n");
for(i=0;i<n;i++) { for(j=0;j<n;j++) PhyML_Printf("%7.3f ",U[i*n+j]); PhyML_Printf("\n"); }
PhyML_Printf("\n");
PhyML_Printf("\n. V\n");
for(i=0;i<n;i++) { for(j=0;j<n;j++) PhyML_Printf("%7.3f ",V[i*n+j]); PhyML_Printf("\n"); }
PhyML_Printf("\n");
PhyML_Printf("\n. Eigen\n");
for(i=0;i<n;i++) PhyML_Printf("%E ",expt[i]);
PhyML_Printf("\n");
PhyML_Printf("\n. Pij\n");
for(i=0;i<n;i++) { For (j,n) PhyML_Printf("%f ",Pij[pos+mod->ns*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;i<ns;i++) */
/* { */
/* PhyML_Printf("\n"); */
/* for(j=0;j<ns;j++) */
/* { */
/* PhyML_Printf("%12f ",Pij[i*ns+j]); */
/* } */
/* } */
/* Exit("\n"); */
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void PMat_Zero_Br_Len(t_mod *mod, int pos, phydbl *Pij)
{
int n = mod->ns;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;i<dim;i++) imbd[i] = log(mod->eigen->e_val[i]);
for(i=0;i<ns;i++) imbd[i] = POW(1. - mod->eigen->e_val[i] * sigsq / mu, -mu*mu / sigsq);
for(i=0;i<ns;i++) for(k=0;k<ns;k++) uexpt[i*ns+k] = mod->eigen->r_e_vect[i*ns+k] * imbd[k];

/* Multiply them by the scaling factor */
for(i=0;i<dim;i++) imbd[i] *= scaling_fact;
for(i=0;i<ns;i++) for(k=0;k<ns;k++) Pij[ns*i+k] = .0;

for(i=0;i<dim;i++) imbd[i] *= -scale;
for(i=0;i<dim;i++) imbd[i] += 1.0;
for(i=0;i<dim;i++) imbd[i] = POW(imbd[i],-shape);

for(i=0;i<dim;i++) for(k=0;k<dim;k++) uexpt[i*dim+k] = mod->eigen->r_e_vect[i*dim+k] * imbd[k];

for(i=0;i<dim;i++) for(k=0;k<dim;k++) Pij[dim*i+k] = .0;

for(i=0;i<dim;i++)
for(i=0;i<ns;i++)
{
for(j=0;j<dim;j++)
for(j=0;j<ns;j++)
{
for(k=0;k<dim;k++)
for(k=0;k<ns;k++)
{
Pij[dim*i+j] += (uexpt[i*dim+k] * mod->eigen->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;i<ns;++i)
{
for(j=0;j<ns;++j)
{
tPij[ns*i+j] = Pij[ns*j+i];
}
}
}

/* printf("\n. Pmat"); */
/* for(i=0;i<dim;i++) */
/* PhyML_Printf("\n. MGF Pmat"); */
/* for(i=0;i<ns;i++) */
/* { */
/* printf("\n"); */
/* for(j=0;j<dim;j++) */
/* { */
/* printf("%12f ",Pij[i*dim+j]); */
/* } */
/* PhyML_Printf("\n"); */
/* for(j=0;j<ns;j++) */
/* { */
/* PhyML_Printf("%12f ",Pij[i*ns+j]); */
/* } */
/* } */
/* Exit("\n"); */

}

//////////////////////////////////////////////////////////////
Expand Down
4 changes: 2 additions & 2 deletions src/models.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ the GNU public licence. See http://www.opensource.org for details.
void PMat(phydbl l, t_mod *mod, int pos, phydbl *Pij, phydbl *tPij);
void PMat_K80(phydbl l,phydbl kappa, int pos, phydbl *Pij);
void PMat_TN93(phydbl l, t_mod *mod, int pos, phydbl *Pij);
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);
void PMat_Zero_Br_Len(t_mod *mod, int pos, phydbl *Pij);
void PMat_Gamma(phydbl l, t_mod *mod, int pos, phydbl *Pij);
int GetDaa (phydbl *daa, phydbl *pi, char *file_name);
Expand All @@ -42,7 +42,7 @@ phydbl Get_Lambda_F84(phydbl *pi, phydbl *kappa);
int Update_Eigen(t_mod *mod);
int Update_RAS(t_mod *mod);
int Update_Efrq(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 Update_Boundaries(t_mod *mod);
void Update_Qmat_TN93(phydbl kappa, phydbl lambda, phydbl *pi, phydbl *qmat);

Expand Down
Loading

0 comments on commit efd0397

Please sign in to comment.