1 /*
2
3 PhyML: a program that computes maximum likelihood phylogenies from
4 DNA or AA homologous sequences.
5
6 Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8 All parts of the source except where indicated are distributed under
9 the GNU public licence. See http://www.opensource.org for details.
10
11 */
12
13 /* Routines that implement Etheridge and Barton's model of continuous-space
14 coalescent.
15 */
16
17 #include "phyrex.h"
18
19 //////////////////////////////////////////////////////////////
20 //////////////////////////////////////////////////////////////
21
PHYREX_Main(int argc,char * argv[])22 int PHYREX_Main(int argc, char *argv[])
23 {
24 #if (defined PHYREXSIM)
25 PHYREX_Main_Simulate(argc,argv);
26 #elif (defined PHYREX)
27 option *io;
28 io = Get_Input(argc,argv);
29 Free(io);
30 #endif
31 return(0);
32 }
33
34 //////////////////////////////////////////////////////////////
35 //////////////////////////////////////////////////////////////
36
37 #if (defined PHYREX)
PHYREX_XML(char * xml_filename)38 void PHYREX_XML(char *xml_filename)
39 {
40 FILE *fp_xml_in;
41 xml_node *xnd,*xroot;
42 t_tree *mixt_tree,*tree;
43 phydbl *res;
44 int seed;
45 char *dum_string;
46
47 mixt_tree = XML_Process_Base(xml_filename);
48 assert(mixt_tree);
49
50 mixt_tree->rates = RATES_Make_Rate_Struct(mixt_tree->n_otu);
51 RATES_Init_Rate_Struct(mixt_tree->rates,NULL,mixt_tree->n_otu);
52
53 mixt_tree->times = TIMES_Make_Time_Struct(mixt_tree->n_otu);
54 TIMES_Init_Time_Struct(mixt_tree->times,NULL,mixt_tree->n_otu);
55
56 mixt_tree->mmod = PHYREX_Make_Migrep_Model(mixt_tree->n_otu,2);
57 mixt_tree->mmod->n_dim = 2;
58 PHYREX_Set_Default_Migrep_Mod(mixt_tree->n_otu,mixt_tree->mmod);
59
60 tree = mixt_tree;
61 do
62 {
63 // All rate stuctures point to the same object
64 tree->rates = mixt_tree->rates;
65 tree->times = mixt_tree->times;
66 tree = tree->next;
67 }
68 while(tree);
69
70
71 fp_xml_in = fopen(xml_filename,"r");
72 if(!fp_xml_in)
73 {
74 PhyML_Fprintf(stderr,"\n. Could not find the XML file '%s'.\n",xml_filename);
75 Exit("\n");
76 }
77
78 /* xroot = XML_Load_File(fp_xml_in); */
79 xroot = mixt_tree->xml_root;
80
81 if(xroot == NULL)
82 {
83 PhyML_Fprintf(stderr,"\n. Encountered an issue while loading the XML file.\n");
84 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
85 }
86
87 xnd = XML_Search_Node_Name("phyrex",NO,xroot);
88
89 if(xnd == NULL)
90 {
91 PhyML_Fprintf(stderr,"\n. Cound not find the \"root\" of the XML file (it should have \'phyrex\' as tag name).\n");
92 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
93 }
94
95 dum_string = XML_Get_Attribute_Value(xnd,"mcmc.chain.len");
96 if(dum_string != NULL) mixt_tree->io->mcmc->chain_len = (int)String_To_Dbl(dum_string);
97
98 dum_string = XML_Get_Attribute_Value(xnd,"mcmc.sample.every");
99 if(dum_string != NULL) mixt_tree->io->mcmc->sample_interval = (int)String_To_Dbl(dum_string);
100
101 dum_string = XML_Get_Attribute_Value(xnd,"mcmc.print.every");
102 if(dum_string != NULL) mixt_tree->io->mcmc->print_every = (int)String_To_Dbl(dum_string);
103
104 dum_string = XML_Get_Attribute_Value(xnd,"mcmc.burnin");
105 if(dum_string != NULL) mixt_tree->io->mcmc->chain_len_burnin = (int)String_To_Dbl(dum_string);
106
107
108 dum_string = XML_Get_Attribute_Value(xnd,"ignore.sequences");
109 if(!dum_string) dum_string = XML_Get_Attribute_Value(xnd,"ignore.seq");
110 if(!dum_string) dum_string = XML_Get_Attribute_Value(xnd,"ignore.data");
111
112 if(dum_string != NULL)
113 {
114 int select = XML_Validate_Attr_Int(dum_string,6,
115 "true","yes","y",
116 "false","no","n");
117 if(select < 3) mixt_tree->eval_alnL = NO;
118 else mixt_tree->eval_alnL = YES;
119 }
120
121 dum_string = XML_Get_Attribute_Value(xnd,"mutmap");
122 if(dum_string != NULL)
123 {
124 int select = XML_Validate_Attr_Int(dum_string,6,
125 "true","yes","y",
126 "false","no","n");
127 if(select < 3) mixt_tree->io->mutmap = YES;
128 else mixt_tree->io->mutmap = NO;
129 }
130
131
132 // Looking for XML node with rate-across-lineage info
133 xnd = XML_Search_Node_Name("lineagerates",YES,xroot);
134
135 if(xnd == NULL)
136 {
137 PhyML_Fprintf(stdout,"\n. The model of rate variation across lineages is not specified.");
138 PhyML_Fprintf(stdout,"\n. Using the geometric Brownian model (see Guindon, 2012, Syst. Biol.).\n");
139 mixt_tree->rates->model = GUINDON;
140 mixt_tree->mod->gamma_mgf_bl = YES;
141 strcpy(mixt_tree->rates->model_name,"geometric Brownian");
142 }
143 else
144 {
145 char *model_name;
146 model_name = XML_Get_Attribute_Value(xnd,"model");
147
148 if(model_name == NULL)
149 {
150 PhyML_Fprintf(stderr,"\n. Please specify a model of rate variation across lineages,");
151 PhyML_Fprintf(stderr,"\n. e.g., <lineagerates model=\"geometricbrownian\"/>.");
152 PhyML_Fprintf(stderr,"\n. See the manual for more options.");
153 assert(FALSE);
154 }
155 else
156 {
157 if(!strcmp(model_name,"geometricbrownian"))
158 {
159 mixt_tree->rates->model = GUINDON;
160 mixt_tree->mod->gamma_mgf_bl = YES;
161 strcpy(mixt_tree->rates->model_name,"geometric Brownian");
162 }
163 else if(!strcmp(model_name,"geometric"))
164 {
165 mixt_tree->rates->model = GUINDON;
166 mixt_tree->mod->gamma_mgf_bl = YES;
167 strcpy(mixt_tree->rates->model_name,"geometric Brownian");
168 }
169 else if(!strcmp(model_name,"brownian"))
170 {
171 mixt_tree->rates->model = GUINDON;
172 mixt_tree->mod->gamma_mgf_bl = YES;
173 strcpy(mixt_tree->rates->model_name,"geometric Brownian");
174 }
175 else if(!strcmp(model_name,"geo"))
176 {
177 mixt_tree->rates->model = GUINDON;
178 mixt_tree->mod->gamma_mgf_bl = YES;
179 strcpy(mixt_tree->rates->model_name,"geometric Brownian");
180 }
181 else if(!strcmp(model_name,"lognormal"))
182 {
183 mixt_tree->rates->model = LOGNORMAL;
184 mixt_tree->mod->gamma_mgf_bl = NO;
185 strcpy(mixt_tree->rates->model_name,"lognormal (uncorrelated)");
186 }
187 else if(!strcmp(model_name,"normal"))
188 {
189 mixt_tree->rates->model = LOGNORMAL;
190 mixt_tree->mod->gamma_mgf_bl = NO;
191 strcpy(mixt_tree->rates->model_name,"lognormal (uncorrelated)");
192 }
193 else if(!strcmp(model_name,"strictclock"))
194 {
195 mixt_tree->rates->model = STRICTCLOCK;
196 mixt_tree->mod->gamma_mgf_bl = NO;
197 strcpy(mixt_tree->rates->model_name,"strict clock");
198 }
199 else if(!strcmp(model_name,"clock"))
200 {
201 mixt_tree->rates->model = STRICTCLOCK;
202 mixt_tree->mod->gamma_mgf_bl = NO;
203 strcpy(mixt_tree->rates->model_name,"strict clock");
204 }
205 else
206 {
207 assert(FALSE);
208 }
209 }
210 }
211
212 // Spatial model
213 xnd = XML_Search_Node_Name("spatialmodel",YES,xroot);
214
215 if(xnd != NULL)
216 {
217 char *modname;
218 modname = XML_Get_Attribute_Value(xnd,"name");
219 XML_Set_Attribute_Value(xnd,"name",modname);
220
221 if(modname != NULL)
222 {
223 if(!strcmp(xnd->attr->value,"slfv")) mixt_tree->mmod->id = SLFV_GAUSSIAN;
224 else if(!strcmp(xnd->attr->value,"rw")) mixt_tree->mmod->id = RW;
225 else if(!strcmp(xnd->attr->value,"rrw")) mixt_tree->mmod->id = RRW;
226 else assert(FALSE);
227 }
228 }
229
230
231 // Looking for XML node with rate-across-lineage info
232 xnd = XML_Search_Node_Name("clockrate",YES,xroot);
233
234 if(xnd != NULL)
235 {
236 char *clock_r;
237 clock_r = XML_Get_Attribute_Value(xnd,"value");
238 if(clock_r == NULL) clock_r = XML_Get_Attribute_Value(xnd,"clock.val");
239 if(clock_r == NULL) clock_r = XML_Get_Attribute_Value(xnd,"val");
240 if(clock_r != NULL)
241 {
242 mixt_tree->rates->clock_r = String_To_Dbl(clock_r);
243 }
244
245 char *opt_clock;
246 opt_clock = XML_Get_Attribute_Value(xnd,"optimise.clock");
247 if(opt_clock == NULL) opt_clock = XML_Get_Attribute_Value(xnd,"optimize.clock");
248 if(opt_clock == NULL) opt_clock = XML_Get_Attribute_Value(xnd,"optimize.rate");
249 if(opt_clock == NULL) opt_clock = XML_Get_Attribute_Value(xnd,"opt.clock");
250
251 if(opt_clock != NULL)
252 {
253 int select = XML_Validate_Attr_Int(opt_clock,6,
254 "true","yes","y",
255 "false","no","n");
256 if(select < 3) mixt_tree->mod->s_opt->opt_clock_r = YES;
257 else mixt_tree->mod->s_opt->opt_clock_r = NO;
258 }
259 }
260
261
262 // Looking for calibration info
263 xnd = XML_Search_Node_Name("calibration",YES,xroot);
264
265 if(xnd == NULL)
266 {
267 PhyML_Fprintf(stderr,"\n. No calibration information seems to be provided.");
268 PhyML_Fprintf(stderr,"\n. Please amend your XML file. \n");
269 assert(FALSE);
270 }
271 else
272 {
273 assert(xnd->child);
274 if(XML_Search_Node_Name("upper",NO,xnd->child) == NULL && XML_Search_Node_Name("lower",NO,xnd->child) == NULL)
275 {
276 PhyML_Fprintf(stderr,"\n. There is no calibration information provided. \n");
277 PhyML_Fprintf(stderr,"\n. Please check your data. \n");
278 assert(FALSE);
279 }
280 }
281
282
283
284 // Looking for coordinate file
285 xnd = XML_Search_Node_Name("coordinates",YES,xroot);
286
287 if(xnd == NULL)
288 {
289 PhyML_Fprintf(stderr,"\n. No spatial information (i.e., coordinates) seems to be provided.");
290 PhyML_Fprintf(stderr,"\n. Please amend your XML file. \n");
291 Exit("\n");
292 }
293 else
294 {
295 char *coord_file;
296 coord_file = XML_Get_Attribute_Value(xnd,"file.name");
297
298 strcpy(mixt_tree->io->in_coord_file,coord_file);
299 mixt_tree->io->fp_in_coord = Openfile(mixt_tree->io->in_coord_file,READ);
300 }
301
302 seed = (mixt_tree->io->r_seed < 0)?(time(NULL)):(mixt_tree->io->r_seed);
303 srand(seed);
304 mixt_tree->io->r_seed = seed;
305
306
307 MIXT_Check_Model_Validity(mixt_tree);
308 MIXT_Init_Model(mixt_tree);
309 Print_Data_Structure(NO,stdout,mixt_tree);
310 tree = MIXT_Starting_Tree(mixt_tree);
311 if(mixt_tree->io->in_tree < 2) Add_Root(tree->a_edges[0],tree);
312 Copy_Tree(tree,mixt_tree);
313 Free_Tree(tree);
314 MIXT_Connect_Cseqs_To_Nodes(mixt_tree);
315 MIXT_Init_T_Beg(mixt_tree);
316 MIXT_Make_Tree_For_Lk(mixt_tree);
317 MIXT_Make_Tree_For_Pars(mixt_tree);
318 MIXT_Make_Spr(mixt_tree);
319 MIXT_Chain_All(mixt_tree);
320 MIXT_Check_Edge_Lens_In_All_Elem(mixt_tree);
321 MIXT_Turn_Branches_OnOff_In_All_Elem(ON,mixt_tree);
322 MIXT_Check_Invar_Struct_In_Each_Partition_Elem(mixt_tree);
323 MIXT_Check_RAS_Struct_In_Each_Partition_Elem(mixt_tree);
324
325
326 XML_Read_Calibration(xroot,mixt_tree);
327 MIXT_Chain_Cal(mixt_tree);
328
329
330 if(TIMES_Calibrations_Apply_To_Tips_Only(mixt_tree) == YES &&
331 mixt_tree->mod->s_opt->opt_topo == NO)
332 {
333 TIMES_Randomize_Tip_Times_Given_Calibrations(mixt_tree); // Topology is unchanged
334 TIMES_Bl_To_Times(mixt_tree);
335 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree);
336 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree);
337 }
338 else
339 {
340 TIMES_Randomize_Tree_With_Time_Constraints(mixt_tree->times->a_cal[0],mixt_tree);
341 }
342
343 MIXT_Propagate_Tree_Update(mixt_tree);
344
345 /* Create ldsks and connect tree tips to them */
346 /* once tip dates have been set properly (in */
347 /* TIMES_Randomize_Tree_With_Time_Constraints) */
348 PHYREX_Make_And_Connect_Tip_Disks(mixt_tree);
349
350
351 /* Read spatial coordinates */
352 PHYREX_Read_Tip_Coordinates(mixt_tree);
353 PHYREX_Init_Migrep_Mod(mixt_tree->mmod,2,
354 mixt_tree->mmod->lim_do->lonlat[0],
355 mixt_tree->mmod->lim_do->lonlat[1],
356 mixt_tree->mmod->lim_up->lonlat[0],
357 mixt_tree->mmod->lim_up->lonlat[1]);
358
359
360
361 /* Initialize parameters of migrep model */
362 mixt_tree->mmod->lbda = 1.;
363 mixt_tree->mmod->mu = 0.8;
364
365 if(mixt_tree->mmod->id == SLFV_GAUSSIAN || mixt_tree->mmod->id == SLFV_UNIFORM)
366 mixt_tree->mmod->rad = 0.1*((mixt_tree->mmod->lim_up->lonlat[0]-mixt_tree->mmod->lim_do->lonlat[0])+
367 (mixt_tree->mmod->lim_up->lonlat[1]-mixt_tree->mmod->lim_do->lonlat[1]));
368 else if(mixt_tree->mmod->id == RW || mixt_tree->mmod->id == RRW)
369 mixt_tree->mmod->rad = 2.0*((mixt_tree->mmod->lim_up->lonlat[0]-mixt_tree->mmod->lim_do->lonlat[0])+
370 (mixt_tree->mmod->lim_up->lonlat[1]-mixt_tree->mmod->lim_do->lonlat[1]));
371 else assert(FALSE);
372
373 mixt_tree->mmod->sigsq = PHYREX_Update_Sigsq(mixt_tree);
374
375
376 /* Random genealogy or user-defined tree */
377 switch(mixt_tree->io->in_tree)
378 {
379 case 0 : case 1 :
380 {
381 PHYREX_Simulate_Backward_Core(mixt_tree->young_disk,YES,mixt_tree);
382 PHYREX_Ldsk_To_Tree(mixt_tree);
383 break;
384 }
385 case 2:
386 {
387 PHYREX_Tree_To_Ldsk(mixt_tree);
388 break;
389 }
390 }
391
392 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree);
393 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree);
394
395 MIXT_Set_Ignore_Root(YES,mixt_tree);
396 MIXT_Set_Bl_From_Rt(YES,mixt_tree);
397
398 PHYREX_Oldest_Sampled_Disk(mixt_tree);
399
400 if(mixt_tree->mmod->id == RW || mixt_tree->mmod->id == RRW)
401 {
402 Make_Contrasts(mixt_tree);
403 mixt_tree->mmod->max_num_of_intervals = 3*tree->n_otu-2;
404 }
405
406 assert(PHYREX_Check_Struct(mixt_tree));
407 PHYREX_Lk(mixt_tree);
408 Set_Update_Eigen(YES,mixt_tree->mod);
409 Lk(NULL,mixt_tree);
410 Set_Update_Eigen(NO,mixt_tree->mod);
411 PhyML_Printf("\n. Init lnPr(seq|phylo): %f lnPr(coor|phylo): %f",mixt_tree->c_lnL,mixt_tree->mmod->c_lnL);
412 PhyML_Printf("\n. Random seed: %d",mixt_tree->io->r_seed);
413
414 res = PHYREX_MCMC(mixt_tree);
415 Free(res);
416
417 // Cleaning up...
418 RATES_Free_Rates(mixt_tree->rates);
419 RATES_Free_Rates(mixt_tree->extra_tree->rates);
420 TIMES_Free_Times(mixt_tree->times);
421 TIMES_Free_Times(mixt_tree->extra_tree->times);
422 MCMC_Free_MCMC(mixt_tree->mcmc);
423 MCMC_Free_MCMC(mixt_tree->extra_tree->mcmc);
424 Free_Mmod(mixt_tree->mmod);
425 Free_Spr_List_One_Edge(mixt_tree);
426 Free_Tree_Pars(mixt_tree);
427 Free_Tree_Lk(mixt_tree);
428
429 if(mixt_tree->io->fp_out_trees) fclose(mixt_tree->io->fp_out_trees);
430 if(mixt_tree->io->fp_out_tree) fclose(mixt_tree->io->fp_out_tree);
431 if(mixt_tree->io->fp_out_stats) fclose(mixt_tree->io->fp_out_stats);
432 if(mixt_tree->io->fp_out_json_trace) fclose(mixt_tree->io->fp_out_json_trace);
433 Free_Input(mixt_tree->io);
434
435
436 tree = mixt_tree;
437 do
438 {
439 Free_Calign(tree->data);
440 tree = tree->next_mixt;
441 }
442 while(tree);
443
444 tree = mixt_tree;
445 do
446 {
447 Free_Optimiz(tree->mod->s_opt);
448 tree = tree->next;
449 }
450 while(tree);
451
452
453 Free_Model_Complete(mixt_tree->mod);
454 Free_Model_Basic(mixt_tree->mod);
455 Free_Tree(mixt_tree->extra_tree);
456 Free_Tree(mixt_tree);
457 Free(res);
458 XML_Free_XML_Tree(xroot);
459 fclose(fp_xml_in);
460 }
461 #endif
462
463 //////////////////////////////////////////////////////////////
464 //////////////////////////////////////////////////////////////
465
PHYREX_Main_Simulate(int argc,char * argv[])466 int PHYREX_Main_Simulate(int argc, char *argv[])
467 {
468 t_tree *tree;
469 int seed,i,modid;
470 char *s;
471 t_dsk *disk;
472 int n_sites,n_otus;
473 phydbl width,height;
474 phydbl lbda, rad,mu;
475
476 s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
477
478 n_otus = (int)atoi(argv[1]);
479 n_sites = 1;
480 width = (phydbl)atof(argv[2]);
481 height = (phydbl)atof(argv[3]);
482 lbda = (phydbl)atof(argv[4]);
483 rad = (phydbl)atof(argv[5]);
484 mu = (phydbl)atof(argv[6]);
485 seed = (int)atoi(argv[7]);
486 modid = (int)atoi(argv[8]);
487
488 printf("\n. seed: %d",seed);
489 srand(seed);
490
491 tree = PHYREX_Simulate(n_otus,n_sites,width,height,lbda,rad,mu,seed,modid);
492 /* tree = PHYREX_Simulate_Independent_Loci(n_otus,500,20.,20.,seed); */
493
494 disk = tree->young_disk;
495 for(i=0;i<disk->n_ldsk_a;i++) Free_Ldisk(disk->ldsk_a[i]);
496 while(disk->prev)
497 {
498 disk = disk->prev;
499 if(disk->next->ldsk != NULL) Free_Ldisk(disk->next->ldsk);
500 Free_Disk(disk->next);
501 }
502
503 /* Root */
504 Free_Ldisk(disk->ldsk);
505 Free_Disk(disk);
506
507 RATES_Free_Rates(tree->rates);
508 /* MCMC_Free_MCMC(tree->mcmc); */
509 Free_Mmod(tree->mmod);
510 Free_Spr_List_One_Edge(tree);
511 Free_Spr_List_All_Edge(tree);
512 Free_Tree_Pars(tree);
513 Free_Tree_Lk(tree);
514 Free_Input(tree->io);
515 Free_Optimiz(tree->mod->s_opt);
516 Free_Model_Complete(tree->mod);
517 Free_Model_Basic(tree->mod);
518 Free_Calign(tree->data);
519 Free_Tree(tree);
520 Free(s);
521
522 return 0;
523 }
524
525 //////////////////////////////////////////////////////////////
526 //////////////////////////////////////////////////////////////
527
528
529 //////////////////////////////////////////////////////////////
530 //////////////////////////////////////////////////////////////
531 /* Test whether coord is in disk. Will actually only works in disk */
532 /* is a rectangle... */
533
PHYREX_Is_In_Disk(t_geo_coord * coord,t_dsk * disk,t_phyrex_mod * mmod)534 int PHYREX_Is_In_Disk(t_geo_coord *coord, t_dsk *disk, t_phyrex_mod *mmod)
535 {
536 int i;
537
538 assert(disk->centr->dim);
539
540 if(mmod->id == SLFV_UNIFORM)
541 {
542 for(i=0;i<disk->centr->dim;i++)
543 {
544 if(fabs(coord->lonlat[i] - disk->centr->lonlat[i]) > disk->mmod->rad + 1.E-20)
545 {
546 return(NO);
547 }
548 }
549 return(YES);
550 }
551 else if(mmod->id == SLFV_GAUSSIAN)
552 {
553 return(YES);
554 }
555
556 return(-1);
557 }
558
559 //////////////////////////////////////////////////////////////
560 //////////////////////////////////////////////////////////////
561
PHYREX_Lk(t_tree * tree)562 phydbl PHYREX_Lk(t_tree *tree)
563 {
564 switch(tree->mmod->id)
565 {
566 case SLFV_GAUSSIAN : case SLFV_UNIFORM :
567 {
568 return(SLFV_Lk_Gaussian(tree));
569 break;
570 }
571 case RW :
572 {
573 return(RW_Lk(tree));
574 break;
575 }
576 case RRW :
577 {
578 return(RRW_Lk(tree));
579 break;
580 }
581 default : assert(FALSE);
582 }
583 }
584
585 //////////////////////////////////////////////////////////////
586 //////////////////////////////////////////////////////////////
587
588
PHYREX_Lk_Core(t_dsk * disk,t_tree * tree)589 phydbl PHYREX_Lk_Core(t_dsk *disk, t_tree *tree)
590 {
591 switch(tree->mmod->id)
592 {
593 case SLFV_GAUSSIAN : case SLFV_UNIFORM :
594 {
595 return(SLFV_Lk_Gaussian_Core(disk,tree));
596 break;
597 }
598 case RW :
599 {
600 return(RW_Lk(tree));
601 break;
602 }
603 case RRW :
604 {
605 return(RRW_Lk(tree));
606 break;
607 }
608 default : assert(FALSE);
609 }
610 }
611
612 //////////////////////////////////////////////////////////////
613 //////////////////////////////////////////////////////////////
614 // Warning: the calculation below does not incoporate time information
615 // since things get messy when considering serial samples
PHYREX_Lk_Range(t_dsk * young,t_dsk * old,t_tree * tree)616 phydbl PHYREX_Lk_Range(t_dsk *young, t_dsk *old, t_tree *tree)
617 {
618 switch(tree->mmod->id)
619 {
620 case SLFV_GAUSSIAN : case SLFV_UNIFORM :
621 {
622 return(SLFV_Lk_Gaussian_Range(young,old,tree));
623 break;
624 }
625 case RW :
626 {
627 return(RW_Lk(tree));
628 break;
629 }
630 case RRW :
631 {
632 return(RRW_Lk(tree));
633 break;
634 }
635 default : assert(FALSE);
636 }
637
638 return(-1.);
639 }
640
641 //////////////////////////////////////////////////////////////
642 //////////////////////////////////////////////////////////////
643
644 #if (defined PHYREX)
PHYREX_MCMC(t_tree * tree)645 phydbl *PHYREX_MCMC(t_tree *tree)
646 {
647 t_mcmc *mcmc;
648 int move,i,n_vars;
649 phydbl u;
650 FILE *fp_tree,*fp_stats;
651 phydbl *res;
652 t_dsk *disk;
653
654 fp_tree = tree->io->fp_out_tree;
655 fp_stats = tree->io->fp_out_stats;
656
657 if(tree->io->mcmc == NULL)
658 {
659 mcmc = MCMC_Make_MCMC_Struct();
660 tree->mcmc = mcmc;
661 }
662 else
663 {
664 tree->mcmc = tree->io->mcmc;
665 mcmc = tree->mcmc;
666 }
667
668
669 mcmc->io = NULL;
670 mcmc->is = NO;
671 mcmc->run = 0;
672 mcmc->randomize = YES;
673 mcmc->norm_freq = 1E+3;
674 mcmc->max_tune = 1.E+20;
675 mcmc->is_burnin = YES;
676 mcmc->nd_t_digits = 1;
677 mcmc->max_lag = 1000;
678 mcmc->sample_size = mcmc->chain_len/mcmc->sample_interval;
679 mcmc->sample_num = 0;
680 disk = NULL;
681
682 PhyML_Fprintf(fp_tree,"#NEXUS");
683 PhyML_Fprintf(fp_tree,"\n\nBegin taxa;");
684 PhyML_Fprintf(fp_tree,"\n\tDimensions ntax=%d;",tree->n_otu);
685 PhyML_Fprintf(fp_tree,"\n\tTaxlabels");
686 for(int i=0;i<tree->n_otu;++i)
687 {
688 PhyML_Fprintf(fp_tree,"\n\t\t'%s'",tree->a_nodes[i]->name);
689 }
690 PhyML_Fprintf(fp_tree,"\n\t;");
691 PhyML_Fprintf(fp_tree,"\n\tend;");
692 PhyML_Fprintf(fp_tree,"\n\n");
693 PhyML_Fprintf(fp_tree,"\nBegin trees;");
694 PhyML_Fprintf(fp_tree,"\n\tTranslate");
695 for(int i=0;i<tree->n_otu;++i)
696 {
697 PhyML_Fprintf(fp_tree,"\n\t%d '%s'",i+1,tree->a_nodes[i]->name);
698 if(i<tree->n_otu-1) PhyML_Fprintf(fp_tree,",");
699 }
700 PhyML_Fprintf(fp_tree,"\n;");
701
702
703 MCMC_Complete_MCMC(mcmc,tree);
704
705 n_vars = 12;
706
707 res = (phydbl *)mCalloc(tree->mcmc->chain_len / tree->mcmc->sample_interval * n_vars,sizeof(phydbl));
708
709
710 PHYREX_Update_Sigsq(tree);
711
712 MIXT_Set_Bl_From_Rt(YES,tree);
713
714 if(XML_Search_Node_Attribute_Value("add","true",NO,tree->xml_root) == NULL)
715 {
716 PhyML_Fprintf(fp_stats,"\n# before rand glnL: %f alnL: %f",tree->mmod->c_lnL,tree->c_lnL);
717 PhyML_Fprintf(fp_stats,"\n# ninter: %d",PHYREX_Total_Number_Of_Intervals(tree));
718 PhyML_Fprintf(fp_stats,"\n# ncoal: %d",PHYREX_Total_Number_Of_Coal_Disks(tree));
719 PhyML_Fprintf(fp_stats,"\n# nhits: %d",PHYREX_Total_Number_Of_Hit_Disks(tree));
720 PhyML_Fprintf(fp_stats,"\n# true lbda: %f",tree->mmod->lbda);
721 PhyML_Fprintf(fp_stats,"\n# true mu: %f",tree->mmod->mu);
722 PhyML_Fprintf(fp_stats,"\n# true rad: %f",SLFV_Update_Radius(tree));
723 PhyML_Fprintf(fp_stats,"\n# true sigsq: %f",tree->mmod->sigsq);
724 PhyML_Fprintf(fp_stats,"\n# true neigh. size: %f",SLFV_Neighborhood_Size(tree));
725 PhyML_Fprintf(fp_stats,"\n# fst-based estimate of neighborhood size: %f",SLFV_Neighborhood_Size_Regression(tree));
726 PhyML_Fprintf(fp_stats,"\n# nucleotide diversity: %f",Nucleotide_Diversity(tree->data));
727 PhyML_Fprintf(fp_stats,"\n# length of a generation: %G time units",SLFV_Generation_Length(tree));
728 PhyML_Fprintf(fp_stats,"\n# clock rate: %G subst. per time unit",tree->rates->clock_r);
729 PhyML_Fprintf(fp_stats,"\n# after rand glnL: %f alnL: %f",tree->mmod->c_lnL,tree->c_lnL);
730 PhyML_Fprintf(fp_stats,"\n# ninter: %d",PHYREX_Total_Number_Of_Intervals(tree));
731 PhyML_Fprintf(fp_stats,"\n# ncoal: %d",PHYREX_Total_Number_Of_Coal_Disks(tree));
732 PhyML_Fprintf(fp_stats,"\n# nhits: %d",PHYREX_Total_Number_Of_Hit_Disks(tree));
733 PhyML_Fprintf(fp_stats,"\n# start lbda: %f",tree->mmod->lbda);
734 PhyML_Fprintf(fp_stats,"\n# start mu: %f",tree->mmod->mu);
735 PhyML_Fprintf(fp_stats,"\n# start rad: %f",tree->mmod->rad);
736 PhyML_Fprintf(fp_stats,"\n# dist. in tree: ");
737 for(int i=0;i<tree->n_otu-1;++i) for(int j=i+1;j<tree->n_otu;++j) PhyML_Fprintf(fp_stats,"%G ",PHYREX_Dist_Between_Two_Ldsk(tree->a_nodes[i]->ldsk,tree->a_nodes[j]->ldsk,tree));
738 PhyML_Fprintf(fp_stats,"\n# dist. in space: ");
739 for(int i=0;i<tree->n_otu-1;++i) for(int j=i+1;j<tree->n_otu;++j) PhyML_Fprintf(fp_stats,"%G ",Euclidean_Dist(tree->a_nodes[i]->ldsk->coord,tree->a_nodes[j]->ldsk->coord));
740 PhyML_Printf("\n");
741
742 fflush(NULL);
743
744 PhyML_Fprintf(fp_stats,"\n");
745 PhyML_Fprintf(fp_stats,"%s\t","sample");
746 PhyML_Fprintf(fp_stats,"%s\t","lnP");
747 PhyML_Fprintf(fp_stats,"%s\t","alnL");
748 PhyML_Fprintf(fp_stats,"%s\t","glnL");
749 PhyML_Fprintf(fp_stats,"%s\t","clock");
750 PhyML_Fprintf(fp_stats,"%s\t","evolrate");
751 PhyML_Fprintf(fp_stats,"%s\t","lbda");
752 PhyML_Fprintf(fp_stats,"%s\t","mu");
753 PhyML_Fprintf(fp_stats,"%s\t","rad");
754 PhyML_Fprintf(fp_stats,"%s\t","sigsq");
755 PhyML_Fprintf(fp_stats,"%s\t","neff");
756 PhyML_Fprintf(fp_stats,"%s\t","neigh");
757 PhyML_Fprintf(fp_stats,"%s\t","rhoe");
758 PhyML_Fprintf(fp_stats,"%s\t","realsigsqroot");
759 PhyML_Fprintf(fp_stats,"%s\t","realsigsqtips");
760 PhyML_Fprintf(fp_stats,"%s\t","realsigsqtipsbis");
761 PhyML_Fprintf(fp_stats,"%s\t","realsigsqtipster");
762 PhyML_Fprintf(fp_stats,"%s\t","dispdist");
763 PhyML_Fprintf(fp_stats,"%s\t","nInt");
764 PhyML_Fprintf(fp_stats,"%s\t","nCoal");
765 PhyML_Fprintf(fp_stats,"%s\t","nHit");
766 PhyML_Fprintf(fp_stats,"%s\t","rootTime");
767 PhyML_Fprintf(fp_stats,"%s\t","rootLon");
768 PhyML_Fprintf(fp_stats,"%s\t","rootLat");
769 for(int i=0;i<Scalar_Len(tree->mod->kappa);++i) PhyML_Fprintf(fp_stats,"tstv%d\t",i);
770 if(tree->mod->ras->free_mixt_rates == NO) PhyML_Fprintf(fp_stats,"alpha\t");
771 else
772 {
773 for(int i=0;i<tree->mod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"p(%d)\t",i+1);
774 for(int i=0;i<tree->mod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"rr(%d)\t",i+1);
775 }
776
777 PhyML_Fprintf(fp_stats,"%s\t","MeanBr");
778 PhyML_Fprintf(fp_stats,"%s\t","TreeLen");
779 PhyML_Fprintf(fp_stats,"%s\t","accLbda");
780 PhyML_Fprintf(fp_stats,"%s\t","accMu");
781 PhyML_Fprintf(fp_stats,"%s\t","accRad");
782 PhyML_Fprintf(fp_stats,"%s\t","accInDelDisk");
783 PhyML_Fprintf(fp_stats,"%s\t","accInDelHit");
784 PhyML_Fprintf(fp_stats,"%s\t","accScaleTime");
785 PhyML_Fprintf(fp_stats,"%s\t","accSPR");
786 PhyML_Fprintf(fp_stats,"%s\t","accSPRlocal");
787 PhyML_Fprintf(fp_stats,"%s\t","accPath");
788 PhyML_Fprintf(fp_stats,"%s\t","accIndelDiskSerial");
789 PhyML_Fprintf(fp_stats,"%s\t","accIndelHitSerial");
790 PhyML_Fprintf(fp_stats,"%s\t","accLdskGivenDisk");
791 PhyML_Fprintf(fp_stats,"%s\t","accDiskGivenLdsk");
792 PhyML_Fprintf(fp_stats,"%s\t","accDiskAndLdsk");
793 PhyML_Fprintf(fp_stats,"%s\t","accLdskMulti");
794 PhyML_Fprintf(fp_stats,"%s\t","accDiskMulti");
795 PhyML_Fprintf(fp_stats,"%s\t","accMoveDiskUD");
796 PhyML_Fprintf(fp_stats,"%s\t","accAddRemoveJump");
797 PhyML_Fprintf(fp_stats,"%s\t","accLdskTipToRoot");
798 PhyML_Fprintf(fp_stats,"%s\t","accSigsqScale");
799 PhyML_Fprintf(fp_stats,"%s\t","tuneLbda");
800 PhyML_Fprintf(fp_stats,"%s\t","tuneRad");
801 PhyML_Fprintf(fp_stats,"%s\t","tuneMu");
802 PhyML_Fprintf(fp_stats,"%s\t","tuneIndelDisk");
803 PhyML_Fprintf(fp_stats,"%s\t","tuneIndelHit");
804 PhyML_Fprintf(fp_stats,"%s\t","tuneLdskGivenDisk");
805 PhyML_Fprintf(fp_stats,"%s\t","tuneIndelDiskSerial");
806 PhyML_Fprintf(fp_stats,"%s\t","tuneIndelHitSerial");
807 for(int i=0;i<2*tree->n_otu-2;++i) PhyML_Fprintf(fp_stats,"s%d\t",i);
808 }
809
810
811 for(i=0;i<mcmc->n_moves;i++) tree->mcmc->start_ess[i] = YES;
812
813 PHYREX_Lk(tree);
814 Set_Update_Eigen(YES,tree->mod);
815 Lk(NULL,tree);
816 Set_Update_Eigen(NO,tree->mod);
817 RATES_Lk_Rates(tree);
818
819 if(isnan(tree->c_lnL) || isinf(tree->c_lnL))
820 {
821 PhyML_Fprintf(stderr,"\n. Cannot compute sequence log-likelihood. Aborting.");
822 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
823 }
824
825 if(isnan(tree->mmod->c_lnL) || isinf(tree->mmod->c_lnL))
826 {
827 PhyML_Fprintf(stderr,"\n. Cannot compute location log-likelihood. Aborting.");
828 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
829 }
830
831
832 Set_Both_Sides(NO,tree);
833 mcmc->always_yes = NO;
834 move = -1;
835 do
836 {
837 MIXT_Propagate_Tree_Update(tree);
838 assert(PHYREX_Check_Struct(tree));
839
840 if(mcmc->run > mcmc->chain_len_burnin)
841 for(i=0;i<mcmc->n_moves;i++) tree->mcmc->adjust_tuning[i] = NO;
842 else
843 {
844 for(i=0;i<mcmc->n_moves;i++) tree->mcmc->adjust_tuning[i] = YES;
845 tree->mcmc->is_burnin = NO;
846 }
847
848 u = Uni();
849
850 for(move=0;move<tree->mcmc->n_moves;move++) if(tree->mcmc->move_weight[move] > u-1.E-10) break;
851
852 tree->mcmc->move_idx = move;
853
854 assert(!(move == tree->mcmc->n_moves));
855
856 if(!(tree->mcmc->run%tree->mcmc->print_every))
857 {
858 PhyML_Fprintf(stdout,"\n. %10d %30s %20f %20f %20f",
859 tree->mcmc->run,
860 tree->mcmc->move_name[move],
861 tree->mmod->c_lnL,
862 tree->c_lnL,
863 tree->mmod->c_lnL+tree->c_lnL+tree->rates->c_lnL_rates);
864 if(tree->numerical_warning == YES) PhyML_Fprintf(stdout," -- WARNING: numerical precision issue detected...");
865 }
866
867
868 if(!(tree->mcmc->run%tree->mcmc->sample_interval))
869 {
870 disk = tree->young_disk;
871 while(disk->prev) disk = disk->prev;
872
873 PhyML_Fprintf(fp_stats,"\n");
874 PhyML_Fprintf(fp_stats,"%6d\t",tree->mcmc->run);
875 PhyML_Fprintf(fp_stats,"%g\t",tree->c_lnL+tree->mmod->c_lnL+tree->rates->c_lnL_rates);
876 PhyML_Fprintf(fp_stats,"%g\t",tree->c_lnL);
877 PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->c_lnL);
878 PhyML_Fprintf(fp_stats,"%g\t",tree->rates->clock_r);
879 PhyML_Fprintf(fp_stats,"%g\t",RATES_Realized_Substitution_Rate(tree));
880 PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->lbda);
881 PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->mu);
882 PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->rad);
883 PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->sigsq);
884 PhyML_Fprintf(fp_stats,"%g\t",tree->times->scaled_pop_size);
885 PhyML_Fprintf(fp_stats,"%g\t",SLFV_Neighborhood_Size(tree));
886 PhyML_Fprintf(fp_stats,"%g\t",SLFV_Effective_Density(tree));
887 PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Root_To_Tip_Realized_Sigsq(tree));
888 PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Tip_To_Root_Realized_Sigsq(tree));
889 PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Tip_To_Root_Realized_Bis_Sigsq(tree));
890 PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Tip_To_Root_Realized_Ter_Sigsq(tree));
891 PhyML_Fprintf(fp_stats,"%g\t",PHYREX_Realized_Dispersal_Dist(tree));
892 PhyML_Fprintf(fp_stats,"%d\t",PHYREX_Total_Number_Of_Intervals(tree));
893 PhyML_Fprintf(fp_stats,"%d\t",PHYREX_Total_Number_Of_Coal_Disks(tree));
894 PhyML_Fprintf(fp_stats,"%d\t",PHYREX_Total_Number_Of_Hit_Disks(tree));
895 PhyML_Fprintf(fp_stats,"%g\t",disk->time);
896 PhyML_Fprintf(fp_stats,"%g\t",disk->ldsk->coord->lonlat[0]);
897 PhyML_Fprintf(fp_stats,"%g\t",disk->ldsk->coord->lonlat[1]);
898 Output_Scalar_Dbl(tree->mod->kappa,"\t",fp_stats);
899 if(tree->mod->ras->free_mixt_rates == NO) PhyML_Fprintf(fp_stats,"%g\t",tree->mod->ras->alpha->v);
900 else
901 {
902 for(int i=0;i<tree->mod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mod->ras->gamma_r_proba->v[i]);
903 for(int i=0;i<tree->mod->ras->n_catg;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mod->ras->gamma_rr->v[i]);
904 }
905 PhyML_Fprintf(fp_stats,"%g\t",RATES_Get_Mean_Rate_In_Subtree(tree->n_root,tree));
906 PhyML_Fprintf(fp_stats,"%g\t",Tree_Length(tree));
907 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_lbda]);
908 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_mu]);
909 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_rad]);
910 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_indel_disk]);
911 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_indel_hit]);
912 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_scale_times]);
913 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_spr]);
914 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_spr_local]);
915 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_traj]);
916 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_indel_disk_serial]);
917 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_indel_hit_serial]);
918 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_ldsk_given_disk]);
919 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_disk_given_ldsk]);
920 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_ldsk_and_disk]);
921 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_ldsk_multi]);
922 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_disk_multi]);
923 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_move_disk_ud]);
924 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_add_remove_jump]);
925 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_ldsk_tip_to_root]);
926 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->acc_rate[tree->mcmc->num_move_phyrex_sigsq_scale]);
927 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_lbda]);
928 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_rad]);
929 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_mu]);
930 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_indel_disk]);
931 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_indel_hit]);
932 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_ldsk_given_disk]);
933 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_indel_disk_serial]);
934 PhyML_Fprintf(fp_stats,"%g\t",tree->mcmc->tune_move[tree->mcmc->num_move_phyrex_indel_hit_serial]);
935 for(int i=0;i<2*tree->n_otu-2;++i) PhyML_Fprintf(fp_stats,"%g\t",tree->mmod->sigsq_scale[i]);
936
937 /* res[0 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = tree->mmod->lbda; */
938 /* res[1 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = tree->mmod->mu; */
939 /* res[2 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Update_Sigsq(tree); */
940 /* res[3 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Neighborhood_Size(tree); */
941 /* res[4 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = tree->mmod->rad; */
942 /* res[5 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Total_Number_Of_Intervals(tree); */
943 /* res[6 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Total_Number_Of_Coal_Disks(tree); */
944 /* res[7 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Total_Number_Of_Hit_Disks(tree); */
945 /* res[8 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Effective_Density(tree); */
946 /* res[9 * tree->mcmc->chain_len / tree->mcmc->sample_interval + tree->mcmc->run / tree->mcmc->sample_interval] = PHYREX_Coalescence_Rate(tree); */
947
948 /* MCMC_Copy_To_New_Param_Val(tree->mcmc,tree); */
949
950 /* for(i=0;i<tree->mcmc->n_moves;i++) if(tree->mcmc->start_ess[i] == YES) MCMC_Update_Effective_Sample_Size(i,tree->mcmc,tree); */
951 /* for(i=0;i<tree->mcmc->n_moves;i++) MCMC_Update_Mode(i,tree->mcmc,tree); */
952
953
954 PHYREX_Ldsk_To_Tree(tree);
955 TIMES_Time_To_Bl(tree);
956 tree->bl_ndigits = 3;
957 tree->write_tax_names = NO;
958 PHYREX_Label_Nodes_With_Locations(tree);
959 PHYREX_Label_Edges(tree);
960 char *s = Write_Tree(tree);
961 PhyML_Fprintf(fp_tree,"\ntree %d [&lnP=%f,precision={1.e-1,1e-02,1e-01}] = [&R] %s",tree->mcmc->run,tree->c_lnL,s);
962 PhyML_Fprintf(fp_tree,"\nend;");
963 fseek(fp_tree,-5,SEEK_END);
964 tree->write_tax_names = YES;
965 Free(s);
966
967 fflush(NULL);
968
969 tree->mcmc->sample_num++;
970 }
971
972 if(!strcmp(tree->mcmc->move_name[move],"phyrex_lbda"))
973 MCMC_PHYREX_Lbda(tree);
974
975 if(!strcmp(tree->mcmc->move_name[move],"phyrex_mu"))
976 MCMC_PHYREX_Mu(tree);
977
978 if(!strcmp(tree->mcmc->move_name[move],"phyrex_rad"))
979 MCMC_PHYREX_Radius(tree);
980
981 if(!strcmp(tree->mcmc->move_name[move],"phyrex_sigsq"))
982 MCMC_PHYREX_Sigsq(tree);
983
984 if(!strcmp(tree->mcmc->move_name[move],"phyrex_neff"))
985 MCMC_PHYREX_Neff(tree);
986
987 if(!strcmp(tree->mcmc->move_name[move],"phyrex_indel_disk"))
988 MCMC_PHYREX_Indel_Disk(tree);
989
990 if(!strcmp(tree->mcmc->move_name[move],"phyrex_indel_hit"))
991 MCMC_PHYREX_Indel_Hit(tree);
992
993 if(!strcmp(tree->mcmc->move_name[move],"phyrex_move_disk_ud"))
994 MCMC_PHYREX_Move_Disk_Updown(tree);
995
996 if(!strcmp(tree->mcmc->move_name[move],"phyrex_swap_disk"))
997 MCMC_PHYREX_Swap_Disk(tree);
998
999 if(!strcmp(tree->mcmc->move_name[move],"phyrex_scale_times"))
1000 MCMC_PHYREX_Scale_Times(tree);
1001
1002 if(!strcmp(tree->mcmc->move_name[move],"phyrex_spr"))
1003 MCMC_PHYREX_Prune_Regraft(tree);
1004
1005 if(!strcmp(tree->mcmc->move_name[move],"phyrex_spr_local"))
1006 MCMC_PHYREX_Prune_Regraft_Local(tree);
1007
1008 if(!strcmp(tree->mcmc->move_name[move],"phyrex_traj"))
1009 MCMC_PHYREX_Lineage_Traj(tree);
1010
1011 if(!strcmp(tree->mcmc->move_name[move],"phyrex_disk_multi"))
1012 MCMC_PHYREX_Disk_Multi(tree);
1013
1014 if(!strcmp(tree->mcmc->move_name[move],"phyrex_ldsk_multi"))
1015 MCMC_PHYREX_Ldsk_Multi(tree);
1016
1017 if(!strcmp(tree->mcmc->move_name[move],"phyrex_ldsk_and_disk"))
1018 MCMC_PHYREX_Ldsk_And_Disk(tree);
1019
1020 if(!strcmp(tree->mcmc->move_name[move],"phyrex_ldsk_tip_to_root"))
1021 MCMC_PHYREX_Ldsk_Tip_To_Root(tree);
1022
1023 if(!strcmp(tree->mcmc->move_name[move],"phyrex_ldsk_given_disk"))
1024 MCMC_PHYREX_Ldsk_Given_Disk(tree);
1025
1026 if(!strcmp(tree->mcmc->move_name[move],"phyrex_disk_given_ldsk"))
1027 MCMC_PHYREX_Disk_Given_Ldsk(tree);
1028
1029 if(!strcmp(tree->mcmc->move_name[move],"phyrex_indel_disk_serial"))
1030 MCMC_PHYREX_Indel_Disk_Serial(tree);
1031
1032 if(!strcmp(tree->mcmc->move_name[move],"phyrex_indel_hit_serial"))
1033 MCMC_PHYREX_Indel_Hit_Serial(tree);
1034
1035 if(!strcmp(tree->mcmc->move_name[move],"phyrex_add_remove_jump"))
1036 MCMC_PHYREX_Add_Remove_Jump(tree);
1037
1038 if(!strcmp(tree->mcmc->move_name[move],"phyrex_sigsq_scale"))
1039 MCMC_PHYREX_Sigsq_Scale(tree);
1040
1041 if(!strcmp(tree->mcmc->move_name[move],"kappa"))
1042 MCMC_Kappa(tree);
1043
1044 if(!strcmp(tree->mcmc->move_name[move],"rr"))
1045 MCMC_RR(tree);
1046
1047 if(!strcmp(tree->mcmc->move_name[move],"ras"))
1048 MCMC_Rate_Across_Sites(tree);
1049
1050 if(!strcmp(tree->mcmc->move_name[move],"br_rate"))
1051 MCMC_Rates_All(tree);
1052
1053 if(!strcmp(tree->mcmc->move_name[move],"tree_rates"))
1054 MCMC_Tree_Rates(tree);
1055
1056 if(!strcmp(tree->mcmc->move_name[move],"clock"))
1057 MCMC_Clock_R(tree);
1058
1059
1060 /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
1061 /* PHYREX_Lk(tree); */
1062 /* Lk(NULL,tree); */
1063
1064 /* PhyML_Printf("\n. Move: %s tree->mmod->c_lnL: %f kappa: %f", */
1065 /* tree->mcmc->move_name[move], */
1066 /* tree->mmod->c_lnL, */
1067 /* tree->mod->kappa->v); */
1068
1069 if(tree->mmod->c_lnL < UNLIKELY || tree->c_lnL < UNLIKELY)
1070 {
1071 PhyML_Printf("\n. Move: %s tree->mmod->c_lnL: %f tree->c_lnL: %f",
1072 tree->mcmc->move_name[move],
1073 tree->mmod->c_lnL,tree->c_lnL);
1074 assert(FALSE);
1075 }
1076
1077 if(tree->mmod->safe_phyrex == YES)
1078 {
1079 /* phydbl c_lnL = tree->c_lnL; */
1080 /* Lk(NULL,tree); */
1081 /* if(Are_Equal(c_lnL,tree->c_lnL,1.E-5) == NO) */
1082 /* { */
1083 /* PhyML_Fprintf(stderr,"\n. Problem detected with move %s",tree->mcmc->move_name[move]); */
1084 /* PhyML_Fprintf(stderr,"\n. c_lnL: %f -> %f",c_lnL,tree->c_lnL); */
1085 /* Generic_Exit(__FILE__,__LINE__,__FUNCTION__); */
1086 /* } */
1087
1088 phydbl g_lnL = tree->mmod->c_lnL;
1089 PHYREX_Lk(tree);
1090 if(Are_Equal(g_lnL,tree->mmod->c_lnL,1.E-5) == NO)
1091 {
1092 PhyML_Fprintf(stderr,"\n. Problem detected with move %s. Iteration %d",tree->mcmc->move_name[move],tree->mcmc->run);
1093 PhyML_Fprintf(stderr,"\n. g_lnL: %f -> %f [%g]",g_lnL,tree->mmod->c_lnL,g_lnL-tree->mmod->c_lnL);
1094 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1095 }
1096
1097 /* phydbl r_lnL = tree->rates->c_lnL_rates; */
1098 /* RATES_Lk_Rates(tree); */
1099 /* if(Are_Equal(r_lnL,tree->rates->c_lnL_rates,1.E-5) == NO) */
1100 /* { */
1101 /* PhyML_Fprintf(stderr,"\n. Problem detected with move %s",tree->mcmc->move_name[move]); */
1102 /* PhyML_Fprintf(stderr,"\n. r_lnL: %f -> %f [%g]",r_lnL,tree->rates->c_lnL_rates,r_lnL-tree->rates->c_lnL_rates); */
1103 /* Generic_Exit(__FILE__,__LINE__,__FUNCTION__); */
1104 /* } */
1105 }
1106
1107 tree->mcmc->run++;
1108 MCMC_Get_Acc_Rates(tree->mcmc);
1109
1110
1111
1112 (void)signal(SIGINT,MCMC_Terminate);
1113 }
1114 while(tree->mcmc->run < tree->mcmc->chain_len);
1115
1116 PhyML_Fprintf(stdout,"\n. The analysis completed !");
1117
1118 fclose(fp_tree);
1119 fclose(fp_stats);
1120 /* fclose(fp_summary); */
1121
1122 return(res);
1123 }
1124 #endif
1125
1126 //////////////////////////////////////////////////////////////
1127 //////////////////////////////////////////////////////////////
1128
PHYREX_Wrap_Lk(t_edge * b,t_tree * tree,supert_tree * stree)1129 phydbl PHYREX_Wrap_Lk(t_edge *b, t_tree *tree, supert_tree *stree)
1130 {
1131 return PHYREX_Lk(tree);
1132 }
1133
1134 /*////////////////////////////////////////////////////////////
1135 ////////////////////////////////////////////////////////////*/
1136
PHYREX_Remove_Disk(t_dsk * disk)1137 void PHYREX_Remove_Disk(t_dsk *disk)
1138 {
1139 t_dsk *prev;
1140 t_dsk *next;
1141
1142 prev = disk->prev;
1143 next = disk->next;
1144
1145 assert(next != NULL);
1146
1147 if(prev != NULL) prev->next = next;
1148 next->prev = prev;
1149 }
1150
1151 /*////////////////////////////////////////////////////////////
1152 ////////////////////////////////////////////////////////////*/
1153
1154 /* Insert disk event based on its time. Insertion above the root is permitted
1155 */
PHYREX_Insert_Disk(t_dsk * ins,t_tree * tree)1156 void PHYREX_Insert_Disk(t_dsk *ins, t_tree *tree)
1157 {
1158 t_dsk *disk;
1159
1160 assert(ins != NULL);
1161
1162 disk = tree->young_disk;
1163 while(disk->prev != NULL && disk->prev->time > ins->time) disk = disk->prev;
1164
1165 ins->prev = disk->prev;
1166 ins->next = disk;
1167 disk->prev = ins;
1168 if(ins->prev != NULL) ins->prev->next = ins;
1169 }
1170
1171 /*////////////////////////////////////////////////////////////
1172 ////////////////////////////////////////////////////////////*/
1173
PHYREX_Prev_Coal_Lindisk(t_ldsk * t)1174 t_ldsk *PHYREX_Prev_Coal_Lindisk(t_ldsk *t)
1175 {
1176 if(t == NULL) return NULL;
1177
1178 if(t->n_next > 1)
1179 {
1180 return t;
1181 }
1182 else
1183 {
1184 return PHYREX_Prev_Coal_Lindisk(t->prev);
1185 }
1186 }
1187
1188 //////////////////////////////////////////////////////////////
1189 //////////////////////////////////////////////////////////////
1190
PHYREX_Next_Coal_Lindisk(t_ldsk * t)1191 t_ldsk *PHYREX_Next_Coal_Lindisk(t_ldsk *t)
1192 {
1193 assert(!(t == NULL));
1194
1195 if(t->n_next > 1 || t->next == NULL) return t;
1196 else
1197 {
1198 if(t->n_next > 1) // Should have t->is_coal = YES
1199 {
1200 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1201 Warn_And_Exit("");
1202 }
1203 return PHYREX_Next_Coal_Lindisk(t->next[0]);
1204 }
1205 }
1206
1207 /*////////////////////////////////////////////////////////////
1208 ////////////////////////////////////////////////////////////*/
1209
1210 /* Generate a new trajectory, including disk event centers, between
1211 'y_ldsk' a ``young'' lindisk event and 'o_ldsk' an old one. 'y_ldsk
1212 and 'o_ldsk' remain unaffected. No disk events should be present
1213 between y_ldsk and o_ldsk: we need to generate some first. n_cur_disk
1214 is the current number of disks between y_ldsk and o_ldsk.
1215 */
PHYREX_One_New_Traj(t_ldsk * y_ldsk,t_ldsk * o_ldsk,int dir_o_y,t_dsk * xtra_dsk,int n_cur_disk,t_tree * tree)1216 int PHYREX_One_New_Traj(t_ldsk *y_ldsk, t_ldsk *o_ldsk, int dir_o_y, t_dsk *xtra_dsk, int n_cur_disk, t_tree *tree)
1217 {
1218 t_phyrex_mod *mmod;
1219 t_dsk *disk,**disk_new;
1220 int i,n,K;
1221 int min_n_disk,n_new_disk,n_disk;
1222
1223 mmod = tree->mmod;
1224 disk = NULL;
1225 disk_new = NULL;
1226 K = 2;
1227
1228 /* printf("\n# New traj from %s to %s",y_ldsk->coord->id,o_ldsk->coord->id); */
1229 /* fflush(NULL); */
1230
1231 /* Minimum number of disks between y_ldsk and o_ldsk */
1232 min_n_disk = 0;
1233 for(i=0;i<mmod->n_dim;i++)
1234 {
1235 /* PhyML_Printf("\n# y_ldsk %s : %f o_ldsk->disk->centr: %f rad: %f", */
1236 /* y_ldsk->coord->id, */
1237 /* y_ldsk->coord->lonlat[i], */
1238 /* o_ldsk->disk->centr->lonlat[i], */
1239 /* mmod->rad); */
1240 if(y_ldsk->coord->lonlat[i] < o_ldsk->disk->centr->lonlat[i])
1241 {
1242 n_disk = 0;
1243 while(y_ldsk->coord->lonlat[i] + (2*n_disk-1)*mmod->rad < o_ldsk->disk->centr->lonlat[i] - mmod->rad) n_disk++;
1244 if(n_disk > min_n_disk) min_n_disk = n_disk;
1245 }
1246 else
1247 {
1248 n_disk = 0;
1249 while(y_ldsk->coord->lonlat[i] - (2*n_disk-1)*mmod->rad > o_ldsk->disk->centr->lonlat[i] + mmod->rad) n_disk++;
1250 if(n_disk > min_n_disk) min_n_disk = n_disk;
1251 }
1252 /* printf(" -- min_n_disk: %d",min_n_disk); */
1253 }
1254
1255 /* printf("\n# min_n_disk: %d cur_n_disk: %d",min_n_disk,n_cur_disk); */
1256 /* fflush(NULL); */
1257
1258 /* How many disks along the new path between y_ldsk and o_ldsk */
1259 n_new_disk = Rand_Int(n_cur_disk-K,n_cur_disk+K);
1260 if(n_new_disk < min_n_disk) n_new_disk = min_n_disk;
1261
1262 if(xtra_dsk != NULL) n_new_disk++;
1263
1264 /* printf("\n# Add n_new_disk: %d",n_new_disk); fflush(NULL); */
1265
1266 if(n_new_disk > 0)
1267 {
1268 /* Make new disks to create a new path between ldsk_left and ldsk_up */
1269 disk_new = (t_dsk **)mCalloc(n_new_disk,sizeof(t_dsk *));
1270 for(i=0;i<n_new_disk-1;i++) disk_new[i] = PHYREX_Make_Disk_Event(mmod->n_dim,tree->n_otu);
1271 if(xtra_dsk != NULL) disk_new[n_new_disk-1] = xtra_dsk;
1272 else disk_new[n_new_disk-1] = PHYREX_Make_Disk_Event(mmod->n_dim,tree->n_otu);
1273
1274 for(i=0;i<n_new_disk;i++) PHYREX_Init_Disk_Event(disk_new[i],mmod->n_dim,mmod);
1275
1276 /* Times of these new disks. If xtra_dsk != NULL, then make sure you do not */
1277 /* reset the time of that disk */
1278 n = (xtra_dsk != NULL) ? (n_new_disk-1) : (n_new_disk);
1279 for(i=0;i<n;i++)
1280 disk_new[i]->time =
1281 Uni()*(y_ldsk->disk->time - o_ldsk->disk->time) + o_ldsk->disk->time;
1282
1283 /* Insert these events */
1284 for(i=0;i<n_new_disk;i++)
1285 {
1286 assert(!tree->young_disk->next);
1287 disk = tree->young_disk;
1288 while(disk->time > disk_new[i]->time) disk = disk->prev;
1289 PHYREX_Insert_Disk(disk_new[i],tree);
1290 }
1291
1292 /* for(i=0;i<n_new_disk;i++) */
1293 /* { */
1294 /* printf("\n> disk_new: %f [%s]",disk_new[i]->time,disk_new[i]->id); fflush(NULL); */
1295 /* } */
1296
1297 /* Add new lindisks to the new disk events */
1298 for(i=0;i<n_new_disk;i++)
1299 {
1300 disk_new[i]->ldsk = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
1301 PHYREX_Init_Lindisk_Node(disk_new[i]->ldsk,disk_new[i],tree->mmod->n_dim);
1302 PHYREX_Make_Lindisk_Next(disk_new[i]->ldsk);
1303 /* printf("\n# Add ldsk %s to %s",disk_new[i]->ldsk->coord->id,disk_new[i]->id); fflush(NULL); */
1304 }
1305
1306 /* Connect them */
1307 PHYREX_Connect_Ldsk_Given_Disk(disk_new,n_new_disk,y_ldsk,o_ldsk,dir_o_y);
1308
1309 Free(disk_new);
1310 }
1311 else
1312 {
1313 o_ldsk->next[dir_o_y] = y_ldsk;
1314 y_ldsk->prev = o_ldsk;
1315 }
1316
1317 /* Generate a trajectory */
1318 PHYREX_One_New_Traj_Given_Disk(y_ldsk,o_ldsk,tree);
1319
1320 return(n_new_disk);
1321 }
1322
1323 /*////////////////////////////////////////////////////////////
1324 ////////////////////////////////////////////////////////////*/
1325
1326 /* Generate a new trajectory, including disk event centers, between
1327 'y_ldsk' a ``young'' lindisk event and 'o_ldsk' an old one. 'y_ldsk
1328 and 'o_ldsk' remain unaffected. Disk events between these two ldsk
1329 should already be set.
1330 */
PHYREX_One_New_Traj_Given_Disk(t_ldsk * y_ldsk,t_ldsk * o_ldsk,t_tree * tree)1331 void PHYREX_One_New_Traj_Given_Disk(t_ldsk *y_ldsk, t_ldsk *o_ldsk, t_tree *tree)
1332 {
1333 int n_disk_btw;
1334 t_ldsk *ldsk;
1335 phydbl min, max;
1336 phydbl *min_disk_coord, *max_disk_coord;
1337 int i,k;
1338 phydbl rad;
1339
1340
1341 /* Number of disks between y_ldsk and o_ldsk */
1342 ldsk = y_ldsk;
1343 n_disk_btw = 0;
1344 while(ldsk->prev != o_ldsk)
1345 {
1346 n_disk_btw++;
1347 ldsk = ldsk->prev;
1348 }
1349
1350 if(n_disk_btw == 0) return;
1351
1352 ldsk = y_ldsk;
1353 rad = ldsk->disk->mmod->rad;
1354 k = 0;
1355 while(ldsk != o_ldsk)
1356 {
1357 if(!ldsk->disk->next) /* Don't change location at tip node */
1358 {
1359 ldsk = ldsk->prev;
1360 continue;
1361 }
1362
1363 PHYREX_Store_Geo_Coord(ldsk->coord);
1364 PHYREX_Store_Geo_Coord(ldsk->disk->centr);
1365
1366 for(i=0;i<tree->mmod->n_dim;i++)
1367 {
1368 min =
1369 MAX(ldsk->disk->mmod->lim_do->lonlat[i],
1370 MAX(ldsk->coord->lonlat[i] - 2.*rad,
1371 o_ldsk->disk->centr->lonlat[i] - rad*(2.*(n_disk_btw-k)-1.)));
1372
1373 max =
1374 MIN(ldsk->disk->mmod->lim_up->lonlat[i],
1375 MIN(ldsk->coord->lonlat[i] + 2.*rad,
1376 o_ldsk->disk->centr->lonlat[i] + rad*(2.*(n_disk_btw-k)-1.)));
1377
1378 assert(!(max < min));
1379
1380 /* New coordinates for the lindisk */
1381 ldsk->coord->lonlat[i] = Uni()*(max - min) + min;
1382 }
1383
1384 /* New coordinate for the centre of the corresponding disk event */
1385 PHYREX_Get_Min_Max_Disk_Given_Ldsk(ldsk->disk,&min_disk_coord,&max_disk_coord,tree);
1386 for(i=0;i<tree->mmod->n_dim;i++) ldsk->disk->centr->lonlat[i] = Uni()*(max_disk_coord[i] - min_disk_coord[i]) + min_disk_coord[i];
1387 Free(min_disk_coord);
1388 Free(max_disk_coord);
1389
1390 ldsk = ldsk->prev;
1391 k++;
1392 }
1393 }
1394
1395 /*////////////////////////////////////////////////////////////
1396 ////////////////////////////////////////////////////////////*/
1397
PHYREX_Uniform_Path_Density(t_ldsk * y_ldsk,t_ldsk * o_ldsk,t_tree * tree)1398 phydbl PHYREX_Uniform_Path_Density(t_ldsk *y_ldsk, t_ldsk *o_ldsk, t_tree *tree)
1399 {
1400 int n_disk_btw;
1401 t_ldsk *ldsk;
1402 phydbl min, max;
1403 int i,k;
1404 phydbl rad;
1405 phydbl log_dens;
1406 phydbl *min_disk_coord, *max_disk_coord;
1407
1408 if(y_ldsk == o_ldsk) return .0;
1409
1410 /* Number of disks between y_ldsk and o_ldsk */
1411 ldsk = y_ldsk;
1412 n_disk_btw = 0;
1413 while(ldsk->prev != o_ldsk)
1414 {
1415 n_disk_btw++;
1416 ldsk = ldsk->prev;
1417 }
1418
1419 if(n_disk_btw == 0) return .0;
1420
1421 log_dens = 0.0;
1422 ldsk = y_ldsk;
1423 rad = ldsk->disk->mmod->rad;
1424 k = 0;
1425 while(ldsk != o_ldsk)
1426 {
1427 if(!ldsk->disk->next)
1428 {
1429 ldsk = ldsk->prev;
1430 continue;
1431 }
1432
1433 for(i=0;i<ldsk->disk->mmod->n_dim;i++)
1434 {
1435 min =
1436 MAX(ldsk->disk->mmod->lim_do->lonlat[i],
1437 MAX(ldsk->coord->lonlat[i] - 2.*rad,
1438 o_ldsk->disk->centr->lonlat[i] - rad*(2.*(n_disk_btw-k)-1.)));
1439
1440 max =
1441 MIN(ldsk->disk->mmod->lim_up->lonlat[i],
1442 MIN(ldsk->coord->lonlat[i] + 2.*rad,
1443 o_ldsk->disk->centr->lonlat[i] + rad*(2.*(n_disk_btw-k)-1.)));
1444
1445 assert(!(max < min));
1446
1447 log_dens -= log(max - min);
1448 }
1449
1450 PHYREX_Get_Min_Max_Disk_Given_Ldsk(ldsk->disk,&min_disk_coord,&max_disk_coord,tree);
1451 for(i=0;i<tree->mmod->n_dim;i++) log_dens -= log(max_disk_coord[i] - min_disk_coord[i]);
1452 Free(min_disk_coord);
1453 Free(max_disk_coord);
1454
1455 ldsk = ldsk->prev;
1456 k++;
1457 }
1458
1459 return(log_dens);
1460 }
1461
1462 /*////////////////////////////////////////////////////////////
1463 ////////////////////////////////////////////////////////////*/
1464
1465 /* Return the index of the 'next' element of 'old' that should be
1466 used in order to reach 'young'.
1467 */
PHYREX_Get_Next_Direction(t_ldsk * young,t_ldsk * old)1468 int PHYREX_Get_Next_Direction(t_ldsk *young, t_ldsk *old)
1469 {
1470 if(young->disk->time < old->disk->time)
1471 {
1472 PhyML_Printf("\n. young (%s) @ time %f; old (%s) @ time %f",
1473 young->coord->id,young->disk->time,
1474 old->coord->id,old->disk->time);
1475 assert(FALSE);
1476 }
1477
1478 assert(!(young == NULL));
1479
1480 if(young->prev == old)
1481 {
1482 int i;
1483 for(i=0;i<old->n_next;i++) if(old->next[i] == young) return i;
1484 }
1485 else
1486 {
1487 return PHYREX_Get_Next_Direction(young->prev,old);
1488 }
1489 return(-1);
1490 }
1491
1492 /*////////////////////////////////////////////////////////////
1493 ////////////////////////////////////////////////////////////*/
1494
PHYREX_Update_Lindisk_List_Range(t_dsk * young,t_dsk * old,t_tree * tree)1495 void PHYREX_Update_Lindisk_List_Range(t_dsk *young, t_dsk *old, t_tree *tree)
1496 {
1497 t_dsk *disk;
1498
1499 disk = young;
1500 do
1501 {
1502 assert(disk);
1503 PHYREX_Update_Lindisk_List_Core(disk,tree);
1504 if(disk == old) break;
1505 disk = disk->prev;
1506 }
1507 while(1);
1508 }
1509
1510 /*////////////////////////////////////////////////////////////
1511 ////////////////////////////////////////////////////////////*/
1512
PHYREX_Update_Lindisk_List(t_tree * tree)1513 void PHYREX_Update_Lindisk_List(t_tree *tree)
1514 {
1515 PHYREX_Update_Lindisk_List_Pre(tree->young_disk->prev,tree);
1516 }
1517
1518 /*////////////////////////////////////////////////////////////
1519 ////////////////////////////////////////////////////////////*/
1520
PHYREX_Update_Lindisk_List_Pre(t_dsk * disk,t_tree * tree)1521 void PHYREX_Update_Lindisk_List_Pre(t_dsk *disk, t_tree *tree)
1522 {
1523 if(!disk) return;
1524 else
1525 {
1526 PHYREX_Update_Lindisk_List_Core(disk,tree);
1527 PHYREX_Update_Lindisk_List_Pre(disk->prev,tree);
1528 }
1529 }
1530
1531 /*////////////////////////////////////////////////////////////
1532 ////////////////////////////////////////////////////////////*/
1533
PHYREX_Update_Lindisk_List_Core(t_dsk * disk,t_tree * tree)1534 void PHYREX_Update_Lindisk_List_Core(t_dsk *disk, t_tree *tree)
1535 {
1536 int i;
1537
1538 if(!disk) return;
1539 if(!disk->next) return;
1540
1541 assert(disk->ldsk_a);
1542
1543 // Set ldsk_a[i] to NULL if it does not point to a tip node
1544 for(i=0;i<tree->n_otu;++i)
1545 if(disk->ldsk_a[i] &&
1546 !(disk->ldsk_a[i]->nd != NULL &&
1547 disk->ldsk_a[i]->nd->tax == YES &&
1548 disk->age_fixed == YES &&
1549 disk->ldsk_a[i]->disk == disk))
1550 disk->ldsk_a[i] = NULL;
1551
1552 disk->n_ldsk_a = 0;
1553 for(i=0;i<tree->n_otu;++i)
1554 if(disk->ldsk_a[i] != NULL)
1555 disk->n_ldsk_a++;
1556
1557 // Make sure the tip nodes are all at the top of ldsk_a
1558 for(i=0;i<disk->n_ldsk_a;++i) assert(disk->ldsk_a[i] != NULL);
1559 for(i=disk->n_ldsk_a;i<tree->n_otu;++i) assert(disk->ldsk_a[i] == NULL);
1560
1561 for(i=0;i<disk->next->n_ldsk_a;++i)
1562 {
1563 // disk->next->ldsk_a[i] does not coalesce or jump on disk->next
1564 // --> add it disk->ldsk_a array
1565 if((disk->next->ldsk_a[i]->prev != NULL &&
1566 disk->next->ldsk_a[i]->prev != disk->next->ldsk) ||
1567 (disk->next->ldsk_a[i]->prev == NULL))
1568 {
1569 disk->ldsk_a[disk->n_ldsk_a] = disk->next->ldsk_a[i];
1570 disk->n_ldsk_a++;
1571 }
1572 }
1573
1574 // A jump or coalescence has occurred on disk->next
1575 // --> add the lineage to disk->ldsk_a array
1576 if(disk->next->ldsk != NULL)
1577 {
1578 disk->ldsk_a[disk->n_ldsk_a] = disk->next->ldsk;
1579 disk->n_ldsk_a++;
1580 }
1581
1582 /* PhyML_Printf("\n"); */
1583 /* for(i=0;i<disk->n_ldsk_a;++i) PhyML_Printf("\n. Disk %s [%12f] ldsk_a: %s (%3d)",disk->id,disk->time,disk->ldsk_a[i]->coord->id,disk->ldsk_a[i]->nd ? disk->ldsk_a[i]->nd->tax : -1); */
1584 /* if(disk->ldsk != NULL) PhyML_Printf("\n* Has %s on it (next: %s @ %12f %s)", */
1585 /* disk->ldsk->coord->id, */
1586 /* disk->ldsk->next[0]->coord->id, */
1587 /* disk->ldsk->next[0]->disk->time, */
1588 /* disk->ldsk->next[0]->prev->coord->id); */
1589 /* PhyML_Printf("\n. n_ldsk_a: %d",disk->n_ldsk_a); */
1590
1591 if(disk->n_ldsk_a == 0 || disk->n_ldsk_a > tree->n_otu)
1592 {
1593 PhyML_Fprintf(stderr,"\n. disk: %s (%p,%d) time: %f next: %s (%f,%d,%c,%d) prev: %s (%f,%d) disk->n_ldsk_a: %d coord: %s n_otu: %d",
1594 disk->id,
1595 disk,
1596 disk->age_fixed,
1597 disk->time,
1598 disk->next ? disk->next->id : "??",
1599 disk->next ? disk->next->time : 0.0,
1600 disk->next ? disk->next->n_ldsk_a : -1,
1601 disk->next->ldsk ? 'y' : 'n',
1602 disk->next->ldsk ? disk->next->ldsk->n_next : -1,
1603 disk->prev ? disk->prev->id : "??",
1604 disk->prev ? disk->prev->time : 0.0,
1605 disk->prev ? disk->prev->n_ldsk_a : -1,
1606 disk->n_ldsk_a,
1607 disk->ldsk?disk->ldsk->coord->id:"??",
1608 tree->n_otu);
1609 assert(FALSE);
1610 }
1611 }
1612
1613 /*////////////////////////////////////////////////////////////
1614 ////////////////////////////////////////////////////////////*/
1615
1616 /* Connect all the ldsk between y_ldsk (young ldsk) and o_ldsk (old ldsk).
1617 The disk between y_ldsk and o_ldsk should have all been set already
1618 Note: the disks in **disk are sorted in ascending order of their
1619 times
1620 */
1621
PHYREX_Connect_Ldsk_Given_Disk(t_dsk ** disk,int n_disk,t_ldsk * y_ldsk,t_ldsk * o_ldsk,int dir_o_y)1622 void PHYREX_Connect_Ldsk_Given_Disk(t_dsk **disk, int n_disk, t_ldsk *y_ldsk, t_ldsk *o_ldsk, int dir_o_y)
1623 {
1624 int i,j;
1625 t_dsk *disk_tmp;
1626
1627 /* Sort these events by ascending order of their times */
1628 for(i=0;i<n_disk-1;i++)
1629 {
1630 for(j=i+1;j<n_disk;j++)
1631 {
1632 if(disk[j]->time > disk[i]->time)
1633 {
1634 disk_tmp = disk[i];
1635 disk[i] = disk[j];
1636 disk[j] = disk_tmp;
1637 }
1638 }
1639 }
1640
1641
1642 for(i=0;i<n_disk;i++)
1643 {
1644 if(!i)
1645 {
1646 disk[i]->ldsk->next[0] = y_ldsk;
1647 y_ldsk->prev = disk[i]->ldsk;
1648 /* printf("\n. connect %s to %s",disk[i]->ldsk->coord->id,y_ldsk->coord->id); fflush(NULL); */
1649 }
1650 else
1651 {
1652 disk[i]->ldsk->next[0] = disk[i-1]->ldsk;
1653 disk[i-1]->ldsk->prev = disk[i]->ldsk;
1654 /* printf("\n. connect %s to %s",disk[i]->ldsk->coord->id,disk[i-1]->ldsk->coord->id); fflush(NULL); */
1655 }
1656 }
1657
1658 /* printf("\n. connect %s next dir: %d [%d] to %s",o_ldsk->coord->id,dir_o_y,o_ldsk->n_next,disk[n_disk-1]->ldsk->coord->id); fflush(NULL); */
1659 o_ldsk->next[dir_o_y] = disk[n_disk-1]->ldsk;
1660 disk[n_disk-1]->ldsk->prev = o_ldsk;
1661
1662 }
1663
1664 /*////////////////////////////////////////////////////////////
1665 ////////////////////////////////////////////////////////////*/
1666
PHYREX_Print_Struct(char sign,t_tree * tree)1667 void PHYREX_Print_Struct(char sign, t_tree *tree)
1668 {
1669 t_dsk *disk;
1670 int i;
1671 t_ldsk *ldisk;
1672
1673 PHYREX_Update_Lindisk_List(tree);
1674
1675 disk = tree->young_disk;
1676 /* while(disk->prev) disk = disk->prev; */
1677 do
1678 {
1679 PhyML_Printf("\n%c Disk: %s @ %7.3f has %3s on it is_coal? %2d rad: %f age fixed? %d #out: %d",
1680 sign,
1681 disk->id,
1682 disk->time,disk->ldsk?disk->ldsk->coord->id:NULL,
1683 disk->ldsk?disk->ldsk->n_next:-1,
1684 tree->mmod->rad,
1685 disk->age_fixed,
1686 PHYREX_Number_Of_Outgoing_Ldsks(disk));
1687 /* for(j=0;j<tree->mmod->n_dim;j++) PhyML_Printf(" %f",disk->centr->lonlat[j]); */
1688 /* fflush(NULL); */
1689
1690
1691 for(i=0;i<disk->n_ldsk_a;i++)
1692 {
1693 ldisk = disk->ldsk_a[i];
1694
1695 PhyML_Printf("\n%c ldisk: %s%c prev: %s",
1696 sign,
1697 ldisk->coord->id,
1698 (ldisk->nd && ldisk->nd->tax == YES && ldisk->disk == disk) ? '*' : ' ',
1699 ldisk->prev ? ldisk->prev->coord->id : NULL);
1700 fflush(NULL);
1701
1702 /* for(j=0;j<tree->mmod->n_dim;j++) */
1703 /* { */
1704 /* PhyML_Printf(" %f",ldisk->coord->lonlat[j]); */
1705 /* fflush(NULL); */
1706
1707 /* if(FABS(ldisk->coord->lonlat[j] - ldisk->disk->centr->lonlat[j]) > 2.*tree->mmod->rad && */
1708 /* ldisk->disk->ldsk == ldisk) PhyML_Printf(" ! "); */
1709
1710 /* if(ldisk->prev) */
1711 /* { */
1712 /* if(ldisk->coord->lonlat[j] < ldisk->prev->disk->centr->lonlat[j] - tree->mmod->rad) PhyML_Printf(" #a "); */
1713 /* if(ldisk->coord->lonlat[j] > ldisk->prev->disk->centr->lonlat[j] + tree->mmod->rad) PhyML_Printf(" #b "); */
1714 /* } */
1715
1716 /* if(ldisk->next) */
1717 /* { */
1718 /* if(ldisk->coord->lonlat[j] < ldisk->disk->centr->lonlat[j] - tree->mmod->rad) PhyML_Printf(" $a "); */
1719 /* if(ldisk->coord->lonlat[j] > ldisk->disk->centr->lonlat[j] + tree->mmod->rad) PhyML_Printf(" $b "); */
1720 /* } */
1721 }
1722
1723
1724 disk = disk->prev;
1725 }
1726 while(disk);
1727 PhyML_Printf("\n. End of struct");
1728 }
1729
1730 /*////////////////////////////////////////////////////////////
1731 ////////////////////////////////////////////////////////////*/
1732
PHYREX_Check_Struct(t_tree * tree)1733 int PHYREX_Check_Struct(t_tree *tree)
1734 {
1735 int i;
1736 t_ldsk *ldisk;
1737 t_dsk *disk;
1738
1739 disk = tree->young_disk;
1740
1741 do
1742 {
1743 if(disk->n_ldsk_a == 0)
1744 {
1745 PhyML_Printf("\n. disk %s time %f has n_ldsk_a=0 (next has %d prev has %d)",
1746 disk->id,
1747 disk->time,
1748 disk->next ? disk->next->n_ldsk_a : -1,
1749 disk->prev ? disk->prev->n_ldsk_a : -1);
1750 assert(FALSE);
1751 return 0;
1752 }
1753
1754 disk = disk->prev;
1755 }
1756 while(disk);
1757
1758 // Check times
1759 for(i=0;i<tree->n_otu;++i)
1760 {
1761 ldisk = tree->a_nodes[i]->ldsk;
1762 assert(ldisk);
1763 do
1764 {
1765 if(ldisk->prev->disk->time > ldisk->disk->time)
1766 {
1767 PhyML_Printf("\n. ldisk->id: %s ldisk->prev->id: %s ldsk->disk->time: %f ldsk->prev->disk->time: %f diff: %g ldisk->prev->disk: %s ldisk->disk: %s",
1768 ldisk->coord->id,
1769 ldisk->prev->coord->id,
1770 ldisk->disk->time,
1771 ldisk->prev->disk->time,
1772 ldisk->prev->disk->time-ldisk->disk->time,
1773 ldisk->prev->disk->id,
1774 ldisk->disk->id);
1775 assert(FALSE);
1776 return 0;
1777 }
1778 ldisk = ldisk->prev;
1779 }
1780 while(ldisk->prev);
1781 }
1782 return 1;
1783 }
1784
1785 /*////////////////////////////////////////////////////////////
1786 ////////////////////////////////////////////////////////////*/
1787
PHYREX_Store_Geo_Coord(t_geo_coord * t)1788 void PHYREX_Store_Geo_Coord(t_geo_coord *t)
1789 {
1790 int i;
1791
1792 for(i=0;i<t->dim;i++) t->cpy->lonlat[i] = t->lonlat[i];
1793 t->cpy->dim = t->dim;
1794 strcpy(t->cpy->id,t->id);
1795 }
1796
1797 /*////////////////////////////////////////////////////////////
1798 ////////////////////////////////////////////////////////////*/
1799
1800
PHYREX_Restore_Geo_Coord(t_geo_coord * t)1801 void PHYREX_Restore_Geo_Coord(t_geo_coord *t)
1802 {
1803 int i;
1804
1805 for(i=0;i<t->dim;i++) t->lonlat[i] = t->cpy->lonlat[i];
1806 t->dim = t->cpy->dim;
1807 strcpy(t->id,t->cpy->id);
1808 }
1809
1810 /*////////////////////////////////////////////////////////////
1811 ////////////////////////////////////////////////////////////*/
1812
PHYREX_Total_Number_Of_Intervals(t_tree * tree)1813 int PHYREX_Total_Number_Of_Intervals(t_tree *tree)
1814 {
1815 t_dsk *disk;
1816 int n_intervals;
1817
1818 assert(!(tree->young_disk->next));
1819
1820 disk = tree->young_disk;
1821 n_intervals = 0;
1822 while(disk->prev)
1823 {
1824 n_intervals++;
1825 disk = disk->prev;
1826 }
1827 return(n_intervals);
1828
1829 }
1830
1831 /*////////////////////////////////////////////////////////////
1832 ////////////////////////////////////////////////////////////*/
1833
PHYREX_Number_Of_Intervals_Range(t_dsk * young,t_dsk * old,t_tree * tree)1834 int PHYREX_Number_Of_Intervals_Range(t_dsk *young, t_dsk *old, t_tree *tree)
1835 {
1836 t_dsk *disk;
1837 int n_intervals;
1838
1839 disk = young;
1840 n_intervals = 0;
1841 do
1842 {
1843 n_intervals++;
1844 disk = disk->prev;
1845 assert(disk);
1846 }
1847 while(disk != old);
1848
1849 return(n_intervals);
1850 }
1851
1852
1853 /*////////////////////////////////////////////////////////////
1854 ////////////////////////////////////////////////////////////*/
1855
PHYREX_Total_Number_Of_Hit_Disks(t_tree * tree)1856 int PHYREX_Total_Number_Of_Hit_Disks(t_tree *tree)
1857 {
1858 t_dsk *disk;
1859 int n_hit_disks;
1860
1861 assert(!(tree->young_disk->next));
1862
1863 disk = tree->young_disk;
1864 n_hit_disks = 0;
1865 while(disk)
1866 {
1867 if(disk->ldsk) n_hit_disks++;
1868 disk = disk->prev;
1869 }
1870 return(n_hit_disks);
1871
1872 }
1873
1874 /*////////////////////////////////////////////////////////////
1875 ////////////////////////////////////////////////////////////*/
1876
PHYREX_Total_Number_Of_Coal_Disks(t_tree * tree)1877 int PHYREX_Total_Number_Of_Coal_Disks(t_tree *tree)
1878 {
1879 t_dsk *disk;
1880 int n_coal_disks;
1881
1882 assert(!(tree->young_disk->next));
1883
1884 disk = tree->young_disk;
1885 n_coal_disks = 0;
1886 while(disk)
1887 {
1888 if(disk->ldsk && disk->ldsk->n_next > 1) n_coal_disks++;
1889 disk = disk->prev;
1890 }
1891
1892 return(n_coal_disks);
1893 }
1894
1895 /*////////////////////////////////////////////////////////////
1896 ////////////////////////////////////////////////////////////*/
1897
PHYREX_Log_Dunif_Rectangle_Overlap(t_ldsk * ldsk,t_dsk * disk,t_phyrex_mod * mmod)1898 phydbl PHYREX_Log_Dunif_Rectangle_Overlap(t_ldsk *ldsk, t_dsk *disk, t_phyrex_mod *mmod)
1899 {
1900 phydbl up, down, left, rght;
1901 phydbl log_dens;
1902
1903 log_dens = 0.0;
1904
1905 up = MIN(disk->centr->lonlat[0] + mmod->rad, mmod->lim_up->lonlat[0]);
1906 down = MAX(disk->centr->lonlat[0] - mmod->rad, mmod->lim_do->lonlat[0]);
1907 rght = MIN(disk->centr->lonlat[1] + mmod->rad, mmod->lim_up->lonlat[1]);
1908 left = MAX(disk->centr->lonlat[1] - mmod->rad, mmod->lim_do->lonlat[1]);
1909
1910
1911 if(ldsk->coord->lonlat[0] < down) return UNLIKELY;
1912 if(ldsk->coord->lonlat[0] > up) return UNLIKELY;
1913 if(ldsk->coord->lonlat[1] < left) return UNLIKELY;
1914 if(ldsk->coord->lonlat[1] > rght) return UNLIKELY;
1915
1916 log_dens = -log(up-down)-log(rght-left);
1917
1918 return(log_dens);
1919 }
1920
1921 /*////////////////////////////////////////////////////////////
1922 ////////////////////////////////////////////////////////////*/
1923 /* Samples uniformly within a rectangle (with truncation for border) */
1924 /* and returns the corresponding log density */
PHYREX_Runif_Rectangle_Overlap(t_ldsk * ldsk,t_dsk * disk,t_phyrex_mod * mmod)1925 phydbl PHYREX_Runif_Rectangle_Overlap(t_ldsk *ldsk, t_dsk *disk, t_phyrex_mod *mmod)
1926 {
1927 phydbl up, down, left, rght;
1928
1929 up = MIN(disk->centr->lonlat[0] + mmod->rad, mmod->lim_up->lonlat[0]);
1930 down = MAX(disk->centr->lonlat[0] - mmod->rad, mmod->lim_do->lonlat[0]);
1931 rght = MIN(disk->centr->lonlat[1] + mmod->rad, mmod->lim_up->lonlat[1]);
1932 left = MAX(disk->centr->lonlat[1] - mmod->rad, mmod->lim_do->lonlat[1]);
1933
1934 ldsk->coord->lonlat[0] = Uni()*(up - down) + down;
1935 ldsk->coord->lonlat[1] = Uni()*(rght - left) + left;
1936
1937 /* printf("\n. disk %s (%f %f) rad: %f up: %f down: %f rght: %f left: %f", */
1938 /* disk->id, */
1939 /* disk->centr->lonlat[0], */
1940 /* disk->centr->lonlat[1], */
1941 /* mmod->rad, */
1942 /* up,down,rght,left); */
1943 return(log(up-down)+log(rght-left));
1944 }
1945
1946 /*////////////////////////////////////////////////////////////
1947 ////////////////////////////////////////////////////////////*/
1948
PHYREX_Rnorm_Trunc(t_ldsk * ldsk,t_dsk * disk,t_phyrex_mod * mmod)1949 phydbl PHYREX_Rnorm_Trunc(t_ldsk *ldsk, t_dsk *disk, t_phyrex_mod *mmod)
1950 {
1951 int i,err;
1952
1953 err = NO;
1954 for(i=0;i<mmod->n_dim;++i)
1955 {
1956 ldsk->coord->lonlat[i] = Rnorm_Trunc(disk->centr->lonlat[i],
1957 mmod->rad,
1958 mmod->lim_do->lonlat[i],
1959 mmod->lim_up->lonlat[i],&err);
1960 assert(err != YES);
1961 }
1962
1963
1964 return(0.0);
1965 }
1966
1967 /*////////////////////////////////////////////////////////////
1968 ////////////////////////////////////////////////////////////*/
1969
PHYREX_Wrap_Prior_Radius(t_edge * e,t_tree * tree,supert_tree * st)1970 phydbl PHYREX_Wrap_Prior_Radius(t_edge *e, t_tree *tree, supert_tree *st)
1971 {
1972 return PHYREX_LnPrior_Radius(tree);
1973 }
1974
1975 /*////////////////////////////////////////////////////////////
1976 ////////////////////////////////////////////////////////////*/
1977
PHYREX_LnPrior_Lbda(t_tree * tree)1978 phydbl PHYREX_LnPrior_Lbda(t_tree *tree)
1979 {
1980 if(tree->mmod->id == RW) return(0.0);
1981 if(tree->mmod->id == RRW) return(0.0);
1982
1983 if(tree->mmod->lbda < tree->mmod->min_lbda) return UNLIKELY;
1984 if(tree->mmod->lbda > tree->mmod->max_lbda) return UNLIKELY;
1985
1986 /* tree->mmod->c_ln_prior_lbda = */
1987 /* log(tree->mmod->prior_param_lbda) - */
1988 /* tree->mmod->prior_param_lbda*tree->mmod->lbda; */
1989
1990 /* tree->mmod->c_ln_prior_lbda -= log(exp(-tree->mmod->prior_param_lbda*tree->mmod->min_lbda)- */
1991 /* exp(-tree->mmod->prior_param_lbda*tree->mmod->max_lbda)); */
1992
1993 tree->mmod->c_ln_prior_lbda = -log(tree->mmod->max_lbda - tree->mmod->min_lbda);;
1994
1995 return(tree->mmod->c_ln_prior_lbda);
1996 }
1997
1998 /*////////////////////////////////////////////////////////////
1999 ////////////////////////////////////////////////////////////*/
2000
PHYREX_LnPrior_Mu(t_tree * tree)2001 phydbl PHYREX_LnPrior_Mu(t_tree *tree)
2002 {
2003 if(tree->mmod->id == RW) return(0.0);
2004 if(tree->mmod->id == RRW) return(0.0);
2005
2006 if(tree->mmod->mu < tree->mmod->min_mu) return UNLIKELY;
2007 if(tree->mmod->mu > tree->mmod->max_mu) return UNLIKELY;
2008
2009 tree->mmod->c_ln_prior_mu = -log(tree->mmod->max_mu - tree->mmod->min_mu);
2010
2011 /* tree->mmod->c_ln_prior_mu = -2.*log(tree->mmod->mu); */
2012
2013 return(tree->mmod->c_ln_prior_mu);
2014 }
2015
2016 /*////////////////////////////////////////////////////////////
2017 ////////////////////////////////////////////////////////////*/
2018
PHYREX_LnPrior_Radius(t_tree * tree)2019 phydbl PHYREX_LnPrior_Radius(t_tree *tree)
2020 {
2021 if(tree->mmod->id == RW) return(0.0);
2022 if(tree->mmod->id == RRW) return(0.0);
2023
2024 if(tree->mmod->rad < tree->mmod->min_rad) return UNLIKELY;
2025 if(tree->mmod->rad > tree->mmod->max_rad) return UNLIKELY;
2026
2027 tree->mmod->c_ln_prior_rad =
2028 log(tree->mmod->prior_param_rad) -
2029 tree->mmod->prior_param_rad*tree->mmod->rad;
2030
2031 tree->mmod->c_ln_prior_rad -= log(exp(-tree->mmod->prior_param_rad*tree->mmod->min_rad)-
2032 exp(-tree->mmod->prior_param_rad*tree->mmod->max_rad));
2033
2034 return(tree->mmod->c_ln_prior_rad);
2035 }
2036
2037 /*////////////////////////////////////////////////////////////
2038 ////////////////////////////////////////////////////////////*/
2039
PHYREX_LnPrior_Sigsq(t_tree * tree)2040 phydbl PHYREX_LnPrior_Sigsq(t_tree *tree)
2041 {
2042 return(0.0);
2043 }
2044
2045 /*////////////////////////////////////////////////////////////
2046 ////////////////////////////////////////////////////////////*/
2047
PHYREX_Initial_Ldsk_Pos(t_tree * tree)2048 void PHYREX_Initial_Ldsk_Pos(t_tree *tree)
2049 {
2050 t_dsk *disk;
2051 int i,j;
2052 phydbl mean;
2053
2054 disk = tree->young_disk->prev;
2055
2056 do
2057 {
2058 if(disk->ldsk)
2059 {
2060 for(i=0;i<tree->mmod->n_dim;i++)
2061 {
2062 mean = 0.0;
2063 for(j=0;j<disk->ldsk->n_next;j++) mean += disk->ldsk->next[j]->coord->lonlat[i];
2064 disk->ldsk->coord->lonlat[i] = mean / (phydbl)disk->ldsk->n_next;
2065 }
2066 }
2067 disk = disk->prev;
2068 }
2069 while(disk);
2070
2071
2072 }
2073
2074 /*////////////////////////////////////////////////////////////
2075 ////////////////////////////////////////////////////////////*/
2076
PHYREX_Min_Radius(t_tree * tree)2077 phydbl PHYREX_Min_Radius(t_tree *tree)
2078 {
2079 phydbl ori_rad, min_rad;
2080
2081 ori_rad = tree->mmod->rad;
2082 tree->mmod->rad = tree->mmod->max_rad;
2083 do
2084 {
2085 PHYREX_Lk(tree);
2086 tree->mmod->rad -= 1.0;
2087 }
2088 while(tree->mmod->c_lnL > UNLIKELY + 0.1);
2089
2090 min_rad = tree->mmod->rad + 2.0;
2091 tree->mmod->rad = ori_rad;
2092 PHYREX_Lk(tree);
2093 return(min_rad);
2094 }
2095
2096 /*////////////////////////////////////////////////////////////
2097 ////////////////////////////////////////////////////////////*/
2098 /* Get the minimum and maximum values a ldsk can take, given the position of the disk centre */
PHYREX_Get_Min_Max_Ldsk_Given_Disk(t_ldsk * ldsk,phydbl ** min,phydbl ** max,t_tree * tree)2099 void PHYREX_Get_Min_Max_Ldsk_Given_Disk(t_ldsk *ldsk, phydbl **min, phydbl **max, t_tree *tree)
2100 {
2101 phydbl *loc_min,*loc_max;
2102 int i;
2103
2104 if(!ldsk->disk->next) return;
2105
2106 loc_min = (phydbl *)mCalloc(tree->mmod->n_dim, sizeof(phydbl));
2107 loc_max = (phydbl *)mCalloc(tree->mmod->n_dim, sizeof(phydbl));
2108
2109 for(i=0;i<tree->mmod->n_dim;i++)
2110 {
2111 loc_min[i] = ldsk->disk->centr->lonlat[i] - tree->mmod->rad;
2112 loc_max[i] = ldsk->disk->centr->lonlat[i] + tree->mmod->rad;
2113
2114 if(ldsk->prev)
2115 {
2116 loc_min[i] = MAX(loc_min[i],ldsk->prev->disk->centr->lonlat[i] - tree->mmod->rad);
2117 loc_max[i] = MIN(loc_max[i],ldsk->prev->disk->centr->lonlat[i] + tree->mmod->rad);
2118 }
2119
2120 loc_min[i] = MAX(tree->mmod->lim_do->lonlat[i],loc_min[i]);
2121 loc_max[i] = MIN(tree->mmod->lim_up->lonlat[i],loc_max[i]);
2122 }
2123
2124 (*min) = loc_min;
2125 (*max) = loc_max;
2126 }
2127
2128 /*////////////////////////////////////////////////////////////
2129 ////////////////////////////////////////////////////////////*/
2130 /* Get the minimum and maximum values a disk centre can take, given the position of the ldsk */
PHYREX_Get_Min_Max_Disk_Given_Ldsk(t_dsk * disk,phydbl ** min,phydbl ** max,t_tree * tree)2131 void PHYREX_Get_Min_Max_Disk_Given_Ldsk(t_dsk *disk, phydbl **min, phydbl **max, t_tree *tree)
2132 {
2133 phydbl *loc_min,*loc_max;
2134 int i,j;
2135 phydbl tmp_min, tmp_max;
2136
2137 if(!disk->next) return;
2138
2139 loc_min = (phydbl *)mCalloc(tree->mmod->n_dim,sizeof(phydbl));
2140 loc_max = (phydbl *)mCalloc(tree->mmod->n_dim,sizeof(phydbl));
2141
2142 if(!disk->ldsk || tree->mmod->id == SLFV_GAUSSIAN)
2143 {
2144 for(i=0;i<tree->mmod->n_dim;i++)
2145 {
2146 loc_min[i] = tree->mmod->lim_do->lonlat[i];
2147 loc_max[i] = tree->mmod->lim_up->lonlat[i];
2148 }
2149 }
2150 else
2151 {
2152 for(i=0;i<tree->mmod->n_dim;i++)
2153 {
2154 tmp_min = +INFINITY;
2155 tmp_max = -INFINITY;
2156 for(j=0;j<disk->ldsk->n_next;j++)
2157 {
2158 if(disk->ldsk->next[j]->coord->lonlat[i] < tmp_min) tmp_min = disk->ldsk->next[j]->coord->lonlat[i];
2159 if(disk->ldsk->next[j]->coord->lonlat[i] > tmp_max) tmp_max = disk->ldsk->next[j]->coord->lonlat[i];
2160 }
2161
2162 if(disk->ldsk->coord->lonlat[i] < tmp_min) tmp_min = disk->ldsk->coord->lonlat[i];
2163 if(disk->ldsk->coord->lonlat[i] > tmp_max) tmp_max = disk->ldsk->coord->lonlat[i];
2164
2165 loc_min[i] = MAX(tree->mmod->lim_do->lonlat[i],
2166 tmp_max - tree->mmod->rad);
2167
2168 loc_max[i] = MIN(tree->mmod->lim_up->lonlat[i],
2169 tmp_min + tree->mmod->rad);
2170
2171 assert(!(loc_max[i] < loc_min[i]));
2172 }
2173 }
2174
2175 (*min) = loc_min;
2176 (*max) = loc_max;
2177 }
2178
2179 /*////////////////////////////////////////////////////////////
2180 ////////////////////////////////////////////////////////////*/
2181
PHYREX_Update_Disk_Ldsk_Subtree(t_ldsk * root_ldsk,t_tree * tree)2182 void PHYREX_Update_Disk_Ldsk_Subtree(t_ldsk *root_ldsk, t_tree *tree)
2183 {
2184 int i;
2185 for(i=0;i<root_ldsk->n_next;i++) PHYREX_Update_Disk_Ldsk_Subtree_Pre(root_ldsk,root_ldsk->next[i],root_ldsk,tree);
2186 }
2187
2188 /*////////////////////////////////////////////////////////////
2189 ////////////////////////////////////////////////////////////*/
2190
PHYREX_Update_Disk_Ldsk_Subtree_Pre(t_ldsk * old_ldsk,t_ldsk * young_ldsk,t_ldsk * root_ldsk,t_tree * tree)2191 void PHYREX_Update_Disk_Ldsk_Subtree_Pre(t_ldsk *old_ldsk, t_ldsk *young_ldsk, t_ldsk *root_ldsk, t_tree *tree)
2192 {
2193 if(!young_ldsk->disk->next)
2194 {
2195 PHYREX_One_New_Traj_Given_Disk(young_ldsk,root_ldsk,tree);
2196 return;
2197 }
2198 else
2199 {
2200 int i;
2201 PHYREX_Update_Disk_Ldsk_Subtree_Pre(young_ldsk,young_ldsk->next[0],root_ldsk,tree);
2202 if(young_ldsk->n_next > 1)
2203 {
2204 for(i=1;i<young_ldsk->n_next;i++) PHYREX_Update_Disk_Ldsk_Subtree_Pre(young_ldsk,young_ldsk->next[i],young_ldsk,tree);
2205 }
2206 }
2207 }
2208
2209 /*////////////////////////////////////////////////////////////
2210 ////////////////////////////////////////////////////////////*/
2211
PHYREX_Restore_Disk_Ldsk_Subtree(t_ldsk * root_ldsk,t_tree * tree)2212 void PHYREX_Restore_Disk_Ldsk_Subtree(t_ldsk *root_ldsk, t_tree *tree)
2213 {
2214 int i;
2215 for(i=0;i<root_ldsk->n_next;i++) PHYREX_Restore_Disk_Ldsk_Subtree_Pre(root_ldsk,root_ldsk->next[i],tree);
2216 }
2217
2218 /*////////////////////////////////////////////////////////////
2219 ////////////////////////////////////////////////////////////*/
2220
PHYREX_Restore_Disk_Ldsk_Subtree_Pre(t_ldsk * old_ldsk,t_ldsk * young_ldsk,t_tree * tree)2221 void PHYREX_Restore_Disk_Ldsk_Subtree_Pre(t_ldsk *old_ldsk, t_ldsk *young_ldsk, t_tree *tree)
2222 {
2223
2224 if(!young_ldsk->disk->next) return;
2225 else
2226 {
2227 int i;
2228
2229 PHYREX_Restore_Geo_Coord(young_ldsk->coord);
2230 PHYREX_Restore_Geo_Coord(young_ldsk->disk->centr);
2231
2232 for(i=0;i<young_ldsk->n_next;i++)
2233 {
2234 PHYREX_Restore_Disk_Ldsk_Subtree_Pre(young_ldsk,young_ldsk->next[i],tree);
2235 }
2236 }
2237 }
2238
2239 /*////////////////////////////////////////////////////////////
2240 ////////////////////////////////////////////////////////////*/
2241
PHYREX_Proposal_Disk_Ldsk_Subtree(t_ldsk * root_ldsk,phydbl * logdens,t_tree * tree)2242 void PHYREX_Proposal_Disk_Ldsk_Subtree(t_ldsk *root_ldsk, phydbl *logdens, t_tree *tree)
2243 {
2244 int i;
2245 (*logdens) = 0.0;
2246 for(i=0;i<root_ldsk->n_next;i++) PHYREX_Proposal_Disk_Ldsk_Subtree_Pre(root_ldsk,root_ldsk->next[i],root_ldsk,logdens,tree);
2247 }
2248
2249 /*////////////////////////////////////////////////////////////
2250 ////////////////////////////////////////////////////////////*/
2251
PHYREX_Proposal_Disk_Ldsk_Subtree_Pre(t_ldsk * old_ldsk,t_ldsk * young_ldsk,t_ldsk * root_ldsk,phydbl * logdens,t_tree * tree)2252 void PHYREX_Proposal_Disk_Ldsk_Subtree_Pre(t_ldsk *old_ldsk, t_ldsk *young_ldsk, t_ldsk *root_ldsk, phydbl *logdens, t_tree *tree)
2253 {
2254
2255 if(!young_ldsk->disk->next)
2256 {
2257 (*logdens) += PHYREX_Uniform_Path_Density(young_ldsk,root_ldsk,tree);
2258 return;
2259 }
2260 else
2261 {
2262 int i;
2263 for(i=0;i<young_ldsk->n_next;i++)
2264 {
2265 PHYREX_Proposal_Disk_Ldsk_Subtree_Pre(young_ldsk,young_ldsk->next[i],root_ldsk,logdens,tree);
2266 }
2267 }
2268 }
2269
2270 /*////////////////////////////////////////////////////////////
2271 ////////////////////////////////////////////////////////////*/
2272 /* Update the tree structure given the whole set of ldsk events */
2273 /* Coalescent events involving multiple lineages are resolved using */
2274 /* very short internal edges. Tip nodes in the tree are always connected */
2275 /* to the corresponding ldsks. */
PHYREX_Ldsk_To_Tree(t_tree * tree)2276 void PHYREX_Ldsk_To_Tree(t_tree *tree)
2277 {
2278 int i,j;
2279 t_dsk *disk;
2280 t_ldsk *root_ldsk;
2281
2282 /* Reset */
2283 for(i=0;i<2*tree->n_otu-1;++i)
2284 {
2285 for(j=0;j<3;++j)
2286 {
2287 tree->a_nodes[i]->v[j] = NULL;
2288 tree->a_nodes[i]->b[j] = NULL;
2289 }
2290 }
2291
2292 // Erase all connections to internal nodes
2293 disk = tree->young_disk;
2294 do
2295 {
2296 if(disk->ldsk)
2297 {
2298 disk->ldsk->nd = NULL;
2299 assert(disk->age_fixed == NO);
2300 }
2301 disk = disk->prev;
2302 }
2303 while(disk->prev);
2304
2305
2306 // Make sure oldest disk has a ldsk on it
2307 assert(disk->ldsk);
2308 root_ldsk = disk->ldsk;
2309
2310 if(tree->n_root == NULL) tree->n_root = tree->a_nodes[2*tree->n_otu-2];
2311 assert(tree->n_root);
2312
2313 i = 2*tree->n_otu-3;
2314 tree->num_curr_branch_available = 0;
2315 PHYREX_Ldsk_To_Tree_Post(tree->n_root,root_ldsk,&i,tree);
2316
2317
2318 for(i=0;i<tree->n_otu;++i) assert(tree->a_nodes[i]->v[0]);
2319
2320 for(i=0;i<3;++i)
2321 if(tree->n_root->v[2]->v[i] == tree->n_root)
2322 {
2323 tree->n_root->v[2]->v[i] = tree->n_root->v[1];
2324 break;
2325 }
2326
2327 for(i=0;i<3;++i)
2328 if(tree->n_root->v[1]->v[i] == tree->n_root)
2329 {
2330 tree->n_root->v[1]->v[i] = tree->n_root->v[2];
2331 break;
2332 }
2333
2334 Connect_Edges_To_Nodes_Serial(tree);
2335
2336
2337 tree->e_root = NULL;
2338 for(i=0;i<2*tree->n_otu-3;++i)
2339 {
2340 if((tree->a_edges[i]->left == tree->n_root->v[1] && tree->a_edges[i]->rght == tree->n_root->v[2]) ||
2341 (tree->a_edges[i]->left == tree->n_root->v[2] && tree->a_edges[i]->rght == tree->n_root->v[1]))
2342 {
2343 tree->e_root = tree->a_edges[i];
2344 break;
2345 }
2346 }
2347 assert(!(tree->e_root == NULL));
2348
2349 tree->n_root->b[1] = tree->a_edges[2*tree->n_otu-3];
2350 tree->n_root->b[2] = tree->a_edges[2*tree->n_otu-2];
2351
2352 tree->n_root->b[1]->left = tree->n_root;
2353 tree->n_root->b[1]->rght = tree->n_root->v[1];
2354
2355 tree->n_root->b[2]->left = tree->n_root;
2356 tree->n_root->b[2]->rght = tree->n_root->v[2];
2357
2358 Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
2359 Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
2360
2361 MIXT_Propagate_Tree_Update(tree);
2362 }
2363
2364 /*////////////////////////////////////////////////////////////
2365 ////////////////////////////////////////////////////////////*/
2366
PHYREX_Ldsk_To_Tree_Post(t_node * a,t_ldsk * ldsk,int * available,t_tree * tree)2367 void PHYREX_Ldsk_To_Tree_Post(t_node *a, t_ldsk *ldsk, int *available, t_tree *tree)
2368 {
2369 assert(ldsk);
2370 assert(a);
2371
2372 ldsk->nd = a;
2373 a->ldsk = ldsk;
2374 tree->times->nd_t[a->num] = ldsk->disk->time;
2375
2376 if(!ldsk->next) return; // Tip node
2377 else
2378 {
2379 t_node *parent,*son;
2380 int idx_next;
2381 t_ldsk *t;
2382
2383 parent = a;
2384 parent->v[1] = NULL;
2385 parent->v[2] = NULL;
2386 idx_next = 0;
2387 do
2388 {
2389 t = ldsk->next[idx_next];
2390
2391 // Descend along that lineage as long as
2392 // one has not reached a tip (t->next == NULL)
2393 // or a coalescent event (t->n_next > 1)
2394 while(t->next && t->n_next <= 1) t = t->next[0];
2395
2396
2397 if(t->nd == NULL) // t->disk is not a sampled disk
2398 {
2399 assert(t->disk->age_fixed == NO);
2400 son = tree->a_nodes[*available];
2401 (*available) = (*available)-1;
2402 }
2403 else
2404 {
2405 assert(t->disk->age_fixed == YES);
2406 son = t->nd;
2407 }
2408
2409 PHYREX_Ldsk_To_Tree_Post(son,t,available,tree);
2410
2411 // Resolve multifurcation
2412 if(parent->v[2] != NULL && idx_next >= 2)
2413 {
2414 t_node *new_parent;
2415
2416 new_parent = tree->a_nodes[*available];
2417 new_parent->ldsk = ldsk;
2418
2419 (*available) = (*available)-1;
2420
2421 new_parent->v[0] = parent;
2422 new_parent->v[1] = parent->v[2];
2423 new_parent->v[2] = son;
2424
2425 parent->v[2] = new_parent;
2426 son->v[0] = new_parent;
2427 new_parent->v[1]->v[0] = new_parent;
2428
2429 /* PhyML_Printf("\n[] connect %d to %d",parent->num,new_parent->num); */
2430 /* PhyML_Printf("\n[] connect %d to %d",new_parent->num,new_parent->v[1]->num); */
2431 /* PhyML_Printf("\n[] connect %d to %d",new_parent->num,new_parent->v[2]->num); */
2432
2433 tree->times->nd_t[new_parent->num] = ldsk->disk->time;
2434
2435 parent = new_parent;
2436 }
2437 else
2438 {
2439 son->v[0] = parent;
2440 if(!parent->v[1]) parent->v[1] = son;
2441 else parent->v[2] = son;
2442
2443 /* printf("\n[] connect %d to %d",parent->num,son->num); */
2444 }
2445
2446 idx_next++;
2447 }
2448 while(idx_next != ldsk->n_next);
2449 }
2450 }
2451
2452 /*////////////////////////////////////////////////////////////
2453 ////////////////////////////////////////////////////////////*/
2454 // Make sure PHYREX_Make_And_Connect_Tip_Disks was called beforehand
PHYREX_Tree_To_Ldsk(t_tree * tree)2455 void PHYREX_Tree_To_Ldsk(t_tree *tree)
2456 {
2457 t_dsk *a_disk,*disk;
2458 t_node *n;
2459
2460 assert(tree->n_root);
2461 assert(tree->young_disk);
2462
2463 // Initialise root disk
2464 a_disk = PHYREX_Make_Disk_Event(tree->mmod->n_dim,tree->n_otu);
2465 PHYREX_Init_Disk_Event(a_disk,tree->mmod->n_dim,NULL);
2466
2467 a_disk->prev = NULL; // last (i.e., oldest) disk
2468
2469 a_disk->ldsk = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
2470 PHYREX_Init_Lindisk_Node(a_disk->ldsk,a_disk,tree->mmod->n_dim);
2471
2472 tree->n_root->ldsk = a_disk->ldsk;
2473 a_disk->ldsk->nd = tree->n_root;
2474
2475 // Initialize centre of event on the root disk
2476 a_disk->centr->lonlat[0] = Uni()*(tree->mmod->lim_up->lonlat[0]-tree->mmod->lim_do->lonlat[0])+tree->mmod->lim_do->lonlat[0];
2477 a_disk->centr->lonlat[1] = Uni()*(tree->mmod->lim_up->lonlat[1]-tree->mmod->lim_do->lonlat[1])+tree->mmod->lim_do->lonlat[1];
2478
2479 /* Its location */
2480 switch(tree->mmod->id)
2481 {
2482 case SLFV_UNIFORM:
2483 {
2484 PHYREX_Runif_Rectangle_Overlap(a_disk->ldsk,a_disk,tree->mmod);
2485 break;
2486 }
2487 case SLFV_GAUSSIAN:
2488 {
2489 PHYREX_Rnorm_Trunc(a_disk->ldsk,a_disk,tree->mmod);
2490 break;
2491 }
2492 }
2493
2494 a_disk->ldsk->nd = tree->n_root;
2495 a_disk->time = tree->times->nd_t[tree->n_root->num];
2496
2497
2498 /* Inflate_Times_To_Get_Reasonnable_Edge_Lengths(1.E-3,tree); */
2499 Get_Node_Ranks_From_Times(tree);
2500
2501 PHYREX_Tree_To_Ldsk_Post(tree->n_root,tree->n_root->v[1],a_disk,tree);
2502 PHYREX_Tree_To_Ldsk_Post(tree->n_root,tree->n_root->v[2],a_disk,tree);
2503
2504 // Create a doubly-chained list of disks
2505 disk = a_disk;
2506 disk->time = tree->times->nd_t[tree->n_root->num];
2507 n = tree->n_root;
2508 while(n->rk_next)
2509 {
2510 // Only jump to next disk if n and n->rk_next are on distinct disks
2511 // If n->ldsk == n->rk_next->ldsk then n->ldsk->disk and n->rk_next->ldsk->disk
2512 // are the same disk (multiple merger)
2513 if((n->ldsk->disk != n->rk_next->ldsk->disk) && (n->ldsk != n->rk_next->ldsk))
2514 {
2515 disk->next = n->rk_next->ldsk->disk;
2516 disk->next->prev = disk;
2517 disk->next->time = tree->times->nd_t[n->rk_next->num];
2518 disk = disk->next;
2519 }
2520
2521 n = n->rk_next;
2522 }
2523
2524
2525 // Fill in ldsk_a arrays throughout the tree
2526 PHYREX_Update_Lindisk_List(tree);
2527
2528 assert(tree->young_disk->prev);
2529
2530 }
2531
2532 /*////////////////////////////////////////////////////////////
2533 ////////////////////////////////////////////////////////////*/
2534
PHYREX_Tree_To_Ldsk_Post(t_node * a,t_node * d,t_dsk * a_disk,t_tree * tree)2535 void PHYREX_Tree_To_Ldsk_Post(t_node *a, t_node *d, t_dsk *a_disk, t_tree *tree)
2536 {
2537 int i;
2538
2539 assert(a);
2540 assert(d);
2541 assert(a_disk);
2542
2543
2544 if(d->tax)
2545 {
2546 assert(d->ldsk);
2547 PHYREX_Make_Lindisk_Next(a_disk->ldsk);
2548 d->ldsk->prev = a_disk->ldsk;
2549 a_disk->ldsk->next[a_disk->ldsk->n_next-1] = d->ldsk;
2550 a_disk->ldsk->next[a_disk->ldsk->n_next-1]->nd = d;
2551 return;
2552 }
2553 else
2554 {
2555 if(tree->times->nd_t[d->num] > tree->times->nd_t[a->num])
2556 {
2557 t_dsk *d_disk;
2558
2559 PHYREX_Make_Lindisk_Next(a_disk->ldsk);
2560
2561 // Make and initialize descendent disk
2562 d_disk = PHYREX_Make_Disk_Event(tree->mmod->n_dim,tree->n_otu);
2563 assert(d_disk);
2564 PHYREX_Init_Disk_Event(d_disk,tree->mmod->n_dim,NULL);
2565
2566 d_disk->time = tree->times->nd_t[d->num];
2567
2568 d_disk->ldsk = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
2569 PHYREX_Init_Lindisk_Node(d_disk->ldsk,d_disk,tree->mmod->n_dim);
2570
2571 // Initialize centre of event on the root disk
2572 d_disk->centr->lonlat[0] = Uni()*(tree->mmod->lim_up->lonlat[0]-tree->mmod->lim_do->lonlat[0])+tree->mmod->lim_do->lonlat[0];
2573 d_disk->centr->lonlat[1] = Uni()*(tree->mmod->lim_up->lonlat[1]-tree->mmod->lim_do->lonlat[1])+tree->mmod->lim_do->lonlat[1];
2574
2575 /* Its location */
2576 switch(tree->mmod->id)
2577 {
2578 case SLFV_UNIFORM:
2579 {
2580 PHYREX_Runif_Rectangle_Overlap(d_disk->ldsk,d_disk,tree->mmod);
2581 break;
2582 }
2583 case SLFV_GAUSSIAN:
2584 {
2585 PHYREX_Rnorm_Trunc(d_disk->ldsk,d_disk,tree->mmod);
2586 break;
2587 }
2588 }
2589
2590 d_disk->ldsk->nd = d;
2591 d->ldsk = d_disk->ldsk;
2592
2593 a_disk->ldsk->next[a_disk->ldsk->n_next-1] = d_disk->ldsk;
2594 d_disk->ldsk->prev = a_disk->ldsk;
2595
2596 for(i=0;i<3;++i)
2597 {
2598 if(d->v[i] != a && d->b[i] != tree->e_root)
2599 {
2600 PHYREX_Tree_To_Ldsk_Post(d,d->v[i],d_disk,tree);
2601 }
2602 }
2603
2604 }
2605 else // time[d] == time[a] -> add lineage to a_disk instead of creating d_disk
2606 {
2607 d->ldsk = a_disk->ldsk;
2608
2609 for(i=0;i<3;++i)
2610 {
2611 if(d->v[i] != a && d->b[i] != tree->e_root)
2612 {
2613 PHYREX_Tree_To_Ldsk_Post(d,d->v[i],a_disk,tree);
2614 }
2615 }
2616 }
2617 }
2618 }
2619
2620 /*////////////////////////////////////////////////////////////
2621 ////////////////////////////////////////////////////////////*/
2622
PHYREX_Simulate_Disk_And_Node_Times(t_tree * tree)2623 void PHYREX_Simulate_Disk_And_Node_Times(t_tree *tree)
2624 {
2625 t_dsk *disk;
2626
2627 disk = tree->young_disk;
2628 assert(disk->age_fixed == YES); // age of youngest disk should be fixed
2629
2630 do
2631 {
2632 // disk->prev is not a sample --> simulate its age by sampling in exp distribution
2633 if(disk->prev->age_fixed == NO)
2634 {
2635 disk->prev->time = disk->time - Rexp(tree->mmod->lbda);
2636
2637 // set time of internal node sitting on disk->prev
2638 if(disk->prev->ldsk != NULL)
2639 tree->times->nd_t[disk->prev->ldsk->nd->num] =
2640 disk->prev->time;
2641 }
2642 PhyML_Printf("\n. Simulate times disk %s time: %f",disk->id,disk->time);
2643 disk = disk->prev;
2644 }
2645 while(disk->prev);
2646 }
2647
2648 /*////////////////////////////////////////////////////////////
2649 ////////////////////////////////////////////////////////////*/
2650 // Connections between tip nodes and corresponding ldsks (plus the
2651 // associated disks) should be made early on and never modified
2652 // after that.
PHYREX_Make_And_Connect_Tip_Disks(t_tree * tree)2653 void PHYREX_Make_And_Connect_Tip_Disks(t_tree *tree)
2654 {
2655 t_dsk *disk,*new_disk;
2656 t_node *n;
2657
2658 Get_Node_Ranks_From_Tip_Times(tree);
2659
2660 // Find most recent tip
2661 n = tree->a_nodes[0];
2662 while(n->rk_next) n = n->rk_next;
2663
2664 // Create most recent disk
2665 disk = PHYREX_Make_Disk_Event(tree->mmod->n_dim,tree->n_otu);
2666 PHYREX_Init_Disk_Event(disk,tree->mmod->n_dim,NULL);
2667 disk->age_fixed = YES; // Sample disks have their ages fixed
2668 disk->time = tree->times->nd_t[n->num]; // Set time of youngest disk
2669 PhyML_Fprintf(stdout,"\n. Youngest sampled disk time set to %f (disk id: %s)",disk->time,disk->id);
2670
2671 // ldsk_a[0] is connected to youngest tip
2672 disk->ldsk_a[0] = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
2673 PHYREX_Init_Lindisk_Node(disk->ldsk_a[0],disk,tree->mmod->n_dim);
2674 disk->n_ldsk_a = 1;
2675 disk->ldsk_a[0]->nd = n;
2676 n->ldsk = disk->ldsk_a[0];
2677 n->ldsk->disk = disk;
2678
2679
2680 // Set pointer to young_disk here and not elsewhere!
2681 tree->young_disk = disk;
2682 new_disk = NULL;
2683 do
2684 {
2685 assert(n->rk_prev);
2686
2687 // n and n->rk_prev have distinct time stamps -> they should be on two distinct disks
2688 if(Are_Equal(tree->times->nd_t[n->num],tree->times->nd_t[n->rk_prev->num],SMALL) == NO)
2689 {
2690 new_disk = PHYREX_Make_Disk_Event(tree->mmod->n_dim,tree->n_otu);
2691 PHYREX_Init_Disk_Event(new_disk,tree->mmod->n_dim,NULL);
2692 new_disk->age_fixed = YES;
2693
2694 new_disk->ldsk_a[0] = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
2695 PHYREX_Init_Lindisk_Node(new_disk->ldsk_a[0],new_disk,tree->mmod->n_dim);
2696 new_disk->ldsk_a[0]->nd = n->rk_prev;
2697 new_disk->ldsk_a[0]->disk = new_disk;
2698 new_disk->n_ldsk_a = 1;
2699
2700 new_disk->next = disk;
2701 disk->prev = new_disk;
2702
2703 disk = new_disk;
2704 }
2705 // n and n->rk_prev have the same time stamp -> they sit on the same disk
2706 else
2707 {
2708 disk->ldsk_a[disk->n_ldsk_a] = PHYREX_Make_Lindisk_Node(tree->mmod->n_dim);
2709 PHYREX_Init_Lindisk_Node(disk->ldsk_a[disk->n_ldsk_a],disk,tree->mmod->n_dim);
2710 disk->ldsk_a[disk->n_ldsk_a]->nd = n->rk_prev;
2711 disk->ldsk_a[disk->n_ldsk_a]->disk = disk;
2712 disk->n_ldsk_a++;
2713 }
2714
2715 // Set sampled disk time
2716 disk->time = tree->times->nd_t[n->rk_prev->num];
2717 PhyML_Fprintf(stdout,"\n. Set sampled disk (id: %s) time to %15f",disk->id,disk->time);
2718
2719 n->rk_prev->ldsk = disk->ldsk_a[disk->n_ldsk_a-1];
2720 n = n->rk_prev;
2721 }
2722 while(n->rk_prev);
2723 }
2724
2725 /*////////////////////////////////////////////////////////////
2726 ////////////////////////////////////////////////////////////*/
2727
PHYREX_Remove_Lindisk_Next(t_ldsk * ldsk,t_ldsk * rm)2728 void PHYREX_Remove_Lindisk_Next(t_ldsk *ldsk, t_ldsk *rm)
2729 {
2730 t_ldsk **new_next;
2731 int i,pos;
2732
2733 new_next = (t_ldsk **)mCalloc(ldsk->n_next-1+NEXT_BLOCK_SIZE,sizeof(t_ldsk *));
2734
2735 pos = 0;
2736 for(i=0;i<ldsk->n_next;i++)
2737 {
2738 if(ldsk->next[i] != rm)
2739 {
2740 new_next[pos] = ldsk->next[i];
2741 pos++;
2742 }
2743 }
2744
2745 ldsk->n_next--;
2746 Free(ldsk->next);
2747 ldsk->next = new_next;
2748 /* printf("\n. remove next for ldsk %s n_next set to %d",ldsk->coord->id,ldsk->n_next); */
2749 /* fflush(NULL); */
2750 }
2751
2752 /*////////////////////////////////////////////////////////////
2753 ////////////////////////////////////////////////////////////*/
2754 /* Returns the vector of average pairwise distances between ldsk on each disk */
PHYREX_Mean_Pairwise_Distance_Between_Lineage_Locations(t_tree * tree)2755 phydbl *PHYREX_Mean_Pairwise_Distance_Between_Lineage_Locations(t_tree *tree)
2756 {
2757 phydbl *dist;
2758 int block,n_disks,i,j, k;
2759 t_dsk *disk;
2760
2761 PHYREX_Update_Lindisk_List(tree);
2762
2763 dist = NULL;
2764 block = 100;
2765 disk = tree->young_disk;
2766 n_disks = 0;
2767 do
2768 {
2769 if(!n_disks) dist = (phydbl *)mCalloc(block,sizeof(phydbl));
2770 else if(!(n_disks%block)) dist = (phydbl *)mRealloc(dist,n_disks+block,sizeof(phydbl));
2771
2772 dist[n_disks] = 0.0;
2773 for(i=0;i<disk->n_ldsk_a-1;i++)
2774 {
2775 for(j=i+1;j<disk->n_ldsk_a;j++)
2776 {
2777 for(k=0;k<tree->mmod->n_dim;k++)
2778 {
2779 dist[n_disks] += FABS(disk->ldsk_a[i]->coord->lonlat[k] - disk->ldsk_a[j]->coord->lonlat[k]);
2780 printf("\n * %d %f %f %f",
2781 k,
2782 disk->ldsk_a[i]->coord->lonlat[k],
2783 disk->ldsk_a[j]->coord->lonlat[k],
2784 tree->mmod->lim_up->lonlat[k]);
2785 }
2786 }
2787 }
2788 dist[n_disks] /= (phydbl)(disk->n_ldsk_a * (disk->n_ldsk_a-1) / 2.);
2789
2790 printf("\n %d %f %f",disk->n_ldsk_a,disk->time,dist[n_disks]);
2791
2792 n_disks++;
2793 disk = disk->prev;
2794 }
2795 while(disk->prev);
2796
2797 return(dist);
2798
2799 }
2800
2801 /*////////////////////////////////////////////////////////////
2802 ////////////////////////////////////////////////////////////*/
2803
PHYREX_Random_Select_Time_Between_Jumps(t_tree * tree)2804 phydbl PHYREX_Random_Select_Time_Between_Jumps(t_tree *tree)
2805 {
2806 t_dsk *disk,**valid_disks;
2807 int n_valid_disks,block,select;
2808 phydbl time;
2809
2810 valid_disks = NULL;
2811 disk = NULL;
2812 block = 100;
2813
2814 assert(!(tree->young_disk->next));
2815
2816 disk = tree->young_disk->prev;
2817 n_valid_disks = 0;
2818 do
2819 {
2820 if(disk->ldsk != NULL && disk->prev != NULL)
2821 {
2822 if(!n_valid_disks) valid_disks = (t_dsk **)mCalloc(block,sizeof(t_dsk *));
2823 else if(!(n_valid_disks%block)) valid_disks = (t_dsk **)mRealloc(valid_disks,n_valid_disks+block,sizeof(t_dsk *));
2824 valid_disks[n_valid_disks] = disk;
2825 n_valid_disks++;
2826 }
2827 disk = disk->prev;
2828 }
2829 while(disk->prev);
2830
2831 if(!n_valid_disks) return -1.0;
2832
2833 select = Rand_Int(0,n_valid_disks-1);
2834
2835 time = valid_disks[select]->time - valid_disks[select]->ldsk->prev->disk->time;
2836
2837 Free(valid_disks);
2838
2839 return(time);
2840 }
2841
2842 /*////////////////////////////////////////////////////////////
2843 ////////////////////////////////////////////////////////////*/
2844
PHYREX_Is_In_Ldscape(t_ldsk * ldsk,t_phyrex_mod * mmod)2845 int PHYREX_Is_In_Ldscape(t_ldsk *ldsk, t_phyrex_mod *mmod)
2846 {
2847 int j;
2848 for(j=0;j<mmod->n_dim;j++)
2849 if(ldsk->coord->lonlat[j] > mmod->lim_up->lonlat[j] ||
2850 ldsk->coord->lonlat[j] < mmod->lim_do->lonlat[j]) return NO;
2851 return YES;
2852 }
2853
2854 /*////////////////////////////////////////////////////////////
2855 ////////////////////////////////////////////////////////////*/
2856
PHYREX_Mean_Time_Between_Events(t_tree * tree)2857 phydbl PHYREX_Mean_Time_Between_Events(t_tree *tree)
2858 {
2859 int n_inter;
2860 phydbl T;
2861
2862 n_inter = PHYREX_Total_Number_Of_Intervals(tree);
2863 T = -tree->times->nd_t[tree->n_root->num];
2864 return((phydbl)(T/n_inter));
2865 }
2866
2867 /*////////////////////////////////////////////////////////////
2868 ////////////////////////////////////////////////////////////*/
2869
PHYREX_Rand_Pairs_Coal_Times_Dist(t_tree * tree)2870 void PHYREX_Rand_Pairs_Coal_Times_Dist(t_tree *tree)
2871 {
2872 t_node *anc;
2873 phydbl dist;
2874 int i, j;
2875
2876 i = Rand_Int(0,tree->n_otu-1);
2877 j = Rand_Int(0,tree->n_otu-1);
2878
2879 if(i == j) PhyML_Printf("\nxxWxx 0.0 0.0");
2880 else
2881 {
2882
2883 anc = Find_Lca_Pair_Of_Nodes(tree->a_nodes[i],tree->a_nodes[j],tree);
2884 if(anc == NULL)
2885 {
2886 PhyML_Fprintf(stderr,"\n. %s",Write_Tree(tree));
2887 PhyML_Fprintf(stderr,"\n. %s %s",tree->a_nodes[i]->name,tree->a_nodes[j]->name);
2888 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2889 }
2890
2891 PhyML_Printf("\nxxWxx %12f",tree->times->nd_t[anc->num]);
2892 dist = Euclidean_Dist(tree->a_nodes[i]->ldsk->coord,tree->a_nodes[j]->ldsk->coord);
2893 PhyML_Printf(" %f",dist);
2894 }
2895 }
2896
2897 /*////////////////////////////////////////////////////////////
2898 ////////////////////////////////////////////////////////////*/
2899
PHYREX_Read_Tip_Coordinates(t_tree * tree)2900 void PHYREX_Read_Tip_Coordinates(t_tree *tree)
2901 {
2902 char *s;
2903 FILE *fp;
2904 int i,*done,found_sw,found_ne;
2905 phydbl sw_lon, sw_lat, ne_lon, ne_lat;
2906 t_node *n;
2907
2908 assert(tree->young_disk);
2909
2910 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2911 fp = tree->io->fp_in_coord;
2912 done = (int *)mCalloc(tree->n_otu,sizeof(int));
2913
2914 found_sw = NO;
2915 found_ne = NO;
2916
2917 for(i=0;i<tree->n_otu;i++) done[i] = NO;
2918
2919 do
2920 {
2921 if(fscanf(fp,"%s",s) == EOF) break;
2922 for(i=0;i<strlen(s);++i) if(s[i] == '#') break; /* skip comment */
2923 if(i != strlen(s)) continue;
2924
2925 for(i=0;i<tree->n_otu;i++) if(strstr(tree->a_nodes[i]->name,s)) break;
2926
2927 if(i != tree->n_otu) /* Found a match */
2928 {
2929 assert(tree->a_nodes[i]->ldsk);
2930 // First column: latitude, second one: longitude
2931 if(fscanf(fp,"%lf",&(tree->a_nodes[i]->ldsk->coord->lonlat[1])) == EOF) break;
2932 if(fscanf(fp,"%lf",&(tree->a_nodes[i]->ldsk->coord->lonlat[0])) == EOF) break;
2933 done[i] = YES;
2934 }
2935 else
2936 {
2937 if(!strcmp(s,"|SouthWest|") || !strcmp(s,"|southwest|") || !strcmp(s,"|Southwest|"))
2938 {
2939 found_sw = YES;
2940 if(fscanf(fp,"%lf",&(sw_lat)) == EOF) break;
2941 if(fscanf(fp,"%lf",&(sw_lon)) == EOF) break;
2942 }
2943 else if(!strcmp(s,"|NorthEast|") || !strcmp(s,"|northeast|") || !strcmp(s,"|Northeast|"))
2944 {
2945 found_ne = YES;
2946 if(fscanf(fp,"%lf",&(ne_lat)) == EOF) break;
2947 if(fscanf(fp,"%lf",&(ne_lon)) == EOF) break;
2948 }
2949 else /* Haven't found any match but still need to skip long and lat for unsampled location */
2950 {
2951 phydbl dum;
2952 if(fscanf(fp,"%lf",&dum) == EOF) break;
2953 if(fscanf(fp,"%lf",&dum) == EOF) break;
2954 }
2955 }
2956 }
2957 while(1);
2958
2959 if(found_ne == NO)
2960 {
2961 PhyML_Fprintf(stderr,"\n. Could not find coordinates for northernmost point.");
2962 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2963 }
2964
2965 if(found_sw == NO)
2966 {
2967 PhyML_Fprintf(stderr,"\n. Could not find coordinates for southernmost point.");
2968 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2969 }
2970
2971 for(i=0;i<tree->n_otu;i++)
2972 if(done[i] == NO)
2973 {
2974 PhyML_Fprintf(stderr,"\n. Could not find coordinates for '%s'.",tree->a_nodes[i]->name);
2975 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2976 }
2977
2978 n = NULL;
2979 for(i=0;i<tree->n_otu;i++)
2980 {
2981 n = tree->a_nodes[i];
2982
2983 PhyML_Printf("\n. Coordinates of '%-50s': %12f\t %12f",
2984 tree->a_nodes[i]->name,
2985 n->ldsk->coord->lonlat[0],
2986 n->ldsk->coord->lonlat[1]);
2987 }
2988
2989 tree->mmod->lim_up->lonlat[0] = ne_lon;
2990 tree->mmod->lim_up->lonlat[1] = ne_lat;
2991
2992 tree->mmod->lim_do->lonlat[0] = sw_lon;
2993 tree->mmod->lim_do->lonlat[1] = sw_lat;
2994
2995
2996 Free(s);
2997 Free(done);
2998
2999 }
3000
3001 /*////////////////////////////////////////////////////////////
3002 ////////////////////////////////////////////////////////////*/
3003
PHYREX_Tree_Height(t_tree * tree)3004 phydbl PHYREX_Tree_Height(t_tree *tree)
3005 {
3006 t_dsk *disk;
3007
3008 disk = tree->young_disk;
3009 while(disk && disk->prev) disk = disk->prev;
3010
3011 return(disk->time);
3012 }
3013
3014 /*////////////////////////////////////////////////////////////
3015 ////////////////////////////////////////////////////////////*/
3016
PHYREX_Random_Insert_Ldsk_In_Next_List(t_ldsk * ins,t_ldsk * where)3017 int PHYREX_Random_Insert_Ldsk_In_Next_List(t_ldsk *ins, t_ldsk *where)
3018 {
3019 int size, pos, i, *rk;
3020 t_ldsk **next_cpy;
3021
3022 size = where->n_next;
3023
3024 rk = (int *)mCalloc(size,sizeof(int));
3025 next_cpy = (t_ldsk **)mCalloc(size,sizeof(t_ldsk *));
3026
3027 for(i=0;i<size;i++) next_cpy[i] = where->next[i];
3028
3029 pos = Rand_Int(0,size);
3030
3031 for(i=0;i<size;i++)
3032 {
3033 if(i < pos) rk[i] = i;
3034 else rk[i] = i+1;
3035 }
3036
3037 PHYREX_Make_Lindisk_Next(where);
3038
3039 for(i=0;i<size;i++) where->next[rk[i]] = next_cpy[i];
3040
3041 where->next[pos]= ins;
3042
3043
3044 Free(rk);
3045 Free(next_cpy);
3046
3047 return pos;
3048 }
3049
3050 /*////////////////////////////////////////////////////////////
3051 ////////////////////////////////////////////////////////////*/
3052
PHYREX_Insert_Ldsk_In_Next_List(t_ldsk * ins,int pos,t_ldsk * where)3053 void PHYREX_Insert_Ldsk_In_Next_List(t_ldsk *ins, int pos, t_ldsk *where)
3054 {
3055 int size, i, *rk;
3056 t_ldsk **next_cpy;
3057
3058 size = where->n_next;
3059
3060 rk = (int *)mCalloc(size,sizeof(int));
3061 next_cpy = (t_ldsk **)mCalloc(size,sizeof(t_ldsk *));
3062
3063 for(i=0;i<size;i++) next_cpy[i] = where->next[i];
3064
3065 for(i=0;i<size;i++)
3066 {
3067 if(i < pos) rk[i] = i;
3068 else rk[i] = i+1;
3069 }
3070
3071 PHYREX_Make_Lindisk_Next(where);
3072
3073 for(i=0;i<size;i++) where->next[rk[i]] = next_cpy[i];
3074
3075 where->next[pos]= ins;
3076
3077 Free(rk);
3078 Free(next_cpy);
3079 }
3080
3081 /*////////////////////////////////////////////////////////////
3082 ////////////////////////////////////////////////////////////*/
3083 /* Remove path between young ldsk and old ldsk. If young->prev == old, */
3084 /* then young->prev set to NULL and old->next[dir_to_young] set */
3085 /* to NULL as well. */
3086
PHYREX_Remove_Path(t_ldsk * young,t_ldsk * old,int * pos_old,t_tree * tree)3087 t_ldsk *PHYREX_Remove_Path(t_ldsk *young, t_ldsk *old, int *pos_old, t_tree *tree)
3088 {
3089 t_ldsk *path,*ldsk;
3090 int dir_old_young,jumps;
3091
3092 path = NULL;
3093
3094 dir_old_young = PHYREX_Get_Next_Direction(young,old);
3095 assert(dir_old_young >= 0);
3096 *pos_old = dir_old_young;
3097 PHYREX_Remove_Lindisk_Next(old,old->next[dir_old_young]);
3098
3099
3100 jumps = 1;
3101 ldsk = young->prev;
3102 while(ldsk != old)
3103 {
3104 if(jumps == 1) path = ldsk;
3105 PHYREX_Remove_Disk(ldsk->disk);
3106 ldsk = ldsk->prev;
3107 jumps++;
3108 }
3109
3110 if(jumps == 1)
3111 path = NULL;
3112 else
3113 {
3114 /* Set end of path to NULL */
3115 ldsk = path;
3116 while(ldsk->prev != old) { ldsk = ldsk->prev; assert(ldsk); }
3117 ldsk->prev = NULL;
3118 }
3119
3120
3121 /* path == NULL if young->prev = old, otherwise path points to the first ldsk */
3122 /* after jump away from young towards the past (i.e., towards old) */
3123 return(path);
3124 }
3125
3126 /*////////////////////////////////////////////////////////////
3127 ////////////////////////////////////////////////////////////*/
3128
PHYREX_Insert_Path(t_ldsk * young,t_ldsk * old,t_ldsk * path,int pos,t_tree * tree)3129 void PHYREX_Insert_Path(t_ldsk *young, t_ldsk *old, t_ldsk *path, int pos, t_tree *tree)
3130 {
3131 t_ldsk *ldsk;
3132
3133 assert(path != young);
3134
3135 if(path == NULL)
3136 {
3137 young->prev = old;
3138 PHYREX_Insert_Ldsk_In_Next_List(young,pos,old);
3139 }
3140 else
3141 {
3142 /* Attach path to the young ldsk */
3143 path->next[0] = young;
3144 young->prev = path;
3145
3146 ldsk = path;
3147 do
3148 {
3149 PHYREX_Insert_Disk(ldsk->disk,tree);
3150 ldsk = ldsk->prev;
3151 }
3152 while(ldsk);
3153
3154 /* Get to the end of path */
3155 ldsk = path;
3156 while(ldsk->prev != NULL) { ldsk = ldsk->prev; assert(ldsk); }
3157
3158 /* Attach it to old ldsk (both ways)*/
3159 PHYREX_Insert_Ldsk_In_Next_List(ldsk,pos,old);
3160 ldsk->prev = old;
3161 }
3162
3163 }
3164
3165 /*////////////////////////////////////////////////////////////
3166 ////////////////////////////////////////////////////////////*/
3167
PHYREX_Time_Tree_Length(t_tree * tree)3168 phydbl PHYREX_Time_Tree_Length(t_tree *tree)
3169 {
3170 phydbl len;
3171 int i;
3172 t_dsk *disk;
3173
3174 disk = tree->young_disk;
3175 while(disk->prev) disk = disk->prev;
3176
3177 len = 0.0;
3178 for(i=0;i<disk->ldsk->n_next;i++) PHYREX_Time_Tree_Length_Pre(disk->ldsk,disk->ldsk->next[i],&len,tree);
3179
3180 assert(!(isnan(len) || isinf(len)));
3181
3182 return(len);
3183 }
3184
3185 /*////////////////////////////////////////////////////////////
3186 ////////////////////////////////////////////////////////////*/
3187
PHYREX_Time_Tree_Length_Pre(t_ldsk * a,t_ldsk * d,phydbl * len,t_tree * tree)3188 void PHYREX_Time_Tree_Length_Pre(t_ldsk *a, t_ldsk *d, phydbl *len, t_tree *tree)
3189 {
3190 int i;
3191
3192 (*len) += FABS(a->disk->time - d->disk->time);
3193
3194 if(d->disk->next == NULL) return;
3195 else
3196 for(i=0;i<d->n_next;i++) PHYREX_Time_Tree_Length_Pre(d,d->next[i],len,tree);
3197 }
3198
3199 /*////////////////////////////////////////////////////////////
3200 ////////////////////////////////////////////////////////////*/
3201
PHYREX_Is_On_Path(t_ldsk * target,t_ldsk * young,t_ldsk * old)3202 int PHYREX_Is_On_Path(t_ldsk *target, t_ldsk *young, t_ldsk *old)
3203 {
3204 t_ldsk *ldsk;
3205
3206 if(target == young || target == old) return NO;
3207
3208 assert(!(young->disk->time < old->disk->time));
3209
3210 ldsk = young->prev;
3211 while(ldsk != old)
3212 {
3213 if(ldsk == target) return YES;
3214 ldsk = ldsk->prev;
3215 assert(!(ldsk == NULL));
3216 }
3217 return NO;
3218 }
3219
3220 /*////////////////////////////////////////////////////////////
3221 ////////////////////////////////////////////////////////////*/
3222
PHYREX_Path_Len(t_ldsk * young,t_ldsk * old)3223 int PHYREX_Path_Len(t_ldsk *young, t_ldsk *old)
3224 {
3225 t_ldsk *ldsk;
3226 int len;
3227
3228 len = 0;
3229 ldsk = young;
3230 while(ldsk != old)
3231 {
3232 len++;
3233 ldsk = ldsk->prev;
3234 }
3235 len++;
3236 return len;
3237 }
3238
3239 /*////////////////////////////////////////////////////////////
3240 ////////////////////////////////////////////////////////////*/
3241
PHYREX_Print_Disk_Lk(t_tree * tree)3242 void PHYREX_Print_Disk_Lk(t_tree *tree)
3243 {
3244 t_dsk *disk;
3245
3246 PHYREX_Update_Lindisk_List(tree);
3247
3248 disk = tree->young_disk->prev;
3249
3250 do
3251 {
3252 PhyML_Printf("\n. Disk: %s time: %12f lk: %12f cumlk: %12f",
3253 disk->id,
3254 disk->time,
3255 PHYREX_Lk_Core(disk,tree),
3256 -1.);
3257
3258 disk = disk->prev;
3259 }
3260 while(disk->prev);
3261
3262 }
3263
3264 /*////////////////////////////////////////////////////////////
3265 ////////////////////////////////////////////////////////////*/
3266
PHYREX_Find_Lca_Pair_Of_Ldsk(t_ldsk * n1,t_ldsk * n2,t_tree * tree)3267 t_ldsk *PHYREX_Find_Lca_Pair_Of_Ldsk(t_ldsk *n1, t_ldsk *n2, t_tree *tree)
3268 {
3269 t_ldsk **list1, **list2, *lca;
3270 int len1, len2;
3271
3272 assert(n1);
3273 assert(n2);
3274
3275 if(n1 == n2) return(n1);
3276
3277 PHYREX_Get_List_Of_Ancestors(n1,&list1,&len1,tree);
3278 PHYREX_Get_List_Of_Ancestors(n2,&list2,&len2,tree);
3279
3280 len1 = len1-1;
3281 len2 = len2-1;
3282
3283 assert(list1[len1] == list2[len2]);
3284
3285 do
3286 {
3287 if((len1 < 0 || len2 < 0) || (list1[len1] != list2[len2])) break;
3288 len1--;
3289 len2--;
3290 }
3291 while(len1 && len2);
3292
3293
3294 lca = list1[len1+1];
3295
3296 Free(list1);
3297 Free(list2);
3298
3299 assert(lca);
3300
3301 return(lca);
3302
3303 }
3304
3305 /*////////////////////////////////////////////////////////////
3306 ////////////////////////////////////////////////////////////*/
3307
PHYREX_Get_List_Of_Ancestors(t_ldsk * start,t_ldsk *** list,int * len,t_tree * tree)3308 void PHYREX_Get_List_Of_Ancestors(t_ldsk *start, t_ldsk ***list, int *len, t_tree *tree)
3309 {
3310 int block,i;
3311 t_ldsk *ldsk;
3312
3313 assert(start);
3314 block = 100;
3315
3316 *list = (t_ldsk **)mCalloc(block,sizeof(t_ldsk *));
3317
3318 ldsk = start;
3319 i = 0;
3320 do
3321 {
3322 (*list)[i] = ldsk;
3323 if(!(i%block)) *list = (t_ldsk **)mRealloc(*list,i+block,sizeof(t_ldsk *));
3324 i++;
3325 ldsk = ldsk->prev;
3326 }
3327 while(ldsk);
3328
3329 (*len) = i;
3330 }
3331
3332 /*////////////////////////////////////////////////////////////
3333 ////////////////////////////////////////////////////////////*/
3334
PHYREX_Dist_To_Lca(t_ldsk * d,t_ldsk * lca)3335 phydbl PHYREX_Dist_To_Lca(t_ldsk *d, t_ldsk *lca)
3336 {
3337 return(fabs(d->disk->time - lca->disk->time));
3338 }
3339
3340 /*////////////////////////////////////////////////////////////
3341 ////////////////////////////////////////////////////////////*/
3342
PHYREX_Dist_Between_Two_Ldsk(t_ldsk * n1,t_ldsk * n2,t_tree * tree)3343 phydbl PHYREX_Dist_Between_Two_Ldsk(t_ldsk *n1, t_ldsk *n2, t_tree *tree)
3344 {
3345 t_ldsk *lca;
3346
3347 lca = PHYREX_Find_Lca_Pair_Of_Ldsk(n1,n2,tree);
3348
3349 return(PHYREX_Dist_To_Lca(n1,lca)+PHYREX_Dist_To_Lca(n2,lca));
3350 }
3351
3352
3353 /*////////////////////////////////////////////////////////////
3354 ////////////////////////////////////////////////////////////*/
3355
PHYREX_Print_MultiTypeTree_Config_File(int n_sites,char * filename,t_tree * tree)3356 void PHYREX_Print_MultiTypeTree_Config_File(int n_sites, char *filename, t_tree *tree)
3357 {
3358 int i, j, n_demes;
3359 char *s,**deme_names;
3360 FILE *fp;
3361
3362
3363 fp = Openfile(filename,WRITE);
3364 assert(fp);
3365
3366 deme_names = (char **)mCalloc(n_sites,sizeof(char *));
3367
3368 n_demes = 0;
3369 for(i=0;i<tree->n_otu;i++)
3370 {
3371 s = strrchr(tree->a_nodes[i]->ldsk->coord->id,'_');
3372 for(j=0;j<n_demes;j++) if(!strcmp(s+1,deme_names[j])) break;
3373 if(j == n_demes)
3374 {
3375 deme_names[n_demes] = (char *)mCalloc(strlen(s+1)+1,sizeof(char));
3376 strcpy(deme_names[n_demes],s+1);
3377 n_demes++;
3378 }
3379 }
3380
3381 // n_demes is the number of non-empty sampling sites
3382
3383
3384 PhyML_Fprintf(fp,"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>");
3385 PhyML_Fprintf(fp,"\n<beast beautitemplate='MultiTypeTree' beautistatus='' namespace=\"beast.core:beast.evolution.alignment:beast.evolution.tree.coalescent:beast.core.util:beast.evolution.nuc:beast.evolution.operators:beast.evolution.sitemodel:beast.evolution.substitutionmodel:beast.evolution.likelihood\" version=\"2.0\">");
3386
3387 PhyML_Fprintf(fp,"\n<data id=\"data\" name=\"alignment\">");
3388
3389
3390 for(i=0;i<tree->n_otu;i++)
3391 {
3392 PhyML_Fprintf(fp,"\n<sequence id=\"%s\" taxon=\"%s\" totalcount=\"4\" value=\"%s\"/>",
3393 tree->a_nodes[i]->ldsk->coord->id,
3394 /* tree->a_nodes[i]->coord->id, */
3395 tree->a_nodes[i]->name,
3396 tree->a_nodes[i]->c_seq->state);
3397 }
3398
3399 PhyML_Fprintf(fp,"\n</data>");
3400
3401 PhyML_Fprintf(fp,"\n<map name=\"Uniform\" >beast.math.distributions.Uniform</map>");
3402 PhyML_Fprintf(fp,"\n<map name=\"Exponential\" >beast.math.distributions.Exponential</map>");
3403 PhyML_Fprintf(fp,"\n<map name=\"LogNormal\" >beast.math.distributions.LogNormalDistributionModel</map>");
3404 PhyML_Fprintf(fp,"\n<map name=\"Normal\" >beast.math.distributions.Normal</map>");
3405 PhyML_Fprintf(fp,"\n<map name=\"Beta\" >beast.math.distributions.Beta</map>");
3406 PhyML_Fprintf(fp,"\n<map name=\"Gamma\" >beast.math.distributions.Gamma</map>");
3407 PhyML_Fprintf(fp,"\n<map name=\"LaplaceDistribution\" >beast.math.distributions.LaplaceDistribution</map>");
3408 PhyML_Fprintf(fp,"\n<map name=\"prior\" >beast.math.distributions.Prior</map>");
3409 PhyML_Fprintf(fp,"\n<map name=\"InverseGamma\" >beast.math.distributions.InverseGamma</map>");
3410 PhyML_Fprintf(fp,"\n<map name=\"OneOnX\" >beast.math.distributions.OneOnX</map>");
3411
3412
3413 PhyML_Fprintf(fp,"\n<run id=\"mcmc\" spec=\"MCMC\" chainLength=\"1000000000\">");
3414 PhyML_Fprintf(fp,"\n<state id=\"state\" storeEvery=\"10000\">");
3415 PhyML_Fprintf(fp,"\n<stateNode id=\"Tree.t:data\" spec=\"beast.evolution.tree.StructuredCoalescentMultiTypeTree\">");
3416 PhyML_Fprintf(fp,"\n<migrationModel id=\"migModelInit.t:data\" spec=\"beast.evolution.tree.MigrationModel\">");
3417
3418
3419 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3420 For(i,n_demes*(n_demes-1)) strcat(s,"1.0 ");
3421 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.0\" dimension=\"%d\" estimate=\"false\" name=\"rateMatrix\">%s</parameter>",n_demes*(n_demes-1),s);
3422 Free(s);
3423
3424 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3425 for(i=0;i<n_demes;i++) strcat(s,"1.0 ");
3426 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.01\" dimension=\"%d\" estimate=\"false\" name=\"popSizes\">%s</parameter>",n_demes,s);
3427 Free(s);
3428
3429 PhyML_Fprintf(fp,"\n</migrationModel>");
3430
3431 PhyML_Fprintf(fp,"\n<typeTrait id=\"typeTraitSet.t:data\" spec=\"beast.evolution.tree.TraitSet\" traitname=\"type\" value=\"");
3432
3433 for(i=0;i<tree->n_otu;i++)
3434 {
3435 s = strrchr(tree->a_nodes[i]->ldsk->coord->id,'_');
3436 PhyML_Fprintf(fp,"%s=%s",
3437 tree->a_nodes[i]->ldsk->coord->id,
3438 s+1);
3439
3440 if(i < tree->n_otu-1) PhyML_Fprintf(fp,",");
3441 else PhyML_Fprintf(fp,"\">");
3442 }
3443
3444 PhyML_Fprintf(fp,"\n<taxa id=\"TaxonSet.0\" spec=\"TaxonSet\">");
3445 PhyML_Fprintf(fp,"\n<alignment idref=\"data\"/>");
3446 PhyML_Fprintf(fp,"\n</taxa>");
3447 PhyML_Fprintf(fp,"\n</typeTrait>");
3448 PhyML_Fprintf(fp,"\n<taxonset idref=\"TaxonSet.0\"/>");
3449 PhyML_Fprintf(fp,"\n</stateNode>");
3450 PhyML_Fprintf(fp,"\n<parameter id=\"kappa.s:data\" lower=\"0.0\" name=\"stateNode\">2.0</parameter>");
3451
3452 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3453 for(i=0;i<n_demes;i++) strcat(s,"1.0 ");
3454 PhyML_Fprintf(fp,"\n<parameter id=\"popSizes.t:data\" dimension=\"%d\" name=\"stateNode\">%s</parameter>",n_demes,s);
3455 Free(s);
3456
3457 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3458 For(i,n_demes*(n_demes-1)) strcat(s,"1.0 ");
3459 PhyML_Fprintf(fp,"\n<parameter id=\"rateMatrix.t:data\" dimension=\"%d\" name=\"stateNode\">%s</parameter>",n_demes*(n_demes-1),s);
3460 Free(s);
3461
3462
3463 PhyML_Fprintf(fp,"\n<parameter id=\"freqParameter.s:data\" dimension=\"4\" lower=\"0.0\" name=\"stateNode\" upper=\"1.0\">0.25</parameter>");
3464 PhyML_Fprintf(fp,"\n</state>");
3465
3466
3467 PhyML_Fprintf(fp,"\n<distribution id=\"posterior\" spec=\"util.CompoundDistribution\">");
3468 PhyML_Fprintf(fp,"\n<distribution id=\"prior\" spec=\"util.CompoundDistribution\">");
3469 PhyML_Fprintf(fp,"\n<prior id=\"KappaPrior.s:data\" name=\"distribution\" x=\"@kappa.s:data\">");
3470 PhyML_Fprintf(fp,"\n<LogNormal id=\"LogNormalDistributionModel.0\" name=\"distr\">");
3471 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.02\" estimate=\"false\" name=\"M\">1.0</parameter>");
3472 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.03\" estimate=\"false\" name=\"S\">1.25</parameter>");
3473 PhyML_Fprintf(fp,"\n</LogNormal>");
3474 PhyML_Fprintf(fp,"\n</prior>");
3475
3476 PhyML_Fprintf(fp,"\n<prior id=\"popSizesPrior.t:data\" name=\"distribution\" x=\"@popSizes.t:data\">");
3477 PhyML_Fprintf(fp,"\n<LogNormal id=\"LogNormalDistributionModel.01\" name=\"distr\">");
3478 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.04\" estimate=\"false\" name=\"M\">1.0</parameter>");
3479 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.05\" estimate=\"false\" lower=\"0.0\" name=\"S\" upper=\"5.0\">1.25</parameter>");
3480 PhyML_Fprintf(fp,"\n</LogNormal>");
3481 PhyML_Fprintf(fp,"\n</prior>");
3482
3483 PhyML_Fprintf(fp,"\n<prior id=\"rateMatrixPrior.t:data\" name=\"distribution\" x=\"@rateMatrix.t:data\">");
3484 PhyML_Fprintf(fp,"\n<LogNormal id=\"LogNormalDistributionModel.02\" name=\"distr\">");
3485 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.06\" estimate=\"false\" name=\"M\">1.0</parameter>");
3486 PhyML_Fprintf(fp,"\n<parameter id=\"RealParameter.07\" estimate=\"false\" lower=\"0.0\" name=\"S\" upper=\"5.0\">1.25</parameter>");
3487 PhyML_Fprintf(fp,"\n</LogNormal>");
3488 PhyML_Fprintf(fp,"\n</prior>");
3489
3490 PhyML_Fprintf(fp,"\n<distribution id=\"structuredCoalescent.t:data\" spec=\"multitypetree.distributions.StructuredCoalescentTreeDensity\" multiTypeTree=\"@Tree.t:data\">");
3491 PhyML_Fprintf(fp,"\n<migrationModel id=\"migModel.t:data\" spec=\"beast.evolution.tree.MigrationModel\" popSizes=\"@popSizes.t:data\" rateMatrix=\"@rateMatrix.t:data\">");
3492 PhyML_Fprintf(fp,"\n</migrationModel>");
3493 PhyML_Fprintf(fp,"\n</distribution>");
3494
3495 PhyML_Fprintf(fp,"\n<distribution id=\"likelihood\" spec=\"util.CompoundDistribution\">");
3496 PhyML_Fprintf(fp,"\n<distribution id=\"treeLikelihood.data\" spec=\"TreeLikelihood\" data=\"@data\" tree=\"@Tree.t:data\">");
3497 PhyML_Fprintf(fp,"\n<siteModel id=\"SiteModel.s:data\" spec=\"SiteModel\">");
3498 PhyML_Fprintf(fp,"\n<parameter id=\"mutationRate.s:data\" estimate=\"false\" name=\"mutationRate\">1.0</parameter>");
3499 PhyML_Fprintf(fp,"\n<parameter id=\"gammaShape.s:data\" estimate=\"false\" name=\"shape\">1.0</parameter>");
3500 PhyML_Fprintf(fp,"\n<parameter id=\"proportionInvariant.s:data\" estimate=\"false\" lower=\"0.0\" name=\"proportionInvariant\" upper=\"1.0\">0.0</parameter>");
3501 PhyML_Fprintf(fp,"\n<substModel id=\"hky.s:data\" spec=\"HKY\" kappa=\"@kappa.s:data\">");
3502 PhyML_Fprintf(fp,"\n<frequencies id=\"estimatedFreqs.s:data\" spec=\"Frequencies\" frequencies=\"@freqParameter.s:data\"/>");
3503 PhyML_Fprintf(fp,"\n</substModel>");
3504 PhyML_Fprintf(fp,"\n</siteModel>");
3505 PhyML_Fprintf(fp,"\n<branchRateModel id=\"StrictClock.c:data\" spec=\"beast.evolution.branchratemodel.StrictClockModel\">");
3506 PhyML_Fprintf(fp,"\n<parameter id=\"clockRate.c:data\" estimate=\"false\" name=\"clock.rate\">%G</parameter>",1.0);
3507 PhyML_Fprintf(fp,"\n</branchRateModel>");
3508 PhyML_Fprintf(fp,"\n</distribution>");
3509 PhyML_Fprintf(fp,"\n</distribution>");
3510 PhyML_Fprintf(fp,"\n</distribution>");
3511 PhyML_Fprintf(fp,"\n</distribution>");
3512 PhyML_Fprintf(fp,"\n");
3513 PhyML_Fprintf(fp,"\n<operator id=\"STX.t:data\" spec=\"multitypetree.operators.TypedSubtreeExchange\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" weight=\"10.0\"/>");
3514 PhyML_Fprintf(fp,"\n<operator id=\"TWB.t:data\" spec=\"multitypetree.operators.TypedWilsonBalding\" alpha=\"0.2\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" weight=\"10.0\"/>");
3515 PhyML_Fprintf(fp,"\n<operator id=\"NR.t:data\" spec=\"multitypetree.operators.NodeRetype\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" weight=\"10.0\"/>");
3516 PhyML_Fprintf(fp,"\n<operator id=\"NSR1.t:data\" spec=\"multitypetree.operators.NodeShiftRetype\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" rootOnly=\"true\" weight=\"10.0\"/>");
3517 PhyML_Fprintf(fp,"\n<operator id=\"NSR2.t:data\" spec=\"multitypetree.operators.NodeShiftRetype\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" noRoot=\"true\" weight=\"10.0\"/>");
3518 PhyML_Fprintf(fp,"<operator id=\"MTU.t:data\" spec=\"multitypetree.operators.MultiTypeUniform\" includeRoot=\"true\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" weight=\"10.0\"/>\n");
3519 PhyML_Fprintf(fp,"<operator id=\"MTTS.t:data\" spec=\"multitypetree.operators.MultiTypeTreeScale\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" scaleFactor=\"0.98\" useOldTreeScaler=\"true\" weight=\"10.0\"/>\n");
3520 PhyML_Fprintf(fp,"\n<operator id=\"MTTUpDown.t:data\" spec=\"multitypetree.operators.MultiTypeTreeScale\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\" scaleFactor=\"0.98\" useOldTreeScaler=\"true\" weight=\"10.0\">");
3521 PhyML_Fprintf(fp,"\n<parameter idref=\"popSizes.t:data\"/>");
3522 PhyML_Fprintf(fp,"\n</operator>");
3523 PhyML_Fprintf(fp,"\n<operator id=\"KappaScaler.s:data\" spec=\"ScaleOperator\" parameter=\"@kappa.s:data\" scaleFactor=\"0.5\" weight=\"0.1\"/>");
3524 PhyML_Fprintf(fp,"\n<operator id=\"popSizesScaler.t:data\" spec=\"ScaleOperator\" parameter=\"@popSizes.t:data\" scaleFactor=\"0.8\" weight=\"1.0\"/>");
3525 PhyML_Fprintf(fp,"\n<operator id=\"rateMatrixScaler.t:data\" spec=\"ScaleOperator\" parameter=\"@rateMatrix.t:data\" scaleFactor=\"0.8\" weight=\"1.0\"/>");
3526 PhyML_Fprintf(fp,"\n<operator id=\"FrequenciesExchanger.s:data\" spec=\"DeltaExchangeOperator\" delta=\"0.01\" weight=\"0.1\">");
3527 PhyML_Fprintf(fp,"\n<parameter idref=\"freqParameter.s:data\"/>");
3528 PhyML_Fprintf(fp,"\n</operator>");
3529 PhyML_Fprintf(fp,"\n");
3530 PhyML_Fprintf(fp,"\n<logger id=\"tracelog\" fileName=\"$(filebase).log\" logEvery=\"10000\">");
3531 PhyML_Fprintf(fp,"\n<log idref=\"likelihood\"/>");
3532 PhyML_Fprintf(fp,"\n<log idref=\"prior\"/>");
3533 PhyML_Fprintf(fp,"\n<log idref=\"treeLikelihood.data\"/>");
3534 PhyML_Fprintf(fp,"\n<log id=\"treeHeight.t:data\" spec=\"beast.evolution.tree.TreeHeightLogger\" tree=\"@Tree.t:data\"/>");
3535 /* PhyML_Fprintf(fp,"\n<log id=\"treeLength.t:data\" spec=\"multitypetree.util.TreeLengthLogger\" tree=\"@Tree.t:data\"/>"); */
3536 /* PhyML_Fprintf(fp,"\n<log id=\"changeCounts.t:data\" spec=\"multitypetree.util.TypeChangeCounts\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\"/>"); */
3537 /* PhyML_Fprintf(fp,"\n<log id=\"rootTypeLogger.t:data\" spec=\"multitypetree.util.TreeRootTypeLogger\" multiTypeTree=\"@Tree.t:data\"/>"); */
3538 PhyML_Fprintf(fp,"\n<log id=\"migModelLogger.t:data\" spec=\"multitypetree.util.MigrationModelLogger\" migrationModel=\"@migModel.t:data\" multiTypeTree=\"@Tree.t:data\"/>");
3539 PhyML_Fprintf(fp,"\n</logger>");
3540 PhyML_Fprintf(fp,"\n");
3541 PhyML_Fprintf(fp,"\n<logger id=\"screenlog\" logEvery=\"50000\">");
3542 PhyML_Fprintf(fp,"\n<log id=\"ESS.0\" spec=\"util.ESS\" arg=\"@posterior\"/>");
3543 PhyML_Fprintf(fp,"\n<log idref=\"likelihood\"/>");
3544 PhyML_Fprintf(fp,"\n</logger>");
3545 PhyML_Fprintf(fp,"\n");
3546 PhyML_Fprintf(fp,"\n</run>");
3547 PhyML_Fprintf(fp,"\n</beast>");
3548
3549 for(i=0;i<n_demes;i++) Free(deme_names[i]);
3550 Free(deme_names);
3551
3552 fclose(fp);
3553 }
3554
3555 /*////////////////////////////////////////////////////////////
3556 ////////////////////////////////////////////////////////////*/
3557
PHYREX_Number_Of_Sampled_Demes(t_tree * tree)3558 int PHYREX_Number_Of_Sampled_Demes(t_tree *tree)
3559 {
3560 int i,j,n_demes;
3561 t_dsk *disk;
3562 char **deme_list;
3563
3564 deme_list = (char **)mCalloc(tree->n_otu,sizeof(char *));
3565
3566 disk = tree->young_disk;
3567
3568 n_demes = 0;
3569 for(i=0;i<tree->n_otu;i++)
3570 {
3571 for(j=0;j<n_demes;j++)
3572 {
3573 if(deme_list[j] != NULL && !strcmp(strstr(disk->ldsk_a[i]->coord->id,"_deme"),deme_list[j]))
3574 {
3575 break;
3576 }
3577 }
3578
3579 if(j == n_demes)
3580 {
3581 deme_list[j] = strstr(disk->ldsk_a[i]->coord->id,"_deme");
3582 /* printf("\n. deme_list[%d]: %s",j,deme_list[j]); */
3583 n_demes++;
3584 }
3585 }
3586
3587 Free(deme_list);
3588
3589
3590 return(n_demes);
3591 }
3592
3593 /*////////////////////////////////////////////////////////////
3594 ////////////////////////////////////////////////////////////*/
3595 /* Returns the number of ldsks going from disk towards the past */
PHYREX_Number_Of_Outgoing_Ldsks(t_dsk * disk)3596 int PHYREX_Number_Of_Outgoing_Ldsks(t_dsk *disk)
3597 {
3598 int i;
3599
3600 if(disk->ldsk == NULL) return disk->n_ldsk_a;
3601 else
3602 {
3603 int n_out;
3604
3605 n_out = 0;
3606 for(i=0;i<disk->n_ldsk_a;++i)
3607 {
3608 if(disk->ldsk_a[i] && disk->ldsk_a[i]->prev != disk->ldsk)
3609 {
3610 n_out++;
3611 }
3612 }
3613 return(n_out+1); /* +1 so as to count disk->ldsk in */
3614 }
3615 }
3616
3617 /*////////////////////////////////////////////////////////////
3618 ////////////////////////////////////////////////////////////*/
3619 /* Select uniformly at random a lineage going "out of" disk (towards the past) */
PHYREX_Random_Select_Outgoing_Ldsk(t_dsk * disk)3620 t_ldsk *PHYREX_Random_Select_Outgoing_Ldsk(t_dsk *disk)
3621 {
3622 int i,*permut,n_ldsk_a_out;
3623 t_ldsk **ldsk_a_out,*target_ldsk;
3624
3625 ldsk_a_out = (t_ldsk **)mCalloc(disk->n_ldsk_a,sizeof(t_ldsk *));
3626
3627 n_ldsk_a_out = 0;
3628 if(disk->ldsk != NULL)
3629 {
3630 ldsk_a_out[0] = disk->ldsk;
3631 n_ldsk_a_out = 1;
3632 }
3633
3634 for(i=0;i<disk->n_ldsk_a;++i)
3635 {
3636 if(disk->ldsk_a[i] && disk->ldsk_a[i]->prev != disk->ldsk)
3637 {
3638 ldsk_a_out[n_ldsk_a_out] = disk->ldsk_a[i];
3639 n_ldsk_a_out++;
3640 }
3641 }
3642
3643 permut = Permutate(n_ldsk_a_out);
3644 target_ldsk = ldsk_a_out[permut[0]];
3645
3646 Free(permut);
3647 Free(ldsk_a_out);
3648
3649 return(target_ldsk);
3650 }
3651
3652 /*////////////////////////////////////////////////////////////
3653 ////////////////////////////////////////////////////////////*/
3654 /* Remove and free all disks and ldsks except for sampled disks.
3655 Connect the sampled disks together */
3656
PHYREX_Strip_And_Reconnect_Tree(t_tree * tree)3657 void PHYREX_Strip_And_Reconnect_Tree(t_tree *tree)
3658 {
3659 t_dsk *disk,*prev;
3660 int i,orig_size;
3661
3662 disk = tree->young_disk;
3663 do
3664 {
3665 if(disk->age_fixed == YES)
3666 {
3667 orig_size = disk->n_ldsk_a;
3668 for(i=0;i<orig_size;++i)
3669 {
3670 if(disk->ldsk_a[i] &&
3671 disk->ldsk_a[i]->disk != disk)
3672 {
3673 disk->ldsk_a[i] = NULL;
3674 disk->n_ldsk_a--;
3675 }
3676 }
3677 for(i=0;i<disk->n_ldsk_a;++i) disk->ldsk_a[i]->prev = NULL;
3678 }
3679 disk = disk->prev;
3680 }
3681 while(disk);
3682
3683 disk = tree->young_disk;
3684 do
3685 {
3686 prev = disk->prev;
3687 if(disk->age_fixed == NO)
3688 {
3689 PHYREX_Remove_Disk(disk);
3690 Free_Ldisk(disk->ldsk);
3691 Free_Disk(disk);
3692 }
3693 disk = prev;
3694 }
3695 while(disk);
3696
3697 }
3698
3699 /*////////////////////////////////////////////////////////////
3700 ////////////////////////////////////////////////////////////*/
3701
PHYREX_Scale_All(phydbl scale,t_dsk * start_disk,t_tree * tree)3702 int PHYREX_Scale_All(phydbl scale, t_dsk *start_disk, t_tree *tree)
3703 {
3704 t_dsk *disk;
3705 t_dsk **sorted_disk;
3706 int n_disk,n_disks_scaled,n_nodes_scaled,sorted,i;
3707
3708 n_disk = 0;
3709 n_disks_scaled = 0;
3710 n_nodes_scaled = 0;
3711
3712 disk = start_disk->prev;
3713 assert(disk);
3714
3715 do
3716 {
3717 if(disk->age_fixed == NO)
3718 {
3719 disk->time = disk->time * scale + start_disk->time * (1.-scale);
3720 n_disks_scaled++;
3721 if(disk->ldsk && disk->ldsk->n_next > 1) n_nodes_scaled++;
3722 }
3723
3724 n_disk++;
3725 disk = disk->prev;
3726 }
3727 while(disk);
3728
3729 sorted_disk = (t_dsk **)mCalloc(n_disk,sizeof(t_dsk *));
3730 disk = start_disk->prev;
3731 n_disk = 0;
3732 do
3733 {
3734 sorted_disk[n_disk] = disk;
3735 n_disk++;
3736 disk = disk->prev;
3737 }
3738 while(disk);
3739
3740
3741 do
3742 {
3743 sorted = YES;
3744 for(i=0;i<n_disk-1;++i)
3745 {
3746 if(sorted_disk[i]->time < sorted_disk[i+1]->time)
3747 {
3748 sorted = NO;
3749 disk = sorted_disk[i];
3750 sorted_disk[i] = sorted_disk[i+1];
3751 sorted_disk[i+1] = disk;
3752 }
3753 }
3754 }
3755 while(sorted == NO);
3756
3757
3758 for(i=0;i<n_disk;++i)
3759 {
3760 PHYREX_Remove_Disk(sorted_disk[i]);
3761 PHYREX_Insert_Disk(sorted_disk[i],tree);
3762 }
3763
3764 Free(sorted_disk);
3765
3766 if(tree->mmod->id == SLFV_GAUSSIAN || tree->mmod->id == SLFV_UNIFORM) return(n_disks_scaled);
3767 else if(tree->mmod->id == RW || tree->mmod->id == RRW) return(n_nodes_scaled);
3768 else return(-1);
3769 }
3770
3771 /*////////////////////////////////////////////////////////////
3772 ////////////////////////////////////////////////////////////*/
3773
PHYREX_Oldest_Sampled_Disk(t_tree * tree)3774 void PHYREX_Oldest_Sampled_Disk(t_tree *tree)
3775 {
3776 t_dsk *disk;
3777
3778 disk = tree->young_disk;
3779 do
3780 {
3781 if(disk->age_fixed == YES) tree->old_samp_disk = disk;
3782 disk = disk->prev;
3783 }
3784 while(disk);
3785 }
3786
3787 /*////////////////////////////////////////////////////////////
3788 ////////////////////////////////////////////////////////////*/
3789
PHYREX_Time_Of_Descendants(t_ldsk * ldsk,t_tree * tree)3790 phydbl PHYREX_Time_Of_Descendants(t_ldsk *ldsk, t_tree *tree)
3791 {
3792 if(ldsk->disk->age_fixed == YES || ldsk->disk == tree->young_disk) return ldsk->disk->time;
3793 else
3794 {
3795 int i;
3796 phydbl t,min;
3797
3798 min = tree->young_disk->time;
3799
3800 for(i=0;i<ldsk->n_next;++i)
3801 {
3802 t = PHYREX_Time_Of_Descendants(ldsk->next[i],tree);
3803 if(t < min) min = t;
3804 }
3805 return(min);
3806 }
3807 return(1.);
3808 }
3809
3810 /*////////////////////////////////////////////////////////////
3811 ////////////////////////////////////////////////////////////*/
3812
PHYREX_Time_Of_Prev_Sampled_Disk(t_dsk * disk,t_tree * tree)3813 phydbl PHYREX_Time_Of_Prev_Sampled_Disk(t_dsk *disk, t_tree *tree)
3814 {
3815
3816 if(disk->prev == NULL) return -INFINITY;
3817 else
3818 {
3819 disk = disk->prev;
3820 while(disk && disk->age_fixed == NO) disk = disk->prev;
3821 if(!disk) return -INFINITY;
3822 else return(disk->time);
3823 }
3824 }
3825
3826 /*////////////////////////////////////////////////////////////
3827 ////////////////////////////////////////////////////////////*/
3828
PHYREX_Time_Of_Next_Sampled_Disk(t_dsk * disk,t_tree * tree)3829 phydbl PHYREX_Time_Of_Next_Sampled_Disk(t_dsk *disk, t_tree *tree)
3830 {
3831
3832 if(disk->next == NULL) return tree->young_disk->time;
3833 else
3834 {
3835 disk = disk->next;
3836 while(disk && disk->age_fixed == NO) disk = disk->next;
3837 if(!disk) return tree->young_disk->time;
3838 else return(disk->time);
3839 }
3840 }
3841
3842 /*////////////////////////////////////////////////////////////
3843 ////////////////////////////////////////////////////////////*/
3844
PHYREX_Lk_Time_Component(t_tree * tree)3845 phydbl PHYREX_Lk_Time_Component(t_tree *tree)
3846 {
3847 phydbl lnL;
3848 t_dsk *disk;
3849 phydbl dt;
3850 int n_evt;
3851
3852 n_evt = 0;
3853 dt = 0.0;
3854 lnL = 0.0;
3855 disk = tree->young_disk->prev;
3856 do
3857 {
3858
3859 if(disk->age_fixed == NO)
3860 {
3861 dt += fabs(disk->next->time - disk->time);
3862 n_evt++;
3863 }
3864
3865 disk = disk->prev;
3866
3867 }
3868 while(disk);
3869
3870 lnL += n_evt * log(tree->mmod->lbda) - tree->mmod->lbda * dt;
3871
3872 return(lnL);
3873 }
3874
3875 /*////////////////////////////////////////////////////////////
3876 ////////////////////////////////////////////////////////////*/
3877
PHYREX_Find_Ldsk_From_Id(char * id,t_ldsk * root)3878 t_ldsk *PHYREX_Find_Ldsk_From_Id(char *id, t_ldsk *root)
3879 {
3880 t_ldsk *ldsk;
3881
3882 ldsk = NULL;
3883
3884 if(!strcmp(root->coord->id,id)) return(root);
3885
3886 if(root->n_next == 0) return(NULL);
3887
3888 for(int i=0;i<root->n_next;++i)
3889 {
3890 ldsk = PHYREX_Find_Ldsk_From_Id(id,root->next[i]);
3891 if(ldsk != NULL) break;
3892 }
3893
3894 return(ldsk);
3895 }
3896
3897
3898 /*////////////////////////////////////////////////////////////
3899 ////////////////////////////////////////////////////////////*/
3900
PHYREX_Mean_Next_Loc(t_ldsk * ldsk,t_tree * tree)3901 t_geo_coord *PHYREX_Mean_Next_Loc(t_ldsk *ldsk, t_tree *tree)
3902 {
3903 t_geo_coord *mean;
3904 int i,j;
3905
3906 mean = GEO_Make_Geo_Coord(tree->mmod->n_dim);
3907
3908 assert(ldsk->n_next > 0);
3909
3910 for(i=0;i<ldsk->n_next;++i)
3911 {
3912 for(j=0;j<tree->mmod->n_dim;++j)
3913 {
3914 mean->lonlat[j] += ldsk->next[i]->coord->lonlat[j];
3915 }
3916 }
3917
3918 for(j=0;j<tree->mmod->n_dim;++j)
3919 {
3920 mean->lonlat[j] /= (phydbl)ldsk->n_next;
3921 }
3922
3923 return(mean);
3924 }
3925
3926 /*////////////////////////////////////////////////////////////
3927 ////////////////////////////////////////////////////////////*/
3928 // realized sigsq = (1/d) E(\sum_x D_x^2), where $D_x$ is the
3929 // square euclidean distance along the x-axis between the
3930 // location of a lineage at time $t$ and time $t+1$.
PHYREX_Root_To_Tip_Realized_Sigsq(t_tree * tree)3931 phydbl PHYREX_Root_To_Tip_Realized_Sigsq(t_tree *tree)
3932 {
3933 t_dsk *disk,*root_disk;
3934 int i;
3935 phydbl res,*sigsq;
3936 phydbl t_root, t_tip;
3937
3938 sigsq = (phydbl *)mCalloc(tree->young_disk->n_ldsk_a,sizeof(phydbl));
3939
3940 disk = tree->young_disk;
3941 while(disk->prev) disk = disk->prev;
3942 root_disk = disk;
3943
3944 t_root = root_disk->time;
3945 t_tip = tree->young_disk->time;
3946 assert(t_tip - t_root > 0.0);
3947
3948 disk = tree->young_disk;
3949 for(i=0;i<disk->n_ldsk_a;++i)
3950 sigsq[i] = pow(Euclidean_Dist(root_disk->ldsk->coord,disk->ldsk_a[i]->coord) / (t_tip - t_root),2);
3951
3952 res = Quantile(sigsq,disk->n_ldsk_a,0.5);
3953
3954 Free(sigsq);
3955
3956 return(res/tree->mmod->n_dim);
3957 }
3958
3959 /*////////////////////////////////////////////////////////////
3960 ////////////////////////////////////////////////////////////*/
3961 // Mean Haversine distance (in km) per year (measured on path between root and tips on the most recent disk)
PHYREX_Realized_Dispersal_Dist(t_tree * tree)3962 phydbl PHYREX_Realized_Dispersal_Dist(t_tree *tree)
3963 {
3964 t_dsk *disk,*root_disk;
3965 phydbl dist,t_tip,t_root;
3966 int i;
3967
3968 disk = tree->young_disk;
3969 while(disk->prev) disk = disk->prev;
3970 root_disk = disk;
3971
3972 t_root = root_disk->time;
3973 t_tip = tree->young_disk->time;
3974 assert(t_tip - t_root > 0.0);
3975 \
3976 disk = tree->young_disk;
3977 dist = 0.0;
3978 for(i=0;i<disk->n_ldsk_a;++i)
3979 {
3980 dist += Haversine_Distance(root_disk->ldsk->coord->lonlat[0],
3981 root_disk->ldsk->coord->lonlat[1],
3982 disk->ldsk_a[i]->coord->lonlat[0],
3983 disk->ldsk_a[i]->coord->lonlat[1])/(t_tip - t_root);
3984
3985 }
3986
3987 return(dist/disk->n_ldsk_a);
3988 }
3989
3990 /*////////////////////////////////////////////////////////////
3991 ////////////////////////////////////////////////////////////*/
3992
PHYREX_Tip_To_Root_Realized_Sigsq(t_tree * tree)3993 phydbl PHYREX_Tip_To_Root_Realized_Sigsq(t_tree *tree)
3994 {
3995 t_dsk *disk;
3996 phydbl *sigsq,res;
3997 int i;
3998
3999
4000 sigsq = (phydbl *)mCalloc(tree->young_disk->n_ldsk_a,sizeof(phydbl));
4001
4002 disk = tree->young_disk;
4003 for(i=0;i<disk->n_ldsk_a;++i)
4004 {
4005 sigsq[i] =
4006 pow(Euclidean_Dist(disk->ldsk_a[i]->coord,
4007 disk->ldsk_a[i]->prev->coord)/fabs(disk->time - disk->ldsk_a[i]->prev->disk->time),2);
4008
4009 }
4010
4011 res = Quantile(sigsq,disk->n_ldsk_a,0.5);
4012
4013 Free(sigsq);
4014
4015 return(res/tree->mmod->n_dim);
4016 }
4017
4018 /*////////////////////////////////////////////////////////////
4019 ////////////////////////////////////////////////////////////*/
4020
PHYREX_Tip_To_Root_Realized_Bis_Sigsq(t_tree * tree)4021 phydbl PHYREX_Tip_To_Root_Realized_Bis_Sigsq(t_tree *tree)
4022 {
4023 t_dsk *disk;
4024 phydbl sumdist,sumt;
4025 int i;
4026
4027 sumdist = 0.0;
4028 sumt = 0.0;
4029 disk = tree->young_disk;
4030 for(i=0;i<disk->n_ldsk_a;++i)
4031 {
4032 sumdist += pow(Euclidean_Dist(disk->ldsk_a[i]->coord,
4033 disk->ldsk_a[i]->prev->coord),2);
4034 sumt += tree->mmod->n_dim * fabs(disk->time - disk->ldsk_a[i]->prev->disk->time);
4035 }
4036
4037 return(sumdist/sumt);
4038 }
4039
4040 /*////////////////////////////////////////////////////////////
4041 ////////////////////////////////////////////////////////////*/
4042
PHYREX_Tip_To_Root_Realized_Ter_Sigsq(t_tree * tree)4043 phydbl PHYREX_Tip_To_Root_Realized_Ter_Sigsq(t_tree *tree)
4044 {
4045 t_ldsk *ldsk;
4046 phydbl mean,dist;
4047 int i;
4048
4049 mean = 0.0;
4050 for(i=0;i<tree->young_disk->n_ldsk_a;++i)
4051 {
4052 ldsk = tree->young_disk->ldsk_a[i];
4053 while(ldsk && ldsk->prev && fabs(ldsk->prev->disk->time-tree->young_disk->ldsk_a[i]->disk->time) < 1.0)
4054 {
4055 ldsk = ldsk->prev;
4056 if(ldsk == NULL) return(-1.);
4057 }
4058 dist = Euclidean_Dist(ldsk->coord,tree->young_disk->ldsk_a[i]->coord);
4059 mean += pow(dist,2);
4060
4061 /* PhyML_Printf("\n. %3d (%12f %12f) (%12f %12f) dist: %12f mean: %12f", */
4062 /* i, */
4063 /* tree->young_disk->ldsk_a[i]->coord->lonlat[0], */
4064 /* tree->young_disk->ldsk_a[i]->coord->lonlat[1], */
4065 /* ldsk->coord->lonlat[0], */
4066 /* ldsk->coord->lonlat[1], */
4067 /* Euclidean_Dist(ldsk->coord,tree->young_disk->ldsk_a[i]->coord), */
4068 /* mean); */
4069
4070 }
4071
4072 return(mean/(tree->young_disk->n_ldsk_a*tree->mmod->n_dim));
4073 }
4074
4075 /*////////////////////////////////////////////////////////////
4076 ////////////////////////////////////////////////////////////*/
4077
PHYREX_Label_Nodes_With_Locations(t_tree * tree)4078 void PHYREX_Label_Nodes_With_Locations(t_tree *tree)
4079 {
4080 t_dsk *disk;
4081 t_node *n;
4082 phydbl lon,lat;
4083
4084
4085 lon = lat = -1.;
4086
4087 for(int i=0;i<tree->n_otu;++i)
4088 {
4089 if(tree->a_nodes[i]->label == NULL)
4090 {
4091 tree->a_nodes[i]->label = Make_Label();
4092 tree->a_nodes[i]->label->next = Make_Label();
4093 }
4094
4095 lat = tree->a_nodes[i]->ldsk->coord->lonlat[1];
4096 lon = tree->a_nodes[i]->ldsk->coord->lonlat[0];
4097
4098 sprintf(tree->a_nodes[i]->label->key,"&location");
4099 sprintf(tree->a_nodes[i]->label->val,"{%f,%f}",lat,lon);
4100 sprintf(tree->a_nodes[i]->label->next->key,"location");
4101 sprintf(tree->a_nodes[i]->label->next->val,"{%f,%f}",lat,lon);
4102 }
4103
4104 disk = tree->young_disk->prev;
4105
4106 do
4107 {
4108 if(disk->ldsk && disk->ldsk->nd != NULL)
4109 {
4110 n = disk->ldsk->nd;
4111 if(n->label == NULL)
4112 {
4113 n->label = Make_Label();
4114 n->label->next = Make_Label();
4115 }
4116
4117 lat = disk->ldsk->coord->lonlat[1];
4118 lon = disk->ldsk->coord->lonlat[0];
4119
4120 sprintf(n->label->key,"&location");
4121 sprintf(n->label->val,"{%f,%f}",lat,lon);
4122 sprintf(n->label->next->key,"location");
4123 sprintf(n->label->next->val,"{%f,%f}",lat,lon);
4124
4125 /* Print same label on all internal nodes with exactly the */
4126 /* same coalescence time. */
4127 for(int i=tree->n_otu;i<2*tree->n_otu-1;++i)
4128 {
4129 if(tree->a_nodes[i] != n &&
4130 Are_Equal(tree->times->nd_t[i],
4131 tree->times->nd_t[n->num],
4132 1.E-10) == YES)
4133 {
4134 n = tree->a_nodes[i];
4135 if(n->label == NULL)
4136 {
4137 n->label = Make_Label();
4138 n->label->next = Make_Label();
4139 }
4140
4141 lat = disk->ldsk->coord->lonlat[1];
4142 lon = disk->ldsk->coord->lonlat[0];
4143
4144 sprintf(n->label->key,"&location");
4145 sprintf(n->label->val,"{%f,%f}",lat,lon);
4146 sprintf(n->label->next->key,"location");
4147 sprintf(n->label->next->val,"{%f,%f}",lat,lon);
4148 }
4149 }
4150 }
4151
4152 disk = disk->prev;
4153 }
4154 while(disk);
4155 }
4156
4157 /*////////////////////////////////////////////////////////////
4158 ////////////////////////////////////////////////////////////*/
4159
PHYREX_Label_Edges(t_tree * tree)4160 void PHYREX_Label_Edges(t_tree *tree)
4161 {
4162 for(int i=0;i<2*tree->n_otu-1;++i)
4163 {
4164 if(tree->a_edges[i]->label == NULL)
4165 {
4166 tree->a_edges[i]->label = Make_Label();
4167 tree->a_edges[i]->label->next = Make_Label();
4168 }
4169
4170 sprintf(tree->a_edges[i]->label->key,"&rate");
4171 sprintf(tree->a_edges[i]->label->val,"0.0");
4172
4173 sprintf(tree->a_edges[i]->label->next->key,"location.rate");
4174 sprintf(tree->a_edges[i]->label->next->val,"0.0");
4175 }
4176
4177
4178 }
4179
4180 /*////////////////////////////////////////////////////////////
4181 ////////////////////////////////////////////////////////////*/
4182
PHYREX_Remove_All_Disks_Except_Coal_And_Tips(t_tree * tree)4183 void PHYREX_Remove_All_Disks_Except_Coal_And_Tips(t_tree *tree)
4184 {
4185 t_dsk *disk,*dum;
4186
4187 disk = tree->young_disk;
4188 dum = NULL;
4189
4190 do
4191 {
4192 dum = disk->prev;
4193
4194 if(disk->age_fixed == NO)
4195 {
4196 if(disk->ldsk == NULL)
4197 {
4198 PHYREX_Remove_Disk(disk);
4199 Free_Disk(disk);
4200 }
4201 else if(disk->ldsk != NULL && disk->ldsk->nd == NULL)
4202 {
4203 t_ldsk *old_ldsk,*young_ldsk;
4204 int dir_old_young;
4205
4206 assert(disk->ldsk->n_next == 1);
4207
4208 old_ldsk = disk->ldsk->prev;
4209 young_ldsk = disk->ldsk->next[0];
4210
4211 dir_old_young = PHYREX_Get_Next_Direction(disk->ldsk,old_ldsk);
4212 assert(dir_old_young != -1);
4213
4214 /* New connections between old_ldsk and young_ldsk */
4215 old_ldsk->next[dir_old_young] = young_ldsk;
4216 young_ldsk->prev = old_ldsk;
4217
4218 PHYREX_Remove_Disk(disk);
4219
4220 Free_Ldisk(disk->ldsk);
4221 Free_Disk(disk);
4222
4223 }
4224 }
4225
4226 disk = dum;
4227 }
4228 while(disk);
4229 }
4230
4231 /*////////////////////////////////////////////////////////////
4232 ////////////////////////////////////////////////////////////*/
4233
PHYREX_Path_Logdensity(t_ldsk * young,t_ldsk * old,phydbl sd,t_tree * tree)4234 phydbl PHYREX_Path_Logdensity(t_ldsk *young, t_ldsk *old, phydbl sd, t_tree *tree)
4235 {
4236 switch(tree->mmod->id)
4237 {
4238 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4239 {
4240 return(SLFV_Path_Logdensity(young,old,sd,tree));
4241 break;
4242 }
4243 case RW : case RRW :
4244 {
4245 return(0.0);
4246 break;
4247 }
4248 default : assert(FALSE);
4249 }
4250 return(-1.);
4251 }
4252
4253 /*////////////////////////////////////////////////////////////
4254 ////////////////////////////////////////////////////////////*/
4255
PHYREX_Sample_Path(t_ldsk * young,t_ldsk * old,phydbl sd,phydbl * global_hr,t_tree * tree)4256 void PHYREX_Sample_Path(t_ldsk *young, t_ldsk *old, phydbl sd, phydbl *global_hr, t_tree *tree)
4257 {
4258 switch(tree->mmod->id)
4259 {
4260 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4261 {
4262 SLFV_Sample_Path(young,old,sd,global_hr,tree);
4263 break;
4264 }
4265
4266 case RW : case RRW :
4267 {
4268 assert(FALSE);
4269 break;
4270 }
4271 default : assert(FALSE);
4272 }
4273 }
4274
4275 /*////////////////////////////////////////////////////////////
4276 ////////////////////////////////////////////////////////////*/
4277
PHYREX_Generate_Path(t_ldsk * young,t_ldsk * old,phydbl n_evt,phydbl sd,t_tree * tree)4278 t_ldsk *PHYREX_Generate_Path(t_ldsk *young, t_ldsk *old, phydbl n_evt, phydbl sd, t_tree *tree)
4279 {
4280 switch(tree->mmod->id)
4281 {
4282 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4283 {
4284 return(SLFV_Generate_Path(young,old,n_evt,sd,tree));
4285 break;
4286 }
4287 case RW : case RRW :
4288 {
4289 return(SLFV_Generate_Path(young,old,0,sd,tree));
4290 break;
4291 }
4292 default : assert(FALSE);
4293 }
4294 return(NULL);
4295 }
4296
4297 /*////////////////////////////////////////////////////////////
4298 ////////////////////////////////////////////////////////////*/
4299
PHYREX_Simulate(int n_otu,int n_sites,phydbl w,phydbl h,phydbl lbda,phydbl rad,phydbl mu,int r_seed,int modid)4300 t_tree *PHYREX_Simulate(int n_otu, int n_sites, phydbl w, phydbl h, phydbl lbda, phydbl rad, phydbl mu, int r_seed, int modid)
4301 {
4302 switch(modid)
4303 {
4304 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4305 {
4306 return(SLFV_Simulate(n_otu,n_sites,w,h,lbda,rad,mu,r_seed));
4307 break;
4308 }
4309
4310 case RW : case RRW :
4311 {
4312 assert(FALSE);
4313 break;
4314 }
4315 default : assert(FALSE);
4316 }
4317 return(NULL);
4318 }
4319
4320 /*////////////////////////////////////////////////////////////
4321 ////////////////////////////////////////////////////////////*/
4322
PHYREX_Simulate_Backward_Core(t_dsk * disk,int avoid_multiple_mergers,t_tree * tree)4323 void PHYREX_Simulate_Backward_Core(t_dsk *disk,int avoid_multiple_mergers, t_tree *tree)
4324 {
4325 switch(tree->mmod->id)
4326 {
4327 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4328 {
4329 SLFV_Simulate_Backward_Core(disk,avoid_multiple_mergers,tree);
4330 break;
4331 }
4332
4333 case RW : case RRW :
4334 {
4335 SLFV_Simulate_Backward_Core(disk,avoid_multiple_mergers,tree);
4336 break;
4337 }
4338
4339 default : assert(FALSE);
4340
4341 }
4342 }
4343
4344 /*////////////////////////////////////////////////////////////
4345 ////////////////////////////////////////////////////////////*/
4346
PHYREX_Update_Sigsq(t_tree * tree)4347 phydbl PHYREX_Update_Sigsq(t_tree *tree)
4348 {
4349 switch(tree->mmod->id)
4350 {
4351 case SLFV_UNIFORM : case SLFV_GAUSSIAN :
4352 {
4353 return(SLFV_Update_Sigsq(tree));
4354 break;
4355 }
4356
4357 case RW : case RRW :
4358 {
4359 return(tree->mmod->sigsq);
4360 break;
4361 }
4362
4363 default : assert(FALSE);
4364
4365 }
4366
4367 return(-1.);
4368 }
4369
4370 /*////////////////////////////////////////////////////////////
4371 ////////////////////////////////////////////////////////////*/
4372 /*////////////////////////////////////////////////////////////
4373 ////////////////////////////////////////////////////////////*/
4374 /*////////////////////////////////////////////////////////////
4375 ////////////////////////////////////////////////////////////*/
4376
4377