Skip to content

Commit

Permalink
Joint sampling of ancestral locations in RRW models
Browse files Browse the repository at this point in the history
  • Loading branch information
@stephaneguindon
stephaneguindon committed Sep 2, 2024
1 parent e43719c commit b984202
Show file tree
Hide file tree
Showing 10 changed files with 480 additions and 426 deletions.
121 changes: 9 additions & 112 deletions src/ibm.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,57 +113,10 @@ void IBM_Integrated_Location_Down(phydbl son_a, phydbl son_b, phydbl son_mu_down
phydbl son_logrem, phydbl bro_logrem,
phydbl *mean, phydbl *var, phydbl *logrem)
{
phydbl m,v,logr;
int err;

err = 0;
v = 0.0;

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(err == YES)
{
PhyML_Printf("\n. son_mu_down: %f son_b: %f son_a: %f bro_mu_down: %f bro_b: %f bro_a: %f son_var_down: %f son_var: %f bro_var_down: %f bro_var: %f",
son_mu_down,
son_b,
son_a,
bro_mu_down,
bro_b,
bro_a,
son_var_down,
son_var,
bro_var_down,
bro_var);
assert(false);
}

if((son_var + son_var_down > 1.E-7) && (bro_var + bro_var_down > 1.E-7)) // Standard case
{
v = pow(son_a,2)/(son_var_down + son_var) + pow(bro_a,2)/(bro_var_down + bro_var);
v = 1./v;

m = (son_a*(son_mu_down-son_b)/(son_var_down + son_var) + bro_a*(bro_mu_down-bro_b)/(bro_var_down + bro_var)) * v;
}
else if(son_var + son_var_down > 1.E-7) // Null variance along d - v2
{
m = (bro_mu_down-bro_b)/bro_a;
}
else if(bro_var + bro_var_down > 1.E-7) // Null variance along d - v1
{
m = (son_mu_down-son_b)/son_a;
}
else
{
m = (son_mu_down-son_b)/son_a;
}

*mean = m;
*var = v;
*logrem = logr;

/* assert(v>0.0); */
RW_Integrated_Lk_Down(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, bro_logrem,
mean, var, logrem);
}

