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
14 /* Routines for molecular clock trees and molecular dating */
15
16
17 #include "tiporder.h"
18
19 //////////////////////////////////////////////////////////////
20 //////////////////////////////////////////////////////////////
21
22
23 /* int TIPO_main(int argc, char **argv) */
24 /* { */
25 /* t_tree **list_tree,*ref_tree,*tree; */
26 /* FILE *fp_ref_tree,*fp_list_tree,*fp_coord,*ps_tree; */
27 /* int i,j,k; */
28 /* int n_trees; */
29 /* option *ref_io,*list_io; */
30 /* char **name_table; */
31 /* int r_seed; */
32
33 /* r_seed = time(NULL); */
34 /* srand(r_seed); */
35
36
37 /* ref_io = (option *)Make_Input(); */
38 /* list_io = (option *)Make_Input(); */
39
40 /* fp_ref_tree = (FILE *)fopen(argv[1],"r"); */
41 /* fp_list_tree = (FILE *)fopen(argv[2],"r"); */
42 /* fp_coord = (FILE *)fopen(argv[3],"r"); */
43
44 /* if(!fp_ref_tree) */
45 /* { */
46 /* PhyML_Printf("\n. Can't find %s",argv[1]); */
47 /* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
48 /* Warn_And_Exit(""); */
49 /* } */
50
51 /* if(!fp_list_tree) */
52 /* { */
53 /* PhyML_Printf("\n. Can't find %s",argv[2]); */
54 /* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
55 /* Warn_And_Exit(""); */
56 /* } */
57
58 /* ref_io->fp_in_tree = fp_ref_tree; */
59 /* Read_Tree_File(ref_io); */
60 /* fclose(ref_io->fp_in_tree); */
61 /* ref_tree = ref_io->treelist->tree[0]; */
62 /* ref_tree->io = ref_io; */
63
64
65 /* list_io->fp_in_tree = fp_list_tree; */
66 /* Read_Tree_File(list_io); */
67 /* fclose(list_io->fp_in_tree); */
68 /* list_tree = list_io->treelist->tree; */
69 /* n_trees = list_io->treelist->list_size; */
70 /* PhyML_Printf("\n. Read %d trees\n",n_trees); */
71
72 /* for(i=0;i<n_trees;i++) list_tree[i]->io = list_io; */
73
74 /* name_table = (char **)mCalloc(ref_tree->n_otu,sizeof(char **)); */
75 /* for(i=0;i<ref_tree->n_otu;i++) name_table[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char)); */
76
77
78 /* /\* Sort translation table such that tree->a_nodes[i]->name == tree->io->short_tax_name[i] for all i *\/ */
79 /* TIPO_Sort_Translation_Table(ref_tree); */
80
81 /* ref_tree->io->z_scores = (phydbl *)mCalloc(ref_tree->n_otu,sizeof(phydbl)); */
82
83 /* /\* TIPO_Read_Taxa_Zscores(fp_coord,ref_tree); *\/ */
84
85 /* for(i=0;i<ref_tree->n_otu;i++) ref_tree->io->z_scores[i] = TIPO_Read_One_Taxon_Zscore(fp_coord,ref_tree->a_nodes[i]->name,ref_tree); */
86
87 /* TIPO_Normalize_Zscores(ref_tree); */
88
89 /* /\* Find matching tips *\/ */
90 /* for(i=0;i<n_trees;i++) */
91 /* { */
92 /* TIPO_Sort_Translation_Table(list_tree[i]); */
93
94 /* for(j=0;j<ref_tree->n_otu;j++) */
95 /* { */
96 /* for(k=0;k<ref_tree->n_otu;k++) */
97 /* { */
98 /* if(!strcmp(ref_io->long_tax_names[j],list_io->long_tax_names[k])) */
99 /* { */
100 /* list_tree[i]->a_nodes[k]->ext_node = ref_tree->a_nodes[j]; */
101 /* break; */
102 /* } */
103 /* } */
104 /* if(k == ref_tree->n_otu) */
105 /* { */
106 /* PhyML_Printf("\n. Could not find matching tips for \"%s\" (tree %d)",ref_tree->a_nodes[j]->name,i); */
107 /* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
108 /* Warn_And_Exit(""); */
109 /* } */
110 /* } */
111 /* } */
112
113 /* PhyML_Printf("\n. Getting ancestors"); fflush(NULL); */
114 /* Update_Ancestors(ref_tree->n_root,ref_tree->n_root->v[2],ref_tree); */
115 /* Update_Ancestors(ref_tree->n_root,ref_tree->n_root->v[1],ref_tree); */
116
117 /* for(i=0;i<n_trees;i++) */
118 /* { */
119 /* Update_Ancestors(list_tree[i]->n_root,list_tree[i]->n_root->v[2],list_tree[i]); */
120 /* Update_Ancestors(list_tree[i]->n_root,list_tree[i]->n_root->v[1],list_tree[i]); */
121 /* list_tree[i]->n_root->anc = NULL; */
122 /* } */
123
124 /* PhyML_Printf("\n. Getting bipartitions"); fflush(NULL); */
125
126 /* Free_Bip(ref_tree); */
127 /* Alloc_Bip(ref_tree); */
128 /* Get_Bip(ref_tree->a_nodes[0], */
129 /* ref_tree->a_nodes[0]->v[0], */
130 /* ref_tree); */
131
132 /* for(i=0;i<n_trees;i++) */
133 /* { */
134 /* if(!(i%10)) */
135 /* PhyML_Printf("\n. Getting bipartition for tree %d",i); */
136 /* Free_Bip(list_tree[i]); */
137 /* Alloc_Bip(list_tree[i]); */
138 /* Get_Bip(list_tree[i]->a_nodes[0], */
139 /* list_tree[i]->a_nodes[0]->v[0], */
140 /* list_tree[i]); */
141 /* } */
142
143
144 /* PhyML_Printf("\n. Getting tip ranks"); fflush(NULL); */
145 /* /\* TIPO_Get_Tips_Y_Rank(ref_tree); *\/ */
146
147 /* TIPO_Get_Tips_Y_Rank_From_Zscores(ref_tree); */
148
149 /* /\* PhyML_Printf("\n. Minimizing"); fflush(NULL); *\/ */
150 /* /\* TIPO_Minimize_Tip_Order_Score(n_trees,list_tree,ref_tree); *\/ */
151
152 /* PhyML_Printf("\n. N_OTU = %d",ref_tree->n_otu); */
153 /* TIPO_Untangle_Tree(ref_tree); */
154 /* PhyML_Printf("\n ** ORIGINAL %f",ref_tree->tip_order_score); */
155
156 /* /\* i = 0; *\/ */
157 /* /\* do *\/ */
158 /* /\* { *\/ */
159 /* /\* /\\* TIPO_Get_Tips_Y_Rank_From_Zscores(ref_tree); *\\/ *\/ */
160 /* /\* TIPO_Randomize_Tip_Y_Ranks(ref_tree); *\/ */
161 /* /\* TIPO_Untangle_Tree(ref_tree); *\/ */
162 /* /\* PhyML_Printf("\n ** %f",ref_tree->tip_order_score); *\/ */
163 /* /\* i++; *\/ */
164 /* /\* }while(i < 5000); *\/ */
165
166 /* Test_Node_Table_Consistency(ref_tree); */
167
168 /* ref_tree->ps_tree = DR_Make_Tdraw_Struct(ref_tree); */
169 /* DR_Get_Tree_Coord(ref_tree); */
170 /* for(j=0;j<ref_tree->n_otu;j++) */
171 /* { */
172 /* ref_tree->ps_tree->ycoord[j] = */
173 /* (ref_tree->a_nodes[j]->y_rank/ref_tree->n_otu)* */
174 /* ref_tree->ps_tree->page_height; */
175 /* } */
176
177 /* ps_tree = (FILE *)fopen("order_tree.ps","w"); */
178
179
180 /* list_io->z_scores = (phydbl *)mCalloc(ref_tree->n_otu,sizeof(phydbl)); */
181
182 /* DR_Print_Postscript_Header(1,ps_tree); */
183 /* for(i=0;i<n_trees;i++) */
184 /* { */
185 /* tree = list_tree[i]; */
186
187 /* Test_Node_Table_Consistency(tree); */
188 /* tree->rates = RATES_Make_Rate_Struct(tree->n_otu); */
189 /* RATES_Init_Rate_Struct(tree->rates,tree->n_otu); */
190 /* TIMES_Least_Square_Node_Times(tree->e_root,tree); */
191 /* TIMES_Adjust_Node_Times(tree); */
192 /* RATES_Update_Cur_Bl(tree); */
193
194 /* tree->ps_tree = DR_Make_Tdraw_Struct(tree); */
195 /* DR_Init_Tdraw_Struct(tree->ps_tree); */
196 /* DR_Get_Tree_Box_Width(tree->ps_tree,tree); */
197 /* Dist_To_Root(tree); */
198 /* tree->ps_tree->max_dist_to_root = DR_Get_Max_Dist_To_Root(tree); */
199
200 /* for(j=0;j<ref_tree->n_otu;j++) tree->io->z_scores[j] = ref_tree->io->z_scores[tree->a_nodes[j]->ext_node->num]; */
201 /* TIPO_Get_Tips_Y_Rank_From_Zscores(tree); */
202 /* TIPO_Untangle_Tree(tree); */
203 /* for(j=0;j<ref_tree->n_otu;j++) tree->ps_tree->ycoord[j] = (tree->a_nodes[j]->y_rank/tree->n_otu)*tree->ps_tree->page_height; */
204
205 /* /\* for(j=0;j<ref_tree->n_otu;j++) tree->ps_tree->ycoord[j] = ref_tree->ps_tree->ycoord[tree->a_nodes[j]->ext_node->num]; *\/ */
206 /* for(j=0;j<ref_tree->n_otu;j++) list_io->z_scores[j] = ref_io->z_scores[tree->a_nodes[j]->ext_node->num]; */
207
208 /* DR_Get_Y_Coord(YES,tree->ps_tree,tree); */
209 /* DR_Get_X_Coord( NO,tree->ps_tree,tree); */
210
211 /* if(!i) DR_Print_Tree_Postscript(1,YES,ps_tree,tree); */
212 /* else DR_Print_Tree_Postscript(1, NO,ps_tree,tree); */
213 /* } */
214 /* DR_Print_Postscript_EOF(ps_tree); */
215 /* fclose(ps_tree); */
216
217
218 /* ps_tree = (FILE *)fopen("ref_tree.ps","w"); */
219
220 /* DR_Print_Postscript_Header(1,ps_tree); */
221 /* tree = ref_tree; */
222 /* tree->rates = RATES_Make_Rate_Struct(tree->n_otu); */
223 /* RATES_Init_Rate_Struct(tree->rates,tree->n_otu); */
224 /* TIMES_Least_Square_Node_Times(tree->e_root,tree); */
225 /* TIMES_Adjust_Node_Times(tree); */
226 /* RATES_Update_Cur_Bl(tree); */
227 /* DR_Init_Tdraw_Struct(tree->ps_tree); */
228 /* DR_Get_Tree_Box_Width(tree->ps_tree,tree); */
229 /* Dist_To_Root(tree->n_root,tree); */
230 /* tree->ps_tree->max_dist_to_root = DR_Get_Max_Dist_To_Root(tree); */
231
232 /* TIPO_Get_Tips_Y_Rank_From_Zscores(tree); */
233 /* /\* TIPO_Untangle_Tree(tree); *\/ */
234 /* for(j=0;j<tree->n_otu;j++) tree->ps_tree->ycoord[j] = (tree->a_nodes[j]->y_rank/tree->n_otu)*tree->ps_tree->page_height; */
235
236 /* DR_Get_Y_Coord(YES,tree->ps_tree,tree); */
237 /* DR_Get_X_Coord( NO,tree->ps_tree,tree); */
238 /* DR_Print_Tree_Postscript(1,YES,ps_tree,tree); */
239 /* DR_Print_Postscript_EOF(ps_tree); */
240 /* fclose(ps_tree); */
241
242
243 /* PhyML_Printf("\n"); */
244
245 /* fclose(fp_ref_tree); */
246 /* fclose(fp_list_tree); */
247 /* fclose(fp_coord); */
248 /* } */
249
TIPO_main(int argc,char ** argv)250 int TIPO_main(int argc, char **argv)
251 {
252 t_tree *tree;
253 option *io;
254 FILE *fp_tree_file, *fp_coord_file;
255 int i;
256
257 /* Rprintf("%s\n",tree_file_name[0]); */
258 /* Rprintf("%s\n",coord_file_name[0]); */
259
260
261 srand(time(NULL)); rand();
262
263
264 fp_tree_file = (FILE *)fopen(argv[1],"r");
265 fp_coord_file = (FILE *)fopen(argv[2],"r");
266
267 io = (option *)Make_Input();
268 io->fp_in_tree = fp_tree_file;
269 /* Read_Tree_File(io); */
270 tree = io->treelist->tree[0];
271 tree->io = io;
272
273 tree->io->z_scores = (phydbl *)mCalloc(tree->n_otu,sizeof(phydbl));
274
275 for(i=0;i<tree->n_otu;i++) tree->io->z_scores[i] = TIPO_Read_One_Taxon_Zscore(fp_coord_file,tree->a_nodes[i]->name,1,tree);
276 /* TIPO_Normalize_Zscores(tree); */
277 Free_Bip(tree);
278 Alloc_Bip(tree);
279 Get_Bip(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
280 TIPO_Get_Tips_Y_Rank_From_Zscores(tree);
281 /* TIPO_Get_Tips_Y_Rank(tree); */
282
283
284 tree->geo_mig_sd = 0.1;
285 Generic_Brent_Lk(&(tree->geo_mig_sd),
286 1.E-3,1.E+2,1.E-6,
287 100,NO,
288 &Wrap_Geo_Lk,
289 NULL,tree,NULL,NO);
290
291 PhyML_Printf("\n. sd=%f",tree->geo_mig_sd);
292
293
294 fclose(fp_tree_file);
295 fclose(fp_coord_file);
296
297 return 0;
298 }
299
300 //////////////////////////////////////////////////////////////
301 //////////////////////////////////////////////////////////////
302
303 /* Z_scores have already been recorder here */
TIPO_Get_Min_Number_Of_Tip_Permut(t_tree * tree)304 void TIPO_Get_Min_Number_Of_Tip_Permut(t_tree *tree)
305 {
306 Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
307 Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
308
309 Free_Bip(tree);
310 Alloc_Bip(tree);
311 Get_Bip(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
312
313 TIPO_Get_Tips_Y_Rank_From_Zscores(tree);
314
315 TIPO_Untangle_Tree(tree);
316
317 }
318
319 //////////////////////////////////////////////////////////////
320 //////////////////////////////////////////////////////////////
321
322
TIPO_Get_Tips_Y_Rank(t_tree * tree)323 void TIPO_Get_Tips_Y_Rank(t_tree *tree)
324 {
325 phydbl curr_rank;
326
327 curr_rank = .0;
328 TIPO_Get_Tips_Y_Rank_Pre(tree->n_root,tree->n_root->v[2],&curr_rank,tree);
329 TIPO_Get_Tips_Y_Rank_Pre(tree->n_root,tree->n_root->v[1],&curr_rank,tree);
330
331 if(curr_rank != tree->n_otu)
332 {
333 PhyML_Printf("\n. tree->n_otu = %d curr_rank = %d",tree->n_otu,curr_rank);
334 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
335 Warn_And_Exit("");
336 }
337 }
338
339 //////////////////////////////////////////////////////////////
340 //////////////////////////////////////////////////////////////
341
342
TIPO_Get_Tips_Y_Rank_Pre(t_node * a,t_node * d,phydbl * curr_rank,t_tree * tree)343 void TIPO_Get_Tips_Y_Rank_Pre(t_node *a, t_node *d, phydbl *curr_rank, t_tree *tree)
344 {
345 if(d->tax)
346 {
347 d->y_rank = *curr_rank;
348 *curr_rank += 1.;
349 return;
350 }
351 else
352 {
353 int i;
354 for(i=0;i<3;i++)
355 {
356 if(d->v[i] != a && d->b[i] != tree->e_root)
357 {
358 TIPO_Get_Tips_Y_Rank_Pre(d,d->v[i],curr_rank,tree);
359 }
360 }
361 }
362 }
363
364 //////////////////////////////////////////////////////////////
365 //////////////////////////////////////////////////////////////
366
367
TIPO_Get_All_Y_Rank(t_tree * tree)368 void TIPO_Get_All_Y_Rank(t_tree *tree)
369 {
370 tree->sum_y_dist_sq = .0;
371 tree->sum_y_dist = .0;
372 TIPO_Get_All_Y_Rank_Pre(tree->n_root,tree->n_root->v[2],tree);
373 TIPO_Get_All_Y_Rank_Pre(tree->n_root,tree->n_root->v[1],tree);
374 tree->n_root->y_rank = (tree->n_root->v[2]->y_rank+tree->n_root->v[1]->y_rank)/2.;
375 tree->n_root->y_rank_min = MIN(tree->n_root->v[2]->y_rank_min,tree->n_root->v[1]->y_rank_min);
376 tree->n_root->y_rank_max = MAX(tree->n_root->v[2]->y_rank_max,tree->n_root->v[1]->y_rank_max);
377 }
378
379 //////////////////////////////////////////////////////////////
380 //////////////////////////////////////////////////////////////
381
382
TIPO_Get_All_Y_Rank_Pre(t_node * a,t_node * d,t_tree * tree)383 void TIPO_Get_All_Y_Rank_Pre(t_node *a, t_node *d, t_tree *tree)
384 {
385 if(d->tax)
386 {
387 d->y_rank_min = d->y_rank;
388 d->y_rank_max = d->y_rank;
389 return;
390 }
391 else
392 {
393 int i;
394 int dir1,dir2;
395 phydbl v1,v2;
396
397 for(i=0;i<3;i++)
398 {
399 if(d->v[i] != a && d->b[i] != tree->e_root)
400 {
401 TIPO_Get_All_Y_Rank_Pre(d,d->v[i],tree);
402 }
403 }
404
405 dir1 = dir2 = -1;
406 for(i=0;i<3;i++)
407 {
408 if(d->v[i] != a && d->b[i] != tree->e_root)
409 {
410 if(dir1 < 0) dir1 = i;
411 else dir2 = i;
412 }
413 }
414
415 v1 = d->v[dir1]->y_rank;
416 v2 = d->v[dir2]->y_rank;
417
418 d->y_rank = (v1+v2)/2.;
419 tree->sum_y_dist_sq += POW(v1-v2,2);
420 tree->sum_y_dist += FABS(v1-v2);
421 d->y_rank_min = MIN(d->v[dir1]->y_rank_min,d->v[dir2]->y_rank_min);
422 d->y_rank_max = MAX(d->v[dir1]->y_rank_max,d->v[dir2]->y_rank_max);
423 }
424 }
425
426 //////////////////////////////////////////////////////////////
427 //////////////////////////////////////////////////////////////
428
429
TIPO_Swap_One_Node(t_node * d,t_tree * tree)430 void TIPO_Swap_One_Node(t_node *d, t_tree *tree)
431 {
432 if(d->tax) return;
433 else
434 {
435 int i;
436 int dir1, dir2;
437 t_node *tmp_n;
438 t_edge *tmp_e;
439
440 if(d != tree->n_root)
441 {
442 dir1 = dir2 = -1;
443 for(i=0;i<3;i++)
444 {
445 if((d->v[i] != d->anc) && (d->b[i] != tree->e_root))
446 {
447 if(dir1 < 0) dir1 = i;
448 else dir2 = i;
449 }
450 }
451 }
452 else
453 {
454 dir1 = 0;
455 dir2 = 1;
456 }
457
458 tmp_n = d->v[dir2];
459 d->v[dir2] = d->v[dir1];
460 d->v[dir1] = tmp_n;
461
462 tmp_e = d->b[dir2];
463 d->b[dir2] = d->b[dir1];
464 d->b[dir1] = tmp_e;
465 }
466 }
467
468 //////////////////////////////////////////////////////////////
469 //////////////////////////////////////////////////////////////
470
471
TIPO_Minimize_Tip_Order_Score(int n_trees,t_tree ** list_tree,t_tree * ref_tree)472 void TIPO_Minimize_Tip_Order_Score(int n_trees, t_tree **list_tree, t_tree *ref_tree)
473 {
474 int i,j;
475 phydbl score,min_score,old_min_score;
476 phydbl diff,eps;
477 t_node **node_table;
478 int swapped;
479 t_node *tmp;
480
481 eps = 1.E-3;
482 old_min_score = min_score = score = INT_MAX;
483 /* TIPO_Print_Tip_Ordered(ref_tree); */
484
485 do
486 {
487 for(i=ref_tree->n_otu;i<2*ref_tree->n_otu-1;i++)
488 {
489 TIPO_Swap_One_Node(ref_tree->a_nodes[i],ref_tree);
490 TIPO_Get_Tips_Y_Rank(ref_tree);
491 score = (phydbl)TIPO_Untangle_Tree_List(n_trees,list_tree,ref_tree);
492 if(score == -1)
493 {
494 return;
495 }
496
497 if(score < min_score)
498 {
499 min_score = score;
500 PhyML_Printf("\n- Score = %f",score);
501 }
502 else
503 {
504 TIPO_Swap_One_Node(ref_tree->a_nodes[i],ref_tree);
505 TIPO_Get_Tips_Y_Rank(ref_tree);
506 /* PhyML_Printf("\n+ Score = %f",score); */
507 }
508 }
509 diff = fabs(old_min_score - min_score);
510 old_min_score = min_score;
511 }while(diff > eps);
512
513 PhyML_Printf("\n");
514
515 node_table = (t_node **)mCalloc(ref_tree->n_otu,sizeof(t_node *));
516
517
518 for(i=0;i<ref_tree->n_otu;i++)
519 {
520 for(j=0;j<ref_tree->n_otu;j++)
521 {
522 if(!strcmp(ref_tree->io->short_tax_names[i],ref_tree->a_nodes[j]->name))
523 {
524 Free(ref_tree->a_nodes[j]->name);
525 ref_tree->a_nodes[j]->name = (char *)mCalloc((int)strlen(ref_tree->io->long_tax_names[i])+1,sizeof(char));
526 strcpy(ref_tree->a_nodes[j]->name,ref_tree->io->long_tax_names[i]);
527 break;
528 }
529 }
530 }
531
532 for(i=0;i<ref_tree->n_otu;i++) node_table[i] = ref_tree->a_nodes[i];
533
534 /* bubble sort of conflict nodes according to their y_rank */
535 do
536 {
537 swapped = NO;
538 for(i=0;i<ref_tree->n_otu-1;i++)
539 {
540 if(node_table[i]->y_rank > node_table[i+1]->y_rank)
541 {
542 swapped = YES;
543 tmp = node_table[i];
544 node_table[i] = node_table[i+1];
545 node_table[i+1] = tmp;
546 }
547 }
548 }while(swapped == YES);
549
550 for(i=0;i<ref_tree->n_otu;i++)
551 {
552 PhyML_Printf("\n%s",node_table[i]->name,node_table[i]->y_rank);
553 }
554
555
556
557 Free(node_table);
558
559 }
560
561 //////////////////////////////////////////////////////////////
562 //////////////////////////////////////////////////////////////
563
564
TIPO_Print_Tip_Ordered(t_tree * tree)565 void TIPO_Print_Tip_Ordered(t_tree *tree)
566 {
567 TIPO_Print_Tip_Ordered_Pre(tree->n_root,tree->n_root->v[2],tree);
568 TIPO_Print_Tip_Ordered_Pre(tree->n_root,tree->n_root->v[1],tree);
569 }
570
571 //////////////////////////////////////////////////////////////
572 //////////////////////////////////////////////////////////////
573
574
TIPO_Print_Tip_Ordered_Pre(t_node * a,t_node * d,t_tree * tree)575 void TIPO_Print_Tip_Ordered_Pre(t_node *a, t_node *d, t_tree *tree)
576 {
577
578 if(d->tax)
579 {
580 PhyML_Printf("\n. %f \"%s\"",d->y_rank,d->name);
581 }
582 else
583 {
584 int i,dir1,dir2;
585
586 dir1 = dir2 = -1;
587 for(i=0;i<3;i++)
588 {
589 if((d->v[i] != a) && (d->b[i] != tree->e_root))
590 {
591 if(dir1 < 0) dir1 = i;
592 else dir2 = i;
593 }
594 }
595 if(d->v[dir1]->y_rank < d->v[dir2]->y_rank)
596 {
597 TIPO_Print_Tip_Ordered_Pre(d,d->v[dir1],tree);
598 TIPO_Print_Tip_Ordered_Pre(d,d->v[dir2],tree);
599 }
600 else
601 {
602 TIPO_Print_Tip_Ordered_Pre(d,d->v[dir2],tree);
603 TIPO_Print_Tip_Ordered_Pre(d,d->v[dir1],tree);
604 }
605 }
606 }
607
608 //////////////////////////////////////////////////////////////
609 //////////////////////////////////////////////////////////////
610
611
TIPO_Untangle_Tree_List(int n_trees,t_tree ** list_tree,t_tree * ref_tree)612 int TIPO_Untangle_Tree_List(int n_trees, t_tree **list_tree, t_tree *ref_tree)
613 {
614 int i,j;
615 int tree_score,score;
616
617 score = 0;
618 for(i=0;i<n_trees;i++)
619 {
620 /* PhyML_Printf("\n. Untangling tree %3d",i); */
621 for(j=0;j<ref_tree->n_otu;j++) list_tree[i]->a_nodes[j]->y_rank = list_tree[i]->a_nodes[j]->ext_node->y_rank;
622 tree_score = TIPO_Untangle_Tree(list_tree[i]);
623 /* PhyML_Printf(" score = %3d",tree_score); */
624 score += tree_score;
625 if(tree_score < 0)
626 {
627 return -1;
628 }
629 }
630
631 return score;
632 }
633
634 //////////////////////////////////////////////////////////////
635 //////////////////////////////////////////////////////////////
636
637
TIPO_Untangle_Tree(t_tree * tree)638 phydbl TIPO_Untangle_Tree(t_tree *tree)
639 {
640 int conflict;
641 int n_trials;
642 t_node **node_table;
643 int i,swapped;
644 t_node *tmp;
645
646 node_table = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *));
647
648 for(i=0;i<tree->n_otu;i++) node_table[i] = tree->a_nodes[i];
649 for(i=0;i<tree->n_otu;i++) tree->a_nodes[i]->y_rank_ori = tree->a_nodes[i]->y_rank;
650
651
652 /* bubble sort of nodes according to their y_rank */
653 do
654 {
655 swapped = NO;
656 for(i=0;i<tree->n_otu-1;i++)
657 {
658 if(node_table[i]->y_rank > node_table[i+1]->y_rank)
659 {
660 swapped = YES;
661 tmp = node_table[i];
662 node_table[i] = node_table[i+1];
663 node_table[i+1] = tmp;
664 }
665 }
666 }
667 while(swapped == YES);
668
669
670 /* Work out the y_rank values for every internal node given the external node ranks */
671 TIPO_Get_All_Y_Rank(tree);
672
673 tree->tip_order_score = .0;
674
675 n_trials = 0;
676 do
677 {
678 conflict= NO;
679 /* Recusrssive untangling of the tree */
680 TIPO_Untangle_Node(tree->n_root,tree->n_root->v[2],node_table,&conflict,tree);
681 TIPO_Untangle_Node(tree->n_root,tree->n_root->v[1],node_table,&conflict,tree);
682 n_trials++;
683
684
685 if(n_trials > 2) /* We should have been able to untangle the tree after just one tree traversal */
686 {
687 int i;
688 FILE *ps_tree;
689
690 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
691 ps_tree = (FILE *)fopen("failed_tree.ps","w");
692
693 Test_Node_Table_Consistency(tree);
694 tree->rates = RATES_Make_Rate_Struct(tree->n_otu);
695 RATES_Init_Rate_Struct(tree->rates,tree->io->rates,tree->n_otu);
696 tree->times = TIMES_Make_Time_Struct(tree->n_otu);
697 TIMES_Init_Time_Struct(tree->times,tree->io->times,tree->n_otu);
698 TIMES_Least_Square_Node_Times(tree->e_root,tree);
699 TIMES_Adjust_Node_Times(tree);
700 RATES_Update_Cur_Bl(tree);
701
702 DR_Print_Postscript_Header(1,ps_tree);
703 tree->ps_tree = DR_Make_Tdraw_Struct(tree);
704 DR_Init_Tdraw_Struct(tree->ps_tree);
705 DR_Get_Tree_Box_Width(tree->ps_tree,tree);
706 Dist_To_Root(tree);
707 tree->ps_tree->max_dist_to_root = DR_Get_Max_Dist_To_Root(tree);
708 for(i=0;i<tree->n_otu;i++) tree->ps_tree->ycoord[i] = tree->a_nodes[i]->y_rank * (int)(tree->ps_tree->page_height / (tree->n_otu));
709 DR_Get_X_Coord(NO,tree->ps_tree,tree);
710 DR_Get_Y_Coord(YES,tree->ps_tree,tree);
711 DR_Print_Tree_Postscript(1,NO,ps_tree,tree);
712 DR_Print_Postscript_EOF(ps_tree);
713 fclose(ps_tree);
714 Warn_And_Exit("");
715 }
716 }while(conflict == YES);
717
718 Free(node_table);
719 return tree->tip_order_score;
720 }
721
722 //////////////////////////////////////////////////////////////
723 //////////////////////////////////////////////////////////////
724
725
TIPO_Untangle_Node(t_node * a,t_node * d,t_node ** node_table,int * conflict,t_tree * tree)726 void TIPO_Untangle_Node(t_node *a, t_node *d, t_node **node_table, int *conflict, t_tree *tree)
727 {
728
729 if(d->tax) return;
730 else
731 {
732 int i,j;
733 int d_a;
734 int beg,end;
735 phydbl min,max;
736 t_node *lca;
737 t_node **conflict_tips, **anc_conflict;
738 int n_conflicts;
739 phydbl eps,tmp_rank;
740 t_node *tmp_node;
741 int n_moved;
742 int n_anc_conflicts;
743
744 anc_conflict = NULL;
745 d_a = -1;
746
747 /* It is a post order traversal */
748 for(i=0;i<3;i++)
749 {
750 if((d->v[i] != d->anc) && (d->b[i] != tree->e_root))
751 {
752 TIPO_Untangle_Node(d,d->v[i],node_table,conflict,tree);
753 }
754 }
755
756
757 lca = NULL;
758 eps = 1.E-10;
759
760 /* Find direction fron node d ((d)escendant) to a ((a)ncestor) */
761 for(i=0;i<3;i++)
762 {
763 if((d->v[i] == d->anc) || (d->b[i] == tree->e_root))
764 {
765 d_a = i;
766 break;
767 }
768 }
769
770
771 /* y_rank_min is the minimum rank among all the ranks of the tips that
772 can be reached when going from a to d */
773 min = d->y_rank_min;
774 max = d->y_rank_max;
775
776 /* Get the list of tip nodes which ranks are between d->y_rank_min and d->y_rank_max */
777 n_conflicts = 0;
778 for(i=0;i<tree->n_otu;i++)
779 {
780 if((node_table[i]->y_rank > min - eps) && (node_table[i]->y_rank < max + eps))
781 {
782 n_conflicts++;
783 }
784 }
785
786 conflict_tips = NULL;
787 for(i=0;i<tree->n_otu;i++)
788 {
789 if(node_table[i]->y_rank > min - eps)
790 {
791 conflict_tips = node_table+i;
792 break;
793 }
794 }
795
796
797 beg = 0;
798 end = n_conflicts;
799 n_moved = 0;
800 n_anc_conflicts = 0;
801 do
802 {
803 for(i=beg;i<end;i++)
804 {
805 For(j,d->bip_size[d_a]) if(conflict_tips[i] == d->bip_node[d_a][j]) break;
806
807 if(j == d->bip_size[d_a])
808 /* conflict_tips[i] does not belong to the list of descendant of node d. It is
809 therefore responsible for a conflict */
810 {
811 *conflict = YES;
812
813 /* printf("\n. Moving %d with rank %f",conflict_tips[i]->num,conflict_tips[i]->y_rank); */
814
815 n_moved++;
816
817 /* Move from conflict_tips[i] towards the root as long as the rank of the node lca is between min and max */
818 lca = conflict_tips[i];
819 while(lca->y_rank_min > min && lca->y_rank_max < max) lca = lca->anc;
820
821
822 if(lca->y_rank_min > max && lca->y_rank_max > max)
823 {
824 PhyML_Printf("\n. lca = %d (%d)",lca->num,lca==tree->n_root);
825 PhyML_Printf("\n. Err in file %s at line %d",__FILE__,__LINE__);
826 Warn_And_Exit("");
827 }
828
829 if(lca->y_rank_min < min && lca->y_rank_max < min)
830 {
831 PhyML_Printf("\n. lca = %d (root ? %d) (tip ? %d)",lca->num,lca==tree->n_root,lca->tax);
832 PhyML_Printf("\n. lca->y_rank_min = %f lca->y_rank_max = %f",lca->y_rank_min,lca->y_rank_max);
833 PhyML_Printf("\n. min=%f max=%f",min,max);
834 PhyML_Printf("\n. Err in file %s at line %d",__FILE__,__LINE__);
835 Warn_And_Exit("");
836 }
837 if(lca->tax)
838 {
839 PhyML_Printf("\n. lca (%d) cannot be a tip.",lca->num);
840 PhyML_Printf("\n. lca->anc->y_rank=%f",lca->anc->y_rank);
841 PhyML_Printf("\n. lca->y_rank=%f",lca->y_rank);
842 PhyML_Printf("\n. lca->y_rank_min=%f lca=>y_rank_max=%f",lca->y_rank_min,lca->y_rank_max);
843 PhyML_Printf("\n. min=%f max=%f",min,max);
844 PhyML_Printf("\n. %p %p",lca,conflict_tips[i]);
845 PhyML_Printf("\n. Err in file %s at line %d",__FILE__,__LINE__);
846 Warn_And_Exit("");
847 }
848
849 /* Have you found lca previously ? */
850 for(j=0;j<n_anc_conflicts;j++) if(anc_conflict[j] == lca) break;
851 if(j == n_anc_conflicts) /* if no, then update the tree score and the list of ancestral nodes at the origin of conflicts */
852 {
853 /* tree->tip_order_score+=1.; */
854 /* tree->tip_order_score+=(lca->y_rank_max-lca->y_rank_min); */
855 n_anc_conflicts++;
856 anc_conflict = (t_node **)realloc(anc_conflict,n_anc_conflicts*sizeof(t_node *));
857 anc_conflict[n_anc_conflicts-1] = lca;
858 }
859
860 /* PhyML_Printf("\n. Detected conflict for ``%s'' (rank:%f min=%f max=%f lca=%f)", */
861 /* conflict_tips[i]->name, */
862 /* conflict_tips[i]->y_rank, */
863 /* min,max,lca->y_rank); */
864
865 /* Solve the conflict by shifting tip nodes to the left or to the right */
866 if(lca->y_rank > d->y_rank)
867 {
868 end--;
869 /* max-=1.; */
870 max = conflict_tips[end-1]->y_rank;
871 /* printf("\n. max=%f %f",max,conflict_tips[end-1]->y_rank); */
872
873 for(j=i;j<n_conflicts-1;j++)
874 {
875 /* PhyML_Printf("\n+ Moved (%d,%d) from (%f,%f)", */
876 /* conflict_tips[j]->num,conflict_tips[j+1]->num, */
877 /* conflict_tips[j]->y_rank,conflict_tips[j+1]->y_rank); */
878
879 tmp_rank = conflict_tips[j]->y_rank;
880 conflict_tips[j]->y_rank = conflict_tips[j+1]->y_rank;
881 conflict_tips[j+1]->y_rank = tmp_rank;
882
883 tmp_node = conflict_tips[j];
884 conflict_tips[j] = conflict_tips[j+1];
885 conflict_tips[j+1] = tmp_node;
886
887 tree->tip_order_score += fabs(conflict_tips[j]->y_rank - conflict_tips[j+1]->y_rank)/n_conflicts;
888
889 /* PhyML_Printf(" to (%d,%d) (%f,%f)", */
890 /* conflict_tips[j]->num,conflict_tips[j+1]->num, */
891 /* conflict_tips[j]->y_rank,conflict_tips[j+1]->y_rank); */
892 }
893 }
894 else
895 {
896 beg++;
897 /* min+=1.; */
898 min = conflict_tips[beg]->y_rank;
899 /* printf("\n. min=%f %f",min,conflict_tips[beg]->y_rank); */
900
901 for(j=i;j>0;j--)
902 {
903 /* PhyML_Printf("\n- Moved (%d,%d) from (%f,%f)", */
904 /* conflict_tips[j]->num,conflict_tips[j-1]->num, */
905 /* conflict_tips[j]->y_rank,conflict_tips[j-1]->y_rank); */
906
907
908 tmp_rank = conflict_tips[j]->y_rank;
909 conflict_tips[j]->y_rank = conflict_tips[j-1]->y_rank;
910 conflict_tips[j-1]->y_rank = tmp_rank;
911
912 tmp_node = conflict_tips[j];
913 conflict_tips[j] = conflict_tips[j-1];
914 conflict_tips[j-1] = tmp_node;
915
916 tree->tip_order_score += fabs(conflict_tips[j]->y_rank - conflict_tips[j+1]->y_rank)/n_conflicts;
917
918 /* PhyML_Printf(" to (%d,%d) (%f,%f)", */
919 /* conflict_tips[j]->num,conflict_tips[j-1]->num, */
920 /* conflict_tips[j]->y_rank,conflict_tips[j-1]->y_rank); */
921 }
922 }
923
924 /* printf("\n"); */
925 /* for(j=0;j<n_conflicts;j++) printf("%.0f ",conflict_tips[j]->y_rank); */
926 /* printf("\n. min=%f max=%f",min,max); */
927 /* printf("\n. Node %d has now rank %f",c_node->num,c_node->y_rank); */
928
929
930 /* Update internal nodes y ranks */
931 TIPO_Get_All_Y_Rank(tree);
932
933 break;
934 }
935 }
936 }while(n_moved + d->bip_size[d_a] != n_conflicts);
937
938
939 for(i=0;i<tree->n_otu;i++)
940 {
941 if((tree->a_nodes[i]->y_rank > min - eps) && (tree->a_nodes[i]->y_rank < max + eps))
942 {
943 For(j,d->bip_size[d_a])
944 {
945 if(tree->a_nodes[i] == d->bip_node[d_a][j]) break;
946 }
947 if(j == d->bip_size[d_a])
948 {
949 printf("\n. Conflict remaining for node %d (%d)",d->num,a->num);
950 PhyML_Printf("\n. Err in file %s at line %d",__FILE__,__LINE__);
951 Warn_And_Exit("");
952 }
953 }
954 }
955
956 if(anc_conflict) Free(anc_conflict);
957
958 return;
959 }
960 }
961
962 //////////////////////////////////////////////////////////////
963 //////////////////////////////////////////////////////////////
964
965
TIPO_Check_Tip_Ranks(t_tree * tree)966 int TIPO_Check_Tip_Ranks(t_tree *tree)
967 {
968 int i,j;
969 phydbl eps;
970
971 eps = 1.E-6;
972
973 for(i=0;i<tree->n_otu-1;i++)
974 {
975 for(j=i+1;j<tree->n_otu;j++)
976 {
977 if(fabs(tree->a_nodes[i]->y_rank - tree->a_nodes[j]->y_rank) < eps)
978 {
979 return 0;
980 }
981 }
982 }
983 return 1;
984 }
985
986 //////////////////////////////////////////////////////////////
987 //////////////////////////////////////////////////////////////
988
989
TIPO_Read_Taxa_Zscores(FILE * fp_coord,t_tree * tree)990 void TIPO_Read_Taxa_Zscores(FILE *fp_coord, t_tree *tree)
991 {
992 char *name,*line;
993 phydbl z;
994 int i;
995
996 name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
997 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
998
999 if(!fgets(line,T_MAX_LINE,fp_coord))
1000 {
1001 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1002 Warn_And_Exit("");
1003 }
1004 Free(line);
1005
1006
1007 do
1008 {
1009 if(fscanf(fp_coord,"%s\t%lf\n",name,&z) == EOF) break;
1010 PhyML_Printf("\n. Read %s. Z-score: %f",name,z);
1011
1012 for(i=0;i<tree->n_otu;i++) if(!strcmp(tree->io->long_tax_names[i],name)) break;
1013
1014 if(i == tree->n_otu)
1015 {
1016 PhyML_Printf("\n. Could not find taxon '%s' in coordinate file.",name);
1017 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1018 Warn_And_Exit("");
1019 }
1020
1021 tree->io->z_scores[i] = z;
1022
1023 }while(1);
1024
1025 Free(name);
1026 }
1027
1028 //////////////////////////////////////////////////////////////
1029 //////////////////////////////////////////////////////////////
1030
1031
TIPO_Read_Taxa_Coordinates(FILE * fp_coord,t_tree * tree)1032 void TIPO_Read_Taxa_Coordinates(FILE *fp_coord, t_tree *tree)
1033 {
1034 char *name,*line;
1035 phydbl lon, lat;
1036 int i;
1037
1038 name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1039 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1040
1041 if(!fgets(line,T_MAX_LINE,fp_coord))
1042 {
1043 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1044 Warn_And_Exit("");
1045 }
1046 Free(line);
1047
1048 tree->io->lat = (phydbl *)mCalloc(tree->n_otu,sizeof(phydbl));
1049 tree->io->lon = (phydbl *)mCalloc(tree->n_otu,sizeof(phydbl));
1050
1051 do
1052 {
1053 if(fscanf(fp_coord,"%s\t%lf\t%lf\n",name,&lat,&lon) == EOF) break;
1054 PhyML_Printf("\n. Read %s %f %f",name,lat,lon);
1055
1056 for(i=0;i<tree->n_otu;i++) if(!strcmp(tree->io->long_tax_names[i],name)) break;
1057
1058 if(i == tree->n_otu)
1059 {
1060 PhyML_Printf("\n. Could not find taxon '%s' in coordinate file.",name);
1061 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1062 Warn_And_Exit("");
1063 }
1064
1065 tree->io->lat[i] = lat;
1066 tree->io->lon[i] = lon;
1067
1068 }while(1);
1069
1070 Free(name);
1071 }
1072
1073 //////////////////////////////////////////////////////////////
1074 //////////////////////////////////////////////////////////////
1075
1076
TIPO_Get_Tips_Y_Rank_From_Zscores(t_tree * tree)1077 void TIPO_Get_Tips_Y_Rank_From_Zscores(t_tree *tree)
1078 {
1079 int i;
1080
1081 for(i=0;i<tree->n_otu;i++) tree->a_nodes[i]->y_rank = .0;
1082
1083 /* Randomization in order to avoid ties */
1084 for(i=0;i<tree->n_otu;i++) tree->io->z_scores[i] += Rnorm(0.0,0.001);
1085
1086 for(i=0;i<tree->n_otu;i++) tree->a_nodes[i]->y_rank = tree->io->z_scores[i];
1087
1088 /* for(i=0;i<tree->n_otu;i++) tree->a_nodes[i]->y_rank = .0; */
1089
1090 /* for(i=0;i<tree->n_otu-1;i++) */
1091 /* { */
1092 /* for(j=i+1;j<tree->n_otu;j++) */
1093 /* { */
1094 /* if(tree->io->z_scores[i] > tree->io->z_scores[j]) */
1095 /* { */
1096 /* tree->a_nodes[i]->y_rank += 1.0; */
1097 /* } */
1098 /* else */
1099 /* if(tree->io->z_scores[i] < tree->io->z_scores[j]) */
1100 /* { */
1101 /* tree->a_nodes[j]->y_rank += 1.0; */
1102 /* } */
1103 /* else */
1104 /* { */
1105 /* PhyML_Printf("\n. Ties not allowed.\n"); */
1106 /* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
1107 /* Warn_And_Exit(""); */
1108 /* } */
1109 /* } */
1110 /* } */
1111
1112 /* for(i=0;i<tree->n_otu;i++) printf("- %f\n",tree->a_nodes[i]->y_rank); */
1113
1114 }
1115
1116 //////////////////////////////////////////////////////////////
1117 //////////////////////////////////////////////////////////////
1118
1119 /* Sort translation table such that tree->a_nodes[i]->name == tree->io->short_tax_name[i] for all i */
TIPO_Sort_Translation_Table(t_tree * tree)1120 void TIPO_Sort_Translation_Table(t_tree *tree)
1121 {
1122 int i,j;
1123 char *s;
1124
1125
1126 Test_Node_Table_Consistency(tree);
1127
1128 for(i=0;i<tree->n_otu-1;i++)
1129 {
1130 for(j=i+1;j<tree->n_otu;j++)
1131 {
1132 if(!strcmp(tree->a_nodes[i]->name,tree->io->short_tax_names[j]))
1133 {
1134 s = tree->io->short_tax_names[i];
1135 tree->io->short_tax_names[i] = tree->io->short_tax_names[j];
1136 tree->io->short_tax_names[j] = s;
1137
1138 s = tree->io->long_tax_names[i];
1139 tree->io->long_tax_names[i] = tree->io->long_tax_names[j];
1140 tree->io->long_tax_names[j] = s;
1141
1142 break;
1143 }
1144 }
1145 }
1146 }
1147
1148 //////////////////////////////////////////////////////////////
1149 //////////////////////////////////////////////////////////////
1150
1151
TIPO_Randomize_Tip_Y_Ranks(t_tree * tree)1152 void TIPO_Randomize_Tip_Y_Ranks(t_tree *tree)
1153 {
1154
1155 int i;
1156 phydbl rk_tmp;
1157 int rnd_node_num;
1158
1159 for(i=0;i<tree->n_otu;i++) tree->a_nodes[i]->y_rank_ori = tree->a_nodes[i]->y_rank;
1160
1161 for(i=0;i<tree->n_otu;i++)
1162 {
1163 rnd_node_num = Rand_Int(0,tree->n_otu-1);
1164
1165 rk_tmp = tree->a_nodes[rnd_node_num]->y_rank;
1166 tree->a_nodes[rnd_node_num]->y_rank = tree->a_nodes[i]->y_rank;
1167 tree->a_nodes[i]->y_rank = rk_tmp;
1168 }
1169 }
1170
1171 //////////////////////////////////////////////////////////////
1172 //////////////////////////////////////////////////////////////
1173
1174
TIPO_Read_One_Taxon_Zscore(FILE * fp_coord,char * seqname_qry,int col,t_tree * tree)1175 phydbl TIPO_Read_One_Taxon_Zscore(FILE *fp_coord, char *seqname_qry, int col, t_tree *tree)
1176 {
1177 char *seqname, *place;
1178 phydbl lat;
1179
1180 seqname = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1181 place = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1182
1183 rewind(fp_coord);
1184
1185 /* skip first line */
1186 /* if(!fgets(line,T_MAX_LINE,fp_coord)) */
1187 /* { */
1188 /* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
1189 /* Warn_And_Exit(""); */
1190 /* } */
1191 /* Free(line); */
1192
1193 do
1194 {
1195 /* if(fscanf(fp_coord,"%s\t%s\t%lf\t%lf\t%d\n",seqname,place,&lat,&lon,&year) == EOF) */
1196 /* if(fscanf(fp_coord,"%s\t%s\t%lf\t%lf\n",seqname,place,&lat,&lon) == EOF) */
1197 /* if(fscanf(fp_coord,"%s\t%lf\t%lf\n",seqname,&lat,&lon) == EOF) */
1198 if(fscanf(fp_coord,"%s %lf\n",seqname,&lat) == EOF)
1199 {
1200 PhyML_Printf("\n. Could not find sequence '%s' in coordinate file",seqname_qry);
1201 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1202 Warn_And_Exit("");
1203 }
1204
1205
1206 if(!strcmp(seqname,seqname_qry)) break;
1207
1208 }while(1);
1209
1210 /* PhyML_Printf("\n. Found %20s %s @ %10.2f %10.2f. Recording %10.2f",seqname,place,lat,lon,(col==1)?lat:lon); */
1211
1212 Free(seqname);
1213 Free(place);
1214
1215 /* if(col == 1) return lat; */
1216 /* else if(col == 2) return lon; */
1217 /* else return -1.; */
1218
1219 return lat;
1220 }
1221
1222 //////////////////////////////////////////////////////////////
1223 //////////////////////////////////////////////////////////////
1224
1225
TIPO_Normalize_Zscores(t_tree * tree)1226 void TIPO_Normalize_Zscores(t_tree *tree)
1227 {
1228 int i;
1229 phydbl min_z,max_z;
1230 phydbl eps;
1231
1232 eps = 1.E-10;
1233
1234 min_z = FLT_MAX;
1235 for(i=0;i<tree->n_otu;i++)
1236 {
1237 if(tree->io->z_scores[i] < min_z)
1238 {
1239 min_z = tree->io->z_scores[i];
1240 }
1241 }
1242
1243 max_z = -FLT_MAX;
1244 for(i=0;i<tree->n_otu;i++)
1245 {
1246 if(tree->io->z_scores[i] > max_z)
1247 {
1248 max_z = tree->io->z_scores[i];
1249 }
1250 }
1251
1252 for(i=0;i<tree->n_otu;i++) tree->io->z_scores[i] = (tree->io->z_scores[i] - min_z)/(max_z-min_z+eps);
1253
1254 }
1255
1256 //////////////////////////////////////////////////////////////
1257 //////////////////////////////////////////////////////////////
1258
1259
TIPO_Lk(t_tree * tree)1260 phydbl TIPO_Lk(t_tree *tree)
1261 {
1262 tree->geo_lnL = 0.0;
1263 TIPO_Lk_Post(tree->n_root,tree->n_root->v[2],tree);
1264 TIPO_Lk_Post(tree->n_root,tree->n_root->v[1],tree);
1265 TIPO_Lk_Core(NULL,tree->n_root,tree);
1266 return(tree->geo_lnL);
1267 }
1268
1269 //////////////////////////////////////////////////////////////
1270 //////////////////////////////////////////////////////////////
1271
1272
TIPO_Lk_Post(t_node * a,t_node * d,t_tree * tree)1273 phydbl TIPO_Lk_Post(t_node *a, t_node *d, t_tree *tree)
1274 {
1275 if(!d->tax)
1276 {
1277 int i;
1278
1279 for(i=0;i<3;i++)
1280 {
1281 if(d->v[i] != a && d->b[i] != tree->e_root)
1282 {
1283 TIPO_Lk_Post(d,d->v[i],tree);
1284 }
1285 }
1286 TIPO_Lk_Core(a,d,tree);
1287 }
1288
1289 return .0;
1290 }
1291
1292 //////////////////////////////////////////////////////////////
1293 //////////////////////////////////////////////////////////////
1294
1295
TIPO_Lk_Core(t_node * a,t_node * d,t_tree * tree)1296 phydbl TIPO_Lk_Core(t_node *a, t_node *d, t_tree *tree)
1297 {
1298
1299 int i,j;
1300 int d_v1,d_v2,v1_d,v2_d;
1301 t_node *v1, *v2;
1302 phydbl dist,dens,min_dist;
1303
1304
1305 if(d->tax)
1306 {
1307 PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1308 Warn_And_Exit("");
1309 }
1310
1311 if(d == tree->n_root)
1312 {
1313 d_v1 = 0;
1314 d_v2 = 1;
1315 }
1316 else
1317 {
1318 d_v1 = d_v2 = -1;
1319 for(i=0;i<3;i++)
1320 {
1321 if(d->v[i] != a && d->b[i] != tree->e_root)
1322 {
1323 if(d_v1 < 0) d_v1 = i;
1324 else d_v2 = i;
1325 }
1326 }
1327 }
1328
1329 v1 = d->v[d_v1];
1330 v2 = d->v[d_v2];
1331
1332 v1_d = v2_d = -1;
1333 if(d == tree->n_root)
1334 {
1335 for(i=0;i<3;i++)
1336 {
1337 if(v1->b[i] == tree->e_root) v1_d = i;
1338 if(v2->b[i] == tree->e_root) v2_d = i;
1339 }
1340 }
1341 else
1342 {
1343 for(i=0;i<3;i++)
1344 {
1345 if(v1->v[i] == d) v1_d = i;
1346 if(v2->v[i] == d) v2_d = i;
1347 }
1348 }
1349
1350
1351 dens = 0.0;
1352 min_dist = FLT_MAX;
1353 For(i,v1->bip_size[v1_d])
1354 {
1355 For(j,v2->bip_size[v2_d])
1356 {
1357 dist = fabs(v1->bip_node[v1_d][i]->y_rank -
1358 v2->bip_node[v2_d][j]->y_rank);
1359
1360 if(dist < min_dist) min_dist = dist;
1361
1362 dens += Dnorm(dist,0.0,tree->geo_mig_sd);
1363 /* printf("\n. dist=%f dens=%f %f %f", */
1364 /* dist,Dnorm(dist,0.0,tree->geo_mig_sd), */
1365 /* v1->bip_node[v1_d][i]->y_rank, */
1366 /* v2->bip_node[v2_d][j]->y_rank); */
1367 }
1368 }
1369
1370 /* printf("\n. min_dist=%f dens=%f", */
1371 /* min_dist,Dnorm(dist,0.0,tree->geo_mig_sd)); */
1372
1373 /* dens = Dnorm(min_dist,0.0,tree->geo_mig_sd); */
1374 tree->geo_lnL += log(dens);
1375
1376 return tree->geo_lnL;
1377 }
1378
1379 //////////////////////////////////////////////////////////////
1380 //////////////////////////////////////////////////////////////
1381
1382 //////////////////////////////////////////////////////////////
1383 //////////////////////////////////////////////////////////////
1384
1385 //////////////////////////////////////////////////////////////
1386 //////////////////////////////////////////////////////////////
1387
1388 //////////////////////////////////////////////////////////////
1389 //////////////////////////////////////////////////////////////
1390
1391