3
PhyML: a program that computes maximum likelihood phylogenies from
4
DNA or AA homologous sequences.
6
Copyright (C) Stephane Guindon. Oct 2003 onward.
8
All parts of the source except where indicated are distributed under
9
the GNU public licence. See http://www.opensource.org for details.
16
/* Tree parser function. We need to pass a pointer to the string of characters
17
since this string might be freed and then re-allocated by that function (i.e.,
18
its address in memory might change)
20
t_tree *Read_Tree(char **s_tree)
23
int i,n_ext,n_int,n_otu;
31
For(i,(int)strlen((*s_tree))) if((*s_tree)[i] == ',') n_otu++;
34
tree = Make_Tree_From_Scratch(n_otu,NULL);
35
subs = Sub_Trees((*s_tree),°ree);
36
Clean_Multifurcation(subs,degree,3);
40
/* Unroot_Tree(subs); */
42
/* root_node = tree->a_nodes[n_otu]; */
43
root_node = tree->a_nodes[2*n_otu-2];
44
root_node->num = 2*n_otu-2;
45
tree->n_root = root_node;
50
root_node = tree->a_nodes[n_otu];
51
root_node->num = n_otu;
55
if(degree > 3) /* Multifurcation at the root. Need to re-assemble the subtrees
56
since Clean_Multifurcation added sets of prevhesis and
57
the corresponding NULL edges */
62
For(i,degree) len += (strlen(subs[i])+1);
65
(*s_tree) = (char *)mCalloc(len,sizeof(char));
67
(*s_tree)[0] = '('; (*s_tree)[1] = '\0';
70
strcat((*s_tree),subs[i]);
71
strcat((*s_tree),",\0");
74
sprintf((*s_tree)+strlen((*s_tree))-1,"%s",");\0");
76
For(i,NODE_DEG_MAX) Free(subs[i]);
79
subs = Sub_Trees((*s_tree),°ree);
84
tree->has_branch_lengths = 0;
85
tree->num_curr_branch_available = 0;
86
For(i,degree) R_rtree((*s_tree),subs[i],root_node,tree,&n_int,&n_ext);
88
for(i=degree;i<NODE_DEG_MAX;i++) Free(subs[i]);
93
tree->e_root = tree->a_edges[tree->num_curr_branch_available];
95
tree->n_root->v[2] = tree->n_root->v[0];
96
tree->n_root->v[0] = NULL;
98
tree->n_root->l[2] = tree->n_root->l[0];
100
For(i,3) if(tree->n_root->v[2]->v[i] == tree->n_root) { tree->n_root->v[2]->v[i] = tree->n_root->v[1]; break; }
101
For(i,3) if(tree->n_root->v[1]->v[i] == tree->n_root) { tree->n_root->v[1]->v[i] = tree->n_root->v[2]; break; }
103
Connect_One_Edge_To_Two_Nodes(tree->n_root->v[2],
108
tree->e_root->l->v = tree->n_root->l[2] + tree->n_root->l[1];
109
if(tree->e_root->l->v > 0.0)
110
tree->n_root_pos = tree->n_root->l[2] / tree->e_root->l->v;
112
tree->n_root_pos = .5;
118
//////////////////////////////////////////////////////////////
119
//////////////////////////////////////////////////////////////
121
/* 'a' in t_node a stands for ancestor. 'd' stands for descendant */
122
void R_rtree(char *s_tree_a, char *s_tree_d, t_node *a, t_tree *tree, int *n_int, int *n_ext)
126
int n_otu = tree->n_otu;
128
if(strstr(s_tree_a," "))
130
PhyML_Printf("\n== [%s]",s_tree_a);
131
Warn_And_Exit("\n== Err: the tree must not contain a ' ' character\n");
134
if(s_tree_d[0] == '(')
142
if((*n_int + n_otu) == (2*n_otu-1))
144
PhyML_Printf("\n== The number of internal nodes in the tree exceeds the number of taxa minus one.");
145
PhyML_Printf("\n== There probably is a formating problem in the input tree.");
146
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
150
d = tree->a_nodes[n_otu+*n_int];
151
d->num = n_otu+*n_int;
153
b = tree->a_edges[tree->num_curr_branch_available];
155
Read_Branch_Label(s_tree_d,s_tree_a,tree->a_edges[tree->num_curr_branch_available]);
156
Read_Branch_Length(s_tree_d,s_tree_a,tree);
163
d->l[0]=tree->a_edges[tree->num_curr_branch_available]->l->v;
164
a->l[i]=tree->a_edges[tree->num_curr_branch_available]->l->v;
170
if(a != tree->n_root)
172
Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[tree->num_curr_branch_available],tree);
175
subs=Sub_Trees(s_tree_d,°ree);
179
Clean_Multifurcation(subs,degree,2);
183
s_tree_d = (char *)mCalloc(strlen(subs[0])+strlen(subs[1])+5,sizeof(char));
185
strcat(s_tree_d,"(");
186
strcat(s_tree_d,subs[0]);
187
strcat(s_tree_d,",");
188
strcat(s_tree_d,subs[1]);
189
strcat(s_tree_d,")");
192
strcat(s_tree_d,"#");
193
strcat(s_tree_d,b->labels[i]);
196
For(i,NODE_DEG_MAX) Free(subs[i]);
199
subs=Sub_Trees(s_tree_d,°ree);
202
R_rtree(s_tree_d,subs[0],d,tree,n_int,n_ext);
203
R_rtree(s_tree_d,subs[1],d,tree,n_int,n_ext);
205
for(i=2;i<NODE_DEG_MAX;i++) Free(subs[i]);
213
d = tree->a_nodes[*n_ext];
216
Read_Node_Name(d,s_tree_d,tree);
217
Read_Branch_Label(s_tree_d,s_tree_a,tree->a_edges[tree->num_curr_branch_available]);
218
Read_Branch_Length(s_tree_d,s_tree_a,tree);
225
d->l[0]=tree->a_edges[tree->num_curr_branch_available]->l->v;
226
a->l[i]=tree->a_edges[tree->num_curr_branch_available]->l->v;
232
if(a != tree->n_root)
234
Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[tree->num_curr_branch_available],tree);
245
//////////////////////////////////////////////////////////////
246
//////////////////////////////////////////////////////////////
249
void Read_Branch_Label(char *s_d, char *s_a, t_edge *b)
255
sub_tp = (char *)mCalloc(3+(int)strlen(s_d)+1,sizeof(char));
256
/* sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
262
p = strstr(s_a,sub_tp);
270
p = strstr(s_a,sub_tp);
276
if(!(b->n_labels%BLOCK_LABELS)) Make_New_Edge_Label(b);
280
while(p[posp] != '#') posp++;
286
b->labels[b->n_labels-1][posl] = p[posp];
291
b->labels[b->n_labels-1][posl] = '\0';
293
if(!(b->n_labels%BLOCK_LABELS)) Make_New_Edge_Label(b);
298
while((p[posp] != ':') &&
302
b->labels[b->n_labels-1][posl] = '\0';
307
/* if(b->n_labels == 1) */
308
/* PhyML_Printf("\n. Read label '%s' on t_edge %3d.",b->labels[0],b->num); */
311
/* PhyML_Printf("\n. Read labels "); */
312
/* For(i,b->n_labels) PhyML_Printf("'%s' ",b->labels[i]); */
313
/* PhyML_Printf("on t_edge %3d.",b->num); */
316
if(!strcmp(b->labels[0],"NULL"))
323
/* PhyML_Printf("\n. No label found on %s",s_d); */
328
//////////////////////////////////////////////////////////////
329
//////////////////////////////////////////////////////////////
331
void Read_Branch_Length(char *s_d, char *s_a, t_tree *tree)
338
b = tree->a_edges[tree->num_curr_branch_available];
340
/* sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
341
sub_tp = (char *)mCalloc(10+strlen(s_d)+1,sizeof(char));
346
strcat(s_d,b->labels[i]);
353
p = strstr(s_a,sub_tp);
361
p = strstr(s_a,sub_tp);
367
b->l->v = atof((char *)p+(int)strlen(sub_tp));
368
tree->has_branch_lengths = YES;
380
//////////////////////////////////////////////////////////////
381
//////////////////////////////////////////////////////////////
383
void Read_Node_Name(t_node *d, char *s_tree_d, t_tree *tree)
387
if(!tree->a_edges[tree->num_curr_branch_available]->n_labels)
389
d->name = (char *)mCalloc(strlen(s_tree_d)+1,sizeof(char ));
390
strcpy(d->name,s_tree_d);
397
d->name = (char *)realloc(d->name,(i+1)*sizeof(char ));
398
d->name[i] = s_tree_d[i];
401
while(s_tree_d[i] != '#');
404
d->ori_name = d->name;
407
//////////////////////////////////////////////////////////////
408
//////////////////////////////////////////////////////////////
411
void Clean_Multifurcation(char **subtrees, int current_deg, int end_deg)
414
if(current_deg <= end_deg) return;
420
/* s_tmp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
421
s_tmp = (char *)mCalloc(10+
422
(int)strlen(subtrees[0])+1+
423
(int)strlen(subtrees[1])+1,
427
strcat(s_tmp,subtrees[0]);
429
strcat(s_tmp,subtrees[1]);
430
strcat(s_tmp,")#NULL\0"); /* Add the label 'NULL' to identify a non-existing edge */
434
for(i=1;i<current_deg-1;i++) strcpy(subtrees[i],subtrees[i+1]);
436
Clean_Multifurcation(subtrees,current_deg-1,end_deg);
440
//////////////////////////////////////////////////////////////
441
//////////////////////////////////////////////////////////////
444
char **Sub_Trees(char *tree, int *degree)
450
if(tree[0] != '(') {*degree = 1; return NULL;}
452
subs=(char **)mCalloc(NODE_DEG_MAX,sizeof(char *));
454
For(i,NODE_DEG_MAX) subs[i]=(char *)mCalloc(strlen(tree)+1,sizeof(char));
462
if(tree[posend] != '(')
464
while((tree[posend] != ',' ) &&
465
(tree[posend] != ':' ) &&
466
(tree[posend] != '#' ) &&
467
(tree[posend] != ')' ))
473
else posend=Next_Par(tree,posend);
475
while((tree[posend+1] != ',') &&
476
(tree[posend+1] != ':') &&
477
(tree[posend+1] != '#') &&
478
(tree[posend+1] != ')')) {posend++;}
481
strncpy(subs[(*degree)],tree+posbeg,posend-posbeg+1);
482
/* strcat(subs[(*degree)],"\0"); */
483
subs[(*degree)][posend-posbeg+1]='\0'; /* Thanks to Jean-Baka Domelevo-Entfellner */
486
while((tree[posend] != ',') &&
487
(tree[posend] != ')')) {posend++;}
492
if((*degree) == NODE_DEG_MAX)
495
PhyML_Printf("\n. Subtree %d : %s\n",i+1,subs[i]);
497
PhyML_Printf("\n. The degree of a t_node cannot be greater than %d\n",NODE_DEG_MAX);
501
while(tree[posend-1] != ')');
507
//////////////////////////////////////////////////////////////
508
//////////////////////////////////////////////////////////////
511
int Next_Par(char *s, int pos)
517
while(*(s+curr) != ')')
519
if(*(s+curr) == '(') curr=Next_Par(s,curr);
526
//////////////////////////////////////////////////////////////
527
//////////////////////////////////////////////////////////////
530
void Print_Tree(FILE *fp, t_tree *tree)
535
s_tree = (char *)Write_Tree(tree,NO);
537
if(OUTPUT_TREE_FORMAT == NEWICK) PhyML_Fprintf(fp,"%s\n",s_tree);
538
else if(OUTPUT_TREE_FORMAT == NEXUS)
540
PhyML_Fprintf(fp,"#NEXUS\n");
541
PhyML_Fprintf(fp,"BEGIN TREES;\n");
542
PhyML_Fprintf(fp,"\tTRANSLATE\n");
543
For(i,tree->n_otu) PhyML_Fprintf(fp,"\t%3d\t%s,\n",i+1,tree->a_nodes[i]->name);
544
PhyML_Fprintf(fp,"\tUTREE PAUP_1=\n");
545
PhyML_Fprintf(fp,"%s\n",s_tree);
546
PhyML_Fprintf(fp,"ENDBLOCK;");
551
//////////////////////////////////////////////////////////////
552
//////////////////////////////////////////////////////////////
555
char *Write_Tree(t_tree *tree, int custom)
561
init_len = 3*(int)T_MAX_NAME;
564
s=(char *)mCalloc(init_len,sizeof(char));
565
available = init_len;
567
s=(char *)mCalloc(T_MAX_LINE,sizeof(char));
578
while((!tree->a_nodes[tree->n_otu+i]->v[0]) ||
579
(!tree->a_nodes[tree->n_otu+i]->v[1]) ||
580
(!tree->a_nodes[tree->n_otu+i]->v[2])) i++;
582
R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[0],&available,&s,tree);
583
R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[1],&available,&s,tree);
584
R_wtree(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[2],&available,&s,tree);
588
R_wtree(tree->n_root,tree->n_root->v[2],&available,&s,tree);
589
R_wtree(tree->n_root,tree->n_root->v[1],&available,&s,tree);
601
while((!tree->a_nodes[tree->n_otu+i]->v[0]) ||
602
(!tree->a_nodes[tree->n_otu+i]->v[1]) ||
603
(!tree->a_nodes[tree->n_otu+i]->v[2])) i++;
605
R_wtree_Custom(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[0],&available,&s,&pos,tree);
606
R_wtree_Custom(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[1],&available,&s,&pos,tree);
607
R_wtree_Custom(tree->a_nodes[tree->n_otu+i],tree->a_nodes[tree->n_otu+i]->v[2],&available,&s,&pos,tree);
611
R_wtree_Custom(tree->n_root,tree->n_root->v[2],&available,&s,&pos,tree);
612
R_wtree_Custom(tree->n_root,tree->n_root->v[1],&available,&s,&pos,tree);
617
s[(int)strlen(s)-1]=')';
618
s[(int)strlen(s)]=';';
623
//////////////////////////////////////////////////////////////
624
//////////////////////////////////////////////////////////////
627
void R_wtree(t_node *pere, t_node *fils, int *available, char **s_tree, t_tree *tree)
632
#if !(defined PHYTIME || defined SERGEII)
636
format = (char *)mCalloc(100,sizeof(char));
638
sprintf(format,"%%.%df",tree->bl_ndigits);
643
/* printf("\n- Writing on %p",*s_tree); */
644
/* ori_len = *pos; */
646
last_len = (int)strlen(*s_tree);
648
if(OUTPUT_TREE_FORMAT == NEWICK)
650
if(tree->write_tax_names == YES)
652
if(tree->io && tree->io->long_tax_names)
654
strcat(*s_tree,tree->io->long_tax_names[fils->num]);
658
strcat(*s_tree,fils->name);
661
else if(tree->write_tax_names == NO)
663
sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
666
else if(OUTPUT_TREE_FORMAT == NEXUS)
668
sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
672
PhyML_Printf("\n== Unknown tree format.");
673
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
674
PhyML_Printf("\n== s=%s\n",*s_tree);
677
if((fils->b) && (fils->b[0]) && (tree->write_br_lens == YES))
679
(*s_tree)[(int)strlen(*s_tree)] = ':';
681
#if !(defined PHYTIME || defined SERGEII)
684
if(tree->is_mixt_tree == NO)
686
mean_len = fils->b[0]->l->v;
688
else mean_len = MIXT_Get_Mean_Edge_Len(fils->b[0],tree);
689
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
693
if(pere == tree->n_root)
695
phydbl root_pos = (fils == tree->n_root->v[2])?(tree->n_root_pos):(1.-tree->n_root_pos);
696
if(tree->is_mixt_tree == NO)
698
mean_len = tree->e_root->l->v;
700
else mean_len = MIXT_Get_Mean_Edge_Len(tree->e_root,tree);
701
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len) * root_pos);
705
if(tree->is_mixt_tree == NO)
707
mean_len = fils->b[0]->l->v;
710
else mean_len = MIXT_Get_Mean_Edge_Len(fils->b[0],tree);
711
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
717
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,fils->b[0]->l->v));
721
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,tree->rates->cur_l[fils->num]));
726
/* strcat(*s_tree,","); */
727
(*s_tree)[(int)strlen(*s_tree)] = ',';
731
(*available) -= ((int)strlen(*s_tree) - last_len);
733
/* printf("\n0 Available = %d [%d %d]",(*available),(int)strlen(*s_tree),last_len); */
734
/* printf("\n0 %s [%d,%d]",*s_tree,(int)(int)strlen(*s_tree),*available); */
738
PhyML_Printf("\n== s=%s\n",*s_tree);
739
PhyML_Printf("\n== len=%d\n",(int)strlen(*s_tree));
740
PhyML_Printf("\n== The sequence names in your input file might be too long.");
741
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
745
if(*available < (int)T_MAX_NAME)
747
(*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
748
For(i,3*(int)T_MAX_NAME) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
749
(*available) = 3*(int)T_MAX_NAME;
750
/* printf("\n. ++ 0 Available = %d",(*available)); */
757
(*s_tree)[(int)strlen(*s_tree)]='(';
762
/* printf("\n1 Available = %d [%d %d]",(*available),(int)strlen(*s_tree),last_len); */
763
/* printf("\n1 %s [%d,%d]",*s_tree,(int)(int)strlen(*s_tree),*available); */
765
if(*available < (int)T_MAX_NAME)
767
(*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
768
For(i,3*(int)T_MAX_NAME) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
769
(*available) = 3*(int)T_MAX_NAME;
770
/* printf("\n. ++ 1 Available = %d",(*available)); */
773
/* (*available)--; */
775
/* if(*available < (int)T_MAX_NAME/2) */
777
/* (*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char)); */
778
/* (*available) = (int)T_MAX_NAME; */
786
if((fils->v[i] != pere) && (fils->b[i] != tree->e_root))
787
R_wtree(fils,fils->v[i],available,s_tree,tree);
795
if(fils->v[i] != pere)
796
R_wtree(fils,fils->v[i],available,s_tree,tree);
803
PhyML_Printf("\n== fils=%p root=%p root->v[2]=%p root->v[1]=%p",fils,tree->n_root,tree->n_root->v[2],tree->n_root->v[1]);
804
PhyML_Printf("\n== tree->e_root=%p fils->b[0]=%p fils->b[1]=%p fils->b[2]=%p",tree->e_root,fils->b[0],fils->b[1],fils->b[2]);
805
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
809
last_len = (int)strlen(*s_tree);
811
(*s_tree)[last_len-1] = ')';
813
if((fils->b) && (tree->write_br_lens == YES))
815
if(tree->print_boot_val)
817
sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->b[p]->bip_score);
819
else if(tree->print_alrt_val)
821
sprintf(*s_tree+(int)strlen(*s_tree),"%f",fils->b[p]->ratio_test);
826
(*s_tree)[(int)strlen(*s_tree)] = ':';
828
#if !(defined PHYTIME || defined SERGEII)
831
if(tree->is_mixt_tree == NO)
833
mean_len = fils->b[p]->l->v;
835
else mean_len = MIXT_Get_Mean_Edge_Len(fils->b[p],tree);
836
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
840
if(pere == tree->n_root)
842
phydbl root_pos = (fils == tree->n_root->v[2])?(tree->n_root_pos):(1.-tree->n_root_pos);
844
if(tree->is_mixt_tree == NO)
846
mean_len = tree->e_root->l->v;
848
else mean_len = MIXT_Get_Mean_Edge_Len(tree->e_root,tree);
849
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len) * root_pos);
853
if(tree->is_mixt_tree == NO)
855
mean_len = fils->b[p]->l->v;
857
else mean_len = MIXT_Get_Mean_Edge_Len(fils->b[p],tree);
858
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,mean_len));
864
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,fils->b[p]->l->v));
868
sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,tree->rates->cur_l[fils->num]));
873
/* strcat(*s_tree,","); */
874
(*s_tree)[(int)strlen(*s_tree)] = ',';
877
(*available) -= ((int)strlen(*s_tree) - last_len);
879
/* printf("\n2 Available = %d [%d %d]",(*available),(int)strlen(*s_tree),last_len); */
880
/* printf("\n2 %s [%d,%d]",*s_tree,(int)(int)strlen(*s_tree),*available); */
884
PhyML_Printf("\n== available = %d",*available);
885
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
889
if(*available < (int)T_MAX_NAME)
891
(*s_tree) = (char *)mRealloc(*s_tree,(int)strlen(*s_tree)+3*(int)T_MAX_NAME,sizeof(char));
892
For(i,3*(int)T_MAX_NAME) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
893
(*available) = 3*(int)T_MAX_NAME;
894
/* printf("\n. ++ 2 Available = %d",(*available)); */
898
/* if(*available < (int)T_MAX_NAME/2) */
900
/* (*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char)); */
901
/* (*available) = (int)T_MAX_NAME; */
908
//////////////////////////////////////////////////////////////
909
//////////////////////////////////////////////////////////////
912
void R_wtree_Custom(t_node *pere, t_node *fils, int *available, char **s_tree, int *pos, t_tree *tree)
917
format = (char *)mCalloc(100,sizeof(char));
919
sprintf(format,"%%.%df",tree->bl_ndigits);
920
/* strcpy(format,"%f"); */
925
/* printf("\n- Writing on %p",*s_tree); */
928
if(OUTPUT_TREE_FORMAT == NEWICK)
930
if(tree->write_tax_names == YES)
932
if(tree->io && tree->io->long_tax_names)
934
strcat(*s_tree,tree->io->long_tax_names[fils->num]);
935
(*pos) += (int)strlen(tree->io->long_tax_names[fils->num]);
939
strcat(*s_tree,fils->name);
940
(*pos) += (int)strlen(fils->name);
943
else if(tree->write_tax_names == NO)
945
(*pos) += sprintf(*s_tree+*pos,"%d",fils->num);
948
else if(OUTPUT_TREE_FORMAT == NEXUS)
950
(*pos) += sprintf(*s_tree+*pos,"%d",fils->num+1);
954
PhyML_Printf("\n== Unknown tree format.");
955
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
956
PhyML_Printf("\n== s=%s\n",*s_tree);
959
if((fils->b) && (fils->b[0]) && (tree->write_br_lens == YES))
964
#if !(defined PHYTIME || defined SERGEII)
967
(*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
971
if(pere == tree->n_root)
973
phydbl root_pos = (fils == tree->n_root->v[2])?(tree->n_root_pos):(1.-tree->n_root_pos);
974
(*pos) += sprintf(*s_tree+*pos,format,tree->e_root->l->v * root_pos);
978
(*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
984
(*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
988
(*pos) += sprintf(*s_tree+*pos,format,tree->rates->cur_l[fils->num]);
993
if(tree->write_labels)
995
if(fils->b[0]->n_labels < 10)
996
For(i,fils->b[0]->n_labels)
998
(*pos) += sprintf(*s_tree+*pos,"::%s",fils->b[0]->labels[i]);
1002
(*pos) += sprintf(*s_tree+*pos,"::%d_labels",fils->b[0]->n_labels);
1007
strcat(*s_tree,",");
1010
(*available) = (*available) - (*pos - ori_len);
1014
PhyML_Printf("\n== s=%s\n",*s_tree);
1015
PhyML_Printf("\n== len=%d\n",strlen(*s_tree));
1016
PhyML_Printf("\n== The sequence names in your input file might be too long.");
1017
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1021
if(*available < (int)T_MAX_NAME)
1023
(*s_tree) = (char *)mRealloc(*s_tree,*pos+3*(int)T_MAX_NAME,sizeof(char));
1024
For(i,3*(int)T_MAX_NAME) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
1025
(*available) = 3*(int)T_MAX_NAME;
1027
/* printf(" %s [%d,%d]",*s_tree,(int)strlen(*s_tree),*available); */
1032
(*s_tree)[(*pos)]='(';
1033
(*s_tree)[(*pos)+1]='\0';
1037
if(*available < (int)T_MAX_NAME/2)
1039
(*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char));
1040
(*available) = (int)T_MAX_NAME;
1047
if((fils->v[i] != pere) && (fils->b[i] != tree->e_root))
1048
R_wtree_Custom(fils,fils->v[i],available,s_tree,pos,tree);
1056
if(fils->v[i] != pere)
1057
R_wtree_Custom(fils,fils->v[i],available,s_tree,pos,tree);
1066
PhyML_Printf("\n== fils=%p root=%p root->v[2]=%p root->v[1]=%p",fils,tree->n_root,tree->n_root->v[2],tree->n_root->v[1]);
1067
PhyML_Printf("\n== tree->e_root=%p fils->b[0]=%p fils->b[1]=%p fils->b[2]=%p",tree->e_root,fils->b[0],fils->b[1],fils->b[2]);
1068
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1072
/* printf("\n+ Writing on %p",*s_tree); */
1073
(*s_tree)[(*pos)-1] = ')';
1074
(*s_tree)[(*pos)] = '\0';
1076
if((fils->b) && (tree->write_br_lens == YES))
1078
if(tree->print_boot_val)
1080
(*pos) += sprintf(*s_tree+*pos,"%d",fils->b[p]->bip_score);
1082
else if(tree->print_alrt_val)
1084
(*pos) += sprintf(*s_tree+*pos,"%f",fils->b[p]->ratio_test);
1089
strcat(*s_tree,":");
1092
#if !(defined PHYTIME || defined SERGEII)
1095
(*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1099
if(pere == tree->n_root)
1101
phydbl root_pos = (fils == tree->n_root->v[2])?(tree->n_root_pos):(1.-tree->n_root_pos);
1102
(*pos) += sprintf(*s_tree+*pos,format,tree->e_root->l->v * root_pos);
1106
(*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1112
(*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1116
(*pos) += sprintf(*s_tree+*pos,format,tree->rates->cur_l[fils->num]);
1122
if((tree->write_labels) && (fils->b[p]->labels != NULL))
1124
if(fils->b[p]->n_labels < 10)
1125
For(i,fils->b[p]->n_labels)
1127
(*pos) += sprintf(*s_tree+*pos,"::%s",fils->b[p]->labels[i]);
1131
(*pos) += sprintf(*s_tree+*pos,"::%d_labels",fils->b[p]->n_labels);
1135
strcat(*s_tree,",");
1137
(*available) = (*available) - (*pos - ori_len);
1141
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1145
if(*available < (int)T_MAX_NAME)
1147
(*s_tree) = (char *)mRealloc(*s_tree,*pos+3*(int)T_MAX_NAME,sizeof(char));
1148
For(i,3*(int)T_MAX_NAME) (*s_tree)[(int)strlen(*s_tree)+i] = '\0';
1149
(*available) = 3*(int)T_MAX_NAME;
1151
/* printf(" %s [%d,%d]",*s_tree,(int)strlen(*s_tree),*available); */
1157
//////////////////////////////////////////////////////////////
1158
//////////////////////////////////////////////////////////////
1160
void Detect_Align_File_Format(option *io)
1165
fgetpos(io->fp_in_align,&curr_pos);
1169
while((c=fgetc(io->fp_in_align)) != EOF)
1171
if(errno) io->data_file_format = PHYLIP;
1174
char s[10],t[6]="NEXUS";
1175
if(!fgets(s,6,io->fp_in_align))
1177
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1182
fsetpos(io->fp_in_align,&curr_pos);
1183
io->data_file_format = NEXUS;
1189
fsetpos(io->fp_in_align,&curr_pos);
1192
//////////////////////////////////////////////////////////////
1193
//////////////////////////////////////////////////////////////
1195
void Detect_Tree_File_Format(option *io)
1200
fgetpos(io->fp_in_tree,&curr_pos);
1204
while((c=fgetc(io->fp_in_tree)) != EOF)
1208
io->tree_file_format = PHYLIP;
1209
PhyML_Printf("\n. Detected PHYLIP tree file format.");
1213
char s[10],t[6]="NEXUS";
1214
if(!fgets(s,6,io->fp_in_tree))
1216
PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1221
fsetpos(io->fp_in_tree,&curr_pos);
1222
io->tree_file_format = NEXUS;
1223
PhyML_Printf("\n. Detected NEXUS tree file format.");
1229
fsetpos(io->fp_in_tree,&curr_pos);
1232
//////////////////////////////////////////////////////////////
1233
//////////////////////////////////////////////////////////////
1235
align **Get_Seq(option *io)
1239
if(!io->fp_in_align)
1241
PhyML_Printf("\n== Filehandle to '%s' seems to be closed.",io->in_align_file);
1245
Detect_Align_File_Format(io);
1247
switch(io->data_file_format)
1251
io->data = Get_Seq_Phylip(io);
1256
io->nex_com_list = Make_Nexus_Com();
1257
Init_Nexus_Format(io->nex_com_list);
1258
Get_Nexus_Data(io->fp_in_align,io);
1264
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1272
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1277
Post_Process_Data(io);
1284
//////////////////////////////////////////////////////////////
1285
//////////////////////////////////////////////////////////////
1287
void Post_Process_Data(option *io)
1291
For(i,io->data[0]->len)
1295
if((io->data[j]->state[i] == '?') || (io->data[j]->state[i] == '-')) io->data[j]->state[i] = 'X';
1296
if((io->datatype == NT) && (io->data[j]->state[i] == 'N')) io->data[j]->state[i] = 'X';
1297
if(io->data[j]->state[i] == 'U') io->data[j]->state[i] = 'T';
1301
For(i,io->n_otu) io->data[i]->len = io->data[0]->len;
1304
/* align **Get_Seq_Nexus(option *io) */
1306
/* char *s,*ori_s; */
1308
/* int in_comment; */
1309
/* nexcom *curr_com; */
1310
/* nexparm *curr_parm; */
1311
/* int nxt_token_t,cur_token_t; */
1313
/* s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
1314
/* token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char)); */
1317
/* in_comment = NO; */
1318
/* curr_com = NULL; */
1319
/* curr_parm = NULL; */
1320
/* nxt_token_t = NEXUS_COM; */
1321
/* cur_token_t = -1; */
1323
/* while(fgets(s,T_MAX_LINE,io->fp_in_align)) */
1327
/* Get_Token(&s,token); */
1329
/* /\* PhyML_Printf("\n. Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); *\/ */
1331
/* if(token[0] == '\0') break; */
1333
/* if(token[0] == ';') */
1335
/* curr_com = NULL; */
1336
/* curr_parm = NULL; */
1337
/* nxt_token_t = NEXUS_COM; */
1338
/* cur_token_t = -1; */
1339
/* break; /\* End of command *\/ */
1342
/* if(nxt_token_t == NEXUS_EQUAL) */
1344
/* cur_token_t = NEXUS_VALUE; */
1345
/* nxt_token_t = NEXUS_PARM; */
1349
/* if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE)) */
1351
/* Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list); */
1354
/* nxt_token_t = curr_com->nxt_token_t; */
1355
/* cur_token_t = curr_com->cur_token_t; */
1357
/* if(cur_token_t != NEXUS_VALUE) continue; */
1360
/* if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE)) */
1362
/* Find_Nexus_Parm(token,&curr_parm,curr_com); */
1365
/* nxt_token_t = curr_parm->nxt_token_t; */
1366
/* cur_token_t = curr_parm->cur_token_t; */
1368
/* if(cur_token_t != NEXUS_VALUE) continue; */
1371
/* if(cur_token_t == NEXUS_VALUE) */
1373
/* if((curr_parm->fp)(token,curr_parm,io)) /\* Read in parameter value *\/ */
1375
/* nxt_token_t = NEXUS_PARM; */
1376
/* cur_token_t = -1; */
1380
/* while(strlen(token) > 0); */
1386
/* return io->data; */
1389
/* /\*********************************************************\/ */
1391
void Get_Nexus_Data(FILE *fp, option *io)
1396
int nxt_token_t,cur_token_t;
1398
token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char));
1402
nxt_token_t = NEXUS_COM;
1407
if(!Get_Token(fp,token)) break;
1409
/* PhyML_Printf("\n+ Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); */
1415
nxt_token_t = NEXUS_COM;
1419
if(nxt_token_t == NEXUS_EQUAL)
1421
cur_token_t = NEXUS_VALUE;
1422
nxt_token_t = NEXUS_PARM;
1426
if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE))
1428
Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list);
1431
nxt_token_t = curr_com->nxt_token_t;
1432
cur_token_t = curr_com->cur_token_t;
1434
if(cur_token_t != NEXUS_VALUE) continue;
1437
if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE))
1439
Find_Nexus_Parm(token,&curr_parm,curr_com);
1442
nxt_token_t = curr_parm->nxt_token_t;
1443
cur_token_t = curr_parm->cur_token_t;
1445
if(cur_token_t != NEXUS_VALUE) continue;
1448
if(cur_token_t == NEXUS_VALUE)
1450
if((curr_parm->fp)(token,curr_parm,io)) /* Read in parameter value */
1452
nxt_token_t = NEXUS_PARM;
1457
while(strlen(token) > 0);
1462
//////////////////////////////////////////////////////////////
1463
//////////////////////////////////////////////////////////////
1466
int Get_Token(FILE *fp, char *token)
1471
while(c == ' ' || c == '\t' || c == '\n')
1474
if(c == EOF) return 0;
1484
if(c == EOF) return 0;
1488
/* c = fgetc(fp); */
1489
if(c == EOF) return 0;
1500
if(c == EOF) return 0;
1504
if(c == '#') { *token = c; token++; }
1505
else if(c == ';') { *token = c; token++; }
1506
else if(c == ',') { *token = c; token++; }
1507
else if(c == '.') { *token = c; token++; }
1508
else if(c == '=') { *token = c; token++; }
1509
else if(c == '(') { *token = c; token++; }
1510
else if(c == ')') { *token = c; token++; }
1511
else if(c == '{') { *token = c; token++; }
1512
else if(c == '}') { *token = c; token++; }
1513
else if(c == '?') { *token = c; token++; }
1514
else if(c == '-') { *token = c; token++; }
1517
while(isgraph(c) && c != ';' && c != '-' && c != ',' && c != '=')
1521
if(c == EOF) return 0;
1524
fseek(fp,-1*sizeof(char),SEEK_CUR);
1532
/* void Get_Token(char *line, char *token) */
1534
/* while(**line == ' ' || **line == '\t') (*line)++; */
1537
/* if(**line == '"') */
1539
/* do { *token = **line; (*line)++; token++; } while(**line != '"'); */
1540
/* *token = **line; */
1542
/* *(token+1) = '\0'; */
1546
/* if(**line == '[') */
1548
/* int in_comment; */
1550
/* in_comment = 1; */
1554
/* if(**line == '[') */
1558
/* else if(**line == ']') in_comment--; */
1560
/* while(in_comment); */
1566
/* if(**line == '#') {*token = **line; (*line)++; token++; } */
1567
/* else if(**line == ';') {*token = **line; (*line)++; token++; } */
1568
/* else if(**line == ',') {*token = **line; (*line)++; token++; } */
1569
/* else if(**line == '.') {*token = **line; (*line)++; token++; } */
1570
/* else if(**line == '=') {*token = **line; (*line)++; token++; } */
1571
/* else if(**line == '(') {*token = **line; (*line)++; token++; } */
1572
/* else if(**line == ')') {*token = **line; (*line)++; token++; } */
1573
/* else if(**line == '{') {*token = **line; (*line)++; token++; } */
1574
/* else if(**line == '}') {*token = **line; (*line)++; token++; } */
1575
/* else if(**line == '?') {*token = **line; (*line)++; token++; } */
1576
/* else if(**line == '-') {*token = **line; (*line)++; token++; } */
1579
/* while(isgraph(**line) && **line != ';' && **line != '=' && **line != ',') */
1581
/* *(token++) = **line; */
1585
/* *token = '\0'; */
1588
//////////////////////////////////////////////////////////////
1589
//////////////////////////////////////////////////////////////
1591
align **Get_Seq_Phylip(option *io)
1593
Read_Ntax_Len_Phylip(io->fp_in_align,&io->n_otu,&io->init_len);
1595
if(io->n_otu > N_MAX_OTU)
1597
PhyML_Printf("\n== The number of taxa should not exceed %d",N_MAX_OTU);
1601
if(io->interleaved) io->data = Read_Seq_Interleaved(io);
1602
else io->data = Read_Seq_Sequential(io);
1607
//////////////////////////////////////////////////////////////
1608
//////////////////////////////////////////////////////////////
1610
void Read_Ntax_Len_Phylip(FILE *fp ,int *n_otu, int *n_tax)
1615
line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1620
if(fscanf(fp,"%s",line) == EOF)
1623
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1628
if(strcmp(line,"\n") && strcmp(line,"\r") && strcmp(line,"\t"))
1630
sscanf(line,"%d",n_otu);
1631
if(*n_otu <= 0) Warn_And_Exit("\n. The number of taxa cannot be negative.\n");
1633
if(!fscanf(fp,"%s",line)) Exit("\n");
1634
sscanf(line,"%d",n_tax);
1635
if(*n_tax <= 0) Warn_And_Exit("\n. The sequence length cannot be negative.\n");
1644
//////////////////////////////////////////////////////////////
1645
//////////////////////////////////////////////////////////////
1647
align **Read_Seq_Sequential(option *io)
1655
format = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1656
line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1657
data = (align **)mCalloc(io->n_otu,sizeof(align *));
1659
/* while((c=fgetc(in))!='\n'); */
1660
/* while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1662
sprintf(format, "%%%ds", T_MAX_NAME);
1666
data[i] = (align *)mCalloc(1,sizeof(align));
1667
data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1668
data[i]->state = (char *)mCalloc(io->init_len*io->state_len+1,sizeof(char));
1670
data[i]->is_ambigu = NULL;
1673
if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1675
Check_Sequence_Name(data[i]->name);
1677
while(data[i]->len < io->init_len * io->state_len) Read_One_Line_Seq(&data,i,io->fp_in_align);
1679
if(data[i]->len != io->init_len * io->state_len)
1681
PhyML_Printf("\n== Err: Problem with species %s's sequence (check the format).\n",data[i]->name);
1682
PhyML_Printf("\n== Observed sequence length: %d, expected length: %d\n",data[i]->len, io->init_len * io->state_len);
1687
For(i,io->n_otu) data[i]->state[data[i]->len] = '\0';
1689
Restrict_To_Coding_Position(data,io);
1697
//////////////////////////////////////////////////////////////
1698
//////////////////////////////////////////////////////////////
1700
align **Read_Seq_Interleaved(option *io)
1702
int i,end,num_block;
1709
line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1710
format = (char *)mCalloc(T_MAX_NAME, sizeof(char));
1711
data = (align **)mCalloc(io->n_otu,sizeof(align *));
1713
/* while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1715
sprintf(format, "%%%ds", T_MAX_NAME);
1720
data[i] = (align *)mCalloc(1,sizeof(align));
1721
data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1722
data[i]->state = (char *)mCalloc(io->init_len*io->state_len+1,sizeof(char));
1725
data[i]->is_ambigu = NULL;
1727
if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1729
Check_Sequence_Name(data[i]->name);
1731
if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1734
if((i != io->n_otu) && (i != io->n_otu-1))
1736
PhyML_Printf("\n== i:%d n_otu:%d",i,io->n_otu);
1737
PhyML_Printf("\n== Err.: problem with species %s's sequence.\n",data[i]->name);
1738
PhyML_Printf("\n== Observed sequence length: %d, expected length: %d\n",data[i]->len, io->init_len * io->state_len);
1745
if(data[0]->len == io->init_len * io->state_len) end = 1;
1747
/* if(end) printf("\n. finished yet '%c'\n",fgetc(io->fp_in_align)); */
1759
if(!fgets(line,T_MAX_LINE,io->fp_in_align)) break;
1761
if(line[0] != 13 && line[0] != 10)
1763
PhyML_Printf("\n== Err.: one or more missing sequences in block %d.\n",num_block-1);
1767
For(i,io->n_otu) if(data[i]->len != io->init_len * io->state_len) break;
1769
if(i == io->n_otu) break;
1773
/* Skip the taxon name, if any, in this interleaved block */
1774
fgetpos(io->fp_in_align,&curr_pos);
1775
if(!fscanf(io->fp_in_align,format,line))
1777
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1780
if(line && strcmp(line,data[i]->name)) fsetpos(io->fp_in_align,&curr_pos);
1782
if(data[i]->len > io->init_len * io->state_len)
1784
PhyML_Printf("\n== Observed sequence length=%d expected length=%d.\n",data[i]->len,io->init_len * io->state_len);
1785
PhyML_Printf("\n== Err.: Problem with species %s's sequence.\n",data[i]->name);
1788
else if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1791
if((i != io->n_otu) && (i != io->n_otu-1))
1793
PhyML_Printf("\n== Err.: Problem with species %s's sequence.\n",data[i]->name);
1794
PhyML_Printf("\n== Observed sequence length: %d, expected length: %d.\n",data[i]->len, io->init_len * io->state_len);
1803
For(i,io->n_otu) data[i]->state[data[i]->len] = '\0';
1807
if(data[i]->len != io->init_len * io->state_len)
1809
PhyML_Printf("\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);
1814
Restrict_To_Coding_Position(data,io);
1822
//////////////////////////////////////////////////////////////
1823
//////////////////////////////////////////////////////////////
1826
int Read_One_Line_Seq(align ***data, int num_otu, FILE *in)
1833
/* if((c == EOF) || (c == '\n') || (c == '\r')) break; */
1835
if((c == 13) || (c == 10))
1837
/* PhyML_Printf("[%d %d]\n",c,nchar); fflush(NULL); */
1845
/* PhyML_Printf("break\n"); */
1851
/* PhyML_Printf("EOL\n"); */
1854
else if((c == ' ') || (c == '\t') || (c == 32))
1856
/* PhyML_Printf("[%d]",c); */
1866
c = (*data)[0]->state[(*data)[num_otu]->len];
1868
Warn_And_Exit("\n== Err: Symbol \".\" should not appear in the first sequence\n");
1870
(*data)[num_otu]->state[(*data)[num_otu]->len]=c;
1871
(*data)[num_otu]->len++;
1872
c = (char)fgetc(in);
1873
/* PhyML_Printf("[%c %d]",c,c); */
1877
/* printf("\n. Exit nchar: %d [%d]\n",nchar,c==EOF); */
1878
if(c == EOF) return 0;
1882
//////////////////////////////////////////////////////////////
1883
//////////////////////////////////////////////////////////////
1885
t_tree *Read_Tree_File(option *io)
1891
PhyML_Printf("\n== Filehandle to '%s' seems to be closed.",io->in_tree_file);
1896
Detect_Tree_File_Format(io);
1898
io->treelist->list_size = 0;
1900
switch(io->tree_file_format)
1906
io->treelist->tree = (t_tree **)realloc(io->treelist->tree,(io->treelist->list_size+1)*sizeof(t_tree *));
1907
io->tree = Read_Tree_File_Phylip(io->fp_in_tree);
1908
if(!io->tree) break;
1909
if(io->treelist->list_size > 1) PhyML_Printf("\n. Reading tree %d",io->treelist->list_size+1);
1910
io->treelist->tree[io->treelist->list_size] = io->tree;
1911
io->treelist->list_size++;
1917
io->nex_com_list = Make_Nexus_Com();
1918
Init_Nexus_Format(io->nex_com_list);
1919
Get_Nexus_Data(io->fp_in_tree,io);
1925
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1931
if(!io->long_tax_names)
1935
tree = io->treelist->tree[0];
1937
io->long_tax_names = (char **)mCalloc(tree->n_otu,sizeof(char *));
1938
io->short_tax_names = (char **)mCalloc(tree->n_otu,sizeof(char *));
1942
io->long_tax_names[i] = (char *)mCalloc(strlen(tree->a_nodes[i]->name)+1,sizeof(char));
1943
io->short_tax_names[i] = (char *)mCalloc(strlen(tree->a_nodes[i]->name)+1,sizeof(char));
1944
strcpy(io->long_tax_names[i],tree->a_nodes[i]->name);
1945
strcpy(io->short_tax_names[i],tree->a_nodes[i]->name);
1951
//////////////////////////////////////////////////////////////
1952
//////////////////////////////////////////////////////////////
1954
char *Return_Tree_String_Phylip(FILE *fp_input_tree)
1961
if(fp_input_tree == NULL)
1963
PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1969
c=fgetc(fp_input_tree);
1971
while((c != '(') && (c != EOF));
1974
if(c==EOF) return NULL;
1976
line = (char *)mCalloc(1,sizeof(char));
1983
if((c == ' ') || (c == '\n'))
1985
c=fgetc(fp_input_tree);
1986
if(c == EOF || c == ';') break;
1992
Skip_Comment(fp_input_tree);
1993
c = fgetc(fp_input_tree);
1994
if(c == EOF || c == ';') break;
1997
line = (char *)mRealloc(line,i+2,sizeof(char));
2001
c=fgetc(fp_input_tree);
2002
if(c==EOF || c==';') break;
2005
if(open>maxopen) maxopen = open;
2010
/* if(maxopen == 1) return NULL; */
2014
//////////////////////////////////////////////////////////////
2015
//////////////////////////////////////////////////////////////
2018
t_tree *Read_Tree_File_Phylip(FILE *fp_input_tree)
2023
line = Return_Tree_String_Phylip(fp_input_tree);
2024
tree = Read_Tree(&line);
2030
//////////////////////////////////////////////////////////////
2031
//////////////////////////////////////////////////////////////
2033
void Print_Site(calign *cdata, int num, int n_otu, char *sep, int stepsize, FILE *fp)
2036
PhyML_Fprintf(fp,"\n");
2039
PhyML_Fprintf(fp,"%20s ",cdata->c_seq[i]->name);
2041
PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[num+j]);
2042
PhyML_Fprintf(fp,"%s",sep);
2044
PhyML_Fprintf(fp,"%s",sep);
2047
//////////////////////////////////////////////////////////////
2048
//////////////////////////////////////////////////////////////
2050
void Print_Site_Lk(t_tree *tree, FILE *fp)
2057
if(!tree->io->print_site_lnl)
2059
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
2063
if(!tree->io->print_trace)
2065
s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2067
PhyML_Fprintf(fp,"Note : P(D|M) is the probability of site D given the model M (i.e., the site likelihood)\n");
2068
if(tree->mod->ras->n_catg > 1 || tree->mod->ras->invar)
2069
PhyML_Fprintf(fp,"P(D|M,rr[x]) is the probability of site D given the model M and the relative rate\nof evolution rr[x], where x is the class of rate to be considered.\nWe have P(D|M) = \\sum_x P(x) x P(D|M,rr[x]).\n");
2070
PhyML_Fprintf(fp,"\n\n");
2073
PhyML_Fprintf(fp, "%-7s",s);
2075
sprintf(s,"Pattern");
2076
PhyML_Fprintf(fp, "%-9s",s);
2078
sprintf(s,"P(D|M)");
2079
PhyML_Fprintf(fp,"%-16s",s);
2081
if(tree->mod->ras->n_catg > 1)
2083
For(catg,tree->mod->ras->n_catg)
2085
sprintf(s,"P(D|M,rr[%d]=%5.4f)",catg+1,tree->mod->ras->gamma_rr->v[catg]);
2086
PhyML_Fprintf(fp,"%-22s",s);
2089
sprintf(s,"Posterior mean");
2090
PhyML_Fprintf(fp,"%-22s",s);
2094
if(tree->mod->ras->invar)
2096
sprintf(s,"P(D|M,rr[0]=0)");
2097
PhyML_Fprintf(fp,"%-16s",s);
2099
PhyML_Fprintf(fp,"\n");
2101
For(site,tree->data->init_len)
2105
PhyML_Fprintf(fp,"%-7d",site+1);
2106
PhyML_Fprintf(fp,"%-9d",tree->data->sitepatt[site]);
2108
PhyML_Fprintf(fp,"%-16g",tree->cur_site_lk[tree->data->sitepatt[site]]);
2109
if(tree->mod->ras->n_catg > 1)
2111
For(catg,tree->mod->ras->n_catg)
2112
PhyML_Fprintf(fp,"%-22g",(phydbl)EXP(tree->log_site_lk_cat[catg][tree->data->sitepatt[site]]));
2115
For(catg,tree->mod->ras->n_catg)
2117
tree->mod->ras->gamma_rr->v[catg] *
2118
EXP(tree->log_site_lk_cat[catg][tree->data->sitepatt[site]]) *
2119
tree->mod->ras->gamma_r_proba->v[catg];
2120
postmean /= tree->cur_site_lk[tree->data->sitepatt[site]];
2122
PhyML_Fprintf(fp,"%-22g",postmean);
2124
if(tree->mod->ras->invar)
2126
if((phydbl)tree->data->invar[tree->data->sitepatt[site]] > -0.5)
2127
PhyML_Fprintf(fp,"%-16g",tree->mod->e_frq->pi->v[tree->data->invar[tree->data->sitepatt[site]]]);
2129
PhyML_Fprintf(fp,"%-16g",0.0);
2134
PhyML_Fprintf(fp,"\n");
2140
For(site,tree->data->init_len)
2141
PhyML_Fprintf(fp,"%.2f\t",LOG(tree->cur_site_lk[tree->data->sitepatt[site]]));
2142
PhyML_Fprintf(fp,"\n");
2146
//////////////////////////////////////////////////////////////
2147
//////////////////////////////////////////////////////////////
2149
void Print_Seq(FILE *fp, align **data, int n_otu)
2153
PhyML_Fprintf(fp,"%d\t%d\n",n_otu,data[0]->len);
2158
if(j<(int)strlen(data[i]->name))
2159
fputc(data[i]->name[j],fp);
2162
/* PhyML_Printf("%10d ",i); */
2165
PhyML_Fprintf(fp,"%c",data[i]->state[j]);
2167
PhyML_Fprintf(fp,"\n");
2171
//////////////////////////////////////////////////////////////
2172
//////////////////////////////////////////////////////////////
2175
void Print_CSeq(FILE *fp, int compressed, calign *cdata)
2180
n_otu = cdata->n_otu;
2181
if(cdata->format == 0)
2183
PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
2187
PhyML_Fprintf(fp,"#NEXUS\n");
2188
PhyML_Fprintf(fp,"begin data\n");
2189
PhyML_Fprintf(fp,"dimensions ntax=%d nchar=%d;\n",n_otu,cdata->init_len);
2190
PhyML_Fprintf(fp,"format sequential datatype=dna;\n");
2191
PhyML_Fprintf(fp,"matrix\n");
2198
if(j<(int)strlen(cdata->c_seq[i]->name))
2199
fputc(cdata->c_seq[i]->name[j],fp);
2203
if(compressed == YES) /* Print out compressed sequences */
2204
PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
2205
else /* Print out uncompressed sequences */
2207
For(j,cdata->init_len)
2209
PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
2212
PhyML_Fprintf(fp,"\n");
2214
PhyML_Fprintf(fp,"\n");
2216
if(cdata->format == 1)
2218
PhyML_Fprintf(fp,";\n");
2219
PhyML_Fprintf(fp,"END;\n");
2223
/* PhyML_Printf("\t"); */
2224
/* For(j,cdata->crunch_len) */
2225
/* PhyML_Printf("%.0f ",cdata->wght[j]); */
2226
/* PhyML_Printf("\n"); */
2229
//////////////////////////////////////////////////////////////
2230
//////////////////////////////////////////////////////////////
2233
void Print_CSeq_Select(FILE *fp, int compressed, calign *cdata, t_tree *tree)
2244
if(tree->rates->nd_t[i] < tree->rates->time_slice_lims[slice] + eps)
2247
PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
2249
n_otu = cdata->n_otu;
2253
if(tree->rates->nd_t[i] < tree->rates->time_slice_lims[slice] + eps)
2257
if(j<(int)strlen(cdata->c_seq[i]->name))
2258
fputc(cdata->c_seq[i]->name[j],fp);
2262
if(compressed == YES) /* Print out compressed sequences */
2263
PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
2264
else /* Print out uncompressed sequences */
2266
For(j,cdata->init_len)
2268
PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
2271
PhyML_Fprintf(fp,"\n");
2275
if(cdata->format == 1)
2277
PhyML_Fprintf(fp,";\n");
2278
PhyML_Fprintf(fp,"END;\n");
2282
/* PhyML_Printf("\t"); */
2283
/* For(j,cdata->crunch_len) */
2284
/* PhyML_Printf("%.0f ",cdata->wght[j]); */
2285
/* PhyML_Printf("\n"); */
2288
//////////////////////////////////////////////////////////////
2289
//////////////////////////////////////////////////////////////
2291
void Print_Dist(matrix *mat)
2297
PhyML_Printf("%s ",mat->name[i]);
2300
PhyML_Printf("%9.6f ",mat->dist[i][j]);
2305
//////////////////////////////////////////////////////////////
2306
//////////////////////////////////////////////////////////////
2309
void Print_Node(t_node *a, t_node *d, t_tree *tree)
2314
For(i,3) if(a->v[i] == d) {dir = i; break;}
2315
PhyML_Printf("Node nums: %3d %3d (dir:%3d) (anc:%3d) ta:%6.3f td:%6.3f;",
2316
a->num,d->num,dir,a->anc?a->anc->num:(-1),
2317
tree->rates?tree->rates->nd_t[a->num]:-1.,
2318
tree->rates?tree->rates->nd_t[d->num]:-1.);
2319
PhyML_Printf("Node names = '%10s' '%10s' ; ",a->name,d->name);
2320
For(i,3) if(a->v[i] == d)
2324
PhyML_Printf("Branch num = %3d%c (%3d %3d) %f",
2325
a->b[i]->num,a->b[i]==tree->e_root?'*':' ',a->b[i]->left->num,
2326
a->b[i]->rght->num,a->b[i]->l->v);
2327
if(a->b[i]->left->tax) PhyML_Printf(" WARNING LEFT->TAX!");
2336
if(d->v[i] != a && d->b[i] != tree->e_root) Print_Node(d,d->v[i],tree);
2339
//////////////////////////////////////////////////////////////
2340
//////////////////////////////////////////////////////////////
2342
void Print_Model(t_mod *mod)
2346
PhyML_Printf("\n. name=%s",mod->modelname->s);
2347
PhyML_Printf("\n. string=%s",mod->custom_mod_string);
2348
PhyML_Printf("\n. mod_num=%d",mod->mod_num);
2349
PhyML_Printf("\n. ns=%d",mod->ns);
2350
PhyML_Printf("\n. n_catg=%d",mod->ras->n_catg);
2351
PhyML_Printf("\n. kappa=%f",mod->kappa->v);
2352
PhyML_Printf("\n. alpha=%f",mod->ras->alpha->v);
2353
PhyML_Printf("\n. lambda=%f",mod->lambda->v);
2354
PhyML_Printf("\n. pinvar=%f",mod->ras->pinvar->v);
2355
PhyML_Printf("\n. br_len_multiplier=%f",mod->br_len_multiplier->v);
2356
PhyML_Printf("\n. whichmodel=%d",mod->whichmodel);
2357
PhyML_Printf("\n. update_eigen=%d",mod->update_eigen);
2358
PhyML_Printf("\n. bootstrap=%d",mod->bootstrap);
2359
PhyML_Printf("\n. n_diff_rr=%d",mod->r_mat->n_diff_rr);
2360
PhyML_Printf("\n. invar=%d",mod->ras->invar);
2361
PhyML_Printf("\n. use_m4mod=%d",mod->use_m4mod);
2362
PhyML_Printf("\n. gamma_median=%d",mod->ras->gamma_median);
2363
PhyML_Printf("\n. state_len=%d",mod->io->state_len);
2364
PhyML_Printf("\n. log_l=%d",mod->log_l);
2365
PhyML_Printf("\n. l_min=%f",mod->l_min);
2366
PhyML_Printf("\n. l_max=%f",mod->l_max);
2367
PhyML_Printf("\n. free_mixt_rates=%d",mod->ras->free_mixt_rates);
2368
PhyML_Printf("\n. gamma_mgf_bl=%d",mod->gamma_mgf_bl);
2370
PhyML_Printf("\n. Pi\n");
2371
For(i,mod->ns) PhyML_Printf(" %f ",mod->e_frq->pi->v[i]);
2373
For(i,mod->ns) PhyML_Printf(" %f ",mod->e_frq->pi_unscaled->v[i]);
2375
PhyML_Printf("\n. Rates\n");
2376
For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_r_proba->v[i]);
2378
For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_r_proba_unscaled->v[i]);
2380
For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_rr->v[i]);
2382
For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_rr_unscaled->v[i]);
2386
PhyML_Printf("\n. Qmat \n");
2387
if(mod->whichmodel == CUSTOM)
2390
For(i,6) {PhyML_Printf(" %12f ",mod->r_mat->rr->v[i]); fflush(NULL);}
2391
For(i,6) {PhyML_Printf(" %12f ",mod->r_mat->rr_val->v[i]); fflush(NULL);}
2392
For(i,6) {PhyML_Printf(" %12d ",mod->r_mat->rr_num->v[i]); fflush(NULL);}
2393
For(i,6) {PhyML_Printf(" %12d ",mod->r_mat->n_rr_per_cat->v[i]); fflush(NULL);}
2399
PhyML_Printf("%8.5f ",mod->r_mat->qmat->v[i*4+j]);
2403
PhyML_Printf("\n. Freqs");
2405
For(i,mod->ns) PhyML_Printf(" %12f ",mod->user_b_freq->v[i]);
2407
For(i,mod->ns) PhyML_Printf(" %12f ",mod->e_frq->pi->v[i]);
2409
For(i,mod->ns) PhyML_Printf(" %12f ",mod->e_frq->pi_unscaled->v[i]);
2411
PhyML_Printf("\n. Eigen\n");
2412
For(i,2*mod->ns) PhyML_Printf(" %f ",mod->eigen->space[i]);
2414
For(i,2*mod->ns) PhyML_Printf(" %f ",mod->eigen->space_int[i]);
2416
For(i,mod->ns) PhyML_Printf(" %f ",mod->eigen->e_val[i]);
2418
For(i,mod->ns) PhyML_Printf(" %f ",mod->eigen->e_val_im[i]);
2420
For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect[i]);
2422
For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect[i]);
2424
For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect_im[i]);
2426
For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->l_e_vect[i]);
2428
For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->q[i]);
2431
PhyML_Printf("\n. Pij");
2432
For(k,mod->ras->n_catg)
2434
PMat(0.01*mod->ras->gamma_rr->v[k],mod,mod->ns*mod->ns*k,mod->Pij_rr->v);
2435
PhyML_Printf("\n. l=%f\n",0.01*mod->ras->gamma_rr->v[k]);
2440
PhyML_Printf("%8.5f ",mod->Pij_rr->v[k*mod->ns*mod->ns+i*mod->ns+j]);
2451
//////////////////////////////////////////////////////////////
2452
//////////////////////////////////////////////////////////////
2454
void Print_Mat(matrix *mat)
2458
PhyML_Printf("%d",mat->n_otu);
2465
if(j>=(int)strlen(mat->name[i])) putchar(' ');
2466
else putchar(mat->name[i][j]);
2472
if(mat->dist[i][j] < .0)
2473
PhyML_Printf("%12s",s);
2475
PhyML_Printf("%12f",mat->dist[i][j]);
2481
//////////////////////////////////////////////////////////////
2482
//////////////////////////////////////////////////////////////
2484
FILE *Openfile(char *filename, int mode)
2486
/* mode = 0 -> read */
2487
/* mode = 1 -> write */
2488
/* mode = 2 -> append */
2494
/* s = (char *)mCalloc(T_MAX_FILE,sizeof(char)); */
2496
/* strcpy(s,filename); */
2506
while(!(fp = (FILE *)fopen(s,"r")) && ++open_test<10)
2508
PhyML_Printf("\n. Can't open file '%s', enter a new name : ",s);
2515
fp = (FILE *)fopen(s,"w");
2520
fp = (FILE *)fopen(s,"a");
2533
//////////////////////////////////////////////////////////////
2534
//////////////////////////////////////////////////////////////
2536
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)
2542
/* For(i,2*tree->n_otu-3) fprintf(fp_out,"\n. * Edge %3d: %f",i,tree->a_edges[i]->l); */
2544
if((!n_data_set) || (!num_tree))
2547
Print_Banner_Small(fp_out);
2550
PhyML_Fprintf(fp_out,"\n. Sequence filename: \t\t\t%s", Basename(io->in_align_file));
2551
PhyML_Fprintf(fp_out,"\n. Data set: \t\t\t\t#%d",n_data_set);
2553
if(io->mod->s_opt->random_input_tree)
2554
PhyML_Fprintf(fp_out,"\n. Random init tree: \t\t\t#%d",num_tree+1);
2555
else if(io->n_trees > 1)
2556
PhyML_Fprintf(fp_out,"\n. Starting tree number: \t\t\t#%d",num_tree+1);
2558
if(io->mod->s_opt->opt_topo)
2560
if(io->mod->s_opt->topo_search == NNI_MOVE) PhyML_Fprintf(fp_out,"\n. Tree topology search : \t\tNNIs");
2561
else if(io->mod->s_opt->topo_search == SPR_MOVE) PhyML_Fprintf(fp_out,"\n. Tree topology search : \t\tSPRs");
2562
else if(io->mod->s_opt->topo_search == BEST_OF_NNI_AND_SPR) PhyML_Fprintf(fp_out,"\n. Tree topology search : \t\tBest of NNIs and SPRs");
2566
PhyML_Fprintf(fp_out,"\n. Tree topology: \t\t\tfixed");
2570
/* was after Sequence file ; moved here FLT */
2571
s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2572
if(io->in_tree == 2)
2574
strcat(strcat(strcat(s,"user tree ("),io->in_tree_file),")");
2578
if(!io->mod->s_opt->random_input_tree)
2580
if(io->in_tree == 0)
2582
if(io->in_tree == 1)
2583
strcat(s,"parsimony");
2587
strcat(s,"random tree");
2591
PhyML_Fprintf(fp_out,"\n. Initial tree: \t\t\t%s",s);
2594
if(tree->io->datatype == NT)
2596
PhyML_Fprintf(fp_out,"\n. Model of nucleotides substitution: \t%s",tree->mod->modelname->s);
2597
if(io->mod->whichmodel == CUSTOM)
2598
PhyML_Fprintf(fp_out," (%s)",io->mod->custom_mod_string);
2600
else if(tree->io->datatype == AA)
2602
PhyML_Fprintf(fp_out,"\n. Model of amino acids substitution: \t%s",tree->mod->modelname->s);
2603
if(io->mod->whichmodel == CUSTOMAA) PhyML_Fprintf(fp_out," (%s)",tree->mod->aa_rate_mat_file->s);
2607
fprintf(fp_out,"\n. Substitution model: \t\t\t%s",tree->mod->modelname->s);
2611
PhyML_Fprintf(fp_out,"\n. Number of taxa: \t\t\t%d",tree->n_otu);/*added FLT*/
2613
PhyML_Fprintf(fp_out,"\n. Log-likelihood: \t\t\t%.5f",tree->c_lnL);/*was last ; moved here FLT*/
2615
Unconstraint_Lk(tree);
2616
PhyML_Fprintf(fp_out,"\n. Unconstrained likelihood: \t\t%.5f",tree->unconstraint_lk);
2618
PhyML_Fprintf(fp_out,"\n. Parsimony: \t\t\t\t%d",tree->c_pars);
2620
PhyML_Fprintf(fp_out,"\n. Tree size: \t\t\t\t%.5f",Get_Tree_Size(tree));
2622
/* if(tree->mod->ras->n_catg > 1 && tree->mod->ras->free_mixt_rates == NO) */
2623
if(tree->mod->ras->free_mixt_rates == NO)
2625
PhyML_Fprintf(fp_out,"\n. Discrete gamma model: \t\t%s","Yes");
2626
PhyML_Fprintf(fp_out,"\n - Number of classes: \t\t\t%d",tree->mod->ras->n_catg);
2627
PhyML_Fprintf(fp_out,"\n - Gamma shape parameter: \t\t%.3f",tree->mod->ras->alpha->v);
2628
For(i,tree->mod->ras->n_catg)
2630
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]);
2633
else if(tree->mod->ras->free_mixt_rates == YES)
2635
PhyML_Fprintf(fp_out,"\n. FreeRate model: \t\t\t%s","Yes");
2636
PhyML_Fprintf(fp_out,"\n - Number of classes: \t\t\t%d",tree->mod->ras->n_catg);
2637
For(i,tree->mod->ras->n_catg)
2639
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]);
2643
if(tree->mod->ras->invar) PhyML_Fprintf(fp_out,"\n. Proportion of invariant: \t\t%.3f",tree->mod->ras->pinvar->v);
2645
if(tree->mod->gamma_mgf_bl == YES) PhyML_Fprintf(fp_out,"\n. Variance of branch lengths: \t\t%f",tree->mod->l_var_sigma);
2647
/*was before Discrete gamma model ; moved here FLT*/
2648
if((tree->mod->whichmodel == K80) ||
2649
(tree->mod->whichmodel == HKY85) ||
2650
(tree->mod->whichmodel == F84))
2651
PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio: \t%.3f",tree->mod->kappa->v);
2652
else if(tree->mod->whichmodel == TN93)
2654
PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio for purines: \t\t%.3f",
2655
tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
2656
PhyML_Fprintf(fp_out,"\n. Transition/transversion ratio for pyrimidines: \t%.3f",
2657
tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
2660
if(tree->io->datatype == NT)
2662
PhyML_Fprintf(fp_out,"\n. Nucleotides frequencies:");
2663
PhyML_Fprintf(fp_out,"\n - f(A)=%8.5f",tree->mod->e_frq->pi->v[0]);
2664
PhyML_Fprintf(fp_out,"\n - f(C)=%8.5f",tree->mod->e_frq->pi->v[1]);
2665
PhyML_Fprintf(fp_out,"\n - f(G)=%8.5f",tree->mod->e_frq->pi->v[2]);
2666
PhyML_Fprintf(fp_out,"\n - f(T)=%8.5f",tree->mod->e_frq->pi->v[3]);
2669
/*****************************************/
2670
if((tree->mod->whichmodel == GTR) ||
2671
(tree->mod->whichmodel == CUSTOM))
2675
Update_Qmat_GTR(tree->mod->r_mat->rr->v,
2676
tree->mod->r_mat->rr_val->v,
2677
tree->mod->r_mat->rr_num->v,
2678
tree->mod->e_frq->pi->v,
2679
tree->mod->r_mat->qmat->v);
2681
PhyML_Fprintf(fp_out,"\n");
2682
PhyML_Fprintf(fp_out,". GTR relative rate parameters : \n");
2683
PhyML_Fprintf(fp_out," A <-> C %8.5f\n", tree->mod->r_mat->rr->v[0]);
2684
PhyML_Fprintf(fp_out," A <-> G %8.5f\n", tree->mod->r_mat->rr->v[1]);
2685
PhyML_Fprintf(fp_out," A <-> T %8.5f\n", tree->mod->r_mat->rr->v[2]);
2686
PhyML_Fprintf(fp_out," C <-> G %8.5f\n", tree->mod->r_mat->rr->v[3]);
2687
PhyML_Fprintf(fp_out," C <-> T %8.5f\n", tree->mod->r_mat->rr->v[4]);
2688
PhyML_Fprintf(fp_out," G <-> T %8.5f\n",tree->mod->r_mat->rr->v[5]);
2691
PhyML_Fprintf(fp_out,"\n. Instantaneous rate matrix : ");
2692
PhyML_Fprintf(fp_out,"\n [A---------C---------G---------T------]\n");
2695
PhyML_Fprintf(fp_out," ");
2697
PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->r_mat->qmat->v[i*4+j]);
2698
PhyML_Fprintf(fp_out,"\n");
2700
PhyML_Fprintf(fp_out,"\n");
2702
/*****************************************/
2705
if(io->ratio_test == 1)
2707
PhyML_Fprintf(fp_out,". aLRT statistics to test branches");
2709
else if(io->ratio_test == 2)
2711
PhyML_Fprintf(fp_out,". aLRT branch supports (cubic approximation, mixture of Chi2s distribution)");
2715
PhyML_Fprintf(fp_out,"\n");
2716
PhyML_Fprintf(fp_out,"\n. Run ID:\t\t\t\t%s", (io->append_run_ID) ? (io->run_id_string): ("none"));
2717
PhyML_Fprintf(fp_out,"\n. Random seed:\t\t\t\t%d", io->r_seed);
2718
PhyML_Fprintf(fp_out,"\n. Subtree patterns aliasing:\t\t%s",io->do_alias_subpatt?"yes":"no");
2719
PhyML_Fprintf(fp_out,"\n. Version:\t\t\t\t%s", VERSION);
2721
hour = div(t_end-t_beg,3600);
2722
min = div(t_end-t_beg,60 );
2724
min.quot -= hour.quot*60;
2726
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));
2729
if(add_citation == YES)
2731
PhyML_Fprintf(fp_out,"\n\n");
2732
PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2733
PhyML_Fprintf(fp_out," Suggested citations:\n");
2734
PhyML_Fprintf(fp_out," S. Guindon, JF. Dufayard, V. Lefort, M. Anisimova, W. Hordijk, O. Gascuel\n");
2735
PhyML_Fprintf(fp_out," \"New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0.\"\n");
2736
PhyML_Fprintf(fp_out," Systematic Biology. 2010. 59(3):307-321.\n");
2737
PhyML_Fprintf(fp_out,"\n");
2738
PhyML_Fprintf(fp_out," S. Guindon & O. Gascuel\n");
2739
PhyML_Fprintf(fp_out," \"A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood\"\n");
2740
PhyML_Fprintf(fp_out," Systematic Biology. 2003. 52(5):696-704.\n");
2741
PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2745
PhyML_Fprintf(fp_out,"\n\n");
2746
PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2747
PhyML_Fprintf(fp_out,"\n");
2752
//////////////////////////////////////////////////////////////
2753
//////////////////////////////////////////////////////////////
2755
/*FLT wrote this function*/
2756
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)
2764
PhyML_Fprintf(fp_out,". Sequence file : [%s]\n\n", Basename(io->in_align_file));
2766
if((tree->io->datatype == NT) || (tree->io->datatype == AA))
2768
(tree->io->datatype == NT)?
2769
(PhyML_Fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",io->mod->modelname->s)):
2770
(PhyML_Fprintf(fp_out,". Model of amino acids substitution : %s\n\n",io->mod->modelname->s));
2773
s = (char *)mCalloc(100,sizeof(char));
2777
case 0: { strcpy(s,"BioNJ"); break; }
2778
case 1: { strcpy(s,"parsimony"); break; }
2779
case 2: { strcpy(s,"user tree (");
2780
strcat(s,io->in_tree_file);
2781
strcat(s,")"); break; }
2784
PhyML_Fprintf(fp_out,". Initial tree : [%s]\n\n",s);
2788
PhyML_Fprintf(fp_out,"\n");
2791
PhyML_Fprintf(fp_out, ". Data\t");
2793
PhyML_Fprintf(fp_out,"Nb of \t");
2795
PhyML_Fprintf(fp_out,"Likelihood\t");
2797
PhyML_Fprintf(fp_out, "Discrete \t");
2799
if(tree->mod->ras->n_catg > 1)
2800
PhyML_Fprintf(fp_out, "Number of \tGamma shape\t");
2802
PhyML_Fprintf(fp_out,"Proportion of\t");
2804
if(tree->mod->whichmodel <= 6)
2805
PhyML_Fprintf(fp_out,"Transition/ \t");
2807
PhyML_Fprintf(fp_out,"Nucleotides frequencies \t");
2809
if((tree->mod->whichmodel == GTR) ||
2810
(tree->mod->whichmodel == CUSTOM))
2811
PhyML_Fprintf(fp_out,"Instantaneous rate matrix \t");
2813
/* PhyML_Fprintf(fp_out,"Time\t");*/
2815
PhyML_Fprintf(fp_out, "\n");
2819
PhyML_Fprintf(fp_out, " set\t");
2821
PhyML_Fprintf(fp_out,"taxa\t");
2823
PhyML_Fprintf(fp_out,"loglk \t");
2825
PhyML_Fprintf(fp_out, "gamma model\t");
2827
if(tree->mod->ras->n_catg > 1)
2828
PhyML_Fprintf(fp_out, "categories\tparameter \t");
2830
PhyML_Fprintf(fp_out,"invariant \t");
2832
if(tree->mod->whichmodel <= 6)
2833
PhyML_Fprintf(fp_out,"transversion\t");
2835
PhyML_Fprintf(fp_out,"f(A) f(C) f(G) f(T) \t");
2837
if((tree->mod->whichmodel == GTR) ||
2838
(tree->mod->whichmodel == CUSTOM))
2839
PhyML_Fprintf(fp_out,"[A---------C---------G---------T------]\t");
2841
/* PhyML_PhyML_Fprintf(fp_out,"used\t");*/
2843
PhyML_Fprintf(fp_out, "\n");
2847
if(tree->mod->whichmodel == TN93)
2849
PhyML_Fprintf(fp_out," \t \t \t \t");
2850
if(tree->mod->ras->n_catg > 1) PhyML_Fprintf(fp_out," \t \t");
2851
PhyML_Fprintf(fp_out," \t");
2852
PhyML_Fprintf(fp_out,"purines pyrimid.\t");
2853
PhyML_Fprintf(fp_out, "\n");
2856
PhyML_Fprintf(fp_out, "\n");
2862
PhyML_Fprintf(fp_out," #%d\t",n_data_set);
2864
PhyML_Fprintf(fp_out,"%d \t",tree->n_otu);
2866
PhyML_Fprintf(fp_out,"%.5f\t",tree->c_lnL);
2868
PhyML_Fprintf(fp_out,"%s \t",
2869
(tree->mod->ras->n_catg>1)?("Yes"):("No "));
2870
if(tree->mod->ras->n_catg > 1)
2872
PhyML_Fprintf(fp_out,"%d \t",tree->mod->ras->n_catg);
2873
PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->ras->alpha->v);
2876
/*if(tree->mod->ras->invar)*/
2877
PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->ras->pinvar->v);
2879
if(tree->mod->whichmodel <= 5)
2881
PhyML_Fprintf(fp_out,"%.3f \t",tree->mod->kappa->v);
2883
else if(tree->mod->whichmodel == TN93)
2885
PhyML_Fprintf(fp_out,"%.3f ",
2886
tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
2887
PhyML_Fprintf(fp_out,"%.3f\t",
2888
tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
2892
if(tree->io->datatype == NT)
2894
PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[0]);
2895
PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[1]);
2896
PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->e_frq->pi->v[2]);
2897
PhyML_Fprintf(fp_out,"%8.5f\t",tree->mod->e_frq->pi->v[3]);
2900
hour = div(t_end-t_beg,3600);
2901
min = div(t_end-t_beg,60 );
2903
min.quot -= hour.quot*60;
2905
PhyML_Fprintf(fp_out,"%dh%dm%ds\t", hour.quot,min.quot,(int)(t_end-t_beg)%60);
2906
if(t_end-t_beg > 60)
2907
PhyML_Fprintf(fp_out,". -> %d seconds\t",(int)(t_end-t_beg));
2910
/*****************************************/
2911
if((tree->mod->whichmodel == GTR) || (tree->mod->whichmodel == CUSTOM))
2919
PhyML_Fprintf(fp_out," \t \t \t \t");
2920
if(tree->mod->ras->n_catg > 1) PhyML_Fprintf(fp_out," \t \t");
2921
PhyML_Fprintf(fp_out," \t \t");
2924
PhyML_Fprintf(fp_out,"%8.5f ",tree->mod->r_mat->qmat->v[i*4+j]);
2925
if (i<3) PhyML_Fprintf(fp_out,"\n");
2928
/*****************************************/
2930
PhyML_Fprintf(fp_out, "\n\n");
2933
//////////////////////////////////////////////////////////////
2934
//////////////////////////////////////////////////////////////
2937
void Print_Freq(t_tree *tree)
2940
switch(tree->io->datatype)
2944
PhyML_Printf("A : %f\n",tree->mod->e_frq->pi->v[0]);
2945
PhyML_Printf("C : %f\n",tree->mod->e_frq->pi->v[1]);
2946
PhyML_Printf("G : %f\n",tree->mod->e_frq->pi->v[2]);
2947
PhyML_Printf("T : %f\n",tree->mod->e_frq->pi->v[3]);
2949
/* PhyML_Printf("U : %f\n",tree->mod->e_frq->pi->v[4]); */
2950
/* PhyML_Printf("M : %f\n",tree->mod->e_frq->pi->v[5]); */
2951
/* PhyML_Printf("R : %f\n",tree->mod->e_frq->pi->v[6]); */
2952
/* PhyML_Printf("W : %f\n",tree->mod->e_frq->pi->v[7]); */
2953
/* PhyML_Printf("S : %f\n",tree->mod->e_frq->pi->v[8]); */
2954
/* PhyML_Printf("Y : %f\n",tree->mod->e_frq->pi->v[9]); */
2955
/* PhyML_Printf("K : %f\n",tree->mod->e_frq->pi->v[10]); */
2956
/* PhyML_Printf("B : %f\n",tree->mod->e_frq->pi->v[11]); */
2957
/* PhyML_Printf("D : %f\n",tree->mod->e_frq->pi->v[12]); */
2958
/* PhyML_Printf("H : %f\n",tree->mod->e_frq->pi->v[13]); */
2959
/* PhyML_Printf("V : %f\n",tree->mod->e_frq->pi->v[14]); */
2960
/* PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[15]); */
2965
PhyML_Printf("A : %f\n",tree->mod->e_frq->pi->v[0]);
2966
PhyML_Printf("R : %f\n",tree->mod->e_frq->pi->v[1]);
2967
PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[2]);
2968
PhyML_Printf("D : %f\n",tree->mod->e_frq->pi->v[3]);
2969
PhyML_Printf("C : %f\n",tree->mod->e_frq->pi->v[4]);
2970
PhyML_Printf("Q : %f\n",tree->mod->e_frq->pi->v[5]);
2971
PhyML_Printf("E : %f\n",tree->mod->e_frq->pi->v[6]);
2972
PhyML_Printf("G : %f\n",tree->mod->e_frq->pi->v[7]);
2973
PhyML_Printf("H : %f\n",tree->mod->e_frq->pi->v[8]);
2974
PhyML_Printf("I : %f\n",tree->mod->e_frq->pi->v[9]);
2975
PhyML_Printf("L : %f\n",tree->mod->e_frq->pi->v[10]);
2976
PhyML_Printf("K : %f\n",tree->mod->e_frq->pi->v[11]);
2977
PhyML_Printf("M : %f\n",tree->mod->e_frq->pi->v[12]);
2978
PhyML_Printf("F : %f\n",tree->mod->e_frq->pi->v[13]);
2979
PhyML_Printf("P : %f\n",tree->mod->e_frq->pi->v[14]);
2980
PhyML_Printf("S : %f\n",tree->mod->e_frq->pi->v[15]);
2981
PhyML_Printf("T : %f\n",tree->mod->e_frq->pi->v[16]);
2982
PhyML_Printf("W : %f\n",tree->mod->e_frq->pi->v[17]);
2983
PhyML_Printf("Y : %f\n",tree->mod->e_frq->pi->v[18]);
2984
PhyML_Printf("V : %f\n",tree->mod->e_frq->pi->v[19]);
2986
PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[20]);
2993
//////////////////////////////////////////////////////////////
2994
//////////////////////////////////////////////////////////////
2997
void Print_Settings(option *io)
3002
s = (char *)mCalloc(100,sizeof(char));
3004
PhyML_Printf("\n\n\n");
3005
PhyML_Printf("\n\n");
3007
PhyML_Printf(" .......................... \n");
3008
PhyML_Printf(" ooooooooooooooooooooooooooooo CURRENT SETTINGS ooooooooooooooooooooooooooooooooooo\n");
3009
PhyML_Printf(" .......................... \n");
3011
PhyML_Printf("\n . Sequence filename:\t\t\t\t %s", Basename(io->in_align_file));
3013
if(io->datatype == NT) strcpy(s,"dna");
3014
else if(io->datatype == AA) strcpy(s,"aa");
3015
else strcpy(s,"generic");
3017
PhyML_Printf("\n . Data type:\t\t\t\t\t %s",s);
3018
PhyML_Printf("\n . Alphabet size:\t\t\t\t %d",io->mod->ns);
3020
PhyML_Printf("\n . Sequence format:\t\t\t\t %s", io->interleaved ? "interleaved": "sequential");
3021
PhyML_Printf("\n . Number of data sets:\t\t\t\t %d", io->n_data_sets);
3023
PhyML_Printf("\n . Nb of bootstrapped data sets:\t\t\t %d", io->mod->bootstrap);
3025
if (io->mod->bootstrap > 0)
3026
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t no");
3029
if(io->ratio_test == 1)
3030
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (aLRT statistics)");
3031
else if(io->ratio_test == 2)
3032
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (Chi2-based parametric branch supports)");
3033
else if(io->ratio_test == 3)
3034
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (Minimum of SH-like and Chi2-based branch supports)");
3035
else if(io->ratio_test == 4)
3036
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (SH-like branch supports)");
3037
else if(io->ratio_test == 5)
3038
PhyML_Printf("\n . Compute approximate likelihood ratio test:\t yes (aBayes branch supports)");
3041
PhyML_Printf("\n . Model name:\t\t\t\t\t %s", io->mod->modelname->s);
3043
if(io->datatype == AA && io->mod->whichmodel == CUSTOMAA) PhyML_Printf(" (%s)",io->mod->aa_rate_mat_file->s);
3045
if (io->datatype == NT)
3047
if ((io->mod->whichmodel == K80) ||
3048
(io->mod->whichmodel == HKY85)||
3049
(io->mod->whichmodel == F84) ||
3050
(io->mod->whichmodel == TN93))
3052
if (io->mod->s_opt->opt_kappa)
3053
PhyML_Printf("\n . Ts/tv ratio:\t\t\t\t\t estimated");
3055
PhyML_Printf("\n . Ts/tv ratio:\t\t\t\t\t %f", io->mod->kappa->v);
3059
if (io->mod->s_opt->opt_pinvar)
3060
PhyML_Printf("\n . Proportion of invariable sites:\t\t estimated");
3062
PhyML_Printf("\n . Proportion of invariable sites:\t\t %f", io->mod->ras->pinvar->v);
3065
PhyML_Printf("\n . Number of subst. rate categs:\t\t\t %d", io->mod->ras->n_catg);
3066
if(io->mod->ras->n_catg > 1)
3068
if(io->mod->ras->free_mixt_rates == NO)
3070
if(io->mod->s_opt->opt_alpha)
3071
PhyML_Printf("\n . Gamma distribution parameter:\t\t\t estimated");
3073
PhyML_Printf("\n . Gamma distribution parameter:\t\t\t %f", io->mod->ras->alpha->v);
3074
PhyML_Printf("\n . 'Middle' of each rate class:\t\t\t %s",(io->mod->ras->gamma_median)?("median"):("mean"));
3079
if(io->datatype == AA)
3080
PhyML_Printf("\n . Amino acid equilibrium frequencies:\t\t %s", (io->mod->s_opt->opt_state_freq) ? ("empirical"):("model"));
3081
else if(io->datatype == NT)
3083
if((io->mod->whichmodel != JC69) &&
3084
(io->mod->whichmodel != K80) &&
3085
(io->mod->whichmodel != F81))
3087
if(!io->mod->s_opt->user_state_freq)
3089
PhyML_Printf("\n . Nucleotide equilibrium frequencies:\t\t %s", (io->mod->s_opt->opt_state_freq) ? ("ML"):("empirical"));
3093
PhyML_Printf("\n . Nucleotide equilibrium frequencies:\t\t %s","user-defined");
3098
PhyML_Printf("\n . Optimise tree topology:\t\t\t %s", (io->mod->s_opt->opt_topo) ? "yes": "no");
3102
case 0: { strcpy(s,"BioNJ"); break; }
3103
case 1: { strcpy(s,"parsimony"); break; }
3104
case 2: { strcpy(s,"user tree (");
3105
strcat(s,Basename(io->in_tree_file));
3106
strcat(s,")"); break; }
3109
if(io->mod->s_opt->opt_topo)
3111
if(io->mod->s_opt->topo_search == NNI_MOVE) PhyML_Printf("\n . Tree topology search:\t\t\t\t NNIs");
3112
else if(io->mod->s_opt->topo_search == SPR_MOVE) PhyML_Printf("\n . Tree topology search:\t\t\t\t SPRs");
3113
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");
3117
PhyML_Printf("\n . Starting tree:\t\t\t\t %s",s);
3119
PhyML_Printf("\n . Add random input tree:\t\t\t %s", (io->mod->s_opt->random_input_tree) ? "yes": "no");
3120
if(io->mod->s_opt->random_input_tree)
3121
PhyML_Printf("\n . Number of random starting trees:\t\t %d", io->mod->s_opt->n_rand_starts);
3124
if(!io->mod->s_opt->random_input_tree)
3125
PhyML_Printf("\n . Evaluated tree:\t\t\t\t \"%s\"",s);
3127
PhyML_Printf("\n . Optimise branch lengths:\t\t\t %s", (io->mod->s_opt->opt_bl) ? "yes": "no");
3130
if(io->mod->s_opt->opt_alpha ||
3131
io->mod->s_opt->opt_kappa ||
3132
io->mod->s_opt->opt_lambda ||
3133
io->mod->s_opt->opt_pinvar ||
3134
io->mod->s_opt->opt_rr) answer = 1;
3136
PhyML_Printf("\n . Optimise substitution model parameters:\t %s", (answer) ? "yes": "no");
3138
PhyML_Printf("\n . Run ID:\t\t\t\t\t %s", (io->append_run_ID) ? (io->run_id_string): ("none"));
3139
PhyML_Printf("\n . Random seed:\t\t\t\t\t %d", io->r_seed);
3140
PhyML_Printf("\n . Subtree patterns aliasing:\t\t\t %s",io->do_alias_subpatt?"yes":"no");
3141
PhyML_Printf("\n . Version:\t\t\t\t\t %s", VERSION);
3144
PhyML_Printf("\n\n oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3146
PhyML_Printf("\n\n");
3152
void Print_Banner(FILE *fp)
3154
PhyML_Fprintf(fp,"\n");
3155
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3156
PhyML_Fprintf(fp," \n");
3157
PhyML_Fprintf(fp," --- PhyML %s --- \n",VERSION);
3158
PhyML_Fprintf(fp," \n");
3159
PhyML_Fprintf(fp," A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood \n");
3160
PhyML_Fprintf(fp," Stephane Guindon & Olivier Gascuel \n");
3161
PhyML_Fprintf(fp," \n");
3162
PhyML_Fprintf(fp," http://www.atgc-montpellier.fr/phyml \n");
3163
PhyML_Fprintf(fp," \n");
3164
PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier II \n");
3165
PhyML_Fprintf(fp," \n");
3166
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3169
//////////////////////////////////////////////////////////////
3170
//////////////////////////////////////////////////////////////
3173
void Print_Banner_Small(FILE *fp)
3175
PhyML_Fprintf(fp,"\n");
3176
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3177
PhyML_Fprintf(fp," --- PhyML %s --- \n",VERSION);
3178
PhyML_Fprintf(fp," http://www.atgc-montpellier.fr/phyml \n");
3179
PhyML_Fprintf(fp," Copyright CNRS - Universite Montpellier II \n");
3180
PhyML_Fprintf(fp," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3183
//////////////////////////////////////////////////////////////
3184
//////////////////////////////////////////////////////////////
3188
void Print_Data_Set_Number(option *io, FILE *fp)
3190
PhyML_Fprintf(fp,"\n");
3191
PhyML_Fprintf(fp," \n");
3192
PhyML_Fprintf(fp," [ Data set number %3d ] \n",io->curr_gt+1);
3193
PhyML_Fprintf(fp," \n");
3195
//////////////////////////////////////////////////////////////
3196
//////////////////////////////////////////////////////////////
3198
void Print_Lk(t_tree *tree, char *string)
3200
time(&(tree->t_current));
3201
PhyML_Printf("\n. (%5d sec) [%15.4f] %s",
3202
(int)(tree->t_current-tree->t_beg),tree->c_lnL,
3209
//////////////////////////////////////////////////////////////
3210
//////////////////////////////////////////////////////////////
3213
void Print_Pars(t_tree *tree)
3215
time(&(tree->t_current));
3216
PhyML_Printf("\n. (%5d sec) [%5d]",(int)(tree->t_current-tree->t_beg),tree->c_pars);
3222
//////////////////////////////////////////////////////////////
3223
//////////////////////////////////////////////////////////////
3225
void Print_Lk_And_Pars(t_tree *tree)
3227
time(&(tree->t_current));
3229
PhyML_Printf("\n. (%5d sec) [%15.4f] [%5d]",
3230
(int)(tree->t_current-tree->t_beg),
3231
tree->c_lnL,tree->c_pars);
3238
//////////////////////////////////////////////////////////////
3239
//////////////////////////////////////////////////////////////
3241
void Read_Qmat(phydbl *daa, phydbl *pi, FILE *fp)
3253
/* if(!fscanf(fp,"%lf",&(daa[i*20+j]))) Exit("\n"); */
3254
if(!fscanf(fp,"%lf",&val))
3256
PhyML_Printf("\n== Rate matrix file does not appear to have a proper format. Please refer to the documentation.");
3259
daa[i*20+j] = (phydbl)val;
3260
daa[j*20+i] = daa[i*20+j];
3268
if(!fscanf(fp,"%lf",&val)) Exit("\n");
3269
pi[i] = (phydbl)val;
3272
For(i,20) sum += pi[i];
3273
if(FABS(sum - 1.) > 1.E-06)
3275
PhyML_Printf("\n. Sum of amino-acid frequencies: %f",sum);
3276
PhyML_Printf("\n. Scaling amino-acid frequencies...\n");
3277
For(i,20) pi[i] /= sum;
3281
//////////////////////////////////////////////////////////////
3282
//////////////////////////////////////////////////////////////
3285
void Print_Qmat_AA(phydbl *daa, phydbl *pi)
3294
PhyML_Printf("daa[%2d*20+%2d] = %10f; ",i,j,daa[i*20+j]);
3296
if(!(cpt%4)) PhyML_Printf("\n");
3300
PhyML_Printf("\n\n");
3301
PhyML_Printf("for (i=0; i<naa; i++) for (j=0; j<i; j++) daa[j*naa+i] = daa[i*naa+j];\n\n");
3302
For(i,20) PhyML_Printf("pi[%d] = %f; ",i,pi[i]);
3304
PhyML_Printf("Ala\tArg\tAsn\tAsp\tCys\tGln\tGlu\tGly\tHis\tIle\tLeu\tLys\tMet\tPhe\tPro\tSer\tThr\tTrp\tTyr\tVal\n");
3308
//////////////////////////////////////////////////////////////
3309
//////////////////////////////////////////////////////////////
3311
void Print_Square_Matrix_Generic(int n, phydbl *mat)
3318
PhyML_Printf("[%3d]",i);
3321
PhyML_Printf("%12.5f ",mat[i*n+j]);
3328
//////////////////////////////////////////////////////////////
3329
//////////////////////////////////////////////////////////////
3332
void Print_Diversity(FILE *fp, t_tree *tree)
3335
Print_Diversity_Pre(tree->a_nodes[0],
3336
tree->a_nodes[0]->v[0],
3337
tree->a_nodes[0]->b[0],
3341
/* mean_div_left = .0; */
3344
/* mean_div_left += (k+1) * tree->a_edges[j]->div_post_pred_left[k]; */
3346
/* mean_div_rght = .0; */
3347
/* For(k,ns) mean_div_rght += (k+1) * tree->a_edges[j]->div_post_pred_rght[k]; */
3349
/* mean_div_left /= (phydbl)tree->data->init_len; */
3350
/* mean_div_rght /= (phydbl)tree->data->init_len; */
3352
/* PhyML_Fprintf(fp,"%4d 0 %f\n",j,mean_div_left); */
3353
/* PhyML_Fprintf(fp,"%4d 1 %f\n",j,mean_div_rght); */
3356
/* mean_div_left = .0; */
3357
/* For(k,ns) mean_div_left += tree->a_edges[j]->div_post_pred_left[k]; */
3359
/* mean_div_rght = .0; */
3362
/* mean_div_rght += tree->a_edges[j]->div_post_pred_rght[k]; */
3365
/* if((mean_div_left != tree->data->init_len) || (mean_div_rght != tree->data->init_len)) */
3367
/* PhyML_Printf("\n. mean_div_left = %f mean_div_rght = %f init_len = %d", */
3368
/* mean_div_left,mean_div_rght,tree->data->init_len); */
3369
/* PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
3370
/* Warn_And_Exit(""); */
3374
//////////////////////////////////////////////////////////////
3375
//////////////////////////////////////////////////////////////
3378
void Print_Diversity_Pre(t_node *a, t_node *d, t_edge *b, FILE *fp, t_tree *tree)
3388
if(tree->io->datatype == NT) ns = 4;
3389
else if(tree->io->datatype == AA) ns = 20;
3391
if(d == b->left) For(k,ns) PhyML_Fprintf(fp,"%4d 0 %2d %4d\n",b->num,k,b->div_post_pred_left[k]);
3392
else For(k,ns) PhyML_Fprintf(fp,"%4d 1 %2d %4d\n",b->num,k,b->div_post_pred_rght[k]);
3394
For(k,3) if(d->v[k] != a) Print_Diversity_Pre(d,d->v[k],d->b[k],fp,tree);
3399
//////////////////////////////////////////////////////////////
3400
//////////////////////////////////////////////////////////////
3402
t_tree *Read_User_Tree(calign *cdata, t_mod *mod, option *io)
3407
PhyML_Printf("\n. Reading tree..."); fflush(NULL);
3408
if(io->n_trees == 1) rewind(io->fp_in_tree);
3409
tree = Read_Tree_File_Phylip(io->fp_in_tree);
3410
if(!tree) Exit("\n. Input tree not found...");
3411
/* Add branch lengths if necessary */
3412
if(!tree->has_branch_lengths) Add_BioNJ_Branch_Lengths(tree,cdata,mod);
3416
//////////////////////////////////////////////////////////////
3417
//////////////////////////////////////////////////////////////
3420
void Print_Time_Info(time_t t_beg, time_t t_end)
3424
hour = div(t_end-t_beg,3600);
3425
min = div(t_end-t_beg,60 );
3426
min.quot -= hour.quot*60;
3428
PhyML_Printf("\n\n. Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
3429
PhyML_Printf("\noooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3432
//////////////////////////////////////////////////////////////
3433
//////////////////////////////////////////////////////////////
3435
void PhyML_Printf(char *format, ...)
3440
if(Global_myRank == 0)
3442
va_start (ptr, format);
3443
vprintf (format, ptr);
3447
va_start (ptr, format);
3448
vprintf (format, ptr);
3455
//////////////////////////////////////////////////////////////
3456
//////////////////////////////////////////////////////////////
3459
void PhyML_Fprintf(FILE *fp, char *format, ...)
3464
if(Global_myRank == 0)
3466
va_start (ptr, format);
3467
vfprintf (fp,format, ptr);
3471
va_start (ptr, format);
3472
vfprintf (fp,format, ptr);
3479
//////////////////////////////////////////////////////////////
3480
//////////////////////////////////////////////////////////////
3482
//////////////////////////////////////////////////////////////
3483
//////////////////////////////////////////////////////////////
3485
void Read_Clade_Priors(char *file_name, t_tree *tree)
3493
phydbl prior_low,prior_up;
3497
PhyML_Printf("\n. Reading prior on node ages.\n");
3499
line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3500
s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3502
clade_list = (char **)mCalloc(tree->n_otu,sizeof(char *));
3503
For(i,tree->n_otu) clade_list[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
3505
fp = Openfile(file_name,0);
3510
if(!fgets(line,T_MAX_LINE,fp)) break;
3518
while(line[pos] == ' ') pos++;
3520
while((line[pos] != ' ') && (line[pos] != '\n') && line[pos] != '#')
3528
/* PhyML_Printf("\n. s = %s\n",s); */
3530
if(line[pos] == '\n' || line[pos] == '#') break;
3535
strcpy(clade_list[clade_size],s);
3544
if(line[pos] != '#' && line[pos] != '\n')
3547
/* sscanf(line+pos,"%lf %lf",&prior_up,&prior_low); */
3548
sscanf(line+pos,"%lf %lf",&val1,&val2);
3549
prior_up = (phydbl)val1;
3550
prior_low = (phydbl)val2;
3552
if(!strcmp("@root@",clade_list[0])) node_num = tree->n_root->num;
3553
else node_num = Find_Clade(clade_list, clade_size, tree);
3561
PhyML_Printf("\n. .................................................................");
3562
PhyML_Printf("\n. WARNING: could not find any clade in the tree referred to with the following taxon names:");
3563
For(i,clade_size) PhyML_Printf("\n. \"%s\"",clade_list[i]);
3564
PhyML_Printf("\n. .................................................................");
3569
tree->rates->t_has_prior[node_num] = YES;
3570
tree->rates->t_prior_min[node_num] = MIN(prior_low,prior_up);
3571
tree->rates->t_prior_max[node_num] = MAX(prior_low,prior_up);
3574
/* Sergei to add stuff here re calibration object */
3577
if(FABS(prior_low - prior_up) < 1.E-6 && tree->a_nodes[node_num]->tax == YES)
3578
tree->rates->nd_t[node_num] = prior_low;
3581
PhyML_Printf("\n. [%3d]..................................................................",n_clade_priors);
3582
PhyML_Printf("\n. Node %4d matches the clade referred to with the following taxon names:",node_num);
3583
For(i,clade_size) PhyML_Printf("\n. - \"%s\"",clade_list[i]);
3584
PhyML_Printf("\n. Lower bound set to: %15f time units.",MIN(prior_low,prior_up));
3585
PhyML_Printf("\n. Upper bound set to: %15f time units.",MAX(prior_low,prior_up));
3586
PhyML_Printf("\n. .......................................................................");
3593
PhyML_Printf("\n. Read prior information on %d %s.",n_clade_priors,n_clade_priors > 1 ? "clades":"clade");
3597
PhyML_Printf("\n. PhyTime could not find any prior on node age.");
3598
PhyML_Printf("\n. This is likely due to a problem in the calibration ");
3599
PhyML_Printf("\n. file format. Make sure, for instance, that there is ");
3600
PhyML_Printf("\n. a blank character between the end of the last name");
3601
PhyML_Printf("\n. of each clade and the character `|'. Otherwise, ");
3602
PhyML_Printf("\n. please refer to the example file.\n");
3603
PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3607
For(i,tree->n_otu) Free(clade_list[i]);
3614
//////////////////////////////////////////////////////////////
3615
//////////////////////////////////////////////////////////////
3617
option *Get_Input(int argc, char **argv)
3627
io = (option *)Make_Input();
3628
mod = (t_mod *)Make_Model_Basic();
3629
s_opt = (t_opt *)Make_Optimiz();
3630
m4mod = (m4 *)M4_Make_Light();
3633
Set_Defaults_Input(io);
3634
Set_Defaults_Model(mod);
3635
Set_Defaults_Optimiz(s_opt);
3637
io->mod->m4mod = m4mod;
3645
rv = Read_Command_Line(io,argc,argv);
3646
#elif (defined PHYTIME || defined SERGEII)
3647
rv = Read_Command_Line(io,argc,argv);
3655
Launch_Interface(io);
3664
rv = Read_Command_Line(io,argc,argv);
3672
//////////////////////////////////////////////////////////////
3673
//////////////////////////////////////////////////////////////
3675
int Set_Whichmodel(int select)
3800
PhyML_Printf("\n== Model number %d is unknown. Please use a valid model name",select);
3810
//////////////////////////////////////////////////////////////
3811
//////////////////////////////////////////////////////////////
3813
void Print_Data_Structure(int final, FILE *fp, t_tree *mixt_tree)
3815
int n_partition_elem;
3817
t_tree *tree,*cpy_mixt_tree;
3818
int c,cc,cc_efrq,cc_rmat,cc_lens;
3820
int *link_efrq,*link_rmat,*link_lens;
3821
phydbl r_mat_weight_sum, e_frq_weight_sum;
3824
PhyML_Fprintf(fp,"\n. Starting tree: %s",
3825
mixt_tree->io->in_tree == 2?mixt_tree->io->in_tree_file:"BioNJ");
3827
PhyML_Fprintf(fp,"\n. Tree topology search: %s",
3828
mixt_tree->io->mod->s_opt->opt_topo==YES?
3829
mixt_tree->io->mod->s_opt->topo_search==SPR_MOVE?"spr":
3830
mixt_tree->io->mod->s_opt->topo_search==NNI_MOVE?"nni":
3833
cpy_mixt_tree = mixt_tree;
3835
n_partition_elem = 1;
3839
tree = tree->next_mixt;
3845
s = (char *)mCalloc(2,sizeof(char));
3851
s = (char *)mRealloc(s,(int)(strlen(s)+strlen(tree->io->in_align_file)+2+2),sizeof(char));
3852
strcat(s,tree->io->in_align_file);
3854
tree = tree->next_mixt;
3858
s[(int)strlen(s)-2]=' ';
3859
s[(int)strlen(s)-1]='\0';
3861
if(final == NO) PhyML_Fprintf(fp,"\n. Processing %d data %s (%s)",n_partition_elem,n_partition_elem>1?"sets":"set",s);
3862
if(final == YES) PhyML_Fprintf(fp,"\n. Processed %d data %s (%s)",n_partition_elem,n_partition_elem>1?"sets":"set",s);
3866
PhyML_Fprintf(fp,"\n. Final log-likelihood: %f",mixt_tree->c_lnL);
3872
PhyML_Fprintf(fp,"\n\n");
3873
PhyML_Fprintf(fp,"\n _______________________________________________________________________ ");
3874
PhyML_Fprintf(fp,"\n| |");
3875
PhyML_Fprintf(fp,"\n| %40s (partition element %2d) |",mixt_tree->io->in_align_file,mixt_tree->dp);
3876
PhyML_Fprintf(fp,"\n|_______________________________________________________________________|");
3877
PhyML_Fprintf(fp,"\n");
3879
PhyML_Fprintf(fp,"\n. Number of rate classes:\t\t%12d",mixt_tree->mod->ras->n_catg);
3880
if(mixt_tree->mod->ras->n_catg > 1)
3882
PhyML_Fprintf(fp,"\n. Model of rate variation:\t\t%12s",
3883
mixt_tree->mod->ras->free_mixt_rates?"FreeRates":
3884
mixt_tree->mod->ras->invar?"Gamma+Inv":"Gamma");
3885
if(mixt_tree->mod->ras->free_mixt_rates == NO)
3887
PhyML_Fprintf(fp,"\n. Gamma shape parameter value:\t\t%12.2f",mixt_tree->mod->ras->alpha->v);
3888
PhyML_Fprintf(fp,"\n Optimise: \t\t\t\t%12s",mixt_tree->mod->s_opt->opt_alpha==YES?"yes":"no");
3890
if(mixt_tree->mod->ras->invar == YES)
3892
PhyML_Fprintf(fp,"\n. Proportion of invariable sites:\t%12.2f",mixt_tree->mod->ras->pinvar->v);
3893
PhyML_Fprintf(fp,"\n Optimise: \t\t\t\t%12s",mixt_tree->mod->s_opt->opt_pinvar==YES?"yes":"no");
3897
r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->r_mat_weight);
3898
e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->e_frq_weight);
3903
if(tree->is_mixt_tree) tree = tree->next;
3905
PhyML_Fprintf(fp,"\n");
3906
PhyML_Fprintf(fp,"\n. Mixture class %d",class+1);
3908
if(mixt_tree->mod->ras->n_catg > 1)
3910
PhyML_Fprintf(fp,"\n Relative substitution rate:\t%12f",mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number]);
3911
PhyML_Fprintf(fp,"\n Relative rate freq.:\t\t%12f",mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number]);
3912
PhyML_Fprintf(fp,"\n Rate class number:\t\t%12d",tree->mod->ras->parent_class_number);
3915
PhyML_Fprintf(fp,"\n Substitution model:\t\t%12s",tree->mod->modelname->s);
3917
if(tree->mod->whichmodel == CUSTOM)
3918
PhyML_Fprintf(fp,"\n Substitution model code:\t%12s",tree->mod->custom_mod_string->s);
3920
if(tree->mod->whichmodel == CUSTOMAA)
3921
PhyML_Fprintf(fp,"\n Rate matrix file name:\t%12s",tree->mod->aa_rate_mat_file->s);
3923
if(tree->mod->whichmodel == K80 ||
3924
tree->mod->whichmodel == HKY85 ||
3925
tree->mod->whichmodel == TN93)
3927
PhyML_Fprintf(fp,"\n Value of the ts/tv raqtio:\t%12f",tree->mod->kappa->v);
3928
PhyML_Fprintf(fp,"\n Optimise ts/tv ratio:\t%12s",tree->mod->s_opt->opt_kappa?"yes":"no");
3930
else if(tree->mod->whichmodel == GTR ||
3931
tree->mod->whichmodel == CUSTOM)
3933
PhyML_Fprintf(fp,"\n Optimise subst. rates:\t%12s",tree->mod->s_opt->opt_rr?"yes":"no");
3936
PhyML_Fprintf(fp,"\n Subst. rate A<->C:\t\t%12.2f",tree->mod->r_mat->rr->v[0]);
3937
PhyML_Fprintf(fp,"\n Subst. rate A<->G:\t\t%12.2f",tree->mod->r_mat->rr->v[1]);
3938
PhyML_Fprintf(fp,"\n Subst. rate A<->T:\t\t%12.2f",tree->mod->r_mat->rr->v[2]);
3939
PhyML_Fprintf(fp,"\n Subst. rate C<->G:\t\t%12.2f",tree->mod->r_mat->rr->v[3]);
3940
PhyML_Fprintf(fp,"\n Subst. rate C<->T:\t\t%12.2f",tree->mod->r_mat->rr->v[4]);
3941
PhyML_Fprintf(fp,"\n Subst. rate G<->T:\t\t%12.2f",tree->mod->r_mat->rr->v[5]);
3946
PhyML_Fprintf(fp,"\n Rate matrix weight:\t\t%12f",tree->mod->r_mat_weight->v / r_mat_weight_sum);
3950
if(tree->io->datatype == NT &&
3951
tree->mod->whichmodel != JC69 &&
3952
tree->mod->whichmodel != K80)
3954
PhyML_Fprintf(fp,"\n Optimise nucletide freq.:\t%12s",tree->mod->s_opt->opt_state_freq?"yes":"no");
3957
PhyML_Fprintf(fp,"\n Freq(A):\t\t\t%12.2f",tree->mod->e_frq->pi->v[0]);
3958
PhyML_Fprintf(fp,"\n Freq(C):\t\t\t%12.2f",tree->mod->e_frq->pi->v[1]);
3959
PhyML_Fprintf(fp,"\n Freq(G):\t\t\t%12.2f",tree->mod->e_frq->pi->v[2]);
3960
PhyML_Fprintf(fp,"\n Freq(T):\t\t\t%12.2f",tree->mod->e_frq->pi->v[3]);
3963
else if(tree->io->datatype == AA)
3967
s = (char *)mCalloc(50,sizeof(char));
3969
if(tree->mod->s_opt->opt_state_freq == YES)
3971
strcpy(s,"Empirical");
3978
PhyML_Fprintf(fp,"\n Amino-acid freq.:\t\t%12s",s);
3984
PhyML_Fprintf(fp,"\n Equ. freq. weight:\t\t%12f",tree->mod->e_frq_weight->v / e_frq_weight_sum);
3990
if(tree && tree->is_mixt_tree == YES) break;
3994
mixt_tree = mixt_tree->next_mixt;
3995
if(!mixt_tree) break;
4000
tree = cpy_mixt_tree;
4004
if(tree->is_mixt_tree) tree = tree->next;
4010
link_efrq = (int *)mCalloc(c,sizeof(int));
4011
link_lens = (int *)mCalloc(c,sizeof(int));
4012
link_rmat = (int *)mCalloc(c,sizeof(int));
4014
PhyML_Fprintf(fp,"\n");
4015
PhyML_Fprintf(fp,"\n");
4016
PhyML_Fprintf(fp,"\n _______________________________________________________________________ ");
4017
PhyML_Fprintf(fp,"\n| |");
4018
PhyML_Fprintf(fp,"\n| Model summary table |");
4019
PhyML_Fprintf(fp,"\n|_______________________________________________________________________|");
4020
PhyML_Fprintf(fp,"\n");
4021
PhyML_Fprintf(fp,"\n");
4022
/* PhyML_Fprintf(fp," "); */
4023
PhyML_Fprintf(fp," ------------------");
4024
tree = cpy_mixt_tree;
4028
if(tree->is_mixt_tree) tree = tree->next;
4029
PhyML_Fprintf(fp,"---");
4034
param = (char *)mCalloc(30,sizeof(char));
4035
PhyML_Fprintf(fp,"\n");
4036
strcpy(param,"Partition element ");
4037
PhyML_Fprintf(fp," %18s",param);
4038
tree = cpy_mixt_tree;
4042
if(tree->is_mixt_tree) tree = tree->next;
4043
PhyML_Fprintf(fp,"%2d ",tree->mixt_tree->dp);
4049
PhyML_Fprintf(fp,"\n");
4050
PhyML_Fprintf(fp," ------------------");
4051
tree = cpy_mixt_tree;
4055
if(tree->is_mixt_tree) tree = tree->next;
4056
PhyML_Fprintf(fp,"---");
4062
tree = cpy_mixt_tree;
4066
if(tree->is_mixt_tree) tree = tree->next;
4075
mixt_tree = cpy_mixt_tree;
4082
if(mixt_tree->is_mixt_tree) mixt_tree = mixt_tree->next;
4084
if(link_efrq[cc] < 0)
4086
link_efrq[cc] = cc_efrq;
4087
tree = mixt_tree->next;
4093
if(tree->is_mixt_tree) tree = tree->next;
4095
if(mixt_tree->mod->e_frq == tree->mod->e_frq) link_efrq[c] = cc_efrq;
4105
if(link_lens[cc] < 0)
4107
link_lens[cc] = cc_lens;
4108
tree = mixt_tree->next;
4114
if(tree->is_mixt_tree) tree = tree->next;
4116
if(mixt_tree->a_edges[0]->l == tree->a_edges[0]->l) link_lens[c] = cc_lens;
4126
if(link_rmat[cc] < 0)
4128
link_rmat[cc] = cc_rmat;
4129
tree = mixt_tree->next;
4135
if(tree->is_mixt_tree) tree = tree->next;
4137
if(mixt_tree->mod->r_mat == tree->mod->r_mat &&
4138
mixt_tree->mod->whichmodel == tree->mod->whichmodel &&
4139
!strcmp(mixt_tree->mod->custom_mod_string->s,
4140
tree->mod->custom_mod_string->s) &&
4141
!strcmp(mixt_tree->mod->aa_rate_mat_file->s,
4142
tree->mod->aa_rate_mat_file->s)) link_rmat[c] = cc_rmat;
4153
mixt_tree = mixt_tree->next;
4157
PhyML_Fprintf(fp,"\n");
4158
strcpy(param,"State frequencies ");
4159
PhyML_Fprintf(fp," %18s",param);
4161
tree = cpy_mixt_tree;
4165
if(tree->is_mixt_tree) tree = tree->next;
4166
PhyML_Fprintf(fp,"%2c ",link_efrq[c]);
4173
PhyML_Fprintf(fp,"\n");
4174
strcpy(param,"Branch lengths ");
4175
PhyML_Fprintf(fp," %18s",param);
4177
tree = cpy_mixt_tree;
4181
if(tree->is_mixt_tree) tree = tree->next;
4182
PhyML_Fprintf(fp,"%2c ",link_lens[c]);
4188
PhyML_Fprintf(fp,"\n");
4189
strcpy(param,"Rate matrix ");
4190
PhyML_Fprintf(fp," %18s",param);
4192
tree = cpy_mixt_tree;
4196
if(tree->is_mixt_tree) tree = tree->next;
4197
PhyML_Fprintf(fp,"%2c ",link_rmat[c]);
4203
PhyML_Fprintf(fp,"\n");
4204
PhyML_Fprintf(fp," ------------------");
4205
tree = cpy_mixt_tree;
4209
if(tree->is_mixt_tree) tree = tree->next;
4210
PhyML_Fprintf(fp,"---");
4214
PhyML_Fprintf(fp,"\n");
4218
tree = cpy_mixt_tree;
4222
PhyML_Fprintf(fp,"\n");
4223
PhyML_Fprintf(fp,"\n. Tree estimated from data partition %d",c++);
4224
s = Write_Tree(tree->next,NO); /*! tree->next is not a mixt_tree so edge lengths
4225
are not averaged over when writing the tree out. */
4226
PhyML_Fprintf(fp,"\n %s",s);
4228
tree = tree->next_mixt;
4238
mixt_tree = cpy_mixt_tree;
4241
//////////////////////////////////////////////////////////////
4242
//////////////////////////////////////////////////////////////
4244
void PhyML_XML(char *xml_filename)
4247
xml_node *root,*p_elem,*m_elem,*parent,*instance;
4251
t_tree *tree,*mixt_tree,*root_tree;
4254
int i,j,n_components;
4257
scalar_dbl **lens,**lens_var,**ori_lens,**lens_old,**lens_var_old,**ori_lens_old,**ori_lens_var,**ori_lens_var_old;
4267
fp = fopen(xml_filename,"r");
4270
PhyML_Printf("\n== Could not find the XML file '%s'.\n",xml_filename);
4275
root = XML_Load_File(fp);
4279
PhyML_Printf("\n== Encountered an issue while loading the XML file.\n");
4283
class_num = (int *)mCalloc(N_MAX_MIXT_CLASSES,sizeof(int));
4284
For(i,N_MAX_MIXT_CLASSES) class_num[i] = i;
4286
component = (char *)mCalloc(T_MAX_NAME,sizeof(char));
4298
lens_var_old = NULL;
4304
ori_lens_old = NULL;
4305
ori_lens_var = NULL;
4306
ori_lens_var_old = NULL;
4309
// Make sure there are no duplicates in node's IDs
4310
XML_Check_Duplicate_ID(root);
4313
XML_Count_Number_Of_Node_With_Name("topology",&count,root);
4317
PhyML_Printf("\n== There should not more than one 'topology' node.");
4318
PhyML_Printf("\n== Found %d. Please fix your XML file",count);
4323
PhyML_Printf("\n== There should be at least one 'topology' node.");
4324
PhyML_Printf("\n== Found none. Please fix your XML file");
4328
p_elem = XML_Search_Node_Name("phyml",NO,p_elem);
4332
outputfile = XML_Get_Attribute_Value(p_elem,"output.file");
4336
PhyML_Printf("\n== The 'outputfile' attribute in 'phyml' tag is mandatory.");
4337
PhyML_Printf("\n== Please amend your XML file accordingly.");
4341
io = (option *)Make_Input();
4342
Set_Defaults_Input(io);
4344
s = XML_Get_Attribute_Value(p_elem,"run.id");
4347
io->append_run_ID = YES;
4348
strcpy(io->run_id_string,s);
4351
strcpy(io->out_tree_file,outputfile);
4352
if(io->append_run_ID) { strcat(io->out_tree_file,"_"); strcat(io->out_tree_file,io->run_id_string); }
4353
strcat(io->out_tree_file,"_phyml_tree");
4354
strcpy(io->out_stats_file,outputfile);
4355
if(io->append_run_ID) { strcat(io->out_stats_file,"_"); strcat(io->out_stats_file,io->run_id_string); }
4356
strcat(io->out_stats_file,"_phyml_stats");
4357
io->fp_out_tree = Openfile(io->out_tree_file,1);
4358
io->fp_out_stats = Openfile(io->out_stats_file,1);
4360
s = XML_Get_Attribute_Value(p_elem,"print.trace");
4364
select = XML_Validate_Attr_Int(s,6,
4370
io->print_trace = YES;
4371
strcpy(io->out_trace_file,outputfile);
4372
if(io->append_run_ID) { strcat(io->out_trace_file,"_"); strcat(io->out_trace_file,io->run_id_string); }
4373
strcat(io->out_trace_file,"_phyml_trace");
4374
io->fp_out_trace = Openfile(io->out_trace_file,1);
4379
s = XML_Get_Attribute_Value(p_elem,"branch.test");
4382
if(!strcmp(s,"aLRT"))
4384
io->ratio_test = ALRTSTAT;
4386
else if(!strcmp(s,"aBayes"))
4388
io->ratio_test = ABAYES;
4390
else if(!strcmp(s,"SH"))
4392
io->ratio_test = SH;
4394
else if(!strcmp(s,"no") || !strcmp(s,"none"))
4400
PhyML_Printf("\n== '%s' is not a valid option for 'branch.test'.",s);
4401
PhyML_Printf("\n== Please use 'aLRT' or 'aBayes' or 'SH'.");
4407
/*! Read all partitionelem nodes and mixturelem nodes in each of them
4411
p_elem = XML_Search_Node_Name("partitionelem",YES,p_elem);
4413
if(p_elem == NULL) break;
4415
buff = (option *)Make_Input();
4416
Set_Defaults_Input(buff);
4418
io->next->prev = io;
4430
/*! Set the datatype (required when compressing data)
4433
dt = XML_Get_Attribute_Value(p_elem,"data.type");
4436
PhyML_Printf("\n== Please specify the type of data ('aa' or 'nt') for partition element '%s'",
4437
XML_Get_Attribute_Value(p_elem,"id"));
4438
PhyML_Printf("\n== Syntax: 'data.type=\"aa\"' or 'data.type=\"nt\"'");
4442
select = XML_Validate_Attr_Int(dt,2,"aa","nt");
4457
PhyML_Printf("\n== Unknown data type. Must be either 'aa' or 'nt'.");
4462
char *format = NULL;
4463
format = XML_Get_Attribute_Value(p_elem,"format");
4466
if(!strcmp(format,"interleave") ||
4467
!strcmp(format,"interleaved"))
4469
io->interleaved = YES;
4471
else if(!strcmp(format,"sequential"))
4473
io->interleaved = NO;
4478
/*! Attach a model to this io struct
4480
io->mod = (t_mod *)Make_Model_Basic();
4481
Set_Defaults_Model(io->mod);
4482
io->mod->ras->n_catg = 1;
4486
/*! Attach an optimization structure to this model
4488
iomod->s_opt = (t_opt *)Make_Optimiz();
4489
Set_Defaults_Optimiz(iomod->s_opt);
4491
iomod->s_opt->opt_kappa = NO;
4492
iomod->s_opt->opt_lambda = NO;
4493
iomod->s_opt->opt_rr = NO;
4497
alignment = XML_Get_Attribute_Value(p_elem,"file.name");
4501
PhyML_Printf("\n== 'file.name' tag is mandatory. Please amend your");
4502
PhyML_Printf("\n== XML file accordingly.");
4506
strcpy(io->in_align_file,alignment);
4508
/*! Open pointer to alignment
4510
io->fp_in_align = Openfile(io->in_align_file,0);
4512
/*! Load sequence file
4514
io->data = Get_Seq(io);
4516
/*! Close pointer to alignment
4518
fclose(io->fp_in_align);
4520
/*! Compress alignment
4522
io->cdata = Compact_Data(io->data,io);
4524
/*! Free uncompressed alignment
4526
Free_Seq(io->data,io->n_otu);
4528
/*! Create new mixture tree
4530
buff = (t_tree *)Make_Tree_From_Scratch(io->cdata->n_otu,io->cdata);
4534
mixt_tree->next_mixt = buff;
4535
mixt_tree->next_mixt->prev_mixt = mixt_tree;
4536
mixt_tree = mixt_tree->next_mixt;
4537
mixt_tree->dp = mixt_tree->prev_mixt->dp+1;
4539
else mixt_tree = buff;
4541
/*! Connect mixt_tree to io struct
4545
/*! Connect mixt_tree to model struct
4547
mixt_tree->mod = iomod;
4549
/*! mixt_tree is a mixture tree
4551
mixt_tree->is_mixt_tree = YES;
4553
/*! mixt_tree is a mixture tree
4555
mixt_tree->mod->is_mixt_mod = YES;
4557
/*! Connect mixt_tree to compressed data
4559
mixt_tree->data = io->cdata;
4561
/*! Set total number of patterns
4563
mixt_tree->n_pattern = io->cdata->crunch_len;
4565
/*! Remove branch lengths from mixt_tree */
4566
For(i,2*mixt_tree->n_otu-1)
4568
Free_Scalar_Dbl(mixt_tree->a_edges[i]->l);
4569
Free_Scalar_Dbl(mixt_tree->a_edges[i]->l_old);
4570
Free_Scalar_Dbl(mixt_tree->a_edges[i]->l_var);
4571
Free_Scalar_Dbl(mixt_tree->a_edges[i]->l_var_old);
4574
/*! Connect last tree of the mixture for the
4575
previous partition element to the next mixture tree
4579
tree->next = mixt_tree;
4580
mixt_tree->prev = tree;
4583
/*! Do the same for the model
4591
if(!root_tree) root_tree = mixt_tree;
4593
/*! Process all the mixtureelem tags in this partition element
4603
m_elem = XML_Search_Node_Name("mixtureelem",YES,m_elem);
4604
if(m_elem == NULL) break;
4606
if(!strcmp(m_elem->name,"mixtureelem"))
4610
/*! Rewind tree and model when processing a new mixtureelem node
4612
if(first_m_elem > 1)
4614
while(tree->prev && tree->prev->is_mixt_tree == NO) { tree = tree->prev; } // tree = tree->prev->next;
4615
while(mod->prev && mod->prev->is_mixt_mod == NO) { mod = mod->prev; } // mod = mod->prev->next;
4618
/*! Read and process model components
4621
list = XML_Get_Attribute_Value(m_elem,"list");
4624
For(i,(int)strlen(list)) if(list[i] == ',') j++;
4626
if(j != n_components && first_m_elem > 1)
4628
PhyML_Printf("\n== Discrepancy in the number of elements in nodes 'mixtureelem' partitionelem id '%s'",p_elem->id);
4629
PhyML_Printf("\n== Check 'mixturelem' node with list '%s'",list);
4635
component[0] = '\0';
4638
if(list[j] == ',' || j == (int)strlen(list))
4640
/*! Reading a new component
4643
if(first_m_elem == YES) // Only true when processing the first mixtureelem node
4650
this_tree = (t_tree *)Make_Tree_From_Scratch(io->cdata->n_otu,io->cdata);
4652
/*! Update the number of mixtures */
4653
iomod->n_mixt_classes++;
4657
tree->next = this_tree;
4658
tree->next->prev = tree;
4662
mixt_tree->next = this_tree;
4663
mixt_tree->next->prev = mixt_tree;
4667
tree->mixt_tree = mixt_tree;
4670
/*! Create a new model
4672
this_mod = (t_mod *)Make_Model_Basic();
4673
Set_Defaults_Model(this_mod);
4674
this_mod->ras->n_catg = 1;
4680
mod->next = this_mod;
4681
mod->next->prev = mod;
4685
this_mod->prev = iomod;
4689
if(!iomod->next) iomod->next = mod;
4692
mod->s_opt = (t_opt *)Make_Optimiz();
4693
Set_Defaults_Optimiz(mod->s_opt);
4695
mod->s_opt->opt_alpha = NO;
4696
mod->s_opt->opt_pinvar = NO;
4698
tree->data = io->cdata;
4699
tree->n_pattern = io->cdata->crunch_len;
4703
if(tree->n_pattern != tree->prev->n_pattern)
4705
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
4710
/*! Read a component
4712
component[i] = '\0';
4713
if(j != (int)strlen(list)-1) i = 0;
4715
/*! Find which node this ID corresponds to
4717
instance = XML_Search_Node_ID(component,YES,root);
4721
PhyML_Printf("\n== Could not find a node with id:'%s'.",component);
4722
PhyML_Printf("\n== Problem with 'mixtureelem' node, list '%s'.",list);
4726
if(!instance->parent)
4728
PhyML_Printf("\n== Node '%s' with id:'%s' has no parent.",instance->name,component);
4732
parent = instance->parent;
4734
////////////////////////////////////////
4735
// SUBSTITUTION MODEL //
4736
////////////////////////////////////////
4738
if(!strcmp(parent->name,"ratematrices"))
4740
/* ! First time we process this 'instance' node which has this 'ratematrices' parent */
4741
if(instance->ds->obj == NULL)
4743
Make_Ratematrice_From_XML_Node(instance, io, mod);
4747
/*! Connect the data structure n->ds to mod->r_mat */
4748
ds->obj = (t_rmat *)(mod->r_mat);
4750
/*! Create and connect the data structure n->ds->next to mod->kappa */
4751
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4753
ds->obj = (scalar_dbl *)(mod->kappa);
4755
/*! Create and connect the data structure n->ds->next to mod->s_opt->opt_kappa */
4756
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4758
ds->obj = (int *)(&mod->s_opt->opt_kappa);
4760
/*! Create and connect the data structure n->ds->next to mod->s_opt->opt_rr */
4761
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4763
ds->obj = (int *)(&mod->s_opt->opt_rr);
4765
/*! Create and connect the data structure n->ds->next to mod->whichmodel */
4766
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4768
ds->obj = (int *)(&mod->whichmodel);
4770
/*! Create and connect the data structure n->ds->next to mod->modelname */
4771
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4773
ds->obj = (t_string *)(mod->modelname);
4775
/*! Create and connect the data structure n->ds->next to mod->ns */
4776
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4778
ds->obj = (int *)(&mod->ns);
4780
/*! Create and connect the data structure n->ds->next to mod->modelname */
4781
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4783
ds->obj = (t_string *)(mod->custom_mod_string);
4787
/*! Connect to already extisting r_mat & kappa structs. */
4792
mod->r_mat = (t_rmat *)ds->obj;
4795
Free_Scalar_Dbl(mod->kappa);
4796
mod->kappa = (scalar_dbl *)ds->obj;
4799
mod->s_opt->opt_kappa = *((int *)ds->obj);
4802
mod->s_opt->opt_rr = *((int *)ds->obj);
4805
mod->whichmodel = *((int *)ds->obj);
4808
Free_String(mod->modelname);
4809
mod->modelname = (t_string *)ds->obj;
4812
mod->ns = *((int *)ds->obj);
4815
Free_String(mod->custom_mod_string);
4816
mod->custom_mod_string = (t_string *)ds->obj;
4820
////////////////////////////////////////
4822
////////////////////////////////////////
4824
else if(!strcmp(parent->name,"equfreqs"))
4827
// If n->ds == NULL, the corrresponding node data structure, n->ds, has not
4828
// been initialized. If not, do nothing.
4829
if(instance->ds->obj == NULL)
4831
Make_Efrq_From_XML_Node(instance,io,mod);
4836
ds->obj = (t_efrq *)mod->e_frq;
4838
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4840
ds->obj = (int *)(&mod->s_opt->opt_state_freq);
4842
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4844
ds->obj = (int *)(&mod->s_opt->user_state_freq);
4846
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4848
ds->obj = (vect_dbl *)(mod->user_b_freq);
4852
// Connect the data structure n->ds to mod->e_frq
4855
mod->e_frq = (t_efrq *)ds->obj;
4858
mod->s_opt->opt_state_freq = *((int *)ds->obj);
4861
mod->s_opt->user_state_freq = *((int *)ds->obj);
4864
Free_Vect_Dbl(mod->user_b_freq);
4865
mod->user_b_freq = (vect_dbl *)ds->obj;
4869
//////////////////////////////////////////
4871
//////////////////////////////////////////
4873
else if(!strcmp(parent->name,"topology"))
4875
if(parent->ds->obj == NULL) Make_Topology_From_XML_Node(instance,io,iomod);
4880
ds->obj = (int *)(& buff);
4883
//////////////////////////////////////////
4885
//////////////////////////////////////////
4887
else if(!strcmp(parent->name,"siterates"))
4889
char *rate_value = NULL;
4890
/* scalar_dbl *r; */
4893
/*! First time we process this 'siterates' node, check that its format is valid.
4894
and process it afterwards.
4896
if(parent->ds->obj == NULL)
4900
Make_RAS_From_XML_Node(parent,iomod);
4904
ds->obj = (t_ras *)iomod->ras;
4906
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4908
ds->obj = (int *)(&iomod->s_opt->opt_alpha);
4910
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4912
ds->obj = (int *)(&iomod->s_opt->opt_free_mixt_rates);
4914
else /*! Connect ras struct to already defined one. Same for opt_alpha & opt_free_mixt_rates */
4916
if(iomod->ras != (t_ras *)parent->ds->obj) Free_RAS(iomod->ras);
4917
iomod->ras = (t_ras *)parent->ds->obj;
4918
iomod->s_opt->opt_alpha = *((int *)parent->ds->next->obj);
4919
iomod->s_opt->opt_free_mixt_rates = *((int *)parent->ds->next->next->obj);
4922
rate_value = XML_Get_Attribute_Value(instance,"init.value");
4925
if(rate_value) val = String_To_Dbl(rate_value);
4927
if(instance->ds->obj == NULL)
4929
instance->ds->obj = (phydbl *)(&val);
4930
instance->ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4931
instance->ds->next->obj = (int *)(class_num + class_number);
4933
iomod->ras->gamma_rr->v[class_number] = val;
4934
iomod->ras->gamma_rr_unscaled->v[class_number] = val;
4936
iomod->ras->init_rr = NO;
4938
if(Are_Equal(val,0.0,1E-20) == NO) class_number++;
4942
/*! Note: ras is already connected to the relevant t_ds stucture. No need
4943
to connect ras->gamma_rr or ras->p_invar */
4946
if(Are_Equal(val,0.0,1E-20))
4948
mod->ras->invar = YES;
4952
mod->ras->parent_class_number = *((int *)instance->ds->next->obj);
4955
xml_node *orig_w = NULL;
4956
orig_w = XML_Search_Node_Attribute_Value("appliesto",instance->id,YES,instance->parent);
4961
weight = XML_Get_Attribute_Value(orig_w,"value");
4962
if(mod->ras->invar == YES)
4964
iomod->ras->pinvar->v = String_To_Dbl(weight);
4969
class = *((int *)instance->ds->next->obj);
4970
iomod->ras->gamma_r_proba->v[class] = String_To_Dbl(weight);
4971
iomod->ras->init_r_proba = NO;
4976
//////////////////////////////////////////////
4977
// BRANCH LENGTHS //
4978
//////////////////////////////////////////////
4980
else if(!strcmp(parent->name,"branchlengths"))
4985
n_otu = tree->n_otu;
4987
if(instance->ds->obj == NULL)
4991
ori_lens = (scalar_dbl **)mCalloc(2*tree->n_otu-1,sizeof(scalar_dbl *));
4992
ori_lens_old = (scalar_dbl **)mCalloc(2*tree->n_otu-1,sizeof(scalar_dbl *));
4994
ori_lens_var = (scalar_dbl **)mCalloc(2*tree->n_otu-1,sizeof(scalar_dbl *));
4995
ori_lens_var_old = (scalar_dbl **)mCalloc(2*tree->n_otu-1,sizeof(scalar_dbl *));
4998
lens_old = ori_lens_old;
5000
lens_var = ori_lens_var;
5001
lens_var_old = ori_lens_var_old;
5003
lens_size = 2*tree->n_otu-1;
5007
ori_lens = (scalar_dbl **)mRealloc(ori_lens,2*tree->n_otu-1+lens_size,sizeof(scalar_dbl *));
5008
ori_lens_old = (scalar_dbl **)mRealloc(ori_lens_old,2*tree->n_otu-1+lens_size,sizeof(scalar_dbl *));
5010
ori_lens_var = (scalar_dbl **)mRealloc(ori_lens_var,2*tree->n_otu-1+lens_size,sizeof(scalar_dbl *));
5011
ori_lens_var_old = (scalar_dbl **)mRealloc(ori_lens_var_old,2*tree->n_otu-1+lens_size,sizeof(scalar_dbl *));
5013
lens = ori_lens + lens_size;;
5014
lens_old = ori_lens_old + lens_size;;
5016
lens_var = ori_lens_var + lens_size;;
5017
lens_var_old = ori_lens_var_old + lens_size;;
5019
lens_size += 2*tree->n_otu-1;
5022
For(i,2*tree->n_otu-1)
5024
lens[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5025
Init_Scalar_Dbl(lens[i]);
5027
lens_old[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5028
Init_Scalar_Dbl(lens_old[i]);
5030
lens_var[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5031
Init_Scalar_Dbl(lens_var[i]);
5033
lens_var_old[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5034
Init_Scalar_Dbl(lens_var_old[i]);
5036
Free_Scalar_Dbl(tree->a_edges[i]->l);
5037
Free_Scalar_Dbl(tree->a_edges[i]->l_old);
5039
Free_Scalar_Dbl(tree->a_edges[i]->l_var);
5040
Free_Scalar_Dbl(tree->a_edges[i]->l_var_old);
5043
tree->prev->a_edges[i]->l == mixt_tree->a_edges[i]->l &&
5044
tree->prev->is_mixt_tree == NO)
5046
PhyML_Printf("\n== %p %p",tree->a_edges[i]->l,mixt_tree->a_edges[i]->l);
5047
PhyML_Printf("\n== Only one set of edge lengths is allowed ");
5048
PhyML_Printf("\n== in a 'partitionelem'. Please fix your XML file.");
5053
char *opt_bl = NULL;
5054
opt_bl = XML_Get_Attribute_Value(instance,"optimise.lens");
5058
if(!strcmp(opt_bl,"yes") || !strcmp(opt_bl,"true"))
5060
iomod->s_opt->opt_bl = YES;
5064
iomod->s_opt->opt_bl = NO;
5070
ds->obj = (scalar_dbl **)lens;
5072
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
5074
ds->obj = (scalar_dbl **)lens_old;
5076
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
5078
ds->obj = (scalar_dbl **)lens_var;
5080
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
5082
ds->obj = (scalar_dbl **)lens_var_old;
5084
ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
5086
ds->obj = (int *)(&iomod->s_opt->opt_bl);
5090
For(i,2*tree->n_otu-1)
5092
Free_Scalar_Dbl(tree->a_edges[i]->l);
5093
Free_Scalar_Dbl(tree->a_edges[i]->l_old);
5094
Free_Scalar_Dbl(tree->a_edges[i]->l_var);
5095
Free_Scalar_Dbl(tree->a_edges[i]->l_var_old);
5100
lens = (scalar_dbl **)ds->obj;
5103
lens_old = (scalar_dbl **)ds->obj;
5106
lens_var = (scalar_dbl **)ds->obj;
5109
lens_var_old = (scalar_dbl **)ds->obj;
5112
iomod->s_opt->opt_bl = *((int *)ds->obj);
5115
if(n_otu != tree->n_otu)
5117
PhyML_Printf("\n== All the data sets should display the same number of sequences.");
5118
PhyML_Printf("\n== Found at least one data set with %d sequences and one with %d sequences.",n_otu,tree->n_otu);
5122
For(i,2*tree->n_otu-1)
5124
tree->a_edges[i]->l = lens[i];
5125
mixt_tree->a_edges[i]->l = lens[i];
5126
tree->a_edges[i]->l_old = lens_old[i];
5127
mixt_tree->a_edges[i]->l_old = lens_old[i];
5129
tree->a_edges[i]->l_var = lens_var[i];
5130
mixt_tree->a_edges[i]->l_var = lens_var[i];
5131
tree->a_edges[i]->l_var_old = lens_var_old[i];
5132
mixt_tree->a_edges[i]->l_var_old = lens_var_old[i];
5136
///////////////////////////////////////////////
5137
///////////////////////////////////////////////
5138
///////////////////////////////////////////////
5140
if(first_m_elem > 1) // Done with this component, move to the next tree and model
5142
if(tree->next) tree = tree->next;
5143
if(mod->next) mod = mod->next;
5147
else if(list[j] != ' ')
5149
component[i] = list[j];
5153
if(j == (int)strlen(list)+1) break;
5155
} // end of mixtureelem processing
5156
} // end of partitionelem processing
5163
if(ori_lens) Free(ori_lens);
5164
if(ori_lens_old) Free(ori_lens_old);
5165
if(ori_lens_var) Free(ori_lens_var);
5166
if(ori_lens_var_old) Free(ori_lens_var_old);
5168
while(io->prev != NULL) io = io->prev;
5169
while(mixt_tree->prev != NULL) mixt_tree = mixt_tree->prev;
5171
/*! Finish making the models */
5172
mod = mixt_tree->mod;
5175
Make_Model_Complete(mod);
5180
Check_Taxa_Sets(mixt_tree);
5182
Check_Mandatory_XML_Node(root,"phyml");
5183
Check_Mandatory_XML_Node(root,"topology");
5184
Check_Mandatory_XML_Node(root,"branchlengths");
5185
Check_Mandatory_XML_Node(root,"ratematrices");
5186
Check_Mandatory_XML_Node(root,"equfreqs");
5187
Check_Mandatory_XML_Node(root,"siterates");
5188
Check_Mandatory_XML_Node(root,"partitionelem");
5189
Check_Mandatory_XML_Node(root,"mixtureelem");
5191
if(!io->mod->s_opt->random_input_tree) io->mod->s_opt->n_rand_starts = 1;
5193
char *most_likely_tree = NULL;
5194
phydbl best_lnL = UNLIKELY;
5196
For(num_rand_tree,io->mod->s_opt->n_rand_starts)
5198
/*! Initialize the models */
5199
mod = mixt_tree->mod;
5202
Init_Model(mod->io->cdata,mod,io);
5207
/* printf("\n%f %f %f %f", */
5208
/* mixt_tree->mod->ras->gamma_r_proba->v[0], */
5209
/* mixt_tree->mod->ras->gamma_r_proba->v[1], */
5210
/* mixt_tree->mod->ras->gamma_r_proba->v[2], */
5211
/* mixt_tree->mod->ras->gamma_r_proba->v[3]); */
5214
Print_Data_Structure(NO,stdout,mixt_tree);
5216
t_tree *bionj_tree = NULL;
5217
switch(mixt_tree->io->in_tree)
5219
case 2: /*! user-defined input tree */
5221
if(!mixt_tree->io->fp_in_tree)
5223
PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
5228
Copy the user tree to all the tree structures
5230
tree = Read_User_Tree(mixt_tree->io->cdata,
5239
/*! Build a BioNJ tree from the analysis of
5240
the first partition element
5242
tree = Dist_And_BioNJ(mixt_tree->data,
5245
bionj_tree = (t_tree *)Make_Tree_From_Scratch(mixt_tree->n_otu,mixt_tree->data);
5246
Copy_Tree(tree,bionj_tree);
5250
tree = (t_tree *)Make_Tree_From_Scratch(mixt_tree->n_otu,mixt_tree->data);
5251
Copy_Tree(bionj_tree,tree);
5258
Copy_Tree(tree,mixt_tree);
5260
if(bionj_tree) Free_Tree(bionj_tree);
5262
if(io->mod->s_opt->random_input_tree)
5264
PhyML_Printf("\n\n. [%3d/%3d]",num_rand_tree+1,io->mod->s_opt->n_rand_starts);
5265
Random_Tree(mixt_tree);
5271
if(tree != mixt_tree) Copy_Tree(mixt_tree,tree);
5272
Connect_CSeqs_To_Nodes(tree->data,tree);
5278
/*! Initialize t_beg in each mixture tree */
5282
time(&(tree->t_beg));
5283
tree = tree->next_mixt;
5287
Prepare_Tree_For_Lk(mixt_tree);
5289
MIXT_Chain_All(mixt_tree);
5291
/*! Check that all the edges in a mixt_tree at pointing
5292
to a single set of lengths
5297
MIXT_Check_Single_Edge_Lens(tree);
5298
tree = tree->next_mixt;
5303
/*! Turn all branches to ON state */
5307
MIXT_Turn_Branches_OnOff(ON,tree);
5308
tree = tree->next_mixt;
5315
X 2) ALLOW MORE THAN ONE VECTOR OF BRANCH LENGTHS -> NEED TO
5316
LOOK CAREFULY AT WHAT IT MEANS FOR SPR, SIMULTANEOUS NNI MOVES
5317
PRUNE AND REGRAFT...
5318
X 3) Branch lengths in Prune and Regarft -> make sure you don't apply
5319
change several times
5320
4) Make sure you can't have the same branch lengths with distinct
5322
5) Finish rewritting Opt_Free_Param
5323
X 6) Make sure you can have only one set of branch lengths per partition
5324
7) BIONJ starting tree
5325
X 8) When iomod->ras->invar = YES, check that only one mod has mod->ras->invar = YES;
5326
9) Reflect changes on simu.c to spr.c, especially regarding invariants
5327
10) Rough_SPR to replace SPR_Shuffle and first step in nni (refining tree)
5328
11) Tree corresponding to invar class is never really used (except for testing
5329
whether tree->mod->ras->invar==YES). Do we need to use memory for this?
5330
X 12) Check that mixt_tree->mod->ras->n_catg corresponds to the same number of
5331
trees in the mixture (do that for each mixt_tree).
5332
13) Change tree->a_edges to tree->a_edges (t is for type a is for array)
5333
14) Change tree->a_nodes to tree->a_nodes (t is for type a is for array)
5334
X 15) tstv should be shared! Check the following:
5335
// Connect the data structure n->ds to mod->r_mat
5336
mod->r_mat = (t_rmat *)instance->ds->obj;
5337
mod->kappa = (scalar_dbl *)instance->ds->next->obj;
5338
16) Replace s_opt struct by opt flag for each param?
5339
X 17) Check the sequences and tree have the same number of otus
5340
X 18) Avoid transformation to lower case of attribute values when processing XML.
5341
X 19) Break mixture: break nodes and branches too.
5343
X 21) Set_Alias_Subpatt
5344
X 22) Check that all partition elements have the same set of taxa.
5345
23) What if ratematrices are not specified ? Default is ok?
5346
X 24) Set upper and lower bounds on Br_Len_Brent
5350
MIXT_Check_Invar_Struct_In_Each_Partition_Elem(mixt_tree);
5351
MIXT_Check_RAS_Struct_In_Each_Partition_Elem(mixt_tree);
5353
Set_Both_Sides(YES,mixt_tree);
5355
if(mixt_tree->mod->s_opt->opt_topo)
5357
if(mixt_tree->mod->s_opt->topo_search == NNI_MOVE) Simu_Loop(mixt_tree);
5358
else if(mixt_tree->mod->s_opt->topo_search == SPR_MOVE) Speed_Spr_Loop(mixt_tree);
5359
else Best_Of_NNI_And_SPR(mixt_tree);
5363
if(mixt_tree->mod->s_opt->opt_subst_param ||
5364
mixt_tree->mod->s_opt->opt_bl)
5367
Round_Optimize(mixt_tree,mixt_tree->data,ROUND_MAX);
5376
PhyML_Printf("\n\n. Log-likelihood = %f",mixt_tree->c_lnL);
5378
if((num_rand_tree == io->mod->s_opt->n_rand_starts-1) && (io->mod->s_opt->random_input_tree))
5381
io->mod->s_opt->random_input_tree = NO;
5384
Br_Len_Involving_Invar(mixt_tree);
5385
Rescale_Br_Len_Multiplier_Tree(mixt_tree);
5387
/*! Print the tree estimated using the current random (or BioNJ) starting tree */
5388
if(io->mod->s_opt->n_rand_starts > 1)
5390
Print_Tree(io->fp_out_trees,mixt_tree);
5394
if(mixt_tree->c_lnL > best_lnL)
5396
best_lnL = mixt_tree->c_lnL;
5397
if(most_likely_tree) Free(most_likely_tree);
5398
if(io->ratio_test) aLRT(mixt_tree);
5399
most_likely_tree = Write_Tree(mixt_tree,NO);
5400
mixt_tree->lock_topo = NO;
5403
/*! Initialize t_end in each mixture tree */
5407
time(&(tree->t_current));
5408
tree = tree->next_mixt;
5412
Print_Data_Structure(YES,mixt_tree->io->fp_out_stats,mixt_tree);
5414
Free_Spr_List(mixt_tree);
5415
Free_Tree_Pars(mixt_tree);
5416
Free_Tree_Lk(mixt_tree);
5417
Free_Triplet(mixt_tree->triplet_struct);
5420
/*! Print the most likely tree in the output file */
5421
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));
5422
PhyML_Fprintf(mixt_tree->io->fp_out_tree,"%s\n",most_likely_tree);
5425
/*! Bootstrap analysis */
5426
MIXT_Bootstrap(most_likely_tree,root);
5428
while(io->prev != NULL) io = io->prev;
5430
Free(most_likely_tree);
5435
Free_Cseq(tree->data);
5436
tree = tree->next_mixt;
5443
Free_Optimiz(tree->mod->s_opt);
5448
Free_Model_Complete(mixt_tree->mod);
5449
Free_Model_Basic(mixt_tree->mod);
5450
Free_Tree(mixt_tree);
5452
if(io->fp_out_trees) fclose(io->fp_out_trees);
5453
if(io->fp_out_tree) fclose(io->fp_out_tree);
5454
if(io->fp_out_stats) fclose(io->fp_out_stats);
5456
XML_Free_XML_Tree(root);
5463
//////////////////////////////////////////////////////////////
5464
//////////////////////////////////////////////////////////////
5466
/*! Check that the same nodes in the different mixt_trees are
5467
connected to the same taxa
5469
void Check_Taxa_Sets(t_tree *mixt_tree)
5481
if(strcmp(tree->a_nodes[i]->name,tree->next->a_nodes[i]->name))
5483
PhyML_Printf("\n== There seems to be a problem in one (or more) of your");
5484
PhyML_Printf("\n== sequence alignment. PhyML could not match taxon");
5485
PhyML_Printf("\n== '%s' found in file '%s' with any of the taxa",tree->a_nodes[i]->name,tree->io->in_align_file);
5486
PhyML_Printf("\n== listed in file '%s'.",tree->next->io->in_align_file);
5497
//////////////////////////////////////////////////////////////
5498
//////////////////////////////////////////////////////////////
5500
void Make_Ratematrice_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5505
model = XML_Get_Attribute_Value(instance,"model");
5509
PhyML_Printf("\n== Poorly formated XML file.");
5510
PhyML_Printf("\n== Attribute 'model' is mandatory in a <ratematrix> node.");
5514
select = XML_Validate_Attr_Int(model,26,
5545
if(io->datatype != NT)
5547
PhyML_Printf("\n== Data type and selected model are incompatible");
5554
if(io->datatype != AA)
5556
PhyML_Printf("\n== Data type and selected model are incompatible");
5561
io->mod->ns = mod->ns;
5563
mod->r_mat = (t_rmat *)Make_Rmat(mod->ns);
5564
Init_Rmat(mod->r_mat);
5566
/*! Set model number & name */
5567
mod->whichmodel = Set_Whichmodel(select);
5568
Set_Model_Name(mod);
5570
if(mod->whichmodel == K80 ||
5571
mod->whichmodel == HKY85 ||
5572
mod->whichmodel == TN93)
5574
char *tstv,*opt_tstv;
5576
tstv = XML_Get_Attribute_Value(instance,"tstv");
5580
mod->s_opt->opt_kappa = NO;
5581
mod->kappa->v = String_To_Dbl(tstv);
5585
mod->s_opt->opt_kappa = YES;
5588
opt_tstv = XML_Get_Attribute_Value(instance,"optimise.tstv");
5592
if(!strcmp(opt_tstv,"true") || !strcmp(opt_tstv,"yes"))
5594
mod->s_opt->opt_kappa = YES;
5598
mod->s_opt->opt_kappa = NO;
5604
mod->s_opt->opt_kappa = NO;
5607
if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM)
5611
opt_rr = XML_Get_Attribute_Value(instance,"optimise.rr");
5615
if(!strcmp(opt_rr,"yes") || !strcmp(opt_rr,"true"))
5617
mod->s_opt->opt_rr = YES;
5622
/*! Custom model for nucleotide sequences. Read the corresponding
5624
if(mod->whichmodel == CUSTOM)
5626
char *model_code = NULL;
5628
model_code = XML_Get_Attribute_Value(instance,"model.code");
5632
PhyML_Printf("\n== No valid 'model.code' attribute could be found.\n");
5633
PhyML_Printf("\n== Please fix your XML file.\n");
5638
strcpy(mod->custom_mod_string->s,model_code);
5643
/*! Custom model for amino-acids. Read in the rate matrix file */
5644
if(mod->whichmodel == CUSTOMAA)
5648
r_mat_file = XML_Get_Attribute_Value(instance,"ratematrix.file");
5652
PhyML_Printf("\n== No valid 'ratematrix.file' attribute could be found.");
5653
PhyML_Printf("\n== Please fix your XML file.\n");
5658
mod->fp_aa_rate_mat = Openfile(r_mat_file,0);
5659
strcpy(mod->aa_rate_mat_file->s,r_mat_file);
5667
buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
5668
if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
5670
mod->s_opt->opt_rmat_weight = YES;
5674
mod->s_opt->opt_rmat_weight = NO;
5678
//////////////////////////////////////////////////////////////
5679
//////////////////////////////////////////////////////////////
5681
void Make_Efrq_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5685
mod->e_frq = (t_efrq *)Make_Efrq(mod->ns);
5686
Init_Efrq(mod->e_frq);
5688
buff = XML_Get_Attribute_Value(instance,"optimise.freqs");
5692
if(!strcmp(buff,"yes") || !strcmp(buff,"true"))
5694
if(io->datatype == AA)
5696
PhyML_Printf("\n== Option 'optimise.freqs' set to 'yes' (or 'true')");
5697
PhyML_Printf("\n== is not allowed with amino-acid data.");
5700
mod->s_opt->opt_state_freq = YES;
5705
buff = XML_Get_Attribute_Value(instance,"aa.freqs");
5709
if(!strcmp(buff,"empirical"))
5711
if(io->datatype == AA)
5713
mod->s_opt->opt_state_freq = YES;
5715
else if(io->datatype == NT)
5717
mod->s_opt->opt_state_freq = NO;
5724
buff = XML_Get_Attribute_Value(instance,"base.freqs");
5728
if(io->datatype == AA)
5730
PhyML_Printf("\n== Option 'base.freqs' is not allowed with amino-acid data.");
5736
sscanf(buff,"%lf,%lf,%lf,%lf",&A,&C,&G,&T);
5737
mod->user_b_freq->v[0] = (phydbl)A;
5738
mod->user_b_freq->v[1] = (phydbl)C;
5739
mod->user_b_freq->v[2] = (phydbl)G;
5740
mod->user_b_freq->v[3] = (phydbl)T;
5741
mod->s_opt->user_state_freq = YES;
5742
mod->s_opt->opt_state_freq = NO;
5747
buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
5748
if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
5750
mod->s_opt->opt_efrq_weight = YES;
5754
mod->s_opt->opt_efrq_weight = NO;
5759
//////////////////////////////////////////////////////////////
5760
//////////////////////////////////////////////////////////////
5762
void Make_Topology_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5767
init_tree = XML_Get_Attribute_Value(instance,"init.tree");
5771
PhyML_Printf("\n== Attribute 'init.tree=bionj|user|random' is mandatory");
5772
PhyML_Printf("\n== Please amend your XML file accordingly.");
5776
if(!strcmp(init_tree,"user") ||
5777
!strcmp(init_tree,"User"))
5779
char *starting_tree = NULL;
5780
starting_tree = XML_Get_Attribute_Value(instance,"file.name");
5782
if(!Filexists(starting_tree))
5784
PhyML_Printf("\n== The tree file '%s' could not be found.",starting_tree);
5789
strcpy(io->in_tree_file,starting_tree);
5791
io->fp_in_tree = Openfile(io->in_tree_file,0);
5794
else if(!strcmp(init_tree,"random") ||
5795
!strcmp(init_tree,"Random"))
5797
char *n_rand_starts = NULL;
5799
io->mod->s_opt->random_input_tree = YES;
5801
n_rand_starts = XML_Get_Attribute_Value(instance,"n.rand.starts");
5805
mod->s_opt->n_rand_starts = atoi(n_rand_starts);
5806
if(mod->s_opt->n_rand_starts < 1) Exit("\n== Number of random starting trees must be > 0.\n\n");
5809
strcpy(io->out_trees_file,io->in_align_file);
5810
strcat(io->out_trees_file,"_phyml_rand_trees");
5811
if(io->append_run_ID) { strcat(io->out_trees_file,"_"); strcat(io->out_trees_file,io->run_id_string); }
5812
io->fp_out_trees = Openfile(io->out_trees_file,1);
5814
else if(!strcmp(init_tree,"parsimony") ||
5815
!strcmp(init_tree,"Parsimony"))
5820
// Estimate tree topology
5821
char *optimise = NULL;
5823
optimise = XML_Get_Attribute_Value(instance,"optimise.tree");
5829
select = XML_Validate_Attr_Int(optimise,6,
5833
if(select > 2) io->mod->s_opt->opt_topo = NO;
5839
search = XML_Get_Attribute_Value(instance,"search");
5840
select = XML_Validate_Attr_Int(search,4,"spr","nni","best","none");
5846
io->mod->s_opt->topo_search = SPR_MOVE;
5847
io->mod->s_opt->opt_topo = YES;
5852
io->mod->s_opt->topo_search = NNI_MOVE;
5853
io->mod->s_opt->opt_topo = YES;
5858
io->mod->s_opt->topo_search = BEST_OF_NNI_AND_SPR;
5859
io->mod->s_opt->opt_topo = YES;
5864
io->mod->s_opt->opt_topo = NO;
5869
PhyML_Printf("\n== Topology search option '%s' is not valid.",search);
5879
//////////////////////////////////////////////////////////////
5880
//////////////////////////////////////////////////////////////
5882
void Make_RAS_From_XML_Node(xml_node *parent, t_mod *mod)
5888
mod->ras->n_catg = 0;
5890
XML_Check_Siterates_Node(parent);
5892
w = XML_Search_Node_Name("weights",YES,parent);
5895
family = XML_Get_Attribute_Value(w,"family");
5896
select = XML_Validate_Attr_Int(family,3,"gamma","gamma+inv","freerates");
5899
case 0: // Gamma model
5901
char *alpha,*alpha_opt;
5903
mod->s_opt->opt_pinvar = NO;
5904
mod->ras->invar = NO;
5906
alpha = XML_Get_Attribute_Value(w,"alpha");
5910
if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") ||
5911
!strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5913
mod->s_opt->opt_alpha = YES;
5917
mod->s_opt->opt_alpha = NO;
5918
mod->ras->alpha->v = String_To_Dbl(alpha);
5922
alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5926
if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5928
mod->s_opt->opt_alpha = YES;
5932
mod->s_opt->opt_alpha = NO;
5936
mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
5938
Make_RAS_Complete(mod->ras);
5942
case 1: // Gamma+Inv model
5944
char *alpha,*pinv,*alpha_opt,*pinv_opt;
5946
mod->ras->invar = YES;
5947
mod->s_opt->opt_pinvar = YES;
5949
alpha = XML_Get_Attribute_Value(w,"alpha");
5953
if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") ||
5954
!strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5956
mod->s_opt->opt_alpha = YES;
5960
mod->s_opt->opt_alpha = NO;
5961
mod->ras->alpha->v = String_To_Dbl(alpha);;
5965
alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5969
if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5971
mod->s_opt->opt_alpha = YES;
5975
mod->s_opt->opt_alpha = NO;
5980
pinv = XML_Get_Attribute_Value(w,"pinv");
5984
if(!strcmp(pinv,"estimate") || !strcmp(pinv,"estimated") ||
5985
!strcmp(pinv,"optimise") || !strcmp(pinv,"optimised"))
5987
mod->s_opt->opt_pinvar = YES;
5991
mod->s_opt->opt_pinvar = NO;
5992
mod->ras->pinvar->v = String_To_Dbl(pinv);;
5996
pinv_opt = XML_Get_Attribute_Value(w,"optimise.pinv");
6000
if(!strcmp(pinv_opt,"yes") || !strcmp(pinv_opt,"true"))
6002
mod->s_opt->opt_pinvar = YES;
6006
mod->s_opt->opt_pinvar = NO;
6010
mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
6013
case 2: // FreeRate model
6015
char *opt_freerates;
6017
mod->ras->free_mixt_rates = YES;
6019
mod->s_opt->opt_free_mixt_rates = YES;
6021
opt_freerates = XML_Get_Attribute_Value(w,"optimise.freerates");
6025
if(!strcmp(opt_freerates,"yes") || !strcmp(opt_freerates,"true"))
6027
mod->s_opt->opt_free_mixt_rates = YES;
6031
mod->s_opt->opt_free_mixt_rates = NO;
6035
mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
6040
PhyML_Printf("\n== family: %s",family);
6041
PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
6047
int nc = XML_Get_Number_Of_Classes_Siterates(parent);
6049
if(w && nc != mod->ras->n_catg)
6051
PhyML_Printf("\n== <siterates> component '%s'. The number of classes ",parent->id);
6052
PhyML_Printf("\n== specified in the <weight> component should be ");
6053
PhyML_Printf("\n== the same as that of <instance> nodes. Please fix");
6054
PhyML_Printf("\n== your XML file accordingly.");
6058
if(!w) mod->ras->n_catg = nc;
6060
Make_RAS_Complete(mod->ras);
6064
//////////////////////////////////////////////////////////////
6065
//////////////////////////////////////////////////////////////
6066
//////////////////////////////////////////////////////////////
6067
//////////////////////////////////////////////////////////////
6068
//////////////////////////////////////////////////////////////
6069
//////////////////////////////////////////////////////////////
added
removed
GDE/PHYML20130708/phyml/src/io.c