//////////////////////////////////////////////////////////////
Expand All @@ -172,69 +125,13 @@ void IBM_Integrated_Location_Down(phydbl son_a, phydbl son_b, phydbl son_mu_down
void IBM_Integrated_Location_Up(phydbl dad_mu_up, phydbl dad_var_up, phydbl dad_logrem_up,
phydbl son_a, phydbl son_b, phydbl son_var,
phydbl bro_a, phydbl bro_b, phydbl bro_mu_down, phydbl bro_var_down, phydbl bro_var, phydbl bro_logrem_down,


phydbl *mean, phydbl *var, phydbl *logrem)
{
phydbl m,v,logr;
int err;
{
RW_Integrated_Lk_Up(dad_mu_up, dad_var_up, dad_logrem_up,
son_a, son_b, son_var,
bro_a, bro_b, bro_mu_down, bro_var_down, bro_var, bro_logrem_down,
mean, var, logrem);

v = logr = 0.0;
err = 0;

logr = dad_logrem_up + bro_logrem_down;
logr -= log(fabs(bro_a));
logr += Log_Dnorm((bro_mu_down-bro_b)/bro_a,dad_mu_up,sqrt((bro_var_down+bro_var)/(bro_a*bro_a)+dad_var_up),&err);

if(err == YES)
{
PhyML_Printf("\n. bro_mu_down: %f bro_b: %f bro_a: %f dad_mu_up: %f bro_var_down: %f bro_var: %f dad_var_up: %f",
bro_mu_down,
bro_b,
bro_a,
dad_mu_up,
bro_var_down,
bro_var,
dad_var_up);
assert(false);
}


if((bro_var_down + bro_var) > 1.E-7 && dad_var_up > 1.E-7) // Standard case
{

v += dad_var_up * (bro_var_down + bro_var);
v /= bro_a * bro_a * dad_var_up + bro_var_down + bro_var;

m = v * (bro_a * (bro_mu_down - bro_b) / (bro_var_down + bro_var) + dad_mu_up / dad_var_up);
m = son_a * m + son_b;

v = son_a * son_a * v + son_var;
}
else if(dad_var_up > 1.E-7) // Null variance along bro
{
v = son_var;
m = son_a * bro_mu_down / bro_a + son_b;
}
else if((bro_var_down + bro_var) > 1.E-7) // Null variance along dad
{
v = son_var;
m = son_a * dad_mu_up + son_b;
}
else
{
v = son_var;
m = son_a * bro_mu_down / bro_a + son_b;
/* m = son_b; */
}


*mean = m;
*var = v;
*logrem = logr;

/* !!!!!!!!!!!!!!!!!!!! */
/* assert(v>0.0); */
}

//////////////////////////////////////////////////////////////
Expand Down
78 changes: 36 additions & 42 deletions src/io.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6249,7 +6249,6 @@ void Collect_Edge_Support_Values(t_tree *tree)
/*////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////*/


#if defined (PHYREX)
void PHYREX_Output_Tree_Structure(FILE *fp, t_tree *tree)
{
Expand Down Expand Up @@ -6989,36 +6988,20 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
/* tree->a_nodes[i] == tree->n_root->v[2]) ? */
/* '*':' '); */

if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES)
{
PhyML_Fprintf(fp_stats,"root_VelocLon\t");
PhyML_Fprintf(fp_stats,"root_VelocLat\t");

for(int i=0;i<tree->n_otu;++i)
{
PhyML_Fprintf(fp_stats,"%s_VelocLon\t",tree->a_nodes[i]->name);
PhyML_Fprintf(fp_stats,"%s_VelocLat\t",tree->a_nodes[i]->name);
}
for(int i=tree->n_otu;i<2*tree->n_otu-1;++i)
{
PhyML_Fprintf(fp_stats,"%s_VelocLon\t","anc");
PhyML_Fprintf(fp_stats,"%s_VelocLat\t","anc");
}

/* for(int i=0;i<tree->n_otu-1;++i) */
/* { */
/* PhyML_Fprintf(fp_stats,"%d_IntVelocLon\t",tree->a_nodes[tree->n_otu+i]->num); */
/* PhyML_Fprintf(fp_stats,"%d_IntVelocLat\t",tree->a_nodes[tree->n_otu+i]->num); */
/* } */
PhyML_Fprintf(fp_stats, "root_VelocLon\t");
PhyML_Fprintf(fp_stats, "root_VelocLat\t");

for (int i = 0; i < tree->n_otu; ++i)
{
PhyML_Fprintf(fp_stats, "%s_VelocLon\t", tree->a_nodes[i]->name);
PhyML_Fprintf(fp_stats, "%s_VelocLat\t", tree->a_nodes[i]->name);
}
for (int i = tree->n_otu; i < 2 * tree->n_otu - 1; ++i)
{
PhyML_Fprintf(fp_stats, "%s_VelocLon\t", "anc");
PhyML_Fprintf(fp_stats, "%s_VelocLat\t", "anc");
}

/* for(int i=0;i<tree->n_otu;++i) */
/* { */
/* PhyML_Fprintf(fp_stats,"%s_Speed\t",tree->a_nodes[i]->name); */
/* } */
}


/* for(int i=0;i<tree->n_otu;++i) */
/* { */
/* PhyML_Fprintf(fp_stats,"%s_DispLon\t",tree->a_nodes[i]->name); */
Expand Down Expand Up @@ -7090,7 +7073,7 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)

PhyML_Fprintf(fp_stats,"%.2f\t",tree->n_root->ldsk->disk->time);

if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations(tree);
if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations_Joint(tree);
if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_All_Node_Locations(tree);

PhyML_Fprintf(fp_stats,"%g\t",tree->n_root->ldsk->coord->lonlat[0]);
Expand Down Expand Up @@ -7170,16 +7153,7 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
/* for(int i=0;i<2*tree->n_otu-2;++i) PhyML_Fprintf(fp_stats,"%f\t",tree->rates->br_r[tree->a_nodes[i]->num]); */

if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES)
{
/* PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Degrees_To_Km(tree->n_root->ldsk->veloc->deriv[0],tree)); */
/* PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Degrees_To_Km(tree->n_root->ldsk->veloc->deriv[1],tree)); */

/* for(int i=0;i<tree->n_otu-1;++i) */
/* { */
/* PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Degrees_To_Km(tree->a_nodes[i+tree->n_otu]->ldsk->veloc->deriv[0])); */
/* PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Degrees_To_Km(tree->a_nodes[i+tree->n_otu]->ldsk->veloc->deriv[1])); */
/* } */

{
PhyML_Fprintf(fp_stats,"%g\t",tree->n_root->ldsk->veloc->deriv[0]);
PhyML_Fprintf(fp_stats,"%g\t",tree->n_root->ldsk->veloc->deriv[1]);

Expand All @@ -7190,6 +7164,21 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
PhyML_Fprintf(fp_stats,"%g\t",tree->a_nodes[i]->ldsk->veloc->deriv[1]);
}
}
else if(RRW_Is_Rw(tree->mmod) == YES)
{
RRW_Tip_Velocities(tree);

PhyML_Fprintf(fp_stats,"%g\t",-1.);
PhyML_Fprintf(fp_stats,"%g\t",-1);

/* for(int i=0;i<tree->n_otu;++i) */
for(int i=0;i<2*tree->n_otu-1;++i)
{
PhyML_Fprintf(fp_stats,"%g\t",tree->a_nodes[i]->ldsk->veloc->deriv[0]);
PhyML_Fprintf(fp_stats,"%g\t",tree->a_nodes[i]->ldsk->veloc->deriv[1]);
}

}


/* for(int i=0;i<tree->n_otu;++i) */
Expand Down Expand Up @@ -7263,7 +7252,7 @@ void PHYREX_Print_MCMC_Tree(t_tree *tree)

