Code review for IBM

stephaneguindon · Sep 12, 2024 · 188952e · 188952e
1 parent 2be8bfc
commit 188952e
Show file tree

Hide file tree

Showing 13 changed files with 454 additions and 359 deletions.
diff --git a/src/free.c b/src/free.c
@@ -1408,6 +1408,8 @@ void Free_Mmod(t_phyrex_mod *mmod)
   Free(mmod->sigsq_scale);
   Free(mmod->sigsq);
   Free(mmod->ou_mu);
+  Free(mmod->rw_root_mean);
+  Free(mmod->rw_root_var);
   Free(mmod);
 }
 

diff --git a/src/init.c b/src/init.c
@@ -3560,6 +3560,8 @@ void PHYREX_Init_Migrep_Mod(t_phyrex_mod *t, int n_dim, phydbl min_lat, phydbl m
   t->rad = 0.01*((max_lat-min_lat)+(max_lon-min_lon));
   t->prior_param_rad = 1./(0.1*((max_lat-min_lat)+(max_lon-min_lon)));
 
+  t->do_location_sampling = YES;
+
 }
 
 //////////////////////////////////////////////////////////////
@@ -3577,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.;
@@ -3641,9 +3643,12 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t)
   t->rw_prior_mean    = 1.0;
   t->rw_prior_sd      = 1.0;
   t->rw_prior_distrib = EXPONENTIAL_PRIOR;
-
-  t->rw_root_mean     = 0.0;
-  t->rw_root_var      = 1.0;
+
+  t->rw_root_mean[LOCATION]     = 0.0;
+  t->rw_root_var[LOCATION]      = 1.0;
+
+  t->rw_root_mean[VELOCITY]     = 0.0;
+  t->rw_root_var[VELOCITY]      = 1.0;
 
   t->integrateAncestralLocations = YES;
 

diff --git a/src/io.c b/src/io.c
@@ -6998,8 +6998,8 @@ void PHYREX_Print_MCMC_Stats(t_tree *tree)
           }
           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");
+            PhyML_Fprintf(fp_stats, "%d_%s_VelocLon%s\t", i,"anc",(tree->a_nodes[i] == tree->n_root->v[1] || tree->a_nodes[i] == tree->n_root->v[2]) ? "*" : "");
+            PhyML_Fprintf(fp_stats, "%d_%s_VelocLat%s\t", i,"anc",(tree->a_nodes[i] == tree->n_root->v[1] || tree->a_nodes[i] == tree->n_root->v[2]) ? "*" : "");
           }
 
           /* for(int i=0;i<tree->n_otu;++i) */
@@ -7074,7 +7074,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_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 && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_Node_Locations_Joint(tree);
 
       PhyML_Fprintf(fp_stats,"%g\t",tree->n_root->ldsk->coord->lonlat[0]);
       PhyML_Fprintf(fp_stats,"%g\t",tree->n_root->ldsk->coord->lonlat[1]);
@@ -7253,7 +7253,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_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 && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_Node_Locations_Joint(tree); */
       Record_Br_Len(tree);
       TIMES_Time_To_Bl(tree);
       tree->bl_ndigits = 3;
@@ -7311,7 +7311,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_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 && tree->mmod->integrateAncestralLocations == YES) VELOC_Sample_Node_Locations_Joint(tree);
 
       if(VELOC_Is_Integrated_Velocity(tree->mmod) == YES) sprintf(s,"%13f",VELOC_Mean_Speed(tree));
       else if(RRW_Is_Rw(tree->mmod) == YES) sprintf(s,"%13f",PHYREX_Realized_Dispersal_Dist(tree->mmod->dist_type,tree)/PHYREX_Time_Tree_Length(tree));      

diff --git a/src/make.c b/src/make.c
@@ -1741,6 +1741,10 @@ t_phyrex_mod *PHYREX_Make_Migrep_Model(int n_otu,int dim)
   t->sigsq = (phydbl *)mCalloc(dim,sizeof(phydbl));
   t->ou_mu = (phydbl *)mCalloc(dim,sizeof(phydbl));
   t->sigsq_scale = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
+
+  t->rw_root_mean = (phydbl *)mCalloc(2,sizeof(phydbl));
+  t->rw_root_var = (phydbl *)mCalloc(2,sizeof(phydbl));
+
   return(t);
 }