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 #include "io.h"
14 #include "assert.h"
15
16 #ifdef BEAGLE
17 #include "libhmsbeagle/beagle.h"
18 #endif
19
20
21 /* Tree parser function. We need to pass a pointer to the string of characters
22 since this string might be freed and then re-allocated by that function (i.e.,
23 its address in memory might change)
24 */
Read_Tree(char ** s_tree)25 t_tree *Read_Tree(char **s_tree)
26 {
27 char **subs;
28 int i,n_ext,n_int,n_otu;
29 t_tree *tree;
30 int degree,len;
31 t_node *root_node;
32
33 n_int = n_ext = 0;
34
35 n_otu=0;
36 for(i=0;i<(int)strlen((*s_tree));++i)
37 {
38 if((*s_tree)[i] == '[') do ++i; while((*s_tree)[i] != ']'); // Skip labels
39 if((*s_tree)[i] == ',') n_otu++;
40 }
41 n_otu+=1;
42
43 tree = Make_Tree_From_Scratch(n_otu,NULL);
44 subs = Sub_Trees((*s_tree),°ree);
45 Clean_Multifurcation(subs,degree,3);
46
47 if(degree == 2)
48 {
49 /* Unroot_Tree(subs); */
50 /* degree = 3; */
51 /* root_node = tree->a_nodes[n_otu]; */
52
53 root_node = tree->a_nodes[2*n_otu-2];
54 root_node->num = 2*n_otu-2;
55 tree->n_root = root_node;
56 n_int -= 1;
57 }
58 else
59 {
60 root_node = tree->a_nodes[n_otu];
61 root_node->num = n_otu;
62 tree->n_root = NULL;
63 }
64
65 if(degree > 3) /* Multifurcation at the root. Need to re-assemble the subtrees
66 since Clean_Multifurcation added sets of parenthesis and
67 the corresponding NULL edges */
68 {
69 degree = 3;
70 Free((*s_tree));
71 len = 0;
72 for(i=0;i<degree;i++) len += (strlen(subs[i])+1);
73 len += 5;
74
75 (*s_tree) = (char *)mCalloc(len,sizeof(char));
76
77 (*s_tree)[0] = '('; (*s_tree)[1] = '\0';
78 for(i=0;i<degree;i++)
79 {
80 strcat((*s_tree),subs[i]);
81 strcat((*s_tree),",\0");
82 }
83
84 sprintf((*s_tree)+strlen((*s_tree))-1,"%s",");\0");
85
86 i = 0;
87 while(subs[i] != NULL) Free(subs[i++]);
88 Free(subs);
89
90 subs = Sub_Trees((*s_tree),°ree);
91 }
92
93 root_node->tax = 0;
94
95 tree->has_branch_lengths = 0;
96 tree->num_curr_branch_available = tree->n_otu;
97 for(i=0;i<degree;i++) R_rtree((*s_tree),subs[i],root_node,tree,&n_int,&n_ext);
98
99 i = degree;
100 while(subs[i] != NULL) Free(subs[i++]);
101 Free(subs);
102
103
104 if(tree->n_root)
105 {
106
107 if(tree->n_root->v[1]->tax == NO && tree->n_root->v[2]->tax == NO)
108 {
109 tree->e_root = tree->a_edges[tree->num_curr_branch_available];
110 tree->n_root->b[1] = tree->a_edges[tree->num_curr_branch_available+1];
111 tree->n_root->b[2] = tree->a_edges[tree->num_curr_branch_available+2];
112 }
113 else
114 {
115 for(i=0;i<tree->n_otu;++i) if(tree->a_edges[i]->left == NULL && tree->a_edges[i]->rght == NULL) break;
116 assert(i != tree->n_otu);
117 tree->e_root = tree->a_edges[i];
118
119 tree->n_root->b[1] = tree->a_edges[tree->num_curr_branch_available];
120 tree->n_root->b[2] = tree->a_edges[tree->num_curr_branch_available+1];
121 }
122
123 tree->n_root->v[2]->v[0] = tree->n_root->v[1];
124 tree->n_root->v[1]->v[0] = tree->n_root->v[2];
125
126 // Reading edge lengths
127 subs = Sub_Trees((*s_tree),°ree);
128 Read_Branch_Length(subs[0],(*s_tree),tree->n_root->b[1],tree);
129 Read_Branch_Length(subs[1],(*s_tree),tree->n_root->b[2],tree);
130 Free(subs);
131
132 Connect_One_Edge_To_Two_Nodes(tree->n_root->v[2],
133 tree->n_root->v[1],
134 tree->e_root,
135 tree);
136
137 tree->e_root->l->v = tree->n_root->b[2]->l->v + tree->n_root->b[1]->l->v;
138 if(tree->e_root->l->v > 0.0)
139 tree->n_root_pos = tree->n_root->b[2]->l->v / tree->e_root->l->v;
140 else
141 tree->n_root_pos = .5;
142
143
144 Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
145 Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
146
147 }
148
149 return tree;
150 }
151
152 //////////////////////////////////////////////////////////////
153 //////////////////////////////////////////////////////////////
154
155 /* 'a' in t_node a stands for ancestor. 'd' stands for descendant */
R_rtree(char * s_tree_a,char * s_tree_d,t_node * a,t_tree * tree,int * n_int,int * n_ext)156 void R_rtree(char *s_tree_a, char *s_tree_d, t_node *a, t_tree *tree, int *n_int, int *n_ext)
157 {
158 int i;
159 t_node *d;
160 int n_otu = tree->n_otu;
161
162
163 if(strstr(s_tree_a," "))
164 {
165 PhyML_Fprintf(stderr,"\n. [%s]",s_tree_a);
166 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
167 }
168
169 // Internal edge
170 if(s_tree_d[0] == '(')
171 {
172 char **subs;
173 int degree;
174 t_edge *b;
175
176 (*n_int)+=1;
177
178 if((*n_int + n_otu) == (2*n_otu-1))
179 {
180 PhyML_Fprintf(stderr,"\n. The number of internal nodes in the tree exceeds the number of taxa minus one.");
181 PhyML_Fprintf(stderr,"\n. There probably is a formating problem in the input tree.");
182 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
183 }
184
185 d = tree->a_nodes[n_otu+*n_int];
186 d->num = n_otu + *n_int;
187 d->tax = 0;
188 b = tree->a_edges[tree->num_curr_branch_available];
189
190 Read_Branch_Support(s_tree_d,s_tree_a,b,tree);
191 Read_Branch_Length(s_tree_d,s_tree_a,b,tree);
192 Read_Node_Label(s_tree_d,s_tree_a,d);
193
194 if(tree->n_root && a == tree->n_root)
195 {
196 if(!a->v[1]) a->v[1] = d;
197 else a->v[2] = d;
198 }
199 else
200 {
201 for(i=0;i<3;i++)
202 {
203 if(!a->v[i])
204 {
205 a->v[i] = d;
206 break;
207 }
208 }
209 }
210 d->v[0]=a;
211
212 if(!(tree->n_root && a == tree->n_root)) Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[tree->num_curr_branch_available],tree);
213
214 subs=Sub_Trees(s_tree_d,°ree);
215
216
217 if(degree < 2)
218 {
219 PhyML_Fprintf(stderr,"\n. A problem was detected in the following subtree:");
220 PhyML_Fprintf(stderr,"\n. %s",s_tree_d);
221 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
222 }
223
224 if(degree > 2)
225 {
226 Clean_Multifurcation(subs,degree,2);
227
228 Free(s_tree_d);
229
230 s_tree_d = (char *)mCalloc(strlen(subs[0])+strlen(subs[1])+5,sizeof(char));
231
232 i = 0;
233 while(subs[i] != NULL) Free(subs[i++]);
234 Free(subs);
235
236 subs=Sub_Trees(s_tree_d,°ree);
237 }
238
239 R_rtree(s_tree_d,subs[0],d,tree,n_int,n_ext);
240 R_rtree(s_tree_d,subs[1],d,tree,n_int,n_ext);
241
242 i = 2;
243 while(subs[i] != NULL) Free(subs[i++]);
244 Free(subs);
245 }
246
247 // External edge
248 else
249 {
250 int i;
251
252 d = tree->a_nodes[*n_ext];
253 d->tax = 1;
254
255 Read_Node_Name(d,s_tree_d,tree->a_edges[*n_ext],tree);
256 Read_Branch_Length(s_tree_d,s_tree_a,tree->a_edges[*n_ext],tree);
257 Read_Node_Label(s_tree_d,s_tree_a,d);
258
259 if(tree->n_root && a == tree->n_root)
260 {
261 if(!a->v[1]) a->v[1] = d;
262 else a->v[2] = d;
263 }
264 else
265 {
266 for(i=0;i<3;i++)
267 {
268 if(!a->v[i])
269 {
270 a->v[i]=d;
271 break;
272 }
273 }
274 }
275 d->v[0]=a;
276
277 if(!(tree->n_root && a == tree->n_root))
278 {
279 Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[*n_ext],tree);
280 }
281
282 d->num=*n_ext;
283 (*n_ext)+=1;
284 }
285
286 Free(s_tree_d);
287
288 }
289
290 //////////////////////////////////////////////////////////////
291 //////////////////////////////////////////////////////////////
292
Read_Node_Name(t_node * d,char * s_tree_d,t_edge * b,t_tree * tree)293 void Read_Node_Name(t_node *d, char *s_tree_d, t_edge *b, t_tree *tree)
294 {
295 d->name = (char *)mCalloc(strlen(s_tree_d)+1,sizeof(char ));
296 strcpy(d->name,s_tree_d);
297 d->ori_name = d->name;
298 }
299
300 //////////////////////////////////////////////////////////////
301 //////////////////////////////////////////////////////////////
302
Read_Branch_Support(char * s_d,char * s_a,t_edge * b,t_tree * tree)303 void Read_Branch_Support(char *s_d, char *s_a, t_edge *b, t_tree *tree)
304 {
305 char *sub_tp;
306 char *p;
307
308 if(s_d[0] != '(') return; // b is an external edge -> no edge support on it
309
310 sub_tp = (char *)mCalloc(10+strlen(s_d)+1,sizeof(char));
311
312 sub_tp[0] = '(';
313 sub_tp[1] = '\0';
314 strcat(sub_tp,s_d);
315 p = strstr(s_a,sub_tp);
316
317 if(!p)
318 {
319 sub_tp[0] = ',';
320 sub_tp[1] = '\0';
321 strcat(sub_tp,s_d);
322 p = strstr(s_a,sub_tp);
323 }
324
325
326 if(p)
327 {
328 b->support_val = atof((char *)p+(int)strlen(sub_tp));
329
330 p = p + strlen(sub_tp);
331 if(p[0] == '[') { while(p[0] != ']') p++; p++; } // Skip label
332 b->support_val = atof((char *)p);
333
334 /* PhyML_Printf("\n. READ SUPPORT for s_d: %s b: %f",s_d,b->support_val); */
335
336 }
337
338 Free(sub_tp);
339 }
340
341 //////////////////////////////////////////////////////////////
342 //////////////////////////////////////////////////////////////
343
Read_Branch_Length(char * s_d,char * s_a,t_edge * b,t_tree * tree)344 void Read_Branch_Length(char *s_d, char *s_a, t_edge *b, t_tree *tree)
345 {
346 char *sub_tp;
347 char *p;
348
349 sub_tp = (char *)mCalloc(10+strlen(s_d)+1,sizeof(char));
350
351 sub_tp[0] = '(';
352 sub_tp[1] = '\0';
353 strcat(sub_tp,s_d);
354 p = strstr(s_a,sub_tp);
355
356 if(!p)
357 {
358 sub_tp[0] = ',';
359 sub_tp[1] = '\0';
360 strcat(sub_tp,s_d);
361 p = strstr(s_a,sub_tp);
362 }
363
364
365 if(p)
366 {
367 p = p + strlen(sub_tp);
368 while(p[0] != ':' && p[0] != '\0') p++;
369 if(p[0] == ':')
370 {
371 p++;
372 b->l->v = atof((char *)p);
373 /* PhyML_Printf("\n. READ LENGTH for s_d: %s b: %f",s_d,b->l->v); */
374 tree->has_branch_lengths = YES;
375 b->does_exist = YES;
376 }
377 else
378 {
379 b->l->v = -1.;
380 }
381 }
382 else
383 {
384 b->l->v = -1.;
385 }
386
387 Free(sub_tp);
388 }
389
390 //////////////////////////////////////////////////////////////
391 //////////////////////////////////////////////////////////////
392
Read_Node_Label(char * s_d,char * s_a,t_node * n)393 void Read_Node_Label(char *s_d, char *s_a, t_node *n)
394 {
395 char *sub_tp,*s_lab;
396 char *p;
397
398 sub_tp = (char *)mCalloc(3+(int)strlen(s_d)+1,sizeof(char));
399 s_lab = (char *)mCalloc(strlen(s_a)+1,sizeof(char));
400
401
402 sub_tp[0] = '(';
403 sub_tp[1] = '\0';
404 strcat(sub_tp,s_d);
405 strcat(sub_tp,"[\0");
406 p = strstr(s_a,sub_tp);
407
408 if(!p)
409 {
410 sub_tp[0] = ',';
411 sub_tp[1] = '\0';
412 strcat(sub_tp,s_d);
413 strcat(sub_tp,"[\0");
414 p = strstr(s_a,sub_tp);
415 }
416
417 if(p)
418 {
419 p = p+strlen(sub_tp)-1;
420 assert(p[0]=='[');
421
422 s_lab[0]='[';
423 if(sscanf(p,"[%[^]]]",s_lab+1) != 1)
424 {
425 PhyML_Fprintf(stderr,"\n. Label is in wrong format. A proper label should");
426 PhyML_Fprintf(stderr,"\n. look as follows: \"[xxx={yyy},xxxx={yy},...]\"");
427 assert(FALSE);
428 }
429
430 s_lab[strlen(s_lab)]=']';
431 s_lab[strlen(s_lab)]='\0';
432
433 /* PhyML_Printf("\n. READ LABEL for s_d: %s b: %s",s_d,s_lab); */
434
435 n->label = Read_Labels(s_lab);
436 }
437
438 Free(sub_tp);
439 Free(s_lab);
440 }
441
442 //////////////////////////////////////////////////////////////
443 //////////////////////////////////////////////////////////////
444
Clean_Multifurcation(char ** subtrees,int current_deg,int end_deg)445 void Clean_Multifurcation(char **subtrees, int current_deg, int end_deg)
446 {
447
448 if(current_deg <= end_deg) return;
449 else
450 {
451 char *s_tmp;
452 int i;
453
454 /* s_tmp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
455 s_tmp = (char *)mCalloc(10+
456 (int)strlen(subtrees[0])+1+
457 (int)strlen(subtrees[1])+1,
458 sizeof(char));
459
460 strcat(s_tmp,"(\0");
461 strcat(s_tmp,subtrees[0]);
462 strcat(s_tmp,",\0");
463 strcat(s_tmp,subtrees[1]);
464 strcat(s_tmp,")#NULL\0"); /* Add the label 'NULL' to identify a non-existing edge */
465 Free(subtrees[0]);
466 subtrees[0] = s_tmp;
467
468 for(i=1;i<current_deg-1;i++) strcpy(subtrees[i],subtrees[i+1]);
469
470 Clean_Multifurcation(subtrees,current_deg-1,end_deg);
471 }
472 }
473
474 //////////////////////////////////////////////////////////////
475 //////////////////////////////////////////////////////////////
476
Sub_Trees(char * tree,int * degree)477 char **Sub_Trees(char *tree, int *degree)
478 {
479 char **subs;
480 int posbeg,posend;
481 int i;
482
483 if(tree[0] != '(')
484 {
485 *degree = 1;
486 return NULL;
487 }
488
489 posbeg=posend=1;
490 (*degree)=0;
491 do
492 {
493 posbeg = posend;
494 if(tree[posend] != '(')
495 {
496 while((tree[posend] != ',' ) &&
497 (tree[posend] != ':' ) &&
498 (tree[posend] != '#' ) &&
499 (tree[posend] != ')' ) &&
500 (tree[posend] != '[' ))
501 {
502 posend++ ;
503 }
504 posend -= 1;
505 }
506 else posend = Next_Par(tree,posend);
507
508 if(*degree == 0)
509 subs = (char **)mCalloc(1,sizeof(char *));
510 else
511 subs = (char **)mRealloc(subs,*degree+1,sizeof(char *));
512
513 subs[(*degree)] = (char *)mCalloc(strlen(tree)+1,sizeof(char));
514 strncpy(subs[(*degree)],tree+posbeg,posend-posbeg+1);
515 subs[(*degree)][posend-posbeg+1]='\0'; /* Thanks to Jean-Baka Domelevo-Entfellner */
516
517 posend += 1;
518 if(tree[posend] == '[') { while(tree[posend] != ']') posend++; posend++; } /* Skip label */
519 while((tree[posend] != ',') && (tree[posend] != ')')) posend++;
520 posend++;
521
522 (*degree)++;
523 if((*degree) == NODE_DEG_MAX)
524 {
525 for(i=0;i<(*degree);++i) PhyML_Fprintf(stderr,"\n. Subtree %d : %s\n",i+1,subs[i]);
526 PhyML_Fprintf(stderr,"\n. The degree of a t_node cannot be greater than %d\n",NODE_DEG_MAX);
527 Warn_And_Exit("\n");
528 }
529 }
530 while(tree[posend-1] != ')');
531
532 subs = (char **)mRealloc(subs,*degree+1,sizeof(char *));
533 subs[(*degree)] = NULL;
534
535 return subs;
536 }
537
538 //////////////////////////////////////////////////////////////
539 //////////////////////////////////////////////////////////////
540
541
Next_Par(char * s,int pos)542 int Next_Par(char *s, int pos)
543 {
544 int curr;
545
546 curr=pos+1;
547
548 while(*(s+curr) != ')')
549 {
550 if(*(s+curr) == '(') curr=Next_Par(s,curr);
551 curr++;
552 }
553
554 return curr;
555 }
556
557 //////////////////////////////////////////////////////////////
558 //////////////////////////////////////////////////////////////
559
Print_Tree(FILE * fp,t_tree * tree)560 void Print_Tree(FILE *fp, t_tree *tree)
561 {
562 char *s_tree;
563 int i;
564
565 s_tree = (char *)Write_Tree(tree);
566
567 if(OUTPUT_TREE_FORMAT == NEWICK) PhyML_Fprintf(fp,"%s\n",s_tree);
568 else if(OUTPUT_TREE_FORMAT == NEXUS)
569 {
570 PhyML_Fprintf(fp,"#NEXUS\n");
571 PhyML_Fprintf(fp,"BEGIN TREES;\n");
572 PhyML_Fprintf(fp,"\tTRANSLATE\n");
573 for(i=0;i<tree->n_otu;++i) PhyML_Fprintf(fp,"\t%3d\t%s,\n",i+1,tree->a_nodes[i]->name);
574 PhyML_Fprintf(fp,"\tUTREE PAUP_1=\n");
575 PhyML_Fprintf(fp,"%s\n",s_tree);
576 PhyML_Fprintf(fp,"ENDBLOCK;");
577 }
578 Free(s_tree);
579 }
580
581 //////////////////////////////////////////////////////////////
582 //////////////////////////////////////////////////////////////
583
Write_Tree(t_tree * tree)584 char *Write_Tree(t_tree *tree)
585 {
586 char *s;
587 int i,available;
588
589 #ifndef MPI
590 int init_len;
591 init_len = 3*(int)T_MAX_NAME;
592 s=(char *)mCalloc(init_len,sizeof(char));
593 available = init_len;
594 #elif defined MPI
595 s=(char *)mCalloc(T_MAX_LINE,sizeof(char));
596 #endif
597
598 i = -1;
599 s[0]='(';
600
601 if(tree->n_root == NULL)
602 {
603 i = 0;
604 while((!tree->a_nodes[tree->n_otu+i]->v[0]) ||
605 (!tree->a_nodes[tree->n_otu+i]->v[1]) ||
606 (!tree->a_nodes[tree->n_otu+i]->v[2])) i++;
607
608 R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[0],tree->a_nodes[tree->n_otu+i]->b[0],&available,&s,tree);
609 R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[1],tree->a_nodes[tree->n_otu+i]->b[1],&available,&s,tree);
610 R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[2],tree->a_nodes[tree->n_otu+i]->b[2],&available,&s,tree);
611 }
612 else
613 {
614 R_wtree(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],&available,&s,tree);
615 R_wtree(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],&available,&s,tree);
616 }
617
618 s[(int)strlen(s)-1]=')';
619
620 if(tree->n_root != NULL) Print_Labels(NULL,s+(int)strlen(s),tree->n_root->label);
621
622 if(tree->io && tree->io->print_node_num == YES)
623 {
624 if(!tree->n_root)
625 sprintf(s+(int)strlen(s),"%d",tree->a_nodes[tree->n_otu+i]->num);
626 else
627 sprintf(s+(int)strlen(s),"%d",tree->n_root->num);
628 }
629 s[(int)strlen(s)]=';';
630
631 return s;
632 }
633
634 //////////////////////////////////////////////////////////////
635 //////////////////////////////////////////////////////////////
636
R_wtree(t_node * pere,t_node * fils,t_edge * b,int * available,char ** s_tree,t_tree * tree)637 void R_wtree(t_node *pere, t_node *fils, t_edge *b, int *available, char **s_tree, t_tree *tree)
638 {
639 int i,p;
640 char *format;
641 int last_len;
642 phydbl mean_len;
643
644 format = (char *)mCalloc(100,sizeof(char));
645
646 sprintf(format,"%%.%df",tree->bl_ndigits);
647
648 p = -1;
649 if(fils->tax == YES) // Tip node
650 {
651 last_len = (int)strlen(*s_tree);
652
653 if(OUTPUT_TREE_FORMAT == NEWICK)
654 {
655 if(tree->write_tax_names == YES)
656 {
657 if(tree->io && tree->io->long_tax_names)
658 {
659 strcat(*s_tree,tree->io->long_tax_names[fils->num]);
660 }
661 else
662 {
663 if(fils->name && strlen(fils->name) > 0)
664 strcat(*s_tree,fils->name);
665 else
666 sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
667 }
668 }
669 else if(tree->write_tax_names == NO)
670 {
671 sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
672 }
673 }
674 else if(OUTPUT_TREE_FORMAT == NEXUS)
675 {
676 sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
677 }
678 else
679 {
680 PhyML_Printf("\n. Unknown tree format.");
681 PhyML_Printf("\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
682 PhyML_Printf("\n. s=%s\n",*s_tree);
683 }
684
685 if((fils->b) && (fils->b[0]) && (tree->write_br_lens == YES))
686 {
687 Print_Labels(NULL,*s_tree+(int)strlen(*s_tree),fils->label);
688
689 (*s_tree)[(int)strlen(*s_tree)] = ':';
690
691 Print_Labels(NULL,*s_tree+(int)strlen(*s_tree),b->label);
692
693 if(tree->is_mixt_tree == NO) mean_len = b->l->v;
694 else mean_len = MIXT_Get_Mean_Edge_Len(b,tree);
695 sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
696 }
697
698 (*s_tree)[(int)strlen(*s_tree)] = ',';
699
700
701 #ifndef MPI
702 (*available) -= ((int)strlen(*s_tree) - last_len);
703
704 /* printf("\n0 Available = %d [%d %d]",(*available),(int)strlen(*s_tree),last_len); */
705 /* printf("\n0 %s [%d,%d]",*s_tree,(int)(int)strlen(*s_tree),*available); */
706
707 if(*available < 0)
708 {
709 PhyML_Fprintf(stderr,"\n. s=%s\n",*s_tree);
710 PhyML_Fprintf(stderr,"\n. len=%d\n",(int)strlen(*s_tree));
711 PhyML_Fprintf(stderr,"\n. The sequence names in your input file might be too long.");
712 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
713 Warn_And_Exit("");
714 }
715
716 if(*available < (int)T_MAX_NAME)
717 {
718 (*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
719 for(i=0;i<3*(int)T_MAX_NAME;++i) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
720 (*available) = 3*(int)T_MAX_NAME;
721 }
722 #endif
723
724 }
725 else // Internal node
726 {
727 (*s_tree)[(int)strlen(*s_tree)]='(';
728
729 #ifndef MPI
730 (*available) -= 1;
731
732
733 if(*available < (int)T_MAX_NAME)
734 {
735 (*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
736 for(i=0;i<3*(int)T_MAX_NAME;++i) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
737 (*available) = 3*(int)T_MAX_NAME;
738 }
739 #endif
740
741
742 for(i=0;i<3;i++)
743 {
744 if((fils->v[i] != pere) && (fils->b[i] != tree->e_root))
745 R_wtree(fils,fils->v[i],fils->b[i],available,s_tree,tree);
746 else p=i;
747 }
748
749 if(p < 0)
750 {
751 PhyML_Fprintf(stderr,"\n. pere: %d fils=%d root=%d root->v[2]=%d root->v[1]=%d",pere->num,fils->num,tree->n_root->num,tree->n_root->v[2]->num,tree->n_root->v[1]->num);
752 PhyML_Fprintf(stderr,"\n. fils=%d root=%d root->v[2]=%d root->v[1]=%d",fils->num,tree->n_root->num,tree->n_root->v[2]->num,tree->n_root->v[1]->num);
753 PhyML_Fprintf(stderr,"\n. tree->e_root=%d fils->b[0]=%d fils->b[1]=%d fils->b[2]=%d",tree->e_root->num,fils->b[0]->num,fils->b[1]->num,fils->b[2]->num);
754 assert(false);
755 }
756
757 last_len = (int)strlen(*s_tree);
758
759 (*s_tree)[last_len-1] = ')';
760
761
762 if((fils->b) && (tree->write_br_lens == YES))
763 {
764 Print_Labels(NULL,*s_tree+(int)strlen(*s_tree),fils->label);
765
766 if(tree->io && tree->io->print_support_val == YES)
767 {
768 if(tree->io->do_boot == YES)
769 {
770 sprintf(*s_tree+(int)strlen(*s_tree),"%.0f",fils->b[p]->support_val);
771 }
772 else
773 {
774 sprintf(*s_tree+(int)strlen(*s_tree),"%f",fils->b[p]->support_val);
775 }
776 }
777 if(tree->io && tree->io->print_node_num == YES)
778 {
779 sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num);
780 }
781
782
783 fflush(NULL);
784
785 (*s_tree)[(int)strlen(*s_tree)] = ':';
786
787 Print_Labels(NULL,*s_tree+(int)strlen(*s_tree),b->label);
788
789 if(tree->is_mixt_tree == NO) mean_len = b->l->v;
790 else mean_len = MIXT_Get_Mean_Edge_Len(b,tree);
791 sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
792 }
793
794 (*s_tree)[(int)strlen(*s_tree)] = ',';
795
796 #ifndef MPI
797 (*available) -= ((int)strlen(*s_tree) - last_len);
798
799
800 if(*available < 0)
801 {
802 PhyML_Fprintf(stderr,"\n. available = %d",*available);
803 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
804 Warn_And_Exit("");
805 }
806
807 if(*available < (int)T_MAX_NAME)
808 {
809 (*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
810 for(i=0;i<3*(int)T_MAX_NAME;++i) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
811 (*available) = 3*(int)T_MAX_NAME;
812 }
813 #endif
814
815 }
816
817 Free(format);
818 }
819
820 //////////////////////////////////////////////////////////////
821 //////////////////////////////////////////////////////////////
822
Detect_Align_File_Format(option * io)823 void Detect_Align_File_Format(option *io)
824 {
825 int c;
826 fpos_t curr_pos;
827
828 fgetpos(io->fp_in_align,&curr_pos);
829
830 errno = 0;
831
832 while((c=fgetc(io->fp_in_align)) != EOF)
833 {
834 if(errno) io->data_file_format = PHYLIP;
835 else if(c == '#')
836 {
837 char s[10],t[6]="NEXUS";
838 if(!fgets(s,6,io->fp_in_align))
839 {
840 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
841 Exit("\n");
842 }
843 if(!strcmp(t,s))
844 {
845 fsetpos(io->fp_in_align,&curr_pos);
846 io->data_file_format = NEXUS;
847 return;
848 }
849 }
850 }
851 fsetpos(io->fp_in_align,&curr_pos);
852 }
853
854 //////////////////////////////////////////////////////////////
855 //////////////////////////////////////////////////////////////
856
Detect_Tree_File_Format(option * io)857 void Detect_Tree_File_Format(option *io)
858 {
859 int c;
860 fpos_t curr_pos;
861
862 fgetpos(io->fp_in_tree,&curr_pos);
863
864 errno = 0;
865
866 while((c=fgetc(io->fp_in_tree)) != EOF)
867 {
868 if(errno)
869 {
870 io->tree_file_format = PHYLIP;
871 PhyML_Printf("\n. Detected PHYLIP tree file format.");
872 }
873 else if(c == '#')
874 {
875 char s[10],t[6]="NEXUS";
876 if(!fgets(s,6,io->fp_in_tree))
877 {
878 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
879 Warn_And_Exit("");
880 }
881 if(!strcmp(t,s))
882 {
883 fsetpos(io->fp_in_tree,&curr_pos);
884 io->tree_file_format = NEXUS;
885 PhyML_Printf("\n. Detected NEXUS tree file format.");
886 return;
887 }
888 }
889 }
890
891 fsetpos(io->fp_in_tree,&curr_pos);
892 }
893
894 //////////////////////////////////////////////////////////////
895 //////////////////////////////////////////////////////////////
896
Get_Seq(option * io)897 align **Get_Seq(option *io)
898 {
899 io->data = NULL;
900
901 if(!io->fp_in_align)
902 {
903 PhyML_Fprintf(stderr,"\n. Filehandle to '%s' seems to be closed.",io->in_align_file);
904 Exit("\n");
905 }
906
907 Detect_Align_File_Format(io);
908
909 switch(io->data_file_format)
910 {
911 case PHYLIP:
912 {
913 io->data = Get_Seq_Phylip(io);
914 break;
915 }
916 case NEXUS:
917 {
918 io->nex_com_list = Make_Nexus_Com();
919 Init_Nexus_Format(io->nex_com_list);
920 Get_Nexus_Data(io->fp_in_align,io);
921 Free_Nexus(io);
922 break;
923 }
924 default:
925 {
926 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
927 Exit("\n");
928 break;
929 }
930 }
931
932 if(!io->data)
933 {
934 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
935 Exit("\n");
936 }
937 else
938 {
939 Post_Process_Data(io);
940 }
941
942 if(io->n_otu < 3)
943 {
944 PhyML_Fprintf(stderr,"\n. PhyML needs at least three sequences to perform an analysis.");
945 assert(FALSE);
946 }
947
948 return io->data;
949 }
950
951 //////////////////////////////////////////////////////////////
952 //////////////////////////////////////////////////////////////
953
Post_Process_Data(option * io)954 void Post_Process_Data(option *io)
955 {
956 int i,j,swap;
957 align *data_buff;
958
959 for(i=0;i<io->data[0]->len;++i)
960 {
961 for(j=0;j<io->n_otu;++j)
962 {
963 if((io->data[j]->state[i] == '*') || (io->data[j]->state[i] == '?') || (io->data[j]->state[i] == '-')) io->data[j]->state[i] = 'X';
964 if((io->datatype == NT) && (io->data[j]->state[i] == 'N')) io->data[j]->state[i] = 'X';
965 if(io->data[j]->state[i] == 'U') io->data[j]->state[i] = 'T';
966 }
967 }
968
969 for(i=0;i<io->n_otu;++i) io->data[i]->len = io->data[0]->len;
970
971 /* Sequences are ordered alphabetically */
972 data_buff = NULL;
973 swap = TRUE;
974 /* swap = FALSE; */
975 while(swap == TRUE)
976 {
977 swap = FALSE;
978 for(i=0;i<io->n_otu-1;i++)
979 {
980 for(j=i+1;j<io->n_otu;j++)
981 {
982 if(strcmp(io->data[i]->name,io->data[j]->name) < 0)
983 {
984 swap = TRUE;
985 data_buff = io->data[i];
986 io->data[i] = io->data[j];
987 io->data[j] = data_buff;
988 }
989 }
990 }
991 }
992 }
993
994 //////////////////////////////////////////////////////////////
995 //////////////////////////////////////////////////////////////
996
997 /* align **Get_Seq_Nexus(option *io) */
998 /* { */
999 /* char *s,*ori_s; */
1000 /* char *token; */
1001 /* int in_comment; */
1002 /* nexcom *curr_com; */
1003 /* nexparm *curr_parm; */
1004 /* int nxt_token_t,cur_token_t; */
1005
1006 /* s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
1007 /* token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char)); */
1008
1009 /* ori_s = s; */
1010 /* in_comment = NO; */
1011 /* curr_com = NULL; */
1012 /* curr_parm = NULL; */
1013 /* nxt_token_t = NEXUS_COM; */
1014 /* cur_token_t = -1; */
1015
1016 /* while(fgets(s,T_MAX_LINE,io->fp_in_align)) */
1017 /* { */
1018 /* do */
1019 /* { */
1020 /* Get_Token(&s,token); */
1021
1022 /* /\* PhyML_Printf("\n. Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); *\/ */
1023
1024 /* if(token[0] == '\0') break; */
1025
1026 /* if(token[0] == ';') */
1027 /* { */
1028 /* curr_com = NULL; */
1029 /* curr_parm = NULL; */
1030 /* nxt_token_t = NEXUS_COM; */
1031 /* cur_token_t = -1; */
1032 /* break; /\* End of command *\/ */
1033 /* } */
1034
1035 /* if(nxt_token_t == NEXUS_EQUAL) */
1036 /* { */
1037 /* cur_token_t = NEXUS_VALUE; */
1038 /* nxt_token_t = NEXUS_PARM; */
1039 /* continue; */
1040 /* } */
1041
1042 /* if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE)) */
1043 /* { */
1044 /* Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list); */
1045 /* if(curr_com) */
1046 /* { */
1047 /* nxt_token_t = curr_com->nxt_token_t; */
1048 /* cur_token_t = curr_com->cur_token_t; */
1049 /* } */
1050 /* if(cur_token_t != NEXUS_VALUE) continue; */
1051 /* } */
1052
1053 /* if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE)) */
1054 /* { */
1055 /* Find_Nexus_Parm(token,&curr_parm,curr_com); */
1056 /* if(curr_parm) */
1057 /* { */
1058 /* nxt_token_t = curr_parm->nxt_token_t; */
1059 /* cur_token_t = curr_parm->cur_token_t; */
1060 /* } */
1061 /* if(cur_token_t != NEXUS_VALUE) continue; */
1062 /* } */
1063
1064 /* if(cur_token_t == NEXUS_VALUE) */
1065 /* { */
1066 /* if((curr_parm->fp)(token,curr_parm,io)) /\* Read in parameter value *\/ */
1067 /* { */
1068 /* nxt_token_t = NEXUS_PARM; */
1069 /* cur_token_t = -1; */
1070 /* } */
1071 /* } */
1072 /* } */
1073 /* while(strlen(token) > 0); */
1074 /* } */
1075
1076 /* Free(ori_s); */
1077 /* Free(token); */
1078
1079 /* return io->data; */
1080 /* } */
1081
1082 /* /\*********************************************************\/ */
1083
Get_Nexus_Data(FILE * fp,option * io)1084 void Get_Nexus_Data(FILE *fp, option *io)
1085 {
1086 char *token;
1087 nexcom *curr_com;
1088 nexparm *curr_parm;
1089 int nxt_token_t,cur_token_t;
1090
1091 token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char));
1092
1093 curr_com = NULL;
1094 curr_parm = NULL;
1095 nxt_token_t = NEXUS_COM;
1096 cur_token_t = -1;
1097
1098 do
1099 {
1100 if(!Get_Token(fp,token)) break;
1101
1102 /* PhyML_Printf("\n+ Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); */
1103
1104 if(token[0] == ';')
1105 {
1106 curr_com = NULL;
1107 curr_parm = NULL;
1108 nxt_token_t = NEXUS_COM;
1109 cur_token_t = -1;
1110 }
1111
1112 if(nxt_token_t == NEXUS_EQUAL)
1113 {
1114 cur_token_t = NEXUS_VALUE;
1115 nxt_token_t = NEXUS_PARM;
1116 continue;
1117 }
1118
1119 if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE))
1120 {
1121 Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list);
1122 if(curr_com)
1123 {
1124 nxt_token_t = curr_com->nxt_token_t;
1125 cur_token_t = curr_com->cur_token_t;
1126 }
1127 if(cur_token_t != NEXUS_VALUE) continue;
1128 }
1129
1130 if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE))
1131 {
1132 Find_Nexus_Parm(token,&curr_parm,curr_com);
1133 if(curr_parm)
1134 {
1135 nxt_token_t = curr_parm->nxt_token_t;
1136 cur_token_t = curr_parm->cur_token_t;
1137 }
1138 if(cur_token_t != NEXUS_VALUE) continue;
1139 }
1140
1141 if(cur_token_t == NEXUS_VALUE)
1142 {
1143 if((curr_parm->fp)(token,curr_parm,io)) /* Read in parameter value */
1144 {
1145 nxt_token_t = NEXUS_PARM;
1146 cur_token_t = -1;
1147 }
1148 }
1149 }
1150 while(strlen(token) > 0);
1151
1152 Free(token);
1153 }
1154
1155 //////////////////////////////////////////////////////////////
1156 //////////////////////////////////////////////////////////////
1157
Get_Token(FILE * fp,char * token)1158 int Get_Token(FILE *fp, char *token)
1159 {
1160 char c;
1161
1162 c = ' ';
1163
1164 while(c == ' ' || c == '\t' || c == '\n' || c == '\r')
1165 {
1166 c = fgetc(fp);
1167 if(c == EOF) return 0;
1168 }
1169
1170 if(c == '"')
1171 {
1172 do
1173 {
1174 *token = c;
1175 token++;
1176 c = fgetc(fp);
1177 if(c == EOF) return 0;
1178 }
1179 while(c != '"');
1180 *token = c;
1181 /* c = fgetc(fp); */
1182 if(c == EOF) return 0;
1183 *(token+1) = '\0';
1184 return 1;
1185 }
1186
1187 if(c == '[')
1188 {
1189 Skip_Comment(fp);
1190 c = fgetc(fp);
1191 *token = c;
1192 token++;
1193 if(c == EOF) return 0;
1194 return 1;
1195 }
1196
1197 if(c == '#') { *token = c; token++; }
1198 else if(c == ';') { *token = c; token++; }
1199 else if(c == ',') { *token = c; token++; }
1200 else if(c == '.') { *token = c; token++; }
1201 else if(c == '=') { *token = c; token++; }
1202 else if(c == '(') { *token = c; token++; }
1203 else if(c == ')') { *token = c; token++; }
1204 else if(c == '{') { *token = c; token++; }
1205 else if(c == '}') { *token = c; token++; }
1206 else if(c == '?') { *token = c; token++; }
1207 else if(c == '-') { *token = c; token++; }
1208 else
1209 {
1210 while(isgraph(c) && c != ';' && c != '-' && c != ',' && c != '=')
1211 {
1212 *(token++) = c;
1213 c = fgetc(fp);
1214 if(c == EOF) return 0;
1215 }
1216
1217 fseek(fp,-1*sizeof(char),SEEK_CUR);
1218
1219 }
1220 *token = '\0';
1221 return 1;
1222 }
1223
1224
1225 /* void Get_Token(char *line, char *token) */
1226 /* { */
1227 /* while(**line == ' ' || **line == '\t') (*line)++; */
1228
1229
1230 /* if(**line == '"') */
1231 /* { */
1232 /* do { *token = **line; (*line)++; token++; } while(**line != '"'); */
1233 /* *token = **line; */
1234 /* (*line)++; */
1235 /* *(token+1) = '\0'; */
1236 /* return; */
1237 /* } */
1238
1239 /* if(**line == '[') */
1240 /* { */
1241 /* int in_comment; */
1242
1243 /* in_comment = 1; */
1244 /* do */
1245 /* { */
1246 /* (*line)++; */
1247 /* if(**line == '[') */
1248 /* { */
1249 /* in_comment++; */
1250 /* } */
1251 /* else if(**line == ']') in_comment--; */
1252 /* } */
1253 /* while(in_comment); */
1254 /* (*line)++; */
1255 /* return; */
1256 /* } */
1257
1258
1259 /* if(**line == '#') {*token = **line; (*line)++; token++; } */
1260 /* else if(**line == ';') {*token = **line; (*line)++; token++; } */
1261 /* else if(**line == ',') {*token = **line; (*line)++; token++; } */
1262 /* else if(**line == '.') {*token = **line; (*line)++; token++; } */
1263 /* else if(**line == '=') {*token = **line; (*line)++; token++; } */
1264 /* else if(**line == '(') {*token = **line; (*line)++; token++; } */
1265 /* else if(**line == ')') {*token = **line; (*line)++; token++; } */
1266 /* else if(**line == '{') {*token = **line; (*line)++; token++; } */
1267 /* else if(**line == '}') {*token = **line; (*line)++; token++; } */
1268 /* else if(**line == '?') {*token = **line; (*line)++; token++; } */
1269 /* else if(**line == '-') {*token = **line; (*line)++; token++; } */
1270 /* else */
1271 /* { */
1272 /* while(isgraph(**line) && **line != ';' && **line != '=' && **line != ',') */
1273 /* { */
1274 /* *(token++) = **line; */
1275 /* (*line)++; */
1276 /* } */
1277 /* } */
1278 /* *token = '\0'; */
1279 /* } */
1280
1281 //////////////////////////////////////////////////////////////
1282 //////////////////////////////////////////////////////////////
1283
Get_Seq_Phylip(option * io)1284 align **Get_Seq_Phylip(option *io)
1285 {
1286 Read_Ntax_Len_Phylip(io->fp_in_align,&io->n_otu,&io->init_len);
1287
1288 if(io->n_otu > N_MAX_OTU)
1289 {
1290 PhyML_Fprintf(stderr,"\n. The number of taxa should not exceed %d",N_MAX_OTU);
1291 assert(FALSE);
1292 }
1293
1294 if(io->interleaved == YES) io->data = Read_Seq_Interleaved(io);
1295 else io->data = Read_Seq_Sequential(io);
1296
1297
1298 return io->data;
1299 }
1300
1301 //////////////////////////////////////////////////////////////
1302 //////////////////////////////////////////////////////////////
1303
Read_Ntax_Len_Phylip(FILE * fp,int * n_otu,int * n_tax)1304 void Read_Ntax_Len_Phylip(FILE *fp ,int *n_otu, int *n_tax)
1305 {
1306 char *line;
1307 int readok;
1308
1309 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1310
1311 readok = 0;
1312 do
1313 {
1314 if(fscanf(fp,"%s",line) == EOF)
1315 {
1316 Free(line);
1317 PhyML_Fprintf(stderr,"\n. PhyML can't read in this alignment.");
1318 PhyML_Fprintf(stderr,"\n. Could it be that sequence file is empty?");
1319 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
1320 Exit("\n");
1321 }
1322 else
1323 {
1324 if(strcmp(line,"\n") && strcmp(line,"\r") && strcmp(line,"\t"))
1325 {
1326 sscanf(line,"%d",n_otu);
1327 if(*n_otu <= 0) Warn_And_Exit("\n. The number of taxa cannot be negative.\n");
1328
1329 if(!fscanf(fp,"%s",line)) Exit("\n");
1330 sscanf(line,"%d",n_tax);
1331 if(*n_tax <= 0) Warn_And_Exit("\n. The sequence length cannot be negative.\n");
1332 else readok = 1;
1333 }
1334 }
1335 }while(!readok);
1336
1337 Free(line);
1338 }
1339
1340 //////////////////////////////////////////////////////////////
1341 //////////////////////////////////////////////////////////////
1342
Read_Seq_Sequential(option * io)1343 align **Read_Seq_Sequential(option *io)
1344 {
1345 int i;
1346 char *line;
1347 align **data;
1348 /* char c; */
1349 char *format;
1350
1351 format = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1352 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1353 data = (align **)mCalloc(io->n_otu,sizeof(align *));
1354
1355 /* while((c=fgetc(in))!='\n'); */
1356 /* while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1357
1358 sprintf(format, "%%%ds", T_MAX_NAME);
1359
1360
1361 for(i=0;i<io->n_otu;i++)
1362 {
1363 data[i] = (align *)mCalloc(1,sizeof(align));
1364 data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1365 data[i]->state = (char *)mCalloc(io->init_len*io->state_len+1,sizeof(char));
1366
1367 data[i]->is_ambigu = NULL;
1368 data[i]->len = 0;
1369
1370 if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1371
1372 Check_Sequence_Name(data[i]->name);
1373
1374 while(data[i]->len < io->init_len * io->state_len) assert(Read_One_Line_Seq(&data,i,io->fp_in_align));
1375
1376 if(data[i]->len != io->init_len * io->state_len)
1377 {
1378 PhyML_Fprintf(stderr,"\n. Err. Problem with species %s's sequence (check the format).\n",data[i]->name);
1379 PhyML_Fprintf(stderr,"\n. Observed sequence length: %d, expected length: %d\n",data[i]->len, io->init_len * io->state_len);
1380 Warn_And_Exit("");
1381 }
1382 }
1383
1384 for(i=0;i<io->n_otu;i++) data[i]->state[data[i]->len] = '\0';
1385
1386 Restrict_To_Coding_Position(data,io);
1387
1388 Free(format);
1389 Free(line);
1390
1391 return data;
1392 }
1393
1394 //////////////////////////////////////////////////////////////
1395 //////////////////////////////////////////////////////////////
1396
Read_Seq_Interleaved(option * io)1397 align **Read_Seq_Interleaved(option *io)
1398 {
1399 int i,end,num_block;
1400 char *line;
1401 align **data;
1402 /* char c; */
1403 char *format;
1404 fpos_t curr_pos;
1405
1406 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1407 format = (char *)mCalloc(T_MAX_NAME, sizeof(char));
1408 data = (align **)mCalloc(io->n_otu,sizeof(align *));
1409
1410 /* while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1411
1412 sprintf(format, "%%%ds", T_MAX_NAME);
1413
1414 end = 0;
1415 for(i=0;i<io->n_otu;i++)
1416 {
1417 data[i] = (align *)mCalloc(1,sizeof(align));
1418 data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1419 data[i]->state = (char *)mCalloc(io->init_len*io->state_len+1,sizeof(char));
1420
1421 data[i]->len = 0;
1422 data[i]->is_ambigu = NULL;
1423
1424 if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1425
1426 Check_Sequence_Name(data[i]->name);
1427
1428 if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1429 {
1430 end = 1;
1431 if((i != io->n_otu) && (i != io->n_otu-1))
1432 {
1433 PhyML_Fprintf(stderr,"\n. Err.: problem with species %s's sequence.\n",data[i]->name);
1434 PhyML_Fprintf(stderr,"\n. Observed sequence length: %d, expected length: %d\n",data[i]->len, io->init_len * io->state_len);
1435 Exit("");
1436 }
1437 break;
1438 }
1439 }
1440
1441 if(data[0]->len == io->init_len * io->state_len) end = 1;
1442
1443 /* if(end) printf("\n. finished yet '%c'\n",fgetc(io->fp_in_align)); */
1444 if(!end)
1445 {
1446 end = 0;
1447 num_block = 1;
1448 do
1449 {
1450 num_block++;
1451
1452 /* interblock */
1453 if(!fgets(line,T_MAX_LINE,io->fp_in_align)) break;
1454
1455 if(line[0] != 13 && line[0] != 10)
1456 {
1457 PhyML_Fprintf(stderr,"\n. Err.: one or more missing sequences in block %d.\n",num_block-1);
1458 Exit("");
1459 }
1460
1461 for(i=0;i<io->n_otu;i++) if(data[i]->len != io->init_len * io->state_len) break;
1462
1463 if(i == io->n_otu) break;
1464
1465 for(i=0;i<io->n_otu;i++)
1466 {
1467 /* Skip the taxon name, if any, in this interleaved block */
1468 fgetpos(io->fp_in_align,&curr_pos);
1469 if(!fscanf(io->fp_in_align,format,line))
1470 {
1471 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
1472 Warn_And_Exit("");
1473 }
1474 if(line && strcmp(line,data[i]->name)) fsetpos(io->fp_in_align,&curr_pos);
1475
1476 if(data[i]->len > io->init_len * io->state_len)
1477 {
1478 PhyML_Fprintf(stderr,"\n. Observed sequence length=%d expected length=%d.\n",data[i]->len,io->init_len * io->state_len);
1479 PhyML_Fprintf(stderr,"\n. Err.: Problem with species %s's sequence.\n",data[i]->name);
1480 Exit("");
1481 }
1482 else if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1483 {
1484 end = 1;
1485 if((i != io->n_otu) && (i != io->n_otu-1))
1486 {
1487 PhyML_Fprintf(stderr,"\n. Err.: Problem with species %s's sequence.\n",data[i]->name);
1488 PhyML_Fprintf(stderr,"\n. Observed sequence length: %d, expected length: %d.\n",data[i]->len, io->init_len * io->state_len);
1489 Exit("");
1490 }
1491 break;
1492 }
1493 }
1494 }while(!end);
1495 }
1496
1497 for(i=0;i<io->n_otu;i++) data[i]->state[data[i]->len] = '\0';
1498
1499 for(i=0;i<io->n_otu;i++)
1500 {
1501 if(data[i]->len != io->init_len * io->state_len)
1502 {
1503 PhyML_Fprintf(stderr,"\n. Check sequence '%s' length (expected length: %d, observed length: %d) [OTU %d].\n",data[i]->name,io->init_len,data[i]->len,i+1);
1504 Exit("");
1505 }
1506 }
1507
1508 Restrict_To_Coding_Position(data,io);
1509
1510 Free(format);
1511 Free(line);
1512
1513 return data;
1514 }
1515
1516 //////////////////////////////////////////////////////////////
1517 //////////////////////////////////////////////////////////////
1518
Read_One_Line_Seq(align *** data,int num_otu,FILE * in)1519 int Read_One_Line_Seq(align ***data, int num_otu, FILE *in)
1520 {
1521 char c = ' ';
1522 int nchar = 0;
1523
1524 while(1)
1525 {
1526 /* if((c == EOF) || (c == '\n') || (c == '\r')) break; */
1527
1528 if((c == 13) || (c == 10))
1529 {
1530 /* PhyML_Printf("[%d %d]\n",c,nchar); fflush(NULL); */
1531 if(!nchar)
1532 {
1533 c=(char)fgetc(in);
1534 continue;
1535 }
1536 else
1537 {
1538 /* PhyML_Printf("break\n"); */
1539 break;
1540 }
1541 }
1542 else if(c == EOF)
1543 {
1544 /* PhyML_Printf("EOL\n"); */
1545 break;
1546 }
1547 else if((c == ' ') || (c == '\t') || (c == 32))
1548 {
1549 /* PhyML_Printf("[%d]",c); */
1550 c=(char)fgetc(in);
1551 continue;
1552 }
1553
1554 nchar++;
1555 Uppercase(&c);
1556
1557 if(c == '.')
1558 {
1559 c = (*data)[0]->state[(*data)[num_otu]->len];
1560 if(!num_otu)
1561 Warn_And_Exit("\n. Err: Symbol \".\" should not appear in the first sequence\n");
1562 }
1563 (*data)[num_otu]->state[(*data)[num_otu]->len]=c;
1564 (*data)[num_otu]->len++;
1565 c = (char)fgetc(in);
1566 /* PhyML_Printf("[%c %d]",c,c); fflush(NULL); */
1567 if(c == ';') break;
1568 }
1569
1570 /* printf("\n. Exit nchar: %d [%d]\n",nchar,c==EOF); */
1571 if(c == EOF) return 0;
1572 else return 1;
1573 }
1574
1575 //////////////////////////////////////////////////////////////
1576 //////////////////////////////////////////////////////////////
1577
1578
Read_Io_Weights(option * io)1579 scalar_dbl *Read_Io_Weights(option *io)
1580 {
1581 scalar_dbl *w,*ori,*prev = NULL;
1582 double val;
1583
1584 assert(io->weight_file);
1585
1586 io->fp_weight_file = Openfile(io->weight_file,READ);
1587
1588 w = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1589 ori = w;
1590
1591 do
1592 {
1593 if(fscanf(io->fp_weight_file,"%lf,",&val) == EOF) break;
1594 w->v = (phydbl)val;
1595 w->next = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1596 prev = w;
1597 w = w->next;
1598 }
1599 while(1);
1600
1601 /* Remove the last allocated and empty element of the list */
1602 if(prev !=NULL && prev->next != NULL)
1603 {
1604 Free(prev->next);
1605 prev->next=NULL;
1606 }
1607
1608 fclose(io->fp_weight_file);
1609
1610 return(ori);
1611 }
1612
1613 //////////////////////////////////////////////////////////////
1614 //////////////////////////////////////////////////////////////
1615
Return_Tree_String_Phylip(FILE * fp_input_tree)1616 char *Return_Tree_String_Phylip(FILE *fp_input_tree)
1617 {
1618 char *line;
1619 int i;
1620 char c;
1621 int open,maxopen;
1622
1623 if(fp_input_tree == NULL)
1624 {
1625 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
1626 Warn_And_Exit("");
1627 }
1628
1629 do
1630 {
1631 c=fgetc(fp_input_tree);
1632 }
1633 while((c != '(') && (c != EOF));
1634
1635
1636 if(c==EOF) return NULL;
1637
1638 line = (char *)mCalloc(1,sizeof(char));
1639 open = 1;
1640 maxopen = open;
1641
1642 i=0;
1643 for(;;)
1644 {
1645 if((c == ' ') || (c == '\n'))
1646 {
1647 c=fgetc(fp_input_tree);
1648 if(c == EOF || c == ';') break;
1649 else continue;
1650 }
1651
1652 /* if(c == '[') */
1653 /* { */
1654 /* Skip_Comment(fp_input_tree); */
1655 /* c = fgetc(fp_input_tree); */
1656 /* if(c == EOF || c == ';') break; */
1657 /* } */
1658
1659 line = (char *)mRealloc(line,i+2,sizeof(char));
1660
1661 line[i]=c;
1662 i++;
1663 c=fgetc(fp_input_tree);
1664 if(c==EOF || c==';') break;
1665 if(c=='(') open++;
1666 if(c==')') open--;
1667 if(open>maxopen) maxopen = open;
1668 }
1669 line[i] = '\0';
1670
1671 /* if(maxopen == 1) return NULL; */
1672 return line;
1673 }
1674
1675 //////////////////////////////////////////////////////////////
1676 //////////////////////////////////////////////////////////////
1677
Read_Tree_File_Phylip(FILE * fp_input_tree)1678 t_tree *Read_Tree_File_Phylip(FILE *fp_input_tree)
1679 {
1680 char *line;
1681 t_tree *tree;
1682
1683 line = Return_Tree_String_Phylip(fp_input_tree);
1684 tree = Read_Tree(&line);
1685 Free(line);
1686
1687 return tree;
1688 }
1689
1690 //////////////////////////////////////////////////////////////
1691 //////////////////////////////////////////////////////////////
1692
Print_Site(calign * cdata,int num,int n_otu,char * sep,int stepsize,FILE * fp)1693 void Print_Site(calign *cdata, int num, int n_otu, char *sep, int stepsize, FILE *fp)
1694 {
1695 int i,j;
1696 PhyML_Fprintf(fp,"\n");
1697 for(i=0;i<n_otu;i++)
1698 {
1699 PhyML_Fprintf(fp,"%20s ",cdata->c_seq[i]->name);
1700 for(j=0;j<stepsize;j++)
1701 PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[num+j]);
1702 PhyML_Fprintf(fp,"%s",sep);
1703 }
1704 PhyML_Fprintf(fp,"%s",sep);
1705 }
1706
1707 //////////////////////////////////////////////////////////////
1708 //////////////////////////////////////////////////////////////
1709
Print_Site_Lk(t_tree * tree,FILE * fp)1710 void Print_Site_Lk(t_tree *tree, FILE *fp)
1711 {
1712 int site;
1713 int catg;
1714 char *s;
1715 phydbl postmean,sum;
1716
1717 assert(fp);
1718 rewind(fp);
1719
1720 if(tree->is_mixt_tree == YES)
1721 {
1722 MIXT_Print_Site_Lk(tree,fp);
1723 return;
1724 }
1725
1726 assert(tree->io->print_site_lnl == YES);
1727
1728 if(!tree->io->print_trace)
1729 {
1730 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1731
1732 PhyML_Fprintf(fp,"Note : P(D|M) is the probability of site D given the model M (i.e., the site likelihood)\n");
1733 if(tree->mod->ras->n_catg > 1 || tree->mod->ras->invar)
1734 {
1735 PhyML_Fprintf(fp,"P*(D|M,rr[x]) is the scaled probability of site D given the model M and the relative rate\n");
1736 PhyML_Fprintf(fp,"of evolution rr[x], where x is the class of rate to be considered.\n");
1737 PhyML_Fprintf(fp,"The actual conditional probability is given by P*(D|M,rr[x])/2^F, where\n");
1738 PhyML_Fprintf(fp,"F is the scaling factor (see column 'Scaler').\n");
1739 PhyML_Fprintf(fp,"For invariant sites, P(D|M,rr[0]=0) is the actual conditional probability\n");
1740 PhyML_Fprintf(fp,"(i.e., it is not scaled).\n");
1741 }
1742
1743 PhyML_Fprintf(fp,"\n\n");
1744
1745 sprintf(s,"Site");
1746 PhyML_Fprintf(fp, "%-12s",s);
1747
1748 sprintf(s,"P(D|M)");
1749 PhyML_Fprintf(fp,"%-15s",s);
1750
1751 sprintf(s,"Scaler");
1752 PhyML_Fprintf(fp,"%-7s",s);
1753
1754 sprintf(s,"Pattern");
1755 PhyML_Fprintf(fp, "%-9s",s);
1756
1757 if(tree->mod->ras->n_catg > 1)
1758 {
1759 for(catg=0;catg<tree->mod->ras->n_catg;catg++)
1760 {
1761 sprintf(s,"P*(D|M,rr[%d]=%5.4f)",catg+1,tree->mod->ras->gamma_rr->v[catg]);
1762 PhyML_Fprintf(fp,"%-23s",s);
1763 }
1764
1765 sprintf(s,"Posterior mean");
1766 PhyML_Fprintf(fp,"%-22s",s);
1767 }
1768
1769
1770 if(tree->mod->ras->invar)
1771 {
1772 sprintf(s,"P(D|M,rr[0]=0)");
1773 PhyML_Fprintf(fp,"%-16s",s);
1774 }
1775
1776 sprintf(s,"NDistinctStates");
1777 PhyML_Fprintf(fp,"%-16s",s);
1778
1779 PhyML_Fprintf(fp,"\n");
1780
1781 Init_Ui_Tips(tree);
1782
1783 for(site=0;site<tree->data->init_len;site++)
1784 {
1785 PhyML_Fprintf(fp,"%-12d",site+1);
1786
1787 PhyML_Fprintf(fp,"%-15g",tree->cur_site_lk[tree->data->sitepatt[site]]);
1788
1789 PhyML_Fprintf(fp,"%-7d",tree->fact_sum_scale[tree->data->sitepatt[site]]);
1790
1791 PhyML_Fprintf(fp,"%-9d",tree->data->sitepatt[site]);
1792
1793 if(tree->mod->ras->n_catg > 1)
1794 {
1795 for(catg=0;catg<tree->mod->ras->n_catg;catg++)
1796 {
1797 PhyML_Fprintf(fp,"%-23g",tree->unscaled_site_lk_cat[tree->data->sitepatt[site]*tree->mod->ras->n_catg + catg]);
1798 }
1799
1800 postmean = .0;
1801 for(catg=0;catg<tree->mod->ras->n_catg;catg++)
1802 postmean +=
1803 tree->mod->ras->gamma_rr->v[catg] *
1804 tree->unscaled_site_lk_cat[tree->data->sitepatt[site]*tree->mod->ras->n_catg + catg] *
1805 tree->mod->ras->gamma_r_proba->v[catg];
1806
1807 sum = .0;
1808 for(catg=0;catg<tree->mod->ras->n_catg;catg++)
1809 {
1810 sum +=
1811 tree->unscaled_site_lk_cat[tree->data->sitepatt[site]*tree->mod->ras->n_catg + catg] *
1812 tree->mod->ras->gamma_r_proba->v[catg];
1813 }
1814
1815 postmean /= sum;
1816
1817 PhyML_Fprintf(fp,"%-22g",postmean);
1818 }
1819
1820 if(tree->mod->ras->invar)
1821 {
1822 if((phydbl)tree->data->invar[tree->data->sitepatt[site]] > -0.5)
1823 PhyML_Fprintf(fp,"%-16g",tree->mod->e_frq->pi->v[tree->data->invar[tree->data->sitepatt[site]]]);
1824 else
1825 PhyML_Fprintf(fp,"%-16g",0.0);
1826 }
1827
1828
1829 PhyML_Fprintf(fp,"%-16d",Number_Of_Diff_States_One_Site(tree->data->sitepatt[site],tree));
1830
1831 PhyML_Fprintf(fp,"\n");
1832 }
1833 Free(s);
1834
1835 }
1836 else
1837 {
1838 for(site=0;site<tree->data->init_len;site++)
1839 PhyML_Fprintf(fp,"%.2f\t",log(tree->cur_site_lk[tree->data->sitepatt[site]]));
1840 PhyML_Fprintf(fp,"\n");
1841 }
1842 }
1843
1844 //////////////////////////////////////////////////////////////
1845 //////////////////////////////////////////////////////////////
1846
Print_Seq(FILE * fp,align ** data,int n_otu)1847 void Print_Seq(FILE *fp, align **data, int n_otu)
1848 {
1849 int i,j;
1850
1851 PhyML_Fprintf(fp,"%d\t%d\n",n_otu,data[0]->len);
1852 for(i=0;i<n_otu;i++)
1853 {
1854 PhyML_Fprintf(fp,"%s\t",data[i]->name);
1855 /* for(j=0;j<20;j++) */
1856 /* { */
1857 /* if(j<(int)strlen(data[i]->name)) */
1858 /* fputc(data[i]->name[j],fp); */
1859 /* else fputc(' ',fp); */
1860 /* } */
1861 /* PhyML_Printf("%10d ",i); */
1862 For(j,data[i]->len)
1863 {
1864 PhyML_Fprintf(fp,"%c",data[i]->state[j]);
1865 }
1866 PhyML_Fprintf(fp,"\n");
1867 }
1868 }
1869
1870 //////////////////////////////////////////////////////////////
1871 //////////////////////////////////////////////////////////////
1872
Print_CSeq(FILE * fp,int compressed,calign * cdata,t_tree * tree)1873 void Print_CSeq(FILE *fp, int compressed, calign *cdata, t_tree *tree)
1874 {
1875 int i,j;
1876 int n_otu;
1877
1878 n_otu = cdata->n_otu;
1879 if(cdata->format == PHYLIP)
1880 {
1881 PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
1882 }
1883 else if(cdata->format == NEXUS)
1884 {
1885 PhyML_Fprintf(fp,"#NEXUS\n");
1886 PhyML_Fprintf(fp,"begin data\n");
1887 PhyML_Fprintf(fp,"dimensions ntax=%d nchar=%d;\n",n_otu,cdata->init_len);
1888 PhyML_Fprintf(fp,"format sequential datatype=dna;\n");
1889 PhyML_Fprintf(fp,"matrix\n");
1890 }
1891 else
1892 {
1893 PhyML_Fprintf(fp,"This sample file is to be processed by IBDSim (http://www1.montpellier.inra.fr/CBGP/software/ibdsim/)");
1894 }
1895
1896 if(cdata->format == PHYLIP || cdata->format == NEXUS)
1897 {
1898 for(i=0;i<n_otu;i++)
1899 {
1900 for(j=0;j<50;j++)
1901 {
1902 if(j<(int)strlen(cdata->c_seq[i]->name))
1903 fputc(cdata->c_seq[i]->name[j],fp);
1904 else fputc(' ',fp);
1905 }
1906
1907 if(compressed == YES) /* Print out compressed sequences */
1908 PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
1909 else /* Print out uncompressed sequences */
1910 {
1911 for(j=0;j<cdata->init_len;j++)
1912 {
1913 PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
1914 }
1915 }
1916 PhyML_Fprintf(fp,"\n");
1917 }
1918 PhyML_Fprintf(fp,"\n");
1919
1920 if(cdata->format == NEXUS)
1921 {
1922 PhyML_Fprintf(fp,";\n");
1923 PhyML_Fprintf(fp,"END;\n");
1924 }
1925 }
1926 else if(cdata->format == IBDSIM)
1927 {
1928 for(i=0;i<cdata->init_len;i++) PhyML_Fprintf(fp,"\nlocus %6d",i);
1929 for(i=0;i<n_otu;i++)
1930 {
1931 PhyML_Fprintf(fp,"\npop");
1932 PhyML_Fprintf(fp,"%12f %12f , ",
1933 tree->a_nodes[i]->coord->lonlat[0],
1934 tree->a_nodes[i]->coord->lonlat[1]);
1935 for(j=0;j<cdata->init_len;j++)
1936 {
1937 switch(tree->a_nodes[i]->c_seq->state[j])
1938 {
1939 case 'A' :
1940 {
1941 PhyML_Fprintf(fp,"001 ");
1942 break;
1943 }
1944 case 'C' :
1945 {
1946 PhyML_Fprintf(fp,"002 ");
1947 break;
1948 }
1949 case 'G' :
1950 {
1951 PhyML_Fprintf(fp,"003 ");
1952 break;
1953 }
1954 case 'T' :
1955 {
1956 PhyML_Fprintf(fp,"004 ");
1957 break;
1958 }
1959 }
1960 }
1961 }
1962 }
1963 /* PhyML_Printf("\t"); */
1964 /* for(j=0;j<cdata->crunch_len;j++) */
1965 /* PhyML_Printf("%.0f ",cdata->wght[j]); */
1966 /* PhyML_Printf("\n"); */
1967 }
1968
1969 //////////////////////////////////////////////////////////////
1970 //////////////////////////////////////////////////////////////
1971
Print_CSeq_Select(FILE * fp,int compressed,calign * cdata,t_tree * tree)1972 void Print_CSeq_Select(FILE *fp, int compressed, calign *cdata, t_tree *tree)
1973 {
1974 int i,j;
1975 int n_otu;
1976 phydbl eps;
1977
1978 int slice = 14;
1979
1980 eps = 1.E-6;
1981 n_otu = 0;
1982 for(i=0;i<cdata->n_otu;i++)
1983 if(tree->times->nd_t[i] < tree->times->time_slice_lims[slice] + eps)
1984 n_otu++;
1985
1986 PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
1987
1988 n_otu = cdata->n_otu;
1989
1990 for(i=0;i<n_otu;i++)
1991 {
1992 if(tree->times->nd_t[i] < tree->times->time_slice_lims[slice] + eps)
1993 {
1994 for(j=0;j<50;j++)
1995 {
1996 if(j<(int)strlen(cdata->c_seq[i]->name))
1997 fputc(cdata->c_seq[i]->name[j],fp);
1998 else fputc(' ',fp);
1999 }
2000
2001 if(compressed == YES) /* Print out compressed sequences */
2002 PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
2003 else /* Print out uncompressed sequences */
2004 {
2005 for(j=0;j<cdata->init_len;j++)
2006 {
2007 PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
2008 }
2009 }
2010 PhyML_Fprintf(fp,"\n");
2011 }
2012 }
2013
2014 if(cdata->format == 1)
2015 {
2016 PhyML_Fprintf(fp,";\n");
2017 PhyML_Fprintf(fp,"END;\n");
2018 }
2019
2020
2021 /* PhyML_Printf("\t"); */
2022 /* for(j=0;j<cdata->crunch_len;j++) */
2023 /* PhyML_Printf("%.0f ",cdata->wght[j]); */
2024 /* PhyML_Printf("\n"); */
2025 }
2026
2027 //////////////////////////////////////////////////////////////
2028 //////////////////////////////////////////////////////////////
2029
Print_Dist(matrix * mat)2030 void Print_Dist(matrix *mat)
2031 {
2032 int i,j;
2033
2034 for(i=0;i<mat->n_otu;i++)
2035 {
2036 PhyML_Printf("%s ",mat->name[i]);
2037
2038 for(j=0;j<mat->n_otu;j++)
2039 PhyML_Printf("%9.6f ",mat->dist[i][j]);
2040 PhyML_Printf("\n");
2041 }
2042 }
2043
2044 //////////////////////////////////////////////////////////////
2045 //////////////////////////////////////////////////////////////
2046
Print_Node(t_node * a,t_node * d,t_tree * tree)2047 void Print_Node(t_node *a, t_node *d, t_tree *tree)
2048 {
2049 int i;
2050 int dir;
2051 dir = -1;
2052 for(i=0;i<3;i++) if(a->v[i] == d) {dir = i; break;}
2053 PhyML_Printf("Node nums: %3d %3d (dir:%3d) (a->anc:%3d) (d->anc:%3d) ta:%8.4f td:%8.4f t_min:%6.2f t_max:%6.2f",
2054 a->num,d->num,dir,a->anc?a->anc->num:(-1),
2055 d->anc?d->anc->num:(-1),
2056 tree->rates?tree->times->nd_t[a->num]:-1.,
2057 tree->rates?tree->times->nd_t[d->num]:-1.,
2058 tree->rates?tree->times->t_prior_min[a->num]:-1.,
2059 tree->rates?tree->times->t_prior_max[a->num]:-1.);
2060
2061 PhyML_Printf(" names = '%10s' '%10s' ; ",a->name,d->name);
2062 for(i=0;i<3;i++) if(a->v[i] == d)
2063 {
2064 if(a->b[i])
2065 {
2066 PhyML_Printf("Branch num = %3d%c (%3d %3d) %f",
2067 a->b[i]->num,a->b[i]==tree->e_root?'*':' ',a->b[i]->left->num,
2068 a->b[i]->rght->num,a->b[i]->l->v);
2069 if(a->b[i]->left->tax) PhyML_Printf(" WARNING LEFT->TAX!");
2070 break;
2071 }
2072 }
2073 PhyML_Printf("\n");
2074
2075 if(d->tax) return;
2076 else
2077 for(i=0;i<3;i++)
2078 if(d->v[i] != a && d->b[i] != tree->e_root) Print_Node(d,d->v[i],tree);
2079 }
2080
2081 //////////////////////////////////////////////////////////////
2082 //////////////////////////////////////////////////////////////
2083
Print_Node_Brief(t_node * a,t_node * d,t_tree * tree,FILE * fp)2084 void Print_Node_Brief(t_node *a, t_node *d, t_tree *tree, FILE *fp)
2085 {
2086 int i;
2087 int dir;
2088
2089 dir = -1;
2090 for(i=0;i<3;i++) if(a->v[i] == d) {dir = i; break;}
2091
2092 PhyML_Fprintf(fp,"\n");
2093 PhyML_Fprintf(fp,"Node nums: %3d %3d (dir:%3d)",
2094 a->num,d->num,dir);
2095
2096 PhyML_Fprintf(fp,"\tnames = '%10s' '%10s' ; ",a->name,d->name);
2097 for(i=0;i<3;i++) if(a->v[i] == d)
2098 {
2099 if(a->b[i])
2100 {
2101 PhyML_Fprintf(fp,"Branch num = %3d%c (%3d %3d) length:%10f",
2102 a->b[i]->num,a->b[i]==tree->e_root?'*':' ',a->b[i]->left->num,
2103 a->b[i]->rght->num,a->b[i]->l->v);
2104 if(a->b[i]->left->tax) PhyML_Printf(" WARNING LEFT->TAX!");
2105 break;
2106 }
2107 }
2108
2109 if(d->tax) return;
2110 else
2111 for(i=0;i<3;i++)
2112 if(d->v[i] != a && d->b[i] != tree->e_root) Print_Node_Brief(d,d->v[i],tree,fp);
2113 }
2114
2115 //////////////////////////////////////////////////////////////
2116 //////////////////////////////////////////////////////////////
2117
Print_Model(t_mod * mod)2118 void Print_Model(t_mod *mod)
2119 {
2120 int i,j,k;
2121
2122 PhyML_Printf("\n. name=%s",mod->modelname->s);
2123 PhyML_Printf("\n. string=%s",mod->custom_mod_string);
2124 PhyML_Printf("\n. mod_num=%d",mod->mod_num);
2125 PhyML_Printf("\n. ns=%d",mod->ns);
2126 PhyML_Printf("\n. n_catg=%d",mod->ras->n_catg);
2127 PhyML_Printf("\n. kappa=%f",mod->kappa->v);
2128 PhyML_Printf("\n. alpha=%f",mod->ras->alpha->v);
2129 PhyML_Printf("\n. lambda=%f",mod->lambda->v);
2130 PhyML_Printf("\n. pinvar=%f",mod->ras->pinvar->v);
2131 PhyML_Printf("\n. br_len_mult=%f",mod->br_len_mult->v);
2132 PhyML_Printf("\n. whichmodel=%d",mod->whichmodel);
2133 PhyML_Printf("\n. update_eigen=%d",mod->update_eigen);
2134 PhyML_Printf("\n. bootstrap=%d",mod->io->n_boot_replicates);
2135 PhyML_Printf("\n. n_diff_rr=%d",mod->r_mat->n_diff_rr);
2136 PhyML_Printf("\n. invar=%d",mod->ras->invar);
2137 PhyML_Printf("\n. use_m4mod=%d",mod->use_m4mod);
2138 PhyML_Printf("\n. gamma_median=%d",mod->ras->gamma_median);
2139 PhyML_Printf("\n. state_len=%d",mod->io->state_len);
2140 PhyML_Printf("\n. log_l=%d",mod->log_l);
2141 PhyML_Printf("\n. l_min=%f",mod->l_min);
2142 PhyML_Printf("\n. l_max=%f",mod->l_max);
2143 PhyML_Printf("\n. free_mixt_rates=%d",mod->ras->free_mixt_rates);
2144 PhyML_Printf("\n. gamma_mgf_bl=%d",mod->gamma_mgf_bl);
2145
2146 PhyML_Printf("\n. Pi\n");
2147 for(i=0;i<mod->ns;i++) PhyML_Printf(" %f ",mod->e_frq->pi->v[i]);
2148 PhyML_Printf("\n");
2149 for(i=0;i<mod->ns;i++) PhyML_Printf(" %f ",mod->e_frq->pi_unscaled->v[i]);
2150
2151 PhyML_Printf("\n. Rates\n");
2152 for(i=0;i<mod->ras->n_catg;i++) PhyML_Printf(" %f ",mod->ras->gamma_r_proba->v[i]);
2153 PhyML_Printf("\n");
2154 for(i=0;i<mod->ras->n_catg;i++) PhyML_Printf(" %f ",mod->ras->gamma_r_proba_unscaled->v[i]);
2155 PhyML_Printf("\n");
2156 for(i=0;i<mod->ras->n_catg;i++) PhyML_Printf(" %f ",mod->ras->gamma_rr->v[i]);
2157 PhyML_Printf("\n");
2158 for(i=0;i<mod->ras->n_catg;i++) PhyML_Printf(" %f ",mod->ras->gamma_rr_unscaled->v[i]);
2159
2160
2161
2162 PhyML_Printf("\n. Qmat \n");
2163 if(mod->whichmodel == CUSTOM)
2164 {
2165 fflush(NULL);
2166 for(i=0;i<6;i++) {PhyML_Printf(" %12f ",mod->r_mat->rr->v[i]); fflush(NULL);}
2167 for(i=0;i<6;i++) {PhyML_Printf(" %12f ",mod->r_mat->rr_val->v[i]); fflush(NULL);}
2168 for(i=0;i<6;i++) {PhyML_Printf(" %12d ",mod->r_mat->rr_num->v[i]); fflush(NULL);}
2169 for(i=0;i<6;i++) {PhyML_Printf(" %12d ",mod->r_mat->n_rr_per_cat->v[i]); fflush(NULL);}
2170 }
2171 for(i=0;i<mod->ns;i++)
2172 {
2173 PhyML_Printf(" ");
2174 for(j=0;j<4;j++)
2175 PhyML_Printf("%8.5f ",mod->r_mat->qmat->v[i*4+j]);
2176 PhyML_Printf("\n");
2177 }
2178
2179 PhyML_Printf("\n. Freqs");
2180 PhyML_Printf("\n");
2181 for(i=0;i<mod->ns;i++) PhyML_Printf(" %12f ",mod->e_frq->user_b_freq->v[i]);
2182 PhyML_Printf("\n");
2183 for(i=0;i<mod->ns;i++) PhyML_Printf(" %12f ",mod->e_frq->pi->v[i]);
2184 PhyML_Printf("\n");
2185 for(i=0;i<mod->ns;i++) PhyML_Printf(" %12f ",mod->e_frq->pi_unscaled->v[i]);
2186
2187 PhyML_Printf("\n. Eigen\n");
2188 For(i,2*mod->ns) PhyML_Printf(" %f ",mod->eigen->space[i]);
2189 /* PhyML_Printf("\n"); */
2190 /* For(i,2*mod->ns) PhyML_Printf(" %f ",mod->eigen->space_int[i]); */
2191 PhyML_Printf("\n");
2192 for(i=0;i<mod->ns;i++) PhyML_Printf(" %f ",mod->eigen->e_val[i]);
2193 PhyML_Printf("\n");
2194 for(i=0;i<mod->ns;i++) PhyML_Printf(" %f ",mod->eigen->e_val_im[i]);
2195 PhyML_Printf("\n");
2196 For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect[i]);
2197 PhyML_Printf("\n");
2198 For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect_im[i]);
2199 PhyML_Printf("\n");
2200 For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->l_e_vect[i]);
2201 PhyML_Printf("\n");
2202 For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->q[i]);
2203 PhyML_Printf("\n");
2204
2205 PhyML_Printf("\n. Pij");
2206 for(k=0;k<mod->ras->n_catg;k++)
2207 {
2208 PMat(0.01*mod->ras->gamma_rr->v[k],mod,mod->ns*mod->ns*k,mod->Pij_rr->v,NULL);
2209 PhyML_Printf("\n. l=%f\n",0.01*mod->ras->gamma_rr->v[k]);
2210 for(i=0;i<mod->ns;i++)
2211 {
2212 PhyML_Printf(" ");
2213 for(j=0;j<mod->ns;j++)
2214 PhyML_Printf("%8.5f ",mod->Pij_rr->v[k*mod->ns*mod->ns+i*mod->ns+j]);
2215 PhyML_Printf("\n");
2216 }
2217 }
2218
2219 PhyML_Printf("\n");
2220
2221 fflush(NULL);
2222 }
2223
2224 //////////////////////////////////////////////////////////////
2225 //////////////////////////////////////////////////////////////
2226
Print_Mat(matrix * mat)2227 void Print_Mat(matrix *mat)
2228 {
2229 int i,j;
2230
2231 PhyML_Printf("%d",mat->n_otu);
2232 PhyML_Printf("\n");
2233
2234 for(i=0;i<mat->n_otu;i++)
2235 {
2236 for(j=0;j<13;j++)
2237 {
2238 if(j>=(int)strlen(mat->name[i])) putchar(' ');
2239 else putchar(mat->name[i][j]);
2240 }
2241
2242 for(j=0;j<mat->n_otu;j++)
2243 {
2244 char s[2]="-";
2245 if(mat->dist[i][j] < .0)
2246 PhyML_Printf("%12s",s);
2247 else
2248 PhyML_Printf("%12f",mat->dist[i][j]);
2249 }
2250 PhyML_Printf("\n");
2251 }
2252 }
2253
2254 //////////////////////////////////////////////////////////////
2255 //////////////////////////////////////////////////////////////
2256
Openfile(char * filename,int mode)2257 FILE *Openfile(char *filename, int mode)
2258 {
2259 FILE *fp;
2260 char *s;
2261 int open_test=0;
2262
2263 s = filename;
2264
2265 fp = NULL;
2266
2267 switch(mode)
2268 {
2269 case READ :
2270 {
2271 while(!(fp = (FILE *)fopen(s,"r")) && ++open_test<10)
2272 {
2273 PhyML_Printf("\n. Can't open file '%s', enter a new name : ",s);
2274 Getstring_Stdin(s);
2275 }
2276 break;
2277 }
2278 case WRITE :
2279 {
2280 fp = (FILE *)fopen(s,"w");
2281 break;
2282 }
2283 case APPEND :
2284 {
2285 fp = (FILE *)fopen(s,"a");
2286 break;
2287 }
2288 case READWRITE :
2289 {
2290 fp = (FILE *)fopen(s,"w+");
2291 break;
2292 }
2293
2294 default : break;
2295
2296 }
2297
2298 /* Free(s); */
2299
2300 return fp;
2301 }
2302
2303 //////////////////////////////////////////////////////////////
2304 //////////////////////////////////////////////////////////////
2305
Print_Fp_Out(FILE * fp_out,time_t t_beg,time_t t_end,t_tree * tree,option * io,int n_data_set,int num_tree,int add_citation,int precision)2306 void Print_Fp_Out(FILE *fp_out, time_t t_beg, time_t t_end, t_tree *tree, option *io, int n_data_set, int num_tree, int add_citation, int precision)
2307 {
2308 char *s;
2309 char format[8];
2310 div_t hour,min;
2311 int i, j;
2312
2313 if (precision > 0) snprintf(format,8,"%%.%huf",(unsigned short)precision);
2314
2315 if(n_data_set == 1)
2316 {
2317 rewind(fp_out);
2318 Print_Banner_Small(fp_out);
2319 }
2320
2321 PhyML_Fprintf(fp_out,"\n. Sequence filename: \t\t\t%s", Basename(io->in_align_file));
2322 PhyML_Fprintf(fp_out,"\n. Data set: \t\t\t\t#%d",n_data_set);
2323
2324 if(io->mod->s_opt->random_input_tree) PhyML_Fprintf(fp_out,"\n. Random init tree: \t\t\t#%d",num_tree+1);
2325 else if(io->n_trees > 1) PhyML_Fprintf(fp_out,"\n. Starting tree number: \t\t#%d",num_tree+1);
2326
2327 /* if(io->mod->s_opt->opt_topo) */
2328 /* PhyML_Fprintf(fp_out,"\n. Tree topology search: \t\tSPRs"); */
2329 /* else */
2330 /* PhyML_Fprintf(fp_out,"\n. Tree topology: \t\t\tfixed"); */
2331
2332 /* was after Sequence file ; moved here FLT */
2333 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2334 if(io->in_tree == 2)
2335 {
2336 strcat(strcat(strcat(s,"user tree ("),io->in_tree_file),")");
2337 }
2338 else
2339 {
2340 if(!io->mod->s_opt->random_input_tree)
2341 {
2342 if(io->in_tree == 0) strcat(s,"BioNJ");
2343 if(io->in_tree == 1) strcat(s,"parsimony");
2344 }
2345 else
2346 strcat(s,"random tree");
2347 }
2348
2349 PhyML_Fprintf(fp_out,"\n. Initial tree: \t\t\t%s",s);
2350 Free(s);
2351
2352 if(tree->io->datatype == NT)
2353 {
2354 PhyML_Fprintf(fp_out,"\n. Model of nucleotides substitution: \t%s",tree->mod->modelname->s);
2355 if(io->mod->whichmodel == CUSTOM)
2356 PhyML_Fprintf(fp_out," (%s)",io->mod->custom_mod_string);
2357 }
2358 else if(tree->io->datatype == AA)
2359 {
2360 PhyML_Fprintf(fp_out,"\n. Model of amino acids substitution: \t%s",tree->mod->modelname->s);
2361 if(io->mod->whichmodel == CUSTOMAA) PhyML_Fprintf(fp_out," (%s)",tree->mod->aa_rate_mat_file->s);
2362 }
2363 else
2364 {
2365 fprintf(fp_out,"\n. Substitution model: \t\t\t%s",tree->mod->modelname->s);
2366 }
2367
2368 PhyML_Fprintf(fp_out,"\n. Number of taxa: \t\t\t%d",tree->n_otu);/*added FLT*/
2369
2370 PhyML_Fprintf(fp_out,"\n. Log-likelihood: \t\t\t%.5f",tree->c_lnL);/*was last ; moved here FLT*/
2371
2372 Unconstraint_Lk(tree);
2373 PhyML_Fprintf(fp_out,"\n. Unconstrained log-likelihood: \t%.5f",tree->unconstraint_lk);
2374
2375 Composite_Lk(tree);
2376 PhyML_Fprintf(fp_out,"\n. Composite log-likelihood: \t\t%.5f",tree->composite_lk);
2377
2378 PhyML_Fprintf(fp_out,"\n. Parsimony: \t\t\t\t%d",tree->c_pars);
2379
2380 PhyML_Fprintf(fp_out,"\n. Tree size: \t\t\t\t%.5f",Get_Tree_Size(tree));
2381
2382 /* if(tree->mod->ras->n_catg > 1 && tree->mod->ras->free_mixt_rates == NO) */
2383 if(tree->mod->ras->free_mixt_rates == NO)
2384 {
2385 PhyML_Fprintf(fp_out,"\n. Discrete gamma model: \t\t%s","Yes");
2386 PhyML_Fprintf(fp_out,"\n - Number of classes: \t\t\t%d",tree->mod->ras->n_catg);
2387 PhyML_Fprintf(fp_out,"\n - Gamma shape parameter: \t\t%.3f",tree->mod->ras->alpha->v);
2388 for(i=0;i<tree->mod->ras->n_catg;i++)
2389 {
2390 PhyML_Fprintf(fp_out,"\n - Relative rate in class %d: \t\t%.5f [freq=%4f] \t\t",i+1,tree->mod->ras->gamma_rr->v[i],tree->mod->ras->gamma_r_proba->v[i]);
2391 }
2392 }
2393 else if(tree->mod->ras->free_mixt_rates == YES)
2394 {
2395 int *rk;
2396 rk = Ranks(tree->mod->ras->gamma_rr->v,tree->mod->ras->n_catg);
2397 PhyML_Fprintf(fp_out,"\n. FreeRate model: \t\t\t%s","Yes");
2398 PhyML_Fprintf(fp_out,"\n - Number of classes: \t\t\t%d",tree->mod->ras->n_catg);
2399 for(i=0;i<tree->mod->ras->n_catg;i++)
2400 {
2401 PhyML_Fprintf(fp_out,"\n - Relative rate in class %d: \t\t%.5f [freq=%4f] \t\t",i+1,tree->mod->ras->gamma_rr->v[rk[i]],tree->mod->ras->gamma_r_proba->v[rk[i]]);
2402 }
2403 Free(rk);
2404 }
2405
2406 if(tree->mod->ras->invar) PhyML_Fprintf(fp_out,"\n. Proportion of invariant: \t\t%.3f",tree->mod->ras->pinvar->v);
2407
2408 if(tree->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp_out,"\n. Variance of branch lengths: \t\t%f",tree->mod->l_var_sigma);
2409
2410
2411 /*was before Discrete gamma model ; moved here FLT*/
2412 if((tree->mod->whichmodel == K80) ||
2413 (tree->mod->whichmodel == HKY85) ||
2414 (tree->mod->whichmodel == F84))
2415 {
2416 PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio: \t");
2417 if (precision > 0)
2418 PhyML_Fprintf(fp_out,format,tree->mod->kappa->v);
2419 else
2420 PhyML_Fprintf(fp_out,"%f",tree->mod->kappa->v);
2421 }
2422 else if(tree->mod->whichmodel == TN93)
2423 {
2424 PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio for purines: \t\t");
2425 if (precision > 0)
2426 PhyML_Fprintf(fp_out,format,tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
2427 else
2428 PhyML_Fprintf(fp_out,"%f",tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
2429
2430 PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio for pyrimidines: \t");
2431 if (precision > 0)
2432 PhyML_Fprintf(fp_out,format,tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
2433 else
2434 PhyML_Fprintf(fp_out,"%f",tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
2435 }
2436
2437 if(tree->io->datatype == NT)
2438 {
2439 PhyML_Fprintf(fp_out,"\n. Nucleotides frequencies:");
2440 if (precision > 0)
2441 {
2442 PhyML_Fprintf(fp_out,"\n - f(A)= ");
2443 PhyML_Fprintf(fp_out,format,tree->mod->e_frq->pi->v[0]);
2444 PhyML_Fprintf(fp_out,"\n - f(C)= ");
2445 PhyML_Fprintf(fp_out,format,tree->mod->e_frq->pi->v[1]);
2446 PhyML_Fprintf(fp_out,"\n - f(G)= ");
2447 PhyML_Fprintf(fp_out,format,tree->mod->e_frq->pi->v[2]);
2448 PhyML_Fprintf(fp_out,"\n - f(T)= ");
2449 PhyML_Fprintf(fp_out,format,tree->mod->e_frq->pi->v[3]);
2450 }
2451 else
2452 {
2453 PhyML_Fprintf(fp_out,"\n - f(A)= %8.5f",tree->mod->e_frq->pi->v[0]);
2454 PhyML_Fprintf(fp_out,"\n - f(C)= %8.5f",tree->mod->e_frq->pi->v[1]);
2455 PhyML_Fprintf(fp_out,"\n - f(G)= %8.5f",tree->mod->e_frq->pi->v[2]);
2456 PhyML_Fprintf(fp_out,"\n - f(T)= %8.5f",tree->mod->e_frq->pi->v[3]);
2457 }
2458 }
2459
2460 /*****************************************/
2461 if((tree->mod->whichmodel == GTR) ||
2462 (tree->mod->whichmodel == CUSTOM))
2463 {
2464 Update_Qmat_GTR(tree->mod->r_mat->rr->v,
2465 tree->mod->r_mat->rr_val->v,
2466 tree->mod->r_mat->rr_num->v,
2467 tree->mod->e_frq->pi->v,
2468 tree->mod->r_mat->qmat->v,
2469 tree->mod->s_opt->opt_rr);
2470
2471 PhyML_Fprintf(fp_out,"\n");
2472 PhyML_Fprintf(fp_out,". GTR relative rate parameters :");
2473 if (precision > 0)
2474 {
2475 PhyML_Fprintf(fp_out,"\n A <-> C ");
2476 if (tree->mod->r_mat->rr->v[0] < 10)
2477 PhyML_Fprintf(fp_out," ");
2478 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[0]);
2479 PhyML_Fprintf(fp_out,"\n A <-> G ");
2480 if (tree->mod->r_mat->rr->v[1] < 10)
2481 PhyML_Fprintf(fp_out," ");
2482 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[1]);
2483 PhyML_Fprintf(fp_out,"\n A <-> T ");
2484 if (tree->mod->r_mat->rr->v[2] < 10)
2485 PhyML_Fprintf(fp_out," ");
2486 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[2]);
2487 PhyML_Fprintf(fp_out,"\n C <-> G ");
2488 if (tree->mod->r_mat->rr->v[3] < 10)
2489 PhyML_Fprintf(fp_out," ");
2490 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[3]);
2491 PhyML_Fprintf(fp_out,"\n C <-> T ");
2492 if (tree->mod->r_mat->rr->v[4] < 10)
2493 PhyML_Fprintf(fp_out," ");
2494 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[4]);
2495 PhyML_Fprintf(fp_out,"\n G <-> T ");
2496 if (tree->mod->r_mat->rr->v[5] < 10)
2497 PhyML_Fprintf(fp_out," ");
2498 PhyML_Fprintf(fp_out,format, tree->mod->r_mat->rr->v[5]);
2499 }
2500 else
2501 {
2502 PhyML_Fprintf(fp_out,"\n A <-> C %8.5f", tree->mod->r_mat->rr->v[0]);
2503 PhyML_Fprintf(fp_out,"\n A <-> G %8.5f", tree->mod->r_mat->rr->v[1]);
2504 PhyML_Fprintf(fp_out,"\n A <-> T %8.5f", tree->mod->r_mat->rr->v[2]);
2505 PhyML_Fprintf(fp_out,"\n C <-> G %8.5f", tree->mod->r_mat->rr->v[3]);
2506 PhyML_Fprintf(fp_out,"\n C <-> T %8.5f", tree->mod->r_mat->rr->v[4]);
2507 PhyML_Fprintf(fp_out,"\n G <-> T %8.5f", tree->mod->r_mat->rr->v[5]);
2508 }
2509
2510 PhyML_Fprintf(fp_out,"\n. Instantaneous rate matrix : ");
2511 if (precision > 0)
2512 {
2513 PhyML_Fprintf(fp_out,"\n [A");
2514 for(i=0;i<precision+4;i++)
2515 PhyML_Fprintf(fp_out,"-");
2516 PhyML_Fprintf(fp_out,"C");
2517 for(i=0;i<precision+4;i++)
2518 PhyML_Fprintf(fp_out,"-");
2519 PhyML_Fprintf(fp_out,"G");
2520 for(i=0;i<precision+4;i++)
2521 PhyML_Fprintf(fp_out,"-");
2522 PhyML_Fprintf(fp_out,"T");
2523 for(i=0;i<precision+1;i++)
2524 PhyML_Fprintf(fp_out,"-");
2525 PhyML_Fprintf(fp_out,"]\n");
2526 for(i=0;i<4;i++)
2527 {
2528 for(j=0;j<4;j++)
2529 {
2530 if (i == j)
2531 PhyML_Fprintf(fp_out," ");
2532 else
2533 PhyML_Fprintf(fp_out," ");
2534 PhyML_Fprintf(fp_out,format,tree->mod->r_mat->qmat->v[i*4+j]);
2535 }
2536 PhyML_Fprintf(fp_out,"\n");
2537 }
2538 }
2539 else
2540 {
2541 PhyML_Fprintf(fp_out,"\n [A---------C---------G---------T------]\n");
2542 for(i=0;i<4;i++)
2543 {
2544 PhyML_Fprintf(fp_out," ");
2545 for(j=0;j<4;j++)
2546 PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->r_mat->qmat->v[i*4+j]);
2547 PhyML_Fprintf(fp_out,"\n");
2548 }
2549 }
2550 //PhyML_Fprintf(fp_out,"\n");
2551 }
2552
2553 /*****************************************/
2554 if(io->ratio_test == 1)
2555 {
2556 PhyML_Fprintf(fp_out,". aLRT statistics to test branches");
2557 }
2558 else if(io->ratio_test == 2)
2559 {
2560 PhyML_Fprintf(fp_out,". aLRT branch supports (cubic approximation, mixture of Chi2s distribution)");
2561 }
2562
2563 PhyML_Fprintf(fp_out,"\n");
2564 PhyML_Fprintf(fp_out,"\n. Run ID:\t\t\t\t%s", (io->append_run_ID) ? (io->run_id_string): ("none"));
2565 PhyML_Fprintf(fp_out,"\n. Random seed:\t\t\t\t%d", io->r_seed);
2566 PhyML_Fprintf(fp_out,"\n. Subtree patterns aliasing:\t\t%s",io->do_alias_subpatt?"yes":"no");
2567 PhyML_Fprintf(fp_out,"\n. Version:\t\t\t\t%s", VERSION);
2568
2569 hour = div(t_end-t_beg,3600);
2570 min = div(t_end-t_beg,60 );
2571 min.quot -= hour.quot*60;
2572
2573 PhyML_Fprintf(fp_out,"\n. Time used:\t\t\t\t%dh%dm%ds (%d seconds)", hour.quot,min.quot,(int)(t_end-t_beg)%60,(int)(t_end-t_beg));
2574
2575 if(add_citation == YES)
2576 {
2577 PhyML_Fprintf(fp_out,"\n\n");
2578 PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2579 PhyML_Fprintf(fp_out," Suggested citations:\n");
2580 PhyML_Fprintf(fp_out," S. Guindon, JF. Dufayard, V. Lefort, M. Anisimova, W. Hordijk, O. Gascuel\n");
2581 PhyML_Fprintf(fp_out," \"New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0.\"\n");
2582 PhyML_Fprintf(fp_out," Systematic Biology. 2010. 59(3):307-321.\n");
2583 PhyML_Fprintf(fp_out,"\n");
2584 PhyML_Fprintf(fp_out," S. Guindon & O. Gascuel\n");
2585 PhyML_Fprintf(fp_out," \"A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood\"\n");
2586 PhyML_Fprintf(fp_out," Systematic Biology. 2003. 52(5):696-704.\n");
2587 PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2588 }
2589 else
2590 {
2591 PhyML_Fprintf(fp_out,"\n\n");
2592 PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2593 PhyML_Fprintf(fp_out,"\n");
2594 }
2595 }
2596
2597 //////////////////////////////////////////////////////////////
2598 //////////////////////////////////////////////////////////////
2599
2600 /*FLT wrote this function*/
Print_Fp_Out_Lines(FILE * fp_out,time_t t_beg,time_t t_end,t_tree * tree,option * io,int n_data_set)2601 void Print_Fp_Out_Lines(FILE *fp_out, time_t t_beg, time_t t_end, t_tree *tree, option *io, int n_data_set)
2602 {
2603 char *s;
2604 /*div_t hour,min;*/
2605
2606 if (n_data_set==1)
2607 {
2608
2609 PhyML_Fprintf(fp_out,". Sequence file : [%s]\n\n", Basename(io->in_align_file));
2610
2611 if((tree->io->datatype == NT) || (tree->io->datatype == AA))
2612 {
2613 (tree->io->datatype == NT)?
2614 (PhyML_Fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",io->mod->modelname->s)):
2615 (PhyML_Fprintf(fp_out,". Model of amino acids substitution : %s\n\n",io->mod->modelname->s));
2616 }
2617
2618 s = (char *)mCalloc(100,sizeof(char));
2619
2620 switch(io->in_tree)
2621 {
2622 case 0: { strcpy(s,"BioNJ"); break; }
2623 case 1: { strcpy(s,"parsimony"); break; }
2624 case 2: { strcpy(s,"user tree (");
2625 strcat(s,io->in_tree_file);
2626 strcat(s,")"); break; }
2627 }
2628
2629 PhyML_Fprintf(fp_out,". Initial tree : [%s]\n\n",s);
2630
2631 Free(s);
2632
2633 PhyML_Fprintf(fp_out,"\n");
2634
2635 /*headline 1*/
2636 PhyML_Fprintf(fp_out, ". Data\t");
2637
2638 PhyML_Fprintf(fp_out,"Nb of \t");
2639
2640 PhyML_Fprintf(fp_out,"Likelihood\t");
2641
2642 PhyML_Fprintf(fp_out, "Discrete \t");
2643
2644 if(tree->mod->ras->n_catg > 1)
2645 PhyML_Fprintf(fp_out, "Number of \tGamma shape\t");
2646
2647 PhyML_Fprintf(fp_out,"Proportion of\t");
2648
2649 if(tree->mod->whichmodel <= 6)
2650 PhyML_Fprintf(fp_out,"Transition/ \t");
2651
2652 PhyML_Fprintf(fp_out,"Nucleotides frequencies \t");
2653
2654 if((tree->mod->whichmodel == GTR) ||
2655 (tree->mod->whichmodel == CUSTOM))
2656 PhyML_Fprintf(fp_out,"Instantaneous rate matrix \t");
2657
2658 /* PhyML_Fprintf(fp_out,"Time\t");*/
2659
2660 PhyML_Fprintf(fp_out, "\n");
2661
2662
2663 /*headline 2*/
2664 PhyML_Fprintf(fp_out, " set\t");
2665
2666 PhyML_Fprintf(fp_out,"taxa\t");
2667
2668 PhyML_Fprintf(fp_out,"loglk \t");
2669
2670 PhyML_Fprintf(fp_out, "gamma model\t");
2671
2672 if(tree->mod->ras->n_catg > 1)
2673 PhyML_Fprintf(fp_out, "categories\tparameter \t");
2674
2675 PhyML_Fprintf(fp_out,"invariant \t");
2676
2677 if(tree->mod->whichmodel <= 6)
2678 PhyML_Fprintf(fp_out,"transversion\t");
2679
2680 PhyML_Fprintf(fp_out,"f(A) f(C) f(G) f(T) \t");
2681
2682 if((tree->mod->whichmodel == GTR) ||
2683 (tree->mod->whichmodel == CUSTOM))
2684 PhyML_Fprintf(fp_out,"[A---------C---------G---------T------]\t");
2685
2686 /* PhyML_PhyML_Fprintf(fp_out,"used\t");*/
2687
2688 PhyML_Fprintf(fp_out, "\n");
2689
2690
2691 /*headline 3*/
2692 if(tree->mod->whichmodel == TN93)
2693 {
2694 PhyML_Fprintf(fp_out," \t \t \t \t");
2695 if(tree->mod->ras->n_catg > 1) PhyML_Fprintf(fp_out," \t \t");
2696 PhyML_Fprintf(fp_out," \t");
2697 PhyML_Fprintf(fp_out,"purines pyrimid.\t");
2698 PhyML_Fprintf(fp_out, "\n");
2699 }
2700
2701 PhyML_Fprintf(fp_out, "\n");
2702 }
2703
2704
2705 /*line items*/
2706
2707 PhyML_Fprintf(fp_out," #%d\t",n_data_set);
2708
2709 PhyML_Fprintf(fp_out,"%d \t",tree->n_otu);
2710
2711 PhyML_Fprintf(fp_out,"%.5f\t",tree->c_lnL);
2712
2713 PhyML_Fprintf(fp_out,"%s \t",
2714 (tree->mod->ras->n_catg>1)?("Yes"):("No "));
2715 if(tree->mod->ras->n_catg > 1)
2716 {
2717 PhyML_Fprintf(fp_out,"%d \t",tree->mod->ras->n_catg);
2718 PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->ras->alpha->v);
2719 }
2720
2721 /*if(tree->mod->ras->invar)*/
2722 PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->ras->pinvar->v);
2723
2724 if(tree->mod->whichmodel <= 5)
2725 {
2726 PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->kappa->v);
2727 }
2728 else if(tree->mod->whichmodel == TN93)
2729 {
2730 PhyML_Fprintf(fp_out,"%.3f ",
2731 tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
2732 PhyML_Fprintf(fp_out,"%.3f\t",
2733 tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
2734 }
2735
2736
2737 if(tree->io->datatype == NT)
2738 {
2739 PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[0]);
2740 PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[1]);
2741 PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[2]);
2742 PhyML_Fprintf(fp_out,"%8.5f\t",tree->mod->e_frq->pi->v[3]);
2743 }
2744 /*
2745 hour = div(t_end-t_beg,3600);
2746 min = div(t_end-t_beg,60 );
2747
2748 min.quot -= hour.quot*60;
2749
2750 PhyML_Fprintf(fp_out,"%dh%dm%ds\t", hour.quot,min.quot,(int)(t_end-t_beg)%60);
2751 if(t_end-t_beg > 60)
2752 PhyML_Fprintf(fp_out,". -> %d seconds\t",(int)(t_end-t_beg));
2753 */
2754
2755 /*****************************************/
2756 if((tree->mod->whichmodel == GTR) || (tree->mod->whichmodel == CUSTOM))
2757 {
2758 int i,j;
2759
2760 for(i=0;i<4;i++)
2761 {
2762 if (i!=0) {
2763 /*format*/
2764 PhyML_Fprintf(fp_out," \t \t \t \t");
2765 if(tree->mod->ras->n_catg > 1) PhyML_Fprintf(fp_out," \t \t");
2766 PhyML_Fprintf(fp_out," \t \t");
2767 }
2768 for(j=0;j<4;j++)
2769 PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->r_mat->qmat->v[i*4+j]);
2770 if (i<3) PhyML_Fprintf(fp_out,"\n");
2771 }
2772 }
2773 /*****************************************/
2774
2775 PhyML_Fprintf(fp_out, "\n\n");
2776 }
2777
2778 //////////////////////////////////////////////////////////////
2779 //////////////////////////////////////////////////////////////
2780
2781
Print_Freq(t_tree * tree)2782 void Print_Freq(t_tree *tree)
2783 {
2784
2785 switch(tree->io->datatype)
2786 {
2787 case NT:
2788 {
2789 PhyML_Printf("A : %f\n",tree->mod->e_frq->pi->v[0]);
2790 PhyML_Printf("C : %f\n",tree->mod->e_frq->pi->v[1]);
2791 PhyML_Printf("G : %f\n",tree->mod->e_frq->pi->v[2]);
2792 PhyML_Printf("T : %f\n",tree->mod->e_frq->pi->v[3]);
2793 break;
2794 }
2795 case AA:
2796 {
2797 PhyML_Printf("A : %f\n",tree->mod->e_frq->pi->v[0]);
2798 PhyML_Printf("R : %f\n",tree->mod->e_frq->pi->v[1]);
2799 PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[2]);
2800 PhyML_Printf("D : %f\n",tree->mod->e_frq->pi->v[3]);
2801 PhyML_Printf("C : %f\n",tree->mod->e_frq->pi->v[4]);
2802 PhyML_Printf("Q : %f\n",tree->mod->e_frq->pi->v[5]);
2803 PhyML_Printf("E : %f\n",tree->mod->e_frq->pi->v[6]);
2804 PhyML_Printf("G : %f\n",tree->mod->e_frq->pi->v[7]);
2805 PhyML_Printf("H : %f\n",tree->mod->e_frq->pi->v[8]);
2806 PhyML_Printf("I : %f\n",tree->mod->e_frq->pi->v[9]);
2807 PhyML_Printf("L : %f\n",tree->mod->e_frq->pi->v[10]);
2808 PhyML_Printf("K : %f\n",tree->mod->e_frq->pi->v[11]);
2809 PhyML_Printf("M : %f\n",tree->mod->e_frq->pi->v[12]);
2810 PhyML_Printf("F : %f\n",tree->mod->e_frq->pi->v[13]);
2811 PhyML_Printf("P : %f\n",tree->mod->e_frq->pi->v[14]);
2812 PhyML_Printf("S : %f\n",tree->mod->e_frq->pi->v[15]);
2813 PhyML_Printf("T : %f\n",tree->mod->e_frq->pi->v[16]);
2814 PhyML_Printf("W : %f\n",tree->mod->e_frq->pi->v[17]);
2815 PhyML_Printf("Y : %f\n",tree->mod->e_frq->pi->v[18]);
2816 PhyML_Printf("V : %f\n",tree->mod->e_frq->pi->v[19]);
2817 PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[20]);
2818 break;
2819 }
2820 default : {break;}
2821 }
2822 }
2823
2824 //////////////////////////////////////////////////////////////
2825 //////////////////////////////////////////////////////////////
2826
2827
Print_Settings(option * io)2828 void Print_Settings(option *io)
2829 {
2830 int answer;
2831 char *s;
2832
2833 s = (char *)mCalloc(100,sizeof(char));
2834
2835 PhyML_Printf("\n\n\n");
2836 PhyML_Printf("\n\n");
2837
2838
2839 PhyML_Printf("\n ////////////////////////////////////.\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\");
2840 PhyML_Printf("\n \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\.//////////////////////////////////////////\n");
2841
2842 PhyML_Printf("\n . Sequence filename:\t\t\t\t %s", Basename(io->in_align_file));
2843
2844 if(io->datatype == NT) strcpy(s,"dna");
2845 else if(io->datatype == AA) strcpy(s,"aa");
2846 else strcpy(s,"generic");
2847
2848 PhyML_Printf("\n . Data type:\t\t\t\t\t %s",s);
2849 PhyML_Printf("\n . Alphabet size:\t\t\t\t %d",io->mod->ns);
2850
2851 PhyML_Printf("\n . Sequence format:\t\t\t\t %s", io->interleaved ? "interleaved": "sequential");
2852 PhyML_Printf("\n . Number of data sets:\t\t\t\t %d", io->n_data_sets);
2853
2854 PhyML_Printf("\n . Nb of bootstrapped data sets:\t\t\t %d", io->n_boot_replicates);
2855
2856 if (io->n_boot_replicates > 0)
2857 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t no");
2858 else
2859 {
2860 if(io->ratio_test == 1)
2861 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (aLRT statistics)");
2862 else if(io->ratio_test == 2)
2863 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (Chi2-based parametric branch supports)");
2864 else if(io->ratio_test == 3)
2865 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (Minimum of SH-like and Chi2-based branch supports)");
2866 else if(io->ratio_test == 4)
2867 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (SH-like branch supports)");
2868 else if(io->ratio_test == 5)
2869 PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (aBayes branch supports)");
2870 }
2871
2872 PhyML_Printf("\n . Model name:\t\t\t\t\t %s", io->mod->modelname->s);
2873
2874 if(io->datatype == AA && io->mod->whichmodel == CUSTOMAA) PhyML_Printf(" (%s)",io->mod->aa_rate_mat_file->s);
2875
2876 if (io->datatype == NT)
2877 {
2878 if ((io->mod->whichmodel == K80) ||
2879 (io->mod->whichmodel == HKY85)||
2880 (io->mod->whichmodel == F84) ||
2881 (io->mod->whichmodel == TN93))
2882 {
2883 if(io->mod->s_opt && io->mod->s_opt->opt_kappa)
2884 PhyML_Printf("\n . Ts/tv ratio:\t\t\t\t\t estimated");
2885 else
2886 PhyML_Printf("\n . Ts/tv ratio:\t\t\t\t\t %f", io->mod->kappa->v);
2887 }
2888 }
2889
2890 if(io->mod->s_opt && io->mod->s_opt->opt_pinvar)
2891 PhyML_Printf("\n . Proportion of invariable sites:\t\t estimated");
2892 else
2893 PhyML_Printf("\n . Proportion of invariable sites:\t\t %f", io->mod->ras->pinvar->v);
2894
2895
2896 if(io->mod->ras->free_mixt_rates == NO)
2897 PhyML_Printf("\n . RAS model:\t\t\t\t\t discrete Gamma");
2898 else
2899 PhyML_Printf("\n . RAS model:\t\t\t\t\t FreeRate");
2900
2901 PhyML_Printf("\n . Number of subst. rate catgs:\t\t\t %d", io->mod->ras->n_catg);
2902
2903 if(io->mod->ras->n_catg > 1)
2904 {
2905
2906 if(io->mod->ras->free_mixt_rates == NO)
2907 {
2908 if(io->mod->s_opt && io->mod->s_opt->opt_alpha)
2909 PhyML_Printf("\n . Gamma distribution parameter:\t\t\t estimated");
2910 else
2911 PhyML_Printf("\n . Gamma distribution parameter:\t\t\t %f", io->mod->ras->alpha->v);
2912 PhyML_Printf("\n . 'Middle' of each rate class:\t\t\t %s",(io->mod->ras->gamma_median)?("median"):("mean"));
2913 }
2914 }
2915
2916
2917 if(io->datatype == AA)
2918 PhyML_Printf("\n . Amino acid equilibrium frequencies:\t\t %s", (io->mod->s_opt->opt_state_freq) ? ("empirical"):("model"));
2919 else if(io->datatype == NT)
2920 {
2921 if((io->mod->whichmodel != JC69) &&
2922 (io->mod->whichmodel != K80) &&
2923 (io->mod->whichmodel != F81))
2924 {
2925 if(io->mod->s_opt && !io->mod->e_frq->user_state_freq)
2926 {
2927 PhyML_Printf("\n . Nucleotide equilibrium frequencies:\t\t %s", (io->mod->s_opt->opt_state_freq) ? ("ML"):("empirical"));
2928 }
2929 else
2930 {
2931 PhyML_Printf("\n . Nucleotide equilibrium frequencies:\t\t %s","user-defined");
2932 }
2933 }
2934 }
2935
2936 PhyML_Printf("\n . Optimise tree topology:\t\t\t %s", (io->mod->s_opt && io->mod->s_opt->opt_topo) ? "yes": "no");
2937
2938 switch(io->in_tree)
2939 {
2940 case 0: { strcpy(s,"BioNJ"); break; }
2941 case 1: { strcpy(s,"parsimony"); break; }
2942 case 2: { strcpy(s,"user tree (");
2943 strcat(s,Basename(io->in_tree_file));
2944 strcat(s,")"); break; }
2945 }
2946
2947 if(io->mod->s_opt && io->mod->s_opt->opt_topo)
2948 {
2949 /* if(io->mod->s_opt->topo_search == NNI_MOVE) PhyML_Printf("\n . Tree topology search:\t\t\t\t NNIs"); */
2950 /* else if(io->mod->s_opt->topo_search == SPR_MOVE) PhyML_Printf("\n . Tree topology search:\t\t\t\t SPRs"); */
2951 /* else if(io->mod->s_opt->topo_search == BEST_OF_NNI_AND_SPR) PhyML_Printf("\n . Tree topology search:\t\t\t\t Best of NNIs and SPRs"); */
2952
2953 PhyML_Printf("\n . Starting tree:\t\t\t\t %s",s);
2954
2955 PhyML_Printf("\n . Add random input tree:\t\t\t %s", (io->mod->s_opt->random_input_tree) ? "yes": "no");
2956 if(io->mod->s_opt->random_input_tree)
2957 PhyML_Printf("\n . Number of random starting trees:\t\t %d", io->mod->s_opt->n_rand_starts);
2958 }
2959 else
2960 if(io->mod->s_opt && !io->mod->s_opt->random_input_tree)
2961 PhyML_Printf("\n . Evaluated tree:\t\t\t\t \"%s\"",s);
2962
2963 PhyML_Printf("\n . Optimise branch lengths:\t\t\t %s", (io->mod->s_opt && io->mod->s_opt->opt_bl) ? "yes": "no");
2964
2965 PhyML_Printf("\n . Minimum length of an edge:\t\t\t %g",io->mod->l_min);
2966
2967 answer = 0;
2968 if(io->mod->s_opt &&
2969 (io->mod->s_opt->opt_alpha ||
2970 io->mod->s_opt->opt_kappa ||
2971 io->mod->s_opt->opt_lambda ||
2972 io->mod->s_opt->opt_pinvar ||
2973 io->mod->s_opt->opt_rr)) answer = 1;
2974
2975 PhyML_Printf("\n . Optimise substitution model parameters:\t %s", (answer) ? "yes": "no");
2976
2977 PhyML_Printf("\n . Run ID:\t\t\t\t\t %s", (io->append_run_ID) ? (io->run_id_string): ("none"));
2978 PhyML_Printf("\n . Random seed:\t\t\t\t\t %d", io->r_seed);
2979 PhyML_Printf("\n . Subtree patterns aliasing:\t\t\t %s",io->do_alias_subpatt?"yes":"no");
2980 PhyML_Printf("\n . Version:\t\t\t\t\t %s", VERSION);
2981 PhyML_Printf("\n . Byte alignment:\t\t\t\t %d",BYTE_ALIGN);
2982 PhyML_Printf("\n . AVX enabled:\t\t\t\t\t %s",
2983 #if defined(__AVX__)
2984 "yes"
2985 #else
2986 "no"
2987 #endif
2988 );
2989 PhyML_Printf("\n . SSE enabled:\t\t\t\t\t %s",
2990 #if defined(__SSE3__)
2991 "yes"
2992 #else
2993 "no"
2994 #endif
2995 );
2996
2997
2998 /* PhyML_Printf("\n\n \u205C \u205C \u205C \u205C \u205C \u205C \u205C \u205C \u205C \u205C \u205C \u205C \n"); */
2999
3000 PhyML_Printf("\n");
3001 PhyML_Printf("\n ////////////////////////////////////.\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\");
3002 PhyML_Printf("\n \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\.//////////////////////////////////////////\n");
3003
3004
3005 PhyML_Printf("\n\n");
3006 fflush(NULL);
3007
3008 Free(s);
3009 }
3010
Print_Banner(FILE * fp)3011 void Print_Banner(FILE *fp)
3012 {
3013 PhyML_Fprintf(fp,"\n");
3014 PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3015 PhyML_Fprintf(fp," \n");
3016 PhyML_Fprintf(fp," --- PhyML %s --- \n",VERSION);
3017 PhyML_Fprintf(fp," \n");
3018 PhyML_Fprintf(fp," A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood \n");
3019 PhyML_Fprintf(fp," Stephane Guindon & Olivier Gascuel \n");
3020 PhyML_Fprintf(fp," \n");
3021 PhyML_Fprintf(fp," http://www.atgc-montpellier.fr/phyml \n");
3022 PhyML_Fprintf(fp," \n");
3023 PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier \n");
3024 PhyML_Fprintf(fp," \n");
3025 PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3026 }
3027
3028 //////////////////////////////////////////////////////////////
3029 //////////////////////////////////////////////////////////////
3030
3031
Print_Banner_Small(FILE * fp)3032 void Print_Banner_Small(FILE *fp)
3033 {
3034 PhyML_Fprintf(fp,"\n");
3035 PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3036 PhyML_Fprintf(fp," --- PhyML %s --- \n",VERSION);
3037 PhyML_Fprintf(fp," http://www.atgc-montpellier.fr/phyml \n");
3038 PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier \n");
3039 PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3040 }
3041
3042 //////////////////////////////////////////////////////////////
3043 //////////////////////////////////////////////////////////////
3044
3045
3046
Print_Data_Set_Number(option * io,FILE * fp)3047 void Print_Data_Set_Number(option *io, FILE *fp)
3048 {
3049 PhyML_Fprintf(fp,"\n");
3050 PhyML_Fprintf(fp," \n");
3051 PhyML_Fprintf(fp," [ Data set number %3d ] \n",io->curr_gt+1);
3052 PhyML_Fprintf(fp," \n");
3053 }
3054 //////////////////////////////////////////////////////////////
3055 //////////////////////////////////////////////////////////////
3056
Print_Lk(t_tree * tree,char * string)3057 void Print_Lk(t_tree *tree, char *string)
3058 {
3059 t_tree *loc_tree;
3060
3061 loc_tree = tree;
3062 /*! Rewind back to the first mixt_tree */
3063 while(loc_tree->prev) loc_tree = loc_tree->prev;
3064
3065 time(&(loc_tree->t_current));
3066 PhyML_Printf("\n. (%5d sec) [%15.4f] %s",
3067 (int)(loc_tree->t_current-loc_tree->t_beg),Get_Lk(tree),
3068 string);
3069 #ifndef QUIET
3070 fflush(NULL);
3071 #endif
3072 }
3073
3074 //////////////////////////////////////////////////////////////
3075 //////////////////////////////////////////////////////////////
3076
Print_List(t_ll * list)3077 void Print_List(t_ll *list)
3078 {
3079 t_ll *ll;
3080 ll = list->head;
3081 do
3082 {
3083 t_node *n;
3084 n = (t_node *)ll->v;
3085 PhyML_Printf("\n. list elem: %p next: %p prev: %p [%d] %p %p",(void *)ll,(void *)ll->next,(void *)ll->prev,n->num,(void *)ll->head,(void *)ll->tail);
3086 ll = ll->next;
3087 }
3088 while(ll != NULL);
3089 }
3090
3091 //////////////////////////////////////////////////////////////
3092 //////////////////////////////////////////////////////////////
3093
Print_Pars(t_tree * tree)3094 void Print_Pars(t_tree *tree)
3095 {
3096 time(&(tree->t_current));
3097 PhyML_Printf("\n. (%5d sec) [%5d]",(int)(tree->t_current-tree->t_beg),tree->c_pars);
3098 #ifndef QUIET
3099 fflush(NULL);
3100 #endif
3101 }
3102
3103 //////////////////////////////////////////////////////////////
3104 //////////////////////////////////////////////////////////////
3105
Print_Lk_And_Pars(t_tree * tree)3106 void Print_Lk_And_Pars(t_tree *tree)
3107 {
3108 time(&(tree->t_current));
3109
3110 PhyML_Printf("\n. (%5d sec) [%15.4f] [%5d]",
3111 (int)(tree->t_current-tree->t_beg),
3112 tree->c_lnL,tree->c_pars);
3113
3114 #ifndef QUIET
3115 fflush(NULL);
3116 #endif
3117 }
3118
3119 //////////////////////////////////////////////////////////////
3120 //////////////////////////////////////////////////////////////
3121
Read_Qmat(phydbl * daa,phydbl * pi,FILE * fp)3122 void Read_Qmat(phydbl *daa, phydbl *pi, FILE *fp)
3123 {
3124 int i,j;
3125 phydbl sum;
3126 double val;
3127
3128 assert(fp);
3129 rewind(fp);
3130
3131 for(i=1;i<20;i++)
3132 {
3133 for(j=0;j<19;j++)
3134 {
3135 /* if(!fscanf(fp,"%lf",&(daa[i*20+j]))) Exit("\n"); */
3136 if(!fscanf(fp,"%lf",&val))
3137 {
3138 PhyML_Fprintf(stderr,"\n. Rate matrix file does not appear to have a proper format. Please refer to the documentation.");
3139 Exit("\n");
3140 }
3141 daa[i*20+j] = (phydbl)val;
3142 daa[j*20+i] = daa[i*20+j];
3143 if(j == i-1) break;
3144 }
3145 }
3146
3147
3148 for(i=0;i<20;i++)
3149 {
3150 if(!fscanf(fp,"%lf",&val)) Exit("\n");
3151 pi[i] = (phydbl)val;
3152 }
3153 sum = .0;
3154 for(i=0;i<20;i++) sum += pi[i];
3155 if(FABS(sum - 1.) > 1.E-06)
3156 {
3157 PhyML_Printf("\n. Sum of amino-acid frequencies: %f",sum);
3158 PhyML_Printf("\n. Scaling amino-acid frequencies...\n");
3159 for(i=0;i<20;i++) pi[i] /= sum;
3160 }
3161 }
3162
3163 //////////////////////////////////////////////////////////////
3164 //////////////////////////////////////////////////////////////
3165
3166
Print_Qmat_AA(phydbl * daa,phydbl * pi)3167 void Print_Qmat_AA(phydbl *daa, phydbl *pi)
3168 {
3169 int i,j,cpt;
3170
3171 cpt = 0;
3172 for(i=0;i<20;i++)
3173 {
3174 for(j=0;j<i;j++)
3175 {
3176 PhyML_Printf("daa[%2d*20+%2d] = %10f; ",i,j,daa[i*20+j]);
3177 cpt++;
3178 if(!(cpt%4)) PhyML_Printf("\n");
3179 }
3180 }
3181
3182 PhyML_Printf("\n\n");
3183 PhyML_Printf("for (i=0; i<naa; i++) for (j=0; j<i; j++) daa[j*naa+i] = daa[i*naa+j];\n\n");
3184 for(i=0;i<20;i++) PhyML_Printf("pi[%d] = %f; ",i,pi[i]);
3185 PhyML_Printf("\n");
3186 PhyML_Printf("Ala\tArg\tAsn\tAsp\tCys\tGln\tGlu\tGly\tHis\tIle\tLeu\tLys\tMet\tPhe\tPro\tSer\tThr\tTrp\tTyr\tVal\n");
3187 }
3188
3189
3190 //////////////////////////////////////////////////////////////
3191 //////////////////////////////////////////////////////////////
3192
Print_Square_Matrix_Generic(int n,phydbl * mat)3193 void Print_Square_Matrix_Generic(int n, phydbl *mat)
3194 {
3195 int i,j;
3196
3197 PhyML_Printf("\n");
3198 for(i=0;i<n;i++)
3199 {
3200 PhyML_Printf("[%3d]",i);
3201 for(j=0;j<n;j++)
3202 {
3203 PhyML_Printf("%12.5f ",mat[i*n+j]);
3204 }
3205 PhyML_Printf("\n");
3206 }
3207 PhyML_Printf("\n");
3208 }
3209
3210 //////////////////////////////////////////////////////////////
3211 //////////////////////////////////////////////////////////////
3212
3213
Print_Diversity(FILE * fp,t_tree * tree)3214 void Print_Diversity(FILE *fp, t_tree *tree)
3215 {
3216
3217 Print_Diversity_Pre(tree->a_nodes[0],
3218 tree->a_nodes[0]->v[0],
3219 tree->a_nodes[0]->b[0],
3220 fp,
3221 tree);
3222
3223 /* mean_div_left = .0; */
3224 /* for(k=0;k<ns;k++) */
3225 /* { */
3226 /* mean_div_left += (k+1) * tree->a_edges[j]->div_post_pred_left[k]; */
3227 /* } */
3228 /* mean_div_rght = .0; */
3229 /* for(k=0;k<ns;k++) mean_div_rght += (k+1) * tree->a_edges[j]->div_post_pred_rght[k]; */
3230
3231 /* mean_div_left /= (phydbl)tree->data->init_len; */
3232 /* mean_div_rght /= (phydbl)tree->data->init_len; */
3233
3234 /* PhyML_Fprintf(fp,"%4d 0 %f\n",j,mean_div_left); */
3235 /* PhyML_Fprintf(fp,"%4d 1 %f\n",j,mean_div_rght); */
3236
3237
3238 /* mean_div_left = .0; */
3239 /* for(k=0;k<ns;k++) mean_div_left += tree->a_edges[j]->div_post_pred_left[k]; */
3240
3241 /* mean_div_rght = .0; */
3242 /* for(k=0;k<ns;k++) */
3243 /* { */
3244 /* mean_div_rght += tree->a_edges[j]->div_post_pred_rght[k]; */
3245 /* } */
3246
3247 /* if((mean_div_left != tree->data->init_len) || (mean_div_rght != tree->data->init_len)) */
3248 /* { */
3249 /* PhyML_Printf("\n. mean_div_left = %f mean_div_rght = %f init_len = %d", */
3250 /* mean_div_left,mean_div_rght,tree->data->init_len); */
3251 /* PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
3252 /* Warn_And_Exit(""); */
3253 /* } */
3254 }
3255
3256 //////////////////////////////////////////////////////////////
3257 //////////////////////////////////////////////////////////////
3258
Print_Diversity_Pre(t_node * a,t_node * d,t_edge * b,FILE * fp,t_tree * tree)3259 void Print_Diversity_Pre(t_node *a, t_node *d, t_edge *b, FILE *fp, t_tree *tree)
3260 {
3261 int k,ns;
3262
3263 ns = -1;
3264
3265 if(d->tax) return;
3266 else
3267 {
3268
3269 if(tree->io->datatype == NT) ns = 4;
3270 else if(tree->io->datatype == AA) ns = 20;
3271
3272 if(d == b->left) for(k=0;k<ns;k++) PhyML_Fprintf(fp,"%4d 0 %2d %4d\n",b->num,k,b->div_post_pred_left[k]);
3273 else for(k=0;k<ns;k++) PhyML_Fprintf(fp,"%4d 1 %2d %4d\n",b->num,k,b->div_post_pred_rght[k]);
3274
3275 for(k=0;k<3;k++) if(d->v[k] != a) Print_Diversity_Pre(d,d->v[k],d->b[k],fp,tree);
3276 }
3277
3278 }
3279
3280 //////////////////////////////////////////////////////////////
3281 //////////////////////////////////////////////////////////////
3282
Print_Tip_Partials(t_tree * tree,t_node * d)3283 void Print_Tip_Partials(t_tree* tree, t_node* d)
3284 {
3285 if(!d->tax)
3286 {
3287 fprintf(stdout,"Node %d is not a Taxa\n",d->num);fflush(stdout);
3288 return;
3289 }
3290
3291 assert(d->b[0]->rght == d);
3292 assert(d->b[0]->rght->tax);
3293 fprintf(stdout,"Taxa/Node %d\n",d->num);
3294 int site, i;
3295 for(site=0;site<tree->n_pattern;++site)
3296 {
3297 fprintf(stdout,"[%d: ",site);
3298 for(i=0;i<tree->mod->ns;++i)
3299 {
3300 int tip_partial = d->b[0]->p_lk_tip_r[site*tree->mod->ns + i];
3301 fprintf(stdout,"%d",tip_partial);fflush(stdout);
3302 }
3303 fprintf(stdout,"] ");fflush(stdout);
3304 }
3305 fprintf(stdout,"\n");fflush(stdout);
3306 }
3307
3308 //////////////////////////////////////////////////////////////
3309 //////////////////////////////////////////////////////////////
3310
Print_Edge_PMats(t_tree * tree,t_edge * b)3311 void Print_Edge_PMats(t_tree* tree, t_edge* b)
3312 {
3313 phydbl *Pij;
3314 #ifdef BEAGLE
3315 Pij = (phydbl*)malloc(tree->mod->ns * tree->mod->ns * tree->mod->ras->n_catg * sizeof(phydbl)); if (NULL==Pij) Warn_And_Exit(__PRETTY_FUNCTION__);
3316 int ret = beagleGetTransitionMatrix(tree->b_inst, b->Pij_rr_idx, Pij);
3317 if(ret<0){
3318 fprintf(stderr, "beagleGetTransitionMatrix() on instance %i failed:%i\n\n",tree->b_inst,ret);
3319 Free(Pij);
3320 Exit("");
3321 }
3322 #else
3323 Pij = b->Pij_rr;
3324 #endif
3325 fprintf(stdout,"\nflattened P-Matrices (for each rate category) state*state*num_rates[%d*%d*%d] for branch num:%i\n",tree->mod->ns,tree->mod->ns,tree->mod->ras->n_catg, b->num);
3326 int i;
3327 for(i=0;i<tree->mod->ns * tree->mod->ns * tree->mod->ras->n_catg;++i){
3328 fprintf(stdout,"%f,",Pij[i]);
3329 fflush(stdout);
3330 }
3331 fprintf(stdout,"\n");
3332 fflush(stdout);
3333 #ifdef BEAGLE
3334 Free(Pij);
3335 #endif
3336 }
3337
3338 //////////////////////////////////////////////////////////////
3339 //////////////////////////////////////////////////////////////
3340
Print_All_Edge_PMats(t_tree * tree)3341 void Print_All_Edge_PMats(t_tree* tree)
3342 {
3343 int i;
3344 for(i=0;i < 2*tree->n_otu-3;++i) Print_Edge_PMats(tree, tree->a_edges[i]);
3345 }
3346
3347 //////////////////////////////////////////////////////////////
3348 //////////////////////////////////////////////////////////////
3349
Print_Edge_Likelihoods(t_tree * tree,t_edge * b,bool scientific)3350 void Print_Edge_Likelihoods(t_tree* tree, t_edge* b, bool scientific/*Print in scientific notation?*/)
3351 {
3352 int catg, site, j;
3353 char* fmt = scientific ? "[%d,%d,%d]%e ":"[%d,%d,%d]%f "; //rate category, site, state, likelihood
3354
3355 phydbl* lk_left = b->p_lk_left;
3356 phydbl* lk_right = b->p_lk_rght;
3357
3358 fprintf(stdout,"\n");fflush(stdout);
3359 if(NULL!=lk_left)//not a tip?
3360 {
3361 fprintf(stdout,"Partial Likelihoods on LEFT subtree of Branch %d [rate,site,state]:\n",b->num);
3362 for(catg=0;catg<tree->mod->ras->n_catg;++catg)
3363 for(site=0;site<tree->n_pattern;++site)
3364 for(j=0;j<tree->mod->ns;++j)
3365 #ifdef BEAGLE
3366 fprintf(stdout,fmt,catg,site,j,lk_left[catg*tree->n_pattern*tree->mod->ns + site*tree->mod->ns + j]);
3367 #else
3368 fprintf(stdout,fmt,catg,site,j,lk_left[catg*tree->mod->ns + site*tree->mod->ras->n_catg*tree->mod->ns + j]);
3369 #endif
3370 fflush(stdout);
3371 }
3372 else //is a tip
3373 {
3374 fprintf(stdout,"Likelihoods on LEFT tip of Branch %d [site,state]:\n",b->num);
3375 for(site=0;site<tree->n_pattern;++site)
3376 for(j=0;j<tree->mod->ns;++j)
3377 fprintf(stdout,"[%d,%d]%.1f ",site,j,b->p_lk_tip_l[site*tree->mod->ns + j]);
3378 fflush(stdout);
3379 }
3380 fprintf(stdout,"\n");fflush(stdout);
3381 if(NULL!=lk_right)//not a tip?
3382 {
3383 fprintf(stdout,"Partial Likelihoods on RIGHT subtree of Branch %d [rate,site,state]:\n",b->num);
3384 for(catg=0;catg<tree->mod->ras->n_catg;++catg)
3385 for(site=0;site<tree->n_pattern;++site)
3386 for(j=0;j<tree->mod->ns;++j)
3387 #ifdef BEAGLE
3388 fprintf(stdout,fmt,catg,site,j,lk_right[catg*tree->n_pattern*tree->mod->ns + site*tree->mod->ns + j]);
3389 #else
3390 fprintf(stdout,fmt,catg,site,j,lk_right[catg*tree->mod->ns + site*tree->mod->ras->n_catg*tree->mod->ns + j]);
3391 #endif
3392 fflush(stdout);
3393 }
3394 else //is a tip
3395 {
3396 fprintf(stdout,"Likelihoods on RIGHT tip of Branch %d [site,state]:\n",b->num);
3397 for(site=0;site<tree->n_pattern;++site)
3398 for(j=0;j<tree->mod->ns;++j)
3399 fprintf(stdout,"[%d,%d]%.1f ",site,j,b->p_lk_tip_r[site*tree->mod->ns + j]);
3400 fflush(stdout);
3401 }
3402 }
3403
3404 //////////////////////////////////////////////////////////////
3405 //////////////////////////////////////////////////////////////
3406
Print_All_Edge_Likelihoods(t_tree * tree)3407 void Print_All_Edge_Likelihoods(t_tree* tree)
3408 {
3409 int i;
3410 for(i=0;i < 2*tree->n_otu-3; ++i) Print_Edge_Likelihoods(tree, tree->a_edges[i],false);
3411 fflush(NULL);
3412 }
3413
3414 //////////////////////////////////////////////////////////////
3415 //////////////////////////////////////////////////////////////
3416
Dump_Arr_S(short int * arr,int num)3417 void Dump_Arr_S(short int* arr, int num)
3418 {
3419 int i;
3420
3421 if(NULL==arr)
3422 {
3423 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
3424 Exit("\n. Trying to print NULL array");
3425 return;
3426 }
3427 fprintf(stdout,"[");
3428 for(i=0;i<num;++i)
3429 {
3430 fprintf(stdout,"%d,",arr[i]);
3431 fflush(stdout);
3432 }
3433 fprintf(stdout,"]\n");fflush(stdout);
3434 }
3435
3436 //////////////////////////////////////////////////////////////
3437 //////////////////////////////////////////////////////////////
3438
Dump_Arr_D(phydbl * arr,int num)3439 void Dump_Arr_D(phydbl* arr, int num)
3440 {
3441 int i;
3442
3443 if(NULL==arr)
3444 {
3445 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
3446 Exit("\n. Trying to print NULL array");
3447 return;
3448 }
3449 fprintf(stdout,"[");
3450 for(i=0;i<num;++i)
3451 {
3452 fprintf(stdout,"%g,",arr[i]);
3453 fflush(stdout);
3454 }
3455 fprintf(stdout,"]\n");fflush(stdout);
3456 }
3457
3458 //////////////////////////////////////////////////////////////
3459 //////////////////////////////////////////////////////////////
3460
Dump_Arr_I(int * arr,int num)3461 void Dump_Arr_I(int* arr, int num)
3462 {
3463 int i;
3464
3465 if(NULL==arr)
3466 {
3467 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
3468 Exit("\n. Trying to print NULL array");
3469 return;
3470 }
3471 fprintf(stdout,"[");
3472 for(i=0;i<num;++i)
3473 {
3474 fprintf(stdout,"%d,",arr[i]);
3475 fflush(stdout);
3476 }
3477 fprintf(stdout,"]\n");fflush(stdout);
3478 }
3479
3480 //////////////////////////////////////////////////////////////
3481 //////////////////////////////////////////////////////////////
3482
Print_Tree_Structure(t_tree * tree)3483 void Print_Tree_Structure(t_tree* tree)
3484 {
3485 int i;
3486 PhyML_Fprintf(stdout,"\n. n_otu: %d",tree->n_otu);
3487
3488 for(i=0; i<2*tree->n_otu-1; ++i)
3489 {
3490 if(tree->a_edges[i])
3491 PhyML_Fprintf(stdout,"\n. Edge %p %3d, Length: %f LeftNode %3d [%s], RightNode %3d [%s]",
3492 tree->a_edges[i],
3493 tree->a_edges[i]->num,
3494 tree->a_edges[i]->l->v,
3495 tree->a_edges[i]->left ? tree->a_edges[i]->left->num : -1,
3496 tree->a_edges[i]->left ? tree->a_edges[i]->left->tax ? tree->a_edges[i]->left->name : "" : "",
3497 tree->a_edges[i]->rght ? tree->a_edges[i]->rght->num : -1,
3498 tree->a_edges[i]->left ? tree->a_edges[i]->rght->tax ? tree->a_edges[i]->rght->name : "" : "");
3499 else PhyML_Fprintf(stdout,"\n. NULL");
3500 }
3501
3502 for(i=0; i<2*tree->n_otu-1; ++i)
3503 {
3504 if(tree->a_nodes[i])
3505 PhyML_Fprintf(stdout,"\n. Node %p %3d v0: %3d v1: %3d v2: %3d b0: %3d b1: %3d b2: %3d",
3506 tree->a_nodes[i],
3507 tree->a_nodes[i]->num,
3508 tree->a_nodes[i]->v[0] ? tree->a_nodes[i]->v[0]->num : -1,
3509 tree->a_nodes[i]->v[1] ? tree->a_nodes[i]->v[1]->num : -1,
3510 tree->a_nodes[i]->v[2] ? tree->a_nodes[i]->v[2]->num : -1,
3511 tree->a_nodes[i]->v[0] ? tree->a_nodes[i]->b[0]->num : -1,
3512 tree->a_nodes[i]->v[1] ? tree->a_nodes[i]->b[1]->num : -1,
3513 tree->a_nodes[i]->v[2] ? tree->a_nodes[i]->b[2]->num : -1);
3514 else PhyML_Fprintf(stdout,"\n. NULL");
3515 }
3516 }
3517 //////////////////////////////////////////////////////////////
3518 //////////////////////////////////////////////////////////////
3519
Read_User_Tree(calign * cdata,t_mod * mod,option * io)3520 t_tree *Read_User_Tree(calign *cdata, t_mod *mod, option *io)
3521 {
3522 t_tree *tree;
3523
3524 PhyML_Printf("\n. Reading tree..."); fflush(NULL);
3525 if(io->n_trees == 1) rewind(io->fp_in_tree);
3526 tree = Read_Tree_File_Phylip(io->fp_in_tree);
3527 /* fclose(io->fp_in_tree); */
3528 /* io->fp_in_tree = NULL; */
3529 if(tree == NULL) Exit("\n. Input tree not found...");
3530 /* Add branch lengths if necessary */
3531 if(tree->has_branch_lengths == NO) Add_BioNJ_Branch_Lengths(tree,cdata,mod,NULL);
3532 return tree;
3533 }
3534
3535 //////////////////////////////////////////////////////////////
3536 //////////////////////////////////////////////////////////////
3537
3538
Print_Time_Info(time_t t_beg,time_t t_end)3539 void Print_Time_Info(time_t t_beg, time_t t_end)
3540 {
3541 div_t hour,min;
3542
3543 hour = div(t_end-t_beg,3600);
3544 min = div(t_end-t_beg,60 );
3545 min.quot -= hour.quot*60;
3546
3547 PhyML_Printf("\n\n. Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
3548 PhyML_Printf("\noooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3549 }
3550
3551 //////////////////////////////////////////////////////////////
3552 //////////////////////////////////////////////////////////////
3553
Print_Time_Info_Brief(time_t t_beg,time_t t_end)3554 void Print_Time_Info_Brief(time_t t_beg, time_t t_end)
3555 {
3556 div_t hour,min;
3557
3558 hour = div(t_end-t_beg,3600);
3559 min = div(t_end-t_beg,60 );
3560 min.quot -= hour.quot*60;
3561
3562 PhyML_Printf("\n. Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
3563 }
3564
3565 //////////////////////////////////////////////////////////////
3566 //////////////////////////////////////////////////////////////
3567
3568
PhyML_Printf(char * format,...)3569 void PhyML_Printf(char *format, ...)
3570 {
3571 va_list ptr;
3572
3573 #ifdef MPI
3574 if(Global_myRank == 0)
3575 {
3576 va_start (ptr, format);
3577 vprintf (format, ptr);
3578 va_end(ptr);
3579 }
3580 #else
3581 va_start (ptr, format);
3582 vprintf (format, ptr);
3583 va_end(ptr);
3584 #endif
3585
3586 fflush (NULL);
3587 }
3588
3589 //////////////////////////////////////////////////////////////
3590 //////////////////////////////////////////////////////////////
3591
PhyML_Fprintf(FILE * fp,char * format,...)3592 void PhyML_Fprintf(FILE *fp, char *format, ...)
3593 {
3594 va_list ptr;
3595
3596 #ifdef MPI
3597 if(Global_myRank == 0)
3598 {
3599 va_start (ptr, format);
3600 vfprintf (fp,format, ptr);
3601 va_end(ptr);
3602 }
3603 #else
3604 va_start (ptr, format);
3605 vfprintf (fp,format, ptr);
3606 va_end(ptr);
3607 #endif
3608
3609 fflush (NULL);
3610 }
3611
3612 //////////////////////////////////////////////////////////////
3613 //////////////////////////////////////////////////////////////
3614
PhyML_Fscanf(FILE * fp,char * format,...)3615 int PhyML_Fscanf(FILE *fp, char *format, ...)
3616 {
3617 va_list ptr;
3618 int rv;
3619
3620 rv = -1;
3621
3622 #ifdef MPI
3623 if(Global_myRank == 0)
3624 {
3625 va_start (ptr, format);
3626 rv = vfscanf(fp,format,ptr);
3627 if(!rv) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
3628 va_end(ptr);
3629 }
3630 #else
3631 va_start (ptr, format);
3632 rv = vfscanf(fp,format,ptr);
3633 if(!rv) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
3634 va_end(ptr);
3635 #endif
3636
3637 return(rv);
3638 }
3639
3640
3641 //////////////////////////////////////////////////////////////
3642 //////////////////////////////////////////////////////////////
3643
Read_Clade_Priors(char * file_name,t_tree * tree)3644 void Read_Clade_Priors(char *file_name, t_tree *tree)
3645 {
3646 FILE *fp;
3647 char *s,*line;
3648 int n_clade_priors;
3649 int clade_size;
3650 char **clade_list;
3651 int i,pos;
3652 phydbl prior_low,prior_up;
3653 int node_num;
3654
3655 PhyML_Printf("\n");
3656 PhyML_Printf("\n. Reading prior on node ages.\n");
3657
3658 line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3659 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3660
3661 clade_list = (char **)mCalloc(tree->n_otu,sizeof(char *));
3662 for(i=0;i<tree->n_otu;i++) clade_list[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
3663
3664 fp = Openfile(file_name,0);
3665
3666 n_clade_priors = 0;
3667 do
3668 {
3669 if(!fgets(line,T_MAX_LINE,fp)) break;
3670
3671 clade_size = 0;
3672 pos = 0;
3673 do
3674 {
3675 i = 0;
3676
3677 while(line[pos] == ' ') pos++;
3678
3679 while((line[pos] != ' ') && (line[pos] != '\n') && line[pos] != '#')
3680 {
3681 s[i] = line[pos];
3682 i++;
3683 pos++;
3684 }
3685 s[i] = '\0';
3686
3687 /* PhyML_Printf("\n. s = %s\n",s); */
3688
3689 if(line[pos] == '\n' || line[pos] == '#') break;
3690 pos++;
3691
3692 if(strcmp(s,"|"))
3693 {
3694 strcpy(clade_list[clade_size],s);
3695 clade_size++;
3696 }
3697 else
3698 break;
3699 }
3700 while(1);
3701
3702
3703 if(line[pos] != '#' && line[pos] != '\n')
3704 {
3705 phydbl val1, val2;
3706 /* sscanf(line+pos,"%lf %lf",&prior_up,&prior_low); */
3707 sscanf(line+pos,"%lf %lf",&val1,&val2);
3708 prior_up = (phydbl)val1;
3709 prior_low = (phydbl)val2;
3710 node_num = -1;
3711 if(!strcmp("@root@",clade_list[0])) node_num = tree->n_root->num;
3712 else node_num = Find_Clade(clade_list, clade_size, tree);
3713
3714 n_clade_priors++;
3715
3716 if(node_num < 0)
3717 {
3718 PhyML_Printf("\n");
3719 PhyML_Printf("\n");
3720 PhyML_Printf("\n. .................................................................");
3721 PhyML_Printf("\n. WARNING: could not find any clade in the tree referred to with the following taxon names:");
3722 for(i=0;i<clade_size;i++) PhyML_Printf("\n. \"%s\"",clade_list[i]);
3723 PhyML_Printf("\n. .................................................................");
3724 }
3725 else
3726 {
3727 tree->times->t_has_prior[node_num] = YES;
3728
3729 tree->times->t_prior_min[node_num] = -MAX(prior_low,prior_up);
3730 tree->times->t_prior_max[node_num] = -MIN(prior_low,prior_up);
3731
3732 if(fabs(prior_low - prior_up) < 1.E-6 && tree->a_nodes[node_num]->tax == YES)
3733 tree->times->nd_t[node_num] = prior_low;
3734
3735 PhyML_Printf("\n");
3736 PhyML_Printf("\n. [%3d]..................................................................",n_clade_priors);
3737 PhyML_Printf("\n. Node %4d matches the clade referred to with the following taxon names:",node_num);
3738 for(i=0;i<clade_size;i++) PhyML_Printf("\n. - \"%s\"",clade_list[i]);
3739 PhyML_Printf("\n. Lower bound set to: %15f time units.",MIN(prior_low,prior_up));
3740 PhyML_Printf("\n. Upper bound set to: %15f time units.",MAX(prior_low,prior_up));
3741 PhyML_Printf("\n. .......................................................................");
3742 }
3743 }
3744 }
3745 while(1);
3746
3747
3748 PhyML_Printf("\n. Read prior information on %d %s.",n_clade_priors,n_clade_priors > 1 ? "clades":"clade");
3749
3750 if(!n_clade_priors)
3751 {
3752 PhyML_Fprintf(stderr,"\n. PhyTime could not find any prior on node age.");
3753 PhyML_Fprintf(stderr,"\n. This is likely due to a problem in the calibration ");
3754 PhyML_Fprintf(stderr,"\n. file format. Make sure, for instance, that there is ");
3755 PhyML_Fprintf(stderr,"\n. a blank character between the end of the last name");
3756 PhyML_Fprintf(stderr,"\n. of each clade and the character `|'. Otherwise, ");
3757 PhyML_Fprintf(stderr,"\n. please refer to the example file.\n");
3758 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
3759 Warn_And_Exit("");
3760 }
3761
3762 for(i=0;i<tree->n_otu;i++) Free(clade_list[i]);
3763 Free(clade_list);
3764 Free(line);
3765 Free(s);
3766 fclose(fp);
3767 }
3768
3769 //////////////////////////////////////////////////////////////
3770 //////////////////////////////////////////////////////////////
3771
Get_Input(int argc,char ** argv)3772 option *Get_Input(int argc, char **argv)
3773 {
3774 option *io;
3775 t_mod *mod;
3776 t_opt *s_opt;
3777 int rv;
3778
3779 rv = 1;
3780
3781 io = (option *)Make_Input();
3782 mod = (t_mod *)Make_Model_Basic();
3783 s_opt = (t_opt *)Make_Optimiz();
3784
3785 Set_Defaults_Input(io);
3786 Set_Defaults_Model(mod);
3787 Set_Defaults_Optimiz(s_opt);
3788
3789 io->mod = mod;
3790 mod->io = io;
3791 mod->s_opt = s_opt;
3792
3793 #ifdef MPI
3794 rv = Read_Command_Line(io,argc,argv);
3795 #elif (defined PHYTIME || defined INVITEE || defined PHYREX || defined TEST)
3796 rv = Read_Command_Line(io,argc,argv);
3797 #else
3798 switch (argc)
3799 {
3800 case 1:
3801 {
3802 Launch_Interface(io);
3803 break;
3804 }
3805 default:
3806 {
3807 rv = Read_Command_Line(io,argc,argv);
3808 }
3809 }
3810 #endif
3811
3812
3813 if(rv) return io;
3814 else return NULL;
3815 }
3816
3817 //////////////////////////////////////////////////////////////
3818 //////////////////////////////////////////////////////////////
3819
Set_Whichmodel(int select)3820 int Set_Whichmodel(int select)
3821 {
3822 int wm;
3823
3824 wm = -1;
3825
3826 switch(select)
3827 {
3828 case 1:
3829 {
3830 wm = JC69;
3831 break;
3832 }
3833 case 2:
3834 {
3835 wm = K80;
3836 break;
3837 }
3838 case 3:
3839 {
3840 wm = F81;
3841 break;
3842 }
3843 case 4:
3844 {
3845 wm = HKY85;
3846 break;
3847 }
3848 case 5:
3849 {
3850 wm = F84;
3851 break;
3852 }
3853 case 6:
3854 {
3855 wm = TN93;
3856 break;
3857 }
3858 case 7:
3859 {
3860 wm = GTR;
3861 break;
3862 }
3863 case 8:
3864 {
3865 wm = CUSTOM;
3866 break;
3867 }
3868 case 11:
3869 {
3870 wm = WAG;
3871 break;
3872 }
3873 case 12:
3874 {
3875 wm = DAYHOFF;
3876 break;
3877 }
3878 case 13:
3879 {
3880 wm = JTT;
3881 break;
3882 }
3883 case 14:
3884 {
3885 wm = BLOSUM62;
3886 break;
3887 }
3888 case 15:
3889 {
3890 wm = MTREV;
3891 break;
3892 }
3893 case 16:
3894 {
3895 wm = RTREV;
3896 break;
3897 }
3898 case 17:
3899 {
3900 wm = CPREV;
3901 break;
3902 }
3903 case 18:
3904 {
3905 wm = DCMUT;
3906 break;
3907 }
3908 case 19:
3909 {
3910 wm = VT;
3911 break;
3912 }
3913 case 20:
3914 {
3915 wm = MTMAM;
3916 break;
3917 }
3918 case 21:
3919 {
3920 wm = MTART;
3921 break;
3922 }
3923 case 22:
3924 {
3925 wm = HIVW;
3926 break;
3927 }
3928 case 23:
3929 {
3930 wm = HIVB;
3931 break;
3932 }
3933 case 24:
3934 {
3935 wm = FLU;
3936 break;
3937 }
3938 case 25:
3939 {
3940 wm = CUSTOMAA;
3941 break;
3942 }
3943 case 26:
3944 {
3945 wm = LG;
3946 break;
3947 }
3948 case 27:
3949 {
3950 wm = AB;
3951 break;
3952 }
3953 default:
3954 {
3955 PhyML_Fprintf(stderr,"\n. Model number %d is unknown. Please use a valid model name",select);
3956 Exit("\n");
3957 break;
3958 }
3959 }
3960
3961 return wm;
3962
3963 }
3964
3965 //////////////////////////////////////////////////////////////
3966 //////////////////////////////////////////////////////////////
3967
Print_Data_Structure(int final,FILE * fp,t_tree * mixt_tree)3968 void Print_Data_Structure(int final, FILE *fp, t_tree *mixt_tree)
3969 {
3970 int n_partition_elem;
3971 char *s;
3972 t_tree *tree,*cpy_mixt_tree;
3973 int c,cc,cc_efrq,cc_rmat,cc_lens;
3974 char *param;
3975 int *link_efrq,*link_rmat,*link_lens;
3976 phydbl r_mat_weight_sum, e_frq_weight_sum;
3977
3978
3979 PhyML_Fprintf(fp,"\n. Starting tree: %s",
3980 mixt_tree->io->in_tree == 2?mixt_tree->io->in_tree_file:"BioNJ");
3981
3982 cpy_mixt_tree = mixt_tree;
3983
3984 n_partition_elem = 1;
3985 tree = mixt_tree;
3986 do
3987 {
3988 tree = tree->next_mixt;
3989 if(!tree) break;
3990 n_partition_elem++;
3991 }
3992 while(1);
3993
3994 s = (char *)mCalloc(2,sizeof(char));
3995 s[0] = ' ';
3996 s[1] = '\0';
3997 tree = mixt_tree;
3998 do
3999 {
4000 s = (char *)mRealloc(s,(int)(strlen(s)+strlen(tree->io->in_align_file)+2+2),sizeof(char));
4001 strcat(s,tree->io->in_align_file);
4002 strcat(s,", ");
4003 tree = tree->next_mixt;
4004 if(!tree) break;
4005 }
4006 while(1);
4007 s[(int)strlen(s)-2]=' ';
4008 s[(int)strlen(s)-1]='\0';
4009
4010 if(final == NO) PhyML_Fprintf(fp,"\n\n. Processing %d data %s (%s)",n_partition_elem,n_partition_elem>1?"sets":"set",s);
4011 if(final == YES) PhyML_Fprintf(fp,"\n\n. Processed %d data %s (%s)",n_partition_elem,n_partition_elem>1?"sets":"set",s);
4012 Free(s);
4013
4014 if(final == YES)
4015 PhyML_Fprintf(fp,"\n\n. Final log-likelihood: %f",mixt_tree->c_lnL);
4016
4017 r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->r_mat_weight);
4018 e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->e_frq_weight);
4019
4020 do
4021 {
4022 int class = 0;
4023
4024 PhyML_Fprintf(fp,"\n\n");
4025 PhyML_Fprintf(fp,"\n _______________________________________________________________________ ");
4026 PhyML_Fprintf(fp,"\n| |");
4027 PhyML_Fprintf(fp,"\n| %40s (partition element %2d) |",mixt_tree->io->in_align_file,mixt_tree->dp);
4028 PhyML_Fprintf(fp,"\n|_______________________________________________________________________|");
4029 PhyML_Fprintf(fp,"\n");
4030
4031 PhyML_Fprintf(fp,"\n. Number of rate classes:\t\t%20d",mixt_tree->mod->ras->n_catg+(mixt_tree->mod->ras->invar ?1:0));
4032 if(mixt_tree->mod->ras->n_catg > 1)
4033 {
4034 PhyML_Fprintf(fp,"\n. Model of rate variation:\t\t%20s",
4035 mixt_tree->mod->ras->free_mixt_rates?"FreeRates":
4036 mixt_tree->mod->ras->invar?"Gamma+Inv":"Gamma");
4037 if(mixt_tree->mod->ras->free_mixt_rates == NO)
4038 {
4039 PhyML_Fprintf(fp,"\n. Gamma shape parameter value:\t\t%20.2f",mixt_tree->mod->ras->alpha->v);
4040 PhyML_Fprintf(fp,"\n Optimise: \t\t\t\t%20s",mixt_tree->mod->s_opt->opt_alpha==YES?"yes":"no");
4041 }
4042 if(mixt_tree->mod->ras->invar == YES)
4043 {
4044 PhyML_Fprintf(fp,"\n. Proportion of invariable sites:\t%20.2f",mixt_tree->mod->ras->pinvar->v);
4045 PhyML_Fprintf(fp,"\n Optimise: \t\t\t\t%20s",mixt_tree->mod->s_opt->opt_pinvar==YES?"yes":"no");
4046 }
4047 }
4048 PhyML_Fprintf(fp,"\n. Relative average rate:\t\t%20f",mixt_tree->mod->br_len_mult->v);
4049
4050
4051 tree = mixt_tree;
4052 do
4053 {
4054 if(tree->is_mixt_tree) tree = tree->next;
4055
4056 PhyML_Fprintf(fp,"\n");
4057 PhyML_Fprintf(fp,"\n. Mixture class %d",class+1);
4058
4059 if(mixt_tree->mod->ras->n_catg > 1)
4060 {
4061 if(tree->mod->ras->invar == NO)
4062 {
4063 PhyML_Fprintf(fp,"\n Relative substitution rate:\t%20f",mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number]);
4064 PhyML_Fprintf(fp,"\n Rel. rate freq. (> 0 rates):\t%20f",mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number]);
4065 PhyML_Fprintf(fp,"\n Rate class number:\t\t%20d",tree->mod->ras->parent_class_number);
4066 }
4067 else
4068 {
4069 PhyML_Fprintf(fp,"\n Relative substitution rate:\t%20f",0.0);
4070 PhyML_Fprintf(fp,"\n Relative rate freq.:\t\t%20f",mixt_tree->mod->ras->pinvar->v);
4071 }
4072 }
4073
4074 PhyML_Fprintf(fp,"\n Substitution model:\t\t%20s",tree->mod->modelname->s);
4075
4076 if(tree->mod->whichmodel == CUSTOM)
4077 PhyML_Fprintf(fp,"\n Substitution model code:\t%20s",tree->mod->custom_mod_string->s);
4078
4079 if(tree->mod->whichmodel == CUSTOMAA)
4080 PhyML_Fprintf(fp,"\n Rate matrix file name:\t%20s",tree->mod->aa_rate_mat_file->s);
4081
4082 if(tree->mod->whichmodel == K80 ||
4083 tree->mod->whichmodel == HKY85 ||
4084 tree->mod->whichmodel == TN93)
4085 {
4086 PhyML_Fprintf(fp,"\n Value of the ts/tv ratio:\t%20f",tree->mod->kappa->v);
4087 PhyML_Fprintf(fp,"\n Optimise ts/tv ratio:\t%20s",tree->mod->s_opt->opt_kappa?"yes":"no");
4088 }
4089 else if(tree->mod->whichmodel == GTR ||
4090 tree->mod->whichmodel == CUSTOM)
4091 {
4092 PhyML_Fprintf(fp,"\n Optimise subst. rates:\t%20s",tree->mod->s_opt->opt_rr?"yes":"no");
4093 if(final == YES)
4094 {
4095 PhyML_Fprintf(fp,"\n Subst. rate A<->C:\t\t%20.2f",tree->mod->r_mat->rr->v[0]);
4096 PhyML_Fprintf(fp,"\n Subst. rate A<->G:\t\t%20.2f",tree->mod->r_mat->rr->v[1]);
4097 PhyML_Fprintf(fp,"\n Subst. rate A<->T:\t\t%20.2f",tree->mod->r_mat->rr->v[2]);
4098 PhyML_Fprintf(fp,"\n Subst. rate C<->G:\t\t%20.2f",tree->mod->r_mat->rr->v[3]);
4099 PhyML_Fprintf(fp,"\n Subst. rate C<->T:\t\t%20.2f",tree->mod->r_mat->rr->v[4]);
4100 PhyML_Fprintf(fp,"\n Subst. rate G<->T:\t\t%20.2f",tree->mod->r_mat->rr->v[5]);
4101 }
4102 }
4103
4104 PhyML_Fprintf(fp,"\n Rate matrix weight:\t\t%20f",tree->mod->r_mat_weight->v / r_mat_weight_sum);
4105
4106 if(tree->io->datatype == NT &&
4107 tree->mod->whichmodel != JC69 &&
4108 tree->mod->whichmodel != K80)
4109 {
4110 PhyML_Fprintf(fp,"\n Optimise nucleotide freq.:\t%20s",tree->mod->s_opt->opt_state_freq?"yes":"no");
4111 if(final == YES)
4112 {
4113 PhyML_Fprintf(fp,"\n Freq(A):\t\t\t%20.2f",tree->mod->e_frq->pi->v[0]);
4114 PhyML_Fprintf(fp,"\n Freq(C):\t\t\t%20.2f",tree->mod->e_frq->pi->v[1]);
4115 PhyML_Fprintf(fp,"\n Freq(G):\t\t\t%20.2f",tree->mod->e_frq->pi->v[2]);
4116 PhyML_Fprintf(fp,"\n Freq(T):\t\t\t%20.2f",tree->mod->e_frq->pi->v[3]);
4117 }
4118 }
4119 else if(tree->io->datatype == AA)
4120 {
4121 char *s;
4122
4123 s = (char *)mCalloc(50,sizeof(char));
4124
4125 if(tree->mod->s_opt->opt_state_freq == YES)
4126 {
4127 strcpy(s,"Empirical");
4128 }
4129 else
4130 {
4131 strcpy(s,"Model");
4132 }
4133
4134 PhyML_Fprintf(fp,"\n Amino-acid freq.:\t\t%20s",s);
4135
4136 Free(s);
4137 }
4138
4139
4140 PhyML_Fprintf(fp,"\n Equ. freq. weight:\t\t%20f",tree->mod->e_frq_weight->v / e_frq_weight_sum);
4141
4142 class++;
4143
4144 tree = tree->next;
4145
4146 if(tree && tree->is_mixt_tree == YES) break;
4147 }
4148 while(tree);
4149
4150 mixt_tree = mixt_tree->next_mixt;
4151 if(!mixt_tree) break;
4152 }
4153 while(1);
4154
4155
4156 tree = cpy_mixt_tree;
4157 c = 0;
4158 do
4159 {
4160 if(tree->is_mixt_tree) tree = tree->next;
4161 c++;
4162 tree = tree->next;
4163 }
4164 while(tree);
4165
4166 link_efrq = (int *)mCalloc(c,sizeof(int));
4167 link_lens = (int *)mCalloc(c,sizeof(int));
4168 link_rmat = (int *)mCalloc(c,sizeof(int));
4169
4170 PhyML_Fprintf(fp,"\n");
4171 PhyML_Fprintf(fp,"\n");
4172 PhyML_Fprintf(fp,"\n _______________________________________________________________________ ");
4173 PhyML_Fprintf(fp,"\n| |");
4174 PhyML_Fprintf(fp,"\n| Model summary table |");
4175 PhyML_Fprintf(fp,"\n|_______________________________________________________________________|");
4176 PhyML_Fprintf(fp,"\n");
4177 PhyML_Fprintf(fp,"\n");
4178 /* PhyML_Fprintf(fp," "); */
4179 PhyML_Fprintf(fp," ------------------");
4180 tree = cpy_mixt_tree;
4181 c = 0;
4182 do
4183 {
4184 if(tree->is_mixt_tree) tree = tree->next;
4185 PhyML_Fprintf(fp,"---");
4186 tree = tree->next;
4187 }
4188 while(tree);
4189
4190 param = (char *)mCalloc(30,sizeof(char));
4191 PhyML_Fprintf(fp,"\n");
4192 strcpy(param,"Partition element ");
4193 PhyML_Fprintf(fp," %18s",param);
4194 tree = cpy_mixt_tree;
4195 c = 0;
4196 do
4197 {
4198 if(tree->is_mixt_tree) tree = tree->next;
4199 PhyML_Fprintf(fp,"%2d ",tree->mixt_tree->dp);
4200 tree = tree->next;
4201 }
4202 while(tree);
4203
4204
4205 PhyML_Fprintf(fp,"\n");
4206 PhyML_Fprintf(fp," ------------------");
4207 tree = cpy_mixt_tree;
4208 c = 0;
4209 do
4210 {
4211 if(tree->is_mixt_tree) tree = tree->next;
4212 PhyML_Fprintf(fp,"---");
4213 tree = tree->next;
4214 }
4215 while(tree);
4216
4217
4218 tree = cpy_mixt_tree;
4219 c = 0;
4220 do
4221 {
4222 if(tree->is_mixt_tree) tree = tree->next;
4223 link_rmat[c] = -1;
4224 link_lens[c] = -1;
4225 link_efrq[c] = -1;
4226 tree = tree->next;
4227 c++;
4228 }
4229 while(tree);
4230
4231 mixt_tree = cpy_mixt_tree;
4232 cc_efrq = 97;
4233 cc_rmat = 97;
4234 cc_lens = 97;
4235 cc = 0;
4236 do
4237 {
4238 if(mixt_tree->is_mixt_tree) mixt_tree = mixt_tree->next;
4239
4240 if(link_efrq[cc] < 0)
4241 {
4242 link_efrq[cc] = cc_efrq;
4243 tree = mixt_tree->next;
4244 c = cc+1;
4245 if(tree)
4246 {
4247 do
4248 {
4249 if(tree->is_mixt_tree) tree = tree->next;
4250
4251 if(mixt_tree->mod->e_frq == tree->mod->e_frq) link_efrq[c] = cc_efrq;
4252
4253 tree = tree->next;
4254 c++;
4255 }
4256 while(tree);
4257 }
4258 cc_efrq++;
4259 }
4260
4261 if(link_lens[cc] < 0)
4262 {
4263 link_lens[cc] = cc_lens;
4264 tree = mixt_tree->next;
4265 c = cc+1;
4266 if(tree)
4267 {
4268 do
4269 {
4270 if(tree->is_mixt_tree) tree = tree->next;
4271
4272 if(mixt_tree->a_edges[0]->l == tree->a_edges[0]->l) link_lens[c] = cc_lens;
4273
4274 tree = tree->next;
4275 c++;
4276 }
4277 while(tree);
4278 }
4279 cc_lens++;
4280 }
4281
4282 if(link_rmat[cc] < 0)
4283 {
4284 link_rmat[cc] = cc_rmat;
4285 tree = mixt_tree->next;
4286 c = cc+1;
4287 if(tree)
4288 {
4289 do
4290 {
4291 if(tree->is_mixt_tree) tree = tree->next;
4292
4293 if(mixt_tree->mod->r_mat == tree->mod->r_mat &&
4294 mixt_tree->mod->whichmodel == tree->mod->whichmodel &&
4295 !strcmp(mixt_tree->mod->custom_mod_string->s,
4296 tree->mod->custom_mod_string->s) &&
4297 !strcmp(mixt_tree->mod->aa_rate_mat_file->s,
4298 tree->mod->aa_rate_mat_file->s)) link_rmat[c] = cc_rmat;
4299
4300 tree = tree->next;
4301 c++;
4302 }
4303 while(tree);
4304 }
4305 cc_rmat++;
4306 }
4307
4308 cc++;
4309 mixt_tree = mixt_tree->next;
4310 }
4311 while(mixt_tree);
4312
4313 PhyML_Fprintf(fp,"\n");
4314 strcpy(param,"State frequencies ");
4315 PhyML_Fprintf(fp," %18s",param);
4316
4317 tree = cpy_mixt_tree;
4318 c = 0;
4319 do
4320 {
4321 if(tree->is_mixt_tree) tree = tree->next;
4322 PhyML_Fprintf(fp,"%2c ",link_efrq[c]);
4323 tree = tree->next;
4324 c++;
4325 }
4326 while(tree);
4327
4328
4329 PhyML_Fprintf(fp,"\n");
4330 strcpy(param,"Branch lengths ");
4331 PhyML_Fprintf(fp," %18s",param);
4332
4333 tree = cpy_mixt_tree;
4334 c = 0;
4335 do
4336 {
4337 if(tree->is_mixt_tree) tree = tree->next;
4338 PhyML_Fprintf(fp,"%2c ",link_lens[c]);
4339 tree = tree->next;
4340 c++;
4341 }
4342 while(tree);
4343
4344 PhyML_Fprintf(fp,"\n");
4345 strcpy(param,"Rate matrix ");
4346 PhyML_Fprintf(fp," %18s",param);
4347
4348 tree = cpy_mixt_tree;
4349 c = 0;
4350 do
4351 {
4352 if(tree->is_mixt_tree) tree = tree->next;
4353 PhyML_Fprintf(fp,"%2c ",link_rmat[c]);
4354 tree = tree->next;
4355 c++;
4356 }
4357 while(tree);
4358
4359 PhyML_Fprintf(fp,"\n");
4360 PhyML_Fprintf(fp," ------------------");
4361 tree = cpy_mixt_tree;
4362 c = 0;
4363 do
4364 {
4365 if(tree->is_mixt_tree) tree = tree->next;
4366 PhyML_Fprintf(fp,"---");
4367 tree = tree->next;
4368 }
4369 while(tree);
4370 PhyML_Fprintf(fp,"\n");
4371
4372 if(final == YES)
4373 {
4374 tree = cpy_mixt_tree;
4375 c = 0;
4376 do
4377 {
4378 PhyML_Fprintf(fp,"\n");
4379 PhyML_Fprintf(fp,"\n. Tree estimated from data partition %d",c++);
4380 Br_Len_Involving_Invar(tree->next);
4381 Rescale_Br_Len_Multiplier_Tree(tree->next);
4382 s = Write_Tree(tree->next); /*! tree->next is not a mixt_tree so edge lengths
4383 are not averaged over when writing the tree out. */
4384 PhyML_Fprintf(fp,"\n %s",s);
4385 Br_Len_Not_Involving_Invar(tree->next);
4386 Unscale_Br_Len_Multiplier_Tree(tree->next);
4387 Free(s);
4388 tree = tree->next_mixt;
4389 }
4390 while(tree);
4391 }
4392
4393 Free(param);
4394 Free(link_efrq);
4395 Free(link_rmat);
4396 Free(link_lens);
4397
4398 mixt_tree = cpy_mixt_tree;
4399 }
4400
4401 //////////////////////////////////////////////////////////////
4402 //////////////////////////////////////////////////////////////
4403
PhyML_XML(char * xml_filename)4404 option *PhyML_XML(char *xml_filename)
4405 {
4406 char *most_likely_tree;
4407 phydbl best_lnL;
4408 int num_rand_tree;
4409 t_tree *mixt_tree,*tree;
4410 option *io,*dum;
4411 xml_node *xml_root;
4412
4413 mixt_tree = XML_Process_Base(xml_filename);
4414 io = mixt_tree->io;
4415 most_likely_tree = NULL;
4416 xml_root = mixt_tree->xml_root;
4417 best_lnL = UNLIKELY;
4418 dum = NULL;
4419
4420 dum = (option *)mCalloc(1,sizeof(option));
4421 dum->use_xml = YES;
4422
4423 for(num_rand_tree=0;num_rand_tree<io->mod->s_opt->n_rand_starts;num_rand_tree++)
4424 {
4425 MIXT_Check_Model_Validity(mixt_tree);
4426 MIXT_Init_Model(mixt_tree);
4427 Print_Data_Structure(NO,stdout,mixt_tree);
4428 tree = MIXT_Starting_Tree(mixt_tree);
4429 Copy_Tree(tree,mixt_tree);
4430 Free_Tree(tree);
4431
4432 if(mixt_tree->io->mod->s_opt->random_input_tree)
4433 {
4434 PhyML_Printf("\n\n. [%3d/%3d]",num_rand_tree+1,mixt_tree->io->mod->s_opt->n_rand_starts);
4435 Random_Tree(mixt_tree);
4436 }
4437
4438 MIXT_Connect_Cseqs_To_Nodes(mixt_tree);
4439 MIXT_Init_T_Beg(mixt_tree);
4440
4441 MIXT_Make_Tree_For_Pars(mixt_tree);
4442 MIXT_Make_Tree_For_Lk(mixt_tree);
4443 MIXT_Make_Spr(mixt_tree);
4444
4445 MIXT_Chain_All(mixt_tree);
4446 MIXT_Check_Edge_Lens_In_All_Elem(mixt_tree);
4447 MIXT_Turn_Branches_OnOff_In_All_Elem(ON,mixt_tree);
4448 MIXT_Check_Invar_Struct_In_Each_Partition_Elem(mixt_tree);
4449 MIXT_Check_RAS_Struct_In_Each_Partition_Elem(mixt_tree);
4450 Br_Len_Not_Involving_Invar(mixt_tree);
4451 Unscale_Br_Len_Multiplier_Tree(mixt_tree);
4452 Set_Both_Sides(YES,mixt_tree);
4453
4454 Set_Update_Eigen(YES,mixt_tree->mod);
4455 Lk(NULL,mixt_tree);
4456 Set_Update_Eigen(NO,mixt_tree->mod);
4457
4458
4459 if(mixt_tree->mod->s_opt->opt_topo)
4460 {
4461 Global_Spr_Search(mixt_tree);
4462 }
4463 else
4464 {
4465 Round_Optimize(mixt_tree,ROUND_MAX);
4466 }
4467
4468
4469 PhyML_Printf("\n\n. Log-likelihood = %f",mixt_tree->c_lnL);
4470
4471 if((num_rand_tree == io->mod->s_opt->n_rand_starts-1) && (io->mod->s_opt->random_input_tree))
4472 {
4473 num_rand_tree--;
4474 io->mod->s_opt->random_input_tree = NO;
4475 }
4476
4477 Br_Len_Involving_Invar(mixt_tree);
4478 Rescale_Br_Len_Multiplier_Tree(mixt_tree);
4479
4480 /*! Print the tree estimated using the current random (or BioNJ) starting tree */
4481 if(io->mod->s_opt->n_rand_starts > 1)
4482 {
4483 Print_Tree(io->fp_out_trees,mixt_tree);
4484 fflush(NULL);
4485 }
4486
4487 if(mixt_tree->c_lnL > best_lnL)
4488 {
4489 best_lnL = mixt_tree->c_lnL;
4490 if(most_likely_tree) Free(most_likely_tree);
4491 if(io->ratio_test) aLRT(mixt_tree);
4492 most_likely_tree = Write_Tree(mixt_tree);
4493 mixt_tree->lock_topo = NO;
4494 }
4495
4496 tree = mixt_tree;
4497 do
4498 {
4499 if(tree->io->print_site_lnl == YES) Print_Site_Lk(tree,tree->io->fp_out_lk);
4500 tree = tree->next_mixt;
4501 }
4502 while(tree);
4503
4504 MIXT_Init_T_End(mixt_tree);
4505 Print_Data_Structure(YES,mixt_tree->io->fp_out_stats,mixt_tree);
4506
4507 Free_Spr_List_One_Edge(mixt_tree);
4508 Free_Tree_Pars(mixt_tree);
4509 Free_Tree_Lk(mixt_tree);
4510 }
4511
4512 /*! Print the most likely tree in the output file */
4513 if(!mixt_tree->io->quiet) PhyML_Printf("\n\n. Printing the most likely tree in file '%s'...\n", Basename(mixt_tree->io->out_tree_file));
4514 PhyML_Fprintf(mixt_tree->io->fp_out_tree,"%s\n",most_likely_tree);
4515
4516 /*! Bootstrap analysis */
4517 MIXT_Bootstrap(most_likely_tree,xml_root);
4518
4519 while(io->prev != NULL) io = io->prev;
4520
4521 Free(most_likely_tree);
4522
4523 tree = mixt_tree;
4524 do
4525 {
4526 Free_Calign(tree->data);
4527 tree = tree->next_mixt;
4528 }
4529 while(tree);
4530
4531 tree = mixt_tree;
4532 do
4533 {
4534 Free_Optimiz(tree->mod->s_opt);
4535 tree = tree->next;
4536 }
4537 while(tree);
4538
4539 Free_Model_Complete(mixt_tree->mod);
4540 Free_Model_Basic(mixt_tree->mod);
4541 Free_Tree(mixt_tree);
4542
4543 if(io->fp_out_trees) fclose(io->fp_out_trees);
4544 if(io->fp_out_tree) fclose(io->fp_out_tree);
4545 if(io->fp_out_stats) fclose(io->fp_out_stats);
4546 if(io->fp_out_json_trace) fclose(io->fp_out_json_trace);
4547 Free_Input(io);
4548
4549 XML_Free_XML_Tree(xml_root);
4550
4551 return(dum);
4552 }
4553
4554 //////////////////////////////////////////////////////////////
4555 //////////////////////////////////////////////////////////////
4556
4557 /*! Check that the same nodes in the different mixt_trees are
4558 connected to the same taxa
4559 */
Check_Taxa_Sets(t_tree * mixt_tree)4560 void Check_Taxa_Sets(t_tree *mixt_tree)
4561 {
4562 t_tree *tree;
4563 int i;
4564
4565 tree = mixt_tree;
4566 do
4567 {
4568 if(tree->next)
4569 {
4570 for(i=0;i<tree->n_otu;i++)
4571 {
4572 if(strcmp(tree->a_nodes[i]->name,tree->next->a_nodes[i]->name))
4573 {
4574 PhyML_Fprintf(stderr,"\n. There seems to be a problem in one (or more) of your");
4575 PhyML_Fprintf(stderr,"\n. sequence alignments. PhyML could not match taxon");
4576 PhyML_Fprintf(stderr,"\n. '%s' found in file '%s' with any of the taxa",tree->a_nodes[i]->name,tree->io->in_align_file);
4577 PhyML_Fprintf(stderr,"\n. listed in file '%s'.",tree->next->io->in_align_file);
4578 Exit("\n");
4579 }
4580 }
4581 }
4582 tree = tree->next;
4583 }
4584 while(tree);
4585
4586 }
4587
4588 //////////////////////////////////////////////////////////////
4589 //////////////////////////////////////////////////////////////
4590
Make_Ratematrix_From_XML_Node(xml_node * instance,option * io,t_mod * mod)4591 void Make_Ratematrix_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
4592 {
4593 char *model = NULL;
4594 int select;
4595
4596 model = XML_Get_Attribute_Value(instance,"model");
4597
4598 if(model == NULL)
4599 {
4600 PhyML_Fprintf(stderr,"\n. Poorly formated XML file.");
4601 PhyML_Fprintf(stderr,"\n. Attribute 'model' is mandatory in a <ratematrix> node.");
4602 Exit("\n");
4603 }
4604
4605 select = XML_Validate_Attr_Int(model,27,
4606 "xxxxx", //0
4607 "JC69", //1
4608 "K80", //2
4609 "F81", //3
4610 "HKY85", //4
4611 "F84", //5
4612 "TN93", //6
4613 "GTR", //7
4614 "CUSTOM", //8
4615 "xxxxx", //9
4616 "xxxxx", //10
4617 "WAG", //11
4618 "DAYHOFF", //12
4619 "JTT", //13
4620 "BLOSUM62", //14
4621 "MTREV", //15
4622 "RTREV", //16
4623 "CPREV", //17
4624 "DCMUT", //18
4625 "VT", //19
4626 "MTMAM", //20
4627 "MTART", //21
4628 "HIVW", //22
4629 "HIVB", //23
4630 "FLU", //24
4631 "CUSTOMAA", //25
4632 "LG", //26
4633 "AB" ); //27
4634
4635 if(select < 9)
4636 {
4637 if(io->datatype != NT)
4638 {
4639 PhyML_Fprintf(stderr,"\n. Data type and selected model are incompatible");
4640 Exit("\n");
4641 }
4642 }
4643 else
4644 {
4645 if(io->datatype != AA)
4646 {
4647 PhyML_Fprintf(stderr,"\n. Data type and selected model are incompatible");
4648 Exit("\n");
4649 }
4650 }
4651
4652 mod->r_mat = (t_rmat *)Make_Rmat(mod->ns);
4653 Init_Rmat(mod->r_mat);
4654
4655 /*! Set model number & name */
4656 mod->whichmodel = Set_Whichmodel(select);
4657 Set_Model_Name(mod);
4658
4659 if(mod->whichmodel == K80 ||
4660 mod->whichmodel == HKY85 ||
4661 mod->whichmodel == TN93)
4662 {
4663 char *tstv,*opt_tstv;
4664
4665 tstv = XML_Get_Attribute_Value(instance,"tstv");
4666
4667 if(tstv)
4668 {
4669 mod->s_opt->opt_kappa = NO;
4670 mod->kappa->v = String_To_Dbl(tstv);
4671 }
4672 else
4673 {
4674 mod->s_opt->opt_kappa = YES;
4675 }
4676
4677 opt_tstv = XML_Get_Attribute_Value(instance,"optimise.tstv");
4678
4679 if(opt_tstv)
4680 {
4681 if(!strcmp(opt_tstv,"true") || !strcmp(opt_tstv,"yes"))
4682 {
4683 mod->s_opt->opt_kappa = YES;
4684 mod->s_opt->opt_subst_param = YES;
4685 }
4686 else
4687 {
4688 mod->s_opt->opt_kappa = NO;
4689 mod->s_opt->opt_subst_param = NO;
4690 }
4691 }
4692 }
4693 else
4694 {
4695 mod->s_opt->opt_kappa = NO;
4696 }
4697
4698 if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM)
4699 {
4700 char *opt_rr;
4701
4702 opt_rr = XML_Get_Attribute_Value(instance,"optimise.rr");
4703
4704
4705 if(opt_rr)
4706 {
4707 if(!strcmp(opt_rr,"yes") || !strcmp(opt_rr,"true"))
4708 {
4709 mod->s_opt->opt_rr = YES;
4710 mod->s_opt->opt_subst_param = YES;
4711 }
4712 else
4713 {
4714 mod->s_opt->opt_rr = NO;
4715 mod->s_opt->opt_subst_param = NO;
4716 }
4717 }
4718 }
4719
4720 if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM)
4721 {
4722 xml_node *rr;
4723 char *val;
4724 phydbl v;
4725
4726 rr = XML_Search_Node_Name("rr",YES,instance);
4727
4728 if(rr != NULL)
4729 {
4730 mod->r_mat = (t_rmat *)Make_Rmat(mod->ns);
4731 Init_Rmat(mod->r_mat);
4732 Make_Custom_Model(mod);
4733
4734
4735 // A<->C
4736 val = XML_Get_Attribute_Value(rr,"AC");
4737 if(val == NULL)
4738 {
4739 PhyML_Printf("\n. Please specify the relative rate of substitution between A and T");
4740 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4741 }
4742 else
4743 {
4744 v = strtod(val,NULL);
4745 if(v != 0 && v > 0.0)
4746 {
4747 mod->r_mat->rr->v[0] = v;
4748 mod->r_mat->rr_val->v[0] = log(v);
4749 }
4750 else
4751 {
4752 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4753 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4754 }
4755 }
4756
4757 // A<->G
4758 val = XML_Get_Attribute_Value(rr,"AG");
4759 if(val == NULL)
4760 {
4761 PhyML_Printf("\n. Please specify the relative rate of substitution between A and T");
4762 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4763 }
4764 else
4765 {
4766 v = strtod(val,NULL);
4767 if(v != 0 && v > 0.0)
4768 {
4769 mod->r_mat->rr->v[1] = v;
4770 mod->r_mat->rr_val->v[1] = log(v);
4771 }
4772 else
4773 {
4774 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4775 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4776 }
4777 }
4778
4779 // A<->T
4780 val = XML_Get_Attribute_Value(rr,"AT");
4781 if(val == NULL)
4782 {
4783 PhyML_Printf("\n. Please specify the relative rate of substitution between A and T");
4784 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4785 }
4786 else
4787 {
4788 v = strtod(val,NULL);
4789 if(v != 0 && v > 0.0)
4790 {
4791 mod->r_mat->rr->v[2] = v;
4792 mod->r_mat->rr_val->v[2] = log(v);
4793 }
4794 else
4795 {
4796 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4797 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4798 }
4799 }
4800
4801 // C<->G
4802 val = XML_Get_Attribute_Value(rr,"CG");
4803 if(val == NULL)
4804 {
4805 PhyML_Printf("\n. Please specify the relative rate of substitution between A and T");
4806 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4807 }
4808 else
4809 {
4810 v = strtod(val,NULL);
4811 if(v != 0 && v > 0.0)
4812 {
4813 mod->r_mat->rr->v[3] = v;
4814 mod->r_mat->rr_val->v[3] = log(v);
4815 }
4816 else
4817 {
4818 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4819 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4820 }
4821 }
4822
4823
4824 // C<->T
4825 val = XML_Get_Attribute_Value(rr,"CT");
4826 if(val == NULL)
4827 {
4828 PhyML_Printf("\n. Please specify the relative rate of substitution between A and T");
4829 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4830 }
4831 else
4832 {
4833 v = strtod(val,NULL);
4834 if(v != 0 && v > 0.0)
4835 {
4836 mod->r_mat->rr->v[4] = v;
4837 mod->r_mat->rr_val->v[4] = log(v);
4838 }
4839 else
4840 {
4841 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4842 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4843 }
4844 }
4845
4846 // G<->T
4847 val = XML_Get_Attribute_Value(rr,"GT");
4848
4849 if(val != NULL)
4850 {
4851 v = strtod(val,NULL);
4852 if(v != 0 && v > 0.0)
4853 {
4854 mod->r_mat->rr->v[5] = v;
4855 mod->r_mat->rr_val->v[5] = log(v);
4856 }
4857 else
4858 {
4859 PhyML_Printf("\n. Invalid relative rate parameter value: '%s'.\n", val);
4860 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
4861 }
4862 }
4863 }
4864 }
4865
4866
4867 /*! Custom model for nucleotide sequences. Read the corresponding
4868 code. */
4869 if(mod->whichmodel == CUSTOM)
4870 {
4871 char *model_code = NULL;
4872
4873 model_code = XML_Get_Attribute_Value(instance,"model.code");
4874
4875 if(!model_code)
4876 {
4877 PhyML_Fprintf(stderr,"\n. No valid 'model.code' attribute could be found.\n");
4878 PhyML_Fprintf(stderr,"\n. Please fix your XML file.\n");
4879 Exit("\n");
4880 }
4881 else
4882 {
4883 strcpy(mod->custom_mod_string->s,model_code);
4884 }
4885 }
4886
4887
4888 /*! Custom model for amino-acids. Read in the rate matrix file */
4889 if(mod->whichmodel == CUSTOMAA)
4890 {
4891 char *r_mat_file;
4892
4893 r_mat_file = XML_Get_Attribute_Value(instance,"ratematrix.file");
4894
4895 if(!r_mat_file)
4896 {
4897 PhyML_Fprintf(stderr,"\n. No valid 'ratematrix.file' attribute could be found.");
4898 PhyML_Fprintf(stderr,"\n. Please fix your XML file.\n");
4899 Exit("\n");
4900 }
4901 else
4902 {
4903 strcpy(mod->aa_rate_mat_file->s,r_mat_file);
4904 }
4905
4906 /* Free(r_mat_file); */
4907 }
4908
4909 char *buff;
4910
4911 buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
4912 if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
4913 {
4914 mod->s_opt->opt_rmat_weight = YES;
4915 }
4916 else
4917 {
4918 mod->s_opt->opt_rmat_weight = NO;
4919 }
4920 }
4921
4922 //////////////////////////////////////////////////////////////
4923 //////////////////////////////////////////////////////////////
4924
Make_Efrq_From_XML_Node(xml_node * instance,option * io,t_mod * mod)4925 void Make_Efrq_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
4926 {
4927 char *buff;
4928
4929 mod->e_frq = (t_efrq *)Make_Efrq(mod->ns);
4930 Init_Efrq(NULL,mod->e_frq);
4931
4932 buff = XML_Get_Attribute_Value(instance,"optimise.freqs");
4933
4934 if(buff)
4935 {
4936 if(!strcmp(buff,"yes") || !strcmp(buff,"true"))
4937 {
4938 if(io->datatype == AA)
4939 {
4940 PhyML_Fprintf(stderr,"\n. Option 'optimise.freqs' set to 'yes' (or 'true')");
4941 PhyML_Fprintf(stderr,"\n. is not allowed with amino-acid data.");
4942 Exit("\n");
4943 }
4944 mod->s_opt->opt_state_freq = YES;
4945 }
4946 }
4947
4948 buff = XML_Get_Attribute_Value(instance,"aa.freqs");
4949
4950 if(buff)
4951 {
4952 if(!strcmp(buff,"empirical"))
4953 {
4954 if(io->datatype == AA)
4955 {
4956 mod->s_opt->opt_state_freq = YES;
4957 mod->e_frq->empirical_state_freq = YES;
4958 }
4959 else if(io->datatype == NT)
4960 {
4961 mod->s_opt->opt_state_freq = NO;
4962 }
4963 }
4964 }
4965
4966
4967 buff = XML_Get_Attribute_Value(instance,"base.freqs");
4968
4969 if(buff)
4970 {
4971 if(io->datatype == AA)
4972 {
4973 PhyML_Fprintf(stderr,"\n. Option 'base.freqs' is not allowed with amino-acid data.");
4974 Exit("\n");
4975 }
4976 else
4977 {
4978 phydbl A,C,G,T;
4979 sscanf(buff,"%lf,%lf,%lf,%lf",&A,&C,&G,&T);
4980 mod->e_frq->user_b_freq->v[0] = (phydbl)A;
4981 mod->e_frq->user_b_freq->v[1] = (phydbl)C;
4982 mod->e_frq->user_b_freq->v[2] = (phydbl)G;
4983 mod->e_frq->user_b_freq->v[3] = (phydbl)T;
4984 mod->e_frq->user_state_freq = YES;
4985 mod->s_opt->opt_state_freq = NO;
4986 }
4987 }
4988
4989
4990 buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
4991 if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
4992 {
4993 mod->s_opt->opt_efrq_weight = YES;
4994 }
4995 else
4996 {
4997 mod->s_opt->opt_efrq_weight = NO;
4998 }
4999 }
5000
5001 //////////////////////////////////////////////////////////////
5002 //////////////////////////////////////////////////////////////
5003
Make_Topology_From_XML_Node(xml_node * instance,option * io,t_mod * mod)5004 void Make_Topology_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5005 {
5006 // Starting tree
5007 char *init_tree;
5008
5009 init_tree = XML_Get_Attribute_Value(instance,"init.tree");
5010
5011 if(!init_tree)
5012 {
5013 PhyML_Fprintf(stderr,"\n. Attribute 'init.tree=bionj|user|random' is mandatory");
5014 PhyML_Fprintf(stderr,"\n. Please amend your XML file accordingly.");
5015 Exit("\n");
5016 }
5017
5018 if(!strcmp(init_tree,"user") ||
5019 !strcmp(init_tree,"User"))
5020 {
5021 char *starting_tree = NULL;
5022 starting_tree = XML_Get_Attribute_Value(instance,"file.name");
5023
5024 if(!Filexists(starting_tree))
5025 {
5026 PhyML_Fprintf(stderr,"\n. The tree file '%s' could not be found.",starting_tree);
5027 Exit("\n");
5028 }
5029 else
5030 {
5031 strcpy(io->in_tree_file,starting_tree);
5032 io->in_tree = 2;
5033 io->fp_in_tree = Openfile(io->in_tree_file,0);
5034 }
5035 }
5036 else if(!strcmp(init_tree,"random") ||
5037 !strcmp(init_tree,"Random"))
5038 {
5039 char *n_rand_starts = NULL;
5040
5041 io->mod->s_opt->random_input_tree = YES;
5042
5043 n_rand_starts = XML_Get_Attribute_Value(instance,"n.rand.starts");
5044
5045 if(n_rand_starts)
5046 {
5047 mod->s_opt->n_rand_starts = atoi(n_rand_starts);
5048 if(mod->s_opt->n_rand_starts < 1) Exit("\n. Number of random starting trees must be > 0.\n\n");
5049 }
5050
5051 strcpy(io->out_trees_file,io->in_align_file);
5052 strcat(io->out_trees_file,"_phyml_rand_trees");
5053 if(io->append_run_ID) { strcat(io->out_trees_file,"_"); strcat(io->out_trees_file,io->run_id_string); }
5054 strcat(io->out_trees_file,".txt");
5055 io->fp_out_trees = Openfile(io->out_trees_file,1);
5056 }
5057 else if(!strcmp(init_tree,"parsimony") ||
5058 !strcmp(init_tree,"Parsimony"))
5059 {
5060 io->in_tree = 1;
5061 }
5062
5063 // Estimate tree topology
5064 char *optimise = NULL;
5065
5066 optimise = XML_Get_Attribute_Value(instance,"optimise.tree");
5067
5068 if(optimise)
5069 {
5070 int select;
5071
5072 select = XML_Validate_Attr_Int(optimise,6,
5073 "true","yes","y",
5074 "false","no","n");
5075
5076 if(select > 2) io->mod->s_opt->opt_topo = NO;
5077 else
5078 {
5079 char *search;
5080 int select;
5081
5082 search = XML_Get_Attribute_Value(instance,"search");
5083
5084 if(search == NULL)
5085 {
5086 io->mod->s_opt->topo_search = SPR_MOVE;
5087 io->mod->s_opt->opt_topo = YES;
5088 }
5089 else
5090 {
5091 select = XML_Validate_Attr_Int(search,4,"spr","nni","best","none");
5092
5093 switch(select)
5094 {
5095 case 0:
5096 {
5097 io->mod->s_opt->topo_search = SPR_MOVE;
5098 io->mod->s_opt->opt_topo = YES;
5099 break;
5100 }
5101 case 1:
5102 {
5103 io->mod->s_opt->topo_search = NNI_MOVE;
5104 io->mod->s_opt->opt_topo = YES;
5105 break;
5106 }
5107 case 2:
5108 {
5109 io->mod->s_opt->topo_search = BEST_OF_NNI_AND_SPR;
5110 io->mod->s_opt->opt_topo = YES;
5111 break;
5112 }
5113 case 3:
5114 {
5115 io->mod->s_opt->opt_topo = NO;
5116 break;
5117 }
5118 default:
5119 {
5120 PhyML_Fprintf(stderr,"\n. Topology search option '%s' is not valid.",search);
5121 Exit("\n");
5122 break;
5123 }
5124 }
5125 }
5126 }
5127 }
5128 }
5129
5130
5131 //////////////////////////////////////////////////////////////
5132 //////////////////////////////////////////////////////////////
5133
Make_RAS_From_XML_Node(xml_node * parent,t_mod * mod)5134 void Make_RAS_From_XML_Node(xml_node *parent, t_mod *mod)
5135 {
5136 xml_node *w;
5137 char *family;
5138 int select;
5139
5140 family = NULL;
5141 mod->ras->n_catg = 0;
5142
5143 XML_Check_Siterates_Node(parent);
5144
5145 w = XML_Search_Node_Name("weights",YES,parent);
5146 if(w)
5147 {
5148 family = XML_Get_Attribute_Value(w,"family");
5149 if(family == NULL) select = -1;
5150 else select = XML_Validate_Attr_Int(family,3,"gamma","gamma+inv","freerates");
5151 switch(select)
5152 {
5153 case 0:// Gamma model
5154 {
5155 char *alpha,*alpha_opt;
5156
5157 mod->s_opt->opt_pinvar = NO;
5158 mod->ras->invar = NO;
5159
5160 alpha = XML_Get_Attribute_Value(w,"alpha");
5161
5162 if(alpha)
5163 {
5164 if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") ||
5165 !strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5166 {
5167 mod->s_opt->opt_alpha = YES;
5168 }
5169 else
5170 {
5171 mod->s_opt->opt_alpha = NO;
5172 mod->ras->alpha->v = String_To_Dbl(alpha);
5173 }
5174 }
5175
5176 alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5177
5178 if(alpha_opt)
5179 {
5180 if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5181 {
5182 mod->s_opt->opt_alpha = YES;
5183 }
5184 else
5185 {
5186 mod->s_opt->opt_alpha = NO;
5187 }
5188 }
5189
5190 mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
5191
5192 Make_RAS_Complete(mod->ras);
5193
5194 break;
5195 }
5196 case 1: // Gamma+Inv model
5197 {
5198 char *alpha,*pinv,*alpha_opt,*pinv_opt;
5199
5200 mod->ras->invar = YES;
5201 mod->s_opt->opt_pinvar = YES;
5202
5203 alpha = XML_Get_Attribute_Value(w,"alpha");
5204
5205 if(alpha)
5206 {
5207 if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") ||
5208 !strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5209 {
5210 mod->s_opt->opt_alpha = YES;
5211 }
5212 else
5213 {
5214 mod->s_opt->opt_alpha = NO;
5215 mod->ras->alpha->v = String_To_Dbl(alpha);;
5216 }
5217 }
5218
5219 alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5220
5221 if(alpha_opt)
5222 {
5223 if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5224 {
5225 mod->s_opt->opt_alpha = YES;
5226 }
5227 else
5228 {
5229 mod->s_opt->opt_alpha = NO;
5230 }
5231 }
5232
5233
5234 pinv = XML_Get_Attribute_Value(w,"pinv");
5235
5236 if(pinv)
5237 {
5238 if(!strcmp(pinv,"estimate") || !strcmp(pinv,"estimated") ||
5239 !strcmp(pinv,"optimise") || !strcmp(pinv,"optimised"))
5240 {
5241 mod->s_opt->opt_pinvar = YES;
5242 }
5243 else
5244 {
5245 mod->s_opt->opt_pinvar = NO;
5246 mod->ras->pinvar->v = String_To_Dbl(pinv);;
5247 }
5248 }
5249
5250 pinv_opt = XML_Get_Attribute_Value(w,"optimise.pinv");
5251
5252 if(pinv_opt)
5253 {
5254 if(!strcmp(pinv_opt,"yes") || !strcmp(pinv_opt,"true"))
5255 {
5256 mod->s_opt->opt_pinvar = YES;
5257 }
5258 else
5259 {
5260 mod->s_opt->opt_pinvar = NO;
5261 }
5262 }
5263
5264 mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
5265 break;
5266 }
5267 case 2: // FreeRate model
5268 {
5269 char *opt_freerates;
5270
5271 mod->ras->free_mixt_rates = YES;
5272
5273 mod->s_opt->opt_free_mixt_rates = YES;
5274
5275 opt_freerates = XML_Get_Attribute_Value(w,"optimise.freerates");
5276
5277 if(opt_freerates)
5278 {
5279 if(!strcmp(opt_freerates,"yes") || !strcmp(opt_freerates,"true"))
5280 {
5281 mod->s_opt->opt_free_mixt_rates = YES;
5282 }
5283 else
5284 {
5285 mod->s_opt->opt_free_mixt_rates = NO;
5286 }
5287 }
5288
5289 mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
5290 break;
5291 }
5292 default:
5293 {
5294 if(family != NULL) PhyML_Printf("\n. family: %s",family);
5295 else
5296 {
5297 PhyML_Printf("\n. Please specify a model family (\"gamma\", \"gamma+inv\" or \"freerate\")");
5298 PhyML_Printf("\n. for every 'weights' element in every 'siterate' element. Note that " );
5299 PhyML_Printf("\n. a \"gamma\" or a \"freerate\" model with a single rate class amounts");
5300 PhyML_Printf("\n. to no variation of rates across sites, if this is the model you'd");
5301 PhyML_Printf("\n. like to implement...");
5302 }
5303
5304 PhyML_Printf("\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
5305 Exit("\n");
5306 }
5307 }
5308 }
5309
5310 int nc = XML_Get_Number_Of_Classes_Siterates(parent);
5311
5312 if(w && nc != mod->ras->n_catg)
5313 {
5314 PhyML_Printf("\n. <siterates> component '%s'. The number of classes ",parent->id);
5315 PhyML_Printf("\n. specified in the <weight> component should be ");
5316 PhyML_Printf("\n. the same as that of <instance> nodes. Please fix");
5317 PhyML_Printf("\n. your XML file accordingly.");
5318 Exit("\n");
5319 }
5320
5321 if(!w) mod->ras->n_catg = nc;
5322
5323 Make_RAS_Complete(mod->ras);
5324
5325 }
5326
5327 /*////////////////////////////////////////////////////////////
5328 ////////////////////////////////////////////////////////////*/
5329
Generic_Exit(const char * file,int line,const char * function)5330 void Generic_Exit(const char *file, int line, const char *function)
5331 {
5332 PhyML_Fprintf(stderr,"\n. Err. in file '%s' (line %d)",file,line);
5333 if(function != NULL) PhyML_Printf(", function '%s'",function);
5334
5335 # if defined(PHYTIME)
5336 PhyML_Fprintf(stderr,"\n. PhyTime finished prematurely.");
5337 # elif defined(PHYREX)
5338 PhyML_Fprintf(stderr,"\n. PhyREX finished prematurely.");
5339 # elif defined(PHYML)
5340 PhyML_Fprintf(stderr,"\n. PhyML finished prematurely.");
5341 #else
5342 PhyML_Fprintf(stderr,"\n. The execution finished prematurely.");
5343 #endif
5344
5345 Exit("\n");
5346 }
5347
5348 /*////////////////////////////////////////////////////////////
5349 ////////////////////////////////////////////////////////////*/
5350
JSON_Write_Object(json_o * obj,FILE * where)5351 void JSON_Write_Object(json_o *obj, FILE *where)
5352 {
5353 json_kv *kv;
5354
5355 assert(obj);
5356 assert(where);
5357
5358 kv = obj->kv;
5359 assert(kv);
5360
5361 PhyML_Fprintf(where,"{");
5362
5363 do
5364 {
5365 PhyML_Fprintf(where,"\"%s\":",kv->key);
5366 if(kv->value != NULL) PhyML_Fprintf(where,"\"%s\"",kv->value);
5367 else if(kv->array != NULL) JSON_Write_Array(kv->array,where);
5368 else if(kv->object != NULL) JSON_Write_Object(kv->object,where);
5369 kv = kv->next;
5370 if(kv) PhyML_Fprintf(where,",");
5371 }
5372 while(kv);
5373
5374 PhyML_Fprintf(where,"}");
5375 }
5376
5377 /*////////////////////////////////////////////////////////////
5378 ////////////////////////////////////////////////////////////*/
5379
JSON_Write_Array(json_a * array,FILE * where)5380 void JSON_Write_Array(json_a *array, FILE *where)
5381 {
5382 json_o *o;
5383
5384 assert(where);
5385 assert(array);
5386
5387 o = array->object;
5388
5389 assert(o);
5390
5391 PhyML_Fprintf(where,"[");
5392
5393 do
5394 {
5395 JSON_Write_Object(o,where);
5396 o = o->next;
5397 if(o) PhyML_Fprintf(where,",");
5398 }
5399 while(o);
5400
5401 PhyML_Fprintf(where,"]\n");
5402 }
5403
5404 /*////////////////////////////////////////////////////////////
5405 ////////////////////////////////////////////////////////////*/
5406
JSON_Write_All(json_a * array,FILE * where)5407 void JSON_Write_All(json_a *array, FILE *where)
5408 {
5409 JSON_Write_Array(array,where);
5410 }
5411
5412 /*////////////////////////////////////////////////////////////
5413 ////////////////////////////////////////////////////////////*/
5414
JSON_Tree_To_Object(t_tree * mixt_tree)5415 json_o *JSON_Tree_To_Object(t_tree *mixt_tree)
5416 {
5417 t_tree *tree;
5418 json_kv *state,*kv;
5419 json_o *ret,*o;
5420 json_a *a;
5421 int string_len;
5422 char *s;
5423 time_t t_end;
5424
5425 string_len = 20;
5426
5427 ret = (json_o *)mCalloc(1,sizeof(json_o));
5428 ret->next = NULL;
5429
5430 state = (json_kv *)mCalloc(1,sizeof(json_kv));
5431 ret->kv = state;
5432
5433 state->key = (char *)mCalloc(string_len,sizeof(char));
5434 strcpy(state->key,"state");
5435
5436 state->value = NULL; state->array = NULL; state->object = NULL; state->next = NULL;
5437
5438 state->array = (json_a *)mCalloc(1,sizeof(json_a));
5439 a = state->array;
5440
5441 a->object = (json_o *)mCalloc(1,sizeof(json_o));
5442 o = a->object;
5443
5444 // State num
5445 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5446 kv = o->kv;
5447 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5448 kv->key = (char *)mCalloc(string_len,sizeof(char));
5449 strcpy(kv->key,"type");
5450 kv->value = (char *)mCalloc(string_len,sizeof(char));
5451 strcpy(kv->value,"t_num");
5452
5453 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5454 kv = kv->next;
5455
5456 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5457 kv->key = (char *)mCalloc(string_len,sizeof(char));
5458 strcpy(kv->key,"id");
5459 kv->value = (char *)mCalloc(string_len,sizeof(char));
5460 strcpy(kv->value,"state_num");
5461
5462 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5463 kv = kv->next;
5464
5465 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5466 kv->key = (char *)mCalloc(string_len,sizeof(char));
5467 strcpy(kv->key,"value");
5468 kv->value = (char *)mCalloc(string_len,sizeof(char));
5469 sprintf(kv->value,"%d",mixt_tree->json_num);
5470 mixt_tree->json_num++;
5471
5472 kv->next = NULL;
5473
5474
5475
5476
5477 // Time
5478 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5479 o = o->next;
5480
5481 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5482 kv = o->kv;
5483
5484 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5485 kv->key = (char *)mCalloc(string_len,sizeof(char));
5486 strcpy(kv->key,"type");
5487 kv->value = (char *)mCalloc(string_len,sizeof(char));
5488 strcpy(kv->value,"t_time");
5489
5490 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5491 kv = kv->next;
5492 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5493 kv->key = (char *)mCalloc(string_len,sizeof(char));
5494 strcpy(kv->key,"id");
5495 kv->value = (char *)mCalloc(string_len,sizeof(char));
5496 strcpy(kv->value,"time");
5497
5498
5499 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5500 kv = kv->next;
5501 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5502 kv->key = (char *)mCalloc(string_len,sizeof(char));
5503 strcpy(kv->key,"value");
5504 kv->value = (char *)mCalloc(string_len,sizeof(char));
5505 time(&t_end);
5506 sprintf(kv->value,"%d",(int)(t_end-mixt_tree->t_beg));
5507
5508 kv->next = NULL;
5509
5510
5511 // Tree
5512 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5513 o = o->next;
5514
5515 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5516 kv = o->kv;
5517
5518 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5519 kv->key = (char *)mCalloc(string_len,sizeof(char));
5520 strcpy(kv->key,"type");
5521 kv->value = (char *)mCalloc(string_len,sizeof(char));
5522 strcpy(kv->value,"t_tree");
5523
5524 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5525 kv = kv->next;
5526 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5527 kv->key = (char *)mCalloc(string_len,sizeof(char));
5528 strcpy(kv->key,"id");
5529 kv->value = (char *)mCalloc(string_len,sizeof(char));
5530 strcpy(kv->value,"tree");
5531
5532
5533 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5534 kv = kv->next;
5535 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5536 kv->key = (char *)mCalloc(string_len,sizeof(char));
5537 strcpy(kv->key,"value");
5538 s = Write_Tree(mixt_tree);
5539 kv->value = (char *)mCalloc((int)strlen(s)+1,sizeof(char));
5540 sprintf(kv->value,"%s",s);
5541 Free(s);
5542
5543 kv->next = NULL;
5544
5545
5546
5547 // Likelihood
5548 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5549 o = o->next;
5550
5551 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5552 kv = o->kv;
5553
5554 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5555 kv->key = (char *)mCalloc(string_len,sizeof(char));
5556 strcpy(kv->key,"type");
5557 kv->value = (char *)mCalloc(string_len,sizeof(char));
5558 strcpy(kv->value,"t_param");
5559
5560 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5561 kv = kv->next;
5562 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5563 kv->key = (char *)mCalloc(string_len,sizeof(char));
5564 strcpy(kv->key,"id");
5565 kv->value = (char *)mCalloc(string_len,sizeof(char));
5566 strcpy(kv->value,"likelihood");
5567
5568
5569 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5570 kv = kv->next;
5571 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5572 kv->key = (char *)mCalloc(string_len,sizeof(char));
5573 strcpy(kv->key,"value");
5574 kv->value = (char *)mCalloc(string_len,sizeof(char));
5575 sprintf(kv->value,"%G",mixt_tree->c_lnL);
5576
5577 kv->next = NULL;
5578
5579
5580
5581 // TsTv
5582 {
5583 int n_tstv,i;
5584 scalar_dbl **tstv;
5585
5586 n_tstv = 0;
5587 tstv = NULL;
5588 tree = mixt_tree;
5589
5590 do
5591 {
5592 if(tree->is_mixt_tree == YES) tree = tree->next;
5593
5594 for(i=0;i<n_tstv;++i) if(tree->mod->kappa == tstv[i]) break;
5595
5596 if(i == n_tstv)
5597 {
5598 if(!tstv) tstv = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
5599 else tstv = (scalar_dbl **)mRealloc(tstv,n_tstv+1,sizeof(scalar_dbl *));
5600 tstv[n_tstv] = tree->mod->kappa;
5601 n_tstv++;
5602
5603 if(tree->mod->s_opt->opt_kappa == YES)
5604 {
5605 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5606 o = o->next;
5607
5608 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5609 kv = o->kv;
5610
5611 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5612 kv->key = (char *)mCalloc(string_len,sizeof(char));
5613 strcpy(kv->key,"type");
5614 kv->value = (char *)mCalloc(string_len,sizeof(char));
5615 strcpy(kv->value,"t_param");
5616
5617 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5618 kv = kv->next;
5619 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5620 kv->key = (char *)mCalloc(string_len,sizeof(char));
5621 strcpy(kv->key,"id");
5622 kv->value = (char *)mCalloc(string_len,sizeof(char));
5623 strcpy(kv->value,"tstv");
5624 sprintf(kv->value+strlen(kv->value),"%d",n_tstv);
5625
5626 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5627 kv = kv->next;
5628 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5629 kv->key = (char *)mCalloc(string_len,sizeof(char));
5630 strcpy(kv->key,"value");
5631 kv->value = (char *)mCalloc(string_len,sizeof(char));
5632 sprintf(kv->value,"%G",tree->mod->kappa->v);
5633
5634 kv->next = NULL;
5635
5636 }
5637 }
5638 tree = tree->next;
5639 }
5640 while(tree);
5641
5642 if(tstv) Free(tstv);
5643 }
5644
5645
5646
5647 // Alpha
5648 {
5649 int n_alpha,i;
5650 scalar_dbl **alpha;
5651
5652 n_alpha = 0;
5653 alpha = NULL;
5654 tree = mixt_tree;
5655
5656 do
5657 {
5658
5659 for(i=0;i<n_alpha;i++) if(tree->mod->ras->alpha == alpha[i]) break;
5660
5661 if(i == n_alpha)
5662 {
5663 if(!alpha) alpha = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
5664 else alpha = (scalar_dbl **)mRealloc(alpha,n_alpha+1,sizeof(scalar_dbl *));
5665 alpha[n_alpha] = tree->mod->ras->alpha;
5666 n_alpha++;
5667
5668 if(tree->mod->s_opt->opt_alpha == YES)
5669 {
5670 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5671 o = o->next;
5672
5673 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5674 kv = o->kv;
5675
5676 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5677 kv->key = (char *)mCalloc(string_len,sizeof(char));
5678 strcpy(kv->key,"type");
5679 kv->value = (char *)mCalloc(string_len,sizeof(char));
5680 strcpy(kv->value,"t_param");
5681
5682 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5683 kv = kv->next;
5684 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5685 kv->key = (char *)mCalloc(string_len,sizeof(char));
5686 strcpy(kv->key,"id");
5687 kv->value = (char *)mCalloc(string_len,sizeof(char));
5688 strcpy(kv->value,"alpha");
5689 sprintf(kv->value+strlen(kv->value),"%d",n_alpha);
5690
5691 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5692 kv = kv->next;
5693 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5694 kv->key = (char *)mCalloc(string_len,sizeof(char));
5695 strcpy(kv->key,"value");
5696 kv->value = (char *)mCalloc(string_len,sizeof(char));
5697 sprintf(kv->value,"%G",tree->mod->ras->alpha->v);
5698
5699 kv->next = NULL;
5700
5701 }
5702 }
5703 tree = tree->next_mixt;
5704 }
5705 while(tree);
5706
5707 if(alpha) Free(alpha);
5708 }
5709
5710
5711
5712 // Tree size
5713 {
5714 tree = mixt_tree;
5715 int tree_num = 1;
5716 do
5717 {
5718 o->next = (json_o *)mCalloc(1,sizeof(json_o));
5719 o = o->next;
5720
5721 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5722 kv = o->kv;
5723
5724 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5725 kv->key = (char *)mCalloc(string_len,sizeof(char));
5726 strcpy(kv->key,"type");
5727 kv->value = (char *)mCalloc(string_len,sizeof(char));
5728 strcpy(kv->value,"t_param");
5729
5730 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5731 kv = kv->next;
5732 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5733 kv->key = (char *)mCalloc(string_len,sizeof(char));
5734 strcpy(kv->key,"id");
5735 kv->value = (char *)mCalloc(string_len,sizeof(char));
5736 strcpy(kv->value,"tree_size");
5737 sprintf(kv->value+strlen(kv->value),"%d",tree_num);
5738
5739 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5740 kv = kv->next;
5741 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5742 kv->key = (char *)mCalloc(string_len,sizeof(char));
5743 strcpy(kv->key,"value");
5744 kv->value = (char *)mCalloc(string_len,sizeof(char));
5745 sprintf(kv->value,"%G",Get_Tree_Size(tree));
5746
5747 kv->next = NULL;
5748
5749 tree = tree->next_mixt;
5750 }
5751 while(tree);
5752
5753 }
5754
5755 return(ret);
5756 }
5757
5758 /*////////////////////////////////////////////////////////////
5759 ////////////////////////////////////////////////////////////*/
5760
5761
JSON_Tree_To_Object_Light(t_tree * mixt_tree)5762 json_o *JSON_Tree_To_Object_Light(t_tree *mixt_tree)
5763 {
5764 t_tree *tree;
5765 json_kv *state,*kv;
5766 json_o *ret,*o;
5767 int string_len;
5768 char *s;
5769 time_t t_end;
5770
5771 string_len = 20;
5772
5773 ret = (json_o *)mCalloc(1,sizeof(json_o));
5774 ret->next = NULL;
5775
5776 state = (json_kv *)mCalloc(1,sizeof(json_kv));
5777 ret->kv = state;
5778
5779 state->key = (char *)mCalloc(string_len,sizeof(char));
5780 strcpy(state->key,"state");
5781
5782 state->value = NULL; state->array = NULL; state->object = NULL; state->next = NULL;
5783
5784 state->object = (json_o *)mCalloc(1,sizeof(json_o));
5785 o = state->object;
5786
5787
5788 o->kv = (json_kv *)mCalloc(1,sizeof(json_kv));
5789 kv = o->kv;
5790
5791
5792 // State num
5793 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5794 kv->key = (char *)mCalloc(string_len,sizeof(char));
5795 strcpy(kv->key,"state_num");
5796 kv->value = (char *)mCalloc(string_len,sizeof(char));
5797 sprintf(kv->value,"%d",mixt_tree->json_num);
5798
5799
5800
5801 // Time
5802 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5803 kv = kv->next;
5804 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5805 kv->key = (char *)mCalloc(string_len,sizeof(char));
5806 strcpy(kv->key,"time");
5807 kv->value = (char *)mCalloc(string_len,sizeof(char));
5808 time(&t_end);
5809 sprintf(kv->value,"%d",(int)(t_end-mixt_tree->t_beg));
5810
5811
5812
5813
5814 // Tree
5815 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5816 kv = kv->next;
5817
5818 kv->value = NULL; kv->array = NULL; kv->object = NULL; kv->next = NULL;
5819 kv->key = (char *)mCalloc(string_len,sizeof(char));
5820 strcpy(kv->key,"tree");
5821 s = Write_Tree(mixt_tree);
5822 kv->value = (char *)mCalloc((int)strlen(s)+1,sizeof(char));
5823 sprintf(kv->value,"%s",s);
5824 Free(s);
5825
5826
5827
5828
5829
5830 // Likelihood
5831 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5832 kv = kv->next;
5833 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5834 kv->key = (char *)mCalloc(string_len,sizeof(char));
5835 strcpy(kv->key,"likelihood");
5836 kv->value = (char *)mCalloc(string_len,sizeof(char));
5837 sprintf(kv->value,"%G",mixt_tree->c_lnL);
5838
5839
5840
5841
5842 // TsTv
5843 {
5844 int n_tstv,i;
5845 scalar_dbl **tstv;
5846
5847 n_tstv = 0;
5848 tstv = NULL;
5849 tree = mixt_tree;
5850
5851 do
5852 {
5853 if(tree->is_mixt_tree == YES) tree = tree->next;
5854
5855 for(i=0;i<n_tstv;++i) if(tree->mod->kappa == tstv[i]) break;
5856
5857 if(i == n_tstv)
5858 {
5859 if(!tstv) tstv = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
5860 else tstv = (scalar_dbl **)mRealloc(tstv,n_tstv+1,sizeof(scalar_dbl *));
5861 tstv[n_tstv] = tree->mod->kappa;
5862 n_tstv++;
5863
5864 if(tree->mod->s_opt->opt_kappa == YES)
5865 {
5866
5867 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5868 kv = kv->next;
5869 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5870 kv->key = (char *)mCalloc(string_len,sizeof(char));
5871 strcpy(kv->key,"tstv");
5872 sprintf(kv->key+strlen(kv->key),"%d",n_tstv);
5873 kv->value = (char *)mCalloc(string_len,sizeof(char));
5874 sprintf(kv->value,"%G",tree->mod->kappa->v);
5875 }
5876 }
5877 tree = tree->next;
5878 }
5879 while(tree);
5880
5881 if(tstv) Free(tstv);
5882 }
5883
5884
5885
5886 // Alpha
5887 {
5888 int n_alpha,i;
5889 scalar_dbl **alpha;
5890
5891 n_alpha = 0;
5892 alpha = NULL;
5893 tree = mixt_tree;
5894
5895 do
5896 {
5897
5898 for(i=0;i<n_alpha;i++) if(tree->mod->ras->alpha == alpha[i]) break;
5899
5900 if(i == n_alpha)
5901 {
5902 if(!alpha) alpha = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
5903 else alpha = (scalar_dbl **)mRealloc(alpha,n_alpha+1,sizeof(scalar_dbl *));
5904 alpha[n_alpha] = tree->mod->ras->alpha;
5905 n_alpha++;
5906
5907 if(tree->mod->s_opt->opt_alpha == YES)
5908 {
5909
5910 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5911 kv = kv->next;
5912 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5913 kv->key = (char *)mCalloc(string_len,sizeof(char));
5914 strcpy(kv->key,"alpha");
5915 sprintf(kv->key+strlen(kv->key),"%d",n_alpha);
5916 kv->value = (char *)mCalloc(string_len,sizeof(char));
5917 sprintf(kv->value,"%G",tree->mod->ras->alpha->v);
5918 }
5919 }
5920 tree = tree->next_mixt;
5921 }
5922 while(tree);
5923
5924 if(alpha) Free(alpha);
5925 }
5926
5927
5928
5929 // Tree size
5930 {
5931 tree = mixt_tree;
5932 int tree_num = 1;
5933 do
5934 {
5935
5936 kv->next = (json_kv *)mCalloc(1,sizeof(json_kv));
5937 kv = kv->next;
5938
5939 kv->value = NULL; kv->array = NULL; kv->object = NULL; o->next = NULL;
5940 kv->key = (char *)mCalloc(string_len,sizeof(char));
5941 strcpy(kv->key,"tree_size");
5942 sprintf(kv->key+strlen(kv->key),"%d",tree_num);
5943 kv->value = (char *)mCalloc(string_len,sizeof(char));
5944 sprintf(kv->value,"%G",Get_Tree_Size(tree));
5945
5946 tree = tree->next_mixt;
5947 }
5948 while(tree);
5949
5950 }
5951
5952 kv->next = NULL;
5953
5954 return(ret);
5955 }
5956
5957 /*////////////////////////////////////////////////////////////
5958 ////////////////////////////////////////////////////////////*/
5959
JSON_Tree_Io(t_tree * tree,FILE * where)5960 void JSON_Tree_Io(t_tree *tree, FILE *where)
5961 {
5962 // Append
5963 json_o *o;
5964 fpos_t pos;
5965 char c;
5966
5967 fgetpos(where,&pos);
5968
5969 rewind(where);
5970 c = fgetc(where);
5971
5972 if(c != '[')
5973 {
5974 PhyML_Fprintf(where,"[");
5975 }
5976 else
5977 {
5978 fsetpos(where,&pos);
5979 fseek(where,-1,SEEK_CUR);
5980 PhyML_Fprintf(where,",");
5981 }
5982
5983 PhyML_Fprintf(where,"\n");
5984 /* o = JSON_Tree_To_Object(tree); */
5985 o = JSON_Tree_To_Object_Light(tree);
5986 JSON_Write_Object(o,where);
5987 JSON_Free_Object(o);
5988 PhyML_Fprintf(where,"]");
5989 fflush(where);
5990
5991
5992 // Print out latest tree + info
5993 /* json_o *o; */
5994
5995 /* rewind(where); */
5996 /* /\* PhyML_Fprintf(where,"["); *\/ */
5997 /* o = JSON_Tree_To_Object(tree); */
5998 /* JSON_Write_Object(o,where); */
5999 /* JSON_Free_Object(o); */
6000 /* /\* PhyML_Fprintf(where,"]"); *\/ */
6001 /* fflush(where); */
6002
6003
6004 }
6005
6006 /*////////////////////////////////////////////////////////////
6007 ////////////////////////////////////////////////////////////*/
6008
Print_Lk_Given_Edge_Recurr(t_node * a,t_node * d,t_edge * b,t_tree * tree)6009 void Print_Lk_Given_Edge_Recurr(t_node *a, t_node *d, t_edge *b, t_tree *tree)
6010 {
6011 PhyML_Printf("\n___ Edge %3d (left=%3d rght=%3d) lnL=%f",
6012 b->num,
6013 b->left->num,
6014 b->rght->num,
6015 Lk(b,tree));
6016
6017 if(d->tax) return;
6018 else
6019 {
6020 int i;
6021 for(i=0;i<3;i++)
6022 if(d->v[i] != a)
6023 Print_Lk_Given_Edge_Recurr(d,d->v[i],d->b[i],tree);
6024 }
6025 }
6026
6027 /*////////////////////////////////////////////////////////////
6028 ////////////////////////////////////////////////////////////*/
6029
Collect_Edge_Support_Values(t_tree * tree)6030 void Collect_Edge_Support_Values(t_tree *tree)
6031 {
6032 int i;
6033
6034 for(i=0;i<2*tree->n_otu-3;++i)
6035 {
6036 if(tree->io->do_boot == YES)
6037 {
6038 tree->a_edges[i]->support_val = tree->a_edges[i]->bip_score;
6039 }
6040 else if(tree->io->do_tbe == YES)
6041 {
6042 int pminus1=MIN(tree->a_edges[i]->left->bip_size[tree->a_edges[i]->l_r], tree->a_edges[i]->rght->bip_size[tree->a_edges[i]->r_l])-1;
6043 tree->a_edges[i]->support_val = 1-((tree->a_edges[i]->tdist_score)/(tree->io->n_boot_replicates*1.0))/pminus1;
6044 }
6045 else if(tree->io->do_alrt == YES)
6046 {
6047 tree->a_edges[i]->support_val = tree->a_edges[i]->ratio_test;
6048 }
6049 }
6050 }
6051
6052 /*////////////////////////////////////////////////////////////
6053 ////////////////////////////////////////////////////////////*/
6054
6055
6056 #if defined (PHYREX)
PHYREX_Output_Tree_Structure(FILE * fp,t_tree * tree)6057 void PHYREX_Output_Tree_Structure(FILE *fp, t_tree *tree)
6058 {
6059 char *s;
6060 s = PHYREX_Print_Tree_Structure(tree);
6061 PhyML_Fprintf(fp,"%s",s);
6062 Free(s);
6063 }
6064 #endif
6065
6066 /*////////////////////////////////////////////////////////////
6067 ////////////////////////////////////////////////////////////*/
6068
6069 #if defined (PHYREX)
PHYREX_Print_Tree_Structure(t_tree * tree)6070 char *PHYREX_Print_Tree_Structure(t_tree *tree)
6071 {
6072 t_dsk *disk;
6073 char *s,*buff;
6074 FILE *fp;
6075 fpos_t pos;
6076
6077 fp = tmpfile();
6078 assert(fp);
6079
6080 buff = (char *)mCalloc(T_MAX_LINE,sizeof(char));
6081
6082 disk = tree->young_disk;
6083 while(disk->prev != NULL) disk = disk->prev;
6084
6085 assert(disk->ldsk);
6086 assert(disk->ldsk->n_next > 0);
6087
6088 fgetpos(fp,&pos);
6089 PhyML_Fprintf(fp,"begin\n");
6090 fsetpos(fp,&pos);
6091 if(fgets(buff,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6092 s = (char *)mCalloc((int)strlen(buff)+1,sizeof(char));
6093 sprintf(s+strlen(s),"%s",buff);
6094
6095 fgetpos(fp,&pos);
6096 PhyML_Fprintf(fp,"%d %d\n",tree->n_otu,tree->mmod->n_dim);
6097 fsetpos(fp,&pos);
6098 if(fgets(buff,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6099 s = (char *)mRealloc(s,(int)(strlen(s)+strlen(buff)+1),sizeof(char));
6100 sprintf(s+strlen(s),"%s",buff);
6101
6102
6103 do
6104 {
6105 fgetpos(fp,&pos);
6106 PhyML_Fprintf(fp,"%s ",disk->id);
6107 PhyML_Fprintf(fp,"%f ",disk->time);
6108 PhyML_Fprintf(fp,"%f %f ",disk->centr->lonlat[0],disk->centr->lonlat[1]);
6109
6110 if(disk->ldsk != NULL)
6111 {
6112 s = (char *)mRealloc(s,(int)strlen(s)+9,sizeof(char));
6113 PhyML_Fprintf(fp,"*%s ",disk->ldsk->coord->id);
6114 PhyML_Fprintf(fp,"%f %f ",disk->ldsk->coord->lonlat[0],disk->ldsk->coord->lonlat[1]);
6115 }
6116 for(int i=0;i<disk->n_ldsk_a;++i)
6117 {
6118 PhyML_Fprintf(fp,"%s ",disk->ldsk_a[i]->coord->id);
6119 }
6120 PhyML_Fprintf(fp,"\n");
6121
6122 fsetpos(fp,&pos);
6123 if(fgets(buff,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6124 s = (char *)mRealloc(s,(int)(strlen(s)+strlen(buff)+1),sizeof(char));
6125 sprintf(s+strlen(s),"%s",buff);
6126
6127 disk = disk->next;
6128 }
6129 while(disk);
6130
6131 fgetpos(fp,&pos);
6132 PhyML_Fprintf(fp,"end");
6133 fsetpos(fp,&pos);
6134 if(fgets(buff,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6135 s = (char *)mRealloc(s,(int)(strlen(s)+strlen(buff)+1),sizeof(char));
6136 sprintf(s+strlen(s),"%s",buff);
6137
6138 fclose(fp);
6139
6140 Free(buff);
6141
6142 return(s);
6143 }
6144 #endif
6145
6146 /*////////////////////////////////////////////////////////////
6147 ////////////////////////////////////////////////////////////*/
6148
6149 #if defined (PHYREX)
PHYREX_Check_Point(FILE * fp,t_tree * tree)6150 void PHYREX_Check_Point(FILE *fp, t_tree *tree)
6151 {
6152 xml_node *n;
6153 char *s;
6154
6155 if(XML_Search_Node_Attribute_Value("add","true",NO,tree->xml_root) == NULL)
6156 {
6157 XML_Add_Attribute(tree->xml_root,"add","true");
6158 }
6159
6160 n = XML_Search_Node_Name("slfv",YES,tree->xml_root);
6161 if(n == NULL)
6162 {
6163 n = XML_Add_Node(tree->xml_root,"slfv");
6164 XML_Set_Node_Id(n,"SLFV1");
6165 }
6166
6167 s = PHYREX_Print_Tree_Structure(tree);
6168 XML_Set_Node_Value(n,s);
6169 XML_Update_XML_Struct_Given_Model_Params(tree);
6170 /* XML_Update_XML_Struct_Given_MCMC_Params(tree); */
6171 XML_Write_XML_Graph(fp,tree->xml_root);
6172 }
6173 #endif
6174
6175 /*////////////////////////////////////////////////////////////
6176 ////////////////////////////////////////////////////////////*/
6177
6178 #if defined (PHYREX)
PHYREX_Input_Tree_Structure(FILE * fp)6179 void PHYREX_Input_Tree_Structure(FILE *fp)
6180 {
6181 char *tk,*s,delim[4]="\n\r ";
6182 int n_otu,n_dim,idx;
6183 t_dsk *disk;
6184 t_ldsk *new_ldsk,*root_ldsk,*ldsk;
6185
6186 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
6187
6188 do
6189 {
6190 if(fgets(s,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6191 }
6192 while(strstr(s,"begin") == NULL);
6193
6194 if(fgets(s,T_MAX_LINE,fp) == NULL) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
6195
6196 tk = strtok(s,delim);
6197 n_otu = strtod(tk,NULL);
6198
6199 tk = strtok(NULL,delim);
6200 n_dim = strtod(tk,NULL);
6201
6202
6203 root_ldsk = NULL;
6204 new_ldsk = NULL;
6205 while(fgets(s,T_MAX_LINE,fp) != NULL)
6206 {
6207 tk = strtok(s,delim);
6208 if(strcmp(tk,"end") && tk != NULL)
6209 {
6210 disk = PHYREX_Make_Disk_Event(n_dim,n_otu);
6211 strcpy(disk->id,tk);
6212
6213 tk = strtok(NULL,delim);
6214 disk->time = strtold(tk,NULL);
6215
6216 tk = strtok(NULL,delim);
6217 disk->centr->lonlat[0] = strtold(tk,NULL);
6218 tk = strtok(NULL,delim);
6219 disk->centr->lonlat[1] = strtold(tk,NULL);
6220
6221
6222 tk = strtok(NULL,delim);
6223 if(tk[0] == '*') // disk has ldsk on it
6224 {
6225 if(root_ldsk == NULL)
6226 {
6227 disk->ldsk = PHYREX_Make_Lindisk_Node(n_dim);
6228 disk->ldsk->disk = disk;
6229 root_ldsk = disk->ldsk;
6230 }
6231 else
6232 {
6233 disk->ldsk = PHYREX_Find_Ldsk_From_Id(tk+1,root_ldsk);
6234 assert(disk->ldsk != NULL);
6235 }
6236
6237 strcpy(disk->ldsk->coord->id,tk+1);
6238
6239 tk = strtok(NULL,delim);
6240 disk->ldsk->coord->lonlat[0] = strtold(tk,NULL);
6241 tk = strtok(NULL,delim);
6242 disk->ldsk->coord->lonlat[1] = strtold(tk,NULL);
6243
6244
6245 do
6246 {
6247 tk = strtok(NULL,delim);
6248 if(tk == NULL) break;
6249 new_ldsk = PHYREX_Make_Lindisk_Node(n_dim);
6250 strcpy(new_ldsk->coord->id,tk);
6251 PHYREX_Make_Lindisk_Next(disk->ldsk);
6252 disk->ldsk->next[disk->ldsk->n_next-1] = new_ldsk;
6253 disk->ldsk_a[disk->ldsk->n_next-1] = new_ldsk;
6254 new_ldsk->prev = disk->ldsk;
6255 }
6256 while(tk);
6257 }
6258 else
6259 {
6260 idx = 0;
6261 do
6262 {
6263 tk = strtok(NULL,delim);
6264 if(tk == NULL) break;
6265 ldsk = PHYREX_Find_Ldsk_From_Id(tk,root_ldsk);
6266 assert(ldsk);
6267 disk->ldsk_a[idx] = ldsk;
6268 idx++;
6269 }
6270 while(tk);
6271
6272 }
6273 }
6274 }
6275
6276 Free(s);
6277 }
6278 #endif
6279
6280 /*////////////////////////////////////////////////////////////
6281 ////////////////////////////////////////////////////////////*/
6282
Output_Scalar_Dbl(scalar_dbl * t,char * sep,FILE * fp)6283 void Output_Scalar_Dbl(scalar_dbl *t, char *sep, FILE *fp)
6284 {
6285 scalar_dbl *l;
6286 l = t;
6287 do
6288 {
6289 PhyML_Fprintf(fp,"%g%s",l->v,sep);
6290 l = l->next;
6291 }
6292 while(l);
6293 }
6294
6295 /*////////////////////////////////////////////////////////////
6296 ////////////////////////////////////////////////////////////*/
6297
6298 // s should look like that: "xxx={yyy},xxxx={yy},..."
Read_Labels(char * s)6299 t_label *Read_Labels(char *s)
6300 {
6301 t_label *lab,*init_lab;
6302 char *s_cpy,*s_big,*s_small;
6303 char *label;
6304 char *key,*val;
6305
6306 if(s[0] != '[' || s[(int)strlen(s)-1] != ']')
6307 {
6308 PhyML_Fprintf(stderr,"\n. Label is in wrong format. A proper label should");
6309 PhyML_Fprintf(stderr,"\n. look as follows: \"[xxx={yyy},xxxx={yy},...]\"");
6310 assert(FALSE);
6311 }
6312
6313 lab = Make_Label();
6314 init_lab = lab;
6315
6316 s_cpy = (char *)mCalloc((int)strlen(s)-1,sizeof(char));
6317 strncpy(s_cpy,s+1,(strlen(s)-2)*sizeof(char));
6318 s_cpy[strlen(s)-2]='\0';
6319
6320
6321 label = strtok_r(s_cpy,",",&s_big);
6322
6323 while(label != NULL)
6324 {
6325
6326 key=strtok_r(label,"=",&s_small);
6327 val=strtok_r(NULL,"=",&s_small);
6328
6329 Free(lab->key);
6330 lab->key = (char *)mCalloc(strlen(key)+1,sizeof(char));
6331 strcpy(lab->key,key);
6332
6333 Free(lab->val);
6334 lab->val = (char *)mCalloc(strlen(val)+1,sizeof(char));
6335 strcpy(lab->val,val);
6336
6337 label = strtok_r(NULL,",",&s_big);
6338
6339 if(label != NULL)
6340 {
6341 lab->sep=',';
6342 lab->next = Make_Label();
6343 lab = lab->next;
6344 }
6345 }
6346
6347 lab = init_lab;
6348
6349 return(lab);
6350 }
6351
6352 /*////////////////////////////////////////////////////////////
6353 ////////////////////////////////////////////////////////////*/
6354
Print_Labels(FILE * fp_where,char * s_where,t_label * label)6355 void Print_Labels(FILE *fp_where, char *s_where, t_label *label)
6356 {
6357 if(label == NULL) return;
6358 else
6359 {
6360 t_label *lab;
6361
6362 lab = label;
6363
6364 if(fp_where != NULL)
6365 {
6366 PhyML_Fprintf(fp_where,"[");
6367 }
6368 else
6369 {
6370 sprintf(s_where+(int)strlen(s_where),"[");
6371 }
6372
6373
6374 while(lab)
6375 {
6376 if(fp_where != NULL)
6377 {
6378 PhyML_Fprintf(fp_where,"%s=%s",lab->key,lab->val);
6379 if(lab->next != NULL) PhyML_Fprintf(fp_where,",");
6380 }
6381 else
6382 {
6383 sprintf(s_where+(int)strlen(s_where),"%s=%s",lab->key,lab->val);
6384 if(lab->next != NULL) sprintf(s_where+(int)strlen(s_where),",");
6385 }
6386
6387 lab = lab->next;
6388 }
6389
6390 if(fp_where != NULL)
6391 {
6392 PhyML_Fprintf(fp_where,"]");
6393 }
6394 else
6395 {
6396 sprintf(s_where+(int)strlen(s_where),"]");
6397 }
6398 }
6399 }
6400
6401
6402 /*////////////////////////////////////////////////////////////
6403 ////////////////////////////////////////////////////////////*/
6404