if(!(tree->mcmc->run%tree->mcmc->sample_interval) && tree->mcmc->sample_interval > 0)
{
/* if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations(tree); */
/* if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations_Joint(tree); */
/* if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_All_Node_Locations(tree); */
Record_Br_Len(tree);
TIMES_Time_To_Bl(tree);
Expand All @@ -7272,6 +7261,11 @@ void PHYREX_Print_MCMC_Tree(t_tree *tree)
tree->write_tax_names = NO;
Label_Nodes_With_Locations(tree);
if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES) Label_Nodes_With_Velocities(tree);
else if(RRW_Is_Rw(tree->mmod) == YES)
{
RRW_Tip_Velocities(tree);
Label_Nodes_With_Velocities(tree);
}
Label_Edges(tree);
char *s = Write_Tree(tree);
Free_All_Node_Labels(tree);
Expand Down Expand Up @@ -7316,7 +7310,7 @@ void PHYREX_Print_MCMC_Summary(t_tree *tree)

if(!(tree->mcmc->run%tree->mcmc->print_every) && tree->mcmc->print_every > 0)
{
if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations(tree);
if(RRW_Is_Rw(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) RRW_Sample_Node_Locations_Joint(tree);
if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_All_Node_Locations(tree);

if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES) sprintf(s,"%13f",VELOC_Mean_Speed(tree));
Expand Down
82 changes: 8 additions & 74 deletions src/iou.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,40 +178,10 @@ void IOU_Integrated_Location_Down(phydbl son_a, phydbl son_b, phydbl son_mu_down
phydbl son_logrem, phydbl bro_logrem,
phydbl *mean, phydbl *var, phydbl *logrem)
{
phydbl m,v,logr;
int err;

err = 0;
v = 0.0;

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 > 1.E-7) && (bro_var + bro_var_down > 1.E-7)) // Standard case
{
v = pow(son_a,2)/(son_var_down + son_var) + pow(bro_a,2)/(bro_var_down + bro_var);
v = 1/v;

m = (son_a*(son_mu_down-son_b)/(son_var_down + son_var) + bro_a*(bro_mu_down-bro_b)/(bro_var_down + bro_var)) * v;
}
else if(son_var + son_var_down > 1.E-7) // Null variance along d - v2
{
m = (bro_mu_down-bro_b)/bro_a;
}
else if(bro_var + bro_var_down > 1.E-7) // Null variance along d - v1
{
m = (son_mu_down-son_b)/son_a;
}
else
{
m = (son_mu_down-son_b)/son_a;
}

*mean = m;
*var = v;
*logrem = logr;

RW_Integrated_Lk_Down(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, bro_logrem,
mean, var, logrem);
}

//////////////////////////////////////////////////////////////
Expand All @@ -222,46 +192,10 @@ void IOU_Integrated_Location_Up(phydbl dad_mu_up, phydbl dad_var_up, phydbl dad_
phydbl bro_a, phydbl bro_b, phydbl bro_mu_down, phydbl bro_var_down, phydbl bro_var, phydbl bro_logrem_down,
phydbl *mean, phydbl *var, phydbl *logrem)
{
phydbl m,v,logr;
int err;

v = logr = 0.0;
err = 0;

logr = dad_logrem_up + bro_logrem_down;
logr -= log(fabs(bro_a));
logr += Log_Dnorm((bro_mu_down-bro_b)/bro_a,dad_mu_up,sqrt((bro_var_down+bro_var)/(bro_a*bro_a)+dad_var_up),&err);

if((bro_var_down + bro_var) > 1.E-7 && dad_var_up > 1.E-7) // Standard case
{

v += dad_var_up * (bro_var_down + bro_var);
v /= bro_a * bro_a * dad_var_up + bro_var_down + bro_var;

m = v * (bro_a * (bro_mu_down - bro_b) / (bro_var_down + bro_var) + dad_mu_up / dad_var_up);
m = son_a * m + son_b;

v = son_a * son_a * v + son_var;
}
else if(dad_var_up > 1.E-7) // Null variance along bro
{
v = son_var;
m = son_a * (bro_mu_down - bro_b) / bro_a + son_b;
}
else if((bro_var_down + bro_var) > 1.E-7) // Null variance along dad
{
v = son_var;
m = son_a * dad_mu_up + son_b;
}
else
{
v = son_var;
m = son_b;
}

*mean = m;
*var = v;
*logrem = logr;
RW_Integrated_Lk_Up(dad_mu_up, dad_var_up, dad_logrem_up,
son_a, son_b, son_var,
bro_a, bro_b, bro_mu_down, bro_var_down, bro_var, bro_logrem_down,
mean, var, logrem);
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
3 changes: 1 addition & 2 deletions src/mcmc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14569,8 +14569,7 @@ void MCMC_Crossvalidate_Locations(t_tree *tree)
(void)signal(SIGINT,MCMC_Terminate);
}
while(mcmc->run < mcmc->chain_len);



PHYREX_Restore_Geo_Coord(tree->a_nodes[i]->ldsk->coord);
}
}
Expand Down
Loading

0 comments on commit b984202

Please # to comment.