source: branches/nameserver/GDE/PHYML20130708/phyml/src/io.c

Last change on this file was 10307, checked in by aboeckma, 12 years ago

added most recent version of phyml
File size: 169.1 KB
Line 
1/*
2
3PhyML:  a program that  computes maximum likelihood phylogenies from
4DNA or AA homologous sequences.
5
6Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8All parts of the source except where indicated are distributed under
9the GNU public licence. See http://www.opensource.org for details.
10
11*/
12
13#include "io.h"
14
15
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)
19*/
20t_tree *Read_Tree(char **s_tree)
21{
22  char **subs;
23  int i,n_ext,n_int,n_otu;
24  t_tree *tree;
25  int degree,len;
26  t_node *root_node;
27 
28  n_int = n_ext = 0;
29 
30  n_otu=0;
31  For(i,(int)strlen((*s_tree))) if((*s_tree)[i] == ',') n_otu++;
32  n_otu+=1;
33
34  tree = Make_Tree_From_Scratch(n_otu,NULL);
35  subs = Sub_Trees((*s_tree),&degree);
36  Clean_Multifurcation(subs,degree,3);
37
38  if(degree == 2) 
39    {
40      /* Unroot_Tree(subs); */
41      /* degree = 3; */
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;
46      n_int         -= 1;
47    }
48  else
49    {     
50      root_node      = tree->a_nodes[n_otu];
51      root_node->num = n_otu;
52      tree->n_root   = NULL;
53   }
54 
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 */
58    {
59      degree = 3;
60      Free((*s_tree));
61      len = 0;
62      For(i,degree) len += (strlen(subs[i])+1);
63      len += 5;
64
65      (*s_tree) = (char *)mCalloc(len,sizeof(char));
66
67      (*s_tree)[0] = '('; (*s_tree)[1] = '\0';
68      For(i,degree) 
69        {
70          strcat((*s_tree),subs[i]);
71          strcat((*s_tree),",\0");
72        }
73
74      sprintf((*s_tree)+strlen((*s_tree))-1,"%s",");\0");
75     
76      For(i,NODE_DEG_MAX) Free(subs[i]);
77      Free(subs);
78
79      subs = Sub_Trees((*s_tree),&degree);
80    }
81
82  root_node->tax = 0;
83
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);
87
88  for(i=degree;i<NODE_DEG_MAX;i++) Free(subs[i]);
89  Free(subs);
90
91  if(tree->n_root)
92    {
93      tree->e_root = tree->a_edges[tree->num_curr_branch_available];
94           
95      tree->n_root->v[2] = tree->n_root->v[0];
96      tree->n_root->v[0] = NULL;
97
98      tree->n_root->l[2] = tree->n_root->l[0];
99
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; }
102
103      Connect_One_Edge_To_Two_Nodes(tree->n_root->v[2],
104                                    tree->n_root->v[1],
105                                    tree->e_root,
106                                    tree);
107
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;
111      else
112        tree->n_root_pos = .5;
113    }
114 
115  return tree;
116}
117
118//////////////////////////////////////////////////////////////
119//////////////////////////////////////////////////////////////
120
121/* 'a' in t_node a stands for ancestor. 'd' stands for descendant */ 
122void R_rtree(char *s_tree_a, char *s_tree_d, t_node *a, t_tree *tree, int *n_int, int *n_ext)
123{
124  int i;
125  t_node *d;
126  int n_otu = tree->n_otu;
127
128  if(strstr(s_tree_a," ")) 
129    {
130      PhyML_Printf("\n== [%s]",s_tree_a);
131      Warn_And_Exit("\n== Err: the tree must not contain a ' ' character\n");
132    }
133
134  if(s_tree_d[0] == '(')
135    {
136      char **subs;
137      int degree;
138      t_edge *b;
139
140      (*n_int)+=1;
141
142      if((*n_int + n_otu) == (2*n_otu-1))
143        {
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__);
147          Exit("\n");     
148        }
149
150      d      = tree->a_nodes[n_otu+*n_int];
151      d->num = n_otu+*n_int;
152      d->tax = 0;
153      b      = tree->a_edges[tree->num_curr_branch_available];
154
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);     
157
158      For(i,3)
159        {
160          if(!a->v[i])
161            {
162              a->v[i]=d;
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;
165              break;
166            }
167        }
168      d->v[0]=a;
169
170      if(a != tree->n_root)
171        {
172          Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[tree->num_curr_branch_available],tree);
173        }
174     
175      subs=Sub_Trees(s_tree_d,&degree);
176
177      if(degree > 2)
178        {
179          Clean_Multifurcation(subs,degree,2);
180
181          Free(s_tree_d);
182
183          s_tree_d = (char *)mCalloc(strlen(subs[0])+strlen(subs[1])+5,sizeof(char));
184
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,")");
190          For(i,b->n_labels)
191            {
192              strcat(s_tree_d,"#");
193              strcat(s_tree_d,b->labels[i]);
194            }
195         
196          For(i,NODE_DEG_MAX) Free(subs[i]);
197          Free(subs);
198
199          subs=Sub_Trees(s_tree_d,&degree);
200        }
201
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);
204 
205      for(i=2;i<NODE_DEG_MAX;i++) Free(subs[i]);
206      Free(subs);
207    }
208
209  else
210    {
211      int i;
212
213      d      = tree->a_nodes[*n_ext];
214      d->tax = 1;
215
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);
219     
220      For(i,3)
221        {
222         if(!a->v[i])
223           {
224             a->v[i]=d;
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;
227             break;
228           }
229        }
230      d->v[0]=a;
231
232      if(a != tree->n_root)
233        {
234          Connect_One_Edge_To_Two_Nodes(a,d,tree->a_edges[tree->num_curr_branch_available],tree);
235        }
236
237      d->num=*n_ext;
238      (*n_ext)+=1;
239    }
240 
241  Free(s_tree_d);
242
243}
244
245//////////////////////////////////////////////////////////////
246//////////////////////////////////////////////////////////////
247
248
249void Read_Branch_Label(char *s_d, char *s_a, t_edge *b)
250{
251  char *sub_tp;
252  char *p;
253  int posp,posl;
254
255  sub_tp = (char *)mCalloc(3+(int)strlen(s_d)+1,sizeof(char));
256  /* sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
257
258  sub_tp[0] = '(';
259  sub_tp[1] = '\0';
260  strcat(sub_tp,s_d);
261  strcat(sub_tp,"#");
262  p = strstr(s_a,sub_tp);
263 
264  if(!p)
265    {
266      sub_tp[0] = ',';
267      sub_tp[1] = '\0';
268      strcat(sub_tp,s_d);
269      strcat(sub_tp,"#");
270      p = strstr(s_a,sub_tp);
271    }
272
273  b->n_labels = 0;
274  if(p)
275    {
276      if(!(b->n_labels%BLOCK_LABELS)) Make_New_Edge_Label(b);
277      b->n_labels++;
278
279      posp = strlen(s_d);
280      while(p[posp] != '#') posp++;
281      posp++;
282
283      posl = 0;
284      do 
285        {
286          b->labels[b->n_labels-1][posl] = p[posp];
287          posl++;
288          posp++;
289          if(p[posp] == '#') 
290            { 
291              b->labels[b->n_labels-1][posl] = '\0';
292              b->n_labels++;
293              if(!(b->n_labels%BLOCK_LABELS)) Make_New_Edge_Label(b);
294              posp++;
295              posl=0;
296            }
297        }
298      while((p[posp] != ':') &&
299            (p[posp] != ',') &&
300            (p[posp] != '('));
301
302      b->labels[b->n_labels-1][posl] = '\0';
303    }
304
305  if(p)
306    {
307      /* if(b->n_labels == 1) */
308      /*        PhyML_Printf("\n. Read label '%s' on t_edge %3d.",b->labels[0],b->num); */
309      /* else */
310      /*        { */
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); */
314      /*        } */
315
316      if(!strcmp(b->labels[0],"NULL"))
317        {
318          b->does_exist = NO;
319        }
320    }
321  /* else */
322  /*   { */
323  /*     PhyML_Printf("\n. No label found on %s",s_d); */
324  /*   } */
325  Free(sub_tp);
326}
327
328//////////////////////////////////////////////////////////////
329//////////////////////////////////////////////////////////////
330
331void Read_Branch_Length(char *s_d, char *s_a, t_tree *tree)
332{
333  char *sub_tp;
334  char *p;
335  t_edge *b;
336  int i;
337
338  b = tree->a_edges[tree->num_curr_branch_available];
339
340  /* sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
341  sub_tp = (char *)mCalloc(10+strlen(s_d)+1,sizeof(char));
342
343  For(i,b->n_labels)
344    {
345      strcat(s_d,"#");
346      strcat(s_d,b->labels[i]);
347    }
348
349  sub_tp[0] = '(';
350  sub_tp[1] = '\0';
351  strcat(sub_tp,s_d);
352  strcat(sub_tp,":");
353  p = strstr(s_a,sub_tp);
354
355  if(!p)
356    {
357      sub_tp[0] = ',';
358      sub_tp[1] = '\0';
359      strcat(sub_tp,s_d);
360      strcat(sub_tp,":");
361      p = strstr(s_a,sub_tp);
362    }
363
364
365  if(p)
366    {
367      b->l->v = atof((char *)p+(int)strlen(sub_tp));
368      tree->has_branch_lengths = YES;
369      b->does_exist = YES;
370    }     
371  else
372    {
373      b->l->v = -1.;
374    }
375
376
377  Free(sub_tp);
378}
379
380//////////////////////////////////////////////////////////////
381//////////////////////////////////////////////////////////////
382
383void Read_Node_Name(t_node *d, char *s_tree_d, t_tree *tree)
384{
385  int i;
386 
387  if(!tree->a_edges[tree->num_curr_branch_available]->n_labels)
388    {
389      d->name = (char *)mCalloc(strlen(s_tree_d)+1,sizeof(char ));
390      strcpy(d->name,s_tree_d);
391    }
392  else
393    {
394      i = 0;
395      do
396        {
397          d->name = (char *)realloc(d->name,(i+1)*sizeof(char ));
398          d->name[i] = s_tree_d[i];
399          i++;
400        }
401      while(s_tree_d[i] != '#');
402      d->name[i] = '\0';
403    }
404  d->ori_name = d->name;
405
406}
407//////////////////////////////////////////////////////////////
408//////////////////////////////////////////////////////////////
409
410
411void Clean_Multifurcation(char **subtrees, int current_deg, int end_deg)
412{
413
414  if(current_deg <= end_deg) return;
415  else
416    {
417      char *s_tmp;
418      int i;
419
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,
424                              sizeof(char));
425     
426      strcat(s_tmp,"(\0");
427      strcat(s_tmp,subtrees[0]);
428      strcat(s_tmp,",\0");
429      strcat(s_tmp,subtrees[1]);
430      strcat(s_tmp,")#NULL\0"); /* Add the label 'NULL' to identify a non-existing edge */
431      Free(subtrees[0]);
432      subtrees[0] = s_tmp;
433
434      for(i=1;i<current_deg-1;i++) strcpy(subtrees[i],subtrees[i+1]);
435
436      Clean_Multifurcation(subtrees,current_deg-1,end_deg);
437    }
438}
439
440//////////////////////////////////////////////////////////////
441//////////////////////////////////////////////////////////////
442
443
444char **Sub_Trees(char *tree, int *degree)
445{
446  char **subs;
447  int posbeg,posend;
448  int i;
449
450  if(tree[0] != '(') {*degree = 1; return NULL;}
451
452  subs=(char **)mCalloc(NODE_DEG_MAX,sizeof(char *));
453
454  For(i,NODE_DEG_MAX) subs[i]=(char *)mCalloc(strlen(tree)+1,sizeof(char));
455
456
457  posbeg=posend=1;
458  (*degree)=0;
459  do
460    {
461      posbeg = posend;
462      if(tree[posend] != '(')
463        {
464          while((tree[posend] != ',' ) &&
465                (tree[posend] != ':' ) &&
466                (tree[posend] != '#' ) &&
467                (tree[posend] != ')' )) 
468            {
469              posend++ ;
470            }
471          posend -= 1;
472        }
473      else posend=Next_Par(tree,posend);
474
475      while((tree[posend+1] != ',') &&
476            (tree[posend+1] != ':') &&
477            (tree[posend+1] != '#') &&
478            (tree[posend+1] != ')')) {posend++;}
479
480
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 */
484
485      posend += 1;
486      while((tree[posend] != ',') &&
487            (tree[posend] != ')')) {posend++;}
488      posend+=1;
489
490
491      (*degree)++;
492      if((*degree) == NODE_DEG_MAX)
493        {
494          For(i,(*degree))
495            PhyML_Printf("\n. Subtree %d : %s\n",i+1,subs[i]);
496
497          PhyML_Printf("\n. The degree of a t_node cannot be greater than %d\n",NODE_DEG_MAX);
498          Warn_And_Exit("\n");
499        }
500    }
501  while(tree[posend-1] != ')');
502
503  return subs;
504}
505
506
507//////////////////////////////////////////////////////////////
508//////////////////////////////////////////////////////////////
509
510
511int Next_Par(char *s, int pos)
512{
513  int curr;
514
515  curr=pos+1;
516
517  while(*(s+curr) != ')')
518    {
519      if(*(s+curr) == '(') curr=Next_Par(s,curr);
520      curr++;
521    }
522
523  return curr;
524}
525
526//////////////////////////////////////////////////////////////
527//////////////////////////////////////////////////////////////
528
529
530void Print_Tree(FILE *fp, t_tree *tree)
531{
532  char *s_tree;
533  int i;
534
535  s_tree = (char *)Write_Tree(tree,NO);
536
537  if(OUTPUT_TREE_FORMAT == NEWICK) PhyML_Fprintf(fp,"%s\n",s_tree);
538  else if(OUTPUT_TREE_FORMAT == NEXUS)
539    {
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;");
547    }
548  Free(s_tree);
549}
550
551//////////////////////////////////////////////////////////////
552//////////////////////////////////////////////////////////////
553
554
555char *Write_Tree(t_tree *tree, int custom)
556{
557  char *s;
558  int i,available;
559  int init_len;
560
561  init_len = 3*(int)T_MAX_NAME;
562
563#ifndef MPI
564  s=(char *)mCalloc(init_len,sizeof(char));
565  available = init_len;
566#elif defined MPI
567  s=(char *)mCalloc(T_MAX_LINE,sizeof(char));
568#endif
569
570 
571  s[0]='(';
572
573  if(custom == NO)
574    {
575      if(!tree->n_root)
576        {
577          i = 0;
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++;
581         
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);
585        }
586      else
587        {
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);
590        }
591    }
592  else
593    {
594      int pos;
595
596      pos = 1;
597
598      if(!tree->n_root)
599        {
600          i = 0;
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++;
604         
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);
608        }
609      else
610        {
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);
613        }
614
615    }
616
617  s[(int)strlen(s)-1]=')';
618  s[(int)strlen(s)]=';';
619
620  return s;     
621}
622
623//////////////////////////////////////////////////////////////
624//////////////////////////////////////////////////////////////
625
626
627void R_wtree(t_node *pere, t_node *fils, int *available, char **s_tree, t_tree *tree)
628{
629  int i,p;
630  char *format;
631  int last_len;
632#if !(defined PHYTIME || defined SERGEII)
633  phydbl mean_len;
634#endif
635
636  format = (char *)mCalloc(100,sizeof(char));
637
638  sprintf(format,"%%.%df",tree->bl_ndigits);
639
640  p = -1;
641  if(fils->tax)
642    {
643      /* printf("\n- Writing on %p",*s_tree); */
644      /* ori_len = *pos; */
645
646      last_len = (int)strlen(*s_tree);
647
648      if(OUTPUT_TREE_FORMAT == NEWICK)
649        {
650          if(tree->write_tax_names == YES)
651            {
652              if(tree->io && tree->io->long_tax_names) 
653                {
654                  strcat(*s_tree,tree->io->long_tax_names[fils->num]);
655                }
656              else
657                {
658                  strcat(*s_tree,fils->name);
659                }
660            }
661          else if(tree->write_tax_names == NO)
662            {
663              sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
664            }
665        }
666      else if(OUTPUT_TREE_FORMAT == NEXUS)
667        {
668          sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->num+1);
669        }
670      else
671        {
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);
675        }
676
677      if((fils->b) && (fils->b[0]) && (tree->write_br_lens == YES))
678        {
679          (*s_tree)[(int)strlen(*s_tree)] = ':';
680
681#if !(defined PHYTIME || defined SERGEII)
682          if(!tree->n_root)
683            {
684              if(tree->is_mixt_tree == NO) 
685                {
686                  mean_len = fils->b[0]->l->v;
687                }
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));
690            }
691          else
692            {
693              if(pere == tree->n_root)
694                {
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) 
697                    {
698                      mean_len = tree->e_root->l->v;
699                    }
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);
702                }
703              else
704                {
705                  if(tree->is_mixt_tree == NO) 
706                    {
707                      mean_len = fils->b[0]->l->v;
708                    }
709
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));
712                }
713            }
714#else
715          if(!tree->n_root)
716            {
717              sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,fils->b[0]->l->v));
718            }
719          else
720            {
721              sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,tree->rates->cur_l[fils->num]));
722            }
723#endif
724        }
725
726      /* strcat(*s_tree,","); */
727      (*s_tree)[(int)strlen(*s_tree)] = ',';
728
729
730#ifndef MPI     
731      (*available) -= ((int)strlen(*s_tree) - last_len);
732
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); */
735
736      if(*available < 0)
737        {
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__);
742          Warn_And_Exit("");
743        }
744
745      if(*available < (int)T_MAX_NAME)
746        {
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)); */
751        }
752#endif
753
754    }
755  else
756    {
757      (*s_tree)[(int)strlen(*s_tree)]='(';
758
759#ifndef MPI
760      (*available) -= 1;
761
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); */
764
765      if(*available < (int)T_MAX_NAME)
766        {
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)); */
771        }
772#endif
773      /* (*available)--; */
774
775      /* if(*available < (int)T_MAX_NAME/2) */
776      /*        { */
777      /*          (*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char)); */
778      /*          (*available) = (int)T_MAX_NAME; */
779      /*        } */
780
781
782      if(tree->n_root)
783        {
784          For(i,3)
785            {
786              if((fils->v[i] != pere) && (fils->b[i] != tree->e_root))
787                R_wtree(fils,fils->v[i],available,s_tree,tree);
788              else p=i;
789            }
790        }
791      else
792        {
793          For(i,3)
794            {
795              if(fils->v[i] != pere)
796                R_wtree(fils,fils->v[i],available,s_tree,tree);
797              else p=i;
798            }
799        }
800     
801      if(p < 0)
802        {
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__);
806          Exit("\n");
807        }
808
809      last_len = (int)strlen(*s_tree);
810
811      (*s_tree)[last_len-1] = ')';
812
813      if((fils->b) && (tree->write_br_lens == YES))
814        {
815          if(tree->print_boot_val)
816            {
817              sprintf(*s_tree+(int)strlen(*s_tree),"%d",fils->b[p]->bip_score);
818            }
819          else if(tree->print_alrt_val)
820            {
821              sprintf(*s_tree+(int)strlen(*s_tree),"%f",fils->b[p]->ratio_test);
822            }
823         
824          fflush(NULL);
825
826          (*s_tree)[(int)strlen(*s_tree)] = ':';
827
828#if !(defined PHYTIME || defined SERGEII)
829          if(!tree->n_root)
830            {
831              if(tree->is_mixt_tree == NO) 
832                {
833                  mean_len = fils->b[p]->l->v;
834                }
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));
837            }
838          else
839            {
840              if(pere == tree->n_root)
841                {
842                  phydbl root_pos = (fils == tree->n_root->v[2])?(tree->n_root_pos):(1.-tree->n_root_pos);
843
844                  if(tree->is_mixt_tree == NO) 
845                    {
846                      mean_len = tree->e_root->l->v;
847                    }
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);
850                }
851              else
852                {
853                  if(tree->is_mixt_tree == NO) 
854                    {
855                      mean_len = fils->b[p]->l->v;
856                    }
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));
859                }
860            }
861#else
862          if(!tree->n_root)
863            {
864              sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,fils->b[p]->l->v));
865            }
866          else
867            {
868              sprintf(*s_tree+(int)strlen(*s_tree),format,MAX(0.0,tree->rates->cur_l[fils->num]));
869            }
870#endif   
871        }
872
873      /* strcat(*s_tree,","); */
874      (*s_tree)[(int)strlen(*s_tree)] = ',';
875     
876#ifndef MPI
877      (*available) -= ((int)strlen(*s_tree) - last_len);
878     
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); */
881
882      if(*available < 0)
883        {
884          PhyML_Printf("\n== available = %d",*available);
885          PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
886          Warn_And_Exit("");
887        }
888
889      if(*available < (int)T_MAX_NAME)
890        {
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)); */
895        }
896#endif
897
898      /* if(*available < (int)T_MAX_NAME/2) */
899      /*        { */
900      /*          (*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char)); */
901      /*          (*available) = (int)T_MAX_NAME; */
902      /*        } */
903    }
904
905  Free(format);
906}
907
908//////////////////////////////////////////////////////////////
909//////////////////////////////////////////////////////////////
910
911
912void R_wtree_Custom(t_node *pere, t_node *fils, int *available, char **s_tree, int *pos, t_tree *tree)
913{
914  int i,p,ori_len;
915  char *format;
916
917  format = (char *)mCalloc(100,sizeof(char));
918
919  sprintf(format,"%%.%df",tree->bl_ndigits);
920  /* strcpy(format,"%f"); */
921
922  p = -1;
923  if(fils->tax)
924    {
925/*       printf("\n- Writing on %p",*s_tree); */
926      ori_len = *pos;
927
928      if(OUTPUT_TREE_FORMAT == NEWICK)
929        {
930          if(tree->write_tax_names == YES)
931            {
932              if(tree->io && tree->io->long_tax_names) 
933                {
934                  strcat(*s_tree,tree->io->long_tax_names[fils->num]);
935                  (*pos) += (int)strlen(tree->io->long_tax_names[fils->num]);
936                }
937              else
938                {
939                  strcat(*s_tree,fils->name);
940                  (*pos) += (int)strlen(fils->name);
941                }         
942            }
943          else if(tree->write_tax_names == NO)
944            {
945              (*pos) += sprintf(*s_tree+*pos,"%d",fils->num);
946            }
947        }
948      else if(OUTPUT_TREE_FORMAT == NEXUS)
949        {
950          (*pos) += sprintf(*s_tree+*pos,"%d",fils->num+1);
951        }
952      else
953        {
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);
957        }
958
959      if((fils->b) && (fils->b[0]) && (tree->write_br_lens == YES))
960        {
961          strcat(*s_tree,":");
962          (*pos)++;
963 
964#if !(defined PHYTIME || defined SERGEII)
965          if(!tree->n_root)
966            {
967              (*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
968            }
969          else
970            {
971              if(pere == tree->n_root)
972                {
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);
975                }
976              else
977                {
978                  (*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
979                }
980            }           
981#else
982          if(!tree->n_root)
983            {
984              (*pos) += sprintf(*s_tree+*pos,format,fils->b[0]->l->v);
985            }
986          else
987            {
988              (*pos) += sprintf(*s_tree+*pos,format,tree->rates->cur_l[fils->num]);
989            }
990#endif
991        }
992     
993      if(tree->write_labels)
994        {
995          if(fils->b[0]->n_labels < 10)
996            For(i,fils->b[0]->n_labels) 
997              {
998                (*pos) += sprintf(*s_tree+*pos,"::%s",fils->b[0]->labels[i]);
999              }
1000          else
1001            {
1002              (*pos) += sprintf(*s_tree+*pos,"::%d_labels",fils->b[0]->n_labels);
1003            }
1004        }
1005 
1006
1007      strcat(*s_tree,",");
1008      (*pos)++;
1009
1010      (*available) = (*available) - (*pos - ori_len);
1011
1012      if(*available < 0)
1013        {
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__);
1018          Warn_And_Exit("");
1019        }
1020
1021      if(*available < (int)T_MAX_NAME)
1022        {
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;
1026        }
1027/*       printf(" %s [%d,%d]",*s_tree,(int)strlen(*s_tree),*available); */
1028    }
1029  else
1030    {
1031
1032      (*s_tree)[(*pos)]='(';
1033      (*s_tree)[(*pos)+1]='\0';
1034      (*pos)++;
1035      (*available)--;
1036
1037      if(*available < (int)T_MAX_NAME/2)
1038        {
1039          (*s_tree) = (char *)mRealloc(*s_tree,*pos+(int)T_MAX_NAME,sizeof(char));
1040          (*available) = (int)T_MAX_NAME;
1041        }
1042
1043      if(tree->n_root)
1044        {
1045          For(i,3)
1046            {
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);
1049              else p=i;
1050            }
1051        }
1052      else
1053        {
1054          For(i,3)
1055            {
1056              if(fils->v[i] != pere)
1057                R_wtree_Custom(fils,fils->v[i],available,s_tree,pos,tree);
1058              else p=i;
1059            }
1060        }
1061 
1062      ori_len = *pos;
1063     
1064      if(p < 0)
1065        {
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__);
1069          Warn_And_Exit("");
1070        }
1071
1072/*       printf("\n+ Writing on %p",*s_tree); */
1073      (*s_tree)[(*pos)-1] = ')';
1074      (*s_tree)[(*pos)]   = '\0';
1075
1076      if((fils->b) && (tree->write_br_lens == YES))
1077        {
1078          if(tree->print_boot_val)
1079            {
1080              (*pos) += sprintf(*s_tree+*pos,"%d",fils->b[p]->bip_score);
1081            }
1082          else if(tree->print_alrt_val)
1083            {
1084              (*pos) += sprintf(*s_tree+*pos,"%f",fils->b[p]->ratio_test);
1085            }
1086         
1087          fflush(NULL);
1088
1089          strcat(*s_tree,":");
1090          (*pos)++;
1091
1092#if !(defined PHYTIME || defined SERGEII)
1093          if(!tree->n_root)
1094            {
1095              (*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1096            }
1097          else
1098            {
1099              if(pere == tree->n_root)
1100                {
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);
1103                }
1104              else
1105                {
1106                  (*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1107                }
1108            }
1109#else
1110          if(!tree->n_root)
1111            {
1112              (*pos) += sprintf(*s_tree+*pos,format,fils->b[p]->l->v);
1113            }
1114          else
1115            {
1116              (*pos) += sprintf(*s_tree+*pos,format,tree->rates->cur_l[fils->num]);
1117            }
1118#endif
1119         
1120        }
1121
1122      if((tree->write_labels) && (fils->b[p]->labels != NULL))
1123        {
1124          if(fils->b[p]->n_labels < 10)
1125            For(i,fils->b[p]->n_labels) 
1126              {
1127                (*pos) += sprintf(*s_tree+*pos,"::%s",fils->b[p]->labels[i]);
1128              }
1129          else
1130            {
1131              (*pos) += sprintf(*s_tree+*pos,"::%d_labels",fils->b[p]->n_labels);
1132            }
1133        }
1134
1135      strcat(*s_tree,",");
1136      (*pos)++;
1137      (*available) = (*available) - (*pos - ori_len);
1138
1139      if(*available < 0)
1140        {
1141          PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1142          Warn_And_Exit("");
1143        }
1144
1145      if(*available < (int)T_MAX_NAME)
1146        {
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;
1150        }
1151/*       printf(" %s [%d,%d]",*s_tree,(int)strlen(*s_tree),*available); */
1152    }
1153
1154  Free(format);
1155}
1156
1157//////////////////////////////////////////////////////////////
1158//////////////////////////////////////////////////////////////
1159
1160void Detect_Align_File_Format(option *io)
1161{
1162  int c;
1163  fpos_t curr_pos;
1164 
1165  fgetpos(io->fp_in_align,&curr_pos);
1166 
1167  errno = 0;
1168
1169  while((c=fgetc(io->fp_in_align)) != EOF)
1170    {
1171      if(errno) io->data_file_format = PHYLIP;
1172      else if(c == '#')
1173        {
1174          char s[10],t[6]="NEXUS";
1175          if(!fgets(s,6,io->fp_in_align))
1176            {     
1177              PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1178              Exit("\n");
1179            }
1180          if(!strcmp(t,s)) 
1181            {
1182              fsetpos(io->fp_in_align,&curr_pos);
1183              io->data_file_format = NEXUS;
1184              return;
1185            }
1186        }
1187    }
1188 
1189  fsetpos(io->fp_in_align,&curr_pos);
1190}
1191
1192//////////////////////////////////////////////////////////////
1193//////////////////////////////////////////////////////////////
1194
1195void Detect_Tree_File_Format(option *io)
1196{
1197  int c;
1198  fpos_t curr_pos;
1199 
1200  fgetpos(io->fp_in_tree,&curr_pos);
1201
1202  errno = 0;
1203
1204  while((c=fgetc(io->fp_in_tree)) != EOF)
1205    {
1206      if(errno) 
1207        {
1208          io->tree_file_format = PHYLIP;
1209          PhyML_Printf("\n. Detected PHYLIP tree file format.");
1210        }
1211      else if(c == '#')
1212        {
1213          char s[10],t[6]="NEXUS";
1214          if(!fgets(s,6,io->fp_in_tree))
1215            {
1216              PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1217              Warn_And_Exit("");
1218            }
1219          if(!strcmp(t,s))
1220            {
1221              fsetpos(io->fp_in_tree,&curr_pos);
1222              io->tree_file_format = NEXUS;
1223              PhyML_Printf("\n. Detected NEXUS tree file format.");
1224              return;
1225            }
1226        }
1227    }
1228 
1229  fsetpos(io->fp_in_tree,&curr_pos);
1230}
1231
1232//////////////////////////////////////////////////////////////
1233//////////////////////////////////////////////////////////////
1234
1235align **Get_Seq(option *io)
1236{
1237  io->data = NULL;
1238
1239  if(!io->fp_in_align)
1240    {
1241      PhyML_Printf("\n== Filehandle to '%s' seems to be closed.",io->in_align_file);
1242      Exit("\n");
1243    }
1244
1245  Detect_Align_File_Format(io);
1246
1247  switch(io->data_file_format)
1248    {
1249    case PHYLIP: 
1250      {
1251        io->data = Get_Seq_Phylip(io);
1252        break;
1253      }
1254    case NEXUS:
1255      {
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);
1259        Free_Nexus(io);
1260        break;
1261      }
1262    default:
1263      {
1264        PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1265        Exit("\n");
1266        break;
1267      }
1268    }
1269
1270  if(!io->data)
1271    {
1272      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1273      Exit("\n");
1274    }
1275  else
1276    {
1277      Post_Process_Data(io);
1278    }
1279
1280
1281  return io->data;
1282}
1283
1284//////////////////////////////////////////////////////////////
1285//////////////////////////////////////////////////////////////
1286
1287void Post_Process_Data(option *io)
1288{     
1289  int i,j;
1290 
1291  For(i,io->data[0]->len)
1292    {
1293      For(j,io->n_otu)
1294        {
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';
1298        }
1299    }
1300 
1301  For(i,io->n_otu) io->data[i]->len = io->data[0]->len;
1302}
1303
1304/* align **Get_Seq_Nexus(option *io) */
1305/* { */
1306/*   char *s,*ori_s; */
1307/*   char *token; */
1308/*   int in_comment; */
1309/*   nexcom *curr_com; */
1310/*   nexparm *curr_parm; */
1311/*   int nxt_token_t,cur_token_t; */
1312
1313/*   s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); */
1314/*   token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char)); */
1315         
1316/*   ori_s      = s; */
1317/*   in_comment = NO; */
1318/*   curr_com   = NULL; */
1319/*   curr_parm  = NULL; */
1320/*   nxt_token_t = NEXUS_COM;  */
1321/*   cur_token_t = -1;  */
1322
1323/*   while(fgets(s,T_MAX_LINE,io->fp_in_align)) */
1324/*     {       */
1325/*       do */
1326/*      {          */
1327/*        Get_Token(&s,token);     */
1328
1329/* /\*    PhyML_Printf("\n. Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); *\/ */
1330
1331/*        if(token[0] == '\0') break; */
1332
1333/*        if(token[0] == ';')  */
1334/*          { */
1335/*            curr_com   = NULL; */
1336/*            curr_parm  = NULL; */
1337/*            nxt_token_t = NEXUS_COM; */
1338/*            cur_token_t = -1; */
1339/*            break; /\* End of command *\/  */
1340/*          } */
1341
1342/*        if(nxt_token_t == NEXUS_EQUAL)  */
1343/*          { */
1344/*            cur_token_t = NEXUS_VALUE; */
1345/*            nxt_token_t = NEXUS_PARM; */
1346/*            continue; */
1347/*          } */
1348
1349/*        if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE))  */
1350/*          { */
1351/*            Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list); */
1352/*            if(curr_com)  */
1353/*              { */
1354/*                nxt_token_t = curr_com->nxt_token_t; */
1355/*                cur_token_t = curr_com->cur_token_t; */
1356/*              } */
1357/*            if(cur_token_t != NEXUS_VALUE) continue; */
1358/*          } */
1359
1360/*        if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE))  */
1361/*          { */
1362/*            Find_Nexus_Parm(token,&curr_parm,curr_com); */
1363/*            if(curr_parm)  */
1364/*              { */
1365/*                nxt_token_t = curr_parm->nxt_token_t; */
1366/*                cur_token_t = curr_parm->cur_token_t; */
1367/*              } */
1368/*            if(cur_token_t != NEXUS_VALUE) continue; */
1369/*          } */
1370
1371/*        if(cur_token_t == NEXUS_VALUE) */
1372/*          { */
1373/*            if((curr_parm->fp)(token,curr_parm,io))  /\* Read in parameter value *\/ */
1374/*              { */
1375/*                nxt_token_t = NEXUS_PARM; */
1376/*                cur_token_t = -1; */
1377/*              } */
1378/*          } */
1379/*      } */
1380/*       while(strlen(token) > 0); */
1381/*     } */
1382
1383/*   Free(ori_s); */
1384/*   Free(token); */
1385
1386/*   return io->data; */
1387/* } */
1388
1389/* /\*********************************************************\/ */
1390
1391void Get_Nexus_Data(FILE *fp, option *io)
1392{
1393  char *token;
1394  nexcom *curr_com;
1395  nexparm *curr_parm;
1396  int nxt_token_t,cur_token_t;
1397
1398  token = (char *)mCalloc(T_MAX_TOKEN,sizeof(char));
1399         
1400  curr_com   = NULL;
1401  curr_parm  = NULL;
1402  nxt_token_t = NEXUS_COM; 
1403  cur_token_t = -1; 
1404
1405  do
1406    {
1407      if(!Get_Token(fp,token)) break;
1408
1409/*       PhyML_Printf("\n+ Token: '%s' next_token=%d cur_token=%d",token,nxt_token_t,cur_token_t); */
1410
1411      if(token[0] == ';') 
1412        {
1413          curr_com    = NULL;
1414          curr_parm   = NULL;
1415          nxt_token_t = NEXUS_COM;
1416          cur_token_t = -1;
1417        }
1418     
1419      if(nxt_token_t == NEXUS_EQUAL) 
1420        {
1421          cur_token_t = NEXUS_VALUE;
1422          nxt_token_t = NEXUS_PARM;
1423          continue;
1424        }
1425     
1426      if((nxt_token_t == NEXUS_COM) && (cur_token_t != NEXUS_VALUE)) 
1427        {
1428          Find_Nexus_Com(token,&curr_com,&curr_parm,io->nex_com_list);
1429          if(curr_com) 
1430            {
1431              nxt_token_t = curr_com->nxt_token_t;
1432              cur_token_t = curr_com->cur_token_t;
1433            }
1434          if(cur_token_t != NEXUS_VALUE) continue;
1435        }
1436     
1437      if((nxt_token_t == NEXUS_PARM) && (cur_token_t != NEXUS_VALUE)) 
1438        {
1439          Find_Nexus_Parm(token,&curr_parm,curr_com);
1440          if(curr_parm) 
1441            {
1442              nxt_token_t = curr_parm->nxt_token_t;
1443              cur_token_t = curr_parm->cur_token_t;
1444            }
1445          if(cur_token_t != NEXUS_VALUE) continue;
1446        }
1447     
1448      if(cur_token_t == NEXUS_VALUE)
1449        {
1450          if((curr_parm->fp)(token,curr_parm,io))  /* Read in parameter value */
1451            {
1452              nxt_token_t = NEXUS_PARM;
1453              cur_token_t = -1;
1454            }
1455        }
1456    }
1457  while(strlen(token) > 0);
1458 
1459  Free(token);
1460}
1461
1462//////////////////////////////////////////////////////////////
1463//////////////////////////////////////////////////////////////
1464
1465
1466int Get_Token(FILE *fp, char *token)
1467{
1468  char c;
1469 
1470  c = ' ';
1471  while(c == ' ' || c == '\t' || c == '\n') 
1472    {
1473      c = fgetc(fp);
1474      if(c == EOF) return 0;
1475    }
1476
1477  if(c == '"')
1478    {
1479      do
1480        {
1481          *token = c;
1482          token++;
1483          c = fgetc(fp);
1484          if(c == EOF) return 0;
1485        }
1486      while(c != '"');
1487      *token = c;
1488      /* c = fgetc(fp); */
1489      if(c == EOF) return 0;
1490      *(token+1) = '\0';
1491      return 1;
1492    }
1493
1494  if(c == '[')
1495    {
1496      Skip_Comment(fp);
1497      c = fgetc(fp);
1498      *token = c;
1499      token++;
1500      if(c == EOF) return 0;
1501      return 1;
1502    }
1503
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++; }
1515  else
1516    {
1517      while(isgraph(c) && c != ';' && c != '-' && c != ',' && c != '=')
1518        {
1519          *(token++) = c;
1520          c = fgetc(fp);
1521          if(c == EOF) return 0;
1522        }
1523
1524      fseek(fp,-1*sizeof(char),SEEK_CUR);
1525
1526    }
1527  *token = '\0';
1528  return 1;
1529}
1530
1531
1532/* void Get_Token(char *line, char *token) */
1533/* { */
1534/*   while(**line == ' ' || **line == '\t') (*line)++; */
1535
1536
1537/*   if(**line == '"')  */
1538/*     { */
1539/*       do { *token = **line; (*line)++; token++; } while(**line != '"'); */
1540/*       *token = **line; */
1541/*       (*line)++; */
1542/*       *(token+1) = '\0'; */
1543/*       return; */
1544/*     } */
1545
1546/*   if(**line == '[')  */
1547/*     { */
1548/*       int in_comment; */
1549
1550/*       in_comment = 1; */
1551/*       do  */
1552/*      {  */
1553/*        (*line)++;  */
1554/*        if(**line == '[')  */
1555/*          { */
1556/*            in_comment++; */
1557/*          } */
1558/*        else if(**line == ']') in_comment--;     */
1559/*      } */
1560/*       while(in_comment); */
1561/*       (*line)++; */
1562/*       return; */
1563/*     } */
1564
1565
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++; } */
1577/*   else */
1578/*     { */
1579/*       while(isgraph(**line) && **line != ';' && **line != '=' && **line != ',')  */
1580/*      { */
1581/*        *(token++) = **line; */
1582/*        (*line)++;  */
1583/*      } */
1584/*     } */
1585/*   *token = '\0'; */
1586/* } */
1587
1588//////////////////////////////////////////////////////////////
1589//////////////////////////////////////////////////////////////
1590
1591align **Get_Seq_Phylip(option *io)
1592{
1593  Read_Ntax_Len_Phylip(io->fp_in_align,&io->n_otu,&io->init_len);
1594
1595  if(io->n_otu > N_MAX_OTU)
1596    {
1597      PhyML_Printf("\n== The number of taxa should not exceed %d",N_MAX_OTU);
1598      Exit("\n");
1599    }
1600
1601  if(io->interleaved) io->data = Read_Seq_Interleaved(io);
1602  else                io->data = Read_Seq_Sequential(io);
1603
1604  return io->data;
1605}
1606
1607//////////////////////////////////////////////////////////////
1608//////////////////////////////////////////////////////////////
1609
1610void Read_Ntax_Len_Phylip(FILE *fp ,int *n_otu, int *n_tax)
1611{
1612  char *line;
1613  int readok;
1614 
1615  line = (char *)mCalloc(T_MAX_LINE,sizeof(char)); 
1616
1617  readok = 0;
1618  do
1619    {
1620      if(fscanf(fp,"%s",line) == EOF)
1621        {
1622          Free(line); 
1623          PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1624          Exit("\n");
1625        }
1626      else
1627        {
1628          if(strcmp(line,"\n") && strcmp(line,"\r") && strcmp(line,"\t"))
1629            {
1630              sscanf(line,"%d",n_otu);
1631              if(*n_otu <= 0) Warn_And_Exit("\n. The number of taxa cannot be negative.\n");
1632
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");
1636              else readok = 1;
1637            }
1638        }
1639    }while(!readok);
1640 
1641  Free(line);
1642}
1643
1644//////////////////////////////////////////////////////////////
1645//////////////////////////////////////////////////////////////
1646
1647align **Read_Seq_Sequential(option *io)
1648{
1649  int i;
1650  char *line;
1651  align **data;
1652/*   char c; */
1653  char *format;
1654
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 *));
1658
1659/*   while((c=fgetc(in))!='\n'); */
1660 /*  while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1661
1662  sprintf(format, "%%%ds", T_MAX_NAME);
1663
1664  For(i,io->n_otu)
1665    {
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));
1669
1670      data[i]->is_ambigu = NULL;
1671      data[i]->len = 0;
1672
1673      if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1674
1675      Check_Sequence_Name(data[i]->name);
1676
1677      while(data[i]->len < io->init_len * io->state_len) Read_One_Line_Seq(&data,i,io->fp_in_align);
1678
1679      if(data[i]->len != io->init_len * io->state_len)
1680        {
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);
1683          Warn_And_Exit("");
1684        }
1685    }
1686
1687  For(i,io->n_otu) data[i]->state[data[i]->len] = '\0';
1688 
1689  Restrict_To_Coding_Position(data,io);
1690
1691  Free(format);
1692  Free(line);
1693 
1694  return data;
1695}
1696
1697//////////////////////////////////////////////////////////////
1698//////////////////////////////////////////////////////////////
1699
1700align **Read_Seq_Interleaved(option *io)
1701{
1702  int i,end,num_block;
1703  char *line;
1704  align **data;
1705/*   char c; */
1706  char *format;
1707  fpos_t curr_pos;
1708
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 *));
1712
1713/*   while(((c=fgetc(io->fp_in_align))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); */
1714
1715  sprintf(format, "%%%ds", T_MAX_NAME);
1716
1717  end = 0;
1718  For(i,io->n_otu)
1719    {
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));
1723
1724      data[i]->len       = 0;
1725      data[i]->is_ambigu = NULL;
1726     
1727      if(!fscanf(io->fp_in_align,format,data[i]->name)) Exit("\n");
1728
1729      Check_Sequence_Name(data[i]->name);
1730
1731      if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1732        {
1733          end = 1;
1734          if((i != io->n_otu) && (i != io->n_otu-1))
1735            {
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);
1739              Exit("");
1740            }
1741          break;
1742        }
1743    }
1744 
1745  if(data[0]->len == io->init_len * io->state_len) end = 1;
1746
1747/*   if(end) printf("\n. finished yet '%c'\n",fgetc(io->fp_in_align)); */
1748  if(!end)
1749    {
1750
1751      end = 0;
1752
1753      num_block = 1;
1754      do
1755        {
1756          num_block++;
1757
1758          /* interblock */
1759          if(!fgets(line,T_MAX_LINE,io->fp_in_align)) break;
1760
1761          if(line[0] != 13 && line[0] != 10)
1762            {
1763              PhyML_Printf("\n== Err.: one or more missing sequences in block %d.\n",num_block-1);
1764              Exit("");
1765            }
1766
1767          For(i,io->n_otu) if(data[i]->len != io->init_len * io->state_len) break;
1768
1769          if(i == io->n_otu) break;
1770
1771          For(i,io->n_otu)
1772            {
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)) 
1776                {
1777                  PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1778                  Warn_And_Exit("");
1779                }
1780              if(line && strcmp(line,data[i]->name)) fsetpos(io->fp_in_align,&curr_pos);                 
1781
1782              if(data[i]->len > io->init_len * io->state_len)
1783                {
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);
1786                  Exit("");
1787                }
1788              else if(!Read_One_Line_Seq(&data,i,io->fp_in_align))
1789                {
1790                  end = 1;
1791                  if((i != io->n_otu) && (i != io->n_otu-1))
1792                    {
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);
1795                      Exit("");
1796                    }
1797                  break;
1798                }
1799            }
1800        }while(!end);
1801    }
1802
1803  For(i,io->n_otu) data[i]->state[data[i]->len] = '\0';
1804
1805  For(i,io->n_otu)
1806    {
1807      if(data[i]->len != io->init_len * io->state_len)
1808        {
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);
1810          Exit("");
1811        }
1812    }
1813
1814  Restrict_To_Coding_Position(data,io);
1815
1816  Free(format);
1817  Free(line);
1818
1819  return data;
1820}
1821
1822//////////////////////////////////////////////////////////////
1823//////////////////////////////////////////////////////////////
1824
1825
1826int Read_One_Line_Seq(align ***data, int num_otu, FILE *in)
1827{
1828  char c = ' ';
1829  int nchar = 0;
1830 
1831  while(1)
1832    {
1833/*       if((c == EOF) || (c == '\n') || (c == '\r')) break; */
1834     
1835      if((c == 13) || (c == 10)) 
1836        {
1837          /* PhyML_Printf("[%d %d]\n",c,nchar); fflush(NULL); */
1838          if(!nchar)
1839            {
1840              c=(char)fgetc(in);
1841              continue;
1842            }
1843          else 
1844            { 
1845              /* PhyML_Printf("break\n"); */
1846              break; 
1847            }
1848        }
1849      else if(c == EOF)
1850        {
1851/*        PhyML_Printf("EOL\n"); */
1852          break;
1853        }
1854      else if((c == ' ') || (c == '\t') || (c == 32)) 
1855        {
1856/*        PhyML_Printf("[%d]",c); */
1857          c=(char)fgetc(in); 
1858          continue;
1859        }
1860
1861      nchar++;
1862      Uppercase(&c);
1863     
1864      if(c == '.')
1865        {
1866          c = (*data)[0]->state[(*data)[num_otu]->len];
1867          if(!num_otu)
1868            Warn_And_Exit("\n== Err: Symbol \".\" should not appear in the first sequence\n");
1869        }
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); */
1874      if(c == ';') break;
1875    }
1876
1877  /* printf("\n. Exit nchar: %d [%d]\n",nchar,c==EOF); */
1878  if(c == EOF) return 0;
1879  else return 1; 
1880}
1881
1882//////////////////////////////////////////////////////////////
1883//////////////////////////////////////////////////////////////
1884
1885t_tree *Read_Tree_File(option *io)
1886{
1887  t_tree *tree;
1888
1889  if(!io->fp_in_tree)
1890    {
1891      PhyML_Printf("\n== Filehandle to '%s' seems to be closed.",io->in_tree_file);
1892      Exit("\n");
1893    }
1894
1895
1896  Detect_Tree_File_Format(io);
1897
1898  io->treelist->list_size = 0;
1899
1900  switch(io->tree_file_format)
1901    {
1902    case PHYLIP: 
1903      {
1904        do
1905          {
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++;
1912          }while(io->tree);
1913        break;
1914      }
1915    case NEXUS:
1916      {
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);
1920        Free_Nexus(io);
1921        break;
1922      }
1923    default:
1924      {
1925        PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1926        Warn_And_Exit("");
1927        break;
1928      }
1929    }
1930 
1931  if(!io->long_tax_names)
1932    {
1933      int i;
1934
1935      tree = io->treelist->tree[0];
1936
1937      io->long_tax_names  = (char **)mCalloc(tree->n_otu,sizeof(char *));
1938      io->short_tax_names = (char **)mCalloc(tree->n_otu,sizeof(char *));
1939
1940      For(i,tree->n_otu)
1941        {
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);
1946        }
1947    }
1948  return NULL;
1949}
1950
1951//////////////////////////////////////////////////////////////
1952//////////////////////////////////////////////////////////////
1953
1954char *Return_Tree_String_Phylip(FILE *fp_input_tree)
1955{
1956  char *line;
1957  int i;
1958  char c;
1959  int open,maxopen;
1960
1961  if(fp_input_tree == NULL)
1962    {
1963      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1964      Warn_And_Exit("");     
1965    }
1966 
1967  do
1968    {
1969      c=fgetc(fp_input_tree);
1970    }
1971  while((c != '(') && (c != EOF));
1972 
1973 
1974  if(c==EOF) return NULL;
1975
1976  line = (char *)mCalloc(1,sizeof(char));
1977  open = 1;
1978  maxopen = open;
1979
1980  i=0;
1981  for(;;)
1982    {
1983      if((c == ' ') || (c == '\n'))
1984        {
1985          c=fgetc(fp_input_tree);
1986          if(c == EOF || c == ';') break;
1987          else continue;
1988        }
1989     
1990      if(c == '[')
1991        {
1992          Skip_Comment(fp_input_tree);
1993          c = fgetc(fp_input_tree);
1994          if(c == EOF || c == ';') break;
1995        }
1996
1997      line = (char *)mRealloc(line,i+2,sizeof(char));
1998
1999      line[i]=c;
2000      i++;
2001      c=fgetc(fp_input_tree);
2002      if(c==EOF || c==';') break;
2003      if(c=='(') open++;
2004      if(c==')') open--;
2005      if(open>maxopen) maxopen = open;
2006    }
2007  line[i] = '\0';
2008 
2009
2010  /* if(maxopen == 1) return NULL; */
2011  return line;
2012}
2013
2014//////////////////////////////////////////////////////////////
2015//////////////////////////////////////////////////////////////
2016
2017
2018t_tree *Read_Tree_File_Phylip(FILE *fp_input_tree)
2019{
2020  char *line;
2021  t_tree *tree;
2022
2023  line = Return_Tree_String_Phylip(fp_input_tree);
2024  tree = Read_Tree(&line);
2025  Free(line);
2026 
2027  return tree;
2028}
2029
2030//////////////////////////////////////////////////////////////
2031//////////////////////////////////////////////////////////////
2032
2033void Print_Site(calign *cdata, int num, int n_otu, char *sep, int stepsize, FILE *fp)
2034{
2035  int i,j;
2036  PhyML_Fprintf(fp,"\n");
2037  For(i,n_otu)
2038    {
2039      PhyML_Fprintf(fp,"%20s ",cdata->c_seq[i]->name);
2040      For(j,stepsize)
2041        PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[num+j]);
2042      PhyML_Fprintf(fp,"%s",sep);
2043    }
2044  PhyML_Fprintf(fp,"%s",sep);
2045}
2046
2047//////////////////////////////////////////////////////////////
2048//////////////////////////////////////////////////////////////
2049
2050void Print_Site_Lk(t_tree *tree, FILE *fp)
2051{
2052  int site;
2053  int catg;
2054  char *s;
2055  phydbl postmean;
2056
2057  if(!tree->io->print_site_lnl)
2058    {
2059      PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
2060      Exit("");
2061    }
2062
2063  if(!tree->io->print_trace)
2064    {
2065      s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2066     
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");
2071     
2072      sprintf(s,"Site");
2073      PhyML_Fprintf(fp, "%-7s",s);
2074
2075      sprintf(s,"Pattern");
2076      PhyML_Fprintf(fp, "%-9s",s);
2077     
2078      sprintf(s,"P(D|M)");
2079      PhyML_Fprintf(fp,"%-16s",s);
2080     
2081      if(tree->mod->ras->n_catg > 1)
2082        {
2083          For(catg,tree->mod->ras->n_catg)
2084            {
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);
2087            }
2088         
2089          sprintf(s,"Posterior mean");
2090          PhyML_Fprintf(fp,"%-22s",s);
2091        }
2092     
2093     
2094      if(tree->mod->ras->invar)
2095        {
2096          sprintf(s,"P(D|M,rr[0]=0)");
2097          PhyML_Fprintf(fp,"%-16s",s);
2098        }
2099      PhyML_Fprintf(fp,"\n");
2100     
2101      For(site,tree->data->init_len)
2102        {
2103         
2104
2105          PhyML_Fprintf(fp,"%-7d",site+1);
2106          PhyML_Fprintf(fp,"%-9d",tree->data->sitepatt[site]);
2107
2108          PhyML_Fprintf(fp,"%-16g",tree->cur_site_lk[tree->data->sitepatt[site]]);     
2109          if(tree->mod->ras->n_catg > 1)
2110            {
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]]));
2113
2114              postmean = .0;
2115              For(catg,tree->mod->ras->n_catg) 
2116                postmean += 
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]];
2121
2122              PhyML_Fprintf(fp,"%-22g",postmean);
2123            }
2124          if(tree->mod->ras->invar)
2125            {
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]]]);
2128              else
2129                PhyML_Fprintf(fp,"%-16g",0.0);
2130            }
2131
2132         
2133
2134          PhyML_Fprintf(fp,"\n");
2135        }
2136      Free(s);
2137    }
2138  else
2139    {
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");
2143    }
2144}
2145
2146//////////////////////////////////////////////////////////////
2147//////////////////////////////////////////////////////////////
2148
2149void Print_Seq(FILE *fp, align **data, int n_otu)
2150{
2151  int i,j;
2152
2153  PhyML_Fprintf(fp,"%d\t%d\n",n_otu,data[0]->len);
2154  For(i,n_otu)
2155    {
2156      For(j,20)
2157        {
2158          if(j<(int)strlen(data[i]->name))
2159            fputc(data[i]->name[j],fp);
2160          else fputc(' ',fp);
2161        }
2162/*       PhyML_Printf("%10d  ",i); */
2163      For(j,data[i]->len)
2164        {
2165          PhyML_Fprintf(fp,"%c",data[i]->state[j]);
2166        }
2167      PhyML_Fprintf(fp,"\n");
2168    }
2169}
2170
2171//////////////////////////////////////////////////////////////
2172//////////////////////////////////////////////////////////////
2173
2174
2175void Print_CSeq(FILE *fp, int compressed, calign *cdata)
2176{
2177  int i,j;
2178  int n_otu;
2179 
2180  n_otu = cdata->n_otu;
2181  if(cdata->format == 0)
2182    {
2183      PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
2184    }
2185  else
2186    {
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");
2192    }
2193 
2194  For(i,n_otu)
2195    {
2196      For(j,50)
2197        {
2198          if(j<(int)strlen(cdata->c_seq[i]->name))
2199            fputc(cdata->c_seq[i]->name[j],fp);
2200          else fputc(' ',fp);
2201        }
2202     
2203      if(compressed == YES) /* Print out compressed sequences */
2204        PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
2205      else /* Print out uncompressed sequences */
2206        {
2207          For(j,cdata->init_len)
2208            {
2209              PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
2210            }
2211        }
2212      PhyML_Fprintf(fp,"\n");
2213    }
2214  PhyML_Fprintf(fp,"\n");
2215
2216  if(cdata->format == 1)
2217    {
2218      PhyML_Fprintf(fp,";\n");
2219      PhyML_Fprintf(fp,"END;\n");
2220    }
2221
2222
2223/*   PhyML_Printf("\t"); */
2224/*   For(j,cdata->crunch_len) */
2225/*     PhyML_Printf("%.0f ",cdata->wght[j]); */
2226/*   PhyML_Printf("\n"); */
2227}
2228
2229//////////////////////////////////////////////////////////////
2230//////////////////////////////////////////////////////////////
2231
2232
2233void Print_CSeq_Select(FILE *fp, int compressed, calign *cdata, t_tree *tree)
2234{
2235  int i,j;
2236  int n_otu;
2237  phydbl eps;
2238
2239  int slice = 14;
2240
2241  eps = 1.E-6;
2242  n_otu = 0;
2243  For(i,cdata->n_otu)
2244    if(tree->rates->nd_t[i] < tree->rates->time_slice_lims[slice] + eps)
2245      n_otu++;
2246 
2247  PhyML_Fprintf(fp,"%d\t%d\n",n_otu,cdata->init_len);
2248
2249  n_otu = cdata->n_otu;
2250
2251  For(i,n_otu)
2252    {
2253      if(tree->rates->nd_t[i] < tree->rates->time_slice_lims[slice] + eps)
2254        {
2255          For(j,50)
2256            {
2257              if(j<(int)strlen(cdata->c_seq[i]->name))
2258                fputc(cdata->c_seq[i]->name[j],fp);
2259              else fputc(' ',fp);
2260            }
2261         
2262          if(compressed == YES) /* Print out compressed sequences */
2263            PhyML_Fprintf(fp,"%s",cdata->c_seq[i]->state);
2264          else /* Print out uncompressed sequences */
2265            {
2266              For(j,cdata->init_len)
2267                {
2268                  PhyML_Fprintf(fp,"%c",cdata->c_seq[i]->state[cdata->sitepatt[j]]);
2269                }
2270            }
2271          PhyML_Fprintf(fp,"\n");
2272        }
2273    }
2274
2275  if(cdata->format == 1)
2276    {
2277      PhyML_Fprintf(fp,";\n");
2278      PhyML_Fprintf(fp,"END;\n");
2279    }
2280
2281
2282/*   PhyML_Printf("\t"); */
2283/*   For(j,cdata->crunch_len) */
2284/*     PhyML_Printf("%.0f ",cdata->wght[j]); */
2285/*   PhyML_Printf("\n"); */
2286}
2287
2288//////////////////////////////////////////////////////////////
2289//////////////////////////////////////////////////////////////
2290
2291void Print_Dist(matrix *mat)
2292{
2293  int i,j;
2294
2295  For(i,mat->n_otu)
2296    {
2297      PhyML_Printf("%s ",mat->name[i]);
2298
2299      For(j,mat->n_otu)
2300        PhyML_Printf("%9.6f ",mat->dist[i][j]);
2301      PhyML_Printf("\n");
2302    }
2303}
2304
2305//////////////////////////////////////////////////////////////
2306//////////////////////////////////////////////////////////////
2307
2308
2309void Print_Node(t_node *a, t_node *d, t_tree *tree)
2310{
2311  int i;
2312  int dir;
2313  dir = -1;
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)
2321    {
2322      if(a->b[i])
2323        {
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!");
2328          break;
2329        }
2330    }
2331  PhyML_Printf("\n");
2332
2333  if(d->tax) return;
2334  else
2335    For(i,3)
2336      if(d->v[i] != a && d->b[i] != tree->e_root) Print_Node(d,d->v[i],tree);
2337}
2338
2339//////////////////////////////////////////////////////////////
2340//////////////////////////////////////////////////////////////
2341
2342void Print_Model(t_mod *mod)
2343{
2344  int i,j,k;
2345
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);
2369 
2370  PhyML_Printf("\n. Pi\n");
2371  For(i,mod->ns) PhyML_Printf(" %f ",mod->e_frq->pi->v[i]);
2372  PhyML_Printf("\n");
2373  For(i,mod->ns) PhyML_Printf(" %f ",mod->e_frq->pi_unscaled->v[i]);
2374 
2375  PhyML_Printf("\n. Rates\n");
2376  For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_r_proba->v[i]);
2377  PhyML_Printf("\n");
2378  For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_r_proba_unscaled->v[i]);
2379  PhyML_Printf("\n");
2380  For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_rr->v[i]);
2381  PhyML_Printf("\n");
2382  For(i,mod->ras->n_catg) PhyML_Printf(" %f ",mod->ras->gamma_rr_unscaled->v[i]);
2383 
2384 
2385 
2386  PhyML_Printf("\n. Qmat \n");
2387  if(mod->whichmodel == CUSTOM)
2388    {
2389      fflush(NULL);
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);}
2394    }
2395  For(i,mod->ns)
2396    {
2397      PhyML_Printf("  ");
2398      For(j,4)
2399        PhyML_Printf("%8.5f  ",mod->r_mat->qmat->v[i*4+j]);
2400      PhyML_Printf("\n");
2401    }
2402
2403  PhyML_Printf("\n. Freqs");
2404  PhyML_Printf("\n");
2405  For(i,mod->ns) PhyML_Printf(" %12f ",mod->user_b_freq->v[i]);
2406  PhyML_Printf("\n");
2407  For(i,mod->ns) PhyML_Printf(" %12f ",mod->e_frq->pi->v[i]);
2408  PhyML_Printf("\n");
2409  For(i,mod->ns) PhyML_Printf(" %12f ",mod->e_frq->pi_unscaled->v[i]);
2410
2411  PhyML_Printf("\n. Eigen\n");
2412  For(i,2*mod->ns)       PhyML_Printf(" %f ",mod->eigen->space[i]);     
2413  PhyML_Printf("\n");
2414  For(i,2*mod->ns)       PhyML_Printf(" %f ",mod->eigen->space_int[i]); 
2415  PhyML_Printf("\n");
2416  For(i,mod->ns)         PhyML_Printf(" %f ",mod->eigen->e_val[i]);     
2417  PhyML_Printf("\n");
2418  For(i,mod->ns)         PhyML_Printf(" %f ",mod->eigen->e_val_im[i]);   
2419  PhyML_Printf("\n");
2420  For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect[i]);   
2421  PhyML_Printf("\n");
2422  For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect[i]);   
2423  PhyML_Printf("\n");
2424  For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->r_e_vect_im[i]);
2425  PhyML_Printf("\n");
2426  For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->l_e_vect[i]);   
2427  PhyML_Printf("\n");
2428  For(i,mod->ns*mod->ns) PhyML_Printf(" %f ",mod->eigen->q[i]);         
2429  PhyML_Printf("\n");
2430 
2431  PhyML_Printf("\n. Pij");
2432  For(k,mod->ras->n_catg)
2433    {
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]);
2436      For(i,mod->ns)
2437        {
2438          PhyML_Printf("  ");
2439          For(j,mod->ns)
2440            PhyML_Printf("%8.5f  ",mod->Pij_rr->v[k*mod->ns*mod->ns+i*mod->ns+j]);
2441          PhyML_Printf("\n");
2442        }
2443    }
2444 
2445
2446  PhyML_Printf("\n");
2447
2448  fflush(NULL);
2449}
2450
2451//////////////////////////////////////////////////////////////
2452//////////////////////////////////////////////////////////////
2453
2454void Print_Mat(matrix *mat)
2455{
2456  int i,j;
2457
2458  PhyML_Printf("%d",mat->n_otu);
2459  PhyML_Printf("\n");
2460
2461  For(i,mat->n_otu)
2462    {
2463      For(j,13)
2464        {
2465          if(j>=(int)strlen(mat->name[i])) putchar(' ');
2466          else putchar(mat->name[i][j]);
2467        }
2468
2469      For(j,mat->n_otu)
2470        {
2471          char s[2]="-";
2472          if(mat->dist[i][j] < .0)
2473            PhyML_Printf("%12s",s);
2474          else
2475            PhyML_Printf("%12f",mat->dist[i][j]);
2476        }
2477      PhyML_Printf("\n");
2478    }
2479}
2480
2481//////////////////////////////////////////////////////////////
2482//////////////////////////////////////////////////////////////
2483
2484FILE *Openfile(char *filename, int mode)
2485{
2486  /* mode = 0 -> read */
2487  /* mode = 1 -> write */
2488  /* mode = 2 -> append */
2489
2490  FILE *fp;
2491  char *s;
2492  int open_test=0;
2493
2494/*   s = (char *)mCalloc(T_MAX_FILE,sizeof(char)); */
2495
2496/*   strcpy(s,filename); */
2497
2498  s = filename;
2499
2500  fp = NULL;
2501
2502  switch(mode)
2503    {
2504    case 0 :
2505      {
2506        while(!(fp = (FILE *)fopen(s,"r")) && ++open_test<10)
2507          {
2508            PhyML_Printf("\n. Can't open file '%s', enter a new name : ",s);
2509            Getstring_Stdin(s);
2510          }
2511        break;
2512      }
2513    case 1 :
2514      {
2515        fp = (FILE *)fopen(s,"w");
2516        break;
2517      }
2518    case 2 :
2519      {
2520        fp = (FILE *)fopen(s,"a");
2521        break;
2522      }
2523
2524    default : break;
2525
2526    }
2527
2528/*   Free(s); */
2529
2530  return fp;
2531}
2532
2533//////////////////////////////////////////////////////////////
2534//////////////////////////////////////////////////////////////
2535
2536void 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)
2537{
2538  char *s;
2539  div_t hour,min;
2540  int i;
2541
2542/*   For(i,2*tree->n_otu-3) fprintf(fp_out,"\n. * Edge %3d: %f",i,tree->a_edges[i]->l); */
2543 
2544  if((!n_data_set) || (!num_tree))
2545    {
2546      rewind(fp_out);
2547      Print_Banner_Small(fp_out);
2548    }
2549
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);
2552
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);
2557 
2558  if(io->mod->s_opt->opt_topo)
2559    {
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");
2563    }
2564  else
2565    {
2566      PhyML_Fprintf(fp_out,"\n. Tree topology: \t\t\tfixed");
2567    }
2568
2569
2570  /* was after Sequence file ; moved here FLT */
2571  s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2572  if(io->in_tree == 2)
2573    {
2574      strcat(strcat(strcat(s,"user tree ("),io->in_tree_file),")");
2575    }
2576  else
2577    {
2578      if(!io->mod->s_opt->random_input_tree)
2579        {
2580          if(io->in_tree == 0)
2581            strcat(s,"BioNJ");
2582          if(io->in_tree == 1)
2583            strcat(s,"parsimony");
2584        }
2585      else
2586        {
2587          strcat(s,"random tree");
2588        }
2589    }
2590
2591  PhyML_Fprintf(fp_out,"\n. Initial tree: \t\t\t%s",s);
2592  Free(s);
2593
2594  if(tree->io->datatype == NT)
2595    {
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);
2599    }
2600  else if(tree->io->datatype == AA)
2601    {
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);
2604    }
2605  else
2606    {
2607      fprintf(fp_out,"\n. Substitution model: \t\t\t%s",tree->mod->modelname->s);
2608    }
2609
2610
2611  PhyML_Fprintf(fp_out,"\n. Number of taxa: \t\t\t%d",tree->n_otu);/*added FLT*/
2612
2613  PhyML_Fprintf(fp_out,"\n. Log-likelihood: \t\t\t%.5f",tree->c_lnL);/*was last ; moved here FLT*/
2614
2615  Unconstraint_Lk(tree);
2616  PhyML_Fprintf(fp_out,"\n. Unconstrained likelihood: \t\t%.5f",tree->unconstraint_lk);
2617
2618  PhyML_Fprintf(fp_out,"\n. Parsimony: \t\t\t\t%d",tree->c_pars);
2619
2620  PhyML_Fprintf(fp_out,"\n. Tree size: \t\t\t\t%.5f",Get_Tree_Size(tree));
2621
2622  /* if(tree->mod->ras->n_catg > 1 && tree->mod->ras->free_mixt_rates == NO) */
2623  if(tree->mod->ras->free_mixt_rates == NO)
2624    {
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)
2629        {
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]);
2631        }
2632    }
2633  else if(tree->mod->ras->free_mixt_rates == YES)
2634    {
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)
2638        {
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]);
2640        }
2641    }
2642
2643  if(tree->mod->ras->invar) PhyML_Fprintf(fp_out,"\n. Proportion of invariant: \t\t%.3f",tree->mod->ras->pinvar->v);
2644
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);
2646
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)
2653    {
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));
2658    }
2659
2660  if(tree->io->datatype == NT)
2661    {
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]);
2667    }
2668
2669  /*****************************************/
2670  if((tree->mod->whichmodel == GTR) ||
2671     (tree->mod->whichmodel == CUSTOM))
2672    {
2673      int i,j;
2674
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);
2680
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]);
2689
2690
2691      PhyML_Fprintf(fp_out,"\n. Instantaneous rate matrix : ");
2692      PhyML_Fprintf(fp_out,"\n  [A---------C---------G---------T------]\n");
2693      For(i,4)
2694        {
2695          PhyML_Fprintf(fp_out,"  ");
2696          For(j,4)
2697            PhyML_Fprintf(fp_out,"%8.5f  ",tree->mod->r_mat->qmat->v[i*4+j]);
2698          PhyML_Fprintf(fp_out,"\n");
2699        }
2700      PhyML_Fprintf(fp_out,"\n");
2701    }
2702  /*****************************************/
2703
2704
2705  if(io->ratio_test == 1)
2706    {
2707      PhyML_Fprintf(fp_out,". aLRT statistics to test branches");
2708    }
2709  else if(io->ratio_test == 2)
2710    {
2711      PhyML_Fprintf(fp_out,". aLRT branch supports (cubic approximation, mixture of Chi2s distribution)");
2712    }
2713
2714
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);
2720
2721  hour = div(t_end-t_beg,3600);
2722  min  = div(t_end-t_beg,60  );
2723
2724  min.quot -= hour.quot*60;
2725
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));
2727
2728 
2729  if(add_citation == YES)
2730    { 
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");
2742    }
2743  else
2744    {
2745      PhyML_Fprintf(fp_out,"\n\n");
2746      PhyML_Fprintf(fp_out," oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
2747      PhyML_Fprintf(fp_out,"\n");
2748
2749    }
2750}
2751
2752//////////////////////////////////////////////////////////////
2753//////////////////////////////////////////////////////////////
2754
2755/*FLT wrote this function*/
2756void Print_Fp_Out_Lines(FILE *fp_out, time_t t_beg, time_t t_end, t_tree *tree, option *io, int n_data_set)
2757{
2758  char *s;
2759  /*div_t hour,min;*/
2760
2761  if (n_data_set==1)
2762      {
2763
2764        PhyML_Fprintf(fp_out,". Sequence file : [%s]\n\n", Basename(io->in_align_file));
2765
2766        if((tree->io->datatype == NT) || (tree->io->datatype == AA))
2767          {
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));
2771          }
2772
2773        s = (char *)mCalloc(100,sizeof(char));
2774
2775        switch(io->in_tree)
2776          {
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; }
2782          }
2783
2784        PhyML_Fprintf(fp_out,". Initial tree : [%s]\n\n",s);
2785
2786        Free(s);
2787
2788        PhyML_Fprintf(fp_out,"\n");
2789
2790        /*headline 1*/
2791        PhyML_Fprintf(fp_out, ". Data\t");
2792
2793        PhyML_Fprintf(fp_out,"Nb of \t");
2794
2795        PhyML_Fprintf(fp_out,"Likelihood\t");
2796
2797        PhyML_Fprintf(fp_out, "Discrete   \t");
2798
2799        if(tree->mod->ras->n_catg > 1)
2800          PhyML_Fprintf(fp_out, "Number of \tGamma shape\t");
2801
2802        PhyML_Fprintf(fp_out,"Proportion of\t");
2803
2804        if(tree->mod->whichmodel <= 6)
2805          PhyML_Fprintf(fp_out,"Transition/ \t");
2806
2807        PhyML_Fprintf(fp_out,"Nucleotides frequencies               \t");
2808
2809        if((tree->mod->whichmodel == GTR) ||
2810           (tree->mod->whichmodel == CUSTOM))
2811          PhyML_Fprintf(fp_out,"Instantaneous rate matrix              \t");
2812
2813        /*    PhyML_Fprintf(fp_out,"Time\t");*/
2814
2815        PhyML_Fprintf(fp_out, "\n");
2816
2817
2818        /*headline 2*/
2819        PhyML_Fprintf(fp_out, "  set\t");
2820
2821        PhyML_Fprintf(fp_out,"taxa\t");
2822
2823        PhyML_Fprintf(fp_out,"loglk     \t");
2824
2825        PhyML_Fprintf(fp_out, "gamma model\t");
2826
2827        if(tree->mod->ras->n_catg > 1)
2828          PhyML_Fprintf(fp_out, "categories\tparameter  \t");
2829
2830        PhyML_Fprintf(fp_out,"invariant    \t");
2831
2832        if(tree->mod->whichmodel <= 6)
2833          PhyML_Fprintf(fp_out,"transversion\t");
2834
2835        PhyML_Fprintf(fp_out,"f(A)      f(C)      f(G)      f(T)    \t");
2836
2837        if((tree->mod->whichmodel == GTR) ||
2838           (tree->mod->whichmodel == CUSTOM))
2839          PhyML_Fprintf(fp_out,"[A---------C---------G---------T------]\t");
2840
2841        /*    PhyML_PhyML_Fprintf(fp_out,"used\t");*/
2842
2843        PhyML_Fprintf(fp_out, "\n");
2844
2845
2846        /*headline 3*/
2847        if(tree->mod->whichmodel == TN93)
2848          {
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");
2854          }
2855
2856          PhyML_Fprintf(fp_out, "\n");
2857      }
2858
2859
2860  /*line items*/
2861
2862  PhyML_Fprintf(fp_out,"  #%d\t",n_data_set);
2863
2864  PhyML_Fprintf(fp_out,"%d   \t",tree->n_otu);
2865
2866  PhyML_Fprintf(fp_out,"%.5f\t",tree->c_lnL);
2867
2868  PhyML_Fprintf(fp_out,"%s        \t",
2869          (tree->mod->ras->n_catg>1)?("Yes"):("No "));
2870  if(tree->mod->ras->n_catg > 1)
2871    {
2872      PhyML_Fprintf(fp_out,"%d        \t",tree->mod->ras->n_catg);
2873      PhyML_Fprintf(fp_out,"%.3f    \t",tree->mod->ras->alpha->v);
2874    }
2875
2876  /*if(tree->mod->ras->invar)*/
2877    PhyML_Fprintf(fp_out,"%.3f    \t",tree->mod->ras->pinvar->v);
2878
2879  if(tree->mod->whichmodel <= 5)
2880    {
2881      PhyML_Fprintf(fp_out,"%.3f     \t",tree->mod->kappa->v);
2882    }
2883  else if(tree->mod->whichmodel == TN93)
2884    {
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));
2889    }
2890
2891
2892  if(tree->io->datatype == NT)
2893    {
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]);
2898    }
2899  /*
2900  hour = div(t_end-t_beg,3600);
2901  min  = div(t_end-t_beg,60  );
2902
2903  min.quot -= hour.quot*60;
2904
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));
2908  */
2909
2910  /*****************************************/
2911  if((tree->mod->whichmodel == GTR) || (tree->mod->whichmodel == CUSTOM))
2912    {
2913      int i,j;
2914
2915      For(i,4)
2916        {
2917          if (i!=0) {
2918            /*format*/
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");
2922          }
2923          For(j,4)
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");
2926        }
2927    }
2928  /*****************************************/
2929
2930  PhyML_Fprintf(fp_out, "\n\n");
2931}
2932
2933//////////////////////////////////////////////////////////////
2934//////////////////////////////////////////////////////////////
2935
2936
2937void Print_Freq(t_tree *tree)
2938{
2939
2940  switch(tree->io->datatype)
2941    {
2942    case NT:
2943      {
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]);
2948
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]); */
2961        break;
2962      }
2963    case AA:
2964      {
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]);
2985
2986        PhyML_Printf("N : %f\n",tree->mod->e_frq->pi->v[20]);
2987        break;
2988      }
2989    default : {break;}
2990    }
2991}
2992
2993//////////////////////////////////////////////////////////////
2994//////////////////////////////////////////////////////////////
2995
2996
2997void Print_Settings(option *io)
2998{
2999  int answer;
3000  char *s;
3001
3002  s = (char *)mCalloc(100,sizeof(char));
3003 
3004  PhyML_Printf("\n\n\n");
3005  PhyML_Printf("\n\n");
3006
3007  PhyML_Printf("                                 ..........................                                      \n");
3008  PhyML_Printf(" ooooooooooooooooooooooooooooo        CURRENT SETTINGS        ooooooooooooooooooooooooooooooooooo\n");
3009  PhyML_Printf("                                 ..........................                                      \n");
3010
3011  PhyML_Printf("\n                . Sequence filename:\t\t\t\t %s", Basename(io->in_align_file));
3012
3013  if(io->datatype == NT) strcpy(s,"dna");
3014  else if(io->datatype == AA) strcpy(s,"aa");
3015  else strcpy(s,"generic");
3016
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);
3019
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);
3022
3023  PhyML_Printf("\n                . Nb of bootstrapped data sets:\t\t\t %d", io->mod->bootstrap);
3024
3025  if (io->mod->bootstrap > 0)
3026    PhyML_Printf("\n                . Compute approximate likelihood ratio test:\t no");
3027  else
3028    {
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)");
3039    }
3040
3041  PhyML_Printf("\n                . Model name:\t\t\t\t\t %s", io->mod->modelname->s);
3042
3043  if(io->datatype == AA && io->mod->whichmodel == CUSTOMAA) PhyML_Printf(" (%s)",io->mod->aa_rate_mat_file->s);
3044
3045  if (io->datatype == NT)
3046    {
3047      if ((io->mod->whichmodel == K80)  ||
3048          (io->mod->whichmodel == HKY85)||
3049          (io->mod->whichmodel == F84)  ||
3050          (io->mod->whichmodel == TN93))
3051        {
3052          if (io->mod->s_opt->opt_kappa)
3053            PhyML_Printf("\n                . Ts/tv ratio:\t\t\t\t\t estimated");
3054          else
3055            PhyML_Printf("\n                . Ts/tv ratio:\t\t\t\t\t %f", io->mod->kappa->v);
3056        }
3057    }
3058
3059  if (io->mod->s_opt->opt_pinvar)
3060    PhyML_Printf("\n                . Proportion of invariable sites:\t\t estimated");
3061  else
3062    PhyML_Printf("\n                . Proportion of invariable sites:\t\t %f", io->mod->ras->pinvar->v);
3063
3064
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)
3067    {
3068      if(io->mod->ras->free_mixt_rates == NO)
3069        {
3070          if(io->mod->s_opt->opt_alpha)
3071            PhyML_Printf("\n                . Gamma distribution parameter:\t\t\t estimated");
3072          else
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"));
3075        }
3076    }
3077   
3078 
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)
3082    {
3083      if((io->mod->whichmodel != JC69) &&
3084         (io->mod->whichmodel != K80)  &&
3085         (io->mod->whichmodel != F81))
3086        {
3087          if(!io->mod->s_opt->user_state_freq)
3088            {
3089              PhyML_Printf("\n                . Nucleotide equilibrium frequencies:\t\t %s", (io->mod->s_opt->opt_state_freq) ? ("ML"):("empirical"));
3090            }
3091          else
3092            {
3093              PhyML_Printf("\n                . Nucleotide equilibrium frequencies:\t\t %s","user-defined");
3094            }
3095        }
3096    }
3097
3098  PhyML_Printf("\n                . Optimise tree topology:\t\t\t %s", (io->mod->s_opt->opt_topo) ? "yes": "no");
3099
3100  switch(io->in_tree)
3101    {
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; }
3107    }
3108
3109  if(io->mod->s_opt->opt_topo)
3110    {
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");
3114
3115
3116
3117      PhyML_Printf("\n                . Starting tree:\t\t\t\t %s",s);
3118
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);     
3122    }
3123  else
3124    if(!io->mod->s_opt->random_input_tree)
3125      PhyML_Printf("\n                . Evaluated tree:\t\t\t\t \"%s\"",s);
3126
3127  PhyML_Printf("\n                . Optimise branch lengths:\t\t\t %s", (io->mod->s_opt->opt_bl) ? "yes": "no");
3128
3129  answer = 0;
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;
3135 
3136  PhyML_Printf("\n                . Optimise substitution model parameters:\t %s", (answer) ? "yes": "no");
3137
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);
3142
3143
3144  PhyML_Printf("\n\n oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo\n");
3145
3146  PhyML_Printf("\n\n");
3147  fflush(NULL);
3148 
3149  Free(s);
3150}
3151
3152void Print_Banner(FILE *fp)
3153{
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");
3167}
3168
3169//////////////////////////////////////////////////////////////
3170//////////////////////////////////////////////////////////////
3171
3172
3173void Print_Banner_Small(FILE *fp)
3174{
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");
3181}
3182
3183//////////////////////////////////////////////////////////////
3184//////////////////////////////////////////////////////////////
3185
3186
3187
3188void Print_Data_Set_Number(option *io, FILE *fp)
3189{
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");
3194}
3195//////////////////////////////////////////////////////////////
3196//////////////////////////////////////////////////////////////
3197
3198void Print_Lk(t_tree *tree, char *string)
3199{
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,
3203               string);
3204#ifndef QUIET
3205  fflush(NULL);
3206#endif
3207}
3208
3209//////////////////////////////////////////////////////////////
3210//////////////////////////////////////////////////////////////
3211
3212
3213void Print_Pars(t_tree *tree)
3214{
3215  time(&(tree->t_current));
3216  PhyML_Printf("\n. (%5d sec) [%5d]",(int)(tree->t_current-tree->t_beg),tree->c_pars);
3217#ifndef QUIET
3218  fflush(NULL);
3219#endif
3220}
3221
3222//////////////////////////////////////////////////////////////
3223//////////////////////////////////////////////////////////////
3224
3225void Print_Lk_And_Pars(t_tree *tree)
3226{       
3227  time(&(tree->t_current));
3228
3229  PhyML_Printf("\n. (%5d sec) [%15.4f] [%5d]",
3230         (int)(tree->t_current-tree->t_beg),
3231         tree->c_lnL,tree->c_pars);
3232
3233#ifndef QUIET
3234  fflush(NULL);
3235#endif
3236}
3237
3238//////////////////////////////////////////////////////////////
3239//////////////////////////////////////////////////////////////
3240
3241void Read_Qmat(phydbl *daa, phydbl *pi, FILE *fp)
3242{
3243  int i,j;
3244  phydbl sum;
3245  double val;
3246
3247  rewind(fp);
3248
3249  for(i=1;i<20;i++)
3250    {
3251      For(j,19)
3252        {
3253/*        if(!fscanf(fp,"%lf",&(daa[i*20+j]))) Exit("\n"); */
3254          if(!fscanf(fp,"%lf",&val)) 
3255            {
3256              PhyML_Printf("\n== Rate matrix file does not appear to have a proper format. Please refer to the documentation.");
3257              Exit("\n");
3258            }
3259          daa[i*20+j] = (phydbl)val;
3260          daa[j*20+i] = daa[i*20+j];
3261          if(j == i-1) break; 
3262        }
3263    }
3264
3265
3266  For(i,20) 
3267    { 
3268      if(!fscanf(fp,"%lf",&val)) Exit("\n");
3269      pi[i] = (phydbl)val;
3270    }
3271  sum = .0;
3272  For(i,20) sum += pi[i];
3273  if(FABS(sum - 1.) > 1.E-06)
3274    {
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;
3278    }
3279}
3280
3281//////////////////////////////////////////////////////////////
3282//////////////////////////////////////////////////////////////
3283
3284
3285void Print_Qmat_AA(phydbl *daa, phydbl *pi)
3286{
3287  int i,j,cpt;
3288
3289  cpt = 0;
3290  For(i,20)
3291    {
3292      for(j=0;j<i;j++)
3293        {
3294          PhyML_Printf("daa[%2d*20+%2d] = %10f;  ",i,j,daa[i*20+j]);
3295          cpt++;
3296          if(!(cpt%4)) PhyML_Printf("\n");
3297        }
3298    }
3299
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]);
3303  PhyML_Printf("\n");
3304  PhyML_Printf("Ala\tArg\tAsn\tAsp\tCys\tGln\tGlu\tGly\tHis\tIle\tLeu\tLys\tMet\tPhe\tPro\tSer\tThr\tTrp\tTyr\tVal\n");
3305}
3306
3307
3308//////////////////////////////////////////////////////////////
3309//////////////////////////////////////////////////////////////
3310
3311void Print_Square_Matrix_Generic(int n, phydbl *mat)
3312{
3313  int i,j;
3314
3315  PhyML_Printf("\n");
3316  For(i,n)
3317    {
3318      PhyML_Printf("[%3d]",i);
3319      For(j,n)
3320        {
3321          PhyML_Printf("%12.5f ",mat[i*n+j]);
3322        }
3323      PhyML_Printf("\n");
3324    }
3325  PhyML_Printf("\n");
3326}
3327
3328//////////////////////////////////////////////////////////////
3329//////////////////////////////////////////////////////////////
3330
3331
3332void Print_Diversity(FILE *fp, t_tree *tree)
3333{
3334 
3335  Print_Diversity_Pre(tree->a_nodes[0],
3336                      tree->a_nodes[0]->v[0],
3337                      tree->a_nodes[0]->b[0],
3338                      fp,
3339                      tree);
3340
3341/*       mean_div_left = .0; */
3342/*       For(k,ns)  */
3343/*      { */
3344/*        mean_div_left += (k+1) * tree->a_edges[j]->div_post_pred_left[k]; */
3345/*      } */
3346/*       mean_div_rght = .0; */
3347/*       For(k,ns) mean_div_rght += (k+1) * tree->a_edges[j]->div_post_pred_rght[k]; */
3348
3349/*       mean_div_left /= (phydbl)tree->data->init_len; */
3350/*       mean_div_rght /= (phydbl)tree->data->init_len; */
3351
3352/*       PhyML_Fprintf(fp,"%4d 0 %f\n",j,mean_div_left); */
3353/*       PhyML_Fprintf(fp,"%4d 1 %f\n",j,mean_div_rght); */
3354
3355
3356/*       mean_div_left = .0; */
3357/*       For(k,ns) mean_div_left += tree->a_edges[j]->div_post_pred_left[k]; */
3358
3359/*       mean_div_rght = .0; */
3360/*       For(k,ns)  */
3361/*      { */
3362/*        mean_div_rght += tree->a_edges[j]->div_post_pred_rght[k]; */
3363/*      } */
3364
3365/*       if((mean_div_left != tree->data->init_len) || (mean_div_rght != tree->data->init_len)) */
3366/*      { */
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(""); */
3371/*      } */
3372}
3373
3374//////////////////////////////////////////////////////////////
3375//////////////////////////////////////////////////////////////
3376
3377
3378void Print_Diversity_Pre(t_node *a, t_node *d, t_edge *b, FILE *fp, t_tree *tree)
3379{
3380  int k,ns;
3381
3382  ns = -1;
3383
3384  if(d->tax) return;
3385  else
3386    {
3387
3388      if(tree->io->datatype == NT)      ns = 4;
3389      else if(tree->io->datatype == AA) ns = 20;
3390
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]);
3393
3394      For(k,3) if(d->v[k] != a) Print_Diversity_Pre(d,d->v[k],d->b[k],fp,tree);
3395    }
3396
3397}
3398
3399//////////////////////////////////////////////////////////////
3400//////////////////////////////////////////////////////////////
3401
3402t_tree *Read_User_Tree(calign *cdata, t_mod *mod, option *io)
3403{
3404  t_tree *tree;
3405
3406 
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); 
3413  return tree;
3414}
3415
3416//////////////////////////////////////////////////////////////
3417//////////////////////////////////////////////////////////////
3418
3419
3420void Print_Time_Info(time_t t_beg, time_t t_end)
3421{
3422  div_t hour,min;
3423
3424  hour = div(t_end-t_beg,3600);
3425  min  = div(t_end-t_beg,60  );
3426  min.quot -= hour.quot*60;
3427
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");
3430}
3431
3432//////////////////////////////////////////////////////////////
3433//////////////////////////////////////////////////////////////
3434
3435void PhyML_Printf(char *format, ...)
3436{
3437  va_list ptr;
3438
3439 #ifdef MPI
3440  if(Global_myRank == 0)
3441    {
3442      va_start (ptr, format);
3443      vprintf (format, ptr);
3444      va_end(ptr);
3445    }
3446#else
3447      va_start (ptr, format);
3448      vprintf (format, ptr);
3449      va_end(ptr);
3450#endif
3451 
3452  fflush (NULL);
3453}
3454
3455//////////////////////////////////////////////////////////////
3456//////////////////////////////////////////////////////////////
3457
3458
3459void PhyML_Fprintf(FILE *fp, char *format, ...)
3460{
3461  va_list ptr;
3462
3463#ifdef MPI
3464  if(Global_myRank == 0)
3465    {
3466      va_start (ptr, format);
3467      vfprintf (fp,format, ptr);
3468      va_end(ptr);
3469    }
3470#else
3471      va_start (ptr, format);
3472      vfprintf (fp,format, ptr);
3473      va_end(ptr);
3474#endif
3475 
3476  fflush (NULL);
3477}
3478
3479//////////////////////////////////////////////////////////////
3480//////////////////////////////////////////////////////////////
3481
3482//////////////////////////////////////////////////////////////
3483//////////////////////////////////////////////////////////////
3484
3485void Read_Clade_Priors(char *file_name, t_tree *tree)
3486{
3487  FILE *fp;
3488  char *s,*line;
3489  int n_clade_priors;
3490  int clade_size;
3491  char **clade_list;
3492  int i,pos;
3493  phydbl prior_low,prior_up;
3494  int node_num;
3495
3496  PhyML_Printf("\n");
3497  PhyML_Printf("\n. Reading prior on node ages.\n");
3498
3499  line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3500  s    = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3501
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));
3504
3505  fp = Openfile(file_name,0);
3506 
3507  n_clade_priors = 0; 
3508  do
3509    {
3510      if(!fgets(line,T_MAX_LINE,fp)) break;
3511
3512      clade_size = 0;
3513      pos = 0;
3514      do
3515        {
3516          i = 0;
3517
3518          while(line[pos] == ' ') pos++;
3519
3520          while((line[pos] != ' ') && (line[pos] != '\n') && line[pos] != '#')
3521            {
3522              s[i] = line[pos];
3523              i++;
3524              pos++;
3525            }
3526          s[i] = '\0';
3527
3528          /* PhyML_Printf("\n. s = %s\n",s); */
3529         
3530          if(line[pos] == '\n' || line[pos] == '#') break;
3531          pos++;
3532
3533          if(strcmp(s,"|"))
3534            {
3535              strcpy(clade_list[clade_size],s);
3536              clade_size++;
3537            }
3538          else 
3539            break;         
3540        }
3541      while(1);
3542
3543     
3544      if(line[pos] != '#' && line[pos] != '\n')
3545        {
3546          phydbl val1, val2;
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;
3551          node_num = -1;
3552          if(!strcmp("@root@",clade_list[0])) node_num = tree->n_root->num;
3553          else node_num = Find_Clade(clade_list, clade_size, tree);
3554
3555          n_clade_priors++; 
3556
3557          if(node_num < 0)
3558            {
3559              PhyML_Printf("\n");
3560              PhyML_Printf("\n");
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. .................................................................");
3565              /* sleep(3); */
3566            }
3567          else
3568            {         
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);
3572
3573
3574              /* Sergei to add stuff here re calibration object */
3575
3576
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;
3579
3580              PhyML_Printf("\n");
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. .......................................................................");
3587            }
3588        }
3589    }
3590  while(1);
3591 
3592 
3593  PhyML_Printf("\n. Read prior information on %d %s.",n_clade_priors,n_clade_priors > 1 ? "clades":"clade");
3594
3595  if(!n_clade_priors)
3596    {
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__);
3604      Warn_And_Exit("");
3605    }
3606
3607  For(i,tree->n_otu) Free(clade_list[i]);
3608  Free(clade_list);
3609  Free(line);
3610  Free(s);
3611  fclose(fp);
3612}
3613
3614//////////////////////////////////////////////////////////////
3615//////////////////////////////////////////////////////////////
3616
3617option *Get_Input(int argc, char **argv)
3618{
3619  option *io;
3620  t_mod *mod;
3621  t_opt *s_opt;
3622  m4 *m4mod;
3623  int rv;
3624 
3625  rv = 1;
3626 
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();
3631
3632
3633  Set_Defaults_Input(io);
3634  Set_Defaults_Model(mod);
3635  Set_Defaults_Optimiz(s_opt);
3636  io->mod        = mod;
3637  io->mod->m4mod = m4mod;
3638  mod->io        = io;
3639  mod->s_opt     = s_opt;
3640
3641
3642
3643
3644#ifdef MPI
3645  rv = Read_Command_Line(io,argc,argv);
3646#elif (defined PHYTIME || defined SERGEII)
3647  rv = Read_Command_Line(io,argc,argv);
3648#else
3649  putchar('\n');
3650
3651  switch (argc)
3652    {
3653    case 1:
3654      {
3655        Launch_Interface(io);
3656        break;
3657      }
3658      /*
3659        case 2:
3660        Usage();
3661        break;
3662      */
3663    default:
3664      rv = Read_Command_Line(io,argc,argv);
3665    }
3666#endif
3667 
3668  if(rv) return io;
3669  else   return NULL;
3670}
3671
3672//////////////////////////////////////////////////////////////
3673//////////////////////////////////////////////////////////////
3674
3675int Set_Whichmodel(int select)
3676{
3677  int wm;
3678
3679  wm = -1;
3680
3681  switch(select)
3682    {
3683    case 1:
3684      {
3685        wm = JC69;
3686        break;
3687      }
3688    case 2:
3689      {
3690        wm = K80;
3691        break;
3692      }
3693    case 3:
3694      {
3695        wm = F81;
3696        break;
3697      }
3698    case 4:
3699      {
3700        wm = HKY85;
3701        break;
3702      }
3703    case 5:
3704      {
3705        wm = F84;
3706        break;
3707      }
3708    case 6:
3709      {
3710        wm = TN93;
3711        break;
3712      }
3713    case 7:
3714      {
3715        wm = GTR;
3716        break;
3717      }
3718    case 8:
3719      {
3720        wm = CUSTOM;
3721        break;
3722      }
3723    case 11:
3724      {
3725        wm = WAG;
3726        break;
3727      }
3728    case 12:
3729      {
3730        wm = DAYHOFF;
3731        break;
3732      }
3733    case 13:
3734      {
3735        wm = JTT;
3736        break;
3737      }
3738    case 14:
3739      {
3740        wm = BLOSUM62;
3741        break;
3742      }
3743    case 15:
3744      {
3745        wm = MTREV;
3746        break;
3747      }
3748    case 16:
3749      {
3750        wm = RTREV;
3751        break;
3752      }
3753    case 17:
3754      {
3755        wm = CPREV;
3756        break;
3757      }
3758    case 18:
3759      {
3760        wm = DCMUT;
3761        break;
3762      }
3763    case 19:
3764      {
3765        wm = VT;
3766        break;
3767      }
3768    case 20:
3769      {
3770        wm = MTMAM;
3771        break;
3772      }
3773    case 21:
3774      {
3775        wm = MTART;
3776        break;
3777      }
3778    case 22:
3779      {
3780        wm = HIVW;
3781        break;
3782      }
3783    case 23:
3784      {
3785        wm = HIVB;
3786        break;
3787      }
3788    case 24:
3789      {
3790        wm = CUSTOMAA;
3791        break;
3792      }
3793    case 25:
3794      {
3795        wm = LG;
3796        break;
3797      }
3798    default:
3799      {
3800        PhyML_Printf("\n== Model number %d is unknown. Please use a valid model name",select);
3801        Exit("\n");
3802        break;
3803      }
3804    }
3805
3806  return wm;
3807
3808}
3809
3810//////////////////////////////////////////////////////////////
3811//////////////////////////////////////////////////////////////
3812
3813void Print_Data_Structure(int final, FILE *fp, t_tree *mixt_tree)
3814{
3815  int n_partition_elem;
3816  char *s;
3817  t_tree *tree,*cpy_mixt_tree;
3818  int c,cc,cc_efrq,cc_rmat,cc_lens;
3819  char *param;
3820  int *link_efrq,*link_rmat,*link_lens;
3821  phydbl r_mat_weight_sum, e_frq_weight_sum;
3822 
3823 
3824  PhyML_Fprintf(fp,"\n. Starting tree: %s",
3825               mixt_tree->io->in_tree == 2?mixt_tree->io->in_tree_file:"BioNJ");
3826
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":
3831               "spr+nni":"no");
3832
3833  cpy_mixt_tree = mixt_tree;
3834 
3835  n_partition_elem = 1;
3836  tree = mixt_tree;
3837  do
3838    {
3839      tree = tree->next_mixt;
3840      if(!tree) break;
3841      n_partition_elem++;
3842    }
3843  while(1);
3844 
3845  s = (char *)mCalloc(2,sizeof(char));
3846  s[0] = ' ';
3847  s[1] = '\0';
3848  tree = mixt_tree;
3849  do
3850    {
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);
3853      strcat(s,", ");
3854      tree = tree->next_mixt;
3855      if(!tree) break;     
3856    }
3857  while(1);
3858  s[(int)strlen(s)-2]=' ';
3859  s[(int)strlen(s)-1]='\0';
3860
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);
3863  Free(s);
3864
3865  if(final == YES)
3866    PhyML_Fprintf(fp,"\n. Final log-likelihood: %f",mixt_tree->c_lnL);
3867
3868  do
3869    {
3870      int class = 0;
3871     
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");
3878     
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)
3881        {
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)
3886            {
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");
3889            }
3890          if(mixt_tree->mod->ras->invar == YES)
3891            {
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");
3894            }
3895        }
3896
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);
3899         
3900      tree = mixt_tree;
3901      do
3902        {
3903          if(tree->is_mixt_tree) tree = tree->next;
3904
3905          PhyML_Fprintf(fp,"\n");
3906          PhyML_Fprintf(fp,"\n. Mixture class %d",class+1);
3907
3908          if(mixt_tree->mod->ras->n_catg > 1)
3909            {
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);
3913            }
3914
3915          PhyML_Fprintf(fp,"\n   Substitution model:\t\t%12s",tree->mod->modelname->s);
3916         
3917          if(tree->mod->whichmodel == CUSTOM)
3918            PhyML_Fprintf(fp,"\n   Substitution model code:\t%12s",tree->mod->custom_mod_string->s);
3919         
3920          if(tree->mod->whichmodel == CUSTOMAA)
3921            PhyML_Fprintf(fp,"\n   Rate matrix file name:\t%12s",tree->mod->aa_rate_mat_file->s);
3922
3923          if(tree->mod->whichmodel == K80 ||
3924             tree->mod->whichmodel == HKY85 ||
3925             tree->mod->whichmodel == TN93)
3926            {
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");
3929            }
3930          else if(tree->mod->whichmodel == GTR ||
3931                  tree->mod->whichmodel == CUSTOM)         
3932            {
3933              PhyML_Fprintf(fp,"\n   Optimise subst. rates:\t%12s",tree->mod->s_opt->opt_rr?"yes":"no");
3934              if(final == YES)
3935                {
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]);
3942                }
3943            }
3944
3945
3946          PhyML_Fprintf(fp,"\n   Rate matrix weight:\t\t%12f",tree->mod->r_mat_weight->v  / r_mat_weight_sum);
3947
3948
3949
3950          if(tree->io->datatype == NT && 
3951             tree->mod->whichmodel != JC69 && 
3952             tree->mod->whichmodel != K80)
3953            {
3954              PhyML_Fprintf(fp,"\n   Optimise nucletide freq.:\t%12s",tree->mod->s_opt->opt_state_freq?"yes":"no");
3955              if(final == YES)
3956                {
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]);                 
3961                }
3962            }
3963          else if(tree->io->datatype == AA)
3964            {
3965              char *s;
3966             
3967              s = (char *)mCalloc(50,sizeof(char));
3968
3969              if(tree->mod->s_opt->opt_state_freq == YES)
3970                {
3971                  strcpy(s,"Empirical");
3972                }
3973              else
3974                {
3975                  strcpy(s,"Model");
3976                }
3977
3978              PhyML_Fprintf(fp,"\n   Amino-acid freq.:\t\t%12s",s);
3979
3980              Free(s);
3981            }
3982
3983
3984          PhyML_Fprintf(fp,"\n   Equ. freq. weight:\t\t%12f",tree->mod->e_frq_weight->v / e_frq_weight_sum);
3985
3986          class++;
3987         
3988          tree = tree->next;
3989
3990          if(tree && tree->is_mixt_tree == YES) break;
3991        }
3992      while(tree);
3993
3994      mixt_tree = mixt_tree->next_mixt;
3995      if(!mixt_tree) break;
3996    }
3997  while(1); 
3998
3999
4000  tree = cpy_mixt_tree;
4001  c = 0;
4002  do
4003    {
4004      if(tree->is_mixt_tree) tree = tree->next;
4005      c++;
4006      tree = tree->next;
4007    }
4008  while(tree);
4009
4010  link_efrq = (int *)mCalloc(c,sizeof(int));
4011  link_lens = (int *)mCalloc(c,sizeof(int));
4012  link_rmat = (int *)mCalloc(c,sizeof(int));
4013
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;
4025  c = 0;
4026  do
4027    {
4028      if(tree->is_mixt_tree) tree = tree->next;
4029      PhyML_Fprintf(fp,"---");
4030      tree = tree->next;
4031    }
4032  while(tree);
4033
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;
4039  c = 0;
4040  do
4041    {
4042      if(tree->is_mixt_tree) tree = tree->next;
4043      PhyML_Fprintf(fp,"%2d ",tree->mixt_tree->dp);
4044      tree = tree->next;
4045    }
4046  while(tree);
4047
4048
4049  PhyML_Fprintf(fp,"\n");
4050  PhyML_Fprintf(fp,"  ------------------");
4051  tree = cpy_mixt_tree;
4052  c = 0;
4053  do
4054    {
4055      if(tree->is_mixt_tree) tree = tree->next;
4056      PhyML_Fprintf(fp,"---");
4057      tree = tree->next;
4058    }
4059  while(tree);
4060
4061 
4062  tree = cpy_mixt_tree;
4063  c    = 0;
4064  do
4065    {
4066      if(tree->is_mixt_tree) tree = tree->next;
4067      link_rmat[c] = -1;
4068      link_lens[c] = -1;
4069      link_efrq[c] = -1;
4070      tree = tree->next;
4071      c++;
4072    }
4073  while(tree);
4074
4075  mixt_tree = cpy_mixt_tree;
4076  cc_efrq   = 97; 
4077  cc_rmat   = 97; 
4078  cc_lens   = 97; 
4079  cc        = 0;
4080  do
4081    {
4082      if(mixt_tree->is_mixt_tree) mixt_tree = mixt_tree->next;
4083
4084      if(link_efrq[cc] < 0)
4085        {
4086          link_efrq[cc] = cc_efrq;
4087          tree = mixt_tree->next;
4088          c    = cc+1;
4089          if(tree)
4090            {
4091              do
4092                {
4093                  if(tree->is_mixt_tree) tree = tree->next;
4094
4095                  if(mixt_tree->mod->e_frq == tree->mod->e_frq) link_efrq[c] = cc_efrq;
4096
4097                  tree = tree->next;
4098                  c++;
4099                }
4100              while(tree);
4101            }
4102          cc_efrq++;
4103        }
4104
4105      if(link_lens[cc] < 0)
4106        {
4107          link_lens[cc] = cc_lens;
4108          tree = mixt_tree->next;
4109          c    = cc+1;
4110          if(tree)
4111            {
4112              do
4113                {
4114                  if(tree->is_mixt_tree) tree = tree->next;
4115
4116                  if(mixt_tree->a_edges[0]->l == tree->a_edges[0]->l) link_lens[c] = cc_lens;
4117                 
4118                  tree = tree->next;
4119                  c++;
4120                }
4121              while(tree);
4122            }
4123          cc_lens++;
4124        }
4125
4126      if(link_rmat[cc] < 0)
4127        {
4128          link_rmat[cc] = cc_rmat;
4129          tree = mixt_tree->next;
4130          c    = cc+1;
4131          if(tree)
4132            {
4133              do
4134                {
4135                  if(tree->is_mixt_tree) tree = tree->next;
4136
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;
4143
4144                  tree = tree->next;
4145                  c++;
4146                }
4147              while(tree);
4148            }
4149          cc_rmat++;
4150        }
4151     
4152      cc++;
4153      mixt_tree = mixt_tree->next;
4154    }
4155  while(mixt_tree);
4156   
4157  PhyML_Fprintf(fp,"\n");
4158  strcpy(param,"State frequencies ");
4159  PhyML_Fprintf(fp,"  %18s",param);
4160
4161  tree = cpy_mixt_tree;
4162  c    = 0;
4163  do
4164    {
4165      if(tree->is_mixt_tree) tree = tree->next;
4166      PhyML_Fprintf(fp,"%2c ",link_efrq[c]);
4167      tree = tree->next;
4168      c++;
4169    }
4170  while(tree);
4171
4172
4173  PhyML_Fprintf(fp,"\n");
4174  strcpy(param,"Branch lengths ");
4175  PhyML_Fprintf(fp,"  %18s",param);
4176
4177  tree = cpy_mixt_tree;
4178  c    = 0;
4179  do
4180    {
4181      if(tree->is_mixt_tree) tree = tree->next;
4182      PhyML_Fprintf(fp,"%2c ",link_lens[c]);
4183      tree = tree->next;
4184      c++;
4185    }
4186  while(tree);
4187
4188  PhyML_Fprintf(fp,"\n");
4189  strcpy(param,"Rate matrix ");
4190  PhyML_Fprintf(fp,"  %18s",param);
4191
4192  tree = cpy_mixt_tree;
4193  c    = 0;
4194  do
4195    {
4196      if(tree->is_mixt_tree) tree = tree->next;
4197      PhyML_Fprintf(fp,"%2c ",link_rmat[c]);
4198      tree = tree->next;
4199      c++;
4200    }
4201  while(tree);
4202
4203  PhyML_Fprintf(fp,"\n");
4204  PhyML_Fprintf(fp,"  ------------------");
4205  tree = cpy_mixt_tree;
4206  c = 0;
4207  do
4208    {
4209      if(tree->is_mixt_tree) tree = tree->next;
4210      PhyML_Fprintf(fp,"---");
4211      tree = tree->next;
4212    }
4213  while(tree);
4214  PhyML_Fprintf(fp,"\n");
4215 
4216  if(final == YES)
4217    {
4218      tree = cpy_mixt_tree;
4219      c = 0;
4220      do
4221        {
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);
4227          Free(s);
4228          tree = tree->next_mixt;
4229        }
4230      while(tree);
4231    }
4232
4233  Free(param);
4234  Free(link_efrq);
4235  Free(link_rmat);
4236  Free(link_lens);
4237
4238  mixt_tree = cpy_mixt_tree;
4239}
4240
4241//////////////////////////////////////////////////////////////
4242//////////////////////////////////////////////////////////////
4243
4244void PhyML_XML(char *xml_filename)
4245{
4246  FILE *fp;
4247  xml_node *root,*p_elem,*m_elem,*parent,*instance;
4248  option *io;
4249  void *buff;
4250  t_mod *mod,*iomod;
4251  t_tree *tree,*mixt_tree,*root_tree;
4252  int select;
4253  char *component;
4254  int i,j,n_components;
4255  int first_m_elem;
4256  int class_number;
4257  scalar_dbl **lens,**lens_var,**ori_lens,**lens_old,**lens_var_old,**ori_lens_old,**ori_lens_var,**ori_lens_var_old;
4258  t_ds *ds;
4259  char *outputfile;
4260  char *alignment;
4261  char *s;
4262  int lens_size;
4263  int first;
4264  int *class_num;
4265
4266
4267  fp = fopen(xml_filename,"r");
4268  if(!fp)
4269    {
4270      PhyML_Printf("\n== Could not find the XML file '%s'.\n",xml_filename);
4271      Exit("\n");
4272    }
4273
4274
4275  root = XML_Load_File(fp); 
4276 
4277  if(!root)
4278    {
4279      PhyML_Printf("\n== Encountered an issue while loading the XML file.\n");
4280      Exit("\n");
4281    }
4282
4283  class_num = (int *)mCalloc(N_MAX_MIXT_CLASSES,sizeof(int));
4284  For(i,N_MAX_MIXT_CLASSES) class_num[i] = i;
4285
4286  component = (char *)mCalloc(T_MAX_NAME,sizeof(char));
4287
4288  m_elem           = NULL;
4289  p_elem           = root;
4290  io               = NULL;
4291  mixt_tree        = NULL;
4292  root_tree        = NULL;
4293  mod              = NULL;
4294  tree             = NULL;
4295  lens             = NULL;
4296  lens_var         = NULL;
4297  lens_old         = NULL;
4298  lens_var_old     = NULL;
4299  select           = -1;
4300  class_number     = -1;
4301  ds               = NULL;
4302  lens_size        = 0;
4303  ori_lens         = NULL;
4304  ori_lens_old     = NULL;
4305  ori_lens_var     = NULL;
4306  ori_lens_var_old = NULL;
4307  first            = YES;
4308
4309  // Make sure there are no duplicates in node's IDs
4310  XML_Check_Duplicate_ID(root);
4311
4312  int count = 0;
4313  XML_Count_Number_Of_Node_With_Name("topology",&count,root);
4314 
4315  if(count > 1)
4316    {
4317      PhyML_Printf("\n== There should not more than one 'topology' node.");
4318      PhyML_Printf("\n== Found %d. Please fix your XML file",count);
4319      Exit("\n");
4320    }
4321  else if(count < 1)
4322    {
4323      PhyML_Printf("\n== There should be at least one 'topology' node.");
4324      PhyML_Printf("\n== Found none. Please fix your XML file");
4325      Exit("\n");
4326    }
4327 
4328  p_elem = XML_Search_Node_Name("phyml",NO,p_elem);
4329 
4330  /*! Input file
4331   */
4332  outputfile = XML_Get_Attribute_Value(p_elem,"output.file"); 
4333
4334  if(!outputfile)
4335    {
4336      PhyML_Printf("\n== The 'outputfile' attribute in 'phyml' tag is mandatory.");
4337      PhyML_Printf("\n== Please amend your XML file accordingly.");
4338      Exit("\n");
4339    }
4340
4341  io = (option *)Make_Input();
4342  Set_Defaults_Input(io);
4343
4344  s = XML_Get_Attribute_Value(p_elem,"run.id");
4345  if(s)
4346    {
4347      io->append_run_ID = YES;
4348      strcpy(io->run_id_string,s);
4349    }
4350
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);         
4359
4360  s = XML_Get_Attribute_Value(p_elem,"print.trace");
4361 
4362  if(s)
4363    {
4364      select = XML_Validate_Attr_Int(s,6,
4365                                     "true","yes","y",
4366                                     "false","no","n");
4367     
4368      if(select < 3) 
4369        {
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);
4375        }
4376    }
4377
4378
4379  s = XML_Get_Attribute_Value(p_elem,"branch.test");
4380  if(s)
4381    {
4382      if(!strcmp(s,"aLRT"))
4383        {
4384          io->ratio_test = ALRTSTAT;
4385        }
4386      else if(!strcmp(s,"aBayes"))
4387        {
4388          io->ratio_test = ABAYES;
4389        }
4390      else if(!strcmp(s,"SH"))
4391        {
4392          io->ratio_test = SH;
4393        }
4394      else if(!strcmp(s,"no") || !strcmp(s,"none"))
4395        {
4396          io->ratio_test = 0;
4397        }
4398      else
4399        {
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'.");
4402          Exit("\n");
4403        }     
4404    }
4405   
4406
4407  /*! Read all partitionelem nodes and mixturelem nodes in each of them
4408   */
4409  do
4410    {
4411      p_elem = XML_Search_Node_Name("partitionelem",YES,p_elem);
4412
4413      if(p_elem == NULL) break;
4414
4415      buff = (option *)Make_Input();
4416      Set_Defaults_Input(buff);
4417      io->next = buff;
4418      io->next->prev = io;
4419     
4420      io = io->next;           
4421      if(first == YES) 
4422        {
4423          io = io->prev;
4424          io->next = NULL;
4425          Free_Input(buff);
4426          first = NO;
4427        }
4428
4429
4430      /*! Set the datatype (required when compressing data)
4431       */
4432      char *dt = NULL;
4433      dt = XML_Get_Attribute_Value(p_elem,"data.type");
4434      if(!dt)
4435        {
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\"'");
4439          Exit("\n");       
4440        }
4441
4442      select = XML_Validate_Attr_Int(dt,2,"aa","nt");
4443      switch(select)
4444        {
4445        case 0: 
4446          {
4447            io->datatype = AA;
4448            break;
4449          }
4450        case 1:
4451          {
4452            io->datatype = NT;
4453            break;
4454          }
4455        default:
4456          {
4457            PhyML_Printf("\n== Unknown data type. Must be either 'aa' or 'nt'.");
4458            Exit("\n");
4459          }
4460        }
4461     
4462      char *format = NULL;
4463      format = XML_Get_Attribute_Value(p_elem,"format");
4464      if(format)
4465        {
4466          if(!strcmp(format,"interleave") ||
4467             !strcmp(format,"interleaved"))
4468            {
4469              io->interleaved = YES;
4470            }
4471          else if(!strcmp(format,"sequential"))
4472            {
4473              io->interleaved = NO;             
4474            }
4475        }
4476
4477
4478      /*! Attach a model to this io struct
4479       */
4480      io->mod = (t_mod *)Make_Model_Basic();
4481      Set_Defaults_Model(io->mod);
4482      io->mod->ras->n_catg = 1;
4483      io->mod->io = io;     
4484      iomod = io->mod;
4485
4486      /*! Attach an optimization structure to this model
4487       */
4488      iomod->s_opt = (t_opt *)Make_Optimiz();
4489      Set_Defaults_Optimiz(iomod->s_opt);
4490
4491      iomod->s_opt->opt_kappa   = NO;
4492      iomod->s_opt->opt_lambda  = NO;
4493      iomod->s_opt->opt_rr      = NO;
4494
4495      /*! Input file
4496       */
4497      alignment = XML_Get_Attribute_Value(p_elem,"file.name"); 
4498
4499      if(!alignment)
4500        {
4501          PhyML_Printf("\n== 'file.name' tag is mandatory. Please amend your");
4502          PhyML_Printf("\n== XML file accordingly.");
4503          Exit("\n");
4504        }
4505
4506      strcpy(io->in_align_file,alignment);
4507     
4508      /*! Open pointer to alignment
4509       */
4510      io->fp_in_align = Openfile(io->in_align_file,0);
4511
4512      /*! Load sequence file
4513       */
4514      io->data  = Get_Seq(io);
4515
4516      /*! Close pointer to alignment
4517       */
4518      fclose(io->fp_in_align);
4519
4520      /*! Compress alignment
4521       */
4522      io->cdata = Compact_Data(io->data,io);
4523     
4524      /*! Free uncompressed alignment
4525       */
4526      Free_Seq(io->data,io->n_otu);
4527
4528      /*! Create new mixture tree
4529       */
4530      buff = (t_tree *)Make_Tree_From_Scratch(io->cdata->n_otu,io->cdata);
4531
4532      if(mixt_tree)
4533        {
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;
4538        }
4539      else mixt_tree = buff;
4540     
4541      /*! Connect mixt_tree to io struct
4542       */
4543      mixt_tree->io = io;
4544
4545      /*! Connect mixt_tree to model struct
4546       */
4547      mixt_tree->mod = iomod;
4548
4549      /*! mixt_tree is a mixture tree
4550       */
4551      mixt_tree->is_mixt_tree = YES;
4552
4553      /*! mixt_tree is a mixture tree
4554       */
4555      mixt_tree->mod->is_mixt_mod = YES;
4556
4557      /*! Connect mixt_tree to compressed data
4558       */
4559      mixt_tree->data = io->cdata;
4560
4561      /*! Set total number of patterns
4562       */
4563      mixt_tree->n_pattern = io->cdata->crunch_len;
4564
4565      /*! Remove branch lengths from mixt_tree */
4566      For(i,2*mixt_tree->n_otu-1)
4567        {
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);
4572        }
4573
4574      /*! Connect last tree of the mixture for the
4575        previous partition element to the next mixture tree
4576      */
4577      if(tree) 
4578        {
4579          tree->next = mixt_tree;
4580          mixt_tree->prev = tree;
4581        }
4582
4583      /*! Do the same for the model
4584       */
4585      if(mod) 
4586        {
4587          mod->next = iomod;
4588          iomod->prev = mod;
4589        }
4590
4591      if(!root_tree) root_tree = mixt_tree;
4592
4593      /*! Process all the mixtureelem tags in this partition element
4594       */
4595      n_components  = 0;
4596      m_elem        = p_elem;
4597      first_m_elem  = 0;
4598      mod           = NULL;
4599      tree          = NULL;
4600      class_number  = 0;
4601      do
4602        {
4603          m_elem = XML_Search_Node_Name("mixtureelem",YES,m_elem);
4604          if(m_elem == NULL) break;
4605         
4606          if(!strcmp(m_elem->name,"mixtureelem"))
4607            {
4608              first_m_elem++;
4609             
4610              /*! Rewind tree and model when processing a new mixtureelem node
4611               */
4612              if(first_m_elem > 1) 
4613                {
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;
4616                }
4617
4618              /*! Read and process model components
4619               */
4620              char *list;
4621              list = XML_Get_Attribute_Value(m_elem,"list");
4622
4623              j = 0;
4624              For(i,(int)strlen(list)) if(list[i] == ',') j++;
4625             
4626              if(j != n_components && first_m_elem > 1)
4627                {
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);
4630                  Exit("\n");
4631                }
4632              n_components = j;
4633
4634              i = j = 0;
4635              component[0] = '\0';           
4636              while(1)
4637              {
4638                if(list[j] == ',' || j == (int)strlen(list))
4639                  {
4640                    /*! Reading a new component
4641                     */
4642
4643                    if(first_m_elem == YES) // Only true when processing the first mixtureelem node
4644                      {
4645                        t_tree *this_tree;
4646                        t_mod *this_mod;
4647
4648                        /*! Create new tree
4649                         */
4650                        this_tree = (t_tree *)Make_Tree_From_Scratch(io->cdata->n_otu,io->cdata);
4651                                               
4652                        /*! Update the number of mixtures */
4653                        iomod->n_mixt_classes++;
4654
4655                        if(tree)
4656                          {
4657                            tree->next = this_tree;
4658                            tree->next->prev = tree;
4659                          }
4660                        else 
4661                          {
4662                            mixt_tree->next = this_tree;
4663                            mixt_tree->next->prev = mixt_tree;
4664                          }
4665
4666                        tree = this_tree;
4667                        tree->mixt_tree = mixt_tree;
4668
4669
4670                        /*! Create a new model
4671                         */
4672                        this_mod = (t_mod *)Make_Model_Basic();
4673                        Set_Defaults_Model(this_mod);
4674                        this_mod->ras->n_catg = 1;
4675                       
4676                       
4677
4678                        if(mod)
4679                          {
4680                            mod->next = this_mod;
4681                            mod->next->prev = mod;
4682                          }
4683                        else
4684                          {
4685                            this_mod->prev = iomod;
4686                          }
4687
4688                        mod = this_mod;
4689                        if(!iomod->next) iomod->next = mod;
4690                        mod->io = io;
4691
4692                        mod->s_opt = (t_opt *)Make_Optimiz();
4693                        Set_Defaults_Optimiz(mod->s_opt);
4694                       
4695                        mod->s_opt->opt_alpha  = NO;
4696                        mod->s_opt->opt_pinvar = NO;
4697                       
4698                        tree->data      = io->cdata;
4699                        tree->n_pattern = io->cdata->crunch_len;
4700                        tree->io        = io;
4701                        tree->mod       = mod;
4702
4703                        if(tree->n_pattern != tree->prev->n_pattern)
4704                          {
4705                            PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
4706                            Warn_And_Exit("");
4707                          }
4708                      }
4709
4710                    /*! Read a component
4711                     */
4712                    component[i] = '\0';                   
4713                    if(j != (int)strlen(list)-1) i = 0;
4714                   
4715                    /*! Find which node this ID corresponds to
4716                     */
4717                    instance = XML_Search_Node_ID(component,YES,root);
4718                 
4719                    if(!instance)
4720                      {
4721                        PhyML_Printf("\n== Could not find a node with id:'%s'.",component);
4722                        PhyML_Printf("\n== Problem with 'mixtureelem' node, list '%s'.",list);
4723                        Exit("\n");
4724                      }
4725
4726                    if(!instance->parent)
4727                      {
4728                        PhyML_Printf("\n== Node '%s' with id:'%s' has no parent.",instance->name,component);
4729                        Exit("\n");
4730                      }
4731                     
4732                    parent = instance->parent;
4733
4734                    ////////////////////////////////////////
4735                    //        SUBSTITUTION MODEL          //
4736                    ////////////////////////////////////////
4737
4738                    if(!strcmp(parent->name,"ratematrices"))
4739                      {
4740                        /* ! First time we process this 'instance' node which has this 'ratematrices' parent */
4741                        if(instance->ds->obj == NULL)
4742                          {
4743                            Make_Ratematrice_From_XML_Node(instance, io, mod);
4744
4745                            ds = instance->ds;
4746
4747                            /*! Connect the data structure n->ds to mod->r_mat */
4748                            ds->obj = (t_rmat *)(mod->r_mat);
4749                           
4750                            /*! Create and connect the data structure n->ds->next to mod->kappa */
4751                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4752                            ds = ds->next;
4753                            ds->obj = (scalar_dbl *)(mod->kappa);
4754
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));
4757                            ds = ds->next;
4758                            ds->obj = (int *)(&mod->s_opt->opt_kappa);
4759
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));
4762                            ds = ds->next;
4763                            ds->obj = (int *)(&mod->s_opt->opt_rr);
4764
4765                            /*! Create and connect the data structure n->ds->next to mod->whichmodel */
4766                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4767                            ds = ds->next;
4768                            ds->obj = (int *)(&mod->whichmodel);
4769
4770                            /*! Create and connect the data structure n->ds->next to mod->modelname */
4771                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4772                            ds = ds->next;
4773                            ds->obj = (t_string *)(mod->modelname);
4774
4775                            /*! Create and connect the data structure n->ds->next to mod->ns */
4776                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4777                            ds = ds->next;
4778                            ds->obj = (int *)(&mod->ns);
4779
4780                            /*! Create and connect the data structure n->ds->next to mod->modelname */
4781                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4782                            ds = ds->next;
4783                            ds->obj = (t_string *)(mod->custom_mod_string);
4784                          }
4785                        else
4786                          {
4787                            /*! Connect to already extisting r_mat & kappa structs. */
4788                            t_ds *ds;
4789
4790                            ds = instance->ds;
4791                            Free(mod->r_mat);
4792                            mod->r_mat             = (t_rmat *)ds->obj;
4793
4794                            ds = ds->next;
4795                            Free_Scalar_Dbl(mod->kappa);
4796                            mod->kappa             = (scalar_dbl *)ds->obj;
4797
4798                            ds = ds->next;
4799                            mod->s_opt->opt_kappa  = *((int *)ds->obj);
4800
4801                            ds = ds->next;
4802                            mod->s_opt->opt_rr     = *((int *)ds->obj);
4803
4804                            ds = ds->next;
4805                            mod->whichmodel        = *((int *)ds->obj);
4806
4807                            ds = ds->next;
4808                            Free_String(mod->modelname);
4809                            mod->modelname         = (t_string *)ds->obj;
4810
4811                            ds = ds->next;
4812                            mod->ns                = *((int *)ds->obj);
4813
4814                            ds = ds->next;
4815                            Free_String(mod->custom_mod_string);
4816                            mod->custom_mod_string = (t_string *)ds->obj;
4817                          }
4818                      }
4819
4820                    ////////////////////////////////////////
4821                    //           STATE FREQS              //
4822                    ////////////////////////////////////////
4823                       
4824                    else if(!strcmp(parent->name,"equfreqs"))
4825                      {                       
4826
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) 
4830                          {
4831                            Make_Efrq_From_XML_Node(instance,io,mod);
4832                           
4833                            t_ds *ds;
4834
4835                            ds = instance->ds;
4836                            ds->obj = (t_efrq *)mod->e_frq;
4837                           
4838                            ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4839                            ds = ds->next;
4840                            ds->obj = (int *)(&mod->s_opt->opt_state_freq);
4841
4842                            ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4843                            ds = ds->next;
4844                            ds->obj = (int *)(&mod->s_opt->user_state_freq);
4845
4846                            ds->next = (t_ds *)mCalloc(1,sizeof(t_ds));
4847                            ds = ds->next;
4848                            ds->obj = (vect_dbl *)(mod->user_b_freq);                           
4849                          }
4850                        else
4851                          {
4852                            // Connect the data structure n->ds to mod->e_frq
4853                           
4854                            ds = instance->ds;
4855                            mod->e_frq = (t_efrq *)ds->obj;
4856
4857                            ds = ds->next;
4858                            mod->s_opt->opt_state_freq = *((int *)ds->obj);
4859
4860                            ds = ds->next;
4861                            mod->s_opt->user_state_freq = *((int *)ds->obj);
4862
4863                            ds = ds->next;
4864                            Free_Vect_Dbl(mod->user_b_freq);
4865                            mod->user_b_freq = (vect_dbl *)ds->obj;
4866                          }
4867                      }
4868                       
4869                    //////////////////////////////////////////
4870                    //             TOPOLOGY                 //
4871                    //////////////////////////////////////////
4872                   
4873                    else if(!strcmp(parent->name,"topology"))
4874                      {
4875                        if(parent->ds->obj == NULL) Make_Topology_From_XML_Node(instance,io,iomod);
4876
4877                        ds = parent->ds;
4878                       
4879                        int buff;
4880                        ds->obj = (int *)(& buff);
4881                      }
4882
4883                    //////////////////////////////////////////
4884                    //                RAS                   //
4885                    //////////////////////////////////////////
4886
4887                    else if(!strcmp(parent->name,"siterates"))
4888                      {
4889                        char *rate_value = NULL;                       
4890                        /* scalar_dbl *r; */
4891                        phydbl val;
4892
4893                        /*! First time we process this 'siterates' node, check that its format is valid.
4894                          and process it afterwards.
4895                        */
4896                        if(parent->ds->obj == NULL)
4897                          {
4898                            class_number = 0;
4899
4900                            Make_RAS_From_XML_Node(parent,iomod);
4901
4902                            ds = parent->ds;
4903
4904                            ds->obj = (t_ras *)iomod->ras;
4905
4906                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4907                            ds = ds->next;
4908                            ds->obj = (int *)(&iomod->s_opt->opt_alpha);
4909
4910                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));
4911                            ds = ds->next;
4912                            ds->obj = (int *)(&iomod->s_opt->opt_free_mixt_rates);
4913                          }
4914                        else /*! Connect ras struct to already defined one. Same for opt_alpha & opt_free_mixt_rates */
4915                          {
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);
4920                          }
4921                       
4922                        rate_value = XML_Get_Attribute_Value(instance,"init.value");
4923
4924                        val = 1.;
4925                        if(rate_value) val = String_To_Dbl(rate_value);
4926                       
4927                        if(instance->ds->obj == NULL) 
4928                          {
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);
4932
4933                            iomod->ras->gamma_rr->v[class_number] = val;
4934                            iomod->ras->gamma_rr_unscaled->v[class_number] = val;
4935
4936                            iomod->ras->init_rr = NO;
4937
4938                            if(Are_Equal(val,0.0,1E-20) == NO) class_number++;
4939                          }
4940                                             
4941
4942                        /*! Note: ras is already connected to the relevant t_ds stucture. No need
4943                          to connect ras->gamma_rr or ras->p_invar */
4944
4945                        /*! Invariant */
4946                        if(Are_Equal(val,0.0,1E-20)) 
4947                          {
4948                            mod->ras->invar = YES;
4949                          }
4950                        else
4951                          {
4952                            mod->ras->parent_class_number = *((int *)instance->ds->next->obj);
4953                          }
4954                       
4955                        xml_node *orig_w = NULL;
4956                        orig_w = XML_Search_Node_Attribute_Value("appliesto",instance->id,YES,instance->parent);
4957                       
4958                        if(orig_w)
4959                          {
4960                            char *weight;
4961                            weight = XML_Get_Attribute_Value(orig_w,"value");
4962                            if(mod->ras->invar == YES)
4963                              {
4964                                iomod->ras->pinvar->v = String_To_Dbl(weight);
4965                              }
4966                            else
4967                              {
4968                                int class;
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;
4972                              }
4973                          }
4974                      }
4975
4976                    //////////////////////////////////////////////
4977                    //           BRANCH LENGTHS                 //
4978                    //////////////////////////////////////////////
4979
4980                    else if(!strcmp(parent->name,"branchlengths"))
4981                      {
4982                        int i;
4983                        int n_otu;
4984       
4985                        n_otu = tree->n_otu;
4986                       
4987                        if(instance->ds->obj == NULL)
4988                          {
4989                            if(!lens)                               
4990                              {                               
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 *));
4993
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 *));
4996
4997                                lens         = ori_lens;
4998                                lens_old     = ori_lens_old;
4999
5000                                lens_var     = ori_lens_var;
5001                                lens_var_old = ori_lens_var_old;
5002
5003                                lens_size = 2*tree->n_otu-1;
5004                              }
5005                            else
5006                              {
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 *));
5009
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 *));
5012
5013                                lens         = ori_lens     + lens_size;;
5014                                lens_old     = ori_lens_old + lens_size;;
5015
5016                                lens_var     = ori_lens_var     + lens_size;;
5017                                lens_var_old = ori_lens_var_old + lens_size;;
5018
5019                                lens_size += 2*tree->n_otu-1;
5020                              }
5021
5022                            For(i,2*tree->n_otu-1)
5023                              {
5024                                lens[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5025                                Init_Scalar_Dbl(lens[i]);
5026
5027                                lens_old[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5028                                Init_Scalar_Dbl(lens_old[i]);
5029
5030                                lens_var[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5031                                Init_Scalar_Dbl(lens_var[i]);
5032
5033                                lens_var_old[i] = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
5034                                Init_Scalar_Dbl(lens_var_old[i]);
5035
5036                                Free_Scalar_Dbl(tree->a_edges[i]->l);
5037                                Free_Scalar_Dbl(tree->a_edges[i]->l_old);
5038
5039                                Free_Scalar_Dbl(tree->a_edges[i]->l_var);
5040                                Free_Scalar_Dbl(tree->a_edges[i]->l_var_old);
5041
5042                                if(tree->prev && 
5043                                   tree->prev->a_edges[i]->l == mixt_tree->a_edges[i]->l &&
5044                                   tree->prev->is_mixt_tree == NO)
5045                                  {
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.");
5049                                    Exit("\n");
5050                                  }
5051                              }
5052                                                   
5053                            char *opt_bl = NULL;
5054                            opt_bl = XML_Get_Attribute_Value(instance,"optimise.lens");
5055                           
5056                            if(opt_bl)
5057                              {
5058                                if(!strcmp(opt_bl,"yes") || !strcmp(opt_bl,"true"))
5059                                  {
5060                                    iomod->s_opt->opt_bl = YES;
5061                                  }
5062                                else
5063                                  {
5064                                    iomod->s_opt->opt_bl = NO;
5065                                  }
5066                              }
5067
5068                            ds = instance->ds;
5069
5070                            ds->obj       = (scalar_dbl **)lens;
5071
5072                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));                           
5073                            ds            = ds->next;
5074                            ds->obj       = (scalar_dbl **)lens_old;
5075
5076                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));                           
5077                            ds            = ds->next;
5078                            ds->obj       = (scalar_dbl **)lens_var;
5079
5080                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));                           
5081                            ds            = ds->next;
5082                            ds->obj       = (scalar_dbl **)lens_var_old;
5083
5084                            ds->next      = (t_ds *)mCalloc(1,sizeof(t_ds));                           
5085                            ds            = ds->next;
5086                            ds->obj       = (int *)(&iomod->s_opt->opt_bl);
5087                          }
5088                        else
5089                          {
5090                            For(i,2*tree->n_otu-1) 
5091                              {
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);
5096                              }
5097
5098                            ds = instance->ds;
5099
5100                            lens     = (scalar_dbl **)ds->obj;
5101
5102                            ds = ds->next;
5103                            lens_old = (scalar_dbl **)ds->obj;
5104
5105                            ds = ds->next;
5106                            lens_var = (scalar_dbl **)ds->obj;
5107
5108                            ds = ds->next;
5109                            lens_var_old = (scalar_dbl **)ds->obj;
5110
5111                            ds = ds->next;
5112                            iomod->s_opt->opt_bl = *((int *)ds->obj);
5113                          }
5114                       
5115                        if(n_otu != tree->n_otu)
5116                          {
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);
5119                            Exit("\n");
5120                          }
5121                       
5122                        For(i,2*tree->n_otu-1) 
5123                          {
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];
5128
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];
5133                          }
5134                      }
5135
5136                    ///////////////////////////////////////////////
5137                    ///////////////////////////////////////////////
5138                    ///////////////////////////////////////////////
5139
5140                    if(first_m_elem > 1) // Done with this component, move to the next tree and model
5141                      {
5142                        if(tree->next) tree = tree->next;
5143                        if(mod->next)   mod  = mod->next;
5144                      }
5145
5146                  }
5147              else if(list[j] != ' ')
5148                  {
5149                    component[i] = list[j];
5150                    i++;
5151                  }
5152                j++;
5153                if(j == (int)strlen(list)+1) break;
5154
5155              } // end of mixtureelem processing
5156            } // end of partitionelem processing         
5157        }
5158      while(1);
5159    }
5160  while(1);
5161
5162
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);
5167 
5168  while(io->prev != NULL) io = io->prev;
5169  while(mixt_tree->prev != NULL) mixt_tree = mixt_tree->prev;
5170
5171  /*! Finish making the models */
5172  mod = mixt_tree->mod;
5173  do
5174    {
5175      Make_Model_Complete(mod);     
5176      mod = mod->next;
5177    }
5178  while(mod);
5179 
5180  Check_Taxa_Sets(mixt_tree);
5181
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");
5190 
5191  if(!io->mod->s_opt->random_input_tree) io->mod->s_opt->n_rand_starts = 1;
5192
5193  char *most_likely_tree = NULL;
5194  phydbl best_lnL = UNLIKELY;
5195  int num_rand_tree;
5196  For(num_rand_tree,io->mod->s_opt->n_rand_starts)
5197    {     
5198      /*! Initialize the models */
5199      mod = mixt_tree->mod;
5200      do
5201        {
5202          Init_Model(mod->io->cdata,mod,io);
5203          mod = mod->next;
5204        }
5205      while(mod);
5206
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]); */
5212      /* Exit("\n"); */
5213
5214      Print_Data_Structure(NO,stdout,mixt_tree);
5215
5216      t_tree *bionj_tree = NULL;     
5217      switch(mixt_tree->io->in_tree)
5218        {     
5219        case 2: /*! user-defined input tree */
5220          {
5221            if(!mixt_tree->io->fp_in_tree)
5222              {
5223                PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
5224                Exit("\n");
5225              }
5226           
5227            /*!
5228              Copy the user tree to all the tree structures
5229            */
5230            tree = Read_User_Tree(mixt_tree->io->cdata,
5231                                  mixt_tree->mod,
5232                                  mixt_tree->io);
5233            break;
5234          }
5235        case 1: case 0:         
5236          {
5237            if(!bionj_tree)
5238              {
5239                /*! Build a BioNJ tree from the analysis of
5240                  the first partition element
5241                */
5242                tree = Dist_And_BioNJ(mixt_tree->data,
5243                                      mixt_tree->mod,
5244                                      mixt_tree->io);
5245                bionj_tree = (t_tree *)Make_Tree_From_Scratch(mixt_tree->n_otu,mixt_tree->data);
5246                Copy_Tree(tree,bionj_tree); 
5247              }
5248            else
5249              {
5250                tree = (t_tree *)Make_Tree_From_Scratch(mixt_tree->n_otu,mixt_tree->data);
5251                Copy_Tree(bionj_tree,tree); 
5252              }
5253            break;
5254          }
5255        }
5256     
5257
5258      Copy_Tree(tree,mixt_tree);
5259      Free_Tree(tree);
5260      if(bionj_tree) Free_Tree(bionj_tree);
5261     
5262      if(io->mod->s_opt->random_input_tree) 
5263        {
5264          PhyML_Printf("\n\n. [%3d/%3d]",num_rand_tree+1,io->mod->s_opt->n_rand_starts);
5265          Random_Tree(mixt_tree);
5266        }
5267           
5268      tree = mixt_tree;
5269      do
5270        {
5271          if(tree != mixt_tree) Copy_Tree(mixt_tree,tree);
5272          Connect_CSeqs_To_Nodes(tree->data,tree);
5273          tree = tree->next;
5274        }
5275      while(tree);
5276     
5277     
5278      /*! Initialize t_beg in each mixture tree */
5279      tree = mixt_tree;
5280      do
5281        {
5282          time(&(tree->t_beg));
5283          tree = tree->next_mixt;
5284        }
5285      while(tree);
5286
5287      Prepare_Tree_For_Lk(mixt_tree);
5288
5289      MIXT_Chain_All(mixt_tree);
5290     
5291      /*! Check that all the edges in a mixt_tree at pointing
5292        to a single set of lengths
5293      */
5294      tree = mixt_tree;
5295      do
5296        {
5297          MIXT_Check_Single_Edge_Lens(tree);
5298          tree = tree->next_mixt;
5299        }
5300      while(tree);
5301                 
5302
5303      /*! Turn all branches to ON state */
5304      tree = mixt_tree;
5305      do
5306        {
5307          MIXT_Turn_Branches_OnOff(ON,tree);
5308          tree = tree->next_mixt;
5309        }
5310      while(tree);
5311
5312      /* TO DO
5313         
5314         1) REMOVE ROOT
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
5321         data types
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.
5342         X 20) Check alrt.c
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
5347      */
5348     
5349     
5350      MIXT_Check_Invar_Struct_In_Each_Partition_Elem(mixt_tree);
5351      MIXT_Check_RAS_Struct_In_Each_Partition_Elem(mixt_tree);
5352     
5353      Set_Both_Sides(YES,mixt_tree);     
5354
5355      if(mixt_tree->mod->s_opt->opt_topo)
5356        {
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);
5360        }
5361      else
5362        {
5363          if(mixt_tree->mod->s_opt->opt_subst_param || 
5364             mixt_tree->mod->s_opt->opt_bl)                       
5365            {
5366             
5367              Round_Optimize(mixt_tree,mixt_tree->data,ROUND_MAX);
5368            }
5369          else                                                   
5370            {
5371              Lk(NULL,mixt_tree);
5372            }
5373        }
5374
5375     
5376      PhyML_Printf("\n\n. Log-likelihood = %f",mixt_tree->c_lnL);
5377
5378      if((num_rand_tree == io->mod->s_opt->n_rand_starts-1) && (io->mod->s_opt->random_input_tree))
5379        {
5380          num_rand_tree--;
5381          io->mod->s_opt->random_input_tree = NO;
5382        }
5383     
5384      Br_Len_Involving_Invar(mixt_tree);
5385      Rescale_Br_Len_Multiplier_Tree(mixt_tree);
5386     
5387      /*! Print the tree estimated using the current random (or BioNJ) starting tree */
5388      if(io->mod->s_opt->n_rand_starts > 1)
5389        {
5390          Print_Tree(io->fp_out_trees,mixt_tree);
5391          fflush(NULL);
5392        }
5393     
5394      if(mixt_tree->c_lnL > best_lnL)
5395        {
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;
5401        }
5402
5403      /*! Initialize t_end in each mixture tree */
5404      tree = mixt_tree;
5405      do
5406        {
5407          time(&(tree->t_current));
5408          tree = tree->next_mixt;
5409        }
5410      while(tree);
5411     
5412      Print_Data_Structure(YES,mixt_tree->io->fp_out_stats,mixt_tree);
5413     
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);
5418    }
5419
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);
5423  Free(component);
5424
5425  /*! Bootstrap analysis */
5426  MIXT_Bootstrap(most_likely_tree,root);
5427
5428  while(io->prev != NULL) io = io->prev;
5429
5430  Free(most_likely_tree);
5431
5432  tree = mixt_tree;
5433  do
5434    {
5435      Free_Cseq(tree->data);
5436      tree = tree->next_mixt;
5437    }
5438  while(tree);
5439
5440  tree = mixt_tree;
5441  do
5442    {
5443      Free_Optimiz(tree->mod->s_opt);
5444      tree = tree->next;
5445    }
5446  while(tree);
5447
5448  Free_Model_Complete(mixt_tree->mod);
5449  Free_Model_Basic(mixt_tree->mod);
5450  Free_Tree(mixt_tree);
5451
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);
5455  Free_Input(io);
5456  XML_Free_XML_Tree(root);
5457
5458  Free(class_num);
5459
5460  fclose(fp);
5461}
5462
5463//////////////////////////////////////////////////////////////
5464//////////////////////////////////////////////////////////////
5465
5466/*! Check that the same nodes in the different mixt_trees are
5467  connected to the same taxa
5468*/
5469void Check_Taxa_Sets(t_tree *mixt_tree)
5470{
5471  t_tree *tree;
5472  int i;
5473
5474  tree = mixt_tree;
5475  do
5476    {
5477      if(tree->next)
5478        {         
5479          For(i,tree->n_otu)
5480            {
5481              if(strcmp(tree->a_nodes[i]->name,tree->next->a_nodes[i]->name))
5482                {
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);             
5487                  Exit("\n");
5488                }
5489            }
5490        }
5491      tree = tree->next;
5492    }
5493  while(tree);
5494
5495}
5496
5497//////////////////////////////////////////////////////////////
5498//////////////////////////////////////////////////////////////
5499
5500void Make_Ratematrice_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5501{
5502  char *model = NULL;
5503  int select;
5504
5505  model = XML_Get_Attribute_Value(instance,"model"); 
5506 
5507  if(model == NULL)
5508    {
5509      PhyML_Printf("\n== Poorly formated XML file.");
5510      PhyML_Printf("\n== Attribute 'model' is mandatory in a <ratematrix> node.");
5511      Exit("\n");
5512    }
5513 
5514  select = XML_Validate_Attr_Int(model,26,
5515                                 "xxxxx",    //0
5516                                 "JC69",     //1
5517                                 "K80",      //2
5518                                 "F81",      //3
5519                                 "HKY85",    //4
5520                                 "F84",      //5
5521                                 "TN93",     //6
5522                                 "GTR",      //7
5523                                 "CUSTOM",   //8
5524                                 "xxxxx",    //9
5525                                 "xxxxx",    //10
5526                                 "WAG",      //11
5527                                 "DAYHOFF",  //12
5528                                 "JTT",      //13
5529                                 "BLOSUM62", //14
5530                                 "MTREV",    //15
5531                                 "RTREV",    //16
5532                                 "CPREV",    //17
5533                                 "DCMUT",    //18
5534                                 "VT",       //19
5535                                 "MTMAM",    //20
5536                                 "MTART",    //21
5537                                 "HIVW",     //22
5538                                 "HIVB",     //23
5539                                 "CUSTOMAA", //24
5540                                 "LG");      //25
5541 
5542  if(select < 9)
5543    {
5544      mod->ns = 4;
5545      if(io->datatype != NT)
5546        {
5547          PhyML_Printf("\n== Data type and selected model are incompatible");
5548          Exit("\n");
5549        }
5550    }
5551  else
5552    {
5553      mod->ns = 20;
5554      if(io->datatype != AA)
5555        {
5556          PhyML_Printf("\n== Data type and selected model are incompatible");
5557          Exit("\n");
5558        }
5559    }
5560 
5561  io->mod->ns = mod->ns;
5562 
5563  mod->r_mat = (t_rmat *)Make_Rmat(mod->ns);
5564  Init_Rmat(mod->r_mat);
5565                 
5566  /*! Set model number & name */
5567  mod->whichmodel = Set_Whichmodel(select);
5568  Set_Model_Name(mod);
5569 
5570  if(mod->whichmodel == K80   || 
5571     mod->whichmodel == HKY85 || 
5572     mod->whichmodel == TN93)
5573    {
5574      char *tstv,*opt_tstv;
5575     
5576      tstv = XML_Get_Attribute_Value(instance,"tstv");
5577     
5578      if(tstv)
5579        {
5580          mod->s_opt->opt_kappa = NO;
5581          mod->kappa->v = String_To_Dbl(tstv);
5582        }
5583      else
5584        {
5585          mod->s_opt->opt_kappa = YES;
5586        }
5587     
5588      opt_tstv = XML_Get_Attribute_Value(instance,"optimise.tstv");
5589     
5590      if(opt_tstv)
5591        {
5592          if(!strcmp(opt_tstv,"true") || !strcmp(opt_tstv,"yes"))
5593            {
5594              mod->s_opt->opt_kappa = YES;
5595            }
5596          else
5597            {
5598              mod->s_opt->opt_kappa = NO;
5599            }
5600        }
5601    }
5602  else
5603    {
5604      mod->s_opt->opt_kappa = NO;
5605    }
5606 
5607  if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM)
5608    {
5609      char *opt_rr;
5610     
5611      opt_rr = XML_Get_Attribute_Value(instance,"optimise.rr");
5612     
5613      if(opt_rr)
5614        {
5615          if(!strcmp(opt_rr,"yes") || !strcmp(opt_rr,"true"))
5616            {
5617              mod->s_opt->opt_rr = YES;
5618            }
5619        }
5620    }
5621 
5622  /*! Custom model for nucleotide sequences. Read the corresponding
5623    code. */
5624  if(mod->whichmodel == CUSTOM)
5625    {
5626      char *model_code = NULL;
5627     
5628      model_code = XML_Get_Attribute_Value(instance,"model.code"); 
5629     
5630      if(!model_code)
5631        {
5632          PhyML_Printf("\n== No valid 'model.code' attribute could be found.\n");
5633          PhyML_Printf("\n== Please fix your XML file.\n");
5634          Exit("\n");                                   
5635        }
5636      else
5637        {
5638          strcpy(mod->custom_mod_string->s,model_code);
5639        }
5640    }
5641 
5642 
5643  /*! Custom model for amino-acids. Read in the rate matrix file */
5644  if(mod->whichmodel == CUSTOMAA)
5645    {
5646      char *r_mat_file;
5647     
5648      r_mat_file = XML_Get_Attribute_Value(instance,"ratematrix.file"); 
5649     
5650      if(!r_mat_file)
5651        {
5652          PhyML_Printf("\n== No valid 'ratematrix.file' attribute could be found.");
5653          PhyML_Printf("\n== Please fix your XML file.\n");
5654          Exit("\n");
5655        }
5656      else
5657        {
5658          mod->fp_aa_rate_mat = Openfile(r_mat_file,0);
5659          strcpy(mod->aa_rate_mat_file->s,r_mat_file);
5660        }
5661     
5662      Free(r_mat_file);
5663    }
5664
5665  char *buff;
5666 
5667  buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
5668  if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
5669    {
5670      mod->s_opt->opt_rmat_weight = YES;
5671    }
5672  else
5673    {
5674      mod->s_opt->opt_rmat_weight = NO;
5675    }
5676}
5677
5678//////////////////////////////////////////////////////////////
5679//////////////////////////////////////////////////////////////
5680
5681void Make_Efrq_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5682{
5683  char *buff;
5684 
5685  mod->e_frq = (t_efrq *)Make_Efrq(mod->ns);
5686  Init_Efrq(mod->e_frq);
5687 
5688  buff = XML_Get_Attribute_Value(instance,"optimise.freqs");
5689 
5690  if(buff)
5691    {
5692      if(!strcmp(buff,"yes") || !strcmp(buff,"true"))
5693        {
5694          if(io->datatype == AA)
5695            {
5696              PhyML_Printf("\n== Option 'optimise.freqs' set to 'yes' (or 'true')");
5697              PhyML_Printf("\n== is not allowed with amino-acid data.");
5698              Exit("\n");
5699            }
5700          mod->s_opt->opt_state_freq = YES;
5701        }
5702      Free(buff);
5703    }
5704 
5705  buff = XML_Get_Attribute_Value(instance,"aa.freqs");
5706 
5707  if(buff)
5708    {
5709      if(!strcmp(buff,"empirical"))
5710        {
5711          if(io->datatype == AA)
5712            {
5713              mod->s_opt->opt_state_freq = YES;
5714            }
5715          else if(io->datatype == NT)
5716            {
5717              mod->s_opt->opt_state_freq = NO;
5718            }
5719        }
5720      Free(buff);
5721    }
5722 
5723 
5724  buff = XML_Get_Attribute_Value(instance,"base.freqs");
5725 
5726  if(buff)
5727    {
5728      if(io->datatype == AA)
5729        {
5730          PhyML_Printf("\n== Option 'base.freqs' is not allowed with amino-acid data.");
5731          Exit("\n");
5732        }
5733      else
5734        {
5735          phydbl A,C,G,T;
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;
5743        }
5744    } 
5745
5746 
5747  buff = XML_Get_Attribute_Value(instance->parent,"optimise.weights");
5748  if(buff && (!strcmp(buff,"yes") || !strcmp(buff,"true")))
5749    {
5750      mod->s_opt->opt_efrq_weight = YES;
5751    }
5752  else
5753    {
5754      mod->s_opt->opt_efrq_weight = NO;
5755    }
5756
5757}
5758
5759//////////////////////////////////////////////////////////////
5760//////////////////////////////////////////////////////////////
5761
5762void Make_Topology_From_XML_Node(xml_node *instance, option *io, t_mod *mod)
5763{
5764  // Starting tree
5765  char *init_tree;
5766 
5767  init_tree = XML_Get_Attribute_Value(instance,"init.tree");
5768 
5769  if(!init_tree)
5770    {
5771      PhyML_Printf("\n== Attribute 'init.tree=bionj|user|random' is mandatory");
5772      PhyML_Printf("\n== Please amend your XML file accordingly.");
5773      Exit("\n");
5774    }
5775 
5776  if(!strcmp(init_tree,"user") || 
5777     !strcmp(init_tree,"User"))
5778    {
5779      char *starting_tree = NULL;
5780      starting_tree = XML_Get_Attribute_Value(instance,"file.name");
5781     
5782      if(!Filexists(starting_tree))
5783        {
5784          PhyML_Printf("\n== The tree file '%s' could not be found.",starting_tree);
5785          Exit("\n");
5786        }
5787      else
5788        {
5789          strcpy(io->in_tree_file,starting_tree);
5790          io->in_tree = 2;
5791          io->fp_in_tree = Openfile(io->in_tree_file,0);
5792        }
5793    }
5794  else if(!strcmp(init_tree,"random") || 
5795          !strcmp(init_tree,"Random"))
5796    {
5797      char *n_rand_starts = NULL;
5798     
5799      io->mod->s_opt->random_input_tree = YES;
5800     
5801      n_rand_starts = XML_Get_Attribute_Value(instance,"n.rand.starts");
5802     
5803      if(n_rand_starts)
5804        {
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");
5807        }
5808     
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);
5813    }
5814  else if(!strcmp(init_tree,"parsimony") || 
5815          !strcmp(init_tree,"Parsimony"))
5816    {
5817      io->in_tree = 1;
5818    }
5819 
5820  // Estimate tree topology
5821  char *optimise = NULL;
5822 
5823  optimise = XML_Get_Attribute_Value(instance,"optimise.tree");
5824                           
5825  if(optimise)
5826    {
5827      int select;
5828     
5829      select = XML_Validate_Attr_Int(optimise,6,
5830                                     "true","yes","y",
5831                                     "false","no","n");
5832     
5833      if(select > 2) io->mod->s_opt->opt_topo = NO;
5834      else
5835        {
5836          char *search;
5837          int select;
5838         
5839          search = XML_Get_Attribute_Value(instance,"search");                             
5840          select = XML_Validate_Attr_Int(search,4,"spr","nni","best","none");
5841         
5842          switch(select)
5843            {
5844            case 0:
5845              {
5846                io->mod->s_opt->topo_search = SPR_MOVE;
5847                io->mod->s_opt->opt_topo    = YES;
5848                break;
5849              }
5850            case 1:
5851              {
5852                io->mod->s_opt->topo_search = NNI_MOVE;
5853                io->mod->s_opt->opt_topo    = YES;
5854                break;
5855              }
5856            case 2:
5857              {
5858                io->mod->s_opt->topo_search = BEST_OF_NNI_AND_SPR;
5859                io->mod->s_opt->opt_topo    = YES;
5860                break;
5861              }
5862            case 3:
5863              {
5864                io->mod->s_opt->opt_topo    = NO;
5865                break;
5866              }
5867            default:
5868              {
5869                PhyML_Printf("\n== Topology search option '%s' is not valid.",search);
5870                Exit("\n");
5871                break;
5872              }
5873            }
5874        }
5875    } 
5876}
5877
5878
5879//////////////////////////////////////////////////////////////
5880//////////////////////////////////////////////////////////////
5881
5882void Make_RAS_From_XML_Node(xml_node *parent, t_mod *mod)
5883{
5884  xml_node *w;
5885  char *family;
5886  int select;
5887
5888  mod->ras->n_catg = 0;
5889 
5890  XML_Check_Siterates_Node(parent);
5891                           
5892  w = XML_Search_Node_Name("weights",YES,parent);
5893  if(w)
5894    {
5895      family = XML_Get_Attribute_Value(w,"family");
5896      select = XML_Validate_Attr_Int(family,3,"gamma","gamma+inv","freerates");
5897      switch(select)
5898        {
5899        case 0: // Gamma model
5900          {
5901            char *alpha,*alpha_opt;
5902           
5903            mod->s_opt->opt_pinvar = NO;
5904            mod->ras->invar        = NO;
5905           
5906            alpha = XML_Get_Attribute_Value(w,"alpha");
5907           
5908            if(alpha)
5909              {
5910                if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") || 
5911                   !strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5912                  {
5913                    mod->s_opt->opt_alpha = YES;
5914                  }
5915                else
5916                  {                                       
5917                    mod->s_opt->opt_alpha = NO;
5918                    mod->ras->alpha->v = String_To_Dbl(alpha);
5919                  }
5920              }
5921           
5922            alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5923           
5924            if(alpha_opt)
5925              {
5926                if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5927                  {
5928                    mod->s_opt->opt_alpha = YES;
5929                  }
5930                else
5931                  {                                       
5932                    mod->s_opt->opt_alpha = NO;
5933                  }
5934              }
5935           
5936            mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
5937           
5938            Make_RAS_Complete(mod->ras);
5939           
5940            break;
5941          }
5942        case 1: // Gamma+Inv model
5943          {
5944            char *alpha,*pinv,*alpha_opt,*pinv_opt;
5945           
5946            mod->ras->invar        = YES;
5947            mod->s_opt->opt_pinvar = YES;
5948           
5949            alpha = XML_Get_Attribute_Value(w,"alpha");
5950           
5951            if(alpha)
5952              {
5953                if(!strcmp(alpha,"estimate") || !strcmp(alpha,"estimated") || 
5954                   !strcmp(alpha,"optimise") || !strcmp(alpha,"optimised"))
5955                  {
5956                    mod->s_opt->opt_alpha = YES;
5957                  }
5958                else
5959                  {
5960                    mod->s_opt->opt_alpha = NO;
5961                    mod->ras->alpha->v = String_To_Dbl(alpha);;
5962                  }
5963              }
5964           
5965            alpha_opt = XML_Get_Attribute_Value(w,"optimise.alpha");
5966           
5967            if(alpha_opt)
5968              {
5969                if(!strcmp(alpha_opt,"yes") || !strcmp(alpha_opt,"true"))
5970                  {
5971                    mod->s_opt->opt_alpha = YES;
5972                  }
5973                else
5974                  {                                       
5975                    mod->s_opt->opt_alpha = NO;
5976                  }
5977              }
5978           
5979           
5980            pinv = XML_Get_Attribute_Value(w,"pinv");
5981           
5982            if(pinv)
5983              {
5984                if(!strcmp(pinv,"estimate") || !strcmp(pinv,"estimated") || 
5985                   !strcmp(pinv,"optimise") || !strcmp(pinv,"optimised"))
5986                  {
5987                    mod->s_opt->opt_pinvar = YES;
5988                  }
5989                else
5990                  {
5991                    mod->s_opt->opt_pinvar = NO;
5992                    mod->ras->pinvar->v = String_To_Dbl(pinv);;
5993                  }
5994              }
5995           
5996            pinv_opt = XML_Get_Attribute_Value(w,"optimise.pinv");
5997           
5998            if(pinv_opt)
5999              {
6000                if(!strcmp(pinv_opt,"yes") || !strcmp(pinv_opt,"true"))
6001                  {
6002                    mod->s_opt->opt_pinvar = YES;
6003                  }
6004                else
6005                  {                                       
6006                    mod->s_opt->opt_pinvar = NO;
6007                  }
6008              }
6009           
6010            mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
6011            break;
6012          }
6013        case 2: // FreeRate model
6014          {
6015            char *opt_freerates;
6016
6017            mod->ras->free_mixt_rates = YES;
6018           
6019            mod->s_opt->opt_free_mixt_rates = YES;
6020
6021            opt_freerates = XML_Get_Attribute_Value(w,"optimise.freerates");
6022                       
6023            if(opt_freerates)
6024              {
6025                if(!strcmp(opt_freerates,"yes") || !strcmp(opt_freerates,"true"))
6026                  {
6027                    mod->s_opt->opt_free_mixt_rates = YES;
6028                  }
6029                else 
6030                  {                                       
6031                    mod->s_opt->opt_free_mixt_rates = NO;
6032                  }
6033              }
6034
6035            mod->ras->n_catg = XML_Get_Number_Of_Classes_Siterates(parent);
6036            break;
6037          }
6038        default:
6039          {
6040            PhyML_Printf("\n== family: %s",family);
6041            PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
6042            Exit("\n");
6043          }
6044        }
6045    }
6046 
6047  int nc = XML_Get_Number_Of_Classes_Siterates(parent);
6048 
6049  if(w && nc != mod->ras->n_catg)
6050    {
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.");
6055      Exit("\n");
6056    }
6057 
6058  if(!w) mod->ras->n_catg = nc;
6059 
6060  Make_RAS_Complete(mod->ras);
6061 
6062}
6063
6064//////////////////////////////////////////////////////////////
6065//////////////////////////////////////////////////////////////
6066//////////////////////////////////////////////////////////////
6067//////////////////////////////////////////////////////////////
6068//////////////////////////////////////////////////////////////
6069//////////////////////////////////////////////////////////////
6070
Note: See TracBrowser for help on using the repository browser.