source: tags/arb-6.0/GDE/RAxML/bipartitionList.c

Last change on this file was 10409, checked in by aboeckma, 11 years ago

Updated raxml to current version
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 114.8 KB
Line 
1/*  RAxML-HPC, a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright March 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 *  stamatak@ics.forth.gr
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 * 
27 *  Alexandros Stamatakis: "An Efficient Program for phylogenetic Inference Using Simulated Annealing".
28 *  Proceedings of IPDPS2005,  Denver, Colorado, April 2005.
29 * 
30 *  AND
31 *
32 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
33 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
34 */
35
36
37#ifndef WIN32 
38#include <sys/times.h>
39#include <sys/types.h>
40#include <sys/time.h>
41#include <unistd.h> 
42#endif
43
44#include <limits.h>
45#include <math.h>
46#include <time.h>
47#include <stdlib.h>
48#include <stdio.h>
49#include <ctype.h>
50#include <string.h>
51#include <stdint.h>
52#include "axml.h"
53
54#ifdef __SIM_SSE3
55
56#include <xmmintrin.h>
57#include <pmmintrin.h>
58
59#endif
60
61#ifdef _USE_PTHREADS
62#include <pthread.h>
63#endif
64
65#ifdef _WAYNE_MPI
66#include <mpi.h>
67extern int processID;
68extern int processes;
69#endif
70
71#define _NEW_MRE
72
73extern FILE *INFILE;
74extern char run_id[128];
75extern char workdir[1024];
76extern char bootStrapFile[1024];
77extern char tree_file[1024];
78extern char infoFileName[1024];
79extern char resultFileName[1024];
80extern char verboseSplitsFileName[1024];
81
82extern double masterTime;
83
84extern const unsigned int mask32[32];
85
86extern volatile branchInfo      **branchInfos;
87extern volatile int NumberOfThreads;
88extern volatile int NumberOfJobs;
89
90
91static void mre(hashtable *h, boolean icp, entry*** sbi, int* len, int which, int n, unsigned int vectorLength, boolean sortp, tree *tr, boolean bootStopping);
92
93
94entry *initEntry(void)
95{
96  entry *e = (entry*)rax_malloc(sizeof(entry));
97
98  e->bitVector     = (unsigned int*)NULL;
99  e->treeVector    = (unsigned int*)NULL;
100  e->supportVector = (int*)NULL;
101  e->bipNumber  = 0;
102  e->bipNumber2 = 0;
103  e->supportFromTreeset[0] = 0;
104  e->supportFromTreeset[1] = 0;
105  e->next       = (entry*)NULL;
106
107  return e;
108} 
109
110
111
112
113hashtable *initHashTable(hashNumberType n)
114{
115  /*
116     init with primes
117   
118     static const hashNumberType initTable[] = {53, 97, 193, 389, 769, 1543, 3079, 6151, 12289, 24593, 49157, 98317,
119     196613, 393241, 786433, 1572869, 3145739, 6291469, 12582917, 25165843,
120     50331653, 100663319, 201326611, 402653189, 805306457, 1610612741};
121  */
122
123  /* init with powers of two */
124
125  static const  hashNumberType initTable[] = {64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384,
126                                              32768, 65536, 131072, 262144, 524288, 1048576, 2097152,
127                                              4194304, 8388608, 16777216, 33554432, 67108864, 134217728,
128                                              268435456, 536870912, 1073741824, 2147483648U};
129 
130  hashtable *h = (hashtable*)rax_malloc(sizeof(hashtable));
131 
132  hashNumberType
133    tableSize,
134    i,
135    primeTableLength = sizeof(initTable)/sizeof(initTable[0]),
136    maxSize = (hashNumberType)-1;   
137
138  assert(n <= maxSize);
139
140  i = 0;
141
142  while(initTable[i] < n && i < primeTableLength)
143    i++;
144
145  assert(i < primeTableLength);
146
147  tableSize = initTable[i];
148
149  /* printf("Hash table init with size %u\n", tableSize); */
150
151  h->table = (entry**)rax_calloc(tableSize, sizeof(entry*));
152  h->tableSize = tableSize; 
153  h->entryCount = 0; 
154
155  return h;
156}
157
158
159
160
161void freeHashTable(hashtable *h)
162{
163  hashNumberType
164    i,
165    entryCount = 0;
166   
167
168  for(i = 0; i < h->tableSize; i++)
169    {
170      if(h->table[i] != NULL)
171        {
172          entry *e = h->table[i];
173          entry *previous;       
174
175          do
176            {
177              previous = e;
178              e = e->next;
179
180              if(previous->bitVector)
181                rax_free(previous->bitVector);
182
183              if(previous->treeVector)
184                rax_free(previous->treeVector);
185
186              if(previous->supportVector)
187                rax_free(previous->supportVector);
188             
189              rax_free(previous);             
190              entryCount++;
191            }
192          while(e != NULL);       
193        }
194
195    }
196
197  assert(entryCount == h->entryCount);
198 
199  rax_free(h->table);
200}
201
202
203
204void cleanupHashTable(hashtable *h, int state)
205{
206  hashNumberType
207    k,
208    entryCount = 0,
209    removeCount = 0;
210 
211  assert(state == 1 || state == 0);
212
213  for(k = 0, entryCount = 0; k < h->tableSize; k++)         
214    {     
215      if(h->table[k] != NULL)
216        {
217          entry *e = h->table[k];
218          entry *start     = (entry*)NULL;
219          entry *lastValid = (entry*)NULL;
220                 
221          do
222            {                           
223              if(state == 0)
224                {
225                  e->treeVector[0] = e->treeVector[0] & 2;     
226                  assert(!(e->treeVector[0] & 1));
227                }
228              else
229                {
230                  e->treeVector[0] = e->treeVector[0] & 1;
231                  assert(!(e->treeVector[0] & 2));
232                }
233             
234              if(e->treeVector[0] != 0)
235                {
236                  if(!start)
237                    start = e;
238                  lastValid = e;
239                  e = e->next;
240                }         
241              else
242                {
243                  entry *remove = e;
244                  e = e->next;
245                 
246                  removeCount++;
247
248                  if(lastValid)                             
249                    lastValid->next = remove->next;
250
251                  if(remove->bitVector)
252                    rax_free(remove->bitVector);
253                  if(remove->treeVector)
254                    rax_free(remove->treeVector);
255                  if(remove->supportVector)
256                    rax_free(remove->supportVector);
257                  rax_free(remove);             
258                }
259             
260              entryCount++;                 
261            }
262          while(e != NULL);     
263
264          if(!start)
265            {
266              assert(!lastValid);
267              h->table[k] = NULL;
268            }
269          else
270            {
271              h->table[k] = start;
272            }           
273        }   
274    }
275
276  assert(entryCount ==  h->entryCount); 
277
278  h->entryCount -= removeCount;
279}
280
281
282
283
284
285
286
287
288
289
290
291unsigned int **initBitVector(tree *tr, unsigned int *vectorLength)
292{
293  unsigned int **bitVectors = (unsigned int **)rax_malloc(sizeof(unsigned int*) * 2 * tr->mxtips);
294  int i;
295
296  if(tr->mxtips % MASK_LENGTH == 0)
297    *vectorLength = tr->mxtips / MASK_LENGTH;
298  else
299    *vectorLength = 1 + (tr->mxtips / MASK_LENGTH); 
300 
301  for(i = 1; i <= tr->mxtips; i++)
302    {
303      bitVectors[i] = (unsigned int *)rax_calloc(*vectorLength, sizeof(unsigned int));
304      bitVectors[i][(i - 1) / MASK_LENGTH] |= mask32[(i - 1) % MASK_LENGTH];
305    }
306 
307  for(i = tr->mxtips + 1; i < 2 * tr->mxtips; i++) 
308    bitVectors[i] = (unsigned int *)rax_malloc(sizeof(unsigned int) * *vectorLength);
309
310  return bitVectors;
311}
312
313void freeBitVectors(unsigned int **v, int n)
314{
315  int i;
316
317  for(i = 1; i < n; i++)
318    rax_free(v[i]);
319}
320
321
322
323
324
325static void newviewBipartitions(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength)
326{
327  if(isTip(p->number, numsp))
328    return;
329  {
330    nodeptr
331      q = p->next->back, 
332      r = p->next->next->back;
333    unsigned int       
334      *vector = bitVectors[p->number],
335      *left  = bitVectors[q->number],
336      *right = bitVectors[r->number];
337    unsigned 
338      int i;           
339
340    while(!p->x)
341      { 
342        if(!p->x)
343          getxnode(p);
344      }
345
346    p->hash = q->hash ^ r->hash;
347
348    if(isTip(q->number, numsp) && isTip(r->number, numsp))
349      {         
350        for(i = 0; i < vectorLength; i++)
351          vector[i] = left[i] | right[i];               
352      }
353    else
354      { 
355        if(isTip(q->number, numsp) || isTip(r->number, numsp))
356          {
357            if(isTip(r->number, numsp))
358              { 
359                nodeptr tmp = r;
360                r = q;
361                q = tmp;
362              }   
363                   
364            while(!r->x)
365              {
366                if(!r->x)
367                  newviewBipartitions(bitVectors, r, numsp, vectorLength);
368              }   
369
370            for(i = 0; i < vectorLength; i++)
371              vector[i] = left[i] | right[i];           
372          }
373        else
374          {         
375            while((!r->x) || (!q->x))
376              {
377                if(!q->x)
378                  newviewBipartitions(bitVectors, q, numsp, vectorLength);
379                if(!r->x)
380                  newviewBipartitions(bitVectors, r, numsp, vectorLength);
381              }                                   
382
383            for(i = 0; i < vectorLength; i++)
384              vector[i] = left[i] | right[i];   
385          }
386
387      }     
388  }     
389}
390
391
392/* compute bit-vectors representing bipartitions/splits for a multi-furcating tree */
393
394static void newviewBipartitionsMultifurcating(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength)
395{
396  if(isTip(p->number, numsp))
397    return;
398  {
399    nodeptr
400      q,
401      firstDescendant;
402   
403    unsigned int       
404      *vector = bitVectors[p->number];
405   
406    unsigned 
407      int i;           
408   
409    int
410      x_set = 0,
411      number;
412   
413    /* Set the directional x token to the correct element in the
414       cyclically linked list representing an inner node */
415
416    q = p->next;
417   
418    if(p->x)
419      x_set++;
420
421    p->x = 1;   
422    number = p->number;
423
424    while(q != p)
425      {
426        if(q->x) 
427          x_set++;
428        q->x = 0;
429        assert(q->number == number);
430        q = q->next;
431      }
432   
433    assert(x_set == 1);
434
435    /* get the first connecting branch of the node to initialize the
436       bipartition hash value and the bit vector representing this split.
437     */
438
439    firstDescendant = p->next->back;     
440   
441    /* if this is not a tip, we first need to recursively
442       compute the hash values and bit-vectors for the subtree rooted at
443       firstDescendant */
444
445    if(!isTip(firstDescendant->number, numsp))
446      { 
447        if(!firstDescendant->x)
448          newviewBipartitionsMultifurcating(bitVectors, firstDescendant, numsp, vectorLength);
449      }
450       
451    /* initialize the bit vector of the current split by the bit vector of the first descandant */
452   
453    for(i = 0; i < vectorLength; i++)
454      vector[i] = bitVectors[firstDescendant->number][i];
455       
456    /* initialize the hash key of the current split by the hash key of the first descandant */
457
458    p->hash = firstDescendant->hash;
459   
460    /* handle all other descendants of this inner node potentially representing a multi-furcation */
461
462    q = p->next->next;
463
464    while(q != p)
465      {
466        /* update the has key by xoring with the current hash with the hash of this descendant */
467       
468        p->hash = p->hash ^ q->back->hash;
469       
470        if(!isTip(q->back->number, numsp))
471          {       
472            if(!q->back->x)
473              newviewBipartitionsMultifurcating(bitVectors, q->back, numsp, vectorLength);
474          }
475               
476        /* update the bit-vector representing the current split by applying a bitwise with the bit vector of the descendant */
477       
478        for(i = 0; i < vectorLength; i++)
479          vector[i] = bitVectors[q->back->number][i] | vector[i];               
480       
481        q = q->next;
482      }     
483  }     
484}
485
486
487static void insertHash(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, int bipNumber, hashNumberType position)
488{
489  entry *e = initEntry();
490
491  e->bipNumber = bipNumber; 
492  /*e->bitVector = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int)); */
493
494  e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
495  memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
496 
497  memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
498 
499  if(h->table[position] != NULL)
500    {
501      e->next = h->table[position];
502      h->table[position] = e;           
503    }
504  else
505    h->table[position] = e;
506
507  h->entryCount =  h->entryCount + 1;
508}
509
510
511
512static int countHash(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, hashNumberType position)
513{ 
514  if(h->table[position] == NULL)         
515    return -1;
516  {
517    entry *e = h->table[position];     
518
519    do
520      { 
521        unsigned int i;
522
523        for(i = 0; i < vectorLength; i++)
524          if(bitVector[i] != e->bitVector[i])
525            goto NEXT;
526           
527        return (e->bipNumber);   
528      NEXT:
529        e = e->next;
530      }
531    while(e != (entry*)NULL); 
532     
533    return -1;   
534  }
535
536}
537
538static void insertHashAll(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, int treeNumber,  hashNumberType position)
539{   
540  if(h->table[position] != NULL)
541    {
542      entry *e = h->table[position];     
543
544      do
545        {       
546          unsigned int i;
547         
548          for(i = 0; i < vectorLength; i++)
549            if(bitVector[i] != e->bitVector[i])
550              break;
551         
552          if(i == vectorLength)
553            {
554              if(treeNumber == 0)
555                e->bipNumber =  e->bipNumber  + 1;
556              else
557                e->bipNumber2 = e->bipNumber2 + 1;
558              return;
559            }
560         
561          e = e->next;   
562        }
563      while(e != (entry*)NULL); 
564
565      e = initEntry(); 
566 
567      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int)); */
568      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
569      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
570
571
572      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
573
574      if(treeNumber == 0)       
575        e->bipNumber  = 1;             
576      else               
577        e->bipNumber2 = 1;
578       
579      e->next = h->table[position];
580      h->table[position] = e;             
581    }
582  else
583    {
584      entry *e = initEntry(); 
585 
586      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int)); */
587
588      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
589      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
590
591      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
592
593      if(treeNumber == 0)       
594        e->bipNumber  = 1;             
595      else   
596        e->bipNumber2 = 1;     
597
598      h->table[position] = e;
599    }
600
601  h->entryCount =  h->entryCount + 1;
602}
603
604
605
606static void insertHashBootstop(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, int treeNumber, int treeVectorLength, hashNumberType position)
607{   
608  if(h->table[position] != NULL)
609    {
610      entry *e = h->table[position];     
611
612      do
613        {       
614          unsigned int i;
615         
616          for(i = 0; i < vectorLength; i++)
617            if(bitVector[i] != e->bitVector[i])
618              break;
619         
620          if(i == vectorLength)
621            {
622              e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
623              return;
624            }
625         
626          e = e->next;
627        }
628      while(e != (entry*)NULL); 
629
630      e = initEntry(); 
631
632      e->bipNumber = h->entryCount;
633       
634      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int));*/
635      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
636      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
637
638
639      e->treeVector = (unsigned int*)rax_calloc(treeVectorLength, sizeof(unsigned int));
640     
641      e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
642      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
643     
644      e->next = h->table[position];
645      h->table[position] = e;         
646    }
647  else
648    {
649      entry *e = initEntry(); 
650
651      e->bipNumber = h->entryCount;
652
653      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int));*/
654
655      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
656      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
657
658      e->treeVector = (unsigned int*)rax_calloc(treeVectorLength, sizeof(unsigned int));
659
660      e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
661      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);     
662
663      h->table[position] = e;
664    }
665
666  h->entryCount =  h->entryCount + 1;
667}
668
669static void insertHashRF(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, int treeNumber, int treeVectorLength, hashNumberType position, int support, 
670                         boolean computeWRF)
671{     
672  if(h->table[position] != NULL)
673    {
674      entry *e = h->table[position];     
675
676      do
677        {       
678          unsigned int i;
679         
680          for(i = 0; i < vectorLength; i++)
681            if(bitVector[i] != e->bitVector[i])
682              break;
683         
684          if(i == vectorLength)
685            {
686              e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
687              if(computeWRF)
688                {
689                  e->supportVector[treeNumber] = support;
690                 
691                  assert(0 <= treeNumber && treeNumber < treeVectorLength * MASK_LENGTH);
692                }
693              return;
694            }
695         
696          e = e->next;
697        }
698      while(e != (entry*)NULL); 
699
700      e = initEntry(); 
701       
702      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int));*/
703      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
704      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
705
706
707      e->treeVector = (unsigned int*)rax_calloc(treeVectorLength, sizeof(unsigned int));
708      if(computeWRF)
709        e->supportVector = (int*)rax_calloc(treeVectorLength * MASK_LENGTH, sizeof(int));
710
711      e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
712      if(computeWRF)
713        {
714          e->supportVector[treeNumber] = support;
715         
716          assert(0 <= treeNumber && treeNumber < treeVectorLength * MASK_LENGTH);
717        }
718
719      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
720     
721      e->next = h->table[position];
722      h->table[position] = e;         
723    }
724  else
725    {
726      entry *e = initEntry(); 
727       
728      /*e->bitVector  = (unsigned int*)rax_calloc(vectorLength, sizeof(unsigned int)); */
729
730      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
731      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
732
733      e->treeVector = (unsigned int*)rax_calloc(treeVectorLength, sizeof(unsigned int));
734      if(computeWRF)   
735        e->supportVector = (int*)rax_calloc(treeVectorLength * MASK_LENGTH, sizeof(int));
736
737
738      e->treeVector[treeNumber / MASK_LENGTH] |= mask32[treeNumber % MASK_LENGTH];
739      if(computeWRF)
740        {
741          e->supportVector[treeNumber] = support;
742         
743          assert(0 <= treeNumber && treeNumber < treeVectorLength * MASK_LENGTH);
744        }
745
746      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);     
747
748      h->table[position] = e;
749    }
750
751  h->entryCount =  h->entryCount + 1;
752}
753
754
755
756void bitVectorInitravSpecial(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength, hashtable *h, int treeNumber, int function, branchInfo *bInf, 
757                             int *countBranches, int treeVectorLength, boolean traverseOnly, boolean computeWRF)
758{
759  if(isTip(p->number, numsp))
760    return;
761  else
762    {
763      nodeptr q = p->next;         
764
765      do 
766        {
767          bitVectorInitravSpecial(bitVectors, q->back, numsp, vectorLength, h, treeNumber, function, bInf, countBranches, treeVectorLength, traverseOnly, computeWRF);
768          q = q->next;
769        }
770      while(q != p);
771           
772      newviewBipartitions(bitVectors, p, numsp, vectorLength);
773     
774      assert(p->x);
775
776      if(traverseOnly)
777        {
778          if(!(isTip(p->back->number, numsp)))
779            *countBranches =  *countBranches + 1;
780          return;
781        }
782
783      if(!(isTip(p->back->number, numsp)))
784        {
785          unsigned int *toInsert  = bitVectors[p->number];
786          hashNumberType position = p->hash % h->tableSize;
787         
788          assert(!(toInsert[0] & 1));   
789
790          switch(function)
791            {
792            case BIPARTITIONS_ALL:           
793              insertHashAll(toInsert, h, vectorLength, treeNumber, position);
794              *countBranches =  *countBranches + 1;     
795              break;
796            case GET_BIPARTITIONS_BEST:             
797              insertHash(toInsert, h, vectorLength, *countBranches, position);       
798             
799              p->bInf            = &bInf[*countBranches];
800              p->back->bInf      = &bInf[*countBranches];       
801              p->bInf->support   = 0;           
802              p->bInf->oP = p;
803              p->bInf->oQ = p->back;
804             
805              *countBranches =  *countBranches + 1;             
806              break;
807            case DRAW_BIPARTITIONS_BEST:             
808              {
809                int found = countHash(toInsert, h, vectorLength, position);
810                if(found >= 0)
811                  bInf[found].support =  bInf[found].support + 1;
812                *countBranches =  *countBranches + 1;
813              }       
814              break;
815            case BIPARTITIONS_BOOTSTOP:       
816              insertHashBootstop(toInsert, h, vectorLength, treeNumber, treeVectorLength, position);
817              *countBranches =  *countBranches + 1;
818              break;
819            case BIPARTITIONS_RF:
820              if(computeWRF)
821                assert(p->support == p->back->support);
822              insertHashRF(toInsert, h, vectorLength, treeNumber, treeVectorLength, position, p->support, computeWRF);
823              *countBranches =  *countBranches + 1;
824              break;       
825            default:
826              assert(0);
827            }             
828        }
829     
830    }
831}
832
833static void linkBipartitions(nodeptr p, tree *tr, branchInfo *bInf, int *counter, int numberOfTrees)
834{
835  if(isTip(p->number, tr->mxtips))   
836    {
837      assert(p->bInf == (branchInfo*) NULL && p->back->bInf == (branchInfo*) NULL);     
838      return;
839    }
840  else
841    {
842      nodeptr q;         
843     
844      q = p->next;
845
846      while(q != p)
847        {
848          linkBipartitions(q->back, tr, bInf, counter, numberOfTrees); 
849          q = q->next;
850        }
851     
852      if(!(isTip(p->back->number, tr->mxtips)))
853        {
854          double support;
855
856          p->bInf       = &bInf[*counter];
857          p->back->bInf = &bInf[*counter]; 
858
859          support = ((double)(p->bInf->support)) / ((double) (numberOfTrees));
860          p->bInf->support = (int)(0.5 + support * 100.0);                       
861
862          assert(p->bInf->oP == p);
863          assert(p->bInf->oQ == p->back);
864         
865          *counter = *counter + 1;
866        }
867
868
869      return;
870    }
871}
872
873
874static int readSingleTree(tree *tr, char *fileName, analdef *adef, boolean readBranches, boolean readNodeLabels, boolean completeTree)
875{ 
876  FILE 
877    *f = myfopen(fileName, "r");
878
879  int 
880    numberOfTaxa,
881    ch,
882    trees = 0;
883
884  while((ch = fgetc(f)) != EOF)
885    if(ch == ';')
886      trees++;
887   
888  assert(trees == 1);
889
890  printBothOpen("\n\nFound 1 tree in File %s\n\n", fileName);
891
892  rewind(f);
893
894  treeReadLen(f, tr, readBranches, readNodeLabels, TRUE, adef, completeTree, FALSE);
895 
896  numberOfTaxa = tr->ntips;
897
898  fclose(f);
899
900  return numberOfTaxa;
901}
902
903/*************** function for drawing bipartitions on a bifurcating tree ***********/
904
905void calcBipartitions(tree *tr, analdef *adef, char *bestTreeFileName, char *bootStrapFileName)
906{ 
907  branchInfo
908    *bInf;
909  unsigned int vLength;
910
911  int 
912    numberOfTaxa = 0,
913    branchCounter = 0,
914    counter = 0, 
915    numberOfTrees = 0,
916    i;
917
918  unsigned int 
919    **bitVectors = initBitVector(tr, &vLength);
920 
921  hashtable *h = 
922    initHashTable(tr->mxtips * 10);
923
924  FILE 
925    *treeFile; 
926 
927  numberOfTaxa = readSingleTree(tr, bestTreeFileName, adef, FALSE, FALSE, TRUE);   
928 
929  bInf = (branchInfo*)rax_malloc(sizeof(branchInfo) * (tr->mxtips - 3));
930
931  bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, GET_BIPARTITIONS_BEST, bInf, &branchCounter, 0, FALSE, FALSE);   
932
933  if(numberOfTaxa != tr->mxtips)
934    {
935      printBothOpen("The number of taxa in the reference tree file \"%s\" is %d and\n",  bestTreeFileName, numberOfTaxa);
936      printBothOpen("is not equal to the number of taxa in the bootstrap tree file \"%s\" which is %d.\n", bootStrapFileName, tr->mxtips);
937      printBothOpen("RAxML will exit now with an error ....\n\n");
938    }
939 
940  assert((int)h->entryCount == (tr->mxtips - 3)); 
941  assert(branchCounter == (tr->mxtips - 3));
942 
943  treeFile = getNumberOfTrees(tr, bootStrapFileName, adef);
944
945  numberOfTrees = tr->numberOfTrees;
946 
947  for(i = 0; i < numberOfTrees; i++)
948    {               
949      int 
950        bCount = 0;
951     
952      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);
953      assert(tr->ntips == tr->mxtips);
954     
955      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, DRAW_BIPARTITIONS_BEST, bInf, &bCount, 0, FALSE, FALSE);     
956      assert(bCount == tr->mxtips - 3);     
957    }     
958 
959  fclose(treeFile);
960   
961  readSingleTree(tr, bestTreeFileName, adef, TRUE, FALSE, TRUE); 
962   
963  linkBipartitions(tr->nodep[1]->back, tr, bInf, &counter, numberOfTrees);
964
965  assert(counter == branchCounter);
966
967  printBipartitionResult(tr, adef, TRUE, FALSE);   
968
969  freeBitVectors(bitVectors, 2 * tr->mxtips);
970  rax_free(bitVectors);
971  freeHashTable(h);
972  rax_free(h); 
973
974  rax_free(bInf);
975}
976
977
978/****** functions for IC computation ***********************************************/
979
980static void insertHash_IC(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, hashNumberType position)
981{   
982  if(h->table[position] != NULL)
983    {
984      entry
985        *e = h->table[position];     
986
987      do
988        {       
989          unsigned int 
990            i;
991         
992          for(i = 0; i < vectorLength; i++)
993            if(bitVector[i] != e->bitVector[i])
994              break;
995         
996          if(i == vectorLength)
997            {       
998              e->bipNumber =    e->bipNumber  + 1;           
999              return;
1000            }
1001         
1002          e = e->next;   
1003        }
1004      while(e != (entry*)NULL); 
1005
1006      e = initEntry(); 
1007       
1008      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
1009      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
1010      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
1011   
1012      e->bipNumber  = 1;       
1013       
1014      e->next = h->table[position];
1015      h->table[position] = e;             
1016    }
1017  else
1018    {
1019      entry
1020        *e = initEntry(); 
1021 
1022      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));
1023      memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
1024      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
1025
1026     
1027      e->bipNumber  = 1;                   
1028
1029      h->table[position] = e;
1030    }
1031
1032  h->entryCount =  h->entryCount + 1;
1033}
1034
1035
1036
1037static unsigned int findHash_IC(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, hashNumberType position)
1038{
1039  if(h->table[position] == NULL)         
1040    return 0;
1041  {
1042    entry *e = h->table[position];     
1043
1044    do
1045      { 
1046        unsigned int 
1047          i;
1048
1049        for(i = 0; i < vectorLength; i++)
1050          if(bitVector[i] != e->bitVector[i])
1051            goto NEXT;
1052           
1053        return e->bipNumber;     
1054      NEXT:
1055        e = e->next;
1056      }
1057    while(e != (entry*)NULL); 
1058     
1059    return 0;   
1060  }
1061}
1062
1063static boolean compatibleIC(unsigned int *A, unsigned int *C, unsigned int bvlen)
1064{
1065  unsigned int 
1066    i;
1067 
1068  for(i = 0; i < bvlen; i++)       
1069    if(A[i] & C[i])
1070      break;
1071         
1072  if(i == bvlen)
1073    return TRUE;
1074 
1075  for(i = 0; i < bvlen; i++)
1076    if(A[i] & ~C[i])
1077      break;
1078   
1079  if(i == bvlen) 
1080    return TRUE; 
1081 
1082  /*
1083    not required -> ask Andre
1084  for(i = 0; i < bvlen; i++)
1085    if(~A[i] & ~C[i])
1086      break;
1087   
1088  if(i == bvlen) 
1089  return TRUE; */
1090
1091 
1092  for(i = 0; i < bvlen; i++)
1093    if(~A[i] & C[i])
1094      break;
1095 
1096  if(i == bvlen)
1097    return TRUE; 
1098  else
1099    return FALSE;
1100}
1101
1102static int sortByBipNumber(const void *a, const void *b)
1103{       
1104  int       
1105    ca = ((*((entry **)a))->bipNumber),
1106    cb = ((*((entry **)b))->bipNumber);
1107 
1108  if (ca == cb) 
1109    return 0;
1110 
1111  return ((ca < cb)?1:-1);
1112}
1113
1114static int sortByTreeSet(const void *a, const void *b)
1115{       
1116  int       
1117    ca = ((*((entry **)a))->supportFromTreeset)[0],
1118    cb = ((*((entry **)b))->supportFromTreeset)[0];
1119 
1120  if (ca == cb) 
1121    return 0;
1122 
1123  return ((ca < cb)?1:-1);
1124}
1125
1126
1127static unsigned int countIncompatibleBipartitions(unsigned int *toInsert, hashtable *h,  hashNumberType vectorLength, unsigned int *maxima, unsigned int *maxCounter, boolean bipNumber, 
1128                                                  unsigned int numberOfTrees, unsigned int **maximaBitVectors)
1129{
1130  unsigned int   
1131    threshold = numberOfTrees / 20,
1132    max = 0,
1133    entryVectorSize = h->entryCount,
1134    entryVectorElements = 0,
1135    k,
1136    uncompatible = 0;
1137
1138  entry
1139    **entryVector = (entry**)rax_malloc(sizeof(entry*) * entryVectorSize);
1140
1141  for(k = 0; k < entryVectorSize; k++)
1142    entryVector[k] = (entry*)NULL;
1143
1144  for(k = 0; k < h->tableSize; k++)         
1145    {   
1146      if(h->table[k] != NULL)
1147        {
1148          entry
1149            *e = h->table[k];             
1150
1151          do
1152            {                     
1153              if(!compatibleIC(toInsert, e->bitVector, vectorLength))
1154                {
1155                  unsigned int
1156                    support;
1157
1158                  if(bipNumber)
1159                    support = e->bipNumber;
1160                  else
1161                    support = e->supportFromTreeset[0];         
1162
1163                  if(support > max)
1164                    max = support;                             
1165                 
1166                  uncompatible += support;
1167
1168                  entryVector[entryVectorElements] = e;
1169                  entryVectorElements++;
1170                  assert(entryVectorElements < entryVectorSize);
1171                }
1172                 
1173              e = e->next;
1174            }
1175          while(e != NULL);
1176        }
1177    }
1178
1179  if(bipNumber)
1180    {
1181      qsort(entryVector, entryVectorElements, sizeof(entry *), sortByBipNumber);
1182      assert(max == entryVector[0]->bipNumber);
1183    }
1184  else
1185    {
1186      qsort(entryVector, entryVectorElements, sizeof(entry *), sortByTreeSet);
1187      assert(max == entryVector[0]->supportFromTreeset[0]);
1188    }
1189
1190  for(k = 0; k < entryVectorElements; k++)
1191    {     
1192      unsigned int 
1193        j,
1194        support;
1195
1196      boolean
1197        uncompat = TRUE;
1198
1199      entry
1200        *referenceEntry = entryVector[k];           
1201     
1202      if(bipNumber)
1203        support = entryVector[k]->bipNumber;
1204      else
1205        support = entryVector[k]->supportFromTreeset[0];
1206     
1207      if(k > 0)
1208        {
1209          if(support > threshold)
1210            {
1211              for(j = 0; j < k; j++)
1212                {
1213                  entry
1214                    *checkEntry = entryVector[j];
1215                 
1216                  if(compatibleIC(checkEntry->bitVector, referenceEntry->bitVector, vectorLength)) 
1217                    {
1218                      uncompat = FALSE;
1219                      break;
1220                    }
1221                }
1222            }
1223          else
1224            uncompat = FALSE;
1225        }
1226
1227      if(uncompat)
1228        {
1229          maximaBitVectors[*maxCounter] = entryVector[k]->bitVector;           
1230          maxima[*maxCounter] = support;                         
1231          *maxCounter = *maxCounter + 1;         
1232        }
1233    }
1234
1235  rax_free(entryVector);
1236
1237  return uncompatible;
1238}
1239
1240static double computeIC_Value(unsigned int supportedBips, unsigned int *maxima, unsigned int numberOfTrees, unsigned int maxCounter, boolean computeIC_All, boolean warnNegativeIC)
1241{
1242  unsigned int 
1243    loopLength,
1244    i,
1245    totalBipsAll = supportedBips;
1246
1247  double 
1248    ic,
1249    n = 1 + (double)maxCounter,
1250    supportFreq;
1251 
1252  boolean
1253    negativeIC = FALSE;   
1254 
1255  if(computeIC_All)
1256    {
1257      loopLength = maxCounter;
1258      n = 1 + (double)maxCounter;
1259    }
1260  else
1261    {
1262      loopLength = 1;
1263      n = 2.0;
1264    }
1265
1266  // should never enter this function when the bip is supported by 100%
1267
1268  assert(supportedBips < numberOfTrees);
1269
1270  // must be larger than 0 in this case
1271
1272  //assert(maxCounter > 0);
1273
1274  // figure out if the competing bipartition is higher support
1275  // can happen for MRE and when drawing IC values on best ML tree
1276
1277  if(maxima[0] > supportedBips)
1278    {
1279      negativeIC = TRUE;
1280     
1281      if(warnNegativeIC)
1282        {
1283          printBothOpen("\nMax conflicting bipartition frequency: %d is larger than frequency of the included bipartition: %d\n", maxima[0], supportedBips);
1284          printBothOpen("This is interesting, but not unexpected when computing extended Majority Rule consensus trees.\n");
1285          printBothOpen("Please send an email with the input files and command line\n");
1286          printBothOpen("to Alexandros.Stamatakis@gmail.com.\n");
1287          printBothOpen("Thank you :-)\n\n");   
1288        }
1289    }
1290
1291  for(i = 0; i < loopLength; i++)
1292    totalBipsAll += maxima[i]; 
1293
1294  //neither support for the actual bipartition, nor for the conflicting ones
1295  //I am not sure that this will happen, but anyway
1296  if(totalBipsAll == 0)
1297    return 0.0;
1298 
1299  supportFreq = (double)supportedBips / (double)totalBipsAll;
1300 
1301  if(supportedBips == 0)
1302    ic = log(n);
1303  else   
1304    ic = log(n) + supportFreq * log(supportFreq);
1305
1306  for(i = 0; i < loopLength; i++)
1307    {
1308      assert(maxima[i] > 0);
1309
1310      supportFreq =  (double)maxima[i] / (double)totalBipsAll;
1311     
1312      if(maxima[i] != 0)
1313        ic += supportFreq * log(supportFreq);
1314    }
1315 
1316  ic /= log(n);
1317 
1318  if(negativeIC)
1319    return (-ic);
1320  else
1321    return ic;
1322}
1323
1324static void printSplit(FILE *f, FILE *v, unsigned int *bitVector, tree *tr, double support, double ic, unsigned int frequency)
1325{
1326  int 
1327    i,
1328    countLeftTaxa = 0,
1329    countRightTaxa = 0,
1330    taxa = 0,
1331    totalTaxa = 0;
1332
1333  for(i = 0; i < tr->mxtips; i++)       
1334    if((bitVector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]))         
1335      countLeftTaxa++;
1336
1337  countRightTaxa = tr->mxtips - countLeftTaxa;
1338
1339  fprintf(f, "((");
1340
1341  for(i = 0; i < tr->mxtips; i++)   
1342    {
1343      if((bitVector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]))       
1344        {
1345          fprintf(v, "*");
1346          fprintf(f, "%s", tr->nameList[i+1]);           
1347          taxa++;
1348          totalTaxa++;
1349          if(taxa < countLeftTaxa)
1350            fprintf(f, ", ");
1351        }   
1352      else
1353        fprintf(v, "-");
1354    } 
1355
1356  fprintf(v, "\t%u/%f/%f\n", frequency, support * 100.0, ic);
1357
1358  fprintf(f, "),(");
1359
1360  taxa = 0;
1361
1362  for(i = 0; i < tr->mxtips; i++)   
1363    {
1364      if(!(bitVector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]))       
1365        {
1366          fprintf(f, "%s", tr->nameList[i+1]);   
1367          taxa++;
1368          totalTaxa++;
1369          if(taxa < countRightTaxa)
1370            fprintf(f, ", ");
1371        }   
1372    }
1373
1374  assert(totalTaxa == tr->mxtips);
1375
1376  fprintf(f, "));\n");
1377}
1378
1379static void printFullySupportedSplit(tree *tr, unsigned int *bitVector, unsigned int numberOfTrees)
1380{
1381  FILE   
1382    *v = myfopen(verboseSplitsFileName, "a");
1383
1384  int 
1385    i;
1386
1387  fprintf(v, "partition: \n");
1388
1389  for(i = 0; i < tr->mxtips; i++)   
1390    {
1391      if((bitVector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]))               
1392        fprintf(v, "*");               
1393      else
1394        fprintf(v, "-");
1395    } 
1396
1397  fprintf(v, "\t%u/%f/%f\n\n\n", numberOfTrees, 100.0, 1.0);
1398
1399  fclose(v);
1400}
1401
1402static void printVerboseTaxonNames(tree *tr)
1403{
1404  FILE 
1405    *f = myfopen(verboseSplitsFileName, "w");
1406
1407  int 
1408    i;
1409
1410  fprintf(f, "\n");
1411
1412  for(i = 1; i <= tr->mxtips; i++)
1413    fprintf(f, "%s \n", tr->nameList[i]);
1414
1415  fprintf(f, "\n");
1416
1417  fclose(f);
1418 
1419}
1420
1421static void printVerboseIC(tree *tr, unsigned int supportedBips, unsigned int *toInsert, unsigned int maxCounter, unsigned int *maxima, 
1422                           unsigned int **maximaBitVectors, unsigned int numberOfTrees, int counter, double ic)
1423{
1424  unsigned int 
1425    i;
1426
1427  double 
1428    support = (double)supportedBips / (double)numberOfTrees;
1429
1430  FILE 
1431    *f,
1432    *v = myfopen(verboseSplitsFileName, "a");
1433
1434  char 
1435    fileName[1024],
1436    id[64];
1437
1438  sprintf(id, "%d", counter);
1439  strcpy(fileName, workdir);
1440  strcat(fileName, "RAxML_verboseIC.");
1441  strcat(fileName, run_id);
1442  strcat(fileName, ".");
1443  strcat(fileName, id);
1444
1445  f = myfopen(fileName, "w");
1446
1447  printBothOpen("Support for split number %d in tree: %f\n", counter, support);
1448
1449  fprintf(v, "partition: \n");
1450
1451  printSplit(f, v, toInsert, tr, support, ic, supportedBips); 
1452
1453  for(i = 0; i < maxCounter; i++)
1454    {
1455      support = (double)maxima[i] / (double)numberOfTrees;
1456
1457      printBothOpen("Support for conflicting split number %u: %f\n", i, support);
1458
1459      printSplit(f, v, maximaBitVectors[i], tr, support, ic, maxima[i]);       
1460    }
1461 
1462  printBothOpen("All Newick-formatted splits for this bipartition have been written to file %s\n", fileName);
1463  printBothOpen("\n\n");
1464
1465  fprintf(v, "\n\n");
1466
1467  fclose(f);
1468  fclose(v);
1469}
1470
1471
1472
1473
1474static void bitVectorInitravIC(tree *tr, unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength, hashtable *h, int treeNumber, int function, branchInfo *bInf, 
1475                               int *countBranches, int treeVectorLength, boolean traverseOnly, boolean computeWRF, double *tc, double *tcAll, boolean verboseIC)
1476{
1477  if(isTip(p->number, numsp))
1478    return;
1479  else
1480    {
1481      nodeptr q = p->next;         
1482
1483      do 
1484        {
1485          bitVectorInitravIC(tr, bitVectors, q->back, numsp, vectorLength, h, treeNumber, function, bInf, countBranches, treeVectorLength, traverseOnly, computeWRF, tc, tcAll, verboseIC);
1486          q = q->next;
1487        }
1488      while(q != p);
1489           
1490      newviewBipartitions(bitVectors, p, numsp, vectorLength);
1491     
1492      assert(p->x);
1493
1494      if(traverseOnly)
1495        {
1496          if(!(isTip(p->back->number, numsp)))
1497            *countBranches =  *countBranches + 1;
1498          return;
1499        }
1500
1501      if(!(isTip(p->back->number, numsp)))
1502        {
1503          unsigned int *toInsert  = bitVectors[p->number];
1504          hashNumberType position = p->hash % h->tableSize;
1505         
1506          assert(!(toInsert[0] & 1));   
1507
1508          switch(function)
1509            {     
1510            case GATHER_BIPARTITIONS_IC:
1511              insertHash_IC(toInsert, h, vectorLength, position);
1512              *countBranches =  *countBranches + 1;
1513              break;
1514            case FIND_BIPARTITIONS_IC:
1515              {
1516                unsigned int
1517                  maxCounter = 0,
1518                  *maxima = (unsigned int *)rax_calloc(h->entryCount, sizeof(unsigned int)),
1519                  **maximaBitVectors = (unsigned int **)rax_calloc(h->entryCount, sizeof(unsigned int *)),
1520                  numberOfTrees = (unsigned int)treeVectorLength,                 
1521                  supportedBips;
1522               
1523                double           
1524                  ic,
1525                  icAll;               
1526
1527                //obtain the support for the bipartitions in the tree
1528                supportedBips = findHash_IC(toInsert, h, vectorLength, position);
1529
1530                if(supportedBips == numberOfTrees)
1531                  {
1532                    ic = 1.0;
1533                    icAll = 1.0;
1534
1535                    if(verboseIC)
1536                      printFullySupportedSplit(tr, toInsert, numberOfTrees);
1537                  }
1538                else
1539                  {
1540                    //find all incompatible bipartitions in the hash table and also
1541                    //get the conflicting bipartition with maximum support
1542                    countIncompatibleBipartitions(toInsert, h, vectorLength, maxima, &maxCounter, TRUE, numberOfTrees, maximaBitVectors);
1543                               
1544                    //this number must be smaller or equal to the total number of trees           
1545
1546                    assert(supportedBips + maxima[0] <= numberOfTrees);     
1547
1548                    //now compute the IC score for this bipartition
1549                   
1550                    ic    = computeIC_Value(supportedBips, maxima, numberOfTrees, maxCounter, FALSE, FALSE);
1551                    icAll = computeIC_Value(supportedBips, maxima, numberOfTrees, maxCounter, TRUE, FALSE);     
1552                   
1553                    if(verboseIC)                     
1554                      printVerboseIC(tr, supportedBips, toInsert, maxCounter, maxima, maximaBitVectors, numberOfTrees, *countBranches, ic);
1555                             
1556                  }
1557
1558                //printf("%d %d %d %d IC %f\n", supportedBips, unsupportedBips, max, treeVectorLength, ic);
1559               
1560                p->bInf            = &bInf[*countBranches];
1561                p->back->bInf      = &bInf[*countBranches];                 
1562                p->bInf->oP = p;
1563                p->bInf->oQ = p->back;
1564
1565                p->bInf->ic    = ic;
1566                p->bInf->icAll = icAll;
1567
1568                //increment tc
1569                *tc    += ic;
1570                *tcAll += icAll;
1571
1572       
1573               
1574                rax_free(maxima);
1575                rax_free(maximaBitVectors);
1576              }
1577              *countBranches =  *countBranches + 1;
1578              break;
1579            default:
1580              assert(0);
1581            }             
1582        }
1583     
1584    }
1585}
1586
1587
1588
1589
1590void calcBipartitions_IC(tree *tr, analdef *adef, char *bestTreeFileName, char *bootStrapFileName)
1591{ 
1592  branchInfo
1593    *bInf;
1594 
1595  unsigned int 
1596    vLength,
1597    **bitVectors = initBitVector(tr, &vLength);
1598 
1599  int 
1600    bCount = 0,
1601    numberOfTaxa = 0,
1602    numberOfTrees = 0,
1603    i;
1604     
1605  hashtable *h = 
1606    initHashTable(tr->mxtips * 10);
1607
1608  FILE 
1609    *treeFile; 
1610 
1611  double
1612    rtc = 0.0,
1613    rtcAll = 0.0,
1614    tc = 0.0,
1615    tcAll = 0.0;
1616
1617  numberOfTaxa = readSingleTree(tr, bestTreeFileName, adef, FALSE, FALSE, TRUE);   
1618 
1619  bInf = (branchInfo*)rax_malloc(sizeof(branchInfo) * (tr->mxtips - 3));
1620
1621  if(adef->verboseIC)
1622    printVerboseTaxonNames(tr);
1623
1624  if(numberOfTaxa != tr->mxtips)
1625    {
1626      printBothOpen("The number of taxa in the reference tree file \"%s\" is %d and\n",  bestTreeFileName, numberOfTaxa);
1627      printBothOpen("is not equal to the number of taxa in the bootstrap tree file \"%s\" which is %d.\n", bootStrapFileName, tr->mxtips);
1628      printBothOpen("RAxML will exit now with an error ....\n\n");
1629    } 
1630 
1631  treeFile = getNumberOfTrees(tr, bootStrapFileName, adef);
1632
1633  numberOfTrees = tr->numberOfTrees;
1634 
1635  for(i = 0; i < numberOfTrees; i++)
1636    {                     
1637      bCount = 0;
1638     
1639      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);
1640      assert(tr->ntips == tr->mxtips);
1641     
1642      bitVectorInitravIC(tr, bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, GATHER_BIPARTITIONS_IC, (branchInfo*)NULL, &bCount, 0, FALSE, FALSE, &tc, &tcAll, FALSE);     
1643      assert(bCount == tr->mxtips - 3);     
1644    }     
1645 
1646  fclose(treeFile);
1647   
1648  readSingleTree(tr, bestTreeFileName, adef, TRUE, FALSE, TRUE); 
1649   
1650  bCount = 0;
1651  tc = 0.0;
1652
1653  bitVectorInitravIC(tr, bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, FIND_BIPARTITIONS_IC, bInf, &bCount, numberOfTrees, FALSE, FALSE, &tc, &tcAll, adef->verboseIC); 
1654
1655  rtc = tc / (double)(tr->mxtips - 3);
1656
1657  assert(bCount == tr->mxtips - 3); 
1658  assert(tc <= (double)(tr->mxtips - 3));
1659 
1660  printBothOpen("Tree certainty for this tree: %f\n", tc);
1661  printBothOpen("Relative tree certainty for this tree: %f\n\n", rtc);
1662
1663  rtcAll = tcAll / (double)(tr->mxtips - 3);
1664
1665  printBothOpen("Tree certainty including all conflicting bipartitions (TC-All) for this tree: %f\n", tcAll);
1666  printBothOpen("Relative tree certainty including all conflicting bipartitions (TC-All) for this tree: %f\n\n", rtcAll);
1667
1668  if(adef->verboseIC)
1669    printBothOpen("Verbose PHYLIP-style formatted bipartition information written to file: %s\n\n",  verboseSplitsFileName);
1670
1671  printBipartitionResult(tr, adef, TRUE, TRUE);   
1672
1673  freeBitVectors(bitVectors, 2 * tr->mxtips);
1674  rax_free(bitVectors);
1675  freeHashTable(h);
1676  rax_free(h); 
1677
1678  rax_free(bInf);
1679}
1680
1681
1682
1683
1684/*******************/
1685
1686static double testFreq(double *vect1, double *vect2, int n);
1687
1688
1689
1690void compareBips(tree *tr, char *bootStrapFileName, analdef *adef)
1691{
1692  int 
1693    numberOfTreesAll = 0, 
1694    numberOfTreesStop = 0,
1695    i; 
1696  unsigned int k, entryCount;
1697  double *vect1, *vect2, p, avg1 = 0.0, avg2 = 0.0, scaleAll, scaleStop;
1698  int 
1699    bipAll = 0,
1700    bipStop = 0;
1701  char bipFileName[1024];
1702  FILE 
1703    *outf,
1704    *treeFile;
1705 
1706  unsigned 
1707    int vLength;
1708 
1709  unsigned int **bitVectors = initBitVector(tr, &vLength);
1710  hashtable *h = initHashTable(tr->mxtips * 100);   
1711  unsigned long int c1 = 0;
1712  unsigned long int c2 = 0;   
1713
1714
1715  /*********************************************************************************************************/
1716 
1717  treeFile = getNumberOfTrees(tr, bootStrapFileName, adef); 
1718  numberOfTreesAll = tr->numberOfTrees;
1719             
1720  for(i = 0; i < numberOfTreesAll; i++)
1721    { 
1722      int 
1723        bCounter = 0;
1724     
1725      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);               
1726      assert(tr->mxtips == tr->ntips); 
1727
1728      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, BIPARTITIONS_ALL, (branchInfo*)NULL, &bCounter, 0, FALSE, FALSE);
1729      assert(bCounter == tr->mxtips - 3);     
1730    }
1731         
1732  fclose(treeFile);
1733
1734
1735  /*********************************************************************************************************/   
1736
1737  treeFile = getNumberOfTrees(tr, tree_file, adef);
1738 
1739  numberOfTreesStop = tr->numberOfTrees;   
1740       
1741  for(i = 0; i < numberOfTreesStop; i++)
1742    {             
1743      int 
1744        bCounter = 0;
1745
1746
1747      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);     
1748      assert(tr->mxtips == tr->ntips);
1749     
1750      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 1, BIPARTITIONS_ALL, (branchInfo*)NULL, &bCounter, 0, FALSE, FALSE);
1751      assert(bCounter == tr->mxtips - 3);     
1752    }
1753         
1754  fclose(treeFile); 
1755
1756  /***************************************************************************************************/
1757     
1758  vect1 = (double *)rax_malloc(sizeof(double) * h->entryCount);
1759  vect2 = (double *)rax_malloc(sizeof(double) * h->entryCount);
1760
1761  strcpy(bipFileName,         workdir); 
1762  strcat(bipFileName,         "RAxML_bipartitionFrequencies.");
1763  strcat(bipFileName,         run_id);
1764
1765  outf = myfopen(bipFileName, "wb");
1766
1767
1768  scaleAll  = 1.0 / (double)numberOfTreesAll;
1769  scaleStop = 1.0 / (double)numberOfTreesStop;
1770
1771  for(k = 0, entryCount = 0; k < h->tableSize; k++)         
1772    {
1773     
1774      if(h->table[k] != NULL)
1775        {
1776          entry *e = h->table[k];
1777
1778          do
1779            {
1780              c1 += e->bipNumber;
1781              c2 += e->bipNumber2;
1782              vect1[entryCount] = ((double)e->bipNumber) * scaleAll;
1783              if(vect1[entryCount] > 0)
1784                bipAll++;
1785              vect2[entryCount] = ((double)e->bipNumber2) * scaleStop;
1786              if(vect2[entryCount] > 0)
1787                bipStop++;
1788              fprintf(outf, "%f %f\n", vect1[entryCount], vect2[entryCount]);
1789              entryCount++;
1790              e = e->next;
1791            }
1792          while(e != NULL);
1793        }     
1794    }
1795 
1796  printBothOpen("%ld %ld\n", c1, c2);
1797
1798  assert(entryCount == h->entryCount);
1799
1800  fclose(outf);
1801
1802  p = testFreq(vect1, vect2, h->entryCount);
1803
1804  for(k = 0; k < h->entryCount; k++)
1805    {
1806      avg1 += vect1[k];
1807      avg2 += vect2[k];
1808    }
1809
1810  avg1 /= ((double)h->entryCount);
1811  avg2 /= ((double)h->entryCount);
1812 
1813 
1814  printBothOpen("Average [%s]: %1.40f [%s]: %1.40f\n", bootStrapFileName, avg1, tree_file, avg2);
1815  printBothOpen("Pearson: %f Bipartitions in [%s]: %d Bipartitions in [%s]: %d Total Bipartitions: %d\n", p, bootStrapFileName, bipAll, tree_file, bipStop, h->entryCount);
1816
1817  printBothOpen("\nFile containing pair-wise bipartition frequencies written to %s\n\n", bipFileName);
1818 
1819  freeBitVectors(bitVectors, 2 * tr->mxtips);
1820  rax_free(bitVectors);
1821  freeHashTable(h);
1822  rax_free(h);
1823
1824  rax_free(vect1);
1825  rax_free(vect2);
1826
1827  exit(0);
1828}
1829
1830/*************************************************************************************************************/
1831
1832void computeRF(tree *tr, char *bootStrapFileName, analdef *adef)
1833{
1834  int     
1835    treeVectorLength, 
1836    numberOfTrees = 0, 
1837    i,
1838    j, 
1839    *rfMat,
1840    *wrfMat,
1841    *wrf2Mat;
1842
1843  unsigned int
1844    vLength; 
1845
1846  unsigned int 
1847    k, 
1848    entryCount,   
1849    **bitVectors = initBitVector(tr, &vLength);
1850 
1851  char rfFileName[1024];
1852
1853  boolean computeWRF = FALSE;
1854
1855  double 
1856    maxRF, 
1857    avgRF,
1858    avgWRF,
1859    avgWRF2;
1860
1861  FILE 
1862    *outf,
1863    *treeFile = getNumberOfTrees(tr, bootStrapFileName, adef);
1864
1865  hashtable *h = (hashtable *)NULL;   
1866 
1867  numberOfTrees = tr->numberOfTrees;
1868 
1869  h = initHashTable(tr->mxtips * 2 * numberOfTrees); 
1870
1871  if(numberOfTrees % MASK_LENGTH == 0)
1872    treeVectorLength = numberOfTrees / MASK_LENGTH;
1873  else
1874    treeVectorLength = 1 + (numberOfTrees / MASK_LENGTH);
1875
1876 
1877  rfMat = (int*)rax_calloc(numberOfTrees * numberOfTrees, sizeof(int));
1878  wrfMat = (int*)rax_calloc(numberOfTrees * numberOfTrees, sizeof(int));
1879  wrf2Mat = (int*)rax_calloc(numberOfTrees * numberOfTrees, sizeof(int));
1880
1881  for(i = 0; i < numberOfTrees; i++)
1882    { 
1883      int 
1884        bCounter = 0,     
1885        lcount   = 0;
1886     
1887      lcount = treeReadLen(treeFile, tr, FALSE, TRUE, TRUE, adef, TRUE, FALSE); 
1888
1889     
1890     
1891      assert(tr->mxtips == tr->ntips); 
1892     
1893      if(i == 0)
1894        {
1895          assert(lcount == 0 || lcount == tr->mxtips - 3);
1896          if(lcount == 0)
1897            computeWRF = FALSE;
1898          else
1899            computeWRF = TRUE;
1900        }
1901     
1902      if(computeWRF)
1903        assert(lcount == tr->mxtips - 3);       
1904
1905      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, i, BIPARTITIONS_RF, (branchInfo *)NULL,
1906                              &bCounter, treeVectorLength, FALSE, computeWRF);
1907     
1908      assert(bCounter == tr->mxtips - 3);         
1909    }
1910         
1911  fclose(treeFile); 
1912
1913  for(k = 0, entryCount = 0; k < h->tableSize; k++)         
1914    {   
1915      if(h->table[k] != NULL)
1916        {
1917          entry *e = h->table[k];
1918
1919          do
1920            {
1921              unsigned int *vector = e->treeVector;
1922                             
1923              if(!computeWRF)
1924                {       
1925                  i = 0;
1926                  while(i < numberOfTrees)
1927                    {
1928                      if(vector[i / MASK_LENGTH] == 0)
1929                        i += MASK_LENGTH;
1930                      else
1931                        {                   
1932                          if((vector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]) > 0)
1933                            {                   
1934                              int *r = &rfMat[i * numberOfTrees];
1935
1936                              for(j = 0; j < numberOfTrees; j++)
1937                                if((vector[j / MASK_LENGTH] & mask32[j % MASK_LENGTH]) == 0)
1938                                  r[j]++;                                                   
1939                            }
1940                          i++;
1941                        }
1942                    }                                                     
1943                }
1944              else
1945                {
1946                  int *supportVector = e->supportVector;
1947
1948                  i = 0;
1949
1950                  while(i < numberOfTrees)
1951                    {
1952                      if(vector[i / MASK_LENGTH] == 0)
1953                        i += MASK_LENGTH;
1954                      else
1955                        {                   
1956                          if((vector[i / MASK_LENGTH] & mask32[i % MASK_LENGTH]) > 0)
1957                            {                   
1958                              int 
1959                                *r = &rfMat[i * numberOfTrees],
1960                                *w   = &wrfMat[i * numberOfTrees],
1961                                *w2  = &wrf2Mat[i * numberOfTrees],
1962                                support = supportVector[i];
1963                             
1964                             
1965
1966                              for(j = 0; j < numberOfTrees; j++)
1967                                {
1968                                  if((vector[j / MASK_LENGTH] & mask32[j % MASK_LENGTH]) == 0)
1969                                    {
1970                                      r[j]++;                                               
1971                                      w[j] += ABS(support - supportVector[j]);
1972                                      w2[j] += ABS(support - supportVector[j]);
1973                                    }
1974                                  else
1975                                    {
1976                                      w2[j] += ABS(support - supportVector[j]);
1977                                    }
1978
1979                                }
1980                            }
1981                          i++;
1982                        }
1983                    }                     
1984                }     
1985             
1986              entryCount++;
1987              e = e->next;
1988            }
1989          while(e != NULL);
1990        }
1991
1992     
1993    }
1994  assert(entryCount == h->entryCount); 
1995
1996
1997  strcpy(rfFileName,         workdir); 
1998  strcat(rfFileName,         "RAxML_RF-Distances.");
1999  strcat(rfFileName,         run_id);
2000
2001  outf = myfopen(rfFileName, "wb");
2002 
2003  maxRF  = ((double)(2 * (tr->mxtips - 3)));
2004  avgRF  = 0.0;
2005  avgWRF = 0.0;
2006  avgWRF2 = 0.0;
2007
2008  if(!computeWRF)
2009    {
2010      for(i = 0; i < numberOfTrees; i++)   
2011        for(j = i + 1; j < numberOfTrees; j++)
2012          rfMat[i * numberOfTrees + j] += rfMat[j * numberOfTrees + i];
2013    }
2014  else
2015    {
2016      for(i = 0; i < numberOfTrees; i++)   
2017        for(j = i + 1; j < numberOfTrees; j++)
2018          {
2019            rfMat[i * numberOfTrees + j] += rfMat[j * numberOfTrees + i];
2020            wrfMat[i * numberOfTrees + j] += wrfMat[j * numberOfTrees + i];
2021            wrf2Mat[i * numberOfTrees + j] += wrf2Mat[j * numberOfTrees + i];
2022          }
2023    }
2024
2025  for(i = 0; i < numberOfTrees; i++)
2026    for(j = i + 1; j < numberOfTrees; j++)
2027      {
2028        int    rf = rfMat[i * numberOfTrees + j];
2029        double rrf = (double)rf / maxRF;
2030        if(computeWRF)
2031          {
2032            double wrf = wrfMat[i * numberOfTrees + j] / 100.0;
2033            double rwrf = wrf / maxRF;
2034            double wrf2 = wrf2Mat[i * numberOfTrees + j] / 100.0;
2035            double rwrf2 = wrf2 / maxRF;
2036       
2037            fprintf(outf, "%d %d: %d %f, %f %f, %f %f\n", i, j, rf, rrf, wrf, rwrf, wrf2, rwrf2);
2038            avgWRF  += rwrf;
2039            avgWRF2 += rwrf2; 
2040          }
2041        else
2042          fprintf(outf, "%d %d: %d %f\n", i, j, rf, rrf);
2043       
2044        avgRF += rrf;
2045      }
2046 
2047  fclose(outf);
2048
2049 
2050  printBothOpen("\n\nAverage relative RF in this set: %f\n", avgRF / ((double)((numberOfTrees * numberOfTrees - numberOfTrees) / 2)));
2051  if(computeWRF)
2052    {
2053      printBothOpen("\n\nAverage relative WRF in this set: %f\n", avgWRF / ((double)((numberOfTrees * numberOfTrees - numberOfTrees) / 2)));
2054      printBothOpen("\n\nAverage relative WRF2 in this set: %f\n", avgWRF2 / ((double)((numberOfTrees * numberOfTrees - numberOfTrees) / 2)));
2055      printBothOpen("\nFile containing all %d pair-wise RF, WRF and WRF2 distances written to file %s\n\n", (numberOfTrees * numberOfTrees - numberOfTrees) / 2, rfFileName);
2056    }   
2057  else
2058    printBothOpen("\nFile containing all %d pair-wise RF distances written to file %s\n\n", (numberOfTrees * numberOfTrees - numberOfTrees) / 2, rfFileName);
2059
2060  rax_free(rfMat);
2061  rax_free(wrfMat);
2062  rax_free(wrf2Mat);
2063  freeBitVectors(bitVectors, 2 * tr->mxtips);
2064  rax_free(bitVectors);
2065  freeHashTable(h);
2066  rax_free(h); 
2067
2068  exit(0);
2069}
2070
2071/********************plausibility checker **********************************/
2072
2073/* function to extract the bit mask for the taxa that are present in the small tree */
2074
2075static void setupMask(unsigned int *smallTreeMask, nodeptr p, int numsp)
2076{
2077  if(isTip(p->number, numsp))
2078    smallTreeMask[(p->number - 1) / MASK_LENGTH] |= mask32[(p->number - 1) % MASK_LENGTH];
2079  else
2080    {   
2081      nodeptr
2082        q = p->next;
2083
2084      /* I had to change this function to account for mult-furcating trees.
2085         In this case an inner node can have more than 3 cyclically linked
2086         elements, because there might be more than 3 outgoing branches
2087         from an inner node */
2088
2089      while(q != p)
2090        {
2091          setupMask(smallTreeMask, q->back, numsp);
2092          q = q->next;
2093        }
2094     
2095      //old code below
2096      //setupMask(smallTreeMask, p->next->back, numsp);   
2097      //setupMask(smallTreeMask, p->next->next->back, numsp);     
2098    }
2099}
2100
2101/* we can not use the default hash numbers generated e.g., in the RF code, based on the tree shape.
2102   we need to compute a hash on the large/long vector that has as many bits as the big tree has taxa */
2103
2104static hashNumberType oat_hash(unsigned char *p, int len)
2105{
2106  unsigned int 
2107    h = 0;
2108  int 
2109    i;
2110 
2111  for(i = 0; i < len; i++) 
2112    {
2113      h += p[i];
2114      h += ( h << 10 );
2115      h ^= ( h >> 6 );
2116    }
2117 
2118  h += ( h << 3 );
2119  h ^= ( h >> 11 );
2120  h += ( h << 15 );
2121 
2122  return h;
2123}
2124
2125/* function that re-hashes bipartitions from the large tree into the new hash table */
2126
2127static void insertHashPlausibility(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, hashNumberType position)
2128{     
2129  if(h->table[position] != NULL)
2130    {
2131      entry
2132        *e = h->table[position];     
2133
2134      do
2135        {       
2136          unsigned int 
2137            i;
2138         
2139          for(i = 0; i < vectorLength; i++)
2140            if(bitVector[i] != e->bitVector[i])
2141              break;
2142         
2143          if(i == vectorLength)                     
2144            return;                 
2145         
2146          e = e->next;
2147        }
2148      while(e != (entry*)NULL); 
2149
2150      e = initEntry(); 
2151           
2152      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));               
2153      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
2154     
2155      e->next = h->table[position];
2156      h->table[position] = e;         
2157    }
2158  else
2159    {
2160      entry
2161        *e = initEntry(); 
2162             
2163      e->bitVector = (unsigned int*)rax_malloc(vectorLength * sizeof(unsigned int));     
2164      memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);     
2165
2166      h->table[position] = e;
2167    }
2168
2169  h->entryCount =  h->entryCount + 1;
2170}
2171
2172/* this function is called while we parse the small trees and extract the bipartitions, it will just look
2173   if the bipartition stored in bitVector is already in the hastable, if it is in there this means that
2174   this bipartition is also present in the big tree */
2175
2176static int findHash(unsigned int *bitVector, hashtable *h, unsigned int vectorLength, hashNumberType position)
2177{ 
2178  if(h->table[position] == NULL)         
2179    return 0;
2180  {
2181    entry *e = h->table[position];     
2182
2183    do
2184      { 
2185        unsigned int i;
2186
2187        for(i = 0; i < vectorLength; i++)
2188          if(bitVector[i] != e->bitVector[i])
2189            goto NEXT;
2190           
2191        return 1;       
2192      NEXT:
2193        e = e->next;
2194      }
2195    while(e != (entry*)NULL); 
2196     
2197    return 0;   
2198  }
2199}
2200
2201/* this function actually traverses the small tree, generates the bit vectors for all
2202   non-trivial bipartitions and simultaneously counts how many bipartitions (already stored in the has table) are shared with the big tree
2203*/
2204
2205static int bitVectorTraversePlausibility(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength, hashtable *h,
2206                                         int *countBranches, int firstTaxon, tree *tr, boolean multifurcating)
2207{
2208
2209  /* trivial bipartition */
2210
2211  if(isTip(p->number, numsp))
2212    return 0;
2213  else
2214    {
2215      int 
2216        found = 0;
2217
2218      nodeptr
2219        q = p->next;         
2220
2221      /* recursively descend into the tree and get the bips of all subtrees first */
2222
2223      do 
2224        {
2225          found = found + bitVectorTraversePlausibility(bitVectors, q->back, numsp, vectorLength, h, countBranches, firstTaxon, tr, multifurcating);
2226          q = q->next;
2227        }
2228      while(q != p);
2229           
2230      /* compute the bipartition induced by the current branch p, p->back,
2231         here we invoke two different functions, depending on whether we are dealing with
2232         a multi-furcating or bifurcating tree.
2233      */
2234
2235      if(multifurcating)
2236        newviewBipartitionsMultifurcating(bitVectors, p, numsp, vectorLength);
2237      else
2238        newviewBipartitions(bitVectors, p, numsp, vectorLength);
2239     
2240      assert(p->x);     
2241
2242      /* if p->back does not lead to a tip this is an inner branch that induces a non-trivial bipartition.
2243         in this case we need to lookup if the induced bipartition is already contained in the hash table
2244      */
2245
2246      if(!(isTip(p->back->number, numsp)))
2247        {       
2248          /* this is the bit vector to insert into the hash table */
2249          unsigned int 
2250            *toInsert = bitVectors[p->number];
2251         
2252          /* compute the hash number on that bit vector */
2253          hashNumberType
2254            position = oat_hash((unsigned char *)toInsert, sizeof(unsigned int) * vectorLength) % h->tableSize;         
2255
2256          /* each bipartition can be stored in two forms (the two bit-wise complements
2257             we always canonically store that version of the bit-vector that does not contain the
2258             first taxon of the small tree, we use an assertion to make sure that all is correct */
2259
2260          assert(!(toInsert[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH])); 
2261                     
2262          /* increment the branch counter to assure that all inner branches are traversed */
2263         
2264          *countBranches =  *countBranches + 1; 
2265                       
2266          /* now look up this bipartition in the hash table, If it is present the number of
2267             shared bipartitions between the small and the big tree is incremented by 1 */
2268           
2269          found = found + findHash(toInsert, h, vectorLength, position);       
2270        }
2271      return found;
2272    }
2273}
2274
2275#define _ONLY_BIFURCATING_TREES
2276
2277#ifdef _ONLY_BIFURCATING_TREES
2278
2279void plausibilityChecker(tree *tr, analdef *adef)
2280{
2281  FILE 
2282    *treeFile,
2283    *rfFile;
2284 
2285  char 
2286    rfFileName[1024];
2287 
2288  /* init has table for big reference tree */
2289 
2290  hashtable
2291    *h      = initHashTable(tr->mxtips * 2 * 2);
2292 
2293
2294  /* init the bit vectors we need for computing and storing bipartitions during
2295     the tree traversal */
2296  unsigned int 
2297    vLength, 
2298    **bitVectors = initBitVector(tr, &vLength);
2299   
2300  int
2301    branchCounter = 0,
2302    i;
2303
2304  double 
2305    avgRF = 0.0;
2306
2307  /* set up an output file name */
2308
2309  strcpy(rfFileName,         workdir); 
2310  strcat(rfFileName,         "RAxML_RF-Distances.");
2311  strcat(rfFileName,         run_id);
2312
2313  rfFile = myfopen(rfFileName, "wb"); 
2314
2315  assert(adef->mode ==  PLAUSIBILITY_CHECKER);
2316
2317  /* open the big reference tree file and parse it */
2318
2319  treeFile = myfopen(tree_file, "r");
2320
2321  printBothOpen("Parsing reference tree %s\n", tree_file);
2322
2323  treeReadLen(treeFile, tr, FALSE, TRUE, TRUE, adef, TRUE, FALSE);
2324
2325  assert(tr->mxtips == tr->ntips);
2326
2327  printBothOpen("The reference tree has %d tips\n", tr->ntips);
2328
2329  fclose(treeFile);
2330 
2331  /* extract all induced bipartitions from the big tree and store them in the hastable */
2332 
2333  bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, BIPARTITIONS_RF, (branchInfo *)NULL,
2334                          &branchCounter, 1, FALSE, FALSE);
2335     
2336  assert(branchCounter == tr->mxtips - 3);   
2337 
2338  /* now see how many small trees we have */
2339
2340  treeFile = getNumberOfTrees(tr, bootStrapFile, adef); 
2341 
2342  /* loop over all small trees */
2343
2344  for(i = 0; i < tr->numberOfTrees;  i++)
2345    {
2346      unsigned int
2347        entryCount = 0,
2348        k,
2349        j,     
2350        *masked    = (unsigned int *)rax_calloc(vLength, sizeof(unsigned int)),
2351        *smallTreeMask = (unsigned int *)rax_calloc(vLength, sizeof(unsigned int));
2352     
2353      int   
2354        bCounter = 0, 
2355        bips,
2356        firstTaxon,
2357        taxa = 0;
2358     
2359      /* allocate a has table for re-hashing the bipartitions of the big tree */
2360
2361      hashtable
2362        *rehash = initHashTable(tr->mxtips * 2 * 2);
2363
2364      double
2365        rf,
2366        maxRF;
2367     
2368      /* parse the small tree */
2369
2370      treeReadLen(treeFile, tr, FALSE, TRUE, TRUE, adef, TRUE, FALSE);
2371      printBothOpen("Small tree %d has %d tips\n", i, tr->ntips);
2372
2373      /* compute the maximum RF distance for computing the relative RF distance later-on */
2374
2375      maxRF = ((double)(2 * (tr->ntips - 3)));
2376
2377      /* now set up a bit mask where only the bits are set to one for those
2378         taxa that are actually present in the small tree we just read */
2379         
2380      setupMask(smallTreeMask, tr->start, tr->mxtips);
2381      setupMask(smallTreeMask, tr->start->back, tr->mxtips);
2382
2383      /* now get the index of the first taxon of the small tree.
2384         we will use this to unambiguously store the bipartitions
2385      */
2386
2387      firstTaxon = tr->start->number;
2388
2389      /* make sure that this bit vector is set up correctly, i.e., that
2390         it contains as many non-zero bits as there are taxa in this small tree
2391      */
2392
2393      for(j = 0; j < vLength; j++)
2394        taxa += BIT_COUNT(smallTreeMask[j]);
2395      assert(taxa == tr->ntips);
2396
2397      /* now re-hash the big tree by applying the above bit mask */
2398
2399
2400      /* loop over hash table */
2401
2402      for(k = 0, entryCount = 0; k < h->tableSize; k++)     
2403        {   
2404          if(h->table[k] != NULL)
2405            {
2406              entry *e = h->table[k];
2407
2408              /* we resolve collisions by chaining, hence the loop here */
2409
2410              do
2411                {
2412                  unsigned int 
2413                    *bitVector = e->bitVector; 
2414                 
2415                  hashNumberType
2416                    position;
2417
2418                  int 
2419                    count = 0;
2420         
2421                  /* double check that our tree mask contains the first taxon of the small tree */
2422
2423                  assert(smallTreeMask[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH]);
2424
2425                  /* if the first taxon is set then we will re-hash the bit-wise complement of the
2426                     bit vector.
2427                     The count variable is used for a small optimization */
2428
2429                  if(bitVector[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH])                   
2430                    {
2431                      //hash complement
2432                     
2433                      for(j = 0; j < vLength; j++)
2434                        {
2435                          masked[j] = (~bitVector[j]) & smallTreeMask[j];                           
2436                          count += BIT_COUNT(masked[j]);
2437                        }
2438                    }
2439                  else
2440                    {
2441                      //hash this vector
2442                     
2443                      for(j = 0; j < vLength; j++)
2444                        {
2445                          masked[j] = bitVector[j] & smallTreeMask[j]; 
2446                          count += BIT_COUNT(masked[j]);     
2447                        }
2448                    }
2449
2450                  /* note that padding the last bits is not required because they are set to 0 automatically by smallTreeMask */       
2451                 
2452                  /* make sure that we will re-hash  the canonic representation of the bipartition
2453                     where the bit for firstTaxon is set to 0!
2454                   */
2455
2456                  assert(!(masked[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH]));
2457         
2458                  /* only if the masked bipartition of the large tree is a non-trivial bipartition (two or more bits set to 1
2459                     will we re-hash it */
2460
2461                  if(count > 1)
2462                    {
2463                      /* compute hash */
2464                      position = oat_hash((unsigned char *)masked, sizeof(unsigned int) * vLength);
2465                      position = position % rehash->tableSize;
2466                     
2467                      /* re-hash to the new hash table that contains the bips of the large tree, pruned down
2468                         to the taxa contained in the small tree
2469                       */
2470                      insertHashPlausibility(masked, rehash, vLength, position);
2471                    }           
2472                 
2473                  entryCount++;
2474                 
2475                  e = e->next;
2476                }
2477              while(e != NULL);
2478            }
2479        }
2480
2481      /* make sure that we tried to re-hash all bipartitions of the original tree */
2482     
2483      assert(entryCount == (unsigned int)(tr->mxtips - 3));
2484
2485      /* now traverse the small tree and count how many bipartitions it shares
2486         with the corresponding induced tree from the large tree */
2487
2488      bips = bitVectorTraversePlausibility(bitVectors, tr->start->back, tr->mxtips, vLength, rehash, &bCounter, firstTaxon, tr, FALSE);
2489     
2490      /* compute the relative RF */
2491
2492      rf = (double)(2 * ((tr->ntips - 3) - bips)) / maxRF;           
2493
2494      avgRF += rf;
2495
2496      printBothOpen("Relative RF tree %d: %f\n\n", i, rf);
2497
2498      fprintf(rfFile, "%d %f\n", i, rf);
2499      assert(bCounter == tr->ntips - 3);         
2500
2501      /* free masks and hast table for this iteration */
2502
2503      rax_free(smallTreeMask);
2504      rax_free(masked);
2505      freeHashTable(rehash);
2506    }
2507
2508  printBothOpen("Average RF distance %f\n\n", avgRF / (double)tr->numberOfTrees);
2509
2510  printBothOpen("Total execution time: %f secs\n\n", gettime() - masterTime);
2511
2512  printBothOpen("\nFile containing all %d pair-wise RF distances written to file %s\n\n", tr->numberOfTrees, rfFileName);
2513
2514  fclose(treeFile);
2515  fclose(rfFile);   
2516 
2517  freeBitVectors(bitVectors, 2 * tr->mxtips);
2518  rax_free(bitVectors);
2519 
2520  freeHashTable(h);
2521  rax_free(h);
2522}
2523
2524#else
2525
2526void plausibilityChecker(tree *tr, analdef *adef)
2527{
2528  FILE 
2529    *treeFile,
2530    *rfFile;
2531 
2532  tree
2533    *smallTree = (tree *)rax_malloc(sizeof(tree));
2534
2535  char 
2536    rfFileName[1024];
2537 
2538  /* init hash table for big reference tree */
2539 
2540  hashtable
2541    *h      = initHashTable(tr->mxtips * 2 * 2);
2542 
2543  /* init the bit vectors we need for computing and storing bipartitions during
2544     the tree traversal */
2545  unsigned int 
2546    vLength, 
2547    **bitVectors = initBitVector(tr, &vLength);
2548   
2549  int
2550    branchCounter = 0,
2551    i;
2552
2553  double 
2554    avgRF = 0.0;
2555
2556  /* set up an output file name */
2557
2558  strcpy(rfFileName,         workdir); 
2559  strcat(rfFileName,         "RAxML_RF-Distances.");
2560  strcat(rfFileName,         run_id);
2561
2562  rfFile = myfopen(rfFileName, "wb"); 
2563
2564  assert(adef->mode ==  PLAUSIBILITY_CHECKER);
2565
2566  /* open the big reference tree file and parse it */
2567
2568  treeFile = myfopen(tree_file, "r");
2569
2570  printBothOpen("Parsing reference tree %s\n", tree_file);
2571
2572  treeReadLen(treeFile, tr, FALSE, TRUE, TRUE, adef, TRUE, FALSE);
2573
2574  assert(tr->mxtips == tr->ntips);
2575
2576  printBothOpen("The reference tree has %d tips\n", tr->ntips);
2577
2578  fclose(treeFile);
2579 
2580  /* extract all induced bipartitions from the big tree and store them in the hastable */
2581 
2582  bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vLength, h, 0, BIPARTITIONS_RF, (branchInfo *)NULL,
2583                          &branchCounter, 1, FALSE, FALSE);
2584     
2585  assert(branchCounter == tr->mxtips - 3);   
2586 
2587  /* now see how many small trees we have */
2588
2589  treeFile = getNumberOfTrees(tr, bootStrapFile, adef); 
2590 
2591  /* allocate a data structure for parsing the potentially mult-furcating tree */
2592
2593  allocateMultifurcations(tr, smallTree);
2594
2595  /* loop over all small trees */
2596
2597  for(i = 0; i < tr->numberOfTrees;  i++)
2598    {
2599      unsigned int
2600        entryCount = 0,
2601        k,
2602        j,     
2603        *masked    = (unsigned int *)rax_calloc(vLength, sizeof(unsigned int)),
2604        *smallTreeMask = (unsigned int *)rax_calloc(vLength, sizeof(unsigned int));
2605     
2606      int
2607        numberOfSplits = 0,
2608        bCounter = 0, 
2609        bips,
2610        firstTaxon,
2611        taxa = 0;
2612     
2613      /* allocate a has table for re-hashing the bipartitions of the big tree */
2614
2615      hashtable
2616        *rehash = initHashTable(tr->mxtips * 2 * 2);
2617
2618      double
2619        rf,
2620        maxRF;
2621     
2622      /* parse the small tree */             
2623
2624      /*
2625         instead of the standard tree parsing function, we parse a multi-furcating tree here.
2626         the function returns the number of inner branches/splits in the multi-furcating tree which can,
2627         of course be smaller than n-3, where n is the number of taxa in the tree.
2628      */
2629
2630      numberOfSplits = readMultifurcatingTree(treeFile, smallTree, adef);
2631      printBothOpen("Small tree %d has %d tips\n", i, smallTree->ntips);   
2632
2633      /* compute the maximum RF distance for computing the relative RF distance later-on */
2634
2635      /* note that here we need to pay attention, since the RF distance is not normalized
2636         by 2 * (n-3) but we need to account for the fact that the multifurcating small tree
2637         will potentially contain less bipartitions.
2638         Hence the normalization factor is obtained as n-3 + numberOfSplits, where n-3 is the number
2639         of bipartitions of the pruned down large reference tree for which we know that it is
2640         bifurcating/strictly binary */
2641
2642      maxRF = (double)((smallTree->ntips - 3) + numberOfSplits);
2643
2644      /* now set up a bit mask where only the bits are set to one for those
2645         taxa that are actually present in the small tree we just read */
2646         
2647      /* note that I had to apply some small changes to this function to make it work for
2648         multi-furcating trees ! */
2649
2650      setupMask(smallTreeMask, smallTree->start,       smallTree->mxtips);
2651      setupMask(smallTreeMask, smallTree->start->back, smallTree->mxtips);
2652
2653      /* now get the index of the first taxon of the small tree.
2654         we will use this to unambiguously store the bipartitions
2655      */
2656
2657      firstTaxon = smallTree->start->number;
2658
2659      /* make sure that this bit vector is set up correctly, i.e., that
2660         it contains as many non-zero bits as there are taxa in this small tree
2661      */
2662
2663      for(j = 0; j < vLength; j++)
2664        taxa += BIT_COUNT(smallTreeMask[j]);
2665      assert(taxa == smallTree->ntips);
2666
2667      /* now re-hash the big tree by applying the above bit mask */
2668
2669
2670      /* loop over hash table */
2671
2672      for(k = 0, entryCount = 0; k < h->tableSize; k++)     
2673        {   
2674          if(h->table[k] != NULL)
2675            {
2676              entry *e = h->table[k];
2677
2678              /* we resolve collisions by chaining, hence the loop here */
2679
2680              do
2681                {
2682                  unsigned int 
2683                    *bitVector = e->bitVector; 
2684                 
2685                  hashNumberType
2686                    position;
2687
2688                  int 
2689                    count = 0;
2690         
2691                  /* double check that our tree mask contains the first taxon of the small tree */
2692
2693                  assert(smallTreeMask[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH]);
2694
2695                  /* if the first taxon is set then we will re-hash the bit-wise complement of the
2696                     bit vector.
2697                     The count variable is used for a small optimization */
2698
2699                  if(bitVector[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH])                   
2700                    {
2701                      //hash complement
2702                     
2703                      for(j = 0; j < vLength; j++)
2704                        {
2705                          masked[j] = (~bitVector[j]) & smallTreeMask[j];                           
2706                          count += BIT_COUNT(masked[j]);
2707                        }
2708                    }
2709                  else
2710                    {
2711                      //hash this vector
2712                     
2713                      for(j = 0; j < vLength; j++)
2714                        {
2715                          masked[j] = bitVector[j] & smallTreeMask[j]; 
2716                          count += BIT_COUNT(masked[j]);     
2717                        }
2718                    }
2719
2720                  /* note that padding the last bits is not required because they are set to 0 automatically by smallTreeMask */       
2721                 
2722                  /* make sure that we will re-hash  the canonic representation of the bipartition
2723                     where the bit for firstTaxon is set to 0!
2724                   */
2725
2726                  assert(!(masked[(firstTaxon - 1) / MASK_LENGTH] & mask32[(firstTaxon - 1) % MASK_LENGTH]));
2727         
2728                  /* only if the masked bipartition of the large tree is a non-trivial bipartition (two or more bits set to 1
2729                     will we re-hash it */
2730
2731                  if(count > 1)
2732                    {
2733                      /* compute hash */
2734                      position = oat_hash((unsigned char *)masked, sizeof(unsigned int) * vLength);
2735                      position = position % rehash->tableSize;
2736                     
2737                      /* re-hash to the new hash table that contains the bips of the large tree, pruned down
2738                         to the taxa contained in the small tree
2739                       */
2740                      insertHashPlausibility(masked, rehash, vLength, position);
2741                    }           
2742                 
2743                  entryCount++;
2744                 
2745                  e = e->next;
2746                }
2747              while(e != NULL);
2748            }
2749        }
2750
2751      /* make sure that we tried to re-hash all bipartitions of the original tree */
2752     
2753      assert(entryCount == (unsigned int)(tr->mxtips - 3));
2754
2755      /* now traverse the small tree and count how many bipartitions it shares
2756         with the corresponding induced tree from the large tree */
2757     
2758      /* the following function also had to be modified to account for multi-furcating trees ! */
2759     
2760      bips = bitVectorTraversePlausibility(bitVectors, smallTree->start->back, smallTree->mxtips, vLength, rehash, &bCounter, firstTaxon, smallTree, TRUE);
2761     
2762      /* compute the relative RF */
2763
2764      rf = (double)(2 * ((smallTree->ntips - 3) - bips)) / maxRF;           
2765
2766      avgRF += rf;
2767
2768      printBothOpen("Relative RF tree %d: %f\n\n", i, rf);
2769
2770      fprintf(rfFile, "%d %f\n", i, rf);
2771
2772      /* I also modified this assertion, we nee to make sure here that we checked all non-trivial splits/bipartitions
2773         in the multi-furcating tree whech can be less than n - 3 ! */
2774     
2775      assert(bCounter == numberOfSplits);         
2776
2777      /* free masks and hast table for this iteration */
2778
2779      rax_free(smallTreeMask);
2780      rax_free(masked);
2781      freeHashTable(rehash);
2782      rax_free(rehash);
2783    }
2784
2785  printBothOpen("Average RF distance %f\n\n", avgRF / (double)tr->numberOfTrees);
2786
2787  printBothOpen("Total execution time: %f secs\n\n", gettime() - masterTime);
2788
2789  printBothOpen("\nFile containing all %d pair-wise RF distances written to file %s\n\n", tr->numberOfTrees, rfFileName);
2790
2791  fclose(treeFile);
2792  fclose(rfFile);   
2793 
2794  /* free the data structure used for parsing the potentially multi-furcating tree */
2795
2796  freeMultifurcations(smallTree);
2797  rax_free(smallTree);
2798
2799  freeBitVectors(bitVectors, 2 * tr->mxtips);
2800  rax_free(bitVectors);
2801 
2802  freeHashTable(h);
2803  rax_free(h);
2804}
2805
2806#endif
2807
2808/********************************************************/
2809
2810double convergenceCriterion(hashtable *h, int mxtips)
2811{
2812  int     
2813    rf = 0; 
2814
2815  unsigned int 
2816    k = 0, 
2817    entryCount = 0;
2818 
2819  double   
2820    rrf; 
2821
2822  for(k = 0, entryCount = 0; k < h->tableSize; k++)         
2823    {     
2824      if(h->table[k] != NULL)
2825        {
2826          entry *e = h->table[k];
2827
2828          do
2829            {
2830              unsigned int *vector = e->treeVector;         
2831              if(((vector[0] & 1) > 0) + ((vector[0] & 2) > 0) == 1)
2832                rf++;       
2833             
2834              entryCount++;
2835              e = e->next;
2836            }
2837          while(e != NULL);
2838        }     
2839    }
2840
2841  assert(entryCount == h->entryCount); 
2842     
2843  rrf = (double)rf/((double)(2 * (mxtips - 3))); 
2844 
2845  return rrf;
2846}
2847
2848
2849
2850
2851/*************************************************************************************************************/
2852
2853static void permute(unsigned int *perm, unsigned int n, long *seed)
2854{
2855  unsigned int  i, j, k;
2856 
2857  for (i = 0; i < n; i++) 
2858    {
2859      k =  (int)((double)(n - i) * randum(seed));
2860      j        = perm[i];   
2861      perm[i]     = perm[i + k];
2862      perm[i + k] = j; 
2863      /*assert(i + k < n);*/
2864    }
2865}
2866
2867
2868
2869
2870
2871static double testFreq(double *vect1, double *vect2, int n)
2872{
2873  int 
2874    i;
2875 
2876  boolean
2877    allEqual = TRUE;
2878
2879  double
2880    avg1 = 0.0, 
2881    avg2 = 0.0,
2882    sum_xy = 0.0, 
2883    sum_x  = 0.0, 
2884    sum_y  = 0.0,
2885    corr   = 0.0;
2886 
2887  for(i = 0; i < n; i++)
2888    {       
2889      allEqual = allEqual && (vect1[i] == vect2[i]);
2890
2891      avg1 += vect1[i];
2892      avg2 += vect2[i];
2893    }
2894     
2895  avg1 /= ((double)n);
2896  avg2 /= ((double)n); 
2897     
2898  for(i = 0; i < n; i++)
2899    {
2900      sum_xy += ((vect1[i] - avg1) * (vect2[i] - avg2));
2901      sum_x  += ((vect1[i] - avg1) * (vect1[i] - avg1));
2902      sum_y  += ((vect2[i] - avg2) * (vect2[i] - avg2));         
2903    }       
2904
2905  if(allEqual)
2906    return 1.0;
2907
2908  if(sum_x == 0.0 || sum_y == 0.0) 
2909    return 0.0;
2910
2911  corr = sum_xy / (sqrt(sum_x) * sqrt(sum_y));
2912   
2913  /*
2914    #ifndef WIN32
2915    if(isnan(corr))
2916    {
2917    printf("Numerical Error pearson correlation is not a number\n");
2918    assert(0);
2919    }
2920    #endif
2921  */
2922
2923  return corr;
2924}
2925
2926static double frequencyCriterion(int numberOfTrees, hashtable *h, int *countBetter, int bootstopPermutations)
2927{
2928  int 
2929    k, 
2930    l;
2931   
2932  long 
2933    seed = 12345;
2934
2935  double     
2936    result, 
2937    avg = 0, 
2938    *vect1, 
2939    *vect2; 
2940
2941  unsigned int
2942    *perm =  (unsigned int *)rax_malloc(sizeof(unsigned int) * numberOfTrees),
2943    j;
2944
2945  assert(*countBetter == 0);
2946         
2947#ifdef _WAYNE_MPI
2948  seed = seed + 10000 * processID;
2949#endif
2950         
2951  for(j = 0; j < (unsigned int)numberOfTrees; j++)
2952    perm[j] = j;
2953         
2954  for(k = 0; k < bootstopPermutations; k++)
2955    {                                         
2956      unsigned int entryCount = 0;
2957
2958      permute(perm, numberOfTrees, &seed);
2959     
2960     
2961
2962      vect1 = (double *)rax_calloc(h->entryCount, sizeof(double));
2963      vect2 = (double *)rax_calloc(h->entryCount, sizeof(double));           
2964
2965       
2966     
2967      for(j = 0; j < h->tableSize; j++)
2968        {               
2969          if(h->table[j] != NULL)
2970            {           
2971              entry *e = h->table[j];
2972             
2973              do
2974                {
2975                  unsigned int *set = e->treeVector;       
2976                 
2977                  for(l = 0; l < numberOfTrees; l++)
2978                    {                       
2979                      if((set[l / MASK_LENGTH] != 0) && (set[l / MASK_LENGTH] & mask32[l % MASK_LENGTH]))
2980                        {
2981                          if(perm[l] % 2 == 0)
2982                            vect1[entryCount] = vect1[entryCount] + 1.0;
2983                          else                 
2984                            vect2[entryCount] = vect2[entryCount] + 1.0;
2985                        }
2986                    }
2987                  entryCount++;
2988                  e = e->next;
2989                }
2990              while(e != NULL);
2991            }                       
2992        }                       
2993     
2994     
2995     
2996     
2997      assert(entryCount == h->entryCount);
2998     
2999     
3000
3001      result = testFreq(vect1, vect2, entryCount);
3002         
3003         
3004     
3005      if(result >= FC_LOWER)
3006        *countBetter = *countBetter + 1;
3007             
3008      avg += result;
3009             
3010      rax_free(vect1);           
3011      rax_free(vect2);           
3012    }
3013
3014  rax_free(perm);
3015         
3016  avg /= bootstopPermutations;
3017       
3018     
3019
3020  return avg;
3021}
3022
3023
3024
3025
3026static double wcCriterion(int numberOfTrees, hashtable *h, int *countBetter, double *wrf_thresh_avg, double *wrf_avg, tree *tr, unsigned int vectorLength, int bootstopPermutations)
3027{
3028  int 
3029    k, 
3030    l,   
3031    wrf,
3032    mr_thresh = ((double)numberOfTrees/4.0);
3033   
3034  unsigned int 
3035    *perm =  (unsigned int *)rax_malloc(sizeof(unsigned int) * numberOfTrees),
3036    j;
3037
3038  long seed = 12345; 
3039 
3040  double 
3041    wrf_thresh = 0.0,
3042    pct_avg = 0.0;
3043
3044#ifdef _WAYNE_MPI
3045  seed = seed + 10000 * processID;
3046#endif
3047
3048  assert(*countBetter == 0 && *wrf_thresh_avg == 0.0 && *wrf_avg == 0.0);
3049                 
3050  for(j = 0; j < (unsigned int)numberOfTrees; j++)
3051    perm[j] = j;
3052         
3053  for(k = 0; k < bootstopPermutations; k++)
3054    {                                     
3055      int mcnt1 = 0;                     
3056      int mcnt2 = 0;
3057      unsigned int entryCount = 0;
3058      double halfOfConsideredBips = 0.0;
3059
3060      entry ** sortedByKeyA = (entry **)NULL;
3061      entry ** sortedByKeyB = (entry **)NULL;
3062      int lenA, lenB;
3063      boolean ignoreCompatibilityP;
3064
3065      int iA, iB;
3066      wrf = 0;
3067             
3068      permute(perm, numberOfTrees, &seed);     
3069   
3070      for(j = 0; j < h->tableSize; j++)
3071        {               
3072          if(h->table[j] != NULL)
3073            {
3074              entry *e = h->table[j];
3075
3076                do
3077                  {
3078                    int cnt1 = 0;
3079                    int cnt2 = 0;
3080
3081                    unsigned int *set = e->treeVector;
3082
3083                    for(l = 0; l < numberOfTrees; l++)
3084                      {                     
3085                        if((set[l / MASK_LENGTH] != 0) && (set[l / MASK_LENGTH] & mask32[l % MASK_LENGTH]))
3086                          {                         
3087                            if(perm[l] % 2 == 0)
3088                              cnt1++;
3089                            else                       
3090                              cnt2++;
3091                          }                         
3092                      }
3093                   
3094                    switch(tr->bootStopCriterion)
3095                      {
3096                      case MR_STOP:
3097                        if(cnt1 <= mr_thresh)                         
3098                          cnt1 = 0;
3099                       
3100                        if(cnt2 <= mr_thresh)       
3101                          cnt2 = 0;
3102
3103                        if(cnt1 > 0)                         
3104                          mcnt1++;
3105                       
3106                        if(cnt2 > 0)                         
3107                          mcnt2++;
3108
3109                        wrf += ((cnt1 > cnt2) ? cnt1 - cnt2 : cnt2 - cnt1);
3110                        break;
3111                      case MRE_STOP:
3112                      case MRE_IGN_STOP:
3113                        e->supportFromTreeset[0] = cnt1;
3114                        e->supportFromTreeset[1] = cnt2;
3115                        break;
3116                      default:
3117                        assert(0);
3118                      }
3119
3120                    entryCount++;
3121                    e = e->next;
3122                  }
3123                while(e != NULL);
3124            }                           
3125        }       
3126     
3127      assert(entryCount == h->entryCount);
3128
3129      if((tr->bootStopCriterion == MRE_STOP) || (tr->bootStopCriterion == MRE_IGN_STOP))
3130        {
3131       
3132         
3133          if (tr->bootStopCriterion == MRE_IGN_STOP)
3134            ignoreCompatibilityP = TRUE;
3135          else
3136            ignoreCompatibilityP = FALSE;
3137           
3138         
3139         
3140          mre(h, ignoreCompatibilityP, &sortedByKeyA, &lenA, 0, tr->mxtips, vectorLength, TRUE, tr, TRUE);
3141          mre(h, ignoreCompatibilityP, &sortedByKeyB, &lenB, 1, tr->mxtips, vectorLength, TRUE, tr, TRUE);
3142         
3143           
3144          mcnt1 = lenA;
3145          mcnt2 = lenB;
3146         
3147          iA = iB = 0;
3148         
3149          while(iA < mcnt1 || iB < mcnt2)
3150            {       
3151              if( iB == mcnt2 || (iA < mcnt1 && sortedByKeyA[iA] < sortedByKeyB[iB]) )
3152                {
3153                  wrf += sortedByKeyA[iA]->supportFromTreeset[0];
3154                  iA++;
3155                }
3156              else
3157                {
3158                  if( iA == mcnt1 || (iB < mcnt2 && sortedByKeyB[iB] < sortedByKeyA[iA]) )
3159                    {
3160                      wrf += sortedByKeyB[iB]->supportFromTreeset[1];
3161                      iB++;
3162                    }
3163                  else
3164                    {
3165                      int cnt1, cnt2;
3166                     
3167                      assert (sortedByKeyA[iA] == sortedByKeyB[iB]);
3168                     
3169                      cnt1 = sortedByKeyA[iA]->supportFromTreeset[0];
3170                      cnt2 = sortedByKeyB[iB]->supportFromTreeset[1];
3171                     
3172                      wrf += ((cnt1 > cnt2) ? cnt1 - cnt2 : cnt2 - cnt1);
3173                     
3174                      iA++; 
3175                      iB++;
3176                    }
3177                }
3178            }
3179
3180          rax_free(sortedByKeyA);
3181          rax_free(sortedByKeyB);
3182         
3183          assert (iA == mcnt1);
3184          assert (iB == mcnt2);   
3185        }
3186
3187      halfOfConsideredBips = ( ((((double)numberOfTrees/2.0) * (double)mcnt1)) + ((((double)numberOfTrees/2.0) * (double)mcnt2)) );
3188
3189      /*
3190         wrf_thresh is the 'custom' threshold computed for this pair
3191         of majority rules trees (i.e. one of the BS_PERMS splits),
3192         and simply takes into account the resolution of the two trees
3193      */
3194
3195      wrf_thresh = (tr->wcThreshold) * halfOfConsideredBips;   
3196     
3197      /*
3198        we count this random split as 'succeeding' when
3199         the wrf between maj rules trees is exceeded
3200         by its custom threshold
3201      */
3202
3203      if((double)wrf <= wrf_thresh)                             
3204        *countBetter = *countBetter + 1;
3205
3206      /*
3207         here we accumulate outcomes and thresholds, because
3208         we're not going to stop until the avg dist is less
3209         than the avg threshold
3210      */
3211
3212      pct_avg += (double)wrf / halfOfConsideredBips  * 100.0;
3213      *wrf_avg += (double)wrf;
3214      *wrf_thresh_avg += wrf_thresh;
3215    }
3216 
3217  rax_free(perm);
3218
3219  pct_avg /= (double)bootstopPermutations; 
3220  *wrf_avg /= (double)bootstopPermutations; 
3221  *wrf_thresh_avg /= (double)bootstopPermutations;   
3222
3223  /*printf("%d \t\t %f \t\t %d \t\t\t\t %f\n", numberOfTrees, *wrf_avg, *countBetter, *wrf_thresh_avg);               */
3224
3225  return pct_avg; 
3226}         
3227
3228
3229
3230
3231
3232
3233void computeBootStopOnly(tree *tr, char *bootStrapFileName, analdef *adef)
3234{
3235  int numberOfTrees = 0, i;
3236  boolean stop = FALSE;
3237  double avg;
3238  int checkEvery;
3239  int treesAdded = 0;
3240  hashtable *h = initHashTable(tr->mxtips * FC_INIT * 10); 
3241  unsigned int 
3242    treeVectorLength, 
3243    vectorLength;
3244  unsigned int **bitVectors = initBitVector(tr, &vectorLength);   
3245 
3246
3247  FILE 
3248    *treeFile = getNumberOfTrees(tr, bootStrapFileName, adef);
3249
3250  assert((FC_SPACING % 2 == 0) && (FC_THRESHOLD < BOOTSTOP_PERMUTATIONS));
3251   
3252  numberOfTrees = tr->numberOfTrees;
3253 
3254 
3255  printBothOpen("\n\nFound %d trees in File %s\n\n", numberOfTrees, bootStrapFileName);
3256 
3257  assert(sizeof(unsigned char) == 1);
3258 
3259  if(numberOfTrees % MASK_LENGTH == 0)
3260    treeVectorLength = numberOfTrees / MASK_LENGTH;
3261  else
3262    treeVectorLength = 1 + (numberOfTrees / MASK_LENGTH); 
3263 
3264  checkEvery = FC_SPACING;
3265       
3266  switch(tr->bootStopCriterion)
3267    {
3268    case FREQUENCY_STOP:
3269      printBothOpen("# Trees \t Average Pearson Coefficient \t # Permutations: pearson >= %f\n", 
3270                    FC_LOWER);
3271      break;
3272    case MR_STOP:
3273    case MRE_STOP:
3274    case MRE_IGN_STOP:
3275      printBothOpen("# Trees \t Avg WRF in %s \t # Perms: wrf <= %1.2f %s\n","%", 100.0 * tr->wcThreshold, "%");
3276      break;
3277    default:   
3278      assert(0);
3279    }
3280 
3281  for(i = 1; i <= numberOfTrees && !stop; i++)
3282    {                 
3283      int 
3284        bCount = 0;           
3285     
3286
3287      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE); 
3288      assert(tr->mxtips == tr->ntips);
3289     
3290      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, (i - 1), BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
3291                              &bCount, treeVectorLength, FALSE, FALSE);
3292     
3293      assert(bCount == tr->mxtips - 3);
3294                 
3295      treesAdded++;     
3296           
3297      if((i > START_BSTOP_TEST) && (i % checkEvery == 0))
3298        { 
3299          int countBetter = 0;
3300                 
3301          switch(tr->bootStopCriterion)
3302            {
3303            case FREQUENCY_STOP:
3304              avg = frequencyCriterion(i, h, &countBetter, BOOTSTOP_PERMUTATIONS);                       
3305              printBothOpen("%d \t\t\t %f \t\t\t\t %d\n", i, avg, countBetter);
3306                 
3307              stop = (countBetter >= FC_THRESHOLD && avg >= FC_LOWER);           
3308              break;
3309            case MR_STOP:
3310            case MRE_STOP:
3311            case MRE_IGN_STOP:
3312              {
3313                double 
3314                  wrf_thresh_avg = 0.0,
3315                  wrf_avg = 0.0;
3316                avg = wcCriterion(i, h, &countBetter, &wrf_thresh_avg, &wrf_avg, tr, vectorLength, BOOTSTOP_PERMUTATIONS);
3317                printBothOpen("%d \t\t %1.2f \t\t\t %d\n", i, avg, countBetter);               
3318               
3319                stop = (countBetter >= FC_THRESHOLD && wrf_avg <= wrf_thresh_avg);
3320              }
3321              break;
3322            default:
3323              assert(0);
3324            }   
3325        }                 
3326     
3327    }
3328 
3329 
3330
3331  if(stop)             
3332    printBothOpen("Converged after %d replicates\n", treesAdded);           
3333  else   
3334    printBothOpen("Bootstopping test did not converge after %d trees\n", treesAdded);     
3335
3336  fclose(treeFile); 
3337 
3338  freeBitVectors(bitVectors, 2 * tr->mxtips);
3339  rax_free(bitVectors);
3340  freeHashTable(h);
3341  rax_free(h);
3342
3343 
3344
3345 
3346  exit(0);
3347}
3348
3349#ifdef _WAYNE_MPI
3350
3351boolean computeBootStopMPI(tree *tr, char *bootStrapFileName, analdef *adef, double *pearsonAverage)
3352{
3353  boolean
3354    stop = FALSE;
3355
3356  int 
3357    bootStopPermutations = 0,
3358    numberOfTrees = 0, 
3359    i,
3360    countBetter = 0;
3361 
3362  unsigned int
3363    treeVectorLength, 
3364    vectorLength;
3365
3366  double 
3367    avg;
3368 
3369  hashtable
3370    *h = initHashTable(tr->mxtips * FC_INIT * 10); 
3371 
3372  unsigned 
3373    int **bitVectors = initBitVector(tr, &vectorLength);   
3374
3375 
3376  FILE 
3377    *treeFile = getNumberOfTrees(tr, bootStrapFileName, adef);
3378 
3379  numberOfTrees = tr->numberOfTrees;
3380
3381  if(numberOfTrees % 2 != 0)
3382    numberOfTrees--;
3383   
3384  /*printf("\n\nProcess %d Found %d trees in File %s\n\n", processID, numberOfTrees, bootStrapFileName);*/
3385 
3386  assert(sizeof(unsigned char) == 1);
3387 
3388
3389  if(BOOTSTOP_PERMUTATIONS % processes == 0)
3390    bootStopPermutations = BOOTSTOP_PERMUTATIONS / processes;
3391  else
3392    bootStopPermutations = 1 + (BOOTSTOP_PERMUTATIONS / processes);
3393 
3394  /*printf("Perms %d\n",  bootStopPermutations);*/
3395
3396  if(numberOfTrees % MASK_LENGTH == 0)
3397    treeVectorLength = numberOfTrees / MASK_LENGTH;
3398  else
3399    treeVectorLength = 1 + (numberOfTrees / MASK_LENGTH);   
3400 
3401  for(i = 1; i <= numberOfTrees; i++)
3402    {                 
3403      int 
3404        bCount = 0;         
3405     
3406      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE); 
3407      assert(tr->mxtips == tr->ntips);
3408     
3409      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, (i - 1), BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
3410                              &bCount, treeVectorLength, FALSE, FALSE);   
3411      assert(bCount == tr->mxtips - 3);
3412    }                                     
3413             
3414  switch(tr->bootStopCriterion)
3415    {
3416    case FREQUENCY_STOP:
3417      {
3418        double 
3419          allOut[2],
3420          allIn[2];
3421       
3422        avg = frequencyCriterion(numberOfTrees, h, &countBetter, bootStopPermutations);                   
3423       
3424        /*printf("%d \t\t\t %f \t\t\t\t %d\n", numberOfTrees, avg, countBetter);*/
3425       
3426        allOut[0] = (double)countBetter;
3427        allOut[1] = avg;
3428
3429        MPI_Allreduce(allOut, allIn, 2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
3430       
3431        /*printf("%d %f %f\n", processID, allIn[0], allIn[1]);*/
3432
3433        stop = (((int)allIn[0]) >= FC_THRESHOLD && (allIn[1] / ((double)processes)) >= FC_LOWER);
3434
3435        *pearsonAverage = (allIn[1] / ((double)processes));
3436      }
3437      break;
3438    case MR_STOP:
3439    case MRE_STOP:
3440    case MRE_IGN_STOP:
3441      {
3442        double 
3443          allOut[4],
3444          allIn[4];
3445       
3446        double 
3447          wrf_thresh_avg = 0.0,
3448          wrf_avg = 0.0;
3449       
3450        avg = wcCriterion(numberOfTrees, h, &countBetter, &wrf_thresh_avg, &wrf_avg, tr, vectorLength, bootStopPermutations);
3451       
3452        /*printf("%d %1.2f  %d %f %f\n", numberOfTrees, avg, countBetter, wrf_thresh_avg, wrf_avg);*/
3453
3454        allOut[0] = (double)countBetter;
3455        allOut[1] = wrf_thresh_avg;
3456        allOut[2] = wrf_avg;
3457        allOut[3] = avg;
3458
3459        MPI_Allreduce(allOut, allIn, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
3460       
3461        /*printf("%d %f %f %f\n", processID, allIn[0], allIn[1], allIn[2]);*/
3462       
3463        stop = (((int)allIn[0]) >= FC_THRESHOLD && (allIn[2] / ((double)processes)) <= (allIn[1] / ((double)processes)));
3464
3465        *pearsonAverage = (allIn[3] / ((double)processes));
3466      }
3467      break;
3468    default:
3469      assert(0);
3470    }
3471       
3472  fclose(treeFile); 
3473 
3474  freeBitVectors(bitVectors, 2 * tr->mxtips);
3475  rax_free(bitVectors);
3476  freeHashTable(h);
3477  rax_free(h);
3478
3479  return stop;
3480}
3481
3482#endif
3483
3484boolean bootStop(tree *tr, hashtable *h, int numberOfTrees, double *pearsonAverage, unsigned int **bitVectors, int treeVectorLength, unsigned int vectorLength)
3485{
3486  int 
3487    n = numberOfTrees + 1,
3488    bCount = 0;
3489
3490  assert((FC_SPACING % 2 == 0) && (FC_THRESHOLD < BOOTSTOP_PERMUTATIONS));
3491  assert(tr->mxtips == tr->rdta->numsp);
3492
3493  bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, numberOfTrees, BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
3494                          &bCount, treeVectorLength, FALSE, FALSE);
3495  assert(bCount == tr->mxtips - 3); 
3496
3497  if((n > START_BSTOP_TEST) && (n % FC_SPACING == 0))
3498    {     
3499      int countBetter = 0;
3500
3501      switch(tr->bootStopCriterion)
3502        {
3503        case FREQUENCY_STOP:
3504          *pearsonAverage = frequencyCriterion(n, h, &countBetter, BOOTSTOP_PERMUTATIONS);                                     
3505
3506          if(countBetter >= FC_THRESHOLD && *pearsonAverage >= FC_LOWER)
3507            return TRUE;
3508          else
3509            return FALSE;
3510          break;
3511        case MR_STOP:
3512        case MRE_STOP:
3513        case MRE_IGN_STOP:
3514         {         
3515           double 
3516             wrf_thresh_avg = 0.0,
3517             wrf_avg = 0.0;
3518           
3519           *pearsonAverage = wcCriterion(n, h, &countBetter, &wrf_thresh_avg, &wrf_avg, tr, vectorLength, BOOTSTOP_PERMUTATIONS);
3520         
3521           if(countBetter >= FC_THRESHOLD && wrf_avg <= wrf_thresh_avg)
3522             return TRUE;
3523           else
3524             return FALSE;
3525         } 
3526        default:
3527          assert(0);
3528        }
3529    }
3530  else
3531    return FALSE;
3532}
3533
3534
3535
3536
3537/* consensus stuff */
3538
3539boolean compatible(entry* e1, entry* e2, unsigned int bvlen)
3540{
3541  unsigned int i;
3542
3543  unsigned int 
3544    *A = e1->bitVector,
3545    *C = e2->bitVector;
3546 
3547  for(i = 0; i < bvlen; i++)       
3548    if(A[i] & C[i])
3549      break;
3550         
3551  if(i == bvlen)
3552    return TRUE;
3553 
3554  for(i = 0; i < bvlen; i++)
3555    if(A[i] & ~C[i])
3556      break;
3557   
3558  if(i == bvlen) 
3559    return TRUE; 
3560 
3561  for(i = 0; i < bvlen; i++)
3562    if(~A[i] & C[i])
3563      break;
3564 
3565  if(i == bvlen)
3566    return TRUE; 
3567  else
3568    return FALSE;
3569}
3570
3571
3572
3573static int sortByWeight(const void *a, const void *b, int which)
3574{
3575  /* recall, we want to sort descending, instead of ascending */
3576     
3577  int 
3578    ca,
3579    cb;
3580   
3581  ca = ((*((entry **)a))->supportFromTreeset)[which];
3582  cb = ((*((entry **)b))->supportFromTreeset)[which];
3583 
3584  if (ca == cb) 
3585    return 0;
3586 
3587  return ((ca<cb)?1:-1);
3588}
3589
3590static int sortByIndex(const void *a, const void *b)
3591{
3592  if ( (*((entry **)a)) == (*((entry **)b)) ) return 0;
3593  return (( (*((entry **)a)) < (*((entry **)b)) )?-1:1);
3594}
3595
3596static int _sortByWeight0(const void *a, const void *b)
3597{
3598  return sortByWeight(a,b,0);
3599}
3600
3601static int _sortByWeight1(const void *a, const void *b)
3602{
3603  return sortByWeight(a,b,1);
3604}
3605
3606boolean issubset(unsigned int* bipA, unsigned int* bipB, unsigned int vectorLen, unsigned int firstIndex)
3607{
3608  unsigned int 
3609    i;
3610 
3611  for(i = firstIndex; i < vectorLen; i++)
3612    if((bipA[i] & bipB[i]) != bipA[i])   
3613      return FALSE;
3614       
3615  return TRUE; 
3616}
3617
3618
3619
3620
3621
3622#ifdef _NEW_MRE
3623
3624static void mre(hashtable *h, boolean icp, entry*** sbi, int* len, int which, int n, unsigned int vectorLength, boolean sortp, tree *tr, boolean bootStopping)
3625{
3626  entry
3627    **sbw;
3628 
3629  unsigned int   
3630    i = 0,
3631    j = 0;
3632 
3633  sbw = (entry **) rax_calloc(h->entryCount, sizeof(entry *));
3634   
3635  for(i = 0; i < h->tableSize; i++) /* copy hashtable h to list sbw */
3636    {           
3637      if(h->table[i] != NULL)
3638        {
3639          entry
3640            *e = h->table[i];
3641         
3642          do
3643            {
3644              sbw[j] = e;
3645              j++;
3646              e = e->next;
3647            }
3648         
3649          while(e != NULL);
3650        }
3651    }
3652
3653  assert(h->entryCount == j);
3654
3655  if(which == 0)                /* sort the sbw list */
3656    qsort(sbw, h->entryCount, sizeof(entry *), _sortByWeight0);
3657  else
3658    qsort(sbw, h->entryCount, sizeof(entry *), _sortByWeight1);
3659
3660  *sbi = (entry **)rax_calloc(n - 3, sizeof(entry *));
3661
3662  *len = 0;
3663
3664  if(icp == FALSE)
3665    {
3666     
3667#ifdef _USE_PTHREADS
3668      /*
3669        We only deploy the parallel version of MRE when not using it
3670        in conjunction with bootstopping for the time being.
3671        When bootstopping it is probably easier and more efficient to
3672        parallelize over the permutations
3673      */
3674
3675      if(!bootStopping)
3676        {
3677          //printf("Parallel region \n" );
3678
3679          tr->h = h;
3680          NumberOfJobs = tr->h->entryCount;
3681          tr->sectionEnd = MIN(NumberOfJobs, NumberOfThreads * MRE_MIN_AMOUNT_JOBS_PER_THREAD); //NumberOfThreads * MRE_MIN_AMOUNT_JOBS_PER_THREAD;
3682          tr->len = len;
3683          tr->sbi = (*sbi);
3684          tr->maxBips = n - 3; 
3685          tr->recommendedAmountJobs = 1;
3686          tr->bitVectorLength = vectorLength;
3687          tr->sbw = sbw;
3688          tr->entriesOfSection = tr->sbw;
3689          tr->bipStatus = (int*)rax_calloc(tr->sectionEnd, sizeof(int));
3690          tr->bipStatusLen = tr->sectionEnd;
3691          masterBarrier(THREAD_MRE_COMPUTE, tr);
3692        }
3693      else
3694#endif
3695        {
3696          for(i = 0; i < h->entryCount && (*len) < n-3; i++)
3697            {                   
3698              boolean
3699                compatflag = TRUE;
3700             
3701              entry
3702                *currentEntry = sbw[i];             
3703             
3704              assert(*len < n-3);
3705             
3706              if(currentEntry->supportFromTreeset[which] <= ((unsigned int)tr->mr_thresh))
3707                {
3708                  int k;
3709                 
3710                  for(k = (*len); k > 0; k--)
3711                    {
3712                      if( ! compatible((*sbi)[k-1], currentEntry, vectorLength))
3713                        {
3714                          compatflag = FALSE;       
3715                          break;
3716                        }
3717                    }
3718                }
3719             
3720              if(compatflag)
3721                {             
3722                  (*sbi)[*len] = sbw[i];
3723                  (*len)++;
3724                }         
3725            }
3726        }
3727    }
3728  else 
3729    {     
3730      for(i = 0; i < (unsigned int)(n-3); i++)
3731        {
3732          (*sbi)[i] = sbw[i];
3733          (*len)++;
3734        }
3735    }
3736
3737
3738  rax_free(sbw);
3739
3740  if (sortp == TRUE)
3741    qsort(*sbi, (*len), sizeof(entry *), sortByIndex);
3742
3743  return;
3744}
3745
3746
3747/* if we encounter the first bits that are set, then we can determine,
3748   whether bip a is a subset of bip b.  We already know, that A has
3749   more bits set than B and that both bips are compatible to each
3750   other. Thus, if A & B is true (and A contains bits), then A MUST be
3751   a proper subset of B (given the setting). */
3752
3753/* check different versions of this ! */
3754
3755
3756
3757
3758static int sortByAmountTips(const void *a, const void *b)
3759{                               
3760  entry
3761    *A = (*(entry **)a),
3762    *B = (*(entry **)b);
3763 
3764  if((unsigned int)A->amountTips == (unsigned int)B->amountTips)
3765    return 0; 
3766
3767  return (((unsigned int)A->amountTips < (unsigned int)B->amountTips) ?  -1 : 1); 
3768}
3769
3770
3771/******* IC function *******************/
3772
3773
3774static void calculateIC(tree *tr, hashtable *h, unsigned int *bitVector, unsigned int vectorLength, int trees, unsigned int supportedBips, double *ic, double *icAll, boolean verboseIC, int counter)
3775{
3776  unsigned int
3777    maxCounter = 0,
3778    *maxima = (unsigned int *)rax_calloc(h->entryCount, sizeof(unsigned int)),
3779    **maximaBitVectors = (unsigned int **)rax_calloc(h->entryCount, sizeof(unsigned int *)),
3780    numberOfTrees = (unsigned int)trees;
3781
3782  *ic = 0.0,
3783  *icAll = 0.0;
3784
3785  //if the support is 100% we don't need to consider any conflicting bipartitions and can save some time
3786
3787  if(supportedBips == numberOfTrees)
3788    {
3789      *ic = 1.0;
3790      *icAll = 1.0;
3791     
3792      if(verboseIC)
3793        printFullySupportedSplit(tr, bitVector, numberOfTrees);
3794    }
3795  else
3796    {
3797      //search conflicting bipartitions
3798
3799      countIncompatibleBipartitions(bitVector, h, vectorLength, maxima, &maxCounter, FALSE, numberOfTrees, maximaBitVectors);     
3800     
3801      //make sure that the sum of raw supports is not higher than the number of trees
3802
3803      assert(supportedBips + maxima[0] <= numberOfTrees);
3804     
3805      *ic    = computeIC_Value(supportedBips, maxima, numberOfTrees, maxCounter, FALSE, TRUE);
3806      *icAll = computeIC_Value(supportedBips, maxima, numberOfTrees, maxCounter, TRUE, FALSE); 
3807
3808       if(verboseIC)                 
3809         printVerboseIC(tr, supportedBips, bitVector, maxCounter, maxima, maximaBitVectors, numberOfTrees, counter, *ic);
3810    }
3811
3812  //printf("IC %f %f IC-all %f %fmaxima: %u\n", ic, _ic, icAll, _icAll, maxCounter);
3813
3814  rax_free(maxima);
3815  rax_free(maximaBitVectors);
3816}
3817
3818/******* IC function end ***************/
3819
3820static void printBipsRecursive(tree *tr, FILE *outf, int consensusBipLen, entry **consensusBips, int numberOfTrees, 
3821                               int currentBipIdx, List **listOfDirectChildren, int bitVectorLength, int numTips, 
3822                               char **nameList, entry *currentBip, boolean *printed, boolean topLevel, unsigned int *printCounter, hashtable
3823                               *h, boolean computeIC, double *tc, double *tcAll, boolean verboseIC)
3824{
3825  List
3826    *idx; 
3827 
3828  int 
3829    i;
3830 
3831  unsigned int 
3832    *currentBitVector = (unsigned int*)rax_malloc(bitVectorLength * sizeof(unsigned int)); 
3833
3834  /* open bip */
3835  if(*printed)
3836    fprintf(outf, ",");
3837  *printed = FALSE;
3838   
3839  if(!topLevel)
3840    fprintf(outf, "(");   
3841
3842  /* determine tips that are not in sub bipartitions */
3843  for(i = 0; i < bitVectorLength; i++)
3844    {     
3845      idx = listOfDirectChildren[currentBipIdx]; 
3846      currentBitVector[i] = currentBip->bitVector[i];
3847     
3848      while(idx)
3849        {
3850          currentBitVector[i] = currentBitVector[i] & ~ consensusBips[*((int*)idx->value)]->bitVector[i]; 
3851          idx = idx->next;
3852        }
3853    }
3854
3855  /* print out those tips that are direct leafs of the current bip */
3856  for(i = 0; i < numTips; i++)
3857    {
3858    if(currentBitVector[i/MASK_LENGTH] & mask32[i%MASK_LENGTH])
3859      {
3860        if(*printed){fprintf(outf, ",");};
3861        fprintf(outf, "%s", nameList[i+1]);   
3862        *printed = TRUE;
3863      }
3864    }
3865
3866  /* process all sub bips */   
3867  idx = listOfDirectChildren[currentBipIdx]; 
3868  while(idx)
3869    {
3870   
3871      if(*printed)
3872        {
3873          fprintf(outf, ",");
3874          *printed = FALSE;
3875        } 
3876     
3877      printBipsRecursive(tr, outf, consensusBipLen, consensusBips, numberOfTrees, 
3878                         *((int*)idx->value), listOfDirectChildren, bitVectorLength, numTips, nameList, 
3879                         consensusBips[*((int*)idx->value)], printed, FALSE, printCounter, h, computeIC, tc, tcAll, verboseIC);
3880      *printed  = TRUE;
3881      idx = idx->next; 
3882    }
3883
3884  /* close the bipartition */
3885  if(currentBipIdx != consensusBipLen)
3886    {
3887      if(computeIC)
3888        {
3889          double 
3890            ic,
3891            icAll;
3892         
3893          calculateIC(tr, h, currentBip->bitVector, bitVectorLength, numberOfTrees, currentBip->supportFromTreeset[0], &ic, &icAll, verboseIC, *printCounter);
3894
3895          *tc    += ic;
3896          *tcAll += icAll;
3897
3898          fprintf(outf,"):1.0[%1.2f,%1.2f]", ic, icAll);
3899        }
3900      else
3901        {
3902          double 
3903            support = ((double)(currentBip->supportFromTreeset[0])) / ((double) (numberOfTrees));
3904         
3905          int 
3906            branchLabel = (int)(0.5 + support * 100.0);
3907         
3908          fprintf(outf,"):1.0[%d]", branchLabel);
3909        }
3910     
3911      *printCounter = *printCounter + 1;
3912    }
3913  else
3914    fprintf(outf, ");\n");
3915 
3916  rax_free(currentBitVector); 
3917}
3918
3919
3920
3921
3922static void printSortedBips(entry **consensusBips, const int consensusBipLen, const int numTips, const unsigned int vectorLen, 
3923                            const int numberOfTrees, FILE *outf, char **nameList , tree *tr, unsigned int *printCounter, hashtable *h, boolean computeIC, boolean verboseIC)
3924{
3925  int 
3926    i;
3927
3928  double
3929    tc = 0.0,
3930    tcAll = 0.0;
3931
3932  List
3933    **listOfDirectChildren = (List**) rax_calloc(consensusBipLen + 1, sizeof(List*)); /* reserve one more: the last one is the bip with all species */
3934 
3935  boolean
3936    *hasAncestor = (boolean*) rax_calloc(consensusBipLen, sizeof(boolean)),
3937    *printed = (boolean*)rax_calloc(1, sizeof(boolean));
3938 
3939  entry
3940    *topBip; 
3941
3942#ifndef _USE_PTHREADS
3943  List
3944    *elems = (List*)rax_malloc((size_t)consensusBipLen *  sizeof(List));
3945  int 
3946    *intList = (int*)rax_malloc(sizeof(int) * (size_t)consensusBipLen); 
3947#endif
3948
3949  /* sort the consensusBips by the amount of tips they contain */
3950 
3951  for( i = 0; i < consensusBipLen; i++)
3952    consensusBips[i]->amountTips = genericBitCount(consensusBips[i]->bitVector, vectorLen); 
3953
3954  qsort(consensusBips, consensusBipLen, sizeof(entry *), &sortByAmountTips);
3955
3956  /* create an artificial entry for the top */
3957  topBip = (entry *)rax_malloc(sizeof(entry));
3958  topBip->bitVector = rax_malloc(sizeof(unsigned int) * vectorLen); 
3959 
3960  for(i = 1; i < numTips ; i++)
3961    topBip->bitVector[i / MASK_LENGTH] |= mask32[i % MASK_LENGTH]; 
3962
3963 
3964
3965  /* find the parent of each bip (in the tree they represent) and construct some kind of hashtable this way */
3966#ifdef _USE_PTHREADS
3967
3968  //printf("Parallel region 2\n");
3969
3970  NumberOfJobs = consensusBipLen;
3971  tr->consensusBipLen = consensusBipLen; 
3972  tr->consensusBips = consensusBips;
3973  tr->mxtips = numTips; /* don't need this ? */
3974  tr->hasAncestor = hasAncestor; 
3975  tr->listOfDirectChildren = listOfDirectChildren;
3976  tr->bitVectorLength = vectorLen;   
3977  tr->mutexesForHashing = (pthread_mutex_t**) rax_malloc(consensusBipLen * sizeof(pthread_mutex_t*)); 
3978 
3979  for(i = 0; i < consensusBipLen; i++)
3980    {
3981      tr->mutexesForHashing[i] = (pthread_mutex_t*) rax_malloc(sizeof(pthread_mutex_t));
3982      pthread_mutex_init(tr->mutexesForHashing[i], (pthread_mutexattr_t *)NULL);
3983    }
3984 
3985  masterBarrier(THREAD_PREPARE_BIPS_FOR_PRINT, tr);
3986
3987  /* cleanup */
3988  for(i = 0; i < consensusBipLen; i++)
3989    rax_free(tr->mutexesForHashing[i]);
3990  rax_free(tr->mutexesForHashing);
3991
3992  /* restore the old variables - necessary? */
3993 
3994  hasAncestor = tr->hasAncestor;
3995  listOfDirectChildren = tr->listOfDirectChildren; 
3996#else
3997  {
3998    int 
3999      j,
4000      highestId = 0; 
4001
4002    for(i = 0; i < consensusBipLen; i++)
4003      {
4004        entry
4005          *bipA = consensusBips[i]; 
4006       
4007        /* find first index  */
4008        unsigned int 
4009          firstIndex = 0;
4010       
4011        while(firstIndex < vectorLen && bipA->bitVector[firstIndex] == 0 )
4012          firstIndex++;
4013       
4014        for(j = i + 1; j < consensusBipLen; j++)
4015          {                 
4016            if((unsigned int)consensusBips[i]->amountTips < (unsigned int)consensusBips[j]->amountTips
4017               && issubset(consensusBips[i]->bitVector, consensusBips[j]->bitVector, vectorLen, firstIndex))
4018              {               
4019                List
4020                  *elem = &(elems[highestId]); 
4021               
4022                int 
4023                  *nmbr = &(intList[highestId]); 
4024               
4025                highestId++; 
4026               
4027                elem->value = rax_calloc(1, sizeof(int));
4028               
4029                *nmbr = i; 
4030                elem->value = nmbr; 
4031                elem->next = (listOfDirectChildren[j])
4032                  ?listOfDirectChildren[j]
4033                  :NULL;
4034                listOfDirectChildren[j] = elem;
4035                hasAncestor[i] = TRUE;
4036                break;
4037              }
4038          }
4039      }   
4040  }
4041#endif
4042
4043  /****************************************************************/
4044  /* print the bips during a DFS search on the ancestor hashtable */
4045  /****************************************************************/
4046
4047  /* insert these toplevel bips into the last field of the array */
4048  for(i = 0; i < consensusBipLen; i++)
4049    if( ! hasAncestor[i])
4050      {
4051        List *elem  = (List*) rax_malloc(sizeof(List));
4052        /* elem->value = &i;  */
4053        elem->value = rax_calloc(1, sizeof(int)); /* TODO omg this needs refactoring... */
4054        *(int*)elem->value = i;
4055        elem->next = (listOfDirectChildren[consensusBipLen]) 
4056          ? listOfDirectChildren[consensusBipLen]
4057          : NULL;
4058        listOfDirectChildren[consensusBipLen] = elem;       
4059      }
4060 
4061  /* start dfs search at the top level */
4062  printBipsRecursive(tr, outf, 
4063                     consensusBipLen, consensusBips, 
4064                     numberOfTrees, consensusBipLen,
4065                     listOfDirectChildren, vectorLen, 
4066                     numTips, nameList, 
4067                     topBip, printed, TRUE, printCounter, h, computeIC, &tc, &tcAll, verboseIC);
4068
4069  if(computeIC)
4070    {
4071      double
4072        rtcAll = tcAll / (double)(tr->mxtips - 3),
4073        rtc    = tc    / (double)(tr->mxtips - 3);
4074     
4075      printBothOpen("Tree certainty for this tree: %f\n", tc);
4076     
4077      /* Leonida: for consensus trees I also calculate the relative tree certainty by dividing
4078         the tc by the total number of bipartitions of a tree with n (tr->mxtips) taxa
4079         to penalize the consensi for potentially being unresolved */
4080         
4081      printBothOpen("Relative tree certainty for this tree: %f\n\n", rtc);
4082
4083      printBothOpen("Tree certainty including all conflicting bipartitions (TC-All) for this tree: %f\n", tcAll);
4084      printBothOpen("Relative tree certainty including all conflicting bipartitions (TC-All) for this tree: %f\n\n", rtcAll);
4085    }
4086
4087  rax_free(topBip->bitVector);
4088  rax_free(topBip);
4089  rax_free(printed);
4090  rax_free(hasAncestor);
4091
4092#ifndef _USE_PTHREADS
4093  rax_free(elems);
4094  rax_free(intList);
4095#endif
4096
4097
4098  /* here is a bug, when I try to rax_free the memory on the veryBig (55K)
4099     dataset. When rax_freeing the toplevel bips
4100     (listOfDirectChildren[consensusBipLen]), he complains of sth like
4101     a double rax_free. At this point the value of ptr is not 0, however
4102     the memory cannot be accessed. Also got a "bus error" instead of
4103     the described error here. This is very strange, I already have
4104     accessed the stuff I try to rax_free.   */
4105  /* for(i = 0; i < consensusBipLen+1; i++) */
4106  /*   { */
4107  /*     list *ptr = listOfDirectChildren[i]; */
4108  /*     while(ptr){ */
4109  /*    list *n = ptr->next; */
4110  /*    /\* rax_free(ptr);              /\\* TODO pthreads: last one  *\\/ *\/ */
4111  /*    ptr = n; */
4112  /*     } */
4113  /*   } */
4114}
4115
4116
4117
4118
4119void computeConsensusOnly(tree *tr, char *treeSetFileName, analdef *adef, boolean computeIC)
4120{       
4121  hashtable
4122    *h = initHashTable(tr->mxtips * FC_INIT * 10);
4123
4124  hashNumberType
4125    entries = 0;
4126
4127  unsigned int 
4128    printCounter  = 0,
4129    numberOfTrees = 0, 
4130    i, 
4131    j,     
4132    treeVectorLength, 
4133    vectorLength;
4134
4135  int
4136    consensusBipsLen = 0; 
4137
4138  unsigned int   
4139    **bitVectors = initBitVector(tr, &vectorLength);
4140
4141  entry
4142    **consensusBips;
4143
4144  char 
4145    someChar[1024],
4146    consensusFileName[1024];   
4147 
4148  FILE 
4149    *outf,
4150    *treeFile = getNumberOfTrees(tr, treeSetFileName, adef);
4151
4152  numberOfTrees = tr->numberOfTrees; 
4153
4154  tr->mr_thresh = ((double)numberOfTrees / 2.0);   
4155
4156  assert(sizeof(unsigned char) == 1);
4157 
4158  treeVectorLength = GET_BITVECTOR_LENGTH(numberOfTrees);
4159
4160  /* read the trees and process the bipartitions */ 
4161
4162  for(i = 1; i <= numberOfTrees; i++)
4163    {                 
4164      int 
4165        bCount = 0;
4166     
4167      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);               
4168     
4169      assert(tr->mxtips == tr->ntips);
4170     
4171      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, (i - 1), BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
4172                              &bCount, treeVectorLength, FALSE, FALSE);
4173     
4174      assert(bCount == tr->mxtips - 3);                     
4175    }
4176 
4177  fclose(treeFile);
4178 
4179  if(tr->consensusType == MR_CONSENSUS || tr->consensusType == STRICT_CONSENSUS || tr->consensusType == USER_DEFINED)
4180    {
4181      consensusBips = (entry **)rax_calloc(tr->mxtips - 3, sizeof(entry *));
4182      consensusBipsLen = 0;   
4183    }
4184 
4185  for(j = 0; j < (unsigned int)h->tableSize; j++) /* determine support of the bips */
4186    {           
4187      if(h->table[j] != NULL)
4188        {
4189          entry *e = h->table[j];
4190         
4191          do
4192            {   
4193              unsigned int 
4194                cnt = genericBitCount(e->treeVector, treeVectorLength);
4195
4196                if((tr->consensusType == MR_CONSENSUS     && cnt > (unsigned int)tr->mr_thresh) || 
4197                 (tr->consensusType == STRICT_CONSENSUS && cnt == numberOfTrees) ||
4198                 (tr->consensusType ==  USER_DEFINED    && cnt > (numberOfTrees * tr->consensusUserThreshold) / 100))
4199                {
4200                  consensusBips[consensusBipsLen] = e;
4201                  consensusBipsLen++;
4202                }
4203
4204              e->supportFromTreeset[0] = cnt;
4205              e = e->next;
4206              entries++;
4207            }
4208          while(e != NULL);             
4209        }                       
4210    }   
4211 
4212  assert(h->entryCount == entries);
4213 
4214  if(tr->consensusType == MR_CONSENSUS || tr->consensusType == STRICT_CONSENSUS || tr->consensusType == USER_DEFINED)
4215    assert(consensusBipsLen <= (tr->mxtips - 3));
4216   
4217  if(tr->consensusType == MRE_CONSENSUS)   
4218    mre(h, FALSE, &consensusBips, &consensusBipsLen, 0, tr->mxtips, vectorLength, FALSE , tr, FALSE); 
4219
4220  /* printf("Bips NEW %d\n", consensusBipsLen); */
4221
4222  strcpy(consensusFileName,         workdir); 
4223 
4224  switch(tr->consensusType)
4225    {
4226    case MR_CONSENSUS:
4227      if(computeIC)
4228        strcat(consensusFileName,         "RAxML_MajorityRuleConsensusTree_IC.");
4229      else
4230        strcat(consensusFileName,         "RAxML_MajorityRuleConsensusTree.");
4231      break;
4232    case MRE_CONSENSUS:
4233      if(computeIC)
4234        strcat(consensusFileName,         "RAxML_MajorityRuleExtendedConsensusTree_IC.");
4235      else
4236        strcat(consensusFileName,         "RAxML_MajorityRuleExtendedConsensusTree.");
4237      break;
4238    case STRICT_CONSENSUS:
4239      assert(!computeIC);
4240      strcat(consensusFileName,         "RAxML_StrictConsensusTree.");
4241      break;
4242    case USER_DEFINED :       
4243      if(computeIC)
4244        sprintf(someChar,         "RAxML_Threshold-%d-ConsensusTree_IC.", tr->consensusUserThreshold);
4245      else
4246        sprintf(someChar,         "RAxML_Threshold-%d-ConsensusTree.", tr->consensusUserThreshold);
4247      strcat(consensusFileName, someChar);       
4248      break; 
4249    default:
4250      assert(0);
4251    }
4252 
4253  strcat(consensusFileName,         run_id);
4254
4255  outf = myfopen(consensusFileName, "wb");
4256
4257  fprintf(outf, "(%s,", tr->nameList[1]);
4258 
4259  if(computeIC)
4260    {
4261      if(adef->verboseIC)
4262        printVerboseTaxonNames(tr);
4263      printSortedBips(consensusBips, consensusBipsLen, tr->mxtips, vectorLength, numberOfTrees, outf, tr->nameList, tr, &printCounter, h, computeIC, adef->verboseIC);
4264    }
4265  else
4266    printSortedBips(consensusBips, consensusBipsLen, tr->mxtips, vectorLength, numberOfTrees, outf, tr->nameList, tr, &printCounter, h, computeIC, FALSE);
4267
4268  assert(printCounter ==  (unsigned int)consensusBipsLen);
4269
4270  /* ????? fprintf(outf, ");\n"); */
4271 
4272  fclose(outf);
4273 
4274  if(adef->verboseIC && computeIC)
4275    printBothOpen("Verbose PHYLIP-style formatted bipartition information written to file: %s\n\n",  verboseSplitsFileName);
4276 
4277  switch(tr->consensusType)
4278    {
4279    case MR_CONSENSUS:
4280      if(computeIC)     
4281        printBothOpen("RAxML Majority Rule consensus tree with IC values written to file: %s\n\n", consensusFileName);           
4282      else
4283        printBothOpen("RAxML Majority Rule consensus tree written to file: %s\n", consensusFileName);
4284      break;
4285    case MRE_CONSENSUS:
4286      if(computeIC)
4287        printBothOpen("RAxML extended Majority Rule consensus tree with IC values written to file: %s\n", consensusFileName);
4288      else
4289        printBothOpen("RAxML extended Majority Rule consensus tree written to file: %s\n", consensusFileName);
4290      break;
4291    case STRICT_CONSENSUS:
4292      printBothOpen("RAxML strict consensus tree written to file: %s\n", consensusFileName);
4293      break;
4294    case USER_DEFINED: 
4295      if(computeIC)
4296        printBothOpen("RAxML consensus tree with threshold %d with IC values written to file: %s\n", tr->consensusUserThreshold,  consensusFileName);
4297      else
4298        printBothOpen("RAxML consensus tree with threshold %d written to file: %s\n", tr->consensusUserThreshold,  consensusFileName);
4299      break;
4300    default:
4301      assert(0);
4302    }
4303 
4304  freeBitVectors(bitVectors, 2 * tr->mxtips);
4305  rax_free(bitVectors);
4306  freeHashTable(h);
4307  rax_free(h);
4308  rax_free(consensusBips);
4309
4310  exit(0);   
4311}
4312
4313
4314#else
4315
4316
4317
4318
4319static void mre(hashtable *h, boolean icp, entry*** sbi, int* len, int which, int n, unsigned int vectorLength, boolean sortp, tree *tr, boolean bootStopping)
4320{
4321  entry **sbw;
4322  unsigned int 
4323    i = 0,
4324    j = 0,
4325    k = 0; 
4326 
4327  sbw = (entry **) rax_calloc(h->entryCount, sizeof(entry *)); 
4328
4329  for(i = 0; i < h->tableSize; i++)
4330    {           
4331      if(h->table[i] != NULL)
4332        {
4333          entry *e = h->table[i];
4334          do
4335            {
4336              sbw[j] = e;
4337              j++;
4338              e = e->next;
4339            }
4340          while(e != NULL);
4341        }
4342    }
4343
4344  assert(j == h->entryCount);
4345
4346  if(which == 0)   
4347    qsort(sbw, h->entryCount, sizeof(entry *), _sortByWeight0);     
4348  else   
4349    qsort(sbw, h->entryCount, sizeof(entry *), _sortByWeight1);     
4350
4351  /* ***********************************          */
4352  /* SOS SBI is never rax_freed ********************* */
4353  /* ******************************************** */
4354  /**** this will cause problems for repeated invocations */
4355  /**** with the bootstopping MRE VERSION !!!!!!        ***/
4356     
4357
4358
4359  *sbi = (entry **)rax_calloc(n - 3, sizeof(entry *));
4360
4361  *len = 0;
4362
4363  if(icp == FALSE)
4364    {     
4365      for(i = 0; i < h->entryCount && (*len) < n-3; i++)
4366        {       
4367          boolean compatflag = TRUE;
4368
4369          assert(*len < n-3);             
4370       
4371          /*  for(k = 0; k < (unsigned int)(*len); k++)   */
4372          /*if(sbw[i]->supportFromTreeset[which] <= mr_thresh) */
4373            for(k = ((unsigned int)(*len)); k > 0; k--)
4374              {
4375                /*
4376                  k indexes sbi
4377                  j indexes sbw
4378                  need to compare the two
4379                */
4380               
4381                if(!compatible((*sbi)[k-1], sbw[i], vectorLength))             
4382                  {
4383                    compatflag = FALSE;     
4384                    break;
4385                  }
4386              }
4387         
4388          if(compatflag)
4389            {         
4390              (*sbi)[*len] = sbw[i];
4391              (*len)++;
4392            }         
4393        }
4394    }
4395  else 
4396    {     
4397      for(i = 0; i < (unsigned int)(n-3); i++)
4398        {
4399          (*sbi)[i] = sbw[i];
4400          (*len)++;
4401        }
4402    }
4403
4404  rax_free(sbw);
4405
4406  if (sortp == TRUE)
4407    qsort(*sbi, (*len), sizeof(entry *), sortByIndex);   
4408
4409  return;
4410}
4411
4412
4413
4414
4415
4416
4417
4418static void printBip(entry *curBip, entry **consensusBips, const unsigned int consensusBipLen, const int numtips, const unsigned int vectorLen, 
4419                     boolean *processed, tree *tr, FILE *outf, const int numberOfTrees, boolean topLevel, unsigned int *printCounter)
4420{
4421  int
4422    branchLabel,     
4423    printed = 0;
4424
4425  unsigned int 
4426    i,
4427    j;
4428
4429  unsigned int *subBip = (unsigned int *)rax_calloc(vectorLen, sizeof(unsigned int));
4430
4431  double 
4432    support = 0.0;
4433
4434  for(i = 0; i < consensusBipLen; i++)
4435    {
4436      if((!processed[i]) && issubset(consensusBips[i]->bitVector, curBip->bitVector, vectorLen))
4437        {
4438          boolean processThisRound = TRUE;
4439         
4440          for (j = 0; j < consensusBipLen; j++)     
4441            if(j != i && !processed[j] && issubset(consensusBips[i]->bitVector, consensusBips[j]->bitVector, vectorLen))               
4442              processThisRound = FALSE;                   
4443         
4444          if(processThisRound == TRUE)
4445            {
4446              processed[i] = TRUE;
4447
4448              for(j = 0; j < vectorLen; j++)
4449                subBip[j] |= consensusBips[i]->bitVector[j];
4450             
4451              if(printed == 0 && !topLevel)             
4452                fprintf(outf, "(");             
4453              else             
4454                fprintf(outf, ",");
4455               
4456              printBip(consensusBips[i], consensusBips, consensusBipLen, numtips, vectorLen, processed, tr, outf, numberOfTrees, FALSE, printCounter);
4457
4458              printed += 1;
4459            }
4460        }
4461    }   
4462 
4463  for(i = 0; i < ((unsigned int)numtips); i++)
4464    {
4465      if((((curBip->bitVector[i/MASK_LENGTH] & mask32[i%MASK_LENGTH]) > 0) && ((subBip[i/MASK_LENGTH] & mask32[i%MASK_LENGTH]) == 0) ) == TRUE)
4466        {
4467          if(printed == 0 && !topLevel)     
4468            fprintf(outf,"(");     
4469          else     
4470            fprintf(outf,",");
4471           
4472          fprintf(outf,"%s", tr->nameList[i+1]);
4473          printed += 1;
4474        }
4475    }
4476
4477  rax_free(subBip);
4478
4479  support = ((double)(curBip->supportFromTreeset[0])) / ((double) (numberOfTrees));
4480  branchLabel = (int)(0.5 + support * 100.0);
4481 
4482  if(!topLevel)
4483    {
4484      *printCounter = *printCounter + 1;
4485      fprintf(outf,"):1.0[%d]", branchLabel);
4486    }
4487}
4488
4489void computeConsensusOnly(tree *tr, char *treeSetFileName, analdef *adef)
4490{       
4491  hashtable
4492    *h = initHashTable(tr->mxtips * FC_INIT * 10); 
4493
4494  hashNumberType
4495    entries = 0;
4496 
4497  int 
4498    numberOfTrees = 0, 
4499    i, 
4500    j, 
4501    l,
4502    treeVectorLength,     
4503    consensusBipsLen,
4504    mr_thresh;
4505
4506  unsigned int
4507    printCounter = 0,
4508    vectorLength,
4509    **bitVectors = initBitVector(tr, &vectorLength),
4510    *topBip;
4511
4512  entry
4513    topBipE,
4514    **consensusBips;
4515 
4516  boolean 
4517    *processed;
4518
4519  char 
4520    consensusFileName[1024];
4521 
4522  FILE 
4523    *outf,
4524    *treeFile = getNumberOfTrees(tr, treeSetFileName, adef);
4525
4526     
4527  numberOfTrees = tr->numberOfTrees; 
4528 
4529  mr_thresh = ((double)numberOfTrees / 2.0); 
4530 
4531  assert(sizeof(unsigned char) == 1);
4532 
4533  if(numberOfTrees % MASK_LENGTH == 0)
4534    treeVectorLength = numberOfTrees / MASK_LENGTH;
4535  else
4536    treeVectorLength = 1 + (numberOfTrees / MASK_LENGTH); 
4537
4538  for(i = 1; i <= numberOfTrees; i++)
4539    {                 
4540      int 
4541        bCount = 0;
4542     
4543      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);               
4544     
4545      assert(tr->mxtips == tr->ntips);
4546     
4547      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, (i - 1), BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
4548                              &bCount, treeVectorLength, FALSE, FALSE);
4549     
4550      assert(bCount == tr->mxtips - 3);                     
4551    } 
4552
4553  if(tr->consensusType == MR_CONSENSUS || tr->consensusType == STRICT_CONSENSUS)
4554    {
4555      consensusBips = (entry **)rax_calloc(tr->mxtips - 3, sizeof(entry *));
4556      consensusBipsLen = 0;
4557    }
4558
4559  for(j = 0; j < (int)h->tableSize; j++)
4560    {           
4561      if(h->table[j] != NULL)
4562        {
4563          entry *e = h->table[j];
4564         
4565          do
4566            {   
4567              int cnt = 0;
4568             
4569              unsigned int 
4570                *set = e->treeVector;
4571             
4572              for(l = 0; l < numberOfTrees; l++)                                             
4573                if((set[l / MASK_LENGTH] != 0) && (set[l / MASK_LENGTH] & mask32[l % MASK_LENGTH]))                                         
4574                  cnt++;                                               
4575             
4576              if(tr->consensusType == MR_CONSENSUS)
4577                {
4578                  if(cnt > mr_thresh)                         
4579                    {
4580                      consensusBips[consensusBipsLen] = e;
4581                      consensusBipsLen++;
4582                    }
4583                }
4584
4585              if(tr->consensusType == STRICT_CONSENSUS)
4586                {
4587                  if(cnt == numberOfTrees)
4588                    {
4589                      consensusBips[consensusBipsLen] = e;
4590                      consensusBipsLen++;
4591                    }
4592                }
4593
4594              e->supportFromTreeset[0] = cnt;
4595              e = e->next;
4596              entries++;
4597            }
4598          while(e != NULL);             
4599        }                       
4600    }   
4601
4602  fclose(treeFile); 
4603  assert(entries == h->entryCount);
4604 
4605  if(tr->consensusType == MR_CONSENSUS || tr->consensusType == STRICT_CONSENSUS)
4606    assert(consensusBipsLen <= (tr->mxtips - 3));
4607
4608  if(tr->consensusType == MRE_CONSENSUS)   
4609    mre(h, FALSE, &consensusBips, &consensusBipsLen, 0, tr->mxtips, vectorLength, FALSE, tr);   
4610
4611 
4612  /* printf("Bips OLD %d\n", consensusBipsLen); */
4613
4614  processed = (boolean *) rax_calloc(consensusBipsLen, sizeof(boolean));
4615
4616  topBip = (unsigned int *) rax_calloc(vectorLength, sizeof(unsigned int));
4617 
4618  for(i = 1; i < tr->mxtips; i++)
4619    topBip[i / MASK_LENGTH] |= mask32[i % MASK_LENGTH]; 
4620
4621  topBipE.bitVector = topBip;
4622  topBipE.supportFromTreeset[0] = numberOfTrees;
4623
4624  strcpy(consensusFileName,         workdir); 
4625 
4626  switch(tr->consensusType)
4627    {
4628    case MR_CONSENSUS:
4629      strcat(consensusFileName,         "RAxML_MajorityRuleConsensusTree.");
4630      break;
4631    case MRE_CONSENSUS:
4632      strcat(consensusFileName,         "RAxML_MajorityRuleExtendedConsensusTree.");
4633      break;
4634    case STRICT_CONSENSUS:
4635      strcat(consensusFileName,         "RAxML_StrictConsensusTree.");
4636      break;
4637    default:
4638      assert(0);
4639    }
4640 
4641  strcat(consensusFileName,         run_id);
4642
4643  outf = myfopen(consensusFileName, "wb");
4644
4645  fprintf(outf, "(%s", tr->nameList[1]);
4646  printBip(&topBipE, consensusBips, consensusBipsLen, tr->mxtips, vectorLength, processed, tr, outf, numberOfTrees, TRUE, &printCounter); 
4647  fprintf(outf, ");\n");
4648
4649  assert(consensusBipsLen == (int)printCounter);
4650
4651  fclose(outf);
4652 
4653  switch(tr->consensusType)
4654    {
4655    case MR_CONSENSUS:
4656      printBothOpen("RAxML Majority Rule consensus tree written to file: %s\n", consensusFileName);
4657      break;
4658    case MRE_CONSENSUS:
4659      printBothOpen("RAxML extended Majority Rule consensus tree written to file: %s\n", consensusFileName);
4660      break;
4661    case STRICT_CONSENSUS:
4662      printBothOpen("RAxML strict consensus tree written to file: %s\n", consensusFileName);
4663      break;
4664    default:
4665      assert(0);
4666    }
4667 
4668 
4669  rax_free(topBip);
4670  rax_free(processed);
4671
4672  freeBitVectors(bitVectors, 2 * tr->mxtips);
4673  rax_free(bitVectors);
4674  freeHashTable(h);
4675  rax_free(h);
4676  rax_free(consensusBips); 
4677 
4678  exit(0);   
4679}
4680
4681#endif
Note: See TracBrowser for help on using the repository browser.