source: trunk/GDE/RAxML/treeIO.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.4 KB
Line 
1/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright August 2006 by Alexandros Stamatakis
3 *
4 *  Partially derived from
5 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6 * 
7 *  and
8 *
9 *  Programs of the PHYLIP package by Joe Felsenstein.
10 *
11 *  This program is free software; you may redistribute it and/or modify its
12 *  under the terms of the GNU General Public License as published by the Free
13 *  Software Foundation; either version 2 of the License, or (at your option)
14 *  any later version.
15 *
16 *  This program is distributed in the hope that it will be useful, but
17 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19 *  for more details.
20 *
21 *
22 *  For any other enquiries send an Email to Alexandros Stamatakis
23 *  Alexandros.Stamatakis@epfl.ch
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 *
27 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
28 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
29 */
30
31#ifndef WIN32
32#include <sys/times.h>
33#include <sys/types.h>
34#include <sys/time.h>
35#include <unistd.h>
36#endif
37
38#include <math.h>
39#include <time.h>
40#include <stdlib.h>
41#include <stdio.h>
42#include <ctype.h>
43#include <string.h>
44
45#include "axml.h"
46
47extern FILE *INFILE;
48extern char infoFileName[1024];
49extern char tree_file[1024];
50extern char *likelihood_key;
51extern char *ntaxa_key;
52extern char *smoothed_key;
53extern int partCount;
54extern double masterTime;
55
56
57
58
59
60stringHashtable *initStringHashTable(hashNumberType n)
61{
62  /*
63     init with primes
64  */
65   
66  static const hashNumberType initTable[] = {53, 97, 193, 389, 769, 1543, 3079, 6151, 12289, 24593, 49157, 98317,
67                                             196613, 393241, 786433, 1572869, 3145739, 6291469, 12582917, 25165843,
68                                             50331653, 100663319, 201326611, 402653189, 805306457, 1610612741};
69 
70
71  /* init with powers of two
72
73  static const  hashNumberType initTable[] = {64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384,
74                                              32768, 65536, 131072, 262144, 524288, 1048576, 2097152,
75                                              4194304, 8388608, 16777216, 33554432, 67108864, 134217728,
76                                              268435456, 536870912, 1073741824, 2147483648U};
77  */
78 
79  stringHashtable *h = (stringHashtable*)rax_malloc(sizeof(stringHashtable));
80 
81  hashNumberType
82    tableSize,
83    i,
84    primeTableLength = sizeof(initTable)/sizeof(initTable[0]),
85    maxSize = (hashNumberType)-1;   
86
87  assert(n <= maxSize);
88
89  i = 0;
90
91  while(initTable[i] < n && i < primeTableLength)
92    i++;
93
94  assert(i < primeTableLength);
95
96  tableSize = initTable[i]; 
97
98  h->table = (stringEntry**)rax_calloc(tableSize, sizeof(stringEntry*));
99  h->tableSize = tableSize;   
100
101  return h;
102}
103
104
105static hashNumberType  hashString(char *p, hashNumberType tableSize)
106{
107  hashNumberType h = 0;
108 
109  for(; *p; p++)
110    h = 31 * h + *p;
111 
112  return (h % tableSize);
113}
114
115 
116
117void addword(char *s, stringHashtable *h, int nodeNumber)
118{
119  hashNumberType position = hashString(s, h->tableSize);
120  stringEntry *p = h->table[position];
121 
122  for(; p!= NULL; p = p->next)
123    {
124      if(strcmp(s, p->word) == 0)               
125        return;         
126    }
127
128  p = (stringEntry *)rax_malloc(sizeof(stringEntry));
129
130  assert(p);
131 
132  p->nodeNumber = nodeNumber;
133  p->word = (char *)rax_malloc((strlen(s) + 1) * sizeof(char));
134
135  strcpy(p->word, s);
136 
137  p->next =  h->table[position];
138 
139  h->table[position] = p;
140}
141
142int lookupWord(char *s, stringHashtable *h)
143{
144  hashNumberType position = hashString(s, h->tableSize);
145  stringEntry *p = h->table[position];
146 
147  for(; p!= NULL; p = p->next)
148    {
149      if(strcmp(s, p->word) == 0)               
150        return p->nodeNumber;           
151    }
152
153  return -1;
154}
155
156
157int countTips(nodeptr p, int numsp)
158{
159  if(isTip(p->number, numsp)) 
160    return 1;   
161  {
162    nodeptr q;
163    int tips = 0;
164
165    q = p->next;
166    while(q != p)
167      { 
168        tips += countTips(q->back, numsp);
169        q = q->next;
170      } 
171   
172    return tips;
173  }
174}
175
176
177static double getBranchLength(tree *tr, int perGene, nodeptr p)
178{
179  double 
180    z = 0.0,
181    x = 0.0;
182
183  assert(perGene != NO_BRANCHES);
184             
185  if(!tr->multiBranch)
186    {
187      assert(tr->fracchange != -1.0);
188      z = p->z[0];
189      if (z < zmin) 
190        z = zmin;       
191     
192      x = -log(z) * tr->fracchange;           
193    }
194  else
195    {
196      if(perGene == SUMMARIZE_LH)
197        {
198          int 
199            i;
200         
201          double 
202            avgX = 0.0;
203                     
204          for(i = 0; i < tr->numBranches; i++)
205            {
206              assert(tr->partitionContributions[i] != -1.0);
207              assert(tr->fracchanges[i] != -1.0);
208              z = p->z[i];
209              if(z < zmin) 
210                z = zmin;       
211              x = -log(z) * tr->fracchanges[i];
212              avgX += x * tr->partitionContributions[i];
213            }
214
215          x = avgX;
216        }
217      else
218        {       
219          assert(tr->fracchanges[perGene] != -1.0);
220          assert(perGene >= 0 && perGene < tr->numBranches);
221         
222          z = p->z[perGene];
223         
224          if(z < zmin) 
225            z = zmin;           
226         
227          x = -log(z) * tr->fracchanges[perGene];         
228        }
229    }
230
231  return x;
232}
233
234
235 
236
237
238static char *Tree2StringREC(char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, 
239                            boolean printLikelihood, boolean rellTree, boolean finalPrint, int perGene, boolean branchLabelSupport, boolean printSHSupport, boolean printIC, boolean printSHSupports)
240{
241  char  *nameptr;           
242     
243  if(isTip(p->number, tr->rdta->numsp)) 
244    {             
245      if(printNames)
246        {
247          nameptr = tr->nameList[p->number];     
248          sprintf(treestr, "%s", nameptr);
249        }
250      else
251        sprintf(treestr, "%d", p->number);   
252       
253      while (*treestr) treestr++;
254    }
255  else 
256    {                   
257      *treestr++ = '(';
258      treestr = Tree2StringREC(treestr, tr, p->next->back, printBranchLengths, printNames, printLikelihood, rellTree, 
259                               finalPrint, perGene, branchLabelSupport, printSHSupport, printIC, printSHSupports);
260      *treestr++ = ',';
261      treestr = Tree2StringREC(treestr, tr, p->next->next->back, printBranchLengths, printNames, printLikelihood, rellTree, 
262                               finalPrint, perGene, branchLabelSupport, printSHSupport, printIC, printSHSupports);
263      if(p == tr->start->back) 
264        {
265          *treestr++ = ',';
266          treestr = Tree2StringREC(treestr, tr, p->back, printBranchLengths, printNames, printLikelihood, rellTree, 
267                                   finalPrint, perGene, branchLabelSupport, printSHSupport, printIC, printSHSupports);
268        }
269      *treestr++ = ')';                   
270    }
271
272  if(p == tr->start->back) 
273    {           
274      if(printBranchLengths && !rellTree)
275        sprintf(treestr, ":0.0;\n");
276      else
277        sprintf(treestr, ";\n");                       
278    }
279  else 
280    {                   
281      if(rellTree || branchLabelSupport || printSHSupport || printIC || printSHSupports)
282        {               
283          if(( !isTip(p->number, tr->rdta->numsp)) && 
284             ( !isTip(p->back->number, tr->rdta->numsp)))
285            {                         
286              assert(p->bInf != (branchInfo *)NULL);               
287             
288              assert(rellTree + branchLabelSupport + printSHSupport + printSHSupports == 1);
289
290              if(rellTree)
291                {
292                  if(printIC)
293                    sprintf(treestr, "%1.2f:%8.20f", p->bInf->ic, p->z[0]);
294                  else
295                    sprintf(treestr, "%d:%8.20f", p->bInf->support, p->z[0]);
296                }
297             
298              if(branchLabelSupport)
299                {
300                  if(printIC)
301                    sprintf(treestr, ":%8.20f[%1.2f,%1.2f]", p->z[0], p->bInf->ic, p->bInf->icAll);
302                  else             
303                    sprintf(treestr, ":%8.20f[%d]", p->z[0], p->bInf->support);
304                }
305             
306              if(printSHSupport)
307                sprintf(treestr, ":%8.20f[%d]", getBranchLength(tr, perGene, p), p->bInf->support);
308
309              if(printSHSupports)
310                {
311                  int 
312                    model;
313                 
314                  sprintf(treestr, ":%8.20f[", getBranchLength(tr, perGene, p));
315                  while(*treestr) 
316                    treestr++;
317                 
318                  for(model = 0; model < tr->NumberOfModels - 1; model++)                   
319                    {
320                      sprintf(treestr, "%d,", p->bInf->supports[model]);
321                       while(*treestr) 
322                         treestr++;
323                    }
324
325                  sprintf(treestr, "%d]", p->bInf->supports[model]);
326                }
327             
328            }
329          else         
330            {
331              if(rellTree || branchLabelSupport)
332                sprintf(treestr, ":%8.20f", p->z[0]);   
333              if(printSHSupport || printSHSupports)
334                sprintf(treestr, ":%8.20f", getBranchLength(tr, perGene, p));
335            }
336        }
337      else
338        {
339          if(printBranchLengths)           
340            sprintf(treestr, ":%8.20f", getBranchLength(tr, perGene, p));                 
341          else     
342            sprintf(treestr, "%s", "\0");           
343        }     
344    }
345 
346  while (*treestr) treestr++;
347  return  treestr;
348}
349
350
351static void collectSubtrees(tree *tr, nodeptr *subtrees, int *count, int ogn)
352{
353  int i;
354  for(i = tr->mxtips + 1; i <= tr->mxtips + tr->mxtips - 2; i++)
355    {
356      nodeptr p, q;
357      p = tr->nodep[i];
358      if(countTips(p, tr->rdta->numsp) == ogn)
359        {
360          subtrees[*count] = p;
361          *count = *count + 1;
362        }
363      q = p->next;
364      while(q != p)
365        {
366          if(countTips(q, tr->rdta->numsp) == ogn)
367            {
368              subtrees[*count] = q;
369              *count = *count + 1;
370            }
371          q = q->next;
372        }
373    }
374}
375
376static void checkOM(nodeptr p, int *n, int *c, tree *tr)
377{
378  if(isTip(p->number, tr->rdta->numsp))
379    {
380      n[*c] = p->number;
381      *c = *c + 1;     
382    }
383  else
384    {
385      nodeptr q = p->next;
386
387      while(q != p)
388        {
389          checkOM(q->back, n, c, tr);
390          q = q->next;
391        }
392    }
393}
394   
395static char *rootedTreeREC(char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, 
396                           boolean printLikelihood, boolean rellTree, 
397                           boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport)
398{
399  char  *nameptr;
400
401  if(isTip(p->number, tr->rdta->numsp)) 
402    {       
403      if(printNames)
404        {
405          nameptr = tr->nameList[p->number];     
406          sprintf(treestr, "%s", nameptr);
407        }
408      else
409        sprintf(treestr, "%d", p->number);
410     
411      while (*treestr) treestr++;
412    }
413  else 
414    {
415      *treestr++ = '(';
416      treestr = rootedTreeREC(treestr, tr, p->next->back, printBranchLengths, printNames, printLikelihood, 
417                              rellTree, finalPrint, adef, perGene, branchLabelSupport, printSHSupport);
418      *treestr++ = ',';
419      treestr = rootedTreeREC(treestr, tr, p->next->next->back, printBranchLengths, printNames, printLikelihood, 
420                              rellTree, finalPrint, adef, perGene, branchLabelSupport, printSHSupport);     
421      *treestr++ = ')';           
422    }
423
424  if(rellTree || branchLabelSupport || printSHSupport)
425    {           
426      if(( !isTip(p->number, tr->rdta->numsp)) && 
427         ( !isTip(p->back->number, tr->rdta->numsp)))
428        {                             
429          assert(p->bInf != (branchInfo *)NULL);
430             
431          if(rellTree)
432            sprintf(treestr, "%d:%8.20f", p->bInf->support, p->z[0]);     
433          if(branchLabelSupport)
434            sprintf(treestr, ":%8.20f[%d]", p->z[0], p->bInf->support);
435          if(printSHSupport)
436            sprintf(treestr, ":%8.20f[%d]", getBranchLength(tr, perGene, p), p->bInf->support);       
437        }
438      else             
439        {
440          if(rellTree || branchLabelSupport)
441            sprintf(treestr, ":%8.20f", p->z[0]);       
442          if(printSHSupport)
443            sprintf(treestr, ":%8.20f", getBranchLength(tr, perGene, p));
444        }
445    }
446  else
447    {
448      if(printBranchLengths)       
449        sprintf(treestr, ":%8.20f", getBranchLength(tr, perGene, p));             
450      else         
451        sprintf(treestr, "%s", "\0");       
452    }     
453
454  while (*treestr) treestr++;
455  return  treestr;
456}
457
458static char *rootedTree(char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, 
459                        boolean printLikelihood, boolean rellTree, 
460                        boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport)
461{
462  double oldz[NUM_BRANCHES];
463
464  for(int i = 0; i < tr->numBranches; i++)
465    oldz[i] = p->z[i];
466
467  if(rellTree)   
468    {   
469      p->z[0] = p->back->z[0] = oldz[0] * 0.5;
470      /*printf("%f\n",  p->z[0]);*/
471    }
472  else
473    {
474      if(printBranchLengths)
475        {
476          double rz, z;
477          assert(perGene != NO_BRANCHES);
478
479          if(!tr->multiBranch)
480            {
481              assert(tr->fracchange != -1.0);
482              z = -log(p->z[0]) * tr->fracchange;
483              rz = exp(-(z * 0.5)/ tr->fracchange);
484              p->z[0] = p->back->z[0] = rz;
485            }
486          else
487            {
488              if(perGene == SUMMARIZE_LH)
489                {                                               
490                  for(int i = 0; i < tr->numBranches; i++)
491                    {   
492                      assert(tr->fracchanges[i] != -1.0);
493                      z    = -log(p->z[i]) * tr->fracchanges[i];                     
494                      rz   = exp(-(z * 0.5)/ tr->fracchanges[i]);
495                      p->z[i] = p->back->z[i] = rz;                 
496                    }           
497                }                   
498              else
499                {               
500                  assert(tr->fracchanges[perGene] != -1.0);
501                  assert(perGene >= 0 && perGene < tr->numBranches);
502                  z = -log(p->z[perGene]) * tr->fracchanges[perGene];
503                  rz = exp(-(z * 0.5)/ tr->fracchanges[perGene]);
504                  p->z[perGene] = p->back->z[perGene] = rz;                           
505                }
506            }
507        }
508    }
509
510  *treestr = '(';
511  treestr++;
512  treestr = rootedTreeREC(treestr, tr, p,  printBranchLengths, printNames, printLikelihood, rellTree, finalPrint, 
513                          adef, perGene, branchLabelSupport, printSHSupport);
514  *treestr = ',';
515  treestr++;
516  treestr = rootedTreeREC(treestr, tr, p->back,  printBranchLengths, printNames, printLikelihood, rellTree, finalPrint, 
517                          adef, perGene, branchLabelSupport, printSHSupport);
518  sprintf(treestr, ");\n");
519  while (*treestr) treestr++;
520
521
522  for(int i = 0; i < tr->numBranches; i++)
523    p->z[i] = p->back->z[i] = oldz[i]; 
524   
525  return  treestr;
526}
527
528
529
530char *Tree2String(char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, boolean printLikelihood, 
531                  boolean rellTree, 
532                  boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport, boolean printIC, boolean printSHSupports)
533{ 
534
535  if(rellTree)
536    assert(!branchLabelSupport && !printSHSupport);
537
538  if(branchLabelSupport)
539    assert(!rellTree && !printSHSupport);
540
541  if(printSHSupport)
542    assert(!branchLabelSupport && !rellTree);
543
544  if(finalPrint && adef->outgroup)
545    {
546      nodeptr startNode = tr->start;
547
548      if(tr->numberOfOutgroups > 1)
549        {
550          nodeptr root;
551          nodeptr *subtrees = (nodeptr *)rax_malloc(sizeof(nodeptr) * tr->mxtips);
552          int i, k, count = 0;
553          int *nodeNumbers = (int*)rax_malloc(sizeof(int) * tr->numberOfOutgroups);
554          int *foundVector = (int*)rax_malloc(sizeof(int) * tr->numberOfOutgroups);
555          boolean monophyletic = FALSE;
556
557          collectSubtrees(tr, subtrees, &count, tr->numberOfOutgroups);
558
559          /*printf("Found %d subtrees of size  %d\n", count, tr->numberOfOutgroups);*/
560         
561          for(i = 0; (i < count) && (!monophyletic); i++)
562            {
563              int l, sum, nc = 0;
564              for(k = 0; k <  tr->numberOfOutgroups; k++)
565                {
566                  nodeNumbers[k] = -1;
567                  foundVector[k] = 0;
568                }
569
570              checkOM(subtrees[i], nodeNumbers, &nc, tr);             
571             
572              for(l = 0; l < tr->numberOfOutgroups; l++)
573                for(k = 0; k < tr->numberOfOutgroups; k++)
574                  {
575                    if(nodeNumbers[l] == tr->outgroupNums[k])
576                      foundVector[l] = 1;
577                  }
578             
579              sum = 0;
580              for(l = 0; l < tr->numberOfOutgroups; l++)
581                sum += foundVector[l];
582             
583              if(sum == tr->numberOfOutgroups)
584                {                         
585                  root = subtrees[i];
586                  tr->start = root;             
587                  /*printf("outgroups are monphyletic!\n");*/
588                  monophyletic = TRUE;           
589                }
590              else
591                {
592                  if(sum > 0)
593                    {
594                      /*printf("outgroups are NOT monophyletic!\n");*/
595                      monophyletic = FALSE;
596                    }       
597                }       
598            }
599         
600          if(!monophyletic)
601            {
602              printf("WARNING, outgroups are not monophyletic, using first outgroup \"%s\"\n", tr->nameList[tr->outgroupNums[0]]);
603              printf("from the list to root the tree!\n");
604             
605
606              {
607                FILE *infoFile = myfopen(infoFileName, "ab");
608
609                fprintf(infoFile, "\nWARNING, outgroups are not monophyletic, using first outgroup \"%s\"\n", tr->nameList[tr->outgroupNums[0]]);
610                fprintf(infoFile, "from the list to root the tree!\n");
611               
612                fclose(infoFile);
613              }
614
615
616              tr->start = tr->nodep[tr->outgroupNums[0]];
617             
618              rootedTree(treestr, tr, tr->start->back, printBranchLengths, printNames, printLikelihood, rellTree, 
619                         finalPrint, adef, perGene, branchLabelSupport, printSHSupport);
620            }
621          else
622            {       
623              if(isTip(tr->start->number, tr->rdta->numsp))
624                {
625                  printf("Outgroup-Monophyly ERROR; tr->start is a tip \n");
626                  errorExit(-1);
627                }
628              if(isTip(tr->start->back->number, tr->rdta->numsp))
629                {
630                  printf("Outgroup-Monophyly ERROR; tr->start is a tip \n");
631                  errorExit(-1);
632                }
633                     
634              rootedTree(treestr, tr, tr->start->back, printBranchLengths, printNames, printLikelihood, rellTree, 
635                         finalPrint, adef, perGene, branchLabelSupport, printSHSupport);
636            }
637         
638          rax_free(foundVector);
639          rax_free(nodeNumbers);
640          rax_free(subtrees);
641        }
642      else
643        {         
644          tr->start = tr->nodep[tr->outgroupNums[0]];   
645          rootedTree(treestr, tr, tr->start->back, printBranchLengths, printNames, printLikelihood, rellTree, 
646                     finalPrint, adef, perGene, branchLabelSupport, printSHSupports);
647        }     
648
649      tr->start = startNode;
650    }
651  else   
652    Tree2StringREC(treestr, tr, p, printBranchLengths, printNames, printLikelihood, rellTree, 
653                   finalPrint, perGene, branchLabelSupport, printSHSupport, printIC, printSHSupports); 
654   
655 
656  while (*treestr) treestr++;
657 
658  return treestr;
659}
660
661
662void printTreePerGene(tree *tr, analdef *adef, const char *fileName, const char *permission)
663{ 
664  FILE *treeFile;
665  char extendedTreeFileName[1024];
666  char buf[16];
667  int i;
668
669  assert(adef->perGeneBranchLengths);
670     
671  for(i = 0; i < tr->numBranches; i++) 
672    {
673      strcpy(extendedTreeFileName, fileName);
674      sprintf(buf,"%d", i);
675      strcat(extendedTreeFileName, ".PARTITION.");
676      strcat(extendedTreeFileName, buf);
677      /*printf("Partitiuon %d file %s\n", i, extendedTreeFileName);*/
678      Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, TRUE, adef, i, FALSE, FALSE, FALSE, FALSE);
679      treeFile = myfopen(extendedTreeFileName, permission);
680      fprintf(treeFile, "%s", tr->tree_string);
681      fclose(treeFile);
682    } 
683   
684}
685
686
687
688/*=======================================================================*/
689/*                         Read a tree from a file                       */
690/*=======================================================================*/
691
692
693/*  1.0.A  Processing of quotation marks in comment removed
694 */
695
696static int treeFinishCom (FILE *fp, char **strp)
697{
698  int  ch;
699 
700  while ((ch = getc(fp)) != EOF && ch != ']') {
701    if (strp != NULL) *(*strp)++ = ch;    /* save character  */
702    if (ch == '[') {                      /* nested comment; find its end */
703      if ((ch = treeFinishCom(fp, strp)) == EOF)  break;
704      if (strp != NULL) *(*strp)++ = ch;  /* save closing ]  */
705    }
706  }
707 
708  if (strp != NULL) **strp = '\0';        /* terminate string  */
709  return  ch;
710} /* treeFinishCom */
711
712
713static int treeGetCh (FILE *fp)         /* get next nonblank, noncomment character */
714{ /* treeGetCh */
715  int  ch;
716
717  while ((ch = getc(fp)) != EOF) {
718    if (whitechar(ch)) ;
719    else if (ch == '[') {                   /* comment; find its end */
720      if ((ch = treeFinishCom(fp, (char **) NULL)) == EOF)  break;
721    }
722    else  break;
723  }
724 
725  return  ch;
726} /* treeGetCh */
727
728
729static boolean treeLabelEnd (int ch)
730{
731  switch (ch) 
732    {
733    case EOF: 
734    case '\0': 
735    case '\t': 
736    case '\n': 
737    case '\r': 
738    case ' ':
739    case ':': 
740    case ',':   
741    case '(':   
742    case ')': 
743    case ';':
744      return TRUE;
745    default:
746      break;
747    }
748  return FALSE;
749} 
750
751
752static boolean  treeGetLabel (FILE *fp, char *lblPtr, int maxlen)
753{
754  int      ch;
755  boolean  done, quoted, lblfound;
756
757  if (--maxlen < 0) 
758    lblPtr = (char *) NULL; 
759  else 
760    if (lblPtr == NULL) 
761      maxlen = 0;
762
763  ch = getc(fp);
764  done = treeLabelEnd(ch);
765
766  lblfound = ! done;
767  quoted = (ch == '\'');
768  if (quoted && ! done) 
769    {
770      ch = getc(fp); 
771      done = (ch == EOF);
772    }
773
774  while (! done) 
775    {
776      if (quoted) 
777        {
778          if (ch == '\'') 
779            {
780              ch = getc(fp); 
781              if (ch != '\'') 
782                break;
783            }
784        }
785      else 
786        if (treeLabelEnd(ch)) break;     
787
788      if (--maxlen >= 0) *lblPtr++ = ch;
789      ch = getc(fp);
790      if (ch == EOF) break;
791    }
792
793  if (ch != EOF)  (void) ungetc(ch, fp);
794
795  if (lblPtr != NULL) *lblPtr = '\0';
796
797  return lblfound;
798}
799
800
801static boolean  treeFlushLabel (FILE *fp)
802{ 
803  return  treeGetLabel(fp, (char *) NULL, (int) 0);
804} 
805
806
807
808
809static int treeFindTipByLabelString(char  *str, tree *tr, boolean check)                   
810{
811  int 
812    lookup = lookupWord(str, tr->nameHash);
813
814  if(lookup > 0)
815    {
816      if(check)
817        assert(! tr->nodep[lookup]->back);
818      return lookup;
819    }
820  else
821    { 
822      printf("ERROR: Cannot find tree species: %s\n", str);
823      return  0;
824    }
825}
826
827
828int treeFindTipName(FILE *fp, tree *tr, boolean check)
829{
830  char    str[nmlngth+2];
831  int      n;
832
833  if(treeGetLabel(fp, str, nmlngth+2))
834    n = treeFindTipByLabelString(str, tr, check);
835  else
836    n = 0;
837   
838
839  return  n;
840} 
841
842
843
844static void  treeEchoContext (FILE *fp1, FILE *fp2, int n)
845{ /* treeEchoContext */
846  int      ch;
847  boolean  waswhite;
848 
849  waswhite = TRUE;
850 
851  while (n > 0 && ((ch = getc(fp1)) != EOF)) {
852    if (whitechar(ch)) {
853      ch = waswhite ? '\0' : ' ';
854      waswhite = TRUE;
855    }
856    else {
857      waswhite = FALSE;
858    }
859   
860    if (ch > '\0') {putc(ch, fp2); n--;}
861  }
862} /* treeEchoContext */
863
864static boolean treeNeedCh (FILE *fp, int c1, const char *where);
865
866
867static boolean treeProcessLength (FILE *fp, double *dptr, int *branchLabel, boolean storeBranchLabels, tree *tr)
868{
869  int
870    ch;
871 
872  if((ch = treeGetCh(fp)) == EOF) 
873    return FALSE;    /*  Skip comments */
874  (void) ungetc(ch, fp);
875 
876  if(fscanf(fp, "%lf", dptr) != 1) 
877    {
878      printf("ERROR: treeProcessLength: Problem reading branch length\n");
879      treeEchoContext(fp, stdout, 40);
880      printf("\n");
881      return  FALSE;
882    }
883   
884  if((ch = getc(fp)) != EOF)
885    {
886      if(ch == '[')
887        {
888          if(fscanf(fp, "%d", branchLabel) != 1)
889            goto handleError;         
890
891          //printf("Branch label: %d\n", *branchLabel);
892         
893          if((ch = getc(fp)) != ']')
894            {
895            handleError:
896              printf("ERROR: treeProcessLength: Problem reading branch label\n");
897              treeEchoContext(fp, stdout, 40);
898              printf("\n");
899              return FALSE;
900            }       
901
902          if(storeBranchLabels)
903            tr->branchLabelCounter = tr->branchLabelCounter + 1;
904
905        }
906      else
907        (void)ungetc(ch, fp);
908    }
909 
910  return  TRUE;
911}
912
913
914static int treeFlushLen (FILE  *fp, tree *tr)
915{
916  double 
917    dummy; 
918  int 
919    dummyLabel,     
920    ch;
921 
922  ch = treeGetCh(fp);
923 
924  if (ch == ':') 
925    {
926      ch = treeGetCh(fp);
927     
928      ungetc(ch, fp);
929      if(!treeProcessLength(fp, &dummy, &dummyLabel, FALSE, tr)) 
930        return 0;
931      return 1;   
932    }
933 
934 
935 
936  if (ch != EOF) (void) ungetc(ch, fp);
937  return 1;
938} 
939
940
941
942
943
944static boolean treeNeedCh (FILE *fp, int c1, const char *where)
945{
946  int  c2;
947 
948  if ((c2 = treeGetCh(fp)) == c1)  return TRUE;
949 
950  printf("ERROR: Expecting '%c' %s tree; found:", c1, where);
951  if (c2 == EOF) 
952    {
953      printf("End-of-File");
954    }
955  else 
956    {           
957      ungetc(c2, fp);
958      treeEchoContext(fp, stdout, 40);
959    }
960  putchar('\n');
961   
962  printf("RAxML may be expecting to read a tree that contains branch lengths\n");
963
964  return FALSE;
965} 
966
967
968
969static boolean addElementLen (FILE *fp, tree *tr, nodeptr p, boolean readBranchLengths, boolean readNodeLabels, int *lcount, analdef *adef, boolean storeBranchLabels)
970{   
971  nodeptr  q;
972  int      n, ch, fres;
973 
974  if ((ch = treeGetCh(fp)) == '(') 
975    { 
976      n = (tr->nextnode)++;
977      if (n > 2*(tr->mxtips) - 2) 
978        {
979          if (tr->rooted || n > 2*(tr->mxtips) - 1) 
980            {
981              printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
982              printf("       Deepest splitting should be a trifurcation.\n");
983              return FALSE;
984            }
985          else 
986            {
987              if(readNodeLabels)
988                {
989                  printf("The program will exit with an error in the next source code line\n");
990                  printf("You are probably trying to read in rooted trees with a RAxML option \n");
991                  printf("that for some reason expects unrooted binary trees\n\n");
992                }
993
994              assert(!readNodeLabels);
995              tr->rooted = TRUE;
996            }
997        }
998     
999      q = tr->nodep[n];
1000
1001      if (! addElementLen(fp, tr, q->next, readBranchLengths, readNodeLabels, lcount, adef, storeBranchLabels))        return FALSE;
1002      if (! treeNeedCh(fp, ',', "in"))             return FALSE;
1003      if (! addElementLen(fp, tr, q->next->next, readBranchLengths, readNodeLabels, lcount, adef, storeBranchLabels))  return FALSE;
1004      if (! treeNeedCh(fp, ')', "in"))             return FALSE;
1005     
1006      if(readNodeLabels)
1007        {
1008          char label[64];
1009          int support;
1010
1011          if(treeGetLabel (fp, label, 10))
1012            {   
1013              int val = sscanf(label, "%d", &support);
1014     
1015              assert(val == 1);
1016
1017              /*printf("LABEL %s Number %d\n", label, support);*/
1018              p->support = q->support = support;
1019              /*printf("%d %d %d %d\n", p->support, q->support, p->number, q->number);*/
1020              assert(p->number > tr->mxtips && q->number > tr->mxtips);
1021              *lcount = *lcount + 1;
1022            }
1023        }
1024      else     
1025        (void) treeFlushLabel(fp);
1026    }
1027  else 
1028    {   
1029      ungetc(ch, fp);
1030      if ((n = treeFindTipName(fp, tr, TRUE)) <= 0)          return FALSE;
1031      q = tr->nodep[n];
1032      if (tr->start->number > n)  tr->start = q;
1033      (tr->ntips)++;
1034    }
1035 
1036  if(readBranchLengths)
1037    {
1038      double 
1039        branch;
1040     
1041      int 
1042        startCounter = tr->branchLabelCounter,
1043        endCounter,
1044        branchLabel = -1;
1045     
1046      if (! treeNeedCh(fp, ':', "in"))                 return FALSE;
1047      if (! treeProcessLength(fp, &branch, &branchLabel, storeBranchLabels, tr))           
1048        return FALSE;
1049
1050      endCounter = tr->branchLabelCounter;
1051     
1052      /*printf("Branch %8.20f %d\n", branch, tr->numBranches);*/
1053      if(adef->mode == CLASSIFY_ML)
1054        {
1055          double 
1056            x[NUM_BRANCHES];
1057         
1058          assert(tr->NumberOfModels == 1);
1059          assert(adef->useBinaryModelFile);
1060          assert(tr->numBranches == 1);
1061
1062          x[0] = exp(-branch / tr->fracchange);           
1063
1064          hookup(p, q, x, tr->numBranches);
1065        }
1066      else
1067        hookup(p, q, &branch, tr->numBranches);
1068
1069      if(storeBranchLabels && (endCounter > startCounter))
1070        {
1071          assert(!isTip(p->number, tr->mxtips) && !isTip(q->number, tr->mxtips));
1072          assert(branchLabel >= 0);
1073          p->support = q->support = branchLabel;
1074        }
1075    }
1076  else
1077    {
1078      fres = treeFlushLen(fp, tr);
1079      if(!fres) return FALSE;
1080     
1081      hookupDefault(p, q, tr->numBranches);
1082    }
1083  return TRUE;         
1084} 
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096static nodeptr uprootTree (tree *tr, nodeptr p, boolean readBranchLengths, boolean readConstraint)
1097{
1098  nodeptr  q, r, s, start;
1099  int      n;
1100
1101  for(int i = tr->mxtips + 1; i < 2 * tr->mxtips - 1; i++)
1102    assert(i == tr->nodep[i]->number);
1103
1104 
1105 
1106
1107  if(isTip(p->number, tr->mxtips) || p->back) 
1108    {
1109      printf("ERROR: Unable to uproot tree.\n");
1110      printf("       Inappropriate node marked for removal.\n");
1111      assert(0);
1112    }
1113 
1114  assert(p->back == (nodeptr)NULL);
1115 
1116  tr->nextnode = tr->nextnode - 1;
1117
1118  assert(tr->nextnode < 2 * tr->mxtips);
1119 
1120  n = tr->nextnode;               
1121 
1122  assert(tr->nodep[tr->nextnode]);
1123
1124  if (n != tr->mxtips + tr->ntips - 1) 
1125    {
1126      printf("ERROR: Unable to uproot tree.  Inconsistent\n");
1127      printf("       number of tips and nodes for rooted tree.\n");
1128      assert(0);
1129    }
1130
1131  q = p->next->back;                  /* remove p from tree */
1132  r = p->next->next->back;
1133  assert(p->back == (nodeptr)NULL);
1134   
1135  if(readBranchLengths)
1136    {
1137      double b[NUM_BRANCHES];
1138      for(int i = 0; i < tr->numBranches; i++)
1139        {
1140          b[i] = (r->z[i] + q->z[i]);     
1141        }
1142      hookup (q, r, b, tr->numBranches);
1143    }
1144  else   
1145    hookupDefault(q, r, tr->numBranches);   
1146
1147  tr->leftRootNode = p->next->back;
1148  tr->rightRootNode = p->next->next->back;
1149
1150  if(readConstraint && tr->grouped)
1151    {   
1152      if(tr->constraintVector[p->number] != 0)
1153        {
1154          printf("Root node to remove should have top-level grouping of 0\n");
1155          assert(0);
1156        }
1157    } 
1158 
1159  assert(!(isTip(r->number, tr->mxtips) && isTip(q->number, tr->mxtips))); 
1160
1161  assert(p->number > tr->mxtips);
1162
1163  if(tr->ntips > 2 && p->number != n) 
1164    {               
1165      q = tr->nodep[n];            /* transfer last node's conections to p */
1166      r = q->next;
1167      s = q->next->next;
1168           
1169      if(readConstraint && tr->grouped) 
1170        tr->constraintVector[p->number] = tr->constraintVector[q->number];       
1171     
1172      hookup(p,             q->back, q->z, tr->numBranches);   /* move connections to p */
1173      hookup(p->next,       r->back, r->z, tr->numBranches);
1174      hookup(p->next->next, s->back, s->z, tr->numBranches); 
1175
1176      if(q == tr->leftRootNode || q == tr->rightRootNode)
1177        {
1178          if(q == tr->leftRootNode)
1179            {
1180              if(p->back == tr->rightRootNode)
1181                tr->leftRootNode = p;
1182              else
1183                {
1184                   if(p->next->back == tr->rightRootNode)
1185                     tr->leftRootNode = p->next;
1186                   else
1187                     {
1188                       if(p->next->next->back == tr->rightRootNode)
1189                         tr->leftRootNode = p->next->next;
1190                       else
1191                         assert(0);
1192                     }
1193                }
1194            }
1195          else
1196            {
1197              assert(q == tr->rightRootNode);
1198
1199              if(p->back == tr->leftRootNode)
1200                tr->rightRootNode = p;
1201              else
1202                {
1203                   if(p->next->back == tr->leftRootNode)
1204                     tr->rightRootNode = p->next;
1205                   else
1206                     {
1207                       if(p->next->next->back == tr->leftRootNode)
1208                         tr->rightRootNode = p->next->next;
1209                       else
1210                         assert(0);
1211                     }
1212                }
1213            }
1214        }
1215     
1216      q->back = q->next->back = q->next->next->back = (nodeptr) NULL;
1217    }
1218  else   
1219    {
1220      assert(tr->ntips > 2);     
1221      p->back = p->next->back = p->next->next->back = (nodeptr) NULL;
1222    }
1223
1224 
1225 
1226  assert(tr->ntips > 2);
1227 
1228  start = findAnyTip(tr->nodep[tr->mxtips + 1], tr->mxtips);
1229 
1230  assert(isTip(start->number, tr->mxtips));
1231  tr->rooted = FALSE;
1232  return  start;
1233}
1234
1235
1236int treeReadLen (FILE *fp, tree *tr, boolean readBranches, boolean readNodeLabels, boolean topologyOnly, analdef *adef, boolean completeTree, boolean storeBranchLabels)
1237{
1238  nodeptr 
1239    p;
1240 
1241  int 
1242    i, 
1243    ch, 
1244    lcount = 0; 
1245
1246  tr->branchLabelCounter = 0;
1247
1248  for (i = 1; i <= tr->mxtips; i++) 
1249    {
1250      tr->nodep[i]->back = (node *) NULL; 
1251      if(topologyOnly)
1252        tr->nodep[i]->support = -1;
1253    }
1254
1255  for(i = tr->mxtips + 1; i < 2 * tr->mxtips; i++)
1256    {
1257      tr->nodep[i]->back = (nodeptr)NULL;
1258      tr->nodep[i]->next->back = (nodeptr)NULL;
1259      tr->nodep[i]->next->next->back = (nodeptr)NULL;
1260      tr->nodep[i]->number = i;
1261      tr->nodep[i]->next->number = i;
1262      tr->nodep[i]->next->next->number = i;
1263
1264      if(topologyOnly)
1265        {
1266          tr->nodep[i]->support = -2;
1267          tr->nodep[i]->next->support = -2;
1268          tr->nodep[i]->next->next->support = -2;
1269        }
1270    }
1271
1272  if(topologyOnly)
1273    tr->start       = tr->nodep[tr->mxtips];
1274  else
1275    tr->start       = tr->nodep[1];
1276
1277  tr->ntips       = 0;
1278  tr->nextnode    = tr->mxtips + 1;     
1279 
1280  for(i = 0; i < tr->numBranches; i++)
1281    tr->partitionSmoothed[i] = FALSE;
1282 
1283  tr->rooted      = FALSE; 
1284  tr->wasRooted   = FALSE;
1285
1286  p = tr->nodep[(tr->nextnode)++]; 
1287 
1288  while((ch = treeGetCh(fp)) != '(');
1289     
1290  if(!topologyOnly)
1291    {
1292      if(adef->mode != CLASSIFY_ML)
1293        {
1294          if(adef->mode != OPTIMIZE_BR_LEN_SCALER)
1295            assert(readBranches == FALSE && readNodeLabels == FALSE);
1296          else           
1297            assert(readBranches == TRUE && readNodeLabels == FALSE);           
1298        }
1299      else
1300        {
1301          if(adef->useBinaryModelFile)
1302            assert(readBranches == TRUE && readNodeLabels == FALSE);           
1303          else
1304            assert(readBranches == FALSE && readNodeLabels == FALSE);
1305        }
1306    }
1307 
1308       
1309  if (! addElementLen(fp, tr, p, readBranches, readNodeLabels, &lcount, adef, storeBranchLabels))                 
1310    assert(0);
1311  if (! treeNeedCh(fp, ',', "in"))               
1312    assert(0);
1313  if (! addElementLen(fp, tr, p->next, readBranches, readNodeLabels, &lcount, adef, storeBranchLabels))
1314    assert(0);
1315  if (! tr->rooted) 
1316    {
1317      if ((ch = treeGetCh(fp)) == ',') 
1318        { 
1319          if (! addElementLen(fp, tr, p->next->next, readBranches, readNodeLabels, &lcount, adef, storeBranchLabels))
1320            assert(0);     
1321        }
1322      else 
1323        {         
1324          /*  A rooted format */
1325         
1326          tr->rooted = TRUE;
1327          tr->wasRooted     = TRUE;
1328         
1329          if (ch != EOF)  (void) ungetc(ch, fp);
1330        }       
1331    }
1332  else 
1333    {           
1334      p->next->next->back = (nodeptr) NULL;
1335      tr->wasRooted     = TRUE;   
1336    }
1337
1338  if(!tr->rooted && adef->mode == ANCESTRAL_STATES)
1339    {
1340      printf("Error: The ancestral state computation mode requires a rooted tree as input, exiting ....\n");
1341      exit(0);
1342    }
1343
1344  if (! treeNeedCh(fp, ')', "in"))               
1345    assert(0);
1346
1347  if(topologyOnly)
1348    assert(!(tr->rooted && readNodeLabels));
1349
1350  (void) treeFlushLabel(fp);
1351 
1352  if (! treeFlushLen(fp, tr))                         
1353    assert(0);
1354 
1355  if (! treeNeedCh(fp, ';', "at end of"))       
1356    assert(0);
1357 
1358  if (tr->rooted) 
1359    {     
1360      assert(!readNodeLabels);
1361
1362      p->next->next->back = (nodeptr) NULL;     
1363      tr->start = uprootTree(tr, p->next->next, readBranches, FALSE);     
1364
1365       
1366      /*tr->leftRootNode  = p->back;
1367        tr->rightRootNode = p->next->back;   
1368      */
1369
1370      if (! tr->start)                             
1371        {
1372          printf("FATAL ERROR UPROOTING TREE\n");
1373          assert(0);
1374        }   
1375    }
1376  else   
1377    tr->start = findAnyTip(p, tr->rdta->numsp);   
1378 
1379   if(!topologyOnly || adef->mode == CLASSIFY_MP)
1380    {     
1381      assert(tr->ntips <= tr->mxtips);
1382     
1383
1384      if(tr->ntips < tr->mxtips)
1385        {
1386          if(completeTree)
1387            {
1388              printBothOpen("Hello this is your friendly RAxML tree parsing routine\n");
1389              printBothOpen("The RAxML option you are uisng requires to read in only complete trees\n");
1390              printBothOpen("with %d taxa, there is at least one tree with %d taxa though ... exiting\n", tr->mxtips, tr->ntips);
1391              exit(-1);
1392            }
1393          else
1394            {
1395              if(adef->computeDistance)
1396                {
1397                  printBothOpen("Error: pairwise distance computation only allows for complete, i.e., containing all taxa\n");
1398                  printBothOpen("bifurcating starting trees\n");
1399                  exit(-1);
1400                }     
1401              if(adef->mode == CLASSIFY_ML || adef->mode == CLASSIFY_MP)
1402                {       
1403                  printBothOpen("RAxML placement algorithm: You provided a reference tree with %d taxa; alignmnet has %d taxa\n", tr->ntips, tr->mxtips);                 
1404                  printBothOpen("%d query taxa will be placed using %s\n", tr->mxtips - tr->ntips, (adef->mode == CLASSIFY_ML)?"maximum likelihood":"parsimony");
1405                  if(adef->mode == CLASSIFY_ML)
1406                    classifyML(tr, adef);         
1407                  else
1408                    {
1409                      assert(adef->mode == CLASSIFY_MP);
1410                      classifyMP(tr, adef);
1411                    }
1412                }
1413              else
1414                {
1415                  printBothOpen("You provided an incomplete starting tree %d alignmnet has %d taxa\n", tr->ntips, tr->mxtips);   
1416                  makeParsimonyTreeIncomplete(tr, adef);                         
1417                }   
1418            }
1419        }
1420      else
1421        {
1422          if(adef->mode == PARSIMONY_ADDITION)
1423            {
1424              printBothOpen("Error you want to add sequences to a trees via MP stepwise addition, but \n");
1425              printBothOpen("you have provided an input tree that already contains all taxa\n");
1426              exit(-1);
1427            }
1428          if(adef->mode == CLASSIFY_ML || adef->mode == CLASSIFY_MP)
1429            {
1430              printBothOpen("Error you want to place query sequences into a tree using %s, but\n", tr->mxtips - tr->ntips, (adef->mode == CLASSIFY_ML)?"maximum likelihood":"parsimony");
1431              printBothOpen("you have provided an input tree that already contains all taxa\n");
1432              exit(-1);
1433            }
1434        }
1435   
1436      onlyInitrav(tr, tr->start);
1437    }
1438 
1439 
1440  return lcount;
1441}
1442
1443
1444
1445/********************************MULTIFURCATIONS************************************************/
1446
1447
1448static boolean  addElementLenMULT (FILE *fp, tree *tr, nodeptr p, int partitionCounter)
1449{ 
1450  nodeptr  q, r, s;
1451  int      n, ch, fres, rn;
1452  double randomResolution;
1453  int old;
1454   
1455  tr->constraintVector[p->number] = partitionCounter; 
1456
1457  if ((ch = treeGetCh(fp)) == '(') 
1458    {
1459      partCount++;
1460      old = partCount;       
1461     
1462      n = (tr->nextnode)++;
1463      if (n > 2*(tr->mxtips) - 2) 
1464        {
1465          if (tr->rooted || n > 2*(tr->mxtips) - 1) 
1466            {
1467              printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
1468              printf("       Deepest splitting should be a trifurcation.\n");
1469              return FALSE;
1470            }
1471          else 
1472            {
1473              tr->rooted = TRUE;           
1474            }
1475        }
1476      q = tr->nodep[n];
1477      tr->constraintVector[q->number] = partCount;
1478      if (! addElementLenMULT(fp, tr, q->next, old))        return FALSE;
1479      if (! treeNeedCh(fp, ',', "in"))             return FALSE;
1480      if (! addElementLenMULT(fp, tr, q->next->next, old))  return FALSE;
1481                 
1482      hookupDefault(p, q, tr->numBranches);
1483
1484      while((ch = treeGetCh(fp)) == ',')
1485        { 
1486          n = (tr->nextnode)++;
1487          if (n > 2*(tr->mxtips) - 2) 
1488            {
1489              if (tr->rooted || n > 2*(tr->mxtips) - 1) 
1490                {
1491                  printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
1492                  printf("       Deepest splitting should be a trifurcation.\n");
1493                  return FALSE;
1494                }
1495              else 
1496                {
1497                  tr->rooted = TRUE;
1498                }
1499            }
1500          r = tr->nodep[n];
1501          tr->constraintVector[r->number] = partCount;   
1502
1503          rn = randomInt(10000);
1504          if(rn == 0) 
1505            randomResolution = 0;
1506          else 
1507            randomResolution = ((double)rn)/10000.0;
1508                 
1509           if(randomResolution < 0.5)
1510            {       
1511              s = q->next->back;             
1512              r->back = q->next;
1513              q->next->back = r;             
1514              r->next->back = s;
1515              s->back = r->next;             
1516              addElementLenMULT(fp, tr, r->next->next, old);         
1517            }
1518          else
1519            {     
1520              s = q->next->next->back;       
1521              r->back = q->next->next;
1522              q->next->next->back = r;       
1523              r->next->back = s;
1524              s->back = r->next;             
1525              addElementLenMULT(fp, tr, r->next->next, old);         
1526            }                     
1527        }       
1528
1529      if(ch != ')')
1530        {
1531          printf("Missing /) in treeReadLenMULT\n");
1532          exit(-1);             
1533        }
1534       
1535
1536
1537      (void) treeFlushLabel(fp);
1538    }
1539  else 
1540    {                             
1541      ungetc(ch, fp);
1542      if ((n = treeFindTipName(fp, tr, TRUE)) <= 0)          return FALSE;
1543      q = tr->nodep[n];     
1544      tr->constraintVector[q->number] = partitionCounter;
1545
1546      if (tr->start->number > n)  tr->start = q;
1547      (tr->ntips)++;
1548      hookupDefault(p, q, tr->numBranches);
1549    }
1550 
1551  fres = treeFlushLen(fp, tr);
1552  if(!fres) return FALSE;
1553   
1554  return TRUE;         
1555} 
1556
1557
1558
1559
1560
1561boolean treeReadLenMULT (FILE *fp, tree *tr, analdef *adef)
1562{
1563  nodeptr  p, r, s;
1564  int      i, ch, n, rn;
1565  int partitionCounter = 0;
1566  double randomResolution;
1567
1568  srand((unsigned int) time(NULL));
1569 
1570  for(i = 0; i < 2 * tr->mxtips; i++)
1571    tr->constraintVector[i] = -1;
1572
1573  for (i = 1; i <= tr->mxtips; i++) 
1574    tr->nodep[i]->back = (node *) NULL;
1575
1576  for(i = tr->mxtips + 1; i < 2 * tr->mxtips; i++)
1577    {
1578      tr->nodep[i]->back = (nodeptr)NULL;
1579      tr->nodep[i]->next->back = (nodeptr)NULL;
1580      tr->nodep[i]->next->next->back = (nodeptr)NULL;
1581      tr->nodep[i]->number = i;
1582      tr->nodep[i]->next->number = i;
1583      tr->nodep[i]->next->next->number = i;
1584    }
1585
1586
1587  tr->start       = tr->nodep[tr->mxtips];
1588  tr->ntips       = 0;
1589  tr->nextnode    = tr->mxtips + 1;
1590 
1591  for(i = 0; i < tr->numBranches; i++)
1592    tr->partitionSmoothed[i] = FALSE;
1593
1594  tr->rooted      = FALSE;
1595 
1596  p = tr->nodep[(tr->nextnode)++]; 
1597  while((ch = treeGetCh(fp)) != '(');
1598     
1599  if (! addElementLenMULT(fp, tr, p, partitionCounter))                 return FALSE;
1600  if (! treeNeedCh(fp, ',', "in"))                return FALSE;
1601  if (! addElementLenMULT(fp, tr, p->next, partitionCounter))           return FALSE;
1602  if (! tr->rooted) 
1603    {
1604      if ((ch = treeGetCh(fp)) == ',') 
1605        {       
1606          if (! addElementLenMULT(fp, tr, p->next->next, partitionCounter)) return FALSE;
1607
1608          while((ch = treeGetCh(fp)) == ',')
1609            { 
1610              n = (tr->nextnode)++;
1611              assert(n <= 2*(tr->mxtips) - 2);
1612       
1613              r = tr->nodep[n]; 
1614              tr->constraintVector[r->number] = partitionCounter;         
1615             
1616              rn = randomInt(10000);
1617              if(rn == 0) 
1618                randomResolution = 0;
1619              else 
1620                randomResolution = ((double)rn)/10000.0;
1621
1622
1623              if(randomResolution < 0.5)
1624                {       
1625                  s = p->next->next->back;               
1626                  r->back = p->next->next;
1627                  p->next->next->back = r;               
1628                  r->next->back = s;
1629                  s->back = r->next;             
1630                  addElementLenMULT(fp, tr, r->next->next, partitionCounter);   
1631                }
1632              else
1633                {
1634                  s = p->next->back;             
1635                  r->back = p->next;
1636                  p->next->back = r;             
1637                  r->next->back = s;
1638                  s->back = r->next;             
1639                  addElementLenMULT(fp, tr, r->next->next, partitionCounter);
1640                }
1641            }                             
1642
1643          if(ch != ')')
1644            {
1645              printf("Missing /) in treeReadLenMULT\n");
1646              exit(-1);                               
1647            }
1648          else
1649            ungetc(ch, fp);
1650        }
1651      else 
1652        { 
1653          tr->rooted = TRUE;
1654          if (ch != EOF)  (void) ungetc(ch, fp);
1655        }       
1656    }
1657  else 
1658    {
1659      p->next->next->back = (nodeptr) NULL;
1660    }
1661   
1662  if (! treeNeedCh(fp, ')', "in"))                return FALSE;
1663  (void) treeFlushLabel(fp);
1664  if (! treeFlushLen(fp, tr))                         return FALSE;
1665   
1666  if (! treeNeedCh(fp, ';', "at end of"))       return FALSE;
1667 
1668
1669  if (tr->rooted) 
1670    {       
1671      p->next->next->back = (nodeptr) NULL;
1672      tr->start = uprootTree(tr, p->next->next, FALSE, TRUE);
1673      if (! tr->start)                              return FALSE;
1674    }
1675  else 
1676    {     
1677      tr->start = findAnyTip(p, tr->rdta->numsp);
1678    }
1679
1680 
1681 
1682
1683  if(tr->ntips < tr->mxtips)         
1684    makeParsimonyTreeIncomplete(tr, adef);         
1685
1686
1687  if(!adef->rapidBoot)
1688    onlyInitrav(tr, tr->start);
1689  return TRUE; 
1690}
1691
1692
1693
1694
1695
1696
1697void getStartingTree(tree *tr, analdef *adef)
1698{
1699  tr->likelihood = unlikely;
1700 
1701  if(adef->restart) 
1702    {                       
1703      INFILE = myfopen(tree_file, "rb");       
1704                               
1705      if(!adef->grouping)       
1706        {
1707          switch(adef->mode)
1708            {
1709            case ANCESTRAL_STATES:         
1710              assert(!tr->saveMemory);
1711
1712              tr->leftRootNode  = (nodeptr)NULL;
1713              tr->rightRootNode = (nodeptr)NULL;
1714
1715              treeReadLen(INFILE, tr, FALSE, FALSE, FALSE, adef, TRUE, FALSE);
1716
1717              assert(tr->leftRootNode && tr->rightRootNode);
1718              break;
1719            case CLASSIFY_MP:
1720              treeReadLen(INFILE, tr, TRUE, FALSE, TRUE, adef, FALSE, FALSE);
1721              break;
1722            case OPTIMIZE_BR_LEN_SCALER:
1723              treeReadLen(INFILE, tr, TRUE, FALSE, FALSE, adef, TRUE, FALSE);
1724              break;
1725            case CLASSIFY_ML:
1726              if(adef->useBinaryModelFile)
1727                {
1728                  if(tr->saveMemory)                             
1729                    treeReadLen(INFILE, tr, TRUE, FALSE, TRUE, adef, FALSE, FALSE);                           
1730                  else             
1731                    treeReadLen(INFILE, tr, TRUE, FALSE, FALSE, adef, FALSE, FALSE);
1732                }
1733              else
1734                {
1735                  if(tr->saveMemory)                             
1736                    treeReadLen(INFILE, tr, FALSE, FALSE, TRUE, adef, FALSE, FALSE);                           
1737                  else             
1738                    treeReadLen(INFILE, tr, FALSE, FALSE, FALSE, adef, FALSE, FALSE);
1739                }
1740              break;
1741            default:         
1742              if(tr->saveMemory)                                 
1743                treeReadLen(INFILE, tr, FALSE, FALSE, TRUE, adef, FALSE, FALSE);                               
1744              else                 
1745                treeReadLen(INFILE, tr, FALSE, FALSE, FALSE, adef, FALSE, FALSE);
1746              break;
1747            }
1748        }
1749      else
1750        {
1751          assert(adef->mode != ANCESTRAL_STATES);
1752
1753          partCount = 0;
1754          if (! treeReadLenMULT(INFILE, tr, adef))
1755            exit(-1);
1756        }                                                                         
1757
1758      if(adef->mode == PARSIMONY_ADDITION)
1759        return; 
1760
1761      if(adef->mode != CLASSIFY_MP)
1762        {
1763          if(adef->mode == OPTIMIZE_BR_LEN_SCALER)
1764            {
1765              assert(tr->numBranches == tr->NumberOfModels);
1766              scaleBranches(tr, TRUE);
1767              evaluateGenericInitrav(tr, tr->start);                                 
1768            }
1769          else
1770            {
1771              evaluateGenericInitrav(tr, tr->start); 
1772              treeEvaluate(tr, 1);
1773            }
1774        }
1775               
1776      fclose(INFILE);
1777    }
1778  else
1779    { 
1780      assert(adef->mode != PARSIMONY_ADDITION &&
1781             adef->mode != MORPH_CALIBRATOR   &&
1782             adef->mode != ANCESTRAL_STATES   &&
1783             adef->mode != OPTIMIZE_BR_LEN_SCALER);
1784
1785      if(adef->randomStartingTree)       
1786        makeRandomTree(tr, adef);                                         
1787      else
1788        makeParsimonyTree(tr, adef);                                           
1789     
1790      if(adef->startingTreeOnly)
1791        {
1792          printStartingTree(tr, adef, TRUE);
1793          exit(0);
1794        }
1795      else               
1796        printStartingTree(tr, adef, FALSE);                     
1797           
1798     
1799      evaluateGenericInitrav(tr, tr->start);   
1800
1801     
1802     
1803      treeEvaluate(tr, 1);               
1804
1805     
1806     
1807    }         
1808
1809  tr->start = tr->nodep[1];
1810}
1811
1812
1813
1814/****** functions for reading true multi-furcating trees *****/
1815
1816/****************************************************************************/
1817
1818static boolean addMultifurcation (FILE *fp, tree *tr, nodeptr _p, analdef *adef, int *nextnode)
1819{   
1820  nodeptr 
1821    p, 
1822    initial_p;
1823 
1824  int     
1825    n, 
1826    ch, 
1827    fres;
1828 
1829  if ((ch = treeGetCh(fp)) == '(') 
1830    { 
1831      int 
1832        i = 0;     
1833     
1834      initial_p = p = tr->nodep[*nextnode];     
1835      *nextnode = *nextnode + 1;
1836
1837      do
1838        {                 
1839          p->next = tr->nodep[*nextnode];       
1840          *nextnode = *nextnode + 1;
1841          p = p->next;
1842       
1843          addMultifurcation(fp, tr, p, adef, nextnode);   
1844          i++;
1845        } 
1846      while((ch = treeGetCh(fp)) == ',');
1847
1848      ungetc(ch, fp);
1849
1850      p->next = initial_p;
1851           
1852      if (! treeNeedCh(fp, ')', "in"))               
1853        assert(0);                   
1854
1855      treeFlushLabel(fp);
1856    }
1857  else 
1858    {   
1859      ungetc(ch, fp);
1860      if ((n = treeFindTipName(fp, tr, FALSE)) <= 0)         
1861        return FALSE;
1862      p = tr->nodep[n];
1863      initial_p = p;
1864      tr->start = p;
1865      (tr->ntips)++;
1866    }
1867 
1868 
1869  fres = treeFlushLen(fp, tr);
1870  if(!fres) 
1871    return FALSE;
1872     
1873  hookupDefault(initial_p, _p, tr->numBranches);
1874 
1875  return TRUE;         
1876} 
1877
1878static void printMultiFurc(nodeptr p, tree *tr)
1879{ 
1880  if(isTip(p->number, tr->mxtips))   
1881    printf("%s", tr->nameList[p->number]);   
1882  else
1883    {
1884      nodeptr
1885        q = p->next;     
1886     
1887      printf("(");
1888
1889      while(q != p)
1890        {       
1891          printMultiFurc(q->back, 
1892                         tr);
1893          q = q->next;
1894          if(q != p)
1895            printf(",");
1896        }
1897
1898      printf(")");
1899    }
1900}
1901
1902void allocateMultifurcations(tree *tr, tree *smallTree)
1903{ 
1904  int
1905    i,
1906    tips,
1907    inter; 
1908
1909  smallTree->numBranches = tr->numBranches;
1910
1911  smallTree->mxtips = tr->mxtips;
1912
1913  smallTree->nameHash = tr->nameHash;
1914
1915  tips  = tr->mxtips;
1916  inter = tr->mxtips - 1;
1917 
1918  smallTree->nodep = (nodeptr *)rax_malloc((tips + 3 * inter) * sizeof(nodeptr));
1919
1920  smallTree->nodep[0] = (node *) NULL;
1921
1922  for (i = 1; i <= tips; i++)
1923    {
1924      smallTree->nodep[i] = (nodeptr)rax_malloc(sizeof(node));     
1925      memcpy(smallTree->nodep[i], tr->nodep[i], sizeof(node));
1926      smallTree->nodep[i]->back = (node *) NULL;
1927      smallTree->nodep[i]->next = (node *) NULL;
1928    }
1929
1930  for(i = tips + 1; i < tips + 3 * inter; i++)
1931    {     
1932      smallTree->nodep[i] = (nodeptr)rax_malloc(sizeof(node));
1933      smallTree->nodep[i]->number = i;     
1934      smallTree->nodep[i]->back = (node *) NULL;
1935      smallTree->nodep[i]->next = (node *) NULL;
1936    }
1937 
1938}
1939
1940void freeMultifurcations(tree *tr)
1941{
1942  int
1943    i,
1944    tips  = tr->mxtips,
1945    inter = tr->mxtips - 1; 
1946
1947  for (i = 1; i < tips + 3 * inter; i++)
1948    rax_free(tr->nodep[i]);
1949 
1950  rax_free(tr->nodep);
1951}
1952
1953static void relabelInnerNodes(nodeptr p, tree *tr, int n, int *nodeCounter)
1954{
1955  if(isTip(p->number, tr->mxtips))
1956    {
1957      assert(0);
1958    }
1959  else
1960    {
1961      nodeptr
1962        q = p->next;
1963
1964      int 
1965        _n = n;
1966
1967      tr->nodep[p->number]->number = n;
1968      p->x = 1;
1969     
1970      while(q != p)
1971        {
1972          tr->nodep[q->number]->number = n;     
1973          q->x = 0;
1974         
1975          if(!isTip(q->back->number, tr->mxtips))
1976            {
1977              _n++;
1978              *nodeCounter = *nodeCounter + 1;
1979              relabelInnerNodes(q->back, tr, _n, nodeCounter);
1980            }
1981          q = q->next;
1982        } 
1983    }
1984}
1985
1986
1987int readMultifurcatingTree(FILE *fp, tree *tr, analdef *adef)
1988{
1989  nodeptr 
1990    p,
1991    initial_p;
1992 
1993  int   
1994    innerBranches = 0,
1995    nextnode,
1996    i, 
1997    ch,
1998    tips  = tr->mxtips,
1999    inter = tr->mxtips - 1; 
2000 
2001  //clean up before parsing !
2002   
2003  for (i = 1; i < tips + 3 * inter; i++)     
2004    {     
2005      tr->nodep[i]->back = (node *) NULL;
2006      tr->nodep[i]->next = (node *) NULL;   
2007      tr->nodep[i]->x = 0;
2008    }
2009 
2010  for(i = tips + 1; i < tips + 3 * inter; i++)   
2011    tr->nodep[i]->number = i;
2012
2013  tr->ntips = 0;
2014  nextnode  = tr->mxtips + 1;         
2015
2016  while((ch = treeGetCh(fp)) != '(');           
2017
2018  i = 0;
2019
2020  do
2021    {         
2022      if(i == 0)
2023        {
2024          initial_p = p = tr->nodep[nextnode];   
2025          nextnode++;
2026        }
2027      else
2028        {
2029          p->next = tr->nodep[nextnode];                 
2030          p = p->next; 
2031          nextnode++;
2032        }
2033     
2034      addMultifurcation(fp, tr, p, adef, &nextnode);       
2035                   
2036      i++;
2037    } 
2038  while((ch = treeGetCh(fp)) == ',');
2039   
2040  if(i < 3)
2041    {
2042      printBothOpen("You need to provide unrooted input trees!\n");
2043      assert(0);
2044    }
2045 
2046  ungetc(ch, fp);
2047 
2048  p->next = initial_p;
2049 
2050  if (! treeNeedCh(fp, ')', "in"))               
2051    assert(0); 
2052
2053  (void)treeFlushLabel(fp);
2054 
2055  if (! treeFlushLen(fp, tr))                         
2056    assert(0);
2057 
2058  if (! treeNeedCh(fp, ';', "at end of"))       
2059    assert(0);
2060   
2061  //printf("%d tips found, %d inner nodes used start %d maxtips %d\n", tr->ntips, tr->nextnode - tr->mxtips, tr->start->number, tr->mxtips);
2062
2063  assert(isTip(tr->start->number, tr->mxtips));
2064
2065  relabelInnerNodes(tr->start->back, tr, tr->mxtips + 1, &innerBranches);
2066
2067  //printf("Inner branches %d\n", innerBranches);
2068
2069  if(0)
2070    {
2071      printf("(");
2072      printMultiFurc(tr->start, tr);
2073      printf(",");
2074      printMultiFurc(tr->start->back, tr);
2075      printf(");\n");
2076    }
2077
2078  return innerBranches;
2079}
2080
Note: See TracBrowser for help on using the repository browser.